\section{Introduction}
\label{sec:intro}

In the early universe, the gravitational collapse of large density fluctuations may give rise to primordial black holes~\cite{Hawking:1971ei, Carr:1974nx, Carr:1975qj} (PBHs). Among the different mechanisms~\cite{Escriva:2022duf} that have been proposed to account for the presence of those large perturbations is the parametric amplification of quantum vacuum fluctuations during cosmic inflation~\cite{Starobinsky:1979ty, Mukhanov:1981xt, Starobinsky:1982ee, Guth:1982ec, Bardeen:1983qw}. This same mechanism is responsible for the cosmological perturbations observed in the cosmic-microwave-background anisotropies~\cite{Planck:2018vyg, Planck:2018nkj}, and in the large-scale structures of the universe~\cite{SDSS:2005xqv, BOSS:2014hwf, Amendola:2016saw, DES:2022qpf}. Those measurements confirmed that cosmological perturbations are almost scale-invariant, Gaussian and adiabatic, which is consistent with single-field models of inflation~\cite{Martin:2013tda, Planck:2018jri}.

If inflation is driven by a single scalar field, perturbations get large when the field velocity is small, which in general requires to cross very flat regions of the potential function. Typical examples involve inflection-points models~\cite{Choudhury:2013woa, Kawasaki:2016pql, Garcia-Bellido:2017mdw, Ezquiaga:2017fvi, Germani:2017bcs,  Ezquiaga:2018gbw, Bhaumik:2019tvl, Figueroa:2020jkf}, false-vacuum configurations~\cite{Dvali:2003vv, Hamada:2014wna, Bezrukov:2014bra, Kitajima:2019ibn, Geller:2022nkr, Gu:2022pbo, Animali:2022otk}, or step-like features~\cite{Inomata:2021tpx, Cai:2021zsp, Cai:2022erk}. Another possibility is to consider the case where the inflaton climbs up its potential: at some point, the field stops climbing and starts rolling down, hence at the turn-around point its velocity necessarily vanishes. When this happens, inflation must take place (since the kinetic energy of the inflaton vanishes), and we refer to this regime as ``uphill inflation''. 

Interestingly, the production of PBHs in such a scenario was studied early on in \Refa{Yokoyama:1998pt}, in double-well potentials similar to the one depicted in \Fig{fig:DWIpot}. If the width of the double well is of the order of the Planck mass, the velocity inherited by the inflaton from the first phase of inflation taking place at large-field values is such that, after overshooting the potential minimum, the inflaton climbs up the local maximum and its velocity vanishes close to the top, before it rolls down again towards one of the two minima. This triggers a second phase of inflation during which perturbations grow at a high rate and may thus later collapse into PBHs.

The configuration leading to maximal PBH production is the one where the inflaton comes to a rest exactly at the potential maximum. Classically, it would inflate there forever. However, quantum fluctuations are expected to carry the system away from such an unstable equilibrium configuration. This indicates that the backreaction of quantum fluctuations onto the background dynamics plays an important role in this model, and the goal of this paper is to carefully study this effect. 
If large fluctuations are produced, such a backreaction effect is known to be crucial in shaping their statistical properties~\cite{Pattison:2017mbe, Ezquiaga:2019ftu, Panagopoulos:2019ail, Figueroa:2020jkf, Pattison:2021oen, Achucarro:2021pdh, Hooshangi:2021ubn, Ezquiaga:2022qpw, Cohen:2022clv} so this is also why it must be properly accounted for. We do so by using the stochastic-$\delta N$ formalism~\cite{Starobinsky:1982ee, Starobinsky:1986fx, Fujita:2013cna, Vennin:2015hra}, in which small-scales quantum fluctuations shift the background dynamics when they are stretched to large distances, and effectively act as a random noise. 

We find that the abundance of PBHs takes a universal value that is always large and that does not depend on the energy at which uphill inflation occurs. Since it takes place over the last few $e$-folds~of inflation, black holes are produced with very small masses that Hawking evaporate before big-bang nucleosynthesis, but they are so abundant that they drive a phase of PBH domination in the early Universe. The amount of fine tuning required for this scenario to take place is substantially smaller than in alternative models, which makes it attractive. 

In passing, when employing the stochastic-$\delta N$ formalism we propose two technical improvements, which are developed in the context of the present model but that are nonetheless of broader applicability. The first one concerns the implementation of the boundary conditions in the first-passage-time problem, where we find an approximation scheme that proves impressively efficient when the end-of-inflation surface is located in the drift-dominated region of phase space. The second one is a computation of the convolution integrals derived in \Refa{Tada:2021zzj} for the one-point statistics of the coarse-grained curvature perturbations, by means of a sampling algorithm that is not more expensive than a mere first-passage-time calculation. This should be useful for future implementations of the stochastic-$\delta N$ approach in other contexts.

The rest of this paper is organised as follows. In \Sec{sec:double-well inflation}, we present the scenario of uphill inflation at the classical level, both for the background and for linear cosmological perturbations. In \Sec{sec:stochastic formalism} we show how quantum diffusion can be included in the model, and we derive the relevant stochastic equations. In \Sec{sec:Fokker Planck} we solve the first-passage time problem associated to these stochastic equations, which leads us to the abundance of PBHs from uphill inflation in \Sec{sec:BW and PBH}. We present our conclusions in \Sec{sec:Conclusions}, and defer to two appendices some of the technical details that are not essential to the understanding of the main arguments developed in the main text. 
\section{Uphill inflation}
\label{sec:double-well inflation}
We consider the situation where the inflaton climbs up a hilltop potential. As the field slows down, a phase of inflation may take place before the inflaton falls off the local maximum. Such a scenario is referred to as ``uphill inflation''.
\subsection{Double-well realisation}
\label{sec:double:well}
\begin{figure}
\centering
\includegraphics[width=0.50\textwidth]{figs/DWI_pot}
\includegraphics[width=0.48\textwidth]{figs/DWI_Ninf}
			\caption{Left panel: double-well potential~\eqref{eq:pot:dwi} for $\f_{\mathrm{0}}=0.8166M_\usssPl$. Right panel: Number of inflationary $e$-folds~realised between the first two crossings of a local minimum, as a function of $\f_{\mathrm{0}}$.} 
			\label{fig:DWIpot}
\end{figure} 
Although the results presented in this work do not depend on the detailed way uphill inflation is realised, for concreteness we first describe the case where it occurs in a double-well potential of the form
\begin{equation}
\label{eq:pot:dwi}
V(\phi) = M^4\left[\left(\frac{\phi}{M_\usssPl}\right)^2-\left(\frac{\f_{\mathrm{0}}}{M_\usssPl}\right)^2\right]^2\, ,
\end{equation}
which is displayed in the left panel of \Fig{fig:DWIpot}. In this model, inflation proceeds from large-field values, say at decreasing $\phi$ in the direction indicated by the arrows in \Fig{fig:DWIpot}. When the inflaton approaches its local minimum, inflation stops by violation of the slow-roll conditions. Details about this first phase of slow-roll inflation can be found for instance in \Refa{Martin:2013tda} (under the name ``double-well inflation''). What matters is that since slow roll is a dynamical attractor, it quickly erases any dependence on the initial field value and its velocity. The inflaton then crosses the local minimum and starts to climb up the hill. 

Whether or not it overshoots the local maximum depends only on the value of $\phi_0$. The reason is the following. The classical dynamics of the inflaton is driven by the Klein-Gordon and Friedmann equations
\begin{equation}\begin{aligned} 
\ddot{\phi}+3H\dot{\phi}+\frac{\mathrm{d} V}{\mathrm{d}\phi}=0
\quad\quad\text{and}\quad\quad
H^2=\frac{V+\frac{\dot{\phi}^2}{2}}{3M_\usssPl^2}\, ,
\end{aligned}\end{equation}
where $H=\dot{a}/a$ is the Hubble parameter and a dot denotes derivation with respect to cosmic time. In terms of the number of $e$-folds~$N=\ln a$, they can be written  as~\cite{Chowdhury:2019otk} 
\begin{equation}
\label{eq:KG}
\frac{2}{6M_\usssPl^2-\pi^2}\frac{\mathrm{d}\pi}{\mathrm{d} N}+\frac{\pi}{M_\usssPl^2}+\frac{\mathrm{d}\ln V}{\mathrm{d}\phi}=0
\quad\quad\text{and}\quad\quad
H^2=\frac{V}{3M_\usssPl^2-\frac{\pi^2}{2}}
\, ,
\end{equation}
where we have introduced $\pi\equiv\mathrm{d}\phi/\mathrm{d} N$. 
The Klein-Gordon equation for $\phi(N)$ does not involve $M^4$, and since initial conditions are washed out by the preceding dynamical attractor, the only relevant parameter of the (classical version of the) model is indeed $\f_{\mathrm{0}}$. Upon solving \Eq{eq:KG} numerically, we find that the inflaton falls on the other side of the hill if $\f_{\mathrm{0}}<\f_{\mathrm{0}}^\mathrm{cri}$, where
\begin{equation}
\label{eq:phi0:cri}
\f_{\mathrm{0}}^\mathrm{cri}\simeq 0.8166M_\usssPl\, .
\end{equation} 
Otherwise, it turns over and remains on the same side of the hill. 

When the inflaton climbs up its potential, its velocity decays and the first Hubble-flow parameter 
\begin{equation}
\label{eq:eps1:def}
\epsilon_1\equiv -\frac{\dot{H}}{H^2}=\frac{\pi^2}{2M_\usssPl^2}
\end{equation}
may become smaller than one (this is necessarily the case if $\f_{\mathrm{0}}>\f_{\mathrm{0}}^{\mathrm{cri}}$, since then the inflaton's velocity vanishes when it turns over, hence $\epsilon_1=0$ at that point). In that case, a second phase of inflation takes place, the duration of which is displayed in \Fig{fig:DWIpot}. One can check that when $\f_{\mathrm{0}}$ is fine-tuned to values close enough to $\f_{\mathrm{0}}^\mathrm{cri}$, the number of inflationary $e$-folds~in the uphill phase becomes large, and that it diverges at $\f_{\mathrm{0}}=\f_{\mathrm{0}}^\mathrm{cri}$. This is because, when $\f_{\mathrm{0}}=\f_{\mathrm{0}}^\mathrm{cri}$, the velocity inherited from the preceding phase of slow-roll inflation is such that it brings the inflaton to an exact rest at the location of the local maximum, where it inflates forever. 

As noted in \Refa{Yokoyama:1998pt}, this stationary configuration is nonetheless unstable, hence after some time it shall be broken by quantum fluctuations.\footnote{In \Refa{Yokoyama:1998pt}, a slightly different potential than \Eq{eq:pot:dwi} is considered, but it is still of the double-well type and the same considerations apply.} Since those fluctuations develop around a flat point in the potential, they acquire large cosmological amplitudes, which may result in the production of primordial black holes. To properly describe this mechanism, one has to account for the fact that fluctuations destabilise the background, and we will show below that this backreaction effect can be incorporated within the formalism of stochastic inflation. The goal of this work is to apply this formalism to uphill-inflation models, in order to investigate the production of large cosmological perturbations at the end of inflation in such scenarios.
\subsection{Classical phase space}
\label{sec:Classical:Phase:Space}
Let us start by studying the classical dynamics of uphill models. In practice, we focus on the region of the potential that is close to its local maximum, where $V(\phi)$ can be approximated by
\begin{equation}
\label{eq:pot:quad:max}
V(\phi) = V_0 - \frac{m^2}{2}\phi^2\, .
\end{equation}
In the double-well model~\eqref{eq:pot:dwi}, $V_0=M^4(\f_{\mathrm{0}}/M_\usssPl)^4$ and $m=2M^2\f_{\mathrm{0}}/M_\usssPl^2$ but we now consider \Eq{eq:pot:quad:max} as a generic parametrisation of any quadratic local maximum. We simply assume that the initial velocity of the inflaton is close to the one that would make it freeze at the top of the hill, whether this velocity is set by a preceding slow-roll attractor phase as in the model detailed in \Sec{sec:double:well} or by any other mechanism.

Close to the local maximum, the inflaton's velocity becomes negligible, $\pi\ll M_\usssPl$, and the potential is approximately constant, $V(\phi)\simeq V_0\simeq 3M_\usssPl^2H_0^2$. In this limit, \Eq{eq:KG} becomes linear, and it can be solved according to
\begin{equation}
\label{eq:KGsol:phi}
\frac{\phi(N)}{M_\usssPl} = x_\uine^{\left(\nu_0-\frac{3}{2}\right)N}+y_\uine^{\left(-\nu_0-\frac{3}{2}\right)N}
\end{equation}
where $x_\mathrm{in}$ and $y_\mathrm{in}$ are two integration constants and 
\begin{equation}
\label{eq:nu0:def}
\nu_0=\frac{3}{2}\sqrt{1+\frac{4m^2}{9H_0^2}}\, .
\end{equation} 
Since $\nu_0>3/2$, the first branch is unstable while the second branch decays away at late time. Therefore, the critical trajectory that  freezes on the local maximum at late time is such that $x_\mathrm{in}=0$ and $y_\mathrm{in}=\phi(N=0)/M_\usssPl$, and in the following we consider trajectories close enough to this configuration, \textsl{i.e.~} such that $x_\mathrm{in}\ll \phi(N=0)/M_\usssPl$. The parameter $x_\mathrm{in}$ thus quantifies the deviation from the critical trajectory. In the context of the double-well model~\eqref{eq:pot:dwi}, as will be made explicit in \Sec{sec:Conclusions}, this is related to the deviation of $\f_{\mathrm{0}}$ from $\f_{\mathrm{0}}^\mathrm{cri}$. 

Let us note that, in this regime, the three terms in the Klein-Gordon equation are comparable if the curvature of the potential at the end of the first phase of inflation is of order $m^2$. This is true in general, since for the first phase of inflation to end by violation of slow roll around a local minimum,  $m$ and $H$ need to be of the same order of magnitude. Assuming that $H$ does not decrease substantially between the two phases of inflation, $m$ is thus of order $H_0$, hence $\nu_0$ is of order one. In the double-well model described in \Sec{sec:double:well}, this is indeed the case since $m=2\sqrt{3}H_0M_\usssPl/\f_{\mathrm{0}}$, so $\nu_0\simeq 4.5$. Then, $\mathrm{d} V/\mathrm{d}\phi = -m^2\phi$, $3H\dot{\phi}=-3H_0^2(3/2+\nu_0)\phi$ and $\ddot{\phi}=H_0^2(3/2+\nu_0)^2\phi$ are indeed all of the same order. This means that uphill inflation proceeds neither in the slow roll regime ($\ddot{\phi}\ll 3H\dot{\phi},\mathrm{d} V/\mathrm{d}\phi$), nor in the ultra slow-roll regime ($\mathrm{d} V/\mathrm{d}\phi \ll \ddot{\phi}, 3H\dot{\phi}$). It constitutes a ``new'' regime of inflation where the field acceleration, the Hubble friction and the potential gradient all play equally important roles.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{figs/trajectory.pdf}
\caption{Phase-space portrait of uphill inflation in the $\phi,\mathrm{d}\phi/\mathrm{d} N$ plane. The $ x $ and $ y $ directions introduced in \Eq{eq:x:y:def} are displayed in red and yellow respectively. Inflation takes place when $\vert \mathrm{d} \phi/\mathrm{d} N\vert =\sqrt{2\epsilon_1}M_\usssPl<\sqrt{2}M_\usssPl$, \textsl{i.e.~} between the two dashed blue lines. The light-blue dashed rectangle corresponds to $m^2\phi^2/2< V_0/10$, and $\left(\mathrm{d} \phi/\mathrm{d} N\right)^2/2 < 3 M_\usssPl^2/10$, which stands for the region where the first slow-roll parameter is small and the potential is almost constant, and is where the approximation~\eqref{eq:KGsol:phi} applies. This region is named ``de Sitter'' as $V$ and $H$ are almost constant (even though $ \epsilon_2 $ and $ \epsilon_3 $ are not small). The green region is where $ \frac{\delta\nu}{\nu_0} < 0.1 $, see also \Fig{fig:nu-tilde-phase-space-color-map}. A few trajectories are shown, on which coloured dots indicate the point where $\frac{\delta\nu}{\nu_0}$ drops below $0.1$ (\textsl{i.e.~} where the green region is entered). Trajectories are computed in the full potential~\eqref{eq:pot:dwi} where $\f_{\mathrm{0}}$ is given by its critical value~\eqref{eq:phi0:cri} and $M/\Mp=5\times 10^{-4}$.}  
\label{fig:Phase-space}
\end{center}
\end{figure}
The two branches appearing in \Eq{eq:KGsol:phi} also suggest that a more convenient parametrisation of phase may be provided by the variables $x$ and $y$ such that $\phi/M_\usssPl=x+y$ and $\pi/M_\usssPl=(\nu_0-3/2)x-(\nu_0+3/2)y$. These relations can be readily inverted and give
\begin{equation}\begin{aligned} 
\label{eq:x:y:def}
x&=\dfrac{1}{4\nu_0M_\usssPl}\left[\left(3+2\nu_0\right)\phi+2\pi\right],\\
y&=-\dfrac{1}{4\nu_0M_\usssPl}\left[\left(3-2\nu_0\right)\phi+2\pi\right].\\
\end{aligned}\end{equation} 
By inserting \Eq{eq:KGsol:phi} into \Eq{eq:x:y:def}, one obtains $x(N)=x_\mathrm{in} e^{\left(\nu_0-3/2\right)N}$ and $y(N)=x_\mathrm{in} e^{\left(-\nu_0-3/2\right)N}$, \textsl{i.e.~} $x$ and $y$ are the growing and decaying directions respectively. 

A phase-space portrait in the plane $(\phi,\pi)$ is shown in \Fig{fig:Phase-space}, where the $x$ and $y$ directions are displayed in red and yellow respectively. Inflation takes place between the two dashed blue lines, which correspond to $\pi=\pm\sqrt{2}M_\usssPl$ [\textsl{i.e.~} $\epsilon_1=1$, see \Eq{eq:eps1:def}]. The solution~\eqref{eq:KGsol:phi} applies when the field velocity can be neglected in the Friedmann equation $H^2=V/(3M_\usssPl^2-\pi^2/2)$, and the potential is dominated by its constant term $V_0$. This region is delineated by the light-blue dashed rectangle, in which $\pi^2/2 < 3 M_\usssPl^2/10$ and $m^2\phi^2/2< V_0/10$. The $x$ and $y$ axes correspond to ``extremal'' trajectories: if initial conditions are set on the $x$ axis then the system remains on that axis and drifts away from the origin, while if initial conditions are set on the $y$ axis then the system is attracted towards the origin along that same axis. In this sense, the origin is an unstable saddle point.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{figs/trajectory-time.pdf}
\caption{Field velocity $\pi=\mathrm{d}\phi/\mathrm{d} N$ as a function of the number of $e$-folds~$N$ for the four trajectories displayed in \Fig{fig:Phase-space}. As in \Fig{fig:Phase-space}, the blue dashed lines stand for $\epsilon_1=\pm 1$, and the coloured dots indicate the point where each trajectory enters the domain of phase space where $ \delta\nu/\nu_0 < 0.1 $. The origin $N=0$ for the number of $e$-folds~is set to the time when inflation resumes. }  
\label{fig:pi-trajectory}
\end{center}
\end{figure}
Four trajectories are represented in \Fig{fig:Phase-space}, one in each $(x,y)$ quarter plane. The green trajectory ($x>0,y>0$) climbs up the potential from positive field values, turns around before reaching the maximum and returns to the rightmost local minimum. Similarly, the light blue trajectory ($x<0,y<0$) climbs up the potential from negative field values, turns around and returns to the leftmost local minimum. The brown (respectively grey) trajectory crosses over the local maximum leftwards (respectively rightwards). These four trajectories are also displayed in \Fig{fig:pi-trajectory}, where $\pi$ is shown as a function of $N$. 
\subsection{Linear cosmological perturbations}
\label{sec:CPT}
Let us now turn our attention to cosmological perturbations. This is required to compute the amplitude of the noise in the stochastic-inflation formalism. The scalar sector can be described by a single gauge-invariant combination, such as the Mukhanov-Sasaki variable $v(\vec{x},t)$~\cite{Mukhanov:1981xt,Kodama:1984ziu}. At linear order, its Fourier modes $v_k$ obey the  
Mukhanov-Sasaki equation
\begin{equation}
\label{eq:MS:eq}
v_k''+\left(k^2-\frac{Z''}{Z}\right) v_k=0\, ,
\end{equation} 
where a prime denotes derivation with respect to conformal time $\eta$, and $Z\equiv aM_\usssPl\sqrt{2\epsilon_1}$. In terms of the Hubble-flow parameters, defined by $\epsilon_{n+1}= \mathrm{d}\ln\epsilon_n/\mathrm{d} N$, where $\epsilon_1$ was introduced in \Eq{eq:eps1:def}, one has
\begin{equation}
\frac{Z''}{Z} = \left(a H\right)^2\left(2-\epsilon_1+\frac{3}{2}\epsilon_2-\frac{1}{2}\epsilon_1\epsilon_2 + \frac{1}{4}\epsilon_2^2+\frac{1}{2}\epsilon_2\epsilon_3\right).
\end{equation}
Close to the phase-space origin, $H\simeq H_0$ hence $aH\simeq -1/\eta$, and $Z''/Z=(\nu^2-1/4)/\eta^2$ with $ \nu = \frac{3}{2}\sqrt{1+\frac{4}{9}(-\epsilon_1+\frac{3}{2}\epsilon_2-\frac{1}{2}\epsilon_1\epsilon_2 + \frac{1}{4}\epsilon_2^2+\frac{1}{2}\epsilon_2\epsilon_3)} $. The Hubble-flow parameters can be expressed in terms of $\phi$ and $\pi$ (or equivalently in terms of $x$ and $y$), hence they are phase-space functions. Indeed, from differentiating \Eq{eq:eps1:def} with respect to time and using the Klein-Gordon equation~\eqref{eq:KG}, one finds
\begin{equation}\begin{aligned}
\epsilon_2 &= \left(1+\frac{M_\usssPl^2}{\pi}\frac{\mathrm{d}\ln V}{\mathrm{d}\phi}\right)\left(\frac{\pi^2}{M_\usssPl^2}-6\right),\\
\epsilon_3 &= M_\usssPl^2\frac{\mathrm{d}^2\ln V}{\mathrm{d}\phi^2}- \frac{M_\usssPl^4  \frac{\mathrm{d}^2\ln V}{\mathrm{d}\phi^2}}{M_\usssPl^2+\frac{\pi}{\frac{\mathrm{d}\ln V}{\mathrm{d} \phi}}}
+3\frac{M_\usssPl^2}{\pi}\frac{\mathrm{d}\ln V}{\mathrm{d} \phi} + \frac{\pi}{2}\frac{\mathrm{d}\ln V}{\mathrm{d}\phi} +\frac{\pi^2}{M_\usssPl^2}\, .
\end{aligned}\end{equation}
It is interesting to notice that these functions are not well defined at the origin of phase space $(\phi=0,\pi=0)$, since they involve the ratio $(\mathrm{d}\ln V/\mathrm{d}\phi)/\pi$, the value of which depends on the direction along which the origin is approached. However, remarkably, this singular contribution cancels out in the combination giving $\nu$, and one finds
\begin{equation}\begin{aligned}
\label{eq:nu:phi:dphi_dN}
\nu=\frac{3}{2} &\left\lbrace
1+\frac{2}{9M_\usssPl^4}\left[ \pi ^4+2M_\usssPl^2 \frac{\mathrm{d}\ln V}{\mathrm{d}\phi} \pi ^3+M_\usssPl^2  \pi ^2\left(M_\usssPl^2\frac{\frac{\mathrm{d}^2 V}{\mathrm{d}\phi^2}}{V}-7\right)
\right.\right. \\ & \left. \left.
-12M_\usssPl^4  \pi \frac{\mathrm{d}\ln V}{\mathrm{d}\phi}-6M_\usssPl^6\frac{\frac{\mathrm{d}^2 V}{\mathrm{d}\phi^2}}{V} \right]
\right\rbrace^{1/2} .
\end{aligned}\end{equation}
This expression is perfectly regular at the phase-space origin where it reduces to $\nu_0$, which was introduced in \Eq{eq:nu0:def}. The Mukhanov-Sasaki equation is therefore well-behaved everywhere in phase space, even though the higher Hubble-flow parameters are ill-defined in the configuration where the field rests with vanishing velocity at a potential's local extremum.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{figs/colormap.pdf}
\caption{Relative difference between the index $\nu$ and its value at the phase-space origin, $ \delta\nu/\nu_0=(\nu-\nu_0)/\nu_0 $, as a function of $\phi$ and $\pi$. The boundaries of the axes match the dashed-blue rectangle in \Fig{fig:Phase-space}, inside of which $V$ and $H$ are almost constant.}  
\label{fig:nu-tilde-phase-space-color-map}
\end{center}
\end{figure}

The function $\nu(\phi,\pi)$ given in \Eq{eq:nu:phi:dphi_dN} is shown in \Fig{fig:nu-tilde-phase-space-color-map}, where the relative difference $\delta\nu/\nu_0 = (\nu-\nu_0)/\nu_0$ is displayed. The frame of the figure corresponds to the dashed-blue rectangle in \Fig{fig:Phase-space}, within which $H\simeq H_0$ and $V\simeq V_0$ and our approximation applies. One can see that $\nu$ does not vary much within this domain, not more than $\sim 10\%$ except in the upper-left and bottom-right corners. In \Fig{fig:Phase-space}, the region where $\delta\nu/\nu_0<0.1$ is displayed in light green. One can check that as soon as a given trajectory enters that region (this moment is denoted by the coloured dots in \Figs{fig:Phase-space} and~\ref{fig:pi-trajectory}), it stays there until inflation ends. Moreover, once the system is attracted close to the $x$ axis (where we will see that most relevant stochastic effects take place), $\delta\nu/\nu_0$ is subject to tiny variations only, see \Fig{fig:nu-tilde-phase-space-color-map}. 

For these reasons, it is a good approximation to assume that $\nu\simeq \nu_0$ in the phase-space region of interest, and we have checked the validity of this approximation numerically by comparing the solution to \Eq{eq:MS:eq} with and without replacing $\nu=\nu_0$. 
One may be worried that, in practice, the system still inflates after exiting the domain where $H\simeq H_0, V\simeq V_0$ (the dashed blue rectangle in \Fig{fig:Phase-space}). However, this last stage of inflation is found to be very short-lived and subject to negligible stochastic effects, since it corresponds to the phase where the field is classically drifted away from the potential maximum. The loss of accuracy in the calculation of the noise amplitude in that phase has therefore negligible impact.

If $\nu$ is constant and equal to $\nu_0$, the Mukhanov-Sasaki equation~\eqref{eq:MS:eq} can be solved analytically in terms of the Hankel function of the first kind $H^{(1)}_{\nu_0}$, 
\begin{equation}
\label{eq:vk:Hankel}
v_k(\eta) = \frac{\sqrt{-\pi\eta}}{2} e^{i\frac{\pi}{2}\left(\nu_0+\frac{1}{2}\right)} H_{\nu_0}^{(1)}(-k\eta) \,,
\end{equation} 
where the integration constants have been set such as the Bunch-Davies vacuum $v_k(\eta)=e^{-ik\eta}/\sqrt{2k}$ is recovered in the sub-Hubble, asymptotic past. In the super-Hubble regime, $-k\eta\ll 1$, \Eq{eq:vk:Hankel} reduces to
\begin{equation}
\label{eq:v:superH}
v_k(\eta)\simeq -i e^{i\frac{\pi}{2}\left(\nu_0+\frac{1}{2}\right)}\frac{\Gamma(\nu_0)}{\sqrt{2\pi k}} \left(\frac{-k\eta}{2}\right)^{\frac{1}{2}-\nu_0}\, .
\end{equation} 
Let us note that, along the decaying branch $y$ [\textsl{i.e.~} when setting $x_\mathrm{in}=0$ in \Eq{eq:KGsol:phi}], the curvature perturbations, $\zeta=v/z$, quickly grows on super-Hubble scales. This is because $z\propto \eta^{1/2+\nu_0}$ in that case, hence \Eq{eq:v:superH} implies that $\zeta\propto \eta^{-2\nu_0}$.\footnote{On the contrary, on the growing branch $x$, one has $z\propto \eta^{1/2-\nu_0}$, hence $\zeta$ is frozen. This is expected as the $x$-branch is a dynamical attractor.} With $\nu_0\simeq 4.5$ (corresponding to the double-well model described in \Sec{sec:double:well}), this leads to $\zeta\sim \eta^{-9}$, so the growth rate is indeed very large (recall that $\zeta$ is conserved in the presence of a phase-space attractor such as slow roll, and that it grows as $\zeta\sim \eta^{-3}$ in ultra slow-roll inflation). This suggests that large cosmological fluctuations develop in uphill inflation, that could be responsible for the formation of primordial black holes. These considerations however ignore the crucial role quantum diffusion plays in this model, which is the topic of the next section.


\section{Quantum diffusion}
\label{sec:stochastic formalism}

If the inflaton comes to a rest at a local maximum of its potential, it will soon be destabilised by quantum fluctuations, which implies that the classical description given in \Sec{sec:Classical:Phase:Space} breaks down. This can be modelled within the formalism of stochastic inflation, which we now briefly review. 
\subsection{Stochastic inflation}
Stochastic inflation~\cite{Starobinsky:1982ee, Starobinsky:1986fx} is an effective description for the long-wavelength part of quantum fields living on (and possibly sourcing) an inflating cosmological background. In practice, the fields are coarse-grained at the physical scale $R=(\sigma H)^{-1}$, where $\sigma\ll 1$ is the ratio between the Hubble radius and the coarse-graining radius, according to
\begin{equation}\begin{aligned}
\hat{\phi}_{\mathrm{IR}}(\vec{x},N)&=\int\frac{\mathrm{d}\vec{k}}{(2\pi)^{3/2}} W\left(\frac{k}{\sigma a H}\right)\left[ e^{-i \vec{k}\cdot\vec{x}}\phi_{k}(N) \hat{a}_{\vec{k}}+e^{i \vec{k}\cdot\vec{x}}\phi_{k}^*(N) \hat{a}^\dagger_{\vec{k}}\right]\, ,\\
\hat{\pi}_{\mathrm{IR}}(\vec{x},N)&=\int\frac{\mathrm{d}\vec{k}}{(2\pi)^{3/2}} W\left(\frac{k}{\sigma a H}\right) \left[ e^{-i \vec{k}\cdot\vec{x}}\pi_{k}(N) \hat{a}_{\vec{k}}+e^{i \vec{k}\cdot\vec{x}}\pi_{k}^*(N) \hat{a}^\dagger_{\vec{k}}\right]\, .
\end{aligned}\end{equation} 
In this expression, $\hat{a}_{\vec{k}}$ and $\hat{a}_{\vec{k}}^\dagger$ are annihilation and creation operators, hats indicate that we are dealing with quantum operators, and $W$ is a window function that selects modes $k$ with wavelength larger than the coarse-graining scale, \textsl{i.e.~} $W\simeq 1$ for $k\ll \sigma a H$ and $0$ for $k\gg \sigma a H$. On super-Hubble scales, gradient terms can be neglected and the mode functions $\phi_k$ and $\pi_k$ follow the same equations of motion as the classical background. This is the so-called separate-universe approximation~\cite{Salopek:1990jq, Sasaki:1995aw, Wands:2000dp, Pattison:2019hef, Artigas:2021zdk}. The dynamics of $\hat{\phi}_{\mathrm{IR}}$ and $\hat{\pi}_{\mathrm{IR}}$ is thus given by~\cite{Grain:2017dqa}
\begin{equation}\begin{aligned} 
\label{eq:Langevin:interm}
\frac{\mathrm{d}\hat{\phi}_{\mathrm{IR}}}{\mathrm{d} N}&=\hat{\pi}_{\mathrm{IR}}+\hat{\xi}_\phi(N)\, ,\\
\frac{\mathrm{d}\hat{\pi}_{\mathrm{IR}}}{\mathrm{d} N}&=-\left(3-\frac{\hat{\pi}_{\mathrm{IR}}^2}{2M_\usssPl^2}\right)\hat{\pi}_{\mathrm{IR}}- \frac{\frac{\mathrm{d} V}{\mathrm{d}\phi}\left(\hat{\phi}_{\mathrm{IR}}\right)}{H^2\left(\hat{\phi}_{\mathrm{IR}},\hat{\pi}_{\mathrm{IR}}\right)}+\hat{\xi}_\pi(N)\, ,
\end{aligned}\end{equation} 
where the source functions $\hat{\xi}_\phi$ and $\hat{\xi}_\pi$ read
\begin{equation}\begin{aligned} 
\hat{\xi}_\phi(N)&=-\int\frac{\mathrm{d}\vec{k}}{(2\pi)^{3/2}} \left[ e^{-i \vec{k}\cdot\vec{x}}\phi_{k}(N) \hat{a}_{\vec{k}}+e^{i \vec{k}\cdot\vec{x}}\phi_{k}^*(N) \hat{a}^\dagger_{\vec{k}}\right]\frac{\mathrm{d} }{\mathrm{d} N}W\left(\frac{k}{\sigma a H }\right)\, ,\\
\hat{\xi}_\pi(N)&=-\int\frac{\mathrm{d}\vec{k}}{(2\pi)^{3/2}} \left[ e^{-i \vec{k}\cdot\vec{x}}\pi_{k}(N) \hat{a}_{\vec{k}}+e^{i \vec{k}\cdot\vec{x}}\pi_{k}^*(N) \hat{a}^\dagger_{\vec{k}}\right]\frac{\mathrm{d} }{\mathrm{d} N}W\left(\frac{k}{\sigma a H }\right)\, ,
\end{aligned}\end{equation} 
and where the mode functions are to be evaluated in the uniform-$N$ gauge~\cite{Pattison:2019hef}.
The stochastic formalism then consists in treating \Eqs{eq:Langevin:interm} as stochastic, Langevin equations, where $\hat{\xi}_\phi$ and $\hat{\xi}_\pi$ are replaced with stochastic noises, the statistics of which are drawn from their quantum expectation values. The noises $\hat{\xi}_\phi$ and $\hat{\xi}_\pi$ involve small wavelength modes, which are described by standard cosmological perturbation theory. At linear order, they are placed in a Gaussian state, hence the noises are centred Gaussian noises. If $W$ is set to a Heaviside step function, \textsl{i.e.~} $W=1$ if $k<\sigma a H$ and $0$ otherwise, then their covariance is given by
\begin{equation} 
\label{eq:cov:interm}
\left\langle \hat{\xi}_f (N) \hat{\xi}_g^\dagger(N') \right\rangle = \frac{1}{6\pi^2} \frac{\mathrm{d} \left(\sigma a H\right)^3}{\mathrm{d} N} 
\Re \mathrm{e}\,\left[f_{\sigma a H}(N) g^*_{\sigma a H}(N)\right]\delta(N-N')\, .
\end{equation} 
where $f,g=\phi,\pi$. 

As mentioned above, $\phi_k$ and $\pi_k$ correspond to the mode functions of the field fluctuations in the uniform-$N$ gauge [since perturbations to the lapse function have been discarded in \Eq{eq:Langevin:interm}]. In \Sec{sec:CPT}, we have explained how to compute the mode function of the Mukhanov-Sasaki variable, so these quantities now need to be related. The relevant relationship is obtained in \Refa{Pattison:2019hef}, and here we simply recall the result. 

The Mukhanov-Sasaki variable can be identified with the scalar field fluctuation in the spatially-flat gauge, $v_k/a=\phi_k^\mathrm{flat}$. A gauge transformation is then needed to compute $\phi_k$ in the uniform-$N$ gauge, which reads~\cite{Pattison:2019hef}
\begin{equation}
\label{eq:gauge:transform}
\phi_k = \phi_k^\mathrm{flat}+ \phi' \alpha_k\, ,
\end{equation}
where $\alpha$ is solution to the differential equation
\begin{equation}
\label{equa alpha}
3\mathcal{H}\alpha_k' + (3\mathcal{H}' -k^2)\alpha_k = S_k,
\end{equation}
where $\mathcal{H}=aH=a'/a$ is the conformal Hubble parameter, and the source $S_k$ is given by 
\begin{equation}
\label{source}
S_k = \frac{v_k\sqrt{2\epsilon_1}}{2 a\Mp}\mathrm{sign}(\dot\phi)\left[\frac{\mathcal{H}\epsilon_2}{2}-\frac{(v_k/a)'}{v_k/a}\right]\, .
\end{equation}
Noting that \Eq{equa alpha} can be solved as
\begin{equation}
\label{formule alpha}
\alpha_k = \frac{1}{3\mathcal{H}}\int_{\eta_0}^{\eta}S_k(\eta')\exp\left[\frac{k^2}{3}\int_{\eta'}^{\eta}\frac{\mathrm{d}\eta''}{\mathcal{H}(\eta'')}\right]\mathrm{d}\eta'\, ,
\end{equation}
the gauge transformation~\eqref{eq:gauge:transform} can be performed explicitly. In \Eq{formule alpha}, $\eta_0$ is an integration constant that defines the slicing relative to which the expansion is measured, which here we take in the asymptotic future ($\eta_0=0^-$).

In uphill inflation, \Eq{eq:v:superH} indicates that $v_k\propto \eta^{1/2-\nu_0}$ on super-Hubble scales, hence $(v_k/a)'/(v_k/a)=(2\nu_0-3)\mathcal{H}/2$.
Along the decaying branch $y$, one has $\epsilon_1\propto \eta^{3+2\nu_0}$, hence $\epsilon_2=-2\nu_0-3$. This leads to $S_k\propto \eta^2$ in \Eq{source}, hence $\alpha_k \propto \eta^4$ in \Eq{formule alpha}. The gauge correction in \Eq{eq:gauge:transform}, $\phi ' \alpha_k \propto \eta^{\nu_0+9/2}$, is therefore negligible compared to $\phi_k^{\mathrm{flat}}=v_k/a\propto \eta^{3/2-v_0}$ on super-Hubble scales. Along the growing branch $x$,   one has $\epsilon_1\propto \eta^{3-2\nu_0}$, hence $\epsilon_2=2\nu_0+3$. This leads to an exact cancellation in the terms $\mathcal{H} \epsilon_2/2-(v_k/a)'/(v_k/a)$ appearing in the source function given in \Eq{source}. It implies that the gauge correction is even more suppressed in this case, which is expected since the growing branch acts as a dynamical attractor~\cite{Pattison:2019hef}.




We thus conclude that the gauge corrections are negligible on super-Hubble scales, where $\phi_k\simeq v_k/a \simeq 2^{\nu_0-1}H_0\Gamma(\nu_0)k^{-3/2}(-k\eta)^{3/2-\nu_0}/\sqrt{\pi }$ (up to an irrelevant global phase) and $\pi_k \simeq (\nu_0-\frac{3}{2})\phi_k$. By substituting these expressions into \Eq{eq:cov:interm}, the covariance matrix of the noises is given by
\begin{equation}\begin{aligned}
\left\langle 
\begin{pmatrix}
\hat{\xi}_\phi(N)\\
\hat{\xi}_\pi(N)
\end{pmatrix}
\begin{pmatrix}
\hat{\xi}_\phi^\dagger(N') & \hat{\xi}_\pi^\dagger(N')
\end{pmatrix}\right\rangle=2\frac{M_\usssPl^2}{\mu^2}
\begin{pmatrix} 1 &\ & \nu_0-\frac{3}{2}  \\ &\ & \\ \nu_0-\frac{3}{2} &\ & (\nu_0-\frac{3}{2})^2 \end{pmatrix}\delta\left(N-N'\right)\, ,
\end{aligned}\end{equation} 
where 
\begin{equation} 
\label{eq:mu:def}
\mu \equiv \frac{\sqrt{2}\pi^{3/2}}{\Gamma(\nu_0)}\frac{M_\usssPl}{H_0}\left(\frac{\sigma}{2}\right)^{\nu_0-\frac{3}{2}}\, .
\end{equation} 
Since the determinant of the covariance matrix vanishes, there is a single independent noise, and one can write $\xi_\phi=\sqrt{2}M_\usssPl\xi/\mu$ and $\xi_\pi=(\nu_0-3/2)\sqrt{2}M_\usssPl\xi/\mu$, where $\xi$ is a normalised white Gaussian noise, \textsl{i.e.~} such that $\langle \xi(N) \xi(N')\rangle = \delta(N-N')$. At the level of the approximation $\pi\ll M_\usssPl$ performed here, the term proportional to $\hat{\pi}^2$ in \Eq{eq:Langevin:interm} must be neglected, and one obtains 
\begin{equation}\begin{aligned}
\label{Langevin}
    \frac{\mathrm{d}\phi}{\mathrm{d} N} &= \pi + \sqrt{2}\frac{M_\usssPl}{\mu} \xi (N) \, ,\\
    \frac{\mathrm{d}\pi}{\mathrm{d} N} &= -3\pi + \frac{m^2}{H_0^2}\phi + \left(\nu_0-\frac{3}{2}\right)\sqrt{2}\frac{M_\usssPl}{\mu} \xi (N)\, .
\end{aligned}\end{equation}
In this expression, the subscripts ``$\mathrm{IR}$'' have been dropped for notational convenience, and the hats have been removed too since one now deals with classical, random variables. These are the Langevin equations we aim at solving in the rest of this paper\footnote{Formally, when $m^2\to 0$ while keeping $H_0$ fixed, $\nu_0\to 3/2$, $\gamma\to H_0/(2\pi)$ and one recovers the ultra-slow roll Langevin equations~\cite{Pattison:2021oen}. This consistency check is nonetheless purely formal as ultra-slow roll is not a limiting case of our model (where $m$ and $H_0$ are of the same order).}.

These Langevin equations can also be rewritten in the $x$ and $y$ coordinates introduced in \Eq{eq:x:y:def}, and one obtains
\begin{equation}\begin{aligned}
\label{Langevin rotated}
    \frac{\mathrm{d}x}{\mathrm{d} N} &= \left(\nu_0-\frac{3}{2}\right)x + \frac{\sqrt{2}}{\mu}\xi (N) \, ,\\
    \frac{\mathrm{d}y}{\mathrm{d} N} &= -\left(\nu_0+\frac{3}{2}\right)y\, .
\end{aligned}\end{equation}
At the classical level (\textsl{i.e.~} without the stochastic noise $\xi$), these equations decouple, which is the reason why the variables $x$ and $y$ were introduced. We now find that these variables are also convenient at the stochastic level, since quantum diffusion only takes place along the $x$ direction. As a consequence, even if initial conditions are set along the $y$ direction (\textsl{i.e.~} such that the field classically freezes at the potential maximum), quantum diffusion explores the $x$ direction, and necessarily destabilises the system. We thus expect that, at the stochastic level, the number of inflationary $e$-folds~remains finite even when $x_\mathrm{in}=0$, contrary to the classical setting (see the right panel of \Fig{fig:DWIpot}).
\subsection{The stochastic-$\delta N$ formalism}
\label{sec:StochaDeltaN}
Having a formalism at our disposal to describe the backreaction of quantum fluctuations onto the inflationary dynamics, the next step is to investigate how quantum backreaction affects observables, \textsl{i.e.~} the statistics of cosmological fluctuations. In standard cosmological perturbation theory, cosmological fluctuations are placed in a Gaussian state, the two-point functions of which are simply given by squaring the mode functions computed in \Sec{sec:CPT}. In stochastic inflation, the statistics of the curvature perturbation can be obtained using the so-called stochastic-$\delta N$ formalism~\cite{Enqvist:2008kt, Fujita:2013cna, Vennin:2015hra}. It relies on the fact that on super-Hubble scales, the local fluctuation in the amount of expansion $\mathcal{N}$ between an initial spatially-flat hypersurface and a final hypersurface of uniform energy density, is nothing but the curvature perturbation~\cite{Starobinsky:1982ee, Starobinsky:1986fxa, Sasaki:1995aw, Wands:2000dp, Lyth:2004gb}, $\zeta(\vec{x})=\mathcal{N}(\vec{x}) - \langle \mathcal{N} \rangle_{\vec{x}} =\delta\mathcal{N}(\vec{x})$. 

More precisely, when the curvature perturbation is coarse-grained at a scale $R$, the one-point statistics of $\zeta_R$ is given by the distribution function~\cite{Tada:2021zzj}
\begin{equation}
\label{eq:zeta:one:point:pdf}
P\left(\zeta_R\right) = \int \mathrm{d} x_*\mathrm{d} y_*P_{\mathrm{BW}}\left[x_*,y_*\lvert N_{\mathrm{BW}}(R)\right]{P_{\mathrm{FPT}}}_{x_\mathrm{in},y_\mathrm{in} \rightarrow x_*,y_*}\left[\zeta_R-\langle\mathcal{{N}}\rangle\left(x_*,y_*\right)+\langle\mathcal{{N}}\rangle\left(x_\mathrm{in},y_\mathrm{in}\right)\right] .
\end{equation}
In this expression, $N_{\mathrm{BW}}(R)$ corresponds to the number of $e$-folds~elapsed between the time when the scale $R$ exits the coarse-graining scale and the end of inflation. If $H$ is almost constant (which we assume here), it is given by $N_{\mathrm{BW}}(R)\simeq \ln(\sigma H R) $. The value of $x$ and $y$, $N_{\mathrm{BW}}(R)$ $e$-folds~before the end of inflation, are denoted $x_*$ and $y_*$ and they follow the ``backward'' distribution function~\cite{Ando:2020fjm}
\begin{equation}\begin{aligned}
\label{eq:PBW:def}
P_{\mathrm{BW}}\left[x_*,y_*\lvert N_{\mathrm{BW}}(R)\right]= P_{\mathrm{FPT}}\left[N_{\mathrm{BW}}(R)\vert x_*,y_*\right]\dfrac{\int_0^\infty \mathrm{d} N' P\left(x_*,y_* \vert N',x_{\mathrm{in}},y_{\mathrm{in}}\right)}{\int_{N_{\mathrm{BW}}(R)}^\infty \mathrm{d}\mathcal{N} P_{\mathrm{FPT}}\left(\mathcal{N}\vert x_\mathrm{in},y_\mathrm{in}\right)}\, .
\end{aligned}\end{equation} 
Here, $P_{\mathrm{FPT}}(N\vert x,y)$ is the probability that, starting from $x$ and $y$, inflation proceeds for $N$ $e$-folds~before ending. This is the so-called ``first-passage time'' (FPT) probability. Conversely, $P(x,y\vert N,x_{\mathrm{in}},y_{\mathrm{in}})$ is the probability that, starting from $x_\mathrm{in}$ and $y_\mathrm{in}$, after $N$ $e$-folds, the system is at field values $x$ and $y$. Finally, in \Eq{eq:zeta:one:point:pdf}, $\langle \mathcal{N}\rangle(x,y)$ denotes the mean first-passage time starting from $x$ and $y$, and ${P_{\mathrm{FPT}}}_{x_\mathrm{in},y_\mathrm{in} \rightarrow x_*,y_*}$ is the distribution for the time of first passage through $(x_*,y_*)$ (assuming this location is indeed visited at least once), starting from $(x_\mathrm{in},y_\mathrm{in})$. 

Similar expressions can be obtained for the two-point statistics, \textsl{i.e.~} the power spectrum (see \Refa{Ando:2020fjm}), and for the one-point statistics of the density contrast and of the compaction function (see \Refa{Tada:2021zzj}). Although seemingly not straightforward, they require to determine only two functions, $P$ and $P_{\mathrm{FPT}}$, and in the next section we show how this can be done. 

When reconstructed numerically from Langevin simulations, $P(\zeta_R)$ is even simpler to obtain, along the following steps:
\begin{enumerate}
    \item \label{item:1} Pre-compute $\langle \mathcal{N} \rangle (x_*,y_*)$ on a grid of points $(x_*,y_*)$ to interpolate that function.
    \item \label{item:2} From the initial condition $(x_{\mathrm{in}},y_{\mathrm{in}})$, simulate one realisation of the Langevin equations and record the number of inflationary $e$-folds~$\mathcal{N}$ as well as the field values $(x_*,y_*)$ at a time $N_{\mathrm{BW}}$ before inflation ends.
    \item \label{item:3}Evaluate $\zeta_R = \mathcal{N}-N_{\mathrm{BW}}+\langle \mathcal{N}\rangle(x_*,y_*)-\langle \mathcal{N}\rangle(x_{\mathrm{in}},y_{\mathrm{in}})$ and store that value.
    \item Repeat steps \ref{item:2}-\ref{item:3} a large number of times.
    \item reconstruct $P(\zeta_R)$ from the obtained list of values for $\zeta_R$.
\end{enumerate}
This procedure allows one to sample $(x_*,y_*)$ according to the backward probability, hence the first-passage-time distribution is indeed weighted by that probability as required by \Eq{eq:zeta:one:point:pdf}. Note that, strictly speaking, the above steps sample ${P_{\mathrm{FPT}}}_{x_\mathrm{in},y_\mathrm{in} \rightarrow x_*,y_*}$ under the condition that between $(x_*,y_*)$ and the end-of-inflation surface, $N_{\mathrm{BW}}$ $e$-folds~are realised. However, the stochastic process being Markovian, the number of $e$-folds~realised before and after $(x_*,y_*)$ are independent, hence this additional condition is inoffensive. Apart from the pre-computation in step \ref{item:1}, this procedure is not computationally more expensive than the reconstruction of a mere first-passage time distribution.

\section{Field and first-passage-time probabilities}
\label{sec:Fokker Planck}

In this section, we derive the distribution functions $P$ and $P_{\mathrm{FPT}}$, which respectively denote the probability associated to the field values after a given time, and the probability associated to the duration of inflation starting from given field values. It is first convenient to recast the Langevin equations~\eqref{Langevin rotated} in terms of the rescaled variables $ z = \mu x \sqrt{(\nu_0-3/2)/2} $ and $ \widetilde{N} = (\nu_0-3/2 )N $,
\begin{align}
\label{Langevin rotated new variables:z}
    \frac{\mathrm{d}z}{\mathrm{d} \widetilde{N}} &= z + \xi (\widetilde{N})\, ,\\
    \frac{\mathrm{d}y}{\mathrm{d} \widetilde{N}} &= -\frac{\nu_0+\frac{3}{2}}{\nu_0-\frac{3}{2}}\,y\, .
\label{Langevin rotated new variables:y}
\end{align}
The Fokker-Planck equation driving $P(z,y\vert\widetilde{N},z_{\mathrm{in}},y_{\mathrm{in}})$ can then be derived~\cite{risken1989fpe},
\begin{equation}
\label{Fokker-Planck PDF}
\frac{\partial P}{\partial \widetilde{N}} = \mathcal{L}_{\mathrm{FP}}P,
\end{equation}
where $\mathcal{L}_{\mathrm{FP}}$ is the so-called Fokker-Planck operator and in the case of \Eqs{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y} it is given by
\begin{equation}\begin{aligned}
\label{eq:LFP:def}
\mathcal{L}_{\mathrm{FP}} &= \frac{1}{2}\frac{\partial^2}{\partial z^2} + \frac{\nu_0+\frac{3}{2}}{\nu_0-\frac{3}{2}}\, \frac{\partial}{\partial y}y - \frac{\partial}{\partial z}z \, .
\end{aligned}\end{equation}
Since the noise vanishes in the $ y $-direction, $y$ follows a purely classical trajectory and $\mathcal{L}_{\mathrm{FP}}$ does not contain second derivatives with respect to $y$.
The Fokker-Planck equation~\eqref{Fokker-Planck PDF} needs to be solved with the initial condition $P(z,y\vert 0,z_{\mathrm{in}},y_{\mathrm{in}})= \delta(z-z_{\mathrm{in}})\delta(y-y_{\mathrm{in}}) $, enforcing the fact that at initial time, the system starts at $z_{\mathrm{in}}$ and $y_{\mathrm{in}}$. Moreover, an absorbing boundary needs to be placed where inflation ends, \textsl{i.e.~} at $\pi=\pm \sqrt{2}M_\usssPl$. In terms of the rescaled coordinates $z$ and $y$, this condition reads $2\sqrt{2\nu_0-3}z/\mu=(2\nu_0+3)y\pm2\sqrt{2}$, where we thus impose that $P$ vanishes at all times. 

Similarly, $P_{\mathrm{FPT}}(\widetilde{N}\vert z,y)$ obeys the adjoint Fokker-Planck equation~\cite{Pattison:2017mbe}
\begin{equation}
\label{eq:adjoint:FP}
\frac{\partial P_{\mathrm{FPT}}}{\partial \widetilde{N}} = \mathcal{L}_{\mathrm{FP}}^\dagger P_{\mathrm{FPT}}\, ,
\end{equation}
where $ \mathcal{L}_{\mathrm{FP}}^\dagger$ is the adjoint of the operator introduced in \Eq{eq:LFP:def} and reads
\begin{equation}\begin{aligned}
\label{eq:adjoint:FP:operator}
\mathcal{L}_{\mathrm{FP}}^{\dag} &= \frac{1}{2}\frac{\partial^2}{\partial z^2} -  \frac{\nu_0+\frac{3}{2}}{\nu_0-\frac{3}{2}}\, y\frac{\partial}{\partial y} + z\frac{\partial}{\partial z}\, .
\end{aligned}\end{equation}
When initial conditions are located on the end-of-inflation surface, the first-passage time vanishes, hence \Eq{eq:adjoint:FP} needs to be solved with the boundary condition $P_{\mathrm{FPT}}(N\vert z,y)=\delta(N)$ when $2\sqrt{2\nu_0-3}z/\mu=(2\nu_0+3)y\pm2\sqrt{2}$.


\subsection{Field probability}
\label{sec:field:proba}

Although the two Langevin equations~\eqref{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y} are uncoupled, the variables $z$ and $y$ are intertwined through the boundary condition defined on the end-of-inflation surface. This makes analytical treatments of the stochastic problem difficult to carry out. However, in practice, the role played by the coordinate $y$ can be neglected for the following reason.

The end-of-inflation condition can be rewritten as $z=\pmz_{\mathrm{end}}(1\pm y /y_\mathrm{c})$, where $z_{\mathrm{end}} = \mu/\sqrt{\nu_0-3/2}$ and $y_\mathrm{c}=\sqrt{2}/(\nu_0+3/2)$. The correction $y$ carries in the boundary condition can thus be quantified by the ratio $y_{\mathrm{end}}/y_\mathrm{c}$. According to \Eq{Langevin rotated new variables:y}, $y$ decays exponentially in time, $y_{\mathrm{end}}=\yine^{-\widetilde{\mathcal{N}}(\nu_0+3/2)/(\nu_0-3/2)}$. In the double-well model introduced in \Sec{sec:double:well}, $y_{\mathrm{in}}\simeq 0.24M_\usssPl$, and the typical value for $y_{\mathrm{end}}$ can be assessed by replacing $\widetilde{\mathcal{N}}$ by its expectation value in the above expression. From the adjoint Fokker-Planck equation~\eqref{eq:adjoint:FP}, one can show that the mean number of $e$-folds~starting from initial field values $z_{\mathrm{in}}$ and $y_{\mathrm{in}}$, $\langle \widetilde{\mathcal{N}}\rangle(z_{\mathrm{in}},y_{\mathrm{in}})=\intP_{\mathrm{FPT}}(\widetilde{\mathcal{N}}\vert z_{\mathrm{in}},y_{\mathrm{in}}) \widetilde{\mathcal{N}}\mathrm{d}\widetilde{\mathcal{N}}$, obeys~\cite{Vennin:2015hra}
\begin{equation}
\label{mean number equation}
\mathcal{L}_{\mathrm{FP}}^{\dag}\langle \mathcal{\widetilde{N}}\rangle = -1\, ,
\end{equation}
with boundary condition $\langle\widetilde{\mathcal{N}}\rangle=0$ when starting on the end-of-inflation surface. Since our goal is to show that $y$ can be neglected, let us solve this equation in the case where $y_{\mathrm{in}}=0$, where it has solution
\begin{equation}
\label{mean number e-fold}
\mean{\mathcal{\widetilde{N}}}(z_{\mathrm{in}}) = z_{\mathrm{end}}^2\, {}_2 F_2\left(1,1,\frac{3}{2},2,-z_{\mathrm{end}}^2\right)-z_{\mathrm{in}}^2\, {}_2 F_2\left(1,1,\frac{3}{2},2,-z_{\mathrm{in}}^2\right),
\end{equation}
with $ {}_2 F_2 $ the generalised hypergeometric function. This leads to the ratio $y_{\mathrm{end}}/y_\mathrm{c} \sim y_{\mathrm{in}}/y_\uce^{-\langle\widetilde{\mathcal{N}}\rangle(z_{\mathrm{in}})(\nu_0+3/2)/(\nu_0-3/2)}$ displayed in the left panel of \Fig{fig:ratio time}.
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.49\textwidth]{figs/ratio.pdf}
\includegraphics[width=0.49\textwidth]{figs/ratio-characteristic-time-boundary.pdf}
\caption{Left panel: Relative error to the end-of-inflation condition carried by $y$, as a function of $ z_{\mathrm{end}} $ and for a few values of $ z_{\mathrm{in}} $, using the value $\nu_0\simeq 4.5$ of the double-well model~\eqref{eq:pot:dwi}.
Right panel: Comparison between the critical time $\widetilde{N}_\mathrm{c}$ above which the return probability becomes small, and the mean duration of inflation $\langle\widetilde{N}\rangle(z_{\mathrm{in}})$, in the same situation as in the left panel.}
\label{fig:ratio time}
\end{center}
\end{figure}

Large values of $z_{\mathrm{end}}$ correspond to large values of $\mu$, \textsl{i.e.~} to configurations where the noise is subdominant, see \Eq{Langevin}. In that regime, from \Eq{Langevin rotated new variables:z} $z$ decays exponentially with $\widetilde{N}$ and inflation ends after a time $\widetilde{\mathcal{N}}=\ln(z_{\mathrm{end}}/z_{\mathrm{in}})$. This is why the ratio displayed in \Fig{fig:ratio time} asymptotes $y_{\mathrm{end}}/y_\mathrm{c}\simeq (y_{\mathrm{in}}/y_\mathrm{c})(z_{\mathrm{in}}/z_{\mathrm{end}})^{(2\nu_0+3)/(2\nu_0-3)}$. In practice, current constraints on the CMB tensor-to-scalar ratio~\cite{Planck:2018jri} impose the upper bound $H<7\times 10^{13}\mathrm{GeV}$ at the time CMB scales exit the Hubble radius. Uphill inflation occurs later, hence at a lower value of $H_0$, and using the double-well value $\nu_0\simeq 4.5$ in \Eq{eq:mu:def} thus gives $\mu\gg 3\times 10^3 \sigma^3$. With $\sigma$ typically of order $0.1$, this means that $\mu$, hence $z_{\mathrm{end}}$, is a large parameter. From \Fig{fig:ratio time}, one concludes that, provided $z_{\mathrm{in}}$ is not too large (\textsl{i.e.~} provided one starts close enough to the critical trajectory), $y$ plays a negligible role in the end-of-inflation condition, which simply becomes $z=\pmz_{\mathrm{end}}$. In that limit, the dynamics of $z$ fully decouples and the stochastic problem becomes one-dimensional. \\

The absorbing boundary at $z=\pm z_{\mathrm{end}}$ models the fact that, as inflation ends, the noise is turned off (given that modes stop crossing out the Hubble radius). Hence the system is drifted away from the inflating region in phase space, and cannot re-enter it. Now, for sufficiently large $\mu$, the noise becomes subdominant around the end-of-inflation surface anyway. Therefore, as we will now argue, the absorbing boundary can be neglected, since the realisations it prevents from re-entering inflation are so rare that they give negligible contributions to the quantities of physical interest.
More precisely, in the absence of absorbing boundaries the dynamics of $y$ and $z$ decouple, and one has $P(y,z \vert \widetilde{N},y_{\mathrm{in}},z_{\mathrm{in}}) = \delta[y-y_{\mathrm{in}} e^{- \widetilde{N}(2\nu_0+3)/(2\nu_0-3) }] P^{\mathrm{NB}}(z\vert \widetilde{N},z_{\mathrm{in}})$, where ``NB'' stands for ``no boundaries'' and
\begin{equation}
\label{Green PDF}
P^{\mathrm{NB}}\left(z\Big\vert \widetilde{N},z_{\mathrm{in}}\right) = \frac{e^{-\frac{\left(z-\zine^{\widetilde{N}}\right)^2}{e^{2\widetilde{N}}-1}}}{\sqrt{\pi \left(e^{2\widetilde{N}}-1\right)}}\, .
\end{equation}
This Gaussian distribution is indeed a solution of \Eq{Fokker-Planck PDF} such that $P^{\mathrm{NB}}(z\vert 0,z_{\mathrm{in}})=\delta(z-z_{\mathrm{in}})$, but it does not vanish on the end-of-inflation surface.\footnote{In the absence of boundaries, \Eq{Langevin rotated new variables:z} has solution $z(\widetilde{N})=\zine^{\widetilde{N}}+\Delta z(\widetilde{N})$, where $\Delta z(\widetilde{N}) = e^{\widetilde{N}}\int_{0}^{\widetilde{N}} e^{-\bar{N}}\xi(\bar{N})\mathrm{d}\bar{N}$. Since $\Delta z$ is linearly related to the Gaussian noises $\xi(\bar{N})$, it is itself a Gaussian random variable and its first two moments can be computed using $\langle \xi(\bar{N}) \rangle = 0$ and $\langle \xi(\bar{N}_1) \xi(\bar{N}_2)\rangle = \delta(\bar{N}_1-\bar{N}_2)$. One obtains $\langle \Delta z\rangle = 0$ and $\langle \Delta z^2 \rangle = (e^{2 \widetilde{N}}-1)/2$, hence the solution~\eqref{Green PDF}.} Removing the boundary condition implies that some realisations in \Eq{Green PDF} resume inflation (by crossing back $\pm z_{\mathrm{end}})$ while they should not. The occurrence of such an event can be assessed by considering the probability $p_{\text{re-entry}}$ that, starting from $ z_{\mathrm{end}} $ at initial time, $ z < z_{\mathrm{end}} $ at some later time $ \widetilde{N} $ (\textsl{i.e.~} the system inflates again at time $\widetilde{N}$). Using \Eq{Green PDF}, it is given by
\begin{equation}
\label{P boundary}
p_{\text{re-entry}}\left(\widetilde{N}\right) = \int_{-\infty}^{z_{\mathrm{end}}}P^{\mathrm{NB}}\left(z\Big\vert \widetilde{N},z_{\mathrm{end}}\right)\mathrm{d}z = \frac{1}{2}\left(1-\erf\left\{z_{\mathrm{end}}\left[\exp\left(\widetilde{N}\right)-1\right]\right\}\right)\, .
\end{equation}
As a function of $\widetilde{N}$, this varies between $1/2$ and $0$, with a transition that occurs when the argument of the error function is of order one, that is at the characteristic time $\widetilde{N}_\mathrm{c} = \ln(1+1/z_{\mathrm{end}})$. This means that, when $\widetilde{N}$ is a few times larger than $\widetilde{N}_\mathrm{c}$ or more, the probability that inflation (wrongly) resumes becomes negligible. The  characteristic time $\widetilde{N}_\mathrm{c}$ is compared with the mean number of inflationary $e$-folds~$\langle\widetilde{\mathcal{N}}\rangle$ in the right panel of \Fig{fig:ratio time}. One can see that, when $z_{\mathrm{end}}$ is large, $\langle\widetilde{\mathcal{N}}\rangle\gg \widetilde{N}_\mathrm{c}$, hence the return probability can be neglected and \Eq{Green PDF} provides a good approximation to the full solution indeed. In what follows, the results obtained from \Eq{Green PDF} will be compared to full numerical simulations of the Langevin equations~\eqref{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y}, and this will confirm the validity of these approximations. 
\subsection{First-passage-time probability}
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.49\textwidth]{figs/Langevin-pdf-N_simu=10e6-entree-inflation_z_end=10.pdf}
\includegraphics[width=0.49\textwidth]{figs/Langevin-log-pdf-N_simu=10e6-entree-inflation_z_end=10.pdf}
\caption{First-passage time probability as obtained from the simulation of $ 10^{6} $ realisations of the Langevin equations \eqref{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y} (solid blue line), and according to the approximation~\eqref{Green PFPT} (dashed red line). The initial conditions are set at the intersection between the end-of-inflation surface $ \pi = - \pi_{\mathrm{end}} $ (where inflation resumes) and the $ y $-axis, see figure \ref{fig:Phase-space}, \textsl{i.e.~} at $ z_{\mathrm{in}} = 0 $ and $ y_{\mathrm{in}} = \sqrt{2}/(\nu_0+3/2) $. The parameters are set such that $ \sigma = 0.1 $ and $ z_{\mathrm{end}} = 10$. In the stochastic simulations, the error bars (in green) correspond to $ 2\sigma$-estimates of the statistical error using the jackknife method. The right panel shows the distribution in logarithmic scale, to better display the (lower and upper) tails.}
\label{fig:Langevin-pdf-N_simu=10e6-entrÃ©e-inflation_z_end=10}
\end{center}
\end{figure}
The first-passage-time probability can be computed from the solution~\eqref{Green PDF} of the Fokker-Planck equation by introducing the survival probability. This is the probability that, at time $\widetilde{N}$, a given realisation is still inflating. It can be equally written as the probability that the system is still in the inflating region at time $\widetilde{N}$, or as the probability that the first-passage time through the end-of-inflation surface is larger than $\widetilde{N}$. One thus has
\begin{equation}\begin{aligned}
\label{eq:survival:proba}
\int\mathrm{d} y \int_{-z_{\mathrm{end}}(y)}^{z_{\mathrm{end}}(y)}\mathrm{d} z P\left(z,y\Big\vert \widetilde{N},z_{\mathrm{in}},y_{\mathrm{in}}\right) = \int^\infty_{\widetilde{N}}P_{\mathrm{FPT}}\left(\mathcal{\widetilde{N}}\Big\vert z_{\mathrm{in}},y_{\mathrm{in}}\right)\mathrm{d}\mathcal{\widetilde{N}} \, .
\end{aligned}\end{equation}
By differentiating this expression with respect to $ \widetilde{N} $, and using the Fokker-Planck equation \eqref{Fokker-Planck PDF} where $y$ has been dropped, one obtains
\begin{equation}
\label{PFPT function of PDF}
P_{\mathrm{FPT}}(\widetilde{N}\big\vert z_{\mathrm{in}}) = -\int_{-z_{\mathrm{end}}}^{z_{\mathrm{end}}}\widetilde{\mathcal{L}}_{\mathrm{FP}} P(z\big\vert \widetilde{N},z_{\mathrm{in}})\mathrm{d}z
=\left[\left(z - \frac{1}{2}\frac{\partial}{\partial z}\right) P(z\big\vert \widetilde{N},z_{\mathrm{in}}) \right]_{z=-z_{\mathrm{end}}}^{z=z_{\mathrm{end}}}\, .
\end{equation}
Inserting \Eq{Green PDF} into this expression yields
\begin{equation}\begin{aligned}
\label{Green PFPT}
    P_{\mathrm{FPT}}^{\mathrm{NB}}(\widetilde{N}\big\vert z_{\mathrm{in}}) =& \frac{2\zee^{2\widetilde{N}}}{\sqrt{\pi \left(e^{2\widetilde{N}}-1\right)^3}}
    \exp\left(-\frac{z_{\mathrm{in}}^2e^{2\widetilde{N}}+z_{\mathrm{end}}^2}{e^{2\widetilde{N}}-1}\right)
    \cosh\left[\frac{\zinz_{\mathrm{end}} }{\sinh(\widetilde{N})}\right] \\
    &\times \left\lbrace 1-\frac{z_{\mathrm{in}}}{z_{\mathrm{end}}}e^{-\widetilde{N}}\tanh\left[\frac{\zinz_{\mathrm{end}}}{\sinh(\widetilde{N})}\right]\right\rbrace .
\end{aligned}\end{equation}
This formula is compared with numerical simulations of $10^6$ realisations of the Langevin equations~\eqref{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y} in \Fig{fig:Langevin-pdf-N_simu=10e6-entrÃ©e-inflation_z_end=10}, where $ z_{\mathrm{in}} = 0 $, $\sigma=0.1$ and $H_0 $ is set such that $ z_{\mathrm{end}} = 10 $. The agreement is excellent, which confirms the validity of the two approximations performed in \Sec{sec:field:proba} (namely dropping $y$ and neglecting the effect of the absorbing boundaries).

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.49\textwidth]{figs/Langevin-pdf-N_simu=10e5-entree-inflation_z_end=4.pdf}
\includegraphics[width=0.49\textwidth]{figs/Langevin-pdf-N_simu=10e5-y_0=0_z_end=4.pdf}
\caption{Left panel: first-passage time probability in the same situation as in \Fig{fig:Langevin-pdf-N_simu=10e6-entrÃ©e-inflation_z_end=10}, except that now $z_{\mathrm{end}}=4$. In that case, the one-dimensional approximation breaks down at small values of $\widetilde{N}$. Right panel: same as in the left panel, where initial conditions are now taken as $ y_{\mathrm{in}} = z_{\mathrm{in}} = 0 $ (and the one-dimensional approximation becomes exact). 
In both panels, the black lines correspond to the first-passage time probability as reconstructed from the method of poles, see main text and \App{app:FPT}. }
\label{fig:breaking 1d approx}
\end{center}
\end{figure}
For smaller values of $z_{\mathrm{end}}$, the one-dimensional approximation is less reliable, see the discussion around \Fig{fig:ratio time}. This can be seen in \Fig{fig:breaking 1d approx}, which is similar to \Fig{fig:Langevin-pdf-N_simu=10e6-entrÃ©e-inflation_z_end=10} except that $z_{\mathrm{end}}=4$. In that case, at small values of $\widetilde{N}$, the approximation~\eqref{Green PFPT} breaks down, and the full distribution displays a feature that the one-dimensional approximation fails to capture. This can be understood by noting that if inflation is too short, $y$ has not decayed away when it terminates, hence one can no longer discard it.  

Let us also note that, in previous works~\cite{Pattison:2017mbe, Ezquiaga:2019ftu, Pattison:2021oen, Animali:2022otk}, the first-passage time distribution has mostly been obtained by solving the adjoint Fokker-Planck equation~\eqref{eq:adjoint:FP} for the characteristic function. Using the residue theorem, $P_{\mathrm{FPT}}$ can then be written as a sum of decaying exponentials, with exponents that correspond to the poles of the characteristic function. This approach is implemented in details in \App{app:FPT}, where the poles are extracted in the large-$\mu$ (hence large-$z_{\mathrm{end}}$) limit. One obtains an expression that is very similar to \Eq{Green PFPT}, see \Eq{PDF large zend}, which is displayed with the black line in \Fig{fig:breaking 1d approx}. Although it provides a good fit to the full result, the approximation is not as good as \Eq{Green PFPT}, in particular around the maximum of the distribution. This may seem surprising, given that the pole-based approach correctly implements the boundary conditions, contrary to the method employed above. The reason is that, away from the $\widetilde{N}\gg 1$ tail, a large number of poles needs to be summed over, and the individual errors made when approximating each residue in the large-$z_{\mathrm{end}}$ limit accumulate. This is why the method presented here, which consists in neglecting the stochastic noise around the end-of-inflation surface, is in fact better, and in what follows we stick to this approach.

One may object that this method performs well as long as one is interested in a first-passage time through a condition where the stochastic noise is negligible, which is commonly the case around the end-of-inflation surface but not for the generic point $(x_*,y_*)$ which \Eq{eq:zeta:one:point:pdf} integrates over. Although this is true, for Markovian stochastic processes in one dimension (which is the case here), the first-passage time across \emph{any} point can be written in terms of fist-passage times across the end-of-inflation surface only. The reason is that, if $z_*$ lies between $z_{\mathrm{in}}$ and $z_{\mathrm{end}}$, one has $\widetilde{\mathcal{N}}_{z_{\mathrm{in}}\toz_{\mathrm{end}}}=\widetilde{\mathcal{N}}_{z_{\mathrm{in}}\toz_*}+\widetilde{\mathcal{N}}_{z_*\toz_{\mathrm{end}}}$ where $\widetilde{\mathcal{N}}_{z_{\mathrm{in}}\toz_*}$ and $\widetilde{\mathcal{N}}_{z_*\toz_{\mathrm{end}}}$ are independent variables (due to the Markovian property), hence
\begin{equation}\begin{aligned}
\label{eq:Markov:convolution}
{P_{\mathrm{FPT}}}_{z_{\mathrm{in}} \rightarrow z_{\mathrm{end}}}\left(\widetilde{N}\right) = \int \mathrm{d}\widetilde{N}' {P_{\mathrm{FPT}}}_{z_{\mathrm{in}} \rightarrow z_*}\left(\widetilde{N}'\right) {P_{\mathrm{FPT}}}_{z_* \rightarrow z_{\mathrm{end}}}\left(\widetilde{N}-\widetilde{N}'\right)
\end{aligned}\end{equation}
for any $ z_* $. Thus, ${P_{\mathrm{FPT}}}_{z_{\mathrm{in}} \rightarrow z_*}$ can be expressed in terms of ${P_{\mathrm{FPT}}}_{z_{\mathrm{in}} \rightarrow z_{\mathrm{end}}}$ and ${P_{\mathrm{FPT}}}_{z_* \rightarrow z_{\mathrm{end}}}$ only, by performing de-convolution of the above expression. This can be done using dedicated de-convolution methods (for a review, see \textsl{e.g.~} \Refa{2002PASP..114.1051S}); we do not further explore this direction here but that makes this alternative method of broader interest. 

In what follows, ${P_{\mathrm{FPT}}}_{z_{\mathrm{in}} \rightarrow z_*}$ is rather estimated using the same approximation as above. We start by writing
\begin{equation}\begin{aligned}
\label{eq:survival:proba:generic}
\int_{-\infty}^{z_*} P\left(z\big\vert \widetilde{N},z_{\mathrm{in}}\right) \mathrm{d} z = p_{\mathrm{exit}}(z_{\mathrm{in}},z_*) \int_{\widetilde{N}}^\infty {P_{\mathrm{FPT}}}_{z_{\mathrm{in}} \rightarrow z_*}\left(\widetilde{\mathcal{N}}\right) \mathrm{d}\widetilde{\mathcal{N}}+1-p_{\mathrm{exit}}(z_{\mathrm{in}},z_*)
\end{aligned}\end{equation} 
as in \Eq{eq:survival:proba}, where we assume $z_{\mathrm{in}}\leq z_*$ (the case $z_{\mathrm{in}}>z_*$ can be treated along similar lines). In this expression, $p_{\mathrm{exit}}(z_{\mathrm{in}},z_*)$ corresponds to the probability that, starting from $z_{\mathrm{in}}$, the stochastic process crosses $z_*$ at least once (which is required for the first-passage time to be defined). Contrary to the situation discussed around \Eq{eq:survival:proba}, where two absorbing boundaries were located at $\pmz_{\mathrm{end}}$ and the exit probability was equal to one, here we have a single absorbing boundary and it is not guaranteed that it gets ever crossed. Formally, $p_{\mathrm{exit}}$ can be related to $P$ by noting that the exit probability is nothing but the probability that, in the asymptotic future, one has $z>z_*$, so
\begin{equation}\begin{aligned}
\label{eq:p:exit:P}
p_{\mathrm{exit}}(z_{\mathrm{in}},z_*) = 1-\lim_{\widetilde{N}\to\infty}\int_{-\infty}^{z_*} P\left(z\big\vert \widetilde{N},z_{\mathrm{in}}\right)\mathrm{d} z\, ,
\end{aligned}\end{equation} 
which indeed corresponds to \Eq{eq:survival:proba:generic} in the $\widetilde{N}\to\infty$ limit. By differentiating \Eq{eq:survival:proba:generic} with respect to $\widetilde{N}$, one obtains
\begin{equation}\begin{aligned}
\label{eq:PFPT:zin:z:interm}
p_{\mathrm{exit}}(z_{\mathrm{in}},z_*) {P_{\mathrm{FPT}}}_{z_{\mathrm{in}} \rightarrow z_*}\left(\widetilde{N}\right) = 
-\int_{-\infty}^{z_*}\widetilde{\mathcal{L}}_{\mathrm{FP}} P(z\big\vert \widetilde{N},z_{\mathrm{in}})\mathrm{d}z
=\left(z_* - \frac{1}{2}\frac{\partial}{\partial z_*}\right) P(z_*\big\vert \widetilde{N},z_{\mathrm{in}}) \, .
\end{aligned}\end{equation} 
Inserting \Eq{Green PDF} into \Eqs{eq:p:exit:P} and~\eqref{eq:PFPT:zin:z:interm} finally leads to $p_{\mathrm{exit}}^{\mathrm{NB}}(z_{\mathrm{in}},z_*)=[1+\mathrm{erf}(z_{\mathrm{in}})]/2$ and
\begin{equation}\begin{aligned} 
\label{eq:GFPT:gen}
{P_{\mathrm{FPT}}^{\mathrm{NB}}}_{z_{\mathrm{in}} \rightarrow z_*>z_{\mathrm{in}}}\left(\widetilde{N}\right) = \frac{2}{1+\mathrm{erf}(z_{\mathrm{in}})} \frac{z_*e^{\widetilde{N}}-z_{\mathrm{in}}}{\sqrt{\pi \left(e^{2 \widetilde{N}}-1\right)^3}}
e^{\widetilde{N}-\frac{\left(z_*-\zine^{\widetilde{N}}\right)^2}{e^{2\widetilde{N}}-1}}\, .
\end{aligned}\end{equation} 
Let us stress that this approximation is expected to be reliable only when $z_*$ lies in the drift-dominated region of phase space and not otherwise. One can also check that, when $z_{\mathrm{in}}=0$ and $z_*=z_{\mathrm{end}}$, one recovers \Eq{Green PFPT}, due to the $z\to-z$ symmetry of the problem in that case. 
\section{Backward probability and PBH abundance}
\label{sec:BW and PBH}
Having solved the Fokker-Planck and adjoint Fokker-Planck problems, we are now in a position to evaluate the different terms appearing in \Eq{eq:zeta:one:point:pdf}, and to study phenomenological consequences of that formula such as the probability to form primordial black holes in uphill inflation.
\subsection{Backward probability}
Inserting \Eqs{Green PDF} and~\eqref{Green PFPT} into \Eq{eq:PBW:def}, the backward probability can be approximated by\footnote{The denominator in \Eq{eq:PBW:def} can be computed by integrating \Eq{Green PFPT} over $\widetilde{N}$, but a simpler calculation is to use \Eq{eq:survival:proba} and integrate \Eq{Green PDF} over $z$.}
\begin{equation}\begin{aligned}
\label{Green backward}
P_{\mathrm{BW}}^{\mathrm{NB}}(z\big\vert\widetilde{N}) = \frac{\sqrt{\pi}}{2}\frac{P_{\mathrm{FPT}}^{\mathrm{NB}}\left(\widetilde{N}\big\vertz\right)e^{z^2}\erfc\left(\lvert z \rvert\right)}{\erf\left(\frac{z_{\mathrm{end}}}{\sqrt{e^{2\widetilde{N}}-1}}\right)}
\end{aligned}\end{equation}
where we have set $z_{\mathrm{in}}=0$ (\textsl{i.e.~} initial conditions are set on the critical trajectory), and where we recall that $P_{\mathrm{FPT}}^{\mathrm{NB}}$ is given by \Eq{Green PFPT}.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{figs/Langevin-pdf-BW-N_simu=10e5-entree-inflation_z_end=10_N_BW=1.pdf}
\caption{Backward probability $P_{\mathrm{BW}}$ of the field value $z$, $N_{\mathrm{BW}}=1$ $e$-fold~before the end of inflation (which corresponds to $\widetilde{N}_{\mathrm{BW}}=\nu_0-3/2$). We have set $ z_{\mathrm{in}} = 0 $ and $ z_{\mathrm{end}} = 10 $. The blue line is reconstructed from $10^5$ realisations of the Langevin equations \eqref{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y}, with $2\sigma$-statistical error bars (inferred from jackknife resampling) in green. The red curve corresponds to the no-boundary approximation~\eqref{Green backward}, which is slightly more peaked at the origin. This is partly due to the finite width of the kernel used to reconstruct $P_{\mathrm{BW}}$ from stochastic simulations, so to account for this effect we display in yellow the approximation~\eqref{Green backward} when smoothed with the same kernel.}
\label{fig:PBW}
\end{center}
\end{figure}
This expression is compared with the backward distribution reconstructed from $10^{5}$ simulations of the Langevin equations \eqref{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y} in \Fig{fig:PBW}. The agreement is excellent, except close to $z=0$ where the approximation~\eqref{Green backward} is slightly more peaked. However, when reconstructing $P_{\mathrm{BW}}$ from a finite sample of Langevin realisations, one uses a kernel with a finite width\footnote{In practice, we use the kernel density estimate \texttt{scipy.stats.gaussian\_kde} available in python.}, which has the effect of smoothing out the underlying distribution. To account for this effect, we have smoothed the approximation~\eqref{Green backward} with the same kernel, \textsl{i.e.~} we have drawn $10^5$ points from that distribution and used the same reconstruction technique as for the Langevin sample. The result is displayed in yellow in \Fig{fig:PBW}, which indeed shows better agreement.
\subsection{Curvature perturbation}
The one-point distribution function of the curvature perturbation, coarse-grained at the scale $R$, is given by \Eq{eq:zeta:one:point:pdf}. In the case $z_{\mathrm{in}}=0$, due to the symmetry $z\to-z$ the integral can be restricted to the range $0\leq z_*\leqz_{\mathrm{end}}$, and accounting for the rescaling $ \widetilde{N} = (\nu_0-3/2 )N $ one has
\begin{equation}
P\left(\zeta_R\right) = \left(2\nu_0-3\right)\!\int_{0}^{z_{\mathrm{end}}}\mathrm{d}z_*P_{\mathrm{BW}}\left(z_*\big\vert \widetilde{N}_{\mathrm{BW}}\right){P_{\mathrm{FPT}}}_{0 \rightarrow z_*}\left[\left(\nu_0-\frac{3}{2}\right)\zeta_R-\langle\mathcal{\widetilde{N}}\rangle\left(z_*\right)+\langle\mathcal{\widetilde{N}}\rangle\left(z_{\mathrm{in}}\right)\right].
\end{equation}
Inserting \Eqs{mean number e-fold}, \eqref{eq:GFPT:gen} and~\eqref{Green backward} into that expression, one finds
\begin{equation}
\label{Curvature perturbation starting zero}
P^{\mathrm{NB}}\left(\zeta_R\right) = \left(2\nu_0-3\right)\int_{0}^{z_{\mathrm{end}}}\mathrm{d}z_*\frac{P_{\mathrm{FPT}}^{\mathrm{NB}}\left(\widetilde{N}_{\mathrm{BW}}\big\vert z_*\right)e^{z_*^2}\erfc\left(z_*\right)z_* \alpha\left(z_*,\zeta_R\right)^2\exp\left[-\frac{z_*^2}{\alpha\left(z_*,\zeta_R\right)^2-1}\right]}{\erf\left[z_{\mathrm{end}} \left(e^{2\widetilde{N}_{\mathrm{BW}}}-1\right)^{-\frac{1}{2}}\right] \left[\alpha\left(z_*,\zeta_R\right)^2-1\right]^{\frac{3}{2}}},
\end{equation}
where $\widetilde{N}_{\mathrm{BW}} = (\nu_0-3/2) N_{\mathrm{BW}}(R) = (\nu_0-3/2)\ln(\sigma H_0 R)$ and we have introduced $ \alpha(z_*,\zeta_R) \equiv \exp[(\nu_0-3/2)\zeta_R+z_*^2\, {}_2 F_2(1,1,3/2,2,-z_*^2)] $. The integral over $z_*$ needs to be performed numerically, and the result is compared with the distribution reconstructed from $10^7$ realisations of the Langevin equations~\eqref{Langevin rotated new variables:z} and~\eqref{Langevin rotated new variables:y} in \Fig{fig:curvature-perturbation}.
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.49\textwidth]{figs/curvature_perturbation_z_end=10_different_N_BW_2d.pdf}
\includegraphics[width=0.49\textwidth]{figs/log-curvature_perturbation_z_end=10_different_N_BW_2d.pdf}
\caption{One-point distribution function for the curvature perturbation when coarse-grained at the scale $R=(\sigma H_0)^{-1} e^{N_{\mathrm{BW}}}$, for a few values of $N_{\mathrm{BW}}$, assuming $z_{\mathrm{in}}=0$ and $z_{\mathrm{end}}=10$. The solid curves are reconstructed from $10^7$ realisations of the Langevin equations~\eqref{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y}, with error bars obtained by jackknife resampling. The dashed curves stand for the one-dimensional, no-boundary approximation~\eqref{Curvature perturbation starting zero}, and provide a good fit when $N_{\mathrm{BW}}$ is not too small. The right panel uses a logarithmic scale to better display the tails behaviour (the case $ N_{\mathrm{BW}} = 1 $ is not shown on the left panel to avoid cluttering).}
\label{fig:curvature-perturbation}
\end{center}
\end{figure}
The agreement between the full numerical result and the one-dimensional, no-boundary approximation is good as soon as $N_{\mathrm{BW}}(R)$ is not too large. The reason is that, for large $N_{\mathrm{BW}}$, the backward probability $P_{\mathrm{PW}}(z_*\vert N_{\mathrm{BW}})$ peaks at small values of $z_*$ (\textsl{i.e.~} far from the end-of-inflation surface), around which the stochastic dynamics is not drift-dominated. Since the approximation performed in \Eq{eq:GFPT:gen} relied on this condition, it explains why \Eq{Curvature perturbation starting zero} is less reliable at large $N_{\mathrm{BW}}$, especially on the upper tail (see the case $N_{\mathrm{BW}}=1.5$ in the right panel). To solve this issue, one would have to use the deconvolution method alluded to below \Eq{eq:Markov:convolution}. One also recovers the fact that, at large $\zeta_R$, $P(\zeta_R)$ decays exponentially (hence the upper tails are much heavier than in Gaussian statistics).


One may notice that \Eq{Curvature perturbation starting zero} gives a diverging distribution at $ \zeta_R = 0 $, as can be shown analytically by expanding the result around $\zeta_R=0$. This is done in detail in \App{app:diverging density}, where we find that it is of the form $P(\zeta_R)\propto \zeta_R^{-1/2}$. Therefore, although divergent, $P(\zeta_R)$ remains integrable, hence it is well-defined (and so are all its moments). This peaky behaviour can be observed in the right panel, although it is so sharp that it can hardly be resolved (except for $N_{\mathrm{BW}}=1.5$ where it is slightly more pronounced). It also occurs in the distributions sampled by numerical simulations, but it is smoothed away by the reconstruction kernels. 


\subsection{PBH abundance}
\begin{figure}[t]
\begin{center}
\includegraphics[width=0.49\textwidth]{figs/mass-fraction-different_z_end.pdf}
\includegraphics[width=0.49\textwidth]{figs/mass-fraction_z_end=10_different_z_0.pdf}
\caption{Mass fraction of primordial black holes $ \beta_{\mathrm{f}} $ as a function of the mass $ M $ in gram, for $\zeta_\mathrm{c}=1$. The corresponding value for $ N_{\mathrm{BW}} $ is displayed on the top horizontal axis in the right panel and is related to $M$ via \Eq{eq:mass:NBW} where we have set $\xi=1$ and $\sigma = 0.1$. The value of $ H_0 $ is related to  $ z_{\mathrm{end}} = \mu/\sqrt{\nu_0-3/2} $ using \Eq{eq:mu:def}. On the left panel, $ \beta_{\mathrm{f}} $ is evaluated for different values of $ z_{\mathrm{end}} $ using simulations of the Langevin equations (solid lines, with $95\%$ error bars), compared with the low-mass approximations \eqref{eq:beta:NB} (dashed line) and \eqref{eq:beta:lowmass} (dotted line). On the right panel, $ \beta_{\mathrm{f}} $ is computed via numerical simulations for different values of $ z_{\mathrm{in}} $ and with $ z_{\mathrm{end}} = 10 $.}
\label{fig:mass-fraction}
\end{center}
\end{figure}
One of the phenomenological consequences of the heavy tails observed in $P(\zeta_R)$ is that the probability to form primordial black holes (or other extreme objects such as massive halos~\cite{Ezquiaga:2022qpw}) is enhanced. PBHs form when a large density fluctuation re-enters the Hubble radius and collapses into a black hole. State-of-the-art criteria for PBH formation involve the compaction function~\cite{Shibata:1999zs, Harada:2015yda}, with the mass of the resulting object being assessed via critical scaling relations~\cite{Choptuik:1992jv, Evans:1994pj, Niemeyer:1997mt}. In the present work, for simplicity, we assume that a PBH forms when the local coarse-grained curvature perturbation is larger than a given threshold $\zeta_\mathrm{c}$ of order one, and assess the PBH abundance with the Press-Schechter estimate
\begin{equation}
\label{mass fraction}
\beta_{\mathrm{f}}\left(M\right) = \int_{\zeta_\mathrm{c}}^{\infty} P\left(\zeta_R\right)\mathrm{d}\zeta_R\int_{N_{\mathrm{BW}}(R)}^\infty  P_{\mathrm{FPT}}(\mathcal{N}\big\vert z_{\mathrm{in}})  \mathrm{d}\mathcal{N}
\, .
\end{equation}
In this expression, the first integral corresponds to the probability that, in regions of the universe inflating at least $N_{\mathrm{BW}}(R)$ $e$-folds, a given patch of size $R$ collapses into a black hole, while the second integral is the probability that at least $N_{\mathrm{BW}}(R)$ $e$-folds~are inflated. This second term is usually not considered  since it is close to one in models producing long phases of inflation, but in uphill inflation, which typically lasts for a few $e$-folds~only, it must be added explicitly.\footnote{Formally, this term cancels out with the denominator of \Eq{eq:PBW:def}, which appears in $P\left(\zeta_R\right)$. When using the algorithm detailed at the end of \Sec{sec:StochaDeltaN}, it simply means that trajectories inflating less than $N_{\mathrm{BW}}(R)$ $e$-folds~are thrown away, but their number is recorded in order to post-weigh the mass fraction.}


In \Eq{mass fraction}, $M$ corresponds to some fraction $\xi$ of the mass contained in a Hubble patch at the time $R$ re-enters the Hubble radius. Assuming that PBHs form in the radiation era (which is the case for the scales that emerge during the last few $e$-folds~of inflation if reheating is quick enough), one has
\begin{equation}\begin{aligned}
\label{eq:mass:NBW}
M=\xi \frac{M_\usssPl^2}{\sigma^2 H_0} e^{2 N_{\mathrm{BW}}(R)}\, .
\end{aligned}\end{equation} 
When $z_{\mathrm{in}}=0$, in the no-boundary approximation developed above one can insert \Eqs{Curvature perturbation starting zero} and~\eqref{Green PFPT} into \Eq{mass fraction} and express the integral over $\zeta_R$ in terms of the error function, leading to
\begin{equation}\begin{aligned}
\label{eq:beta:NB}
\beta_{\mathrm{f}}^{\mathrm{NB}}\left(M\right) =  \int_0^{z_{\mathrm{end}}}\mathrm{d} z_*\sqrt{\pi}\,P_{\mathrm{FPT}}^{\mathrm{NB}}\left(\widetilde{N}_{\mathrm{BW}}\big\vert z_*\right)e^{z_*^2}\erfc\left(z_*\right) 
\erf\left[\frac{z_* }{\sqrt{\alpha^2\left(z_*,\zeta_\mathrm{c}\right)-1}}\right]  ,
\end{aligned}\end{equation} 
where we recall that $P_{\mathrm{FPT}}^{\mathrm{NB}}$ is given by \Eq{Green PFPT}.

This expression is compared with the mass fraction extracted from numerical integrations of the Langevin equations~\eqref{Langevin rotated new variables:z}-\eqref{Langevin rotated new variables:y}, following the procedure outlined in \Sec{sec:StochaDeltaN},  in the left panel of \Fig{fig:mass-fraction}. One can check that it gives a good fit at low values of $N_{\mathrm{BW}}$, which is expected since in that case the backward probability peaks at values of $z_*$ that fall in the drift-dominated region, as explained above. In fact, \Eq{eq:beta:NB} can be further expanded in the limit $z_{\mathrm{end}} e^{-\widetilde{N}_{\mathrm{BW}}}\gg 1$, where one obtains
\begin{equation}\begin{aligned}
\label{eq:beta:lowmass}
\beta_{\mathrm{f}}\left(M\right) \simeq  \mathrm{erf}\left[\frac{1}{2}e^{-\left(\nu_0-\frac{3}{2}\right)\zeta_\mathrm{c}-\frac{\gamma_{\mathrm{E}}}{2}}\right]
\end{aligned}\end{equation} 
where $\gamma_{\mathrm{E}}\simeq 0.57$ is the Euler-Mascheroni constant. This formula is displayed with the dotted line in the left panel of \Fig{fig:mass-fraction} and correctly captures the mass fraction at low masses. For larger masses it becomes less accurate, but given that the mass fraction does not vary across orders of magnitude in the relevant range of masses, and since our ability to reproduce its detailed shape is anyway limited by the crude Press-Schechter, curvature-perturbation based method we have employed, only the order of magnitude of $\beta_{\mathrm{f}}$ can be safely assessed here. To that end, \Eq{eq:beta:lowmass} provides a reliable estimate.

A striking feature of \Eq{eq:beta:lowmass} is that it neither depends on the mass, nor on $z_{\mathrm{end}}$, \textsl{i.e.~} on the energy scale at which uphill inflation occurs. Usually, the abundance of black holes increases with the energy density during inflation, since this sets the size of cosmological perturbations. Here, the details of the background dynamics are such that this dependence cancels out, and one obtains mass fractions of order $10^{-2}$ (with $\nu_0\simeq 4.5$ as obtained in the double-well model described in \Sec{sec:double:well}). During the radiation era, the abundance of black holes (which behave as pressure-less matter) grows linearly with the scale factor, hence after a few $e$-folds~of radiation, PBHs dominate the energy budget of the universe. This conclusion is, again, independent of the energy scale at which hilltop inflation proceeds.
In this scenario, PBHs are produced with low masses, smaller than $10^9\, \mathrm{g}$, so they Hawking-evaporate before big-bang nucleosynthesis. This implies that the only imprint left from the PBH-dominated phase is a background of gravitational waves, as studied recently in \Refs{Inomata:2020lmk, Papanikolaou:2020qtd, Domenech:2020ssp}. 
\section{Conclusion}
\label{sec:Conclusions}
Let us now summarise our main findings. In this work, we have investigated uphill inflation, where inflation proceeds as the inflaton climbs up a local maximum in its potential energy. If the initial velocity of the inflaton is close to the value classically required to freeze near the local maximum, we have found that inflation proceeds in a new regime, which is neither close to slow roll, fast roll nor ultra-slow roll, and where all three terms in the Klein-Gordon equation are of comparable size. In that regime, although the second and third Hubble-flow functions are singular at the phase-space origin, the Mukhanov-Sasaki effective mass remains regular, and the curvature perturbation $\zeta$ grows as $a^9$, i.e. much faster than in slow roll ($\zeta$ constant) or ultra slow roll ($\zeta\propto a^3$).

This suggests that quantum fluctuations destabilise the freezing of the inflaton at the local maximum by exciting the unstable direction of phase space, and we have studied this mechanism using the stochastic-inflation formalism. Although two-dimensional, the Langevin equations are tractable since only the unstable direction receives a contribution from the noise, while the stable direction remains classical. From the results of \Refa{Tada:2021zzj} we have derived a simple algorithm to extract the statistics of the coarse-grained curvature perturbation $\zeta_R$, which is numerically not more expensive than the computation of a mere first-passage time distribution, and which is fully generic. We have applied this method to the problem at hand and obtained the one-point distribution of $\zeta_R$, as well as the probability for it to exceed the threshold value required to form primordial black holes.

We have also interpreted these results by means of analytical approximations. To that end, we noticed that instead of solving the adjoint Fokker-Planck problem, which requires to impose non-trivial boundary conditions, the first-passage time distributions can be obtained from the solutions of the Fokker-Planck equation directly. The Fokker-Planck equation also comes with absorbing boundary conditions on the end-of-inflation surface, however they can be neglected if one assumes that the stochastic noise is subdominant when inflation ends, which in practice we found to be an excellent approximation. We compared it to the more traditional characteristic-function method, and found that it provides more accurate, and remarkably simpler to derive, results. 

The distribution function of $\zeta_R$ features exponential tails, as is now known to be a ubiquitous consequence of quantum diffusion. The abundance of primordial black holes was found to be independent of the energy scale at which uphill inflation proceeds, and to only involve the curvature of the potential in Hubble units, which is of order unity when uphill inflation is realised in double-well models as explained in \Sec{sec:double:well}. In that case, PBHs are formed with an abundance of order $10^{-2}$, so they quickly dominate the universe content after the first few $e$-folds~of the radiation era. Being formed with ultra-light masses, they evaporate before big-bang nucleosynthesis, but leave a stochastic gravitational-wave background behind~\cite{Inomata:2020lmk, Papanikolaou:2020qtd, Domenech:2020ssp}, possibly as their only imprint. \\

Let us note that, in double-well potentials, the initial velocity required for the field to classically freeze at the local maximum of the potential corresponds to the one inherited from the preceding attractor phase of slow-roll inflation, if the width $\f_{\mathrm{0}}$ matches the critical value $\f_{\mathrm{0}}^\mathrm{cri}$ given in \Eq{eq:phi0:cri}. This parameter thus needs to be sufficiently fine tuned, which may not come as a surprise since models producing PBHs are  generically known to rely on some fine tuning. Before closing this paper, let us try to better assess the amount of fine tuning that is required.

\begin{figure}[t]
\begin{center}
\includegraphics[width=0.7\textwidth]{figs/N_average_function_phi_0.pdf}
\caption{Same as in the right panel of \Fig{fig:DWIpot}, \textsl{i.e.~}~number of inflationary $e$-folds~realised between the first two crossings of a local minimum in the double-well model discussed in \Sec{sec:double:well}, where the mean number of stochastic $e$-folds~is also displayed. On the top horizontal axis, the value of $x_{\mathrm{in}}$ that corresponds to $\f_{\mathrm{0}}$ is displayed, see main text.}
\label{fig:duration inflation stochastic}
\end{center}
\end{figure}
If $\f_{\mathrm{0}} = \f_{\mathrm{0}}^\mathrm{cri} $, initial conditions are set on the unstable branch exactly, \textsl{i.e.~} $x_\mathrm{in}=0$. A deviation of $\f_{\mathrm{0}}$ from its critical value can thus be described by a non-vanishing value for $x_\mathrm{in}$, which simply translates the fact that by tweaking the potential, one may either overshoot the local maximum or turn back before reaching it. In the first case, corresponding to $\f_{\mathrm{0}}<\f_{\mathrm{0}}^\mathrm{cri}$, the effective value of $x_\mathrm{in}$ can be identified by matching the field velocity when crossing the potential maximum. In the second case, corresponding to $\f_{\mathrm{0}}>\f_{\mathrm{0}}^\mathrm{cri}$, it can be done by matching the field value when turning back. This leads to the function $x_\mathrm{in}(\f_{\mathrm{0}})$ used to label the upper horizontal axis in \Fig{fig:duration inflation stochastic}. There, the mean duration of uphill inflation is displayed, and one can check that the classical divergence observed in \Fig{fig:DWIpot}, and reproduced here with the blue curve for convenience, is renormalised away by quantum diffusion. When $\f_{\mathrm{0}}$ deviates more than a few percents from its critical value, the duration of uphill inflation is reduced, and stochastic effects become less important, which explains why the stochastic and classical curves gets closer to each other (the reason why they do not exactly match far from $\f_{\mathrm{0}}^\mathrm{cri}$ is because, in the classical problem, the full Klein-Gordon and Friedmann equations are solved, while in the stochastic problem their linearised version around the phase-space origin is employed).

From \Fig{fig:duration inflation stochastic} one can also see that a few-percent deviation for $\f_{\mathrm{0}}$ corresponds to working with $x_{\mathrm{in}}$ of order a few percents too, which in turn corresponds to $z_{\mathrm{in}}/z_{\mathrm{end}}$ of order a few percents [using that $z_{\mathrm{in}}/z_{\mathrm{end}} = (\nu_0/2-3/4)x_{\mathrm{in}}$]. This is why, in the right panel of \Fig{fig:mass-fraction}, we show the PBH mass fraction obtained from solving the coupled Langevin equations with two non-vanishing values of $z_{\mathrm{in}}/z_{\mathrm{end}}$, namely $0.01$ and $0.1$. We observe that, although the mass fraction decreases with $z_{\mathrm{in}}/z_{\mathrm{end}}$ as expected, even for $z_{\mathrm{in}}/z_{\mathrm{end}}=0.1$ the abundance of black holes is still of a similar order of magnitude. This corresponds to working with $x_{\mathrm{in}}\simeq 0.07$, hence a value of $\f_{\mathrm{0}}$ that is $10\%$ larger than $\f_{\mathrm{0}}\equiv\f_{\mathrm{0}}^\mathrm{cri} $. As a consequence, the PBH production mechanism seems to be subject to much less fine tuning in uphill inflation than in most alternative models. This makes it an attractive scenario for the production of ultra-light PBHs in the early universe. 


\section*{Acknowledgements}
It is a pleasure to thank Louis-Pierre Chaintron for valuable discussions on stochastic processes.



