HIP-2021-31/TH

# Implications of stochastic effects for primordial black hole production in
ultra-slow-roll inflation

Daniel G. Figueroa  , Sami Raatikainen  , Syksy Räsänen  , and Eemeli Tomberg

###### Abstract

We study the impact of stochastic noise on the generation of primordial black
hole (PBH) seeds in ultra-slow-roll (USR) inflation with numerical
simulations. We consider the non-linearity of the system by consistently
taking into account the noise dependence on the inflaton perturbations, while
evolving the perturbations on the coarse-grained background affected by the
noise. We capture in this way the non-Markovian nature of the dynamics, and
demonstrate that non-Markovian effects are subleading. Using the $\Delta N$
formalism, we find the probability distribution $P(\mathcal{R})$ of the
comoving curvature perturbation $\mathcal{R}$. We consider inflationary
potentials that fit the CMB and lead to PBH dark matter with $i)$ asteroid,
$ii)$ solar, or $iii)$ Planck mass, as well as $iv)$ PBHs that form the seeds
of supermassive black holes. We find that stochastic effects enhance the PBH
abundance by a factor of $\mathcal{O}(10)-\mathcal{O}(10^{8})$, depending on
the PBH mass. We also show that the usual approximation, where stochastic
kicks depend only on the Hubble rate, either underestimates or overestimates
the abundance by orders of magnitude, depending on the potential. We evaluate
the gauge dependence of the results, discuss the quantum-to-classical
transition, and highlight open issues of the application of the stochastic
formalism to USR inflation.

##  1 Introduction

##### Stochastic ultra-slow-roll inflation and primordial black holes.

Inflation is the most successful scenario to explain the primordial universe
[1], predicting that the average spatial curvature is negligible, and that the
primordial perturbations observed in the cosmic microwave background (CMB) and
large scale matter perturbations are almost scale-invariant, highly Gaussian,
mostly adiabatic, dominated by a growing scalar mode, and statistically very
homogeneous and isotropic [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16]. In the simplest inflationary models, a scalar field called the inflaton
rolls down its potential slowly, with Hubble friction and potential push
balanced, so that its acceleration is negligible, like for a falling
parachute. This is called slow-roll (SR).

If the potential has a very flat section or a shallow minimum, the potential
push is instead negligible, and the inflaton acceleration is dominated by
friction, leading to an exponentially falling velocity, like for a brick
sliding on a flat surface. Such a regime is called ultra-slow-roll (USR).
Perturbations generated in USR are far from scale-invariant, and can be highly
non-Gaussian. Therefore, the inflaton cannot be in USR when the modes observed
in the CMB are generated. However, as the inflaton perturbations are enhanced
by the fact that its background velocity is small, USR outside the CMB region
could generate large perturbations that seed primordial black holes (PBH) [17,
18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30]. PBHs are a longstanding
dark matter candidate [31, 32, 33, 34, 35, 36, 37, 38], which can constitute
all of the dark matter, or they can act as seeds for supermassive black holes
[39, 40, 41]. Large inhomogeneities on small scales could also affect big bang
nucleosynthesis [42] and recombination [43, 44].

During inflation, perturbations are constantly created on sub-Hubble scales.
Due to accelerated expansion, perturbations are stretched beyond the Hubble
radius, creating a phase space populated by modes with exponentially different
wavelengths. In the leading approximation, perturbations evolve according to
linear perturbation equations around a Friedmann–Lemaître–Robertson–Walker
(FLRW) background [45, 46]. Beyond the leading approximation, modes with
different wavelengths couple. Small wavelength modes can be treated as an
environment for the system of large wavelength modes. On super-Hubble scales,
modes are squeezed and their probability distribution classicalises [47, 48,
49, 50, 51, 52, 53, 54, 55]. Once the originally quantum small wavelength
modes become super-Hubble, they can be coarse-grained and treated as
stochastic noise for the effectively classical long wavelength modes. This
effect is captured by the formalism of stochastic inflation, first introduced
in [56] and further developed in [57, 47, 58, 59, 50, 60, 61, 62, 63, 64, 65,
66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84,
85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97].

In typical SR inflationary scenarios stochastic effects are small, though
there can be exceptions, especially when inflation lasts much longer than 60
e-folds [61, 62, 87]. However, in USR stochastic effects can be expected to be
important. This is due to three reasons. First, as the potential push is
absent and the field velocity decays exponentially, the field evolution is
more sensitive to stochastic kicks, especially if USR lasts long. Second, the
USR regime can produce larger perturbations than standard SR, and hence give
stronger stochastic kicks. Third, the main observable related to PBH abundance
– the prevalence of high density regions – is exponentially sensitive to the
amplitude of the tail of the curvature perturbation distribution, where the
usual Gaussian contribution is highly suppressed and the strong non-
Gaussianities from USR may be important111As pointed out in [87], stochastic
effects can also have a large effect on the CMB even when the USR region of
the potential is far from the CMB region, because some rare modes can remain
in the part of the potential where stochastic effects are important for so
long that their wavelength stretches to CMB scales..

The effect of stochastic kicks on the power spectrum in USR has been studied
in [84, 27, 28, 90, 91, 74, 29, 94, 87, 93, 96, 97], see also [84, 91, 93] for
higher moments. Using the $\Delta N$ formalism [98, 99, 46, 100] it has been
demonstrated that stochastic effects can produce an exponential tail in the
probability distribution $P(\mathcal{R})$ of the curvature perturbation
$\mathcal{R}$, which dominates over the Gaussian tail expected from linear
theory [85, 89]. As PBHs are produced from rare perturbations, the amplitude
of the power spectrum may then not be enough to determine the PBH abundance,
even when the power spectrum correctly describes perturbations close to the
mean.

The calculations in [85, 89] were in the context of SR, and in our recent
Letter [88] – which we refer to as $Paper\,I$ – we presented the first
stochastic calculation in USR inflation. We calculated numerically the
probability distribution $P(N)$ for the number of e-folds $N$ of inflation in
the $\Delta N$ formalism, which gives $P(\mathcal{R})$. We simulated a large
number of realisations of the stochastic dynamics, choosing a realistic
inflationary potential tailored to create perturbations during USR that lead
to PBHs with mass $7\times 10^{-15}M_{\odot}$, roughly the mass of the
asteroid Eros. After our results appeared, the form of $P(N)$ in stochastic
USR, including an exponential tail, was calculated analytically [95]. Our
calculation in $Paper\,I$ included non-Markovian backreaction effects for the
first time, incorporated in the dynamics by solving the stochastic noise self-
consistently for each realisation. We considered the non-linearity of the
system by taking into account that the noise depends on the inflaton
perturbations, while the perturbations evolve on the coarse-grained background
affected by the noise. In [68, 69, 73], such backreaction effects had been
considered recursively. However, to our knowledge, to date no calculation
apart from $Paper\,I$ has included the non-Markovian effects at every
timestep, realisation by realisation. We obtained the full probability
distribution, all the way from the Gaussian peak to the exponential tail for
large $\Delta N$ values where PBHs form (as well as negative values of $\Delta
N$). We found that in the model studied (referred to in the present paper as
the asteroid case), the exponential tail led to an enhancement of the observed
PBH abundance by a factor of $\sim 10^{5}$ over the Gaussian approximation
based on the power spectrum. We also found that the non-trivial evolution of
the perturbations enhances the PBH abundance by orders of magnitude, when
compared to standard stochastic calculations where the mode functions are
fixed to their SR super-Hubble limit, which leads to noise proportional to the
Hubble parameter. Our numerical results showed that such simplifying
assumptions, used in studies of stochastic inflation from the original work
[56] to the most recent papers [95], may be inaccurate in USR. We now extend
this conclusion to other USR potentials.

In the present paper we return to asteroid mass PBHs, and extend our numerical
treatment to other observationally interesting mass scales for PBH dark matter
[39, 40, 41]. We consider PBHs with a mass of $5M_{\odot}$ and Planck scale
relic PBHs with an initial mass of $10^{3}$ kg. We also consider PBHs that act
as seeds for supermassive black holes in the early universe, with an initial
mass of $2000M_{\odot}$. As in $Paper\,I$, we build the inflationary
potentials starting from the scenario where the Standard Model Higgs is the
inflaton [101, 102], and use the chiral Standard Model renormalisation group
running to create a shallow minimum that leads to USR [25] (see [18, 24, 19,
103, 20] for similar constructions.) In the Higgs inflation model, a potential
that produces PBHs with the right abundance (according to the Gaussian
approximation) gives a spectral index that is too red on the CMB scales,
except in the case of Planck mass relics. Hence, as in $Paper\,I$, we adjust
the large field part of the potentials by hand to fit the CMB. Our models are
therefore only inspired by Higgs inflation, not real examples of Higgs
inflation (except in the relic case).

In section 2 we go over the formalism of stochastic inflation. We highlight
open issues and ambiguities in the application of the stochastic equations to
the USR regime. In section 3 we first verify our algorithm by applying it to
toy models discussed in [89]. We apply our numerical machinery to more
realistic potentials that simultaneously fit the CMB constraints on the
observable part of the spectrum and produce PBHs with the right abundance (in
the Gaussian approximation). In section 4 we summarise our findings. In
appendix A we give the details of how we construct our inflationary
potentials. In appendix B we explain how we discretise and solve numerically
the relevant equations. Finally, in appendix C we discuss gauge issues.

##  2 Stochastic $\Delta N$ formalism

###  2.1 Coarse-graining and noise

####  2.1.1 Background equations

We consider stochastic effects in single field inflation. We follow the
notation of [74]. The action is

$\displaystyle
S=\int\mathrm{d}^{4}x\sqrt{-g}\left[\frac{1}{2}M_{{}_{\mathrm{Pl}}}^{2}R-\frac{1}{2}g^{\alpha\beta}\partial_{\alpha}\phi\partial_{\beta}\phi-V(\phi)\right]\
,$ (2.1)

where $M_{{}_{\mathrm{Pl}}}$ is the reduced Planck mass, $R$ is the Ricci
scalar and $\phi$ is the inflaton. We choose units such that
$M_{{}_{\mathrm{Pl}}}=1$. We consider a spatially flat
Friedmann–Robertson–Lemaître–Walker universe with scalar perturbations, so the
metric is

$\displaystyle\mathrm{d}s^{2}=-(1+2A)\mathrm{d}t^{2}+2a(t)\partial_{i}B\mathrm{d}x^{i}\mathrm{d}t+a(t)^{2}[(1-2\psi)\delta_{ij}+2\partial_{i}\partial_{j}E]\mathrm{d}x^{i}\mathrm{d}x^{j}\
.$ (2.2)

We decompose the inflaton into long and short wavelength parts as

$\displaystyle\phi(t,\mathbf{x})$ $\displaystyle=$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{\frac{3}{2}}}\phi_{\mathbf{k}}(t)e^{-i\mathbf{k}\cdot\mathbf{x}}$
(2.3) $\displaystyle=$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{\frac{3}{2}}}\left[1-W\left(\frac{k}{\sigma
aH}\right)\right]\phi_{\mathbf{k}}(t)e^{-i\mathbf{k}\cdot\mathbf{x}}+\int\frac{\mathrm{d}^{3}k}{(2\pi)^{\frac{3}{2}}}W\left(\frac{k}{\sigma
aH}\right)\phi_{\mathbf{k}}(t)e^{-i\mathbf{k}\cdot\mathbf{x}}$
$\displaystyle\equiv$
$\displaystyle\bar{\phi}(t,\mathbf{x})+\delta\phi(t,\mathbf{x})\ ,$

where $k\equiv|\mathbf{k}|$, $H\equiv\dot{a}/a$ is the Hubble rate (dot
denotes derivative with respect to $t$), and $W\left(\frac{k}{\sigma
aH}\right)$ is a window function: $W(x)\simeq 0$ for $x\ll 1$ and $W(x)\simeq
1$ for $x\gg 1$. The parameter $\sigma$ controls the separation scale
$k_{\mathrm{c}}\equiv\sigma aH$, and modes with $k\ll k_{\mathrm{c}}$ are
considered long wavelength, and those with $k\gg k_{\mathrm{c}}$ are
considered short wavelength. In the following, we consider patches of size
$\sim k_{\mathrm{c}}^{-1}$. Within each patch, the value of the long-
wavelength field $\bar{\phi}$ is approximately constant.

As the universe expands, modes are stretched to super-Hubble scales. As a
result, the long-wavelength background over which short-scale modes evolve
changes in time. Coarse-graining the short wavelength environment to get an
effective theory for the long wavelength system leads to quantum noise that
can be treated as classical stochastic noise if the probability distribution
of the perturbations is classical. This requires the wavelength of the modes
to be much larger than the Hubble radius, $\sigma\ll 1$ (we discuss
classicalisation in section 2.2). The value of $\sigma$ gives an offset
between the time a mode exits the Hubble radius and the time it is coarse-
grained and gives a stochastic kick to the background. In SR, modes freeze in
the super-Hubble regime, so the stochastic results are insensitive to the
choice of $\sigma$, as long as it is sufficiently small that modes have
stopped evolving [56] (but not too small [58, 69, 83, 82, 104, 86]), as the
decaying mode dies out and the growing mode stops evolving quickly after the
Hubble exit.

This is not the case in USR. Instead, the decaying mode turns into the growing
mode, leading to a growth of the comoving curvature perturbation $\mathcal{R}$
on super-Hubble scales, and consequently to the generation of large
perturbations that can potentially seed PBHs. Growth of $\mathcal{R}$ on
super-Hubble scales due to violation of SR is well-known [105, 106, 107, 108,
109, 110], see e.g. [25] for discussion of this in USR. The longer USR lasts,
the larger are the scales to which the growth of $\mathcal{R}$ extends. This
makes the results potentially sensitive to the choice of window function and
the value of $\sigma$. Dependence of the noise on the window function has been
discussed in [57], see also [53, 111, 112, 73, 83]. The validity of the choice
of $\sigma$ (more generally, the form of the Langevin equation) would
ultimately have to be checked with a first principle derivation of the
separation between system and environment in quantum field theory. As it is
customary, we choose a step function in momentum space,
$W\left(\frac{k}{\sigma aH}\right)=\theta(k-\sigma aH)$. This corresponds to a
position space window function that takes negative values, and underestimates
the rate at which the two-point function falls as a function of spatial
distance [113]. For our purposes this should not make a difference, as we do
not consider spatial correlations.

We follow the usual heuristic derivation of the stochastic equations. We start
from the FLRW background equations of motion for the action (2.1) [45, 46],

$\displaystyle\frac{\partial\bar{\phi}}{\partial t}$ $\displaystyle=$
$\displaystyle\bar{\pi}_{t}$ (2.4)
$\displaystyle\frac{\partial\bar{\pi}_{t}}{\partial t}$ $\displaystyle=$
$\displaystyle-3H\bar{\pi}_{t}-V_{,\bar{\phi}}$ (2.5) $\displaystyle 3H^{2}$
$\displaystyle=$ $\displaystyle\frac{1}{2}\bar{\pi}_{t}^{2}+V\ ,$ (2.6)

where the subscript $t$ indicates that the field momentum is defined with
respect to the time derivative. We rewrite the $t$ derivatives as derivatives
with respect to the number of e-folds $N\equiv\ln a/a_{*}$, where $*$ refers
to the time when the CMB pivot scale $k_{*}=0.05$ Mpc-1 crosses the Hubble
radius (note that in our convention $N$ grows with time). We then replace the
FLRW field $\bar{\phi}$ with the full field $\phi=\bar{\phi}+\delta\phi$
(where $\bar{\phi}$ is not the FLRW field but the coarse-grained field). We do
the same for the momentum, and denote the momentum defined with respect to the
e-fold derivative by $\pi\equiv\pi_{t}/H$. We keep the short wavelength part
in the time derivatives only, reinterpreted as noise. All in all, we get the
equations

$\displaystyle\bar{\phi}^{\prime}$ $\displaystyle=$
$\displaystyle\bar{\pi}+\xi_{\phi}$ (2.7) $\displaystyle\bar{\pi}^{\prime}$
$\displaystyle=$
$\displaystyle-\left(3+\frac{H^{\prime}}{H}\right)\bar{\pi}-\frac{V_{,\bar{\phi}}}{H^{2}}+\xi_{\pi}$
(2.8) $\displaystyle 2V$ $\displaystyle=$
$\displaystyle(6-\bar{\pi}^{2})H^{2}\ ,$ (2.9)

where prime denotes derivative with respect to $N$, and $\xi_{\phi}$ and
$\xi_{\pi}$ are noise terms. It is important to use $N$ as the time coordinate
when defining the noise [64, 114, 84, 74]. Physics is invariant under a change
of coordinates, but the noise terms are not. In this treatment, they are

$\displaystyle\xi_{\phi}$ $\displaystyle\equiv$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\frac{\mathrm{d}W}{\mathrm{d}N}\phi_{\mathbf{k}}(N)e^{-i\mathbf{k}\cdot\mathbf{x}}=-\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\frac{\mathrm{d}(\sigma
aH)}{\mathrm{d}N}\delta(k-\sigma
aH)\delta\phi_{\mathbf{k}}(N)e^{-i\mathbf{k}\cdot\mathbf{x}}$ (2.10)
$\displaystyle\xi_{\pi}$ $\displaystyle\equiv$
$\displaystyle\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\frac{\mathrm{d}W}{\mathrm{d}N}\phi^{\prime}_{\mathbf{k}}(N)e^{-i\mathbf{k}\cdot\mathbf{x}}=-\int\frac{d^{3}k}{(2\pi)^{\frac{3}{2}}}\frac{\mathrm{d}(\sigma
aH)}{\mathrm{d}N}\delta(k-\sigma
aH)\delta\phi^{\prime}_{\mathbf{k}}(N)e^{-i\mathbf{k}\cdot\mathbf{x}}\,.$
(2.11)

The average over the noise is zero, and its two-point function is obtained by
quantising $\phi$. We write

$\displaystyle\hat{\phi}_{\mathbf{k}}$ $\displaystyle\equiv$
$\displaystyle\phi_{k}\hat{a}_{\mathbf{k}}+\phi^{*}_{k}\hat{a}^{\dagger}_{\mathbf{-k}}$
(2.12) $\displaystyle\hat{\pi}_{\mathbf{k}}$ $\displaystyle\equiv$
$\displaystyle\phi^{\prime}_{k}\hat{a}_{\mathbf{k}}+\phi^{\prime*}_{k}\hat{a}^{\dagger}_{\mathbf{-k}}\
,$ (2.13)

where the creation and annihilation operators obey the usual commutation
relations

$\quantity[\hat{a}_{\mathbf{k}},\hat{a}^{\dagger}_{\mathbf{k^{\prime}}}]=\delta^{(3)}(\mathbf{k}-\mathbf{k^{\prime}})\
,\quad\quantity[\hat{a}_{\mathbf{k}},\hat{a}_{\mathbf{k^{\prime}}}]=\quantity[\hat{a}^{\dagger}_{\mathbf{k}},\hat{a}^{\dagger}_{\mathbf{k^{\prime}}}]=0\
,$ (2.14)

and the vacuum state satisfies $\hat{a}_{\mathbf{k}}\ket{0}=0$.

The resulting quantum noise correlators from (2.10)–(2.14) are:

$\displaystyle\expectationvalue{\hat{\xi}_{\phi}(N)\hat{\xi}_{\phi}(N^{\prime})}{0}$
$\displaystyle=$ $\displaystyle\frac{1}{6\pi^{2}}\frac{\mathrm{d}(\sigma
aH)^{3}}{\mathrm{d}N}|\delta\phi_{k=\sigma aH}|^{2}\delta(N-N^{\prime})$
(2.15)
$\displaystyle\expectationvalue{\hat{\xi}_{\pi}(N)\hat{\xi}_{\pi}(N^{\prime})}{0}$
$\displaystyle=$ $\displaystyle\frac{1}{6\pi^{2}}\frac{\mathrm{d}(\sigma
aH)^{3}}{\mathrm{d}N}|\delta\phi^{\prime}_{k=\sigma
aH}|^{2}\delta(N-N^{\prime})$ (2.16)
$\displaystyle\expectationvalue{\hat{\xi}_{\phi}(N)\hat{\xi}_{\pi}(N^{\prime})}{0}$
$\displaystyle=$ $\displaystyle\frac{1}{6\pi^{2}}\frac{\mathrm{d}(\sigma
aH)^{3}}{\mathrm{d}N}\delta\phi_{k=\sigma aH}\delta\phi^{\prime*}_{k=\sigma
aH}\delta(N-N^{\prime})\ .$ (2.17)

The physical reason for the noise is that, unlike in standard perturbation
theory, the split between the background and perturbations is time-dependent,
and modes cross from small wavelengths into large wavelengths. (As $\sigma\ll
1$, rather than small and large wavelength modes, it would be more accurate to
talk about infrared and far infrared modes.) Choosing the window function to
be a step function offers considerable simplification, as each mode
contributes to the noise only when it crosses the coarse-graining radius
$1/(\sigma aH)$, and (2.7), (2.8) become Langevin equations. For other window
functions the noise involves correlations between different $k$ values.

There are several detailed analyses of the stochastic equations, including
derivations from and comparisons to quantum field theory [57, 59, 47, 61, 62,
63, 64, 65, 68, 66, 67, 72, 73, 115, 75, 76, 77, 78, 84, 79, 80, 81, 82, 104,
86, 74]. However, none with explicit Langevin equations apply to the USR case,
where the results can be sensitive to the choice of the window function, which
for the step function reduces to the choice of $\sigma$. The value of $\sigma$
has to be sufficiently small for the probability distribution of the short
modes to have become classical. Beyond that, different choices correspond to
different splits between the environment and the system, and hence to
different stochastic equations. It is likely that the correct value for the
equations (2.7)–(2.8) and (2.15)–(2.17) cannot be determined heuristically,
but requires a quantum field theory calculation applied to the USR case. A
proper derivation could also modify the form of the Langevin equation beyond
the choice of $\sigma$. Such analysis goes beyond the scope of our paper, and
we simply take a step function and choose a value of $\sigma$ for which the
probability distribution of the modes has become sufficiently classical. We
also consider smaller and larger values of $\sigma$, and find that the results
are rather insensitive to the precise value.

####  2.1.2 Perturbation equation

We also have to consider the evolution of the modes, which determine the
amplitudes of the stochastic kicks. In most explicit studies of stochastic
inflation, the modes are taken to evolve on a non-stochastic FLRW background.
Furthermore, the mode functions are fixed to their asymptotic SR and super-
Hubble limit, $|\delta\phi_{\mathbf{k}}|=H/(\sqrt{2}k^{3/2})$, and the time
derivative of $H$ in the noise amplitude (2.15) is neglected. The momentum
noise then vanishes, and the field noise correlator becomes

$\displaystyle\expectationvalue{\hat{\xi}_{\phi}(N)\hat{\xi}_{\phi}(N^{\prime})}{0}$
$\displaystyle=$
$\displaystyle\left(\frac{H}{2\pi}\right)^{2}\delta(N-N^{\prime})\ ,$ (2.18)

where $H$ is the FLRW Hubble rate. This simplifies the problem considerably,
because there is no need to follow the evolution of the modes, and the
stochastic process becomes Markovian. However, for consistency, the
backreaction of non-trivial mode evolution and change of the environment
should be included, particularly in USR where modes can evolve after Hubble
crossing. This has been done recursively in [68, 69, 73], and consistently at
every timestep in $Paper\,I$.

In standard linear perturbation theory, the short wavelength modes obey linear
perturbation equations over the FLRW background [45, 46]. In the spatially
flat gauge ($\psi=E=0$), these equations reduce to

$\frac{\partial^{2}\delta\phi_{\mathbf{k}}}{\partial
t^{2}}+3H\frac{\partial\delta\phi_{\mathbf{k}}}{\partial
t}+\quantity[\frac{k^{2}}{a^{2}}+V^{\prime\prime}(\bar{\phi})-\frac{1}{a^{3}}\frac{\mathrm{d}}{\mathrm{d}t}\quantity(\frac{a^{3}}{H}\bar{\pi}^{2})]{\delta\phi}_{\mathbf{k}}=0\
,$ (2.19)

where the barred symbols are FLRW background quantities, and
$\delta\phi_{\mathbf{k}}$ is the perturbation around the FLRW field (or the
mode function in the quantum operator (2.12)). We include the change of the
coarse-grained background in the evolution of the small wavelength modes by
assuming that the perturbation $\delta\phi_{\mathbf{k}}$ follows (2.19) with
the FLRW fields replaced by the coarse-grained fields. Switching to
derivatives with respect to the number of e-folds, the resulting equation for
the short wavelength modes is

$\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\delta\phi_{\mathbf{k}}^{\prime\prime}+\left(3+\frac{H^{\prime}}{H}\right)\delta\phi_{\mathbf{k}}^{\prime}+\omega_{k}^{2}\delta\phi_{\mathbf{k}}=0\
,\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ $ (2.20) $\displaystyle{\rm where}$
(2.21) $\displaystyle\leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\omega_{k}^{2}\equiv{k^{2}\over(aH)^{2}}+\bar{\pi}^{2}\left(3+2{H^{\prime}\over
H}-{H^{\prime}\over H^{2}}\right)+2\bar{\pi}{V_{,\bar{\phi}}\over
H^{2}}+{V_{,\bar{\phi}\bar{\phi}}\over H^{2}}\,.\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
\leavevmode\nobreak\ $

As in standard perturbation theory, we assume that the initial conditions in
the deep sub-Hubble regime are given by the Bunch–Davies vacuum (we discuss
the vacuum in more detail in section 2.2),

$\delta\phi_{\mathbf{k}}=\frac{1}{\sqrt{2k}a}\
,\qquad(a\delta\phi_{\mathbf{k}})^{\prime}=-i\frac{k}{H}\delta\phi_{\mathbf{k}}\
.$ (2.22)

####  2.1.3 Remarks

A few comments are in order. In (2.9) we have chosen the first Friedmann
equation not to receive noise corrections. This is usual in stochastic
studies. However, as the continuity equation (i.e. the field equation of
motion) has noise, both Friedmann equations cannot be noise-free. Therefore,
given our choice for the first Friedmann equation, the second Friedmann
equation must have noise. One reason for our choice is that the first
Friedmann equation written in terms of $H$ does not involve time derivatives,
which in the heuristic derivation of the Langevin equation are the origin of
the noise. Still, in order to derive the second order field equation (2.19)
from the perturbed Friedmann equations, it is necessary to use the second
Friedmann equation. If we added noise directly to the Einstein equations
perturbed around the FLRW universe, and derived the second order field
equation by taking derivatives as usual, the resulting stochastic counterpart
to the mode equation (2.20) would involve noise explicitly. Correspondingly,
the noise amplitudes (2.15) and (2.16) would themselves depend on the noise.
We choose to avoid such complications, and take factors of $H^{\prime}$
appearing in the perturbation equation to be related to the coarse-grained
fields, in the same way as in the FLRW case,
$H^{\prime}/H=-\frac{1}{2}\bar{\pi}^{2}$. While the noise correction to
$H^{\prime}/H$ can be one order of magnitude larger than the classical
contribution (which is strongly suppressed in USR), it is always orders of
magnitude smaller than unity. Because $H^{\prime}/H$ only appears in the
combination $3+H^{\prime}/H$, we expect this noise correction to have
negligible impact on the calculation.

In addition to the full non-Markovian calculation, for the sake of comparison
we also do a simplified calculation where the kicks are given by (2.18),
corresponding to the asymptotic SR and long-wavelength behaviour of the mode
functions (and also neglecting a time derivative of $H$). We find that this
approximation, commonly used in the literature, severely underestimates or
overestimates the PBH abundance, depending on the potential. We discuss the
details in section 3.2. We also use a recursive treatment where the modes are
first taken to evolve on the FLRW background, and the noise given by these
modes determines the coarse-grained background for the next iteration, as in
[68, 69, 73]. We find that one iteration is sufficient to reproduce the
results of the full treatment.

As the stochastic formalism goes beyond perturbation theory, the correct
equations cannot be reliably deduced using the heuristic arguments we have
given, but should be derived from quantum field theory. Our prescription
amounts to using the equations for perturbations around homogeneous and
isotropic background also for perturbations around the coarse-grained
background field. There could be couplings between the perturbations and the
modes that form the mean field that our treatment does not capture. Note that
while all the modes counted as part of the background are super-Hubble and
classical, the background is not exactly homogeneous and isotropic. Modes with
$k$ close to the Hubble scale will significantly evolve again during USR,
while modes with smaller $k$ will not, leading to a time-dependent gradient in
the background.

We also have to consider gauge-dependence. The stochastic equations
(2.7)–(2.8) for $\bar{\phi}$ have been formulated in the uniform $N$ gauge,
whereas the equation of motion (2.20) for the perturbations is in the
spatially flat gauge. It was shown in [74] that the correction from the gauge
transformation is small both in SR and deep in USR. However, this is not
guaranteed to be the case during the transition between the two regimes, or
when the potential makes a sizeable contribution during USR. We have verified
numerically that for the inflationary potentials we use to produce asteroid,
solar mass and supermassive PBHs, the correction from the gauge
transformation, computed as the ratio of the correction over the mode
amplitude, is less than $10^{-5}$ for most modes. There is, however, a narrow
spike of modes with very small amplitude at the coarse-graining scale, for
which the gauge correction rises to a maximum of $\sim 10^{-2}$. These modes
exit the coarse-graining scale roughly halfway in the USR region, when the
field has just passed the local minimum in the potential. The contribution of
these modes to the overall dynamics is small. For details, see appendix C. In
the relic case, the situation is different and the gauge transformations have
a large impact also on modes outside the spikes. As that case also suffers
from other technical and conceptual problems, we postpone discussion of this
issue to section 3.2.5.

Note that the uniform $N$ gauge has a residual gauge freedom, which is not
reflected in the Langevin equation. Under such a residual gauge
transformation, the perturbation equations change, while the Langevin equation
does not, which is in principle inconsistent. The system of equations of large
and small wavelengths as a whole should be gauge-invariant, and a derivation
from first principles in quantum field theory should determine its correct
form. Such a derivation goes beyond the scope of this paper.

Finally, defining stochastic equations involves specifying the discretisation,
for example choosing between the Itô and the Stratonovich prescription. This
issue for stochastic inflation is discussed in detail in [116, 86], but it is
unimportant for single field inflation with a sub-Planckian potential, $V\ll
1$, such as ours. Our exact numerical implementation, including the
discretisation details of the evolution equations, are given in appendix B.

###  2.2 Squeezing, freezing, and classicalisation

As noted earlier, in order for the stochastic formalism to apply, the coarse-
graining parameter $\sigma$ has to be small enough that the probability
distribution of the perturbations has become classical by the time they join
the background. As modes cross the Hubble scale, they are squeezed, and the
expectation values of powers of fields and momenta approach those of a
classical stochastic distribution [117, 118, 119, 120, 121, 122, 123].
Squeezing is due to the rapid decline of the decaying mode of the
perturbation, and leads to what has been called “decoherence without
decoherence” [49], since it leads to classicalisation in the sense explained
above, but not to proper decoherence, i.e. suppression of the off-diagonal
elements of the density matrix [117, 124, 54, 52, 125, 126, 127, 128, 55,
129]. Therefore, squeezing captures only part of decoherence and
classicalisation, and the entire process is not fully understood [130, 131,
132, 133, 134, 135, 136]. We will be content to use the degree of squeezing as
a criterion for classicalisation, as in [49]. This approach can be considered
conservative in the sense that any decoherence due to interactions is expected
to only increase the degree of classicality.

To simplify the discussion and make contact with the standard squeezing
formalism, in this section we work with a discrete set of $\mathbf{k}$ modes
instead of a continuum, replacing Fourier transforms such as (2.3) with
Fourier series, and Dirac delta functions such as those in the ladder
commutator (2.14) with Kronecker deltas. This corresponds to working in a
finite box with periodic boundary conditions. The operators
$\hat{\phi}_{\mathbf{k}}$, $\hat{\pi}_{\mathbf{k}}$ mix the $\mathbf{k}$ and
$-\mathbf{k}$ modes and are non-Hermitian. We obtain two independent Hermitian
field operators for each $(\mathbf{k},-\mathbf{k})$ pair by defining the real
and imaginary parts
$\hat{\phi}^{R}_{\mathbf{k}}\equiv\frac{1}{\sqrt{2}}(\hat{\phi}_{\mathbf{k}}+\hat{\phi}_{\mathbf{k}}^{\dagger})$,
$\hat{\phi}^{I}_{\mathbf{k}}\equiv-\frac{i}{\sqrt{2}}(\hat{\phi}_{\mathbf{k}}-\hat{\phi}_{\mathbf{k}}^{\dagger})$,
and similarly for $\hat{\pi}_{\mathbf{k}}$ [121]. They behave identically, and
in the following we concentrate on one of them, dropping the superscript for
brevity.

In the Schrödinger picture, a squeezed state can be written as (we omit a
phase which does not contribute to the expectation values we consider) [49,
123]

$\displaystyle\ket{\psi_{\mathbf{k}}}=\exp[\frac{1}{2}\quantity(s_{\mathbf{k}}^{*}\hat{c}_{\mathbf{k}}^{2}-s_{\mathbf{k}}\hat{c}_{\mathbf{k}}^{\dagger}{}^{2})]\ket{0}\
,$ (2.23)

where $s_{\mathbf{k}}=r_{k}e^{2i\varphi_{\mathbf{k}}}$ is the squeezing
parameter and $\hat{c}^{\dagger}_{\mathbf{k}}$ and $\hat{c}_{\mathbf{k}}$ are
creation and annihilation operators. The amplitude $r_{k}$ indicates how
squeezed the state is, and the phase $\varphi_{\mathbf{k}}$ gives the
squeezing direction in phase space. The $\hat{c}_{\mathbf{k}}$ operators
determine the vacuum state with respect to which squeezing is measured. They
are related in the standard way to the system’s canonical variables
$\hat{q}_{\mathbf{k}}$ and $\hat{p}_{\mathbf{k}}$, which we choose to be

$\displaystyle\hat{q}_{\mathbf{k}}$ $\displaystyle=$
$\displaystyle\frac{1}{\sqrt{2}}\quantity(\hat{c}_{\mathbf{k}}+\hat{c}^{\dagger}_{\mathbf{k}})=\sqrt{k}a\hat{\phi}_{\mathbf{k}}$
(2.24) $\displaystyle\hat{p}_{\mathbf{k}}$ $\displaystyle=$
$\displaystyle-\frac{i}{\sqrt{2}}\quantity(\hat{c}_{\mathbf{k}}-\hat{c}^{\dagger}_{\mathbf{k}})=\frac{Ha^{2}}{\sqrt{k}}\hat{\pi}_{\mathbf{k}}\
.$ (2.25)

Deep inside the Hubble radius the state annihilated by $\hat{c}_{\mathbf{k}}$
is the adiabatic Bunch–Davies vacuum, and the creation and annihilation
operators coincide with $\hat{a}^{\dagger}_{\mathbf{k}}$ and
$\hat{a}_{\mathbf{k}}$ in (2.12) (up to the normalisation difference between
the continuum and discrete case). When $k\lesssim aH$, this is no longer the
case; in fact, there is no adiabatic vacuum at all. The choice (2.24)–(2.25)
gives the closest match to the usual discussion of squeezing for the
Sasaki–Mukhanov variable, which in the spatially flat gauge reduces to
$a\phi_{\mathbf{k}}$.222The Bunch–Davies vacuum does not uniquely fix the
canonical variables. If we add terms that are subleading in $(aH)/k$, they
vanish in the deep sub-Hubble limit where the vacuum is defined, but can be
important for super-Hubble wavelengths. Our choice is consistent with the
Bunch–Davies vacuum, but not determined by it.

From (2.23) and (2.24)–(2.25) we get the correlators [123]

$\displaystyle\Sigma_{{}_{qq}}(\mathbf{k})$ $\displaystyle\equiv$
$\displaystyle\expectationvalue{\hat{q}_{\mathbf{k}}^{2}}{\psi_{\mathbf{k}}}=\frac{1}{2}\quantity[\cosh(2r_{k})-\cos(2\varphi_{k})\sinh(2r_{k})]=ka^{2}|\phi_{k}|^{2}$
(2.26) $\displaystyle\Sigma_{{}_{pp}}(\mathbf{k})$ $\displaystyle\equiv$
$\displaystyle\expectationvalue{\hat{p}_{\mathbf{k}}^{2}}{\psi_{\mathbf{k}}}=\frac{1}{2}\quantity[\cosh(2r_{k})+\cos(2\varphi_{k})\sinh(2r_{k})]=\frac{H^{2}a^{4}}{k}|\phi^{\prime}_{k}|^{2}$
(2.27) $\displaystyle\Sigma_{{}_{qp}}(\mathbf{k})$ $\displaystyle\equiv$
$\displaystyle\frac{1}{2}\expectationvalue{\hat{q}_{\mathbf{k}}\hat{p}_{\mathbf{k}}+\hat{p}_{\mathbf{k}}\hat{q}_{\mathbf{k}}}{\psi_{\mathbf{k}}}=-\frac{1}{2}\sin(2\varphi_{k})\sinh(2r_{k})=a^{3}H\real(\phi_{k}^{*}\phi^{\prime}_{k})\
,$ (2.28)

where we have made contact with the Heisenberg picture formalism by expressing
the correlators in terms of the mode functions $\phi_{k}$ used in (2.12).
These expectation values are generically large if the squeezing amplitude
$r_{k}$ is large. However, for particular values of the squeezing angle
$\varphi_{k}$, some of the expectation values are small even for large
$r_{k}$. The sum
$\Sigma_{{}_{qq}}(\mathbf{k})+\Sigma_{{}_{pp}}(\mathbf{k})=\cosh(2r_{k})$ is
independent of $\varphi_{k}$, so it is always large when the state is very
squeezed.

For the Bunch–Davies vacuum, the mode is initially described by the minimum
uncertainty wave packet, for which $r_{k}=0$. As the mode stretches further
outside the Hubble radius, the phase space probability distribution gets
squeezed and $r_{k}$ grows. Because
$n_{\mathbf{k}}\equiv\expectationvalue{\hat{c}^{\dagger}_{\mathbf{k}}\hat{c}_{\mathbf{k}}}{\psi_{\mathbf{k}}}=\frac{1}{2}[1+\cosh(2r_{k})]$,
the condition $r_{k}\gg 1$ corresponds to large occupation number
$n_{\mathbf{k}}\gg 1$. For large $r_{k}$, the probability distribution covers
a large region in phase space where the expectation values of the even powers
of the fields are large compared to the expectation value of the commutator
$[\hat{q}_{\mathbf{k}},\hat{p}_{\mathbf{k}}]=i$. Therefore, all relevant
expectation values can be well approximated by a classical probability
distribution.

Importantly for the stochastic calculation, squeezing also makes the operators
$\hat{q}_{\mathbf{k}}$ and $\hat{p}_{\mathbf{k}}$ highly correlated, so the
field and momentum kicks become proportional to each other. This can be seen
from the classical Gaussian probability distribution that reproduces the
correlators (2.26)–(2.28), together with
$\langle\hat{q}_{\mathbf{k}}\rangle=\langle\hat{p}_{\mathbf{k}}\rangle=0$.
Dropping the subscript $\mathbf{k}$ to lighten the notation, the distribution
is

$P(q,p)=\frac{1}{2\pi\sqrt{\Sigma_{qq}\Sigma_{pp}-\Sigma_{qp}^{2}}}\exp[-\frac{\Sigma_{pp}q^{2}+\Sigma_{qq}p^{2}-2\Sigma_{qp}qp}{2(\Sigma_{qq}\Sigma_{pp}-\Sigma_{qp}^{2})}]\
.$ (2.29)

From (2.26)–(2.28), the squeezing limit $r_{k}\to\infty$ corresponds to
$\Sigma^{2}_{qp}\to\Sigma_{qq}\Sigma_{pp}$. In this limit, the distribution
becomes

$P(q,p)\to\frac{1}{\sqrt{2\pi\Sigma_{qq}}}e^{-\frac{q^{2}}{2\Sigma_{qq}}}\delta\quantity(p-\mathrm{sign}(\Sigma_{qp})\sqrt{\frac{\Sigma_{pp}}{\Sigma_{qq}}}q)\
.$ (2.30)

The variables $q$ and $p$ become tightly correlated, and the phase space
distribution is squeezed to the line
$p=\mathrm{sign}(\Sigma_{qp})\sqrt{\frac{\Sigma_{pp}}{\Sigma_{qq}}}q$. This
also applies to $\phi_{\mathbf{k}}$ and $\phi^{\prime}_{\mathbf{k}}$ and the
kicks, giving
$\xi_{\pi}=\xi_{\phi}\real\frac{\phi^{\prime}_{k}}{\phi_{k}}$.333To see this,
first note that in the squeezed limit
$\Sigma^{2}_{qp}\to\Sigma_{qq}\Sigma_{pp}$, (2.26)–(2.28) give
$|\real(\phi_{k}^{*}\phi_{k}^{\prime})|=|\phi_{k}^{*}\phi_{k}^{\prime}|$. This
is only possible if
$\real(\phi_{k}^{*}\phi_{k}^{\prime})=\phi_{k}^{*}\phi_{k}^{\prime}$, that is,
$\phi_{k}^{*}\phi_{k}^{\prime}$ is real. Then
$\mathrm{sign}(\Sigma_{qp})\sqrt{\frac{\Sigma_{pp}}{\Sigma_{qq}}}=\frac{\phi_{k}^{*}\phi_{k}^{\prime}}{|\phi_{k}^{*}\phi_{k}^{\prime}|}\frac{|\phi^{\prime}_{k}|}{|\phi_{k}|}=\frac{\phi^{\prime}_{k}}{\phi_{k}}=\real\frac{\phi^{\prime}_{k}}{\phi_{k}}\
.$ The last equality follows from the reality of
$\phi_{k}^{*}\phi_{k}^{\prime}$ by noticing that
$\phi^{\prime}_{k}/\phi_{k}=|\phi^{\prime}_{k}|^{2}/(\phi_{k}^{*}\phi_{k}^{\prime})^{*}\in\mathbb{R}$.
Since the state is not completely squeezed in our practical computations, we
take the real part to get rid of the small imaginary contribution. Therefore
we effectively have only one random variable. In SR, this variable is aligned
almost purely with the $\phi$ direction, as discussed around (2.18). Outside
of SR, the $\pi$ component can, in principle, become significant, though it
remains correlated with the $\phi$ component.

We use the value of $r_{k}$ as a measure of classicalisation, and demand that
modes satisfy $\cosh(2r_{k})>100$ before they are counted as part of the
background. Setting this as both the necessary and sufficient condition would
lead to a time-dependent coarse-graining parameter $\sigma$. We skip such a
complication, as our treatment of classicalisation is in any case simplified.
Instead, we choose a constant value of $\sigma$ small enough that the
condition $\cosh(2r_{k})>100$ is satisfied for all modes before they are
coarse-grained. With the exception of the Planck mass relic dark matter case,
the choice $\sigma=0.01$ is sufficient to guarantee this. Slightly larger
values of $\sigma$ could also be used, up to $\sigma=0.07$ for the asteroid
case, and $\sigma=0.03$ for the solar and the supermassive case. We also
consider smaller values and find that the results remain pretty much the same,
unless very small values of $\sigma$ are considered, so that the amplitude of
the kicks changes significantly between the Hubble exit and coarse-graining.
For the Planck mass relic case we have to take $\sigma=0.35$ for the modes
that produce PBHs to reach the classicalisation threshold. We discuss the
dependence on the value of $\sigma$, and the problems of the relic case, in
section 3.

A possible complication arises from the fact that $r_{k}$ is not always a
monotonic function of time. This is not a problem in SR, where modes only get
more squeezed with time. However, when SR is broken, squeezing can be reduced,
even completely undone. In figure 1 (left) we show $r_{k}$ for a single mode
as a function of time for the asteroid case. The mode exits the Hubble radius
at $N=30.1$, and inflation ends at $N=51.7$. The field is in USR from $N=30.9$
to $N=34.5$. Squeezing grows until a few e-folds after Hubble exit, and the
classicality threshold $\cosh(2r_{k})=100$ is reached at $N=31.2$. However,
after $N=31.4$, the squeezing amplitude $r_{k}$ starts to decrease, and the
mode goes under the threshold at $N=31.7$, before crossing back at $N=32.9$.
The situation does not change if we change the classicality threshold $r_{k}$,
as some other modes will cross back under the new threshold. To get insight
into the origin of this feature, we show in figure 1 (right) the normalised
field perturbation $\delta{\phi}_{k}/\bar{\phi}^{\prime}$, which, like the
power spectrum $\mathcal{P}_{\mathcal{R}}(k)$, freezes to a constant value
during SR on super-Hubble scales. The modes for which squeezing is partly
undone exit the Hubble radius just before or during the USR regime. This is
why the plot does not have a constant part at early times: the mode
transitions from the standard sub-Hubble time evolution directly to USR time
evolution. The plateau on the right is the super-Hubble SR regime after USR.
We see that the decrease in $r_{k}$ is related to (but does not perfectly
correlate with) a sharp dip of $\delta{\phi}_{k}/\bar{\phi}^{\prime}$ followed
by rapid growth. Such behaviour is generic in transitions between SR and USR
[137, 138, 139, 140, 141].

Figure 1: Squeezing parameter (left) and perturbation amplitude
$\delta\phi_{k}/\bar{\phi}^{\prime}$ (right) as a function of e-folds $N$ in
the asteroid case for a mode that exits the Hubble radius shortly before USR.
The dashed vertical line marks the Hubble exit, and the green area indicates
the USR region.

We simply assume that once modes cross our classicality threshold for the
first time, they irreversibly become part of the classical stochastic
background. In our simulations, 1% to 8% of the modes cross back and are
desqueezed, depending on the potential and stochastic realisation. These modes
exit the Hubble radius shortly before USR, and reach the coarse-graining scale
soon after the end of USR, when the background field is still near the feature
in the potential. We checked that disregarding this small range of modes does
not have a noticeable impact on the results. Note that in addition to
squeezing, classicalisation also proceeds through decoherence via
interactions, which there is no reason to expect to be undone by USR. For this
reason, treating classicalisation as irreversible seems justified. Ultimately,
this issue would have to be addressed by properly treating the division into
classical background and quantum perturbations from first principles in
quantum field theory.

###  2.3 Stochastic $\Delta N$ formalism and black hole abundance

Stochastic effects lead to an exponential tail [85, 89, 88] for the
probability distribution $P(\mathcal{R})$ of the comoving curvature
perturbation $\mathcal{R}$. It is therefore not enough to consider the power
spectrum $\mathcal{P}_{\mathcal{R}}(k)$ to calculate the PBH abundance. We
instead use the $\Delta N$ formalism [98, 99, 46, 100] to find the value of
$\mathcal{R}$ on a coarse-grained patch, given by

$\mathcal{R}=N-\bar{N}\equiv\Delta N\,.$ (2.31)

Here $N$ is the number of e-folds in the local patch from an initial
unperturbed hypersurface with constant initial field value
${\bar{\phi}}_{\mathrm{i}}$ to a final hypersurface with constant field value
${\bar{\phi}}_{\mathrm{f}}$ that corresponds roughly to the end of inflation,
and $\bar{N}$ is the mean number of e-folds averaged over stochastic
realisations (i.e. patches). Choosing a final hypersurface with constant field
value ${\bar{\phi}}_{\mathrm{f}}$ corresponds to calculating the curvature
perturbation in the comoving gauge, where field perturbations vanish.

The patch size is determined by the coarse-graining scale
$k_{\mathrm{c}}=\sigma aH$. We start the simulation when the wavenumber of the
mode corresponding to the CMB pivot scale $k_{*}$ equals $100aH$. We have
checked that starting the simulation earlier to accumulate kicks from a longer
period has negligible impact on our results, as the kicks during the SR stage
are small and far removed from the PBH scale. We tested this with 50, 100 and
500 e-folds of evolution before the CMB pivot scale. Our results do not
exclude the possibility that there are effects that grow logarithmically or
otherwise very slowly with $N$, in which case the history of kicks could be
important if inflation lasts very long before the scale $k_{*}$ exits the
Hubble radius. Also, if inflation lasts long, expectation values may not be
representative of the observations made in our universe over a finite volume
[142].

We stop the stochastic kicks when the patch reaches the scale
$k_{\mathrm{PBH}}$ that corresponds to the PBH mass we consider, given below
in (2.32). The reason is that perturbations with wavelengths smaller than the
size of the collapsing region do not affect PBH formation at this scale: they
behave as noise that averages out in the coarse-graining process. (For the
asteroid case we have explicitly checked that continuing the kicks until the
end of inflation and averaging over scales below the PBH scale gives the same
result as simply stopping the kicks when $k_{\mathrm{c}}=k_{\mathrm{PBH}}$.)
Subsequently, we evolve the equations without kicks until
${\bar{\phi}}_{\mathrm{f}}$ is reached. We record the final value of $N$, and
repeat the process many times to find the probability distribution $P(N)$, and
thus $P(\mathcal{R})$.

This approximation does not take into account correlations between different
PBHs. Two overdensities separated by distance $L$ have shared the same kicks
from modes with wavelengths larger than $L$, but have a different history of
kicks from shorter wavelength modes. See [87] for a discussion of how to
obtain the two-point function from the Langevin equations. We consider only
PBH abundance, i.e. the one-point function, not spatial correlations.

Sometimes the uniform density curvature perturbation $\zeta$ is used instead
of $\mathcal{R}$. In SR, these perturbations are equal up to a sign, but
during USR they evolve differently, see e.g. [143]. However, if we have a
period of SR after USR, it doesn’t make a difference which one we use. If USR
is not followed by SR, or the SR period is short (as in the case of relic
PBHs), $\zeta$ and $\mathcal{R}$ may not have had time to equalise, and one
should be careful regarding which curvature perturbation to use and how to
relate it to the density contrast that determines PBH formation.

When a sufficiently overdense perturbation of wavenumber $k$ re-enters the
Hubble radius during the radiation-dominated phase after inflation, almost the
entire Hubble radius collapses into a PBH with mass [25, 144]

$\displaystyle M=\frac{4}{3}\pi\gamma H_{k}^{-3}\rho_{k}\approx 5.6\times
10^{15}\gamma\left(\frac{k}{k_{*}}\right)^{-2}M_{\odot}\,,$ (2.32)

where $M_{\odot}\approx 2\times 10^{33}$ g, $\gamma\approx 0.2$ is a parameter
that characterises the collapse [144], and $H_{k}$ and $\rho_{k}$ are the
background Hubble rate and the energy density at Hubble entry, respectively.
We assume standard radiation-dominated expansion history with the Standard
Model degrees of freedom.

There are various analyses of and refinements to the treatment of the collapse
[144, 145, 146, 147, 148, 149, 22, 150, 151, 152, 153, 154, 155, 156, 157,
158, 159, 160, 161, 162, 163, 164, 165]. We simply assume that a region
collapses if the comoving curvature perturbation exceeds a critical value
$\mathcal{R}_{c}$. We adopt the value $\mathcal{R}_{c}=1$, close to the most
recent estimates. The abundance is extremely sensitive to this choice.
Increasing or decreasing $\mathcal{R}_{c}$ by 10% changes the abundance by
$\sim$ 1–2 orders of magnitude in the asteroid, solar and supermassive cases;
the relic case is more stable. Decreasing to the original estimate
$\mathcal{R}_{c}=3/4$ from [144]444This corresponds to the critical density
perturbation $\delta_{c}=1/3$, equal to the equation of state parameter $w$ in
radiation domination, related to the comoving curvature perturbation by
$\delta=-\frac{2(1+w)}{5+3w}\nabla^{2}\mathcal{R}/(aH)^{2}$, or
$\mathcal{R}_{c}=\frac{9}{4}\delta_{c}$ at Hubble crossing (see e.g. [166]).
would increase the abundance by 3, 3, and 5 orders of magnitude in the
asteroid, solar and supermassive cases, respectively, and by less than an
order of magnitude in the relic case. This issue concerns the correct
modelling of the PBH formation process, given the probability distribution
$P(\mathcal{R})$. Therefore, we do not consider it further, since our aim is
to determine the shape of $P(\mathcal{R})$, and see how it affects the
abundance, given a procedure to calculate the abundance. Furthermore, a tiny
change in the parameters of the inflationary potential leads to a large
difference in the PBH abundance, which can compensate for our simplified
modelling of the collapse.

The fraction of simulations where $\mathcal{R}>\mathcal{R}_{c}$ gives the
initial fraction $\beta$ of the energy density in PBHs,

$\displaystyle\beta=2\int_{\Delta N>\mathcal{R}_{c}}^{\infty}\mathrm{d}NP(N)\
,$ (2.33)

where the factor of 2 comes from the Press–Schechter prescription to solve the
cloud-in-a-cloud problem [167]. Since the energy density of PBHs scales as
non-relativistic matter, their fraction of the energy density grows during
radiation domination, and today it is [25]

$\Omega_{\text{PBH}}\approx 9\times
10^{7}\gamma^{\frac{1}{2}}\beta\left(\frac{M}{M_{\odot}}\right)^{-\frac{1}{2}}\
.$ (2.34)

In order not to overproduce dark matter, the initial abundance $\beta$ has to
be small, as the collapse happens deep in the radiation-dominated era, and
hence the energy density fraction has a lot of time to increase. In the cases
we consider below, $\beta$ varies between $10^{-17}$ and $10^{-8}$.

If the distribution of $\mathcal{R}_{k}$ were Gaussian [168, 149], as is often
assumed in the literature, its variance would be
$\sigma_{\mathcal{R}}^{2}=\int_{k_{\mathrm{IR}}}^{k_{\mathrm{PBH}}}\mathrm{d}(\ln
k)\mathcal{P}_{\mathcal{R}}(k)$, where $\mathcal{P}_{\mathcal{R}}$ is the
curvature power spectrum and ${k_{\mathrm{IR}}}$ is a cutoff corresponding
roughly to the present Hubble radius. (The precise value of $k_{\mathrm{IR}}$
makes no difference.) The initial energy density fraction is then

$\beta=2\int_{\mathcal{R}_{c}}^{\infty}\mathrm{d}\mathcal{R}\frac{1}{\sqrt{2\pi}\sigma_{\mathcal{R}}}e^{-\frac{\mathcal{R}^{2}}{2\sigma_{\mathcal{R}}^{2}}}\approx\frac{\sqrt{2}\sigma_{\mathcal{R}}}{\sqrt{\pi}\mathcal{R}_{c}}e^{-\frac{\mathcal{R}_{c}^{2}}{2\sigma_{\mathcal{R}}^{2}}}\
.$ (2.35)

We build our inflationary potentials so that $\Omega_{\text{PBH}}$ calculated
from the power spectrum in the Gaussian approximation matches the observed
dark matter abundance today. We will see that the Gaussian approximation
strongly underestimates the real abundance.

##  3 Results

###  3.1 Analytical toy potentials

The details of how we discretise the continuum equations and solve them
numerically are given in appendix B. We first test our procedure on two toy
model potentials studied in [89], where $P(N)$ was calculated analytically. We
consider the family of cubic potentials

$\displaystyle V(\phi)$ $\displaystyle=$ $\displaystyle
V_{0}[1+\alpha(\phi-\phi_{0})+\beta(\phi-\phi_{0})^{3}]\ ,$ (3.1)

where $V_{0}>0$, $\alpha\geq 0$, $\beta\geq 0$, and $\phi_{0}$ are constants.
The potential is called flat or tilted if $\alpha=0$ or $\alpha\neq 0$,
respectively.

We use the simplified treatment of noise, where the momentum kicks vanish and
the field kick amplitude is given by (2.18), for two reasons. First, we are
interested in checking whether we can replicate the results of [89], which
used the simplified treatment. Second, our full treatment cannot start
directly in the flat section of the potential, as modes need to be first
evolved before any kicks can be given. The full treatment also requires the
potentials to be classically traversable, which is not the case for some of
the parameter values considered here.

In the flat cubic potential case we set $V_{0}=\frac{6\pi^{2}}{25}$,
$\alpha=0$, and $\beta=0.01$, and vary the starting point of the field so that
$(\phi-\phi_{0})/\Delta\phi_{\text{well}}$ takes the values $-0.5$, $-0.25$,
$0$, $0.25$, $0.5$, where
$\phi_{0}=\Delta\phi_{\text{well}}/2+0.3/\sqrt{\beta}$ and
$\Delta\phi_{\text{well}}=2\left(150\beta\right)^{-1/3}$. The results are
shown in figure 2 (left). In all cases, the exponential fit $P(N)=Ae^{-BN}$ to
the tail of the distribution gives robustly the slope $B=7.4\times 10^{-3}$.
Our results are to be compared to those shown in figure 9 of [89]. The
analytical leading term flat approximation in [89] gives the slope $B=8\times
10^{-3}$, while including higher order terms gives $B=7.4\times 10^{-3}$, in
perfect agreement with our result555The analytical result for the slopes is
only given to leading order in [89], but the numbers including higher order
corrections have been privately communicated to us by the authors..

In the tilted case we set $V_{0}=\frac{3\pi^{2}}{25}$, $\phi_{0}=0.1/\alpha$,
$\beta=10^{-3}$, and $\alpha$ takes the values $0.02$, $0.02666$, $0.03333$
and $0.04$. The results are shown in figure 2 (right). We again fit an
exponential and find that the slope $B$ takes the values $0.026$, $0.042$,
$0.064$ and $0.087$, respectively. These results are to be compared to those
in figure 11 of [89]. The leading term flat approximation in [89] gives the
corresponding slopes $0.022$, $0.037$, $0.057$ and $0.081$, while including
higher order terms5 gives values $0.025$, $0.042$, $0.062$ and $0.086$,
respectively, close to our results, especially given the errors shown in
figure 2.

These results show that our algorithm works as intended, and we confirm the
existence of the exponential tails. We now move to our realistic potentials.

Figure 2: Our results for the probability distribution $P(N)$ for the
potentials used in [89]. For the flat cubic potential (left), the different
colours correspond to different values of
$(\phi-\phi_{0})/\Delta\phi_{\text{well}}$. For the tilted cubic potential
(right), the different colours correspond to different values of the parameter
$\alpha$. The black dashed lines are exponential fits to the numerical
results.

###  3.2 Realistic potentials

####  3.2.1 Potentials for the four cases

The cubic models considered above describe potentials only around the region
where the PBH seeds are generated. In a realistic model, we have to fit not
only the PBH adundance, but also CMB observables, which correspond to larger
inflaton field values. When building the models, the location of the USR
region is determined by the relation (2.32) between the PBH mass and
wavenumber. For longer wavelengths, the Hubble volume at collapse is bigger,
so the PBH mass is larger. Hence, the larger the PBH mass, the closer the USR
region is to the CMB scales, and correspondingly the more difficult it is to
have a potential that leads to USR and produces enough PBHs without spoiling
the CMB predictions. Although the difference in the wavenumber of modes that
collapse to form PBHs and modes in the observed CMB region can be
exponentially large, the difference in field values is moderate. For the
potentials we consider, the field value at the CMB pivot scale is between 3 to
6 times its value at the beginning of USR, which we define as the period when
the second SR parameter $\eta=-\ddot{\bar{\phi}}/(H\dot{\bar{\phi}})$ is
larger than 1.

We look at four cases: PBHs with asteroid mass, solar mass, and Planck mass,
which can be dark matter, as well as PBHs with mass $1800M_{\odot}$, which can
act as seeds for supermassive black holes [39, 40, 41]. Our starting point for
the inflaton potentials is non-minimally coupled Standard Model Higgs
inflation [101, 102], and we use loop corrections in the chiral Standard Model
to create a shallow minimum that leads to USR [25] (see also [18, 24, 19, 24,
103, 20]). This model can fit the PBH abundance and the CMB observables
simultaneously only if the PBHs are Planck mass relics. In the other cases, we
adjust the potentials by hand on the CMB scales, like we did in $Paper\,I$. In
practice, this means changing numerical factors in the potential above a
threshold that lies between the PBH and CMB scales, as described in appendix
A. The resulting potentials for the four cases are shown in figure 3, and the
masses, CMB observables, and PBH parameters are listed in table 1.

Figure 3: The potentials for the asteroid, solar, supermassive, and relic
case. The dashed lines mark the initial field value $\bar{\phi}_{i}$ at the
CMB scale, the final field value $\bar{\phi}_{f}$ near the end of inflation,
and the PBH scale $\bar{\phi}_{\text{PBH}}$ corresponding to the Hubble exit
of the last mode that gives a kick. The red sections indicate the USR region
(small $\phi$) and the kink that interpolates between the low and high
sections of the potential (large $\phi$).

We use the metric formulation of general relativity. In [25], PBHs from USR
Higgs inflation were considered also in the Palatini formulation. This makes
the potential much flatter, and as a consequence the PBH abundance is much
more sensitive to the input parameters, and it is difficult to find the
parameters for which the Gaussian approximation gives even roughly the
observed PBH abundance.

| Asteroid | Solar | Supermassive | Relic
---|---|---|---|---
$M$ | $7.2\times 10^{-15}M_{\odot}$ | 4.7 $M_{\odot}$ | $1.8\times 10^{3}M_{\odot}$ |  $1.4\times 10^{3}$ kg
$n_{s}$ | 0.966 | 0.965 | 0.962 | 0.958
$r$ | 0.01 | 0.05 | $8\times 10^{-3}$ | $5\times 10^{-3}$
$\alpha_{s}$ | $-3\times 10^{-4}$ | $7\times 10^{-4}$ | $-7\times 10^{-4}$ | $-9\times 10^{-4}$
Gaussian |  |  |  |
$\sigma_{\mathcal{R}}^{\text{G}}$ | 0.122 | 0.1739 | 0.118 | 0.149
$\Omega_{\text{PBH}}^{\text{G}}$ | 0.13 | 0.17 | $1.4\times 10^{-5}$ | 0.11
Full |  |  |  |
$\sigma_{\mathcal{R}}$ | 0.123 | 0.1744 | 0.119 | 0.184
$A$ | $1699\pm 61$ | $1403\pm 26$ | $2136\pm 27$ | $225\pm 5$
$B$ | $32.7\pm 1.2$ | $26.7\pm 0.5$ | $41.3\pm 0.5$ | $4.38\pm 0.09$
$\Omega_{\text{PBH}}^{\text{fit}}$ | $(1.6\pm 0.4)\times 10^{4}$ | $1.58\pm 0.06$ | $0.049\pm 0.01$ | $(2.38\pm 0.06)\times 10^{7}$
$\Omega_{\text{PBH}}^{\text{num}}$ | $-$ | 1.55 | $-$ | $2.39\times 10^{7}$
$\Omega_{\text{PBH}}^{\text{fit}}/\Omega_{\text{PBH}}^{\text{G}}$ | $10^{5}$ | $10$ | $10^{3}$ | $10^{8}$
Simplified |  |  |  |
$A$ | $2080\pm 63$ | $998\pm 4$ | $1833\pm 58$ | $3995\pm 57$
$B$ | $40.0\pm 1.2$ | $19.0\pm 0.1$ | $35.4\pm 1.1$ | $77.8\pm 1.1$
$\Omega_{\text{PBH}}$ | $26\pm 12$ | $(1.25\pm 0.01)\times 10^{2}$ | $17\pm 5$ | $(5.5\pm 4.9)\times 10^{-24}$
$\Omega_{\text{PBH}}/\Omega_{\text{PBH}}^{\text{G}}$ | $10^{2}$ | $10^{3}$ | $10^{6}$ | $10^{-23}$
Table 1: Masses, CMB observables and PBH parameters for the four potentials.
Here $M$ is the initial PBH mass, $n_{s}$ is the CMB scalar spectral index,
$r$ is the tensor-to-scalar ratio, and $\alpha_{s}$ is the running of the
scalar spectral index. In the rows labeled “Gaussian”,
$\sigma_{\mathcal{R}}^{\text{G}}$ is the standard deviation of the
perturbations calculated from the power spectrum as explained above (2.35),
and $\Omega_{\text{PBH}}^{\text{G}}$ is the corresponding PBH abundance today.
In the rows labeled “Full”, $\sigma_{\mathcal{R}}$ is the standard deviation
of the Gaussian fit to the full $P(N)$, $A$ and $B$ characterise the
exponential fit $P(N)=e^{A-BN}$ to the tail,
$\Omega_{\text{PBH}}^{\text{fit}}$ is the PBH density fraction today
calculated from the exponential fit, and $\Omega_{\text{PBH}}^{\text{num}}$ is
the density fraction calculated directly from $P(N)$ by numerical integration.
The rows labeled “Simplified” give the corresponding numbers for the
simplified treatment of kicks (2.18). For relics, the loss of mass through
evaporation has been taken into account in $\Omega_{\text{PBH}}$. The errors
are based on a jackknife analysis with 20 subsamples.

####  3.2.2 Asteroid mass dark matter

In the mass window $\sim 10^{-15}$–$10^{-11}M_{\odot}$ PBHs can constitute all
of the dark matter [39, 40, 41]. We consider a value close to the lower end of
this range, $7.2\times 10^{-15}M_{\odot}$, around the mass of the asteroid
Eros. As discussed above, we adjust the renormalisation group improved Higgs
potential by hand in the CMB region to make the spectral index less red, to
fit observations. The tensor-to-scalar ratio is compatible with observations
even without the adjustment.

We have run $1024\times 10^{8}$ simulations of the inflationary dynamics, each
representing a realisation of the solution to (the discretised version of) the
set of equations (2.7), (2.8), (2.9), and (2.20). The main results for this
potential have already been reported in $Paper\,I$. The probability
distribution $P(N)$ is shown in figure 4, and some relevant numbers are given
in table 1. For $|\Delta N|\gtrsim 0.5$ the shape deviates significantly from
a Gaussian, with an exponential tail. The distribution is clearly skewed
towards large $\Delta N$. This can be easily understood. Each kick moves the
field forward or back with equal probability, but in USR the field velocity
decays rapidly. Hence, the longer the field spends in the USR region, the
slower it moves, so being kicked back increases $\Delta N$ more than being
kicked forward decreases it. The effect is sizeable only in the tail of the
distribution; the mean $\bar{N}$ remains essentially unchanged. We have data
points up to $\Delta N=1.09$, but the region $\Delta N\gtrsim 0.95$ is poorly
resolved. The computational cost of resolving further into the tail would be
prohibitive given our computational resources. We fit an exponential
$e^{A-BN}$ to the tail in the range $\Delta N=0.75\ldots 0.95$, where the
statistics are robust. We then calculate the abundance from (2.33) using an
extrapolation of the fit. As shown in table 1, the PBH abundance is enhanced
by a factor of $\sim 10^{5}$ compared to the Gaussian case.

To assess the importance of mode evolution in a background affected by
previous kicks, we have done $256\times 10^{8}$ runs with the simplified
treatment where the momentum kicks vanish and the amplitude of the field kicks
is given by (2.18), as discussed in section 2.1.3. This leads to a much
steeper slope, as shown in figure 4. In this case, we only have points up to
$\Delta N=0.79$. Fitting an exponential to the well-resolved range $\Delta
N=0.6\ldots 0.7$ and extrapolating to $\Delta N\geq 1$ gives
$\Omega_{\text{PBH}}=26$. This is $\sim 10^{2}$ times higher than in the
Gaussian case, but a factor $\sim 10^{3}$ smaller than in the full treatment.
The result is subject to extrapolation errors, but it is clear that the
simplified treatment significantly underestimates the stochastic enhancement
of the PBH abundance.

We also considered a recursive treatment where we first evolve the mode
functions on the FLRW background without noise, and then do a run where these
modes are used to calculate the noise amplitude (2.15). The results are
(within the statistical error) indistinguishable from those of the full
treatment, unless $\sigma$ is so small that the kicks are pushed close to the
end of inflation. While the simplified treatment above shows that it is
crucial to take into account mode evolution when calculating the noise, the
recursive treatment establishes that the change of the coarse-grained
background due to noise has little impact on the mode evolution (except for
small $\sigma$). This is perhaps not surprising, given that the Hubble rate,
which characterises the background evolution, is almost constant (except
during the transition from SR to USR), so the kicks don’t actually perturb the
background much. Thus, while the full treatment is not Markovian, and hence
does not reduce to a Fokker–Planck equation in principle, the non-Markovian
nature of the dynamics turns out to be irrelevant. Furthermore, we have
checked that the difference in $\Omega_{\text{PBH}}$ is only $21\%$ (i.e.
smaller than our 25% jackknife error) if we fix the amplitude of the kicks to
the square root of the variance (which evolves in time), and randomise only
their sign, i.e. whether the field gets kicked forward or backward. This
approximation may simplify analytical calculations of the stochastic effects.

Figure 4: The probability distribution $P(N)$ for the asteroid case. The
bottom label gives the number of e-folds at the end of the simulation, the top
label gives the deviation from the mean. The blue solid line is the full
numerical result. The black dotted line is a Gaussian fit, and the black
dashed line is the exponential fit to the tail. The red dash-dotted line is
the simplified treatment. The vertical line marks the collapse threshold
$\Delta N=1$.

In figure 5, we show the kicks in the full and in the simplified treatment for
a typical run. In the simplified treatment, the typical kick amplitude is
almost constant. In contrast, in the full treatment the kick amplitude varies
by orders of magnitude, with a distinct dip during USR. This structure traces
back to the curvature power spectrum, shown in figure 6, which drops for modes
that exit the Hubble radius just before the start of USR due to a cancellation
between terms in the transition region [137, 138, 139, 140, 141], and peaks
for modes that exit the Hubble radius during USR. Note the offset between the
two figures. The strongest perturbations originate from USR, but give their
kicks later, when they exit the coarse-graining scale, by which time USR is
already over. After the end of USR, the kick amplitude continues to rise, and
the variance is a maximum of 80 times larger in the full treatment than in the
simplified treatment. This should be compared to the classical drift of the
field, which dips during USR when the field slows down, and grows again
afterwards. In practice, the kicks have the strongest effect somewhat after
the end of USR, when the kicks have already grown to large values but the
classical drift is still small. Individual kicks are much smaller than the
classical contribution, which may in part explain why fixing the kick
amplitude to the square root of the variance gives results indistinguishable
from the full treatment: many kicks are required to have an impact.

Decreasing the coarse-graining scale parameter $\sigma$ from $10^{-2}$ to
$10^{-4}$ changes $\Omega_{\text{PBH}}$ by only 4%, much less than our
statistical error of 25%. For $\sigma=10^{-5}$ the abundance grows by a factor
of 3, and for $10^{-6}$ it grows by a factor of 5. Making $\sigma$ smaller
increases the offset between when the mode exits the Hubble radius and when it
delivers its kick. As long as the kick amplitude remains similar, the ensuing
difference in the results is small, but it grows as the kicks are moved closer
to the end of inflation. The agreement between the full and recursive
treatment also degrades in this case. The difference between the full and
recursive treatment grows to an order of magnitude for $\sigma=10^{-7}$. For
$\sigma=10^{-8}$, the kicks are pushed to the very end of inflation, and their
amplitude becomes large. The $P(N)$ distribution in the full treatment changes
drastically, becoming similar to the relic case discussed below in section
3.2.5. We then also encounter the issue that some of the simulations never
finish (as also happens in the relic case). A lower limit for $\sigma$ has
been given from various requirements, including that the noise is given by the
simplified form (2.18) [58, 69], a derivation of the stochastic formalism from
quantum field theory [82, 104], SR remaining an attractor [83], and the
validity of neglecting heavier fields in the multifield case [86]. The
requirement in [83] that the slope of the potential should not change
significantly between the time a mode exits the Hubble radius and gives its
kick agrees with the condition we find for the results to be insensitive to
the value of $\sigma$.

Figure 5: Stochastic kicks as a function of e-folds for the asteroid case. The
black solid line is the classical evolution of the field, the kicks in the
full treatment are in blue, and the kicks in the simplified treatment are in
red. The classical evolution has been scaled by $\differential{N}^{-1}$ and
the kicks by $\differential{N}^{-1/2}$ to make the overall amplitudes in the
plot independent of the time step $\differential{N}$. (See appendix B for
details.) The green area denotes the USR region. The plot shows all the kicks
in a typical run. Figure 6: The power spectrum $\mathcal{P}_{\mathcal{R}}(k)$
for the asteroid case over different scales, evaluated at the end of inflation
when all the super-Hubble scales have frozen. The e-fold number $N$
corresponds to the time the mode exits the Hubble radius. The dashed lines
mark the CMB scale and the end of inflation, and the green area indicates USR.
The last scale that contributes to the PBHs, corresponding to the end of USR,
is also marked.

####  3.2.3 Solar mass dark matter

Observations seem to imply that PBHs around the solar mass cannot form all of
the dark matter, although some of the constraints may be ameliorated by
nonstandard clustering [39, 40, 41]. In any case, such PBHs could contribute a
subdominant component of dark matter. We consider PBHs with mass
$4.7M_{\odot}$. As in the asteroid case, the USR region that seeds PBHs is
produced by renormalisation group running, and we adjust the CMB part of the
potential by hand to make the spectral index less red to fit observations. The
tensor-to-scalar ratio fits the CMB data available up to Sept. 2021 [1] with
or without the adjustment.666As this paper was being finalised, a new combined
analysis of new BICEP/Keck, Planck and WMAP CMB data appeared [169]. The
updated constraint $r<0.036$ excludes our solar mass potential, which has
$r=0.05$. This could be fixed by adjusting the height of the potential in the
same way as for the supermassive case discussed in section 3.2.4. However,
this would require re-running all simulations for the solar mass case, which
would take a lot of computational resources without changing the qualitative
conclusions.

Figure 7: The probability distribution $P(N)$ for the solar case, similarly to
figure 4. Figure 8: Stochastic kicks for solar case, similarly to figure 5.
Figure 9: The power spectrum $\mathcal{P}_{\mathcal{R}}(k)$ for the solar
case, similarly to figure 6.

As the PBH mass is a factor of $\sim 10^{15}$ higher than in the asteroid
case, fewer PBHs are needed today to match the observed dark matter abundance.
However, as heavier PBHs collapse later, they have less time to increase their
fraction of the energy density. As a result, the initial PBH abundance has to
be larger than in the asteroid case, as quantified in (2.34). Therefore, the
variance from the power spectrum is larger, as shown in table 1. The slope of
the full numerical $P(N)$ distribution shown in figure 7 also turns out to be
less steep than in the asteroid case. With $256\times 10^{8}$ simulations, the
distribution reaches up to $\Delta N=1.20$. The mean $\bar{N}$ is again
unchanged from the Gaussian result. We can now numerically integrate the PBH
abundance, with no need to extrapolate to reach the collapse threshold. The
abundance is enhanced by one order of magnitude compared to the Gaussian case.
The tail is well fit by an exponential, which gives the same result for the
abundance as numerical integration of the distribution.

In the simplified treatment, with $256\times 10^{8}$ simulations we reach up
to $\Delta N=1.45$. In contrast to the asteroid case, now the simplified
treatment gives a larger enhancement than the full treatment, by two orders of
magnitude. The enhancement over the Gaussian case is three orders of
magnitude. In figure 8, we show the kicks for a typical run in the full and
simplified treatments. The maximum amplitude of the kicks after USR in the
full treatment is a factor of 5 larger than in the simplified treatment; the
difference is much smaller than in the asteroid case. This enhancement after
USR in the full treatment is not sufficient to overcome the dip during USR,
compared to the steady kicks in the simplified case, which explains why the
simplified case gives overall a larger enhancement. As in the asteroid case,
the recursive treatment gives results indistinguishable from the full
treatment after one iteration.

The calculation is less sensitive to the choice of the coarse-graining scale
than in the asteroid case. Decreasing $\sigma$ from $10^{-2}$ to $10^{-8}$
changes $\Omega_{\text{PBH}}$ by less than 30%. The effect is large only when
the kicks are pushed to the very end of inflation, which corresponds to the
value $\sigma=10^{-16}$.

As can be seen from figure 9, the USR region is now closer to the CMB region
than in the asteroid case: because the PBHs are more massive, their seeds have
longer wavelengths and exit the Hubble radius earlier. In the asteroid case,
USR begins 31 e-folds after the CMB pivot scale exits the Hubble radius, but
in the solar case this interval is 14 e-folds. In terms of field values, the
difference is less drastic: the value of the field changes by
$4M_{{}_{\mathrm{Pl}}}$ between the CMB pivot scale exit and the start of USR
in the asteroid case and by $3M_{{}_{\mathrm{Pl}}}$ in the solar case. The
larger the mass, the smaller the separation, and the more difficult it is to
keep the CMB perturbations unaffected by the USR region. This is even clearer
in the case of supermassive PBHs.

####  3.2.4 Supermassive black hole seeds

Apart from dark matter, PBHs have been suggested as seeds for supermassive
black holes observed at high redshifts. It is not clear how these objects can
form sufficiently rapidly in the early universe [39, 41]. There may be an
astrophysical explanation, but initial PBHs that jumpstart the accretion
process are also a possible solution. We consider initial seeds with mass
$1.8\times 10^{3}M_{\odot}$, which then grow by accretion to
$10^{9}M_{\odot}$, with abundance today $\Omega_{\text{PBH}}=1.4\times
10^{-5}$ in the Gaussian approximation. USR begins 10 e-folds after the exit
of the CMB pivot scale, which is only a difference of $2M_{{}_{\mathrm{Pl}}}$
in the field value. In addition to adjusting the potential in the CMB region
by hand, as in the previous cases, we now have to also adjust the overall
height of the whole potential. Otherwise, the tensor-to-scalar ratio $r$ grows
with the PBH mass as the USR and CMB regions come closer, and exceeds the
observational CMB upper limit $0.036$ [169], or even the previous one
$r<0.061$ [1].777In the Palatini formulation of general relativity mentioned
above, $r$ is suppressed [170, 25], so there would be no such problem.
However, the spectral index would still need to be adjusted by hand.

Figure 10: The probability distribution $P(N)$ for the supermassive case,
similarly to figure 4. Figure 11: Stochastic kicks as a function of e-folds
for the supermassive case, similarly to figure 5. Figure 12: The power
spectrum $\mathcal{P}_{\mathcal{R}}(k)$ for the supermassive case, similarly
to figure 6.

With $1024\times 10^{8}$ simulations, the $P(N)$ distribution for this case
reaches up to $\Delta N=0.83$, and is shown in figure 10. The mean $\bar{N}$
is again unchanged from the Gaussian result. Fitting an exponential to the
well-resolved range $0.6\ldots 0.75$ and extrapolating beyond $\Delta N=1$
gives an enhancement of $10^{3}$ over the Gaussian case. In the simplified
case, $256\times 10^{8}$ simulations reach up to $\Delta N=0.93$, leading to
an extrapolated enhancement of $\sim 10^{6}$ over the Gaussian case. As in the
solar case, the simplified treatment overestimates the enhancement of the PBH
abundance. We show typical kicks in figure 11. As in the solar case, the
maximum amplitude of the kicks after USR is only a factor $5$ larger in the
full treatment compared to the simplified case. Again, the recursive treatment
gives results identical to the full treatment with one iteration.

As in the previous cases, the results depend very weakly on the value of
$\sigma$: decreasing it from $10^{-2}$ to $10^{-8}$ changes
$\Omega_{\text{PBH}}$ only by 70%. Significant change is seen only at
$\sigma=10^{-17}$, when the kicks are delayed until the end of inflation.

In this model, the peak of the power spectrum (shown in figure 12),
$\mathcal{P}_{\mathcal{R}}=0.005$, is located at $k=5400$ Mpc-1, close enough
to the CMB scale that it can be probed by future spectral distortion
observations [171, 172]. The model may already be in tension with the FIRAS
measurements [173, 174, 137, 40], although a full spectral distortion analysis
would be needed to draw definite conclusions. This is an important
consideration for all USR models that generate PBHs in this mass range.

####  3.2.5 Planck mass relic dark matter

There are essentially no observational constraints (apart from the dark matter
relic abundance) on the scenario where dark matter consists of uncharged
Planck mass objects (see [175] for a proposal of how to observe them). We
consider PBHs that form with mass $1.4\times 10^{3}$ kg and evaporate down to
the Planck mass $M_{{}_{\mathrm{Pl}}}=4.3\mu$g. The usual semiclassical
computation of Hawking evaporation breaks down close to the Planck mass, and
it is possible that evaporation stops and leaves a stable relic [176, 177,
178, 179, 180, 181, 182, 183, 184]. In the relic case, the power spectrum
peaks on scales that leave the Hubble radius a little after the end of USR.
For this reason, we have moved the PBH scale forward from the end of USR to
capture the modes with the largest amplitude, as shown in figure 15. Unlike in
the previous cases, there is no need to adjust the potential by hand. The non-
minimally coupled Standard Model Higgs potential with renormalisation group
corrections can produce both PBHs with the observed dark matter abundance (in
the Gaussian approximation) and the observed CMB perturbations [25]. As in the
other cases, the feature that leads to USR is extremely sensitive to the
parameters entering the renormalisation group running.

In the potentials shown in figure 3, a kink in the non-relic cases marks where
the adjusted inflationary plateau is matched to the part of the potential that
has the USR feature, whereas in the relic case, the potential continues
smoothly from the CMB scale down to the USR region. In all four mass cases the
USR region (unlike the CMB region) is roughly at the same field value:
$\bar{\phi}$ at the start of USR varies from $1.06M_{{}_{\mathrm{Pl}}}$ in the
supermassive case to $0.89M_{{}_{\mathrm{Pl}}}$ in the relic case, although
the difference in e-folds (and consequently PBH mass) is large. In the relic
case, USR ends only about 2 e-folds before the end of inflation. This is why
no adjustment is needed: the quantum corrections modify the successful Higgs
inflation scenario only at the very end of inflation. However, this leads to
various problems.

In the relic case, most modes that exit during USR do not have enough time to
reach the coarse-graining scale before the end of inflation, unless we
increase the value of $\sigma$ considerably compared to the other cases we
have explored. The lack of a long SR period after USR also means that the
perturbations $\mathcal{R}$ and $\zeta$ do not necessarily have time to
equalise, which should be taken into account when relating the inflationary
curvature perturbation to the post-inflationary density contrast. Furthermore,
as the modes exit close to the end of inflation, they quickly re-enter the
Hubble radius. Therefore PBHs may collapse during preheating instead of the
radiation-dominated era as we have assumed. This might affect the collapse
process and the evolution of their density fraction after formation, depending
on the equation of state. The duration of preheating is also crucial here, and
depends on the detailed shape of the potential and the couplings of the fields
(including between the different components of the Standard Model Higgs
doublet). In tree-level Higgs inflation (in the metric formulation), it is not
clear whether preheating takes few e-folds or is almost instant [185, 186,
187, 188, 189, 190, 191, 192, 193]. In the Palatini formulation, the field re-
enters the very flat inflationary plateau during preheating, leading to new
periods of inflation, and preheating takes less than one e-fold [194, 195]. In
general, if the inflaton remains close to its potential minimum during
preheating and the potential is monomial, $V(\phi)\propto|\phi|^{p}$ (with
$p\geq 2$), preheating is relatively rapid, taking one to a few e-folds [196],
depending on the couplings. It therefore depends on the details of preheating
whether the universe has had time to transition to radiation domination before
the PBHs form.

Another problem is that in the transition from SR to USR the field velocity
grows so much that inflation ends for 0.3 e-folds. (This problem seems to be
specific to the potential; for example, with the Palatini formulation
discussed above, similar relic potentials do not have this problem.) Then some
modes that have become part of the background start crossing back inside the
coarse-graining radius, and when inflation resumes, they exit for the second
time. Our heuristic treatment of the stochastic formalism does not prescribe
how these modes should be treated. Should they contribute stochastic kicks
when they first exit the coarse-graining radius, and again at second exit?
Such questions go beyond the simple system-environment split we use, and
should be resolved by a proper derivation from quantum field theory. We opt to
skip the second kicks, so we have a kick-free period until all modes that re-
entered the coarse-graining radius have re-exited, after which the kicks
resume. While this raises questions of principle, in practice these kicks are
so subdominant that it makes little quantitative difference how we treat them.
As shown in figure 14, the largest kicks during this period are roughly four
orders of magnitude smaller than the classical evolution. Towards the end of
inflation, the classical contribution increases as the background velocity
$\bar{\pi}$ starts to rise. The kicks also grow proportionally to $\bar{\pi}$,
so their impact relative to the classical contribution does not increase
towards the end. Similar but less drastic growth of the background and kicks
is also seen in the other mass cases.

In order for the modes relevant to Planck mass relics to have enough time to
give their kicks before the end of inflation, we set $\sigma=0.35$, so that
the delay between Hubble exit and giving kicks is short. The side effect of
choosing such a large value of $\sigma$ is that while the modes around the USR
region do classicalise in time according to our criterion,
$\cosh(2r_{k})>100$, the vast majority of modes before the USR region do not.
Therefore, we miss out on about 42 e-folds worth of kicks starting from the
beginning of the simulation until the first classicalised mode crosses the
coarse-graining scale. Using such a large value of $\sigma$ can be naively
expected to significantly overestimate the stochastic kicks, and not to give a
sufficient decoupling between the system and the environment. In any case, as
we have discussed, the treatment of stochastic effects in USR should be put on
more firm footing before any reliable conclusions can be drawn.

Running $256\times 10^{5}$ simulations, the $P(N)$ distribution goes up to
$\Delta N=3.77$. Figure 13 shows that the changes due to the stochastic
effects are qualitatively similar to the previous three cases, but more
extreme. The kicks now distort the distribution everywhere, not only in the
tail. The tail is still well described by an exponential. The mean value is
barely affected: classical FLRW evolution has $\bar{N}=51.31$, and the full
stochastic treatment gives $\bar{N}=51.34$. However, the distribution around
the mean is no longer described by a Gaussian. As in the other cases, the
first iteration of the recursive treatment gives results indistinguishable
from the full calculation.

Numerical integration gives a PBH abundance that is enhanced by a factor of
$10^{8}$ over the Gaussian calculation. The exponential fit gives results in
agreement. However, the wide $P(N)$ distribution leads to serious
interpretation problems. In approximately 6% of the runs inflation does not
end by $N=70$, when we stop the simulation, if the field has not reached the
reached the end-of-simulation value $\bar{\phi}_{\text{f}}$. Instead, the
field gets stuck in the local minimum. The corresponding points would be far
to the right of the points plotted in figure 13. As we do not know how long
the runs would continue, the shape of this far tail is unknown. The estimate
of PBH abundance is thus not reliable, as the number of these unfinished runs
is larger than the number of finished runs with $\Delta N\geq 1$. We could
also have the problem of eternal inflation, i.e. diverging expectation value
of the volume. Any exponential tail with an asymptotic slope $P(N)\propto
e^{-BN}$ where $B\leq 3$ leads to this problem. Although the value $B=4.38\pm
0.09$ we find is safely above this limit, we cannot draw definite conclusions
due to the runs that did not finish, which follow a different distribution.
There is a related problem for any infinitely extended exponential tail. For
example, for $3<B\leq 6$ (including our relic case), the expectation value of
the volume is finite, but the variance diverges, so the expectation value does
not describe a typical member of the ensemble. Iterating this reasoning seems
to lead to a problem of interpretation for all $P(N)$ with any exponential (or
less steep) tail.

Figure 13: The probability distribution $P(N)$ for the relic case, similarly
to figure 4. Figure 14: Stochastic kicks as a function of e-folds for the
relic case, similarly to figure 5. The region where inflation temporarily ends
is marked in violet, and its overlap with the USR region is indicated by the
diagonal lines. Figure 15: The power spectrum $\mathcal{P}_{\mathcal{R}}(k)$
for the relic case, similarly to figure 6. The region where inflation
temporarily ends is marked in violet, and its overlap with the USR region is
indicated by the diagonal lines. Note that near this region, the same modes
are depicted multiple times as they re-enter and again exit the Hubble scale.
This leads to the dip in $\mathcal{P}_{\mathcal{R}}$ as a function of $N$ near
$N=48$. The PBH scale is located somewhat after the end of USR to capture the
peak in $\mathcal{P}_{\mathcal{R}}(k)$.

In the simplified treatment, the distribution is much sharper than the linear
theory Gaussian. Running $256\times 10^{7}$ simulations reaches only up to
$\Delta N=0.26$. Such a drastic narrowing of the classical distribution by
stochastic effects is not unknown. In [197], where the same simplified
treatment is used, stochastic effects significantly narrow $P(N)$, and reduce
$\bar{N}$ by an order of magnitude. For the potential studied there, the field
classically spends many e-folds in a flat region, while the noise kicks it out
fast. In our case the number of e-folds in the USR region is classically
small. The difference between the full and simplified treatment is due to the
difference in the kick amplitudes shown in figure 14. (The interval without
kicks corresponds to the temporary end of inflation.) The kicks in the full
case are at their maximum four orders of magnitude larger than in the
simplified case. Recall that we define USR as the time when the second SR
parameter $\eta=-\ddot{\bar{\phi}}/(H\dot{\bar{\phi}})$ is larger than 1, so
for some time after the transition from SR to USR, we have USR without
inflation. In the simplified case all the runs finish.

Perhaps surprisingly, the results of the full treatment again depend only
weakly on the value of $\sigma$: decreasing it to $10^{-8}$ changes
$\Omega_{\text{PBH}}$ only by 70%. Curiously, while we lose the peak modes by
considering such a small $\sigma$, it also enables earlier modes (which get
enhanced as we approach the end of inflation) to become classicalised and give
their kicks close to the end of inflation instead, giving almost the same
result. However, the fraction of simulations that finish is strongly affected.
For instance, for $\sigma=10^{-4}$, only 10% of the simulations complete, in
the rest the field gets stuck in the minimum of the potential. There does not
appear to be a clear pattern to the relationship between this fraction and the
value of $\sigma$.

In the relic case, our calculation also has a problem with gauge-dependence.
In contrast to the other cases, the effect of the transformation from the
spatially flat gauge to the uniform $N$ gauge is above 10% even away from the
narrow spike mentioned in section 2.1.3. It would be straightforward to take
into account the gauge transformation, but given the other problems with the
relic case, there is little point in correcting only this issue. Furthermore,
as discussed above, a proper derivation of the Langevin equations from first
principles, including a full accounting of the gauge issue, remains to be
done.

In summary, for the relic case we simply present the data obtained from our
simulations, acknowledging that no robust conclusion about the dark matter
relic abundance can be reached without a more detailed treatment.

##  4 Conclusions

Using the $\Delta N$ formalism, we have calculated numerically the probability
distribution $P(N)$ of e-folds of inflation in order to find the stochastic
correction to the abundance of primordial black hole (PBH) seeds produced
during ultra-slow-roll (USR). We have used the same set of stochastic
equations as in $Paper\,I$ [88], which for the first time consistently took
into account both the dependence of the noise that drives the coarse-grained
background on the evolution of the perturbations, and the change in the
evolution of the perturbations due to the change in the background. This leads
to a non-Markovian process.

We test our algorithm by considering USR in some toy model potentials used in
[89], finding results in close agreement. We then build realistic inflationary
potentials using as a starting point Higgs inflation [101, 102], where loop
corrections in the chiral Standard Model provide a feature at the desired PBH
scale [25]. We have adjusted these potentials by hand (except in the relic
case) in the CMB region to give the right spectral index. In the supermassive
case, we have also adjusted the overall height of the potential to give a
sufficiently small tensor-to-scalar ratio. (This would also be needed in the
solar case to fit the latest CMB data [169].) We have considered four cases,
tuning the potentials to have a shallow minimum that leads to USR. This leads
to the generation of roughly the observed abundance of dark matter in the form
of PBHs of asteroid, solar, or Planck mass, or alternatively the seeds of
supermassive black holes with density parameter $10^{-5}$ today, in the
Gaussian approximation. The stochastic effects generate an exponential tail
for $P(N)$, which dominates over the linear theory Gaussian [85, 89, 88]. The
stochastic noise enhances the PBH abundance by roughly 5, 1, 8, and 3 orders
of magnitude in the asteroid, solar, relic, and supermassive case,
respectively.

We compare our results to the usual treatment of stochastic inflation, where
the modes evolve on the FLRW background, and the noise amplitude is fixed to
its slow-roll super-Hubble limit, where it is given by the FLRW Hubble rate.
We also compare to a recursive treatment where the modes first evolve on the
FLRW background, and are then used as the noise for the next iteration [68,
69, 73]. The recursive treatment gives the same results as our full treatment
with just one iteration, while the simplified treatment significantly over-
predicts or under-predicts the stochastic effect. In the simplified treatment,
the enhancement over the Gaussian case is 2, 3 and 6 orders of magnitude in
the asteroid, solar and supermassive cases, respectively. In the relic case,
the simplified treatment predicts a negligible PBH abundance. This shows that
the evolution of the modes that produce the stochastic noise has to be taken
into account, whereas the change in the coarse-grained background in the
evolution of perturbations is not important, at least for the potentials we
consider. In other words, the non-Markovian nature of the process is
unimportant.

We set the coarse-graining scale by the requirement that the mode squeezing
parameter $r_{k}$ satisfies $\cosh(2r_{k})>100$ for all modes when they give
their kicks, so that they are sufficiently classical. The coarse-graining
scale gives a delay between the time when the modes exit the Hubble radius and
the time they join the background and give a kick to the evolution. The
results are rather insensitive to the choice of the coarse-graining scale, as
long as the delay is not so long that the amplitude of the kicks has changed
significantly between Hubble exit and coarse-graining.

We have also considered the gauge-dependence of our results. We have verified
numerically that in the asteroid, solar and supermassive cases the correction
from the gauge transformation between the spatially flat gauge to the uniform
$N$ gauge, computed as the ratio of the gauge correction over the mode
amplitude, is less than $10^{-5}$ for most modes. There is a narrow spike of
modes for which the gauge correction rises to a maximum of $\sim 10^{-2}$.
These correspond to modes with very small amplitudes at the coarse-graining
scale, and their contribution to the overall dynamics is small.

In the case of Planck mass relics the gauge transformations have a large
effect, and this model also suffers from other problems. Hence, no robust
conclusions may be drawn. In this case, USR occurs right before the end of
inflation, which leads to several issues. If we wanted each mode in the
simulation to classicalise before crossing the coarse-graining scale, the
modes that contribute most to PBH formation would not have enough time to give
their kicks. We therefore set a lower coarse-graining scale, meaning we miss
most kicks from the beginning of the simulation in order to preserve the ones
around the USR region. The effect of the kicks is large, as they grow towards
the end, and in $\sim$ 6% of the runs they slow down the field so much that it
gets stuck in the shallow minimum in the USR region, and inflation does not
end by the time we finish the simulation. This percentage is larger than the
fraction of simulations where $\Delta N$ exceeds the PBH formation threshold,
so our calculation of the abundance is not reliable. There is also the problem
that inflation temporarily ends before and during USR, and modes start
crossing back through the coarse-graining scale. It is not clear how to treat
this in the stochastic formalism. Furthermore, as the period of slow-roll
after USR is very short, the relation between the curvature perturbations and
the density contrast may not correspond to the usual one. Also, the assumption
that PBHs form during radiation domination may need to be corrected, as the
modes re-enter the Hubble radius quickly after inflation, so the collapse may
happen during preheating.

Our results show that in models with USR, it is crucial to check the impact of
stochastic effects on the curvature perturbation when considering extreme
perturbations such as PBH seeds. They also indicate that fixing the noise
amplitude to its slow-roll super-Hubble limit can be very inaccurate. The
quantitative impact depends on the details of the potential, and has to be
checked on a case by case basis. On a more technical side, to obtain the full
treatment results for our four potentials, we have run over 2 million CPU
hours of simulations in total. The need for such a considerable amount of
computational time may hinder doing more similar calculations. To reduce the
time needed to obtain $P(N)$, it is more practical to use the recursive
procedure, where the mode functions are only evolved once, and then used as
input for the noise, without having to evolve the modes again. The results of
this recursive treatment for all the four potentials we have considered is
indistinguishable (within the statistical errors) from those of the full
treatment, but the running time is between $1\%$ to $10\%$ of that in the full
treatment, depending on the potential.

As a final note, we remark that the derivation of the stochastic equations we
have used is heuristic, and would need to be worked out from first principles
in quantum field theory. The existing derivations that lead to explicit
Langevin equations assume slow-roll, and do not necessarily apply in USR.

## Acknowledgments

We thank Peter Johansson and Vincent Vennin, Jose María Ezquiaga and Juan
García-Bellido for discussion and correspondence. DGF (ORCID
0000-0002-4005-8915) is supported by a Ramón y Cajal contract with Ref.
RYC-2017-23493. This work was supported by the Generalitat Valenciana grant
PROMETEO/2021/083, the Ministerio de Ciencia e Innovación grant
PID2020-113644GB-I00, the Estonian Research Council grants PRG803, PRG1055,
and MOBTT5, and by the EU through the European Regional Development Fund CoE
program TK133 “The Dark Side of the Universe”. This work used computational
resources provided by the Finnish Grid and Cloud Infrastructure
(urn:nbn:fi:research-infras-2016072533). The authors wish to acknowledge CSC -
IT Center for Science, Finland, for computational resources.

##  Appendix A Inflaton potential

The inflaton potentials used in this work come from the quantum-corrected
potential of Higgs inflation in the Einstein frame, as described in detail in
[25]. In summary, the potential can be written as

$V[\phi(h)]=c_{0}\frac{\lambda_{\mathrm{eff}}(h)}{4}F(h)^{4}\ ,$ (A.1)

where $\phi$ is the Einstein frame canonical field. It is related to the
original Jordan frame Higgs field $h$ through the field transformation

$\frac{\mathrm{d}\phi}{\mathrm{d}h}=\frac{\sqrt{1+\xi
h^{2}+6\xi^{2}h^{2}}}{1+\xi h^{2}}\ ,$ (A.2)

where $\xi$ is a non-minimal coupling to gravity. The constant $c_{0}$ is an
extra scaling factor introduced by hand for the supermassive case to reduce
the tensor-to-scalar ratio $r$. The effective coupling
$\lambda_{\mathrm{eff}}(h)$ contains the one-loop Coleman–Weinberg corrections
evaluated so that all couplings run according to the one-loop beta functions
of the chiral Standard Model. The function $F(h)$ is

$F(h)\equiv\begin{cases}\frac{h}{\sqrt{1+\xi h^{2}}}\,,&h\leq h_{1}\\\
\text{interpolation}\,,&h_{1}\leq h\leq h_{2}\\\
c_{1}\frac{h+c_{2}h_{0}}{\sqrt{c_{3}+\xi(h+c_{2}h_{0})^{2}}}\,,&h\geq h_{2}\
,\end{cases}$ (A.3)

where the first expression is the same as in standard Higgs inflation and is
used near the PBH scale where $h<h_{1}$, and the third expression is modified
by the coefficients $c_{1},c_{2},c_{3}$, so that it produces perturbations
consistent with the observations in the CMB region where $h>h_{2}$. The
reference scale $h_{0}$ is the location of the feature in the potential. The
interpolation is done with an order five polynomial, chosen so that $F$ and
its first and second derivatives are continuous everywhere. The interpolation
gives the kinks in the potentials shown in figure 3.

The couplings featured in the quantum corrections are tuned so that the first
derivative of the potential vanishes at $h_{0}$ and the second derivative is
fixed by the second SR parameter
$\eta=\frac{1}{V(\phi)}\frac{\mathrm{d}^{2}V(\phi)}{\mathrm{d}\phi^{2}}$. In
table 2 we give $h_{0}$ in terms of $\delta_{0}\equiv\frac{1}{\xi h_{0}^{2}}$,
the value of $\eta$ at $h=h_{0}$, and the parameters in $F(h)$ for each of our
example potentials.

In our numerical calculations, we use the full renormalisation group improved
expression for $\lambda_{\mathrm{eff}}(h)$. It can be roughly approximated in
the different cases as

$\displaystyle\lambda_{\mathrm{eff}}(h)$ $\displaystyle\approx 2.556\times
10^{-6}+1.029\times 10^{-5}\ln\frac{F(h)}{F(h_{0})}+2.326\times
10^{-5}\ln^{2}\frac{F(h)}{F(h_{0})}$ (asteroid) (A.4)
$\displaystyle\lambda_{\mathrm{eff}}(h)$ $\displaystyle\approx 2.457\times
10^{-6}+9.799\times 10^{-6}\ln\frac{F(h)}{F(h_{0})}+2.316\times
10^{-5}\ln^{2}\frac{F(h)}{F(h_{0})}$ (solar)
$\displaystyle\lambda_{\mathrm{eff}}(h)$ $\displaystyle\approx 2.464\times
10^{-6}+9.829\times 10^{-6}\ln\frac{F(h)}{F(h_{0})}+2.312\times
10^{-5}\ln^{2}\frac{F(h)}{F(h_{0})}$ (supermassive)
$\displaystyle\lambda_{\mathrm{eff}}(h)$ $\displaystyle\approx 3.795\times
10^{-6}+1.1617\times 10^{-5}\ln\frac{F(h)}{F(h_{0})}+2.329\times
10^{-5}\ln^{2}\frac{F(h)}{F(h_{0})}$ $\displaystyle\text{(relic)}\ .$

| Asteroid | Solar | Supermassive | Relic
---|---|---|---|---
$\delta_{0}$ | 1.48137 | 1.484794 | 1.4878 | 1.5167
$\xi$ | 38.8 | 16.332(54668015058) | 12.75 | 75
$c_{0}$ | 1 | 1 | 0.1 | 1
$c_{1}$ | 0.998 | 0.985 | 0.958 | –
$c_{2}$ | 4.7 | 6.2 | 1.469 | –
$c_{3}$ | 2.5 | 5 | 0.28 | –
$h_{1}/h_{0}$ | 5.5 | 3 | 3 | $\infty$
$h_{2}/h_{0}$ | 7 | 4 | 4 | –
$\eta$ | 0.4141(0445743324537) | 0.141 | 0.1466 | 7.28023
Table 2: The parameter values used to build our example potentials. In the
relic case, the expression for the potential at $h<h_{1}$ applies for all $h$,
since there is no need to adjust the CMB scales. The overall scale factor
$c_{0}$ is only used for the supermassive black hole seeds. For the parameters
with many digits, the full accuracy is needed to produce our results reliably;
ignoring the digits in the parentheses gives Gaussian PBH abundances to within
$10\%$.

##  Appendix B Numerical implementation

In our discretisation procedure, a single time step consists of the classical
evolution of all quantities and a stochastic kick, where the kick is computed
using the newly evolved values. The continuum equations of motion for the
background and the perturbations read

$\displaystyle\tilde{\phi}^{\prime}$ $\displaystyle=\tilde{\pi}$ (B.1)
$\displaystyle\tilde{\pi}^{\prime}$
$\displaystyle=\left(\frac{1}{2}\tilde{\pi}^{2}-3\right)\left(\tilde{\pi}+\frac{\tilde{V}_{,\phi}}{\tilde{V}}\right)\equiv
f(\tilde{\phi},\tilde{\pi})$ (B.2)
$\displaystyle\delta{\tilde{\phi}}_{k}^{\prime}$
$\displaystyle=\delta\tilde{{\pi}}_{k}$ (B.3)
$\displaystyle\delta\tilde{{\pi}}_{k}^{\prime}$
$\displaystyle=\left(\frac{1}{2}\tilde{\pi}^{2}-3\right)\left(\delta\tilde{{\pi}}_{k}+\left[\tilde{\pi}^{2}+\frac{e^{2(\ln{\tilde{k}}-N)}}{\tilde{V}}+\frac{2\tilde{\pi}\tilde{V}_{,\phi}}{\tilde{V}}+\frac{\tilde{V}_{,\phi\phi}}{\tilde{V}}\right]\delta{\tilde{\phi}}_{k}\right)\equiv
g_{k}(N,\tilde{\phi},\tilde{\pi},\delta{\tilde{\phi}}_{k},\delta\tilde{{\pi}}_{k})\
,$ (B.4)

where prime denotes derivative with respect to the number of e-folds, and we
have introduced the dimensionless variables
$\tilde{\phi}\equiv\frac{\phi}{M_{\text{Pl}}}$,
$\tilde{\pi}\equiv\frac{\pi}{M_{\text{Pl}}}$,
$\delta{\tilde{\phi}}_{k}\equiv\sqrt{2H_{0}\tilde{k}^{3}}\delta{\phi}_{k}$,
$\delta\tilde{{\pi}}_{k}\equiv\sqrt{2H_{0}\tilde{k}^{3}}\delta{\pi}_{k}$,
$\tilde{V}\equiv\frac{V}{H_{0}^{2}M_{\text{Pl}}^{2}}$,
$\tilde{V}_{,\phi}\equiv\frac{V_{,\phi}}{H_{0}^{2}M_{\text{Pl}}}$,
$\tilde{V}_{,\phi\phi}\equiv\frac{V_{,\phi\phi}}{H_{0}^{2}}$, and
$\tilde{k}\equiv\frac{k}{H_{0}}$, with $H_{0}$ being the initial value of the
Hubble rate in the beginning of the simulation.

To numerically integrate our equations of motion, we use an explicit
Runge–Kutta method of order 4 with fixed time step $\differential{N}$.
Denoting consecutive time steps by $i$, the discretised solutions are

$\displaystyle\tilde{\phi}_{i+1}$
$\displaystyle=\tilde{\phi}_{i}+\differential{N}\tilde{\pi}_{i}+\frac{1}{6}\differential{N}^{2}(m_{1}+m_{2}+m_{3})$
(B.5) $\displaystyle\tilde{\pi}_{i+1}$
$\displaystyle=\tilde{\pi}_{i}+\frac{1}{6}\differential{N}(m_{1}+2m_{2}+2m_{3}+m_{4})$
(B.6) $\displaystyle\delta{\tilde{\phi}}_{k,i+1}$
$\displaystyle=\delta{\tilde{\phi}}_{k,i}+\differential{N}\delta\tilde{{\pi}}_{k,i}+\frac{1}{6}\differential{N}^{2}(n_{k,1}+n_{k,2}+n_{k,3})$
(B.7) $\displaystyle\delta\tilde{{\pi}}_{k,i+1}$
$\displaystyle=\delta\tilde{{\pi}}_{k,i}+\frac{1}{6}\differential{N}(n_{k,1}+2n_{k,2}+2n_{k,3}+n_{k,4})\
,$ (B.8)

where

$\displaystyle m_{1}$ $\displaystyle=f(\tilde{\phi}_{i},\tilde{\pi}_{i})$
(B.9) $\displaystyle m_{2}$
$\displaystyle=f(\tilde{\phi}_{i}+\tfrac{1}{2}\differential{N}\tilde{\pi}_{i},\tilde{\pi}_{i}+\tfrac{1}{2}\differential{N}m_{1})$
(B.10) $\displaystyle m_{3}$
$\displaystyle=f(\tilde{\phi}_{i}+\tfrac{1}{2}\differential{N}\tilde{\pi}_{i}+\tfrac{1}{4}\differential{N}^{2}m_{1},\tilde{\pi}_{i}+\tfrac{1}{2}\differential{N}m_{2})$
(B.11) $\displaystyle m_{4}$
$\displaystyle=f(\tilde{\phi}_{i}+\differential{N}\tilde{\pi}_{i}+\tfrac{1}{2}\differential{N}^{2}m_{2},\tilde{\pi}_{i}+\differential{N}m_{3})$
(B.12) $\displaystyle n_{k,1}$
$\displaystyle=g_{k}(N_{i},\tilde{\phi}_{i},\tilde{\pi}_{i},\delta{\tilde{\phi}}_{k,i},\delta\tilde{{\pi}}_{k,i})$
(B.13) $\displaystyle n_{k,2}$
$\displaystyle=g_{k}(N_{i}+\tfrac{1}{2}\differential{N},\tilde{\phi}_{i}+\tfrac{1}{2}\differential{N}\tilde{\pi}_{i},\tilde{\pi}_{i}+\tfrac{1}{2}\differential{N}m_{1},$
$\displaystyle\qquad\qquad\delta{\tilde{\phi}}_{k,i}+\tfrac{1}{2}\differential{N}\delta\tilde{{\pi}}_{k,i},\delta\tilde{{\pi}}_{k,i}+\tfrac{1}{2}\differential{N}n_{k,1})$
(B.14) $\displaystyle n_{k,3}$
$\displaystyle=g_{k}(N_{i}+\tfrac{1}{2}\differential{N},\tilde{\phi}_{i}+\tfrac{1}{2}\differential{N}\tilde{\pi}_{i}+\tfrac{1}{4}\differential{N}^{2}m_{1},\tilde{\pi}_{i}+\tfrac{1}{2}\differential{N}m_{2},$
$\displaystyle\qquad\qquad\delta{\tilde{\phi}}_{k,i}+\tfrac{1}{2}\differential{N}\delta\tilde{{\pi}}_{k,i}+\tfrac{1}{4}\differential{N}^{2}n_{k,1},\delta\tilde{{\pi}}_{k,i}+\tfrac{1}{2}\differential{N}n_{k,2})$
(B.15) $\displaystyle n_{k,4}$
$\displaystyle=g_{k}(N_{i}+\differential{N},\tilde{\phi}_{i}+\differential{N}\tilde{\pi}_{i}+\tfrac{1}{2}\differential{N}^{2}m_{2},\tilde{\pi}_{i}+\differential{N}m_{3},$
$\displaystyle\qquad\qquad\delta{\tilde{\phi}}_{k,i}+\differential{N}\delta\tilde{{\pi}}_{k,i}+\tfrac{1}{2}\differential{N}^{2}n_{k,2},\delta\tilde{{\pi}}_{k,i}+\differential{N}n_{k,3})\
,$ (B.16)

and the real and imaginary parts of the perturbation modes are solved
independently.

In the full treatment of the kicks, the stochastic field kicks
$\Delta{\tilde{\phi}_{i}}$ are normally distributed with variance

$\displaystyle\expectationvalue{(\Delta\tilde{\phi}_{i})^{2}}$
$\displaystyle=\evaluated{\zeta\differential{N}\left(1-\frac{1}{2}\tilde{\pi}_{i}^{2}\right)\absolutevalue{\delta{\tilde{\phi}}_{k,i}}^{2}}_{k=\sigma
aH}\ ,$ (B.17)

where $\zeta\equiv\frac{H_{0}^{2}}{4\pi^{2}M_{\text{Pl}}^{2}}$, and the
subscript $k=\sigma aH$ indicates that the mode amplitude is taken at the
coarse-graining scale. The momentum kicks are proportional to the realised
field kicks,

$\displaystyle\Delta{\tilde{\pi}_{i}}=\evaluated{\Delta{\tilde{\phi}_{i}}\frac{\real(\delta\tilde{{\pi}}_{k,i}\delta{\tilde{\phi}}_{k,i}^{*})}{\absolutevalue{\delta{\tilde{\phi}}_{k,i}}^{2}}}_{k=\sigma
aH}\ .$ (B.18)

The proportionality factor is equal to
$\real(\delta\tilde{{\pi}}_{k,i}/\delta\tilde{{\phi}}_{k,i})$ discussed below
(2.30), just written in a more convenient form.

In the recursive treatment, the mode amplitudes in the above expressions are
replaced with ones pre-calculated on a non-stochastic FLRW background, so that
no perturbation modes are evolved during the simulation. In the simplified
treatment, complexity is reduced even further by removing also the pre-
calculated modes, setting the field kick variance to

$\displaystyle\expectationvalue{(\Delta\tilde{\phi}_{i})^{2}}$
$\displaystyle=\zeta\differential{N}\frac{\tilde{V}}{3-\frac{1}{2}\tilde{\pi}_{i}^{2}}\
,$ (B.19)

and removing the momentum kicks entirely.

For the perturbations, we consider a discrete grid of modes with modulus
evenly distributed on a logarithmic scale as
$\ln(k_{i+1})=\ln(k_{i})+\differential{k}$, where $\differential{k}$ is the
mode spacing. The evolution of each mode begins when $k=\alpha aH$, with
$\alpha=100$ (the results are insensitive to making $\alpha$ larger). The
initial conditions are $\delta{\tilde{\phi}}_{k}=\alpha\tilde{H}$ and
$\delta{\tilde{\pi}}_{k}=-\alpha\tilde{H}-i\alpha^{2}\tilde{H}$, where
$\tilde{H}\equiv\frac{H}{H_{0}}$. The longest wavelength mode we consider
corresponds to the CMB pivot scale $k_{*}$, and its evolution starts
immediately at the onset of each simulation. The next mode to begin evolving
is denoted by $k_{\mathrm{next}}$. For each realisation, the code executes the
algorithm below:

Set initial values for $N$, $\bar{\phi}$, $\bar{\pi}$. Set
$k_{\mathrm{next}}=k_{\mathrm{*}}$. Set current kick coefficient to zero.

while  _$\bar{\phi} >{\bar{\phi}}_{\mathrm{f}}$_ do

Evolve $N$, $\bar{\phi}$, $\bar{\pi}$.

for  _all modes k in the simulation_ do

if  _$k >\sigma aH$_ then

Evolve $\delta\phi_{\mathbf{k}}$, $\delta\phi_{\mathbf{k}}^{\prime}$.

else

Evolve $\delta\phi_{\mathbf{k}}$, $\delta\phi_{\mathbf{k}}^{\prime}$ to
$k=\sigma aH$. Set the current kick coefficient from
$\delta\phi_{\mathbf{k}}$, $\delta\phi_{\mathbf{k}}^{\prime}$. Remove mode $k$
from the simulation.

if  _$k_{\mathrm{next}}\leq k_{\mathrm{PBH}}$_ then

if  _$k_{\mathrm{next}}\leq\alpha aH$_ then

Add mode $k=k_{\mathrm{next}}$ to the simulation. Set initial values for
$\delta\phi_{\mathbf{k}}$, $\delta\phi_{\mathbf{k}}^{\prime}$. Evolve
$\delta\phi_{\mathbf{k}}$, $\delta\phi_{\mathbf{k}}^{\prime}$ from $k=\alpha
aH$. Set $k_{\mathrm{next}}=e^{\differential{k}}k_{\mathrm{next}}$.

else

if  _$k_{\mathrm{next}}\leq\sigma aH$_ then

Set the current kick coefficient to zero.

if  _$aH\geq\max(aH)$_ then

Add stochastic kick to $\bar{\phi}$, $\bar{\pi}$ using the current kick
coefficient.

Algorithm 1 Evolution for each run

Initially, the background field is in slow-roll. Although our code generally
runs on a fixed time step, it takes partial steps for increased accuracy
whenever we encounter events that are not synchronised with the time steps,
e.g. when we start or stop evolving modes, and at the end of inflation. The
condition $aH\geq\max(aH)$ guarantees that we do not apply duplicate kicks in
the event that $aH$ momentarily decreases in the middle of the simulation, as
happens for the relic potential.

The values used for the time step $\differential{N}$ and the mode spacing
$\differential{k}$ depend on the potential. The time steps are generally
between $1/256$ and $1/512$, and the mode spacings are between $1/32$ and
$1/256$. These values are the result of extensive testing for each of the
potentials; they have been chosen so that any further halving of the time step
or mode spacing would change the results by less than 1%.

##  Appendix C Gauge transformations

Under a gauge transformation from the spatially flat gauge to the uniform $N$
gauge, the field perturbation and its derivative transform as

$\displaystyle\delta{\phi}_{k}$
$\displaystyle\rightarrow\delta{\phi}_{k}+\phi^{\prime}\alpha_{k}$ (C.1)
$\displaystyle\delta{\phi}_{k}^{\prime}$
$\displaystyle\rightarrow\delta{\phi}_{k}^{\prime}+(\phi^{\prime}\alpha_{k})^{\prime}\
,$ (C.2)

where $\alpha_{k}$ is a gauge transformation parameter defined in [74], and in
this appendix prime exceptionally denotes a derivative with respect to
conformal time.

Figure 16 shows the amplitude of the gauge correction terms relative to the
corresponding perturbation at the coarse-graining scale. In the asteroid,
solar and supermassive cases, the relative correction is less than $10^{-5}$
for most modes. However, there is a single spike in each plot just before or
during USR. For these modes the amplitude of the perturbation (or its
derivative) in the denominator becomes small, driving up the value of the
fraction. The contribution of such modes to the overall dynamics is very
small. For the relic case the relative gauge corrections are much larger,
reaching almost order unity in general, and even larger for the spikes. Note
that we show the modes at the coarse-graining scale, not at Hubble exit. In
the relic case, there are of order unity gauge corrections during USR even at
Hubble crossing. In [74] it was assumed that in USR the potential can be
neglected completely, and that $\sigma\ll 1$. These approximations do not hold
in our case.

Figure 16: Relative gauge transformation corrections to $\delta{\phi}_{k}$
(left column) and $\delta{\phi}_{k}^{\prime}$ (right column) at the coarse-
graining scale for the asteroid, solar, supermassive and relic cases, in order
from top to bottom. The green area indicates the USR region, the violet area
indicates the period where inflation temporarily ends, and its overlap with
the USR region is marked by the diagonal lines.

## References

  * [1] Planck collaboration, Y. Akrami et al., _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10, [1807.06211].
  * [2] A. A. Starobinsky, _Spectrum of relic gravitational radiation and the early state of the universe_ , _JETP Lett._ 30 (1979) 682–685.
  * [3] A. A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett._ 91B (1980) 99–102.
  * [4] D. Kazanas, _Dynamics of the Universe and Spontaneous Symmetry Breaking_ , _Astrophys. J._ 241 (1980) L59–L63.
  * [5] A. H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys. Rev._ D23 (1981) 347–356.
  * [6] K. Sato, _First Order Phase Transition of a Vacuum and Expansion of the Universe_ , _Mon. Not. Roy. Astron. Soc._ 195 (1981) 467–479.
  * [7] V. F. Mukhanov and G. V. Chibisov, _Quantum Fluctuations and a Nonsingular Universe_ , _JETP Lett._ 33 (1981) 532–535.
  * [8] A. D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Phys. Lett._ 108B (1982) 389–393.
  * [9] A. Albrecht and P. J. Steinhardt, _Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking_ , _Phys. Rev. Lett._ 48 (1982) 1220–1223.
  * [10] S. W. Hawking and I. G. Moss, _Supercooled Phase Transitions in the Very Early Universe_ , _Phys. Lett._ 110B (1982) 35–38.
  * [11] G. V. Chibisov and V. F. Mukhanov, _Galaxy formation and phonons_ , _Mon. Not. Roy. Astron. Soc._ 200 (1982) 535–550.
  * [12] S. W. Hawking, _The Development of Irregularities in a Single Bubble Inflationary Universe_ , _Phys. Lett._ 115B (1982) 295.
  * [13] A. H. Guth and S. Y. Pi, _Fluctuations in the New Inflationary Universe_ , _Phys. Rev. Lett._ 49 (1982) 1110–1113.
  * [14] A. A. Starobinsky, _Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations_ , _Phys. Lett._ 117B (1982) 175–178.
  * [15] M. Sasaki, _Large Scale Quantum Fluctuations in the Inflationary Universe_ , _Prog. Theor. Phys._ 76 (1986) 1036.
  * [16] V. F. Mukhanov, _Quantum Theory of Gauge Invariant Cosmological Perturbations_ , _Sov. Phys. JETP_ 67 (1988) 1297–1302.
  * [17] J. Garcia-Bellido and E. Ruiz Morales, _Primordial black holes from single field models of inflation_ , _Phys. Dark Univ._ 18 (2017) 47–54, [1702.03901].
  * [18] J. M. Ezquiaga, J. Garcia-Bellido and E. Ruiz Morales, _Primordial Black Hole production in Critical Higgs Inflation_ , _Phys. Lett._ B776 (2018) 345–349, [1705.04861].
  * [19] K. Kannike, L. Marzola, M. Raidal and H. VeermÃ¤e, _Single Field Double Inflation and Primordial Black Holes_ , _JCAP_ 1709 (2017) 020, [1705.06225].
  * [20] M. P. Hertzberg and M. Yamada, _Primordial Black Holes from Polynomial Potentials in Single Field Inflation_ , _Phys. Rev. D_ 97 (2018) 083509, [1712.09750].
  * [21] C. Germani and T. Prokopec, _On primordial black holes from an inflection point_ , _Phys. Dark Univ._ 18 (2017) 6–10, [1706.04226].
  * [22] H. Motohashi and W. Hu, _Primordial Black Holes and Slow-Roll Violation_ , _Phys. Rev._ D96 (2017) 063503, [1706.06784].
  * [23] Y. Gong and Y. Gong, _Primordial black holes and second order gravitational waves from ultra-slow-roll inflation_ , _JCAP_ 1807 (2018) 007, [1707.09578].
  * [24] G. Ballesteros and M. Taoso, _Primordial black hole dark matter from single field inflation_ , _Phys. Rev._ D97 (2018) 023501, [1709.05565].
  * [25] S. Räsänen and E. Tomberg, _Planck scale black hole dark matter from Higgs inflation_ , _JCAP_ 01 (2019) 038, [1810.12608].
  * [26] M. Drees and Y. Xu, _Overshooting, Critical Higgs Inflation and Second Order Gravitational Wave Signatures_ , _Eur. Phys. J. C_ 81 (2021) 182, [1905.13581].
  * [27] M. Biagetti, G. Franciolini, A. Kehagias and A. Riotto, _Primordial Black Holes from Inflation and Quantum Diffusion_ , _JCAP_ 07 (2018) 032, [1804.07124].
  * [28] J. M. Ezquiaga and J. Garcia-Bellido, _Quantum diffusion beyond slow-roll: implications for primordial black-hole production_ , _JCAP_ 08 (2018) 018, [1805.06731].
  * [29] G. Ballesteros, J. Rey, M. Taoso and A. Urbano, _Stochastic inflationary dynamics beyond slow-roll and consequences for primordial black hole formation_ , _JCAP_ 2008 (2020) 043, [2006.14597].
  * [30] C. Fu, P. Wu and H. Yu, _Primordial black holes and oscillating gravitational waves in slow-roll and slow-climb inflation with an intermediate noninflationary phase_ , _Phys. Rev. D_ 102 (2020) 043527, [2006.03768].
  * [31] G. F. Chapline, _Cosmological effects of primordial black holes_ , _Nature_ 253 (1975) 251–252.
  * [32] A. Dolgov and J. Silk, _Baryon isocurvature fluctuations at small scales and baryonic dark matter_ , _Phys. Rev._ D47 (1993) 4244–4255.
  * [33] K. Jedamzik, _Primordial black hole formation during the QCD epoch_ , _Phys. Rev._ D55 (1997) R5871(R)–R5875(R), [astro-ph/9605152].
  * [34] P. Ivanov, P. Naselsky and I. Novikov, _Inflation and primordial black holes as dark matter_ , _Phys. Rev._ D50 (1994) 7173–7178.
  * [35] J. Garcia-Bellido, A. D. Linde and D. Wands, _Density perturbations and black hole formation in hybrid inflation_ , _Phys. Rev._ D54 (1996) 6040–6058, [astro-ph/9605094].
  * [36] J. Yokoyama, _Formation of MACHO primordial black holes in inflationary cosmology_ , _Astron. Astrophys._ 318 (1997) 673, [astro-ph/9509027].
  * [37] P. Ivanov, _Nonlinear metric perturbations and production of primordial black holes_ , _Phys. Rev._ D57 (1998) 7145–7154, [astro-ph/9708224].
  * [38] D. Blais, C. Kiefer and D. Polarski, _Can primordial black holes be a significant part of dark matter?_ , _Phys. Lett._ B535 (2002) 11–16, [astro-ph/0203520].
  * [39] B. Carr and F. Kuhnel, _Primordial Black Holes as Dark Matter: Recent Developments_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 355–394, [2006.02838].
  * [40] A. M. Green and B. J. Kavanagh, _Primordial Black Holes as a dark matter candidate_ , _J. Phys. G_ 48 (2021) 043001, [2007.10722].
  * [41] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _Constraints on Primordial Black Holes_ , 2002.12778.
  * [42] A. Arbey, J. Auffinger and J. Silk, _Stellar Signatures of Inhomogeneous Big Bang Nucleosynthesis_ , 2006.02446.
  * [43] K. Jedamzik and L. Pogosian, _Relieving the Hubble tension with primordial magnetic fields_ , 2004.09487.
  * [44] M. Rashkovetskyi, J. B. Muñoz, D. J. Eisenstein and C. Dvorkin, _Small-scale Clumping at Recombination and the Hubble Tension_ , 2108.02747.
  * [45] D. Salopek and J. Stewart, _Hamilton-Jacobi theory for general relativity with matter fields_ , _Class. Quant. Grav._ 9 (1992) 1943–1968.
  * [46] D. Wands, K. A. Malik, D. H. Lyth and A. R. Liddle, _A New approach to the evolution of cosmological perturbations on large scales_ , _Phys. Rev._ D62 (2000) 043527, [astro-ph/0003278].
  * [47] S. Habib, _Stochastic inflation: The Quantum phase space approach_ , _Phys. Rev. D_ 46 (1992) 2408–2427, [gr-qc/9208006].
  * [48] E. Calzetta and B. L. Hu, _Quantum fluctuations, decoherence of the mean field, and structure formation in the early universe_ , _Phys. Rev._ D52 (1995) 6770–6788, [gr-qc/9505046].
  * [49] D. Polarski and A. A. Starobinsky, _Semiclassicality and decoherence of cosmological perturbations_ , _Class. Quant. Grav._ 13 (1996) 377–392, [gr-qc/9504030].
  * [50] M. Bellini, H. Casini, R. Montemayor and P. Sisterna, _Stochastic approach to inflation: Classicality conditions_ , _Phys. Rev. D_ 54 (1996) 7172–7180.
  * [51] J. Lesgourgues, D. Polarski and A. A. Starobinsky, _Quantum to classical transition of cosmological perturbations for nonvacuum initial states_ , _Nucl. Phys._ B497 (1997) 479–510, [gr-qc/9611019].
  * [52] C. Kiefer and D. Polarski, _Emergence of classicality for primordial fluctuations: Concepts and analogies_ , _Annalen Phys._ 7 (1998) 137–158, [gr-qc/9805014].
  * [53] H. Casini, R. Montemayor and P. Sisterna, _Stochastic approach to inflation. 2. Classicality, coarse graining and noises_ , _Phys. Rev. D_ 59 (1999) 063512, [gr-qc/9811083].
  * [54] C. Kiefer, D. Polarski and A. A. Starobinsky, _Quantum to classical transition for fluctuations in the early universe_ , _Int. J. Mod. Phys._ D07 (1998) 455–462, [gr-qc/9802003].
  * [55] C. Kiefer and D. Polarski, _Why do cosmological perturbations look classical to us?_ , _Adv. Sci. Lett._ 2 (2009) 164–173, [0810.0087].
  * [56] A. A. Starobinsky, _STOCHASTIC DE SITTER (INFLATIONARY) STAGE IN THE EARLY UNIVERSE_ , _Lect. Notes Phys._ 246 (1986) 107–126.
  * [57] M. Morikawa, _Dissipation and Fluctuation of Quantum Fields in Expanding Universes_ , _Phys. Rev. D_ 42 (1990) 1027–1034.
  * [58] A. A. Starobinsky and J. Yokoyama, _Equilibrium state of a selfinteracting scalar field in the De Sitter background_ , _Phys. Rev. D_ 50 (1994) 6357–6368, [astro-ph/9407016].
  * [59] M. Mijic, _Stochastic dynamics of coarse grained quantum fields in the inflationary universe_ , _Phys. Rev. D_ 49 (1994) 6434–6441, [gr-qc/9401030].
  * [60] J. Martin and M. Musso, _On the reliability of the Langevin perturbative solution in stochastic inflation_ , _Phys. Rev. D_ 73 (2006) 043517, [hep-th/0511292].
  * [61] N. Tsamis and R. Woodard, _Stochastic quantum gravitational inflation_ , _Nucl. Phys. B_ 724 (2005) 295–328, [gr-qc/0505115].
  * [62] R. P. Woodard, _A Leading logarithm approximation for inflationary quantum field theory_ , _Nucl. Phys. B Proc. Suppl._ 148 (2005) 108–119, [astro-ph/0502556].
  * [63] M. van der Meulen and J. Smit, _Classical approximation to quantum cosmological correlations_ , _JCAP_ 11 (2007) 023, [0707.0842].
  * [64] F. Finelli, G. Marozzi, A. Starobinsky, G. Vacca and G. Venturi, _Generation of fluctuations during inflation: Comparison of stochastic and field-theoretic approaches_ , _Phys. Rev. D_ 79 (2009) 044007, [0808.1786].
  * [65] M. Beneke and P. Moch, _On “dynamical mass” generation in Euclidean de Sitter space_ , _Phys. Rev. D_ 87 (2013) 064018, [1212.3058].
  * [66] F. Gautier and J. Serreau, _Infrared dynamics in de Sitter space from Schwinger-Dyson equations_ , _Phys. Lett. B_ 727 (2013) 541–547, [1305.5705].
  * [67] B. Garbrecht, G. Rigopoulos and Y. Zhu, _Infrared correlations in de Sitter space: Field theoretic versus stochastic approach_ , _Phys. Rev. D_ 89 (2014) 063506, [1310.0367].
  * [68] L. Perreault Levasseur, _Lagrangian formulation of stochastic inflation: Langevin equations, one-loop corrections and a proposed recursive approach_ , _Phys. Rev. D_ 88 (2013) 083537, [1304.6408].
  * [69] L. Perreault Levasseur, V. Vennin and R. Brandenberger, _Recursive Stochastic Effects in Valley Hybrid Inflation_ , _Phys. Rev. D_ 88 (2013) 083538, [1307.2575].
  * [70] T. Fujita, M. Kawasaki, Y. Tada and T. Takesako, _A new algorithm for calculating the curvature perturbations in stochastic inflation_ , _JCAP_ 1312 (2013) 036, [1308.4754].
  * [71] T. Fujita, M. Kawasaki and Y. Tada, _Non-perturbative approach for curvature perturbations in stochastic $\delta N$ formalism_, _JCAP_ 10 (2014) 030, [1405.2187].
  * [72] B. Garbrecht, F. Gautier, G. Rigopoulos and Y. Zhu, _Feynman Diagrams for Stochastic Inflation and Quantum Field Theory in de Sitter Space_ , _Phys. Rev. D_ 91 (2015) 063520, [1412.4893].
  * [73] L. Perreault Levasseur and E. McDonough, _Backreaction and Stochastic Effects in Single Field Inflation_ , _Phys. Rev. D_ 91 (2015) 063513, [1409.7399].
  * [74] C. Pattison, V. Vennin, H. Assadullahi and D. Wands, _Stochastic inflation beyond slow roll_ , _JCAP_ 07 (2019) 031, [1905.06300].