# The hand-made tail:
Non-perturbative tails from multifield inflation

Ana Achúcarro  Sebastian Céspedes  Anne-Christine Davis  Gonzalo A. Palma

###### Abstract

It is becoming increasingly clear that large but rare fluctuations of the
primordial curvature field, controlled by the tail of its probability
distribution, could have dramatic effects on the current structure of the
universe —e.g. via primordial black-holes. However, the use of standard
perturbation theory to study the evolution of fluctuations during inflation
fails in providing a reliable description of how non-linear interactions
induce non-Gaussian tails. Here, we use the stochastic inflation formalism to
study the non-perturbative effects from multi-field fluctuations on the
statistical properties of the primordial curvature field. Starting from the
effective action describing multi-field fluctuations, we compute the joint
probability density function and show that enhanced non-Gaussian tails are a
generic feature of slow-roll inflation with additional degrees of freedom.

##  1 Introduction

The reconstruction of our universe’s history relies on the assumption that the
primordial curvature fluctuation $\zeta({\vec{x}})$ (responsible for our
universe’s inhomogeneities) was initially distributed according to a Gaussian
statistics, parametrised by an almost scale invariant power spectrum
$P_{\zeta}(k)$. Although this assumption agrees with every relevant
cosmological observation [1], there are good reasons to suspect that our
primordial universe could not have been perfectly Gaussian. To start with, the
simplest models of cosmic inflation —the theory that explains the origin of
$\zeta({\vec{x}})$— predict tiny, but non-vanishing, levels of non-
Gaussianities [2]. Unfortunately, this minimal prediction will likely remain
out of reach for the next generation of cosmological surveys. On the other
hand, large non-Gaussianity can arise if, during inflation, $\zeta({\vec{x}})$
evolved experiencing large self-interactions and/or interactions with other
relevant degrees of freedom [3, 4, 5, 6, 7]. One way to parametrise the
observable effects of these interactions on the distribution of
$\zeta({\vec{x}})$ is in the form of $n$-point correlation functions
$\langle\zeta(\vec{x}_{1})\cdots\zeta(\vec{x}_{n})\rangle$. The shape of these
$n$-point functions in momentum space can display distinctive signatures,
providing a powerful diagnostic of the types of fields present during
inflation. For example, massive fields with spin can leave oscillatory
features in the primordial bispectrum (the amplitude of the 3-point
correlation function of primordial fluctuations) with a shape determined by
their spin [8, 9, 10, 11].

However, $n$-point correlation functions computed with standard perturbation
theory are inappropriate to assess the occurrence of large statistical
excursions of $\zeta({\vec{x}})$. The prevalence of large statistical
excursions is dictated by the shape of the tail of the probability
distribution function describing the statistics of $\zeta({\vec{x}})$. But
perturbative methods fail to correctly determine the profile of tails. As
emphasised in [12], perturbation theory schematically relates the 1-point
probability distribution of $\zeta$ and connected $n$th-moments
$\langle\zeta^{n}\rangle_{c}$ [which, in turn, are related to connected
$n$-point correlation functions of $\zeta({\vec{x}})$] as

$\displaystyle P(\zeta)$ $\displaystyle\sim$
$\displaystyle\exp\left[-\frac{\zeta^{2}}{2\sigma_{\zeta}^{2}}+\sum_{n=3}^{\infty}\frac{\langle\zeta^{n}\rangle_{c}}{\sigma_{\zeta}^{2n}}\zeta^{n}\right],$
(1.1)

where $\sigma_{\zeta}^{2}$ is the Gaussian variance of the distribution,
determined by the power spectrum $P_{\zeta}(k)$ as $\sigma_{\zeta}^{2}=\int
d^{3}kP_{\zeta}(k)$. In terms of the usual $f_{\rm NL}$ and $g_{\rm NL}$
parameters for the first few terms in the expansion, the previous expression
takes the form

$\displaystyle P(\zeta)$ $\displaystyle\sim$
$\displaystyle\exp\left[-\frac{\zeta^{2}}{2\sigma_{\zeta}^{2}}\left(1+f_{\rm
NL}\zeta+g_{\rm NL}\zeta^{2}+\cdots\right)\right].$ (1.2)

For typical statistical excursions $\zeta\sim\sigma_{\zeta}\ll 1$, the
expansion of the distribution function in terms of moments
$\langle\zeta^{n}\rangle_{c}$ remains under control as long as
$\langle\zeta^{n}\rangle_{c}/\sigma_{\zeta}^{n}\ll 1$, which can be satisfied
in perturbation theory even for values of $f_{\rm NL}$ and $g_{\rm NL}$ of
order 1. On the other hand, for unlikely large statistical excursions
$\zeta\sim 1$ this expansion may fail, particularly in models predicting
$f_{\rm NL}$ and $g_{\rm NL}$ of order 1 (or larger) commonly encountered in
theories of inflation involving sizable non-linear interactions (for instance,
in the form of interactions with other degrees of freedom). In such models,
not only $\langle\zeta^{3}\rangle_{c}$ and $\langle\zeta^{4}\rangle_{c}$, but
all $n$-point moments are expected to contribute corrections of order $1$ on
the tail of the distribution, making the expansion (1.2) useless to study
extreme fluctuations. This failure of perturbation theory to parametrise large
statistical excursions of $\zeta$ in certain models of inflation motivates the
consideration of non-perturbative techniques to study the consequence of non-
linear interactions of $\zeta$ during inflation [13, 14, 15, 16, 12, 17, 18].

As unlikely as they might be, large statistical fluctuations of the primordial
field $\zeta$ can have dramatic effects on the formation of structure in our
universe. More to the point, after inflation, large fluctuations of
$\zeta({\vec{x}})$ can lead to overdense regions of space that inevitably
collapse into primordial black holes (PBHs) (see Refs. [19, 20] for recent
reviews). These black holes could become the seeds of supermassive black holes
at the center of galaxies, and form a substantial part of the dark matter
content of our universe. The abundance and clustering properties of these PBHs
are extremely sensitive to the shape of the tails of the PDF dictating the
distribution of $\zeta({\vec{x}})$ [21, 22, 23, 24]. Thus, to correctly
understand the possible generation of PBHs as a result of inflation, we need a
reliable, non-perturbative approach to reconstruct the non-Gaussian tails of
the primordial fluctuation’s PDF.

The purpose of this work is to quantify precisely the effects of light
isocurvature fluctuations on the probability density function of
$\zeta({\vec{x}})$ by using the non-perturbative approach offered by the
stochastic inflation formalism [25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35,
36, 37, 38, 39, 40, 41]. The authors of [42] have argued that the interaction
between $\zeta$ and a light field $\psi$ can introduce non-Gaussian
corrections that modify the shape of tails of the probability density function
of $\zeta({\vec{x}})$, valid at the end of inflation. Here we confirm this
scenario and we show that the joint PDF describing the statistic of $\zeta$
and $\psi$ in two-field models of inflation with canonical kinetic terms is
given by

$\displaystyle
P(\zeta,\psi)\sim\exp\left[-\frac{\psi^{2}}{2\sigma_{\psi}^{2}}-\frac{1}{2\sigma_{\zeta}^{2}}\left(\zeta-\kappa\frac{\psi^{2}}{2\sigma_{\psi}^{2}}\right)^{2}+\cdots\right],$
(1.3)

where $\kappa$ is related to the strength of the coupling between $\zeta$ and
$\psi$, and the ellipses stand for additional subleading contributions that we
calculate in some specific examples. The non-perturbative nature of (1.3) is
not obvious, but it becomes apparent after integrating over $\psi$ to reveal
the tail of the distribution for $\zeta$, which becomes strongly non-Gaussian

$P(\zeta)\sim\exp(-\zeta/\kappa).$ (1.4)

The dependence on $\kappa$ makes manifest the non-perturbative sensitivity of
tails to non-linear interactions between $\zeta$ and other degrees of freedom.
Similar non-Gaussian tails have been found in other single-field scenarios
where the background is non trivial [43, 44, 45, 46], and quantum diffusion
plays an important role. And in [42], for a sudden, transient coupling between
the curvature field and a light spectator field. Instead, in our calculation,
slow-roll is preserved throughout and all interactions and couplings are
constant in time. In multi-field models of slow-roll inflation with non-
canonical kinetic terms we expect corrections in (1.3) that would change the
details on how (1.4) is obtained, leading to a different profile for the tail.

An important advantage of the analysis presented here, based on the
fluctuations, is that we are able to show that such non-Gaussian tails are a
generic consequence of multifield inflation, and also that they do not require
the interruption of slow roll. The leading non-Gaussian contribution to the
PDF can be traced back to an ever-present quadratic derivative coupling111This
term is always present unless the inflationary trajectory follows a geodesic
in field space $\dot{\zeta}\psi$ between the curvature and the isocurvature
perturbations [47, 48]. In minimal multifield scenarios the coupling $\kappa$
is related to the angular velocity $\Omega$ of the inflationary trajectory in
field space. However, our results apply to any model where this derivative
coupling is present.

We will be particularly interested in the case of a very light or even
massless (ultralight) isocurvature fluctuation [49, 50]. This provides an
alternative inflationary scenario –potentially relevant to string
compactifications– with predictions currently indistinguishable from those of
single-field inflation but where light fields do not need to be stabilized. We
consider the background to be quasi de Sitter during the whole duration of
inflation. The UV completion of such systems in terms of an effective field
space metric and an effective multifield potential has been discussed in [10,
50, 51, 52].

To derive (1.3), we start from the effective action for the perturbations of a
two field model of inflation [47, 48]. From the action, we will coarse grain
the equations of motion to obtain a Fokker-Planck equation for the long
wavelength modes. In order to introduce derivative interactions we write the
equation in phase space. Because the time scale associated with the approach
to equilibrium of the velocity field $v_{\zeta}$ is much shorter than the one
of the other fields, we can integrate out $v_{\zeta}$ directly from the
Fokker-Planck equation. This leads to Eq. (4.48) which involves only $\psi$
and $\zeta$. This equation assumes that the entropy mass of the second field
is light, that the coupling $\Omega^{2}<H^{2}$ and that the curvature power
spectrum is smaller that one, and to our knowledge has not been previously
derived. Surprisingly enough, it will be possible to solve the time dependent
Fokker-Planck equation on a myriad of cases, which among other consequences
show that the derivative coupling $\Omega$, both enhances the Gaussian
variance and modifies the PDF introducing a coupling
$\kappa\sim\Omega^{2}/H^{2}$.

The structure of the paper is as follows. In Section 2 we study the statistics
of primordial curvature perturbations and we show how to integrate out
$v_{\zeta}$. We also review some known results in the case of spectator fields
on fixed de Sitter. In Section 3 we present the linear Fokker-Planck equation
for the curvature perturbation coupled to another light field. Since the
distribution is Gaussian is possible to obtain exact expressions for the
variances of the fields, that as we will show, match known results using
standard techniques. Section 4 contains the main results of this paper, where
we study the full non linear Fokker-Planck equation. After integrating out
$v_{\zeta}$ we will show how it is to obtain (1.3) and under which assumptions
it holds. Finally, in Section 5 we conclude and present different ideas to
explore in the future. There are a series of appendices where we present
technical details of the calculations.

##  2 Statistics of primordial curvature perturbations

Before studying the effects of isocurvature fields on the statistics of
$\zeta$, we first review the use of the stochastic formalism, showing how it
allows a derivation of the probability density function describing the
statistics of single fields in quasi-de Sitter backgrounds.

###  2.1 Primordial curvature perturbation

We start by considering the task of deriving the probability distribution of
the primordial curvature perturbation $\zeta$. First, let us recall that the
canonical quadratic action for $\zeta$ describing its dynamics during
inflation is given by

$S={M_{\mathrm{Pl}}^{2}}\int
d^{4}xa^{3}\epsilon\left[\dot{\zeta}^{2}-\frac{1}{a^{2}}(\nabla\zeta)^{2}\right],$
(2.1)

where $M_{\rm Pl}$ is the reduced Planck mass. In the previous expression,
$a=a(t)$ is the usual scale factor, and $\epsilon$ is the first slow-roll
parameter, determined by the Hubble parameter $H=\dot{a}/a$ as
$\epsilon=-\dot{H}/H^{2}$, and required to be much smaller than 1 throughout
inflation. For simplicity, we will disregard slow-roll corrections and take
both $H$ and $\epsilon$ as constants. Then, the equation of motion for the
long wavelength modes (with wavelengths much larger than the Hubble radius
$H^{-1}$) is given by

$\ddot{\zeta}+3H\dot{\zeta}=0.$ (2.2)

One can modify this equation to quantify the influence of short wavelength
fluctuations on the evolution of $\zeta$ by introducing a source term
representing noise [53]. The resulting equation takes the form

$\ddot{\zeta}+3H\dot{\zeta}=3H\eta_{\zeta},$ (2.3)

where $\eta_{\zeta}=\eta_{\zeta}(t)$ is a time-dependent Gaussian noise with a
two-point correlator given by:

$\langle\eta_{\zeta}(t)\eta_{\zeta}(t^{\prime})\rangle=\frac{H^{3}}{8\pi^{2}\epsilon{M_{\mathrm{Pl}}}^{2}}\delta(t-t^{\prime}).$
(2.4)

Equation (2.3) allows one to obtain a Fokker-Planck equation satisfied by the
probability density function (PDF) $P(\zeta)$ describing the statistics of
long wavelength modes. In order to see this, it is useful to rewrite (2.3) in
terms of the following two first order differential equations

$\dot{\phi}_{i}=\sum_{j}A_{ij}(t)\phi_{j}+f_{i}(t),$ (2.5)

where we have identified $\phi_{1}=\zeta$, and $\phi_{2}=\dot{\zeta}$.
Equation (2.5) is a Langevin equation with a drift matrix $A_{ij}$ and noise
vector $f_{i}$ given by

$A_{ij}=\begin{pmatrix}0&1\\\ 0&-3H\end{pmatrix},\qquad
f_{i}=3H\eta_{\zeta}\begin{pmatrix}0\\\ 1\end{pmatrix}.$ (2.6)

From (2.4) it follows that the noise vector $f_{i}$ must satisfy $\langle
f_{i}(t)f_{j}(t^{\prime})\rangle=D_{ij}\delta(t-t^{\prime})$, where $D_{ij}$
is the diffusion matrix, given by

$D_{ij}=\begin{pmatrix}0&0\\\ 0&D_{\zeta}\end{pmatrix},\qquad
D_{\zeta}\equiv\frac{9H^{5}}{8\epsilon\pi^{2}}.$ (2.7)

In general, one might be interested in computing correlation functions of the
stochastic fields fields $\phi_{i}$ of Eq. (2.5). These can be computed with
the help of a probability density function $P(\phi^{i},t)$ derived from the
associated Fokker-Planck equation [54]. The Fokker-Planck equation is
determined by $A_{ij}$ and $D_{ij}$ as

$\displaystyle\frac{\partial P}{\partial
t}+A_{ij}\frac{\partial}{\partial\phi_{i}}(\phi_{j}P)-\frac{1}{2}D_{ij}\frac{\partial^{2}}{\partial\phi_{i}\partial\phi_{j}}P=0.$
(2.8)

Using (2.6) and (2.7), we re-express the Fokker-Planck equation in terms of
$\phi_{1}=\zeta$ and $\phi_{2}=v_{\zeta}\equiv\dot{\zeta}$ as

$\frac{\partial P}{\partial
t}+\frac{\partial}{\partial\zeta}\left(v_{\zeta}P\right)-3H\frac{\partial}{\partial
v_{\zeta}}(v_{\zeta}P)-\frac{D_{\zeta}}{2}\frac{\partial^{2}}{\partial
v_{\zeta}^{2}}P=0.$ (2.9)

To solve this equation, let us assume a general Gaussian profile of the form

$\displaystyle P(\zeta,v_{\zeta},t)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi\sqrt{\det
S^{-1}}}\exp\left(-\frac{1}{2}\sum_{ij}S^{-1}_{ij}\phi_{i}\phi_{j}\right),$
(2.10)

where $S^{-1}_{ij}$ are the elements of the (symmetric) covariance matrix,
whose inverse $S_{ij}$ is constituted by two-point moments as

$S_{ij}(t)=\begin{pmatrix}\langle\zeta^{2}\rangle(t)&\langle\zeta
v\rangle(t)\\\ \langle\zeta v\rangle(t)&\langle v^{2}\rangle(t)\end{pmatrix}.$
(2.11)

The time dependence of $S$ is determined by (2.9) together with initial
conditions. To determine $S$ we can take the fields $\zeta$ and $v_{\zeta}$ to
be coordinates with Fourier transforms $p$ and $q$ respectively. Then, the
Fourier transformed version of (2.9) is

$\displaystyle\frac{\partial\tilde{P}}{\partial t}-p\frac{\partial}{\partial
q}\tilde{P}+3Hq\frac{\partial}{\partial
q}\tilde{P}+\frac{D_{\zeta}}{2}q^{2}\tilde{P}=0,$ (2.12)

where $\tilde{P}$ represents the Fourier transform of $P$. The ansatz given in
Eq. (2.10) then implies the following form for $\tilde{P}$:

$\tilde{P}(p,q,t)=\exp\left(-\frac{1}{2}\left[S_{\zeta\zeta}p^{2}+2S_{\zeta
v}pq+S_{vv}q^{2}\right]\right).$ (2.13)

Replacing this expression back into (2.12) we get the following set of
equations satisfied by the elements of the matrix $S$:

$\displaystyle\frac{1}{2}\dot{S}_{\zeta\zeta}$ $\displaystyle=S_{\zeta v},$
(2.14) $\displaystyle\dot{S}_{\zeta v}+3HS_{\zeta v}$ $\displaystyle=S_{vv},$
(2.15) $\displaystyle\dot{S}_{vv}+6HS_{vv}$ $\displaystyle=D_{\zeta}.$ (2.16)

To solve these equations, we need to impose initial conditions at a given time
$t_{0}$. For instance, consider an initial Gaussian distribution (2.10) such
that at $t=t_{0}$ the matrix $S$ contains initial values

$S_{ij}(t_{0})=\begin{pmatrix}S^{(0)}_{\zeta\zeta}&S^{(0)}_{\zeta v}\\\
S^{(0)}_{v\zeta}&S^{(0)}_{vv}\end{pmatrix}.$ (2.17)

Solving Eqs. (2.14)-(2.16) with these initial conditions, we then arrive at:

$\displaystyle S_{\zeta\zeta}(t)$ $\displaystyle=$
$\displaystyle-\frac{2}{3H}\left(S_{\zeta
v}^{(0)}-\frac{D_{\zeta}}{9H^{2}}\right)(e^{-3H(t-t_{0})}-1)+\frac{1}{9H^{2}}\left(S_{vv}^{(0)}-\frac{D_{\zeta}}{6H}\right)(e^{-6H(t-t_{0})}-1)$
(2.18)
$\displaystyle+\left(\frac{2}{3H}S_{vv}^{(0)}+\frac{D_{\zeta}}{9H^{2}}\right)(t-t_{0})+S_{\zeta\zeta}^{(0)},$
$\displaystyle S_{\zeta v}(t)$ $\displaystyle=$ $\displaystyle\left(S_{\zeta
v}^{(0)}-\frac{D_{\zeta}}{9H^{2}}\right)e^{-3H(t-t_{0})}-\frac{1}{3H}\left(S_{vv}^{(0)}-\frac{D_{\zeta}}{6H}\right)e^{-6H(t-t_{0})},$
(2.19) $\displaystyle+\frac{1}{3H}S_{vv}^{(0)}+\frac{D_{\zeta}}{18H^{2}},$
$\displaystyle S_{vv}(t)$ $\displaystyle=$
$\displaystyle\left(S_{vv}^{(0)}-\frac{D_{\zeta}}{6H}\right)e^{-6H(t-t_{0})}+\frac{D_{\zeta}}{6H}.$
(2.20)

The initial values $S_{\zeta\zeta}^{(0)}$, $S_{\zeta v}^{(0)}$ and
$S_{vv}^{(0)}$ are the variances associated with long wavelength fluctuations
that have already crossed the horizon prior to $t_{0}$. If we are interested
only in the statistics of those modes that cross the horizon starting at
$t_{0}$, we can set the initial values of $S_{ij}$ to $0$. This corresponds to
a distribution where the position $\zeta$ and rapidity $v_{\zeta}$ of the
fluctuation are exactly localized at the origin of the field phase space
$P(\zeta,v_{\zeta},t_{0})=\delta(\zeta)\delta(v_{\zeta})$. Then, the solutions
take the form:

$\displaystyle S_{\zeta\zeta}(t)$
$\displaystyle=\frac{D_{\zeta}}{54H^{3}}\left(-3+6H(t-t_{0})+4e^{-3H(t-t_{0})}-e^{-6H(t-t_{0})}\right),$
(2.21) $\displaystyle S_{\zeta v}(t)$
$\displaystyle=\frac{D_{\zeta}}{18H^{2}}\left(1-e^{-3H(t-t_{0})}\right)^{2},$
(2.22) $\displaystyle S_{vv}(z)$
$\displaystyle=\frac{D_{\zeta}}{6H}(1-e^{-6H(t-t_{0})}).$ (2.23)

Replacing these expressions back into (2.10) we obtain the desired expression
for the distribution $P$. The coefficients $S_{\zeta\zeta}$, $S_{\zeta v}$ and
$S_{vv}$ depend on time with a characteristic timescale determined by
$H^{-1}$. In the limit $t-t_{0}\gg H^{-1}$ the distribution simplifies to an
asymptotic expression given by

$P(\zeta,v,t)=\frac{1}{2\pi}\left(\frac{54H^{3}}{D_{\zeta}^{2}(t-t_{0})}\right)^{1/2}\exp\left[-\frac{9H^{2}}{2D_{\zeta}(t-t_{0})}\zeta^{2}+\frac{3H}{D_{\zeta}(t-t_{0})}\zeta
v_{\zeta}-\frac{6H}{2D_{\zeta}}v_{\zeta}^{2}\right].$ (2.24)

Notice that the widths associated with $\zeta$ and $v_{\zeta}$ differ in their
time dependence, with $v_{\zeta}$ sharply localized around $0$ (signaling that
$v_{\zeta}$ decays quickly after it becomes super-horizon). In fact, we may
marginalize $v_{\zeta}$ by integrating it from the distribution (2.24), in
which case we obtain

$P(\zeta,t)=\frac{1}{\sqrt{2\pi}\sigma_{\zeta}}\exp\left(-\frac{1}{2\sigma_{\zeta}^{2}}\zeta^{2}\right),$
(2.25)

which is a Gaussian distribution for $\zeta$ with variance
$\sigma_{\zeta}^{2}$ given by

$\sigma_{\zeta}^{2}=\frac{D_{\zeta}}{9H^{2}}(t-t_{0})=\frac{H^{3}}{8\epsilon\pi^{2}}(t-t_{0}).$
(2.26)

Recall that this expression is valid for $t\gg H^{-1}$, provided the initial
condition $S_{ij}=0$ at $t_{0}=0$. The time dependence of the variance
$\sigma_{\zeta}^{2}$ just reflects the fact that as time progresses, more and
more modes populate the long wavelength regime. In this way, at a given time
$t$, the probability distribution $P(\zeta,t)$ describes the statistics of
long-wavelengths that crossed the horizon between $t_{0}$ and $t$.

###  2.2 Integrating out $v_{\zeta}$

Although $\zeta$ and $v_{\zeta}$ had the same status in the treatment leading
to (2.24) we are ultimately interested only in the statistics of $\zeta$
(after all, $v_{\zeta}$ decays quickly on super-horizon scales). This led us
to derive (2.25) after marginalizing $v_{\zeta}$. Alternatively, we can
integrate out $v_{\zeta}$ at an early stage, and obtain a Fokker-Planck
equation only for $\zeta$. For instance, if we neglect the second derivative
of $\zeta$ in (2.5), the system reduces to a single Langevin equation
$\dot{\zeta}=\eta_{\zeta}$, yielding the following Fokker-Planck equation

$\displaystyle\frac{\partial P}{\partial
t}=\frac{D_{\zeta}}{18H^{2}}\frac{\partial^{2}P}{\partial\zeta^{2}},$ (2.27)

whose solution is precisely given by Eq. (2.25). Notice that (2.9) included a
term involving a second derivative with respect to $v_{\zeta}$, and only a
first derivative with respect to $\zeta$. In contrast, Eq. (2.27) contains a
second derivative with respect to $\zeta$. This can be understood as the
effect of integrating out $v_{\zeta}$ over the Fokker-Planck equation (2.9)
with the ideas of Refs. [55, 54], which go as follows: First, by noticing that
the time scale in which $v_{\zeta}$ becomes time independent is given by
$t_{v}=1/3H$, we can rewrite Eq. (2.9) as:

$\displaystyle\frac{\partial}{\partial
v_{\zeta}}\left(v_{\zeta}+\frac{t_{v}D_{\zeta}}{2}\frac{\partial}{\partial
v_{\zeta}}\right)P=t_{v}\left(\frac{\partial}{\partial
t}+v_{\zeta}\frac{\partial}{\partial\zeta}\right)P.$ (2.28)

Here, the terms on the right hand side (RHS) are much smaller than those on
the left hand side (LHS). Indeed if we assume that the fields are given by
their typical values $\zeta\sim\sqrt{\sigma_{\zeta}^{2}}$ then we have that
the term containing the time derivative is of order
$t_{v}\dot{\sigma}^{2}_{\zeta}P\sim\Delta_{\zeta}^{2}P$ while the second term
is of order $1/\sqrt{Ht}P$. On the other hand the terms on the LHS are both of
order $\mathcal{O}(1)\times P$. We can make use of this hierarchy if we expand
the PDF in powers of $Ht_{v}$:

$\displaystyle
P(\zeta,v_{\zeta},t)=P^{(0)}+Ht_{v}P^{(1)}+(Ht_{v})^{2}P^{(2)}+\dots.$ (2.29)

Then, by replacing this expression back into (2.28) we find that at leading
order in $t_{v}$ the equation for the first term $P^{(0)}$ is simply given by

$\displaystyle\frac{\partial}{\partial
v_{\zeta}}\left(v_{\zeta}+\frac{t_{v}2D_{\zeta}}{2}\frac{\partial}{\partial
v_{\zeta}}\right)P^{(0)}=0,$ (2.30)

whose solution can be written as

$\displaystyle
P^{(0)}(\zeta,v_{\zeta};t)=e^{-\frac{6H}{2D_{\zeta}}v_{\zeta}^{2}}\phi_{0}(\zeta,t),$
(2.31)

where $\phi_{0}$ is a function of $\zeta$ and $t$ that can be determined by
considering the equation for $P^{(1)}$:

$\displaystyle\frac{\partial}{\partial
v_{\zeta}}\left(v_{\zeta}+\frac{t_{v}D_{\zeta}}{2}\frac{\partial}{\partial
v_{\zeta}}\right)HP^{(1)}=\left(\frac{\partial\phi_{0}}{\partial
t}+v_{\zeta}\frac{\partial\phi_{0}}{\partial\zeta}\right)e^{-\frac{6H}{2D_{\zeta}}v_{\zeta}^{2}}.$
(2.32)

Now, notice that in this last equation the LHS corresponds to a total
derivative of $v_{\zeta}$, whereas the RHS is proportional to Gaussian
function of $v_{\zeta}$. Hence we can integrate it with respect to $v_{\zeta}$
and obtain the constraint equation $\frac{\partial\phi_{0}}{\partial t}=0$.
This in turns allows us to write a solution for $P_{1}$ given by

$\displaystyle
P^{(1)}(\zeta,v_{\zeta};t)=-\frac{v_{\zeta}}{H}\frac{\partial\phi_{0}}{\partial\zeta}e^{-\frac{6H}{D_{\zeta}}v_{\zeta}^{2}}+\frac{\phi_{1}(\zeta,t)}{H}e^{-\frac{6H}{2D_{\zeta}}v_{\zeta}^{2}},$
(2.33)

which now depends on another function $\phi_{1}(\zeta,t)$. Repeating the same
step we can find, at the next order in $t_{v}$, the following equation for
$P^{(2)}$:

$\displaystyle\frac{\partial}{\partial
v_{\zeta}}\left(v_{\zeta}+\frac{t_{v}D_{\zeta}}{2}\frac{\partial}{\partial
v_{\zeta}}\right)H^{2}P^{(2)}(\zeta,v_{\zeta},\psi;t)=\left(\frac{\partial\phi_{1}}{\partial
t}+v_{\zeta}\frac{\partial\phi_{1}}{\partial\zeta}+v_{\zeta}^{2}\frac{\partial^{2}\phi_{0}}{\partial_{\zeta}^{2}}\right)e^{-\frac{6H}{2D_{\zeta}}v_{\zeta}^{2}}.$
(2.34)

Given that this equation has the same structure as (2.32), we immediately
infer, after integrating over $v_{\zeta}$, that

$\displaystyle\frac{\partial\phi_{1}}{\partial
t}=\frac{D_{\zeta}}{6H}\frac{\partial^{2}\phi_{0}}{\partial_{\zeta}^{2}}.$
(2.35)

Collecting the terms for the PDF we obtain,

$\displaystyle
P(\zeta,v_{\zeta},t)=\Big{[}\phi_{0}(\zeta)+t_{v}v_{\zeta}\frac{\partial\phi_{0}}{\partial\zeta}+t_{v}\phi_{1}(\zeta,t)\Big{]}e^{-\frac{6H}{2D_{\zeta}}v_{\zeta}^{2}}.$
(2.36)

After integrating over $v_{\zeta}$, this result reduces to

$\displaystyle P(\zeta;t)=\phi_{0}(\zeta)+t_{v}\phi_{1}(\zeta,t).$ (2.37)

Finally, using the constraint equations for $\phi_{0}$ and $\phi_{1}$ we find,
up to first order in $t_{v}$, that $P(\zeta;t)$ must satisfy

$\displaystyle\frac{\partial}{\partial
t}P(\zeta,t)=\frac{D_{\zeta}}{18H^{2}}\frac{\partial^{2}P(\zeta,t)}{\partial\zeta^{2}},$
(2.38)

which is the Fokker-Planck equation (2.27) obtained by ignoring the
$\ddot{\zeta}$ term in the Langevin equation. Going beyond second order in
$t_{v}$ does not add any further correction as the equations obtained for the
other terms in the PDF expansion are the same as those in (2.34).

###  2.3 Spectator fields in de Sitter

To complement the previous discussions, we now study the statistics of a light
spectator field on a de Sitter background. The light scalar field $\psi$ has
potential $V(\psi)$ and its equation of motion is given by

$\displaystyle\ddot{\psi}+3H\dot{\psi}+\frac{k^{2}}{a^{2}}\psi+V^{\prime}(\psi)=0.$
(2.39)

We take $V^{\prime\prime}(\psi)\ll H^{2}$ since the field is light. As we did
with $\zeta$, the statistical properties of $\psi$ can be studied by dividing
the field into long and short wavelength modes. In this way, the long
wavelength field $\psi$, satisfies the Langevin equation

$\displaystyle\dot{\psi}=-\frac{V^{\prime}(\psi)}{3H}+\eta_{\psi}(t),$ (2.40)

where $\eta_{\psi}$ is a Gaussian noise representing the effects of the short
wavelength modes. The correlation function of the noise term is given by

$\displaystyle\langle\eta_{\psi}(t)\eta_{\psi}(t^{\prime})\rangle=\frac{H^{3}}{4\pi^{2}}\delta(t-t^{\prime}).$
(2.41)

Notice that we have chosen to disregard the role of $\ddot{\psi}$ in the
Langevin equation, which can be justified with the same arguments given in
Section 2.2. From the Langevin equation (2.40) it is possible to compute the
one point probability distribution function $P(\psi,t)$ by writing the
associated Fokker-Planck equation:

$\displaystyle\frac{\partial P}{\partial
t}=\frac{1}{3H}\frac{\partial}{\partial\psi}\left(V^{\prime}(\psi)P\right)+\frac{H^{3}}{8\pi^{2}}\frac{\partial^{2}P}{\partial\psi^{2}}.$
(2.42)

This equation is highly non linear since the drift
$-\frac{V^{\prime}(\psi)}{3H}$ is an arbitrary function of $\psi$.
Nevertheless it is possible to find an exact solution. This is done by first
noticing that (2.42) has an equilibrium solution

$\displaystyle\lim_{t\to\infty}P(\psi,t)=\exp\left(-\frac{8\pi^{2}V(\psi)}{3H^{4}}\right),$
(2.43)

which is obtained by imposing that $P$ is time independent. To obtain
solutions to (2.42) we can now write $P(\psi,t)$ as

$\displaystyle
P(\psi,t)=\exp\left(-\frac{4\pi^{2}V(\psi)}{3H^{4}}\right)\sum_{n=0}^{\infty}a_{n}\Phi_{n}(\psi)e^{-\Lambda_{n}(t-t_{0})},$
(2.44)

where the coefficients $\Lambda_{n}$ and the functions $\Phi_{n}$ satisfy the
following eigenvalue problem

$\displaystyle\left(-\frac{1}{2}\frac{\partial^{2}}{\partial\psi^{2}}+\frac{1}{2}(v^{\prime}(\psi)^{2}-v^{\prime\prime}(\psi))\right)\Phi_{n}(\psi)=\frac{4\pi^{2}\Lambda_{n}}{H^{3}}\Phi_{n}(\psi),\qquad
v(\phi)\equiv\frac{4\pi^{2}}{3H^{4}}V(\psi).$ (2.45)

The time dependence of $P(\psi,t)$ is controlled by the eigenvalues
$\Lambda_{n}$, which are positive and, for general potentials $V(\psi)$, their
value increase with $n$ (with $\Lambda_{0}=0$). This implies that the decay
rate to the equilibrium distribution is given by $1/\Lambda_{1}$. For
instance, when the potential is quadratic ($V(\psi)=\frac{1}{2}m^{2}\psi^{2}$)
one finds the solution of (2.45) is given by Hermite polynomials with
eigenvalues given by $\Lambda_{n}=\frac{m^{2}}{3H^{2}}\times n$. In this case,
the solution reaches equilibrium for $\Delta N\gg H^{2}/m^{2}\gg 1$.

Using the decomposition (2.44) it is also possible to deduce the statistical
properties of the equal time correlation function $G(R)$, where $R$ is the
distance between two points, found as

$G(R)=N\sum_{n}|A_{n}|^{2}e^{-2\log(RH)\Lambda_{n}/H},$ (2.46)

where the coefficients $A_{n}$ are given by

$A_{n}=N^{-1}\int d\psi\psi
e^{-\frac{4\pi^{2}}{3H^{4}}V(\psi)}\Phi_{n}(\psi).$ (2.47)

From this result, we may conclude that the stochastic approach is valid for a
patch of size $R\sim H^{-1}e^{H/\Lambda_{1}}$. For a quadratic potential this
is of order $R\sim H^{-1}e^{H^{2}/m^{2}}\gg H^{-1}$, which implies that the
statistical average occurs over a large number of Hubble patches. This
quantity has to be compared with the correlation length of the observed
universe, given by $\sim H^{-1}e^{\Delta N}$. This implies that a field fits
inside the observed universe if $\Delta N<H^{2}/m^{2}$.

The present analysis assumed a fixed de Sitter background, but it can be
generalised to the case of quasi-de Sitter backgrounds, as required to study
slow-roll inflation. In this case there are added difficulties. One problem
involves the role of gauge transformations on Hamiltonian constraints
satisfied at the level of the Langevin equations. In the following discussion
we will avoid this issue by assuming that the graviton is decoupled from
scalar fields (i.e. we consider the decoupling limit, in which the mixing with
gravity is negligible for energies larger than $\sqrt{\epsilon}H$).
Furthermore we will assume that the time dependence of the couplings is
negligible over the time scales we will consider (typically an $e$-fold).
Within this regime, the dynamics reduces to study the action for the curvature
perturbation $\zeta$ coupled to an isocurvature field $\psi$ via derivative
couplings. In this way the problem is analogous to studying two coupled
spectator fields, evolving on de Sitter.

##  3 Statistics for two-field inflation

In this section we use the tools introduced in the previous section to derive
the probability density function describing the statistics of fluctuations in
multifield theories. For now, we shall restrict our treatment to the case of
two-field models, and focus on the case of theories with linear interactions.
In Section 4 we consider the role of non-linear interactions.

Our starting point is to consider the two-field action (background plus
perturbations) describing inflation:

$S=S_{\rm EH}-\frac{1}{2}\int
d^{4}x\sqrt{-g}\left[\gamma_{ab}(\phi)\partial\phi^{a}\partial\phi^{b}+V(\phi)\right].$
(3.1)

where $S_{\rm EH}$ is the Einstein-Hilbert action, $\gamma_{ab}(\phi)$ is a
sigma model metric describing the geometry of the scalar field target space
and $V(\phi)$ is the scalar potential driving inflation. Deviations from a
geodesic trajectory are parametrised by the angular velocity $\Omega$. The
action for the curvature field $\zeta$ and the isocurvature field $\psi$ is
obtained by decomposing the fields along tangent and normal directions to the
inflationary trajectory (see [47, 48, 56, 6] for a more detailed explanation).
The quadratic action is given by [49]

$S=\frac{1}{2}\int
d^{4}xa^{3}\left[f_{\zeta}^{2}\left(\dot{\zeta}-\frac{2\Omega}{f_{\zeta}}\psi\right)^{2}-\frac{f_{\zeta}}{a^{2}}(\nabla\zeta)^{2}+\dot{\psi}^{2}+\frac{1}{a^{2}}(\nabla\psi)^{2}+\mu^{2}\psi^{2}\right],$
(3.2)

where $\Omega$ is the coupling between the two fields, and $\mu$ is the so
called entropy mass of $\psi$. In addition, we use $f_{\zeta}^{2}\equiv
2{M_{\mathrm{Pl}}}^{2}|\dot{H}|/H^{2}=2\epsilon{M_{\mathrm{Pl}}}^{2}$. The
equations of motion resulting from the variation of the previous action are

$\displaystyle\ddot{\zeta}+3H\dot{\zeta}+\frac{k^{2}}{a^{2}}\zeta$
$\displaystyle=-\frac{2\Omega}{f_{\zeta}}(\dot{\psi}+3H\psi),$ (3.3)
$\displaystyle\ddot{\psi}+3H\dot{\psi}+\frac{k^{2}}{a^{2}}\psi+m^{2}\psi$
$\displaystyle=2\Omega f_{\zeta}\dot{\zeta},$ (3.4)

where222As discussed in Refs. [7, 10], $\mu$ is the physical mass that
identifies the rest energy of one of the quanta in the spectrum of the theory
on subhorizon scales, whereas $m$ is just a mass parameter entering the
equation of motion (3.4). The quantity $m^{2}$ can be large and negative
without affecting the stability of the system, as long as $\mu^{2}\geq 0$.
$m^{2}=\mu^{2}-4\Omega^{2}$. Notice that the coupling $\Omega$ mixes both
fields, but these equations can still be solved using perturbation theory if
we assume that the coupling satisfies $\Omega/H\ll 1$. In this case we have
that the two point functions for $\zeta$ and $\psi$ at horizon crossing, are
given by

$\displaystyle\Delta_{\zeta}^{2}=\frac{1}{8\pi^{2}{M_{\mathrm{Pl}}}^{2}}\frac{H^{2}}{\epsilon},\qquad\Delta_{\psi}^{2}=\frac{H^{2}}{4\pi^{2}}.$
(3.5)

Since $\psi$ is a light field, it can continue evolving after crossing the
horizon. When the coupling $\Omega\neq 0$, $\psi$ seeds the perturbations of
$\zeta$ which could lead to its correlation functions growing on superhorizon
scales [49]. As we will see all these effects can be properly incorporated by
studying the linearised Langevin equations.

In order to analyse the stochastic dynamics let us note that $\psi$ couples to
$v_{\zeta}$. In order to include these terms it is more convenient to use the
phase space formulation of the Langevin equations. As explained in Section 2
this is achieved by introducing the time derivatives of the fields in the
Langevin equations. For (3.4) these correspond to

$\displaystyle\dot{\zeta}$ $\displaystyle=v_{\zeta},$
$\displaystyle\dot{v}_{\zeta}$
$\displaystyle=-3Hv_{\zeta}-2\frac{\Omega}{f_{\zeta}}v_{\psi}-\frac{6\Omega
H}{f_{\zeta}}\psi+\eta_{\zeta},$ $\displaystyle\dot{\psi}$
$\displaystyle=v_{\psi},$ $\displaystyle\dot{v}_{\psi}$
$\displaystyle=-3Hv_{\psi}-{m^{2}}\psi+2\Omega
f_{\zeta}v_{\zeta}+\eta_{\psi},$ (3.6)

where $\eta_{\zeta}$ and $\eta_{\psi}$ are Gaussian noise terms with
correlation functions given by

$\displaystyle\langle\eta_{\zeta}(t)\eta_{\zeta}(t^{\prime})\rangle=H^{2}\Delta_{\zeta}^{2}\delta(t-t^{\prime}),\qquad\langle\eta_{\psi}(t)\eta_{\psi}(t^{\prime})\rangle=\frac{H^{3}}{4\pi^{2}}\delta(t-t^{\prime}).$
(3.7)

Notice that the Langevin equations are coupled, hence the PDF do not factorise
into $P\propto P_{\zeta}(\zeta,v_{\zeta})P(\psi,v_{\psi})$. This adds some
complications due to the fact that $\zeta$ does not reach equilibrium, and so
it is not possible to a priori make use of the decomposition (2.44) to find
the PDF.

Nevertheless, since the Langevin equations are linear it is possible to find
an exact Gaussian solution. In this case all we need to do is to compute the
covariance matrix, as explained in Appendix C. Since we have the drift and the
noise matrices given by

$\displaystyle A=\left(\begin{array}[]{cccc}0&0&1&0\\\ 0&0&0&1\\\
0&-\frac{6H\Omega}{f_{\zeta}}&-3H&-\frac{2\Omega}{f_{\zeta}}\\\
0&-m^{2}&2f_{\zeta}\Omega&-3H\\\ \end{array}\right),\qquad
D=\left(\begin{array}[]{cccc}0&0&0&0\\\ 0&0&0&0\\\ 0&0&D_{\zeta}&0\\\
0&0&0&\frac{9H^{5}}{4\pi^{2}}\\\ \end{array}\right),$ (3.16)

the covariance matrix is given by

$\displaystyle
C(t)=\int_{0}^{t}\exp({(t-t^{\prime})A})D\exp((t-t^{\prime})A^{t})dt^{\prime}\
,$ (3.17)

where we have assumed initial conditions given by
$P=\delta(\zeta)\delta(v_{\zeta})\delta(\psi)\delta(v_{\psi})$ and we have set
the initial time to zero. Notice that it is not necessary to write the
explicit PDF since the variances are given by

$\displaystyle\sigma^{2}_{\phi_{a}}=C_{\phi_{a}\phi_{a}},$ (3.18)

where $\phi_{a}$ is one of the fields and also the corresponding element on
the diagonal of $C$. Off diagonal elements of $C$ are cross correlations
between different fields. Before writing explicit expressions for $C(t)$ let
us notice that the coupled dynamics imply that there are several time scales
over which the field decays. A useful way of understanding this is by noticing
that the time dependence is encoded in the exponential of the drift matrix
$A$. Since for $\mu\neq 0$, $A$ is diagonalisable it can be written as
$A=UDU^{-1}$ with $D$ a diagonal matrix containing the eigenvalues
$\lambda_{i}$ of $A$. Using this decomposition, the exponential of the drift
matrix can be written as $\exp(At)=U\exp(Dt)U^{-1}$ which implies that the
time dependence will appear in terms containing $e^{\lambda_{i}t}$.

For the drift matrix $A$ given in (3.16) the eigenvalues $\lambda_{i}$ are,
$0$,$-3H$, $-3H(1+\mu^{2}/H^{2})$ and $-\mu^{2}/3H$. Since the integrand will
contain factors of $\exp t(\lambda_{i}+\lambda_{j})$ there are three main
cases. First if $\lambda_{i}=\lambda_{j}=0$ after integrating implies the
appearance of terms linear in $t$, which are due to the variance of $\zeta$
always growing with time. The others cases arise when one of the eigenvalues
is different from zero. When the sum of the eigenvalues is proportional to
$3H(1+\mu^{2}/H^{2})$ the term decays after a time $t\sim 1/(3H)$ analogous to
$v_{\zeta}$ in the single field case. After integrating this term will
generate two pieces, one is constant because it is evaluated at $t^{\prime}=t$
and the other one decays. The third case is when the sum of the eigenvalues is
proportional to $\mu^{2}/H$. These terms decay on a longer time scale dictated
by the isocurvature mass $t_{\psi}\sim H/\mu^{2}$. Notice, moreover, that the
decay depends on the isocurvature mass, not on $m$. As an example let us look
at the variance of $\zeta$. From (3.17), this quantity is found to be given by

$\displaystyle\sigma_{\zeta}^{2}=C_{\zeta\zeta}\approx
Ht\Delta_{\zeta}^{2}\left(1+\frac{36H^{2}\Omega^{2}}{\mu^{4}}\right)-54\Delta_{\zeta}^{2}\frac{H^{4}\Omega^{2}}{\mu^{6}}\left(3-4e^{-\frac{\mu^{2}}{3H}t}+e^{-\frac{2\mu^{2}}{3H}t}\right),$
(3.19)

where we have kept only leading order terms and we have assumed $Ht\gg 1$. As
anticipated, the time dependence appears both linearly in $t$ and in powers of
$\exp(-\mu^{2}t/H)$. The first dependence is due to the modes entering the
comoving horizon as we discussed in Section 1. Notice, however, that there is
a new piece which depends on the mass and the coupling of the extra field. The
second term is due to the superhorizon evolution of the light scalar field
which takes a longer time to decay. At very large times $t\gg t_{\psi}$ this
contribution is suppressed and the time dependence is as usual. However, at
intermediate times $1\ll Ht\leq H^{2}/\mu^{2}$ the time dependence from the
exponential can dominate. To see this, let us expand the variance in small
$\mu$

$\displaystyle\sigma_{\zeta}^{2}=C_{\zeta\zeta}\approx\Delta_{\zeta}^{2}Ht\left(1+\frac{4}{3}(Ht)^{2}\frac{\Omega^{2}}{H^{2}}\right)+\Delta_{\zeta}^{2}\Omega^{2}/H^{2}\times\mathcal{O}((Ht)^{4}\mu^{2}/H^{2}).$
(3.20)

There are now cubic terms in $Ht$, which do not depend on the entropy mass but
only on the factor $\Omega/H$. This effect is due to the light field behaving
as a massless (ultralight) field before it settles into an equilibrium
distribution. This can be seen more directly if we look at the variance of
$\psi$ for $Ht\gg 1$

$\displaystyle\sigma_{\psi}^{2}=C_{\sigma\sigma}\approx\frac{3H^{4}}{8\pi^{2}\mu^{2}}\left(1-e^{-\frac{2\mu^{2}}{3H}t}\right).$
(3.21)

At late times the exponential can be neglected and one recovers the
equilibrium value. Interestingly at intermediate times $1\ll Ht\leq
H^{2}/\mu^{2}$ , expanding in small $\mu$ we find

$\displaystyle\sigma_{\psi}^{2}\approx\frac{H^{3}t}{4\pi^{2}}-\frac{H^{2}t^{2}}{12\pi^{2}}\mu^{2},$
(3.22)

which implies that for small $(Ht)\mu^{2}/H^{2}$ the field behaves as a
massless field. Later, when the second piece becomes of order
$(Ht)\mu^{2}/H^{2}\sim 1$ or larger this expansion stops being valid and we
need to consider the full expression. It is also interesting to rewrite the
variances in terms of the scale invariant power spectrum. This can be done by
writing time in terms of the number of efolds after horizon crossing. If a
mode with wavenumber $k$ leaves the horizon at a time $t_{k}$ we can then
write $H(t_{k}+t_{*})=\log(k/k_{*})$ where $k_{*}$ and $t_{*}$ are the longest
scale and time measured. Since the variance of a field $\phi$ is computed as

$\displaystyle\sigma_{\phi}^{2}=\int\frac{dk}{k}\Delta_{\phi}^{2},$ (3.23)

it is possible to invert this relation and write $\Delta_{\phi}$ as the
logarithmic derivative of the variance with respect to $k$. It is also
convenient to express all quantities in terms of e-folds using $\Delta
N=H(t_{k}+t_{*})$. Finally we find that

$\displaystyle\Delta_{\zeta}^{2}(t)$
$\displaystyle=\frac{H^{2}}{8\pi^{2}{M_{\mathrm{Pl}}}^{2}\epsilon}\left(1+36\frac{H^{2}\Omega^{2}}{\mu^{4}}\left(1-e^{-\frac{\mu^{2}}{3H^{2}}\Delta
N}\right)^{2}\right),$ (3.24) $\displaystyle\Delta_{\psi}^{2}(t)$
$\displaystyle=\frac{H^{2}}{4\pi^{2}}e^{-\frac{2\mu^{2}}{3H^{2}}\Delta N}.$
(3.25)

Notice that at $t=0$ we have $\Delta_{\zeta}^{2}(0)=\Delta_{\zeta}^{2}$, where
$\Delta_{\zeta}^{2}$ is defined in (3.5). That is, the value of
$\Delta_{\zeta}^{2}$ at horizon crossing acts as an initial condition for
$\Delta_{\zeta}^{2}(t)$, which continues to evolve since the light field has
not yet reached equilibrium. Notice that similar formulae were obtained by
using the in-in formalism in [52]. From the last formulae we also see that, in
general, light fields behave as massless fields during a time $1\ll\Delta N\ll
3H^{2}/\mu^{2}$. Expanding again in small $\mu$ we find that

$\displaystyle\Delta_{\zeta}^{2}(t)$
$\displaystyle=\frac{H^{2}}{8\pi^{2}{M_{\mathrm{Pl}}}^{2}\epsilon}\left(1+4(\Delta
N)^{2}\frac{\Omega^{2}}{H^{2}}\right).$ (3.26)

Notice that this is same superhorizon growth described in the case of
ultralight fields [49]; this is the case when the entropy mass is exactly
zero. We plot (3.24) and (3.25) in Fig.1. We see that for curvature field
there is an initial ultralight phase which will last as long as $\mu^{2}\Delta
N/3H^{2}\ll 1$. Notice that if this inequality holds until the end of
inflation then $H^{2}/\mu^{2}\ll\Delta N$, the ultralight phase is all there
is, since the field does not have time to start moving away from a massless
distribution. Another feature which we have previously discussed is that the
coupling changes the final amplitude of the primordial fluctuation. In the
case of the light field we see that the power spectrum eventually decays to
zero

Figure 1: We plot the scale invariant power spectra (3.24) amd (3.25). We set
$\Delta_{\zeta}^{0}=2\times 10^{-9}$. The parameters for the black line are
$\mu=0.2H$, $\Omega=.05$, the orange line $\mu=0.2H$, $\Omega=.2H$ and for the
gray line $\mu=.05H$, $\Omega=.2H$

###  3.1 Probability distribution function

Having understood the time dependence of the variances, now we would like to
derive the PDF of the two-field model. To start with, let us notice that the
time scales associated with the velocity fields $t_{v_{\psi}}$ and
$t_{v_{\zeta}}$ are much smaller than the scale of $\zeta$ and $\psi$. Using
this it is possible to integrate out $v_{\psi}$ and $v_{\zeta}$ from the
Fokker-Planck equation. We will explain in detail how this is done in Section
4.1.1. The resulting Fokker-Planck equation is

$\displaystyle\frac{\partial P}{\partial
t}=-\frac{2\Omega}{f_{\zeta}}\frac{\partial}{\partial\zeta}(\psi
P)+\frac{2\Delta_{\zeta}^{2}\Omega}{H}\frac{\partial^{2}}{\partial\zeta\partial\psi}(\psi
P)+\frac{H\Delta_{\zeta}^{2}}{2}\frac{\partial^{2}P}{\partial\zeta^{2}}+\frac{\partial}{\partial\psi}\left(t_{\psi}^{-1}\psi
P+\frac{D_{\psi}}{2}\frac{\partial P}{\partial\psi}\right),$ (3.27)

where

$\displaystyle D\psi\equiv\frac{H^{3}}{4\pi^{2}},\qquad
t_{\psi}\equiv\frac{3H}{\mu^{2}}.$ (3.28)

This Fokker-Planck equation is still linear, even though it has a mixed noise
term, and it can be solved with the techniques of Appendix C. Further
simplification is possible if we consider the following: As explained in
Appendix C, the noise term is computed by using the two-point function of the
field at horizon crossing. Nevertheless, we have seen that the value of the
variance grows with time on superhorizon scales. On the other hand, a direct
computation of the two point function also shows a superhorizon growth [49,
50]

$\displaystyle\Delta_{\zeta}^{2}(t)=\frac{H^{2}}{8\pi^{2}{M_{\mathrm{Pl}}}^{2}\epsilon}\left(1+36\frac{H^{2}\Omega^{2}}{\mu^{4}}(1-e^{-t/t_{\psi}})^{2}\right),$
(3.29)

which, of course, is the same result we obtained in (3.24). This can be
understood as follows. As it has been previously pointed out in Refs. [49,
50], in the long wavelength limit the equations of motion for $\zeta$ decouple
if they are written in terms of
$v_{\zeta}=\dot{\zeta}-\frac{2\Omega}{f_{\zeta}}\psi$

$\displaystyle\dot{v}_{\zeta}+3Hv_{\zeta}+\frac{k^{2}}{a^{2}}=0,$ (3.30)
$\displaystyle\ddot{\psi}+3H\dot{\psi}+\frac{k^{2}}{a^{2}}\psi+\mu^{2}\psi=2\Omega
f_{\zeta}v_{\zeta}.$ (3.31)

Notice that in the long wavelength limit the first equation admits as a
solution $v_{\zeta}=v_{\zeta}^{0}a^{-3}$. Plugging this solution back into the
second equation implies that the last source term vanishes for $t\gg 3H$,
hence both equations are decoupled. Nonetheless, the fact that $v_{\zeta}$
also depends on $\psi$ explains the superhorizon growth since the curvature
mode will depend on integrals of $\psi$ which does not immediately decay on
superhorizon scales.

We can make use of the fact that the equation of motion of $v_{\zeta}$ is free
if we define it as $v_{\zeta}=\dot{\tilde{\zeta}}$ instead of the usual
definition. In order to do this consistently we also require the noise of the
field $v_{\zeta}$ to include the superhorizon growth. This can be done by
shifting $D_{\zeta}$ to be

$\displaystyle D_{\zeta}\to 9H^{3}\Delta_{\zeta}^{2}(t)\equiv
9H^{4}\Delta_{\zeta}^{2}\left(1+36\frac{H^{2}\Omega^{2}}{\mu^{4}}(1-e^{-t/t\psi})^{2}\right),$
(3.32)

where $\Delta_{\zeta}^{2}(t)$ is the same as that given in (3.24). We defer a
more detailed analysis on how to modify the noise term to Appendix D. Since
the variable that appears in the equations of motion is $v_{\zeta}$ there is
no difference between which variable we use, which implies that the statistics
of $\tilde{\zeta}$ are the same than for $\zeta$. We can check explicitly by
solving the Fokker-Planck equation, which now becomes

$\displaystyle\frac{\partial P}{\partial
t}=+\frac{2\Delta_{\zeta}^{2}(t)\Omega}{H}\frac{\partial^{2}}{\partial\tilde{\zeta}\partial\psi}(\psi
P)+\frac{H\Delta_{\zeta}^{2}(t)}{2}\frac{\partial^{2}P}{\partial\tilde{\zeta}^{2}}+\frac{\partial}{\partial\psi}\left(t_{\psi}^{-1}\psi
P+\frac{D_{\psi}}{2}\frac{\partial P}{\partial\psi}\right).$ (3.33)

This equation can be solved using the techniques described in Appendix C. The
solution is found to be given by

$\displaystyle
P(\zeta,\psi,t)\sim\exp\left(-\frac{\psi^{2}}{\sigma_{\psi}^{2}}-\frac{1}{2\sigma_{\zeta}^{2}-2\kappa^{2}/\sigma_{\psi}^{2}}\left(\tilde{\zeta}-\frac{\kappa}{2\sigma_{\psi}^{2}}\psi\right)^{2}\right),$
(3.34)

where the variances are defined in (3.19) and (3.21), and

$\kappa\equiv\frac{2}{3}t_{\psi}\Delta_{\zeta}^{2}f_{\zeta}\Omega(1-e^{-t/t_{\psi}}),$
(3.35)

is the off-diagonal variance which grows from an initial value $\kappa\sim
Ht\Omega/f_{\zeta}$ in the ultralight phase to a value of $\kappa\sim
Ht_{\psi}\Omega/f_{\zeta}$ when the distribution for $\psi$ reaches
equilibrium. Notice that the combination $\kappa/\sigma_{\psi}^{2}$ is almost
constant and is related to the size of the coupling between the two fields.
Also we have that $\sigma_{\zeta}^{2}\sigma_{\psi}^{2}\gg\kappa^{2}$ and so
the denominator of the second term in the PDF is always positive.

After marginalising over $\psi$ we find that the variance of $\tilde{\zeta}$
is given by (3.19), as anticipated, hence the Gaussian statistics of this
field are independent on how we define $v_{\zeta}$. Clearly the same happens
for $\psi$. Moreover, since all couplings are linear, the variances do not
change but the minima of $\psi$ is displaced. By minimizing the PDF with
respect to $\psi$ we find that the minima $\bar{\psi}$ is at

$\displaystyle\bar{\psi}=\frac{\kappa}{2\sigma_{\zeta}^{2}}\zeta\sim\frac{\Omega}{H}f_{\zeta}\tilde{\zeta},$
(3.36)

which implies that the classical trajectory of $\psi$ is shifted by the
interaction. Of course, this does not mean that the statistical fluctuations
are modified, since they are still simply given by $\sigma^{2}_{\psi}$. From
now on, for simplicity, we will remove the tilde from $\zeta$.

##  4 Non linear interactions in a two-field model

So far we have focused our attention on the second order action (3.2) where it
was possible to find an exact solution of the Fokker-Planck equation. Now we
would like to consider the role of non linear terms and study how they modify
the probability distribution (3.34). As already stated, for a spectator field
there are well known techniques which allow us to find non-perturbative
solutions. However these techniques are not useful to uncover the PDF for the
curvature field $\zeta$ since, as we have discussed, there is no equilibrium
distribution for $\zeta$. Despite this shortcoming, we will be able to uncover
precise non-perturbative effects on the joint distribution $P(\zeta,\psi)$,
based on the Gaussian distribution we found in the previous section. To do so,
our strategy will be to ignore non linear terms in $v_{\zeta}$, while keeping
higher order terms in $\psi$.

Let us start this discussion by writing down the action for perturbations in
the case of a canonical two-field model of inflation [57, 58]. Up to leading
order in slow-roll the action is,

$\displaystyle S$ $\displaystyle=\frac{1}{2}\int
d^{4}xa^{3}\left\\{\left(f_{\zeta}+\frac{\Omega}{H}\psi\right)^{2}\left(\dot{\zeta}^{2}-\frac{(\nabla\zeta)^{2}}{a^{2}}\right)-2\frac{\Omega}{H}(2f_{\zeta}H+\Omega\psi)\psi\dot{\zeta}\right.$
$\displaystyle\hskip
65.44142pt\left.+\dot{\psi}^{2}-\frac{(\nabla\psi)^{2}}{a^{2}}-V(\psi)\right\\}+\dots,$
(4.1)

For simplicity we consider the first few powers of the potential
$V(\psi)=m^{2}\psi^{2}+\frac{\lambda}{3}\psi^{3}+\frac{g}{12}\psi^{4}+\cdots$.
A crucial point is that the action contains non linear interactions between
$\zeta$ and $\psi$, which are due to the non geodesic motion of the background
trajectory. Apart from those appearing explicitly in (4.1), there are no
further interactions between the two fields, which can be understood as
arising from the original canonical kinetic term in the action. Non canonical
kinetic terms will generate higher order interactions between $\dot{\zeta}$
and $\psi$ which we are assuming to be suppressed. Furthermore, notice that we
have not written interactions including gradients coming from gravitational
couplings, as they will be negligible for the stochastic dynamics. Finally, as
in the linear case, it will be convenient to write the action in terms of
$v_{\zeta}=\dot{\zeta}-\frac{2\Omega}{f_{\zeta}}\psi$. Doing this, the action
becomes

$\displaystyle S$ $\displaystyle=\frac{1}{2}\int
d^{4}xa^{3}\left\\{f_{\zeta}^{2}\left(\dot{\zeta}-\frac{2\Omega}{f_{\zeta}}\psi\right)^{2}-f_{\zeta}^{2}\frac{(\nabla\zeta)^{2}}{a^{2}}+\dot{\psi}^{2}-\frac{(\nabla\psi)^{2}}{a^{2}}-\mu^{2}\psi^{2}\right.$
$\displaystyle\hskip
62.59596pt\left.+\frac{6\Omega^{2}}{H}\psi^{2}\left(\dot{\zeta}-\frac{2\Omega}{f_{\zeta}}\psi\right)+\frac{2f_{\zeta}\Omega}{H}\psi\left(\dot{\zeta}-\frac{2\Omega}{f_{\zeta}}\psi\right)^{2}-\frac{\tilde{\lambda}}{3}\psi^{3}+\cdots\right\\}\
,$ (4.2)

where we have kept terms up to cubic order with respect to $\psi$. Higher
interactions are suppressed by further powers of $\Omega/H$ (which we take as
a small parameter) although mixed terms are only up to fourth order and they
can be reincorporated without trouble. Notice that in the same way as the mass
term of $\psi$ becomes the entropic mass $\mu$, other self interaction
couplings are also modified, $\tilde{g}=g-\frac{12\Omega^{3}}{f_{\zeta}H}$ and
$\tilde{\lambda}=\lambda-48\frac{\Omega^{4}}{H^{2}f_{\zeta}^{2}}$.

For simplicity, let us examine the case when there is a large cubic
interaction for $\psi$ but the equation for $\zeta$ can be considered as free.
The equations of motion are

$\displaystyle\frac{d}{dt}v_{\zeta}+3Hv_{\zeta}+\frac{k^{2}}{a^{2}}\zeta$
$\displaystyle=-6\frac{\Omega}{f_{\zeta}}v_{\zeta}\psi-\frac{6\Omega^{2}}{Hf_{\zeta}^{2}}(3H\psi^{2}+2\dot{\psi}\psi)\
,$ (4.3)
$\displaystyle\ddot{\psi}+3H\dot{\psi}+\frac{k^{2}}{a^{2}}\psi+\mu^{2}\psi$
$\displaystyle=-2\Omega f_{\zeta}v_{\zeta}+\frac{2\Omega^{2}}{H}\psi
v_{\zeta}+\frac{\Omega f_{\zeta}}{H}v_{\zeta}^{2}\ ,$ (4.4)

where we have disregarded higher order self interactions of $\psi$ (which can
be included back at any point of our analysis). In order to apply the
stochastic approximation for $\zeta$ we have to demand that higher order
interactions are suppressed. The first term on the RHS of Eq. 4.4 is
suppressed for typical fluctuations, since $\Omega\ll f_{\zeta}$. We will
assume that the third term is at most of the size of the second one. For the
second term on the RHS we have that

$\displaystyle\frac{6\Omega^{2}}{f_{\zeta}^{2}}\frac{\psi^{2}}{Hv_{\zeta}}\ll
1\ ,$ (4.5)

which follows from the fact that $\Omega\ll f_{\zeta}$ and that the variance
of $\zeta$ is much larger than the one for $\psi$ for $\Omega\neq 0$. For
larger fluctuations of $\zeta$ the inequality (4.5) still holds. Since at
leading order Eq. (4.3) is free then $v_{\zeta}$ decays after leaving the
horizon and the LHS of Eq. (4.4) is negligible.

###  4.1 Fokker-Planck equation

In order to study the effect of $v_{\zeta}$ over $\psi$ more systematically we
will analyse the stochastic dynamics of the two-field system. For this, we use
the strategy employed for the linear case, that is, we coarse grain the fields
directly from the equation of motion for the perturbations (4.4). Leading non
linearities come only from long-wavelength modes, with interactions involving
short-wavelength modes being subdominant. In the end, the effects of short
wavelength modes reduce to the same linear noise terms as in the linear
Langevin equations [30]. Moreover, we will consider the couplings to be small
with respect to $H$ so it is possible to treat interactions using perturbation
theory. This implies that the noise terms are as in the linear case considered
in (3.6). Another simplification comes from the fact that there are no
interactions involving $\dot{\psi}$ (as they are gravitationally suppressed).
Because of this, we can neglect all terms with time derivatives of $\psi$
except for the leading friction term. Indeed this is related to the fact that
for typical fluctuations $\ddot{\psi}\ll H^{2}\psi$, since in the long
wavelength limit we have that

$\displaystyle\frac{\dot{\psi}}{H\psi}\sim\frac{2\mu^{2}}{3H^{2}}\ll 1\ ,$
(4.6)

where the last inequality follows from the fact that we are considering light
fields. By the same argument we may ignore the second derivative of $\psi$ in
the first equation. Of course this can be understood as integrating out
$\dot{\psi}$ from the Fokker-Planck equation and the details will be analogous
to those examined in the case of the curvature field in Section 2.2. Finally,
after separating the equations into long- and short-wavelength modes, we find
that the Langevin equations for the long-wavelength mode are

$\displaystyle\frac{dv_{\zeta}}{dt}+3Hv_{\zeta}+\frac{6\Omega}{f_{\zeta}}v_{\zeta}\psi+\frac{18\Omega^{2}}{f_{\zeta}^{2}}\psi^{2}$
$\displaystyle=3H\eta_{\zeta},$ (4.7) $\displaystyle
3H\dot{\psi}+\mu^{2}\psi+2\Omega f_{\zeta}v_{\zeta}-\frac{2\Omega^{2}}{H}\psi
v_{\zeta}-\frac{\Omega f_{\zeta}}{H}v_{\zeta}^{2}$
$\displaystyle=\eta_{\psi}.$ (4.8)

As previously discussed, there are two ways of introducing $\zeta$ to the
Langevin equations. We follow the simpler one, whereby we consider an extra
Langevin equation for the field $\dot{\zeta}=v_{\zeta}$. As we described in
Section 3.1, this means that we need to include a time dependent noise for
$v_{\zeta}$. After considering these steps, we finally find that the
associated Fokker-Planck equation is given by

$\displaystyle\frac{\partial P}{\partial t}=$
$\displaystyle-\frac{\partial}{\partial\zeta}(v_{\zeta}P)+\frac{\partial}{\partial
v_{\zeta}}\left(\left(3Hv_{\zeta}+\frac{18\Omega^{2}}{f_{\zeta}^{2}}\psi^{2}+\frac{6\Omega}{f_{\zeta}}v_{\zeta}\psi\right)P\right)$
$\displaystyle+\frac{\partial}{\partial\psi}\left(\left(\frac{\mu^{2}}{3H}\psi+\frac{2\Omega
f_{\zeta}}{3H}v_{\zeta}-\frac{2\Omega^{2}}{3H^{2}}\psi v_{\zeta}-\frac{\Omega
f_{\zeta}}{3H^{2}}v_{\zeta}^{2}\right)P\right)+\frac{D_{\psi}}{2}\frac{\partial^{2}}{\partial\psi^{2}}P+\frac{9}{2}H^{3}\Delta_{\zeta}^{2}(t)\frac{\partial^{2}P}{\partial
v_{\zeta}^{2}},$ (4.9)

where $\Delta_{\zeta}^{2}(t)$ is the same quantity found in Eq. (3.24),
$D_{\psi}=H^{3}/4\pi^{2}$, and where the variances are given by

$\displaystyle\sigma_{\zeta}^{2}(t)$
$\displaystyle=\Delta_{\zeta}^{2}Ht\left(1+\frac{36H^{2}\Omega^{2}}{\mu^{2}}\right)-54\Delta_{\zeta}^{2}\frac{H^{4}\Omega^{2}}{\mu^{6}}\left(3-4e^{-t/t_{\psi}}+e^{-2t/t_{\psi}}\right),$
(4.10) $\displaystyle\sigma_{\psi}^{2}(t)$
$\displaystyle=\frac{3H^{4}}{8\pi^{2}\mu^{2}}\left(1-e^{-2t/t_{\psi}}\right).$
(4.11)

Before continuing, let us comment on the fact that drift terms including
powers of $v_{\zeta}$ will become $\zeta$ derivatives. This can be understood
as a consequence of the shift symmetry of the curvature mode (see also [59]).
Using this we can deduce that terms including two derivatives of $\zeta$ will
change the variance of the $\zeta$ distribution, and the tail of the
distribution of $\psi$. Since we are interested in the tail of the
distribution of $\zeta$ we can ignore them for now and include them later.
This is achieved by imposing that the quadratic terms in the drift for $\psi$
is larger than the quadratic terms in the drift in $v_{\zeta}$ or,
equivalently, that

$\displaystyle\frac{\Omega}{H}\frac{f_{\zeta}\zeta}{\psi}\gg 1\ ,$ (4.12)

where we have used the fact that $v_{\zeta}\sim H\zeta$. Notice that this is
achieved only for $\Omega$ relatively large, although not necessarily larger
than $H$. A non zero $\Omega$ increases the variance of $\zeta$ making it much
larger than that of $\psi$, which otherwise will be very similar. Besides
that, if $\Omega/H$ is very suppressed, it will make the inequality in (4.12)
impracticable. In what follows we will assume that (4.12) holds, and we will
comment on its effect on the PDF later. Similarly as we saw in Section 2.2,
terms in the drift containing powers of $v_{\zeta}$ transform into derivatives
of $\zeta$, $v_{\zeta}\to\Delta_{\zeta}^{2}\frac{\partial}{\partial_{\zeta}}$,
which we will show is due to the shift symmetry of $\zeta$. Using this
property, we can neglect the term proportional to $v_{\zeta}^{2}$ in the drift
of $\psi$ since is subleading with respect to the term proportional to $\psi
v_{\zeta}$.

####  4.1.1 Adiabatic elimination of $v_{\zeta}$

We have used the variable $v_{\zeta}$ since it was a convenient way of
studying the derivative couplings that appear in the action. Nevertheless, we
are not interested in the statistical properties of $v_{\zeta}$ and moreover
it decays faster than the other fields. It is then useful to eliminate
$v_{\zeta}$ from the Fokker-Planck equation and obtain $P(\zeta,\sigma)$
directly. The way of doing this systematically is called adiabatic elimination
of fast variables [55, 54]. We made a similar computation in Section 2.2,
which we now generalise to include the coupling with another field. The idea
is to expand the probability density function in powers of the time scale of
the fast variable. Then replacing order by order it will be possible to
factorise the dependence on the fast variable from the slow variables. For our
case, it is convenient to write the Fokker-Planck equation as

$\displaystyle P(\zeta,v_{\zeta},\psi;t)$ $\displaystyle=$ $\displaystyle
P^{(0)}(\zeta,v_{\zeta},\psi;t)+(Ht_{v})P^{(1)}(\zeta,v_{\zeta},\psi;t)$
(4.13) $\displaystyle+(Ht_{v})^{2}P^{(2)}(\zeta,v_{\zeta},\psi;t)+\cdots.$

This expansion becomes useful if we write the Fokker-Planck equation as

$\displaystyle\frac{d}{dv_{\zeta}}\left(v_{\zeta}+\frac{3}{2}H^{2}\Delta^{2}_{\zeta}(t)\frac{d}{dv_{\zeta}}\right)P=\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad$
$\displaystyle
t_{v}\left(\frac{d}{dt}+v_{\zeta}\frac{d}{d\zeta}-\frac{\partial}{\partial\psi}(t_{\psi}^{-1}\psi+\frac{2\Omega
f_{\zeta}}{3H}v_{\zeta}-\frac{2\Omega^{2}}{3H^{2}}\psi
v_{\zeta})+\frac{D_{\psi}}{2}\frac{\partial^{2}}{\partial\psi^{2}}\right)P\ ,$
(4.14)

where we have written the mass in terms of the time scale
$t_{\psi}^{-1}=\mu^{2}/3H$. Notice that the operator on the RHS of Eq. (4.14)
is of order $t_{v}/t_{\psi}\ll 1$ whereas the operator on the LHS is of order
one, which is the reason why the approximation is well justified. After
replacing (4.13) into the Fokker-Planck equation this can also be written as a
series in powers of $t_{v}$. At zeroth order in $t_{v}$ we find the following
equation

$\displaystyle\frac{\partial}{\partial
v_{\zeta}}(v_{\zeta}P^{(0)})+\frac{3}{2}H^{2}\Delta^{2}_{\zeta}(t)\frac{\partial^{2}P^{(0)}}{\partial
v_{\zeta}^{2}}=0.$ (4.15)

Notice that this equation only specifies the $v_{\zeta}$ dependence of
$P^{(0)}$. Since this is a first order ODE we can solve it and write $P^{(0)}$
as

$\displaystyle
P^{(0)}(\zeta,v_{\zeta},\psi;t)=e^{-\frac{v_{\zeta}^{2}}{3H^{2}\Delta_{\zeta}^{2}(t)}}\phi_{0}(\zeta,\psi;t)\
,$ (4.16)

which is a Gaussian distribution of $v_{\zeta}$ with variance given by
$(3/2)H^{2}\Delta_{\zeta}^{2}(t)$ in line with the linear solution for $P$. Of
course, had we included higher order terms in $v_{\zeta}$, they would have
entered into the LHS of (4.14), and the solution (4.16) would have had to be
modified accordingly. To figure out what the restrictions are for $P^{(1)}$ we
plug (4.16) into the expansion for $P$. We find that at order $t_{v}$

$\displaystyle\frac{\partial}{\partial
v_{\zeta}}\left(v_{\zeta}+\frac{3}{2}H^{2}\Delta^{2}_{\zeta}(t)\frac{\partial}{\partial
v_{\zeta}}\right)HP^{(1)}=\left(\frac{\partial}{\partial\psi}(t_{\psi}^{-1}\psi\phi_{0})+\frac{D_{\psi}}{2}\frac{\partial^{2}\phi_{0}}{\partial\psi^{2}}-\frac{\partial\phi_{0}}{\partial
t}+\right.$ $\displaystyle\qquad\left.+\left(\frac{2\Omega
f_{\zeta}}{3H}\frac{\partial\phi_{0}}{\partial\psi}-\frac{2\Omega^{2}}{3H^{2}}\frac{\partial}{\partial\psi}(\psi\phi_{0})-\frac{\partial\phi_{0}}{\partial\zeta}\right)v-\frac{\dot{\Delta}_{\zeta}^{2}(t)}{3H^{2}\Delta_{\zeta}^{4}(t)}\phi_{0}v_{\zeta}^{2}\right)e^{-\frac{v_{\zeta}^{2}}{3H^{2}\Delta_{\zeta}^{2}(t)}}\
.$ (4.17)

Which has a similar structure to (4.15) but with a more complicated RHS. Due
to this, it is not possible to immediately solve for $P^{(1)}$. Yet, it is
possible to simplify the last equation by noticing that the LHS is a total
derivative and that the RHS is multiplied by a Gaussian function. After
integrating over $v_{\zeta}$ the LHS vanishes while on the RHS only some terms
with even powers of $v_{\zeta}$ remain. Since these have to add up to zero, we
finally find

$\displaystyle\frac{\partial}{\partial\psi}(t_{\psi}^{-1}\psi\phi_{0})+\frac{D_{\psi}}{2}\frac{\partial^{2}\phi_{0}}{\partial\psi^{2}}-\frac{\partial\phi_{0}}{\partial
t}-\frac{\dot{\Delta}_{\zeta}^{2}(t)}{2\Delta_{\zeta}^{2}(t)}\phi_{0}=0.$
(4.18)

Notice that the first three terms form a Fokker-Planck equation for
$\phi_{0}(\psi,t)$ and hence are of order $1/t_{\psi}$. The last term is of
order $t_{v}^{-1}$ so we can neglect it. We can in principle solve for the
$\psi$ dependence of $\phi_{0}$, which is just given by a Gaussian with
variance $\sigma_{\psi}^{2}$, although the explicit solution will not be
important. Instead, we can now find the solution for $P_{1}$ by using the
above condition. Equation (4.17) then reduces to

$\displaystyle\frac{\partial}{\partial
v_{\zeta}}\left(v_{\zeta}+\frac{3}{2}H^{2}\Delta^{2}_{\zeta}(t)\frac{\partial}{\partial
v_{\zeta}}\right)HP^{(1)}=-\left(-\frac{2\Omega
f_{\zeta}}{3H}\frac{\partial\phi_{0}}{\partial\psi}+\frac{2\Omega^{2}}{3H^{2}}\frac{\partial}{\partial\psi}(\psi\phi_{0})+\frac{\partial\phi_{0}}{\partial\zeta}\right)v_{\zeta}e^{-\frac{v_{\zeta}^{2}}{3H^{2}\Delta_{\zeta}^{2}(t)}},$
(4.19)

whose solution can be written as

$\displaystyle P^{(1)}=-\frac{v_{\zeta}}{H}\left(-\frac{2\Omega
f_{\zeta}}{3H}\frac{\partial\phi_{0}}{\partial\psi}+\frac{2\Omega^{2}}{3H^{2}}\frac{\partial}{\partial_{\psi}}(\psi\phi_{0})+\frac{\partial\phi_{0}}{\partial\zeta}\right)e^{-\frac{v_{\zeta}^{2}}{3H^{2}\Delta_{\zeta}^{2}(t)}}+\frac{\phi_{1}(\zeta,\psi,t)}{H}e^{-\frac{v_{\zeta}^{2}}{3H^{2}\Delta_{\zeta}^{2}(t)}}\
.$ (4.20)

If we now plug this solution into the terms of order $t_{v}^{2}$, we find that

$\displaystyle\frac{\partial}{\partial
v_{\zeta}}\left(v_{\zeta}+\frac{3}{2}H^{2}\Delta^{2}_{\zeta}(t)\frac{\partial}{\partial
v_{\zeta}}\right)H^{2}P^{(2)}=\Bigg{[}\frac{\partial}{\partial\psi}(t_{\psi}^{-1}\psi\phi_{1})+\frac{D_{\psi}}{2}\frac{\partial^{2}\phi_{1}}{\partial\psi^{2}}$
$\displaystyle\qquad-\frac{\partial\phi_{1}}{\partial
t}+(\dots)v_{\zeta}+(\dots)v_{\zeta}^{3}+\bigg{(}-\frac{4\Omega
f_{\zeta}}{3H}\frac{\partial^{2}\phi_{0}}{\partial\psi\partial\zeta}+\frac{4\Omega^{2}}{3H^{2}}\frac{\partial^{2}}{\partial\psi\partial\zeta}(\psi\phi_{0})$
$\displaystyle\qquad+\frac{\partial^{2}\phi_{0}}{\partial\zeta^{2}}-\frac{4\Omega^{4}}{H^{4}}(\phi_{0}+\psi(3\frac{\partial\phi_{0}}{\partial\psi}+\psi{\frac{\partial\phi_{0}}{\partial\psi}})\bigg{)}v_{\zeta}^{2}\Bigg{]}e^{-\frac{v_{\zeta}^{2}}{3H^{2}\Delta_{\zeta}^{2}(t)}}\
.$ (4.21)

This equation is similar to the constraint for $P^{(1)}$, thus we can again
integrate over $v_{\zeta}$. We find that

$\displaystyle\frac{\partial}{\partial\psi}(t_{\psi}^{-1}\psi\phi_{1})+\frac{D_{\psi}}{2}\frac{\partial^{2}\phi_{1}}{\partial\psi^{2}}-\frac{\partial\phi_{1}}{\partial
t}-2\Omega
f_{\zeta}\Delta_{\zeta}^{2}\frac{\partial^{2}\phi_{0}}{\partial\psi\partial\zeta}$
$\displaystyle\qquad\qquad+2\Omega^{2}\Delta_{\zeta}^{2}(t)\frac{\partial^{2}}{\partial\psi\partial\zeta}(\psi\phi_{0})+\frac{3H^{2}\Delta_{\zeta}^{2}(t)}{2}\frac{\partial^{2}\phi_{0}}{\partial\zeta^{2}}=0,$
(4.22)

where we have neglected sub leading terms. Notice that the terms proportional
to $\Omega^{4}$ can be recast as subleading corrections to the noise of
$\zeta$ of the form $\Omega^{2}\zeta^{2}$. These terms have been recently
studied [38, 60, 59] and we will leave their analysis to future work. Instead
of finding a solution for $P_{2}$, let us collect the terms in the PDF up to
order $t_{v}$

$\displaystyle
P(\zeta,v_{\zeta},\psi;t)=\left(\phi_{0}(\zeta,\psi;t)+t_{v}\left(v_{\zeta}\left(\frac{2\Omega^{2}}{H^{2}}\frac{\partial}{\partial_{\psi}}(\psi\phi_{0})+\frac{\partial\phi_{0}}{\partial\zeta}\right)+\phi_{1}(\zeta,\psi;t)\right)\right)e^{-\frac{v_{\zeta}^{2}}{3H^{2}\Delta_{\zeta}^{2}(t)}}\
.$ (4.23)

Since we are not interested in the full distribution, but only on the one on
configuration space, we can integrate out $v_{\zeta}$, which leads to
$P(\zeta,\psi,t)=\sqrt{3\pi
H^{2}\Delta_{\zeta}^{2}(t)}(\phi_{0}+t_{v}\phi_{1})$. Finally collecting all
the ingredients, taking the time derivative of $P(\zeta,\psi,t)$ and using
(4.18) and (4.23) we find that

$\displaystyle\boxed{\frac{\partial P}{\partial
t}=\frac{\partial}{\partial\psi}\left(\frac{\psi
P}{t_{\psi}}+\frac{D_{\psi}}{2}\frac{\partial
P}{\partial\psi}\right)+\frac{2\Omega}{3}\Delta_{\zeta}^{2}(t)\frac{\partial^{2}}{\partial\psi\partial\zeta}\left(\frac{\Omega}{H}\psi
P-f_{\zeta}P\right)+\frac{H\Delta_{\zeta}^{2}(t)}{2}\frac{\partial^{2}P}{\partial^{2}\zeta}+\mathcal{O}(t_{v}),}$
(4.24)

which is a Fokker-Planck equation for $\zeta$ and $\psi$ after integrating out
$v_{\zeta}$. Notice that the couplings between $v_{\zeta}$ and $\psi$ now
translate into mixed derivatives between $\psi$ and $\zeta$. Of course in the
absence of such couplings, the system reduces to two linear Fokker-Planck
equations. We can find an analytical solution if we Fourier transform $\zeta$
and $\psi$ to $p$ and $q$. The resulting equation can be solved by looking for
solutions of the form

$\displaystyle
P(p,k;t)=e^{-\frac{\sigma_{\zeta}^{2}(t)}{2}p^{2}}e^{-M(p,t)k^{2}-L(p,t)k},$
(4.25)

which translates into two independent ODEs for $M$ and $L$:

$\displaystyle\frac{dM(t)}{dt}+\frac{2}{3H^{2}t_{\psi}}(3H^{2}+2ipt_{\psi}\Omega^{2}H\Delta_{\zeta}^{2}(t))M(t)-\frac{D_{\psi}}{2}=0,$
(4.26)
$\displaystyle\frac{dL(t)}{dt}+\frac{1}{3H^{2}t_{\psi}}(3H^{2}+2ipt_{\psi}\Omega^{2}H\Delta_{\zeta}^{2}(t))L(t)-\frac{2f_{\zeta}\Omega\Delta_{\zeta}^{2}(t)}{3}p=0.$
(4.27)

Imposing that at $t=0$, $M=0$, we find

$\displaystyle
M(t)=\frac{1}{2}D_{\psi}e^{-\frac{2t}{t_{\psi}}-\frac{4ip\Omega^{2}\sigma_{\zeta}^{2}(t)}{3H^{2}}}\int_{0}^{t}dt^{\prime}e^{\frac{2t}{t_{\psi}}+\frac{4ip\Omega^{2}\sigma_{\zeta}^{2}(t)}{3H^{2}}}\
.$ (4.28)

We can approximate the integral using a saddle point approximation to find

$\displaystyle M(t)$
$\displaystyle=\frac{3D_{\psi}t_{\psi}}{4(3H+2ipt_{\psi}\Omega^{2}\Delta_{\zeta}^{2}(0))}e^{-\frac{2t}{t_{\psi}}-\frac{4ip\Omega^{2}\sigma_{\zeta}^{2}(t)}{3H^{2}}}\left(-1+e^{\frac{2t}{t_{\psi}}+\frac{4ipt\Omega^{2}\Delta_{\zeta}^{2}(0)}{3H}}\right)\
.$ (4.29)

This approximation works well for early times $1\ll Ht\ll Ht_{\psi}$ times,
when $\psi$ has not reached its equilibrium distribution. At later times a
better approximation is obtained by considering that $M(t)$ is time
independent and we will comment on this later. Following the same method a
straightforward computation shows that $L(t)$ is

$\displaystyle
L(t)=-\frac{i2f_{\zeta}Ht_{\psi}p}{(3H+2ipt_{\psi}\Omega^{2}\Delta_{\zeta}^{2}(0))}e^{-\frac{t}{t_{\psi}}-\frac{2ip\Omega^{2}\sigma_{\zeta}^{2}}{3H^{2}}}\left(1-e^{\frac{t}{t_{\psi}}+\frac{2ip\Omega^{2}\Delta_{\zeta}^{2}(0)}{3H}}\right).$
(4.30)

Using the solution (4.29) and (4.30) we can immediately compute the Fourier
transform in $k$ since the PDF is Gaussian. Computing the Fourier transform in
$p$ is non trivial, since $p$ appears in the exponents. In order to Fourier
transform we expand $M(t)$ and $L(t)$ in powers of $\Delta_{\zeta}^{2}(0)$.
Notice that this requires $\zeta\gg\Delta_{\zeta}^{2}(0)\Omega^{2}/\mu^{2}$
which we will assume to hold true and indeed it does for typical values of
$\zeta$ and small amplitude of the density perturbations. Retaining terms only
up to second order in $p$ and Fourier transforming back, we finally obtain:

$\displaystyle\boxed{P(\zeta,\psi)=\exp\left[-\frac{\psi^{2}}{2\sigma_{\psi}^{2}}-\frac{1}{2\sigma_{\zeta}^{2}}\left(\zeta+\frac{2f_{\zeta}\Omega}{3H}\frac{\sigma_{\zeta}^{2}}{\sigma_{\psi}^{2}}\psi-\frac{\Omega^{2}}{3H^{2}}\frac{\sigma_{\zeta}^{2}}{\sigma_{\psi}^{2}}\psi^{2}\right)^{2}\right].}$
(4.31)

which is valid for $1\ll Ht\ll Ht_{\psi}$.

###  4.2 Non Gaussian tails

In order to understand the effect of the non linear interactions in the PDF it
is useful to ignore first the linear mixing term. This is not well justified
since it means that we are ignoring a $\psi\zeta$ interaction that is
important for large $\zeta$. Nevertheless its addition does not change the
qualitative effect of adding the non linear derivative interaction. As we
discuss in Section 3.1, this is because the effect of the linear mixing term
is shifting the trajectory of $\psi$ but not its variance. Due to this let us
ignore this term for the moment. Furthermore, we see that by doing this, it
becomes possible to integrate over $\psi$ analytically which greatly
simplifies the analysis. After these considerations let us study the following
PDF

$\displaystyle\boxed{P(\zeta,\psi)=\exp\left[-\frac{\psi^{2}}{2\sigma_{\psi}^{2}}-\frac{1}{2\sigma_{\zeta}^{2}}\left(\zeta-\bar{\kappa}\frac{\psi^{2}}{2\sigma_{\psi}^{2}}\right)^{2}\right],}$
(4.32)

where we have defined

$\displaystyle\bar{\kappa}\equiv\frac{\Omega^{2}}{3H^{2}}\sigma_{\zeta}^{2},$
(4.33)

which is a time independent parameter. Notice that this coupling is related to
(3.35) through the relation

$\displaystyle\bar{\kappa}=\frac{\kappa^{2}}{\sigma_{\psi}^{2}},$ (4.34)

by which it should be clear that $\kappa$ is related to the size of the
interactions. We are interested in the distribution for $\zeta$ which we
obtain after integrating over $\psi$. The integral can be done analytically
and expressed in term of Bessel functions:

$\displaystyle
P(\zeta)\propto\begin{cases}\exp\left(\frac{(\zeta\bar{\kappa}-\sigma_{\zeta}^{2})^{2}}{4\bar{\kappa}^{2}\sigma_{\zeta}^{2}}-\frac{\zeta^{2}}{2\sigma_{\zeta}^{2}}\right)K_{1/4}\left(\frac{(\zeta\bar{\kappa}-\sigma_{\zeta}^{2})^{2}}{4\bar{\kappa}^{2}\sigma_{\zeta}^{2}}\right)\hskip
128.0374pt\mathrm{for}\ \zeta<\frac{\sigma^{2}_{\zeta}}{\bar{\kappa}}\\\
\exp\left(\frac{(\zeta\bar{\kappa}-\sigma_{\zeta}^{2})^{2}}{4\bar{\kappa}^{2}\sigma_{\zeta}^{2}}-\frac{\zeta^{2}}{4\sigma_{\zeta}^{2}}\right)\left(I_{-1/4}\left(\frac{(\zeta\bar{\kappa}-\sigma_{\zeta}^{2})^{2}}{4\bar{\kappa}^{2}\sigma_{\zeta}^{2}}\right)+I_{1/4}\left(\frac{(\zeta\bar{\kappa}-\sigma_{\zeta}^{2})^{2}}{4\bar{\kappa}^{2}\sigma_{\zeta}^{2}}\right)\right)\qquad\mathrm{for}\
\zeta>\frac{\sigma^{2}_{\zeta}}{\bar{\kappa}}\end{cases}\\!\\!\\!.$ (4.35)

Figure 2: Toy model of the distribution described in (4.35). The black region
is the when $\zeta<\sigma_{\zeta}^{2}/\bar{\kappa}$ while the orange line is
otherwise. The dotted line is a Gaussian distribution that fits the region in
the left of the plot. The parameters are made up to highlight the fact that
the distribution on the right has a very non Gaussian tail.

To appreciate the distribution we plot it in Fig. 2. From there we see that
for it is a displaced Gaussian around the centre but it becomes strongly non
Gaussian for
$\zeta>\zeta_{\mathrm{cr}}\equiv\frac{\sigma^{2}_{\zeta}}{\bar{\kappa}}$. We
can understand the asymptotic behavior of $P(\zeta)$ as emerging from a change
on the saddle points in (4.32). Let us notice from it, that the shifted
distribution of $\zeta$ has the overall effect of changing the coefficient in
front of $\psi^{2}$, which becomes negative for large values of
$\zeta>\zeta_{\mathrm{cr}}$. This implies that there are three saddle points
for $\psi$, one at $\psi=0$, and other two at

$\displaystyle\bar{\psi}=\pm\sqrt{\frac{2\sigma_{\psi}^{2}}{\bar{\kappa}}}\sqrt{\zeta-\frac{\sigma_{\zeta}^{2}}{\bar{\kappa}}}\
.$ (4.36)

The behaviour for the asymptotics of $P(\zeta)$ are then similar to the Stokes
phenomena, in the sense that for large values of $\zeta$ the saddle point
changes from $0$ to $\bar{\psi}$. If we expand around $0$ we find a Gaussian
distribution for $\zeta$, whereas if we expand (4.32) around $\bar{\psi}$
after integrating over $\psi$, we find that

$\displaystyle P(\zeta)\sim\exp\left(-\frac{\zeta}{\bar{\kappa}}\right),$
(4.37)

valid for $\zeta\gg\zeta_{\mathrm{cr}}$ and that it coincides with the large
$\zeta$ limit of (4.35). This effect is non perturbative in nature since
$\bar{\kappa}\ll 1$ and so typically $\zeta/\bar{\kappa}\geq 1$. Let us now
comment on the regime of validity of $\bar{\kappa}$. As we mentioned before
Eq. (4.29) was only valid for intermediate times $t_{v}\ll t\ll t_{\psi}$,
whereas for later $t\geq t_{\psi}$ it is more accurate to consider that $M(t)$
is time independent. Solving for $M$ and $L$ we find that the PDF is given by,

$\displaystyle
P(\zeta,\psi)=\exp\left(-\frac{\psi^{2}}{2\sigma_{0\psi}^{2}}-\frac{1}{2\sigma_{\zeta}^{2}}\left(\zeta-\bar{\kappa}^{\prime}\frac{\psi^{2}}{2\sigma_{0\psi}^{2}}\right)^{2}\right)\
,\qquad\bar{\kappa}^{\prime}=\frac{1}{3}t_{\psi}H\Delta_{\zeta}^{2}(t)\frac{\Omega^{2}}{H^{2}}\
,$ (4.38)

with $\bar{\kappa}$ replaced by $\bar{\kappa}^{\prime}$ and
$\sigma_{0\psi}^{2}=3H^{4}/(8\pi^{2}\mu^{2})$ the equilibrium distribution for
$\psi$. In absence of $\zeta$, the PDF in (4.38) reduces to the equilibrium
distribution of $\psi$, in accordance with the fact the we are considering the
distribution at late times $t\geq t_{\psi}$. In this sense, (4.38) is the
distribution when the field $\psi$ has settled into its equilibrium
distribution. Of course if the field is ultralight then $t_{\psi}\to\infty$
and the transition between $\bar{\kappa}$ and $\bar{\kappa}^{\prime}$ does not
take place during inflation.

Furthermore, we see that for $t\ll t_{\psi}$,
$\bar{\kappa}<\bar{\kappa}^{\prime}$ but otherwise
$\bar{\kappa}>\bar{\kappa}^{\prime}$ since $\bar{\kappa}$ keeps growing.
Expanding in powers of $t_{\psi}$ we find that

$\displaystyle\bar{\kappa}=\frac{\Omega^{2}}{3H^{2}}Ht\Delta_{\zeta}^{2}(t)\
,\qquad\bar{\kappa}^{\prime}=\frac{\Omega^{2}}{3H^{2}}Ht_{\psi}\Delta_{\zeta}^{2}(t)\
,$ (4.39)

from where we can see that $\bar{\kappa}$ grows until it reaches the
equilibrium value $\bar{\kappa}^{\prime}$. Anyhow since at equilibrium
$\sigma_{\zeta}\approx Ht_{\psi}\Delta_{\zeta}^{2}$, we can write
$\kappa^{\prime}\sim\Omega^{2}/H^{2}\sigma_{\zeta}^{2}$ to see that the
coupling does not change its dependence on the parameters. In any case is
worth mentioning that $\bar{\kappa}^{\prime}$ is a limiting value where the
time dependence has become negligible.

All of this implies that initially the tail of the distribution becomes non
Gaussian at smaller values of $\zeta$ until it settles down on
$\sigma_{\zeta}^{2}/\bar{\kappa}^{\prime}$. When this happens the coupling of
the tail becomes constant

$\displaystyle
P(\zeta)\sim\exp\left(-\frac{\zeta}{\bar{\kappa}}\right)=\exp\left(-\frac{\mu^{2}}{2\Omega^{2}\Delta_{\zeta}^{2}(t)}\zeta\right).$
(4.40)

Still at early times the distribution is more localised, hence the values at
which the tail becomes non-Gaussian are smaller, as can be seen in Fig. 3. In
the end, even though the saddle point changes at larger $\zeta$, larger values
of $\zeta$ are more likely due to the growth of $\sigma_{\zeta}^{2}$. Finally,
let us pay attention to the fact that the tail is typically very suppressed
since we have that

Figure 3: a) The figure is the PDF of (4.35). The orange line is for $N=30$
and the black line is for $N=40$. The other parameters are $\Omega=.2H$ and
$\Delta_{\zeta}=10{-3}$. The dotted lines are the region where
$\zeta\geq\sigma_{\zeta}^{2}/\bar{\kappa}$ b) Both curves are at $N=40$ but
the green line has $\Omega=0.3H$.

$\displaystyle\zeta/\bar{\kappa}\gg\sigma_{\zeta}^{2}/\bar{\kappa}^{2}\gg 1\
,$ (4.41)

since $\bar{\kappa}^{-1}\gg 1$. In order to have a larger effect the power
spectrum has to be several orders of magnitude larger than the CMB values.
Furthermore, let us note that the critical value for $\zeta$ does not depend
on the amplitude of the fluctuations $\Delta_{\zeta}^{2}(t)$, which in the end
implies that for smaller values of $\Delta_{\zeta}^{2}(t)$ the probability for
the regions where the tail changes is very suppressed. In any case the tail we
have found is always larger than a Gaussian tail, which is due to the fact
that $\bar{\kappa}\zeta/\sigma_{\zeta}^{2}>1$ is equivalent to
$\zeta/\bar{\kappa}>\zeta^{2}/\sigma_{\zeta}^{2}$.

To conclude let us stress the point that changing $\Omega$ has a large effect,
as can be seen in Fig. 3b, where the overall effect is flattening and shifting
the PDF.

### Local non-Gaussianity

In order to understand how the tail is related to usual perturbation theory,
we can estimate the size of the non-Gaussianities produced by the interaction
$\sim v_{\zeta}^{2}\psi$ we are considering. From the Langrangian (4.2), we
have that $f_{\mathrm{NL}}$ can be estimated to be of order

$\displaystyle
f_{\mathrm{NL}}\zeta\sim\frac{\mathcal{L}_{3}}{\mathcal{L}_{2}}\sim\frac{6\Omega^{2}}{H^{2}f_{\zeta}^{2}}\frac{\psi^{2}}{\zeta}\sim\frac{\Omega^{2}}{H^{2}f_{\zeta}^{2}}\frac{\sigma_{\psi}^{2}}{\sigma_{\zeta}^{2}}\zeta\
,$ (4.42)

where the last term corresponds to the off diagonal coupling parameter.
Indeed, using (3.35) the last relation can be recast as

$\displaystyle f_{\mathrm{NL}}\zeta\sim\frac{\Omega^{2}}{H^{2}}\zeta\ ,$
(4.43)

which implies the following: Perturbation theory is usually under control when
$f_{\mathrm{NL}}\zeta\ll 1$ which in our case corresponds to expanding around
the Gaussian saddle point. Observables can be computed by expanding them in
powers of the power spectrum since it is always suppressed. On the other hand,
when $f_{\mathrm{NL}}\zeta\sim 1$ perturbation theory fails since the
distribution is non Gaussian, which means that the expansion in powers of the
power spectrum cannot be justified. We can see this explicitly now since
$\zeta_{\mathrm{cr}}\sim 1/f_{\mathrm{NL}}$. Naively one should expect that
for $\zeta\geq\zeta_{\mathrm{cr}}$ perturbation theory fails which we now see
translates into the tail becoming strongly non-Gaussian. Clearly if
$f_{\mathrm{{NL}}}\ll 1$ one might be worried about corrections from other
interactions becoming important for large values of the tail. Nevertheless,
since there is a finite number of interactions in the EFT we are considering
we can still deduce the behaviour of the tail up to certain (large) values of
$\zeta$ when other interactions become important. In order to estimate the
validity of our results we need to include the terms we have been neglecting
so far, which is what we are going to discuss now.

#### Including the linear mixing

In the previous section we ignored the linear mixing term, since we argued it
does not affect the appearance of non-Gaussian tails. When adding it, it turns
out that it is not possible to integrate over $\psi$ analytically but it is
still possible to obtain the saddle points and check how the tail changes for
large $\zeta$. Firstly, let us notice that the point where the saddle point
changes is not modified significantly. To see this we can expand the PDF and
notice that the terms proportional to $\psi^{2}$ in Eq. (4.31) are

$\displaystyle\frac{\psi^{2}}{2\sigma_{\zeta}^{2}}\left(1-\frac{4\Omega^{2}}{3H^{2}}\left(\zeta-\frac{1}{3}\right)\right),$
(4.44)

which implies that for small $\zeta<1/3$ the saddle point does not change even
when $\Omega^{2}/H^{2}$ is very large. Since we are interested in small
$\Omega^{2}/H^{2}$ this does not significantly change our results. Next, if we
ignore suppressed terms we find that the displaced saddle points are now at

$\displaystyle\bar{\psi}=\frac{3f_{\zeta}H}{\Omega}\pm\frac{3f_{\zeta}H^{2}}{\Omega^{2}}\sqrt{-1/2+(1+24\zeta)\frac{\Omega^{2}}{36H^{2}}}.$
(4.45)

The fact that the saddle point now contains a constant piece is reflected in
the tail. Indeed, expanding around $\bar{\psi}$ we find that for large $\zeta$

$\displaystyle
P(\zeta)\sim\exp\left(-\frac{3\zeta\pm\sqrt{6\zeta}}{4\sigma_{\zeta}^{2}}\frac{H^{2}}{\Omega^{2}}\right),$
(4.46)

where the different signs correspond to different saddle points. In general
there is one correct saddle point which can be picked based on the analytical
properties of the PDF, a task which we leave for future work.333See Refs. [61,
62, 63] for related attempts in dealing with complex saddle points. In any
case, if we assume that the fluctuation in $\zeta$ is large we find that the
tail gets a small correction. Finally let us notice that the effects of the
linear mixing term can also be understood as $f_{\mathrm{NL}}\zeta$ becoming
large. In this case it corresponds to the interaction $\mathcal{L}_{3}\sim
2f_{\zeta}\Omega/H\psi v_{\zeta}^{2}$, and so the constant piece of the saddle
point, which appears for smaller $\zeta$, appears as long as
$\Omega/H\zeta\sim 1$. Nevertheless, the tail of the distribution will not
change until $\Omega^{2}/H^{2}\zeta\sim 1$ where the other saddle point
becomes real.

As a side comment let us notice that the coupling does not grow unbounded but
it settles into a constant value at approximately $t\sim t_{\psi}$, given by

$\displaystyle\frac{2f_{\zeta}\Omega}{3H^{2}}\frac{t_{\psi}\Delta_{\zeta}^{2}(t)}{\sigma_{\psi}^{2}}\psi\sim\frac{\Omega
t_{\psi}}{f_{\zeta}}\ ,$ (4.47)

as explained in Eq. (3.35).

###  4.3 Adding more interactions

Having understood the leading interaction effect we can now add the rest of
the terms to the system. In particular, let us consider a general potential
$V(\psi)$. Since the new terms do not add any new shorter time scale than
$t_{v}$ we can eliminate $v_{\zeta}$ from the Fokker-Plank equation following
the same steps we described in Section 4.1.1. This results in the following
Fokker-Plank equation

$\displaystyle\frac{\partial P}{\partial t}$
$\displaystyle=\frac{\partial}{\partial\psi}\left(\frac{V^{\prime}(\psi)}{3H}P+\frac{D_{\psi}}{2}\frac{\partial
P}{\partial\psi}\right)+H\Delta_{\zeta}^{2}(t)\frac{\partial^{2}}{\partial\psi\partial\zeta}\left(\left(\frac{2\Omega^{2}}{H^{2}}\psi+\frac{2f_{\zeta}\Omega}{3H}\right)P\right)$
$\displaystyle+\frac{H\Delta_{\zeta}(t)^{2}}{2}\frac{\partial^{2}P}{\partial\zeta^{2}}+\frac{6\Omega^{2}}{f_{\zeta}^{2}H}\frac{\partial}{\partial\zeta}\left(\psi^{2}P\right)-\frac{\Omega\Delta_{\zeta}^{2}(t)}{f_{\zeta}}\frac{\partial^{2}}{\partial\zeta^{2}}(\psi
P)-\frac{1}{4}\Omega
f_{\zeta}\Delta_{\zeta}^{4}(t)\frac{\partial^{3}P}{\partial^{2}\zeta\partial\psi}\
,$ (4.48)

where we have ignored sub leading terms in the noise, and expanded to cubic
order in $t_{v}$ in order to obtain the last two terms. Before solving the
Fokker-Planck equation let us pay attention to the fact that $\zeta$ appears
only through derivatives of the PDF. This stems from the fact that $\zeta$
posses a shift symmetry, by which the only allowed interaction contains time
derivatives or gradients. As for the stochastic dynamics, this result in the
drift being an explicit function of $v_{\zeta}$. Now, at leading order in
$t_{v}$, we see that this translates into derivatives of $\zeta$. Indeed
allowing higher powers of $v_{\zeta}$ in the drift, translates into higher
derivatives of $\zeta$ in the Fokker-Planck equation [33, 64]. Another
important feature of Eq. (4.48) is that, as expected, in the absence of the
coupling $\Omega$ the equation reduces to uncoupled Fokker-Plank equations
whose solution factorises as $P\sim P_{\psi}(\psi,t)P_{\zeta}(\zeta,t)$.

Finally let us stress that in order to obtain Eq. (4.48) we only need to
assume that $f_{\zeta}^{2}\gg H^{2}$, $\Omega^{2}\ll H^{2}$, and $\mu^{2}\ll
H^{2}$. The assumption about $\Omega$ can be relaxed, but doing that will
modify the value of the two point function which will modify the noise
function.

#### Including higher $\zeta$ derivatives

Let us now study the last two terms in (4.48) and restrict to
$V(\psi)=\mu^{2}\psi^{2}/2$. These contain two derivatives of $\zeta$, and as
expected will modify the variance of the distribution for $\zeta$. As we will
see they restrict the range of the Gaussian fluctuations for $\psi$. Indeed we
can find the solution for the PDF at equilibrium to be

$\displaystyle\log
P(\zeta,\psi)\propto-\frac{\psi^{2}}{2\sigma_{\psi}^{2}}-\frac{\left(\zeta+\frac{2\Omega
f_{\zeta}}{3H}\frac{\sigma_{\zeta}^{2}}{\sigma_{\psi}^{2}}\psi-\frac{\Omega^{2}}{3H^{2}}\frac{\sigma_{\zeta}^{2}}{\sigma_{\psi}^{2}}\psi^{2}\right)^{2}}{\sigma_{\zeta}^{2}(1+\frac{3\Omega
f_{\zeta}}{t_{\psi}H}\frac{\sigma_{\zeta}^{2}}{2\sigma_{\psi}^{2}}\psi)}\ ,$
(4.49)

where we have neglected other terms that appear in the denominator which are
suppressed by additional powers of $\Delta_{\zeta}^{2}$. We see from (4.49)
that the variance gets shifted for larger values of $\psi$. Notice, moreover,
that for typical values of $\psi$ this effect is very suppressed since
$\Omega\ll f_{\zeta}$. In any case, at that point the distribution is not
valid and other terms will become dominant. This effect would modify the tail
of $\psi$, and as long as we consider small fluctuations around the new
trajectory of $\psi$ the description we have given for the distribution in
$\zeta$ remains valid. If we expand in small $\psi$ we find that there are new
terms proportional to $\psi^{2}\zeta^{2}$. These terms would actually change
the saddle point moving it back to $\psi=0$ at larger values of $\zeta$. This
effect is related to the condition we impose in (4.12). In the end we find
that a violation of (4.12) is related to the fact that there is no change in
the saddle points due to the quadratic term in $\psi$. In general we found
that for late times and larger $\Omega^{2}/H^{2}$ the effect of these new
terms is suppressed. The reason behind this is that for late times $t\gtrsim
t_{\psi}$ the variance of $\sigma_{\psi}$ stops growing whereas the variance
of $\zeta$ grows until the end of inflation.

#### General $V(\psi)$

We have focused only on quadratic potentials for $\psi$ but it possible to
study a general potential. To start with, let us neglect the linear coupling
and focus only on the leading order interactions. If we Fourier transform in
$p$ and write $P=\exp(-\sigma_{\zeta}^{2}/2p^{2})\tilde{P}$ we obtain a
Fokker-Planck equation with only derivatives of $\psi$ and with a corrected
drift term

$\displaystyle\frac{\partial\tilde{P}}{\partial
t}=\frac{\partial}{\partial\psi}\left(\left(\frac{V^{\prime}(\psi)}{3H}+ip\frac{2\Omega^{2}\Delta_{\zeta}^{2}(t)}{H}\psi\right)\tilde{P}\right)-\frac{D_{\psi}}{2}\frac{\partial^{2}\tilde{P}}{\partial\psi^{2}}\
.$ (4.50)

We can rewrite this equation using the replacement

$\displaystyle\tilde{P}(\zeta,\psi,t)=e^{-v/2}\sum_{n=0}^{\infty}a_{n}\Phi_{n}e^{\Lambda_{n}t},\qquad
v(\zeta,\psi,t)=\frac{8\pi^{2}}{3H}\left(V(\psi)+3ip\frac{\Omega^{2}}{H^{2}}\Delta_{\zeta}^{2}(t)\psi^{2}\right)\
,$ (4.51)

by which the Fokker-Planck equation becomes an eigenvalue problem of the form

$\displaystyle\frac{1}{2}\left(-\frac{\partial}{\partial\psi}+\frac{\partial
v}{\partial\psi}\right)\left(\frac{\partial}{\partial\psi}+\frac{\partial
v}{\partial\psi}\right)\Phi_{n}(\psi)=\frac{8\pi^{2}\Lambda_{n}}{H^{3}}\Phi_{n}(\psi)\
,$ (4.52)

which is similar to Eq. (2.45) with the particular difference that the
equilibrium distribution $\exp(-v)$ in this case is time independent. Solving
this equation in general lies beyond the scope of this paper, but we can still
comment on the case that $\psi$ reaches the equilibrium PDF. Fourier
transforming in $p$ we find that the PDF is given by

$\displaystyle
P(\zeta,\psi)\sim\exp\left[-\frac{8\pi^{2}V(\psi)}{3H^{2}}-\frac{1}{2\sigma_{\zeta}^{2}}\left(\zeta-\frac{\bar{\kappa}^{\prime}}{2\sigma_{0\psi}^{2}}\psi^{2}\right)^{2}\right]\
,$ (4.53)

which is a generalisation of (4.38), including general potentials. Let us
consider the example $V(\psi)=\mu^{2}\psi^{2}/2+\lambda\psi^{4}/4$ which
allows us to integrate $\psi$ out of the equilibrium distribution. This leads
to

Figure 4: We plot the log of the $\zeta$ distribution given by (4.54) for
three different values of $g$. The black line corresponds to $\lambda=0$, the
red line to $\lambda=10^{-5}$ and the orange line to $\lambda=10^{-2}$. The
dotted lines are the regions where $3\zeta>Ht\frac{\mu^{2}}{\Omega^{2}}$.

$\displaystyle
P(\zeta)\propto\begin{cases}\exp\left(\frac{1}{4\sigma_{\zeta}^{2}}\frac{(\zeta-\sigma_{\zeta}^{2}/\bar{\kappa}^{\prime})^{2}}{1+\frac{6\lambda
t_{\psi}\sigma_{\psi}^{2}\sigma_{\zeta}^{2}}{H\bar{\kappa}^{\prime
2}}}-\frac{\zeta^{2}}{2\sigma_{\zeta}^{2}}\right)K_{1/4}\left(\frac{1}{4\sigma_{\zeta}^{2}}\frac{(\zeta-\sigma_{\zeta}^{2}/\bar{\kappa}^{\prime})^{2}}{1+\frac{6\lambda
t_{\psi}\sigma_{\psi}^{2}\sigma_{\zeta}^{2}}{H\bar{\kappa}^{\prime
2}}^{2}}\right)\\\ \hskip 239.00298pt\mathrm{for}\
\zeta<\frac{\sigma^{2}_{\zeta}}{\bar{\kappa}^{\prime}}\ ,\\\ \\\
\exp\left(\frac{1}{4\sigma_{\zeta}^{2}}\frac{(\zeta-\sigma_{\zeta}^{2}/\bar{\kappa}^{\prime})^{2}}{1+\frac{6\lambda
t_{\psi}\sigma_{\psi}^{2}\sigma_{\zeta}^{2}}{H\bar{\kappa}^{\prime
2}}}-\frac{\zeta^{2}}{2\sigma_{\zeta}^{2}}\right)\left(I_{-1/4}\left(\frac{1}{4\sigma_{\zeta}^{2}}\frac{(\zeta-\sigma_{\zeta}^{2}/\bar{\kappa}^{\prime})^{2}}{1+\frac{6\lambda
t_{\psi}\sigma_{\psi}^{2}\sigma_{\zeta}^{2}}{H\bar{\kappa}^{\prime
2}}}\right)+I_{1/4}\left(\frac{1}{4\sigma_{\zeta}^{2}}\frac{(\zeta-\sigma_{\zeta}^{2}/\bar{\kappa}^{\prime})^{2}}{1+\frac{6\lambda
t_{\psi}\sigma_{\psi}^{2}\sigma_{\zeta}^{2}}{H\bar{\kappa}^{\prime
2}}}\right)\right)\\\ \hskip 239.00298pt\mathrm{for}\
\zeta>\frac{\sigma^{2}_{\zeta}}{\bar{\kappa}^{\prime}}\ .\end{cases}$ (4.54)

This expression is a bit complicated but we can see that a quartic self
interaction reduces the amplitude of the tail. This has an overall effect of
making the distribution Gaussian for larger values of $\zeta$. We plot (4.54)
in Fig. 4 where it is noticeable how the non Gaussian effects diminish by
increasing $\lambda$. This can be shown by expanding around the non trivial
saddle point. In this case after integrating out $\psi$ we find that

$\displaystyle
P(\zeta)\sim\exp\left(-\frac{\zeta^{2}}{\frac{H\bar{\kappa}^{2}}{3t_{\psi}\lambda\sigma_{\zeta}\sigma_{\psi}^{2}}+2\sigma_{\zeta}^{2}}-\frac{\zeta}{\bar{\kappa}+\frac{6t_{\psi}\lambda\sigma_{\psi}^{2}\sigma_{\zeta}^{2}}{H\bar{\kappa}}}\right).$
(4.55)

Notice that the tail contains a quadratic term, whose variance has a
correction which is the inverse of the correction of the coefficient of the
tail. This means that increasing $\lambda$ makes the correction of the
quadratic term smaller while the correction for the linear term becomes
larger. The effect can be understood by estimating the sizes of the non
Gaussianites. Indeed we have that for the quartic interaction
$f_{\mathrm{NL}}^{(\lambda)}$, is of order

$\displaystyle f_{\mathrm{NL}}^{(\lambda)}\
\zeta\sim\frac{\lambda}{H^{2}f_{\zeta}^{2}}\frac{\sigma_{\psi}^{4}}{\sigma_{\zeta}^{2}}\times\zeta.$
(4.56)

If we compare the ratio between $f_{\mathrm{NL}}^{(\lambda)}$ and
$f_{\mathrm{NL}}$ we find that it coincides with the ratio that controls
whether the Gaussian term dominates

$\displaystyle\frac{f_{\mathrm{NL}}^{(\lambda)}}{f_{\mathrm{NL}}}\sim\frac{\lambda}{\Omega^{2}}\sigma_{\psi}^{2}.$
(4.57)

Let us first detail the case in which the ratio is much larger than one. When
this happens, since perturbation theory breaks down for
$f_{\mathrm{NL}}^{(\lambda)}\zeta\sim 1$, the action cannot be trusted anymore
and the effect of the cubic coupling between $\zeta$ and $\Omega$ is not seen.
Clearly at this point other self interactions have to be taken into account
and the PDF at large $\zeta$ might become dominated by other higher order
terms. On the other hand when the ratio in (4.57) becomes smaller than one the
saddle point changes before $f_{\mathrm{NL}}^{(\lambda)}\zeta\sim 1$ and so
the tail becomes non Gaussian. For any non zero $\lambda$ there it will be a
point when $f_{\mathrm{NL}}^{(\lambda)}\zeta\sim 1$. At that stage the
computations of the tail in are not valid. Nevertheless this will produce an
exponentially small effect on the whole PDF.

##  5 Conclusions

We have studied the statistics of large but rare fluctuations within the
multi-field inflation paradigm using the stochastic inflation formalism. In
the simplest class of two-field models, the primordial curvature fluctuation
$\zeta$ interacts with an isocurvature field $\psi$ as a result of turns of
the background trajectory in the target space of scalar fields. This
translates into a derivative coupling proportional to $\Omega$, the rate of
turn of the trajectory, appearing at both linear and non-linear level in the
evolution of perturbations. We found that the non-linear interactions induced
by $\Omega$ imply non-Gaussian deformations affecting the tails of the joint
probability distribution of the perturbations.

By assuming that the evolution of the background is close to de Sitter, we
derived the Fokker-Planck equation that is satisfied by the probability
density function characterising both fields $\zeta$ and $\psi$. When only the
linear evolution of the fields is considered, we find that a non vanishing
$\Omega$ enhances the growth of the variance of $\zeta$. This matches with
results obtained previously using perturbation theory. A particular case of
this scenario is when the entropy mass is exactly zero (the ultralight limit),
studied recently in Ref. [49]. We showed that initially all spectator fields
behave as ultralight fields after horizon crossing but after some time, that
depends on their mass, the fields decay. If the entropy mass is zero, we
recover the exact ultralight case.

On the other hand, when non linearities are taken into account, we find that
after integrating out $\psi$ the tail of the PDF of $\zeta$ becomes non-
Gaussian for values of $\zeta\gtrsim H^{2}/\Omega^{2}$. This can be understood
as the Stokes phenomenon, whereby a Gaussian saddle point leads to non-
Gaussian saddle points for large values of a parameter. Crucially the coupling
makes the exponential tail to be larger than the Gaussian tail. Such a result
has important consequences. For instance the abundance of PBHs formed during
inflation depends strongly on the tail of the PDF. In this way a non-Gaussian
tail implies that the abundance of PBHs can be substantially enhanced in
models with derivative couplings. Another consequence would be a possible
modification on the clustering of galaxies which depend on rare large
fluctuations, whose probability would now be enhanced.

There are several paths along which our work can be expanded. For instance,
here we considered a simple EFT of multi field inflation in which the number
of interaction terms with derivatives couplings is limited. There are other
known examples with a larger number of interaction terms which might become
relevant for large values of $\zeta$. In those cases one might need to resum
the implied corrections to obtain accurate expressions for the tail. Another
interesting topic would be to understand how corrections to the Fokker-Planck
equation are related to resummation of loops. This has been well understood
for light spectator fields on pure de Sitter [28, 65, 38, 20], and a similar
result should follow from the Fokker-Planck equation arising from our two
field model. Also, it would be interesting to compare our results to other
derivations of the Fokker-Planck equations within the multi-field paradigm
[27, 66, 67, 41] and its relation with the Hamilton-Jacobi formalism [68, 50].
On these examples the derivation of the Fokker-Planck equation was done
directly from the background equations, whereas in this article we obtained
the stochastic dynamics directly from perturbations.

Our results suggest that the abundance of PBHs in multifield models can be
much larger than that obtained when the Gaussian approximation is used to
study the production of PBHs through the enhancement of the power spectrum. It
would be interesting to apply our results to models such us those of Refs.
[69, 70] to reassess the production of PBHs. To do this one would need to go
beyond the assumption that $\Omega\ll H$, implying that some of the terms that
in our analysis were suppressed would now become dominant. Finally we see that
stochastic inflation might allow one to go beyond perturbation theory. In that
sense it will be interesting to understand our results in the light of recent
works such as [12, 59] (see [71, 72] for a discussion on the importance of
this for PBHs). For instance it has been suggested that when there is a non
perturbative tail there is an exponential enhancement of the large $N$ point
correlation function [73]. Whether those result apply to the case we study
here, we leave for future work.

## Acknowledgements

We are grateful to Guillermo Ballesteros, Lucas Pinol and Spiros Sypsas for
useful discussions and comments on this work. G.A.P is supported by the
Fondecyt Regular Project No. 1210876 (ANID). The work of SC has been funded by
a Contrato de Atracción de Talento (Modalidad 1) de la Comunidad de Madrid
(Spain), number 2017-T1/TIC-5520 and the IFT Centro de Excelencia Severo Ochoa
Grant SEV-2016-0597. ACD acknowledges partial support from STFC consolidated
grant ST/T000694/1. AA’s work is partially supported by the Netherlands
Organization for Scientific Research (N.W.O), by the Basque Government
(IT-979-16) and by the Spanish Ministry MINECO (FPA2015-64041-C2-1P).

##  Appendix A Langevin equations

In this appendix we give details on how to obtain the Langevin equation (2.3).
Starting from the action

$\displaystyle S=\int
d^{4}xdta^{3}\epsilon\left[\dot{\zeta}^{2}-\frac{1}{a^{2}}(\nabla\zeta)^{2}\right]\
,$ (A.1)

the equations of motion are

$\displaystyle\ddot{\zeta}+3H\dot{\zeta}+\frac{k^{2}}{a^{2}}\zeta=0\ .$ (A.2)

We can split the solution to this equation into long- and short-wavelengths

$\displaystyle\zeta(x)$ $\displaystyle=\zeta_{l}(x)+\zeta_{s}$
$\displaystyle=\zeta_{l}+\int\frac{d^{3}k}{(2\pi)^{3}}\theta(k-\epsilon
a(t)H)\left[a_{k}\zeta_{k}(t)e^{-i\vec{k}\cdot\vec{x}}+a_{k}^{\dagger}\zeta^{*}_{k}e^{i\vec{k}\cdot\vec{x}}\right]\
,$ (A.3)

with $\epsilon\ll 1$ and where $\zeta_{k}$ is the solution of (A.2)

$\displaystyle\zeta_{k}=z\sqrt{\frac{\pi}{4}}H(-\tau)^{3/2}H_{3/2}^{(1)}(-k\tau)\
,$ (A.4)

where $a,a^{\dagger}$ are the corresponding creation and annihilation
operators and $\tau$ is conformal time. Notice that picking a different window
function might change the result of the correlation function (see [74] for a
discussion) although in this work we are only interested in equal-time
correlation functions. Clearly the definition of the long wavelength modes
implies that we neglect the gradient term in the equation of motion. Moreover,
since horizon crossing happens at $k/aH\simeq 1$, it also implies that the
leading piece from the equation of motion is $3H\dot{\zeta}_{l}$. Due to this
we can write the equation of motion as

$\displaystyle\ddot{\zeta}_{l}+3H\dot{\zeta}_{l}=\eta_{\zeta}(t,x)\ ,$ (A.5)

where $\eta_{\zeta}(t,x)$ is the contribution from the short wavelength
perturbations.

###  A.1 Computing the noise terms

In order to compute the noise term let us consider a light scalar field $\phi$
in de Sitter space of mass $m^{2}\ll H^{2}$ . Going back to $\zeta$ can be
done by rescaling our results. Following [25], let us split the field into
long and short wavelengths by writing it as

$\displaystyle\phi(x)$ $\displaystyle=\phi_{l}(x)+\phi_{s}$
$\displaystyle\equiv\bar{\phi}+\int\frac{d^{3}k}{(2\pi)^{3}}\theta(k-\epsilon
a(t)H)\left[a_{k}\phi_{k}(t)e^{-i\vec{k}\cdot\vec{x}}+a_{k}^{\dagger}e^{i\vec{k}\cdot\vec{x}}\right]\
,$ (A.6)

where we have introduced the time dependent cut-off $k_{c}=\epsilon a(t)H$,
with $\epsilon$ a small parameter, and where the mode function $\phi_{k}$ is
given by

$\phi_{k}=\sqrt{\frac{\pi}{4}}H(-\tau)^{3/2}H_{\nu}^{(1)}(-k\tau),$ (A.7)

where $\tau$ is conformal time and $\nu=\sqrt{9/4-m^{2}/H^{2}}$ . The noise
term comes from averaging over the time derivative of the second term. Indeed
defining

$f\equiv\epsilon a(t)H^{2}\int\frac{d^{3}k}{(2\pi)^{3}}\delta(k-\epsilon
a(t)H)\left[a_{k}\phi_{k}(t)e^{-i\vec{k}\cdot\vec{x}}+a_{k}^{\dagger}e^{i\vec{k}\cdot\vec{x}}\right]\
,$ (A.8)

we then have that

$\displaystyle\langle f(x,t)f(x^{\prime},t^{\prime})\rangle$
$\displaystyle=\epsilon^{2}H^{4}\int\frac{d^{3}k}{(2\pi)^{3}}\phi_{k}\phi^{*}_{k}e^{-i\vec{k}\cdot(\vec{x}-\vec{x}^{\prime})}\delta(k-\epsilon
a(t)H)a(t)\delta(k-\epsilon a(t^{\prime})H)a(t^{\prime})$
$\displaystyle=\epsilon^{2}H^{4}\int\frac{k^{2}dk}{2\pi^{2}}\phi_{k}\phi_{k}^{*}\frac{\sin(k|\vec{x}-\vec{x}^{\prime}|)}{k|\vec{x}-\vec{x}^{\prime}|}\delta(k-\epsilon
a(t)H)a(t)\delta(k-\epsilon a(t^{\prime})H)a(t^{\prime}).$

Now for a massless field we have that

$\phi_{k}=\frac{H}{\sqrt{2k}}\left(\tau-\frac{i}{k}\right)e^{-ik\tau},$ (A.10)

so instead we obtain

$\displaystyle\langle f(x,t)f(x^{\prime},t^{\prime})\rangle$
$\displaystyle=\frac{\epsilon^{2}H^{4}}{4\pi^{2}}\frac{\sin(\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|)}{\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|}\frac{a(t)^{2}}{\epsilon
a(t)H}\frac{\delta(t-t^{\prime})}{\epsilon H^{2}a(t)}$
$\displaystyle=\frac{H^{3}}{4\pi^{2}}\frac{\sin(\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|)}{\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|}\delta(t-t^{\prime}),$ (A.11)

which implies that at coincident points we have

$\langle
f(x,t)f(x,t^{\prime})\rangle=\frac{H^{3}}{4\pi^{2}}\delta(t-t^{\prime}).$
(A.12)

In general we should expand $\phi_{k}\phi_{k}^{*}$ in small $\epsilon$ after
replacing the cut-off in the mode functions. At leading order in $\epsilon$ we
find that

$\phi_{k}\phi_{k}^{*}=\frac{\Gamma(\nu)^{2}4^{\nu-1}}{\pi}\frac{\epsilon^{-2\nu}}{Ha^{3}}.$
(A.13)

Replacing into (LABEL:eq:SW_correlation) we get

$\displaystyle\langle f(x,t)f(x^{\prime},t^{\prime})\rangle$
$\displaystyle=\frac{\epsilon^{2}H^{4}}{2\pi^{2}}\frac{\sin(\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|)}{\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|}\epsilon^{2}a(t)^{2}H^{2}\frac{\Gamma(\nu)^{2}4^{\nu-1}}{\pi}\frac{\epsilon^{-2\nu}}{Ha^{3}}\frac{a(t)^{2}\delta(t-t^{\prime})}{\epsilon
H^{2}a(t)}$ $\displaystyle=\frac{H^{3}}{4\pi^{2}}\frac{\sin(\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|)}{\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|}\frac{\Gamma(\nu)^{2}4^{\nu-1/2}}{\pi}\epsilon^{3-2\nu}\delta(t-t^{\prime}),$
(A.14)

Expanding in powers of $m/H$ we get,

$\displaystyle\langle f(x,t)f(x^{\prime},t^{\prime})\rangle$ $\displaystyle=$
$\displaystyle\frac{H^{3}}{4\pi^{2}}\frac{\sin(\epsilon
a(t)H|\vec{x}-\vec{x}^{\prime}|)}{\epsilon a(t)H|\vec{x}-\vec{x}^{\prime}|}$
(A.15)
$\displaystyle\\!\\!\\!\\!\\!\\!\times\left(1+2\left(-2+\gamma_{E}+\log(2\epsilon)\right)\frac{m^{2}}{H^{2}}+\mathcal{O}\left(m^{2}/H^{2}\right)\right)\delta(t-t^{\prime}).\quad$

##  Appendix B Secular growth of $\sigma_{\zeta}^{2}$

In this appendix we will relate the variances of the stochastic fields to the
correlation functions in real space. First let us recall that the field
perturbation in inflation can be written approximately as

$\displaystyle\delta\phi_{k}=\sqrt{\frac{\pi}{4}}H(-\tau)^{3/2}H_{3/2}^{(1)}(-k\tau),$
(B.1)

where $\tau$ is conformal time and $a(\tau)=-1/(H\tau)$. From this it is
possible to write the two point function in momentum space
$\langle\delta\phi_{k}^{2}\rangle=\frac{H^{2}}{2k^{3}}(1+k^{2}\tau^{2})$.
Since we would like to compare this to the result obtained through the
Langevin equation let us Fourier transform the two point function. We have
that

$\displaystyle\langle\delta\phi(\vec{x},t)\delta\phi(\vec{y},t)\rangle$
$\displaystyle=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot|\vec{x}-\vec{y}|}\langle\delta\phi_{k}^{2}\rangle$
$\displaystyle=\frac{H^{2}}{4\pi^{2}}\int
k^{2}dke^{ik|\vec{x}-\vec{y}|}\frac{1}{k^{3}}(1+k^{2}\tau^{2})$
$\displaystyle=\frac{1}{(2\pi)^{2}}\frac{1}{|\vec{x}-\vec{y}|^{2}a(t)^{2}}-\frac{H}{(2\pi)^{2}}\log\left(\frac{|\vec{x}-\vec{y}|}{L}\right),$
(B.2)

where the infrared cut-off $L\equiv a_{0}H_{0}$ is the largest scale during
inflation. This is related to the end of inflation as we can write
$(a_{0}H_{0})^{-1}=-\tau_{0}>0$ where $\tau_{0}<0$ is the time when inflation
ends. Notice also that during eternal inflation $\tau_{0}\to 0$ [75].

Eq. (B.2) can be understood as follows, the first piece corresponds to the
flat space two point function in physical coordinates. As we move deeper into
the bulk this expression dominates. This is expected as we have picked a
vacuum that reproduces the Minkowski vacuum. This term dilutes as we approach
the horizon. The second term does not depend explicitly on time and it
corresponds to the two point function on a scale invariant theory. This is due
to the symmetries of inflation at horizon crossing.

In order to compare to the stochastic result let us compute the correlation
functions at coincident points. This corresponds to the variance of
$\delta\phi$. Since the expression diverges let us introduce a time dependent
cut-off such that

$\displaystyle|\vec{x}-\vec{y}|e^{Ht}\gg\Lambda^{-1}.$ (B.3)

For shorter distances we evaluate the two point at the cut-off, which
schematically implies that

$\displaystyle\langle\delta\phi(t)^{2}\rangle$
$\displaystyle=a\Lambda^{2}+b\log\Lambda+\frac{H^{3}t}{(2\pi)^{2}}=\frac{H^{3}t}{(2\pi)^{2}}+\mathrm{const},$
(B.4)

which is the result we have obtained by solving the Fokker-Planck equation.
Notice that the time dependence comes from the fact that at each time more
modes are included in the region below the cut off. A similar computation
shows that at coincident points the $\dot{\phi}$ correlation function goes as

$\displaystyle\langle\delta\dot{\phi}(\vec{x},t)\delta\dot{\phi}(\vec{x},t^{\prime})\rangle_{\rm{vac}}$
$\displaystyle=\frac{1}{2\pi^{2}}\int_{0}^{\infty}d^{3}k|\delta\dot{\phi}_{k}(\tau)\delta\dot{\phi}_{k}(\tau^{\prime})|$
$\displaystyle=\frac{6H^{4}}{4\pi^{2}}\frac{\tau^{2}\tau^{\prime
2}}{(\tau-\tau^{\prime})^{4}}$
$\displaystyle=\frac{3H^{4}}{32\pi^{2}}\left(\sinh(H|t-t^{\prime}|\right)^{-4},$
(B.5)

where we have pointed out that these are vacuum fluctuations, to distinguish
them from the statistically averaged two point functions. To compute this, let
us note that after smearing the field the noise function is directly
proportional to the smeared speed, as can be seen after taking the time
derivative from (A.6). Then, we have that the correlation function is given by

$\displaystyle\langle\delta\dot{\phi}(\vec{x},t)\delta\dot{\phi}(\vec{x},t^{\prime})\rangle_{\mathrm{av}}$
$\displaystyle=\frac{H^{4}}{8\pi^{2}}\delta(t-t^{\prime}).$ (B.6)

Notice that, while at separate time both decay to zero, at equal times the
correlation function diverges while the statistical average is finite. This is
an effect of smearing out over a region where there are statistical
fluctuations, which in the end translates into the correlation function for
the speed being finite at equal time.

##  Appendix C General solution for linear coefficients

In this appendix we will follow [54] to derive general solutions of the
Fokker-Planck equation. Let us start by considering the following equation,

$\frac{\partial P}{\partial
t}=-A_{ij}\frac{\partial}{\partial\phi_{i}}(\phi^{j}P)+\frac{1}{2}D_{ij}\frac{\partial^{2}}{\partial\phi_{i}\phi_{j}}P\
,$ (C.1)

where both $A_{ij}$ and $D_{ij}$ are $n\times n$ constant matrices and in
addition $D_{ij}$ is symmetric and semipositive definite. Subject to initial
conditions

$P(\phi,0)=\prod_{i}\delta(\phi^{i}-\phi^{i}_{0}).$ (C.2)

The solution of this equation is Gaussian which we will show. First if we
multiply the equation for $\phi^{i}$ and integrate over $\phi$ we find after
integration by parts,

$\partial_{t}\langle\phi^{i}\rangle=A_{kj}\langle\phi^{j}\rangle\ ,$ (C.3)

whose solution is given by

$\langle\phi^{i}\rangle=e^{tA}y_{0}\ ,$ (C.4)

in matrix notation. Now if we insert $\phi^{i}\phi^{j}$ into the Fokker-Planck
equation we find,

$\partial_{t}\langle\phi^{i}\phi^{j}\rangle=A^{i}_{k}\langle\phi^{k}\phi^{j}\rangle+A^{j}_{k}\langle\phi^{k}\phi^{i}\rangle+D^{ij}\
.$ (C.5)

It is more convenient to use covariance matrix
$C_{ij}=\langle\phi_{i}\phi_{j}\rangle-\langle\phi_{i}\rangle\langle\phi_{j}\rangle$.
In which case the above equation reduces to,

$\partial_{t}C=AC+CA^{t}+D\ ,$ (C.6)

in matrix notation. Now if we write the covariance matrix as
$C=e^{tA}\bar{C}e^{tA^{t}}$ then we get

$\partial_{t}\bar{C}=e^{-tA}De^{-tA^{t}}\ ,$ (C.7)

where assuming that $\bar{C}(0)=0$ has a solution given by,

$C(t)=\int_{0}^{t}e^{(t-t^{\prime})A}De^{(t-t^{\prime})A^{t}}dt^{\prime}\ .$
(C.8)

Even though this expression looks abstract it can be easily computed . To
conclude given that we assume that the distribution was Gaussian then it is
fully determined by the covariance matrix, hence we have,

$P(\phi,t)=\frac{1}{(2\pi)^{n/2}}\frac{1}{\sqrt{\det
C}}\exp\left(\frac{1}{2}(\phi-\langle\phi\rangle)^{t}\rangle
C^{-1}(\phi-\langle\phi\rangle)\right)\ .$ (C.9)

Let us check now that the expression for $C$ agrees with what we found before.
We first have that,

$\displaystyle A=\begin{pmatrix}0&1\\\ 0&-3H\end{pmatrix},\qquad
D=\begin{pmatrix}0&0\\\ 0&D_{\zeta}\end{pmatrix}.$ (C.10)

The eigenvalues of $A$ are $0$ and $3H$, and we can compute the exponential of
$A$

$\displaystyle\exp(-tA)=\left(\begin{array}[]{cc}1&\frac{1-e^{-3Ht}}{3H}\\\
0&e^{-3Ht}\\\ \end{array}\right),$ (C.13)

We then have that,

$\displaystyle
e^{(t-t^{\prime})A}De^{(t-t^{\prime})A^{t}}=\left(\begin{array}[]{cc}\frac{D\left(e^{3H\left(t^{\prime}-t\right)}-1\right)^{2}}{9H^{2}}&\frac{De^{6H\left(t^{\prime}-t\right)}\left(e^{3H\left(t-t^{\prime}\right)}-1\right)}{3H}\\\
\frac{De^{6H\left(t^{\prime}-t\right)}\left(e^{3H\left(t-t^{\prime}\right)}-1\right)}{3H}&De^{6H\left(t^{\prime}-t\right)}\\\
\end{array}\right).$ (C.16)

Performing the integral we find ,

$\displaystyle
C(t)=\left(\begin{array}[]{cc}-\frac{D\left(-6Ht+e^{-6Ht}-4e^{-3Ht}+3\right)}{54H^{3}}&\frac{De^{-6Ht}\left(e^{3Ht}-1\right)^{2}}{18H^{2}}\\\
\frac{De^{-6Ht}\left(e^{3Ht}-1\right)^{2}}{18H^{2}}&-\frac{D\left(e^{-6Ht}-1\right)}{6H}\\\
\end{array}\right),$ (C.19)

which coincides with the expression we found before. Taking the limit $t\gg
1/H$ we have

$C(t)\to\left(\begin{array}[]{cc}\frac{Dt}{9H^{2}}&\frac{D}{18H^{2}}\\\
\frac{D}{18H^{2}}&\frac{D}{6H}\\\ \end{array}\right),$ (C.20)

so we find that

$P=\frac{1}{2\pi}\sqrt{\frac{54H^{3}}{D_{\zeta}^{2}t}}\exp\left(-\frac{9H^{2}}{2D_{\zeta}t}\zeta^{2}+\frac{3H}{D_{\zeta}t}v_{\zeta}\zeta-\frac{3H}{D_{\zeta}}v_{\zeta}^{2}\right)\
.$ (C.21)

##  Appendix D Ultralight field

In this appendix we will show how to modify the noise term in the Fokker-
Planck equation to take into account the superhorizon time dependence of the
curvature power spectrum. This can be done in general but in this case we will
focus on the case of an ultralight field [49]. Let us start from the Langevin
Eqs. (3.6), specialised to the case of $\mu^{2}=0$. For convenience we will
define a new variable $\tilde{\zeta}$ such that the Langevin equations are
now,

$\displaystyle\dot{\tilde{\zeta}}$ $\displaystyle=v_{\zeta},$
$\displaystyle\dot{v}_{\zeta}$
$\displaystyle=-3Hv_{\zeta}+\eta_{\tilde{\zeta}}(t),$
$\displaystyle\dot{\psi}$ $\displaystyle=v_{\psi},$ $\displaystyle\dot{\psi}$
$\displaystyle=-3Hv_{\psi}+2\Omega f_{\zeta}v_{\zeta}+\eta_{\psi}(t).$ (D.1)

Notice that in this case $\tilde{\zeta}$ grows outside the horizon due to the
interactions with $\psi$. This implies that the noise term
$\eta_{\tilde{\zeta}}(t)$, which is computed by coarse graining
$\tilde{\zeta}$, should also grow on superhorizon scales. To do so, let us
first recall that for an ultralight field, the two point function for the
curvature mode is given by [49]. ,

$\langle\zeta^{2}\rangle=\frac{1}{2\epsilon
a^{2}}\frac{1}{2k^{3}\tau^{2}}\left(1+\lambda^{2}\left[A_{1}-A_{2}\log(-k\tau)+\log^{2}(-k\tau)\right]\right)\
,$ (D.2)

where $\lambda\equiv\frac{2\Omega}{H}$, and $A_{1}$, $A_{2}$ are given by,

$\displaystyle A_{1}$ $\displaystyle=-\frac{\pi^{2}}{6}+(3-\ln 2)(1-\ln
2)-\gamma_{E}(4-\gamma-2\log 2)\simeq-2.11\ ,$ $\displaystyle A_{2}$
$\displaystyle=4-2\gamma-2\ln 2\simeq 1.46\ .$ (D.3)

Notice that it diverges in the limit $k\tau\to 0$. However, since inflation
will last for a finite amount of time, then there is an natural cut-off for
the power spectrum. A more systematic way to deal with this IR behaviour is to
regularise the growing logs by introducing boundary counterterms as in [76].
Doing so results in the following regularised expression,

$\langle\zeta^{2}\rangle=\frac{1}{2\epsilon
a^{2}}\frac{1}{2k^{3}\tau^{2}}\left(1+\lambda^{2}\left[A_{1}-A_{2}\log(k/\mu_{IR})+\log^{2}(k/\mu_{IR})\right]\right),$
(D.4)

where $\mu_{IR}$ is an infrared cut-off. If we set $\mu_{IR}=\epsilon H$ the
diffusion coefficient for $\zeta$ changes to,

$\displaystyle D_{\zeta}\to
D_{\zeta}\left(1+\lambda^{2}\left[A_{1}-A_{2}\log(a)+\log^{2}(a)\right]\right).$
(D.5)

To compute the covariance matrix we can use the results from Appendix C, which
are compatible with a time dependent difussion matrix. At leading order in
$\lambda^{2}\Delta N^{2}$, and for $t\gg 1/H$, we have that,

$\displaystyle
C=\left(\begin{array}[]{cccc}\frac{Dt(3+H^{2}t^{2}\lambda^{2})}{27H^{2}}&-\frac{\sqrt{2\epsilon}Dt(3+H^{2}t^{2}\lambda^{2})}{81H^{2}}&\frac{D}{18H^{2}}&0\\\
-\frac{\sqrt{2\epsilon}Dt(3+H^{2}t^{2}\lambda^{2})}{81H^{2}}&\frac{H^{3}t}{4\pi^{2}}&0&\frac{H^{3}}{8\pi^{2}}\\\
\frac{D}{18H^{2}}&0&\frac{D}{6H}&0\\\
0&\frac{H^{3}}{8\pi^{2}}&0&\frac{3H^{4}}{8\pi^{2}}\\\ \end{array}\right)\ ,$
(D.10)

where we have kept terms at order $\sqrt{\epsilon}\lambda\Delta N^{2}$. Notice
that this steady state is reached within a couple of efolds. If we drop the
slow roll terms, the PDF is given by

$P=\frac{6\sqrt{3}}{\sqrt{D_{\zeta}^{2}H^{6}t^{4}\lambda}}\exp\left(-\frac{27H^{2}}{2D_{\zeta}t^{3}\lambda^{2}}\tilde{\zeta}^{2}+\frac{9}{D_{\zeta}Ht^{3}\lambda^{2}}\tilde{\zeta}v_{\zeta}-\frac{3H}{D}v_{\zeta}^{2}-\frac{2\pi^{2}}{H^{3}t}\psi^{2}+\frac{4\pi^{2}}{3H^{4}t}v_{\psi}\psi-\frac{4\pi^{2}}{3H^{4}}v_{\psi}^{2}\right)\
,$ (D.11)

It is possible to compute the correlation function directly from the
distribution by integrating over all the fields

$\displaystyle\langle\phi^{a}\phi^{b}\rangle=\int\prod_{i}D\phi^{i}\phi^{a}\phi^{b}P\
,$ (D.12)

doing so we get,

$\displaystyle\langle\tilde{\zeta}^{2}\rangle=\frac{D_{\zeta}t^{3}\lambda^{2}}{27},\qquad\langle\psi^{2}\rangle=\frac{H^{3}t}{4\pi^{2}}\
.$ (D.13)

These are real space correlation functions, to compare them with the power
spectrum we have to use the relation,

$\langle\phi^{2}\rangle=\int d\log k\Delta_{\phi}^{2}.$ (D.14)

Taking derivatives on both sides implies that,

$\frac{d}{d\log k_{*}}(\langle\phi^{2}\rangle)=\Delta_{\phi}^{2},$ (D.15)

where $k_{*}$ is the horizon crossing wavenumber. Since for modes that have
crossed the horizon, $t$ can be written as $t_{*}=\frac{1}{H}\log(k_{*}/H)$,
we have that

$\Delta_{\tilde{\zeta}}^{2}=\frac{D_{\zeta}\Delta
N^{2}\lambda^{2}}{9H}=\frac{H^{2}}{4\pi^{2}}\alpha^{2}\Delta
N^{2},\qquad\Delta_{\psi}^{2}=\frac{H^{2}}{4\pi^{2}}\ ,$ (D.16)

where we have used that $\frac{1}{H}\log(k_{*}/H)=\Delta N$, is the number of
efolds until the end of inflation and where
$D_{\zeta}=9H^{5}/(8\epsilon\pi^{2})$ and $\lambda=\sqrt{2\epsilon\alpha}/H$.
This result coincide with the power spectrum computed in [49]. Finally let us
notice that the faster growth in the variance avoids inflation becoming
eternal. This can be seen by the following argument. In general inflation
becomes eternal if during an interval $t\sim H{-1}$ the quantum fluctuations
$\langle\delta\phi^{2}\rangle^{1/2}$ is larger than the classical change of
the field $\Delta\phi=\dot{\phi}/H$. From (D.13) we have that this is avoid if

$\displaystyle\frac{\dot{\phi}}{H}\sqrt{\frac{D_{\zeta}\lambda^{2}}{27H^{3}}}=\frac{H\lambda}{\sqrt{12}\pi}\leq\frac{\dot{\phi}}{H}\
.$ (D.17)

where we have used that at horizon crossing
$\zeta=-\frac{H}{\dot{\phi}}\delta\phi$. The last inequality implies that the
condition for eternal inflation is more strict than in single fields inflation
(which is that $\sqrt{\Delta_{\zeta}^{2}}\leq 1$). Indeed (D.17) can be
written as,

$\displaystyle\frac{\sqrt{\Delta^{2}_{\zeta}}}{\sqrt{3}\Delta{N}}\ll\sqrt{\Delta_{\zeta}^{2}}\leq
1\ .$ (D.18)

Notice that this assumes that the power spectrum didn’t vary significantly
during the whole inflation. If there is a momentarily increase of the power
spectrum, such that $\lambda\gg 1$ then the last inequality will not hold.

##  Appendix E Computing further corrections to the PDF.

In this appendix we will include higher order corrections to the solution of
the Fokker-Planck equation in (4.31). First, let us solve the following,

$\displaystyle\frac{dP}{dt}$
$\displaystyle=\frac{\partial}{\partial\psi}\left(t_{\psi}^{-1}\psi
P+\frac{D_{\psi}}{2}\frac{\partial
P}{\partial\psi}\right)+\frac{H\Delta_{\zeta}^{2}}{2}\frac{\partial^{2}P}{\partial\zeta^{2}}$
$\displaystyle+H^{2}\Delta_{\zeta}^{2}\frac{\partial^{2}}{\partial_{\zeta}\partial\psi}\left(\left(\frac{2\Omega^{2}}{H^{2}}\psi+\frac{2f\Omega}{3H}\right)P\right)+\frac{6\Omega^{2}}{f_{\zeta}^{2}H}\frac{\partial}{\partial\zeta}(\psi^{2}P)\
,$ (E.1)

Our task will be to add the last term to the PDF (4.31). To simplify we will
look for late time solutions such as $\psi$ has reached its equilibrium
distribution. We can eliminate some of the terms in the equation by Fourier
transforming $\zeta$ to $p$, and look for solutions of the form,

$\displaystyle P(p,\psi,t)\sim\exp(-\sigma_{\zeta}^{2}p^{2}/2)F(p,\psi)$ (E.2)

After replacing into the Fokker-Planck equation (E.1) we obtain,

$\displaystyle 0=$ $\displaystyle
F^{\prime\prime}(k)+\frac{f_{\zeta}^{2}H}{3D_{\psi}t_{\psi}}(6\psi+4ifHt_{\psi}\Delta_{\zeta}^{2}\Omega
p+4it\psi\Delta_{\zeta}^{2}\Omega^{2}\psi p)F^{\prime}(\psi)$
$\displaystyle+\frac{1}{3D_{\psi}t_{\psi}}(36ipt_{\psi}\psi^{2}\Omega^{2}+2f^{2}H(3+2ipt_{\psi}\Delta_{\zeta}^{2}\Omega^{2}))F(\psi).$
(E.3)

Ignoring the linear mixing term the solution is given by,

$\displaystyle F(p,\psi)$
$\displaystyle\propto\exp\left(-\frac{\psi^{2}}{2\sigma_{\psi}^{2}}\left(1+\frac{2ipt_{\psi}\Omega^{2}\Delta_{\zeta}^{2}}{3}-\frac{i}{2}\sqrt{f(p)}\right)\right)\times
H_{n}\left(\sqrt{\frac{\psi^{2}}{2\sigma_{\psi}^{2}}}\left(-\frac{f(p)}{4}\right)^{1/4}\right).$
(E.4) $\displaystyle n$
$\displaystyle=\frac{1}{2}-\frac{-3i+2pt_{\psi}\Delta_{\zeta}^{2}\Omega^{2}}{\sqrt{f(p)}},$
(E.5) $\displaystyle f(p,q)$
$\displaystyle=108\frac{t_{\psi}\Omega^{2}\sigma_{\psi}^{2}}{f_{\zeta}^{2}H}p+(3i-2pt_{\psi}\Delta_{\zeta}^{2}\Omega^{2})^{2},$
(E.6)

where we have discarded the second solution since it grows for large
$\sqrt{p}\psi$. In order to obtain a simplified expression let us notice that
a typical fluctuation of $p\sim 1/\sqrt{\sigma_{\zeta}^{2}}$. Using this we
can deduce that the at leading order $f(p)$ is constant, as the ratio between
the two leading order terms is given by,

$\displaystyle
108t_{\psi}\frac{\Omega^{2}\sigma_{\psi}^{2}}{f_{\zeta}H}\frac{1}{9}=\frac{12t_{\psi}^{2}\Omega^{2}\Delta_{\zeta}^{2}}{\sqrt{\sigma_{\zeta}^{2}}}\ll
1.$ (E.7)

This inequality still holds for larger values of the $\zeta$. If we write
$f_{(}p)$ at leading order we find that the Hermite function reduces to one
and we recover the usual distribution (4.32). When adding the linear mixing
term the distribution is more complicated but still depends on $f(p)$. If we
ignore this term, we find that at leading order the PDF contains further
corrections at order $p^{2}$. The effect of those add up to the quadratic
terms that appeared in the drift for $v_{\zeta}$. In the end this will modify
the value of the tail for very large values of $\zeta$, acting as
exponetentially suppressed corrections as expected.

## References

  * [1] Planck Collaboration, Y. Akrami et al., “Planck 2018 results. IX. Constraints on primordial non-Gaussianity,” Astron. Astrophys. 641 (2020) A9, arXiv:1905.05697 [astro-ph.CO].
  * [2] J. M. Maldacena, “Non-Gaussian features of primordial fluctuations in single field inflationary models,” JHEP 05 (2003) 013, arXiv:astro-ph/0210603.
  * [3] C. Cheung, P. Creminelli, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, “The Effective Field Theory of Inflation,” JHEP 03 (2008) 014, arXiv:0709.0293 [hep-th].
  * [4] X. Chen and Y. Wang, “Quasi-Single Field Inflation and Non-Gaussianities,” JCAP 04 (2010) 027, arXiv:0911.3380 [hep-th].
  * [5] X. Chen and Y. Wang, “Large non-Gaussianities with Intermediate Shapes from Quasi-Single Field Inflation,” Phys. Rev. D 81 (2010) 063511, arXiv:0909.0496 [astro-ph.CO].
  * [6] A. Achucarro, J.-O. Gong, S. Hardeman, G. A. Palma, and S. P. Patil, “Features of heavy physics in the CMB power spectrum,” JCAP 01 (2011) 030, arXiv:1010.3693 [hep-ph].
  * [7] A. Achucarro, V. Atal, S. Cespedes, J.-O. Gong, G. A. Palma, and S. P. Patil, “Heavy fields, reduced speeds of sound and decoupling during inflation,” Phys. Rev. D 86 (2012) 121301, arXiv:1205.0710 [hep-th].
  * [8] N. Arkani-Hamed and J. Maldacena, “Cosmological Collider Physics,” arXiv:1503.08043 [hep-th].
  * [9] H. Lee, D. Baumann, and G. L. Pimentel, “Non-Gaussianity as a Particle Detector,” JHEP 12 (2016) 040, arXiv:1607.03735 [hep-th].
  * [10] A. Achúcarro, S. Céspedes, A.-C. Davis, and G. A. Palma, “Constraints on Holographic Multifield Inflation and Models Based on the Hamilton-Jacobi Formalism,” Phys. Rev. Lett. 122 (2019) no. 19, 191301, arXiv:1809.05341 [hep-th].
  * [11] N. Arkani-Hamed, D. Baumann, H. Lee, and G. L. Pimentel, “The Cosmological Bootstrap: Inflationary Correlators from Symmetries and Singularities,” JHEP 04 (2020) 105, arXiv:1811.00024 [hep-th].
  * [12] M. Celoria, P. Creminelli, G. Tambalo, and V. Yingcharoenrat, “Beyond perturbation theory in inflation,” JCAP 06 (2021) 051, arXiv:2103.09244 [hep-th].
  * [13] R. Flauger, M. Mirbabayi, L. Senatore, and E. Silverstein, “Productive Interactions: heavy particles and non-Gaussianity,” JCAP 10 (2017) 058, arXiv:1606.00513 [hep-th].
  * [14] X. Chen, G. A. Palma, B. Scheihing Hitschfeld, and S. Sypsas, “Reconstructing the Inflationary Landscape with Cosmological Data,” Phys. Rev. Lett. 121 (2018) no. 16, 161302, arXiv:1806.05202 [astro-ph.CO].
  * [15] X. Chen, G. A. Palma, W. Riquelme, B. Scheihing Hitschfeld, and S. Sypsas, “Landscape tomography through primordial non-Gaussianity,” Phys. Rev. D 98 (2018) no. 8, 083528, arXiv:1804.07315 [hep-th].
  * [16] G. A. Palma, B. Scheihing Hitschfeld, and S. Sypsas, “Non-Gaussian CMB and LSS statistics beyond polyspectra,” JCAP 02 (2020) 027, arXiv:1907.05332 [astro-ph.CO].
  * [17] S. Hooshangi, M. H. Namjoo, and M. Noorbala, “Rare Events Are Nonperturbative: Primordial Black Holes From Heavy-Tailed Distributions,” arXiv:2112.04520 [astro-ph.CO].
  * [18] Y.-F. Cai, X.-H. Ma, M. Sasaki, D.-G. Wang, and Z. Zhou, “One Small Step for an Inflaton, One Giant Leap for Inflation: a novel non-Gaussian tail and primordial black holes,” arXiv:2112.13836 [astro-ph.CO].
  * [19] B. Carr and F. Kuhnel, “Primordial Black Holes as Dark Matter: Recent Developments,” Ann. Rev. Nucl. Part. Sci. 70 (2020) 355–394, arXiv:2006.02838 [astro-ph.CO].
  * [20] A. M. Green and B. J. Kavanagh, “Primordial Black Holes as a dark matter candidate,” J. Phys. G 48 (2021) no. 4, 043001, arXiv:2007.10722 [astro-ph.CO].
  * [21] G. Franciolini, A. Kehagias, S. Matarrese, and A. Riotto, “Primordial Black Holes from Inflation and non-Gaussianity,” JCAP 03 (2018) 016, arXiv:1801.09415 [astro-ph.CO].
  * [22] V. Atal and C. Germani, “The role of non-gaussianities in Primordial Black Hole formation,” Phys. Dark Univ. 24 (2019) 100275, arXiv:1811.07857 [astro-ph.CO].
  * [23] I. Musco, V. De Luca, G. Franciolini, and A. Riotto, “Threshold for primordial black holes. II. A simple analytic prescription,” Phys. Rev. D 103 (2021) no. 6, 063538, arXiv:2011.03014 [astro-ph.CO].