# Numerical simulations of inflationary dynamics: slow-roll and beyond

Siddharth S. Bhatt  Swagat S. Mishra  Soumen Basak  and Surya N. Sahoo

###### Abstract

Cosmic inflation is a period of rapid accelerated expansion of space in the
very early universe. During inflation, vacuum quantum fluctuations are
amplified and stretched to cosmological scales which seed the fluctuations in
the cosmic microwave background as well as the large-scale structure of our
universe. Large quantum fluctuations may lead to the formation of primordial
black holes (PBHs) in the post-inflationary universe. Numerical simulations of
the inflationary dynamics are presented here for a single canonical scalar
field minimally coupled to gravity. We spell out the basic equations governing
the inflationary dynamics in terms of cosmic time $t$ and define a set of
dimensionless variables convenient for numerical analysis. We then provide a
link to our simple numerical Python code on GitHub that can be used to
simulate the background dynamics as well as the evolution of linear
perturbations during inflation. The code computes both scalar and tensor power
spectra for a given inflaton potential $V(\phi)$. We discuss a concrete
algorithm to use the code for various purposes, especially for computing the
enhanced scalar power spectrum in the context of PBH formation. We intend to
extend the framework to simulate the dynamics of a number of different
quantities, including the computation of scalar-induced second-order tensor
power spectrum in the revised version of this manuscript in the near future.

##  1 Introduction

Cosmic inflation has emerged as the leading scenario for describing the very
early universe prior to the commencement of the radiative hot Big Bang Phase
[1, 2, 3, 4, 5, 6]. According to the inflationary paradigm, a transient epoch
of at least 60-70 e-folds of rapid accelerated expansion suffices in setting
natural initial conditions for the background space-time in the form of
spatial flatness as well as statistical homogeneity and isotropy on large
angular scales [2, 3, 4, 7]. Additionally, (and more significantly,) quantum
fluctuations during inflation naturally generate a spectrum of almost scale-
invariant initial scalar fluctuations which seed the temperature and
polarisation fluctuations in the Cosmic Microwave Background (CMB) Radiation,
and later, the formation of structure in the universe [8, 9, 10, 11, 7]. In
addition to scalar perturbations, quantum fluctuations during inflation also
create a spectrum of almost scale-invariant tensor perturbations which later
become gravitational waves [12, 13].

The simplest models of inflation comprising of a single scalar field, called
the ‘inflaton’, which is minimally coupled to gravity, makes several distinct
predictions [14] (i.e an almost scale-invariant, nearly Gaussian, and
adiabatic spectrum of scalar fluctuations) most of which have received
spectacular observational confirmation, particularly from the latest CMB
missions [15].

However, as mentioned earlier, inflation also generates tensor perturbations
that later constitute the relic gravitational wave background (GW) which
imprints a distinct signature on the CMB power spectrum in the form of the
B-mode polarization [15]. The amplitude of these relic GWs provides us
information about the inflationary energy scale while their spectrum enables
us to access general properties of the epoch of reheating, being exceedingly
sensitive to the post-inflationary equation of state [13, 16]. The amplitude
of inflationary tensor fluctuations, relative to that of scalar fluctuations,
is usually characterised by the tensor-to-scalar ratio $r$. Different models
of inflation predict different values for $r$ which is sensitive to the
gradient of the inflaton potential $V_{,\phi}(\phi)=\frac{dV(\phi)}{d\phi}$
relative to its height $V(\phi)$. Convex potentials predict large values for
$r$, while concave potentials predict relatively small values of $r$. While
the spectrum of inflationary tensor fluctuations has not yet been observed,
current CMB observations are able to place an upper bound on the tensor-to-
scalar ratio on large angular scales. In particular, the latest CMB
observations of BICEP/Keck [17], combined with those of the PLANCK mission
[15], place the strong upper bound $r\leq 0.036$ (at $95\%$ confidence).

This most recent upper bound on $r$ has important consequences for single
field canonical inflation. In particular, given $r\leq 0.036$, all
monotonically increasing convex potentials, including the whole family of
monomial potentials $V(\phi)\propto\phi^{p}$, are completely ruled out in the
canonical framework. Among these strongly disfavoured models are the simplest
classic inflaton potentials $\frac{1}{2}m^{2}\phi^{2}$ and $\lambda\phi^{4}$.
Instead, the observational upper bound on $r$ appears to favour
asymptotically-flat potentials possessing one or two plateau-like wings; see
[19, 18]. Current observational data lead to a scenario in which the inflaton
$\phi$ slowly rolls down a shallow potential $V(\phi)$ thereby giving rise to
a quasi-de Sitter early stage of near-exponential expansion. A thorough
analysis of the inflationary phase-space dynamics $\\{\phi,\dot{\phi}\\}$ for
plateau potentials shows [20] that a large range of initial conditions leads
to adequate inflation in these models.

However, it is important to stress that the CMB window constitutes only a tiny
part of the observationally available field space between the Hubble-exit of
the largest scales in the sky and the smallest scale at the end of inflation.
Consequently, a substantial period of the inflationary dynamics corresponding
to potentially interesting small-scale primordial physics (which accounts
roughly to the last 40–50 e-folds of accelerated expansion during inflation)
remains observationally unexplored, being inaccessible to the CMB and LSS
observations. Any departure from the slow-roll regime, that might be triggered
by a change in the dynamics of the inflaton field, would lead to interesting
observational consequences on small-scales. In particular, the presence of a
feature at intermediate field values that might lead to large enough
amplification of the small-scale scalar fluctuations, would facilitate the
formation of Primordial Black Holes upon the Hubble re-entry of these modes
during the post-inflationary epochs.

Primordial Black Holes (PBHs) are extremely interesting compact objects which
might have been formed from the collapse of large density fluctuations in the
early universe [21, 22, 23, 24] and they constitute a potential candidate for
dark matter [27, 26, 25, 28, 29, 30]. Seeds for such large fluctuations can be
generated during inflation, as mentioned above. For instance, a feature in the
inflaton potential in the form of a flat inflection point can further slow
down the already slowly rolling inflaton field substantially, leading to an
enhancement of the primordial scalar power $P_{\zeta}$. A number of different
features for enhancing small-scale power during inflation have been proposed
in the recent years [31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44,
45, 46, 47, 48, 49]. Hence PBHs (and the associated induced relic GWs) are
excellent probes of the small-scale primordial physics.

In this preliminary version of our paper, we discuss a simple code developed
for numerical simulations of the inflationary dynamics. We introduce the
relevant dimensionless variables used in our numerical analysis and provide a
link to the code in our GitHub account. We also discuss how to use the code in
various scenarios, which include phase-space analysis of inflationary initial
conditions, inflationary background dynamics and determining scalar and tensor
power spectra both under slow-roll approximation and beyond. The latter case
has important implications for PBH formation and we discuss how to use the
code to simulate the Mukhanov-Sasaki equation mode by mode. We also discuss a
number of important future directions that are to be included in the
forthcoming version of our paper. The primary version of our numerical
framework is quite simple and less compact. It is intended to provide a
pedagogical guideline for researchers who are relatively new to numerical
simulations of inflation. In the forthcoming version of our work, we will
introduce a much more compact numerical framework that we are currently
working on which will incorporate additional new features.

This paper is organised as follows: we begin with a brief introduction of the
inflationary scalar field dynamics in section 2 and quantum fluctuations in
section 3. We then proceed to discuss numerical simulations of the background
dynamics in section 4. This also includes studying the scalar and tensor
fluctuations under the slow-roll approximations. Section 5 is dedicated to
studying the inflationary quantum fluctuations by numerically solving the
Mukhanov-Sasaki equation, and its application to inflaton potentials
possessing a slow-roll violating feature. We also mention a number of future
extensions of our numerical set-up in section 6, before concluding with a
discussion section.

We work in the units $c,\hbar=1$. The reduced Planck mass is defined to be
$m_{p}\equiv 1/\sqrt{8\pi G}=2.43\times 10^{18}~{}{\rm GeV}$. We assume the
background universe to be described by a spatially flat Friedmann-Lemaitre-
Robertson-Walker (FLRW) metric with signature $(-,+,+,+)$.

##  2 Inflationary Dynamics

The Action for a scalar field which couples minimally to gravity has the
following general form

$S[\phi]=\int d^{4}x\,\sqrt{-g}\;{\cal L}(F,\phi),$ (2.1)

where the Lagrangian density ${\cal L}(\phi,F)$ is a function of the field
$\phi$ and the kinetic term

$F=\frac{1}{2}\partial_{\mu}\phi\;\partial^{\mu}\phi.$ (2.2)

Varying (2.1) with respect to $\phi$ results in the equation of motion

$\frac{\partial{\cal
L}}{\partial\phi}-\left(\frac{1}{\sqrt{-g}}\right)\partial_{\mu}\left(\sqrt{-g}\frac{\partial{\cal
L}}{\partial\left(\partial_{\mu}\phi\right)}\right)=0.$ (2.3)

The energy-momentum tensor associated with the scalar field is

$T^{\mu\nu}=\left(\frac{\partial{\cal L}}{\partial
F}\right)\,\left(\partial^{\mu}\phi\;\partial^{\nu}\phi\right)-g^{\mu\nu}\,{\cal
L}~{}.$ (2.4)

Specializing to a spatially flat FRW universe and a homogeneous scalar field,
one gets

$ds^{2}=-dt^{2}+a^{2}(t)\left[dx^{2}+dy^{2}+dz^{2}\right]\,,$ (2.5)
$T^{\mu}_{\;\>\;\nu}=\mathrm{diag}\left(-\rho_{{}_{\phi}},p_{{}_{\phi}},p_{{}_{\phi}},p_{{}_{\phi}}\right),$
(2.6)

where the energy density, $\rho_{{}_{\phi}}$, and pressure, $p_{{}_{\phi}}$,
are given by

$\displaystyle\rho_{{}_{\phi}}$ $\displaystyle=\left(\frac{\partial{\cal
L}}{\partial F}\right)\,(2\,F)-{\cal L}\,,$ (2.7) $\displaystyle
p_{{}_{\phi}}$ $\displaystyle={\cal L}\,,$ (2.8)

and $F=-({\dot{\phi}}^{2}/2)$. The evolution of the scale factor $a(t)$ is
governed by the Friedmann equations:

$\displaystyle\left(\frac{\dot{a}}{a}\right)^{2}\equiv H^{2}$ $\displaystyle=$
$\displaystyle\frac{1}{3m_{p}^{2}}\,\rho_{{}_{\phi}},$ (2.9)
$\displaystyle\frac{\ddot{a}}{a}\equiv\dot{H}+H^{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{6m_{p}^{2}}\,\left(\rho_{{}_{\phi}}+3\,p_{{}_{\phi}}\right),$
(2.10)

where $H\equiv\dot{a}/a$ is the Hubble parameter and $\rho_{{}_{\phi}}$
satisfies the conservation equation

${\dot{\rho}_{{}_{\phi}}}=-3\,H\left(\rho_{{}_{\phi}}+p_{{}_{\phi}}\right)~{}.$
(2.11)

In the standard single field inflationary paradigm, inflation is sourced by a
minimally coupled canonical scalar field $\phi$ with a suitable potential
$V(\phi)$ (see figure 1). For such a canonical scalar field

${\cal L}(F,\phi)=-F-V(\phi),$ (2.12)

Figure 1: This figure schematically depicts a prototype inflation potential
$V(\phi)$, plotted in solid green curve. The ‘CMB Window’ represents field
values corresponding to the Hubble-exit epochs of scales
$k\in\left[0.0005,0.5\right]~{}{\rm Mpc}^{-1}$ that are observable by the
latest CMB missions.

Substituting (2.12) into (2.7) and (2.8), we find

$\displaystyle\rho_{{}_{\phi}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}{\dot{\phi}}^{2}+\;V(\phi),$ $\displaystyle
p_{{}_{\phi}}$ $\displaystyle=$
$\displaystyle\frac{1}{2}{\dot{\phi}}^{2}-\;V(\phi),~{}~{}$ (2.13)

consequently the two Friedmann equations (2.9), (2.10) and the equation (2.11)
become

$\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.14) $\displaystyle\dot{H}\equiv\frac{\ddot{a}}{a}-H^{2}$ $\displaystyle=$
$\displaystyle-\frac{1}{2m_{p}^{2}}\,\dot{\phi}^{2}\,,$ (2.15)
$\displaystyle{\ddot{\phi}}+3\,H{\dot{\phi}}+V_{,\phi}(\phi)$ $\displaystyle=$
$\displaystyle 0\,.$ (2.16)

The epoch of inflation at any time $t<t_{\rm end}$ is conveniently marked by
the number of e-folds before the end of inflation

$N_{e}=\log_{e}{\frac{a_{\rm end}}{a(t)}}=\int_{t}^{t_{\rm
end}}H(t^{\prime})dt^{\prime}\,,$ (2.17)

where $H(t)$ is the Hubble parameter during inflation. $a(t)$ and $a_{\rm
end}$ denote the scale factor at time $t$ and at the end of inflation
respectively. Typically a period of quasi-de Sitter inflation lasting for at
least 60-70 e-folds is required in order to address the problems of the
standard hot Big Bang model. We denote $N_{*}$ as the number of e-folds
(before the end of inflation) when the CMB pivot scale
$k_{*}=(aH)_{*}=0.05~{}\rm Mpc^{-1}$ left the comoving Hubble radius during
inflation. For convenience, we have chosen $N_{*}=60$ for the most part of
this work, although the exact value of $N_{*}$ depends upon the particular
detail of reheating history.

The quasi-de Sitter like phase corresponds to the inflaton field rolling
slowly down the potential $V(\phi)$. This slow-roll regime111It is well known
that the slow-roll phase of the inflation is actually a local attractor for
many different models of inflation, see [50, 20] and the references therein.
of inflation, ensured by the presence of the Hubble friction term in the
equation (2.16), is usually characterised by the first two kinematical Hubble
slow-roll parameters $\epsilon_{H}$, $\eta_{H}$, defined by [7]

$\displaystyle\epsilon_{H}$ $\displaystyle=$
$\displaystyle-\frac{\dot{H}}{H^{2}}=\frac{1}{2m_{p}^{2}}\,\frac{\dot{\phi}^{2}}{H^{2}},$
(2.18) $\displaystyle\eta_{H}$ $\displaystyle=$
$\displaystyle-\frac{\ddot{\phi}}{H\dot{\phi}}=\epsilon_{H}+\frac{1}{2\epsilon_{H}}\,\frac{d\epsilon_{H}}{dN_{e}}~{},$
(2.19)

where the slow-roll regime of inflation corresponds to

$\epsilon_{H},~{}\eta_{H}\ll 1~{}.$ (2.20)

The slow-roll regime is also often characterised by the dynamical potential
slow-roll parameters [7], defined by

$\epsilon_{{}_{V}}=\frac{m_{p}^{2}}{2}\left(\frac{V_{,\phi}}{V}\right)^{2}~{},~{}~{}\eta_{{}_{V}}=m_{p}^{2}\,\left(\frac{V_{,\phi\phi}}{V}\right)~{}.$
(2.21)

For small values of these parameters $\epsilon_{H},\,\eta_{H}\ll 1$, one finds
$\epsilon_{H}\simeq\epsilon_{{}_{V}}$ and
$\eta_{H}\simeq\eta_{{}_{V}}-\epsilon_{{}_{V}}$. Using the definition of
Hubble parameter, $H=\dot{a}/a$, we have
$\ddot{a}/a=\dot{H}+H^{2}=H^{2}(1+\dot{H}/H^{2})$. From the expression for
$\epsilon_{H}$ in (2.18), it is easy to see that

$\frac{\ddot{a}}{a}=\big{(}1-\epsilon_{H}\big{)}\,H^{2}~{}.$ (2.22)

Which implies that the universe accelerates, ${\ddot{a}}>0$, when
$\epsilon_{H}<1$. Using equation (2.14), the expression for $\epsilon_{H}$ in
(2.18) reduces to $\epsilon_{H}\simeq\frac{3}{2}\frac{\dot{\phi}^{2}}{V}$ when
${\dot{\phi}}^{2}\ll V$.

Figure 2: This figure schematically illustrates the evolution of the comoving
Hubble radius $(aH)^{-1}$ with scale factor. During inflation $(aH)^{-1}$
decreases which causes physical scales to exit the Hubble radius. After
inflation ends $(aH)^{-1}$ increases, and physical scales begin to re-enter
the Hubble radius. The CMB pivot scale, as used by the Planck mission, is set
at $k_{*}=0.05~{}{\rm Mpc}^{-1}$ and has been depicted by the dashed blue
line. The ‘CMB Window’ corresponds to the scales
$k\in\left[0.0005,0.5\right]~{}{\rm Mpc}^{-1}$ that are observable by the
latest CMB missions. Note that $(aH)^{-1}\propto a$ during the radiative
regime and $(aH)^{-1}\propto a^{-1}$ during inflation.

##  3 Quantum fluctuations during inflation

In the standard scenario of a minimally coupled single canonical scalar field
$\phi$ as the inflaton, two gauge-independent massless fields, one scalar, and
one transverse traceless tensor, get excited during inflation and receive
quantum fluctuations that are correlated over super-Hubble scales [51] at late
times.

### Scalar fluctuations

The evolution of the scalar degree of freedom, called the curvature
perturbation222Note that the comoving curvature perturbation ${\cal R}$ is
related to the curvature perturbation on uniform-density hypersurfaces,
$\zeta$, and both are equal during slow-roll inflation as well as on super-
Hubble scales, $k\ll aH$, in general (see [7]). $\zeta$ is described by the
following quadratic Action [7, 52]

$\boxed{S_{(2)}[\zeta]=\frac{1}{2}\int{\rm d}\tau{\rm
d}x^{3}\,z^{2}\,\left[\,(\zeta^{\prime})^{2}-(\partial_{i}\zeta)^{2}\,\right]}~{},$
(3.1)

which upon the change of variable

$v\equiv z\,{\zeta}~{},\quad\mbox{with}\quad
z=am_{p}\sqrt{2\epsilon_{H}}=a\,\frac{\dot{\phi}}{H}\,,$ (3.2)

takes the form

$S_{(2)}\left[v\right]=\frac{1}{2}\int{\rm d}\tau{\rm
d}x^{3}\left[\left(v^{\prime}\right)^{2}-\left(\partial_{i}v\right)^{2}+\frac{z^{\prime\prime}}{z}v^{2}\right]~{},$
(3.3)

where the $(^{\prime})$ denotes derivative with respect to conformal time
$\tau=\int\frac{dt}{a(t)}$ ($\simeq\frac{-1}{aH}$ for quasi-de Sitter
expansion). The variable $v$, which itself is a scalar quantum field like
$\zeta$, is called the Mukhanov-Sasaki variable in literature. Its Fourier
modes $v_{k}$ satisfy the famous Mukhanov-Sasaki equation given by [53, 54]

$\boxed{v^{\prime\prime}_{k}\,+\,\left(k^{2}-\frac{z^{\prime\prime}}{z}\right)v_{k}=0}~{},$
(3.4)

where the effective mass term is given by the following exact expression [55]

$\displaystyle\boxed{\frac{z^{\prime\prime}}{z}=(aH)^{2}\left(2-\epsilon_{1}+\frac{3}{2}\epsilon_{2}+\frac{1}{4}\epsilon_{2}^{2}-\frac{1}{2}\epsilon_{1}\epsilon_{2}+\frac{1}{2}\epsilon_{2}\epsilon_{3}\right)}~{},$
(3.5)
$\displaystyle\Rightarrow\boxed{\frac{z^{\prime\prime}}{z}=(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]}~{},$
(3.6)

with $\epsilon_{1}=\epsilon_{H}$ and where

$\epsilon_{n+1}=-\frac{d\ln{\epsilon_{n}}}{dN_{e}}~{}$ (3.7)

are the ‘Hubble flow’ parameters. Given a mode $k$, at sufficiently early
times when it is sub-Hubble i.e $k\gg aH$, we can assume $v$ to be in the
Bunch-Davies vacuum [56] satisfying

$v_{k}\rightarrow\frac{1}{\sqrt{2k}}e^{-ik\tau}~{}.$ (3.8)

During inflation as the comoving Hubble radius falls (see figure 2), modes
start becoming super-Hubble i.e $k\ll aH$ and equation (3.4) dictates that
$|v_{k}|\propto z$ and hence $\zeta_{k}$ approaches a constant value. By
solving the Mukhanov-Sasaki equation we can estimate the dimensionless
primordial power spectrum of $\zeta$ using the following relation [51]

$\boxed{P_{\zeta}\equiv\frac{k^{3}}{2\pi^{2}}\,|{\zeta}_{k}|^{2}\Big{|}_{k\ll
aH}=\frac{k^{3}}{2\pi^{2}}\,\frac{|v_{k}|^{2}}{z^{2}}\Big{|}_{k\ll aH}}~{}.$
(3.9)

During slow-roll inflation, the factor
$\frac{z^{\prime\prime}}{z}=\frac{\nu^{2}-0.25}{\tau^{2}}$ with $\nu\approx
1.5+\epsilon_{H}+\frac{\dot{\epsilon}_{H}}{2H\epsilon_{H}}$. Solving the
Mukhanov-Sasaki equation with suitable Bunch-Davies vacuum conditions leads to
the famous slow-roll approximation formula [7]

$\boxed{P_{\zeta}=\frac{1}{8\pi^{2}}\left(\frac{H}{m_{p}}\right)^{2}\frac{1}{\epsilon_{H}}}~{}.$
(3.10)

Note that one could also directly try to solve for the fourier modes of the
comoving curvature perturbation $\zeta$ (instead of the Mukhanov-Sasaki
variable $v$) which satisfies the equation

$\boxed{{\zeta}^{\prime\prime}_{k}+2\left(\frac{z^{\prime}}{z}\right){\zeta}^{\prime}_{k}+k^{2}{\zeta}_{k}=0}$
(3.11)

and implement the corresponding Bunch-Davies initial conditions for
${\zeta}_{k}$. The friction term in equation (3.11) is given by

$\boxed{\frac{z^{\prime}}{z}=aH\,(1+\epsilon_{H}-\eta_{H})}~{}.$ (3.12)

Before moving forward, we stress that the slow-roll regime of inflation
necessarily requires both the slow-roll parameters to be small i.e
$\epsilon_{H}\ll 1$ and $\eta_{H}\ll 1$. Violation of either of these
conditions invalidates the above analytical treatment. When either of the
slow-roll conditions is violated, which is the situation in the context of
primordial black hole formation, a more accurate determination of $P_{\zeta}$
is provided by solving the Mukhanov-Sasaki equation (3.4) numerically. The
computation of power spectrum when the slow-roll approximation is violated
will be our primary focus.

### Tensor fluctuations

The corresponding quadratic Action for tensor fluctuations is given by [7, 52]

$\boxed{S_{(2)}[\gamma_{ij}]=\frac{1}{2}\int{\rm d}\tau{\rm
d}^{3}x\,a^{2}\left(\frac{m_{p}}{2}\right)^{2}\left[(\gamma_{ij}^{\prime})^{2}-(\partial\gamma_{ij})^{2}\right]}\,.$
(3.13)

The Mukhanov-Sasaki variable for tensor fluctuations are defined as

$\frac{m_{p}}{2}\,a\,\gamma_{ij}\equiv\left(\begin{smallmatrix}h_{+}&h_{\times}&0\\\
h_{\times}&-h_{+}&0\\\ 0&0&0\end{smallmatrix}\right)$ (3.14)

or,

$\frac{m_{p}}{2}\,a\,\gamma_{ij}\equiv\sum_{s=+,\times}\Pi_{ij}^{s}\ h_{s}\,,$
(3.15)

where $\Pi^{+}$ and $\Pi^{\times}$ are the 2 polarization modes, written as

$\Pi^{+}_{ij}=\begin{pmatrix}1&0&0\\\ 0&-1&0\\\
0&0&0\end{pmatrix}\quad\quad{\rm
and}\quad\quad\Pi^{\times}_{ij}=\begin{pmatrix}0&1&0\\\ 1&0&0\\\
0&0&0\end{pmatrix}\,.$ (3.16)

Thus,

$\gamma_{ij}=\frac{2}{m_{p}}\sum_{s=+,\times}\Pi_{ij}^{s}\,\frac{h_{s}}{a}\,.$
(3.17)

The evolution equation for the mode functions (by dropping the ‘s’ subscript
and remembering that it is summed over for 2 polarization states) is given by

$\boxed{h_{k}^{\prime\prime}+\left(k^{2}-\frac{a^{\prime\prime}}{a}\right)h_{k}=0}\,.$
(3.18)

The subsequent computation of tensor power spectrum

$P_{T}(k)\equiv
2\times\frac{k^{3}}{2\pi^{2}}\,\frac{|h_{k}|^{2}}{a^{2}}\Big{|}_{k\ll aH}\,,$
(3.19)

under quasi-de Sitter approximation leads to [7]

$\boxed{P_{T}(k)=\frac{2}{\pi^{2}}\left(\frac{H}{m_{p}}\right)^{2}}\,.$ (3.20)

Note that, unlike the Mukhanov-Sasaki equation (3.4) for scalar fluctuations,
the tensor mode equation (3.18) does not depend upon $z$, rather it depends
only upon the scale factor $a$. Hence, as long as the quasi-de Sitter
approximation is valid, i.e $\epsilon_{H}\ll 1$, power spectrum of tensor
fluctuations does not get affected by an appreciable amount even if slow-roll
is violated. Although, this statement is true only at linear order in
perturbation theory. Tensor fluctuations at second order in perturbation
theory can be induced by large first-order scalar fluctuations [57, 58, 59,
60] (also see [61] and references therein).

###  3.1 Large scale primordial fluctuations

On large cosmological scales which are accessible to CMB observations, the
scalar power spectrum typically takes the form of a power law represented by

$P_{\zeta}(k)=A_{{}_{S}}\left(\frac{k}{k_{*}}\right)^{n_{{}_{S}}-1},$ (3.21)

where $A_{{}_{S}}=P_{\zeta}(k_{*})$ is the amplitude of the scalar power
spectrum at the pivot scale $k=k_{*}$, given by333Note that in general, $k$
may correspond to any observable CMB scale in the range
$k\in\left[0.0005,0.5\right]~{}{\rm Mpc}^{-1}$. However, in order to derive
constraints on the inflationary observables $\\{n_{{}_{S}},r\\}$, we mainly
focus on the CMB pivot scale, namely $k\equiv k_{*}=0.05~{}{\rm Mpc}^{-1}$.

$\boxed{A_{{}_{S}}=\frac{1}{8\pi^{2}}\left(\frac{H}{m_{p}}\right)^{2}\frac{1}{\epsilon_{H}}\,\bigg{|}_{\phi=\phi_{*}}}~{},$
(3.22)

where $\phi_{*}$ is the value of the inflaton field at the epoch of Hubble
exit of the CMB pivot scale $k_{*}$. The scalar spectral tilt $n_{{}_{S}}$, in
the slow-roll regime is given by [7]

$\displaystyle\boxed{n_{{}_{S}}-1\equiv\frac{d\,\mathrm{ln}P_{\zeta}}{d\,\mathrm{ln}k}=2\eta_{H}-4\epsilon_{H}}~{}.$
(3.23)

Similarly the tensor power spectrum, in the slow-roll limit, is represented by

$P_{T}(k)=A_{{}_{T}}\left(\frac{k}{k_{*}}\right)^{n_{{}_{T}}},$ (3.24)

with the amplitude of tensor power spectrum at the CMB pivot scale is given by
[7, 51]

$\boxed{A_{{}_{T}}\equiv
P_{T}(k_{*})=\frac{2}{\pi^{2}}\left(\frac{H}{m_{p}}\right)^{2}\bigg{|}_{\phi=\phi_{*}}}~{},$
(3.25)

and the tensor spectral index (with negligible running) is given by

$\boxed{n_{{}_{T}}=-2\,\epsilon_{H}}~{}.$ (3.26)

The tensor-to-scalar ratio $r$ is defined by

$\boxed{r\equiv\frac{A_{{}_{T}}}{A_{{}_{S}}}=16\,\epsilon_{H}}~{},$ (3.27)

yielding the single field consistency relation

$\boxed{r=-8\,n_{{}_{T}}}~{}.$ (3.28)

Hence the slow-roll parameters $\epsilon_{H}$ and $\eta_{H}$ play an important
role in characterising the power spectra of scalar and tensor fluctuations
during inflation. Before going forward, we briefly discuss the implications of
the latest CMB observations for the slow-roll parameters as well as for other
relevant inflationary observables. In order to relate the CMB observables to
the inflaton potential $V(\phi)$, we work with the potential slow-roll
parameters defined in equation (2.21).

Consider a canonical scalar field minimally coupled to gravity and having the
potential

$V(\phi)=V_{0}\,f\left(\frac{\phi}{m_{p}}\right)~{}.$ (3.29)

The potential slow-roll parameters (2.21) are given by

$\displaystyle\epsilon_{{}_{V}}=\frac{m_{p}^{2}}{2}\left(\frac{f_{,\phi}}{f}\right)^{2}~{},$
(3.30)
$\displaystyle\eta_{{}_{V}}=m_{p}^{2}\left(\frac{f_{,\phi\phi}}{f}\right)~{}.$
(3.31)

In the slow-roll limit $\epsilon_{{}_{V}},\,\eta_{{}_{V}}\ll 1$, the scalar
power spectrum is given by the expression (3.21) with the amplitude of scalar
power at the CMB pivot scale $k\equiv k_{*}=0.05~{}{\rm Mpc}^{-1}$ expressed
as [7]

$A_{{}_{S}}\equiv
P_{\cal\zeta}(k_{*})\simeq\frac{1}{24\pi^{2}}\frac{V_{0}}{m_{p}^{4}}\frac{f\left(\phi_{k}\right)}{\epsilon_{{}_{V}}(\phi_{k})}\bigg{|}_{k=k_{*}}~{},$
(3.32)

and the scalar spectral index (with negligible running) is given by

$n_{{}_{S}}\simeq
1+2\,\eta_{{}_{V}}(\phi_{*})-6\,\epsilon_{{}_{V}}(\phi_{*})~{},$ (3.33)

Similarly the amplitude of tensor power spectrum at the CMB pivot scale is
given by

$A_{{}_{T}}\equiv
P_{T}(k_{*})=\frac{2}{\pi^{2}}\left(\frac{H}{m_{p}}\right)^{2}\bigg{|}_{\phi=\phi_{*}}\simeq\frac{2}{3\pi^{2}}\frac{V_{0}}{m_{p}^{4}}f\left(\phi_{*}\right)~{},$
(3.34)

and the tensor spectral index (3.26) becomes

$n_{{}_{T}}\simeq-2\,\epsilon_{{}_{V}}(\phi_{*})~{},$ (3.35)

and the tensor-to-scalar ratio (3.27) can be written as

$r\simeq 16\,\epsilon_{{}_{V}}(\phi_{*})~{},$ (3.36)

satisfying the single field consistency relation (3.28). From the CMB
observations of Planck 2018 [15], we have

$A_{{}_{S}}=2.1\times 10^{-9}~{},$ (3.37)

while the $2\sigma$ constraint on the scalar spectral index is given by

$n_{{}_{S}}\in[0.957,0.976]~{}.$ (3.38)

Similarly the constraint on the tensor-to-scalar ratio $r$, from the latest
combined observations of Planck 2018 [15] and BICEP/Keck [17], is given by

$r\leq 0.036~{},$ (3.39)

which translates into $A_{{}_{T}}\leq 3.6\times 10^{-2}\,A_{{}_{S}}$. Equation
(3.34) helps place the following upper bound on the inflationary Hubble scale
$H^{\rm inf}$ and the energy scale during inflation $E_{\inf}$

$\displaystyle H^{\rm inf}\leq 4.7\times 10^{13}~{}{\rm GeV}~{},$ (3.40)
$\displaystyle E_{\inf}\equiv\left[\sqrt{3}\,m_{p}\,H^{\rm
inf}\right]^{1/2}\leq 1.4\times 10^{16}~{}{\rm GeV}~{}~{}.$ (3.41)

Similarly the CMB bound on $r$ when combined with (3.36) translates into an
upper bound on the first slow-roll parameter

$\epsilon_{H}\simeq\epsilon_{{}_{V}}\leq 0.00225~{},$ (3.42)

rendering the tensor tilt from equation (3.35) to be negligibly small

$|n_{{}_{T}}|\leq 0.0045~{}.$ (3.43)

Given the upper limit on $\epsilon_{{}_{V}}$, using the CMB bound on
$n_{{}_{S}}$ from (3.38) in (3.33), we infer that the second slow-roll
parameter is negative and obtain interesting upper and lower limits on its
magnitude, given by

$|\eta_{H}|\in[0.0075,0.0215]~{}.$ (3.44)

The EOS $w_{\phi}$ of the inflaton field is given by

$w_{\phi}=\frac{\frac{1}{2}\dot{\phi}^{2}-V(\phi)}{\frac{1}{2}\dot{\phi}^{2}+V(\phi)}\simeq-1+\frac{2}{3}\epsilon_{{}_{V}}(\phi)\,,$
(3.45)

Therefore one finds from (3.42) the following constraint on the inflationary
EOS at the pivot scale

$w_{\phi}\leq-0.9985\,,$ (3.46)

implying that the expansion of the universe during inflation was near
exponential (quasi-de Sitter like).

Figure 3: This figure shows the Starobinsky potential (3.48) with CMB pivot
scale $k_{*}=0.05~{}{\rm Mpc}^{-1}$ by a blue color star as well as the CMB
window $k_{\rm CMB}\in\left[0.0005,0.5\right]~{}{\rm Mpc}^{-1}$ in grey color
shade in the field space. From this figure, it is clear that CMB window
constitutes only a tiny portion of the available field space in between
$\phi_{\rm CMB}$ and end of inflation $\phi_{e}$.

The CMB observations, in the context of single field slow-roll inflationary
paradigm, favours asymptotically-flat potentials (featuring either one or two
plateau wings) with $n_{{}_{S}}\simeq 0.965$ and $r\leq 0.036$. A typical
plateau-potential is demonstrated in figure 1 and this is the standard/vanilla
scenario. Given that power spectrum is almost scale-invariant with slightly
red tilt, i.e $n_{{}_{S}}-1\lesssim 0$, large-scale fluctuations are more
important while nothing drastic is expected to happen on smaller cosmological
scales that are super-Hubble at the end of inflation.

Before proceeding further to discuss small-scale inflationary fluctuations,
let us make our nomenclature concrete (which is consistent with the standard
nomenclature in the inflationary literature).

• Quasi-de Sitter inflation corresponds to the condition $\epsilon_{H}\ll 1$.
• Slow-roll inflation corresponds to $\epsilon_{H},\,\eta_{H}\ll 1$.

This distinction will be important for the rest of the discussions in this
paper. Under either of the aforementioned assumptions, the expression for
conformal time is given by

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

###  3.2 Small-scale primordial fluctuations

As mentioned above, the recent CMB observations support the scenario of the
inflaton field rolling slowly down an asymptotically-flat potential at the
time when the observable CMB scales made their Hubble exit during inflation.
However, the current CMB and LSS observations probe only about 7-8 e-folds of
inflation around the Hubble exit time of the CMB pivot scale. We explicitly
mention that CMB observations probe primordial fluctuations with comoving
scales $k_{\rm CMB}\in\left[0.0005,0.5\right]~{}{\rm Mpc}^{-1}$ (which
includes pivot scale $k_{*}=0.05~{}{\rm Mpc}^{-1}$) corresponding to multipole
$l\in\left[2,2500\right]$ in the angular sky. Additionally, Lyman-$\alpha$
forest observations enforce constraints on the primordial power spectrum upto
$k\simeq{\cal O}(1)~{}{\rm Mpc}^{-1}$ (see [29]).

Hence a large portion of evolution during the inflationary phase that accounts
roughly to about 50 e-folds of expansion, corresponding to scales smaller than
those probed by CMB, remains observationally inaccessible at present.
Consequently the associated dynamics of the inflaton field also remains
unprobed. For example, figure 3 demonstrates that the CMB window constitutes
only a tiny part of the observationally available field space between the
largest scales in the sky and the smallest scale at the end of inflation for
Starobinsky potential [1, 62]

$\boxed{V(\phi)=V_{0}\,\left(1-e^{-\frac{2}{\sqrt{6}}\,\frac{\phi}{m_{p}}}\right)^{2}}\,.$
(3.48)

Any deviation from the quasi-de Sitter expansion and/or departure from the
slow-roll regime $\epsilon_{{}_{H}},\eta_{{}_{H}}\ll 1$ that might be
triggered by a change in the dynamics of the inflaton field, would lead to
interesting observational consequences on small scales. In particular if the
inflaton potential possesses a near inflection point-like broad feature at
some intermediate field values, then the scalar quantum fluctuations
corresponding to scales becoming super-Hubble around the time when the
inflaton rolls past such features, might receive enough amplification to
facilitate the formation of Primordial Black Holes (PBHs) upon their Hubble
re-entry during the radiative epoch.

Figure 4: This figure schematically depicts a prototype plateau potential,
plotted in solid green curve. The ‘CMB Window’ represents field values
corresponding to the Hubble-exit epochs of scales
$k\in\left[0.0005,0.5\right]~{}{\rm Mpc}^{-1}$ that are observable by the
latest CMB missions. The potential exhibits a small-scale feature (shown in
the salmon colour shading) in the form of a flat inflection point-like segment
which results in ultra slow-roll (USR) inflation. After exiting the first
slow-roll phase (SR-I) near the CMB window, the inflaton enters into an USR
phase, during which the second slow-roll condition is violated, namely
$\eta_{H}\simeq+3$. This leads to an enhancement of power spectrum at small-
scales. Later, the inflaton emerges out of the USR to another slow-roll phase
(SR-II) before the end of inflation.

We note that, there are a plethora of possible features that would lead to
deviation from standard scale-invariant power spectrum [63]. However, in this
paper we will only focus on potentials with a tiny local bump/dip like feature
[44] in order to illustrate the efficiency of our numerical analysis. Our code
can be used to simulate a number of different types of features in the
inflaton potential.

PBH formation requires the enhancement of the inflationary power spectrum by
roughly a factor of $10^{7}$ within less than 40 e-folds of expansion (on
scales smaller than the pivot scale $k_{*}$) as illustrated in figure 5.
Therefore the quantity $\Delta\ln{\epsilon_{H}}/\Delta N$, and hence also
$|\eta_{H}|$, can grow to become of order unity, thereby violating the second
slow-roll condition [64]. In fact the second Hubble slow-roll parameter
$|\eta_{H}|$ becomes larger than unity even though $\epsilon_{H}$ itself
remains much smaller than unity. As a result, equation (3.22) can no longer be
trusted to compute the power spectrum and one must determine $P_{\cal R}$ by
numerically integrating the Mukhanov-Sasaki equation (3.4).

We proceed as follows. We first discuss the simulations of inflationary
background dynamics in section 4, where we also discuss how to generate phase-
space portrait $\\{\phi,\dot{\phi}\\}$ during inflation. In section 5.1, we
introduce our numerical scheme for studying quantum fluctuations during slow-
roll inflation. We work with convex as well as asymptotically-flat potentials.
In section 5.2, we apply our numerical scheme to potentials featuring a local
bump/dip like feature that facilitates the amplification of scalar power on
small primordial scales. We demonstrate that slow-roll formula (3.22)
underestimates both the location as well as the height of scalar power
spectrum ${\cal P}_{\zeta}(k)$ for both type of aforementioned features and
hence one must solve the Mukhanov-Sasaki equation (3.4) numerically to
estimate the power spectrum accurately. We also demonstrate that the growth of
the power spectrum obeys the steepest growth bounds discussed in [65, 66, 67,
68].

Figure 5: This figure schematically illustrates the typical amplification of
inflationary primordial power spectrum at smaller length scales required for
PBH formation.

##  4 Numerical analysis of inflationary background dynamics

A complete analysis of the inflationary background dynamics can be obtained
from the evolution of $\phi$, $\dot{\phi}$ and $H$. All of these quantities
can be simulated by numerically solving equations (2.14), (2.15) and (2.16).
The evolution of the scale factor follows directly from $H=\dot{a}/a$. Our
system is defined by the following set of equations (as a function of cosmic
time $t$)

$\displaystyle H^{2}$
$\displaystyle=\frac{1}{3m_{p}^{2}}\left[\frac{1}{2}{\dot{\phi}}^{2}+V(\phi)\right]\,,$
(4.1) $\displaystyle\dot{H}$ $\displaystyle=-\frac{\dot{\phi}^{2}}{2m_{p}}\,,$
(4.2) $\displaystyle\ddot{\phi}$
$\displaystyle=-3H\dot{\phi}-V_{,\phi}(\phi)\,,$ (4.3)

where the functional form of the potential $V(\phi)$ is given by the specific
inflationary model. However, the rest of the algorithm is largely model-
independent. We can re-write the potential as

$V(\phi)=V_{0}\,f(\phi)\,.$ (4.4)

In order to carry out numerical simulations, it is convenient to write down
the dynamical equations in terms of dimensionless variables (which also
ensures that we do not need to worry about keeping track of units).
Furthermore, it is important to re-scale the time variable by a factor $S$
which can be suitably chosen according to the energy scale of the
dynamics444Depending upon the potential, we usually choose the value of $S$ to
be in the range $S\in[10^{-5},\,10^{-3}]$.. Our primary dimensionless
variables are defined as

$\displaystyle T$ $\displaystyle=\left(t\,m_{p}\right)\,S\,,$ (4.5)
$\displaystyle x$ $\displaystyle=\frac{\phi}{m_{p}}\,,$ (4.6) $\displaystyle
y$ $\displaystyle=\left(\frac{\dot{\phi}}{m_{p}^{2}}\right)\,\frac{1}{S}\,,$
(4.7) $\displaystyle z$
$\displaystyle=\left(\frac{H}{m_{p}}\right)\,\frac{1}{S}\,,$ (4.8)
$\displaystyle A$ $\displaystyle=\left(a\,m_{p}\right)\,S\,.$ (4.9)

In terms of these variables, the dynamical equations (to be simulated) take
the form

$\displaystyle\frac{{\rm d}x}{{\rm d}T}$ $\displaystyle=y\,,$ (4.10)
$\displaystyle\frac{{\rm d}y}{{\rm d}T}$
$\displaystyle=-3\,z\,y-\frac{v_{0}}{S^{2}}\,f_{,x}(x)\,,$ (4.11)
$\displaystyle\frac{{\rm d}z}{{\rm d}T}$
$\displaystyle=-\frac{1}{2}\,y^{2}\,,$ (4.12) $\displaystyle\frac{{\rm
d}A}{{\rm d}T}$ $\displaystyle=A\,z\,.$ (4.13)

We can also define the dimensionless potential to be

$\frac{V(\phi)}{m_{p}^{4}}\equiv v_{0}\,f(x)=\frac{V_{0}}{m_{p}^{4}}\,f(x)\,.$
(4.14)

We can solve the aforementioned set of equations with appropriate initial
conditions. In our analysis, we use the odeint function provided in the
scipy.integration package. By incorporating initial conditions
$\\{x_{i},\,y_{i},\,z_{i},\,A_{i}\\}$ for the primary dynamical variables
$\\{x,\,y,\,z,\,A\\}$, we simulate their time evolution during inflation.
Accordingly, we determine the crucial derived/secondary (dimensionless)
dynamical variables from the primary ones by555Note that the observables
$A_{{}_{S}},\,A_{{}_{T}},\,n_{{}_{S}},\,n_{{}_{T}},\,{\rm and}\,r$ are related
to inflationary scalar and tensor fluctuations and the expressions given here
are under slow-roll approximations, i.e $\epsilon_{H},|\eta_{H}|\ll 1$, during
which they can be determined purely from the dynamics of background quantities
such as $H,\,\epsilon_{H},\,\eta_{H}$. Computation of inflationary power
spectra when slow-roll is violated is described in section 5.2.

$\displaystyle N$ $\displaystyle=\log{\frac{A}{A_{i}}}$ (4.15)
$\displaystyle\epsilon_{H}$
$\displaystyle=\frac{1}{2}\,\frac{y^{2}}{z^{2}}~{},$ $\displaystyle\eta_{H}$
$\displaystyle=-\frac{1}{yz}\frac{{\rm d}y}{{\rm d}T}\,,$ (4.16)
$\displaystyle A_{{}_{S}}$
$\displaystyle=\frac{1}{8\pi^{2}}\,\frac{\left(Sz\right)^{2}}{\epsilon_{H}}~{},$
$\displaystyle A_{{}_{T}}$
$\displaystyle=\frac{2}{\pi^{2}}\,\left(Sz\right)^{2}\,,$ (4.17)
$\displaystyle n_{{}_{S}}$ $\displaystyle=1+2\,\eta_{H}-4\,\epsilon_{H}~{},$
$\displaystyle n_{{}_{T}}$ $\displaystyle=-2\,\epsilon_{H}\,,$ (4.18)
$\displaystyle r$ $\displaystyle=16\,\epsilon_{H}\,.$ (4.19)

We define $N_{T}$ to be the number of e-folds of accelerated expansion
realised in between an arbitrary initial time and the end of inflation, which
is marked by $\epsilon_{H}=1$. We then define the more important quantity
$N_{e}=N_{T}-N$ as the number of e-folds before the end of inflation. Note
that $N_{e}=0$ at the end of inflation while $N_{e}>0$ at early times. This
will be our primary time variable against which we will be plotting the
dynamics of different inflationary observables. In order to realise adequate
amount of inflation, i.e. $N_{T}>60$, initial value of scalar field must be
large enough (and is model dependent). For most large field potentials, this
value is of the order $\phi_{i}\lesssim{\cal O}(10)\,m_{p}$.

Composing the code involves typing down the dimensionless equations in the
appropriate syntax and solving them by using an ODE solver. See our
supplementary Python code 666https://github.com/bhattsiddharth/NumDynInflation
for details. For a particular model of interest, we need to input the inflaton
potential in the form $\frac{V(\phi)}{m_{p}^{4}}=v_{0}\,f(x)$. Since the slow-
roll parameters and the duration $N_{T}$ do not strongly depend on $v_{0}$, we
can initially set its value to roughly $v_{0}=10^{-10}$. We can later adjust
the value of $v_{0}$ to yield the correct CMB normalised value of scalar power
spectrum (3.37) at the pivot scale $N=N_{*}$. One can proceed in the following
step-by-step algorithm.

Figure 6: Time evolution of the number of e-folds (scale factor in the
logarithm scale) of expansion of the universe is shown for Starobinsky
potential (3.48). For the most part of inflation, the expansion is almost
exponential (quasi-de Sitter) i.e $a\sim e^{Ht}$, leading to a rapid growth in
the number of e-folds within a small amount of time. While after the end of
inflation, the expansion is decelerated, leading to a much slower growth in
the scale factor.

  1. 1.

After setting the parameters of the potential and defining the function
$f(x)$, we need to incorporate initial conditions for the four primary
variables $\\{x,\,y,\,z,\,A\\}$. We enter appropriate initial conditions
$x_{i},\,y_{i}$ and $A_{i}$ in the following way. $A_{i}$ can be set
arbitrarily in a spatially flat universe, however depending upon the energy
scale of inflation, one can provide an appropriate value. We suggest a typical
$A_{i}=1\times 10^{-3}$, although its precise value does not affect the
dynamics. In regard to the initial value of $x$, we need to ensure that
$x_{i}$ is large enough (or small enough if we are working with symmetry
breaking hilltop type potentials) to yield adequate amount of inflation, i.e
$N_{T}\geq 70$. As mentioned before, the typical value for large field models
is $x_{i}\lesssim{\cal O}(10)$.

Since we will be mostly working with potentials that exhibit slow-roll
behaviour at initial times, and given that slow-roll trajectory is an
attractor in relatively large field models [20], we can safely set $y_{i}=0$,
as long as $x_{i}$ is large enough. One can also incorporate slow-roll initial
conditions from the beginning, namely
$y_{i}=-\frac{v_{0}}{3}\,\frac{f_{,x}}{Sz}$, as is usually done in
practice777Note that for phase-space analysis, we need to incorporate
arbitrary values of $x_{i},\,y_{i}$ (consistent with fixed initial $z_{i}$)
which may be away from the slow-roll trajectory as discuss below.. Finally,
the initial value of $z$ can be incorporated in terms of $x_{i},\,y_{i}$ using
the dimensionless Friedmann equation

$\displaystyle
z_{i}=\sqrt{\frac{1}{6}\,y_{i}^{2}+\frac{1}{3}\,v_{0}\,f(x_{i})\,\frac{1}{S^{2}}}\,.$
(4.20)

  2. 2.

We then proceed to solve the system of equations by taking adequately small
time steps $T$ in the appropriate range $T\in[T_{i}=0,\,T_{f}]$. We then plot
$N$ vs $T$ as given in figure 6 for Starobinsky potential. Typically, $N$
grows linearly with $T$ during near exponential inflation and a substantial
decrease in the rate of growth of $N$ indicates the end of inflation.

  3. 3.

In order to concretely determine the value of $N_{T}$, we plot $\epsilon_{H}$
vs $N$, and note the value of $N$ after which $\epsilon_{H}\geq 1$. By
definition, initially $N=0$. If $N_{T}<70$, then we repeat this step by
increasing the value of $x_{i}$, until we get $N_{T}\geq 70$. (Alternatively,
if inflation has not ended888Note that if one simulates the cosmological
equations in terms of number of e-folds $N$, rather than cosmic time $t$, this
step can usually be avoided by simulating the system from $N=0$ to $N=70$.
However, one has to adjust the value of $x_{i}$ in order to get enough
inflation., i.e $\epsilon_{H}<1$, at the end of our simulation, then either
one can increase the value of $T_{f}$ or decrease $x_{i}$. We suggest the
latter.) We can then define the number of e-folds before the end of inflation
to be $N_{e}=N_{T}-N$.

Figure 7: This figure describes the evolution of inflaton field $\phi$, and
its speed $\dot{\phi}$ in the left panel, while the Hubble parameter $H$ in
the right panel as a function of the number of e-folds before the end of
inflation $N_{e}$ for Starobinsky potential (3.48). Note that during slow-roll
inflation, $\dot{\phi}$ and $H$ are nearly constant, while $\phi$ changes
quite slowly. However, $\phi$ and $H$ begin to change rapidly towards the end
of inflation. After inflation ends, $\phi$ and $\dot{\phi}$ start oscillating
around the minimum of the potential (which is not shown in this figure).
Figure 8: Evolution of the slow-roll parameters $\epsilon_{H}$ and $\eta_{H}$
is shown as a function of the number of e-folds before the end of inflation
$N_{e}$ for Starobinsky potential (3.48). From this plot, it is easy to notice
that at early times when $N_{e}\gg 1$, the slow-roll conditions are satisfied
i.e $\epsilon_{H},\,|\eta_{H}|\ll 1$. However, the slow-roll conditions are
violated towards the end of inflation (marked by $N_{e}=0$ and
$\epsilon_{H}=1$).

  4. 4.

The pivot scale can then be fixed to an appropriate value, for example
$N_{e}=60$, as used in this work. Figure 7 describes the evolution of
$\phi,\,\dot{\phi},\,{\rm and}\,H$, while figure 8 illustrates the dynamics of
slow-roll parameters $\epsilon_{H},\,|\eta_{H}|$ for Starobinsky potential, as
determined from our code. As mentioned before, we usually plot the dynamics of
inflation as a function of $N_{e}$.

  5. 5.

In order to accurately fix the value of $v_{0}$, we need to impose
$A_{S}=2.1\times 10^{-9}$ at the pivot scale $N_{e}=60$. If $A_{S}$ is lower
than expected for the given value of $v_{0}$, then we increase the value of
$v_{0}$ or vice versa until we arrive at the correct value of $A_{S}$, and fix
the corresponding value of $v_{0}$.

Following the aforementioned algorithm, we can easily simulate the
inflationary background dynamics and investigate the evolution of relevant
quantities of our interest. Before going forward, we would like to stress that
many of the aforementioned steps (in the present version of our code) are
rather meant to be carried out manually by the user. While we are already
developing an automated version of this code (which will be presented in the
revised version of our paper), we believe that the present version will help
the user to understand the inflationary dynamics much better.

###  4.1 Phase-space analysis

Phase-space analysis of inflationary dynamics is usually carried out to
determine the set of initial conditions that results in adequate amount of
inflation, and hence it is important to access the generality of initial
conditions for inflation [69, 50, 20]. For a spatially flat background, the
phase-space portrait consists of trajectories of $\\{\phi,\,\dot{\phi}\\}$ for
different initial conditions, with fixed $H_{i}$. The standard algorithm to
generate such a plot is the following.

  1. 1.

The initial energy scale of inflation is kept constant by fixing the value of
initial Hubble parameter in the phase-space portrait simulations. A typical
value often used is $H_{i}\leq m_{p}$ (see [20] for detail). Hence, the user
is expected to incorporate an appropriate value of $z_{i}$.

  2. 2.

One can then input a suitable value of $x_{i}$ and determine the value of
$y_{i}$ for a given potential function $f(x)$ from the dimensionless Friedmann
equation (4.20) as

$y_{i}=\pm\sqrt{6}\,\sqrt{z_{i}^{2}-\frac{1}{3}\,v_{0}\,f(x_{i})\,\frac{1}{S^{2}}}~{}.$
(4.21)

  3. 3.

With these initial conditions, one can then simulate the system of
dimensionless differential equations for $\\{x,\,y,\,z,\,A\\}$ from $T_{i}=0$
till an appropriate $T_{f}$. One can then repeat the same step by
incorporating a number of different values of of $x_{i}$ in order to generate
the phase-space portrait for the given potential. We provide the GitHub link
to our phase-space portrait framework here
999https://github.com/bhattsiddharth/NumDynInflation/blob/main/inf_dyn_phase.py.
The phase-space portraits for Starobinsky potential (3.48), and quadratic
potential $V(\phi)\propto\phi^{2}$ are illustrated in the left and right
panels of figure 9 respectively.

Figure 9: The phase-space portrait $\\{\phi,\,\dot{\phi}\\}$ of the inflaton
field has been illustrated for Starobinsky potential (3.48) in the left panel,
and for quadratic potential $V(\phi)=\frac{1}{2}\,m^{2}\,\phi^{2}$ in the
right panel corresponding to different initial conditions
$\\{\phi_{i},\,\dot{\phi_{i}}\\}$ (plotted in solid black colour) with a fixed
initial scale $H_{i}$. The figure demonstrates that trajectories commencing
from a large class of initial field values (including those with large initial
velocities $\dot{\phi_{i}}$) quickly converge towards the slow-roll attractor
separatrix $\dot{\phi}=-V_{,\phi}/3H\simeq\mathrm{const.}$ (plotted in green
colour) as can be seen from the rapid decline in the inflaton speed until they
meet the green colour curve. After the end of inflation, the inflaton begins
to oscillate around the minimum of the potential.

In order to determine the degree of generality of inflation, we need to define
a measure for the distribution of $\\{\phi_{i},\,\dot{\phi_{i}}\\}$. The
correct choice for the measure might depend on the quantum theory of gravity.
However, a uniform measure is usually considered in the literature. Interested
readers are referred to [69, 20] for further detail.

###  4.2 Quantum fluctuations under the slow-roll approximation

In section 3, we described the expressions for a number of inflationary
observables such as $A_{{}_{S}},\,A_{{}_{T}},\,n_{{}_{S}},\,n_{{}_{T}},\,{\rm
and}\,r$ associated with the scalar and tensor power spectra which can be
determined purely from the dynamics of background quantities such as
$H,\,\epsilon_{H},\,\eta_{H}$ under the slow-roll approximation. Hence they
can be conveniently determined from our background dynamics code as discussed
earlier in section 4. For example, the scalar and tensor power spectra for
Starobinsky inflation have been plotted in figure 10 as a function $N_{e}$.
Similarly, one can plot the spectral indices101010In the standard literature,
one usually plots $r$ vs $n_{{}_{S}}$ for a given inflaton potential for a
range of possible values of $N_{*}\in[50,60]$ which can also be done easily
using our code. $n_{{}_{S}}-1$ and $n_{{}_{T}}$ and determine their values at
the pivot scale $N_{*}$. The spectral indices for Starobinsky potential have
been plotted in figure 11.

Figure 10: The power spectra of scalar and tensor quantum fluctuations
(computed using the slow-roll formulae (3.10) and (3.25) respectively) are
shown for comoving modes exiting the Hubble radius at different number of
e-folds $N_{e}$ before the end of inflation for Starobinsky potential (3.48).
The CMB window (shown in grey shaded region) corresponds to comoving modes in
the range $k_{\rm CMB}\in[0.0005,0.5]~{}{\rm Mpc}^{-1}$ that are being probed
by the current CMB missions. Fluctuations over larger scales (shown in red
shaded region) are outside the observable universe at present and those over
smaller scales remain to be (potentially) probed by a plethora of upcoming
missions, from GW observatories to PBHs.

Figure 11: The scalar and tensor spectral indices $n_{{}_{S}}$ and
$n_{{}_{T}}$ are shown as a function of $N_{e}$ for Starobinsky potential
(3.48) as determined by their slow-roll approximated formulae (3.23) and
(3.26) respectively. Around the pivot scale, they take the approximate values
$n_{{}_{S}}\simeq 0.967\text{ and }n_{{}_{T}}\simeq-0.0004$. Note that we have
plotted $n_{{}_{S}}-1$ (instead of $n_{{}_{S}}$) since it is the correct
scalar spectral index. The tensor-to-scalar ratio is given by
$r\simeq-8\,n_{{}_{T}}$.

##  5 Numerical analysis for quantum fluctuations during inflation

In the previous section we used slow-roll approximated formulae to study the
spectra of inflationary fluctuations in terms of background quantities such as
$H,\,\epsilon_{H},\,\eta_{H}$. Hence, we only had to simulate the background
dynamics for a given potential in order to plot the relevant inflationary
observables. However, if we want to analyze the behaviour of quantum
fluctuations more accurately, especially in situations where one or both the
slow-roll conditions (2.20) are violated, we need to numerically solve the
Mukhanov-Sasaki equation (3.4) corresponding to each comoving scale $k$.

For this purpose, we first rewrite the Mukhanov-Sasaki equation (3.4) in
cosmic time as

$\frac{{\rm d}^{2}v_{k}}{{\rm d}t^{2}}+H\frac{{\rm d}v_{k}}{{\rm
d}t}+\left[\frac{k^{2}}{a^{2}}-\frac{1}{a^{2}}\frac{z^{\prime\prime}}{z}\right]v_{k}=0\,.$
(5.1)

Note that here $z$ is not the dimensionless Hubble parameter used in our
numerical code, rather it is the variable $z=am_{p}\sqrt{2\epsilon_{H}}$ in
the Mukhanov-Sasaki equation (3.4). The effective mass term
$z^{\prime\prime}/z$ in (3.6) can be re-written as

$\frac{z^{\prime\prime}}{z}=a^{2}\left[\frac{5}{2}\frac{\dot{\phi}^{2}}{m_{p}^{2}}+2\frac{\dot{\phi}\ddot{\phi}}{Hm_{p}^{2}}+2H^{2}+\frac{1}{2}\frac{\dot{\phi}^{4}}{H^{2}m_{p}^{4}}-V_{,\phi\phi}(\phi)\right]\,.$
(5.2)

Since $v_{k}$ is a complex valued function, it is convenient to split it into
its real and imaginary parts to study their evolution separately for the
numerical analysis. While both will follow the same evolution equation, they
will be supplied with different initial conditions in the form of the real and
imaginary parts of the Bunch-Davies vacuum (3.8). Writing the Mukhanov-Sasaki
equation for scalar fluctuations in terms of dimensionless variables, we
obtain

$\boxed{\frac{{\rm d}^{2}v_{k}}{{\rm d}T^{2}}+z\,\frac{{\rm d}v_{k}}{{\rm
d}T}+\left[\frac{k^{2}}{A^{2}}-\frac{5}{2}\,y^{2}+2\,\frac{y}{z}\left(3\,z\,y+\frac{v_{0}}{S^{2}}\,f_{,x}\right)-2\,z^{2}-\frac{1}{2}\,\frac{y^{4}}{z^{2}}+\frac{v_{0}}{S^{2}}\,f_{,xx}\right]v_{k}=0}\,.$
(5.3)

Figure 12: Evolution of scalar power $\frac{k^{3}}{2\pi^{2}}|\zeta_{k}|^{2}$
is plotted by numerically solving the Mukhanov-Sasaki equation (5.3) for a
mode exiting the Hubble radius at about 60 e-folds before the end of inflation
(for Starobinsky potential). At early times when the mode is sub-Hubble, i.e.
$k\gg aH$, the power decreases as ${\cal P}_{\zeta}\sim(aH)^{-2}$ as expected.
After the Hubble-exit, the power freezes to a constant in the super-Hubble
regime when $k\ll aH$. We note down its value after the mode-freezing as the
super-Hubble scale power corresponding to that mode. Repeating the procedure
for a range of scales $k$ yields us the power spectrum of scalar fluctuations.
The same numerical analysis can be carried out for tensor fluctuations.

Our primary goal in this section is to numerically solve equation (5.3) for
the Fourier modes $v_{k}$ corresponding to each comoving scale $k$ and plot
the frozen value of the scalar power spectrum of $\zeta_{k}$ given by (3.9)
after the mode becomes super-Hubble. We can conveniently relate a comoving
scale $k$ to its Hubble-exit epoch by $k=aH$. Since we are only interested in
the super-Hubble power spectra, we only need to simulate the system to evolve
$v_{k}$ for a small duration of time around the Hubble-exit of scale $k$
(which should be sufficiently early enough to impose Bunch-Davies initial
conditions and sufficiently late enough for the mode to be frozen outside the
Hubble radius).

In the following, we discuss the algorithm to solve the Mukhanov-Sasaki
equation (5.3) and determine the scalar power spectrum (3.9) numerically. We
also discuss how to solve the corresponding equation (3.18) for the tensor
power spectrum (at linear order in perturbation theory)111111Second-order
tensor fluctuations which are induced by first-order scalar fluctuations will
be discussed in the revised version of our manuscript.. The dimensionless
Mukhanov-Sasaki equation for tensor fluctuations is given by

$\boxed{\frac{{\rm d}^{2}h_{k}}{{\rm d}T^{2}}+z\,\frac{{\rm d}h_{k}}{{\rm
d}T}+\left[\frac{k^{2}}{A^{2}}+\frac{1}{2}\,y^{2}-2\,z^{2}\right]h_{k}=0}\,.$
(5.4)

We explicitly write down the Mukhanov-Sasaki equations for scalar and tensor
fluctuations in terms of dimensionless variables (as used in our code) in the
following way

$\displaystyle v_{k,T}$ $\displaystyle=\frac{{\rm d}v_{k}}{{\rm d}T}\,,$ (5.5)
$\displaystyle\frac{{\rm d}v_{k,T}}{{\rm d}T}$
$\displaystyle=-z\,v_{k,T}-\left[\frac{k^{2}}{A^{2}}-\frac{5}{2}\,y^{2}+2\,\frac{y}{z}\left(3\,z\,y+\frac{v_{0}}{S^{2}}\,f_{,x}\right)-2\,z^{2}-\frac{1}{2}\,\frac{y^{4}}{z^{2}}+\frac{v_{0}}{S^{2}}\,f_{,xx}\right]v_{k}\,;$
(5.6) $\displaystyle h_{k,T}$ $\displaystyle=\frac{{\rm d}h_{k}}{{\rm d}T}\,,$
(5.7) $\displaystyle\frac{{\rm d}h_{k,T}}{{\rm d}T}$
$\displaystyle=-z\,h_{k,T}-\left(\frac{k^{2}}{A^{2}}+\frac{1}{2}\,y^{2}-2\,z^{2}\right)\,h_{k}\,.$
(5.8)

In our numerical set up, we split $v_{k}$ and $h_{k}$ into their real and
imaginary parts and simulate them separately with appropriate Bunch-Davies
initial conditions. We begin with a discussion of numerical simulations of the
Mukhanov-Sasaki equations (5.3) and (5.4) for a purely slow-roll potential
(which we choose to be the Starobinsky potential (3.48) as usual), before
moving forward to discuss the same for a potential with a slow-roll violating
feature. This latter case is of the primary focus of our paper. In particular,
we will illustrate our numerical scheme for the case of a base slow-roll
potential possessing a tiny local bump feature, which was proposed in [44] in
the context of PBH formation.

Figure 13: The super-Hubble power spectra of scalar fluctuations (in green
colour) and tensor fluctuations (in red colour) are plotted for modes exiting
the Hubble radius at different number of e-folds $N_{e}$ before the end of
inflation for Starobinsky potential (3.48). The solid curves represent power
spectra computed under the slow-roll approximation (3.10), while the dotted
curves represent the power computed by numerically solving the Mukhanov-Sasaki
equation (5.3). We conclude that for Starobinsky model, since slow-roll
conditions $\epsilon_{H},\,|\eta_{H}|\ll 1$ are easily satisfied for most part
of inflation, the power spectra computed under the slow-roll approximation
match quite well with their numerically determined counterparts.

###  5.1 Numerical analysis for slow-roll potentials

  1. 1.

As the first step, we numerically solve the background dynamics for a given
potential, determine the values of all relevant parameters of the potential
and the evolution of relevant primary dynamical variables
$\\{x,\,y,\,z,\,A\\}$ as well as the derived quantities such as
$\\{N_{e},\,\epsilon_{H},\,\eta_{H}\\}$ (as discussed in section 4).

  2. 2.

We then proceed to identify different comoving scales $k$. This can be done by
determining their Hubble-exit epochs in the following way. For example, we
plot $aH$ (in log scale) against $N_{e}$ and identify the value of $aH$ at
$N_{e}=N_{*}$ to be the CMB pivot scale $k_{p}$. As mentioned before, we take
$N_{*}=60$ in all our analysis. Similarly, we associate a corresponding value
of $N_{e}$ to each comoving scale $k$ by the value of $aH$ at its Hubble-exit
epoch. This step ensures that we have a one-to-one correspondence between $k$
and $N_{e}$ in our analysis and we can use them interchangeably.

  3. 3.

We intend to impose Bunch-Davies initial conditions for a given mode $v_{k}$
at an epoch when it is sub-Hubble. As it tuns out, for most potentials, the
Bunch-Davies initial conditions can be safely imposed as long as $k\geq
100\,aH$. Hence, rather than simulating the Mukhanov-Sasaki equation for each
mode (making Hubble-exit at the corresponding value of $N_{e}$) all through
the inflationary history (starting from $\phi_{i}>\phi_{*}$), we actually
impose the initial conditions from the background solutions for
$\\{x,\,y,\,z,\,A\\}$ at around $5$ e-folds before the Hubble-exit of that
mode. This step greatly  reduces the running-time of the code. We then
incorporate the initial value of scale factor $A_{i}$ at the same epoch,
namely $A_{i}\exp{(N_{T}-N_{e}-5)}$ and do the same for the initial values of
the field $x_{i}$, and its derivative $y_{i}$. The initial conditions for the
mode functions $v_{k}$ and their derivatives $\dot{v_{k}}$ can then be safely
taken to be of Bunch-Davies type.

  4. 4.

We solve the set of cosmological equations with these initial conditions for a
period of time $T=T_{i}\to T=T_{f}$ such that the mode becomes super-Hubble
and its power ($k^{3}|\zeta_{k}|^{2}/2\pi^{2}$) is frozen to a constant value
(see figure 12), which is typically within $5$ e-folds after Hubble-exit in
the kind of models we are interested in. We note down this frozen value as the
value of the power spectrum of that mode. While we have been discussing about
scalar fluctuations mostly, the same can be done for tensor fluctuations which
we have incorporated in our code.

  5. 5.

We then select another mode that leaves the Hubble radius at some epoch
$N_{e}$ and repeat the procedure until we have collected the frozen super-
Hubble power spectra of a range of scales that we are interested in (see
figure 13).

Figure 14: This is a schematic plot of an asymptotically-flat inflationary
potential with a tiny bump/dip feature (5.9) near intermediate field values
$\phi\simeq\phi_{\rm PBH}$ that leads to an enhancement of the scalar power
spectrum. The full potential asymptotes to the base (slow-roll) potential near
the CMB window $\phi\simeq\phi_{*}$, thus satisfying observational constraints
on large cosmological scales. Slow-roll is violated around the feature, whose
position $\phi_{\rm PBH}$ dictates the range of moving scales $k$ that receive
amplification of power (which accordingly determines the mass and abundance of
formed PBHs). Note that the feature has been greatly exaggerated for
illustration purpose. In most realistic models, both the height and the width
of the feature are too small to be seen (without zooming-in considerably).

From the above numerical analysis of featureless vanilla potentials which
exhibit slow-roll dynamics until close to the end of inflation, we observe
that the power spectra of scalar and tensor fluctuations are nearly scale-
invariant (with small red-tilt) and their behaviour (as obtained from
numerically solving the Mukhanov-Sasaki equation) matches quite well with the
analytical predictions under the slow-roll approximations (see figure 13).

However, for potentials exhibiting a small-scale feature at intermediate field
values $\phi<\phi_{*}$, there might exist a short period of slow-roll
violating phase before the end of inflation during which slow-roll
approximations break down. In particular, as we will see, while the first
slow-roll parameter remains small $\epsilon_{H}\ll 1$, the second slow-roll
parameter might become $\eta_{H}\sim{\cal O}(1)$. Hence, a numerical analysis
of the Mukhanov-Sasaki equation is desired in order to determine the scalar
power spectrum more accurately. This will be the main focus of discussion in
the next subsection.

Figure 15: Evolution of the field value $\phi$ for the KKLT potential with a
tiny bump (5.10) is shown in solid green curve as a function of number of
e-folds $N_{e}$ before the end of inflation. We note that at CMB scales,
$\phi$ is much smaller than the corresponding value in the base model (KKLT)
(shown in dashed black curve). At intermediate scales when the inflaton
evolves across the bump feature in the potential, we gain a lot of extra
e-folds of expansion $\Delta N_{e}\simeq 15$ with little change in the field
value. After crossing the feature, evolution of $\phi$ mimics its
corresponding value in the base model. Figure 16: Evolution of the slow-roll
parameters $\epsilon_{H}$ and $\eta_{H}$ is shown in solid green and solid red
curves respectively for the KKLT potential with a tiny bump (5.10). Both
$\epsilon_{H}$ and $\eta_{H}$ are close to their corresponding values for the
base KKLT potential (shown in dashed curves) at early times near the CMB
window. At $N_{e}\simeq 30$, the value of $\epsilon_{H}$ starts decreasing
rapidly leading to an increase in $\eta_{H}$ from $|\eta_{H}|\ll 1$ to a
higher and positive value $\eta_{H}\simeq+3.3$ (almost USR phase). Thereafter,
the inflaton enters a phase of constant-roll inflation (where
$\eta_{H}\simeq-0.37$) before returning to the final slow-roll phase.

###  5.2 Numerical analysis for potentials with a local bump/dip feature

In order to facilitate PBH formation, we need a large amplification of scalar
power spectrum at smaller scales during inflation. This can be achieved by
introducing a small-scale feature in the potential which leads to a transient
period of slow-roll violating phase (including a short almost-USR phase).
Adequate amplification in the super-Hubble scalar power spectrum results in a
large density contrast in the post-inflationary universe (upon the Hubble-
entry of the corresponding modes) which in turn can collapse to form PBHs.
Usually, such a PBH feature in the potential leads to an increase in the value
of the second slow-roll parameter $\eta_{H}$ from near-zero to a positive
value of $\eta_{H}\sim{\cal O}(1)$.

A number of models with different types of features have been proposed in the
recent literature (as mentioned before) that facilitate the amplification of
scalar power spectrum at small scales. The most common amongst them is an
inflection point-like feature. However, we choose the model proposed in [44]
in which the base inflaton potential $V_{b}(\phi)$ possesses a tiny local bump
or dip $\pm\varepsilon(\phi,\phi_{0})$ at an intermediate field value
$\phi_{0}$ of the form

$V(\phi)=V_{b}(\phi)\,\left[1\pm\varepsilon(\phi,\phi_{0})\right]\,,$ (5.9)

where we assume $V_{b}(\phi)$ to be a symmetric or an anti-symmetric
asymptotically-flat potential in order to satisfy CMB constraints at large
cosmological scales. Such a potential has been schematically illustrated in
figure 14. To be specific, in this paper we choose the base potential to be
the D-brane KKLT potential [70, 71, 72, 73] with a tiny Gaussian bump of the
form [44]

$\boxed{V(\phi)=V_{0}\,\frac{\phi^{2}}{m^{2}+\phi^{2}}\,\left[1+A\,\exp{\left({-\frac{1}{2}\,\frac{(\phi-\phi_{0})^{2}}{\sigma^{2}}}\right)}\right]}~{},$
(5.10)

where $m$ is a mass scale in the KKLT model, while $A$ and $\sigma$ represent
the height and the width of the tiny bump respectively. We use this particular
model to demonstrate our numerical framework because of its simplicity and
efficiency. However, one can choose any model of their interest. Values of all
the parameters appearing in (5.10), which we use in our numerical analysis,
have been explicitly shown in figure 17.

The height and the width of the (bump/dip) feature required to facilitate
adequate amount of power amplification are quite small, and hence the feature
is tiny and local (in contrast to inflection point-like features). This
ensures that the feature does not significantly affect the CMB observables.
However, since we gain a lot of extra e-folds of expansion $\Delta N_{e}\simeq
15$ with little change in the field value (as shown in figure 15) when the
inflaton crosses the feature, the CMB pivot scale gets shifted towards smaller
values as compared to the same for the base potential.

Figure 17: The super-Hubble power spectra of scalar fluctuations are plotted
for modes exiting the Hubble radius at different number of e-folds $N_{e}$
before the end of inflation for KKLT potential with a tiny bump (5.10). The
solid green curve represents the power spectrum computed under the slow-roll
approximation (3.10), while the dotted red curve represents the power computed
by numerically solving the Mukhanov-Sasaki equation (5.3). Since the second
slow-roll condition is violated due to the presence of the bump (leading to
$\eta_{H}>1$), the slow-roll approximation underestimates the value of power
spectrum near the peak (as well as the position of the peak) for modes exiting
the Hubble radius near the epoch when the inflaton crosses the local maximum
around the bump feature. This leads to an incorrect estimation of the mass
fraction as well as the central mass of the PBHs formed in this scenario.

Since slow-roll is violated in these models (as shown in figure 16), we need
to solve the Mukhanov-Sasaki equation numerically in order to accurately
compute the scalar power spectrum. This has been explicitly demonstrated in
figure 17. Note that the inflationary dynamics in such models contains a
number of phases that include an early slow-roll phase SR-I near the CMB
window, a transition T-I from the early SR-I to the subsequent almost ultra
slow-roll phase USR and a transition T-II back to the next slow-roll phase SR-
II after passing through an intermediate constant-roll phase CR (shown in
figure 16).

Figure 18: Evolution of the effective mass term in the Mukhanov-Sasaki
equation is shown here in solid green curve for the KKLT potential with a tiny
bump (5.10). Important transient phases of the scalar field dynamics are also
highlighted following the behaviour of the second slow-roll parameter
$\eta_{H}$ (plotted in dashed blue curve). There is a sharp dip in the
effective mass term when the field transitions from the first slow-roll phase
(SR-I) to an almost ultra slow-roll phase (USR). The inflaton later makes a
transition to a phase of constant-roll inflation with $\eta_{H}\simeq-0.37$,
before reaching a final slow-roll phase until the end of inflation.

The effective mass term $z^{\prime\prime}/z$ in the Mukhanov-Sasaki equation
(3.4), which primarily governs the dynamics of scalar fluctuations, has been
shown in figure 18, and the resultant Hubble-exit behaviour of different modes
is described in figure 19. As the inflaton approaches the PBH feature,
$\eta_{H}$ starts to increase rapidly. This is accompanied by an initial sharp
dip in the effective mass term $z^{\prime\prime}/z$, which then increases to a
higher plateau as $\eta_{H}$ approaches its maximum in the USR type phase. The
modes that leave the Hubble radius slightly before the transition already
start receiving power amplification on super-Hubble scales as shown by the
orange color curve in figure 19.

It is worth noting that the power spectrum exhibits a dip which corresponds to
very narrow range of scales $k\simeq k_{\rm dip}$ that leave the Hubble radius
a few e-folds before the USR phase (shown by the blue color curve in figure
19). The maximum rate of growth observed in this model is consistent with the
steepest growth bound discussed in [65]. Near the USR phase when
$\eta_{H}\gtrsim 3$, $z^{\prime\prime}/z$ saturates to a constant value. Modes
leaving the Hubble radius around this USR epoch receive maximal amplification
in their super-Hubble power spectrum.

Figure 19: This figure demonstrates the horizon exit behaviour of different
modes i.e the evolution of $\sqrt{\cal{P}_{\zeta}}$ for different modes as
they cross the Hubble radius in KKLT model with a Gaussian bump (5.10). The
dot on each curve corresponds to its Hubble-exit epoch. The sharp dip in the
power spectrum (figure 17) corresponds to the mode $k_{\mathrm{dip}}$ (plotted
in blue color) that exits the Hubble radius a few e-folds before the
commencement of the USR phase. The mode $k_{\mathrm{PBH}}$ which exits the
Hubble radius during the USR phase receives a maximal amplification of power
(plotted in green color).

As the field crosses the maximum of the bump feature, $\eta_{H}$ decreases to
a constant negative value. It stays in this constant-roll phase until the
inflaton meets the base potential eventually and approaches the final slow-
roll phase in its dynamics before the end of inflation. The dynamics of scalar
fluctuations in the aforementioned phases are quite rich, and interesting.
However, since the main aim of this paper is to illustrate how to use our
numerical code with an example, we do not discuss these phases and their
impact on the power spectrum (some of which have been explicitly shown in
figures 18, 19, 20), and refer the interested readers to [63] for more detail.

Figure 20: The super-Hubble power spectrum of scalar fluctuations obtained by
numerically solving the Mukhanov-Sasaki equation (5.3) is plotted here for the
KKLT potential with a tiny bump (5.10). At large scales (around the CMB
window), the power spectrum matches that of the base KKLT potential. The power
spectrum (after exhibiting a sharp dip) receives a large amplification at
intermediate scales $k\sim k_{\rm PBH}$ around the USR phase. The maximum rate
of growth observed in this model is consistent with the steepest growth bound
discussed in [65]. After reaching the peak, the power then decreases at a
steady rate during the constant-roll phase before finally asymptoting towards
its base slow-roll value towards the end of inflation.

##  6 Future extension of our numerical framework

In preceding sections, we described the relevant cosmological equations
governing the inflationary dynamics in terms of dimensionless variables. We
also introduced our numerical code (written in terms of cosmic time) that can
easily simulate the inflationary dynamics both at the background level in
section 4 and at linear order in perturbation theory in section 5. However,
with minimal to moderate extension, our numerical code can be used to simulate
a number of different scenarios associated with scalar field dynamics both
during inflation as well as in the post-inflationary universe. We have already
started working on some of these aspects which will appear in the revised
version of our manuscript. In the following we discuss some of the important
future extensions of our code that we intend to include in our revised
version.

  1. 1.

As stressed in section 4, the present version of our numerical code, although
quite fast and neat, contains segments that require the user to carry out a
number of tasks manually. While we believe that the present version will
definitely help a user (who is relatively new to the field) to understand the
inflationary dynamics much better, we are already developing an automated
version of this code that is much more compact and requires substantially less
manual involvement of the user. We also plan to make the code even faster. We
will present the updated version of our code in the revised version of our
paper and refer to it in the same GitHub link
121212https://github.com/bhattsiddharth/NumDynInflation.

  2. 2.

In section 5, we briefly discussed how to solve the evolution equation for the
tensor fluctuations at linear order in perturbation theory. However, first-
order scalar fluctuations induce tensor fluctuations at second order in
perturbation theory which might be significant when slow-roll is violated,
especially in the scenario where scalar power spectrum is largely amplified in
order to source PBH formation. We will incorporate the computation of such
scalar-induced Gravitational Waves in the updated version of our code.

  3. 3.

In our analysis, we used the Mukhanov-Sasaki variable $v_{k}$ with its
corresponding evolution equation (5.1) in order to simulate scalar
fluctuations at linear order. Our code runs quite quickly and generates the
power spectrum without requiring a large computational time. However, one can
study the scalar power spectrum by using other variables proposed in the
literature. We are particularly interested in two such variables. The first
one is the curvature perturbation $\zeta_{k}$ itself using the corresponding
evolution equation (3.11). For example, it was claimed in [74] that numerical
simulations are more stable in terms of $\zeta_{k}$ since it is explicitly
frozen well outside the Hubble radius. Similarly, authors of [75] have
suggested a change in the variable $v_{k}$ in the form
$g_{k}\equiv\frac{v_{k}}{z}\,e^{ik\tau}$ that is supposed to make the
simulations much more stable since it removes the early time
oscillations131313We are thankful to Christian Byrnes for bringing this to our
attention..

In the revised version of our work, we plan on carrying out a comprehensive
numerical analysis to compare both stability as well speed of the numerical
simulations using all three of the aforementioned variables $v_{k}$,
$\zeta_{k}$, and $g_{k}$.

  4. 4.

It is easy to extend our numerical analysis to incorporate the dynamics of
more than one scalar fields during inflation, at least in the background
level. In the future, we are going to provide an updated numerical code to
study both the background dynamics as well as quantum fluctuations in two-
field inflationary dynamics (where the second field might also source
inflation, or act as a spectator field).

  5. 5.

Additionally, our code can be extended to simulate the post-inflationary
dynamics of the inflaton field as well as to study parametric resonance by
simulating the evolution of different Fourier modes of the inflaton
fluctuations. During the post-inflationary oscillations, it is usually
advisable to make a change in the Mukhanov-Sasaki variable of the form
${\tilde{v}_{k}}=a^{1/2}\,v_{k}$ as suggested in [76, 77]. The code can also
be extended to study the dynamics of scalar field dark matter and quintessence
by suitably redefining the dimensionless variables as per the energy scale of
the dynamics.

##  7 Discussion

In the present version of our manuscript, we introduced our numerical approach
to simulate the cosmological equations in order to study the inflationary
dynamics at the level of both background as well as linear-order in
perturbation theory. We provided the link to our open-source GitHub repository
where we have supplied a Python-based simple numerical code to simulate
inflationary dynamics in terms of cosmic time $t$. We explicitly demonstrated
how to use the code to study the inflationary background dynamics in section 4
that includes plotting the phase-space portrait of inflation as well as to
characterise quantum fluctuations during inflation using the simulations of
the background dynamics.

Section 5 was dedicated to study quantum fluctuations during inflation
(without using slow-roll approximated expressions) by numerically solving the
mode function equations of scalar and tensor fluctuations. For a featureless
slow-roll inflaton potential, the difference between the results obtained
numerically and those obtained under slow-roll approximations were negligible
until close to the end of inflation, as expected. We used the Starobinsky
potential (3.48) as an example to illustrate our analysis. Our primary focus
was the numerical evaluation of the scalar power spectrum ${\cal
P}_{\zeta}(k)$ for potentials that exhibit a slow-roll violating feature. In
particular, we used the example of an asymptotically-flat base inflationary
potential that possesses a tiny local bump feature (5.9) to illustrate our
numerical scheme in section 5.2. By suitably choosing the parameters of the
potential (5.10), one can achieve a large enough amplification of the scalar
power spectrum in order to facilitate the formation of PBHs in the post-
inflationary epoch.

In our numerical analysis for the case of potentials with a slow-roll
violating feature, we explicitly chose the parameters in order to amplify the
small-scale scalar power by a factor of $\sim 10^{7}$ with respect to the
corresponding power at large cosmological scales, as is usually assumed in the
literature to facilitate the formation of of PBHs in the subsequent radiation
dominated epoch. However, it is important to stress that a significant growth
in the scalar power spectrum at small-scales might engender the dynamics to
enter into non-perturbative regime. For example, a careful computation of loop
corrections to the two-point scalar fluctuations demonstrates [78, 79] that
contribution from $1$-loop effects becomes of the same order as the tree-level
computation (which we carried out in this work) if ${\cal P}_{\zeta}(k)\sim
10^{-2}$ indicating a breakdown of the perturbative analysis.

Moreover, the mechanism of PBH formation in the context of single field models
of inflation involves additional intricacies that demand for a non-
perturbative analysis of primordial fluctuations. Firstly, a sharp drop in the
classical drift speed of the inflaton due to the presence of the PBH-producing
feature often incites the system to enter into a phase where stochastic
quantum diffusion effects become non-negligible, at times even significant.
More importantly, since PBHs form from rare extreme peaks, and hence are
determined by the tail of the probability distribution function (PDF)
$P[\zeta]$ of the primordial fluctuations, perturbative computations based on
only power-spectrum lead to an inaccurate estimation of the PBH mass fraction.
Consequently, determination of the full primordial PDF becomes crucial, which
can be computed non-perturbatively using the stochastic inflation framework
[80, 81, 82, 83, 84, 85, 86, 88, 87, 89, 90, 91, 92, 94, 95, 93, 96, 97, 98,
99, 100, 101, 102, 103, 104, 105], which usually predicts a non-Gaussian
exponential tail [86, 90]. Tail of the primordial PDF can also be computed by
using semi-classical techniques discussed in [106]. Determining the tail of
the primordial PDF is an important and active topic of research [90, 106, 107,
108, 109, 110, 111] at present, which is beyond the scope of our perturbative
analysis presented here.

Before concluding, let us mention that we also stressed upon various important
future extensions of our numerical scheme in section 6 that will result in
making our code more efficient, and enable us to simulate the scalar field
dynamics in a number of interesting scenarios. This includes updating our code
to make it more compact and automated, as well as extending it to study the
spectrum of scalar-induced gravitational waves, inflaton dynamics in the post-
inflationary epoch, multi-field inflationary dynamics, and even scalar field
models of dark matter and dark energy. We will incorporate most of these
additional features in the revised version of our manuscript. In the meantime,
we welcome constructive comments and suggestions as well as queries from
interested readers which will help us in improving the quality of our
numerical work and presentation in the updated version of the paper.

##  8 Acknowledgements

S.S.M thanks Satadru Bag, Shabbir Shaikh, and Varun Sahni for crucial inputs
during the early stages of development of this code. The authors are grateful
to Parth Bhargava and Sanket Dave for stimulating discussions on various
topics related to numerical dynamics discussed in this paper. S.S.M. is
supported as a postdoctoral Research Associate at the School of Physics and
Astronomy, University of Nottingham by the STFC funded consolidated grant, UK.
S.S.B was supported by the INSPIRE scholarship of the Department of Science
and Technology (DST), Govt. of India during his Master’s thesis work during
which a significant portion of this work was carried out.

## References

  * [1] A. A. Starobinsky, “A New Type of Isotropic Cosmological Models Without Singularity,” Phys. Lett. B 91, 99-102 (1980).
  * [2] A. H. Guth, “The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems,” Phys. Rev. D 23, 347 (1981).
  * [3] A. D. Linde, “A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems”, Phys. Lett. B 108, 389 (1982).
  * [4] A. Albrecht and P. J. Steinhardt, “Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking,” Phys. Rev. Lett. 48, 1220 (1982).
  * [5] A. D. Linde, “Chaotic Inflation,” Phys. Lett. B 129, 177-181 (1983).
  * [6] A.D. Linde, Particle Physics and Inflationary Cosmology, Harwood, Chur, Switzerland (1990).
  * [7] D. Baumann, “TASI Lectures on Inflation”, [arXiv:0907.5424].
  * [8] V. F. Mukhanov and G. V. Chibisov, JETP Lett. 33, 532 (1981).
  * [9] S. W. Hawking, Phys. Lett. B 115, 295 (1982).
  * [10] A. A. Starobinsky, Phys. Lett. B 117, 175 (1982).
  * [11] A. H. Guth and S. -Y. Pi, Phys. Rev. Lett. 49, 1110 (1982).
  * [12] A.A. Starobinsky, JETP Lett., 30 , 682 (1979).
  * [13] V. Sahni, Phys. Rev. D42, 453 (1990).
  * [14] M. Tegmark, “What does inflation really predict?,” JCAP 04 (2005), 001 [arXiv:astro-ph/0410281 [astro-ph]].
  * [15] Y. Akrami et al. [Planck], “Planck 2018 results. X. Constraints on inflation” Astron. Astrophys. 641, A10 (2020) [arXiv:1807.06211 [astro-ph.CO]].
  * [16] S. S. Mishra, V. Sahni and A. A. Starobinsky, “Curing inflationary degeneracies using reheating predictions and relic gravitational waves,” JCAP 05, 075 (2021) [arXiv:2101.00271 [gr-qc]].
  * [17] P. A. R. Ade et al. [BICEP and Keck], “Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season,” Phys. Rev. Lett. 127, no.15, 151301 (2021) [arXiv:2110.00483 [astro-ph.CO]].
  * [18] S. S. Mishra and V. Sahni, “Canonical and Non-canonical Inflation in the light of the recent BICEP/Keck results,” [arXiv:2202.03467 [astro-ph.CO]].
  * [19] R. Kallosh and A. Linde, “BICEP/Keck and cosmological attractors,” JCAP 12, no.12, 008 (2021) [arXiv:2110.10902 [astro-ph.CO]].
  * [20] S. S. Mishra, V. Sahni and A. V. Toporensky, “Initial conditions for Inflation in an FRW Universe,” Phys. Rev. D 98, no.8, 083538 (2018) [arXiv:1801.04948 [gr-qc]].
  * [21] S. Hawking, “Gravitationally collapsed objects of very low mass,” Mon. Not. Roy. Astron. Soc. 152 (1971), 75
  * [22] B. J. Carr and S. W. Hawking, “Black holes in the early Universe,” Mon. Not. Roy. Astron. Soc. 168 (1974), 399-415
  * [23] B. J. Carr, “The Primordial black hole mass spectrum,” Astrophys. J. 201 (1975), 1-19
  * [24] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, “Primordial black holes—perspectives in gravitational wave astronomy,” Class. Quant. Grav. 35 (2018) no.6, 063001 [arXiv:1801.05235 [astro-ph.CO]].
  * [25] G. F. Chapline, “Cosmological effects of primordial black holes,” Nature 253 (1975) no.5489, 251-252
  * [26] P. Meszaros, “Primeval black holes and galaxy formation,” Astron. Astrophys. 38 (1975), 5-13
  * [27] P. Ivanov, P. Naselsky and I. Novikov, “Inflation and primordial black holes as dark matter,” Phys. Rev. D 50 (1994), 7173-7178
  * [28] B. Carr, F. Kuhnel and M. Sandstad, “Primordial Black Holes as Dark Matter,” Phys. Rev. D 94 (2016) no.8, 083504 [arXiv:1607.06077 [astro-ph.CO]].
  * [29] A. M. Green and B. J. Kavanagh, “Primordial Black Holes as a dark matter candidate,” J. Phys. G 48, no.4, 043001 (2021) [arXiv:2007.10722 [astro-ph.CO]].
  * [30] B. Carr and F. Kuhnel, “Primordial Black Holes as Dark Matter: Recent Developments,” Ann. Rev. Nucl. Part. Sci. 70 (2020), 355-394 [arXiv:2006.02838 [astro-ph.CO]].
  * [31] E. Bugaev and P. Klimai, “Large curvature perturbations near horizon crossing in single-field inflation models,” Phys. Rev. D 78 (2008), 063515 [arXiv:0806.4541 [astro-ph]].
  * [32] C. Germani and T. Prokopec, “On primordial black holes from an inflection point,” Phys. Dark Univ. 18 (2017), 6-10 [arXiv:1706.04226 [astro-ph.CO]].
  * [33] J. M. Ezquiaga, J. Garcia-Bellido and E. Ruiz Morales, “Primordial Black Hole production in Critical Higgs Inflation,” Phys. Lett. B 776 (2018), 345-349 doi:10.1016/j.physletb.2017.11.039 [arXiv:1705.04861 [astro-ph.CO]].
  * [34] J. Garcia-Bellido and E. Ruiz Morales, “Primordial black holes from single field models of inflation,” Phys. Dark Univ. 18 (2017), 47-54 [arXiv:1702.03901 [astro-ph.CO]].
  * [35] I. Dalianis, A. Kehagias and G. Tringas, “Primordial black holes from $\alpha$-attractors,” JCAP 01 (2019), 037 [arXiv:1805.09483 [astro-ph.CO]].
  * [36] O. Özsoy, S. Parameswaran, G. Tasinato and I. Zavala, “Mechanisms for Primordial Black Hole Production in String Theory,” JCAP 07 (2018), 005 [arXiv:1803.07626 [hep-th]].
  * [37] M. Cicoli, V. A. Diaz and F. G. Pedro, “Primordial Black Holes from String Inflation,” JCAP 06 (2018), 034 [arXiv:1803.02837 [hep-th]].
  * [38] M. P. Hertzberg and M. Yamada, “Primordial Black Holes from Polynomial Potentials in Single Field Inflation,” Phys. Rev. D 97 (2018) no.8, 083509 [arXiv:1712.09750 [astro-ph.CO]].
  * [39] G. Ballesteros and M. Taoso, “Primordial black hole dark matter from single field inflation,” Phys. Rev. D 97 (2018) no.2, 023501 [arXiv:1709.05565 [hep-ph]].
  * [40] H. Di and Y. Gong, “Primordial black holes and second order gravitational waves from ultra-slow-roll inflation,” JCAP 07 (2018), 007 [arXiv:1707.09578 [astro-ph.CO]].
  * [41] H. Motohashi, S. Mukohyama and M. Oliosi, “Constant Roll and Primordial Black Holes,” JCAP 03 (2020), 002 [arXiv:1910.13235 [gr-qc]].
  * [42] R. Mahbub, “Primordial black hole formation in inflationary $\alpha$-attractor models,” Phys. Rev. D 101 (2020) no.2, 023533 [arXiv:1910.10602 [astro-ph.CO]].
  * [43] N. Bhaumik and R. K. Jain, “Primordial black holes dark matter from inflection point models of inflation and the effects of reheating,” JCAP 01 (2020), 037 [arXiv:1907.04125 [astro-ph.CO]].
  * [44] S. S. Mishra and V. Sahni, “Primordial Black Holes from a tiny bump/dip in the Inflaton potential,” JCAP 04, 007 (2020) [arXiv:1911.00057 [gr-qc]].
  * [45] K. Kefala, G. P. Kodaxis, I. D. Stamou and N. Tetradis, “Features of the inflaton potential and the power spectrum of cosmological perturbations,” Phys. Rev. D 104 (2021) no.2, 023506 [arXiv:2010.12483 [astro-ph.CO]].
  * [46] H. V. Ragavendra, P. Saha, L. Sriramkumar and J. Silk, “Primordial black holes and secondary gravitational waves from ultraslow roll and punctuated inflation,” Phys. Rev. D 103 (2021) no.8, 083510 [arXiv:2008.12202 [astro-ph.CO]].
  * [47] O. Özsoy and Z. Lalak, “Primordial black holes as dark matter and gravitational waves from bumpy axion inflation,” JCAP 01 (2021), 040 [arXiv:2008.07549 [astro-ph.CO]].
  * [48] R. Zheng, J. Shi and T. Qiu, “On Primordial Black Holes and secondary gravitational waves generated from inflation with solo/multi-bumpy potential,” [arXiv:2106.04303 [astro-ph.CO]].
  * [49] K. Inomata, E. McDonough and W. Hu, “Primordial black holes arise when the inflaton falls,” Phys. Rev. D 104 (2021) no.12, 123553 [arXiv:2104.03972 [astro-ph.CO]].
  * [50] R. Brandenberger, “Initial conditions for inflation — A short review,” Int. J. Mod. Phys. D 26 (2016) no.01, 1740002 [arXiv:1601.01918 [hep-th]].
  * [51] D. Baumann, PoS TASI 2017, 009 (2018) [arXiv:1807.03098 [hep-th]].
  * [52] J. M. Maldacena, “Non-Gaussian features of primordial fluctuations in single field inflationary models,” JHEP 05, 013 (2003) [arXiv:astro-ph/0210603 [astro-ph]].
  * [53] M. Sasaki, “Large Scale Quantum Fluctuations in the Inflationary Universe,” Prog. Theor. Phys. 76, 1036 (1986).
  * [54] V. F. Mukhanov, “Quantum Theory of Gauge Invariant Cosmological Perturbations,” Sov. Phys. JETP 67, 1297 (1988) [Zh. Eksp. Teor. Fiz. 94N7, 1 (1988)].
  * [55] H. Motohashi, A. A. Starobinsky and J. Yokoyama, “Inflation with a constant rate of roll,” JCAP 1509, 018 (2015) [arXiv:1411.5021 [astro-ph.CO]].
  * [56] T. S. Bunch and P. C. W. Davies, “Quantum Field Theory in de Sitter Space: Renormalization by Point Splitting,” Proc. Roy. Soc. Lond. A 360, 117 (1978).
  * [57] S. Matarrese, S. Mollerach and M. Bruni, “Second order perturbations of the Einstein-de Sitter universe,” Phys. Rev. D 58 (1998), 043504
  * [58] D. Baumann, P. J. Steinhardt, K. Takahashi and K. Ichiki, “Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations,” Phys. Rev. D 76 (2007), 084019 [arXiv:hep-th/0703290 [hep-th]].
  * [59] N. Bartolo, V. De Luca, G. Franciolini, A. Lewis, M. Peloso and A. Riotto, “Primordial Black Hole Dark Matter: LISA Serendipity,” Phys. Rev. Lett. 122 (2019) no.21, 211301 [arXiv:1810.12218 [astro-ph.CO]].
  * [60] R. g. Cai, S. Pi and M. Sasaki, “Gravitational Waves Induced by non-Gaussian Scalar Perturbations,” Phys. Rev. Lett. 122 (2019) no.20, 201101 [arXiv:1810.11000 [astro-ph.CO]].
  * [61] G. Domènech, “Scalar Induced Gravitational Waves Review,” Universe 7 (2021) no.11, 398 [arXiv:2109.01398 [gr-qc]].
  * [62] B. Whitt, “Fourth Order Gravity as General Relativity Plus Matter,” Phys. Lett. B 145 (1984), 176-178
  * [63] A. Karam, N. Koivunen, E. Tomberg, V. Vaskonen and H. Veermäe, “Anatomy of single-field inflationary models for primordial black holes,” [arXiv:2205.13540 [astro-ph.CO]].
  * [64] H. Motohashi and W. Hu, “Primordial Black Holes and Slow-Roll Violation,” Phys. Rev. D 96 (2017) no.6, 063503 doi:10.1103/PhysRevD.96.063503 [arXiv:1706.06784 [astro-ph.CO]].
  * [65] C. T. Byrnes, P. S. Cole and S. P. Patil, “Steepest growth of the power spectrum and primordial black holes,” JCAP 06, 028 (2019) [arXiv:1811.11158 [astro-ph.CO]].
  * [66] P. Carrilho, K. A. Malik and D. J. Mulryne, “Dissecting the growth of the power spectrum for primordial black holes,” Phys. Rev. D 100, no.10, 103529 (2019) [arXiv:1907.05237 [astro-ph.CO]].
  * [67] O. Özsoy and G. Tasinato, “On the slope of the curvature power spectrum in non-attractor inflation,” JCAP 04, 048 (2020) [arXiv:1912.01061 [astro-ph.CO]].
  * [68] P. S. Cole, A. D. Gow, C. T. Byrnes and S. P. Patil, “Steepest growth re-examined: repercussions for primordial black hole formation,” [arXiv:2204.07573 [astro-ph.CO]].
  * [69] V. A. Belinsky, I. M. Khalatnikov, L. P. Grishchuk and Y. B. Zeldovich, “INFLATIONARY STAGES IN COSMOLOGICAL MODELS WITH A SCALAR FIELD,” Phys. Lett. B 155 (1985), 232-236.
  * [70] S. Kachru, R. Kallosh, A. D. Linde and S. P. Trivedi, “De Sitter vacua in string theory,” Phys. Rev. D 68, 046005 (2003) [arXiv:hep-th/0301240 [hep-th]].
  * [71] S. Kachru, R. Kallosh, A. D. Linde, J. M. Maldacena, L. P. McAllister and S. P. Trivedi, “Towards inflation in string theory,” JCAP 0310, 013 (2003) [hep-th/0308055].
  * [72] R. Kallosh and A. Linde, “CMB targets after the latest Planck data release,” Phys. Rev. D 100, no.12, 123523 (2019) [arXiv:1909.04687 [hep-th]].
  * [73] J. Martin, C. Ringeval and V. Vennin, “Encyclopædia Inflationaris,” Phys. Dark Univ. 5-6 (2014), 75-235 [arXiv:1303.3787 [astro-ph.CO]].
  * [74] W. J. Handley, A. N. Lasenby, H. V. Peiris and M. P. Hobson, “Bayesian inflationary reconstructions from Planck 2018 data,” Phys. Rev. D 100 (2019) no.10, 103511 [arXiv:1908.00906 [astro-ph.CO]].
  * [75] S. Rasanen and E. Tomberg, “Planck scale black hole dark matter from Higgs inflation,” JCAP 01 (2019), 038 [arXiv:1810.12608 [astro-ph.CO]].
  * [76] F. Finelli and R. H. Brandenberger, “Parametric amplification of gravitational fluctuations during reheating,” Phys. Rev. Lett. 82 (1999), 1362-1365 [arXiv:hep-ph/9809490 [hep-ph]].
  * [77] K. Jedamzik, M. Lemoine and J. Martin, “Collapse of Small-Scale Density Perturbations during Preheating in Single Field Inflation,” JCAP 09 (2010), 034 [arXiv:1002.3039 [astro-ph.CO]].
  * [78] K. Inomata, M. Braglia and X. Chen, “Questions on calculation of primordial power spectrum with large spikes: the resonance model case,” [arXiv:2211.02586 [astro-ph.CO]].
  * [79] J. Kristiano and J. Yokoyama, “Ruling Out Primordial Black Hole Formation From Single-Field Inflation,” [arXiv:2211.03395 [hep-th]].
  * [80] A. A. Starobinsky, “STOCHASTIC DE SITTER (INFLATIONARY) STAGE IN THE EARLY UNIVERSE,” Lect. Notes Phys. 246, 107-126 (1986) doi:10.1007/3-540-16452-9_6
  * [81] D. S. Salopek and J. R. Bond, “Stochastic inflation and nonlinear gravity,” Phys. Rev. D 43, 1005-1031 (1991) doi:10.1103/PhysRevD.43.1005
  * [82] A. A. Starobinsky and J. Yokoyama, “Equilibrium state of a selfinteracting scalar field in the De Sitter background,” Phys. Rev. D 50, 6357-6368 (1994) [arXiv:astro-ph/9407016 [astro-ph]].
  * [83] T. Fujita, M. Kawasaki, Y. Tada and T. Takesako, “A new algorithm for calculating the curvature perturbations in stochastic inflation,” JCAP 12, 036 (2013) [arXiv:1308.4754 [astro-ph.CO]].
  * [84] T. Fujita, M. Kawasaki and Y. Tada, “Non-perturbative approach for curvature perturbations in stochastic $\delta N$ formalism,” JCAP 10, 030 (2014) [arXiv:1405.2187 [astro-ph.CO]].
  * [85] V. Vennin and A. A. Starobinsky, Eur. Phys. J. C 75, 413 (2015) [arXiv:1506.04732].
  * [86] C. Pattison, V. Vennin, H. Assadullahi and D. Wands, JCAP 10, 046 (2017) [arXiv:1707.00537].
  * [87] J. M. Ezquiaga and J. García-Bellido, “Quantum diffusion beyond slow-roll: implications for primordial black-hole production,” JCAP 08 (2018), 018 [arXiv:1805.06731 [astro-ph.CO]].
  * [88] M. Biagetti, G. Franciolini, A. Kehagias and A. Riotto, “Primordial Black Holes from Inflation and Quantum Diffusion,” JCAP 07 (2018), 032 [arXiv:1804.07124 [astro-ph.CO]].
  * [89] C. Pattison, V. Vennin, H. Assadullahi and D. Wands, “Stochastic inflation beyond slow roll,” JCAP 07, 031 (2019) [arXiv:1905.06300 [astro-ph.CO]].
  * [90] J. M. Ezquiaga, J. García-Bellido and V. Vennin, “The exponential tail of inflationary fluctuations: consequences for primordial black holes,” JCAP 03 (2020), 029 [arXiv:1912.05399 [astro-ph.CO]].
  * [91] H. Firouzjahi, A. Nassiri-Rad and M. Noorbala, “Stochastic Ultra Slow Roll Inflation,” JCAP 01 (2019), 040 [arXiv:1811.02175 [hep-th]].
  * [92] V. Vennin, “Cosmological Inflation: Theoretical Aspects and Observational Constraints,” (PhD Thesis) tel-01094199.
  * [93] G. Ballesteros, J. Rey, M. Taoso and A. Urbano, “Stochastic inflationary dynamics beyond slow-roll and consequences for primordial black hole formation,” JCAP 08 (2020), 043 [arXiv:2006.14597 [astro-ph.CO]].
  * [94] K. Ando and V. Vennin, “Power spectrum in stochastic inflation,” JCAP 04, 057 (2021) [arXiv:2012.02031 [astro-ph.CO]].
  * [95] A. De and R. Mahbub, “Numerically modeling stochastic inflation in slow-roll and beyond,” Phys. Rev. D 102, no.12, 123509 (2020) [arXiv:2010.12685 [astro-ph.CO]].
  * [96] D. G. Figueroa, S. Raatikainen, S. Rasanen and E. Tomberg, “Non-Gaussian Tail of the Curvature Perturbation in Stochastic Ultraslow-Roll Inflation: Implications for Primordial Black Hole Production,” Phys. Rev. Lett. 127, no.10, 101302 (2021) [arXiv:2012.06551 [astro-ph.CO]].
  * [97] D. Cruces and C. Germani, Phys. Rev. D 105, no.2, 023533 (2022) [arXiv:2107.12735].
  * [98] G. Rigopoulos and A. Wilkins, JCAP 12, no.12, 027 (2021) [arXiv:2107.05317].
  * [99] C. Pattison, V. Vennin, D. Wands and H. Assadullahi, “Ultra-slow-roll inflation with quantum diffusion,” JCAP 04, 080 (2021) [arXiv:2101.05741 [astro-ph.CO]].
  * [100] E. Tomberg, “A numerical approach to stochastic inflation and primordial black holes,” J. Phys. Conf. Ser. 2156, no.1, 012010 (2021) [arXiv:2110.10684 [astro-ph.CO]].
  * [101] D. G. Figueroa, S. Raatikainen, S. Rasanen and E. Tomberg, “Implications of stochastic effects for primordial black hole production in ultra-slow-roll inflation,” [arXiv:2111.07437 [astro-ph.CO]].
  * [102] Y. Tada and V. Vennin, “Statistics of coarse-grained cosmological fields in stochastic inflation,” JCAP 02, no.02, 021 (2022) [arXiv:2111.15280 [astro-ph.CO]].
  * [103] R. Mahbub and A. De, “Smooth coarse-graining and colored noise dynamics in stochastic inflation,” [arXiv:2204.03859 [astro-ph.CO]].
  * [104] N. Ahmadi, M. Noorbala, N. Feyzabadi, F. Eghbalpoor and Z. Ahmadi, “Quantum Diffusion in Sharp Transition to Non-Slow-Roll Phase,” JCAP 08 (2022), 078 [arXiv:2207.10578 [gr-qc]].
  * [105] J. H. P. Jackson, H. Assadullahi, K. Koyama, V. Vennin and D. Wands, “Numerical simulations of stochastic inflation using importance sampling,” JCAP 10 (2022), 067 [arXiv:2206.11234 [astro-ph.CO]].
  * [106] M. Celoria, P. Creminelli, G. Tambalo and V. Yingcharoenrat, “Beyond perturbation theory in inflation,” JCAP 06, 051 (2021) doi:10.1088/1475-7516/2021/06/051 [arXiv:2103.09244 [hep-th]].
  * [107] S. Hooshangi, M. H. Namjoo and M. Noorbala, “Rare events are nonperturbative: Primordial black holes from heavy-tailed distributions,” Phys. Lett. B 834 (2022), 137400 [arXiv:2112.04520 [astro-ph.CO]].
  * [108] Y. F. Cai, X. H. Ma, M. Sasaki, D. G. Wang and Z. Zhou, “One small step for an inflaton, one giant leap for inflation: A novel non-Gaussian tail and primordial black holes,” Phys. Lett. B 834 (2022), 137461 [arXiv:2112.13836 [astro-ph.CO]].
  * [109] J. M. Ezquiaga, J. García-Bellido and V. Vennin, “Could ”El Gordo” be hinting at primordial quantum diffusion?,” [arXiv:2207.06317 [astro-ph.CO]].
  * [110] Y. F. Cai, X. H. Ma, M. Sasaki, D. G. Wang and Z. Zhou, “Highly non-Gaussian tails and primordial black holes from single-field inflation,” [arXiv:2207.11910 [astro-ph.CO]].
  * [111] A. D. Gow, H. Assadullahi, J. H. P. Jackson, K. Koyama, V. Vennin and D. Wands, “Non-perturbative non-Gaussianity and primordial black holes,” [arXiv:2211.08348 [astro-ph.CO]].