section[1.2em]1.2em6pt subsection[5em]2em6pt cr usr

IFT-UAM/CSIC-24-75 DESY-24-071

Non-Gaussian tails without stochastic inflation

Guillermo<EMAIL_ADDRESS>Thomas
<EMAIL_ADDRESS>
Alejandro Pérez<EMAIL_ADDRESS>Mathias
<EMAIL_ADDRESS>and Julián<EMAIL_ADDRESS>
1Departamento de Física Teórica, Universidad Autónoma de Madrid (UAM),
Campus de Cantoblanco, 28049 Madrid, Spain
2Instituto de Física Teórica (IFT) UAM-CSIC, Campus de Cantoblanco, 28049
Madrid, Spain
3Deutsches Elektronen-Synchrotron DESY, Notkestr. 85, 22607 Hamburg, Germany

Abstract

We show, both analytically and numerically, that non-Gaussian tails in the
probability density function of curvature perturbations arise in ultra-slow-
roll inflation from the $\delta N$ formalism, without invoking stochastic
inflation. Previously reported discrepancies between both approaches are a
consequence of not correctly accounting for momentum perturbations. Once they
are taken into account, both approaches agree to an excellent degree. The
shape of the tail depends strongly on the phase space of inflation.

###### Contents

  1. 1 Introduction
  2. 2 Inflationary setup
  3. 3 The $\delta N$ formalism
    1. 3.1 Approximate solutions in classical $\delta N$
    2. 3.2 Full calculation
    3. 3.3 Comparison between the full and the approximate implementations of $\delta N$
  4. 4 Stochastic inflation
  5. 5 Summary and conclusions
  6. A The $\delta N$ formula

##  1 Introduction

There is a broad consensus that current observational bounds allow primordial
black holes (PBHs) with masses between $10^{-16}\,M_{\odot}$ and
$10^{-12}\,M_{\odot}$ to account for all the dark matter of the Universe [1].
The most popular idea for the formation of PBHs that could constitute the dark
matter assumes the existence of very large density fluctuations in the early
Universe (before primordial nucleosynthesis), which collapse under their own
gravitational pull. Indeed, if the collapse is assumed to start during
radiation domination, each fluctuation needs to overcome a threshold for the
relative overdensity of nearly $\mathcal{O}(1)$ in order to become a black
hole. The statistical properties of these fluctuations are the main factor
determining the current PBH abundance and mass distribution. Thinking in terms
of Fourier modes to describe the early Universe inhomogeneities, the collapse
towards a black hole starts to occur as soon as the wavelength of a
fluctuation above the threshold is comparable to the Hubble radius. This
wavelength is the characteristic scale that determines the mass of the
resulting black holes.

The possible origin of these large density fluctuations is uncertain. A
popular way to realize the idea is to play with the dynamics of inflation in
such a way that curvature fluctuations acquire a large variance at adequate
distance scales. Using linear cosmological perturbation theory to relate
density and curvature fluctuations and assuming Gaussianity of the latter, a
simple estimate (see e.g. [2]) shows that a power spectrum of curvature
perturbations $\sim 10^{-2}$ –i.e. about seven orders of magnitude larger than
fluctuations at the scales we are able to probe with the Cosmic Microwave
Background (CMB)– is needed to account for all dark matter. One way to attain
this consists in introducing a period of ultra-slow roll inflation (USR) [3],
during which Fourier modes of curvature fluctuations experience super-Hubble
growth [4]. Concrete models that implement this mechanism are e.g. the ones
proposed in [2] and [5].

It is known that in USR inflation, Fourier modes with small amplitude acquire
some amount of local non-Gaussianity [6]. This already indicates that the
previous Gaussian estimate for the abundance of PBHs in such models is too
simplistic. In addition, several recent studies have highlighted the
appearance in models of USR of non-Gaussian (exponential) tails in the
probability distribution function (PDF) of $\zeta$ for large values [7, 8, 9].
These results have been obtained using mainly two different (but related)
frameworks that go beyond standard cosmological perturbation theory and
several computational techniques: the $\delta N$ formalism and the stochastic
$\delta N$ formalism.111In what follows we will often use the term classical
$\delta N$ formalism to refer to the non-stochastic version of it. The term
classical is used here in contrast to stochastic, and in no way implies that
the fluctuations in the $\delta N$ formalism are not quantum in nature. In
this work we aim to shed further light on the appearance of these tails.
Specifically, we are interested in understanding under which conditions (if
any) the results obtained with these two frameworks differ.

In essence, the $\delta N$ formalism exploits a non-linear relation [10]
between the inflaton fluctuations, $\delta\phi$, and the curvature
perturbation on uniform density slices $\zeta$.222See Appendix A and ref. [11]
for the definition of $\zeta$. Once this non-linear relation is established
and if the PDF of $\delta\phi$ is known, the PDF of $\zeta$ can be easily
obtained. If $\delta\phi$ is assumed to be a Gaussian variable, its non-linear
relation to $\zeta$ makes the PDF of the latter non-Gaussian. In particular,
in the case of USR, an exponentially decaying tail arises for large $\zeta$.
Therefore, within this formalism, the appearance of a non-Gaussian tail in the
PDF of $\zeta$ is just the consequence of a change of variables. The
stochastic $\delta N$ formalism adds another twist to this argument by
incorporating the framework of stochastic inflation [12]. This framework is
used to study the effect of short-wavelength modes on the dynamics of long-
wavelength modes. This backreaction is included via stochastic noise terms
acting on the dynamical equations of the latter. Within this framework (whose
full non-linear implementation is highly costly in computational terms [8]),
the appearance of exponential tails has been presented as a consequence of
quantum diffusion, see e.g. [7]. A natural question to ask is how different
are the tails obtained in these two formalisms [13].

In our work we use two different implementations (models) of USR to carefully
analyze the relevance of the phase space of inflation in order to understand
the practical differences between the two formalisms. In the process of doing
so, we discuss the validity and convenience of using a specific approximation
(the unperturbed-trajectory approximation) that simplifies the treatment of
the fluctuations of the inflaton assuming they follow, at all times, the
unperturbed solution of the background equations of motion.

In the next section we introduce our notation and the two models we use. One
of them (Model 1) is a piecewise implementation of the USR phase and the other
(Model 2), first presented in [2], involves a polynomial potential for the
inflaton. In Section 3 we discuss the application of the classical $\delta N$
formalism to the calculation of the PDF of $\zeta$ assuming Gaussian inflaton
fluctuations, first by using the unperturbed-trajectory approximation and then
by performing the full calculation in both models. We find that the
unperturbed-trajectory approximation works well only for Model 2 and it fails
for Model 1 due to the more complex shape of its phase space. In Section 4 we
discuss the stochastic $\delta N$ formalism and compute the PDF of $\zeta$ for
both models. An analogous approximation to the one we explored for the
classical $\delta N$ formalism can also be used now, but this time it relates
the stochastic terms. Again, we will find a qualitative difference among the
behaviours of both models, which can be understood in terms of a phase space
analysis. We present our conclusions in Section 5. An appendix is devoted to a
more detailed presentation of the $\delta N$ formalism which includes several
technical aspects that we do not discuss in the main body of the text. We use
natural units and set the reduced Planck mass to 1.

##  2 Inflationary setup

Our inflationary setup consists of a single scalar field $\phi$ with an action
given by

$S=\int\mathop{}\\!\mathrm{d}^{4}x\sqrt{-g}\bigg{[}\frac{1}{2}R+\frac{1}{2}(\partial_{\mu}\phi)^{2}-V(\phi)\bigg{]},$
(2.1)

where $g$ denotes the determinant of the spacetime metric, $R$ the Ricci
scalar, and $V$ the potential. We have also set the (reduced) Planck mass to
1. We work with the following sign convention for the FLRW metric

$\mathop{}\\!\mathrm{d}s^{2}=-\mathop{}\\!\mathrm{d}t^{2}+a(t)^{2}\mathop{}\\!\mathrm{d}{\bm{x}}^{2},$
(2.2)

where $a(t)$ denotes the scale factor. With this metric, the equations of
motion for the (spatially homogeneous) field, $\phi$, and its (also
homogeneous) conjugate momentum $\pi$ are

$\phi^{\prime}=\pi\qquad{\rm
and}\qquad\pi^{\prime}=-(3-\epsilon)(\pi+\partial_{\phi}\log V)\,,$ (2.3)

where primes denote derivatives with respect to the number of $e$-folds

$\displaystyle\mathop{}\\!\mathrm{d}N=H\mathop{}\\!\mathrm{d}t\,,$ (2.4)

and we use the following definition for the slow-roll parameters,

$\epsilon=-\frac{H^{\prime}}{H}=\frac{1}{2}\pi^{2}\qquad{\rm
and}\qquad\eta=\epsilon-\frac{\epsilon^{\prime}}{2\epsilon}=-\frac{\ddot{\phi}}{{\dot{\phi}}H}\,,$
(2.5)

where dots indicate derivatives with respect to cosmic time, $t$. In
principle, (2.3) is just a set of equations with no specified initial
conditions. A specific solution to these equations is what we will call the
_unperturbed trajectory_ , denoted by $\bar{\phi}$ and $\bar{\pi}$. By this,
we refer to a solution of (2.3) with initial conditions set in the far past
(i.e. for $\phi$ larger than those corresponding to the USR part of the
potential) and that has reached the SR attractor before entering the USR
phase.333Since the SR part of the potential indeed corresponds to an attractor
in phase space, there is a wide range of initial conditions that, if set
sufficiently far in the past, satisfy this condition. This corresponds to what
is usually called the _background solution_ in standard perturbation theory,
and it is depicted as a black dashed line in the phase spaces in Fig. 3.

Given explicit expressions for the slow-roll parameter $\eta$ as a function of
time and an initial condition for $\epsilon$, the potential can be obtained as
an integral over the number of $e$-folds as

$V(N)=V(N_{\star})\exp\bigg{[}2\int_{N_{\star}}^{N}\,\epsilon(\tilde{N})\,\frac{\eta(\tilde{N})-3}{3-\epsilon(\tilde{N})}\mathop{}\\!\mathrm{d}\tilde{N}\bigg{]}\,.$
(2.6)

In this expression $N_{\star}$ is just some reference $e$-fold value.

We are interested in models of inflation characterized by a potential, $V$,
featuring a non-stationary ($\partial_{\phi}V\neq 0$) inflection point so that
the inflaton experiences an ultra-slow-roll (USR) phase, defined by
$\eta\gtrsim 3/2$,444USR is often defined in the literature as the condition
$\ddot{\phi}+3H\dot{\phi}\simeq 0$. This corresponds to $\eta\simeq 3/2$ (and
$\epsilon\ll 1$), which implies $V^{\prime}\ll\ddot{\phi}$. In this work we
adopt a broader definition, using the term USR to refer to any phase with
$\eta\gtrsim 3/2$, since for these values of $\eta$ the curvature perturbation
grows exponentially, see Eq. (2.13). followed by a constant-roll (CR) phase,
defined by $(\epsilon^{\prime}/\epsilon)^{\prime}\simeq 0$ and negative
$\eta$. Before reaching the inflection point, the inflaton undergoes a
conventional slow-roll (SR) phase, where both $\epsilon$ and $\eta$ are small.
As the field approaches the local maximum of the potential, it decelerates,
and $\eta$ becomes large. Indeed, approximating the potential in a
neighbourhood of its maximum, $\phi_{0}$, by a parabola and assuming the field
decreases as inflation proceeds, if we neglect $\epsilon$ in Eq. (2.3) and if
$\phi^{\prime}_{0}\equiv\phi^{\prime}(N_{0})\neq 0$ at $\phi_{0}$, then

$\displaystyle\phi\simeq\phi_{0}+\frac{\phi^{\prime}_{0}}{\sqrt{9-12\eta_{V}}}\left[-e^{-\eta_{+}(N-N_{0})\,}+e^{-\eta_{-}(N-N_{0})}\right]\,,$
(2.7)

where $\eta_{V}=\partial_{\phi}^{\,2}V/V<0$ is a constant to a good degree of
approximation and, in the limit $\phi\gg\phi_{0}$ we have that
$\eta\simeq\eta_{+}\equiv\frac{3}{2}\left(1+\sqrt{1-4\eta_{V}/3}\right)$,
whereas for $\phi\ll\phi_{0}$,
$\eta\simeq\eta_{-}\equiv\frac{3}{2}\left(1-\sqrt{1-4\,\eta_{V}/3}\right)$.
Therefore, the inflaton climbs up the potential in USR (with $\eta>0$), and
afterwards, it reaches a phase of CR (with $\eta<0$). As noted in [14], the
values of $\eta$ in the two phases are related by the Wands duality [15], such
that their sum is $3$.

Figure 1:  Left panel: Slow-roll parameters $\epsilon$ (blue) and $|\eta|$
(red) for Model 1 (dashed lines) and Model 2 (solid lines). Right panel:
Scalar power spectra corresponding to the two models in the left panel.
$N_{*}$ is the time at which $\eta=3/2$ in the USR-CR transition along the
unperturbed trajectory, see the comment below (2.13).

In what follows, we will consider two different ways of implementing a phase
of USR inflation. We refer to them as Model 1 and Model 2. Their slow-roll
parameters are shown in the left panel of Fig. 1, and the corresponding linear
power spectra of curvature fluctuations, in the right panel.

$\star$ Model 1. The first model consists of a phenomenological
parameterization of $\eta$ as a function of the number of $e$-folds:555Similar
parameterizations have been been used before, see e.g. [6]. The difference
between ours and the one in [6] is on the phase preceding USR. In our case
$\eta$ gradually drops towards zero in the limit $N\rightarrow-\infty$,
whereas the parameterization of [6] drops to zero more abruptly.

$\eta(N)=\frac{\eta_{\rm
USR}}{4}\Big{[}1+\tanh(N-{N}_{1})\Big{]}\Big{[}1+\tanh({N}_{2}-N)\Big{]}+\frac{\eta_{\rm
CR}}{2}\Big{[}1+\tanh{(N-{N}_{2})}\Big{]}.$ (2.8)

For our numerical example we choose $\eta_{\rm USR}=4$, $\eta_{\rm CR}=-1$ and
${N}_{2}-{N}_{1}=2.3$. The USR phase (as defined above, with $\eta>3/2$) sets
in at $N-N_{2}\simeq-2.54$ and ends at $N\simeq N_{2}$. In between these two
instants, there is an interval (from $N-N_{1}\simeq-1.64$ to
$N-N_{2}\simeq-0.79$ ) during which $\eta$ is just above $3$ (and nearly
constant). The subsequent CR phase begins approximately at $N\simeq{N}_{2}+3$.
In the limit $N\rightarrow-\infty$, $\eta$ tends rapidly to zero, as shown in
Fig. 1. The definition of $\eta$ in Eq. (2.5) can be integrated to obtain
$\epsilon$, up to an initial condition which we choose in such a way that
$\epsilon=10^{-5}$ in the far past. The potential can be obtained via Eq.
(2.6), with $V(N_{\star})$ chosen in such a way that the power spectrum is
correctly normalized on CMB scales (again, in the far past). We take
$\mathcal{P}_{\mathcal{\zeta}}(k_{\rm CMB})=2\times 10^{-9}$ at $k_{\rm
CMB}=6.3\times 10^{-17}k_{*}$, where $k_{*}=a(N_{*})H(N_{*})$ and $N_{*}$
denotes the time at which $\eta\simeq 3/2$ after USR, see the comment below
(2.13).

$\star$ Model 2. The second model we study is the one presented in [2, 5], in
which a scalar field $\chi$ with a fourth order polynomial potential couples
to the Ricci scalar $R$ via a term $\xi\sqrt{-g}\chi^{2}R$ in the Lagrangian.
In the Einstein frame

$U(\chi)=\frac{\lambda\chi^{4}}{4!(1+\xi\chi^{2})^{2}}\bigg{[}3+\xi^{2}\chi_{0}^{4}-8(1+c_{3})\frac{\chi_{0}}{\chi}+2(1+c_{2})(3+\xi\chi_{0}^{2})\frac{\chi_{0}^{2}}{\chi^{2}}\bigg{]},$
(2.9)

and $\chi$ has a non-canonical kinetic term. The coefficients of the quartic
polynomial have been rewritten in a convenient way so that the potential
features an inflection point at $\chi_{0}$ for $c_{2}=c_{3}=0$. The relation

$\phi=\sqrt{\frac{1+6\xi}{\xi}}\sinh^{-1}\Big{[}\chi\sqrt{\xi(1+6\xi)}\Big{]}-\sqrt{6}\tanh^{-1}\bigg{[}\frac{\sqrt{6}\xi\chi}{\sqrt{1+\xi(1+6\xi)\chi^{2}}}\bigg{]}\,,$
(2.10)

allows to work with the field $\phi$, which is canonically normalized. We
choose $c_{2}=0.011$, $c_{3}=0.0089$, $\chi_{0}=1$, $\xi=0.325479$ and
$\lambda$ so that the power spectrum is correctly normalized on CMB scales
[5]. Again, we take $\mathcal{P}_{\mathcal{\zeta}}(k_{\rm CMB})=2.1\times
10^{-9}$ at $k_{\rm CMB}=6.3\times 10^{-17}k_{*}$.

Physical properties of the universe such as the CMB temperature fluctuations
are related to inflation via the (linear) comoving curvature perturbation
$\mathcal{R}$,666In flat gauge (see Appendix A), the linear comoving curvature
perturbation reduces to $\mathcal{R}=-\delta\phi/\phi^{\prime}$, where
$\delta\phi$ is the inflaton fluctuation in that gauge. At linear order in
perturbation theory and for single-field inflation, $\mathcal{R}$ and $\zeta$
(which is the variable used in the bulk of this work) agree for frozen super-
Hubble Fourier modes.  see e.g. [16, 17], which is gauge invariant at linear
order and whose linearized Fourier modes obey the equation

$\mathcal{R}_{k}^{\prime\prime}+(3-\epsilon-2\eta)\mathcal{R}_{k}^{\prime}+\frac{k^{2}}{a^{2}H^{2}}\mathcal{R}_{k}=0\,.$
(2.11)

In this equation, $H=\dot{a}/a$ and $k$ denotes the comoving wavenumber of the
Fourier mode $\mathcal{R}_{k}$. At early times during inflation the modes are
assumed to be in the Bunch-Davies vacuum:

$\mathcal{R}_{k}=-\frac{1}{a\phi^{\prime}\sqrt{2k}}\exp\bigg{(}-ik\int\frac{\mathop{}\\!\mathrm{d}N}{aH}\bigg{)}.$
(2.12)

Let us assume that $\epsilon\ll 1$ (so $H^{\prime}\simeq 0$) and $\eta$ is
constant, which holds both during USR and CR. Then, for $k\ll a\,H$,

$\mathcal{R}_{k}(N)\simeq\mathcal{R}_{k}(N_{c})+\mathcal{R}_{k}^{\prime}(N_{c})\int_{N_{c}}^{N}\exp\bigg{\\{}\int_{N_{c}}^{\hat{N}}\Big{[}2\eta(\tilde{N})-3\Big{]}\mathop{}\\!\mathrm{d}\tilde{N}\bigg{\\}}\mathop{}\\!\mathrm{d}\hat{N},$
(2.13)

where $N_{c}$ is some reference time which we may choose to be given by
$k=a(N_{c})H(N_{c})$. In the SR limit, $\eta\simeq 0$ and the second term of
Eq. (2.13) decays exponentially, so
$\mathcal{R}_{k}(N)\simeq\mathcal{R}_{k}(N_{c})$. However, during the USR
phase this solution grows, and the power spectrum is enhanced. During the
subsequent CR phase with $\eta\lesssim 0$ the solution decays again, and the
modes of $\mathcal{R}$ that contribute to the peak of the power spectrum
freeze to their final values roughly at the end of the USR phase, see e.g.
[18]. From (2.11) (also from (2.13)), one can see that the monotony of
$\mathcal{R}$ changes at $\eta=(3-\epsilon)/2\approx 3/2$. The value
$\eta=3/2$ is achieved twice along the unperturbed trajectory: once in the SR-
USR transition, and another time in the USR-CR transition. We denote as
$N_{*}$ the time at which $\eta=3/2$ along the unperturbed trajectory in the
USR-CR transition.

##  3 The $\delta N$ formalism

The $\delta N$ formalism [19] allows to compute the curvature perturbation
$\zeta$ non-linearly on large distances (compared to the Hubble scale) as a
function of the flat-gauge, linear perturbations $\delta\phi$ and $\delta\pi$,
using only the unperturbed equations of motion (2.3). This provides an
explicit relation between the non-linear $\zeta$ and the linear $\delta\phi$
and $\delta\pi$, which can be exploited to compute the probability
distribution function (PDF) of the former provided one knows the PDF of the
latter. We will now summarize the essence of the $\delta N$ formalism and we
refer the reader to Appendix A (and references therein) for technical details.

The $\delta N$ formalism crucially relies on the validity of a _gradient
expansion_ , i.e. an expansion of the non-linear equations of motion in powers
of gradients (in position space), or powers of $k/(aH)$ (in Fourier space).
First, using such an expansion, one can write down a _$\delta N$ formula_
relating $\zeta$ (on large scales) to the local expansion between two
conveniently-chosen hypersurfaces (which are appropriately determined by the
values of $\delta\phi$ and $\delta\pi$ on them). Second, from the gradient
expansion follows the so-called _separate universe approach_ (see [20]), which
states that non-linearly perturbed fields, coarse-grained on large scales, are
described by differential equations that look like those governing the
unperturbed (homogeneous) fields. Invoking the separate universe approach, one
can then relate the aforementioned local expansion to the number of $e$-folds
computed using the unperturbed equations of motion.

In practice, the prescription is the following. We specify an initial
hypersurface, determined by an $e$-fold time $N_{i}$ and unperturbed field and
momentum values $\bar{\phi}(N_{i})$, $\bar{\pi}(N_{i})$.777We recall that
$\bar{\phi}$ and $\bar{\pi}$ are solutions to (2.3) with initial conditions
fixed sufficiently before the USR phase, such that they have converged to the
SR attractor before entering the USR phase of their dynamics (see comment
below (2.5)). Ideally, in a calculation without approximations, the specific
choice of $N_{i}$ would be immaterial. However, in a practical application of
the $\delta N$ formalism such as the one we implement, $N_{i}$ should be
chosen late enough such that the linear Fourier modes of $\zeta$ (or
equivalently $\mathcal{R}$) have frozen and gradients become negligible [21].
At the same time, $\delta\phi$ behaves linearly only at early times when the
kinetic terms dominate the dynamics, so choosing $N_{i}$ too late would
neglect non-linearities and therefore using the linear approximation to
compute the variance of $\zeta$ would become less accurate. As a consequence
of the impossibility of optimizing both aspects simultaneously, a spurious
dependence on $N_{i}$ remains. As we shall discuss in this section, the
details of the phase space of the inflationary model determine how strongly
the PDF of $\delta N$ actually ends up depending on the choice of $N_{i}$. In
Sec. 4, we will discuss the relation between the choice of $N_{i}$ and that of
the coarse-graining parameter $\sigma$ in stochastic inflation. For a strongly
peaked power spectra like the ones arising in the presence of an USR phase,
the time of freezeout of the Fourier modes corresponding to the peak of the
power spectrum provides a choice for $N_{i}$ yielding consistent results with
the stochastic $\delta N$ formalism. Under certain circumstances, the
dependence on $N_{i}$ of the PDF of $\zeta$ may vanish completely, as we will
discuss later.

Second, we perturb $\bar{\phi}(N_{i})$, $\bar{\pi}(N_{i})$ by
$\delta\phi(N_{i})$ and $\delta\pi(N_{i})$ (henceforth $\delta\phi$ and
$\delta\pi$). Third, we compute the difference in the number of $e$-folds
$\delta N$ needed to reach some later stage of inflation given by a fixed
field value $\bar{\phi}(N_{f})$ (e.g. the end of inflation) with respect to
the unperturbed case by evolving the _perturbed initial conditions_
$\bar{\phi}(N_{i})+\delta\phi$, $\bar{\pi}(N_{i})+\delta\pi$ with the
_unperturbed equations of motion_ (2.3). By virtue of the $\delta N$ formula,
this difference $\delta N$ coincides with the (frozen) non-linear perturbation
$\zeta$. We thus reach a non-linear relation between $\delta N=\zeta$ and
$\delta\phi$, $\delta\pi$. Formally:

$\displaystyle\zeta(N_{f})=\delta N(N_{f})=$
$\displaystyle\,\,N\left[\\{\bar{\phi}(N_{i})+\delta\phi,\bar{\pi}(N_{i})+\delta\pi\\}\to\bar{\phi}(N_{f})\right]$
$\displaystyle-N\left[\\{\bar{\phi}(N_{i}),\bar{\pi}(N_{i})\\}\to\bar{\phi}(N_{f})\right]\,,$
(3.1)

where $N[\\{\phi_{1}\,,\pi_{1}\\}\to\phi_{2}]$ denotes the number of $e$-folds
that it takes to reach the field value $\phi_{2}$ evolving the unperturbed
equations of motion (2.3) with the initial conditions
$\\{\phi_{1}\,,\pi_{1}\\}$. The dependence on $N_{f}$ in (3.1) vanishes
because of the conservation of $\zeta$ for $N>N_{i}$ [21]. The residual
dependence on $N_{i}$ will be discussed in the rest of this section.

As explained above, the $\delta N$ formalism provides a non-linear relation
between $\delta N$ and the field perturbations $\delta\phi$, $\delta\pi$ at
some $e$-fold time $N_{i}$, which can then be used to obtain the PDF of the
non-linear perturbation $\zeta=\delta N$. Let us emphasize that the dynamics
or statistics of $\delta\phi$ and $\delta\pi$ are not determined in this
formalism through any action or equation of motion; in this context,
$\delta\phi$ and $\delta\pi$ are just a set of variables whose functional
relation to $\delta N$ we need to know in order to infer the statistical
properties of the latter by sampling them. In general, $\delta N$ will depend
non-linearly on both $\delta\phi$ and $\delta\pi$. Let us now consider for
simplicity that $\delta\phi$ and $\delta\pi$ are not independent of each other
(an assumption whose validity we will discuss later on). In this case, and if
we know the PDF of $\delta\phi$, we can obtain the PDF of $\zeta$ via:

$P(\zeta)=P(\delta\phi)\bigg{|}\frac{\mathop{}\\!\mathrm{d}\delta\phi}{\mathop{}\\!\mathrm{d}\delta
N}\bigg{|}\,.$ (3.2)

It is however important to stress that the $\delta N$ formalism does not
provide any information a priori about the statistical properties of
$\delta\phi$, which is an input that needs to be introduced by hand. If, as it
is often done in the literature [22, 23, 9, 13, 24], the probability density
function (PDF) of $\delta\phi$ is assumed to be Gaussian

$P({\delta\phi})=\frac{1}{\sigma_{\delta\phi}\sqrt{2\pi}}\exp\bigg{(}-\frac{\delta\phi^{2}}{2\sigma_{\delta\phi}^{2}}\bigg{)}\,,$
(3.3)

with variance $\sigma_{\delta\phi}^{2}$, the non-linear relation between
$\delta N$ and $\delta\phi$ given by the $\delta N$ formalism results in a
non-Gaussian PDF for $\zeta$ through Eq. (3.2). We will assume Gaussianity of
$\delta\phi$ in this work for the sake of simplicity, but we emphasize that
obtaining the correct PDF of $\zeta$ requires, in principle, full knowledge of
the statistics of inflation fluctuations. Determining if $\delta\phi$ is
Gaussian away from the maximum of the PDF is, in general, not a simple
task.888For instance, the bispectrum in a model similar to our Model 1 was
discussed in [6] using the in-in formalism. For smooth transitions, the final
bispectrum of $\mathcal{R}$ in that case is of the local form and its
magnitude is controlled by the value of $\eta$ in the CR phase. A consistent
result can be obtained by expanding the $\delta N$ result to second order in
perturbations, see [9]. The use of perturbation theory to compute correlators
around the maximum of the PDF does not seem a practical strategy to find the
statistical properties of large fluctuations because an extrapolation is in
general not well justified. For this reason, some non-perturbative approaches
have recently started to be considered, see in particular [25].

###  3.1 Approximate solutions in classical $\delta N$

In this section, we discuss two different (but related) approximate
implementations of the $\delta N$ formalism to compute $\delta
N(\delta\phi,\delta\pi)$. Both approaches relate $\delta\phi$ and $\delta N$
according to

$\delta\phi=\bar{\phi}(N_{i}-\delta N)-\bar{\phi}(N_{i})\,,$ (3.4)

where $\bar{\phi}$ evolves according to (2.3) and the initial time $N_{i}$ is
chosen after the linear $\zeta$, or equivalently $\mathcal{R}$, has become
constant, as mentioned earlier. The importance of the non-linear relation
between $\delta\phi$ and $\zeta$ is manifest in this expression. As we shall
see in Eq. (3.6), solving Eq. (3.4) for constant $\eta$ and Taylor-expanding
the result around $\delta N=0$, we obtain

$\delta\phi=-\pi\,\delta N+\mathcal{O}(\delta N^{2})\,,$ (3.5)

where the equality has to be evaluated at a time $N_{i}$ after (the linear)
$\mathcal{R}$ has become constant. The leading order in (3.5) is nothing but
the standard linear relation between $\zeta$ and $\delta\phi$. The non-linear
terms $\mathcal{O}(\delta N^{2})$ are responsible for the non-Gaussian
behavior in $\zeta$ for large values of the fluctuations, even if $\delta\phi$
is assumed to be fully Gaussian, as discussed above.

Notice that the formula (3.4) implicitly fixes $\delta\pi$ as a function of
$\delta\phi$ in such a way that $\bar{\pi}(N_{i})+\delta\pi$ lies on the
unperturbed trajectory of the inflaton in phase space. There is, in principle,
no reason why this (implicit) choice of $\delta\pi$ should be the right one.
The prescription has nevertheless been used in the literature, see e.g. [26,
24], to obtain analytical results.

####  3.1.1 Eternal CR approximation

Let us assume that the field $\phi$ moves along its unperturbed trajectory in
phase space, according to Eq. (3.4). In addition, we assume that this
trajectory describes a CR phase. We will call this the eternal CR
approximation. As we will see, this approximation always works well if very
small fluctuations are considered. We will also see (later in Section 3.2)
that if this approximation is used for larger fluctuations it can work well
(or not) depending on the shape of the phase space. Physically, this
corresponds to an inflaton released close to the top of a quadratic potential
in the distant past.

Along a CR trajectory, $\bar{\phi}(N)\sim e^{-\eta N}$ (where $\eta$ is
constant) and the homogeneous perturbations $\delta\phi$ (along the
unperturbed trajectory) obey [13]

$\delta\phi=\left(e^{\eta\,\delta N}-1\right){\bar{\phi}(N_{i})}\,.$ (3.6)

Denoting the linear power spectrum of $\zeta$ as

$\mathcal{P}_{\mathcal{\zeta}}(k)\equiv\frac{k^{3}}{2\pi^{2}}|\mathcal{\zeta}_{k}|^{2}\,,$
(3.7)

evaluated once the modes are super-Hubble and have become constant, and using
that999 We emphasize that $\mathcal{P}_{\mathcal{\zeta}}$ is the power
spectrum of $\zeta$ computed _in linear perturbation theory_ , which coincides
with the linear $\mathcal{P}_{\mathcal{R}}$, and therefore can be computed by
solving (2.11). Correspondingly, $\sigma^{2}_{\zeta}$ is the variance of the
linear (Gaussian) $\zeta$, which does not coincide (in general) with the
variance of the non-linear $\zeta$ (whose full PDF we intend to compute).

$\sigma_{\delta\phi}^{2}=\int\mathcal{P}_{\delta\phi}(k)\,\mathop{}\\!\mathrm{d}\log
k=\bar{\pi}(N_{i})^{2}\int\mathcal{P}_{\zeta}(k)\mathop{}\\!\mathrm{d}\log
k\equiv\bar{\pi}(N_{i})^{2}\sigma_{\zeta}^{2}\,,$ (3.8)

the PDF of $\zeta=\delta N$ turns out to be [13, 9]

$P(\zeta)=\frac{1}{\sqrt{2\pi\,\sigma^{2}_{\zeta}}}\exp\left[-\frac{1}{2\sigma^{2}_{\zeta}\,\eta^{2}}\left(1-e^{\eta\,\delta
N}\right)^{2}+\eta\,\delta N\right]\,.$ (3.9)

In the limit of large $\zeta$, its PDF scales as

$\displaystyle P(\zeta)\sim e^{\eta\zeta}.$ (3.10)

We point out that the dependence on $N_{i}$ of the PDF of $\delta N$ that
might have come from $\bar{\pi}(N_{i})^{2}$ in (3.8) actually vanishes
identically. However, there is in principle a residual dependence on $N_{i}$
that is not explicit in the above expressions. It comes from the definition of
the linear $\sigma^{2}_{\mathcal{\zeta}}$; in particular from the upper limit
of integration.

####  3.1.2 Unperturbed-trajectory approximation

Let us now consider a model of USR inflation where the subsequent CR phase has
a finite duration in the past. Thus, a sufficiently large $\delta\phi$
perturbation makes $\bar{\phi}+\delta\phi$ fall outside of that CR phase and
into the USR part of the unperturbed trajectory (with $\eta$ changing from
negative to positive).101010This approximation and the eternal CR one
discussed above coincide form sufficiently small fluctuations, since we are
are considering models in which USR is followed by a constant CR phase of
finite duration. Again, the solution for $\delta N(\delta\phi)$ can be
obtained by numerically solving (3.4). This _unperturbed-trajectory
approximation_ then leads to a bump in the PDF of $\zeta=\delta N$ located at
a large enough value of $\delta N$ and superimposed on the behaviour described
by Eq. (3.10). The emergence of such a (spurious) bump can be understood from
(3.2). As seen there, the PDF of $\delta N$, $P(\delta N)$ is the product of
$P(\delta\phi)$ and the absolute value of the Jacobian of the change of
variables. While $\bar{\phi}+\delta\phi$ falls inside the CR region (where
$\eta<0$), the latter provides the dominant contribution to $P(\delta N)$ for
large $\delta N$ (cf. Eq. (3.10)), which corresponds to an exponential decay.
Now, if $\delta\phi$ is such that $\bar{\phi}+\delta\phi$ falls into the USR
part of the unperturbed trajectory (as it will eventually happen in the
unperturbed-trajectory approximation), the same functional form for
$P(\delta\phi)$ and the Jacobian still holds, with the sole addition of some
constant factors that ensure the continuity of $\delta N$ across the CR-USR
transition (that in any case do not affect the behaviour of $P(\delta\phi)$ or
the Jacobian). In the USR phase, $\eta>0$, and therefore the Jacobian term
grows as $\exp(\eta\,\delta N)$, while $P(\delta\phi)$ provides a
superexponential suppression for large $\delta N$ (via the term
$\exp\\{-[1-\exp(\eta\,\delta N)]^{2}\\}$ in (3.9)). The competition between
these two terms (dominated by the former initially, and by the latter for
large $\delta N$) yields a bump followed by a violent (superexponential)
cutoff. Moreover, this ultimately induces a dependence on $N_{i}$ in the PDF
of $\delta N$, due to the larger (or smaller) values of $\delta\phi$ needed
for $\bar{\phi}+\delta\phi$ to probe the USR phase of the unperturbed
trajectory if $N_{i}$ is chosen later (or earlier) on in the CR phase. This
behaviour is illustrated in the red curves of Fig. 2 for the two models under
consideration, together with the PDF resulting from the eternal CR
approximation in Eq. (3.6). As we will show shortly, the spurious bump in the
PDF disappears when momentum perturbations are appropriately accounted for.
Interestingly, once these are included, the result for the PDF of $\delta N$
obtained is actually closer to the one obtained by _unrealistically_ assuming
the CR phase of the unperturbed trajectory to extend indefinitely to the past
of $N_{i}$, although significant deviations from it can arise depending on the
specific shape of the phase space attractor close to the local maximum of the
potential, as we shall discuss.

Figure 2:  Left panel: Probability density for $\delta N$ for Model 1. The
black dashed-line corresponds to the eternal CR approximation in Eq. (3.6),
with $\eta$ in the CR phase following USR. The red lines correspond to a
numerical evaluation of Eq. (3.4) for $N_{i}-N_{*}=2$ (dotted), $4$ (dashed),
and $6$ (solid). The blue lines correspond to the solutions obtained by
including the momentum perturbations according to Eq. (3.11) for the same
values of $N_{i}-N_{*}$. The dependence of the PDF on the choice of $N_{i}$ is
evident. Right panel: Same as left panel but for Model 2. The red lines
correspond to $N_{i}-N_{*}=3$ (dotted), $5$ (dashed), and $7$ (solid). The
blue lines correspond to $N_{i}-N_{*}=3$ (solid) and $7$ (dashed).

###  3.2 Full calculation

Let us now consider the full calculation of $\delta N$, consistently
accounting for both field and momentum perturbations without making any
approximation. As we mentioned earlier, $N_{i}$ must be chosen as a time
posterior to the time at which the linear $\mathcal{R}=-\delta\phi/\bar{\pi}$
(or linear $\zeta$) becomes constant. Differentiating this last expression, we
obtain (see also [27]) the following relation between the fluctuations
$\delta\pi$ and $\delta\phi$:

$\delta\pi=\frac{\bar{\pi}^{\prime}(N_{i})}{\bar{\phi}^{\prime}(N_{i})}\delta\phi\,.$
(3.11)

In other words, the freezeout condition for linear perturbations (which is
necessary to apply consistently the $\delta N$ formalism) implies a relation
between $\delta\phi$ and $\delta\pi$, which allows to still use (3.2) to
compute the PDF of (the non-linear) $\zeta=\delta N$ using the PDF of
$\delta\phi$ as sole input. Yet, this relation is not, in general, the one
implicitly imposed assuming (3.4). To compute $\delta N$, one simply needs to
solve Eq. (2.3) scanning through the initial conditions consistently with
(3.11), in order to sample the PDF. Specifically, once we have determined
$N_{i}$ via the peak freezing condition –see the discussion above Eq. (3.1)–,
we perturb $\bar{\phi}(N_{i})$ by some value $\delta\phi$, and
$\bar{\pi}(N_{i})$ by the $\delta\pi$ given by (3.11). The corresponding
$\delta N$ is the difference in the number of realized $e$-folds of inflation
between the shifted trajectory and the unperturbed one. In this way, one
obtains a relation $\delta
N[\delta\phi,\delta\pi(\delta\phi)]=\mathcal{\zeta}(\delta\phi)$, and assuming
$\delta\phi$ to be Gaussian, one can again use (3.2) to obtain the PDF of
$\delta N=\zeta$. The PDF resulting from applying this prescription to the
models under consideration is depicted in the blue curves of Fig. 2.

Figure 3:  Left panel: Phase space for Model 1. For red trajectories, the
field turns back before reaching the local maximum and gets stuck in the local
minimum. The black dot represents the unstable equilibrium point $\phi_{\rm
max}$ at which the field sits in the local maximum with zero velocity. The
black, dashed line represents the unperturbed trajectory, and the blue line
denotes the CR line. We use $\pi_{\rm min}$ to denote the minimum velocity
reached in the unperturbed trajectory, and $\phi_{\rm min}$ denotes the
corresponding field value. Right panel: Same as left panel for Model 2.

It is instructive to analyse the structure of the phase space for both of our
models, depicted in Fig. 3.111111For the model in Eq. (2.8) we reconstruct the
potential as a function of the field using Eq. (2.6), and then solve the
Klein-Gordon equation with different initial conditions using the resulting
$V[N(\phi)]$. The tangent to any trajectory in phase space (in particular, the
unperturbed trajectory) is given by

$\frac{\mathop{}\\!\mathrm{d}\pi}{\mathop{}\\!\mathrm{d}\phi}=\epsilon-\eta\simeq-\eta\,.$
(3.12)

For field values smaller than the local maximum of the potential,121212We
recall that we assume (without loss of generality) that the field decreases as
inflation proceeds. the trajectories of the system are driven to an attractor,
see Fig. 3. In the thin grey trajectories, the field gets past the local
maximum of the potential and reaches its absolute minimum, whereas for the
thin red trajectories it gets stuck in the local minimum, leading to eternal
inflation. Notice how the grey trajectories converge to an attractor for
$\phi<\phi_{\max}$ (the local maximum of the potential), which, as it is to be
expected, at some point overlaps with the unperturbed trajectory (black
dashed). The tangent to this attractor, given by (3.12), reaches a phase of
constant slope (i.e. the attractor is a CR trajectory). This occurs almost
immediately to the left of $\phi_{\max}$ for Model 2 (right panel), but not
for Model 1 (left panel), where the $\eta$ of the attractor varies rapidly
before stabilizing. The existence of a region of non-constant slope of the
attractor in Model 1 close to the local maximum of the potential can be traced
back to $\eta_{V}(\phi)$ varying rapidly near that point.131313In the right
panel, the opposite occurs: $\eta_{V}(\phi)$ is close-to-constant (i.e. the
potential is well described by an inverted parabola) around the local maximum,
yielding a close-to-straight (i.e. CR) attractor right after the local
maximum. The blue line represents the _CR line_ , whose tangent vector is
given by (3.12) evaluated at $N_{i}$ (or equivalently $N>N_{i}$, since $\eta$
is already constant along the unperturbed trajectory at that point). The
points along this line indicate the initial conditions of a perturbed
trajectory $\\{\bar{\phi}(N_{i})+\delta\phi,\bar{\pi}(N_{i})+\delta\pi\\}$
consistent with (3.11). Choosing $N_{i}$ when the modes around the peak of the
power spectrum freeze– toward the beginning of the CR phase– corresponds to
fixing $\\{\bar{\phi}(N_{i}),\bar{\pi}(N_{i})\\}$ in the bottom left corner of
both panels of Fig. 3.

###  3.3 Comparison between the full and the approximate implementations of
$\delta N$

Fig. 3 illustrates the difference between the unperturbed-trajectory
approximation and the full calculation. If momentum perturbations $\delta\pi$
were chosen as in the unperturbed-trajectory approximation (Sec. 3.1), the
field would be shifted at $N_{i}$ to some point in the black dashed line, and
evolve hereon along said line. Instead, accounting for the momentum
perturbations as in (3.11), the field is shifted at $N_{i}$ to some point on
the blue line, and evolves hereon along the corresponding grey trajectory. For
small enough perturbations (that is, for sufficiently small $\delta\phi$, or
$\delta N$), the field shifted at $N_{i}$ remains on the overlapping locus of
the CR attractor and the unperturbed trajectory, and both approaches are
essentially equivalent. This is consistent with Fig. 2, as the PDF for $\delta
N$ obtained through all the different approaches agree sufficiently close to
the peak of the distribution. Deviations from the unperturbed-trajectory
$\delta N$ formalism due to the inclusion of $\delta\pi$ satisfying (3.11)
only appear for fluctuations that are large enough for the trajectory of the
perturbed field to become sensitive to the region of phase space beyond the
turnaround of the unperturbed trajectory, where the CR line (solid blue) and
the unperturbed trajectory (black dashed) coincide no more. Finally, notice
that when momentum perturbations are chosen consistently with (3.11), values
of $\delta N\to\infty$ are sampled for
$\delta\phi\to\bar{\phi}(N_{i})-\phi_{\max}$ (more specifically, for the
$(\delta\phi,\delta\pi)$ at which the blue line crosses from the grey region
to the red region in Fig. 3). Conversely, in the unperturbed-trajectory
approximation, $\delta N\to\infty$ is only achieved as $\delta\phi\to\infty$
along the unperturbed trajectory of the field.

One can notice in Fig. 2 that the PDF for Model 2 obtained by sampling
$\delta\phi$ and $\delta\pi$ consistently with (3.11) is remarkably close to
the one obtained using the eternal CR approximation in Eq. (3.6) for a wide
range of $N_{i}$ values, and only begins to deviate for very late (and
therefore inadequately chosen) values of $N_{i}$, long after the modes
corresponding to the peak of the power spectrum have frozen. The reason for
this behaviour is that the eternal CR approximation effectively corresponds to
shifting the field at $N_{i}$ along the CR line (with tangent vector given by
(3.12) evaluated at $N_{i}$) and assuming it evolves hereon along the same CR
line. This approximation is indeed an accurate description of the actual
dynamics of Model 2, shown in the right panel of Fig. 3, where the CR line
(solid blue) and the phase space attractor (convergence of grey lines) lie on
top of each other. However, this is not necessarily the case for other
implementations of USR, and in fact it does not occur for Model 1, as shown in
the left panel of the same figure. In that case, the field originally shifted
along the CR line does not evolve afterwards along it, but along the grey
curves that intersect it. This translates into notable differences in the
deep-tail regime when comparing the PDF obtained with the eternal CR
approximation and the one obtained by including the momentum perturbations as
per (3.11), as seen in Fig. 2.

##  4 Stochastic inflation

The framework of stochastic inflation [28] splits the inflaton into short and
long-wavelength Fourier modes in order to study the backreaction of the former
on the latter. The modes with small wavelength evolve as free fields,
relatively insensitive to the non-linearities induced by the interaction terms
in the scalar potential. Hence, it is assumed that the evolution of these
short modes is well captured by linear, deterministic, equations of motion.
Their effect on the long modes is then accounted for by stochastic source
terms which are assumed to behave as white noise, due to the underlying
Gaussian statistics assumed for the short modes and the choice of the coarse-
graining window that separates short and long modes (see [29] for a review).
This separation is done introducing a (small) coarse-graining parameter
$\sigma$ such that a mode of wavelength $k$ is considered a long one at any
given time provided its wavenumber is smaller than a quantity

$\displaystyle k_{\sigma}\equiv\sigma aH\,.$ (4.1)

In the presence of an USR phase, $\sigma$ must be chosen small enough in such
a way that the modes around the peak of the power spectrum141414These modes
provide the largest kicks, as it will become clear in (4.7). have had enough
time to freeze out and classicalize when they become part of the coarse-
grained field151515The modes around the peak (which undergo the most
superhorizon enhancement due to USR) must not become part of the coarse-
grained field until they freeze, which occurs slightly after the end of USR.
Thus, if such a mode takes $\Delta N_{\rm freeze}$ $e$-folds to freeze,
$\sigma$ should be, at least, $e^{-\Delta N_{\rm freeze}}$. Later on we will
illustrate this in an example, and explain how it relates to the choice of
$N_{i}$ in classical $\delta N$. [18] (see also [30] for a recent discussion).
However, it has been argued that $\sigma$ cannot be chosen arbitrarily small,
as it would result in the loss of information about the non-linear evolution
of the curvature perturbation [31]. There is therefore a tradeoff between
choosing $\sigma$ too small and too large. In the following, we use
$\sigma=0.1$ and $\sigma=0.01$ as two representative values for the USR phases
of Models 1 and 2 with duration $\Delta N\approx 2.5\approx-\log(0.08)$.

Once the quantum fluctuations of the field freeze out they can be treated as
classical random variables [18]. These fluctuations backreact on the
background trajectory, leading to a system of stochastic differential
equations which reads

$\displaystyle\mathop{}\\!\mathrm{d}\phi$
$\displaystyle=\pi\mathop{}\\!\mathrm{d}N+\mathop{}\\!\mathrm{d}W_{\phi}\,,$
(4.2) $\displaystyle\mathop{}\\!\mathrm{d}\pi$
$\displaystyle=-(3-\epsilon)(\pi+\partial_{\phi}\log
V)\mathop{}\\!\mathrm{d}N+\mathop{}\\!\mathrm{d}W_{\pi}\,,$ (4.3)

where $\phi$ and $\pi$ are now stochastic variables and
$\mathop{}\\!\mathrm{d}W_{\phi,\pi}$ denote Wiener increments, often expressed
in terms of stochastic “kicks”
$\xi_{\phi,\pi}=\mathop{}\\!\mathrm{d}W_{\phi,\pi}/\mathop{}\\!\mathrm{d}N$
satisfying

$\displaystyle\langle\xi_{\phi}(N)\xi_{\phi}(N^{\prime})\rangle=\mathcal{P}_{\delta\phi}(N,k_{\sigma})\delta(N-N^{\prime})\,,$
(4.4)
$\displaystyle\langle\xi_{\pi}(N)\xi_{\pi}(N^{\prime})\rangle=\mathcal{P}_{\delta\pi}(N,k_{\sigma})\delta(N-N^{\prime})\,.$
(4.5)

In what follows we assume that the power spectra appearing in the variance of
these noises is computed solving (2.11) with the coefficients evaluated in the
background trajectory. In principle, a full treatment of stochastic inflation
would involve using the stochastic variables in the coefficients and
recalculating the power spectrum at every time step, as in [8]. As the results
of [8] imply, however, not only does this require plenty of computational
power, but, more importantly, does not seem to have any effect on the final
result for the probability distribution. Using the background quantities for
the coefficients it is therefore a practical working approximation which
should provide excellent accuracy. We will comment on the reasons for this
behaviour at the end of this section.

Despite having two seemingly independent Wiener increments in the equations of
motion, there is actually only one independent kick at each time step. This is
due to the fact that the quantum state of the short-wavelength perturbations
is highly squeezed [32, 8, 27, 13], which forces the kicks to be aligned in a
very specific way, with $\xi_{\pi}=(\delta\pi_{k}/\delta\phi_{k})\xi_{\phi}$.
Once the system enters the last CR phase, the (linear) curvature perturbation
freezes to a time-independent value on super-horizon scales. By
differentiating the linear expression $\mathcal{R}=-\delta\phi/\bar{\pi}$, we
obtain the relation161616Cf. (3.11), where an analogous relation for
field/momentum perturbations holds.

$\dfrac{\xi_{\pi}}{\xi_{\phi}}=\dfrac{\overline{\pi}^{\prime}}{\overline{\phi}^{\prime}}\,,$
(4.6)

where the quantities carrying a bar over are evaluated on the unperturbed
trajectory, and are therefore deterministic. The variance of the Wiener
increment at each time step can be expressed as (for the case of $\phi$ and
analogously for $\pi$):

$\langle\mathop{}\\!\mathrm{d}W_{\phi}^{2}\rangle=\overline{\pi}(N)^{2}\mathcal{P}_{\cal\zeta}(k_{\sigma})\frac{\mathop{}\\!\mathrm{d}\log
k_{\sigma}}{\mathop{}\\!\mathrm{d}N}\mathop{}\\!\mathrm{d}N\,,$ (4.7)

where $\mathcal{P}_{\zeta}$ represents the linear power spectrum of
$\zeta=\mathcal{R}$ after the modes have frozen (recall Footnote 9). Since the
constant parameter $\sigma\ll 1$ is $k$-independent, the field receives a kick
from a single scale at every time step.171717See [30] for a recent proposal on
the topic. By differentiating $k_{\sigma}$, we find
$\mathop{}\\!\mathrm{d}\log
k_{\sigma}/\mathop{}\\!\mathrm{d}N=1-\epsilon\simeq 1$. We remark that the
factor $\overline{\pi}(N)^{2}$ denotes the momentum of the unperturbed
background field (i.e. the solution to eqs. (4.2, 4.3) without the stochastic
noise terms), and it takes into account the time dependence of the power
spectrum
$\mathcal{P}_{\delta\phi}(N,k)=\overline{\pi}(N)^{2}\mathcal{P}_{\zeta}(k)$
where, if $\sigma$ has been chosen adequately, we can assume that
$\mathcal{P}_{\cal\zeta}(k)$ has settled to its frozen value.

The resulting shift in the duration of inflation $\delta N$ induced by the
stochastic kicks can be extracted at late times, once the peak modes have
frozen out (and hence have become coarse-grained), and therefore subsequent
kicks become inefficient. This can be accomplished in practice by turning off
the noise terms in the stochastic equations once they become irrelevant in
comparison to the drift term, and integrating the hereafter unperturbed (i.e.
noiseless) system until inflation ends (that is, until $\epsilon=1$). The
value of $\delta N$ can then be extracted by comparing the duration of
inflation in this realization to that of the unperturbed trajectory. We remark
that it is important to integrate the system until the end of inflation since,
even if the noise is turned off, the field may not be in the unperturbed
trajectory (or even in the phase-space attractor) if integration is cut short.
In the limit in which the linear power spectrum of $\zeta$ is strongly peaked
at some particular scale, the field receives the most significant stochastic
kicks over a reduced amount of time. This corresponds to a single (possibly
non-Gaussian) kick in the extreme case of a Dirac delta spectrum. The latter
case is exactly equivalent to the standard $\delta N$ calculation of the
previous section. This can be checked numerically by manually adjusting the
width of the spectrum. It is therefore reasonable to expect both formalisms to
yield very similar results for inflection-point models leading to highly
peaked spectra.

This discussion also makes clear that there is an approximate one-to-one
correspondence between the parameter $\sigma$ in stochastic inflation and the
parameter $N_{i}$ that denotes the time at which the field is perturbed in the
classical $\delta N$ formalism. Indeed, for both results to match, $N_{i}$
must be chosen roughly as the time at which the peak in the power spectrum is
coarse-grained, $k_{\rm peak}\simeq\sigma a(N_{i})H(N_{i})$. We reiterate,
however, that this correspondence is only valid in the limit in which the
power spectrum is highly peaked.

As in the classical $\delta N$ formalism, one can attempt to simplify the
system of equations by assuming that the velocity of the scalar field adheres
to the background trajectory at every time step, as suggested in [27, 13].
This would indeed be the case if the stochastic kicks never took the field far
from the phase space attractor. As we saw in the previous section for the case
of the classical $\delta N$ formalism, this is a good approximation close to
the peak of the PDF of $\delta N$. Whether it is a good approximation in the
deep tail, however, strongly depends on the morphology of the phase space, and
in particular, the variation of the slope of the attractor for
$\phi<\phi_{\max}$. The corresponding equation implementing this approximation
reads

$\mathop{}\\!\mathrm{d}\phi=\overline{\pi}_{c}(\phi)\mathop{}\\!\mathrm{d}N+\mathop{}\\!\mathrm{d}W_{\phi},$
(4.8)

where we have used $\overline{\pi}_{c}$ to denote the momentum of the
unperturbed trajectory written as a function of the field value (and not as a
function of $N$).181818To be crystal clear, one should solve the background
equation to determine $\overline{\pi}(N)$ and $\overline{\phi}(N)$, invert the
latter to obtain $\overline{\pi}[N(\overline{\phi})]$, and then make the
substitution $\overline{\phi}\to\phi$ in this deterministic equation to obtain
$\overline{\pi}_{c}=\overline{\pi}[N(\phi)]$, which is now a stochastic
variable. This condition ensures that the scalar field stays on the
unperturbed trajectory, denoted with a black dashed line in Fig. 3, at every
time step. Let us emphasize that, in this scheme, the amplitude of the
stochastic kicks in Eq. (4.7) is still determined by $\overline{\pi}(N)$, and
not $\overline{\pi}_{c}(\phi)$. The former is a deterministic quantity,
whereas the latter is a stochastic variable due to its dependence on the
stochastic variable $\phi$. Indeed, if the stochastic term were evaluated
using $\overline{\pi}_{c}(\phi)$ instead of $\overline{\pi}(N)$, the induced
stochastic time shift $\delta N$ would be necessarily Gaussian. This statement
can be made clearer by transforming Eq. (4.8) into a stochastic equation for
the time shift $\delta N$ [27, 13]. This equation can be expressed as

$\mathop{}\\!\mathrm{d}\delta
N=-\frac{\overline{\pi}^{\prime}(N)}{2\overline{\pi}(N)}\mathcal{P}_{\zeta}(k_{\sigma})(1-\epsilon)^{2}\mathop{}\\!\mathrm{d}N+\mathop{}\\!\mathrm{d}W_{\delta
N},$ (4.9)

where the first term follows from Ito’s lemma, and

$\langle\mathop{}\\!\mathrm{d}W_{\delta
N}^{2}\rangle=\frac{\langle\mathop{}\\!\mathrm{d}W_{\phi}^{2}\rangle}{\overline{\pi}_{c}^{2}(\phi)}.$
(4.10)

It is clear from Eq. (4.7) and Eq. (4.9) that replacing the momentum
$\overline{\pi}(N)$ of Eq. (4.7) by $\overline{\pi}_{c}(\phi)$ would entirely
remove the path-dependence from the stochastic process and lead to a Gaussian
distribution for $\delta N$.

Let us now return to the case of interest for us, which is the solution to the
full system of Eqs. (4.2, 4.3). In the stochastic approach, the trajectory in
phase space follows a random walk dictated by these equations, which can be
recast as

$\displaystyle\begin{pmatrix}\phi^{\prime}\\\
\pi^{\prime}\end{pmatrix}=\mathbfcal{D}+\begin{pmatrix}\xi_{\phi}\\\
\xi_{\pi}\end{pmatrix}\,,$ (4.11)

where $\mathbfcal{D}$ is called the drift vector and (see Eqs. (4.6) and
(3.12))

$\begin{pmatrix}\xi_{\phi}\\\
\xi_{\pi}\end{pmatrix}\propto\hat{\xi}\begin{pmatrix}1\\\
-\eta\end{pmatrix}\,,$ (4.12)

where $\hat{\xi}$ is a Gaussian-distributed noise. Eq. (4.12) makes manifest
that the CR line (blue line in Fig. 3) gives the direction of the stochastic
kicks corresponding to the modes that are coarse-grained at $N_{i}$, see Sec.
3.2.191919These modes are the ones providing the _larger_ kicks, as $N_{i}$
corresponds to the freezing-out of the modes around the peak of the power
spectrum.  Also, the drift vector $\mathbfcal{D}$ is tangent to the
trajectories in phase space (grey lines in Fig. 3). If the attractor
trajectory coincides with the CR line (as it happens for Model 2), the field
remains on this line, kicked back and forth by the noise and driven along it
by the drift. The system is therefore effectively described by a single degree
of freedom. In this scenario, the dependence on $\sigma$ is lost; in other
words, it does not matter when specifically the kicks are applied, but only
their accumulated effect, as already noticed in [13]. This can be understood
by analytically solving the single Langevin equation that governs the inflaton
field constrained to drift and diffuse along the superposition of the CR line
and CR attractor with constant $\eta$:

$\phi^{\prime}=-\eta\,\phi+\sqrt{2\epsilon(N)\mathcal{P}_{\mathcal{\zeta}}(k_{\sigma})}\xi_{\phi}\,.$
(4.13)

The solution is a Gaussian $\phi(N)=\bar{\phi}(N)+\delta\phi(N)$ with

$\langle\phi(N)\rangle=\bar{\phi}(N)\,,$ (4.14)

and

$\sigma^{2}_{\delta\phi}=2\epsilon(N)\int_{N_{\rm
IR}}^{N}\mathop{}\\!\mathrm{d}N^{\prime}\,\mathcal{P}_{\zeta}(k_{\sigma}(N^{\prime}))=2\epsilon(N)\int_{\log
k_{\rm IR}}^{\log k_{\sigma}(N)}\mathop{}\\!\mathrm{d}\log
k\,\mathcal{P}_{\zeta}(k)\,.$ (4.15)

If evaluated at $N_{i}$, the Gaussian $\phi(N_{i})$ described by (4.14),
(4.15) coincides with the classical $\delta N$ initial condition, given by
$\bar{\phi}(N_{i})+\delta\phi$, with $\delta\phi$ a Gaussian variable with
variance (3.8). Choosing $N_{i}$ at a time such that the modes in the peak of
$\mathcal{P}_{\mathcal{\zeta}}$ have frozen means that (4.15) has already
received its main contribution from the maximum of its integrand, and
therefore the time dependence through the upper limit of integration is
lost.202020The lower limit of integration, which in practice corresponds to an
infrared cutoff, is the same for the integral in (4.15) and the one in (3.8).
From the point of view of classical $\delta N$, this scenario corresponds to a
realization of the eternal CR case.212121See the discussion in the last
paragraph of Sec. 3, where the dependence of the PDF of $\delta N$ on $N_{i}$
is lost (see Eq. (3.9)) provided it is enough for the peak modes of
$\mathcal{P}_{\mathcal{R}}$ to have frozen. From the point of view of
stochastic $\delta N$, this corresponds to an approximate independence on the
choice of $\sigma$ (Fig. 4, right panels), again provided it is chosen small
enough for the modes of the peak to be coarse-grained once they have frozen.

If the attractor trajectory is not aligned with the CR line (as it happens for
Model 1), the situation is different. While the noise vector kicks the field
along the CR line, the drift vector tends to drive the system toward the
attractor (where the grey lines converge). These two terms are not aligned and
it is not possible to condense the system into a single variable. In general,
the field and the momentum experience a random walk in the region of phase
space between the CR line and the attractor trajectory which is not correctly
captured by the eternal CR approximation (hence the change of slope of the
exponential tail of the PDF of $\delta N$, see Fig. 4, left panel). In
principle, the equivalence between the classical and stochastic $\delta N$
formalisms is not exact, as the argument in (4.13) does not extend to this
case. Nevertheless, the fact that the largest stochastic kicks are clustered
in a short time period (because of the power spectrum having a narrow peak),
with a relatively constant $\eta$, still makes both approaches (classical and
stochastic) very similar, provided $N_{i}$ (in classical $\delta N$) and
$\sigma$ (in stochastic $\delta N$) are chosen consistently. Unlike in the
previous case, there is a dependence of the PDF of $\delta N$ on $N_{i}$ for
(classical $\delta N$) and $\sigma$ (for stochastic $\delta N$). The reason
for this is that the earlier (or later) we chose $N_{i}$, the smaller (or
larger) values of $\delta N$ will begin to probe the region of phase space in
which the CR line and the attractor do not overlap anymore.

Moreover, in the stochastic approach, the field and its momentum can, in
principle, explore a region of parameter space inaccessible to the classical
$\delta N$ formalism, namely, the red regions in Fig. 3. This would be the
case for stochastic realizations of the system in which the field receives
kicks that take it to the region $\phi>\phi_{\max}$, but then it also receives
enough kicks that bring it back to $\phi<\phi_{\max}$ and thus does not get
stuck in the local minimum, as would happen classically. This is a difference
between the two formalism (classical and stochastic). The stochastic $\delta
N$ formalism could in principle describe these events more accurately. These
effects are, however, largely mitigated in the case in which the power
spectrum is very peaked, and thus the kicks occur in a narrow window of time.
This is in fact the case for Model 1, for which both formalisms yield very
close results, as we will now show.

We computed the PDF of $\delta N$ numerically by running $10^{6}$ realizations
of the Langevin system in Eqs. (4.2, 4.3). The results for Model 1 are shown
on the left panel of Fig. 4 for the two representative values $\sigma=10^{-1}$
(orange crosses) and $10^{-2}$ (blue crosses). We have also solved the
simplified equation (4.8) for $\sigma=10^{-2}$ (red dots). The results for
Model 2 are shown on the right panel of the same figure. The black dashed line
represents the eternal CR approximation of Eq. (3.6). The black solid line
corresponds to the result obtained using the classical $\delta N$ formalism
with momentum perturbations included, i.e. by following the procedure outlined
below Eq. (3.11), with $N_{i}-N_{*}=3.5$.222222The precise value of $N_{i}$
that makes the classical and stochastic calculations agree cannot be exactly
determined a priori. However, as explained before, a reasonable estimate for
it can be obtained by taking into account the following. For $\sigma=10^{-2}$,
the mode corresponding to the peak of the power spectrum, $k_{\rm peak}$, is
coarse-grained $-\log(10^{-2})\approx 4.6$ $e$-folds after crossing the
horizon, which, for our model, corresponds to an $N$ such that $N-N_{*}=3.2$
$e$-folds. If the peak were infinitely narrow, we would therefore take
$N_{i}-N_{*}=3.2$. However, since the peak of the spectrum has a certain
width, we need $N_{i}-N_{*}\gtrsim 3.2$ to account for the freezing of modes
close to the peak with $k\lesssim k_{\rm peak}$. Numerically, we found to best
match to occur for $N_{i}-N_{*}=3.5$.

Figure 4:  Left panel: Probability density function for the model of Eq.
(2.8). The black dashed-line denotes the eternal CR approximation. The black
solid line is obtained by using the $\delta N$ formalism including the
momentum perturbations, with $N_{i}-N_{*}=3.5$. The blue and orange crosses
represent $10^{6}$ stochastic realizations of the system of two Langevin
equations with $\sigma=10^{-2}$ and $10^{-1}$, respectively, and the red dots
correspond to realizations of the simplified Langevin equation with
$\sigma=10^{-2}$. Right panel: Same as left panel for the model in Eq. (2.9)
with $\sigma=10^{-2}$.

As shown in Fig. 4, both the stochastic and classical versions of the $\delta
N$ formalism (with momentum perturbations included) lead to extended, non-
Gaussian tails for the PDF of the curvature perturbation. In the case of Model
1, the eternal CR approximation of Eq. (3.6) fails to correctly capture the
tail behaviour due to the complexity of the phase space. Similarly, the
simplified, single-Langevin equation (4.8) deviates from the full result in
the deep-tail regime for the reasons highlighted earlier. On the other hand,
for Model 2 the phase space attractor coincides with the CR line, as shown in
Fig. 3, and therefore all approaches yield essentially the same
result.232323We remark that, due to the steeper decrease of the PDF for this
model, studying the deep-tail regime with the stochastic approach is much more
costly, since probing large values of $\delta N$ requires many more
realizations. This figure also illustrates the importance of the choice of
$\sigma$ in the stochastic approach, which mirrors the choice of $N_{i}$ in
the standard $\delta N$ calculation. Model 2 is less sensitive to this choice
than Model 1 due to the simpler structure of the phase space.

Before closing this section, let us comment on the effects of backreaction we
mentioned earlier. As we stated before, a full treatment of stochastic
inflation would involve recalculating the power spectrum, which appears in the
amplitude of the noise coefficients, at every time step for each realization.
As we already mentioned, this was done in [8], and no deviation was found in
the PDF with respect to the result obtained by computing the spectrum once by
evaluating the coefficients in the background. The reason for this behaviour
is simply that by the time the kicks begin to act on the field during the CR
phase, the modes around the peak of the spectrum are coarse-grained and
frozen. Thus, these modes of the power spectrum remain the same independently
of changes in the background field and velocity. Moreover, the value at which
these modes have frozen was determined by the background evolution in a period
when the noise was negligible (i.e. before the peak of the power spectrum was
coarse-grained). The only way this backreaction could have some effect is for
very rare realizations in which the field is pushed back beyond the local
maximum, and thus a second USR phase begins, momentarily changing the sign of
$\eta$, until the field is driven back to CR so that inflation can end. Even
in this case, further super-Hubble evolution for $\zeta$ is highly unlikely
due to the fact that, at this point, the second mode in Eq. (2.13) has been
decaying since the first USR phase ended, and thus the second USR phase would
need to either occur very shortly after the first or have an extreme duration,
and the probability of this happening is significantly suppressed (as
illustrated by our results).

##  5 Summary and conclusions

We have computed the PDF of the non-linear curvature perturbation on uniform-
density hypersurfaces, $\zeta=\delta N$, in the presence of an USR phase using
the classical and stochastic $\delta N$ formalisms. We have studied the degree
of equivalence of these methods, as well as the validity of several
approximations for both of them. We considered two different models, with two
qualitatively different phase spaces, to understand how the morphology of
phase space affects our conclusions.

In Sec. 3, we applied the classical $\delta N$ formalism, imposing consistent
conditions on field and momentum perturbations as per (3.11), to numerically
compute the PDF of $\zeta$. As illustrated in Fig. 3, Eq. (3.11) is equivalent
to shifting the field at $N_{i}$ along a specific line in phase space that we
refer to as the _CR line_. We addressed the validity of two approximations:

  * •

The _eternal CR approximation_ , which assumes that the field and momentum
perturbations lie along the CR line and do not depart from it throughout their
evolution. This approximation works well for models in which the CR line and
the phase space attractor (approximately) overlap. In our examples, this is
the case for _Model 2_ , but not for _Model 1_. In the cases where this
approximation holds, accurate analytical estimates of the PDF, as in Eq.
(3.9), can be obtained.

  * •

The _unperturbed trajectory approximation_ , which assumes the perturbations
are such that the field and its momentum evolve along the unperturbed
trajectory in phase space. We found that such an approximation has a finite
regime of validity, correspondent to the region where the CR line and the
phase space attractor overlap, but is bound to eventually fail at some value
of $\delta N$. Moreover, we argued that this failure manifests itself in the
PDF of $\delta N$ in the form of a (spurious) bump.

In Sec. 4, we applied the stochastic $\delta N$ formalism, again imposing
consistent conditions on field and momentum noises as per Eq. (4.6). In the
same spirit of the unperturbed-trajectory approximation, we studied the
outcome of simplifying the system of Langevin equations by making the field
drift and diffuse along the unperturbed trajectory, as in (4.8). Similarly to
the classical case, we found the approximation to fail eventually. Also in
Sec. 4, we analysed the equivalence of the stochastic and classical $\delta N$
formalisms. From here, we extracted several important conclusions:

  1. 1.

For the stochastic and classical $\delta N$ formalisms to be equivalent, a
consistency relation needs to hold between $N_{i}$ (which indicates the
initial hypersurface in the classical $\delta N$ formalism) and $\sigma$ (the
coarse-graining parameter in the stochastic $\delta N$ formalism). For a
peaked power spectrum, this relation is $k_{\rm peak}\approx\sigma
a(N_{i})H(N_{i})$, where $k_{\rm peak}$ is the mode corresponding to the
maximum of the power spectrum. The relation is physically meaningful, as both
parameters indicate (in their respective formalisms) the delay from the
horizon crossing of $k_{\rm peak}$ to its actual freezing. From a more
technical point of view, this freezing is a necessary condition for the
separate universe approach (underlying in both classical and stochastic
$\delta N$ formalisms) to hold. The reason why the relation is approximate and
not equal is that the peak of the spectrum is not strictly monochromatic, but
slightly extended in $k$, as discussed in Footnote 22.

  2. 2.

For models with a phase space in which the attractor and the CR line
approximately overlap (as in _Model 2_), the equivalence between classical and
stochastic $\delta N$ formalisms is exact, as the system is effectively
described by dynamics in one dimension (given by the CR line). Moreover, in
this case, both approaches are equivalent to the eternal CR approximation for
classical $\delta N$, where the dependence on the final PDF of $\delta N$
(which has a single-slope exponential tail) on $N_{i}$ (equivalently, on
$\sigma$) is lost, as it can be shown analytically (see (3.9)). This is in
agreement with the results in [13].

  3. 3.

For models with a phase space in which the attractor has non-constant slope,
and thus leaves the CR line, there is still an approximate equivalence between
stochastic and classical $\delta N$ formalisms, ensured by the narrowness of
the peak of the power spectrum (in short, this implies that most stochastic
kicks occur during a short time range where $\eta$ of the attractor is almost
constant). However, in this case, the tail of the PDF is given by a non-
constant-slope exponential, related to the non-constant slope attractor in
phase space. The position of the breaking points between different slopes
depends on the choice of $N_{i}$ (or $\sigma$), therefore in this case the
result is not $N_{i}-$independent.

A corollary of the last two points is that the non-Gaussian tails in the PDF
of the curvature perturbation $\zeta$ can be reproduced without invoking the
stochastic inflation formalism when momentum fluctuations are properly taken
into account. These tails are a consequence of the non-linear relation between
$\zeta$ and $\delta\phi$. The reason why they also arise in stochastic $\delta
N$ is because of the implementation of the $\delta N$ formalism into
stochastic inflation, not because of the stochastic formalism itself. From a
practical point of view, this implies that the non-Gaussian tail of the PDF of
$\zeta$ can be reproduced with a simple procedure that takes only a few lines
of code and a negligible amount of computational power. Moreover, depending on
the structure of the phase space, it may not be necessary to perform any
numerical calculations at all, and the simple expression in Eq. (3.6) can be
used to obtain the PDF analytically, as previously discussed in
[13].242424Note that, in [13], a deviation is found between the results
obtained by using the simplified equation (4.8) and the eternal CR
approximation (3.6). This difference disappears once one considers the full
system (4.2, 4.3) or, alternatively, the classical $\delta N$ formalism with
momentum given by (3.11).

As it has been emphasized, a key property to understand the relation between
the classical and stochastic $\delta N$ formalisms and their possible
simplifications is the constant slope (or lack thereof) of the phase space
attractor for $\phi<\phi_{\max}$ (or equivalently, whether the attractor and
the CR line lie on top of each other or not). As already mentioned, this
property can be traced back to the second derivative of the potential. If it
is close to constant around the local maximum, this maximum is well described
by an inverted parabola, for which the equation of motion of the inflaton
becomes linear (see (2.7)), with an attractor with (asymptotically) constant
$\eta$. This is what happens in _Model 2_. In _Model 1_ , the situation is the
opposite. The second derivative of the potential varies around the local
maximum, meaning that the equation of motion of the inflaton is non-linear,
and the attractor is not characterized by a constant value of $\eta$.
Likewise, the values of $\eta$ during USR and CR are not related by the Wands
duality. In principle, it is to be expected that a single-field potential with
a local minimum - local maximum configuration whose potential can be described
with a small number of free parameters, like the radiative plateau potential
in [2], will be reasonably well approximated by an inverse parabola around its
local maximum, and hence share the properties found for _Model 2_.
Reciprocally, potentials with many free parameters and complicated functional
forms, or reverse-engineered from the evolution of some time-dependent
quantity (as it is indeed the case of _Model 1_), can develop the issues found
in _Model 1_.

Another conclusion of our work is that, unless the unperturbed-trajectory
approximation is imposed (which, as already discussed, is unreliable when
probing the deep tail of the PDF of $\zeta$), the field might get stuck in the
local minimum of the potential. This can happen equivalently in the classical
$\delta N$ formalism ($\delta\phi$ and $\delta\pi$ being large enough for
$(\bar{\phi}(N_{i})+\delta\phi,\bar{\pi}(N_{i})+\delta\pi)$ to be in the red
region of Fig. 3) and the stochastic one (an accumulation of kicks parallel to
the CR line bringing the field to the red region of Fig. 3, where the drift
term carries it to the local minimum252525In this case, there is a small
probability of the noise bringing the field back again to the grey region of
Fig. 3.). Since the potential is designed to lead to a substantial slowdown of
the scalar field, the probability of this happening is not, in principle,
totally negligible. For the models under consideration, we found
$\mathcal{O}(1\%)$ of the trajectories end up trapped in the minimum.

Finally, let us reiterate that the results obtained using both stochastic
inflation and the classical $\delta N$ formalism rely crucially on the
assumption that the field fluctuations $\delta\phi$ follow a Gaussian
distribution. Studies about this assumption in the presence of a USR phase are
scarce, particularly in the deep-tail regime, where these deviations would
become more relevant. We remark, however, that the $\delta N$ formalism might
provide an additional advantage in this context, since it is able to handle
any non-Gaussian distribution for $\delta\phi$, as is clear from the simple
relation in (3.2), whereas the Gaussianity of $\delta\phi$ is essentially
built into the standard formulation of stochastic inflation in terms of
Langevin equations and is therefore much more difficult to incorporate into
the formalism.

## Acknowledgments

The work of GB and APR has been funded by the following grants: 1) Contrato de
Atracción de Talento (Modalidad 1) de la Comunidad de Madrid (Spain),
2017-T1/TIC-5520 and 2021-5A/TIC-20957, 2) PID2021-124704NB-I00 funded by
MCIN/AEI/10.13039/501100011033 and by ERDF A way of making Europe, 3)
CNS2022-135613 MICIU/AEI/10.13039/501100011033 and by the European Union
NextGenerationEU/PRTR, 4) the IFT Centro de Excelencia Severo Ochoa Grant
CEX2020-001007-S, funded by MCIN/AEI/10.13039/501100011033. APR has been
supported by Universidad Autónoma de Madrid with a PhD contract contrato
predoctoral para formación de personal investigador (FPI), call of 2021, and a
grant Programa de ayudas UAM-Santander para la movilidad de jóvenes
investigadores, call of 2023. APR thanks DESY for hospitality during part of
the realization of this work. MP and JR acknowledge support by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy – EXC 2121 “Quantum Universe” - 390833306. MP acknowledges
support by the DFG Emmy Noether Grant No. PI 1933/1-1. We thank S. Céspedes,
J. Gambín Egea and E. Tomberg for discussions.

##  Appendix A The $\delta N$ formula

In this appendix, we discuss in more detail the arguments and assumptions
leading to the $\delta N$ formula (3.1). The starting point is the ADM
decomposition of the metric, which we write following the standard
parameterization used e.g. in [19]

$\mathop{}\\!\mathrm{d}s^{2}=-\alpha\,\mathop{}\\!\mathrm{d}t^{2}+\gamma_{ij}(\mathop{}\\!\mathrm{d}x^{i}+\beta^{i}\,\mathop{}\\!\mathrm{d}t)(\mathop{}\\!\mathrm{d}x^{j}+\beta^{j}\,\mathop{}\\!\mathrm{d}t)\,,$
(A.1)

where $\alpha$ and $\beta^{i}$ are respectively the lapse and shift functions,
and

$\gamma_{ij}=a^{2}\,e^{2\psi}\,\tilde{\gamma}_{ij}\,,\quad\tilde{\gamma}_{ij}=(e^{h})_{ij},\quad\text{Tr}(h)=0\,.$
(A.2)

Given the normal 1-form to $t=$ constant hypersurfaces,
$n_{\mu}=(-\alpha,0,0,0)$, and denoting $D_{i}$ the covariant derivative
associated to the spatial metric $\gamma_{ij}$, the extrinsic curvature on
spatial hypersurfaces can be written as [33]

$K_{ij}=\nabla_{i}n_{j}=\frac{1}{2\alpha}\left(-\dot{\gamma}_{ij}+D_{i}\beta_{j}+D_{j}\beta_{i}\right)=\frac{K}{3}\gamma_{ij}+a^{2}e^{2\psi}\tilde{A}_{ij}\,,$
(A.3)

where $\tilde{A}_{ij}$ is traceless and $K$ is the trace of the extrinsic
curvature, given by

$K=-\frac{1}{\alpha}\left[3(H+\dot{\psi})-D_{i}\beta^{i}\right],\quad
H(t)=\dot{a}(t)/a(t)\,.$ (A.4)

The trace of the extrinsic curvature can be used to define the _generalized
number of $e$-folds_,

$\mathop{}\\!\mathrm{d}\mathcal{N}=-\frac{\alpha\,K}{3}\mathop{}\\!\mathrm{d}t\,,$
(A.5)

which reduces to $\mathop{}\\!\mathrm{d}N=H\,\mathop{}\\!\mathrm{d}t$ for a
FLRW metric.

In the ADM formalism, the quantities $\psi$, $\tilde{A}_{ij}$, $K$ and
$\tilde{\gamma}_{ij}$, together with the inflaton $\phi$, are dynamical
degrees of freedom. Their equations of motion are supplemented with two
constraint equations obtained varying the action (2.1) with respect to the
lapse and the shift [34], known as _Hamiltonian_ and _momentum constraints_
[33, 19]. Instead of expanding the equations in powers of the $\psi$,
$\tilde{A}_{ij}$, $K$, $\tilde{\gamma}_{ij}$ and
$\delta\phi=\phi-\bar{\phi}(t)$ 262626In this appendix, $\bar{\phi}$ denotes
the solution of the scalar field equation when (A.1) is a FLRW metric. –which
would then be assumed to be small fluctuations, as it is done in standard
perturbation theory–, we can expand them in spatial gradients, denoting with
powers of $\epsilon_{g}$ the different orders in the gradient expansion.
Working in a gauge in which $h_{ij}$ is transverse and assuming that the
anisotropic stress is $\mathcal{O}(\epsilon_{g}^{2})$, one can show [26, 19]
that $\tilde{A}_{ij}$, $\dot{h}_{ij}$ and $\partial_{j}\beta^{i}$ are also
$\mathcal{O}(\epsilon_{g}^{2})$. Using this result, and changing the time
variable to the generalized number of $e$-folds, one finds that the
Hamiltonian constraint and the equation of motion for $\phi$ have the same
functional forms as their unperturbed counterparts (for which we assume a FLRW
metric) at lowest order in the gradient expansion:272727Second order gradients
are _not_ negligible in an USR phase, already at linear order in perturbations
[35, 36]. This is consistent with the fact that both $N_{i}$ (in the classical
$\delta N$ formalism) and $\sigma$ (in the stochastic one) have to be chosen
such that modes exiting the horizon during USR do not affect the unperturbed
solutions (via classical perturbations or stochastic kicks) until they freeze
after the USR phase.

$\displaystyle
3\tilde{H}^{2}-\frac{\tilde{H}^{2}(\partial_{\mathcal{N}}\phi)^{2}}{2}-V$
$\displaystyle=\mathcal{O}(\epsilon_{g}^{2})\,,$ (A.6)
$\displaystyle\tilde{H}\,\partial_{\mathcal{N}}(\tilde{H}\,\partial_{\mathcal{N}}\phi)+3\tilde{H}^{2}\partial_{\mathcal{N}}\phi+\partial_{\phi}V$
$\displaystyle=\mathcal{O}(\epsilon_{g}^{2})\,,$ (A.7)

where $\tilde{H}=-K/3$. This non-linear scheme to treat non-linear
cosmological perturbations of large wavelength is referred to as the _separate
universe approach_.

Integrating (A.5) and neglecting $\mathcal{O}(\epsilon_{g}^{2})$ terms, we get
[21] 282828Aside from being mathematically convenient for the argument, the
integral of (A.5) can be physically interpreted (up to
$\mathcal{O}(\epsilon_{g}^{2})$ corrections) as the expansion along a comoving
congruence of geodesics [21, 26].

$\psi(t_{2},\bm{x})-\psi(t_{1},\bm{x})=\mathcal{N}(t_{1}\to
t_{2},\bm{x})-N(t_{1}\to t_{2})+\mathcal{O}(\epsilon_{g}^{2})\,,$ (A.8)

where

$\displaystyle\mathcal{N}(t_{1}\to
t_{2},\bm{x})=\int_{t_{1}}^{t_{2}}\mathop{}\\!\mathrm{d}\mathcal{N}\,,$ (A.9)

and $N(t_{1}\to t_{2})=\int_{t_{1}}^{t_{2}}\mathop{}\\!\mathrm{d}t\,H$
corresponds to the number of $e$-folds among the $t=t_{1}$ and $t=t_{2}$
hypersurfaces in the unperturbed (FLRW) spacetime. No specific gauge has been
introduced so far. We now choose a gauge that interpolates between the _flat
gauge_ ($\partial^{i}h_{ij}=\psi=0$) hypersurface $\Sigma_{\text{f},t_{1}}$ at
$t=t_{1}$ and the _uniform density gauge_ ($\partial^{i}h_{ij}=\delta\rho=0$)
hypersurface $\Sigma_{\text{ud},t_{2}}$ at $t=t_{2}$. Applying this gauge
choice in (A.8), we get

$\left.\psi(t_{2},\bm{x})\right|_{\partial^{i}h_{ij}=\delta\rho=0}=\mathcal{N}(\Sigma_{\text{f},t_{1}}\to\Sigma_{\text{ud},t_{2}},\bm{x})-N(t_{1}\to
t_{2})+\mathcal{O}(\epsilon_{g}^{2})\,.$ (A.10)

The left-hand side is, by definition, the _non-linear perturbation_ 292929The
non-linear $\zeta$ is usually defined in this way [19, 26]. See, however, the
discussion in [36]. An alternative definition is presented in [21], which is
however _not_ gauge invariant and only reduces to the standard definition in
the perfectly adiabatic limit $p=p(\rho)$.  $\zeta(t_{2},x^{i})$. Hence,

$\zeta(t_{2};x^{i})=\mathcal{N}(\Sigma_{\text{f},t_{1}}\to\Sigma_{\text{ud},t_{2}},\bm{x})-N(t_{1}\to
t_{2})+\mathcal{O}(\epsilon_{g}^{2})\,.$ (A.11)

In order for (A.11) to be of practical use, we need a method to compute
$\mathcal{N}(\Sigma_{\text{f},t_{1}}\to\Sigma_{\text{ud},t_{2}};\bm{x})$,
which is provided by the separate universe approach.303030Notice that the
separate universe approach is applicable to (A.11) because both the flat and
uniform-density gauges imply $\partial^{i}h_{ij}=0$, which is necessary to
prove the separate universe approach, see [19]. Thanks to the explicit analogy
between eqs. (A.6) and (A.7) and the ones for the homogeneous $\phi$, the
dependence of the non-linear $\phi(\mathcal{N})$ and $\pi(\mathcal{N})$ on
$\mathcal{N}$ is the same as the unperturbed $\bar{\phi}(N)$ and
$\bar{\pi}(N)$ on $N$. This implies that the functional form of
$\mathcal{N}(\phi,\pi)$ is the same as that of $N(\bar{\phi},\bar{\pi})$. As a
consequence, $\mathcal{N}(\phi,\pi)$ and $N(\bar{\phi},\bar{\pi})$ may differ
_only_ due to their respective boundary conditions.

In Eq. (A.11), the boundary conditions are specified by the gauge choice at
$t=t_{1}$ and $t=t_{2}$. In particular, at $t=t_{1}$, we impose the flat
gauge, therefore
$\phi(t_{1},\bm{x})=\bar{\phi}(t_{1})+\delta\phi_{\text{f}}(t_{1},\bm{x})$,
$\pi(t_{1},\bm{x})=\bar{\pi}(t_{1})+\delta\pi_{\text{f}}(t_{1};\bm{x})$ where
$\delta\phi_{\text{f}}$ and $\delta\pi_{\text{f}}$ are the field and velocity
perturbations computed in the flat gauge. In practice, $\delta\phi_{\text{f}}$
and $\delta\pi_{\text{f}}$ are computed (approximately) by using linear
perturbation theory in the flat gauge. At $t=t_{2}$, we impose the uniform
density gauge, therefore $\delta\rho(t_{2},\bm{x})=0$ for all $\bm{x}$. The
energy density is a function of the field and its velocity, therefore
$\delta\rho[\delta\phi_{\text{ud}}(t_{2},\bm{x}),\delta\pi_{\text{ud}}(t_{2},\bm{x})]=0$.

As argued in Sec. 3, momentum perturbations in the initial hypersurface are
not negligible in general. Regarding the final hypersurface, the velocity
perturbations can be safely neglected.313131See [9] for a full calculation
involving field and momentum perturbations at both initial and final
hypersurfaces. This is because:

  1. 1.

The field has evolved long enough in the presence of a CR attractor for it to
effectively reach the unperturbed trajectory, i.e. the evolution of the field
dynamically erases $\delta\pi_{\text{ud}}(t_{2},\bm{x})$.

  2. 2.

The dependence on $\delta\pi_{\text{ud}}(t_{2},\bm{x})$ is through
$\delta\rho[\delta\phi_{\text{ud}}(t_{2},\bm{x}),\delta\pi_{\text{ud}}(t_{2},\bm{x})]=0$.
Along the unperturbed trajectory, this is suppressed by the slow-rolling of
the inflaton (the main contribution to the energy density is the potential,
i.e. the field value $\delta\phi_{\text{ud}}(t_{2},\bm{x})$).

This implies that the gauge condition in the last hypersurface (on top of
$\partial^{i}h_{ij}$) is

$\delta\rho[\delta\phi_{\text{ud}}(t_{2},\bm{x}),\delta\pi_{\text{ud}}(t_{2},\bm{x})]\simeq\delta\rho[\delta\phi_{\text{ud}}(t_{2},\bm{x})]=0\Rightarrow\delta\phi_{\text{ud}}(t_{2},\bm{x})=0\Rightarrow\phi(t_{2},\bm{x})=\bar{\phi}(t_{2})\,.$
(A.12)

As a consequence, the second hypersurface, which in rigour is a uniform-
density slice, can in practice be taken as a _uniform-field slice_.

Taking the above into account, the separate universe approach provides the
following prescription to compute
$\mathcal{N}(\Sigma_{{\text{f}},t_{1}}\to\Sigma_{\text{ud},t_{2}},\bm{x})$:
solve the _unperturbed_ equations of motion with initial conditions
$\bar{\phi}(t_{1})+\delta\phi_{\text{f}}$,
$\bar{\pi}(t_{1})+\delta\pi_{\text{f}}$, and then compute the number of
$e$-folds until $\phi=\bar{\phi}(t_{2})$ along said solution. The subtrahend
in (A.11) is similarly computed, but using $\bar{\phi}(t_{1})$ and
$\bar{\pi}(t_{1})$ as initial conditions. In summary,

$\displaystyle\zeta(t_{2},\bm{x})=\,$ $\displaystyle
N\left[\\{\bar{\phi}(t_{1})+\delta\phi_{\text{f}}(t_{1},\bm{x}),\bar{\pi}(t_{1})+\delta\pi_{\text{f}}(t_{1},\bm{x})\\}\to\bar{\phi}(t_{2})\right]$
$\displaystyle-N\left[\\{\bar{\phi}(t_{1}),\bar{\pi}(t_{1})\\}\to\bar{\phi}(t_{2})\right]\,,$
(A.13)

where $N$ denotes the functional form of $e$-folds as a function of field
value and velocity derived from the unperturbed equations of motion.

## References

  * [1] A. M. Green and B. J. Kavanagh, “Primordial Black Holes as a dark matter candidate,” J. Phys. G 48 no. 4, (2021) 043001, arXiv:2007.10722 [astro-ph.CO].
  * [2] G. Ballesteros and M. Taoso, “Primordial black hole dark matter from single field inflation,” Phys. Rev. D 97 no. 2, (2018) 023501, arXiv:1709.05565 [hep-ph].
  * [3] W. H. Kinney, “Horizon crossing and inflation with large eta,” Phys. Rev. D 72 (2005) 023515, arXiv:gr-qc/0503017.
  * [4] P. Ivanov, P. Naselsky, and I. Novikov, “Inflation and primordial black holes as dark matter,” Phys. Rev. D 50 (1994) 7173–7178.
  * [5] G. Ballesteros, J. Rey, M. Taoso, and A. Urbano, “Primordial black holes as dark matter and gravitational waves from single-field polynomial inflation,” JCAP 07 (2020) 025, arXiv:2001.08220 [astro-ph.CO].
  * [6] M. Taoso and A. Urbano, “Non-gaussianities for primordial black hole formation,” JCAP 08 (2021) 016, arXiv:2102.03610 [astro-ph.CO].
  * [7] J. M. Ezquiaga, J. García-Bellido, and V. Vennin, “The exponential tail of inflationary fluctuations: consequences for primordial black holes,” JCAP 03 (2020) 029, arXiv:1912.05399 [astro-ph.CO].
  * [8] D. G. Figueroa, S. Raatikainen, S. Rasanen, and E. Tomberg, “Implications of stochastic effects for primordial black hole production in ultra-slow-roll inflation,” JCAP 05 no. 05, (2022) 027, arXiv:2111.07437 [astro-ph.CO].
  * [9] S. Pi and M. Sasaki, “Logarithmic Duality of the Curvature Perturbation,” Phys. Rev. Lett. 131 no. 1, (2023) 011002, arXiv:2211.13932 [astro-ph.CO].
  * [10] V. Atal, J. Garriga, and A. Marcos-Caballero, “Primordial black hole formation with non-Gaussian curvature perturbations,” JCAP 09 (2019) 073, arXiv:1905.13202 [astro-ph.CO].
  * [11] J. M. Maldacena, “Non-Gaussian features of primordial fluctuations in single field inflationary models,” JHEP 05 (2003) 013, arXiv:astro-ph/0210603.
  * [12] A. A. Starobinsky, “STOCHASTIC DE SITTER (INFLATIONARY) STAGE IN THE EARLY UNIVERSE,” Lect. Notes Phys. 246 (1986) 107–126.
  * [13] E. Tomberg, “Stochastic constant-roll inflation and primordial black holes,” Phys. Rev. D 108 no. 4, (2023) 043502, arXiv:2304.10903 [astro-ph.CO].
  * [14] A. Karam, N. Koivunen, E. Tomberg, V. Vaskonen, and H. Veermäe, “Anatomy of single-field inflationary models for primordial black holes,” JCAP 03 (2023) 013, arXiv:2205.13540 [astro-ph.CO].
  * [15] D. Wands, “Duality invariance of cosmological perturbation spectra,” Phys. Rev. D 60 (1999) 023507, arXiv:gr-qc/9809062.
  * [16] D. Baumann, Cosmology.  Cambridge University Press, 7, 2022.
  * [17] G. Ballesteros, J. Beltran Jimenez, and M. Pieroni, “Black hole formation from a general quadratic action for inflationary primordial fluctuations,” JCAP 06 (2019) 016, arXiv:1811.03065 [astro-ph.CO].
  * [18] G. Ballesteros, J. Rey, M. Taoso, and A. Urbano, “Stochastic inflationary dynamics beyond slow-roll and consequences for primordial black hole formation,” JCAP 08 (2020) 043, arXiv:2006.14597 [astro-ph.CO].
  * [19] N. S. Sugiyama, E. Komatsu, and T. Futamase, “$\delta$N formalism,” Phys. Rev. D 87 no. 2, (2013) 023530, arXiv:1208.1073 [gr-qc].
  * [20] D. Wands, K. A. Malik, D. H. Lyth, and A. R. Liddle, “A New approach to the evolution of cosmological perturbations on large scales,” Phys. Rev. D 62 (2000) 043527, arXiv:astro-ph/0003278.
  * [21] D. H. Lyth, K. A. Malik, and M. Sasaki, “A General proof of the conservation of the curvature perturbation,” JCAP 05 (2005) 004, arXiv:astro-ph/0411220.
  * [22] V. Atal, J. Cid, A. Escrivà, and J. Garriga, “PBH in single field inflation: the effect of shape dispersion and non-Gaussianities,” JCAP 05 (2020) 022, arXiv:1908.11357 [astro-ph.CO].
  * [23] M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias, and A. Riotto, “The formation probability of primordial black holes,” Phys. Lett. B 820 (2021) 136602, arXiv:2105.07810 [astro-ph.CO].
  * [24] S. Hooshangi, M. H. Namjoo, and M. Noorbala, “Tail diversity from inflation,” JCAP 09 (2023) 023, arXiv:2305.19257 [astro-ph.CO].
  * [25] M. Celoria, P. Creminelli, G. Tambalo, and V. Yingcharoenrat, “Beyond perturbation theory in inflation,” JCAP 06 (2021) 051, arXiv:2103.09244 [hep-th].
  * [26] J. Garriga, Y. Urakawa, and F. Vernizzi, “$\delta N$ formalism from superpotential and holography,” JCAP 02 (2016) 036, arXiv:1509.07339 [hep-th].
  * [27] E. Tomberg, “Numerical stochastic inflation constrained by frozen noise,” JCAP 04 (2023) 042, arXiv:2210.17441 [astro-ph.CO].
  * [28] A. A. Starobinsky, “Stochastic de Sitter (inflationary) Stage in the Early Universe,” in Field Theory, Quantum Gravity and Strings, H. J. de Vega and N. Sánchez, eds., vol. 246, p. 107.  1986\.
  * [29] V. Vennin, Stochastic inflation and primordial black holes.  PhD thesis, U. Paris-Saclay, 6, 2020.  arXiv:2009.08715 [astro-ph.CO].
  * [30] J. H. P. Jackson, H. Assadullahi, A. D. Gow, K. Koyama, V. Vennin, and D. Wands, “The separate-universe approach and sudden transitions during inflation,” arXiv:2311.03281 [astro-ph.CO].
  * [31] D. G. Figueroa, S. Raatikainen, S. Rasanen, and E. Tomberg, “Non-Gaussian Tail of the Curvature Perturbation in Stochastic Ultraslow-Roll Inflation: Implications for Primordial Black Hole Production,” Phys. Rev. Lett. 127 no. 10, (2021) 101302, arXiv:2012.06551 [astro-ph.CO].
  * [32] J. Grain and V. Vennin, “Canonical transformations and squeezing formalism in cosmology,” JCAP 02 (2020) 022, arXiv:1910.01916 [astro-ph.CO].
  * [33] T. Hamazaki, “Long wavelength limit of evolution of nonlinear cosmological perturbations,” Phys. Rev. D 78 (2008) 103513, arXiv:0811.2366 [astro-ph].
  * [34] D. S. Salopek and J. R. Bond, “Nonlinear evolution of long wavelength metric fluctuations in inflationary models,” Phys. Rev. D 42 (1990) 3936–3962.
  * [35] A. E. Romano, S. Mooij, and M. Sasaki, “Adiabaticity and gravity theory independent conservation laws for cosmological perturbations,” Phys. Lett. B 755 (2016) 464–468, arXiv:1512.05757 [gr-qc].
  * [36] D. Cruces, C. Germani, and A. Palomares, “An update on adiabatic modes in cosmology and $\delta$N formalism,” JCAP 06 (2023) 002, arXiv:2212.05112 [astro-ph.CO].