# Primordial black holes and stochastic inflation beyond slow roll: I - noise
matrix elements

Swagat S. Mishra  Edmund J. Copeland  and Anne M. Green

###### Abstract

Primordial Black Holes (PBHs) may form in the early Universe, from the
gravitational collapse of large density perturbations, generated by large
quantum fluctuations during inflation. Since PBHs form from rare over-
densities, their abundance is sensitive to the tail of the primordial
probability distribution function (PDF) of the perturbations. It is therefore
important to calculate the full PDF of the perturbations, which can be done
non-perturbatively using the ‘stochastic inflation’ framework. In single field
inflation models generating large enough perturbations to produce an
interesting abundance of PBHs requires violation of slow roll. It is therefore
necessary to extend the stochastic inflation formalism beyond slow roll. A
crucial ingredient for this are the stochastic noise matrix elements of the
inflaton potential. We carry out analytical and numerical calculations of
these matrix elements for a potential with a feature which violates slow roll
and produces large, potentially PBH generating, perturbations. We find that
the transition to an ultra slow-roll phase results in the momentum induced
noise terms becoming larger than the field noise whilst each of them falls
exponentially for a few e-folds. The noise terms then start rising with their
original order restored, before approaching constant values which depend on
the nature of the slow roll parameters in the post transition epoch. This will
significantly impact the quantum diffusion of the coarse-grained inflaton
field, and hence the PDF of the perturbations and the PBH mass fraction.

##  1 Introduction

A multitude of cosmological and astrophysical observations indicate that most
of the matter in the Universe is non-baryonic cold dark matter [1, 2, 3].
Primordial black holes (PBHs), black holes that might have been formed from
over-densities in the early Universe [4, 5, 6, 7, 8], are a potential dark
matter candidate [5, 9, 11, 10, 12, 13, 14]. There has been a surge of
interest in PBHs in recent years, both their formation and observational
probes, following the detection of gravitational waves from binary black hole
mergers [15, 16, 17, 18].

‘Cosmic inflation’ (a period of accelerated expansion) has emerged as the
leading scenario for the very early Universe, prior to the commencement of the
radiation-dominated hot Big Bang [19, 20, 21, 22, 23, 24]. A period of at
least $60$-$70$ e-folds of inflation generates natural initial conditions [25,
26, 27, 28, 29, 30]. Furthermore quantum fluctuations of the inflaton field
can generate the density perturbations from which structure forms.
Observations of the anisotropies in the Cosmic Microwave Background (CMB)
radiation [31, 32] provide strong evidence that structure formation on
cosmological scales is seeded by almost scale-invariant, nearly Gaussian,
adiabatic initial density fluctuations, consistent with the predictions of the
simplest single field slow-roll inflation scenario [33, 24].

CMB observations [34, 35] are consistent with the inflaton field, $\phi$,
rolling slowly down an asymptotically flat potential $V(\phi)$ during the
epoch when cosmological scales exit the Hubble radius, $50$-$60$ e-foldings
before the end of inflation. However the scales probed by CMB and large scale
structure (LSS) observations correspond to only $7$-$8$ e-folds of inflation,
and hence a relatively small region of the inflaton potential. On smaller
scales, deviations from slow roll may lead to interesting changes in the
primordial perturbations. In particular, if the scalar perturbations are
sufficiently large on small scales, then PBHs may form when these modes
reenter the Hubble radius during the post-inflationary epoch. PBHs are
therefore a powerful probe of the inflaton potential over the full range of
field values. Large, PBH forming, fluctuations can be generated by a feature
in the inflationary potential, such as a flat inflection point (see Fig. 1).
Such a feature can substantially slow down the already slowly rolling inflaton
field, causing the inflaton to enter into an ultra slow-roll (USR) phase,
which leads to an enhancement of the power spectrum, ${\cal P}_{\zeta}$, of
the primordial curvature perturbation, $\zeta$.

Figure 1: A schematic illustration of a plateau potential (solid green line).
The ‘CMB Window’ represents field values corresponding to cosmological scales
$k\in\left[0.0005,0.5\right]\leavevmode\nobreak\ {\rm Mpc}^{-1}$ that are
probed by CMB observations. The blue star represents the CMB pivot scale
$k_{*}=0.05\,{\rm Mpc}^{-1}$. The potential has a flat inflection-point like
segment (highlighted with pink shading) which results in ultra slow-roll (USR)
inflation. After the first slow-roll phase (SR-I) near the CMB Window, the
inflaton enters into an USR phase. During this transient phase of USR, the
second slow-roll condition (see Eq. (2.7)) is violated, specifically
$\eta_{H}\simeq+3$. This leads to an enhancement in the primordial
perturbations on small scales. Later, the inflaton emerges from the USR into
another slow-roll phase (SR-II) before inflation ends at $\phi_{\rm end}$.

There are several subtleties in calculating the abundance of PBHs formed from
inflation models with a feature in the potential. Firstly, the sharp drop in
the classical drift speed of the inflaton means that the effects of stochastic
quantum diffusion on its motion become non-negligible, and potentially even
dominant. Even more importantly, since PBHs form from the rare extreme peaks
of curvature fluctuations111Note that while we will be working with the PDF of
curvature fluctuations $P[\zeta]$, the criterion for PBH formation is best
specified in terms of the density contrast, $\delta$. Hence, an accurate
computation of the PBH abundance requires the multivariate joint probabilities
of the curvature perturbations [40, 39, 36, 37, 38]. their mass fraction at
formation, $\beta_{\rm PBH}$, is sensitive to the tail of the probability
distribution function (PDF) $P[\zeta]$ of the primordial fluctuations.
Consequently, perturbative computations using the power-spectrum lead to an
inaccurate estimation of the PBH mass fraction. Hence the calculation of the
full primordial PDF, which can be done non-perturbatively using the
‘Stochastic Inflation’ (SI) framework [27, 41, 43, 42, 44, 45, 46, 47], is
extremely important (see Refs. [48, 49, 50, 52, 53, 54, 36, 37, 51]).

There has been a profusion of work on SI in the recent literature in the
context of PBH formation [55, 56, 59, 57, 58, 60, 62, 63, 61, 64, 65, 66, 67,
68, 69, 70, 71, 72, 73], and it has been shown that SI generically predicts a
(highly non-Gaussian) exponential tail in the PDF [55, 58]. However, since PBH
formation usually requires slow-roll violation, it is important to develop the
SI formalism beyond slow-roll [56, 60, 59, 61, 65]. In this context, the
validity of various assumptions of the SI formalism have been scrutinised [60,
66, 67].

The stochastic inflation formalism is an effective treatment of the dynamics
of the long-wavelength (IR) part of the inflaton field coarse-grained on
scales much greater than the Hubble radius i.e. $k\leq\sigma\,aH$, with the
constant $\sigma\ll 1$. In this framework, the evolution of the coarse-grained
inflaton field is governed by two first-order non-linear classical stochastic
differential equations (Langevin equations) which receive constant quantum
kicks from the small scale UV modes that are exiting the Hubble radius due to
the accelerated expansion during inflation. Hence the small-scale fluctuations
constitute classical stochastic noise terms in the Langevin equations denoted
by $\Sigma_{\phi\phi}$, $\Sigma_{\pi\pi}$, and $\Sigma_{\phi\pi}$
corresponding to the inflaton field noise, momentum noise, and the cross-noise
terms (defined in Eq. (3.15)). The SI formalism is generally combined with the
classical $\delta N$ formalism [27, 74, 75, 76, 77] in order to compute
cosmological correlators in this framework. This leads to the emergence of the
stochastic $\delta{\cal N}$ formalism222Note that $N\propto{\rm ln}(a)$, the
number of e-folds during inflation, is a deterministic variable, while ${\cal
N}$ is a stochastic variable as defined in Sec. 3.1. [45, 46, 47, 55].

The PDF $P[\zeta]$ of the primordial curvature perturbation can then be
determined by using the techniques of first-passage time analysis for the
stochastic distribution of the number of e-folds $\cal N$ with fixed boundary
conditions on the coarse-grained inflaton field. A convenient analytic
approach to obtain the distribution of first-passage e-folds (and hence the
PDF $P[\zeta]$) is to solve the corresponding Fokker-Planck equation (FPE)
with the same boundary conditions [55, 58, 68]. In the analysis of stochastic
dynamics, the noise terms in the FPE are assumed [47, 55, 58, 68] to be of de
Sitter-type, i.e. $\Sigma_{\pi\pi},\Sigma_{\phi\pi}\simeq 0$, and the field
noise, $\Sigma_{\phi\phi}=(H/(2\pi))^{2}$, is constant (where $H$ is the
Hubble expansion rate during the de-Sitter type phase). However, in the
context of PBH formation, since slow roll is usually violated, it is important
to compute the stochastic noise matrix elements more accurately, otherwise the
determination of the PDF becomes inaccurate, which in turn leads to an
imprecise estimation of the PBH mass fraction $\beta_{\rm PBH}$.

Our aim is to develop analytical and semi-analytical techniques to estimate
the full PDF $P[\zeta]$ using the stochastic inflation formalism beyond slow-
roll. In the present paper, we carry out a thorough analytical and numerical
computation of the stochastic noise-matrix elements accurately beyond the
slow-roll approximations. In a forthcoming paper [78] we determine the PDF
$P[\zeta]$ by solving the Fokker-Planck equation beyond slow-roll with
appropriate noise matrix elements and discuss the implications for estimating
the mass fraction of PBHs.

In what follows, we begin with a brief introduction to the classical
inflationary dynamics in Sec. 2, with particular focus on the ultra slow-roll
dynamics across a flat segment in the inflaton potential. In Sec. 3, we
discuss the quantum dynamics in the stochastic inflation framework. We
introduce the Langevin equation in Sec. 3.1 and emphasise the importance of
the noise matrix elements in the adjoint Fokker-Planck equation in Sec. 3.2.
Section 4 is dedicated to the computation of the noise matrix elements which
is the primary focus of this work. We numerically compute the noise matrix
elements for a slow-roll potential as well as a potential with a slow-roll
violating feature in Sec. 4.2.1 before proceeding to carry out a thorough
analytical treatment in Sec. 4.2.2 for instantaneous transitions between
different phases during inflation. We discuss the potential implications of
our results for the computation of the PBH mass fraction and spell out a
number of complexities associated with the computation in Sec. 5 before
concluding with a summary of our main results in Sec. 6. Appendix A provides a
derivation of the Mukhanov-Sasaki equation in spatially flat gauge. Appendix B
deals with the analytical solutions of the Mukhanov-Sasaki equation in the
absence of any transition, while Appendix C provides analytical expressions
for the noise matrix elements in the super-Hubble limit. Appendices D and E
are dedicated to the dynamics during instantaneous transitions.

We work in natural units with $c=\hbar=1$ and define the reduced Planck mass
to be $m_{\rm p}\equiv 1/\sqrt{8\pi G}=2.43\times 10^{18}\leavevmode\nobreak\
{\rm GeV}$. We assume the background Universe to be described by a spatially
flat Friedmann-Lemaitre-Robertson-Walker (FLRW) metric with signature
$(-,+,+,+)$. An overdot (.) denotes derivative with respect to cosmic time
$t$, while an overdash ${}^{(^{\prime})}$ denotes derivative with respect to
the conformal time $\tau$.

##  2 Inflationary dynamics beyond slow roll

We focus on the inflationary scenario of a single canonical scalar field
$\phi$ with a self-interaction potential $V(\phi)$ which is minimally coupled
to gravity. The system is described by the action

$S[g_{\mu\nu},\phi]=\int\,{\rm
d}^{4}x\,\sqrt{-g}\,\left[\,\frac{m_{p}^{2}}{2}\,R-\frac{1}{2}\,\partial_{\mu}\phi\,\partial_{\nu}\phi\,g^{\mu\nu}-V(\phi)\right]\leavevmode\nobreak\
,$ (2.1)

where $R$ is the Ricci scalar and $g_{\mu\nu}$ is the metric tensor.
Specializing to the spatially flat FLRW background metric

${\rm d}s^{2}=-{\rm d}t^{2}+a^{2}(t)\left[{\rm d}x^{2}+{\rm d}y^{2}+{\rm
d}z^{2}\right]\,,$ (2.2)

the evolution equations for the scale factor, $a(t)$, and inflaton, $\phi(t)$,
are

$\displaystyle H^{2}\equiv\frac{1}{3m_{p}^{2}}\,\rho_{\phi}$ $\displaystyle=$
$\displaystyle\frac{1}{3m_{p}^{2}}\left[\frac{1}{2}{\dot{\phi}}^{2}+V(\phi)\right]\,,$
(2.3) $\displaystyle\dot{H}\equiv\frac{\ddot{a}}{a}-H^{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2m_{p}^{2}}\,\dot{\phi}^{2}\,,$ (2.4)
$\displaystyle{\ddot{\phi}}+3\,H{\dot{\phi}}+V_{,\phi}(\phi)$ $\displaystyle=$
$\displaystyle 0\,.$ (2.5)

where $H\equiv\dot{a}/a$ and $V_{,\phi}\equiv dV/d\phi$.

The slow-roll regime of inflation is usually characterised by the first two
kinematic Hubble slow-roll parameters $\epsilon_{H}$ and $\eta_{H}$, defined
by

$\displaystyle\epsilon_{H}$ $\displaystyle=$
$\displaystyle-\frac{\dot{H}}{H^{2}}=\frac{1}{2m_{p}^{2}}\,\frac{\dot{\phi}^{2}}{H^{2}}\,,$
(2.6) $\displaystyle\eta_{H}$ $\displaystyle=$
$\displaystyle-\frac{\ddot{\phi}}{H\dot{\phi}}=\epsilon_{H}+\frac{1}{2\epsilon_{H}}\,\frac{{\rm
d}\epsilon_{H}}{{\rm d}N}\,,$ (2.7)

where $N={\rm ln}(a/a_{i})$ is the number of e-folds of expansion during
inflation with $a_{i}$ the initial scale factor at some early epoch during
inflation before the Hubble-exit of CMB scale fluctuations. The slow-roll
conditions correspond to

$\epsilon_{H},\leavevmode\nobreak\ \eta_{H}\ll 1\leavevmode\nobreak\ .$ (2.8)

It follows from the definition of the Hubble parameter, $H$, and
$\epsilon_{H}$ in Eq. (2.6), that the condition for the Universe to
accelerate, ${\ddot{a}}>0$, is $\epsilon_{H}<1$. Before proceeding further, we
remind the reader of the distinction between the quasi-de Sitter (qdS) and
slow-roll (SR) approximations.

  * •

Quasi-de Sitter inflation corresponds to the condition $\epsilon_{H}\ll 1$.

  * •

Slow-roll inflation corresponds to both $\epsilon_{H},\,\eta_{H}\ll 1$.

Hence, one can deviate from the slow-roll regime by having $|\eta_{H}|\geq 1$
while still maintaining the qdS expansion by keeping $\epsilon_{H}\ll 1$,
which is exactly what happens during ultra slow-roll (USR) inflation. This
distinction will be important for the rest of this paper. Under either of the
aforementioned assumptions, the conformal time, $\tau$, is given by

$\displaystyle-\tau\simeq\frac{1}{aH}\,.$ (2.9)

As discussed in Sec. 1, in order to facilitate PBH formation, we need to
significantly amplify the scalar power at small-scales. This can be achieved
with a feature in the inflaton potential, such as an inflection point-like
feature (as shown in Fig. 1) for which $V_{,\phi}\ll 3H\dot{\phi}$. The
following criteria need to be satisfied for an inflationary potential to be
compatible with observations on cosmological scales [34] while also generating
perturbations on smaller scales that are large enough to form an interesting
abundance of PBHs:

  * •

At the CMB pivot scale, $k_{*}=(aH)_{*}=0.05\leavevmode\nobreak\ {\rm
Mpc}^{-1}$, the amplitude of the scalar power spectrum is

${\cal P}_{\zeta}(k_{*})=2.1\times 10^{-9}\,,$ (2.10)

with the scalar spectral index $n_{{}_{S}}$ and tensor-to-scalar ratio $r$
satisfying

$n_{{}_{S}}(k_{*})\in\left[0.957,0.975\right]\leavevmode\nobreak\ ,\quad
r(k_{*})\leq 0.036\leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
at}\leavevmode\nobreak\ \leavevmode\nobreak\ 95\%\leavevmode\nobreak\
\mbox{C.L}\,.$ (2.11)

  * •

A feature in $V(\phi)$ on a smaller scale $k\gg k_{*}$ (closer to the end of
inflation $N_{e}<N_{*}$) to enhance the primordial scalar power spectrum by a
factor of roughly $10^{7}$ with respect to its value at the CMB pivot scale.
Here $N_{e}$ is the number of e-folds before the end of inflation and $N_{*}$
is the value of $N_{e}$ when the CMB pivot scale made its Hubble-exit.
Typically $N_{*}\in[50,60]$ depending on the reheating history of the Universe
(see e.g. Ref. [79]). Throughout this work we take $N_{*}=60$.

  * •

The potential steepens, so that inflation ends. Reheating (and the transition
to the subsequent radiation dominated epoch) then occurs as the field
oscillates around a minimum in the potential.

Given that PBH formation requires the enhancement of the inflationary power
spectrum by a factor of $10^{7}$ within less than 40 e-folds of expansion (see
e.g. Ref. [80]), the quantity $\Delta\ln{\epsilon_{H}}/\Delta N$, and hence
$|\eta_{H}|$, must grow to become of order unity, thereby violating the second
slow-roll condition in Eq. (2.8). In the particular case of a flat plateau
region in the potential at intermediate field values, the inflaton enters a
transient period of ultra slow-roll (USR) (see Refs. [81, 82, 55, 62, 68]).
Since $V_{,\phi}=0$, from the equation of motion Eq. (2.5),
$\ddot{\phi}+3H\dot{\phi}=0\Rightarrow-\ddot{\phi}/(H\,\dot{\phi})=+3$,
leading to

$\eta_{H}=+3\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ ({\rm during\leavevmode\nobreak\
USR}).$ (2.12)

As a consequence, the inflaton speed drops exponentially with the number of
e-folds during this USR phase:

$\dot{\phi}=\dot{\phi}_{\rm en}\,e^{-3\,H\,(t-t_{\rm en})}\propto
e^{-3\,N}\leavevmode\nobreak\ ,$ (2.13)

where $\dot{\phi}_{\rm en}$ is the entry velocity to the USR phase at time
$t_{\rm en}$.

Figure 2: A zoomed-in version of Fig. 1 in order to schematically illustrate
the intermediate flat quantum well feature (highlighted with pink shading) in
the inflaton potential. The height and width of the flat segment are denoted
by $V_{\rm well}$ and $\Delta\phi_{\rm well}$ respectively. After exiting the
first slow-roll phase (SR-I) near the CMB window, the inflaton enters the flat
region at $\phi=\phi_{\rm en}$ at intermediate field values. During this USR
phase, the effects of quantum diffusion might become significant and hence one
should use the stochastic inflation formalism to compute the primordial PDF of
$\zeta$. Later, the inflaton emerges from the USR phase to another slow-roll
phase (SR-II) at $\phi=\phi_{\rm ex}$, before the end of inflation.

Since USR is a transient non-attractor phase, the inflaton dynamics during
this phase are sensitive to the initial conditions, in particular to the speed
$\dot{\phi}_{\rm en}$ with which the inflaton enters the plateau. In this
context, the inflaton potential exhibits three important regimes, namely, the
slow-roll SR-I phase for $\phi>\phi_{\rm en}$ around the CMB scale, the USR
phase at some intermediate field values $\phi_{\rm ex}\leq\phi\leq\phi_{\rm
en}$, succeeded by the final SR-II phase for $\phi<\phi_{\rm ex}$ before the
end of inflation at $\phi=\phi_{\rm end}$. Figure 2 schematically illustrates
the three regimes. The flat regime (flat quantum well)333Note that we refer to
the flat USR regime as the ‘flat quantum well’ because the inflaton dynamics
are usually dominated by stochastic quantum diffusion, as discussed in the
subsequent Sections. is characterised by its width $\Delta\phi_{\rm
well}=\phi_{\rm en}-\phi_{\rm ex}$, and height, $V_{\rm well}$. During this
regime

$\dot{\phi}=\dot{\phi}_{\rm en}-3H\left(\phi-\phi_{\rm
en}\right)\leavevmode\nobreak\ .$ (2.14)

The total number of e-folds of expansion during the USR period up to $\phi$,
where $\phi_{\rm ex}\leq\phi\leq\phi_{\rm en}$ is given by

$N_{\rm USR}(\phi)=\frac{1}{3}\,{\rm log}\left(\frac{\pi_{\rm
en}}{\pi}\right)=\frac{1}{3}\,{\rm log}\left(\frac{\pi_{\rm en}}{\pi_{\rm
en}-3\,\left(\phi-\phi_{\rm en}\right)}\right)\leavevmode\nobreak\ ,$ (2.15)

where

$\pi=\frac{{\rm d}\phi}{{\rm d}N}=\frac{\dot{\phi}}{H}\,.$ (2.16)

In order to amplify the perturbations sufficiently to generate an interesting
abundance of PBHs, the USR phase typically has to last for around $2-3$
e-folds (see Refs. [83, 84]). The above expression for $N_{\rm USR}(\phi)$ can
be used in the ‘non-linear classical $\delta$N formalism’ to determine the PDF
of primordial fluctuations [50].

From the above expressions, it is clear that the dynamics of inflation during
USR is sensitive to the initial conditions $\\{\phi_{\rm en},\pi_{\rm en}\\}$.
Let us define the critical entry velocity $\dot{\phi}_{\rm cr}$ to be the
speed at which the inflaton must enter the flat quantum well in order to come
to a halt at $\phi_{\rm ex}$. From Eqs. (2.14) and (2.16) it follows that

$\dot{\phi}_{\rm cr}=-3\,H\,\Delta\phi_{\rm well}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \pi_{\rm
cr}=-3\,\Delta\phi_{\rm well}\,.$ (2.17)

If $\dot{\phi}_{\rm en}>\dot{\phi}_{\rm cr}$, then the classical speed of the
inflaton is large enough to drive it all the way across the quantum well,
while for $\dot{\phi}_{\rm en}<\dot{\phi}_{\rm cr}$, the inflaton comes to a
halt at some intermediate point $\phi\in\left(\phi_{\rm ex},\phi_{\rm
en}\right)$. Another important constraint comes from requiring inflation to
continue, $\epsilon_{H}<1$, hence from Eq. (2.6) we get

$0\leq\dot{\phi}_{\rm en}<\sqrt{2}\,H\,m_{p}\,.$ (2.18)

In this Section, we have discussed the classical dynamics of the inflaton
field beyond slow roll, with the specific example of ultra slow-roll inflation
across a flat potential well. We now move on to describe the large-scale
quantum dynamics of the inflaton field which is coarse-grained over super-
Hubble scales, using the stochastic inflation formalism. This will enable us
to study the PDF of the primordial fluctuations generated by the quantum
diffusion of the inflaton.

##  3 Quantum dynamics: stochastic inflation formalism

Stochastic inflation is an effective long wavelength IR treatment of inflation
in which the inflaton field is coarse-grained over super-Hubble scales
$k\leq\sigma\,aH$, with the constant $\sigma\ll 1$. On the other hand, the
Hubble-exiting smaller scale UV modes are constantly converted into IR modes
due to the accelerated expansion during inflation. Hence the coarse-grained
inflaton field follows a Langevin-type stochastic differential equation
featuring classical stochastic noise terms sourced by the smaller scale UV
modes, on top of the classical drift terms sourced by the gradient of the
self-interaction potential $V_{,\phi}(\phi)$.

We start with the Hamiltonian equations [85] of the system, Eq. (2.1), for
Heisenberg operators of the inflaton $\hat{\phi}$ and its momentum
$\hat{\pi}_{\phi}$

$\displaystyle\frac{{\rm d}\hat{\phi}}{{\rm d}N}$ $\displaystyle=$
$\displaystyle\hat{\pi}_{\phi}\leavevmode\nobreak\ ,$ (3.1)
$\displaystyle\frac{{\rm d}\hat{\pi}_{\phi}}{{\rm d}N}$ $\displaystyle=$
$\displaystyle-\left(3-\epsilon_{H}\right)\hat{\pi}_{\phi}-\frac{V_{,\phi}}{H^{2}}\leavevmode\nobreak\
,$ (3.2)

where we choose the number of e-folds $N$ as our time evolution variable for
$\phi(N,\vec{x})$ and $\pi_{\phi}(N,\vec{x})$ following Refs. [47, 55].

We split the inflaton $\hat{\phi}(N,\vec{x})$ and its conjugate momentum
$\hat{\pi}_{\phi}(N,\vec{x})$ into the corresponding IR
$\\{\hat{\Phi},\hat{\Pi}\\}$ and UV $\\{\hat{\varphi},\hat{\pi}\\}$ parts:

$\hat{\phi}=\hat{\Phi}+\hat{\varphi}\leavevmode\nobreak\ \leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\
\hat{\pi}_{\phi}=\hat{\Pi}+\hat{\pi}\leavevmode\nobreak\ ,$ (3.3)

where the UV fields are defined as

$\displaystyle\hat{\varphi}(N,\vec{x})$ $\displaystyle=$
$\displaystyle\int\,\frac{{\rm
d}^{3}\vec{k}}{\left(2\pi\right)^{\frac{3}{2}}}\,W\left(\frac{k}{\sigma
aH}\right)\,\left[\phi_{k}(N)\,\hat{a}_{\vec{k}}\,e^{-i\vec{k}.\vec{x}}+\phi_{k}^{*}(N)\,\hat{a}_{\vec{k}}^{\dagger}\,e^{i\vec{k}.\vec{x}}\right]\,,$
(3.4) $\displaystyle\hat{\pi}(N,\vec{x})$ $\displaystyle=$
$\displaystyle\int\,\frac{{\rm
d}^{3}\vec{k}}{\left(2\pi\right)^{\frac{3}{2}}}\,W\left(\frac{k}{\sigma
aH}\right)\,\left[\pi_{k}(N)\,\hat{a}_{\vec{k}}\,e^{-i\vec{k}.\vec{x}}+\pi_{k}^{*}(N)\,\hat{a}_{\vec{k}}^{\dagger}\,e^{i\vec{k}.\vec{x}}\right]\,.$
(3.5)

Here $W\left(k/\sigma aH\right)$ is the ‘window function’ that selects out
modes with $k>\sigma aH$. This coarse-graining guarantees that the IR fields
$\\{\hat{\Phi},\hat{\Pi}\\}$ are composed of modes of super-Hubble wavelengths
and hence can be treated as classical (stochastic) variables. Substituting the
expressions for the UV fields from Eqs. (3.4) and (3.5) into Eqs. (3.1) and
(3.2), the equations for the coarse-grained fields are [68]

$\displaystyle\frac{{\rm d}\hat{\Phi}}{{\rm d}N}$ $\displaystyle=$
$\displaystyle\hat{\Pi}+\hat{\xi}_{\phi}(N)\leavevmode\nobreak\ ,$ (3.6)
$\displaystyle\frac{{\rm d}\hat{\Pi}}{{\rm d}N}$ $\displaystyle=$
$\displaystyle-\left(3-\epsilon_{H}\right)\hat{\Pi}-\frac{V_{,\phi}(\hat{\Phi})}{H^{2}}+\hat{\xi}_{\pi}(N)\leavevmode\nobreak\
,$ (3.7)

where the field and momentum noise operators $\hat{\xi}_{\phi}(N)$ and
$\hat{\xi}_{\pi}(N)$ are given by

$\displaystyle\hat{\xi}_{\phi}(N)$ $\displaystyle=$
$\displaystyle-\int\,\frac{{\rm
d}^{3}\vec{k}}{\left(2\pi\right)^{\frac{3}{2}}}\,\frac{{d}}{{\rm
d}N}W\left(\frac{k}{\sigma
aH}\right)\,\left[\phi_{k}(N)\,\hat{a}_{\vec{k}}\,e^{-i\vec{k}.\vec{x}}+\phi_{k}^{*}(N)\,\hat{a}_{\vec{k}}^{\dagger}\,e^{i\vec{k}.\vec{x}}\right]\leavevmode\nobreak\
,$ (3.8) $\displaystyle\hat{\xi}_{\pi}(N)$ $\displaystyle=$
$\displaystyle-\int\,\frac{{\rm
d}^{3}\vec{k}}{\left(2\pi\right)^{\frac{3}{2}}}\,\frac{{d}}{{\rm
d}N}W\left(\frac{k}{\sigma
aH}\right)\,\left[\pi_{k}(N)\,\hat{a}_{\vec{k}}\,e^{-i\vec{k}.\vec{x}}+\pi_{k}^{*}(N)\,\hat{a}_{\vec{k}}^{\dagger}\,e^{i\vec{k}.\vec{x}}\right]\leavevmode\nobreak\
.$ (3.9)

We assume a window function which imposes a sharp cut off444This assumption
has been questioned, especially because the sharp cut-off window function may
not lead to well-behaved coarse-grained field correlators in the physical
space. See Refs. [86, 87, 72] for more discussion. between the IR and UV
momentum space modes:

$W\left(\frac{k}{\sigma aH}\right)=\Theta\left(\frac{k}{\sigma
aH}-1\right)\,.$ (3.10)

It has the advantage of making the calculation of the noise correlation matrix
elements more tractable.

Physically, the noise terms $\hat{\xi}_{\phi}$ and $\hat{\xi}_{\pi}$ in the
Langevin Eqs. (3.6) and (3.9) are sourced by the constant outflow of UV modes
into the IR modes, i.e. as a UV mode exits the cut-off scale $k=\sigma aH$ to
become part of the IR field on super-Hubble scales, the IR field receives a
‘quantum kick’ whose typical amplitude is given by $\sim\sqrt{\langle
0|\hat{\xi}(N)\hat{\xi}(N^{\prime})|0\rangle}$, where $|0\rangle$ is usually
taken to be the Bunch-Davies vacuum. Given that $\sigma\ll 1$, this happens on
ultra super-Hubble scales, where the UV modes must have already become
classical fluctuations555While the quantum-to-classical transition is still an
open problem, the treatment of UV noise operators as stochastic noise terms is
ensured to be valid as long as the constant mode of $\phi_{k}$ is dominant
over the decaying mode on super-Hubble scales [60]. due to the rapid decline
of the non-commuting parts of the fields $\\{\phi_{k},\pi_{k}\\}$ outside the
Hubble radius [88, 89, 90]. This leads to the classical stochastic description
of the dynamics of the coarse-grained quantum fields $\hat{\Phi}$, $\hat{\Pi}$
as discussed in the following subsection(s).

###  3.1 Langevin equation

The Langevin equations corresponding to Eqs. (3.6) and (3.7) take the compact
form

$\frac{{\rm d}\Phi_{i}}{{\rm d}N}=D_{i}+\xi_{i}(N)\,$ (3.11)

with coarse-grained IR variables $\Phi_{i}=\\{\Phi,\Pi\\}$ and the drift terms

$D_{i}=\left\\{\Pi,-\left(3-\epsilon_{H}\right)\Pi-\frac{V_{,\phi}(\Phi)}{H^{2}}\right\\},$
(3.12)

along with the noise operator terms $\xi_{i}=\\{\xi_{\phi},\xi_{\pi}\\}$.

In this compact notation the expressions for the noise operators, Eqs. (3.8)
and (3.9), become

$\hat{\xi}_{i}(N)=-\int\,\frac{{\rm
d}^{3}\vec{k}}{\left(2\pi\right)^{\frac{3}{2}}}\,\frac{{d}}{{\rm
d}N}W\left(\frac{k}{\sigma
aH}\right)\,\left[\phi_{i_{k}}(N)\,\hat{a}_{\vec{k}}\,e^{-i\vec{k}.\vec{x}}+\phi_{i_{k}}^{*}(N)\,\hat{a}_{\vec{k}}^{\dagger}\,e^{i\vec{k}.\vec{x}}\right]\,,$
(3.13)

with $\phi_{i_{k}}=\\{\phi_{k},\pi_{k}\\}$ being the field and momentum mode
functions respectively. Assuming the sharp $k$-space window function, Eq.
(3.10), it is easy to show that the equal-space noise correlators (auto-
correlators) take the form [85]

$\langle\xi_{i}(N)\,\xi_{j}(N^{\prime})\rangle=\Sigma_{ij}(N)\,\delta_{D}(N-N^{\prime})\,,$
(3.14)

where the noise correlation matrix $\Sigma_{ij}$ has the form

$\Sigma_{ij}(N)=(1-\epsilon_{H})\,\frac{k^{3}}{2\pi^{2}}\,\phi_{i_{k}}(N)\phi_{j_{k}}^{*}(N)\bigg{|}_{k=\sigma
aH}\leavevmode\nobreak\ .$ (3.15)

The stochastic nature of the noise leads to a probabilistic description of the
system $\Phi_{i}=\\{\Phi,\Pi\\}$. One approach to analyse the system is by
solving the Langevin equation, Eq. (3.11), numerically for many (tens of
millions) stochastic realizations and then proceeding to compute different
moments of the physical (stochastic) variables. This method is direct, however
cumbersome, non-analytical and requires significant computational power. See
Refs. [45, 46] for some of the earlier attempts in this direction, while for a
more concrete analysis beyond slow-roll, see Ref. [64], and for state-of-the-
art numerical simulations, relevant for determining the PDF of primordial
fluctuations, see Refs. [65, 69, 70, 73].

There is also an analytically concrete way to study this system, using the
first-passage time analysis which involves making a transition from the
Langevin equations to an equivalent second order partial differential Fokker-
Planck equation (FPE) [91, 92, 41, 44], that describes the time evolution of
the PDF of the stochastic variables $\\{\Phi,\Pi\\}$, subject to appropriate
boundary conditions. Given our primary goal of computing the full PDF
$P[\zeta]$, we take this route following Refs. [55, 58, 68].

The FPE corresponding to the Langevin equation, Eq. (3.11), takes the form

$\frac{\partial}{\partial N}\,P(\Phi_{i};N)={\cal L}_{\rm
FP}(\Phi_{i})\,.\,P(\Phi_{i};N)\,,$ (3.16)

where ${\cal L}_{\rm FP}(\Phi_{i})$ is the second-order Fokker-Planck
differential operator and $P(\Phi_{i};N)$ is the probability density function
of the stochastic process that is related to the probability of finding the
phase-space variables at a given value $\Phi_{i}=\\{\Phi,\Pi\\}$ at some time
$N$. However such a quantity is not of primary concern to us since we are not
interested in studying the phase-space dynamics of the inflaton666This would
have been our primary goal if we were studying the initial conditions for
inflation, or exit from eternal inflation to a SR classical regime [93, 94]..
Rather, we are interested in finding the probability distribution
$P_{\Phi_{i}}({\cal N})$ of the number of e-folds ${\cal N}$. Note the
important difference between our time variable $N$ and the stochastic variable
${\cal N}$. $N$ denotes the background expansion of the Universe, while ${\cal
N}$ is the number of e-folds of expansion obtained from the Langevin equations
with fixed boundary conditions in the IR field space, $\phi_{\rm en}$ and
$\phi_{\rm ex}$. The coarse-grained curvature perturbation $\zeta_{\rm cg}$ is
related to the stochastic number of e-folds ${\cal N}$ via the stochastic
$\delta{\cal N}$ formalism [47, 55, 58, 68]

$\zeta_{\rm cg}\equiv\zeta(\Phi_{i})={\cal N}-\langle{\cal
N}(\Phi_{i})\rangle\,,$ (3.17)

with

$\langle{\cal N}(\Phi_{i})\rangle=\int_{0}^{\infty}\,{\cal
N}\,P_{\Phi_{i}}({\cal N})\,{\rm d}{\cal N}\,,$ (3.18)

where the PDF $P_{\Phi_{i}}({\cal N})$ of the number of e-folds satisfies the
adjoint FPE which we discuss below in Sec. 3.2. Note that ${\cal N}(\Phi_{i})$
and $P_{\Phi_{i}}({\cal N})$ correspond to ${\cal N}(\Phi,\Pi)$ and
$P_{(\Phi,\Pi)}({\cal N})$ respectively.

###  3.2 Adjoint Fokker-Planck equation and first-passage time analysis

The adjoint FPE for the PDF $P_{\Phi_{i}}({\cal N})$ corresponding to the
general Langevin equation, Eq. (3.11), is given by

$\frac{\partial}{\partial{\cal N}}\,P_{\Phi_{i}}({\cal
N})=\left[D_{i}\,\frac{\partial}{\partial\Phi_{i}}+\frac{1}{2}\,\Sigma_{ij}\,\frac{\partial^{2}}{\partial\Phi_{i}\partial\Phi_{j}}\right]P_{\Phi_{i}}({\cal
N})\,.$ (3.19)

Our primary goal is to solve Eq. (3.19), with appropriate boundary conditions
for $\Phi_{i}\equiv\\{\Phi,\Pi\\}$ in order to compute the PDF
$P_{\Phi_{i}}({\cal N})\equiv P_{\Phi,\Pi}({\cal N})$. A physically well-
motivated set of boundary conditions includes an absorbing boundary at smaller
field values $\phi^{(\rm A)}$ closer to the end of inflation and a reflecting
boundary at a larger field value $\phi^{(\rm R)}$ closer to the CMB scale. The
PDF at the boundaries satisfies

  1. 1.

Absorbing boundary at $\phi^{(\rm A)}$

$P_{\Phi=\phi^{(\rm A)},\Pi}({\cal N})=\delta_{D}({\cal
N})\leavevmode\nobreak\ ,$ (3.20)

  2. 2.

Reflecting boundary at $\phi^{(\rm R)}$

$\frac{\partial}{\partial\Phi}P_{\Phi=\phi^{(\rm R)},\Pi}({\cal
N})=0\leavevmode\nobreak\ .$ (3.21)

The absorbing boundary condition ensures that for $\Phi<\phi^{(\rm A)}$, the
dynamics is heavily drift dominated and quantum diffusion effects are
negligible. Similarly, the reflecting boundary condition arises from assuming
that the potential is steep enough in the region $\Phi>\phi^{(\rm R)}$ so that
a freely diffusing inflaton can not climb back to a region of the potential
beyond $\phi^{(\rm R)}$. Both the boundary conditions play a crucial role in
determining the functional form of the PDF, thus affecting the PBH mass
fraction.

A convenient method for determining the PDF, as discussed in Ref. [55],
involves considering the ‘characteristic function’ (CF) $\chi_{\cal
N}(q;\Phi_{i})$, given by 777The subscript ${\cal N}$ in $\chi_{\cal
N}(q;\Phi_{i})$ denotes that the characteristic function is obtained by taking
the Fourier transformation of the PDF with respect to ${\cal N}$, and hence
$\chi_{\cal N}(q;\Phi_{i})$ is in fact independent of ${\cal N}$ and only a
function of $q,\,\Phi$.

$\chi_{\cal N}(q;\Phi_{i})\equiv\langle e^{i\,q\,{\cal
N}}\rangle=\int_{-\infty}^{\infty}e^{i\,q\,{\cal N}}P_{\Phi_{i}}({\cal
N})\,{\rm d}{\cal N}\,,$ (3.22)

which is the Fourier transform of the PDF $P_{\Phi_{i}}({\cal N})$ w.r.t the
dummy variable $q$ (which is a complex number in general). Hence the PDF is
the inverse Fourier transform of the CF:

$P_{\Phi_{i}}({\cal N})=\frac{1}{2\pi}\int_{-\infty}^{\infty}e^{-i\,q\,{\cal
N}}\chi_{\cal N}(q;\Phi_{i})\,{\rm d}q\,.$ (3.23)

Since the PDF satisfies the adjoint FPE, Eq. (3.19), the CF satisfies

$\left[D_{i}\,\frac{\partial}{\partial\Phi_{i}}+\frac{1}{2}\,\Sigma_{ij}\,\frac{\partial^{2}}{\partial\Phi_{i}\partial\Phi_{j}}+iq\right]\chi_{\cal
N}(q;\Phi_{i})=0\,,$ (3.24)

which is a partial differential equation with one less dynamical variable than
the adjoint FPE. The corresponding boundary conditions, Eqs. (3.20) and
(3.21), for the characteristic function are given by

$\chi_{\cal N}(q;\phi^{(\rm A)},\Pi)=1\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\frac{\partial}{\partial\Phi}\,\chi_{\cal N}(q;\phi^{(\rm
R)},\Pi)=0\leavevmode\nobreak\ .$ (3.25)

The characteristic equation, Eq. (3.24), corresponding to a potential
$V(\phi)$ in a general situation is quite difficult to solve. In practice, one
has to make crucial approximations regarding the classical drift $D_{i}$ and
the quantum noise $\xi_{i}$. The most common approximation used in the
literature assumes that the noise matrix elements $\Sigma_{ij}$ in Eq. (3.15)
are of the de Sitter-type, i.e. (see Sec. 4)

$\Sigma_{\phi\phi}=\left(\frac{H}{2\pi}\right)^{2}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\Sigma_{\phi\pi},\,\Sigma_{\pi\pi}\simeq 0\,.$ (3.26)

We now specialise to the case of quantum diffusion across a flat segment of
the inflaton potential, as discussed in Sec. 2 and shown in Fig. 2. It is
helpful to make a change of variables

$f=\frac{\Phi-\phi_{\rm ex}}{\Delta\phi_{\rm well}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ y=\frac{\Pi}{\pi_{\rm cr}}\,,$
(3.27)

where $f$ is the fraction of the flat well which remains to be traversed and
$y$ is the momentum relative to the critical momentum defined in Eq. (2.17),
the initial momentum for which the fields comes to a halt at $\phi_{\rm ex}$.
The CF, Eq. (3.24), then becomes (see Ref. [68])

$\left[\frac{1}{\mu^{2}}\frac{\partial^{2}}{\partial
f^{2}}-3y\left(\frac{\partial}{\partial f}+\frac{\partial}{\partial
y}\right)+iq\right]\chi_{\cal N}(q;f,y)=0\,,$ (3.28)

where

$\mu^{2}\simeq\frac{\Delta\phi^{2}_{\rm well}}{m_{p}^{2}}\,\frac{1}{v_{\rm
well}\,}\,,$ (3.29)

with $v_{\rm well}=V_{\rm well}/m_{p}^{4}$, where $V_{\rm well}$ is the height
of the flat quantum well. The corresponding boundary conditions now become

$\chi_{\cal N}(q;0,y)=1\,,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \frac{\partial}{\partial f}\,\chi_{\cal N}(q;1,y)=0\,,$
(3.30)

Such a system has been solved [68] in two distinct limits, namely

  * •

Free stochastic diffusion for which $\pi_{\rm en}\ll\pi_{\rm cr}\Rightarrow
y_{\rm en}\ll 1$, implying that the classical drift term, Eq. (3.12), can be
safely ignored, in which case the PDF takes the form (see Refs. [55, 58])

$P_{f}({\cal
N})=\sum_{n=0}^{\infty}A_{n}\,\sin{\left[\left(2n+1\right)\,\frac{\pi}{2}\,f\right]}\,{\rm
e}^{-\Lambda_{n}\,{\cal N}}\leavevmode\nobreak\ ,$ (3.31)

where

$A_{n}=\left(2n+1\right)\frac{\pi}{\mu^{2}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\Lambda_{n}=\left(2n+1\right)^{2}\frac{\pi^{2}}{4}\frac{1}{\mu^{2}}\,.$ (3.32)

Figure 3:  Left panel: the full PDF for the flat quantum well as a function of
${\cal N}$ for different values of the initial condition $\Phi=\phi_{i}$,
expressed in terms of $f=\Phi/\Delta\phi_{\rm well}$ (for simplicity, we
assume $\phi_{\rm ex}=0$ here). Right panel: the full PDF as a function of
initial field value $\Phi$ for realizations which have different values of
$\cal N$.

The full PDF as a function of ${\cal N}$ is plotted in the left panel of Fig.
3. In the limit ${\cal N}\gg 1$, the PDF exhibits an ‘exponential tail’ of the
form

$P_{\Phi}({\cal N})\simeq A_{0}\,e^{-\Lambda_{0}{\cal N}}\,.$

It follows from Eq. (3.32) that the amplitude, $A_{0}$ and coefficient of the
exponential, $\Lambda_{0}$, are given by

$A_{0}=\frac{\pi}{\mu^{2}}\,,\leavevmode\nobreak\ \leavevmode\nobreak\
\Lambda_{0}=\frac{1}{4}\,\frac{\pi^{2}}{\mu^{2}}\,,$

(note that $\Lambda_{0}$ is independent of $\Phi$), so that

$P_{\Phi}({\cal
N})\simeq\frac{\pi}{\mu^{2}}\,\cos{\left[\frac{\pi}{2}\left(\frac{\Phi}{\Delta\phi_{\rm
well}}-1\right)\right]}\,e^{-\frac{1}{4}\,\frac{\pi^{2}}{\mu^{2}}\,{\cal
N}}\,.$

In fact, the exponential tail was shown in Ref. [58] to be a universal feature
of the PDF for quantum diffusion across a generic slow-roll potential with
absorbing and reflecting boundary conditions, Eqs. (3.20) and (3.21). Larger
values of $f$ correspond to more quantum diffusion before exiting the flat
quantum well and hence result in more prominent exponential tails. Notice that
the PDFs saturate towards $f=1$ which is a consequence of the reflective
boundary condition given in Eq. (3.21).

The full PDF as a function of initial field value $\Phi$ is plotted in the
right panel of Fig. 3 for realizations which have different values of $\cal
N$. It is clear that for ${\cal N}\rightarrow 0$, starting from
$f=\Phi/\Delta\phi_{\rm well}\simeq 0$ yields a sharply peaked distribution,
in accordance with the absorbing boundary condition given in Eq. (3.20).

  * •

Large classical drift where $\pi_{\rm en}\gg\pi_{\rm cr}\Rightarrow y_{\rm
en}\gg 1$. In this case the dynamics of the inflaton is primarily governed by
its classical drift and hence the PDF is approximately Gaussian even for
${\cal N}\gg 1$ (see Ref. [68]).

However, in cases where the power spectrum is amplified sufficiently to form
an interesting abundance of PBHs, the inflaton typically enters the
intermediate flat USR segment from the CMB scale SR-I phase (see Fig. 2), with
speed of the order $\pi_{\rm en}\simeq\pi_{\rm cr}\Rightarrow y_{\rm en}\simeq
1$. In this case, both classical drift and stochastic diffusion become
important (at least initially during the entry into the USR segment) and hence
the aforementioned approximations will not be valid. Furthermore, the de
Sitter approximations for the noise matrix elements, Eq. (3.26), might
breakdown [95] during the transition into the USR phase. Consequently, it
becomes important to estimate the noise matrix elements more accurately.

We conclude that in order to properly use the stochastic $\delta{\cal N}$
formalism to estimate the abundance of PBHs, one must correctly determine the
PDF $P_{\Phi_{i}}({\cal N})$ from the adjoint FPE Eq. (3.19) with appropriate
boundary conditions. As discussed above, this can be carried out in two
important steps:

  1. 1.

Calculate the noise matrix elements $\Sigma_{ij}$ from Eq. (3.15) accurately
for the transitions between the CMB scale slow roll and subsequent slow-roll
violating phases.

  2. 2.

Determine the form of the PDF $P_{\Phi,\Pi}({\cal N})$, taking into account
the initial momentum with which the inflaton enters the USR segment.

In the rest of this paper, we carry out the first task of accurately computing
the noise matrix elements, first numerically in Sec. 4.2.1 for a potential
with a slow-roll violating feature, and then analytically in Sec. 4.2.2 for
the case of instantaneous transitions between different phases during
inflation. We reserve the second task to an upcoming paper [78].

##  4 Noise matrix elements in stochastic inflation

In this section we calculate the expressions for the noise matrix elements
$\Sigma_{ij}$, i.e. the correlators of the field and momentum noise operators
$\hat{\xi}_{i}=\\{\hat{\xi}_{\phi},\hat{\xi}_{\pi}\\}$. We do this initially
for standard slow-roll inflation, and compare the estimates for $\Sigma_{ij}$
computed using the pure de Sitter approximation to those obtained using the
slow-roll approximations.

The key equations that we use are the following: the definition of the noise
operators, Eq. (3.13), which along with a step-like $k$-space window function,
Eq. (3.10), leads to the noise correlators of Eq. (3.14), with the noise
correlation matrix $\Sigma_{ij}$ being given by Eq. (3.15). It is important to
note that these UV-noise mode functions are to be computed, not at Hubble
crossing, but at $k=\sigma aH$, where they chronologically become part of the
coarse-grained IR field and momentum, and provide quantum kicks. Hence, in
order to compute the elements of the noise matrix $\Sigma_{ij}$, we need to
compute the mode functions $\phi_{i_{k}}=\\{\phi_{k},\pi_{k}\\}$. This can be
done by solving the Mukhanov-Sasaki (MS) equation in terms of conformal time
$\tau$ defined in Eq. (2.9) 888Note that depending upon the situation, the MS
equation, Eq. (4.1), written in terms of the number of e-folds $N\sim{\rm
ln}(a)$ as  $\frac{{\rm d}^{2}v_{k}}{{\rm
d}N^{2}}+\left(1-\epsilon_{H}\right)\frac{{\rm d}v_{k}}{{\rm
d}N}+\left[\left(\frac{k}{aH}\right)^{2}-\frac{z^{\prime\prime}}{z}\right]v_{k}=0\,$
might be more useful. We note that in terms of $N$, the MS equation features a
friction term, and both the terms inside the square bracket evolve with time.
However, in terms of conformal time, $\tau$, it is a simple harmonic
oscillator equation with time dependent mass terms
$(aH)^{-2}z^{\prime\prime}/z$, while the comoving mode frequency $k$ is fixed,
which is why we choose to work with conformal time.

$v_{k}^{\prime\prime}+\left(k^{2}-\frac{z^{\prime\prime}}{z}\right)v_{k}=0\leavevmode\nobreak\
,$ (4.1)

where

$\displaystyle z$ $\displaystyle=$ $\displaystyle
am_{p}\sqrt{2\epsilon_{H}}\leavevmode\nobreak\ ,$ (4.2)
$\displaystyle\frac{z^{\prime\prime}}{z}$ $\displaystyle=$
$\displaystyle(aH)^{2}\left[2+2\epsilon_{H}-3\eta_{H}+2\epsilon_{H}^{2}+\eta_{H}^{2}-3\epsilon_{H}\eta_{H}-\frac{1}{aH}\,\eta^{\prime}_{H}\right]\leavevmode\nobreak\
,$ (4.3)

with appropriate initial conditions. The expressions for the mode functions
$\phi_{i_{k}}$ in the spatially flat gauge999In this work we compute the mode
functions $\\{\phi_{k},\pi_{k}\\}$, and hence the noise matrix elements,
$\Sigma_{ij}$, in the spatially flat gauge, while the Langevin equations are
written in the uniform-$N$ gauge. This introduces small corrections to the
noise terms which we assume to be negligible [60]. We discuss this further in
Sec. 5. are given by (see App. A)

$\phi_{k}=\frac{v_{k}}{a}\leavevmode\nobreak\ ,\leavevmode\nobreak\
\leavevmode\nobreak\ \pi_{k}=\frac{{\rm d}}{{\rm
d}N}\left(\frac{v_{k}}{a}\right)\,.$ (4.4)

From Sec. 3.2, it is clear that in order to accurately compute the noise
matrix elements $\Sigma_{ij}$ appearing in the adjoint FPE, Eq. (3.19), we
need to solve the MS equation as accurately as possible. For slow-roll
inflation, all relevant scales were sub-Hubble at early times, and hence we
impose the Bunch-Davies [96] initial conditions

$\lim_{k\tau\longrightarrow-\infty}v_{k}(\tau)=\frac{1}{\sqrt{2k}}\,e^{-ik\tau}\leavevmode\nobreak\
.$ (4.5)

We introduce a convenient new time variable, $T$, defined as

$T=-k\tau=\frac{k}{aH}\,.$ (4.6)

During quasi-dS expansion, the conformal time $\tau$ runs from $-\infty$ to
$0$, so $T$ runs from $\infty$ to $0$. Modes undergo Hubble-exit at $T\equiv
k/(aH)=1$, and the sub- and super-Hubble regimes correspond to $T\gg 1$ and
$T\ll 1$ respectively.

In terms of $T$ the MS equation, Eq. (4.1), takes the form

$\frac{{\rm d}^{2}v_{k}}{{\rm
d}T^{2}}+\left(\,1-\frac{\nu^{2}-\frac{1}{4}}{T^{2}}\,\right)v_{k}=0\,,$ (4.7)

where

$\nu^{2}=\frac{1}{(aH)^{2}}\,\frac{z^{\prime\prime}}{z}+\frac{1}{4}\,.$ (4.8)

For slow-roll inflation, $\nu^{2}$ is greater than or equal to $9/4$ at early
times and increases monotonically towards the end of inflation. In the limit
where $\nu$ is a constant, the MS Eq. (4.7) can be converted to a Bessel
equation as shown in App. B. In what follows, we start with the computation of
the noise-matrix elements for the case of featureless slow-roll potentials,
before proceeding to discuss the case of potentials possessing a slow-roll
violating feature.

###  4.1 Featureless potentials

In the case of a featureless potential for which slow roll is a good
approximation up until the end of inflation, the effective mass term
$(aH)^{-2}z^{\prime\prime}/z$ in the MS Eq. (4.1) is almost a constant and
evolves monotonically. Hence the MS Eq. (4.7) can be solved analytically by
approximating $\nu$ in Eq. (4.8) to be a constant.

Let us first demonstrate this calculation for the case of the pure de Sitter
limit which is usually employed in the computation of noise matrices in the
stochastic formalism. In the pure dS limit, both $\epsilon_{H},\,\eta_{H}=0$,
leading to $z^{\prime\prime}/z=2a^{2}H^{2}$ and $\nu^{2}=9/4$. Since
$a(\tau)=-1/(H\tau)$ in the pure dS approximation,

$a(T)=\frac{k}{HT}\leavevmode\nobreak\ ,$ (4.9)

and hence the MS Eq. (4.7) reduces to the familiar form

$\frac{{\rm d}^{2}v_{k}}{{\rm
d}T^{2}}+\left(\,1-\frac{2}{T^{2}}\,\right)v_{k}=0\,.$ (4.10)

The general solution of this equation is given by

$v_{k}(T)=\frac{1}{\sqrt{2k}}\,\left[\,\alpha_{k}\,\left(1+\frac{i}{T}\right)\,e^{i\,T}\,+\,\beta_{k}\,\left(1-\frac{i}{T}\right)\,e^{-i\,T}\,\right]\,,$
(4.11)

where the positive and negative frequency Bogolyubov coefficients satisfy the
canonical normalisation (Wronskian) condition

$|\alpha_{k}|^{2}-|\beta_{k}|^{2}=1\leavevmode\nobreak\ .$ (4.12)

Imposing the Bunch-Davies boundary conditions given in Eq. (4.5)

$\lim_{T\gg 1}v_{k}(T)=\frac{1}{\sqrt{2k}}\,e^{iT}\,,$ (4.13)

selects out the positive frequency solution only, i.e $\alpha_{k}=1$,
$\beta_{k}=0$, resulting in the final expression for the MS mode functions

$v_{k}(T)=\frac{1}{\sqrt{2k}}\,\left(1+\frac{i}{T}\right)\,e^{i\,T}\,,$ (4.14)

from which we find the field and momentum mode functions from Eq. (4.4) to be

$\displaystyle\phi_{k}(T)$ $\displaystyle=$
$\displaystyle\frac{H}{\sqrt{2k^{3}}}\,\left(T+i\right)\,e^{i\,T}\,,$ (4.15)
$\displaystyle\pi_{k}(T)$ $\displaystyle=$
$\displaystyle-\frac{H}{\sqrt{2k^{3}}}\,i\,T^{2}\,e^{i\,T}\,.$ (4.16)

Using the above expressions for the mode functions, we derive exact
expressions101010Note that we have dropped the imaginary part of the off-
diagonal terms of the noise correlator matrix since they do not correspond to
classical noise sources [85]. for the noise matrix elements, Eq. (3.15), in
the form (recall they are evaluated at $k=\sigma aH$, hence when $T=\sigma$)

$\displaystyle\Sigma_{\phi\phi}$ $\displaystyle=$
$\displaystyle\left(1+\sigma^{2}\right)\,\left(\frac{H}{2\pi}\right)^{2}\leavevmode\nobreak\
,$ (4.17) $\displaystyle{\rm Re}(\Sigma_{\phi\pi})$ $\displaystyle=$
$\displaystyle-\sigma^{2}\,\left(\frac{H}{2\pi}\right)^{2}\leavevmode\nobreak\
,$ (4.18) $\displaystyle\Sigma_{\pi\pi}$ $\displaystyle=$
$\displaystyle\sigma^{4}\,\left(\frac{H}{2\pi}\right)^{2}\leavevmode\nobreak\
.$ (4.19)

In the stochastic inflation formalism the field and momentum variables are
coarse-grained on ultra-Hubble scales, where $\sigma\ll 1$. For example,
taking $\sigma=0.01$, we get ${\rm
Re}(\Sigma_{\phi\pi})=10^{-4}\,\Sigma_{\phi\phi}$ and
$\Sigma_{\pi\pi}=10^{-8}\,\Sigma_{\phi\phi}$ under the pure de Sitter
approximation. This motivates the usual practice of dropping the momentum-
induced noise terms $\Sigma_{\phi\pi}$ and $\Sigma_{\pi\pi}$ from the adjoint
FPE, Eq. (3.19).

Turning now to slow-roll inflation, even though $\epsilon_{H},\,\eta_{H}\ll
1$, the slow-roll parameters do not exactly vanish unlike in pure dS space.
Nevertheless, as long as the quasi-de Sitter expansion is valid (which is
justified since $\epsilon_{H}\ll 1$), the expression for the scale factor in
terms of $T$ is still given by Eq. (4.9). Under the slow-roll approximations,
the MS equation takes the general form Eq. (4.7) with $\nu\neq 3/2$. In fact
for realistic SR potentials, $\nu$ is roughly equal to $3/2$ and evolves
slowly and monotonically. Assuming $\nu$ to be a constant, and imposing Bunch-
Davies initial conditions, the expression for $v_{k}$ takes the form (see App.
B)

$v_{k}(T)=e^{i\left(\nu+\frac{1}{2}\right)\frac{\pi}{2}}\,\sqrt{\frac{\pi}{2}}\,\frac{1}{\sqrt{2k}}\,\sqrt{T}\,H_{\nu}^{(1)}(T)\,,$
(4.20)

where $H_{\nu}^{(1)}(T)$ is the Hankel function of the first kind. For
$\nu\neq 3/2$, using the expression for the super-Hubble limit of the Hankel
function111111 Expressions for $\Sigma_{ij}$ which are valid for any value of
$\nu$ in the super-Hubble limit are provided in App. C. given in Eq. (B.8), we
obtain expressions for the field and momentum mode functions

$\displaystyle\phi_{k}(T)$ $\displaystyle=$ $\displaystyle
e^{i\left(\nu-\frac{1}{2}\right)\frac{\pi}{2}}\,2^{\nu-\frac{3}{2}}\,\frac{\Gamma(\nu)}{\Gamma(3/2)}\,\frac{H}{\sqrt{2\,k^{3}}}\,T^{-\nu+\frac{3}{2}}\,,$
(4.21) $\displaystyle\pi_{k}(T)$ $\displaystyle=$
$\displaystyle-e^{i\left(\nu-\frac{1}{2}\right)\frac{\pi}{2}}\,2^{\nu-\frac{3}{2}}\,\frac{\Gamma(\nu)}{\Gamma(3/2)}\,\frac{H}{\sqrt{2\,k^{3}}}\,\left(-\nu+\frac{3}{2}\right)\,T^{-\nu+\frac{3}{2}}\,,$
(4.22)

which leads to the following expressions for the noise matrix elements
$\Sigma_{ij}$ on super-Hubble scales

$\displaystyle\Sigma_{\phi\phi}$ $\displaystyle=$ $\displaystyle
2^{2\left(\nu-\frac{3}{2}\right)}\,\left[\frac{\Gamma(\nu)}{\Gamma(3/2)}\right]^{2}\,\left(\frac{H}{2\pi}\right)^{2}\,T^{2\left(-\nu+\frac{3}{2}\right)}\,,$
(4.23) $\displaystyle{\rm Re}\left(\Sigma_{\phi\pi}\right)$ $\displaystyle=$
$\displaystyle-2^{2\left(\nu-\frac{3}{2}\right)}\,\left[\frac{\Gamma(\nu)}{\Gamma(3/2)}\right]^{2}\,\left(\frac{H}{2\pi}\right)^{2}\,\left(-\nu+\frac{3}{2}\right)\,T^{2\left(-\nu+\frac{3}{2}\right)}\,,$
(4.24) $\displaystyle\Sigma_{\pi\pi}$ $\displaystyle=$ $\displaystyle
2^{2\left(\nu-\frac{3}{2}\right)}\,\left[\frac{\Gamma(\nu)}{\Gamma(3/2)}\right]^{2}\,\left(\frac{H}{2\pi}\right)^{2}\,\left(-\nu+\frac{3}{2}\right)^{2}\,T^{2\left(-\nu+\frac{3}{2}\right)}\,.$
(4.25)

Figure 4: The numerically determined noise matrix elements for the slow-roll
D-brane KKLT potential given in Eq. (4.26) for $M=0.5\,m_{p}$, in the absence
of any features with $\sigma=0.01$, as a function of number of e-folds from
the end of inflation $N_{e}$: from top to bottom $|\Sigma_{\phi\phi}|$, $|{\rm
Re}(\Sigma_{\phi\pi})|$ and $|\Sigma_{\pi\pi}|$, blue, green and purple lines
respectively. We find significant differences between the numerical
calculation and the analytical estimation under the assumption of a de Sitter
expansion, Eqs. (4.17)-(4.19), in which case the ratio of
$\Sigma_{\phi\phi}:|{\rm Re(\Sigma_{\phi\pi}})|:\Sigma_{\pi\pi}$ is
$1:10^{-4}:10^{-8}$.

Recalling the definition of $T$ in Eq. (4.6) and the fact super-Hubble scales
correspond to $k=\sigma aH$, hence $T=\sigma$, it follows that the above
expressions demonstrate that all three noise terms scale as
$\Sigma_{ij}\propto\sigma^{2\left(-\nu+3/2\right)}$ on super-Hubble scales.
This is in contrast to the pure dS limit where the three noise terms in Eqs.
(4.17)-(4.19) behave differently, namely, $\Sigma_{\phi\phi}={\rm const.}$,
$\Sigma_{\phi\pi}\propto\sigma^{2}$ and $\Sigma_{\pi\pi}\propto\sigma^{4}$.
Hence, during SR inflation for which $\nu\simeq 3/2$, even though the
momentum-induced noise terms $\Sigma_{\phi\pi}$ and $\Sigma_{\pi\pi}$ are
small compared to the field noise $\Sigma_{\phi\phi}$, they may not be
negligible, depending upon the value of $(\nu-3/2)$. As mentioned previously,
for most slow-roll potentials, $\nu$ evolves slowly and monotonically. The
numerically determined noise matrix elements, $\Sigma_{ij}$, are shown in Fig.
4 for an example asymptotically flat SR potential, which we choose to be the
D-Brane KKLT potential [97, 98, 99, 100] which has the form

$V(\phi)=V_{0}\,\frac{\phi^{2}}{M^{2}+\phi^{2}}\,,$ (4.26)

where $M$ is the mass scale in the KKLT model which we have chosen to be
$M=0.5\,m_{p}$. We have chosen $\sigma=0.01$ as is the standard practice (see
Ref. [85]). We notice that the momentum induced noise terms $\Sigma_{\phi\pi}$
and $\Sigma_{\pi\pi}$ are much higher than their corresponding values in the
pure de Sitter limit. In particular, we find the ratio121212The ratio of the
noise terms is not strongly dependent on the value of the KKLT mass scale,
$M$. Smaller values of $M$ result in smaller values of $\epsilon_{H}$, without
changing the value of $\eta_{H}$ significantly. Therefore, in the quasi-dS
limit $\epsilon_{H}\ll 1$, the value of $\nu^{2}$ from Eq. (4.8) and hence the
noise matrix elements, Eqs. (4.23)-(4.25), do not change significantly. Such
weak dependence of the ratio of noise terms on the parameters of the potential
is also the case for $\alpha$-attractors as well as for a number of other
asymptotically flat potentials (see Ref. [99, 101]). However, in general, if a
change in the parameters of the potential changes the value of $\nu$, it will
change the ratio of the noise terms. of $\Sigma_{\phi\phi}:|{\rm
Re(\Sigma_{\phi\pi}})|:\Sigma_{\pi\pi}$ to be about $1:10^{-2}:10^{-4}$ for
large $N_{e}$ as opposed to the de Sitter analytic estimate of
$1:10^{-4}:10^{-8}$. Additionally, the momentum induced noise terms scale
approximately in the same way as the field noise $\Sigma_{\phi\phi}$ in
accordance with Eqs. (4.23)-(4.25) at early times during inflation when
$\nu\simeq{\rm const.}$. Towards the end of inflation, since $\nu$ starts to
evolve faster, our analytical results based on $\nu\simeq{\rm const.}$ are no
longer applicable.

###  4.2 Potentials with a slow-roll violating feature

Potentials possessing a feature that generates large, PBH-forming,
perturbations, typically exhibit slow-roll violation, during which the quasi-
dS approximation is still valid ($\epsilon_{H}\ll 1$), while $\eta_{H}\geq 1$
(see Ref. [80]). In particular, $\eta_{H}\simeq 3$ during an ultra slow-roll
phase as discussed in Sec. 2. From Eq. (4.3), the expression for the effective
mass term $z^{\prime\prime}/z$ under the quasi-dS approximation becomes

$\frac{1}{(aH)^{2}}\,\frac{z^{\prime\prime}}{z}\simeq
2-3\,\eta_{H}+\eta_{H}^{2}+\tau\,\frac{{\rm d}\eta_{H}}{{\rm d}\tau}\,.$
(4.27)

In this case, the inflationary dynamics undergoes transitions between a number
of phases driven by the behaviour of $\eta_{H}$. In single field models in
which perturbations grow sufficiently to produce an interesting abundance of
PBHs, the inflaton typically undergoes two important transitions (see Ref.
[84]). The first transition T-I occurs from the CMB scale SR-I to a near-USR
phase, followed by a second transition T-II, from the near-USR phase to the
subsequent second slow-roll phase, SR-II, before the end of inflation. For
some class of features (see Refs. [84, 102]), the second transition T-II also
leads to an intermediate constant-roll (CR) phase [103] during which
$\eta_{H}$ is negative, almost constant, and of order unity.

As a specific example, we consider a modified KKLT potential with an
additional tiny Gaussian bump-like feature [104]:

$V_{\rm
b}(\phi)=V_{0}\,\frac{\phi^{2}}{M^{2}+\phi^{2}}\,\left[1+A\,\exp{\left({-\frac{1}{2}\,\frac{(\phi-\phi_{0})^{2}}{\tilde{\sigma}^{2}}}\right)}\right]\leavevmode\nobreak\
,$ (4.28)

where $A$, $\tilde{\sigma}$ and $\phi_{0}$ represent the height, width and
position of the bump respectively. The evolution of $\eta_{H}$ and
$z^{\prime\prime}/z$ for this potential is shown in Fig. 5. Following Ref.
[104], we fix $M=0.5\,m_{p}$, and take the bump parameters to be $A=1.87\times
10^{-3}$, $\tilde{\sigma}=1.993\times 10^{-2}\,m_{p}$ and
$\phi_{0}=2.005\,m_{p}$. These bump parameter values lead to a peak in the
scalar power-spectrum of ${\cal P}_{\zeta}\sim 10^{-2}$ at a $k$ value
corresponding to $\sim 10^{17}\,{\rm g}$ PBHs, i.e. at the lower end of the
asteroid mass window where PBHs can make up all of the dark matter (see e.g.
Refs. [13, 14]).

Figure 5: Evolution of the effective mass term $(1/aH)^{2}z^{\prime\prime}/z$
in the Mukhanov-Sasaki equation Eq. (4.1) (solid green curve) and $\eta_{H}$
(dashed blue curve) for the modified KKLT potential featuring a tiny Gaussian
bump as given in Eq. (4.28). The black-dashed line is
$(1/aH)^{2}z^{\prime\prime}/z$ for a de Sitter expansion (Eq. 4.8), namely
$\nu=3/2$. In the modified KKLT case, $(1/aH)^{2}z^{\prime\prime}/z$ makes a
sharp yet smooth dip around the transition from the CMB scale SR-I to the
subsequent near-USR phase, after which it remains almost constant throughout
the USR and constant-roll (CR) phases (but with $\nu>3/2$), until the inflaton
enters into another slow-roll phase, SR-II, before the end of inflation. The
corresponding noise matrix elements associated with this potential are shown
in Fig. 6.

The inflationary dynamics in this case display the aforementioned three key
phases, namely SR-I, USR and CR with $\eta_{H}$ making sharp (yet smooth)
transitions between them, as shown by the dashed blue curve in Fig. 5.
However, during the second transition from USR to the CR phase ($15\leq
N_{e}\leq 30$), the effective mass term $(aH)^{-2}z^{\prime\prime}/z$ remains
nearly constant131313 The effective mass term remains almost constant during
the second transition because of the upward step-like evolution of $\eta_{H}$
as a function of $N_{e}$. In the quasi-dS approximation, $\epsilon_{H}\ll 1$,
the effective mass term, Eq. (4.27), becomes
$(aH)^{-2}z^{\prime\prime}/z\simeq
2-3\eta_{H}+\eta_{H}^{2}+\tau\eta_{H}^{\prime}$. If $\eta_{H}$ has the form
$\eta_{H}=\frac{3}{2}+C\,\tanh{\left[C\left(N_{e}-\tilde{N_{e}}\right)\right]}\,$
where $C$ is a constant and $\tilde{N_{e}}$ is the value of $N_{e}$ at which
$\eta_{H}=3/2$, i.e. an upward step, then
$(aH)^{-2}z^{\prime\prime}/z=C^{2}-1/4={\rm const.}\,$. Note that the
effective mass term is only constant for an upwards step in $\eta_{H}$, and
not for a downward step, as occurs at the first transition., as emphasized in
Ref. [84]. The evolution of the mode functions (and hence the noise matrix
elements) is determined by $(aH)^{-2}z^{\prime\prime}/z$ through the MS Eq.
(4.1). The expression for the mode functions therefore remains the same in the
subsequent CR phase through the second transition because of the duality first
noticed by Wands (see Ref. [105]). Hence it is only necessary to follow the
evolution through the first transition, T-I, from SR-I to the near-USR phase.

In what follows, we will first describe how to compute the noise matrix
elements $\Sigma_{ij}$ numerically for the potential Eq. (4.28), before
finding accurate analytic solutions for them. Note that we use this particular
model to demonstrate our numerical framework because of its mathematical
simplicity and efficiency. However, the results we present are representative
of models with a broad range of features, including inflection point-like
behaviour. This is because, as shown in Ref. [84], the behaviour of the
effective mass term $z^{\prime\prime}/z$ is similar across these large class
of models, hence our primary conclusions will apply to all of them, and not
just this modified KKLT model.

####  4.2.1 Numerical analysis

In order to numerically compute the noise matrix elements for the potential
Eq. (4.28), our strategy is to split the mode functions $\phi_{k}$, $\pi_{k}$
and $v_{k}$ into their real and imaginary parts (see Ref. [102])

$v_{k}=v^{(R)}_{k}+i\,v^{(I)}_{k}\leavevmode\nobreak\ ,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\phi_{k}=\phi^{(R)}_{k}+i\,\phi^{(I)}_{k}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\pi_{k}=\pi^{(R)}_{k}+i\,\pi^{(I)}_{k}\leavevmode\nobreak\ .$ (4.29)

We impose Bunch-Davies initial conditions Eq. (4.5) deep in the sub-Hubble
regime, $k\gg aH$, to obtain $v^{(R)}_{k}$, $v^{(I)}_{k}$. Using Eq. (4.4), we
then find the real and imaginary parts of $\phi_{k}$ and $\pi_{k}$:

$\displaystyle\phi_{k}\equiv\phi^{(R)}_{k}+i\,\phi^{(I)}_{k}$ $\displaystyle=$
$\displaystyle\frac{1}{a}\,\left[\,v^{(R)}_{k}+i\,v^{(I)}_{k}\,\right]\leavevmode\nobreak\
,$ (4.30) $\displaystyle\pi_{k}\equiv\pi^{(R)}_{k}+i\,\pi^{(I)}_{k}$
$\displaystyle=$
$\displaystyle\frac{1}{a}\,\left[\,\left(\,\frac{v^{\prime(R)}_{k}}{aH}-v^{(R)}_{k}\,\right)+i\,\left(\,\frac{v^{\prime(I)}_{k}}{aH}-v^{(I)}_{k}\,\right)\,\right]\leavevmode\nobreak\
.$ (4.31)

Figure 6: The numerically determined noise matrix elements, $\Sigma_{ij}$,
with $\sigma=0.01$, for the modified KKLT potential with a tiny Gaussian bump,
Eq. (4.28), leading to a realistic smooth transition from SR-I to a near-USR
phase. (Note that the plot shows the behaviour of $\Sigma_{ij}$ only in the
vicinity of the USR regime.) The transition leads to an enhancement of the
momentum induced noise terms, $\Sigma_{\phi\pi}$ and $\Sigma_{\pi\pi}$,
relative to the field noise, $\Sigma_{\phi\phi}$, in the USR epoch.

Substituting Eqs. (4.29)-(4.31) into Eq. (3.15), we derive the following
compact expressions for the noise matrix elements $\Sigma_{ij}$

$\displaystyle\Sigma_{\phi\phi}$ $\displaystyle=$
$\displaystyle\left(1-\epsilon_{H}\right)\frac{k^{3}}{2\pi^{2}}\times\frac{1}{a^{2}}\times\left[\,\left(v^{(R)}_{k}\right)^{2}+\left(v^{(I)}_{k}\right)^{2}\,\right]\Bigg{|}_{k=\sigma
aH}\leavevmode\nobreak\ \leavevmode\nobreak\ $ (4.32) $\displaystyle{\rm
Re}\left(\Sigma_{\pi\phi}\right)$ $\displaystyle=$ $\displaystyle{\rm
Re}\left(\Sigma_{\phi\pi}\right)=\left(1-\epsilon_{H}\right)\frac{k^{3}}{2\pi^{2}}\times\frac{1}{a^{2}}\times\left[v^{(R)}_{k}\left(\frac{v^{\prime(R)}_{k}}{aH}-v^{(R)}_{k}\right)+v^{(I)}_{k}\left(\frac{v^{\prime(I)}_{k}}{aH}-v^{(I)}_{k}\right)\right]\Bigg{|}_{k=\sigma
aH}\leavevmode\nobreak\ \leavevmode\nobreak\ $ (4.33)
$\displaystyle\Sigma_{\pi\pi}$ $\displaystyle=$
$\displaystyle\left(1-\epsilon_{H}\right)\frac{k^{3}}{2\pi^{2}}\times\frac{1}{a^{2}}\times\left[\,\left(\frac{v^{\prime(R)}_{k}}{aH}-v^{(R)}_{k}\right)^{2}+\left(\frac{v^{\prime(I)}_{k}}{aH}-v^{(I)}_{k}\right)^{2}\,\right]\Bigg{|}_{k=\sigma
aH}\,.$ (4.34)

As we mentioned earlier, the imaginary part of the off-diagonal term
$\Sigma_{\pi\phi}$ does not correspond to a stochastic classical noise source
[85], hence we only need consider its real part in Eq. (4.33). The evolution
of the absolute values of $\Sigma_{\phi\phi}$, Re($\Sigma_{\phi\pi}$) and
$\Sigma_{\pi\pi}$ for the potential Eq. (4.28) are plotted in Fig. 6 for
$\sigma=0.01$, while Fig. 7 shows the ratios between the momentum-induced
noise terms and the field noise, $|{\rm
Re}(\Sigma_{\phi\pi})|/|\Sigma_{\phi\phi}|$ and
$|\Sigma_{\pi\pi}|/|\Sigma_{\phi\phi}|$ around the transition epoch. The
transition leads to an enhancement of the momentum induced noise terms
relative to the field noise with $\Sigma_{\pi\pi}>|{\rm
Re}(\Sigma_{\phi\pi})|>\Sigma_{\phi\phi}$. This is followed by a near-
exponential fall of each $\Sigma_{ij}$ during USR, since the slope of
$\Sigma_{ij}$ is almost constant during this epoch. We see that
$|\Sigma_{\pi\pi}|/|\Sigma_{\phi\phi}|\gtrsim 3\times|{\rm
Re}(\Sigma_{\phi\pi})|/|\Sigma_{\phi\phi}|$. At late times the noise matrix
elements begin to rise again and asymptote to constant values, and the
hierarchy between the noise terms gets reversed back to $\Sigma_{\pi\pi}<|{\rm
Re}(\Sigma_{\phi\pi})|<\Sigma_{\phi\phi}$. We also notice that the asymptotic
value of each $\Sigma_{ij}$ at late times is greater than its corresponding
value in the SR-I phase.

Figure 7: The ratios of the momentum-induced noise terms and the field noise,
$|{\rm Re}(\Sigma_{\phi\pi})|/|\Sigma_{\phi\phi}|$ in green and
$|\Sigma_{\pi\pi}|/|\Sigma_{\phi\phi}|$ in purple, with $\sigma=0.01$, for the
potential Eq. (4.28) with a tiny Gaussian bump as a function of $N_{\rm e}$
around the SR-I to USR transition. The transition from SR-I to USR leads to an
enhancement of the momentum induced noise terms, $\Sigma_{\phi\pi}$ and
$\Sigma_{\pi\pi}$, relative to the field noise, $\Sigma_{\phi\phi}$, in the
USR epoch.

From Figs. 6 and 7, we conclude that the noise matrix elements for a potential
with a PBH-forming feature evolve in a more complicated way than for pure de
Sitter or pure slow-roll. We next show that the aforementioned interesting
features of the noise terms across different epochs, such as during SR-I,
immediately after the transition from SR-I to USR, as well as the late time
asymptote, can be understood by making appropriate analytical approximations.
In the following subsection, we compute the noise matrix elements analytically
by assuming the transition T-I from SR-I to the near-USR phase to be
instantaneous. We will also demonstrate that in the quasi-dS limit,
$\epsilon_{H}\ll 1$, the noise terms are completely determined by the second
slow-roll parameter $\eta_{H}$.

####  4.2.2 Analytical treatment for instantaneous transitions

In order to compute $\Sigma_{ij}$ analytically, we consider an approach which
captures the key features of the full numerical evolution, namely solving the
MS Eq. (4.7) under the following assumptions.

  1. 1.

We assume the second slow-roll parameter $\eta_{H}$ to be a piece-wise
constant function which makes an instantaneous (yet finite) transition,
$\eta_{H}:\,\eta_{1}\to\eta_{2}$ at time $\tau=\tau_{1}$, given by

$\eta_{H}(\tau)=\eta_{1}+\left(\eta_{2}-\eta_{1}\right)\,\Theta(\tau-\tau_{1})\leavevmode\nobreak\
,$ (4.35)

where $\Theta$ is the Heaviside step function:

$\Theta(\tau-\tau_{1})=\begin{cases}\leavevmode\nobreak\
0\,,&\tau<\tau_{1}\leavevmode\nobreak\ ,\\\ \leavevmode\nobreak\
1\,,&\tau>\tau_{1}\leavevmode\nobreak\ .\end{cases}$ (4.36)

  2. 2.

The corresponding expression for $\nu$ given in Eq. (4.8) is then determined
using Eq. (4.3), which under the quasi-dS approximation, $\epsilon_{H}\simeq
0$, becomes Eq. (4.27). Using the expression for $\eta_{H}$ from Eq. (4.35) in
Eq. (4.27), we obtain

$\nu^{2}-\frac{1}{4}\equiv\frac{z^{\prime\prime}}{z}\,\tau^{2}={\cal
A}\,\tau\,\delta_{D}(\tau-\tau_{1})+\nu_{1}^{2}-\frac{1}{4}+\left(\nu_{2}^{2}-\nu_{1}^{2}\right)\,\Theta(\tau-\tau_{1})\,,$
(4.37)

where

${\cal A}=\eta_{2}-\eta_{1}\,,\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\
\nu_{1,2}^{2}-\frac{1}{4}=2-3\,\eta_{1,2}+\eta_{1,2}^{2}\,.$ (4.38)

Hence the piece-wise constant $\eta_{H}$ in Eq. (4.35) results in a piece-wise
constant $\nu$ in Eq. (4.37). We notice that the effective mass term
$z^{\prime\prime}/z$ contains a Dirac delta-function arising from the
derivative of the $\Theta$ function in Eq. (4.35). Note that for
$\eta_{2}>\eta_{1}$ (which is the case for the SR-I $\to$ USR transition in
Fig. 5) we have ${\cal A}>0$ and hence the term containing the Dirac delta-
function in Eq. (4.37) is negative (since $\tau<0$ during inflation). This
delta-function dip for an instantaneous transition analytically represents the
observed dip of finite width and depth for potentials with a smooth feature,
as seen in $(1/aH)^{2}z^{\prime\prime}/z$ in Fig. 5 (around $N_{e}\sim 32.5$).

  3. 3.

We impose Bunch-Davies initial conditions, Eq. (4.5), only for modes that
become super-Hubble at early times before the transition i.e. $\tau<\tau_{1}$.

  4. 4.

General solutions to the MS equation in different piece-wise constant $\nu$
regimes are matched during the transition $\tau=\tau_{1}$ by using the Israel
Junction conditions [106, 83, 95]

$\displaystyle v_{k}^{L}(\tau_{1})$ $\displaystyle=$ $\displaystyle
v_{k}^{E}(\tau_{1})\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ ({\rm
Continuity})\leavevmode\nobreak\ ,$ (4.39) $\displaystyle\frac{{\rm d}}{{\rm
d}\tau}v_{k}^{L}\bigg{|}_{\tau_{1}^{+}}-\frac{{\rm d}}{{\rm
d}\tau}v_{k}^{E}\bigg{|}_{\tau_{1}^{-}}$ $\displaystyle=$
$\displaystyle\int_{\tau_{1}^{-}}^{\tau_{1}^{+}}{\rm
d}\tau\frac{z^{\prime\prime}}{z}\,v_{k}^{L}(\tau)\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ ({\rm Differentiability})\,,$ (4.40)

where $v_{k}^{E}(\tau)$ and $v_{k}^{L}(\tau)$ are the mode functions before
and after the transition respectively, represented by

$\displaystyle v_{k}(\tau)\,$
$\displaystyle=\,\begin{cases}\leavevmode\nobreak\
v_{k}^{E}(\tau)\,,&\tau<\tau_{1}\leavevmode\nobreak\ ,\\\ \leavevmode\nobreak\
v_{k}^{L}(\tau)\,,&\tau>\tau_{1}\leavevmode\nobreak\ .\end{cases}$ (4.41)

We would ultimately like to derive expressions for the noise matrix elements
which can be expressed in terms of the mode functions $v_{k}$ in the following
compact form

$\displaystyle\Sigma_{\phi\phi}$ $\displaystyle=$
$\displaystyle\left(\frac{H}{2\pi}\right)^{2}\,T^{2}\,\Big{|}\,\sqrt{2k}\,v_{k}(T)\,\Big{|}^{2}\bigg{|}_{T=\sigma},$
(4.42) $\displaystyle{\rm Re}\left(\Sigma_{\pi\phi}\right)$ $\displaystyle=$
$\displaystyle-\left(\frac{H}{2\pi}\right)^{2}T^{2}{\rm
Re}\left(\sqrt{2k}v_{k}^{*}(T)\left[T\frac{{\rm d}}{{\rm
d}T}\left(\sqrt{2k}v_{k}(T)\right)+\sqrt{2k}v_{k}(T)\right]\right)\bigg{|}_{T=\sigma},$
(4.43) $\displaystyle\Sigma_{\pi\pi}$ $\displaystyle=$
$\displaystyle\left(\frac{H}{2\pi}\right)^{2}\,T^{2}\,\Big{|}\,T\,\frac{{\rm
d}}{{\rm
d}T}\left(\sqrt{2k}\,v_{k}(T)\right)+\sqrt{2k}\,v_{k}(T)\Big{|}^{2}\bigg{|}_{T=\sigma},$
(4.44)

where we take $\sigma=0.01$ as discussed earlier. We start with the
computation of noise matrix elements for an instantaneous transition in the
pure dS limit where $\nu_{1}=\nu_{2}=3/2$, before moving on to a general
transition between constant values of $\nu:\,\nu_{1}\to\nu_{2}$, with
$\nu_{2}>\nu_{1}$.

Case 1: Instantaneous transition in the pure dS limit with
$\nu_{1}=\nu_{2}=\frac{3}{2}$

In the case of an instantaneous transition at $\tau=\tau_{1}$ in the pure dS
limit (first considered in Ref. [107]), we have $\eta_{1}=0$ and $\eta_{2}=3$
and the system makes a transition from a SR to an exact USR phase.
Accordingly, the effective mass term in the MS equation takes the form

$\frac{z^{\prime\prime}}{z}\,\tau^{2}={\cal
A}\,\tau\,\delta_{D}(\tau-\tau_{1})\,+\,2\,,$ (4.45)

where the transition strength is ${\cal A}=3$. The expressions for the mode
functions, obtained by solving Eq. (4.1) are given (in terms of $T=-k\tau$) by

$\displaystyle v_{k}(T)\equiv\begin{cases}\leavevmode\nobreak\
v_{k}^{E}(T)=\frac{1}{\sqrt{2k}}\,\left(1+\frac{i}{T}\right)\,e^{i\,T}\,,&T>T_{1}\leavevmode\nobreak\
,\\\
v_{k}^{L}(T)=\frac{1}{\sqrt{2k}}\,\left[\,\alpha_{k}\,\left(1+\frac{i}{T}\right)\,e^{i\,T}\,+\,\beta_{k}\,\left(1-\frac{i}{T}\right)\,e^{-i\,T}\,\right]\,,&T<T_{1}\,,\end{cases}$
(4.46)

where $\alpha_{k}$ and $\beta_{k}$ are constants of integration (to be
determined from the Israel junction conditions given in Eqs. (4.39) and
(4.40)), while their derivatives are given by

$\displaystyle\frac{{\rm d}v_{k}}{{\rm
d}T}\equiv\begin{cases}\leavevmode\nobreak\ \frac{{\rm d}v_{k}^{E}}{{\rm
d}T}=\frac{1}{\sqrt{2k}}\,\left[-\frac{1}{T}+i\,\left(1-\frac{1}{T^{2}}\right)\right]\,e^{i\,T}\,,&T>T_{1}\leavevmode\nobreak\
,\\\ \frac{{\rm d}v_{k}^{L}}{{\rm
d}T}=\frac{1}{\sqrt{2k}}\left[\alpha_{k}\left[-\frac{1}{T}+i\left(1-\frac{1}{T^{2}}\right)\right]e^{i\,T}+\beta_{k}\left[-\frac{1}{T}-i\left(1-\frac{1}{T^{2}}\right)\right]e^{-i\,T}\right]\,,&T<T_{1}\,,\end{cases}$
(4.47)

where recall $T>T_{1}$ corresponds to the epoch before the transition and
$T<T_{1}$ to the epoch after the transition. Note that we have imposed Bunch-
Davies initial conditions on the mode function $v_{k}^{E}(T)$ before the
transition $T>T_{1}$, in accordance with our third assumption as discussed
above. The corresponding Fourier modes of the field fluctuations are obtained
from Eq. (4.4)

$\displaystyle\phi_{k}(T)\,=\,\begin{cases}\leavevmode\nobreak\
\phi_{k}^{E}(T)=\frac{H}{\sqrt{2k^{3}}}\,\left(T+i\right)\,e^{i\,T}\,,&T>T_{1}\leavevmode\nobreak\
,\\\
\phi_{k}^{L}(T)=\frac{H}{\sqrt{2k^{3}}}\,\left[\,\alpha_{k}\,\left(T+i\right)\,e^{i\,T}\,+\,\beta_{k}\,\left(T-i\right)\,e^{-i\,T}\,\right]\,,&T<T_{1}\leavevmode\nobreak\
.\end{cases}$ (4.48)

as are those of the field momentum fluctuations

$\displaystyle\pi_{k}(T)\,=\,\begin{cases}\leavevmode\nobreak\
\pi_{k}^{E}(T)=-\frac{H}{\sqrt{2k^{3}}}\,i\,T^{2}\,e^{i\,T}\,,&T>T_{1}\leavevmode\nobreak\
,\\\
\pi_{k}^{L}(T)=-\frac{H}{\sqrt{2k^{3}}}\,i\,T^{2}\,\left[\,\alpha_{k}\,e^{i\,T}-\beta_{k}\,e^{-i\,T}\,\right]\,,&T<T_{1}\leavevmode\nobreak\
.\end{cases}$ (4.49)

The Bogolyubov coefficients $\alpha_{k}$ and $\beta_{k}$, determined by
implementing the Israel junction conditions, Eqs. (4.39) and (4.40), are given
by

$\displaystyle a_{1}\,\alpha_{k}\,+\,b_{1}\,\beta_{k}\,$
$\displaystyle=\,d_{1}$ (4.50) $\displaystyle
a_{2}\,\alpha_{k}\,+\,b_{2}\,\beta_{k}\,$
$\displaystyle=\,d_{2}\leavevmode\nobreak\ ,$ (4.51)

which yields

$\displaystyle\alpha_{k}=\frac{d_{1}\,b_{2}-d_{2}\,b_{1}}{a_{1}\,b_{2}-a_{2}\,b_{1}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\beta_{k}=\frac{d_{2}\,a_{1}-d_{1}\,a_{2}}{a_{1}\,b_{2}-a_{2}\,b_{1}}\leavevmode\nobreak\
,$ (4.52)

where

$\displaystyle a_{1}$ $\displaystyle=T_{1}+i$ (4.53) $\displaystyle b_{1}$
$\displaystyle=\left(T_{1}-i\right)\,e^{-i\,2T_{1}}$ (4.54) $\displaystyle
a_{2}$ $\displaystyle=\left(1+{\cal A}\right)\,T_{1}+i\,\left(1+{\cal
A}-T_{1}^{2}\right)$ (4.55) $\displaystyle b_{2}$ $\displaystyle=\left(1+{\cal
A}\right)\,T_{1}-i\,\left(1+{\cal A}-T_{1}^{2}\right)\,e^{-i\,2T_{1}}$ (4.56)
$\displaystyle d_{1}$ $\displaystyle=T_{1}+i$ (4.57) $\displaystyle d_{2}$
$\displaystyle=T_{1}+i\,\left(1-T_{1}^{2}\right)$ (4.58)

With a little bit of algebra, we obtain141414Eqs. (4.59) and (4.60) agree with
the results obtained in Refs. [107, 95] for an inflaton potential consisting
of two linear regimes $V(\phi)\propto\phi$ and slopes $\alpha$, $\beta$ that
are joined at a point $\phi=\phi_{1}$ in which case $T_{1}$ is related to
$\phi_{1}$ and ${\cal A}=3(\alpha-\beta)/\alpha$.

$\displaystyle\alpha_{k}$ $\displaystyle=$ $\displaystyle 1-i\,\frac{{\cal
A}}{2}\,\frac{1}{T_{1}}\,\left(1+\frac{1}{T_{1}^{2}}\right)\,,$ (4.59)
$\displaystyle\beta_{k}$ $\displaystyle=$ $\displaystyle i\,\frac{{\cal
A}}{2}\,\frac{1}{T_{1}}\,\left(1+\frac{i}{T_{1}}\right)^{2}\,e^{i\,2T_{1}}\,,$
(4.60)

From the above expressions, we can compute the noise matrix elements
$\Sigma_{\phi\phi},\,\Sigma_{\phi\pi},\,\Sigma_{\pi\pi}$ in Eqs. (4.42) -
(4.44). The instantaneous transition leads to a mixing between the positive
and negative frequency solutions for the mode functions. Our analytical
results are shown in the left panel of Fig. 8 for ${\cal A}=+3$ 151515 The
noise matrix elements for the case of an instantaneous transition were also
estimated in Ref. [95] under the pure dS approximation for the Starobinsky
model [107]. Using our calculations, with a transition strength in the range
${\cal A}\leq 3$, one can obtain the noise matrix elements for the entire
parameter space of the Starobinsky model.. The asymptotic behaviour of
$\Sigma_{ij}$ at different epochs can be inferred from the scale dependence of
the Bogolyubov coefficients, given by Eqs. (4.59) and (4.60). Note that
different comoving modes contribute to the noise terms at different times.
Immediately after the transition, when $T\lesssim T_{1}\ll 1$, the noise
matrix elements are due to modes joining the coarse-graining scale $T=\sigma$
at this epoch, for which, using Eqs. (4.59) and (4.60),
$\alpha_{k},\,\beta_{k}\simeq-3i/(2T_{1}^{-3})$. Consequently, the noise terms
from Eqs. (4.42)-(4.44) behave as $\Sigma_{ij}\propto T_{1}^{-6}$. Since
$T_{1}=k/(a_{1}H_{1})$, $\Sigma_{ij}\propto k^{-6}\propto e^{6N_{e}}$, and the
ratio $\Sigma_{\phi\phi}:|{\rm Re(\Sigma_{\phi\pi}})|:\Sigma_{\pi\pi}$ is
given by $1:3:9$, to leading order in $T=\sigma$. Following this epoch, the
noise terms begin to rise exponentially, and the hierarchy between the field
and momentum induced terms gets reversed back to $\Sigma_{\phi\phi}>|{\rm
Re(\Sigma_{\phi\pi}})|>\Sigma_{\pi\pi}$. At sufficiently late times, when
$T\ll 1\ll T_{1}$, the noise matrix elements are due to modes for which
$\alpha_{k}\to 1$, while $\beta_{k}\simeq 3i/(2T_{1})\,e^{i2T_{1}}$ decays to
zero with oscillations. Hence the noise matrix elements asymptote to their
corresponding pre-transition (constant) values given by Eqs. (4.17)-(4.19).

Comparing the analytical results for an instantaneous transition in the pure
dS limit, shown in the left panel of Fig. 8, with the numerical results for a
potential with a PBH forming feature shown in Fig. 6, we conclude that the
former fails to capture161616Note, however, that the pure dS-transition is a
good approximation for potentials featuring a ‘flat’ segment (as shown in Fig.
2) rather than a bump. In this case Eqs. (4.46)-(4.60) can be used to compute
the noise terms (see Ref. [95]). the late time asymptotic properties of
$\Sigma_{ij}$. Therefore, in the following, we will compute $\Sigma_{ij}$
relaxing the pure dS approximation.

Case 2: Instantaneous transition between two constant values of $\nu$:
$\nu_{1}\to\nu_{2}$

In the case of an instantaneous transition at $\tau=\tau_{1}$ where $\nu$
makes a jump between the constant values $\nu_{1}\to\nu_{2}$, once again by
solving Eq. (4.1), the expressions for the mode functions are given (in terms
of $T=-k\tau$) by

$\displaystyle v_{k}(T)\,=\,\begin{cases}\leavevmode\nobreak\
\frac{1}{\sqrt{2k}}\,\sqrt{\frac{\pi}{2}}\,e^{i\left(\nu_{1}+\frac{1}{2}\right)\frac{\pi}{2}}\,\sqrt{T}\,H_{\nu_{1}}^{(1)}(T)\,,&T>T_{1}\leavevmode\nobreak\
,\\\ \leavevmode\nobreak\
\sqrt{T}\left[\,C_{1}^{L}\,H_{\nu_{2}}^{(1)}(T)\,+\,C_{2}^{L}\,H_{\nu_{2}}^{(2)}(T)\,\right]\,,&T<T_{1}\leavevmode\nobreak\
,\end{cases}$ (4.61)

and their derivatives are given by

$\displaystyle\frac{{\rm d}v_{k}}{{\rm d}T}=\begin{cases}\leavevmode\nobreak\
\frac{1}{\sqrt{2k}}\sqrt{\frac{\pi}{2}}e^{i\left(\nu_{1}+\frac{1}{2}\right)\frac{\pi}{2}}\frac{1}{\sqrt{T}}\left[\left(\frac{1}{2}-\nu_{1}\right)H_{\nu_{1}}^{(1)}(T)+TH_{\nu_{1}-1}^{(1)}(T)\right]\,,&T>T_{1}\\\
\leavevmode\nobreak\
\frac{1}{\sqrt{T}}\left[C_{1}^{L}\left(\left(\frac{1}{2}-\nu_{2}\right)H_{\nu_{2}}^{(1)}(T)+TH_{\nu_{2}-1}^{(1)}(T)\right)+C_{2}^{L}\left(\left(\frac{1}{2}-\nu_{2}\right)H_{\nu_{2}}^{(2)}(T)+TH_{\nu_{2}-1}^{(2)}(T)\right)\right]\,,&T<T_{1}\end{cases}$
(4.62)

By implementing the Israel junction conditions, Eqs. (4.39) and (4.40), the
constant coefficients of integration $C_{1}^{L}$ and $C_{2}^{L}$ can be shown
to satisfy the algebraic equations

$\displaystyle a_{1}\,C_{1}^{L}\,+\,b_{1}\,C_{2}^{L}\,$
$\displaystyle=\,d_{1}$ (4.63) $\displaystyle
a_{2}\,C_{1}^{L}\,+\,b_{2}\,C_{2}^{L}\,$
$\displaystyle=\,d_{2}\leavevmode\nobreak\ ,$ (4.64)

which yields

$\displaystyle
C_{1}^{L}=\frac{d_{1}\,b_{2}-d_{2}\,b_{1}}{a_{1}\,b_{2}-a_{2}\,b_{1}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
C_{2}^{L}=\frac{d_{2}\,a_{1}-d_{1}\,a_{2}}{a_{1}\,b_{2}-a_{2}\,b_{1}}\leavevmode\nobreak\
,$ (4.65)

where

$\displaystyle a_{1}$ $\displaystyle=H_{\nu_{2}}^{(1)}(T_{1})$ (4.66)
$\displaystyle b_{1}$ $\displaystyle=H_{\nu_{2}}^{(2)}(T_{1})$ (4.67)
$\displaystyle a_{2}$ $\displaystyle=\left(\,\frac{1}{2}-\nu_{2}-{\cal
A}\,\right)\,H_{\nu_{2}}^{(1)}(T_{1})+T_{1}\,H_{\nu_{2}-1}^{(1)}(T_{1})$
(4.68) $\displaystyle b_{2}$ $\displaystyle=\left(\,\frac{1}{2}-\nu_{2}-{\cal
A}\,\right)\,H_{\nu_{2}}^{(2)}(T_{1})+T_{1}\,H_{\nu_{2}-1}^{(2)}(T_{1})$
(4.69) $\displaystyle d_{1}$
$\displaystyle=\frac{1}{\sqrt{2k}}\,\sqrt{\frac{\pi}{2}}\,e^{i\left(\nu_{1}+\frac{1}{2}\right)\frac{\pi}{2}}\,H_{\nu_{1}}^{(1)}(T_{1})$
(4.70) $\displaystyle d_{2}$
$\displaystyle=\frac{1}{\sqrt{2k}}\,\sqrt{\frac{\pi}{2}}\,e^{i\left(\nu_{1}+\frac{1}{2}\right)\frac{\pi}{2}}\,\left[\,\left(\frac{1}{2}-\nu_{1}\right)\,H_{\nu_{1}}^{(1)}(T_{1})+T_{1}\,H_{\nu_{1}-1}^{(1)}(T_{1})\,\right]$
(4.71)

Figure 8: The noise matrix elements $\Sigma_{ij}$ computed analytically using
Eqs. (4.42)-(4.44), for $\sigma=0.01$. Left panel: an instantaneous transition
from SR to USR using the pure de Sitter limit, $\nu_{1}=\nu_{2}=3/2$. Right
panel: an instantaneous transition from a SR phase with $\nu_{1}=1.52$ to a
near-USR phase with $\nu_{2}=1.8$. In both cases the thin dashed lines show
the analytical asymptotes immediately after the transition, $\Sigma_{ij}\sim
e^{2{\cal A}N_{e}}$.

The resulting noise matrix elements, computed using Eqs. (4.42)-(4.44), are
shown in the right panel of Fig. 8. In order to compare our results with the
numerical calculation in Fig. 6, we choose $\nu_{1}=1.52$ and $\nu_{2}=1.8$.
These values correspond to $\eta_{1}=-0.02$ and $\eta_{2}=3.3$ respectively,
to match the values of $\eta_{H}$ during the SR-I and the near-USR epochs for
the modified KKLT potential with a Gaussian bump used for the numerical
calculation in Fig. 6.

As in the pure dS case, immediately after the transition, when $T\lesssim
T_{1}\ll 1$, the noise matrix elements fall nearly-exponentially with
$\Sigma_{ij}\sim e^{2{\cal A}N_{e}}$. The ratio $\Sigma_{\phi\phi}:|{\rm
Re(\Sigma_{\phi\pi}})|:\Sigma_{\pi\pi}$ is approximately $1:{\cal A}:{\cal
A}^{2}$ (where ${\cal A}\equiv\eta_{2}-\eta_{1}=3.32$ from Eq. (4.38)), and
nearly constant. However, following this epoch the noise terms begin to rise
and the hierarchy between the field and momentum induced terms is reversed
back to $\Sigma_{\phi\phi}>|{\rm Re}(\Sigma_{\phi\pi})|>\Sigma_{\pi\pi}$. At
sufficiently late times, $T\ll 1\ll T_{1}$, the coefficient of the negative
frequency solution $C_{2}^{L}$ becomes negligible, and the behaviour of
$\Sigma_{ij}$ can be understood from the constant $\nu$ expressions for the
noise terms, Eqs. (4.23)-(4.25). The late time ratio of noise terms is given
by $\Sigma_{\phi\phi}:|{\rm Re(\Sigma_{\phi\pi}})|:\Sigma_{\pi\pi}\to
1:\left(\nu_{2}-\frac{3}{2}\right):\left(\nu_{2}-\frac{3}{2}\right)^{2}$ i.e.
the values of $\Sigma_{ij}$ are higher than their pre-transition counterparts
in the SR-I phase. This matches the behaviour of the numerically calculated
noise matrix elements for the modified KKLT potential with a Gaussian bump
shown in Fig. 6.

The key results of our analytical calculations for an instantaneous transition
are:

  1. 1.

The expressions for the noise matrix elements in the pre-transition epoch are
given by Eqs. (4.23)-(4.25) with $\nu=\nu_{1}$, resulting in the ratios

$\Sigma_{\phi\phi}:|{\rm Re(\Sigma_{\phi\pi}})|:\Sigma_{\pi\pi}\to
1:\left(\nu_{1}-\frac{3}{2}\right):\left(\nu_{1}-\frac{3}{2}\right)^{2}\,.$
(4.72)

  2. 2.

Immediately after the transition, $\Sigma_{ij}\propto e^{2{\cal A}N_{e}}$, and

$\Sigma_{\phi\phi}:|{\rm Re(\Sigma_{\phi\pi}})|:\Sigma_{\pi\pi}\to 1:{\cal
A}:{\cal A}^{2}\,.$ (4.73)

  3. 3.

At sufficiently late times, the noise terms are again given by
Eqs.(4.23)-(4.25), but with $\nu=\nu_{2}$, yielding the ratios

$\Sigma_{\phi\phi}:|{\rm Re(\Sigma_{\phi\pi}})|:\Sigma_{\pi\pi}\to
1:\left(\nu_{2}-\frac{3}{2}\right):\left(\nu_{2}-\frac{3}{2}\right)^{2}\,.$
(4.74)

Comparing Fig. 8 with Fig. 6, we see that the analytical treatment assuming an
instantaneous transition between two constant values of $\nu$, $\nu_{1}$ and
$\nu_{2}$ (Case 2) captures most of the asymptotic properties of $\Sigma_{ij}$
for a potential with a PBH forming feature. This is in contrast to the pure dS
transition (Case 1) which was not able to capture the late-time asymptote
accurately, due to the assumption that $\nu_{1}=\nu_{2}=3/2$. Furthermore, the
pure dS transition also underestimates the momentum induced noise terms
$\Sigma_{\phi\pi}$ and $\Sigma_{\pi\pi}$ in the SR-I phase, as discussed in
Sec. 4.1.

We conclude this Section by briefly commenting on the degree of correlation
between the field and momentum noise terms, $\xi_{\phi}$ and $\xi_{\pi}$,
which can be quantified in terms of the ratio171717We thank David Wands for
pointing out the relevance of calculating $\gamma$.

$\gamma=\frac{|{\rm
Re}(\Sigma_{\phi\pi})|}{\sqrt{\Sigma_{\phi\phi}\Sigma_{\pi\pi}}}\,,$ (4.75)

where $\gamma$ is related to the determinant of the noise matrix by
$\gamma^{2}=1-{\rm det(\Sigma_{ij})}/(\Sigma_{\phi\phi}\Sigma_{\pi\pi})$. The
noise terms are maximally correlated if $\gamma=1$, while $\gamma=0$ implies
that $\xi_{\phi}$ and $\xi_{\pi}$ are independent (see Ref. [85]). For
featureless potentials, we find that $\gamma\simeq 1$ under the pure dS
approximation, using Eqs. (4.17)-(4.19), as well as under the SR
approximations, using Eqs. (4.23)-(4.25). For potentials with a PBH forming
feature, we also find that $\gamma\simeq 1$ throughout the three asymptotic
regimes181818We find $\gamma\ll 1$ only for brief transient periods when the
noise terms begin to rise after their exponential fall post transition.
described by Eqs. (4.72)-(4.74). This property of maximal correlation between
the noise terms is a direct consequence of the fact that for $\sigma\ll 1$,
the super-Hubble UV mode functions, $\phi_{k}$, are frozen by the time they
join the coarse-graining scale $k=\sigma aH$ (see Refs. [70, 110] for a
detailed discussion on the freezing behaviour of the UV modes). Consequently,
we conclude that quantum diffusion can be assumed to be sourced by a single
random noise term.

##  5 Outstanding steps to calculate the PBH mass fraction

In this Section we discuss what remains to be done to accurately calculate the
PBH mass fraction. In Sec. 5.1 we provide expressions for the primordial
probability distribution of ${\cal N}$, $P_{\Phi,\Pi}({\cal N})$, and the PBH
mass fraction, $\beta$, in terms of quantities which can be calculated using
the Fokker-Plank Eq. (3.19) and the noise terms that we have calculated. In
Sec. 5.2 we discuss the various complexities in calculating the PBH mass
fraction, in the stochastic inflation framework and more generally.

###  5.1 Analytic expressions

In Sec. 4 we accurately calculated the stochastic noise matrix elements
$\Sigma_{ij}$, for a sharp transition from SR to USR, using both analytical
and numerical techniques. Our ultimate aim is to determine the PDF of the
number of e-folds, $P_{\Phi,\Pi}({\cal N})$, by solving the adjoint Fokker-
Planck Eq. (3.19) (using appropriate boundary conditions) and then calculate
the mass fraction of PBHs $\beta_{\rm PBH}$. While this will be the primary
focus of our upcoming paper [78], here we outline the key expressions required
to compute the mass fraction.

The primordial probability distribution of ${\cal N}$, can be written in the
form

$P_{\Phi,\Pi}({\cal
N})=\sum_{n=0}^{\infty}\,B_{n}(\Phi,\Pi)\,e^{-\Lambda_{n}\,{\cal N}}\,,$ (5.1)

where the exact form of $B_{n}(\Phi,\Pi)$ has to be determined by imposing
appropriate boundary conditions and using the correct expressions for the
noise matrix elements $\Sigma_{ij}$ as determined in Sec. 4. Using the above
mathematical form for the primordial PDF, general expressions for the moments
of the number of e-folds, $\langle{\cal N}^{m}(\Phi,\Pi)\rangle$, (see Refs.
[47, 63]) and the mass fraction of PBHs, $\beta(\Phi,\Pi)$, can be derived.

From the stochastic $\delta{\cal N}$ formalism, the coarse-grained comoving
curvature perturbation is given by

$\zeta_{\rm cg}\equiv\zeta_{\rm cg}(\Phi,\Pi)={\cal N}-\langle{\cal
N}(\Phi,\Pi)\rangle\,,$ (5.2)

where the average number of e-folds, $\langle{\cal N}(\Phi,\Pi)\rangle$, can
be determined via

$\langle{\cal N}(\Phi,\Pi)\rangle=\int_{0}^{\infty}\,{\rm d}{\cal N}\,{\cal
N}\,P_{\Phi,\Pi}({\cal
N})=\sum_{n=0}^{\infty}\,B_{n}(\Phi,\Pi)\,\int_{0}^{\infty}\,{\rm d}{\cal
N}\,{\cal N}\,e^{-\Lambda_{n}\,{\cal N}}\leavevmode\nobreak\ ,$

which leads to

$\langle{\cal
N}(\Phi,\Pi)\rangle=\sum_{n=0}^{\infty}\,\frac{B_{n}(\Phi,\Pi)}{\Lambda_{n}^{2}}\,.$
(5.3)

Using the Press-Schechter formalism [108], the PBH mass fraction can be
estimated by integrating the probability distribution of the coarse grained
curvature perturbation, $P(\zeta_{cg})$, above the threshold for PBH
formation, $\zeta_{c}$. The threshold for PBH formation in terms of the number
of e-folds, ${\cal N}_{c}$, is given by

${\cal N}_{c}=\zeta_{c}+\langle{\cal
N}(\Phi,\Pi)\rangle=\zeta_{c}+\sum_{n=0}^{\infty}\,\frac{B_{n}(\Phi,\Pi)}{\Lambda_{n}^{2}}\,.$
(5.4)

Hence the PBH mass fraction is given by

$\beta(\Phi,\Pi)\equiv\int_{\zeta_{c}}^{\infty}\,P(\zeta_{\rm cg})\,{\rm
d}\zeta_{\rm cg}=\int_{\zeta_{c}+\langle{\cal
N}(\Phi,\Pi)\rangle}^{\infty}\,P_{\Phi,\Pi}({\cal N})\,{\rm d}{\cal
N}\leavevmode\nobreak\ ,$ (5.5)

which leads to the final expression for the (non-Gaussian) PDF resulting from
the stochastic inflation formalism

$\beta(\Phi,\Pi)=\sum_{n=0}^{\infty}\,\frac{B_{n}(\Phi,\Pi)}{\Lambda_{n}}\,\exp{\left[-\Lambda_{n}\,\left[\zeta_{c}\,+\,\sum_{m=0}^{\infty}\,\frac{B_{m}(\Phi,\Pi)}{\Lambda_{m}^{2}}\right]\right]}\leavevmode\nobreak\
.$ (5.6)

This can be compared with the Gaussian PDF obtained for typical fluctuations
in the perturbative approach, $\beta^{\rm G}(\Phi,\Pi)$

$\beta^{\rm G}(\Phi,\Pi)\simeq\frac{\sigma_{\rm
cg}}{\sqrt{2\pi}\zeta_{c}}\exp{\left[-\frac{\zeta_{c}^{2}}{2\sigma_{\rm
cg}^{2}}\right]}\,,$ (5.7)

with

$\sigma_{\rm cg}^{2}(\Phi,\Pi)=\int_{k(\Phi,\Pi)}^{k_{e}}\frac{{\rm
d}k}{k}\,{\cal P}_{\zeta}(k)\leavevmode\nobreak\ .$

To calculate the PBH mass fraction using Eq. (5.6) it is necessary to find the
coefficients, $B_{n}(\Phi,\Pi)$, and the exponents, $\Lambda_{n}$, by solving
the adjoint Fokker-Planck Eq. (3.19) using the field and momentum induced
noise terms in Sec. 4. While this task is reserved for our upcoming paper, we
expect that the sharp decline of the noise terms after the transition will
decrease the amount of quantum diffusion of the IR fields across the PBH-
forming feature. Therefore we expect the tail of the PDF to decline more
rapidly than what is usually found using the pure dS approximation without any
transitions. Indeed such behaviour of the PDF was found in Ref. [95] which
focused on a sharp transition in pure dS space using the linear potential
model of Starobinsky [107]. Numerical simulations carried out in Refs. [65,
69, 70] show that the canonical computation based on the pure dS noise terms
without any transition typically leads to inaccurate estimates of the PBH
abundance (over-estimates for some potentials and under-estimates for other
potentials). However, it is important to study the relative contributions of
the noise terms, $\Sigma_{ij}$, the potential, $V(\phi)$, and the boundary
conditions to the PDF separately. An analytical approach is well-suited to
this, and this is one of the primary goals of our upcoming paper.

###  5.2 Complexities

In this subsection we overview the outstanding complexities in accurately
calculating the PBH mass fraction.

  * •

Curvature perturbation vs density contrast: In the previous sub-section, we
outlined how the PBH mass fraction can be calculated using the PDF of the
curvature perturbation obtained from the stochastic inflation framework.
However the criterion for PBH formation is most accurately formulated in terms
of the non-linear density contrast $\delta_{l}$, see e.g. Ref. [109, 71, 40,
39, 38]. An accurate computation of the PDF of the density contrast needs the
knowledge of all the higher order ($n$-point) connected correlators of the
curvature perturbation, $\zeta$. Therefore a high-precision calculation of the
PBH mass fraction requires the joint probabilities rather than the one-point
PDF, $P[\zeta]$ (see Ref. [71, 39] for discussion of this issue).

  * •

Gauge corrections: In this work we compute the mode functions
$\\{\phi_{k},\pi_{k}\\}$, and hence the noise correlators of
$\\{\hat{\xi}_{\phi},\hat{\xi}_{\pi}\\}$, in the spatially flat gauge, however
the Langevin equations are written in the uniform-$N$ gauge. This induces
corrections to the noise terms that could be non-negligible when the slow-roll
approximations are violated [67, 66, 60]. However Refs. [60, 70] showed that
the gauge corrections are negligible for $\sigma\ll 1$.

  * •

Effects of backreaction: We have calculated the noise matrix elements by
treating the mode functions $\\{\phi_{k},\pi_{k}\\}$ as linear perturbations
in a deterministic (non-stochastic) inflationary background, as is the usual
practice in perturbation theory. In stochastic inflation, the noise terms
should in principle be evaluated in the stochastically evolving background of
the coarse-grained IR fields $\\{\Phi,\Pi\\}$. However, numerical simulations
demonstrate that such non-Markovian corrections due to the backreaction
effects of the stochastic IR background are negligible for single field
inflationary potentials with a large class of PBH forming features, such as a
flat segment, an inflection point, or a bump/dip and hence can be safely
ignored (see Refs. [70, 110]).

  * •

Classical $\delta N$ formalism: The PDF of the comoving curvature perturbation
can also be computed using the classical (non-linear) $\delta N$ formalism
developed in Refs. [74, 75, 76]. For potentials with a broad class of
features, the non-perturbative non-Gaussianity induced by stochastic effects
is usually expected to be dominant (see Ref. [70]). The relative significance
of the stochastic effects can be inferred from a classicality criterion,
expressed in terms of the parameter $j_{\rm
cl}=|V(\phi)\eta_{H}/(24\pi^{2}m_{p}^{4}\,\epsilon_{H})|$, obtained from a
saddle-point approximation of the stochastic integrals [47]. For a potential
with $j_{\rm cl}\ll 1$, stochastic effects can safely be ignored (except in
the far tail of the PDF). Hence, it is possible to construct potentials for
which the stochastic effects are specifically negligible by design, while the
classical non-linearities can be significant (see Ref. [50]). In such cases
the classical $\delta N$ formalism can be successfully used to compute the PDF
(see also Refs. [52, 54]).

  * •

Loop corrections: As outlined in Sec. 2, to generate a non-negligible
abundance of PBHs, the power spectrum of the primordial scalar perturbations
on small scales has to be roughly seven orders of magnitude larger than its
measured value on CMB scales, i.e ${\cal P}_{\zeta}(k)\simeq 10^{-2}$.
Therefore it is crucial to ask whether such a large enhancement of power at
smaller scales might induce non-negligible loop corrections to the CMB scale
power spectrum at higher orders in perturbation theory. Recently such
calculations were carried out perturbatively in Refs. [111, 112, 113, 114,
115]. These papers find that the one-loop corrections to the CMB scale power
spectrum can become significant if ${\cal P}_{\zeta}(k)\gtrsim{\cal
O}(10^{-2})$. This appears to rule out the formation of an interesting
abundance of PBHs in single field inflationary models. However this conclusion
is currently the subject of debate. It has been argued in Refs. [117, 116,
118] that the loop corrections are negligible if the transition from USR to
the subsequent attractor phase is smooth enough. Nevertheless we stress that
the amplitude of the small scale power spectrum required to form an
interesting abundance of PBHs depends on the PDF of the perturbations. The
standard value, ${\cal P}_{\zeta}(k)\simeq 10^{-2}$, assumes the PDF is
Gaussian. This amplitude, and therefore the size of the one-loop corrections
to the CMB scale power spectrum, will be different for the non-Gaussian tail
usually generated by stochastic effects.

##  6 Conclusions

PBHs can form due to the gravitational collapse of large fluctuations, in the
non-perturbative tail of the PDF. An accurate calculation of the full PDF of
the perturbations is therefore crucial to calculate their abundance.
Stochastic inflation is a powerful framework for computing the cosmological
correlators non-perturbatively. Using the stochastic $\delta{\cal N}$
formalism, the full PDF can be calculated from the first-passage statistics of
the number of e-folds, ${\cal N}$, during inflation. However to correctly
account for the back-reaction effect of small scale (UV) fluctuations,
$\varphi_{k}$, on the long wavelength coarse-grained (IR) field, $\Phi$, it is
essential to compute the noise matrix elements accurately. Since most single
field inflationary potentials with a PBH-forming feature violate the slow-roll
conditions, a precise calculation of the stochastic noise matrix elements
beyond slow roll is required. In this paper we have done this, both
analytically and numerically.

After a brief overview of single-field inflationary dynamics beyond slow roll
in Sec. 2, we set up the relevant equations underlying the stochastic
inflation formalism in Sec. 3. There are two key steps to using the stochastic
inflation formalism to calculate the full PDF of fluctuations in slow-roll
violating, PBH-producing, models:

  1. 1.

compute the statistics of both field and momentum-induced noise terms
$\\{\xi_{\phi},\xi_{\pi}\\}$,

  2. 2.

set up the Langevin equations (or, the corresponding adjoint Fokker-Planck
equation) without ignoring the inflaton IR momentum $\Pi$,

We have addressed the first issue here and will focus on the second in a
forthcoming paper [78].

In Sec. 4 we computed the matrix elements, $\Sigma_{ij}$, defined in Eq.
(3.15), which characterise the statistics of the field and momentum noise
terms. First, in Sec. 4.1, we derived expressions for $\Sigma_{ij}$ for
featureless potentials where the slow-roll conditions
$\epsilon_{H},\,\eta_{H}\ll 1$ remain valid until almost the end of inflation.
We compared the results of our analytical calculations, Eqs. (4.23)-(4.25),
and numerical calculations (shown in Fig. 4) for the KKLT potential, Eq.
(4.26), with the corresponding estimates under the pure de Sitter
approximation, Eqs. (4.17)-(4.19). We found that the dS approximation
underestimates the momentum induced noise terms, $\Sigma_{\phi\pi}$ and
$\Sigma_{\pi\pi}$, by several orders of magnitude, even for a slow-roll
potential.

In Sec. 4.2, we calculated the noise matrix elements for single field
inflationary potentials with a slow-roll violating, PBH-forming feature. For
the numerical calculations we used the modified KKLT potential featuring a
tiny Gaussian bump, Eq. (4.28), as a proto-typical single-field PBH-forming
potential. This potential has a sharp transition from the CMB scale SR-I phase
to the subsequent near-USR phase, as shown in Fig. 5. Our results, plotted in
Figs. 6 and 7, show that following the transition, $\Sigma_{ij}$ falls
exponentially and the momentum induced noise terms dominate the field noise
with the hierarchy $\Sigma_{\pi\pi}>|{\rm
Re}(\Sigma_{\phi\pi})|>\Sigma_{\phi\phi}$. Subsequently, the noise terms
return back to their original hierarchy, before growing and tending to
constant values.

To understand the asymptotic behaviour of the noise terms, we calculated the
noise matrix elements analytically using several approximations. Firstly we
treated the sharp transition between the SR-I phase and the subsequent near-
USR phase as instantaneous, and assumed the second slow-roll parameter
$\eta_{H}$ to be piece-wise constant. By solving the Mukhanov-Sasaki equation,
Eq. (4.7), analytically for a constant $\eta_{H}$ (and hence a constant
$\nu$), and applying the Israel junction matching conditions across the
transition, we computed the noise matrix elements shown in the right panel of
Fig. 8. We found that the behaviour of the noise terms post transition is
governed by a single parameter, namely the transition strength, ${\cal A}$,
which is defined as the difference between the values of the second slow-roll
parameter $\eta_{H}$ post- and pre-transition as given in Eq. (4.38). This
analytical computation based on an instantaneous transition
$\nu:\nu_{1}\to\nu_{2}$ captures the key features of the noise matrix elements
for potentials with a smooth feature (see Eqs. (4.72)-(4.74)).

We also compared our calculations with those for an instantaneous transition
using the pure dS approximation, i.e. $\nu_{1}=\nu_{2}=3/2$, which was carried
out in Ref. [95] for the Starobinsky model [107], see the left panel of Fig.
8. We found that the dS approximation underestimates the noise terms not only
in the SR-I phase (as mentioned before), but also a long time after the
transition. However, the pure dS-transition estimates are a good approximation
to the behaviour of the noise terms immediately after the transition.
Furthermore, for potentials with a pure ‘flat’ feature (as shown in Fig. 2)
rather than a bump, the dS-transition approximations work quite well.

In our analytical solutions of the MS equation, we focused on a single sharp
transition, T-I. In this case the effective mass term $z^{\prime\prime}/z$
remains almost constant throughout the USR, T-I and CR phases, as can be seen
in Fig. 5. Therefore the expression for the mode functions remains the same
after the second transition, due to the Wands duality as discussed in Sec.
4.2. This is a common characteristic of a broad class of single field
inflationary models with a PBH forming feature (see Ref. [84]). Our analytical
scheme can be extended to situations where the effective mass term
$z^{\prime\prime}/z$ undergoes two or more sharp transitions. We provide the
relevant analytical expressions for the mode functions in this case in App. E.

We conclude that in order to accurately determine the PDF of the curvature
perturbation, $P[\zeta]$, beyond slow roll, one must solve the adjoint Fokker-
Planck Equation (3.19) using the correct asymptotic forms of the noise matrix
elements given in Eqs. (4.72)-(4.74). Our upcoming paper [78] will be
dedicated to developing analytical and semi-analytical techniques to solve the
adjoint FPE with the knowledge of $\Sigma_{ij}$ obtained here. While numerical
simulations of the Langevin equations can be carried out in full generality,
they are often quite time-consuming, and demand large computational resources.
Furthermore, the analytical approach will allow us to calculate the asymptotic
behaviour of the PDF and study the effects of the noise terms, $\Sigma_{ij}$,
the potential, $V(\phi)$, and the boundary conditions on the PDF separately.
It is therefore complementary to the fully numerical simulations of the
Langevin equations discussed in Ref. [65, 69, 70, 73].

## Acknowledgements

We would like to thank David Wands and Antonio Riotto for helpful comments.
SSM, EJC and AMG are supported by a STFC Consolidated Grant [Grant No.
ST/T000732/1], and EJC by a Leverhulme Research Fellowship [RF-2021 312]. For
the purpose of open access, the authors have applied a CC BY public copyright
licence to any Author Accepted Manuscript version arising.

Data Availability Statement: This work is entirely theoretical and has no
associated data.

##  Appendix A Mukhanov-Sasaki equation in spatially flat gauge

In this Appendix we outline the derivation of the Mukhanov-Sasaki equation in
a spatially flat gauge, which we solved for in Sec. 4. Starting from the
action of a canonical scalar field minimally coupled to gravity

$S[g_{\mu\nu},\phi]=\int\,{\rm
d}^{4}x\,\sqrt{-g}\,\left(\,\frac{m_{p}^{2}}{2}\,R-\frac{1}{2}\,\partial_{\mu}\phi\,\partial_{\nu}\phi\,g^{\mu\nu}-V(\phi)\right)\leavevmode\nobreak\
,$

we consider linear field fluctuations
$\phi(t,\vec{x})=\bar{\phi}(t)+\delta\phi(t,\vec{x})$ and the linearly
perturbed ADM metric in the spatially flat gauge

${\rm d}s^{2}=-\alpha^{2}\,{\rm d}t^{2}+a^{2}(t)\,\delta_{ij}\left({\rm
d}x^{i}+\beta^{i}\,{\rm d}t\right)\left({\rm d}x^{j}+\beta^{j}\,{\rm
d}t\right)\,,$

where $\alpha=1+\delta\alpha$ and $\beta^{i}$ are the lapse and shift
functions. Imposing the GR momentum and Hamiltonian constraints, one obtains
(see Ref. [119])

$\delta\alpha=\sqrt{\frac{\epsilon_{H}}{2}}\,\frac{\delta\phi}{m_{p}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\partial_{i}\beta^{i}=-\epsilon_{H}\,\frac{{\rm d}}{{\rm
d}t}\left(\frac{1}{\sqrt{2\epsilon_{H}}}\,\frac{\delta\phi}{m_{p}}\right)\leavevmode\nobreak\
.$

Incorporating the above expressions into the action, expanding around the
background, the quadratic action for $\delta\phi$ fluctuations (in the
spatially flat gauge) is [119, 30]

$S[\delta\phi]=\frac{1}{2}\,\int\,{\rm d}t{\rm
d}^{3}{\vec{x}}\leavevmode\nobreak\
a^{3}\left[\dot{\delta\phi}^{2}-\frac{(\partial_{i}\delta\phi)^{2}}{a^{2}}-\left(\frac{{\rm
d}^{2}V(\phi)}{{\rm
d}\phi^{2}}+2\,\epsilon_{H}\,H^{2}\,\left(2\,\eta_{H}-\epsilon_{H}-3\right)\right)\delta\phi^{2}\right]\,.$
(A.1)

The Euler-Lagrange equation for $\delta\phi$ is given by

$\ddot{\delta\phi}+3\,H\,\dot{\delta\phi}-\frac{\nabla^{2}}{a^{2}}\,\delta\phi-a^{2}H^{2}\left(2-\epsilon_{H}-\frac{1}{H^{2}}\,\frac{{\rm
d}^{2}V(\phi)}{{\rm
d}\phi^{2}}-\frac{2\dot{\bar{\phi}}}{m_{p}^{2}H^{3}}\,\frac{{\rm
d}V(\phi)}{{\rm
d}\phi}-\frac{\dot{\bar{\phi}}^{2}}{m_{p}^{4}H^{4}}\,V\right)\,\delta\phi=0\leavevmode\nobreak\
.$ (A.2)

With the change of variables $v=a\,\delta\phi$, the Fourier modes $v_{k}$
satisfy the Mukhanov-Sasaki equation

$v_{k}^{{}^{\prime\prime}}+\left(k^{2}-\frac{z^{\prime\prime}}{z}\
\right)\,v_{k}=0\leavevmode\nobreak\ ,$

with

$\displaystyle\frac{z^{\prime\prime}}{z}$ $\displaystyle\equiv
a^{2}H^{2}\left(2-\epsilon_{H}-\frac{1}{H^{2}}\,\frac{{\rm d}^{2}V(\phi)}{{\rm
d}\phi^{2}}-\frac{2\dot{\bar{\phi}}}{m_{p}^{2}H^{3}}\,\frac{{\rm
d}V(\phi)}{{\rm
d}\phi}-\frac{\dot{\bar{\phi}}^{2}}{m_{p}^{4}H^{4}}\,V\right)\,,$ (A.3)
$\displaystyle=(aH)^{2}\left[2+2\epsilon_{H}-3\eta_{H}+2\epsilon_{H}^{2}+\eta_{H}^{2}-3\epsilon_{H}\eta_{H}-\frac{1}{aH}\,\eta^{\prime}_{H}\right]\,.$
(A.4)

##  Appendix B Analytical solution of the Mukhanov-Sasaki equation

For the featureless slow-roll potentials that we study in Sec. 4.1, $\nu^{2}$
is greater than or equal to $9/4$ and effectively constant. In this case the
Mukhanov-Sasaki (MS) Eq. (4.1) can be written as a Bessel equation with
constant $\nu$, which can be solved analytically (see Refs. [120, 121]). In
this Appendix we present this solution in terms of both Hankel functions (App.
B.1) and Bessel functions (App. B.2).

For the analytical treatment, assuming
$(aH)^{-2}\,z^{\prime\prime}/z=\nu^{2}-1/4$ to be a constant, the Mukhanov-
Sasaki (MS) Eq. (4.1) can be written as

$\frac{{\rm d}^{2}v_{k}}{{\rm
d}\tau^{2}}+\left[\,k^{2}-\frac{\nu^{2}-\frac{1}{4}}{\tau^{2}}\,\right]v_{k}=0\leavevmode\nobreak\
,$ (B.1)

using the new time variable, $T$, defined in Eq. (4.6)

$T=-k\tau=\frac{k}{a\,H}\,.$

All modes undergo Hubble-exit at $T=1$, with sub (super)-Hubble scales
corresponding to $T\gg(\ll)1$. In terms of this new time variable, the MS
equation takes the form

$\frac{{\rm d}^{2}v_{k}}{{\rm
d}T^{2}}+\left[\,1-\frac{\nu^{2}-\frac{1}{4}}{T^{2}}\,\right]v_{k}=0\,.$

Using the variable redefinition $F=v_{k}/\sqrt{T}$, this equation can be
transformed into the more familiar Bessel equation:

$\frac{{\rm d}^{2}F}{{\rm d}T^{2}}+\frac{1}{T}\,\frac{{\rm d}F}{{\rm
d}T}+\left[\,1-\frac{\nu^{2}}{T^{2}}\,\right]F=0\leavevmode\nobreak\ .$ (B.2)

The general solution to Eq. (B.2) (when $\nu$ is not an integer) can be
written either as a linear combination of Hankel functions of the first and
second kind $\\{H_{\nu}^{(1)}(T),\,H_{\nu}^{(2)}(T)\\}$ or as a linear
combination of positive and negative order ($\pm\nu$) Bessel functions of the
first kind $\\{J_{-\nu}(T),\,J_{\nu}(T)\\}$. The functions are related by
[120]

$H_{\nu}^{(1,2)}(T)=\frac{\pm J_{-\nu}(T)\mp e^{\mp
i\pi\nu}\,J_{\nu}(T)}{i\sin{(\pi\nu)}}\,.$ (B.3)

###  B.1 In terms of Hankel functions

The general solution to the Bessel Eq. (B.2) in terms of the Hankel functions
is given by

$F(T)=C_{1}\,H_{\nu}^{(1)}(T)\,+\,C_{2}\,H_{\nu}^{(2)}(T)\leavevmode\nobreak\
,$ (B.4)

where the coefficients $C_{1}$ and $C_{2}$ are fixed by initial/boundary
conditions. Hence the solution to the MS equation can be written as

$v_{k}(T)=\sqrt{T}\,\left[\,C_{1}\,H_{\nu}^{(1)}(T)\,+\,C_{2}\,H_{\nu}^{(2)}(T)\,\right]\,.$
(B.5)

In the sub-Hubble limit, $T\gg 1$, the Hankel functions take the form

$\displaystyle
H_{\nu}^{(1)}(T)\bigg{|}_{T\to\infty}\simeq\sqrt{\frac{2}{\pi}}\,\frac{1}{\sqrt{T}}\,e^{iT}\,e^{-i\left(\nu+\frac{1}{2}\right)\frac{\pi}{2}}\,,$
(B.6) $\displaystyle
H_{\nu}^{(2)}(T)\bigg{|}_{T\to\infty}\simeq\sqrt{\frac{2}{\pi}}\,\frac{1}{\sqrt{T}}\,e^{-iT}\,e^{i\left(\nu+\frac{1}{2}\right)\frac{\pi}{2}}\,,$
(B.7)

while in the super-Hubble limit, $T\ll 1$, the Hankel functions take the form

$\displaystyle H_{\nu}^{(1)}(T)\bigg{|}_{T\to
0}\simeq\sqrt{\frac{2}{\pi}}\,e^{-i\frac{\pi}{2}}\,2^{\nu-\frac{3}{2}}\,\frac{\Gamma(\nu)}{\Gamma(\frac{3}{2})}\,T^{-\nu}\,,$
(B.8) $\displaystyle H_{\nu}^{(2)}(T)\bigg{|}_{T\to
0}\simeq\leavevmode\nobreak\
-\sqrt{\frac{2}{\pi}}\,e^{-i\frac{\pi}{2}}\,2^{\nu-\frac{3}{2}}\,\frac{\Gamma(\nu)}{\Gamma(\frac{3}{2})}\,T^{-\nu}\,.$
(B.9)

The Bunch-Davies conditions, Eq. (4.5), for the mode functions take the form

$v_{k}(T)\bigg{|}_{T\to\infty}\to\frac{1}{\sqrt{2k}}\,e^{iT}=\sqrt{T}\,C_{1}H_{\nu}^{(1)}(T)\bigg{|}_{T\to\infty}\leavevmode\nobreak\
,$

which yields

$C_{1}=\frac{1}{\sqrt{2k}}\,\sqrt{\frac{\pi}{2}}\,e^{i\left(\nu+\frac{1}{2}\right)\frac{\pi}{2}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
C_{2}=0\leavevmode\nobreak\ ,$

and hence the final expression for the mode functions becomes

$v_{k}(T)=e^{i\left(\nu+\frac{1}{2}\right)\frac{\pi}{2}}\,\sqrt{\frac{\pi}{2}}\,\frac{1}{\sqrt{2k}}\,\sqrt{T}\,H_{\nu}^{(1)}(T)\,.$
(B.10)

###  B.2 In terms of Bessel functions

The general solution to the Bessel equation, Eq. (B.2), in terms of the Bessel
functions of the first kind of order $\pm\nu$ is given by

$F(T)=\sqrt{\frac{\pi}{2}}\,\frac{1}{\sin{(\pi\nu)}}\,\left[\,C_{+}\,J_{-\nu}(T)\,+\,C_{-}\,J_{\nu}(T)\,\right]\leavevmode\nobreak\
,$ (B.11)

where the coefficients $C_{+}$ and $C_{-}$ are again to be fixed by
initial/boundary conditions. Hence the solution to MS Eq. (4.1) can be written
as

$v_{k}(T)=\sqrt{\frac{\pi}{2}}\,\frac{1}{\sin{(\pi\nu)}}\,\sqrt{T}\,\left[\,C_{+}\,J_{-\nu}(T)\,+\,C_{-}\,J_{\nu}(T)\,\right]\,.$
(B.12)

Imposing Bunch-Davies initial conditions, Eq. (4.5), we get

$C_{+}=-i\,\frac{1}{\sqrt{2k}}\,e^{i\left(\nu+\frac{1}{2}\right)\frac{\pi}{2}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ {\rm
and}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
C_{-}=i\,\frac{1}{\sqrt{2k}}\,e^{i\left(\nu+\frac{1}{2}\right)\frac{\pi}{2}}\,e^{-i\pi\nu}\leavevmode\nobreak\
,$

and hence the final expression for the mode functions becomes

$v_{k}(T)=-i\,\sqrt{\frac{\pi}{2}}\,\frac{e^{i\left(\nu+\frac{1}{2}\right)\frac{\pi}{2}}}{\sin{(\pi\nu)}}\,\frac{1}{\sqrt{2k}}\,\sqrt{T}\left[\,J_{-\nu}(T)\,-e^{-i\pi\nu}\,J_{\nu}(T)\,\right]\leavevmode\nobreak\
.$ (B.13)

With the help of Eq. (B.3), we see by equating Eqs. (B.5) and (B.12), that the
relation between the Hankel coefficients $\\{C_{1},\,C_{2}\\}$ and Bessel
coefficients $\\{C_{+},\,C_{-}\\}$ is given by

$C_{1}=i\sqrt{\frac{\pi}{2}}\,\left[\,\frac{C_{+}+e^{-i2\pi\nu}C_{-}}{1-e^{-i2\pi\nu}}\,\right]\,;\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
C_{2}=i\sqrt{\frac{\pi}{2}}\,\left[\,\frac{C_{+}+e^{i\pi\nu}C_{-}}{1-e^{-i2\pi\nu}}\,\right]e^{-i2\pi\nu}\,.$
(B.14)

In this work, we express analytical solutions of the MS equation in terms of
Hankel functions. However, our results can alternatively be easily expressed
in terms of the Bessel functions by using Eqs. (B.3) and (B.14).

##  Appendix C Super-Hubble expansion of the noise matrix elements

The full expressions for the noise matrix elements are given in Eqs.
(4.42)-(4.44). Since we need to evaluate $\Sigma_{ij}$ in the super-Hubble
limit with $T=\sigma\ll 1$, here we provide expressions for the noise terms,
derived from the mode functions, $v_{k}(T)$, in Eq. (4.20), as a series
expansion in $T$, up to ${\cal O}(T^{4})$, for constant $\nu$:

$\displaystyle\Sigma_{\phi\phi}$ $\displaystyle=$ $\displaystyle
2^{2\left(\nu-\frac{3}{2}\right)}\,\left[\frac{\Gamma(\nu)}{\Gamma(3/2)}\right]^{2}\left(\frac{H}{2\pi}\right)^{2}T^{2\left(-\nu+\frac{3}{2}\right)}\left[1+\frac{1}{2(-1+\nu)}\,T^{2}\right.$
(C.1)
$\displaystyle\left.+\frac{\nu-3/2}{8(-2+\nu)(-1+\nu)^{2}}\,T^{4}\right.$
$\displaystyle\left.+\left(2+\frac{T^{2}}{(1-\nu^{2})}\right)\frac{2^{-2\nu}\Gamma[-\nu]}{\Gamma[\nu]}\cos{(\pi\nu)}\,T^{2\nu}+{\cal
O}(T^{4+2\nu})\right]\,,$

$\displaystyle{\rm Re}(\Sigma_{\pi\phi})$ $\displaystyle=$
$\displaystyle-2^{2\left(\nu-\frac{3}{2}\right)}\,\left[\frac{\Gamma(\nu)}{\Gamma(3/2)}\right]^{2}\left(\frac{H}{2\pi}\right)^{2}T^{2\left(-\nu+\frac{3}{2}\right)}\left[\left(\frac{3}{2}-\nu\right)+\frac{5-2\nu}{4(-1+\nu)}\,T^{2}\right.$
(C.2)
$\displaystyle\left.+\frac{-7+2\nu}{16(-2+\nu)(-1+\nu)^{2}}\left(\frac{3}{2}-\nu\right)\,T^{4}\right.$
$\displaystyle\left.-\left(3+\frac{5}{2(-1+\nu^{2})}T^{2}\right)\frac{2^{-2\nu}}{\Gamma(1+\nu)^{2}}\pi\nu\cot{(\pi\nu)}\,T^{2\nu}+{\cal
O}(T^{4+2\nu})\right]\,,$

$\displaystyle\Sigma_{\pi\pi}$ $\displaystyle=$ $\displaystyle
2^{2\left(\nu-\frac{3}{2}\right)}\,\left[\frac{\Gamma(\nu)}{\Gamma(3/2)}\right]^{2}\left(\frac{H}{2\pi}\right)^{2}T^{2\left(-\nu+\frac{3}{2}\right)}\left[\left(-\nu+\frac{3}{2}\right)^{2}+\left(-\nu+\frac{3}{2}\right)\left(\frac{-2\nu+7}{4(-\nu+1)}\right)\,T^{2}\right.$
(C.3)
$\displaystyle\left.+\left(\frac{131-2\nu(83-34\nu+4\nu^{2})}{64(-\nu+2)(-1+\nu)^{2}}\right)\,T^{4}\right.$
$\displaystyle\left.-\left(\frac{-4\nu^{2}+9}{2}+\frac{21+4\nu^{2}}{4(-1+\nu^{2})}\,T^{2}\right)\frac{2^{-2\nu}}{\Gamma[1+\nu]^{2}}\pi\nu\cot{(\pi\nu)}\,T^{2\nu}+{\cal
O}(T^{4+2\nu})\right]\,.$

The above expressions for $\Sigma_{ij}$ are valid for any value of $\nu$, and
accurately reproduce the de Sitter results for $\nu=3/2$. We have verified
that there are no higher order corrections in the dS limit.

##  Appendix D Functional form of $z^{\prime\prime}/z$ during sharp
transitions

In this Appendix we derive the analytic expressions for $z^{\prime\prime}/z$
for the instantaneous transitions that we use in Sec. 4.2.2. Under the quasi-
dS approximation, $\epsilon_{H}\simeq 0$, the effective mass term in the MS
Eq. (4.1) becomes

$\frac{z^{\prime\prime}}{z}\,\tau^{2}\simeq
2-3\,\eta_{H}+\eta_{H}^{2}+\tau\,\frac{{\rm d}\eta_{H}}{{\rm d}\tau}\,,$ (D.1)

where this final form of $z^{\prime\prime}/z$ depends upon the expression for
the second slow-roll parameter $\eta_{H}(\tau)$. In the following, we will
assume that $\eta_{H}$ is piece-wise constant and makes instantaneous, but
finite transitions. We will begin with the simplest case where $\eta_{H}$
makes only one transition and later generalise this to the case of two or more
successive transitions.

###  D.1 Single instantaneous transition

Suppose that during inflation, the second SR parameter makes a sharp
transition from $\eta_{H}=\eta_{1}\to\eta_{2}$ at time $\tau=\tau_{1}$. In
this case $\eta_{H}$ can be written as

$\eta_{H}(\tau)=\eta_{1}+\left(\eta_{2}-\eta_{1}\right)\,\Theta(\tau-\tau_{1})\leavevmode\nobreak\
.$ (D.2)

Using Eq. (D.1), we obtain the following expression for $z^{\prime\prime}/z$

$\frac{z^{\prime\prime}}{z}\,\tau^{2}={\cal
A}\,\tau\,\delta_{D}(\tau-\tau_{1})+\nu_{1}^{2}-\frac{1}{4}+\left(\nu_{2}^{2}-\nu_{1}^{2}\right)\,\Theta(\tau-\tau_{1})\,,$
(D.3)

with ${\cal A}=\eta_{2}-\eta_{1}$ and
$\nu_{1,2}^{2}=(9/4)-3\,\eta_{1,2}+\eta_{1,2}^{2}$. For $\eta_{1}=0$ and
$\eta_{2}=+3$, the system reduces to the case of pure de Sitter SR $\to$ USR
transition, where $\nu_{1}^{2}=\nu_{2}^{2}=9/4$ and ${\cal A}=+3$.

###  D.2 Two successive instantaneous transitions

If the second SR parameter makes another sharp transition from
$\eta_{H}=\eta_{2}\to\eta_{3}$ at time $\tau=\tau_{2}$, then $\eta_{H}$ can be
written as

$\eta_{H}(\tau)=\eta_{1}+\left(\eta_{2}-\eta_{1}\right)\,\Theta(\tau-\tau_{1})+\left(\eta_{3}-\eta_{2}\right)\,\Theta(\tau-\tau_{2})\,.$
(D.4)

Once again using Eq. (D.1), we now obtain

$\frac{z^{\prime\prime}}{z}\,\tau^{2}=\nu_{1}^{2}-\frac{1}{4}+\left(\nu_{2}^{2}-\nu_{1}^{2}\right)\,\Theta(\tau-\tau_{1})+\left(\nu_{3}^{2}-\nu_{2}^{2}\right)\,\Theta(\tau-\tau_{2})+{\cal
A}_{1}\,\tau\,\delta_{D}(\tau-\tau_{1})+{\cal
A}_{2}\,\tau\,\delta_{D}(\tau-\tau_{2})\,,$ (D.5)

with ${\cal A}_{1,2}=\eta_{2,3}-\eta_{1,2}$ and
$\nu_{1,2,3}^{2}=(9/4)-3\,\eta_{1,2,3}+\eta_{1,2,3}^{2}$.

###  D.3 Generalising to multiple instantaneous transitions

If the inflaton potential exhibits a number of tiny features/modulations, then
the second SR parameter $\eta_{H}$ might undergo a number of successive
transitions before the end of inflation. Assuming each transition to be
instantaneous for ease of analytical treatment, we can write the following
general expression for $\eta_{H}$

$\eta_{H}(\tau)=\eta_{1}+\left(\eta_{2}-\eta_{1}\right)\,\Theta(\tau-\tau_{1})+\left(\eta_{3}-\eta_{2}\right)\,\Theta(\tau-\tau_{2})+...+\left(\eta_{n+1}-\eta_{n}\right)\,\Theta(\tau-\tau_{n})\,,$
(D.6)

where ‘$n$’ is the total number of instantaneous transitions occurring at
times $\tau_{1},\tau_{2},....,\tau_{n}$. In this case, the expression for
$z^{\prime\prime}/z$ becomes

$\begin{split}\frac{z^{\prime\prime}}{z}\,\tau^{2}=\nu_{1}^{2}-\frac{1}{4}+\left(\nu_{2}^{2}-\nu_{1}^{2}\right)\,\Theta(\tau-\tau_{1})+\left(\nu_{3}^{2}-\nu_{2}^{2}\right)\,\Theta(\tau-\tau_{2})+...+\left(\nu_{n+1}^{2}-\nu_{n}^{2}\right)\,\Theta(\tau-\tau_{n})\\\
+{\cal A}_{1}\,\tau\,\delta_{D}(\tau-\tau_{1})+{\cal
A}_{2}\,\tau\,\delta_{D}(\tau-\tau_{2})+...+{\cal
A}_{n}\,\tau\,\delta_{D}(\tau-\tau_{n})\,,\end{split}$ (D.7)

where ${\cal A}_{n}=\eta_{n+1}-\eta_{n}$ and
$\nu_{n}^{2}=(9/4)-3\,\eta_{n}+\eta_{n}^{2}$.

In compact notation, we can write the general expressions for $\eta_{H}$ and
$z^{\prime\prime}/z$ for ‘$n$’ successive instantaneous transitions as

$\displaystyle\eta_{H}(\tau)$ $\displaystyle=$
$\displaystyle\eta_{1}\,+\,\sum_{i=1}^{n}\,\left[\,\eta_{i+1}-\eta_{i}\,\right]\,\Theta(\tau-\tau_{i})\,,$
(D.8) $\displaystyle\frac{z^{\prime\prime}}{z}\,\tau^{2}$ $\displaystyle=$
$\displaystyle\nu_{1}^{2}-\frac{1}{4}\,+\,\sum_{i=1}^{n}\,\left[\,\nu_{i+1}^{2}-\nu_{i}^{2}\,\right]\,\Theta(\tau-\tau_{i})\,+\,{\cal
A}_{i}\,\tau\,\delta_{D}(\tau-\tau_{i})\,,$ (D.9)

where ${\cal A}_{i}=\eta_{i+1}-\eta_{i}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\nu_{i}^{2}=(9/4)-3\,\eta_{i}+\eta_{i}^{2}$ .

##  Appendix E Noise matrix elements for two successive instantaneous
transitions

In this Appendix, we present the full calculations for the mode functions in a
closed form for the case of two successive instantaneous transitions during
inflation. This generalises the results presented in Case 1 (Pure dS limit)
and Case 2 (transition between two different values $\nu_{1}\to\nu_{2}$) of
Sec. 4.2.2.

###  E.1 Pure dS limit

For two successive instantaneous transitions SR $\to$ USR $\to$ SR in the pure
dS limit, the effective mass term in the MS equation takes the form

$\frac{z^{\prime\prime}}{z}\,\tau^{2}={\cal
A}\,\tau\,\left[\,\delta_{D}(\tau-\tau_{1})-\delta_{D}(\tau-\tau_{2})\,\right]\,+\,2\,,$
(E.1)

where $\tau_{1}$ is the transition time from SR to USR and $\tau_{2}$ is the
transition time from USR back to SR. ${\cal A}=+3$ is the strength of the
transition, as discussed before

The MS (complex) mode functions generalise those derived in Eq. (4.46) and are
given by

$\displaystyle
v_{k}(T)\equiv\begin{cases}v_{k}^{E}(T)=\frac{1}{\sqrt{2k}}\,\left(1+\frac{i}{T}\right)\,e^{i\,T}\,,&T_{2}<T_{1}<T\leavevmode\nobreak\
,\\\
v_{k}^{I}(T)=\frac{1}{\sqrt{2k}}\,\left[\,C_{1}^{I}\,\left(1+\frac{i}{T}\right)\,e^{i\,T}\,+\,C_{2}^{I}\,\left(1-\frac{i}{T}\right)\,e^{-i\,T}\,\right]\,,&T_{2}<T<T_{1}\leavevmode\nobreak\
,\\\
v_{k}^{L}(T)=\frac{1}{\sqrt{2k}}\,\left[\,C_{1}^{L}\,\left(1+\frac{i}{T}\right)\,e^{i\,T}\,+\,C_{2}^{L}\,\left(1-\frac{i}{T}\right)\,e^{-i\,T}\,\right]\,,&T<T_{2}<T_{1}\leavevmode\nobreak\
,\end{cases}$ (E.2)

and the derivatives of the mode functions generalise those derived in Eq.
(4.47) and are given by

$\displaystyle\frac{{\rm d}v_{k}}{{\rm d}T}\equiv\begin{cases}\frac{{\rm
d}v_{k}^{E}}{{\rm
d}T}=\frac{1}{\sqrt{2k}}\,\left[-\frac{1}{T}+i\,\left(1-\frac{1}{T^{2}}\right)\right]\,e^{i\,T}\,,&T_{2}<T_{1}<T\leavevmode\nobreak\
,\\\ \frac{{\rm d}v_{k}^{I}}{{\rm
d}T}=\frac{1}{\sqrt{2k}}\,\left[\,C_{1}^{I}\,\left[-\frac{1}{T}+i\,\left(1-\frac{1}{T^{2}}\right)\right]\,e^{i\,T}\,+\,C_{2}^{I}\,\left[-\frac{1}{T}-i\,\left(1-\frac{1}{T^{2}}\right)\right]\,e^{-i\,T}\,\right]\,,&T_{2}<T<T_{1}\leavevmode\nobreak\
,\\\ \frac{{\rm d}v_{k}^{L}}{{\rm
d}T}=\frac{1}{\sqrt{2k}}\,\left[\,C_{1}^{L}\,\left[-\frac{1}{T}+i\,\left(1-\frac{1}{T^{2}}\right)\right]\,e^{i\,T}\,+\,C_{2}^{L}\,\left[-\frac{1}{T}-i\,\left(1-\frac{1}{T^{2}}\right)\right]\,e^{-i\,T}\,\right]\,,&T<T_{2}<T_{1}\leavevmode\nobreak\
.\end{cases}$ (E.3)

Here the superscripts ‘$E$’, ‘$I$’ and ‘$L$’ stand for early, intermediate and
late respectively. After implementing the Israel junction matching conditions,
we obtain the final expressions for the mode functions

$\displaystyle
C_{1}^{L}=\frac{d_{3}\,b_{4}-d_{4}\,b_{3}}{a_{3}\,b_{4}-a_{4}\,b_{3}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
C_{2}^{L}=\frac{d_{4}\,a_{3}-d_{3}\,a_{4}}{a_{3}\,b_{4}-a_{4}\,b_{3}}\leavevmode\nobreak\
,$ (E.4)

where

$\displaystyle a_{3}$ $\displaystyle=T_{2}+i\,,$ (E.5) $\displaystyle b_{3}$
$\displaystyle=\left(T_{2}-i\right)\,e^{-i\,2T_{2}}\,,$ (E.6) $\displaystyle
a_{4}$ $\displaystyle=\left({\cal
A}-1\right)\,T_{2}+i\,\left(T_{2}^{2}-1+{\cal A}\right)\,,$ (E.7)
$\displaystyle b_{4}$ $\displaystyle=\ \left({\cal
A}-1\right)\,T_{2}-i\,\left(T_{2}^{2}-1+{\cal A}\right)\,e^{-i\,2T_{2}}\,,$
(E.8) $\displaystyle d_{3}$
$\displaystyle=\left(T_{2}+i\right)\,C_{1}^{I}+\left(T_{2}-i\right)\,e^{-i\,2T_{2}}\,C_{2}^{I}\,,$
(E.9) $\displaystyle d_{4}$
$\displaystyle=\left[T_{2}-i\,\left(T_{2}^{2}-1\right)\right]\,C_{1}^{I}+\left[T_{2}+i\,\left(T_{2}^{2}-1\right)\right]\,C_{2}^{I}\leavevmode\nobreak\
,$ (E.10)

and the intermediate transition coefficients $C_{1}^{I}$ and $C_{2}^{I}$ are
just $\alpha_{k}$ and $\beta_{k}$, respectively, derived previously in Eqs.
(4.59) and (4.60).

###  E.2 $\nu_{1}\to\nu_{2}\to\nu_{3}$

In case of two successive instantaneous transitions from
$\nu_{1}\to\nu_{2}\to\nu_{3}$, the effective mass term in the MS equation
takes the form

$\frac{z^{\prime\prime}}{z}\,\tau^{2}=\nu_{1}^{2}-\frac{1}{4}+\left(\nu_{2}^{2}-\nu_{1}^{2}\right)\,\Theta(\tau-\tau_{1})+\left(\nu_{3}^{2}-\nu_{2}^{2}\right)\,\Theta(\tau-\tau_{2})+{\cal
A}_{1}\,\tau\,\delta_{D}(\tau-\tau_{1})+{\cal
A}_{2}\,\tau\,\delta_{D}(\tau-\tau_{2})\,,$ (E.11)

where $\tau_{1}$ is the transition time from $\nu_{1}\to\nu_{2}$ and
$\tau_{2}$ is the transition time from $\nu_{2}\to\nu_{3}$. ${\cal
A}_{1},\,{\cal A}_{2}$ are the strengths of the first and second transitions
respectively.

The MS (complex) mode functions and their derivatives generalise those derived
in Eqs. (4.61) and (4.62) to yield

$\displaystyle
v_{k}(T)\equiv\begin{cases}v_{k}^{E}(T)=\frac{1}{\sqrt{2k}}\,\sqrt{\frac{\pi}{2}}\,\sqrt{T}\,H_{\nu_{1}}^{(1)}(T)\,e^{i\left(\nu_{1}+\frac{1}{2}\right)\frac{\pi}{2}}\,,&T_{2}<T_{1}<T\leavevmode\nobreak\
,\\\
v_{k}^{I}(T)=\sqrt{T}\,\left[\,C_{1}^{I}\,H_{\nu_{2}}^{(1)}(T)\,+\,C_{2}^{I}\,H_{\nu_{2}}^{(2)}(T)\,\right]\,,&T_{2}<T<T_{1}\leavevmode\nobreak\
,\\\
v_{k}^{L}(T)=\sqrt{T}\,\left[\,C_{1}^{L}\,H_{\nu_{3}}^{(1)}(T)\,+\,C_{2}^{L}\,H_{\nu_{3}}^{(2)}(T)\,\right]\,,&T<T_{2}<T_{1}\leavevmode\nobreak\
,\end{cases}$ (E.12)

$\displaystyle\frac{{\rm d}v_{k}}{{\rm
d}T}\equiv\begin{cases}\frac{1}{\sqrt{2k}}\sqrt{\frac{\pi}{2}}e^{i\left(\nu_{1}+\frac{1}{2}\right)\frac{\pi}{2}}\frac{1}{\sqrt{T}}\left[\left(\frac{1}{2}-\nu_{1}\right)\,H_{\nu_{1}}^{(1)}(T)+T\,H_{\nu_{1}-1}^{(1)}(T)\right]\,,&T_{2}<T_{1}<T\leavevmode\nobreak\
,\\\
\frac{1}{\sqrt{T}}\left[C_{1}^{I}\left(\left(\frac{1}{2}-\nu_{2}\right)H_{\nu_{2}}^{(1)}(T)+T\,H_{\nu_{2}-1}^{(1)}(T)\right)+C_{2}^{I}\left(\left(\frac{1}{2}-\nu_{2}\right)H_{\nu_{2}}^{(2)}(T)+T\,H_{\nu_{2}-1}^{(2)}(T)\right)\right]\,,&T_{2}<T<T_{1}\leavevmode\nobreak\
,\\\
\frac{1}{\sqrt{T}}\left[C_{1}^{L}\left(\left(\frac{1}{2}-\nu_{3}\right)H_{\nu_{3}}^{(1)}(T)+T\,H_{\nu_{3}-1}^{(1)}(T)\right)+C_{2}^{L}\left(\left(\frac{1}{2}-\nu_{3}\right)H_{\nu_{3}}^{(2)}(T)+T\,H_{\nu_{3}-1}^{(2)}(T)\right)\right]\,,&T<T_{2}<T_{1}\end{cases}$
(E.13)

The superscripts ‘$E$’, ‘$I$’ and ‘$L$’ again stand for early, intermediate
and late respectively. After implementing the Israel junction matching
conditions, we find the final expressions for the mode functions

$\displaystyle
C_{1}^{L}=\frac{d_{3}\,b_{4}-d_{4}\,b_{3}}{a_{3}\,b_{4}-a_{4}\,b_{3}}\leavevmode\nobreak\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
C_{2}^{L}=\frac{d_{4}\,a_{3}-d_{3}\,a_{4}}{a_{3}\,b_{4}-a_{4}\,b_{3}}\leavevmode\nobreak\
,$ (E.14)

where

$\displaystyle a_{3}$ $\displaystyle=H_{\nu_{3}}^{(1)}(T_{2})\,,$ (E.15)
$\displaystyle b_{3}$ $\displaystyle=H_{\nu_{3}}^{(2)}(T_{2})\,,$ (E.16)
$\displaystyle a_{4}$ $\displaystyle=\left(\,\frac{1}{2}-\nu_{3}-{\cal
A}_{2}\,\right)\,H_{\nu_{3}}^{(1)}(T_{2})+T_{2}\,H_{\nu_{3}-1}^{(1)}(T_{2})\,,$
(E.17) $\displaystyle b_{4}$ $\displaystyle=\left(\,\frac{1}{2}-\nu_{3}-{\cal
A}_{2}\,\right)\,H_{\nu_{3}}^{(2)}(T_{2})+T_{2}\,H_{\nu_{3}-1}^{(2)}(T_{2})\,,$
(E.18) $\displaystyle d_{3}$