# Numerical simulations of stochastic inflation using importance sampling

Joseph H.P. Jackson, Hooshyar Assadullahi, Kazuya Koyama, Vincent

Vennin, David Wands

► To cite this version:

Joseph H.P. Jackson, Hooshyar Assadullahi, Kazuya Koyama, Vincent Vennin, David Wands. Numerical simulations of stochastic inflation using importance sampling. Journal of Cosmology and Astroparticle Physics, 2022, 10, pp.067. 10.1088/1475-7516/2022/10/067. hal-03712829

HAL Id: hal-03712829

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

Submitted on 9 Oct 2024

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>

Distributed under a Creative Commons Attribution 4.0 International License
# Numerical simulations of stochastic inflation using importance sampling

Joseph H. P. Jackson,ª Hooshyar Assadullahi,ª,b Kazuya Koyama,ª Vincent Vennin, c,d,a David Wandsª

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

**Abstract.** We show how importance sampling can be used to reconstruct the statistics of rare cosmological fluctuations in stochastic inflation. We have developed a publicly available package, PYFPT,<fnref n="1" /> that solves the first-passage time problem of generic one-dimensional Langevin processes. In the stochastic-δN formalism, these are related to the curvature perturbation at the end of inflation. We apply this method to quadratic inflation, where the existence of semi-analytical results allows us to benchmark our approach. We find excellent agreement within the estimated statistical error, both in the drift- and diffusion-dominated regimes. The computation takes at most a few hours on a single CPU, and can reach probability values corresponding to less than one Hubble patch per observable universe at the end of inflation. With direct sampling, this would take more than the age of the universe to simulate even with the best current supercomputers. As an application, we study how the presence of large-field boundaries might affect the tail of the probability distribution. We also find that non-perturbative deviations from Gaussianity are not always of the simple exponential type.

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

<footnote>ªInstitute of Cosmology & Gravitation, University of Portsmouth, Dennis Sciama Building, Burnaby Road, Portsmouth, PO1 3FX, United Kingdom</footnote>

<footnote>ᵇSchool of Mathematics and Physics, University of Portsmouth, Lion Gate Building, Lion Terrace, Portsmouth, PO1 3HF, United Kingdom</footnote>

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

<footnote>ᵈLaboratoire Astroparticule et Cosmologie, CNRS, Université Paris Cité, 75013 Paris, France</footnote>

<footnote n="1">¹https://github.com/JacksOnJ/PyFPT</footnote>
# Contents

<table><tr><td>1 Introduction</td><td>1</td></tr><tr><td>2 Importance sampling</td><td>4</td></tr><tr><td>2.1 Direct sampling of a Langevin equation</td><td>4</td></tr><tr><td>2.2 Importance sampling of a Langevin equation</td><td>5</td></tr><tr><td>2.3 Stochastic inflation</td><td>6</td></tr><tr><td>2.4 Lognormal estimator</td><td>7</td></tr><tr><td>3 Quadratic inflation</td><td>9</td></tr><tr><td>3.1 Diffusion domination: m ≥ M<sub>P1</sub></td><td>11</td></tr><tr><td>3.2 Interim case: m ≤ M<sub>P1</sub></td><td>13</td></tr><tr><td>3.3 Drift domination: m ≪ M<sub>P1</sub></td><td>16</td></tr><tr><td>3.4 Bias optimization</td><td>17</td></tr><tr><td>4 Conclusions</td><td>19</td></tr><tr><td>A Semi-analytical results</td><td>22</td></tr><tr><td>B Edgeworth expansion</td><td>26</td></tr></table>

## 1 Introduction

Inflation is a period of accelerated cosmic expansion in the very early universe [1–6], proposed to explain its observed homogeneity, isotropy and flatness [7–13]. The process of accelerated expansion leads to microscopic quantum fluctuations in light scalar fields growing to macroscopic scales [14–19]. These result in primordial curvature perturbations, $\zeta$, whose imprint can be observed in the Cosmic Microwave Background (CMB) [11, 13, 20] and which act as the seeds of cosmic structure [21]. Measurements of the CMB provide tight constraints on the dynamics of inflation on scales exiting the horizon 50-60 e-folds before the end of inflation [22]. However, the properties of inflation on scales smaller than those observed in the CMB are not as strongly constrained.

While future experiments [23] will be able to constrain inflation by observations of the stochastic gravitational wave background [24], a complementary probe is the possible detection of black holes formed from very large perturbations produced by inflation, known as Primordial Black Holes (PBHs) [25–27]. Not only would the detection of PBHs (or lack thereof) constrain inflation on smaller scales, but PBHs may also explain the origin of supermassive black holes [28], the LIGO–Virgo–KAGRA gravitational wave detections [29–33] and some (if not all) of the dark matter, see Refs [34–37] for recent reviews.
In the standard approach to estimate PBH abundances, perturbation theory is used for canonical single-field inflation where the background evolution of the inflaton field is given by the Klein-Gordon equation

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

In this expression, an over-dot denotes a time derivative, $V(\phi)$ is the field's potential and $H = \dot{a}/a$ is the Hubble rate, with $a$ the scale factor of the universe. It is related to the inflaton field and its velocity by Friedmann's equation $H^2 = [V(\phi) + \dot{\phi}^2/2]/(3M_{Pl}^2)$, with $M_{Pl}$ being the reduced Planck mass.

In linear perturbation theory, the probability density function (PDF) of the primordial curvature perturbation is given by a Gaussian

$$ P(\zeta) \propto \exp\left(-\frac{\zeta^2}{2\sigma_{\zeta}^2}\right), \qquad (1.2) $$

where $\sigma_{\zeta}^2$ is the variance of $\zeta$ on the scale considered. Large density perturbations, resulting from large values of $\zeta$ in the tail of the PDF, collapse to form PBHs. A criterion to form a black hole can be given in terms of the compaction function [38, 39] (or its smoothed version [40]), which is non-linearly related to $\zeta$ [41, 42]. A PBH forms when the compaction exceeds a threshold value, with a mass given by critical scaling [43]. The production of PBHs based on the PDF given in Eq. (1.2) is therefore expected to be suppressed by a Gaussian factor, and given that $\sigma_{\zeta}^2$ is observed to be small on CMB scales [13], a large enhancement in $\sigma_{\zeta}^2$ on small scales is required for a significant abundance of PBHs [35].

One mechanism for producing a peak in $\sigma_{\zeta}^2$, is a period of inflation about an inflection point in the scalar field potential driving inflation, where $dV/d\phi \to 0$, resulting in the power spectrum growing rapidly [44–50]. However, this can also lead to quantum diffusion effects becoming non-negligible [51–56] (see Refs. [57, 58] for a different viewpoint).

The stochastic approach to inflation enables us to study non-perturbative effects where the quantum diffusion can be large. This formalism introduces a coarse-graining scale which separates short and long wavelength modes of the scalar field. Quantum field fluctuations on short wavelengths are swept up into the long wavelength regime, where they are incorporated in the coarse-grained field, $\bar{\phi}$, resulting in a stochastic noise term, $\xi$, in the dynamical equations [16, 59–67]. We thus model inflation as a non-perturbative stochastic process, with the results of linear perturbation theory recovered in the low-diffusion limit [68]. In the slow-roll approximation [neglecting the field acceleration in Eq. (1.1)], this process is described by a first-order Langevin equation [59]

$$ \frac{\partial \bar{\phi}}{\partial N} = -\frac{1}{3H^2(\bar{\phi})} \frac{dV(\bar{\phi})}{d\bar{\phi}} + \frac{H(\bar{\phi})}{2\pi} \xi, \qquad (1.3) $$

where the local Hubble rate $H$ is now a function of the local coarse-grained field,

$$ H^2(\bar{\phi}) = \frac{V(\bar{\phi})}{3M_{Pl}^2}, \qquad (1.4) $$
with $\xi$ being a white Gaussian noise. As the number of e-folds elapsed during inflation, $N = \int H dt$, has been used as the time variable, we are implicitly working in the uniform-$N$ gauge [69] and this choice allows $\zeta$ to be found by use of the $\delta N$ formalism [16, 70–73]. In this approach, the integrated local expansion $\mathcal{N}$ of a homogeneous patch, treated as a separate universe [71, 74–78], is measured from an initially flat hypersurface to a hypersurface of uniform energy density. The curvature perturbation $\zeta$ is then given by the difference between this local expansion and its mean value,

$$ \zeta = \mathcal{N} - \langle \mathcal{N} \rangle \equiv \delta \mathcal{N}, \qquad (1.5) $$

where angle brackets $\langle \cdot \rangle$ denote the ensemble average. Finding the PDF of $\mathcal{N}$ thus corresponds to solving a first-passage time (FPT) problem [68]. This is the stochastic-$\delta N$ formalism [68, 79, 80], which allows $\zeta$ to be calculated beyond perturbation theory. The statistics of the curvature perturbation (and other quantities of cosmological interest such as the density contrast or the compaction function) when coarse-grained at a fixed physical scale can then be reconstructed using backward probabilities [81, 82].

If the inflating domain is bounded in field space, it has been shown that the PDF of $\mathcal{N}$ (and hence of $\zeta$) is given by a sum of decaying exponentials [51]

$$ P_{\phi}(\mathcal{N}) = \sum_{n} a_{n}(\phi) e^{-\Lambda_{n} \mathcal{N}}. \qquad (1.6) $$

Here $\phi$ denotes the initial field configuration (from now on we drop the explicit over-bar notation denoting the coarse-grained field for convenience), on which the decay rates $\Lambda_n$ do not depend. Technically, these decays rates appear as poles of the characteristic function of the PDF [83], hence they will be refereed to as “poles” in what follows. Therefore, while the peak of the PDF may be well approximated by a Gaussian, the far tail is rather exponential. Beyond a very few test cases [51, 83] for which this exponential tail can be calculated analytically, numerical simulations are in general required to reconstruct the PDF of $\mathcal{N}$ [84–89]. However, direct simulations primarily sample the peak of the distribution, corresponding to the most likely realisations, while PBHs form from those rare fluctuations living in the far tail. This implies that billions of simulations need to be run on supercomputers, from which only a tiny fraction is kept to reconstruct the tail. The poor efficiency of direct sampling therefore calls for new approaches.

In this work, we explain how the method of importance sampling [90, 91] can address this issue, and we apply it to stochastic inflation for the first time. Importance sampling deliberately over-samples the rare, large $\zeta$ events, which are then re-weighted to recover the true probability distribution far into the tail. We have developed the publicly available PYFPT package,<fnref n="1" /> which is general and applicable to any one-dimensional Langevin equation. Here we apply it to stochastic inflation.

This paper is organised as follows. In Sec. 2 we introduce the importance sampling method and the associated data-analysis techniques. In Sec. 3 we use both analytical and semi-analytical test cases to illustrate the accuracy of the importance sampling method

<footnote n="1">¹https://github.com/JacksOnJ/PyFPT</footnote>
and present results for the case of slow-roll inflation driven by a scalar field with a quadratic potential. Supporting calculations for this test case are given in Appendices A and B. We draw our conclusions in Sec. 4.

# 2 Importance sampling

## 2.1 Direct sampling of a Langevin equation

Consider a general one-dimensional Langevin equation

$$ \frac{dx}{dt} = D(t, x) + S(t, x)\xi, \qquad (2.1) $$

where $x$ is the stochastic variable to be propagated, $t$ is the time variable, $D(t, x)$ is the deterministic drift and $S(t, x)$ is the amplitude of the stochastic diffusion. $\xi$ is a random white Gaussian noise, normalised such that $\langle \xi(t)\xi(t') \rangle = \delta(t-t')$. Starting from an initial state $x = x_m$ at $t = t_m$, an estimate for the stochastic variable $x = x_{m+1}$ at time $t_{m+1} = t_m + \Delta t_m$ is given, for a sufficiently small interval $\Delta t_m$, by [92]

$$ x_{m+1} = x_m + D(t_m, x_m)\Delta t_m + S(t_m, x_m)\xi_m \sqrt{\Delta t_m}. \qquad (2.2) $$

Here, $\xi_m$ is a random number drawn from a normal distribution with unit variance, meaning that each step is given by a Gaussian with a mean $D(t_m, x_m)\Delta t_m$ and a standard deviation $S(t_m, x_m)\sqrt{\Delta t_m}$. This is the Euler-Maruyama method [92] of solving a stochastic differential equation in the Itô prescription, which we adopt for numerical simplicity. In single-field slow-roll inflation it produces similar results to the Stratonovitch prescription, which has been shown to preserve field-space covariance in the multiple field setup [93]. In a first-passage time problem (FPT), the procedure (2.2) is repeated until a given final condition $x_{\text{end}}$ is reached, and the corresponding elapsed time $t_{\text{FPT}}$ is recorded.

To estimate the PDF for the first-passage times, $P(t_{\text{FPT}})$, multiple completed simulations, known as runs, are required. These runs can then be binned, with our estimate of the PDF of $j$th bin with $t_{\text{FPT}}^{(j)} < t_{\text{FPT}} \le t_{\text{FPT}}^{(j+1)}$ given by

$$ \hat{P}(t_{\text{FPT}}^{(j)}) = \frac{n_j}{(t_{\text{FPT}}^{(j+1)} - t_{\text{FPT}}^{(j)}) n_{\text{total}}}, \qquad (2.3) $$

where $n_j$ is the number of runs in the $j$th bin and $n_{\text{total}} = \sum_j n_j$ is the total number of simulations.

A ‘hat’ is used to indicate that Eq. (2.3) is only a numerical estimate of the true distribution $P(t_{\text{FPT}})$. There are indeed two sources of errors. First, since finite time steps $\Delta t_m$ are used in Eq. (2.2), each run comes with numerical error. Second, since only a finite number of runs is available, Eq. (2.3) is subject to statistical error, which scales as $1/\sqrt{n_{\text{total}}}$ according to the central limit theorem [94], given that the runs are independent. By dividing the full set of runs into $n_{\text{sub}}$ subsets, the estimate (2.3) can be
computed in each subset, and the variance of the results provides an estimate of the error when the sample size is $n_{total}/n_{sub}$. This can then be extrapolated to the full sample size, $n_{total}$, using the scaling mentioned above. This is called jackknife resampling, which in this work we find to be reliable for $n_{sub} > 10$ and in bins where $n_j > 400$.

As mentioned above, this direct-sampling method is efficient at reconstructing the peak of the PDF, and thus provides reliable estimates of e.g. its lowest moments. However, rare events lying in the tail suffer from large statistical error.

## 2.2 Importance sampling of a Langevin equation

To efficiently investigate the tail of the probability distribution, a bias can be introduced into the Langevin equation (2.1), to increase the occurrence of rare realisations of the stochastic process. This is done by modifying the drift term [91]

$$ \frac{dx}{dt} = [D(t, x) + \mathcal{B}(t, x)] + S(t, x)\xi, \qquad (2.4) $$

with the associated Euler-Maruyama step (2.2) becoming

$$ x_{m+1} - x_m = [D(t_m, x_m) + \mathcal{B}(t_m, x_m)]\Delta t_m + S(t_m, x_m)\xi_m \sqrt{\Delta t_m}. \qquad (2.5) $$

$\mathcal{B}$ is known as the bias, as it modifies the mean of the Gaussian distribution for each step. By choosing an appropriate bias, the resulting stochastic process can primarily sample the area of interest of the probability distribution function. Importance sampling is then achieved by calculating the probability associated with the unmodified Langevin equation (2.1), known as the target distribution (T), relative to the distribution given by the modified equation (2.4), known as the sample distribution (S).

To calculate this relative probability, known as the weight $w$, consider a run made of $M$ numerical steps $\mathbf{X} = (x_1, x_2, \dots, x_M)$ starting from some initial value $x_0$ and using the modified Langevin equation (2.4). The statistical weight of this run is defined as

$$ w(\mathbf{X}) = \frac{p_T(\mathbf{X}|x_0)}{p_S(\mathbf{X}|x_0)}, \qquad (2.6) $$

where $p_T$ denotes the probability that the run $\mathbf{X}$ is generated by the target stochastic process (2.1), and $p_S$ is the same probability in the sample process (2.4). Since these processes are Markovian, the steps $x_{m+1} - x_m$ are independent random variables, with Gaussian distributions since they are linearly related to $\xi_m$. The probability associated to the run $\mathbf{X}$ can thus be written as the product of the probabilities for each step,

$$ p_T(\mathbf{X}|x_0) = \prod_{m=0}^{M-1} \frac{1}{S(t_m, x_m)\sqrt{2\pi\Delta t_m}} e^{-\frac{[x_{m+1} - x_m - D(t_m, x_m)\Delta t_m]^2}{2S^2(t_m, x_m)\Delta t_m}}, \qquad (2.7) $$

$$ p_S(\mathbf{X}|x_0) = \prod_{m=0}^{M-1} \frac{1}{S(t_m, x_m)\sqrt{2\pi\Delta t_m}} e^{-\frac{\{x_{m+1} - x_m - [D(t_m, x_m) + \mathcal{B}(t_m, x_m)]\Delta t_m\}^2}{2S^2(t_m, x_m)\Delta t_m}}. \qquad (2.8) $$
By substituting these expressions into Eq. (2.6), one obtains for the weight

$$w(\mathbf{X}) = \exp \left\{ \sum_{m=0}^{M-1} \frac{\mathcal{B}(t_m, x_m)}{S^2(t_m, x_m)} \left[ \frac{\mathcal{B}(t_m, x_m)}{2} \Delta t_m + S(t_m, x_m) \xi_m \sqrt{\Delta t_m} \right] \right\}. \quad (2.9)$$

Note that $\xi_m$ is the noise of the sample process.

When solving Eq. (2.4), one can then update the sum appearing in Eq. (2.9) at each step. When a run $\mathbf{X}$ is complete, the weight is recorded alongside the first-passage time, and the target PDF can be estimated as

$$\hat{P}(t_{\text{FPT}}^{(j)}) = \frac{\sum_{i=0}^{n_j-1} w_j^{(i)}}{(t_{\text{FPT}}^{(j+1)} - t_{\text{FPT}}^{(j)}) n_{\text{total}}}, \quad (2.10)$$

where $w_j^{(i)}$ is weight of the $i$th run belonging to the $j$th bin. Hereafter this will be refereed to as the “naïve” estimate of the PDF. We recover the simple estimate (2.3) of the PDF for a direct simulation, where $\mathcal{B} = 0$, since each run has $w_j^{(i)} = 1$ that case. This method is straightforward to implement and jackknife resampling can still be used to estimate the uncertainty. However we shall see in the following that if there is a large dispersion in the weights associated with different runs within a given bin, then this naïve method may underestimate the statistical error, which will lead us to consider alternative estimates in such cases.

## 2.3 Stochastic inflation

Let us now apply importance sampling to stochastic inflation. Upon comparing Eqs. (1.3) and (2.1), one can see that the drift and diffusion terms are given by

$$D = -\frac{1}{3H^2(\phi)} \frac{dV(\phi)}{d\phi}, \qquad S = \frac{H(\phi)}{2\pi}, \qquad (2.11)$$

if time is labeled by the number of e-folds $N$ and $H^2 = V/(3M_{\text{Pl}}^2)$. To simulate a particular model of slow-roll inflation, the PYFPT package requires the potential function, $V(\phi)$, its derivative, $dV/d\phi$, and the start and end values of the field, $\phi_{\text{in}}$ and $\phi_{\text{end}}$, as inputs.

In practice, we set the bias to only be a function of the field, which we parameterise as

$$\mathcal{B}(\phi) = \mathcal{A}f(\phi), \qquad (2.12)$$

where $\mathcal{A}$ is the bias amplitude and $f(\phi)$ describes its field dependence. The function $f(\phi)$ can be set to optimise the convergence of the importance-sampling procedure, as we will further discuss in Sec. 3.4. Once $f$ is fixed, varying $\mathcal{A}$ leads to sampling different regions of the distribution. Indeed, the mean number of e-folds in the sample distribution reads [68, 95]

$$\langle \mathcal{N} \rangle_{\text{S}}(\phi) = 2 \int_{\phi_{\text{end}}}^{\phi} d\phi_1 \int_{\phi_1}^{\phi_{\text{UV}}} \frac{d\phi_2}{S^2(\phi_2)} \exp \left[ 2 \int_{\phi_1}^{\phi_2} d\phi_3 \frac{D(\phi_3) + \mathcal{B}(\phi_3)}{S^2(\phi_3)} \right]. \quad (2.13)$$
Here, inflation starts from $\phi$ and is allowed to take place between $\phi_{end}$ (where it ends) and $\phi_{UV}$ (which may be infinite and where a reflective boundary is placed [96, 97]). By varying $\mathcal{A}$, one can tune $\langle N \rangle_S$, hence one can choose the typical values of $\mathcal{N}$ that are best sampled. In practice, until we further discuss how the $f$ function can be optimised, we take it to match the noise amplitude, i.e. $f(\phi) = S(\phi)$.

## 2.4 Lognormal estimator

The value of the weights within the $j$th bin, $w_j^{(i)}$, can vary by several orders of magnitude when sampling the far tail of $P(\mathcal{N})$, as illustrated in the left-hand panel of Fig. 1, for the model that will be considered in Sec. 3. In this case, a very large number of simulation runs are needed to prevent the naïve estimator (2.10) from being dominated by only a small number of runs with the highest weights, leading to a systematic error in $\hat{P}(\mathcal{N})$. This greatly reduces the numerical efficiency.

To illustrate this problem, consider the distribution of weights for the bin centred at $\mathcal{N} = 10.02$ in Fig. 1 for quadratic inflation (3.1) with mass $m = 0.001 M_{Pl}$, initial field value $\phi_{in} = \sqrt{42} M_{Pl}$ and bias amplitude $\mathcal{A} = 3$. As the value of $w$ varies by 8 orders of magnitude in this bin, only the runs with the largest, but rare, $w$ values dominate the naïve estimate of the PDF in Eq. (2.10). The majority of the runs contribute very little as they have much smaller weights. If by chance a few more, very large $w$ values are sampled, then an over estimation occurs. Equally, if very few of the rare, but large $w$ values are randomly sampled, then an underestimation occurs. The jackknife resampling used to estimate the uncertainty in the estimate, $\hat{P}(\mathcal{N})$, also suffers a systematic error for similar reasons. This effect can be seen in the right-hand panel of Fig. 1, where most points obtained using the naïve estimator lie below the analytical prediction, and just two points lie above the line.

In practice, one requires at least $n > \langle w \rangle^2 \exp(\sigma_{\ln w}^2)$ runs to overcome this effect [98], where $\sigma_{\ln w}$ is the standard deviation of $\ln w$. Therefore, for the simulation shown in Fig. 1, which already contains $2 \times 10^5$ runs, an order of magnitude more simulations are required for the naïve estimator $\hat{P}(\mathcal{N})$ to converge on the true $P(\mathcal{N})$.

The situation can be improved provided the distribution function of the weights is known, or at least can be approximated. In the present case, with the choice $f(\phi) = S(\phi)$, Eq. (2.9) leads to

$$w(\mathbf{X}) = \exp\left(\frac{\mathcal{A}^2}{2} \sum_{m=0}^{M-1} \Delta t_m + \mathcal{A} \sum_{m=0}^{M-1} \sqrt{\Delta t_m \xi_m}\right) = \exp\left(\frac{\mathcal{A}^2}{2} \mathcal{N}_S + \mathcal{A} \sqrt{\mathcal{N}_S \xi}\right), \quad (2.14)$$

where $\xi$ is a random normal variable with vanishing mean and unit variance, independent of $\mathcal{N}_S$, namely the duration of inflation in the sample process. This leads to the following distribution for the weight

$$P(w) = \frac{1}{w\sqrt{2\pi}} \int_{0}^{\infty} d\mathcal{N}_S \frac{P(\mathcal{N}_S)}{\mathcal{A}\sqrt{\mathcal{N}_S}} \exp\left[-\frac{\left(\ln w - \frac{\mathcal{A}^2\mathcal{N}_S}{2}\right)^2}{2\mathcal{A}^2\mathcal{N}_S}\right]. \quad (2.15)$$
<figure><img src="image_2.png" /><figcaption>Figure 1. Left panel: a 2D histogram of the logarithm of the weights, $log_{10}(w)$, with first-passage times, $N$, for an importance sampling simulation of $2 \times 10^5$ runs for quadratic inflation (3.1) with mass $m = 0.001 M_{\text{Pl}}$, initial field value $\phi_{\text{in}} = \sqrt{42} M_{\text{Pl}}$ and bias amplitude $A = 3$. Right panel: estimations of the PDF, $\hat{P}(N)$, using this data, with the naïve method (2.10), shown in blue and the lognormal method (2.17), shown in orange. The analytical prediction for $P(N)$ in the low-diffusion limit, corresponding to a Gaussian distribution (see Appendix B), is shown as a solid black curve.</figcaption></figure>

This relates the PDF of $w$ to the PDF of $N$, which is precisely the object we are trying to reconstruct. However, it can serve as the basis of an iterative procedure. Indeed, in the limit where diffusion is sub-dominant, $P(N_S)$ is a Dirac distribution centered on its classical value $N_{S,cl} = \int_{\phi}^{\phi_{\text{end}}} d\phi'/[D(\phi') + B(\phi')]$, which can be obtained by either setting the noise to zero in Eq. (2.4) or by performing a saddle-point expansion of Eq. (2.12) in the same limit. In that case, $P(w)$ is nothing but a lognormal distribution, i.e. $\ln w$ is normally distributed with $\langle \ln w \rangle = A^2 N_{S,cl}/2$ and $\sigma_{\ln w} = A^2 N_{S,cl}$.

This result only applies to the low-diffusion limit, but it may be used as a starting point to reconstruct $P(N_S)$, which can then be updated in Eq. (2.15), leading to a new estimate for $P(w)$, so on and so forth.

In practice, PYFPT first establishes whether or not $w$ is lognormally distributed, by applying D'Agostino and Pearson's normality test [99] to $\ln(w)$. This method gives the probability, or $p$-value, that a sample is drawn from an underlying normal distribution. If any $p$-values are smaller than a specified threshold (in practice we shall use 0.5%), the PYFPT package identifies that $w$ does not follow a lognormal distribution and alerts the user. If, however, the $p$-values are all greater than the threshold, then $P(w)$ is assumed to be lognormal. This implies a relationship between the mean of $w$ and the mean of its logarithm, namely

$$ \langle w \rangle = \exp \left( \langle \ln w \rangle + \frac{\sigma_{\ln w}^2}{2} \right). \qquad (2.16) $$

Crucially, the statistical reconstruction of $\langle \ln w \rangle$ is more robust than the one of $\langle w \rangle$. Indeed, in the later, one is dominated by those few occurrences that have a large value of $w$, while all runs with negligible values of $w$ are indistinguishable from $w = 0$ when
evaluating $\langle \widehat{w}_j \rangle = \sum_{i=0}^{n_j-1} w_j^{(i)}$. This is the sampling problem mentioned above, when $w_j^{(i)}$ spans several orders of magnitude. This is however not an issue for $\langle \ln w \rangle$, since $\ln w_j^{(i)}$ does not cover several orders of magnitude. The strategy is therefore the following [98]: estimate $\langle \ln w_j \rangle$ and $\sigma_{\ln w_j}$ from performing ensemble averages in the set of $\ln w_j^{(i)}$ values, and deduce $\langle w_j \rangle$ from Eq. (2.16). One can then replace Eq. (2.10) with

$$ \hat{P}(\mathcal{N}_j) = \frac{n_j \langle \widehat{w}_j \rangle}{n_{\text{total}}(\mathcal{N}_{j+1} - \mathcal{N}_j)}. \qquad (2.17) $$

Likewise, the uncertainty in $P(\mathcal{N}_j) = \hat{P}(\mathcal{N}_j) \pm \hat{\epsilon}_j^{\pm}$ is given by the standard error in estimating $\langle \ln w_j \rangle + \sigma_{\ln w_j}^2/2$ [100], propagated through the exponential of Eq. (2.16)

$$ \hat{\epsilon}_j^{\pm} = \frac{n_j \langle \widehat{w}_j \rangle}{n_{\text{total}}(\mathcal{N}_{j+1} - \mathcal{N}_j)} \left| \exp \left[ \pm \sqrt{\frac{(\hat{\sigma}_{\ln w_j})^2}{n_j} + \frac{(\hat{\sigma}_{\ln w_j})^4}{2n_j - 2}} \right] - 1 \right| \qquad (2.18) $$

Using the lognormal estimators (2.17) and (2.18) leads to more robust estimates in the low-diffusion regime. It greatly reduces the number of simulation runs required to obtain a specified accuracy, as $\hat{P}(\mathcal{N}_j)$ now benefits from all the runs in the jth bin, rather than just the few ones with the largest weights $w_j^{(i)}$. It also prevents the systematic error in the estimation of $P(\mathcal{N})$ from using the naïve estimate, as well as improves the uncertainty estimation, as can be checked explicitly in the right panel of Fig. 1.

## 3 Quadratic inflation

To test our implementation of importance sampling using the PYFPT package, we will compare our numerical method against analytical results. For stochastic slow-roll inflation, only a few analytic solutions of the first-passage time problem are known in the literature, among which are the quantum well, with a constant potential over a finite interval, and inflation driven by a quadratic potential [51]. The quantum well is dominated by quantum diffusion at every stage so the benefit in using importance sampling is not as striking as in quadratic inflation, on which we will thus focus.

The scalar field potential for quadratic inflation is sketched in Fig. 2 and reads

$$ V(\phi) = \frac{m^2 \phi^2}{2}. \qquad (3.1) $$

In the last few e-folds of inflation, in which we will be mainly interested, three regimes can distinguished: diffusion dominates over the drift when $m \gtrsim M_{\text{Pl}}$, there is an intermediate regime when $m \lesssim M_{\text{Pl}}$, and the drift dominates when $m \ll M_{\text{Pl}}$. These three cases will be investigated in Secs. 3.1, 3.2 and 3.3 respectively.

Although quadratic inflation driven by the simple potential (3.1) is not consistent with CMB anisotropies observed on large scales [13], it may still describe the last e-folds of inflation, when the inflaton approaches a simple minimum of the potential with
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Sketch of the potential (3.1) for quadratic inflation. The initial value of the field is given by $\phi_{\text{in}} > \phi_{\text{end}}$ and inflation ends at $\phi_{\text{end}} = \sqrt{2}M_{\text{Pl}}$, when the first slow-roll parameter becomes large, i.e. when $\epsilon_1 = -\dot{H}/H^2 = 1$. We introduce a high-energy UV reflective boundary at $\phi_{\text{UV}}$, shown in dashed grey.</figcaption></figure>

effective mass $m$ at the minimum. The runs are therefore started at an initial field value $\phi_{\text{in}} = \sqrt{42}M_{\text{Pl}}$, which corresponds to $N_{\text{cl}} = 10$ e-folds before the end of inflation in the drift-dominated limit. Obviously, values of $m$ close to or larger than the Planck mass are not allowed (since they would imply that inflation proceeds at super-Planckian energies), so here we consider the cases $m \gtrsim M_{\text{Pl}}$ and $m \lesssim M_{\text{Pl}}$ only to test our code in regimes where stochastic noise dominates. Indeed, there are models for which quantum diffusion dominates even at sub-Planckian energies (see e.g. Refs. [53, 88]), for which this analysis is relevant. Inflation ends by slow-roll violation once the field reaches $\phi_{\text{end}} = \sqrt{2}M_{\text{Pl}}$ and the first-passage times are recorded. Semi-analytical results for the resulting PDF for first-passage times that will be used as benchmarks for the code are presented in Appendix A.

Stochastic diffusion can result in the field climbing arbitrarily far back up the potential. In order to prevent the field from exploring too far up the potential, one may introduce a reflecting boundary at a finite field value, $\phi_{\text{UV}}$. Formally, this boundary can be removed by sending the parameter $\phi_{\text{UV}}$ to infinity, which in the case of single-field quadratic inflation does not lead to any divergence in the FPT statistics [96, 97]. In the following we will numerically investigate the dependence of PDF of first-passage times on $\phi_{\text{UV}}$.

To emphasise the efficiency of PYFPT, only $2 \times 10^5$ runs will be used for each
<figure><img src="image_4.png" /><figcaption>**Figure 3.** First-passage time PDF, $P(\mathcal{N})$, for quadratic inflation (3.1) with $m = M_{\text{Pl}}$ and $\phi_{\text{UV}} = \phi_{\text{in}} = \sqrt{42}M_{\text{Pl}}$. The red, blue and orange data points correspond to bias amplitudes $\mathcal{A} = 0$ (direct sampling), $\mathcal{A} = 0.4$ and $\mathcal{A} = 0.8$ respectively. The solid black curve is given by numerically calculating Eq. (A.19) for the first 50 poles. The dotted gray lines show the maximum $\mathcal{N}$ achieved for the direct and $\mathcal{A} = 0.4$ data sets.</figcaption></figure>

simulation set shown. The longest simulation time, required for drift-dominated cases, was ~ 8 hours using a laptop with a quad-core processor at 1.8 GHz base frequency. The longer simulation time for the sub-dominant diffusion case is due to the requirement that the step size $\Delta N$ must be much smaller than the standard deviation of $\mathcal{N}$. For diffusion-dominated and interim cases however, the simulations only took ~2 minutes and ~8 minutes respectively.

## 3.1 Diffusion domination: $m \gtrsim M_{\text{Pl}}$

As a first test of the importance sampling method we investigate the diffusion-dominated regime, $m \gtrsim M_{\text{Pl}}$. The probability of substantial deviations from the mean first-passage time is large for this case, therefore allowing us to investigate large deviations from a Gaussian distribution in the tail of the PDF, and in particular exponential tails in the presence of a finite boundary at $\phi_{\text{UV}}$ [83]. As will be further discussed below, the naïve estimator (2.10) is used in the figures of this subsection.

In Fig. 3 we show the results of numerical simulations of stochastic inflation for the first-passage times distribution $P(\mathcal{N})$, when $m = M_{\text{Pl}}$ and $\phi_{\text{in}} = \phi_{\text{UV}}$. For a direct simulation using $2 \times 10^5$ runs we are only able to probe the PDF down to $P(\mathcal{N}) \sim 10^{-3}$. Beyond this point there are too few runs in each bin ($n_j < 400$) to be able to reliably estimate the statistical error using jackknife resampling. However by introducing a bias we are able to map the PDF down to $P(\mathcal{N}) \sim 10^{-6}$ with a bias amplitude of $\mathcal{A} = 0.4$,
<figure><img src="image_5.png" /><figcaption>**Figure 4.** First-passage time PDF, $P(N)$, for different locations of the reflective boundary $\phi_{UV}$, in quadratic inflation (3.1) with $m = M_{Pl}$, $\phi_{in} = \sqrt{42}M_{Pl}$ and bias amplitude $A = 0.4$. The semi-analytical, solid black curves only contains the next-to-leading pole $\Lambda_1$ in Eq. (1.6), and the dashed black curve corresponds to the limit $\phi_{UV} \to \infty$, which is obtained by numerically evaluating Eq. (A.28).</figcaption></figure>

while data for $A = 0.8$ reach down to $P(N) \sim 10^{-12}$, using the same total number of runs ($2 \times 10^5$) in each case. This shows the numerical efficiency of the importance sampling method, and that as the bias amplitude $A$ is increased, one probes further into the tail, up to a numerical limit discussed later.

As expected, all three simulations show large deviations from a Gaussian distribution and an exponential tail for $P(N)$. The numerical results accurately reproduce the semi-analytical curve shown in Fig. 3, which corresponds to keeping the first 50 poles in Eq. (1.6), see Appendix A. For $10 \ll N \le 50$, it turns out to be dominated by the next-to-leading order pole $\Lambda_1$, and is well-described by a simple exponential tail. The leading-order pole $\Lambda_0$ is seen further into the tail at $N > 60$, beyond the regime that we can probe with $2 \times 10^5$ runs.

In Fig. 4 we show the effect of varying the position of the reflecting boundary $\phi_{UV}$. Analytically, we find that the distance between the poles in Eq. (1.6) decreases as $\ln^{-2}(\phi_{UV})$, while the leading pole $\Lambda_0$ is independent of $\phi_{UV}$ (see Appendix A). This implies that the exponential tail, i.e. the regime dominated by a single pole, occurs further into the tail with increasing $\phi_{UV}$. This is exactly what is seen in Fig. 4, where our numerical results agree, within the estimated error, with the expected exponential tail (where the first two poles are kept) for sufficiently large $N$. For $\phi_{UV} = 1.2 \phi_{in}$ and $\phi_{UV} = 1.4 \phi_{in}$, one can also see that the neighbourhood of the peak of the distribution
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Left panel: a 2D histogram of the logarithm of the weights, $\log_{10}(w)$, with first-passage times $N$ for the $A = 0.8$ data shown in Fig. 3 before truncation of the bins with less than 400 runs. Right panel: the standard deviation of $\ln(w)$ in each of the $N$ bins after truncation for the two importance samples shown in Fig. 3.</figcaption></figure>

is well reproduced by the $\phi_{UV} \to \infty$ limit, see Appendix A, and the degeneracy between different values of $\phi_{UV}$ is only lifted far in the tail. This is because the reflective boundary only affects the runs that are reflected against it, i.e. those that are sufficiently long to explore the large-field regions of the potential.

Note that we find the D'Agostino and Pearson's normality test to always fail in this regime $m \gtrsim M_{PI}$, where almost all the $p$-values are found to be below the 0.5% threshold. This means that the weights do not have a lognormal distribution, as illustrated in the left panel of Fig. 5. This is because quantum diffusion is not sub-dominant in this case, and it implies that the naïve estimator (2.10) for reconstructing the first-passage times PDF should be used, see Sec. 2.4.

Moreover, we find that the variance of the logarithms of the weights, $\ln(w)$, increases with $N$ and with the bias amplitude $A$, as shown in the right panel of Fig. 5. Therefore, as $A$ is increased to probe further into the tail of $P(N)$, the accuracy of the naïve estimator (2.10) decreases, as discussed in Sec. 2.4. This explains why this method does not allow us to reconstruct arbitrarily far regions of the tail with a fixed number of runs.

## 3.2 Interim case: $m \lesssim M_{PI}$

The interim mass case is phenomenologically interesting, as it produces a PDF, $P(N)$, with a Gaussian peak and non-Gaussian tail. This is therefore a regime where stochastic effects might be expected to predict an increased PBH production without disturbing the predictions of linear perturbation theory. While one can use the third and fourth central moments of $N$ in sub-dominant diffusion regime, calculated using the saddle-point approximation [68], to describe corrections to the Gaussian distribution in the near-tail behaviour using the Edgeworth expansion, see Appendix B for details, the far tail has no analytical prediction. This is because central moments beyond the fourth
<figure><img src="image_7.png" /><figcaption>**Figure 6.** First-passage time PDF, $P(\mathcal{N})$, for two locations of the reflective boundary, $\phi_{UV} = \phi_{in} = \sqrt{42}M_{Pl}$ and $\phi_{UV} \to \infty$, in quadratic inflation (3.1) with interim mass $m = 0.1 M_{Pl}$ and bias amplitude $A = 1.5$. The Gaussian and Edgeworth approximations (see Appendix B) are displayed with solid and dashed black curves, respectively.</figcaption></figure>

are increasingly difficult to calculate. Therefore, $m \lesssim M_{Pl}$ provides a useful test case for PYFPT, as we can see if it reproduces the analytic predictions at the peak of the PDF and for the near-tail before deviating further into the tail.

The far tail can be probed with importance sampling (even in the limit $\phi_{UV} \to \infty$), as diffusion is sub-dominant and the lognormal estimator for the PDF, Eq. (2.17), can be used. We verified that the p-values for the distributions of the logarithms of the weights, $\ln(w)$, in each bin are above the 0.5% threshold.

In Fig. 6 we show the PDF of first-passage times when $m = 0.1 M_{Pl}$ and with a bias amplitude $A = 1.5$, which allows us to probe down to $P(\mathcal{N}) \sim 10^{-10}$ using only $2 \times 10^5$ runs. The PDF is displayed for the two most extreme locations of the reflective boundary, $\phi_{UV} = \phi_{in}$ and $\phi_{UV} \to \infty$ (in the latter case, in practice, we set $\phi_{UV} = 100 \phi_{in}$ and find no simulations ever reach this high boundary). Any other choice for $\phi_{UV}$ leads to results that lie between these two extremes. While there appears to be a slight shift to larger $P(\mathcal{N})$ at large $\mathcal{N}$ for $\phi_{UV} \to \infty$ compared to $\phi_{UV} = \phi_{in}$, the data sets are broadly similar, suggesting that the dependence on $\phi_{UV}$ appears at even larger $\mathcal{N}$ [hence smaller $P(\mathcal{N})$] for $m \lesssim M_{Pl}$. One may also see that the logarithmic slope is not constant even at large $\mathcal{N}$, which means that several poles still contribute to the PDF, and the asymptotic exponential tail has not yet been reached. We finally note that both simulations agree within the estimated error with the Gaussian and Edgeworth approximations until $\mathcal{N} \sim 10.6$ and $\mathcal{N} \sim 11.25$ respectively, as expected.

In Fig. 7, the dependence of the result on the bias amplitude $A$ is explored, in terms
<figure><img src="image_8.png" /><figcaption>**Figure 7.** 2D histograms for the number of simulation runs in each $log_{10}(w)$ and $N$ bin represented as contours for quadratic inflation (3.1) with $m = 0.1 M_{\text{Pl}}$, $\phi_{\text{in}} = \sqrt{42} M_{\text{Pl}}$ and bias amplitudes $A = 1$ (blue), $A = 2$ (orange) and $A = 3$ (purple). The UV boundary used is $\phi_{\text{UV}} \to \infty$. Going outwards, the contours represent where 10, 50 and 90 percent of the runs lie, respectively.</figcaption></figure>

of the distributions of $N$ and of the weights $w$. When increasing $A$, one can see that the size of the contours only mildly grow, the main effect being that the location on which they are centred shifts to larger $N$ and smaller $w$. This confirms that, as discussed around Eq. (2.13), different bias amplitudes allow one to probe different regions of the PDF, since the sample process peaks around different values of $N$. This suggests a simple method for reconstructing $P(N)$ over a wide range, where simulations with different bias amplitudes are combined. This contrasts with the diffusion-domination case where we found that the contours become substantially wider when increasing $A$, see Fig. 5. Therefore, when increasing the bias, one indeed can probe further into the tail even as the number of runs is fixed, but at the expense of mildly diluting the data points into a wider region of $N$.

The multi-bias reconstruction method is employed in Fig. 8, where the PDF of first-passage times is displayed for $m = 0.1 M_{\text{Pl}}$ and several bias amplitudes $A$. We have magnified the error bars by a factor of 10, otherwise they would not be visible by eye, which shows how efficient the reconstruction is. Where they overlap, we see that the simulations are consistent with one another within the estimated errors and they are in agreement with the Edgeworth expansion for small $N$, before the expected deviation in the far tail. This clearly demonstrates the efficiency of the importance sampling method, as we required only $10^6$ runs in total to investigate all the way down to $P(N) \sim 10^{-33}$,
<figure><img src="image_9.png" /><figcaption>**Figure 8.** First-passage time PDF, $P(\mathcal{N})$, for the data shown in Fig. 7, along with $\mathcal{A} = 0$ (red, direct sampling). To improve their visibility, the error bars have been magnified by a factor of 10, i.e. their length corresponds to 10 times the error in $\log_{10}[P(\mathcal{N})]$. The Gaussian and Edgeworth approximations (see Appendix B) are displayed as solid and dashed black curves, respectively.</figcaption></figure>

when a direct sample would have required at least $10^{29}$ simulations to reach this far into the tail. A direct sample with $2 \times 10^5$ runs was only able to investigate down to $P(\mathcal{N}) \sim 10^{-4}$ and $\mathcal{N} = 10.66$, which does not reach the point where there is a significant deviation from Gaussianity at $\mathcal{N} = 10.91$. Note that as in Fig. 6, the exponential tail is still not reached in Fig. 8 (i.e. the logarithmic slope is not yet constant at large $\mathcal{N}$), although the tail is explored down to $P(\mathcal{N}) \sim 10^{-33}$ and strong, non-perturbative deviations from Gaussianity are observed. This was also the case when $\phi_{UV} = \phi_{in}$ was investigated.

Let us also note that the size of the error bars still increases with $\mathcal{A}$ in Fig. 8, and the reason is twofold. First, the mild increase in the dispersion of the weights noticed in Fig. 7 leads to larger error, since $\hat{\epsilon}_j^{\pm}$, grows with $\hat{\sigma}_{\ln w_j}$ in Eq. (2.18). Second, when the bias increases the dispersion in $\mathcal{N}$ also increases, which means that fewer data points are found in each bin, and this also results in a larger error since $\hat{\epsilon}_j^{\pm}$ decreases with $n_j$ (and more precisely as $1/\sqrt{n_j}$ if $n_j$ is large) in Eq. (2.18). However, let us stress that the two effects have different implications for the quality of the reconstruction: spreading the weight leads to increased error bars, while spreading $\mathcal{N}$ simply leads to diluting the information across a wider region of the PDF.

### 3.3 Drift domination: $m \ll M_{Pl}$

Let us now turn to the drift-dominated regime where $m \ll M_{Pl}$. We find that a Gaussian distribution remains an excellent approximation for $P(\mathcal{N})$ up to 15 standard deviations
<figure><img src="image_10.png" /><figcaption>**Figure 9.** First-passage time PDF, $P(\mathcal{N})$, for $\phi_{\text{UV}} = \phi_{\text{in}}$, in quadratic inflation (3.1) with $m = 0.001 M_{\text{Pl}}$, $\phi_{\text{in}} = \sqrt{42} M_{\text{Pl}}$ and bias amplitude $A = 0$ (red, direct sampling), $A = 1$ (blue), $A = 2$ (orange), $A = 3$ (purple) and $A = 4$ (green). The Gaussian approximations (see Appendix B) is displayed with a solid black curve.</figcaption></figure>

away from the mean (for $m = 0.001 M_{\text{Pl}}$ and $\langle \mathcal{N} \rangle = 10$). As diffusion is even more subdominant than in the interim-mass case discussed above, the D'Agostino and Pearson's test is always passed and the lognormal estimator for $P(\mathcal{N})$ is employed. The multi-bias reconstruction method can also be used. As the interim-mass regime gave results with little dependence on the choice of $\phi_{\text{UV}}$, we expect that the drift-dominated case will be approximately independent of the UV cutoff (which we have verified in practice) and therefore only the results obtained for $\phi_{\text{UV}} \to \infty$ are shown.

In Fig. 9 we show the results for a series of simulations with varying bias amplitude $A$. The simulations are consistent with one another when they overlap and the Gaussian approximation within the estimated error. The PyFPT package is able to accurately reproduce $P(\mathcal{N})$ from the peak of the distribution down to $P(\mathcal{N}) \sim 10^{-39}$ at 13.8 standard deviations away from the mean, using only $10^6$ runs (while direct sampling would require at least $10^{39}$ runs, and would take more than the age of the universe to complete). This demonstrates the accuracy of PyFPT and suggests that it can be used for exploring the very far tail of the distribution in the drift-dominated regime.

## 3.4 Bias optimization

So far we have used the bias function $f(\phi) = S(\phi) = H(\phi)/(2\pi)$, and have considered how the region of the PDF being sampled can be selected by hand tuning the bias amplitude $A$. However, in the drift-dominated case, the reconstruction is still numerically
<figure><img src="image_11.png" /><figcaption>Figure 10. Left panel: scatter plot in $log_{10}(w)$ and $N$ for quadratic inflation (3.1) with $m = 0.001 M_{\text{Pl}}$ and different powers of $p$ in the bias function (3.2). Only $2 \times 10^3$ runs are shown for clarity. The value of $A_p$ is tuned such that the mean of the sample distribution, $\langle N \rangle_S$, is 8 standard deviations away from the mean of the target distribution. Right panel: comparison of the first-passage time PDFs, $P(N)$, for the two data sets. The $p = 1$ set uses the lognormal estimator, see Sec. 2.4, while the $p = 3$ set is analysed with the naïve estimator, see Sec. 2.2. Both simulation sets use $4 \times 10^4$ runs.</figcaption></figure>

expensive. This is because, when solving the Langevin equation, one needs to work with a time step $\Delta N$ that is smaller than the standard deviation of $N$, and that standard deviation becomes tiny when quantum diffusion is subdominant (in the present model it grows linearly with $m$). Therefore many steps are required. This prompts us to explore whether further efficiency can be gained by optimising the bias function, $f(\phi)$, for the inflation model under investigation. This is why in this section we consider a more general, power-law form of the bias function

$$ \mathcal{B}(\phi) = A_p \phi^p. \qquad (3.2) $$

The previous bias function for quadratic inflation (3.1) corresponds to the choice $p = 1$, since from Eq. (1.4), $H(\phi) \propto \phi$.

In Fig. 10, we compare the previous bias function ($p = 1$) with the case where $p = 3$, which we will see below is the optimal choice. The left panel displays the scatter in $N$ and $w$ for $2 \times 10^3$ data points, and one can see that the choice $p = 3$ leads to a distribution that is substantially more squeezed, hence values that are more correlated. This implies that, within a given $N$ bin, the dispersion of the weights is smaller. This leads to reduced errors, as confirmed in the right panel of Fig. 10 where the PDFs corresponding to these two values of $p$ are displayed. Note that the typical values of $\ln(w)$ are smaller for $p = 1$ than for $p = 3$ in the left panel. This is because, to (eventually) converge on the same $P(N)$ in each bin, Eq. (2.16) indicates that $\langle \ln w \rangle$ must be smaller if $\sigma_{\ln w}$ is larger.

The dispersion of the weights is so small in the $p = 3$ case that the variance in $\ln w$, $\sigma_{\ln w}$, is difficult to resolve. This means that the lognormal estimate for the error, Eq. (2.18), is difficult to evaluate, which explains why we rather use the naïve estimate
<figure><img src="image_12.png" /><figcaption>**Figure 11.** The modulus of Pearson's correlation coefficient [101], $|r|$, of $log_{10}(w)$ with $N$ for different values of the power $p$ appearing in the bias function given in Eq. (3.2), for quadratic inflation (3.1) with $m = 0.001 M_{\text{P1}}$. The dashed gray line shows $|r| = 1$, corresponding to perfect correlation between $log_{10}(w)$ and $N$.</figcaption></figure>

for $p = 3$. The two estimates anyway converge for vanishingly small $\sigma_{\ln w}$. The right panel of Fig. 10 was produced using only $4 \times 10^4$ runs in each simulation. For $p = 3$ we obtain results that are as good as the lognormal data points shown in the right panel of Fig. 1, for a simulation with $2 \times 10^5$ runs. For $p = 1$, larger error bars are obtained.

This suggests a method for optimising the value of $p$, by measuring how much the $(N, w)$ contours are squeezed. This can be done by computing Pearson's correlation coefficient [101] between $log_{10}(w)$ and $N$, denoted $r$, which vanishes for uncorrelated variables and equals 1 and -1 for fully correlated and anti-correlated variables respectively. The result is displayed in Fig. 11, and confirms that the optimal value is indeed around $p \simeq 3$.

When the bias function is optimised, the different estimators (naive and lognormal) give similar results and can be used equivalently. However, the bias-optimisation programme may not be tractable beyond simple models like the one investigated here, so in general it remains important to choose the estimator appropriately.

# 4 Conclusions

In this work, we have shown how importance sampling can be used to investigate rare cosmological fluctuations in stochastic inflation. We have developed the publicly available PYFPT package, which reconstructs the far tail of the probability distribution of
first-passage times for general one-dimensional Langevin processes. Importance sam-
pling is achieved by introducing a bias into the simulated Langevin equation, to increase
the occurrence of rare realisations of the stochastic process. As by the stochastic-δN
formalism the local duration of inflation in e-folds, N, is directly related to the curvature
perturbation by ζ = N − ⟨N⟩, finding the PDF of ζ (and hence of N) corresponds to
solving a first-passage time problem.

Testing the robustness of importance-sampling methods for stochastic inflation was
our main objective. To this end, quadratic inflation was used as a benchmark model
that provides analytical or semi-analytical comparisons to the numerical results. This
confirmed the accuracy of the PYFPT package within the estimated errors. More pre-
cisely, we found that for diffusion domination (obtained within quadratic inflation by
setting m ≳ MPl), with 2 × 10⁵ simulation runs, importance sampling allows one to
reconstruct the PDF of N down to P(N) ~ 10⁻¹², while direct sampling is only able
to reach P(N) ~ 10⁻³, see Fig. 3. In the interim case where m ≲ MPl, and in the
drift-domination regime (m ≪ MPl), importance sampling allows the PDF down to
P(N) ≲ 10⁻³³ to be simulated, see Figs. 8 and 9 respectively. Randomly realising just
one of these events without importance sampling would take longer than the age of the
universe with current super computers. Here, only a laptop CPU was used.

The dependence of the statistics of rare fluctuations on the UV reflective bound-
ary, φUV, was also investigated. In the diffusion-dominated regime, the first-passage
time distribution carries a significant dependence on φUV, although this dependence is
removed to further into the tail as φUV increases, see Fig. 4. The asymptotic exponential
tail of the distribution, i.e. the part dominated by the leading pole of the characteristic
function, is also found deeper into the tail as φUV increases, and the seemingly expo-
nential tails in Fig. 4 are in fact given by the next-to-leading pole. In the interim and
drift-dominated regimes, no dependence on φUV is observed in the range probed [namely
down to P(N) ~ 10⁻³³, see Fig. 6]. In the interim case, we still observe non-perturbative
deviations from Gaussianity in the tail, although they are not of the simple exponential
type, denoting simultaneous contributions from multiple poles. This also indicates that,
in regimes where quantum diffusion does not dominate, it may not be enough to approx-
imate the non-Gaussian tail with an exponential fit. In the drift-dominated regime, we
find no significant deviation from Gaussianity down to probabilities P(N) ~ 10⁻³⁹, see
Fig. 9.

By tuning the bias, one can select which part of the distribution is to be best
reconstructed. In the interim and drift-dominated cases, this can be used to probe
almost arbitrarily far in the tail<fnref n="2" />. This is because the distributions of the weights, the
relative probability of realising importance sample without the bias, were found to be
given by lognormal distributions. This allows one to design estimators of the PDF that
benefit from all realisations within a given bin, even those with small weight.

We also found that the functional form of the bias could be optimised to the po-
tential investigated and have solved this optimisation problem for quadratic inflation,

<footnote n="2">²We could have investigated as far as P(N) ~ 10⁻⁶⁵, but beyond this point there would be less than
one such Hubble patch in the entire observable universe.</footnote>
where the variance in the weights within each $N$ bin was minimized (see Fig. 10). This suggests that for numerically expensive cases, such as drift domination, further efficiency can be gained by choosing a suitable form for the bias. This also has the additional advantage of making lognormal estimator unnecessary, as the naïve method is sufficient for vanishing weight variance. Of course this requires an optimized form for the bias, which is not in general known, while the lognormal estimator can be employed for all drift dominated cases. Finding the optimal bias function in a generic setup is an open issue, which would potentially allow us to increase the efficiency of importance sampling, and requires further investigation.

In the diffusion-dominated case, the distribution of the weights is more involved and cannot be modeled with a lognormal law. This implies that a more “naïve” estimator must be used, which is dominated by the few runs with largest weight. This leads to larger statistical error, so in practice one is limited in how far into the tail one can go.

The next step is to extend PYFPT to solve multi-dimensional first-passage time problems. This would allow investigations beyond the single-field slow-roll approximation to be done, such as phases of ultra-slow-roll [54–56] or multiple-field setups [84, 85, 93, 96, 97, 102–105]. In principle, the noise amplitude is determined by the field mode equation when integrated along a given realisation, and this leads to non-Markovian effects [88, 106] that we could also incorporate, together with other recently-proposed refinements of the stochastic-inflation equations [89, 107].

Finally we note that to make direct contact with observational constraints such as the abundance of primordial black holes, one still needs to map the distribution of first-crossing times at the end of inflation to the compaction function at a given scale, when primordial curvature perturbations are highly non Gaussian [81, 82, 108, 109]. This important step is left for future work.

## Data Availability Statement

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

## Acknowledgments

We would like to thank Andrew Gow, Chris Pattison, Lucas Pinol, Sébastien Renaux-Petel, Lukas Witkowski and the participants of the Paris-Portsmouth mini workshop for insightful discussions. JJ would like to thank Ian Harry and Coleman Krawczyk for help in developing the code into a package, as well as APC Paris for their hospitality during his Long Term Attachment. This work was supported by the Science and Technology Facilities Council (grant numbers ST/T506345/1, ST/S000550/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.
## A Semi-analytical results

To test the importance-sampling method, it is useful to compare its result with analytical or semi-analytical predictions. In this appendix, we derive such predictions for the simple model studied in this work, namely quadratic inflation (3.1).

Following the approach laid out in Ref. [51], we first introduce the characteristic function

$$ \chi_N(t, \phi) = \int_{-\infty}^{\infty} e^{itN} P_{\phi}(N) dN, \qquad (A.1) $$

where $t$ is a dummy variable. In stochastic slow-roll inflation, the characteristic function obeys the adjoint Fokker-Planck equation [51]

$$ \left( \frac{\partial^2}{\partial \phi^2} - \frac{v'}{v^2} \frac{\partial}{\partial \phi} + \frac{it}{v M_{\text{Pl}}^2} \right) \chi_N(t, \phi) = 0, \qquad (A.2) $$

where the reduced potential

$$ v(\phi) = \frac{V(\phi)}{24\pi^2 M_{\text{Pl}}^4} \qquad (A.3) $$

was introduced to simplify notations. Due to the absorbing boundary at $\phi_{\text{end}}$ and the reflective boundary at $\phi_{\text{UV}}$, the FPT distribution must be such that $P_{\phi_{\text{end}}}(N) = \delta(N)$ and $\frac{\partial}{\partial \phi} P_{\phi_{\text{UV}}}(N) = 0$, hence the characteristic function must satisfy the boundary conditions

$$ \chi_N(t, \phi_{\text{end}}) = 1 \quad \text{and} \quad \frac{\partial \chi_N}{\partial \phi}(t, \phi_{\text{UV}}) = 0. \qquad (A.4) $$

In general, the solution to Eq. (A.2) under the conditions (A.4) can be cast in the form [83]

$$ \chi_N(t, \phi) = \sum_n \frac{a_n(\phi)}{\Lambda_n - it} + g(t, \phi), \qquad (A.5) $$

where $a_n(\phi)$ and $\Lambda_n$ are the $n^{th}$ residual and pole respectively, and $g(t, \phi)$ is a regular function. By inverse Fourier transforming Eq. (A.1), one obtains

$$ P_{\phi}(N) = \frac{1}{2\pi} \int_{-\infty}^{\infty} e^{-itN} \chi_{N}(t, \phi) dt. \qquad (A.6) $$

Computing the integral with the residue theorem, one then finds

$$ P_{\phi}(N) = \sum_{n} a_{n}(\phi) e^{-\Lambda_{n} N}, \qquad (A.7) $$

i.e. $P_{\phi}(N)$ is a sum of decaying exponentials.

Let us now apply this program to quadratic inflation,

$$ v = v_0 \left( \frac{\phi}{M_{\text{Pl}}} \right)^2, \qquad (A.8) $$
see Eq. (3.1), where inflation takes place between $\phi_{\text{end}} = \sqrt{2}M_{\text{Pl}}$ (where inflation ends by slow-roll violation, hence an absorbing wall is set) and $\phi_{\text{UV}}$ (where a reflective wall is set). Substituting the potential function (A.8) into Eq. (A.2), one obtains

$$ \left( \frac{\partial^2}{\partial \phi^2} - \frac{2M_{\text{Pl}}^2}{v_0 \phi^3} \frac{\partial}{\partial \phi} + \frac{it}{v_0 \phi^2} \right) \chi(t, \phi) = 0, \qquad (A.9) $$

the generic solution of which is given by

$$ \chi(t, \phi) = \lambda(t)f(t, \phi) + \mu(t)g(t, \phi), \qquad (A.10) $$

where

$$ \begin{aligned} f(t, \phi) &= [v(\phi)]^{\frac{1+\alpha(t)}{4}} {}_1F_1\left[\frac{-1-\alpha(t)}{4}; 1-\frac{\alpha(t)}{2}; -\frac{1}{v(\phi)}\right], \\ g(t, \phi) &= [v(\phi)]^{\frac{1-\alpha(t)}{4}} {}_1F_1\left[\frac{-1+\alpha(t)}{4}; 1+\frac{\alpha(t)}{2}; -\frac{1}{v(\phi)}\right]. \end{aligned} \qquad (A.11) $$

Here, $\lambda(t)$ and $\mu(t)$ are two integration constants, and to simplify notations we have introduced the parameter

$$ \alpha(t) = \sqrt{1 - 4 \frac{it}{v_0}}. \qquad (A.12) $$

The boundary conditions given in Eq. (A.4) imply that

$$ \mu = \frac{B}{B-A} \quad \text{and} \quad \lambda = -\frac{A}{B-A}, \qquad (A.13) $$

where

$$ A \equiv \frac{f'_{\text{UV}}}{g'_{\text{UV}}} \quad \text{and} \quad B \equiv \frac{f_{\text{end}}}{g_{\text{end}}}. \qquad (A.14) $$

Here, $f_{\text{end}}$ is a shorthand for $f(t, \phi_{\text{end}})$ etc, and a prime denotes a derivation with respect to the field value. The characteristic function can thus be written as

$$ \chi(t, \phi) = \frac{B(t)}{B(t) - A(t)} \frac{f(t, \phi)}{f_{\text{end}}(t)} - \frac{A(t)}{B(t) - A(t)} \frac{g(t, \phi)}{g_{\text{end}}(t)}. \qquad (A.15) $$

## Numerical reconstruction

Two methods can then be used to reconstruct the FPT distribution function. The first one consists in performing the integral in Eq. (A.6) numerically. The second method relies on numerically solving the pole equation

$$ A(t) = B(t) \qquad (A.16) $$

and derive the corresponding values of the poles $\Lambda_n = it_n$. The residues can then be obtained with the formula

$$ a_n(\phi) = \lim_{t \to -i\Lambda_n} (\Lambda_n - it) \frac{B(t) \frac{f(t, \phi)}{f_{\text{end}}(t)} - A(t) \frac{g(t, \phi)}{g_{\text{end}}(t)}}{B(t) - A(t)}, \qquad (A.17) $$
which can be further simplified using L'Hôpital's rule, leading to

$$a_n(\phi) = \lim_{t \to -i\Lambda_n} -i \frac{B(t) \frac{f(t,\phi)}{f_{\text{end}}(t)} - A(t) \frac{g(t,\phi)}{g_{\text{end}}(t)}}{\frac{\partial}{\partial t} [B(t) - A(t)]}. \quad (A.18)$$

The PDF can then be evaluated with Eq. (A.7), i.e.

$$P_{\phi}(\mathcal{N}) = \sum_{n} \lim_{t \to -i\Lambda_n} -i \frac{B(t) \frac{f(t,\phi)}{f_{\text{end}}(t)} - A(t) \frac{g(t,\phi)}{g_{\text{end}}(t)}}{\frac{\partial}{\partial t} [B(t) - A(t)]} e^{-\Lambda_n \mathcal{N}}. \quad (A.19)$$

In practice, only a finite number of poles $\Lambda_n$ for $n \in [0, n_{\text{max}}]$ is extracted with their corresponding residues. One then needs to check that $n_{\text{max}}$ is large enough for the reconstructed PDF to not depend on it. We find that such a number required for an accurate estimation depends not only on $m$ and $\phi_{\text{UV}}$ but also on how far into the tail one is approximating.

For diffusion domination with $m \gtrsim M_{\text{Pl}}$ and a low energy UV cutoff such as $\phi_{\text{UV}} \simeq \phi_{\text{in}}$, including the first 50 poles in Eq. (A.19) gives a well-converged approximation even for the peak of the distribution. Further in the tail, only the leading- and next-to-leading poles are needed. However, for the interim case with $m \lesssim M_{\text{Pl}}$ and high energy UV cutoff $\phi_{\text{UV}} \gg \phi_{\text{in}}$, thousands of poles are required for an accurate estimation even at $P(\mathcal{N}) \sim 10^{-20}$. As we shall now see, the reason is the distance between the poles decreases with $\phi_{\text{UV}}$, hence more poles contribute to the result as $\phi_{\text{UV}}$ increases, even deep in the tail.

## Large-$\phi_{\text{UV}}$ limit

The distance between two consecutive poles depends on the value of $\phi_{\text{UV}}$, and to show this property explicitly let us now study the large-$\phi_{\text{UV}}$ regime. The terms $f'_{\text{UV}}$ and $g'_{\text{UV}}$ can be approximated by expanding the hypergeometric functions in the limit where their last argument is close to zero while being negative, leading to<fnref n="3" />

$$\begin{aligned} f'_{\text{UV}} &= (v_{\text{UV}})^{\frac{1+\alpha}{4}} \frac{1+\alpha}{2} \frac{1}{\phi_{\text{UV}}} \left[ 1 + \mathcal{O}\left(\frac{1}{v_{\text{UV}}}\right) \right], \\ g'_{\text{UV}} &= (v_{\text{UV}})^{\frac{1-\alpha}{4}} \frac{1-\alpha}{2} \frac{1}{\phi_{\text{UV}}} \left[ 1 + \mathcal{O}\left(\frac{1}{v_{\text{UV}}}\right) \right], \end{aligned} \quad (A.20)$$

where $v_{\text{UV}} = v(\phi_{\text{UV}})$. Note that for this expansion to be valid, one needs $v_{\text{UV}}$ to be large, which implies $\phi_{\text{UV}} \gg v_0^{-1/2} \sim M_{\text{Pl}}^2/m$. For sub-Planckian values of $m$ such field values are very large, and presumably out of the validity range of the model, but here the large-$\phi_{\text{UV}}$ limit is worked out only as a formal way to understand how $\phi_{\text{UV}}$ affects the pole locations. These expressions lead to

$$A = \frac{f'_{\text{UV}}}{g'_{\text{UV}}} = \frac{1 + \alpha}{1 - \alpha} (v_{\text{UV}})^{\frac{\alpha}{2}} \left[ 1 + \mathcal{O}\left(\frac{1}{v_{\text{UV}}}\right) \right]. \quad (A.21)$$

<footnote n="3">³This is obtained by combining Eqs. (13.3.15), (13.2.39) and (13.2.13) of Ref. [110].</footnote>
To approximate $B = f_{\text{end}}/g_{\text{end}}$, let us further assume that $v_0 \ll 1$, hence $v_{\text{end}} \ll 1$ with $\phi_{\text{end}} = \sqrt{2}M_{\text{Pl}}$. By expanding the hypergeometric functions in the limit where their last argument is very large negative, one obtains$^4$

$$B = \frac{f_{\text{end}}}{g_{\text{end}}} \approx \frac{\Gamma\left(1 - \frac{\alpha}{2}\right) \Gamma\left(\frac{5+\alpha}{4}\right)}{\Gamma\left(1 + \frac{\alpha}{2}\right) \Gamma\left(\frac{5-\alpha}{4}\right)}. \qquad (\text{A.22})$$

In this regime the pole equation (A.16) reduces to

$$\frac{1+\alpha}{1-\alpha} (v_{UV})^{\frac{\alpha}{2}} = \frac{\Gamma\left(1-\frac{\alpha}{2}\right) \Gamma\left(\frac{5+\alpha}{4}\right)}{\Gamma\left(1+\frac{\alpha}{2}\right) \Gamma\left(\frac{5-\alpha}{4}\right)}. \qquad (\text{A.23})$$

This equation needs to be solved on the negative imaginary axis, i.e. when $t = -i\Lambda$, where $\Lambda$ is real positive. For $0 < \Lambda < v_0/4$, $\alpha$ is real and comprised between 0 and 1, see Eq. (A.12). In that case, it can then be shown that the pole equation has no solution if $v_{UV}$ is of order one or larger (which we assumed). Moreover, if $\Lambda = v_0/4$, then $\alpha = 0$ and the pole equation is trivially satisfied. This thus gives us the leading pole,

$$\Lambda_0 = \frac{v_0}{4}. \qquad (\text{A.24})$$

Note that when $\alpha = 0$, one has $f = g$ and it is clear that $A = B = 1$ exactly. The expression we just derived for the leading pole is therefore valid in general, beyond the large-$\phi_{\text{UV}}$ limit.

Let us now study the location of the first higher poles. When $\Lambda > v_0/4$, $\alpha$ is imaginary, thus both hands of the pole equation (A.23) oscillate. While the left-hand side oscillates in $\alpha$ with a frequency of the order $\ln(v_{UV})$, the right-hand side oscillates at a slower rate with a frequency of order one. If $\ln(v_{UV}) \gg 1$, the first oscillations in the left-hand side occur while the right-hand side is still approximately constant, so one may replace $\alpha \simeq 0$ in Eq. (A.23) except in the term $(v_{UV})^{\alpha/2}$. This leads to

$$(v_{UV})^{\frac{\alpha}{2}} = 1, \qquad (\text{A.25})$$

which, by replacing $t = -i\Lambda$ and substituting in Eq. (A.12), has solutions

$$\Lambda_n|_{n \ll \frac{\ln(v_{UV})}{4\pi}} \approx \frac{v_0}{4} \left[ 1 + \left( \frac{4\pi n}{\ln v_{UV}} \right)^2 \right]. \qquad (\text{A.26})$$

This approximation is valid as long as $|\alpha| \ll 1$, which implies that $n \ll \ln(v_{UV})/(4\pi)$. This is why only the first poles are approached by this formula, which is relevant only when $\ln(v_{UV})$ is large (hence $\phi_{\text{UV}}$ is tremendously large). It however makes explicit that the distance between poles, $\Lambda_{n-1} - \Lambda_n$, decreases with $\phi_{\text{UV}}$. This decrease is only logarithmic but it does show that as $\phi_{\text{UV}}$ increases, more poles contribute to the PDF. It also means that as $\phi_{\text{UV}}$ goes to infinity, the set of poles become continuous, hence the asymptotic tail is not exponential anymore, as we confirm below.

<footnote>$^4$This is obtained by combining Eqs. (13.2.39) and (13.2.23) of Ref. [110].</footnote>
## Infinite-$\phi_{UV}$ limit

The case where $\phi_{UV}$ is infinite is indeed peculiar and needs to be treated separately. When $\phi \to \infty$, $1/v \to 0$ and the hypergeometric functions in Eq. (A.11) tend to a constant. The behaviour of $f$ and $g$ in Eq. (A.11) is then dictated by the prefactors $v^{(1\pm\alpha)/4}$. Whether this diverges or not when $v \to \infty$ depends on the sign of $\Re(1+\alpha)$. We restrict this sign analysis to the case where $t$ is real, recalling that the PDF can be extracted using Eq. (A.6), which only relies on evaluating the characteristic function with real $t$.

Since $\alpha^2 = 1 - 4it/v_0$, one has $\Re(\alpha^2) = 1$ so can write $\alpha^2 = e^{i\theta}/\cos(\theta)$, where $-\pi/2 < \theta < \pi/2$. This gives rise to $\alpha = e^{i\theta/2}/\sqrt{\cos(\theta)}$, meaning $\Re(\alpha) = \sqrt{[1+1/\cos(\theta)]/2}$. When $-\pi/2 < \theta < \pi/2$, one thus has $\Re(\alpha) > 1$. Therefore, when $v$ goes to infinity, $f$ diverges while $g$ asymptotes a constant. In order to satisfy the second boundary condition in Eq. (A.4), one must thus keep the $g$ branch only and set $\lambda(t) = 0$. The remaining integration constant, $\mu(t)$, can then be set in order to satisfy the first boundary condition and this gives rise to

$$ \chi(t, \phi) = \left( \frac{v}{v_{\text{end}}} \right)^{\frac{1-\alpha}{4}} \frac{1 F_1 \left( \frac{-1+\alpha}{4}; 1+\frac{\alpha}{2}; -\frac{1}{v} \right)}{1 F_1 \left( \frac{-1+\alpha}{4}; 1+\frac{\alpha}{2}; -\frac{1}{v_{\text{end}}} \right)}. \qquad (A.27) $$

The PDF can then obtained from Eq. (A.6) (a numerical extraction of the poles is indeed impossible, given that they constitute a continuous spectrum when $\phi_{UV}$ is infinite). On the negative imaginary axis when $it > v_0/4$, the first derivatives of the characteristic function are discontinuous, i.e. there is a branch cut. Together with the fact that there are no poles, we can deform the contour of integration to rewrite Eq. (A.6) as

$$ P_{\phi}(\mathcal{N}) = -2 \int_{v_0/4}^{\infty} \Im[\chi(-iy, \phi)]e^{-y\mathcal{N}}dy. \qquad (A.28) $$

This integral can be computed numerically, but requires very high precision for $v_0 \ll 1$ due to the large argument in the confluent hypergeometric function.

Let us finally note that the far tail of the PDF at $\mathcal{N} \gg 1/v_0$ is determined by the behaviour of the imaginary part of the characteristic function around $y = v_0/4$. Expanding the characteristic function for small $|\alpha|$, one can show that $\Im[\chi(-iy, \phi)] \propto \Im[\alpha(-iy)]$ and $P_{\phi}(\mathcal{N}) = \mathcal{N}^{-3/2} \exp(-v_0\mathcal{N}/4)$. Therefore, as announced above, in the case where $\phi_{UV}$ is infinite the tail is not exponential anymore, but only quasi exponential, due to the presence of a continuous spectrum of poles.

## B Edgeworth expansion

In the main text, in order to describe the non-Gaussian features of the FPT distributions close to their maximum, we compare them to a Gaussian approximation and to an Edgeworth expansion, which parametrises the first deviations from Gaussian statistics. In this appendix, we recall how the Edgeworth expansion [111, 112] is constructed.
Let us consider a random variable X with vanishing mean and unit variance, and let {$X_1, X_2, \cdots, X_n$} be n independent copies of X. We want to determine the distribution function of the normalised sum

$$x = \frac{1}{\sqrt{n}} \sum_{i=1}^{n} X_{n}. \qquad (B.1)$$

Since the $X_i$ are independent, $x$ also has a vanishing mean and unit variance. By virtue of the central-limit theorem, we know that when $n \to \infty$, $x$ becomes a Gaussian random variable, with a vanishing mean and unit variance. This can be seen as the leading-order, Gaussian approximation. Our goal is to go beyond that leading-order result and describe the first non-Gaussian corrections when $n$ is large but finite.

To that end, we introduce the characteristic function of $x$,

$$ \begin{aligned} \chi(t, x) &= \langle \exp(itx) \rangle = \left\langle \exp\left(\frac{it}{\sqrt{n}} \sum_{i=1}^{n} X_i \right) \right\rangle \\ &= \Pi_{i=1}^{n} \left\langle \exp\left(\frac{it}{\sqrt{n}} X_i \right) \right\rangle = \chi^n\left(\frac{t}{\sqrt{n}}, X \right), \end{aligned} \qquad (B.2) $$

where we have used that the $X_i$ are independent, and where we recognise $\chi(t, X)$, the characteristic function of the $X$ variable.

Let us then Taylor expand $\chi(t/\sqrt{n}, X)$ in the limit where $n$ is large:

$$ \begin{aligned} \chi\left(\frac{t}{\sqrt{n}}, X\right) &= \left\langle \exp\left(\frac{it}{\sqrt{n}}X\right) \right\rangle \\ &= 1 - \frac{t^2}{2n} - i \frac{t^3}{6n\sqrt{n}} \langle X^3 \rangle + \frac{t^4}{24n^2} \langle X^4 \rangle + \mathcal{O}\left(\frac{1}{n^2\sqrt{n}}\right), \end{aligned} \qquad (B.3) $$

where we have used that $\langle X \rangle = 0$ and $\langle X^2 \rangle = 1$. Taking this expression to the $n^{th}$ power and further expanding in $n$ leads to

$$ \chi(t, x) = e^{-\frac{t^2}{2}} \left[ 1 - \frac{it^3}{6\sqrt{n}} \langle X^3 \rangle + \frac{t^4}{24n} (\langle X^4 \rangle - 3) - \frac{t^6}{72n} \langle X^3 \rangle^2 \right] + \mathcal{O}\left(\frac{1}{n\sqrt{n}}\right). \quad (B.4) $$

This can be cast in terms of the moments of $x$. Since $\langle X_i \rangle = 0$ and $\langle X_i X_j \rangle = \delta_{ij}$, Eq. (B.1) leads to $\langle x \rangle = 0$, $\langle x^2 \rangle = 1$, $\langle x^3 \rangle = \langle X^3 \rangle / \sqrt{n}$ and $\langle x^4 \rangle = 3 + (\langle X^4 \rangle - 3) / n$. One thus has

$$ \chi(t, x) = e^{-\frac{t^2}{2}} \left[ 1 - \frac{it^3}{6} \langle x^3 \rangle + \frac{t^4}{24} (\langle x^4 \rangle - 3) - \frac{t^6}{72} \langle x^3 \rangle^2 \right] + \mathcal{O}\left(\frac{1}{n\sqrt{n}}\right). \quad (B.5) $$

The PDF can then be obtained by Fourier transforming this expression, see Eq. (A.6), and this leads to

$$ P(x) \approx \frac{e^{-x^2/2}}{\sqrt{2\pi}} \left[ 1 + \frac{\langle x^3 \rangle}{6} He_3(x) + \frac{(\langle x^4 \rangle - 3)}{24} He_4(x) + \frac{\langle x^3 \rangle^2}{72} He_6(x) \right], \quad (B.6) $$
where

$$He_3(x) = x^2 - 1, \qquad (B.7)$$

$$He_4(x) = x^4 - 6x^2 + 3, \qquad (B.8)$$

$$He_6(x) = x^6 - 15x^4 + 45x^2 - 15, \qquad (B.9)$$

are the 3rd, 4th and 6th Hermite polynomials respectively. This expression can be readily generalised to the case where $x$ has a non-trivial mean $\mu$ and standard deviation $\sigma$, by applying the above to $\tilde{x} \equiv (x - \mu)/\sigma$. One thus finds

$$P(x) \approx \frac{e^{-\frac{(x-\mu)^2}{2\sigma^2}}}{\sqrt{2\pi}\sigma^2} \left[ 1 + \frac{\gamma}{6} He_3\left(\frac{x-\mu}{\sigma}\right) + \frac{\kappa}{24} He_4\left(\frac{x-\mu}{\sigma}\right) + \frac{\gamma^2}{72} He_6\left(\frac{x-\mu}{\sigma}\right) \right], \qquad (B.10)$$

where $\gamma$ and $\kappa$ are called “skewness” and “excess kurtosis” respectively, and are given by

$$\gamma = \frac{\langle (x - \mu)^3 \rangle}{\sigma^3} \quad \text{and} \quad \kappa = \frac{\langle (x - \mu)^4 \rangle}{\sigma^4} - 3. \qquad (B.11)$$

For a Gaussian random variable, $\gamma$ and $\kappa$ vanish, and one recovers indeed that $P(x)$ is a Gaussian function. Otherwise, the above expression allows one to parametrise small deviations from Gaussianity. It involves the first few moments of the PDF, which in our case can be computed exactly since the characteristic function is known analytically, using the formula

$$\langle \mathcal{N}^n \rangle = i^{-n} \frac{\partial^n}{\partial t^n} \chi(t, \phi) \Big|_{t=0} . \qquad (B.12)$$

We find that the Edgeworth expansion provides a reliable approximation in the drift-dominated regime close to the maximum of the distribution, and is therefore a useful comparison for the PYFPT package in the near tail of $P(\mathcal{N})$. However, it breaks down in the presence of exponential tails, where higher than 4th order moments yield significant contributions.

## References

[1] A. A. Starobinsky, *A New Type of Isotropic Cosmological Models Without Singularity*, Phys. Lett. B91 (1980) 99-102.

[2] K. Sato, *First Order Phase Transition of a Vacuum and Expansion of the Universe*, Mon. Not. Roy. Astron. Soc. 195 (1981) 467-479.

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

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

[5] A. Albrecht and P. J. Steinhardt, *Cosmology for Grand Unified Theories with Radiatively Induced Symmetry Breaking*, Phys. Rev. Lett. 48 (1982) 1220-1223.
[6] A. D. Linde, *Chaotic Inflation*, Phys.Lett. B129 (1983) 177-181.

[7] D. J. Eisenstein et al., SDSS-III: Massive spectroscopic surveys of the distant universe, the milky way, and extra-solar planetary systems, Astrophys. J. 142 (aug, 2011) 72, [1101.1529].

[8] C. Blake et al., The WiggleZ Dark Energy Survey: mapping the distance-redshift relation with baryon acoustic oscillations, MNRAS 418 (Dec., 2011) 1707-1724, [1108.2635].

[9] K. S. Dawson et al., The baryon oscillation spectroscopic survey of SDSS-III, Astrophys. J. 145 (Dec., 2012) 10, [1208.0022].

[10] S. Alam et al., The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, MNRAS 470 (03, 2017) 2617-2652, [1607.03155].

[11] PLANCK collaboration, P. A. R. Ade et al., Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. 594 (2016) A13, [1502.01589].

[12] D. Saadeh, S. M. Feeney, A. Pontzen, H. V. Peiris and J. D. McEwen, *How isotropic is the Universe?*, Phys. Rev. Lett. 117 (2016) 131302, [1605.07178].

[13] PLANCK collaboration, N. Aghanim et al., Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6, [1807.06209].

[14] V. F. Mukhanov and G. Chibisov, *Quantum Fluctuation and Nonsingular Universe*, JETP Lett. 33 (1981) 532-535.

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

[16] A. A. Starobinsky, *Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations*, Phys.Lett. B117 (1982) 175-178.

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

[18] S. Hawking, *The Development of Irregularities in a Single Bubble Inflationary Universe*, Phys.Lett. B115 (1982) 295.

[19] J. M. Bardeen, P. J. Steinhardt and M. S. Turner, *Spontaneous Creation of Almost Scale Free Density Perturbations in an Inflationary Universe*, Phys.Rev. D28 (1983) 679.

[20] PLANCK collaboration, P. Ade et al., Planck 2015 results. XX. Constraints on inflation, 1502.02114.

[21] D. H. Lyth and A. R. Liddle, *The primordial density perturbation: Cosmology, inflation and the origin of structure*. Cambridge University Press, Apr., 2009.

[22] A. R. Liddle and S. M. Leach, *How long before the end of inflation were observable perturbations produced?*, Phys. Rev. D 68 (2003) 103503, [astro-ph/0305263].

[23] 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, [1505.01834].

[24] N. Christensen, *Stochastic Gravitational Wave Backgrounds*, Rept. Prog. Phys. 82 (2019) 016903, [1811.08797].
[25] S. Hawking, *Gravitationally collapsed objects of very low mass*, Mon. Not. Roy. Astron. Soc. **152** (1971) 75.

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

[27] B. J. Carr, *The Primordial black hole mass spectrum*, Astrophys. J. **201** (1975) 1-19.

[28] J. García-Bellido, *Massive Primordial Black Holes as Dark Matter and their detection with Gravitational Waves*, J. Phys. Conf. Ser. **840** (2017) 012032, [1702.08275].

[29] LIGO SCIENTIFIC, VIRGO collaboration, B. P. Abbott et al., *GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs*, Phys. Rev. X **9** (2019) 031040, [1811.12907].

[30] LIGO SCIENTIFIC, VIRGO collaboration, R. Abbott et al., *GWTC-2.1: Deep Extended Catalog of Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run*, 2108.01045.

[31] LIGO SCIENTIFIC, VIRGO collaboration, R. Abbott et al., *Properties and astrophysical implications of the 150 Msr binary black hole merger GW190521*, Astrophys. J. Lett. **900** (2020) L13, [2009.01190].

[32] S. Clesse and J. Garcia-Bellido, *GW190425 and GW190814: Two candidate mergers of primordial black holes from the QCD epoch*, 2007.06481.

[33] LIGO SCIENTIFIC, VIRGO, KAGRA collaboration, R. Abbott et al., *First joint observation by the underground gravitational-wave detector*, KAGRA, with GEO600, 2203.01270.

[34] B. Carr, F. Kuhnel and M. Sandstad, *Primordial Black Holes as Dark Matter*, Phys. Rev. D **94** (2016) 083504, [1607.06077].

[35] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, *Constraints on Primordial Black Holes*, 2002.12778.

[36] B. Carr and F. Kuhnel, *Primordial Black Holes as Dark Matter: Recent Developments*, Ann. Rev. Nucl. Part. Sci. **70** (2020) 355-394, [2006.02838].

[37] A. M. Green and B. J. Kavanagh, *Primordial Black Holes as a dark matter candidate*, J. Phys. G **48** (2021) 043001, [2007.10722].

[38] M. Shibata and M. Sasaki, *Black hole formation in the Friedmann universe: Formulation and computation in numerical relativity*, Phys. Rev. D **60** (1999) 084002, [gr-qc/9905064].

[39] T. Harada, C.-M. Yoo, T. Nakama and Y. Koga, *Cosmological long-wavelength solutions and primordial black hole formation*, Phys. Rev. D **91** (2015) 084057, [1503.03934].

[40] A. Escrivà, C. Germani and R. K. Sheth, *Universal threshold for primordial black hole formation*, Phys. Rev. D **101** (2020) 044022, [1907.13311].

[41] I. Musco, *Threshold for primordial black holes: Dependence on the shape of the cosmological perturbations*, Phys. Rev. D **100** (2019) 123524, [1809.02127].

[42] S. Young, I. Musco and C. T. Byrnes, *Primordial black hole formation and abundance: contribution from the non-linear relation between the density and curvature perturbation*, JCAP **11** (2019) 012, [1904.00984].
[43] I. Musco, J. C. Miller and A. G. Polnarev, *Primordial black hole formation in the radiative era: Investigation of the critical nature of the collapse*, Class. Quant. Grav. **26** (2009) 235001, [0811.1452].

[44] J. Garcia-Bellido and E. Ruiz Morales, *Primordial black holes from single field models of inflation*, 1702.03901.

[45] J. M. Ezquiaga, J. Garcia-Bellido and E. Ruiz Morales, *Primordial Black Hole production in Critical Higgs Inflation*, Phys. Lett. B **776** (2018) 345-349, [1705.04861].

[46] C. Germani and T. Prokopec, *On primordial black holes from an inflection point*, 1706.04226.

[47] H. Motohashi and W. Hu, *Primordial Black Holes and Slow-roll Violation*, 1706.06784.

[48] G. Ballesteros and M. Taoso, *Primordial black hole dark matter from single field inflation*, Phys. Rev. D **97** (2018) 023501, [1709.05565].

[49] S. Rasanen and E. Tomberg, *Planck scale black hole dark matter from Higgs inflation*, JCAP **01** (2019) 038, [1810.12608].

[50] S. Geller, W. Qin, E. McDonough and D. I. Kaiser, *Primordial Black Holes from Multifield Inflation with Nonminimal Couplings*, 2205.04471.

[51] C. Pattison, V. Vennin, H. Assadullahi and D. Wands, *Quantum diffusion during inflation and primordial black holes*, JCAP **1710** (2017) 046, [1707.00537].

[52] M. Biagetti, G. Franciolini, A. Kehagias and A. Riotto, *Primordial Black Holes from Inflation and Quantum Diffusion*, 1804.07124.

[53] J. M. Ezquiaga and J. García-Bellido, *Quantum diffusion beyond slow-roll: implications for primordial black-hole production*, JCAP **08** (2018) 018, [1805.06731].

[54] H. Firouzjahi, A. Nassiri-Rad and M. Noorbala, *Stochastic Ultra Slow Roll Inflation*, JCAP **1901** (2019) 040, [1811.02175].

[55] G. Ballesteros, J. Rey, M. Taoso and A. Urbano, *Stochastic inflationary dynamics beyond slow-roll and consequences for primordial black hole formation*, JCAP **08** (2020) 043, [2006.14597].

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

[57] T. Prokopec and G. Rigopoulos, *ΔN and the stochastic conveyor belt of Ultra Slow-Roll*, 1910.08487.

[58] G. Rigopoulos and A. Wilkins, *Inflation is always semi-classical: diffusion domination overproduces Primordial Black Holes*, JCAP **12** (2021) 027, [2107.05317].

[59] A. A. Starobinsky, *Stochastic de Sitter (inflationary) stage in the early Universe*, Lect. Notes Phys. **246** (1986) 107-126.

[60] Y. Nambu and M. Sasaki, *Stochastic Stage of an Inflationary Universe Model*, Phys. Lett. **B**205 (1988) 441-446.

[61] Y. Nambu and M. Sasaki, *Stochastic Approach to Chaotic Inflation and the Distribution of Universes*, Phys. Lett. **B**219 (1989) 240-246.

[62] H. E. Kandrup, *Stochastic inflation as a time dependent random walk*, Phys. Rev. **D**39 (1989) 2245.
[63] K.-i. Nakao, Y. Nambu and M. Sasaki, *Stochastic Dynamics of New Inflation*, Prog. Theor. Phys. **80** (1988) 1041.

[64] Y. Nambu, *Stochastic Dynamics of an Inflationary Model and Initial Distribution of Universes*, Prog. Theor. Phys. **81** (1989) 1037.

[65] S. Mollerach, S. Matarrese, A. Ortolan and F. Lucchin, *Stochastic inflation in a simple two field model*, Phys. Rev. D **44** (1991) 1670-1679.

[66] A. D. Linde, D. A. Linde and A. Mezhlumian, *From the Big Bang theory to the theory of a stationary universe*, Phys. Rev. D **49** (1994) 1783-1826, [gr-qc/9306035].

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

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

[69] C. Pattison, V. Vennin, H. Assadullahi and D. Wands, *Stochastic inflation beyond slow roll*, JCAP **1907** (2019) 031, [1905.06300].

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

[71] 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-78, [astro-ph/9507001].

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

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

[74] D. S. Salopek and J. R. Bond, *Nonlinear evolution of long wavelength metric fluctuations in inflationary models*, Phys. Rev. D **42** (1990) 3936-3962.

[75] 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].

[76] D. H. Lyth and D. Wands, *Conserved cosmological perturbations*, Phys. Rev. D **68** (2003) 103515, [astro-ph/0306498].

[77] G. I. Rigopoulos and E. P. S. Shellard, *The separate universe approach and the evolution of nonlinear superhorizon cosmological perturbations*, Phys. Rev. D **68** (2003) 123518, [astro-ph/0306620].

[78] D. H. Lyth and Y. Rodriguez, *The Inflationary prediction for primordial non-Gaussianity*, Phys. Rev. Lett. **95** (2005) 121302, [astro-ph/0504045].

[79] K. Enqvist, S. Nurmi, D. Podolsky and G. Rigopoulos, *On the divergences of inflationary superhorizon perturbations*, JCAP **0804** (2008) 025, [0802.0395].

[80] T. Fujita, M. Kawasaki, Y. Tada and T. Takesako, *A new algorithm for calculating the curvature perturbations in stochastic inflation*, JCAP **12** (2013) 036, [1308.4754].

[81] K. Ando and V. Vennin, *Power spectrum in stochastic inflation*, JCAP **04** (2021) 057, [2012.02031].
[82] Y. Tada and V. Vennin, *Statistics of coarse-grained cosmological fields in stochastic inflation*, JCAP 02 (2022) 021, [2111.15280].

[83] 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, [1912.05399].

[84] J. Martin and V. Vennin, *Stochastic Effects in Hybrid Inflation*, Phys. Rev. D 85 (2012) 043525, [1110.2070].

[85] M. Kawasaki and Y. Tada, *Can massive primordial black holes be produced in mild waterfall hybrid inflation?*, JCAP 1608 (2016) 041, [1512.03515].

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

[87] 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, [2012.06551].

[88] 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, [2111.07437].

[89] R. Mahbub and A. De, *Smooth coarse-graining and colored noise dynamics in stochastic inflation*, 2204.03859.

[90] T. Kloek and H. K. van Dijk, *Bayesian estimates of equation system parameters: An application of integration by monte carlo*, Econometrica 46 (1978) 1-19.

[91] 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-354, [nucl-th/9809075].

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

[93] L. Pinol, S. Renaux-Petel and Y. Tada, *A manifestly covariant theory of multifield stochastic inflation in phase space: solving the discretisation ambiguity in stochastic inflation*, JCAP 04 (2021) 048, [2008.07497].

[94] A. Araújo, *The Central Limit Theorem for Real and Banach Valued Random Variables*. Wiley, New York, 1980.

[95] S. Redner, *A Guide to First-Passage Processes*. Cambridge University Press, 2001, 10.1017/CBO9780511606014.

[96] H. Assadullahi, H. Firouzjahi, M. Noorbala, V. Vennin and D. Wands, *Multiple Fields in Stochastic Inflation*, JCAP 1606 (2016) 043, [1604.04502].

[97] V. Vennin, H. Assadullahi, H. Firouzjahi, M. Noorbala and D. Wands, *Critical Number of Fields in Stochastic Inflation*, Phys. Rev. Lett. 118 (2017) 031301, [1604.06017].

[98] H. Shen, L. D. Brown and H. Zhi, *Efficient estimation of log-normal means with application to pharmacokinetic data*, Statist. Med. 25 (2006) 3023-3038.

[99] R. D'Agostino and E. S. Pearson, *Tests for departure from normality. empirical results for the distributions of b₂ and √b₁*, Biometrika 60 (Dec., 1973) 613.
[100] X.-H. Zhou and S. Gao, *Confidence intervals for the log-normal mean*, *Statist. Med.* **16** (Apr., 1997) 783-790.

[101] D. Freedman, R. Pisani and R. Purves, *Statistics*. W.W. Norton, New York, 1998.

[102] J. Garcia-Bellido, A. D. Linde and D. Wands, *Density perturbations and black hole formation in hybrid inflation*, *Phys. Rev. D* **54** (1996) 6040-6058, [astro-ph/9605094].

[103] S. Clesse and J. García-Bellido, *Massive Primordial Black Holes from Hybrid Inflation as Dark Matter and the seeds of Galaxies*, *Phys. Rev. D92* (2015) 023524, [1501.07565].

[104] M. Noorbala and H. Firouzjahi, *Boundary crossing in stochastic inflation with a critical number of fields*, *Phys. Rev. D* **100** (2019) 083510, [1907.13149].

[105] S. Hooshangi, A. Talebian, M. H. Namjoo and H. Firouzjahi, *Multiple field ultraslow-roll inflation: Primordial black holes from straight bulk and distorted boundary*, *Phys. Rev. D* **105** (2022) 083525, [2201.07258].

[106] L. Perreault Levasseur, V. Vennin and R. Brandenberger, *Recursive Stochastic Effects in Valley Hybrid Inflation*, *Phys. Rev. D* **88** (2013) 083538, [1307.2575].

[107] T. Cohen, D. Green, A. Premkumar and A. Ridgway, *Stochastic Inflation at NNLO*, JHEP 09 (2021) 159, [2106.09728].

[108] M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias and A. Riotto, *The formation probability of primordial black holes*, *Phys. Lett. B* **820** (2021) 136602, [2105.07810].

[109] N. Kitajima, Y. Tada, S. Yokoyama and C.-M. Yoo, *Primordial black holes in peak theory with a non-Gaussian tail*, JCAP 10 (2021) 053, [2109.00791].

[110] “NIST Digital Library of Mathematical Functions.” http://dlmf.nist.gov/20.7E30, Release 1.0.14 of 2016-12-21.

[111] F. Edgeworth, *The law of error I*, Proc. Cambridge Philos. Soc. **20** (1905) 3665.

[112] D. L. Wallace, *Asymptotic approximations to distributions*, The Annals of Mathematical Statistics **29** (Sept., 1958) 635-654.