# Stochastic inflation beyond slow roll: noise modelling and importance sampling

Joseph H.P Jackson, Hooshyar Assadullahi, Andrew D Gow, Kazuya Koyama, Vincent Vennin, David Wands

► To cite this version:

Joseph H.P Jackson, Hooshyar Assadullahi, Andrew D Gow, Kazuya Koyama, Vincent Vennin, et al. Stochastic inflation beyond slow roll: noise modelling and importance sampling. Journal of Cosmology and Astroparticle Physics, 2025, 2025 (04), pp.073. 10.1088/1475-7516/2025/04/073. hal-04751495

HAL Id: hal-04751495

https://hal.science/hal-04751495v1

Submitted on 19 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
PAPER • OPEN ACCESS

# Stochastic inflation beyond slow roll: noise modelling and importance sampling

To cite this article: Joseph H.P. Jackson et al JCAP04(2025)073

View the [article online](#) for updates and enhancements.

You may also like

- A relativistic calculation of super-Hubble suppression of inflation with thermal dissipation
Wolung Lee and Li-Zhi Fang

- A natural model for curved inflation
Quentin Vigneron and Julien Larena

- Trans-Planckian issues for inflationary cosmology
Robert H Brandenberger and Jérôme Martin
# Stochastic inflation beyond slow roll: noise modelling and importance sampling

Joseph H.P. Jackson<sup>1, a</sup> Hooshyar Assadullahi<sup>1, a, b</sup> Andrew D. Gow<sup>1, a</sup>

Kazuya Koyama<sup>1, a</sup> Vincent Vennin<sup>1, c, a</sup> and David Wands<sup>1, a</sup>

E-mail: joseph.jackson@port.ac.uk, hooshyar.assadullahi@port.ac.uk, andrew.gow@port.ac.uk, kazuya.koyama@port.ac.uk, vincent.vennin@ens.fr, david.wands@port.ac.uk

**ABSTRACT:** We simulate the distribution of very rare, large excursions in the primordial density field produced in models of inflation in the very early universe which include a strong enhancement of the power spectrum. The stochastic $\delta N$ formalism is used to identify the probability distribution for the primordial curvature perturbation with the first-passage-time distribution, $P(\delta N)$, and we compare our stochastic results with those obtained in the classical $\delta N$ approach. We extend the PYFPT numerical code to simulate the full 2D phase space, and apply importance sampling which allows very rare fluctuations to be simulated in $O(10)$ minutes on a single CPU, where previous direct simulations required supercomputers. We demonstrate that the stochastic noise due to quantum fluctuations after a sudden transition to ultra-slow roll can be accurately modelled using an analytical Bessel-function ansatz to identify the homogeneous growing mode. The stochastic noise found in this way is a function of the field value only. This enables us to coarse grain the inflation field at the Hubble scale and include non-linear, stochastic evolution on all super-Hubble length scales.

**KEYWORDS:** inflation, physics of the early universe, primordial black holes, Statistical sampling techniques

ARXIV EPRINT: 2410.13683

<figure><img src="image_2.png" /></figure>

<footnote><sup>a</sup>Institute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama Building, Burnaby Road, Portsmouth, PO1 3FX, U.K.</footnote>

<footnote><sup>b</sup>School of Mathematics and Physics, University of Portsmouth, Lion Gate Building, Lion Terrace, Portsmouth, PO1 3HF, U.K.</footnote>

<footnote><sup>c</sup>Laboratoire de Physique de l'École Normale Supérieure, ENS, CNRS, Université PSL, Sorbonne Université, Université Paris Cité, F-75005 Paris, France</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 The separate-universe approach and δN</td><td>3</td></tr><tr><td>2.1 Background phase space and the separate-universe approach</td><td>3</td></tr><tr><td>2.2 Classical δN formalism</td><td>5</td></tr><tr><td>2.3 Stochastic formalism</td><td>7</td></tr><tr><td>3 Noise models</td><td>9</td></tr><tr><td>3.1 Consistency criterion</td><td>9</td></tr><tr><td>3.2 Analytical approximations</td><td>9</td></tr><tr><td>3.3 Numerical matching</td><td>10</td></tr><tr><td>3.4 Semi-analytical approach</td><td>11</td></tr><tr><td>4 Stochastic δN simulations and results</td><td>13</td></tr><tr><td>4.1 Numerical approach</td><td>13</td></tr><tr><td>4.2 Piece-wise linear potential</td><td>16</td></tr><tr><td>4.3 Gaussian bump</td><td>18</td></tr><tr><td>5 Conclusions</td><td>21</td></tr><tr><td>A Noise model accuracy</td><td>24</td></tr></table>

## 1 Introduction

Inflation is a period of accelerated expansion in the very early universe which can explain the observed isotropy, homogeneity and flatness of our Universe [1–6]. Its prediction of an almost scale-invariant spectrum of primordial perturbations is in excellent agreement with current Cosmic Microwave Background (CMB) observations [7, 8]. These observations constrain the distribution of density fluctuations, which originate from quantum vacuum fluctuations during inflation but are stretched to macroscopic scales, thereby constraining the dynamics of inflation on large scales [9–14].

There are fewer constraints on smaller scales during inflation, leading to much interest in the possibility of larger perturbations on these scales [15]. Sufficiently large over-densities can collapse post-inflation to form primordial black holes (PBHs) [16–18]. These are a dark matter candidate [19–22] and could possibly explain several other observational puzzles as well [23]. Large scalar perturbations can also source a stochastic background of gravitational waves at second order in perturbation theory [24, 25] and therefore possibly explain the recent pulsar-timing array detections [26, 27].

Very large and rare perturbations are required to produce a measurable abundance of PBHs [28]. However, such an enhanced power spectrum also brings into question the validity
of linear calculations, and so a non-linear treatment of fluctuations is required to reliably
calculate PBH abundances [29].

Stochastic inflation is an effective-theory approach to modelling non-linear inflationary
dynamics [30, 31]. This is done by mapping the quantum field theory in curved spacetime
onto a stochastic description, splitting the full quantum field into short- and long-wavelength
regimes. While perturbative quantum field theory is used to describe the evolution of short-
scale modes in a local patch, the long-wavelength regime is modelled using a separate-universe
approach, where above the coarse-graining scale the fields obey the evolution equations for a
locally homogeneous and isotropic (Friedmann-Lemaître-Robertson-Walker, FLRW) universe.
As small-scale quantum fluctuations are stretched and cross the coarse-graining scale, they
give a stochastic kick to the long-wavelength fields.

In the $\delta N$ formalism [11, 32–35], the non-linear primordial density perturbation is identified with the fluctuation in the local duration of inflation in different parts of our Universe. In general there is an exponential enhancement in the probability of large perturbations found using the stochastic $\delta N$ formalism compared to linear theory [29, 36–39]. However, computing the amplitude of this enhancement can be difficult. Adopting a brute force approach, billions of stochastic simulations are required to find the few rare events with perturbations large enough to produce PBHs, requiring supercomputers to directly sample the expected distribution [40, 41]. Here for the first time we apply importance sampling [42, 43] to inflationary models which experience a strong enhancement of density perturbations on small scales. With this approach, the events of interest can be directly simulated and their probability found, reducing the number of simulations required by orders of magnitude. What would otherwise require many hours on a supercomputer takes only minutes on a laptop.

We compare the results from our stochastic simulations with those obtained in the classical $\delta N$ approach, where a scalar field perturbation is applied at one single time during inflation and the subsequent non-linear, but classical evolution is followed to the end of inflation. In this case the scalar field fluctuations are computed about a fixed classical background, rather than the stochastically varying local background. Nonetheless, the non-linear evolution subsequently in the long-wavelength limit can lead to an enhancement of the probability of large fluctuations with respect to linear theory [44–49] and in some cases are expected to give a good description of the stochastic simulations [50, 51]. The classical $\delta N$ results can therefore be used to test the stochastic noise model, found using the classical background for a fixed range of modes, and the importance sampling approach proposed here, which is the main focus of this work.

Either formalism uses the separate-universe approach which requires one to identify the
long-wavelength limit (what we will refer to as the homogeneous part) of the linearised field
perturbation with the inhomogeneous field perturbation at a finite scale, usually at, or soon
after, Hubble exit during inflation. Although this is relatively straight-forward during slow
roll, it is non-trivial in models with a sudden transition and non-slow-roll evolution [52, 53],
such as periods of ultra-slow roll [54].

In this paper we first introduce the background equations for the separate-universe
approach, the $\delta N$ formalism and stochastic inflation in section 2. A classical $\delta N$ formula is
also derived, for comparison with the stochastic results. In section 3, different methods for
finding the stochastic noise are discussed. We argue that the noise remains one-dimensional in the two-dimensional phase space, even immediately after the sudden transition. Post-transition we find the noise to be only a function of the field-space position, and not the full trajectory. Numerical results of our stochastic simulations are presented in section 4, along with details about numerical modelling methods. We present our conclusions and possible directions for future work in section 5.

## 2 The separate-universe approach and $\delta N$

### 2.1 Background phase space and the separate-universe approach

In canonical single-field inflation, the universe is dominated by a scalar field $\phi$ with a potential $V(\phi)$ at very early times. The Klein-Gordon equation of motion is

$$ \ddot{\phi} + 3H\dot{\phi} + \frac{dV(\phi)}{d\phi} = 0, \qquad (2.1) $$

where a dot denotes a cosmic time derivative. The Hubble rate, $H$, is given by the Friedmann constraint equation

$$ 3H^2 = V(\phi) + \frac{1}{2}\dot{\phi}^2, \qquad (2.2) $$

where we use natural units, such that $8\pi G = 1$. One way to characterise the dynamics of inflation is to use the Hubble-flow parameters

$$ \epsilon_{i+1} = \frac{d \ln \epsilon_i}{dN} \quad \text{where} \quad \epsilon_0 \equiv \frac{H_{\text{in}}}{H}. \qquad (2.3) $$

In particular, $\epsilon_1 = \dot{\phi}^2/2H^2$, and for inflation to occur one requires $\epsilon_1 < 1$.

We will first consider linear perturbation theory to describe scalar field fluctuations during inflation. In the spatially-flat gauge, which we shall use unless otherwise specified, this gives an equation of motion for each Fourier mode with comoving wavenumber $k$ [55]

$$ \ddot{\delta\phi}_k + 3H\dot{\delta\phi}_k + \left[ \frac{k^2}{a^2} + \frac{d^2V}{d\phi^2} - a^{-3}\frac{d}{dt}\left(\frac{a^3\dot{\phi}^2}{H}\right) \right] \delta\phi_k = 0. \qquad (2.4) $$

If we quantise the field fluctuations and assume they start in the Bunch-Davies vacuum state at early times during inflation, then we have the particular solution on small scales

$$ \delta\phi_k^{\text{BD}} = \frac{e^{-ik\eta}}{a\sqrt{2k}} \quad \text{for} \quad k \gg aH, \qquad (2.5) $$

where the conformal time $\eta \equiv \int dt/a$.

During inflation small, sub-Hubble scales ($k \gg aH$) are stretched up to super-Hubble scales ($k \ll aH$). In the long-wavelength limit, where $k/aH \to 0$, eq. (2.4) reduces to

$$ \ddot{\delta\phi}_h + 3H\dot{\delta\phi}_h + \left[ \frac{d^2V}{d\phi^2} - a^{-3}\frac{d}{dt}\left(\frac{a^3\dot{\phi}^2}{H}\right) \right] \delta\phi_h = 0, \qquad (2.6) $$
which is consistent with the separate-universe approach [34, 35, 56, 57]. Here, $\delta\phi_h$ is the perturbation to the field value in a locally-homogeneous super-Hubble patch, and it follows the homogeneous behaviour on super-Hubble scales. Specifically, one needs to evolve the local Klein-Gordon equation (2.1) using the local Hubble rate obtained from (2.2), while allowing for perturbations of the local proper time with respect to the background cosmic time.<fnref n="1" /> The locally-homogeneous dynamics are expected to give a good approximation of the full inhomogenous dynamics on sufficiently large scales.

The general solution of eq. (2.6), in terms of the number of e-folds $N = \ln(a)$, is given by

$$ \delta\phi_h(N) = C\sqrt{2\epsilon_1(N)} + D\sqrt{2\epsilon_1(N)} \int_N^{N_{\text{end}}} \frac{d\tilde{N}}{2\tilde{a}^3\tilde{H}\tilde{\epsilon}_1}, \quad (2.7) $$

where $C$ and $D$ are constants of integration and $N_{\text{end}}$ is the number of e-folds at the end of inflation. The first term is known as the *growing mode* and the second as the *decaying mode*. By construction the decaying mode vanishes at the end of inflation.

The comoving curvature perturbation during inflation is given by $\mathcal{R} = H\delta\phi/\dot{\phi}$ in linear theory.<fnref n="2" /> Its power spectrum can be related to the modulus-squared of the mode function for the field perturbation during inflation

$$ \mathcal{P}_{\mathcal{R}}(k) = \frac{k^3 |\delta\phi_k|^2}{2\pi^2 2\epsilon_1}. \quad (2.8) $$

We see from eqs. (2.7) and (2.8) that $\mathcal{P}_{\mathcal{R}}$ remains constant in the long-wavelength limit when the decaying mode becomes negligible, corresponding to adiabatic perturbations [52, 59, 60]. More generally it is the amplitude of the growing mode on long wavelengths, $C$ in eq. (2.7), that determines the comoving curvature perturbation at the end of inflation which is then used to characterise the amplitude of primordial density fluctuations after inflation. Thus we are interested in determining the “homogeneous” part $\delta\phi_{k,h}$ of the mode function $\delta\phi_k$ for each wave number $k$ that survives in the long-wavelength limit as $k/aH \to 0$. If the behaviour of $\delta\phi_{k,h}(N)$ is known at a given time $N_*$ (for example, at or soon after Hubble exit, $k = aH$), then the integration constants $C_k$ and $D_k$ can be assessed via continuous matching to the solution (2.7),

$$ \begin{aligned} C_k &= \frac{\delta\phi_{k,h*}}{\sqrt{2\epsilon_{1*}}} - a_*^3 H_* \sqrt{2\epsilon_{1*}} \left( \frac{\epsilon_{2*}}{2} \delta\phi_{k,h*} - \partial_N \delta\phi_{k,h*} \right) \int_{N_*}^{N_{\text{end}}} \frac{d\tilde{N}}{2\tilde{\epsilon}_1 \tilde{a}^3 \tilde{H}}, \\ D_k &= a_*^3 H_* \sqrt{2\epsilon_{1*}} \left( \frac{\epsilon_{2*}}{2} \delta\phi_{k,h*} - \partial_N \delta\phi_{k,h*} \right). \end{aligned} \quad (2.9) $$

The linear power spectrum at the end of inflation is then

$$ \mathcal{P}_{\mathcal{R}}(k) = \frac{k^3 |C_k|^2}{2\pi^2}. \quad (2.10) $$

This relation, by construction, is exact when $C_k$ is the amplitude of the true growing mode. In section 3 we will investigate different approaches to determine the amplitude of the homogeneous parts of linear field perturbations during inflation.

<footnote n="1">¹The momentum constraint in perturbation theory relates the perturbation in the local lapse function to the local field perturbation, which is needed to relate the local proper time in the separate-universe patch to the global time coordinate being used in eq. (2.6) [34, 57, 58].</footnote>

<footnote n="2">²In this paper $\mathcal{R}$ and $\delta\phi$ will only be used to describe the linear free-field fluctuations.</footnote>
## 2.2 Classical $\delta\mathcal{N}$ formalism

To study large perturbations originating from quantum field fluctuations during inflation, we will need to go beyond linear perturbation theory. One such method is the $\delta\mathcal{N}$ formalism [11, 32–35, 56]. Here the curvature perturbation on uniform density hypersurface, is identified with the difference, $\delta\mathcal{N}$, between the local integrated expansion from an initial spatially-flat hypersurface to a final uniform-density hypersurface at the end of inflation, and the background expansion. In linear theory, we can identify this with the linear curvature perturbation at the end of inflation, $\mathcal{R} = H\delta\phi/\dot{\phi} = -\delta\mathcal{N}_{\text{lin}}$, where the end of inflation hypersurface corresponds to a comoving (and hence uniform-field) hypersurface. Here, recall that $\mathcal{R}$ and $\delta\phi$ denote linear fluctuations, and we have introduced the number of e-folds until the end of inflation, $\mathcal{N}$.<fnref n="3" /> As we have seen, we can identify $\mathcal{R}$ with the amplitude, $C$, of the growing mode of the homogeneous solution (2.7) using the long-wavelength limit of linear perturbation theory. More generally, we can calculate the local, non-linear expansion on large scales in the separate-universe approach using the local Klein-Gordon equation (2.1) and local Friedmann equation (2.2) [35, 56].

In classical, single-field inflation the local expansion is determined by the local field value and its time derivative on an initial spatially-flat hypersurface at time $N_*$ at a point $\vec{x}$ in real space, which gives the final non-linear curvature perturbation on uniform-density hypersurfaces [34, 61]

$$ \delta\mathcal{N} = \mathcal{N}(\phi_* + \delta\phi_{h*}, \dot{\phi}_* + \delta\dot{\phi}_{h*}) - \mathcal{N}(\phi_*, \dot{\phi}_*). \quad (2.12) $$

Here $\mathcal{N}(\phi, \dot{\phi})$ is the classical solution that gives the total number of e-folds from the phase-space position $(\phi, \dot{\phi})$ until the end of inflation. Hence we refer to this as the classical $\delta\mathcal{N}$ formalism. Since we are evaluating the expansion in the phase space for homogeneous cosmologies, we need to identify the homogeneous part, $\phi_* + \delta\phi_{h*}$, of the full inhomogeneous field, $\phi_* + \delta\phi_*$, on the initial hypersurface, as well as its proper time derivative, $\dot{\phi}_* + \delta\dot{\phi}_{h*}$.

The $\delta\mathcal{N}$ formalism is most often applied in cases where there is an attractor solution, reducing the 2D phase space to a 1D curve, $\dot{\phi} = \dot{\phi}(\phi)$. Consider the local homogeneous field value, $\phi_* + \delta\phi_{h*}$, at some time $N_*$ during inflation, on a spatially-flat hypersurface, which can be split into a background part and a perturbation. As there is a one-to-one relation between $\phi$ and $N$ on the background trajectory, the local homogeneous field perturbation, $\delta\phi_{h*}$, corresponds to a local adiabatic fluctuation, $\delta\mathcal{N}$, such that

$$ \phi_{\text{bg}}(N_*) + \delta\phi_{h*} = \phi_{\text{bg}}(N_* - \delta\mathcal{N}). \quad (2.13) $$

For monotonic $\phi_{\text{bg}}(N)$ and $\delta\phi_{h*} > 0$ the sign of $\delta\mathcal{N}$ in eq. (2.13) is positive if $\dot{\phi} < 0$ and negative if $\dot{\phi} > 0$. If we consider the field coarse-grained on some large scale, $k_S < aH$, then in linear theory the perturbation, $\delta\phi$, is due to the sum of independent fluctuations on large

$$ \mathcal{N}(N) = \int_{N}^{N_{\text{end}}} d\tilde{N}, \quad (2.11) $$

where we use a calligraphic font to distinguish this from the e-folding time, $N$, during inflation, since they go in opposite directions, $d\mathcal{N} = -dN$.

<footnote n="3">³The number of e-folds until the end of inflation is defined as</footnote>
scales $k_{in} < k < k_S$, where $k_S$ is the smoothing scale. Thus $\delta\phi$ is the sum of many Gaussian variables and is itself a Gaussian, and the linear curvature perturbation, $\mathcal{R} = H\delta\phi/\dot{\phi}$, is also a Gaussian. From the definition of the power spectrum, the coarse-grained curvature perturbation has a variance

$$ \sigma_{\mathcal{R}}^2(k_S) = \int_{k_{in}}^{k_S} \mathcal{P}_{\mathcal{R}}(k) \mathrm{d}\ln k. \qquad (2.14) $$

Since $\mathcal{R}$ is constant on super-Hubble scales for adiabatic perturbations [52], $\sigma_{\mathcal{R}}^2(k_S)$ has no explicit $N_*$ dependence, and conservation of probability gives the probability density function (PDF) for $\delta\mathcal{N}$ as

$$ P(\delta\mathcal{N}) = \frac{e^{-\frac{\mathcal{R}^2}{2\sigma_{\mathcal{R}}^2(k_S)}}}{\sqrt{2\pi}\sigma_{\mathcal{R}}(k_S)} \left| \frac{d\mathcal{R}}{d\delta\mathcal{N}} \right|. \qquad (2.15) $$

Therefore to apply the $\delta\mathcal{N}$ formalism, we need to know $\sigma_{\mathcal{R}}(k_S)$, the linear relation $\delta\phi = \dot{\phi}\mathcal{R}/H$, and the non-linear relation $\delta\mathcal{N} = \delta\mathcal{N}(\delta\phi)$ to determine $P(\delta\mathcal{N})$.

Let us consider, as an example, a phase of constant roll, i.e. constant $\epsilon_2$ [50], which will play an important role in our analysis in section 4. Integrating eq. (2.3) for $i=1$ in this case gives

$$ \epsilon_1 = \epsilon_{1*} e^{\epsilon_2(N-N_*)}. \qquad (2.16) $$

On the other hand, as mentioned below eq. (2.3), $\epsilon_1 = \dot{\phi}^2/(2H^2)$, hence

$$ \frac{d\phi}{dN} = \pm\sqrt{2\epsilon_1}. \qquad (2.17) $$

The choice of sign comes from taking the square root. We set it to match the dynamics of the background, and for explicitness let us consider the case where $\phi$ decreases with time classically, which corresponds to a negative sign. Integrating eq. (2.17) using eq. (2.16) gives the background solution

$$ \phi_{\mathrm{bg}}(N) = \phi_* + \frac{2}{\epsilon_2} \sqrt{2\epsilon_1} \left[ 1 - e^{\frac{\epsilon_2}{2}(N-N_*)} \right]. \qquad (2.18) $$

From eqs. (2.13) and (2.18), one then obtains$^4$

$$ \mathcal{R} = -\frac{\delta\phi_{h*}}{\sqrt{2\epsilon_{1*}}} = -\frac{2}{\epsilon_2} \left(1 - e^{-\frac{\epsilon_2}{2}\delta\mathcal{N}}\right). \qquad (2.19) $$

We can now find $P(\delta\mathcal{N})$ using eq. (2.15), giving [47, 50, 62]

$$ P(\delta\mathcal{N}) = \frac{1}{\sqrt{2\pi}\sigma_{\mathcal{R}}(k_S)} \exp \left[ -\frac{2}{\sigma_{\mathcal{R}}^2(k_S)\epsilon_2^2} \left(1 - e^{-\frac{\epsilon_2}{2}\delta\mathcal{N}}\right)^2 - \frac{\epsilon_2}{2}\delta\mathcal{N} \right], \qquad (2.20) $$

where $\sigma_{\mathcal{R}}(k_S)$ is given by eq. (2.14). The result has no explicit dependence on $N_*$, except through the value of $k_S$.

By applying the classical $\delta\mathcal{N}$ formalism at a point $\vec{x}$, including fluctuations across a range of comoving scales, we can compare it against the stochastic formalism, the main interest of our investigation.

<footnote>$^4$A “bg” subscript has been omitted from $\epsilon_1$, and higher derivatives, for brevity, as there is no possibility for confusion with another variable here.</footnote>
## 2.3 Stochastic formalism

Stochastic inflation is an effective-field-theory approach to modelling inflationary dynamics, where the full quantum-field-theory problem is mapped into a stochastic description. This is done by splitting the full quantum field, $\hat{\phi}$, and its momentum, $\hat{\pi}$, into short- and long-wavelength components, $\hat{\phi} = \phi + \hat{\phi}_s$ and $\hat{\pi} = \pi + \hat{\pi}_s$, where “short” and “long” are defined with respect to a coarse-graining scale $k_\sigma = \sigma aH$. Short-wavelength modes originate as quantum vacuum fluctuations but they impart a stochastic kick to the long-wavelength field as they continuously cross the coarse-graining scale during inflation. Crucially, we assume that the Hubble damping drives the modes into a squeezed state on long wavelengths, which means that the quantum fluctuations can be treated as a classical random field [63, 64]. This allows the backreaction of perturbations to be included in the classical dynamics on long wavelengths. We will model the long-wavelength dynamics using the separate-universe approach, which allows us to use the same non-linear evolution equations as in an unperturbed, background cosmology, apart from the inclusion of stochastic noise.

The stochastic equations of motion for single-field inflation are

$$ 
\begin{aligned}
\frac{\partial \phi}{\partial N} &= \pi_{\phi} + \xi_{\phi}, \\
\frac{\partial \pi_{\phi}}{\partial N} &= -\left(3 - \frac{\pi_{\phi}^2}{2}\right) \left[\pi_{\phi} + \frac{1}{V(\phi)} \frac{dV(\phi)}{d\phi}\right] + \xi_{\pi},
\end{aligned}
\quad (2.21) $$

where $\xi_\phi$ and $\xi_\pi$ are stochastic source terms, and $\phi$ and $\pi$ are now random variables. The covariance of the noise is identified with the quantum vacuum expectation value of the fields they stem from, which gives

$$ 
\begin{aligned}
\langle \xi_{\phi}(N_1) \xi_{\phi}(N_2) \rangle &= \frac{k_{\sigma}^3}{2\pi^2} |\delta\phi_{k_{\sigma},h}(N_1)|^2 \delta_D(N_1 - N_2), \\
\langle \xi_{\pi}(N_1) \xi_{\pi}(N_2) \rangle &= \frac{k_{\sigma}^3}{2\pi^2} |\delta\pi_{k_{\sigma},h}(N_1)|^2 \delta_D(N_1 - N_2), \\
\langle \xi_{\phi}(N_1) \xi_{\pi}(N_2) \rangle &= \frac{k_{\sigma}^3}{2\pi^2} \mathrm{Re} \left[ \delta\phi_{k_{\sigma},h}(N_1) \delta\pi_{k_{\sigma},h}^*(N_1) \right] \delta_D(N_1 - N_2).
\end{aligned}
\quad (2.22) $$

Here, we have considered the case where coarse graining is performed via a step function in Fourier space, with $\epsilon_1 \ll 1$, and $\delta_D$ is the Dirac distribution. A “*” superscript denotes the complex conjugate, and angle brackets $\langle \cdot \rangle$ are used to denote an ensemble average. Since the separate-universe approach is being used in the stochastic equations (2.21), $\delta\phi_{k_\sigma,h}$ and $\delta\pi_{k_\sigma,h}$ in eq. (2.22) correspond to the homogeneous part of the fluctuations [52].

The stochastic system of equations (2.21) can be written in a matrix form. We introduce the covariance matrix

$$ \Xi = \frac{k_{\sigma}^3}{2\pi^2} \begin{pmatrix} |\delta\phi_{k_{\sigma},h}|^2 & \mathrm{Re} \left[ \delta\phi_{k_{\sigma},h} \delta\pi_{k_{\sigma},h}^* \right] \\ \mathrm{Re} \left[ \delta\phi_{k_{\sigma},h} \delta\pi_{k_{\sigma},h}^* \right] & |\delta\pi_{k_{\sigma},h}|^2 \end{pmatrix}. \quad (2.23) $$

If $\Xi$ is non-singular, then its square-root, $S$, is real and symmetric. We can then rewrite eq. (2.21) as

$$ \frac{d}{dN} \begin{pmatrix} \phi \\ \pi_{\phi} \end{pmatrix} = \begin{pmatrix} \pi_{\phi} \\ -\left(3 - \frac{\pi_{\phi}^2}{2}\right) \left[ \pi_{\phi} + \frac{1}{V(\phi)} \frac{dV(\phi)}{d\phi} \right] \end{pmatrix} + \begin{pmatrix} S_{\phi\phi} & S_{\phi\pi} \\ S_{\phi\pi} & S_{\pi\pi} \end{pmatrix} \begin{pmatrix} \xi_1 \\ \xi_2 \end{pmatrix}, \quad (2.24) $$
where $S^2 = \Xi$, with $\xi_1$ and $\xi_2$ being two independent, normally-distributed random variables with unit variance.

The covariance matrix (2.23) defines a quadratic form, the contours of which are ellipses in phase space. Their properties determine how far the stochastic system is from an effective one-dimensional noise. The area of the ellipses is related to the determinant of $\Xi$, which is proportional to the squared Wrońskiskian of the mode functions, which itself is set by the canonical commutator between the field and its momentum at the quantum level. Working with $\phi$ and $\pi_\phi$ as the phase-space coordinates, the area thus decays exponentially as $a^{-6}H^{-2}$. The eccentricity of the ellipses is related to the amount of quantum squeezing. As the mode gets more and more squeezed, the eccentricity approaches one. In this limit the ellipse collapses to its semi-major axis. The angle made between this semi-major axis and the $\phi$ direction is given by

$$ \theta_n \equiv \arctan \left( \frac{\lambda_+ - \Xi_{\phi\phi}}{\Xi_{\pi\phi}} \right), \qquad (2.25) $$

where $\lambda_+$ is the positive eigenvalue of $\Xi$.

If $\delta\phi_{k_\sigma,h}$ is dominated by a single mode in eq. (2.7), then the determinant of $\Xi$ exactly vanishes, and the ellipses are infinitely squeezed. There is then only one degree of freedom in the noise, as $\delta\phi_{k_\sigma,h}$ and $\delta\pi_{k_\sigma,h}$ are fully correlated. The 2D noise of the stochastic system reduces to 1D, although the phase space of the local field remains 2D in general. The angle of the noise then reduces to

$$ \theta_n \equiv \arctan \left( \frac{\delta\pi_{k,h}}{\delta\phi_{k,h}} \right). \qquad (2.26) $$

If $\tan\theta_n = \tan\theta_{bg}$ as well, where the background trajectory is at an angle $\tan\theta_{bg} \equiv \partial_N\pi_\phi/\partial_N\phi = \epsilon_2/2$, then the noise is aligned with the background and the fluctuations are adiabatic. In this case the system given in eq. (2.21) is fully described by a 1D Langevin system.

The non-linear curvature perturbation (2.12) can be found using the $\delta\mathcal{N}$ formalism in stochastic inflation [36]. Although $N$ is used as the time variable in (2.24), by solving the stochastic system of equations for the field, the local integrated expansion from an initial spatially-flat hypersurface to the final uniform-density hypersurface at the end of inflation, $\mathcal{N}$, becomes a stochastic variable. Finding $\mathcal{N}$ and its associated PDF, $P(\mathcal{N})$, then becomes a first-passage-time problem [36]. The full non-linear curvature perturbation of a coarse-grained patch at the end of inflation is thus given by

$$ \delta\mathcal{N} = \mathcal{N} - \langle\mathcal{N}\rangle, \qquad (2.27) $$

with an associated PDF, $P(\delta\mathcal{N})$. While $P(\delta\mathcal{N})$ can be found analytically in certain cases [29, 37–39, 65], in general it needs to be found numerically. It has however been shown that in the presence of quantum diffusion we expect the far tail of $P(\delta\mathcal{N})$ to be a simple exponential [29], potentially giving a significant enhancement of very large curvature fluctuations compared to the approximately Gaussian tail found in perturbation theory.
# 3 Noise models

## 3.1 Consistency criterion

Here we will detail how the homogeneous parts of the fluctuations on a given scale, $\delta\phi_{k_\sigma,h}$ and $\delta\pi_{k_\sigma,h}$ in eq. (2.22), can be calculated. There are a number of approaches one can take even in the context of linear perturbation theory, but a consistency requirement of any approach is that it gives the correct linear power spectrum at the end of inflation [52]. This is to make sure the correct constants of integration $C_k$ and $D_k$ have been chosen in the homogeneous solution (2.7) for each wavenumber $k$. Therefore we require that the fractional difference between the power spectrum at the end of inflation, eq. (2.10), found by identifying the homogeneous part of the perturbations at $N_\sigma$, and the power spectrum (2.8) found by numerically solving the full (inhomogeneous) mode equation (2.4), is less than 1%:

$$ \left| \frac{\mathcal{P}_{\mathcal{R}}(\delta\phi_{k_\sigma,h}, \partial_N\delta\phi_{k_\sigma,h}, N_\sigma) - \mathcal{P}_{\mathcal{R}}(k, N_{\text{end}})}{\mathcal{P}_{\mathcal{R}}(k, N_{\text{end}})} \right| \le 0.01, \qquad (3.1) $$

for all of the modes included in the numerical simulations. Here $N_\sigma$ is the matching time corresponding to $k = \sigma aH$. For the two potentials considered in this paper, the above criterion is shown in appendix A for the noise model used in the simulations, as will be detailed below.

## 3.2 Analytical approximations

The simplest approach, commonly adopted in the $\delta N$ formalism [29, 36, 37, 44, 66], is to approximate the noise as that for a light field in a quasi-de Sitter background,<fnref n="5" /> corresponding to slow-roll inflation with $|\epsilon_i| \ll 1 \forall i \ge 1$. The general solution for the homogeneous part of the field perturbation in the spatially-flat gauge can then be readily identified as [52]

$$ \delta\phi_{k,h} = \frac{iH}{\sqrt{2k^3}} \left[ (\alpha_k - \beta_k) - i(\alpha_k + \beta_k) \frac{(k\eta)^3}{3} \right]. \qquad (3.2) $$

On super-Hubble scales the decaying mode, $\sim (k\eta)^3$, exponentially decays, and is dropped. In particular, if Bunch-Davies initial conditions ($\alpha_k = 1$ and $\beta_k = 0$) are used, then the covariance matrix (2.23) is

$$ \Xi = \frac{H^2}{(2\pi)^2} \begin{pmatrix} 1 & 0 \\ 0 & 0 \end{pmatrix}. \qquad (3.3) $$

Note that in this simple limit the homogeneous field perturbation is constant and hence the noise is independent of the choice of coarse-graining scale, $k_\sigma$. As the noise is only in the $\phi$ direction, only a single random variable is needed to characterise the noise at each time step. The consistency relation (3.1) is obeyed for slow-roll models with $|\epsilon_2| \ll 1$.

This approach can be further refined by requiring $\epsilon_1 \ll 1$ but allowing $|\epsilon_i| \ne 0 \forall i \ge 2$. We will require

$$ \nu^2 = \frac{9}{4} - \epsilon_1 + \frac{3}{2}\epsilon_2 - \frac{1}{2}\epsilon_1\epsilon_2 + \frac{1}{4}\epsilon_2^2 + \frac{1}{2}\epsilon_2\epsilon_3 \qquad (3.4) $$

<footnote n="5">⁵In the de Sitter limit, $\epsilon_1 \to 0$, the field perturbations decouple from the metric and we effectively have a test field in de Sitter spacetime.</footnote>
to be approximately constant. This gives an analytic approximation for the general solution of the mode equation eq. (2.4) of the form [67]

$$ \delta\phi_k = \frac{\sqrt{-\eta}}{a} [A_k J_\nu(-k\eta) + B_k Y_\nu(-k\eta)], \qquad (3.5) $$

where $J_\nu$ and $Y_\nu$ are Bessel functions of the first and second kind. If we again assume Bunch-Davies initial conditions (2.5), then the Bessel coefficients are given by

$$ A_k = \frac{\sqrt{\pi}}{2} e^{i(\nu + \frac{1}{2})\frac{\pi}{2}}, \quad B_k = \frac{\sqrt{\pi}}{2} e^{i(\nu + \frac{3}{2})\frac{\pi}{2}}. \qquad (3.6) $$

The homogeneous component is defined as the part of eq. (3.5) which is a solution of eq. (2.6) and is found by taking the leading order term in a $(-k\eta)$ expansion of $J_\nu$ and $Y_\nu$. The homogeneous growing mode is thus given by

$$ \delta\phi_{k,h} = e^{i(\nu - \frac{1}{2})\frac{\pi}{2}} \frac{H}{2^{\frac{3}{2}-\nu}} \frac{(-k\eta)^{\frac{3}{2}-\nu}}{\sqrt{2}k^3}. \qquad (3.7) $$

The associated correlation matrix (2.23), which is evaluated at $-k\eta = \sigma$, is then

$$ \Xi = \frac{H^2}{4\pi^2} \left[ \frac{\Gamma(\nu)}{\Gamma(\frac{3}{2})} \right]^2 \left(\frac{\sigma}{2}\right)^{3-2\nu} \begin{pmatrix} 1 & \nu - \frac{3}{2} \\ \nu - \frac{3}{2} & (\nu - \frac{3}{2})^2 \end{pmatrix}. \qquad (3.8) $$

While there is noise in both the $\phi$ and $\pi_\phi$ directions, the determinant of this matrix is exactly zero, hence there is only one degree of freedom and the stochastic system is described by a single (1D) noise term, as was discussed near eq. (2.26). One can verify that the consistency relation (3.1) is satisfied for slow-roll inflation, but we will also show that the consistency relation is satisfied in non-slow-roll models after a sudden transition, including, but not limited to, constant-roll regimes. Note that for $\nu \to 3/2$, we recover the result (3.3) for a massless field in de Sitter and the $\sigma$ dependence is lost.

## 3.3 Numerical matching

The analytic approximations used to estimate the noise in section 3.2 rely on $\nu^2$ in eq. (3.4) being constant. This will not be valid when the inflaton experiences a sudden transition. Note that, during sudden transitions, not only scales around and above the Hubble radius are affected, sub-Hubble modes are also altered, and this implies that when they later cross the Hubble radius they may do so with “initial” conditions that depart from the Bunch-Davies profile.

The most general approach to determine the noise is to match the general homogeneous solution (2.7) to the full numerical solution for a given scale, $\delta\phi_k$ and $\partial_N\delta\phi_k$, on sufficiently long wavelengths. The covariance matrix is then

$$ \Xi = \frac{k_{\sigma}^3}{2\pi^2} \left( \begin{array}{c|c} |\delta\phi_{k_{\sigma}}|^2 & \Re[\delta\phi_{k_{\sigma}} \delta\pi_{k_{\sigma}}^*] \\ \hline \Re[\delta\phi_{k_{\sigma}} \delta\pi_{k_{\sigma}}^*] & |\delta\pi_{k_{\sigma}}|^2 \end{array} \right) \Big|_{N=N_{\sigma}}, \qquad (3.9) $$

with $\sigma$ sufficiently small such that the behaviour of $\delta\phi_{k_\sigma}$ has minimal contribution from gradient, $O(k^2)$ terms. This approach can be used for any background dynamics and can
account for non-Bunch-Davies initial conditions but is limited to scales significantly larger than the Hubble scale, $\sigma \ll 1$. As eq. (3.9) has a non-zero determinant, the associated phase space is 2D and two random variables are required to describe the system.

While general, this approach must either accept a finite $k^2$ contamination or use a very small $\sigma$ [40, 41, 51, 68–70], potentially missing non-linear effects associated with stochastic inflation [52]. We find that the consistency relation (3.1) is obeyed, with $C_k$ found using eq. (2.9) with $\delta\phi_{k,h*}$ and $\delta\pi_{k,h*}$ replaced by $\delta\phi_{k*}$ and $\delta\pi_{k*}$, when $\sigma \le 0.01$ in the absence of a sudden transition. An even smaller value is required to be accurate for all modes through a sudden transition.

## 3.4 Semi-analytical approach

The third approach is to try and combine the advantages of the previous two. Namely being able to directly identify the homogeneous component with the analytical approach, and allowing for non-Bunch-Davies initial conditions with the numerical approach, effectively using the Bessel function (3.5) as an ansatz.

This is done by assuming that $\nu$ is locally constant in eq. (2.4), and continuously matching the numerical solution to eq. (3.5) to find the Bessel coefficients

$$ 
\begin{aligned}
A_k &= \frac{a}{\sqrt{-\eta}} \frac{\delta\phi_k \partial_N [Y_\nu(-k\eta)] - \left(\frac{3}{2}\delta\phi_k - \partial_N \delta\phi_k\right) Y_\nu(-k\eta)}{J_\nu(-k\eta) \partial_N [Y_\nu(-k\eta)] - \partial_N [J_\nu(-k\eta)] Y_\nu(-k\eta)} \Big|_{N=N_\sigma}, \\
B_k &= -\frac{a}{\sqrt{-\eta}} \frac{\delta\phi_k \partial_N [J_\nu(-k\eta)] - \left(\frac{3}{2}\delta\phi_k - \partial_N \delta\phi_k\right) J_\nu(-k\eta)}{J_\nu(-k\eta) \partial_N [Y_\nu(-k\eta)] - \partial_N [J_\nu(-k\eta)] Y_\nu(-k\eta)} \Big|_{N=N_\sigma}.
\end{aligned}
\quad (3.10) $$

This method allows for possible time dependence of $\nu$ before Hubble exit, hence for non-trivial sub-Hubble dynamics during the transition.

The homogeneous behaviour can then be identified as before, using the leading terms in a $(-k\eta)$ expansion of eq. (3.5). One then finds

$$ \delta\phi_{k,h} = \frac{\sqrt{-\eta}}{a} \left[ -\frac{B_k}{\pi} \Gamma(\nu) + \frac{A_k \tan(\pi\nu) + B_k}{\Gamma(\nu+1)\tan(\pi\nu)} \left(-\frac{k\eta}{2}\right)^{2\nu} \right] \left(-\frac{k\eta}{2}\right)^{-\nu}, \quad (3.11) $$

where the leading term is the homogeneous growing mode and the second term is the homogeneous decaying mode, as $\eta \to 0$ at late times and $\nu > 0$. Here $\epsilon_1 \ll 1$ was assumed, and $\nu$ is evaluated at the time $k$ crosses the coarse-graining scale (hence it carries an implicit $k$ dependence). If the growing mode dominates over the decaying mode, the covariance matrix is then given by

$$ \Xi = \frac{H^2}{4\pi^2} \frac{4|B_k|^2}{\pi} \left[ \frac{\Gamma(\nu)}{\Gamma(\frac{3}{2})} \right]^2 \left( \frac{\sigma}{2} \right)^{3-2\nu} \begin{pmatrix} 1 & \nu - \frac{3}{2} \\ \nu - \frac{3}{2} & (\nu - \frac{3}{2})^2 \end{pmatrix}. \quad (3.12) $$

When $B_k$ is given by its Bunch-Davies value (3.6), then the above coincides with eq. (3.8). In this limit the determinant of the covariant matrix vanishes, and the stochastic system is subject to a 1D noise term [see again the discussion near eq. (2.26)].

The homogeneous growing mode found in eq. (3.12) obeys the consistency relation (3.1) as long as $\nu$ is approximately constant as the coarse-graining scale is crossed, see appendix A
<figure><img src="image_3.png" /><figcaption>**Figure 1.** The phase-space behaviour for the Gaussian bump potential (4.12), for field values near the sudden transition. Here $\phi_T$ is the field value at the first time when $\epsilon_2 < -3$ on the reference trajectory which starts on the slow-roll attractor at early times. The colour bar shows the value of $\nu^2$ defined in eq. (3.4), at each point in the phase space.</figcaption></figure>

for more details. In practice, we define the period where the constant-$\nu$ approximation is valid by requiring that the consistency relation (3.1) holds.

We expect a phase of approximately constant $\nu$ after a sudden transition for most single-field models of inflation that produce a measurable abundance of PBHs. This can be seen by using the Klein-Gordon equation (2.1) and the Hubble-flow parameters (2.3) to rewrite eq. (3.4) as

$$ \nu^2 = \frac{9}{4} - \epsilon_1 - \frac{1}{2} \frac{\epsilon_1 \epsilon_2^2}{(3 - \epsilon_1)} + (\epsilon_1 - 3) \left[ \frac{1}{V} \frac{d^2 V}{d\phi^2} - \left( \frac{1}{V} \frac{dV}{d\phi} \right)^2 \right]. \quad (3.13) $$

Note that this expression is exact. Post-transition, to grow the power spectrum sufficiently to produce a measurable abundance of PBHs, we have $\epsilon_1 \ll 1$. Then taking eq. (3.13) to zeroth-order in $\epsilon_1$, while requiring$^6 \epsilon_1 \epsilon_2^2 \ll 1$, one obtains

$$ \nu^2 \approx \frac{9}{4} - 3 \left[ \frac{1}{V} \frac{d^2 V}{d\phi^2} - \left( \frac{1}{V} \frac{dV}{d\phi} \right)^2 \right]. \quad (3.14) $$

Therefore post-transition $\nu^2$ is only a function of $\phi$. This can be checked in figure 1, where the colour gradient encoding variations of $\nu^2$ is aligned with the $\phi$ direction, regardless of the phase-space trajectory. In fact, the maximum difference between eq. (3.14) and $\nu^2$ for any of the trajectories shown is 0.05.

Taking the $N$ derivative of eq. (3.14) then gives

$$ |\partial_N \nu^2| \approx \sqrt{18\epsilon_1} \left| \frac{1}{V} \frac{d^3 V}{d\phi^3} - \frac{3}{V^2} \frac{dV}{d\phi} \frac{d^2 V}{d\phi^2} + \left( \frac{1}{V} \frac{dV}{d\phi} \right)^3 \right| . \quad (3.15) $$

<footnote>$^6$In single-field models of interest, the maximum value of $|\epsilon_2|$ obtainable is $\mathcal{O}(10)$. Therefore $\epsilon_1 \ll 10^{-2}$ is required for $\epsilon_1 \epsilon_2^2 \ll 1$, which is easily obtained for most PBH producing models.</footnote>
In the limit where $\epsilon_1 \ll 1$, we can take $\sqrt{18\epsilon_1} < 1$. In the models we will consider, apart from at a sudden transition, the potential derivative terms in the square brackets are at a maximum of order unity. Therefore in the post-transition phase, $|\partial_N \nu^2| < 1$ and $\nu$ is then approximately constant. The duration of this constant $\nu$ phase depends on the potential and the minimum value of $\epsilon_1$ reached.

# 4 Stochastic $\delta N$ simulations and results

In this section we will apply the techniques developed in section 3 to describe the quantum noise used to calculate the full probability distribution for the first-passage time in the stochastic $\delta N$ formalism. We compare this to predictions using the classical $\delta N$ approach in models involving a sudden transition from an initial slow-roll inflation attractor, through an intermediate ultra-slow-roll phase (where $\epsilon_2 \sim -6$), before returning to an attractor solution at late times.

## 4.1 Numerical approach

We simulate the Langevin system given in eq. (2.24) using Euler-Maruyama steps [71]. The $(i+1)^{th}$ step is given by

$$ \phi^{(i+1)} = \phi^{(i)} + \pi^{(i)} \Delta N + \left( S_{\phi\phi}^{(i)} \xi_1^{(i)} + S_{\phi\pi}^{(i)} \xi_2^{(i)} \right) \sqrt{\Delta N}, \\ \pi^{(i+1)} = \pi^{(i)} - \left[ 3 - \frac{(\pi^{(i)})^2}{2} \right] \left[ \pi^{(i)} + \frac{1}{V(\phi^{(i)})} \frac{dV(\phi^{(i)})}{d\phi} \right] \Delta N \quad (4.1) \\ + \left( S_{\pi\phi}^{(i)} \xi_1^{(i)} + S_{\pi\pi}^{(i)} \xi_2^{(i)} \right) \sqrt{\Delta N}, $$

where $\xi_1^{(i)}$ and $\xi_2^{(i)}$ are independent random numbers drawn from a Gaussian distribution and $\Delta N$ is the time step. The system is initiated at $(\phi_{init}, \pi_{init})$. It is then propagated using the above Euler-Maruyama steps until $\epsilon_1(\phi^{(n)}, \pi^{(n)})$ crosses one and the first-passage time $N = n\Delta N$ is recorded, where $n$ is the number of steps required. By running many such simulation runs, $P(\delta N)$ can be sampled numerically.

In principle, the noise covariance matrix, $\Xi$ in eq. (2.23), for each simulation should be calculated for the quantum field within that specific stochastic realisation. This can be done in two ways: either using the exact numerical-matching approach presented in section 3.3, or using the semi-analytical one discussed in section 3.4. While the former greatly increases the computational complexity, requiring supercomputers to solve with current techniques [40, 41], the latter is more straightforward to implement since it only requires to compute the local value of $\nu$ in the background at each time step (the main limitation being that $\nu$ needs to be approximately constant). Although one of these methods would be required for a detailed investigation of the far tail, in what follows we rather evaluate the noise using the unperturbed, classical background. The reason is that, since the main goal of this work is to design a robust modeling of the noise, and to generalise importance-sampling methods to multi-dimensional phase-spaces, it is important for us to be able to compare our results with the classical-$\delta N$ formalism in order to test their validity. This is why we solve stochastic equations in a setup that is otherwise classical: the noise is pre-computed along a reference,
classical trajectory<fnref n="7" /> and after a fixed number of time steps, corresponding to a chosen range of comoving wavenumbers, the noise is turned off and the time to reach the end of inflation classically is recorded. With this approach, each simulation set presented used 10⁵ numerical runs, taking less than 15 minutes on an 8 core 1.6 GHz processor. For more details about the computation and open source code, see the PyFPT (https://github.com/Jacks0nJ/PyFPT) Github repository. The main speedup comes from the technique of importance sampling [43], which we have extended to the 2D phase space.

To probe the rare and large $\delta N$ events associated with PBH production efficiently, importance sampling is utilised ref. [43]. Importance sampling introduces a bias term, $\vec{B}$, at each numerical step in eq. (4.1) to increase the frequency of large first-passage-time events being realised. The induced change to the relative probability for each step is then recorded, with the total relative probability of the biased path compared to the unbiased path being given by the weight, $w$. For an unbiased multi-dimensional Langevin equation, the $(i+1)^{th}$ step is given by

$$ \vec{x}^{(i+1)} - \vec{x}^{(i)} = \vec{D}(\vec{x}^{(i)}, N^{(i)})\Delta N + S(\vec{x}^{(i)}, N^{(i)})\vec{\xi}^{(i)}\sqrt{\Delta N}, \quad (4.2) $$

with drift vector $\vec{D}$ and noise matrix $S$. The biased step is

$$ \vec{x}^{(i+1)} - \vec{x}^{(i)} = \vec{D}(\vec{x}^{(i)}, N^{(i)})\Delta N + \vec{B}(\vec{x}^{(i)}, N^{(i)})\Delta N + S(\vec{x}^{(i)}, N^{(i)})\vec{\xi}^{(i)}\sqrt{\Delta N}. \quad (4.3) $$

The weight of the whole stochastic process for $n$ steps is

$$ w = \exp \left( - \sum_{i=0}^{n-1} \Delta A^{(i)} \right), \quad (4.4) $$

where [42]

$$ \Delta A^{(i)} = \left( \frac{\vec{B}^{(i)}}{2} \Delta N + S^{(i)} \vec{\xi}^{(i)} \sqrt{\Delta N} \right)^{\mathrm{T}} (\Xi^{(i)})^{-1} \vec{B}^{(i)}, \quad (4.5) $$

and $\Xi = S^2$ by definition. PyFPT is able to apply importance sampling to a 2D Langevin system, as long as the correlation matrix is non-singular.

The error in the reconstruction of $P(\delta N)$ is found numerically. When direct sampling is used, corresponding to $\vec{B} = 0$, as each simulation run is independent, jackknife resampling is used. When importance sampling is used, the error in $P(\delta N)$ is found from the error in the estimate of the mean of a lognormal distribution. See ref. [43] for details of both methods.

As already mentioned, when the growing mode dominates over the decaying mode, the noise becomes effectively 1D. In appendix A we check that this a valid approximation even at Hubble crossing with $\sigma = 1$, as long as $\nu$ is approximately constant. Since this condition is also required for the semi-analytical noise description of section 3.4 to apply, we will focus on this regime from now on. In practice, this means that we will start our simulations after a sudden transition to ultra-slow roll, once $\nu^2$ is approximately constant, although PyFPT is applicable beyond this so long as the separate-universe approach applies.

<footnote n="7">⁷In eqs. (2.21)–(2.24) the local time in terms of the number of e-folds, $N$, is used throughout. However, by limiting ourselves to using the pre-computed noise, in practice our simulations are run with the global time $N$.</footnote>
The implementation of importance sampling is more straightforward in 1D. If the noise is given solely by a single growing mode, then the system in general can be written as

$$ 
\begin{aligned}
\frac{d\phi}{dN} &= \pi + S_{\phi\phi}\hat{\xi}, \\
\frac{d\pi}{dN} &= -\left(3 - \frac{\pi^2}{2}\right) \left[\pi + \frac{1}{V(\phi)}\frac{dV(\phi)}{d\phi}\right] + \tan(\theta_n) S_{\phi\phi}\hat{\xi},
\end{aligned}
\qquad (4.6) $$

where $\hat{\xi}$ is a random variable with unit variance and $\tan(\theta_n)$ is given in eq. (2.26). Therefore the system is 2D but with a 1D noise. If $\tan(\theta_n) = \partial_N\pi/\partial_N\phi \equiv \epsilon_2/2$ throughout, then the noise is aligned with the unperturbed background trajectory and the phase space is also 1D, e.g. during slow roll or constant roll. However this is not in general the case in the absence of a dynamical attractor, e.g. during ultra-slow roll, and in this case the full 2D phase space needs to be used to account for stochastic kicks misaligned with the unperturbed background trajectory.

Applying a bias to eq. (4.6) gives for the $(i+1)^{th}$ step

$$ 
\begin{aligned}
\Delta\phi^{(i+1)} &= \pi^{(i)}\Delta N + \mathcal{B}_{\phi}^{(i)}\Delta N + S_{\phi\phi}^{(i)}\hat{\xi}^{(i)}\sqrt{\Delta N}, \\
\Delta\pi^{(i+1)} &= \left(3 - \frac{(\pi^{(i)})^2}{2}\right) \left[\pi^{(i)} + \frac{1}{V(\phi^{(i)})}\frac{dV(\phi^{(i)})}{d\phi}\right]\Delta N + \mathcal{B}_{\pi}^{(i)}\Delta N \\
&\quad + \tan(\theta_n^{(i)}) S_{\phi\phi}^{(i)}\hat{\xi}^{(i)}\sqrt{\Delta N}.
\end{aligned}
\qquad (4.7) $$

As there is only a single noise term, eq. (4.5) does not apply directly, as the determinant of $\Xi = S^2$ is zero and it is thus singular. Instead, the relative probability of the importance sampling step must be the same along each dimension. Therefore we apply the 1D version of eq. (4.5) for both $\phi$ and $\pi$ and equate them, giving the following constraint

$$ 
\begin{aligned}
\frac{\mathcal{B}_{\phi}^{(i)}}{(S_{\phi\phi}^{(i)})^2} \left( \frac{\mathcal{B}_{\phi}^{(i)}}{2} \Delta N + S_{\phi\phi}^{(i)} \xi^{(i)} \sqrt{\Delta N} \right) = \\
\frac{\mathcal{B}_{\pi}^{(i)}}{[\tan(\theta_n^{(i)}) S_{\phi\phi}^{(i)}]^2} \left[ \frac{\mathcal{B}_{\pi}^{(i)}}{2} \Delta N + \tan(\theta_n^{(i)}) S_{\phi\phi}^{(i)} \xi^{(i)} \sqrt{\Delta N} \right],
\end{aligned}
\qquad (4.8) $$

which reduces to

$$ \mathcal{B}_{\phi}^{(i)} = \mathcal{B}_{\pi}^{(i)} \tan(\theta_n^{(i)}). \qquad (4.9) $$

For diffusion-based bias with $\mathcal{B}_{\phi}^{(i)} = A S_{\phi\phi}^{(i)}$, where $A$ is known as the bias amplitude, one then finds

$$ \mathcal{B}_{\phi}^{(i)} = A \tan(\theta_n^{(i)}) S_{\phi\phi}^{(i)}. \qquad (4.10) $$

It is this form of the bias that will be used in the numerical results below.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Left: the Hubble-flow parameters (2.3) and $\nu^2$ (3.4) for the piece-wise linear potential defined in eq. (4.11), near the sudden transition at $N_T$. The quantity $\log_{10}(\epsilon_1) + 5$ is plotted to fit on the same vertical scale as $\epsilon_2$. Right: the associated power spectrum (2.8) where $k_T$ is the scale that crosses the Hubble radius at time $N_T$. The gray dotted vertical line corresponds to when ultra-slow roll ends, here defined as when $\epsilon_2$ crosses $-3$, and the black dashed lines bound the modes used in our stochastic simulations, $k_{min} < k < k_{max}$ (the vertical lines on the left panel label the times at which the scales labeled by the same lines on the right panel cross the Hubble radius).</figcaption></figure>

## 4.2 Piece-wise linear potential

Let us first consider one of the simplest models with a sudden transition from slow-roll to ultra-slow-roll inflation as a test case for PYFPT, the piece-wise linear potential [72]. This toy model has many well-known analytical results [69, 73–76]. The potential is

$$V(\phi) = \begin{cases} V_0 + A_+(\phi - \phi_T) & \text{for } \phi \ge \phi_T, \\ V_0 + A_-(\phi - \phi_T) & \text{for } \phi < \phi_T, \end{cases} \qquad (4.11)$$

where $A_+ > A_- > 0$ and we choose $\phi_{init} \gg \phi_T$ such that the slow-roll attractor, $\dot{\phi} = -A_+/3H$, has been reached before the transition. Immediately after the transition, there is a mismatch between the scalar-field velocity and the potential gradient, $A_-$, resulting in slow roll being broken and giving a transitory phase of ultra-slow roll (as shown in the left-hand panel of figure 2). Eventually the field relaxes back onto a new slow-roll trajectory with $\dot{\phi} = -A_-/3H$. In the following we take $V_0 = 2.8719 \times 10^{-12}$ and $A_+ = 10^{-14}$ to have a scale-invariant power spectrum with an amplitude in agreement with CMB observations on large scales [7], and $A_- = 10^{-3}A_+$, such that the power spectrum peaks at $P_R \sim 5 \times 10^{-3}$. The slow-roll parameters (2.3) and the power spectrum (2.8) are shown in figure 2. The range of modes that left the Hubble radius after the transition, and whose homogeneous behaviour we can identify at Hubble exit using our semi-analytical approach (3.12), is displayed with the black dashed lines. As $\nu$ is approximately constant for the entire post-transition phase, the criterion in eq. (3.1) is always met in that range, see appendix A.

In the left-hand plot of figure 3, the ellipses corresponding to the covariance matrix (2.23) are shown at different coarse-graining scales. The scale shown leaves the Hubble scale immediately after the sudden transition and is soon dominated by the growing mode. This is why it enters a highly-squeezed state, as the eccentricity of the ellipses tends to one. This
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Left: evolution of the ellipse associated with the covariance matrix (2.23) with coarse-graining scale $\sigma$, for the piece-wise linear potential of section 4.2. The mode shown exited the Hubble radius immediately after the transition with $k = 1.1 k_{\text{T}}$. The areas of each ellipse is normalised to $\pi$ to allow visual comparison. Right: the tangent of the phase-space angle for the background, $\partial_N^2\phi/\partial_N\phi$, and for the 1D noise found using eq. (3.12) near the sudden transition at $N_{\text{T}}$. The dotted gray vertical line corresponds to the mode that left the Hubble radius when ultra-slow roll ended and the dashed black lines bound the modes used in the stochastic simulations.</figcaption></figure>

squeezing occurs even quicker for subsequent modes leaving the Hubble radius. Although this provides an illustrative justification for our use of a 1D treatment of the noise after the transition, the squeezing amplitude (i.e. the eccentricity of the ellipse) is not a symplectic invariant quantity and it can be changed by a mere field rescaling, or in general by any canonical transformation [64]. This is why, in practice, we use the criterion given in eq. (3.1) to determine when the growing mode dominates and the noise is effectively 1D.

In the right-hand plot of figure 3, the tangent of the phase-space angle of both the background dynamics, $\partial_N\pi/\partial_N\phi = \epsilon_2/2$, and the 1D noise (3.12) are shown. While at late times the two agree and the phase space is truly 1D, during ultra-slow roll there is a misalignment between the noise and the background, justifying the need to simulate the full 2D phase space.

In figure 4 we investigate the semi-analytical noise model of section 3.4 combined with the numerical-simulation strategy of section 4.1. Around the peak of the first-passage time distribution, see the left-hand plot, we run direct sampling stochastic simulations, either using the 1D noise (4.6) or the full 2D noise (2.24) approaches, for a range of $\sigma$ values. As $\sigma$ is varied, so are $\phi_{\text{init}}/\pi_{\text{init}}$ and $\phi_{\text{end}}$, to correspond to when the first and last mode in the range shown by the vertical dashed lines in figure 2 left the coarse-graining scale. The result is found to be mostly independent of $\sigma$, even for values up to $\sigma = 1$, and to agree with the expectation from linear perturbation theory displayed with the black solid curve. This suggests that for modes that left the Hubble radius post-transition, the resulting noise can be modelled using the growing mode at Hubble-crossing, and this validates the noise modelling proposed in section 3.4. This is also consistent with the fact that the 1D treatment of the noise gives an excellent approximation of the full 2D treatment, see the blue and brown data, in agreement with the discussion around the left panel of figure 3, as well as previous works which used $\sigma = 0.01$ [70, 77].
<figure><img src="image_6.png" /><figcaption>Figure 4. The numerically reconstructed first-passage-time distribution for the piece-wise linear potential (4.11). The semi-analytical noise model described in section 3.4 is used. Left: comparison of the 1D (4.6) and 2D (2.24) models. The value of $\sigma = k/(aH)$ defines when the noise is applied and the start/end field values. The solid black curve is Gaussian with variance given by eq. (2.14). Right: importance sampling using the 1D noise approach (4.7) to find the tail of the distribution for a range of bias amplitudes, given by the different values of $\mathcal{A}$ shown, for $\sigma = 1$. The numerical error bars correspond to one standard deviation, as discussed in section 4.1.</figcaption></figure>

The right-hand plot of figure 4 shows the results of applying importance sampling using the 1D noise approach (4.7) at Hubble crossing, with $\sigma = 1$. As the bias amplitude is increased and the numerics probe further into the tail, each data set overlaps with the previous one, including direct sampling with $\mathcal{A} = 0$, showing that importance sampling is self-consistent. This confirms our ability to probe the tail of the distribution with importance sampling.

Let us also notice that the Gaussian distribution obtained from standard perturbation theory, whose variance is given in eq. (2.14), provides a reliable fit to the whole distribution, both around the maximum and on the far tail. This is because, sufficiently long after the transition, there is a linear relationship between $N$ and $\phi$, and thus between $\delta N$ and $\delta\phi$. Therefore the classical $\delta N$ formalism (2.15) gives a Gaussian distribution, and we find the stochastic simulations not to deviate substantially from that expectation. Let us however stress that, since a fixed range of comoving modes is integrated over, the non-linear curvature perturbation $\delta N$ is coarse-grained at a scale that exits the Hubble radius in the late slow-roll phase, where the linear relationship between $N$ and $\phi$ holds.

## 4.3 Gaussian bump

Now let us consider a more realistic model with a smooth transition from slow roll to ultra-slow roll. One such model is made of a slow-roll potential with the addition of a Gaussian bump [70, 78],

$$V(\phi) = V_0 \frac{\phi^2}{m^2 + \phi^2} \left\{ 1 + K \exp \left[ -\frac{1}{2} \frac{(\phi - \phi_0)^2}{\Sigma^2} \right] \right\}, \quad (4.12)$$

where we use $K = 1.17 \times 10^{-3}$, $m = 0.5$, $\Sigma = 1.59 \times 10^{-2}$ and $\phi_0 = 2.18812$. These parameters are chosen such that CMB constraints [7] are obeyed and the power spectrum peaks at $\mathcal{P}_R \sim 5 \times 10^{-3}$ within the asteroid mass window [79]. In figure 5 the background dynamics and power spectrum are shown, along with the range of $k$ modes integrated over to
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Left: the Hubble-flow parameters (2.3) and $\nu^2$ defined in eq. (3.4) for the Gaussian-bump potential (4.12), near the sudden transition at $N_T$, corresponding to the start of ultra-slow roll with $\epsilon_2 < -3$. Right: the associated power spectrum (2.8) where $k_T$ is the scale that crosses the Hubble radius at time $N_T$. The dotted gray vertical lines frame the modes that left the Hubble radius during ultra-slow roll ($\epsilon_2 < -3$) and the dashed black lines bound the modes integrated over for the noise, as in figure 2.</figcaption></figure>

compute the noise, $k_{min}$ and $k_{max}$. This range is chosen such that the linear power spectrum consistency relation (3.1) is obeyed for $\sigma \in [0.01, 1]$ for the semi-analytical noise model given in eq. (3.12), see appendix A. While excluding half of the modes from the ultra-slow-roll phase, this range still includes >85% of the contribution of the linear power spectrum to $\sigma_R$ when calculated over the full $k$ range shown.

In the left-hand plot of figure 6 we investigate the semi-analytical noise model of section 3.4 combined with the numerical strategy of section 4.1. The variance of the distribution is insensitive to the choice of the $k_{max}$ used when it was sufficiently far away from the power spectrum peak mode. The peak of $P(\delta\mathcal{N})$ and the small near-tail non-Gaussianity are independent (within the error bars) both of the choice of $\sigma$ and whether the 2D or 1D noise model is used. Therefore the 1D importance sampling approach (4.7) is justified.

In the right-hand plot of figure 6 the far tail of $P(\delta\mathcal{N})$ is investigated for two different $\sigma$ choices using 1D importance sampling (4.7). For large $\delta\mathcal{N} \gtrsim 0.1$ there is a clear deviation from the Gaussian distribution with variance given by linear theory eq. (2.14). While the peaks of the distributions agree, the exponential-like far tails deviate. At $\delta\mathcal{N} = 1$, the $\sigma = 1$ and $\sigma = 0.01$ results for $P(\delta\mathcal{N})$ are enhanced by factors of $\sim 10^6$ and $\sim 10^8$ respectively, compared to the Gaussian prediction.

The classical $\delta\mathcal{N}$ prediction (2.20) is in agreement with the results of figure 6 if we choose the value of $\epsilon_2$ appropriately. As $\epsilon_2$ is not exactly constant during this period, some post-hoc fitting is required (the dashed lines). The value of $\epsilon_2$ needed for the $\sigma_1 = 1$ and $\sigma_2 = 0.01$ data are $\epsilon_2 = 0.732$ and $\epsilon_2 = 1.02$, respectively. The e-folding times when these values occurred were $N_1 - N_T = 4.7$ and $N_2 - N_T = 9.3$, respectively, and were separated by $N_1 - N_2 = \ln(\sigma_1/\sigma_2)$. Both values occurred during the constant $\epsilon_2$ phase, suggesting an “effective” $\epsilon_2$ value corresponding to a time shortly after the power spectrum peak scale crossed the coarse-graining scale. If the value of $\epsilon_2$ corresponding to the coarse-graining crossing time of the last mode integrated is used (solid lines), then the classical $\delta\mathcal{N}$ does
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The numerically reconstructed first-passage time distribution for the Gaussian bump potential (4.12). The semi-analytical noise model described in section 3.4 is used. Left: comparison of the 1D (4.6) and 2D (2.24) models. The value of $\sigma = k/(aH)$ defines when the noise is applied and the start/end field values. The solid black curve is Gaussian with variance given by eq. (2.14), i.e. it is the result of standard perturbation theory. Right: importance sampling using the 1D noise approach (4.7) to find the tail of the distribution for a range of bias amplitudes combined, for two $\sigma$ values shown, along with the classical $\delta N$ prediction assuming constant $\epsilon_2$ (2.20), colour coded to match. The solid lines used the $\epsilon_2$ value at coarse-graining crossing for the end of the integrated mode range, and the dashed lines are fitted values. The numerical error bars correspond to one standard deviation, as discussed in section 4.1. By coincidence, the $\sigma = 1$ solid curve and the $\sigma = 0.01$ dashed curve overlap.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Left: the same as the left-hand plot of figure 5 but $k_{max}$ has been shifted to the constant $\epsilon_2$ phase. Right: the same as the right-hand plot of figure 6 but for the mode range shown on the left panel.</figcaption></figure>

not fit the data. This is expected for the solid $\sigma = 0.01$ curve, as $\epsilon_2$ is not constant when the coarse-graining crossing of that scale occurred.

If however the range of modes integrated, depicted in figure 5, is changed such that the last mode integrated left the coarse-graining scale during the constant-$\epsilon_2$ phase, then the classical $\delta N$ prediction (2.20) provides a reasonable fit to the data. This is shown in figure 7 for $\sigma = 1$ data, with both classical $\delta N$ curves giving a reasonable fit. The $\sigma = 0.01$ data still requires a fitted $\epsilon_2$ but the classical $\delta N$ curve with $\epsilon_2$ value corresponding to when the last mode left the coarse-graining scale is a much closer fit than for figure 6.
We see an increased scatter in the importance sampling data in figure 6 compared to figures 4 and 7. This is because the correlation between the weights, $\ln(w)$, and the value of $\delta\mathcal{N}$ is weaker, increasing the error when the lognormal estimator is used. See ref. [43] for details.

# 5 Conclusions

In this paper we have, for the first time, shown how to numerically investigate the far tail of the first-passage-time probability distribution, $P(\delta\mathcal{N})$, using importance sampling in inflation models including a transient period of ultra-slow roll following a sudden transition. We investigated two specific inflaton potentials leading to strongly enhanced power spectra on small scales, $P_R \sim 5 \times 10^{-3}$, as typically studied in single-field models of primordial black hole production from inflation [28]. To do so we have extended the open-source PYFPT code to describe the stochastic dynamics in a 2D phase space, eq. (2.21), incorporating noise that departs from the standard Bunch-Davies profile due to non-trivial sub-Hubble dynamics during the transition. We have applied importance sampling techniques and a 1D noise model, see eq. (4.7). The noise was pre-computed on the classical background for efficiency. With this approach, we have been able to explore the probability distribution down to $P(\delta\mathcal{N}) \lesssim 10^{-10}$ in simulations taking less than 15 minutes to run on a laptop computer.

We have demonstrated that a 1D model for the noise was justified by the squeezing in phase space of the field perturbations on super-Hubble scales even after a sudden transition to ultra-slow roll, as was shown in figure 3. In our simulations we required that keeping only the growing mode of the homogeneous ($k=0$) field perturbation found at coarse graining, eq. (2.10), reproduced within 1% the linear power spectrum at the end of inflation obtained by evolving the full inhomogeneous field perturbation, eq. (2.8), see appendix A. This ensured the 1D noise model was consistent with linear perturbation theory for all the modes used in the simulations and the separate-universe approach was correctly applied [52].

In general it is non-trivial to identify the long-wavelength limit of the field perturbation from the numerical solution of the full mode equation at a finite wavelength at the coarse-graining time. One possibility is simply to wait long enough (long after Hubble exit) until the wavelength is large enough that the spatial gradients are negligible. However we showed that it was also possible to accurately identify the long-wavelength limit by matching the full numerical solution to a Bessel function ansatz soon after, or even at, Hubble exit. We were then able to faithfully reproduce the linear power spectrum even when coarse graining at the Hubble scale, $k=aH$. For a piece-wise linear potential, see section 4.2, our stochastic numerical results matched the classical Gaussian prediction far into the tail of the PDF (which we expect to be an excellent approximation in this linear case), showing the validity of the 1D noise approach.

In section 4.3 we then applied our PYFPT numerical code to a smooth inflaton potential with a small Gaussian bump, leading to a sudden transition to ultra-slow roll, on top of a plateau-type potential, compatible with CMB constraints on large-scale perturbations [78]. The probability distribution, $P(\delta\mathcal{N})$, is well-described by a Gaussian close to the peak of the distribution, but we found strong deviations from Gaussianity and an exponential tail for $P(\delta\mathcal{N}) \lesssim 10^{-5}$. A simple analytic function (2.20) gave an excellent fit to the distribution. This form for the PDF is consistent with the classical $\delta\mathcal{N}$ formalism [30, 33–35, 56] in
constant-roll inflation (where $\epsilon_2$ = constant), in agreement with previous studies [50, 51]. Since $\epsilon_2$ is not exactly constant in this model, there is some ambiguity in the precise value of $\epsilon_2$ used in eq. (2.20) introduced by the time variation, which is seen in figure 6. If however we stop the simulation when the variation of $\epsilon_2$ is very small, then the classical $\delta N$ provides an excellent fit, as seen in figure 7.

Let us also stress that in the Gaussian-bump model, the tail behaviour was found to depend quite significantly on the value of $\sigma$, i.e. on the ratio between the coarse-graining and the Hubble scale. In principle, $\sigma$ should be set such that the coarse-graining scale is the one at which backreaction effectively takes place in the local FLRW patches. There is no clear prescription for where that should be, and this is somehow inherent to the accuracy of the separate-universe picture itself, but values of order one for $\sigma$ seem reasonable (i.e. assuming that backreaction occurs at the Hubble scale). However, to correctly identify the long-wavelength limit and extract the asymptotic growing-mode may require us to wait long enough after Hubble exit such that the spatial gradients are negligible, and this is why values of $\sigma \ll 1$ are usually employed in stochastic inflation. This is indeed required in standard modelling of the noise, but in regimes where $\nu$ is approximately constant, the semi-analytical approach we have proposed was found to work accurately even at values of $\sigma = 1$. This makes it particularly well-suited to cases where a substantial dependence on $\sigma$ is found, as seems to be the case for non-Gaussian tails, where values of $\sigma \ll 1$ might underestimate backreaction effects.

The focus of this work was to propose a noise modelling and to design an importance-sampling approach that allows one to study post-transition phases during inflation, where slow roll is violated. This is why we focused on the statistics of the first-passage time arising from a fixed range of comoving scales, and computed the noise along a fixed background reference trajectory, to simplify the discussion and prevent more subtle stochastic effects to interfere in the analysis. This led us to situations where the stochastic-$\delta N$ results seem to be easily interpreted in terms of classical-$\delta N$ effects only. Although some works suggest that the heavy tails found in the stochastic-$\delta N$ formalism come from purely classical effects, with results from classical- and stochastic-$\delta N$ formalisms in excellent agreement [51], let us stress that this cannot be true in general. The example considered in section 4.2 is a good counter-example: since the classical dynamics gives rise to a linear relationship between $\phi$ and $N$ at late time, a Gaussian distribution is necessarily found. However, the full stochastic-$\delta N$ formalism always yields exponential tails [29, 37], due to each realisation integrating over a different range of modes (here we integrated over the same range), so the two approaches necessarily deviate. In addition, the curvature perturbation is never measured in terms of comoving scales, but rather in terms of physical scales. In the classical approach, going from comoving to physical scales is straightforward (although it can lead to non-trivial effects in more than one-dimensional phase spaces [80]) since there is a one-to-one correspondence between scales and field configurations when they exit the Hubble radius. In the stochastic picture however, this is no longer the case, convolutions across the whole inflating domain need to be performed [81–83]. This leads to drastic changes in observables involving large fluctuations that cannot be accounted for by the classical approach.
To include non-Gaussianities in the tail of $P(\delta\mathcal{N})$ due to stochastic effects beyond the classical $\delta\mathcal{N}$ approximation, one would need to go beyond pre-computing the noise using the classical background trajectory, i.e. one should include backreaction of the stochastic evolution on the noise itself. In principle this would involve solving the mode equation for quantum fluctuations at every time-step in every realisation which would be computationally prohibitive and lead to non-Markovian noise. However we have seen that in the decoupling limit (when $\epsilon_1 \to 0$) during ultra-slow roll, the noise at coarse graining could be described as a function of the local field value, using eq. (3.12) and eq. (3.14), independent of the particular classical phase-space trajectory, see figure 1. Hence the local field value at each step could be used to estimate the noise including the stochastic backreaction in each realisation. This has the advantage of not requiring the mode equations to be solved on each stochastic trajectory separately, greatly speeding up the numerics. It is still an open question if this approach is valid for a fully stochastic background, and we leave the investigation of this for future work. Additionally, how best to find the potentially $k$-dependent amplitude of the growing mode in this approach (the Bessel coefficient $B_k$ in eq. (3.12)) requires further consideration, which we also leave for future work.

Due to the sudden transition from slow roll to ultra-slow roll leading to a discontinuity in the behaviour of the homogeneous part of super-Hubble modes [52], we have only simulated the stochastic evolution after the transition. It remains an unsolved challenge to build a full stochastic model of super-Hubble modes through a sudden transition (i.e. in phases where $\nu$ varies quickly). Nonetheless in this work we have demonstrated the persistence of exponential tails in the probability distribution $P(\delta\mathcal{N})$ incorporating non-slow-roll dynamics and non-adiabatic perturbations immediately after a sudden transition, including the peak of the primordial power spectrum and more than 85% of the total power integrated over all observable scales.

Finally we note that while the diffusion-based bias model (4.10) used for importance sampling was able to probe the far tail of $P(\delta\mathcal{N})$ for the potentials investigated, it might not be sufficient in more general cases. This is due to the quality of the importance sampling data (the correlation between $\ln(w)$ and $\delta\mathcal{N}$, see ref. [43]) depending on the choice of bias and the model investigated. The problem of finding the optimal bias, which gives the highest quality data, has been solved in 1D [77] but is still an important open question in the 2D phase space.

**Data availability statement.** The data from which the results presented here are derived is freely available at https://github.com/Jacks0nJ/PyFPT, along with the code required to produce the data.

## Acknowledgments

The authors would like to thank Parth Bhargava, Laura Iacconi, Swagat Mishra, Shi Pi and Eemeli Tomberg for insightful discussions. This work was supported by the Science and Technology Facilities Council (grant numbers ST/T506345/1 and ST/W001225/1). For the purpose of open access, the authors have applied a Creative Commons Attribution (CC-BY) licence to any Author Accepted Manuscript version arising from this work.
<figure><img src="image_10.png" /><figcaption>**Figure 8.** The percentage difference in the linear power spectrum at the end of inflation, $P_R(k)$, obtained using different models for the noise to determine the amplitude of the homogeneous growing mode, $C_k$ in eq. (2.10), against the true value found by solving the full mode equation (2.4) for $\delta\phi_k$ in eq. (2.8) at the end of inflation. The solid curves correspond to the semi-analytical procedure, detailed in section 3.4, to identify the homogeneous part of the field fluctuations when $k = \sigma aH$, while the dashed curve corresponds to the simple analytical noise for the Bunch-Davies vacuum state for a massless field in de Sitter (3.3) leading to typical field fluctuations $H/2\pi$ at Hubble exit. The homogeneous growing mode is then calculated at the matching time using eq. (2.9). The vertical dashed black lines indicate the range of modes used in our stochastic simulations, corresponding to $k_{\text{min}}$ on the left and $k_{\text{max}}$ on the right. The horizontal dotted gray line is the 1% accuracy threshold. Left: the piece-wise linear potential (4.11). Right: the Gaussian-bump model (4.12).</figcaption></figure>

## A Noise model accuracy

In this appendix, the requirement that the homogeneous perturbations used to calculate the stochastic noise, $\delta\phi_{k_\sigma, h}$ and $\delta\pi_{k_\sigma, h}$ in eq. (2.22), give the correct linear power spectrum at the end of inflation using the separate-universe approach (3.1), is tested for the two potentials investigated. Here the homogeneous perturbations were found using the semi-analytical Bessel-function method of section 3.4, and keeping only the homogeneous growing mode, see eq. (3.12).

The results are shown for the piece-wise linear (4.11) and Gaussian bump (4.12) potentials in the left and right panels respectively of figure 8. For the piece-wise linear potential, the criterion given in eq. (3.1) is always met after the sudden transition. Therefore a range of modes to include the majority of the oscillations in the power spectrum post-transition, see figure 2, is chosen for the simulations. For the Gaussian-bump potential, the criterion is met for $0.01 \le \sigma \le 1$ for the range indicated by the vertical dashed lines. Note that using the standard $H/2\pi$ field fluctuations for a massless field in de Sitter to calculate the noise, i.e. using eq. (3.3), gives a $\gtrsim 10\%$ error close to the transition for the piece-wise linear potential and $\sim 100\%$ error throughout for the Gaussian bump potential, showing the importance of improving upon this simple analytical approximation.

## References

[1] A.A. Starobinsky, *A new type of isotropic cosmological models without singularity*, Phys. Lett. B **91** (1980) 99 [INSPIRE].
[2] K. Sato, First-order phase transition of a vacuum and the expansion of the Universe, Mon. Not. Roy. Astron. Soc. 195 (1981) 467 [INSPIRE].

[3] A.H. Guth, *The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems*, Phys. Rev. D 23 (1981) 347 [INSPIRE].

[4] A.D. Linde, *A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems*, Phys. Lett. B 108 (1982) 389 [INSPIRE].

[5] A. Albrecht and P.J. Steinhardt, *Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking*, Phys. Rev. Lett. 48 (1982) 1220 [INSPIRE].

[6] A.D. Linde, *Chaotic inflation*, Phys. Lett. B 129 (1983) 177 [INSPIRE].

[7] PLANCK collaboration, Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[8] M. Tristram et al., Improved limits on the tensor-to-scalar ratio using BICEP and Planck data, Phys. Rev. D 105 (2022) 083524 [arXiv:2112.07961] [INSPIRE].

[9] V.F. Mukhanov and G.V. Chibisov, *Quantum fluctuations and a nonsingular universe*, JETP Lett. 33 (1981) 532 [INSPIRE].

[10] V.F. Mukhanov and G.V. Chibisov, *Vacuum energy and large scale structure of the universe*, Sov. Phys. JETP 56 (1982) 258 [INSPIRE].

[11] A.A. Starobinsky, *Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations*, Phys. Lett. B 117 (1982) 175 [INSPIRE].

[12] A.H. Guth and S.Y. Pi, *Fluctuations in the New Inflationary Universe*, Phys. Rev. Lett. 49 (1982) 1110 [INSPIRE].

[13] S.W. Hawking, *The Development of Irregularities in a Single Bubble Inflationary Universe*, Phys. Lett. B 115 (1982) 295 [INSPIRE].

[14] J.M. Bardeen, P.J. Steinhardt and M.S. Turner, *Spontaneous creation of almost scale-free density perturbations in an inflationary universe*, Phys. Rev. D 28 (1983) 679 [INSPIRE].

[15] J. Chluba, J. Hamann and S.P. Patil, *Features and New Physical Scales in Primordial Observables: Theory and Observation*, Int. J. Mod. Phys. D 24 (2015) 1530023 [arXiv:1505.01834] [INSPIRE].

[16] Y.B. Zel'dovich and I.D. Novikov, *The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model*, Sov. Astron. 10 (1967) 602 [INSPIRE].

[17] S. Hawking, *Gravitationally collapsed objects of very low mass*, Mon. Not. Roy. Astron. Soc. 152 (1971) 75 [INSPIRE].

[18] B.J. Carr and S.W. Hawking, *Black holes in the early Universe*, Mon. Not. Roy. Astron. Soc. 168 (1974) 399 [INSPIRE].

[19] B. Carr, F. Kühnel and M. Sandstad, *Primordial Black Holes as Dark Matter*, Phys. Rev. D 94 (2016) 083504 [arXiv:1607.06077] [INSPIRE].

[20] B. Carr and F. Kühnel, *Primordial Black Holes as Dark Matter: Recent Developments*, Ann. Rev. Nucl. Part. Sci. 70 (2020) 355 [arXiv:2006.02838] [INSPIRE].

[21] A.M. Green and B.J. Kavanagh, *Primordial Black Holes as a dark matter candidate*, J. Phys. G 48 (2021) 043001 [arXiv:2007.10722] [INSPIRE].
[22] A.M. Green, *Primordial black holes as a dark matter candidate* - a brief overview, Nucl. Phys. B 1003 (2024) 116494 [arXiv:2402.15211] [INSPIRE].

[23] B. Carr, S. Clesse, J. García-Bellido and F. Kühnel, *Cosmic conundra explained by thermal history and primordial black holes*, Phys. Dark Univ. 31 (2021) 100755 [arXiv:1906.08217] [INSPIRE].

[24] K.N. Ananda, C. Clarkson and D. Wands, *The Cosmological gravitational wave background from primordial density perturbations*, Phys. Rev. D 75 (2007) 123518 [gr-qc/0612013] [INSPIRE].

[25] D. Baumann, P.J. Steinhardt, K. Takahashi and K. Ichiki, *Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations*, Phys. Rev. D 76 (2007) 084019 [hep-th/0703290] [INSPIRE].

[26] EPTA and INPTA collaborations, *The second data release from the European Pulsar Timing Array — IV. Implications for massive black holes, dark matter, and the early Universe*, Astron. Astrophys. 685 (2024) A94 [arXiv:2306.16227] [INSPIRE].

[27] NANOGRAV collaboration, *The NANOGlav 15 yr Data Set: Search for Signals from New Physics*, Astrophys. J. Lett. 951 (2023) L11 [Erratum ibid. 971 (2024) L27] [arXiv:2306.16219] [INSPIRE].

[28] LISA COSMOLOGY WORKING GROUP collaboration, *Primordial black holes and their gravitational-wave signatures*, Living Rev. Rel. 28 (2025) 1 [arXiv:2310.19857] [INSPIRE].

[29] C. Pattison, V. Vennin, H. Assadullahi and D. Wands, *Quantum diffusion during inflation and primordial black holes*, JCAP 10 (2017) 046 [arXiv:1707.00537] [INSPIRE].

[30] A.A. Starobinsky, *Stochastic De Sitter (inflationary) stage in the early universe*, Lect. Notes Phys. 246 (1986) 107 [INSPIRE].

[31] A.A. Starobinsky and J. Yokoyama, *Equilibrium state of a selfinteracting scalar field in the De Sitter background*, Phys. Rev. D 50 (1994) 6357 [astro-ph/9407016] [INSPIRE].

[32] A.A. Starobinsky, *Multicomponent de Sitter (Inflationary) Stages and the Generation of Perturbations*, JETP Lett. 42 (1985) 152 [INSPIRE].

[33] M. Sasaki and E.D. Stewart, *A General analytic formula for the spectral index of the density perturbations produced during inflation*, Prog. Theor. Phys. 95 (1996) 71 [astro-ph/9507001] [INSPIRE].

[34] M. Sasaki and T. Tanaka, *Superhorizon scale dynamics of multiscalar inflation*, Prog. Theor. Phys. 99 (1998) 763 [gr-qc/9801017] [INSPIRE].

[35] D.H. Lyth, K.A. Malik and M. Sasaki, *A General proof of the conservation of the curvature perturbation*, JCAP 05 (2005) 004 [astro-ph/0411220] [INSPIRE].

[36] V. Vennin and A.A. Starobinsky, *Correlation Functions in Stochastic Inflation*, Eur. Phys. J. C 75 (2015) 413 [arXiv:1506.04732] [INSPIRE].

[37] J.M. Ezquiaga, J. García-Bellido and V. Vennin, *The exponential tail of inflationary fluctuations: consequences for primordial black holes*, JCAP 03 (2020) 029 [arXiv:1912.05399] [INSPIRE].

[38] C. Animali and V. Vennin, *Primordial black holes from stochastic tunnelling*, JCAP 02 (2023) 043 [arXiv:2210.03812] [INSPIRE].

[39] V. Briaud and V. Vennin, *Uphill inflation*, JCAP 06 (2023) 029 [arXiv:2301.09336] [INSPIRE].

[40] D.G. Figueroa, S. Raatikainen, S. Rasanen and E. Tomberg, *Non-Gaussian Tail of the Curvature Perturbation in Stochastic Ultraslow-Roll Inflation: Implications for Primordial Black Hole Production*, Phys. Rev. Lett. 127 (2021) 101302 [arXiv:2012.06551] [INSPIRE].
[41] D.G. Figueroa, S. Raatikainen, S. Rasanen and E. Tomberg, *Implications of stochastic effects for primordial black hole production in ultra-slow-roll inflation*, JCAP 05 (2022) 027 [arXiv:2111.07437] [INSPIRE].

[42] O. Mazonka, C. Jarzynski and J. Blocki, *Computing probabilities of very rare events for Langevin processes: A New method based on importance sampling*, Nucl. Phys. A 641 (1998) 335 [nucl-th/9809075] [INSPIRE].

[43] J.H.P. Jackson et al., *Numerical simulations of stochastic inflation using importance sampling*, JCAP 10 (2022) 067 [arXiv:2206.11234] [INSPIRE].

[44] M. Biagetti, G. Francolini, A. Kehagias and A. Riotto, *Primordial Black Holes from Inflation and Quantum Diffusion*, JCAP 07 (2018) 032 [arXiv:1804.07124] [INSPIRE].

[45] H. Firouzjahi, A. Nassiri-Rad and M. Noorbala, *Stochastic nonattractor inflation*, Phys. Rev. D 102 (2020) 123504 [arXiv:2009.04680] [INSPIRE].

[46] Y.-F. Cai et al., *Highly non-Gaussian tails and primordial black holes from single-field inflation*, JCAP 12 (2022) 034 [arXiv:2207.11910] [INSPIRE].

[47] S. Pi and M. Sasaki, *Logarithmic Duality of the Curvature Perturbation*, Phys. Rev. Lett. 131 (2023) 011002 [arXiv:2211.13932] [INSPIRE].

[48] S. Hooshangi, M.H. Namjoo and M. Noorbala, *Tail diversity from inflation*, JCAP 09 (2023) 023 [arXiv:2305.19257] [INSPIRE].

[49] X. Wang, Y.-L. Zhang and M. Sasaki, *Enhanced curvature perturbation and primordial black hole formation in two-stage inflation with a break*, JCAP 07 (2024) 076 [arXiv:2404.02492] [INSPIRE].

[50] E. Tomberg, *Stochastic constant-roll inflation and primordial black holes*, Phys. Rev. D 108 (2023) 043502 [arXiv:2304.10903] [INSPIRE].

[51] G. Ballesteros et al., *Non-Gaussian tails without stochastic inflation*, JCAP 11 (2024) 013 [arXiv:2406.02417] [INSPIRE].

[52] J.H.P. Jackson et al., *The separate-universe approach and sudden transitions during inflation*, JCAP 05 (2024) 053 [arXiv:2311.03281] [INSPIRE].

[53] D. Artigas, S. Pi and T. Tanaka, *Extended $\delta N$ formalism*, arXiv:2408.09964 [INSPIRE].

[54] K. Dimopoulos, *Ultra slow-roll inflation demystified*, Phys. Lett. B 775 (2017) 262 [arXiv:1707.05644] [INSPIRE].

[55] B.A. Bassett, S. Tsujikawa and D. Wands, *Inflation dynamics and reheating*, Rev. Mod. Phys. 78 (2006) 537 [astro-ph/0507632] [INSPIRE].

[56] D. Wands, K.A. Malik, D.H. Lyth and A.R. Liddle, *A New approach to the evolution of cosmological perturbations on large scales*, Phys. Rev. D 62 (2000) 043527 [astro-ph/0003278] [INSPIRE].

[57] C. Pattison, V. Vennin, H. Assadullahi and D. Wands, *Stochastic inflation beyond slow roll*, JCAP 07 (2019) 031 [arXiv:1905.06300] [INSPIRE].

[58] H. Kodama and T. Hamazaki, *Evolution of cosmological perturbations in the long wavelength limit*, Phys. Rev. D 57 (1998) 7177 [gr-qc/9712045] [INSPIRE].

[59] C. Gordon, D. Wands, B.A. Bassett and R. Maartens, *Adiabatic and entropy perturbations from inflation*, Phys. Rev. D 63 (2000) 023506 [astro-ph/0009131] [INSPIRE].
[60] A.E. Romano, S. Mooij and M. Sasaki, *Adiabaticity and gravity theory independent conservation laws for cosmological perturbations*, Phys. Lett. B 755 (2016) 464 [arXiv:1512.05757] [INSPIRE].

[61] K.A. Malik and D. Wands, *Cosmological perturbations*, Phys. Rept. 475 (2009) 1 [arXiv:0809.4944] [INSPIRE].

[62] R. Inui et al., *Constant roll and non-Gaussian tail in light of logarithmic duality*, JCAP 02 (2025) 042 [arXiv:2409.13500] [INSPIRE].

[63] D. Polarski and A.A. Starobinsky, *Semiclassicality and decoherence of cosmological perturbations*, Class. Quant. Grav. 13 (1996) 377 [gr-qc/9504030] [INSPIRE].

[64] J. Grain and V. Vennin, *Stochastic inflation in phase space: Is slow roll a stochastic attractor?*, JCAP 05 (2017) 045 [arXiv:1703.00447] [INSPIRE].

[65] C. Pattison, V. Vennin, D. Wands and H. Assadullahi, *Ultra-slow-roll inflation with quantum diffusion*, JCAP 04 (2021) 080 [arXiv:2101.05741] [INSPIRE].

[66] H. Firouzjahi, A. Nassiri-Rad and M. Noorbala, *Stochastic Ultra Slow Roll Inflation*, JCAP 01 (2019) 040 [arXiv:1811.02175] [INSPIRE].

[67] E.D. Stewart and D.H. Lyth, *A More accurate analytic calculation of the spectrum of cosmological perturbations produced during inflation*, Phys. Lett. B 302 (1993) 171 [gr-qc/9302019] [INSPIRE].

[68] A. De and R. Mahbub, *Numerically modeling stochastic inflation in slow-roll and beyond*, Phys. Rev. D 102 (2020) 123509 [arXiv:2010.12685] [INSPIRE].

[69] N. Ahmadi et al., *Quantum diffusion in sharp transition to non-slow-roll phase*, JCAP 08 (2022) 078 [arXiv:2207.10578] [INSPIRE].

[70] S.S. Mishra, E.J. Copeland and A.M. Green, *Primordial black holes and stochastic inflation beyond slow roll*. Part I. Noise matrix elements, JCAP 09 (2023) 005 [arXiv:2303.17375] [INSPIRE].

[71] P.E. Kloeden and E. Platen, *Numerical Solution of Stochastic Differential Equations*, Springer Berlin Heidelberg (1992) [DOI:10.1007/978-3-662-12616-5].

[72] A.A. Starobinsky, *Spectrum of adiabatic perturbations in the universe when there are singularities in the inflation potential*, JETP Lett. 55 (1992) 489 [INSPIRE].

[73] S.M. Leach, M. Sasaki, D. Wands and A.R. Liddle, *Enhancement of superhorizon scale inflationary curvature perturbations*, Phys. Rev. D 64 (2001) 023512 [astro-ph/0101406] [INSPIRE].

[74] J. Martin and L. Sriramkumar, *The scalar bi-spectrum in the Starobinsky model: The equilateral case*, JCAP 01 (2012) 008 [arXiv:1109.5838] [INSPIRE].

[75] J. Martin, L. Sriramkumar and D.K. Hazra, *Sharp inflaton potentials and bi-spectra: Effects of smoothening the discontinuity*, JCAP 09 (2014) 039 [arXiv:1404.6093] [INSPIRE].

[76] S. Pi and J. Wang, *Primordial black hole formation in Starobinsky's linear potential model*, JCAP 06 (2023) 018 [arXiv:2209.14183] [INSPIRE].

[77] E. Tomberg, *Numerical stochastic inflation constrained by frozen noise*, JCAP 04 (2023) 042 [arXiv:2210.17441] [INSPIRE].

[78] S.S. Mishra and V. Sahni, *Primordial Black Holes from a tiny bump/dip in the Inflaton potential*, JCAP 04 (2020) 007 [arXiv:1911.00057] [INSPIRE].
[79] P.S. Cole, A.D. Gow, C.T. Byrnes and S.P. Patil, *Primordial black holes from single-field inflation: a fine-tuning audit*, JCAP 08 (2023) 031 [arXiv:2304.01997] [INSPIRE].

[80] Y. Tada and V. Vennin, *Squeezed bispectrum in the $\delta N$ formalism: local observer effect in field space*, JCAP 02 (2017) 021 [arXiv:1609.08876] [INSPIRE].

[81] Y. Tada and V. Vennin, *Statistics of coarse-grained cosmological fields in stochastic inflation*, JCAP 02 (2022) 021 [arXiv:2111.15280] [INSPIRE].

[82] K. Ando and V. Vennin, *Power spectrum in stochastic inflation*, JCAP 04 (2021) 057 [arXiv:2012.02031] [INSPIRE].

[83] C. Animali and V. Vennin, *Clustering of primordial black holes from quantum diffusion during inflation*, JCAP 08 (2024) 026 [arXiv:2402.08642] [INSPIRE].