singularities in the complex plane which preclude the rotation to Euclidean
time: whether this leads to specific features for the WFU is worth exploring.
The UV completion of the type of models studied in this work contains
additional states: it would be interesting to understand when and how they
modify the predictions for the tail of the wavefunction. Eventually, it would
be valuable to relate our non-perturbative computations of the wavefunction of
the universe, which fully characterizes the state of primordial fluctuations,
to the ongoing developments about how to best extract information from
cosmological data, see for example [30] and references therein.

## Acknowledgements

It is a pleasure to thank C. Duaso Pueyo, M. Mirbabayi, E. Pajer, M. Peloso,
L. Senatore, M. Serone, V. Vennin and G. Villadoro for illuminating
discussions. SRP is supported by the European Research Council under the
European Union’s Horizon 2020 research and innovation programme (grant
agreement No 758792, Starting Grant project GEODESI). V.Y. is supported by the
WPI Research Center Initiative, MEXT, Japan. The work of V.Y. has been
supported by JSPS KAKENHI grant JP22K20367.

##  Appendix A Derivative versus polynomial form of the action

In this appendix we explicitly show that our action (2.7) indeed agrees with
the action (24) of [13], and we demonstrate the practical advantage of using
our form, even for the computation of $n$-point correlation functions. In
doing so, it is important to remember that temporal boundary terms in the
action contribute in general to correlation functions [31, 32, 33, 34]. In
particular, total (temporal) derivative terms depending on $\dot{\pi}$ (by
contrast to the ones involving $\pi$ only) contribute to correlation functions
of $\pi$, which are our object of interest. Hence, we will only keep such
boundary terms in the following discussion.

###  A.1 Action of $\pi$ in two different forms

Our starting point here is the action (2.7):

$\displaystyle S=\int\,\mathrm{d}^{4}x\,a^{3}M_{\rm
Pl}^{2}\dot{H}(t+\pi)(\partial_{\mu}\pi)^{2}\;.$ (A.1)

Expanding $\dot{H}(t+\pi)$ as a power series gives

$\displaystyle S=M_{\rm
Pl}^{2}\sum_{n=0}^{\infty}\int\,\mathrm{d}^{4}x\,\frac{a^{3}}{n!}H^{(n+1)}(t)~{}\pi^{n}\partial_{\mu}\pi\partial^{\mu}\pi\;,$
(A.2)

where $H^{(n)}(t)$ denotes the $n^{\rm th}$ derivative with respect to time
acting on $H(t)$. Here we are interested in the terms with $n\geq 1$. We
rewrite the term $\pi^{n}\partial_{\mu}\pi$ as
$\partial_{\mu}(\pi^{n+1})/(n+1)$, so that the interacting part of the action
above becomes

$\displaystyle S_{\textrm{int}}$ $\displaystyle=M_{\rm Pl}^{2}\sum_{n\geq
1}^{\infty}\int\,\mathrm{d}^{4}x\,\frac{a^{3}}{(n+1)!}H^{(n+1)}(t)\,\partial_{\mu}\pi^{n+1}\partial^{\mu}\pi\;,$
$\displaystyle=-M_{\rm Pl}^{2}\sum_{n\geq
1}^{\infty}\int\,\mathrm{d}^{4}x\,\frac{1}{(n+1)!}\partial_{\mu}\bigg{[}a^{3}H^{(n+1)}(t)\,\partial^{\mu}\pi\bigg{]}\pi^{n+1}+S_{\textrm{boundary}}\;,$
(A.3)

where we have performed an integration by parts in the second line, and we
kept the relevant boundary term

$\displaystyle S_{\textrm{boundary}}=-M_{\rm Pl}^{2}\sum_{n\geq
1}^{\infty}\int\,\mathrm{d}^{4}x\,\partial_{t}\left[\frac{a^{3}}{(n+1)!}H^{(n+1)}(t)\,\pi^{n+1}\dot{\pi}\right]\,.$
(A.4)

We define $\Box\pi\equiv\partial_{\mu}(a^{3}\partial^{\mu}\pi)/a^{3}$. The
action (A.3) then becomes

$\displaystyle S_{\textrm{int}}=-M_{\rm Pl}^{2}\sum_{n\geq
1}^{\infty}\int\,\mathrm{d}^{4}x\,\frac{1}{(n+1)!}\bigg{[}-a^{3}H^{(n+2)}(t)\,\pi^{n+1}\dot{\pi}+a^{3}H^{(n+1)}(t)\,\pi^{n+1}\Box\pi\bigg{]}+S_{\textrm{boundary}}\;.$
(A.5)

Performing an integration by parts on the first term on the RHS of the action
above, we obtain

$\displaystyle S_{\textrm{int}}=-M_{\rm
Pl}^{2}\sum_{n>2}^{\infty}\int\,\mathrm{d}^{4}x\,a^{3}\bigg{\\{}\frac{1}{n!}\bigg{[}3H(t)H^{(n)}(t)+H^{(n+1)}(t)\bigg{]}\pi^{n}+\frac{1}{(n-1)!}H^{(n-1)}(t)\,\pi^{n-1}\Box\pi\bigg{\\}}+S_{\textrm{boundary}}\;.$
(A.6)

where we have changed $n\rightarrow n-2$. Additionally, note that the equation
of motion derived from the quadratic Lagrangian for $\pi$ reads

$\displaystyle\frac{\delta\mathcal{L}_{2}}{\delta\pi}=-2M_{\rm
Pl}^{2}a^{3}\dot{H}(t)\Box\pi+2M_{\rm Pl}^{2}a^{3}\ddot{H}(t)\dot{\pi}\,.$
(A.7)

Hence, to work out the action at first order in $b$, one can neglect the last
term in (A.7) when inserted back in (A.6), which gives

$\displaystyle
S_{\textrm{int}}=\sum_{n>2}^{\infty}\int\,\mathrm{d}^{4}x\,\bigg{\\{}-\frac{a^{3}M_{\rm
Pl}^{2}}{n!}\bigg{[}3H(t)H^{(n)}(t)+H^{(n+1)}(t)\bigg{]}\pi^{n}+f_{n}(\pi)\frac{\delta\mathcal{L}_{2}}{\delta\pi}\bigg{\\}}+\mathcal{O}(b^{2})+S_{\textrm{boundary}}\;,$
(A.8)

where we have defined

$\displaystyle
f_{n}(\pi)\equiv\frac{1}{2(n-1)!}\frac{H^{(n-1)}(t)}{\dot{H}(t)}\,\pi^{n-1}\;.$
(A.9)

We see that the action (A.8) is the same as the action (24) of [13] 222222One
can straightforwardly show that the action (A.8) is also equivalent to the
$n^{\rm th}$-order interacting action obtained in [22], as shown in [13].,
although the authors here did not write the boundary terms, despite their role
in the computation of $n$-point functions, as we will see in (A.20).
Therefore, our non-perturbative form of the action (2.7) is equivalent to
their action.

Although what ultimately matters, in perturbative computations, are
correlation functions of $\pi$, it is interesting to study the structure of
the equation of motion of $\pi$. In particular, the fact that the action (A.8)
involves self-interactions of $\pi$ without derivatives acting on it seems to
suggest that a constant $\pi$ is not a solution to the equation of motion,
implying that $\pi$ is not conserved on super-horizon scales. However, this is
not the case due to the contributions coming from the term in
$\delta\mathcal{L}_{2}/\delta\pi$ to the equation of motion of $\pi$: one can
verify that at each order in $\pi$ the equation of motion derived from (A.8)
admits a constant solution. In contrast, this property is manifestly valid,
and non-perturbatively, from our form (2.7).

###  A.2 Time-independence of the bispectrum

It is instructive to see explicitly the late-time constancy of the 3-point
function of $\pi$ computed from the form (A.8), and to contrast it with the
computation starting from our form (2.7) of the action. We will see that the
constancy of the bispectrum is immediate in the latter case, whereas in the
former, it requires taking into subleading terms in the expansion in large
$\alpha$ as well as cancellations between the contributions from the bulk and
the boundary terms. This provides a non-trivial consistency check of our
computations, and showcases the usefulness of our form of the action.

##### Polynomial form of the action.

There are three types of terms in the Lagrangian (A.8): the bulk terms that
are polynomial in $\pi$, the terms proportional to the linear equation of
motion, and boundary terms. When computing correlation functions in
perturbation theory, terms proportional to the linear equation of motion never
contribute, as the interaction-picture fields are by definition free fields
obeying that equation. We treat the two other types of terms in turn.

The contribution to the 3-point function of $\pi$ at time $\eta$ from the bulk
terms in (A.8), reads

$\displaystyle\langle\pi(\boldsymbol{k}_{1},\eta)\pi(\boldsymbol{k}_{2},\eta)\pi(\boldsymbol{k}_{3},\eta))\rangle^{\prime}_{\textrm{bulk}}$
$\displaystyle=\frac{1}{(4\epsilon_{\star}M_{\rm
Pl}^{2})^{2}2\epsilon_{\star}H^{4}\prod_{i}k_{i}^{3}}\times$
$\displaystyle\hskip
14.22636pt\times\textrm{Im}\left[\int_{-\infty(1-i\epsilon)}^{\eta}\frac{\,\mathrm{d}\eta^{\prime}}{\eta^{\prime
4}}\left(H^{(4)}+3HH^{(3)}\right)(\eta^{\prime})f(\eta^{\prime})f^{*}(\eta)\right]\;,$
(A.10)

with $f(\eta)=e^{ik_{\rm t}\eta}\prod_{i}(1-ik_{i}\eta)$ and $k_{\rm
t}=\sum_{i}k_{i}$. The prime means that we drop the factor
$(2\pi)^{3}\delta^{(3)}(\boldsymbol{k}_{1}+\boldsymbol{k}_{2}+\boldsymbol{k}_{3})$.
All the integrals here can be done analytically, but in terms of incomplete
Gamma functions which are not very illuminating for our purposes. Instead,
remember that our goal is simply to show that the bispectrum of $\pi$ goes to
a constant at late times. Hence, we will only keep track of the relevant
contributions at the asymptotic future, which are simply oscillations in
$\eta^{\pm i\alpha}$, and we check that they cancel in the final result. For
that purpose, it is sufficient to expand $f(\eta)$ at late times:

$\displaystyle
f(\eta)=1+\frac{\eta^{2}}{2}\sum_{i}k_{i}^{2}+\frac{i\eta^{3}}{3}\sum_{i}k_{i}^{3}+{\cal
O}(\eta^{4})\;,$ (A.11)

where we have to go up to third order as one needs

$\displaystyle\textrm{Im}\left[f(\eta^{\prime})f^{*}(\eta)\right]=\frac{1}{3}(\eta^{\prime
3}-\eta^{3})\sum_{i}k_{i}^{3}+\textrm{subleading}\;.$ (A.12)

Here, we use the $i\epsilon$ prescription that is relevant only in the
asymptotic past, so that, for our purposes, taking the imaginary part in
(A.12) effectively applies only to $f(\eta^{\prime})f^{*}(\eta)$. As for the
terms denoted as subleading, they do not lead to contributions that survive at
late times. At first order in $b$, one can consider $H$ constant in the term
$HH^{(3)}$, and with

$\displaystyle\dot{H}(\eta)=-\epsilon_{\star}H_{\star}^{2}\left[1-\frac{6b}{\alpha}\cos\left(\alpha\log\left(\eta/\eta_{\star}\right)-\tilde{\delta}\right)\right]\,,$
(A.13)

one obtains that the terms that survive at late times are

$\displaystyle\langle\pi(\boldsymbol{k}_{1},\eta)$
$\displaystyle\pi(\boldsymbol{k}_{2},\eta)\pi(\boldsymbol{k}_{3},\eta))\rangle^{\prime}_{\textrm{bulk}}\supset-\frac{H_{\star}b\,\alpha^{2}}{(4\epsilon_{\star}M_{\rm
Pl}^{2})^{2}}\frac{\sum_{i}k_{i}^{3}}{\prod_{i}k_{i}^{3}}\times$
$\displaystyle\times\int^{\eta}_{-\infty(1-i\epsilon)}\frac{\,\mathrm{d}\eta^{\prime}}{\eta^{\prime
4}}\left(\eta^{\prime
3}-\eta^{3}\right)\left[\sin\left(\alpha\log\left(\eta^{\prime}/\eta_{\star}\right)-\tilde{\delta}\right)+\frac{3}{\alpha}\cos\left(\alpha\log\left(\eta^{\prime}/\eta_{\star}\right)-\tilde{\delta}\right)\right]$
$\displaystyle=-\frac{3H_{\star}b}{(4\epsilon_{\star}M_{\rm
Pl}^{2})^{2}}\frac{\sum_{i}k_{i}^{3}}{\prod_{i}k_{i}^{3}}\times\sin\left(\alpha\log\left(\eta/\eta_{\star}\right)-\tilde{\delta}\right)\,.$
(A.14)

Note that cancellations between the various contributions—the sine and the
cosine terms, as well as the ones with and without $\eta$-dependence in the
integrand—lead to the amplitude of the result being independent of $\alpha$.
In particular, it was important to take into account the cosine term, despite
the fact that its contribution to the integrand is subdominant in the limit of
large $\alpha$.

Let us now consider the boundary term

$S^{(3)}_{\textrm{boundary}}=-\frac{M_{\rm
Pl}^{2}}{2}\int\,\mathrm{d}t\,\mathrm{d}^{3}x\,\partial_{t}\left(a^{3}\ddot{H}(t)\pi^{2}\dot{\pi}\right)\,.$
(A.15)

The contribution from any such term can be worked out from first principles
using the commutation relation (see e.g. Sec. 3.3 of [34]). Using that the
(linear) conjugate momentum to $\pi$ is $\partial{\cal
L}_{2}/\partial\dot{\pi}=-2a^{3}M_{\rm Pl}^{2}\dot{H}\dot{\pi}$, this gives
here:

$\displaystyle\langle\pi(\boldsymbol{k}_{1},\eta)$
$\displaystyle\pi(\boldsymbol{k}_{2},\eta)\pi(\boldsymbol{k}_{3},\eta))\rangle^{\prime}_{\textrm{boundary}}=-\left(\frac{\ddot{H}}{2\dot{H}}\right)(\eta)\left(P_{\pi}(k_{1},\eta)P_{\pi}(k_{2},\eta)+\textrm{
2 perms.}\right)\;,$ (A.16)

where, at first order in $b$, one can use the standard power spectrum for
$\pi$, $P_{\pi}(k,\eta)=(1+k^{2}\eta^{2})/(4\epsilon_{\star}M_{\rm
Pl}^{2}k^{3})$. Keeping the terms that survive at late times, this gives

$\displaystyle\langle\pi(\boldsymbol{k}_{1},\eta)$
$\displaystyle\pi(\boldsymbol{k}_{2},\eta)\pi(\boldsymbol{k}_{3},\eta))\rangle^{\prime}_{\textrm{boundary}}=\frac{3H_{\star}b}{(4\epsilon_{\star}M_{\rm
Pl}^{2})^{2}}\frac{\sum_{i}k_{i}^{3}}{\prod_{i}k_{i}^{3}}\times\sin\left(\alpha\log\left(\eta/\eta_{\star}\right)-\tilde{\delta}\right)\,,$
(A.17)

which precisely cancel with (A.14). Therefore, we have explicitly shown that
the $3$-point function of $\pi$ goes to a constant at late times, using the
polynomial form of the action accompanied with the necessary boundary terms.

##### Derivative form of the action:

When using the derivative form of the action (2.7), the cubic part reads

$S^{(3)}=\int\,\mathrm{d}^{4}x\,a^{3}M_{\rm
Pl}^{2}\ddot{H}\pi(-\dot{\pi}^{2}+(\partial_{i}\pi)^{2}/a^{2})\,,$ (A.18)

which gives for the time-dependent $3$-point function

$\displaystyle\langle\pi(\boldsymbol{k}_{1},\eta)$
$\displaystyle\pi(\boldsymbol{k}_{2},\eta)\pi(\boldsymbol{k}_{3},\eta))\rangle^{\prime}=\frac{1}{(4\epsilon_{\star}M_{\rm
Pl}^{2})^{2}\epsilon_{\star}H^{2}\prod_{i}k_{i}^{3}}\times$
$\displaystyle\times\textrm{Im}\left[\int_{-\infty(1-i\epsilon)}^{\eta}\frac{\,\mathrm{d}\eta^{\prime}}{\eta^{\prime
2}}\ddot{H}(\eta^{\prime})\left(\boldsymbol{k}_{2}\cdot\boldsymbol{k}_{3}f^{*}(\eta)f(\eta^{\prime})+k_{2}k_{3}f^{*}(\eta)g_{1}(\eta^{\prime})\right)\right]+\textrm{2
perms.}\;,$ (A.19)

with $f(\eta)$ as above, and $g_{1}(\eta)=e^{ik_{\rm
t}\eta}(1-ik_{1}\eta)k_{2}k_{3}\eta^{2}$. Contrary to the polynomial form,
here all the integrals are manifestly convergent, and the late-time constancy
of the $3$-point function is immediate.

Let us also compute the bispectrum starting from (A.19). Taking $\eta\to 0$
and writing the sine function as a sum of two exponentials, this gives, with
$\beta=-\tilde{\delta}-\alpha\log\eta_{\star}$:

$\displaystyle B(k_{1},k_{2},k_{3})=\ $
$\displaystyle\frac{3H_{\star}b}{(4\epsilon_{\star}M_{\rm
Pl}^{2})^{2}\prod_{i}k_{i}^{3}}\times$
$\displaystyle\times\textrm{Im}\left[\sum_{n=0}^{3}\sum_{\sigma=\pm
1}\Gamma(i\sigma\alpha+n-1)\frac{a_{n}}{k_{\rm t}^{i\sigma\alpha+n-1}}\
e^{i\sigma\beta}e^{-\sigma\alpha\pi/2}\right]\;,$ (A.20)

with

$\displaystyle a_{0}$
$\displaystyle=-\frac{1}{2}\sum_{i}k_{i}^{2}\;,\qquad\hskip
66.01059pta_{1}=-\frac{1}{2}k_{\rm t}(\sum_{i}k_{i}^{2})\,,$ (A.21)
$\displaystyle a_{2}$ $\displaystyle=\frac{k_{\rm t}}{2}\bigg{(}-\sum_{i\neq
j}k_{i}^{2}k_{j}+k_{1}k_{2}k_{3}\bigg{)}\;,\quad
a_{3}=-k_{1}k_{2}k_{3}\bigg{(}\sum_{i}k_{i}^{2}/2+\sum_{i<j}k_{i}k_{j}\bigg{)}\,.$

Simplifying, one finds

$B(k_{1},k_{2},k_{3})=\frac{3H_{\star}b\,\alpha^{3}\,\cosh(\alpha\pi/2)}{(4\epsilon_{\star}M_{\rm
Pl}^{2})^{2}\prod_{i}k_{i}^{3}}\,\textrm{Im}\,z\;,$ (A.22)

with

$z=\left[\prod_{i}k_{i}-\frac{i}{\alpha}\left(\prod_{i}k_{i}-\sum_{i\neq
j}k_{i}^{2}k_{j}\right)-\frac{\sum_{i}k_{i}^{3}}{\alpha^{2}}\right]\Gamma(-1-i\alpha)e^{i(\alpha\log
k_{\rm t}-\beta+\pi/2)}\;.$ (A.23)

In the limit of large $\alpha$, with $\cosh(\pi\alpha/2)\to e^{\pi\alpha/2}/2$
and
$\Gamma(-1-i\alpha)\to\sqrt{2\pi}\,e^{-\pi\alpha/2}\alpha^{-3/2}e^{i(\alpha-\alpha\log\alpha+3\pi/4)}[1+13i/(12\alpha)]$
taking into account next-to-leading order (NLO) terms, this reduces to

$B(k_{1},k_{2},k_{3})=\frac{3H_{\star}b\,\alpha^{3/2}\sqrt{2\pi}}{2(4\epsilon_{\star}M_{\rm
Pl}^{2})^{2}\prod_{i}k_{i}^{2}}\left[\sin(\alpha\log k_{\rm
t}+\varphi)+\frac{1}{\alpha}\cos(\alpha\log k_{\rm
t}+\varphi)\left(\sum_{i\neq
j}\frac{k_{i}}{k_{j}}+\frac{1}{12}\right)\right]\;,$ (A.24)

with $\varphi=\alpha-\alpha\log\alpha+5\pi/4-\beta$, and where we have kept
the first subleading term, see [14] for similar result.

##  Appendix B Beyond the decoupling limit: mixing with gravity

In the main text, we have derived the full nonlinear action (2.7), at all
orders in $\pi$ and zeroth-order in the decoupling limit. This action scales
like ${\cal O}(\epsilon)$. Here, we derive first-order corrections to this
action, i.e. up to ${\cal O}(\epsilon^{2})$, also keeping the full nonlinear
structure in $\pi$. From this, we deduce an upper bound on $\bar{\zeta}$ for
the mixing with gravity to be negligible.

We perform the space-time dependent time diffeomorphism $t\rightarrow
t+\pi(t,\boldsymbol{x})$, starting from the unitary gauge action (2.1). We use
the ADM parametrization of the metric:

$\,\mathrm{d}s^{2}=-N^{2}\,\mathrm{d}t^{2}+\hat{g}_{ij}(\,\mathrm{d}x^{i}+N^{i}\,\mathrm{d}t)(\,\mathrm{d}x^{j}+N^{j}\,\mathrm{d}t)\,,$
(B.1)

we neglect tensor modes, and we choose the spatially flat gauge
$\hat{g}_{ij}=a^{2}\delta_{ij}$. One thus obtains the action

$\displaystyle S=\int\,\mathrm{d}^{4}x\;a^{3}M_{\rm Pl}^{2}$
$\displaystyle\left\\{\frac{1}{2N}(E_{ij}^{2}-E^{2})-N(3H^{2}(t+\pi)+\dot{H}(t+\pi))\right.$
$\displaystyle\qquad+\left.\dot{H}(t+\pi)\left[-N^{-1}(1+\dot{\pi})^{2}+2N^{-1}(1+\dot{\pi})N^{i}\partial_{i}\pi+N\left(\partial_{i}\pi\right)^{2}-N^{-1}(N^{i}\partial_{i}\pi)^{2}\right]\right\\}\,,$

where

$E_{ij}=\frac{1}{2}\dot{\hat{g}}_{ij}-N_{(i;j)}=a^{2}H\delta_{ij}-N_{(i;j)}\
,\quad N_{(i;j)}={1\over 2}(N_{i;j}+N_{j;i})\,,$ (B.2)

and ; stands for the covariant derivative with respect to the spatial metric
$\hat{g}_{ij}$, which in this gauge is simply an ordinary derivative.

Performing similar manipulations as in the main text, this reads

$\displaystyle S=\int\,\mathrm{d}^{4}x\,a(t)^{3}M_{\rm Pl}^{2}$
$\displaystyle\left\\{-\dot{H}(t+\pi)(\dot{\pi}^{2}-(\partial_{i}\pi)^{2})-3(H(t+\pi)-H(t))^{2}\right.$
$\displaystyle+\left.\delta
N\left(3N^{-1}H^{2}(t)-3H^{2}(t+\pi)-2N^{-1}HN_{i}^{\,;i}\right)+\frac{1}{2}N^{-1}(N_{(i;j)}^{2}-(N_{i}^{\,;i})^{2})\right.$
$\displaystyle+\left.\dot{H}(t+\pi)\left(\delta
N(N^{-1}(1+\dot{\pi})^{2}-1+(\partial_{i}\pi)^{2})+2(1+\dot{\pi})N^{-1}N^{i}\partial_{i}\pi-N^{-1}(N^{i}\partial_{i}\pi)^{2}\right)\right\\}\;,$

where contractions are made with the spatial metric
$\hat{g}_{ij}=a^{2}\delta_{ij}$, we wrote $N=1+\delta N$, and no approximation
has been made so far. Varying the action with respect to $\delta N$ and
$N_{i}$, we obtain, respectively, the lapse constraint:

$\displaystyle 3(H^{2}(t)-H^{2}(t+\pi))-3H^{2}(t+\pi)\delta N(2+\delta
N)-2H(t)\partial^{i}N_{i}-\frac{1}{2}(N_{(i;j)}^{2}-(\partial^{i}N_{i})^{2})$
(B.3)
$\displaystyle+\dot{H}(t+\pi)[2\dot{\pi}+\dot{\pi}^{2}+(N^{i}\partial_{i}\pi)^{2}-2(1+\dot{\pi})N^{i}\partial_{i}\pi+\left(\partial_{i}\pi\right)^{2}(1+\delta
N)^{2}-\delta N(2+\delta N)]=0\,,$

and shift constraint:

$\displaystyle
2N^{-1}\dot{H}(t+\pi)[1+\dot{\pi}-N^{j}\partial_{j}\pi]\partial_{i}\pi=\partial_{j}\biggl{[}\frac{1}{2}N^{-1}(N_{i}^{\,;j}+N^{j}_{\,;i})-\delta_{i}^{j}(2H(t)\delta
N/N+N_{k}^{\,\,;k})\bigg{]}\,.$ (B.4)

In order to derive the first correction to the decoupling limit action, we see
that it is enough to work out $\delta N$ and $N_{i}$ at first order in
$\epsilon$, which we denote with the superscript (1). From Eqs. (B.3)–(B.4),
and using the Helmoltz decomposition of the shift,
$N_{i}=\partial_{i}\psi+\tilde{N}_{i}$, with $\partial^{i}\tilde{N}_{i}=0$,
one obtains the compact expressions

$\delta
N^{(1)}=\frac{1}{H(t)}\partial^{-2}\partial_{i}X_{i}\,\quad\textrm{and}\quad\tilde{N}^{(1)}_{i}=4a^{2}\partial^{-2}\left(\partial_{i}\partial^{-2}\partial_{j}X_{j}-X_{i}\right)\,,\quad\textrm{with}\quad
X_{i}\equiv-\dot{H}(t+\pi)(1+\dot{\pi})\partial_{i}\pi\,,$ (B.5)

and

$2H(t)\frac{\partial^{2}\psi^{(1)}}{a^{2}}=3\left(H^{2}(t)-H^{2}(t+\pi)\right)-6H^{2}(t+\pi)\delta
N^{(1)}+\dot{H}(t+\pi)[2\dot{\pi}+\dot{\pi}^{2}+(\partial_{i}\pi)^{2}]\,,$
(B.6)

where indices are simply contracted with $\delta_{ij}$ here. One then obtains
the final expression of the action including first-order $\epsilon$
corrections:

$\displaystyle S=\int\,\mathrm{d}^{4}x\,a(t)^{3}M_{\rm Pl}^{2}$
$\displaystyle\left\\{-\dot{H}(t+\pi)(\dot{\pi}^{2}-(\partial_{i}\pi)^{2})-3(H(t+\pi)-H(t))^{2}-3(H^{2}(t+\pi)-H^{2}(t))\delta
N^{(1)}\right.$ (B.7) $\displaystyle+\left.3H^{2}(t)(\delta
N^{(1)})^{2}+\dot{H}(t+\pi)(2\dot{\pi}+\dot{\pi}^{2}+(\partial_{i}\pi)^{2})\delta
N^{(1)}+\frac{1}{2}\tilde{N}^{(1)}_{(i;j)}\tilde{N}^{(1)(i;j)}-2X_{i}\tilde{N}^{(1)i}\right\\}\,,$

where we note that the explicit expression of $\psi^{(1)}$ is actually not
needed for this result because of structural cancellations in the computation.

Let us now use the action (B.7) to identify the regime of validity of the
decoupling limit analysis. For this, in the same spirit as in Sec. 4, we
consider a boundary profile $\bar{\zeta}$ characterized by a typical momentum
scale $k$. This way, the effects of spatial derivatives simply read
$\partial_{i}\sim k_{i}$ and $\partial^{-2}\sim k^{-2}$. We also use the
behaviour of the free mode function (3.14) to deduce the estimates
$\pi\sim\bar{\zeta}/H(1+k\eta)$ and $\dot{\pi}\sim\bar{\zeta}(k\eta)^{2}$. We
thus find $X_{i}\sim\epsilon
H\bar{\zeta}(1+\bar{\zeta}(k\eta)^{2})(1+k\eta)k_{i}$,
$\tilde{N}^{(1)}_{i}\sim a^{2}/k^{2}X_{i}$ and in particular

$\delta
N^{(1)}\sim\epsilon\bar{\zeta}\,[1+\bar{\zeta}(k\eta)^{2}](1+k\eta)\,.$ (B.8)

The validity of the decoupling limit necessitates $\delta N^{(1)}\ll 1$. The
right-hand side of (B.8) grows with $k\eta$, but keep in mind that the
corresponding interactions are shut off deep inside the horizon, due to the
$i\epsilon$ prescription projecting onto the interacting vacuum. However, we
should demand that $\delta N^{(1)}\ll 1$ at the resonance, such that
$|k\eta|\sim\alpha$. This imposes the bound

$|\bar{\zeta}|\ll\frac{1}{\sqrt{\epsilon\alpha^{3}}}\,.$ (B.9)

This turns out to be the most stringent bound on $\bar{\zeta}$ when requiring
the ${\cal O}(\epsilon^{2})$ interactions to be negligible compared to the
leading ${\cal O}(\epsilon)$ part of the action. Checking this is
straightforward, except for the second term in (B.7). The corresponding
polynomial interactions $(H(t+\pi)-H(t))^{2}\sim\epsilon^{2}H^{4}\pi^{2}$ do
not decay outside the Hubble radius, contrary to the leading two-derivative
part of the action. This is simply a manifestation of the fact that $\pi$
acquires a mass beyond the decoupling limit. Requiring that these interactions
are subdominant at the late-time saddle point (4.12) gives the bound
$\epsilon|\bar{\zeta}|\ll 1$, which is indeed less stringent than (B.9).

Importantly, let us highlight that the regime $\alpha^{2}|\bar{\zeta}|\sim 1$,
in which full non-perturbative results are needed and we found qualitative
deviations from perturbation theory, always lies within the regime of validity
(B.9) of the decoupling limit analysis. Notice as well that in order to get
the range of validity for $\bar{\zeta}$ of the decoupling limit it is crucial
to use our non-perturbative form of the action; it is not enough to stick to
the perturbative analysis discussed in [13].

##  Appendix C Correlators from the WFU at one loop

In this section we are going to investigate how Witten diagrams are related to
the correlators at loop level. In the WFU approach, equal-time correlators are
obtained by performing a path integral over field configurations at that given
time. At one-loop level, on top of Witten-diagram loops (needed to obtain the
WFU), there are in principle additional loops originating from ‘averaging’
tree-level Witten diagrams over boundary (late-times) field configurations. In
turn, this step can lead to additional divergences that are not manifest in
the WFU. When estimating the size of loops, we therefore need to take into
account both contributions.

On the other hand, correlators can also be evaluated using different methods
that do not involve the WFU, as for instance the ‘in-in formalism’ [23]. In
this case, one does not make a distinction between WFU and boundary loops.
Notice also that the propagators running inside loops used in the two methods
are different: the in-in uses Wightman functions, while the WFU uses bulk-to-
bulk propagators. As it is perhaps expected, on the WFU side we will obtain
that the two sets of loop diagrams combine to yield back the in-in result: the
effect of boundary loops is to change the boundary conditions of the
propagator. 232323For a similar discussion on the relation between the in-in
and WFU approaches to loop diagrams, see [24].

For concreteness, we consider $\lambda\phi^{4}$ in dS, although our conclusion
appears general. The action is taken to be

$S=\int\,\mathrm{d}^{4}x\,a^{4}(\eta)\left[-\frac{1}{2}(\partial_{\mu}\phi)^{2}-\frac{\lambda}{4!}\phi^{4}\right]\;.$
(C.1)

We furthermore focus on the equal time two-point function
$\braket{\phi(\eta,\boldsymbol{x})\phi(\eta,\boldsymbol{y})}$ at linear order
in $\lambda$. At this order, the WFU at time $\eta_{\rm f}$,
$\Psi[\bar{\phi};\eta_{\rm f}]$, contains a tree-level correction to the four-
point coefficient $\psi_{4}$ and a one-loop correction to the two-point
coefficient $\psi_{2}$. Both are needed to obtain the one-loop correlator at
$\mathcal{O}(\lambda)$.

In the following calculations we will need the bulk-to-bulk and the bulk-to-
boundary propagators for massless fields in dS. First, let us define the wave-
modes in $k$-space that solve the free equation of motion for $\phi$

$\varphi_{\pm}(k,\eta)\equiv(1\mp ik\eta)e^{\pm ik\eta}\;.$ (C.2)

Then, the bulk-to-bulk propagator is obtained as

$\displaystyle G(\eta,\eta^{\prime};\boldsymbol{k})$
$\displaystyle=-\frac{iH^{2}}{2k^{3}}\left[\theta(\eta-\eta^{\prime})\varphi_{+}(k,\eta^{\prime})\varphi_{-}(k,\eta)+\theta(\eta^{\prime}-\eta)\varphi_{+}(k,\eta)\varphi_{-}(k,\eta^{\prime})\right.$
$\displaystyle\;\qquad\qquad\left.-\frac{\bar{\varphi}_{-}}{\bar{\varphi}_{+}}\varphi_{+}(k,\eta)\varphi_{+}(k,\eta^{\prime})\right]\;,$
(C.3)

where $\bar{\varphi}_{\pm}(k)\equiv\varphi_{\pm}(k,\eta_{\rm f})$ and $\theta$
is the Heaviside theta function. By construction, the function
$G(\eta,\eta^{\prime};\boldsymbol{k})$ vanishes at early times (when $\eta$,
$\eta^{\prime}\to-\infty(1-i\epsilon)$) and at late times (when $\eta$,
$\eta^{\prime}\to\eta_{\rm f}$). Requiring that $G$ satisfies the free
equations of motion with a $\delta(\eta-\eta^{\prime})$ source then fixes the
normalization in Eq. (C) (i.e. we are working with canonically-normalized
fields). Note that this normalization changes if we work with $\zeta$ (in this
case the factor $H^{2}$ gets replaced by $P_{\zeta}$). On the other hand, the
bulk-to-boundary propagator is

$K_{\boldsymbol{k}}(\eta)\equiv\frac{\varphi_{+}(k,\eta)}{\bar{\varphi}_{+}(k)}\;.$
(C.4)

This function is obtained as the free solution of the equations of motion
satisfying $K_{\boldsymbol{k}}(\eta_{\rm f})=1$ and vanishing for
$\eta\to-\infty(1-i\epsilon)$.

After introducing these quantities, we are ready to evaluate the Witten
diagrams contributing to the WFU. The tree-level contribution is obtained by
evaluating the interaction part of the action times $i$ on the free modes for
$\phi$ (given by Eq. (3.14), but for $\phi$). Therefore, we obtain

$\displaystyle\log\Psi_{4}[\bar{\phi};\eta_{\rm f}]$
$\displaystyle=-\frac{i\lambda}{4!}\int_{-\infty}^{\eta_{\rm
f}}\frac{\,\mathrm{d}\eta}{(H\eta)^{4}}\int\prod_{i=1}^{4}\left[\frac{\,\mathrm{d}^{3}\boldsymbol{k}_{i}}{(2\pi)^{3}}\,\phi(\eta,\boldsymbol{k}_{i})\right](2\pi)^{3}\,\delta^{(3)}(\boldsymbol{k}_{\rm
t})$
$\displaystyle=\frac{1}{4!}\int\prod_{i=1}^{4}\frac{\,\mathrm{d}^{3}\boldsymbol{k}_{i}}{(2\pi)^{3}}\,(2\pi)^{3}\,\delta^{(3)}(\boldsymbol{k}_{\rm
t})\,\psi_{4}(\boldsymbol{k}_{1},\ldots,\boldsymbol{k}_{4};\eta_{\rm
f})\,\bar{\phi}(\boldsymbol{k}_{1})\ldots\bar{\phi}(\boldsymbol{k}_{i})\;,$
(C.5)

where $\log\Psi_{n}[\bar{\phi};\eta_{\rm f}]$ is the correction to the
exponent of the WFU with $n$ fields and $\boldsymbol{k}_{\rm
t}\equiv\boldsymbol{k}_{1}+\boldsymbol{k}_{2}+\boldsymbol{k}_{3}+\boldsymbol{k}_{4}$.
Thus, written in terms of bulk-to-boundary propagators
$K_{\boldsymbol{k}_{i}}(\eta)$, the coefficient $\psi_{4}$ reads

$\psi_{4}(\boldsymbol{k}_{1},\boldsymbol{k}_{2},\boldsymbol{k}_{3},\boldsymbol{k}_{4};\eta_{\rm
f})=-i\lambda\int_{-\infty}^{\eta_{\rm
f}}\frac{\,\mathrm{d}\eta}{(H\eta)^{4}}\prod_{i=1}^{4}K_{\boldsymbol{k}_{i}}(\eta)\;.$
(C.6)

The one-loop contribution to $\psi_{2}$ can instead be obtained using the
standard Witten rules (see [35] for a derivation in $\rm AdS$)

$\displaystyle\log\Psi_{2}[\bar{\phi};\eta_{\rm f}]$
$\displaystyle=\frac{\lambda}{4}\int_{-\infty}^{\eta_{\rm
f}}\frac{\,\mathrm{d}\eta}{(H\eta)^{4}}\int\frac{\,\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}K_{\boldsymbol{k}}^{2}(\eta)\,\bar{\phi}(\boldsymbol{k})\bar{\phi}(-\boldsymbol{k})\int\frac{\,\mathrm{d}^{3}\boldsymbol{p}}{(2\pi)^{3}}G(\eta,\eta;\boldsymbol{p})$
$\displaystyle\equiv\frac{1}{2!}\int\frac{\,\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}\,\delta\psi_{2}(\boldsymbol{k};\eta_{\rm
f})\,\bar{\phi}(\boldsymbol{k})\bar{\phi}(-\boldsymbol{k})\;,$ (C.7)

where we used the total momentum-conserving delta function to remove one
integral. Therefore, the correction to $\psi_{2}$ is given by

$\delta\psi_{2}(\boldsymbol{k};\eta_{\rm
f})=\frac{\lambda}{2}\int_{-\infty}^{\eta_{\rm
f}}\frac{\,\mathrm{d}\eta}{(H\eta)^{4}}K_{\boldsymbol{k}}^{2}(\eta)\int\frac{\,\mathrm{d}^{3}\boldsymbol{p}}{(2\pi)^{3}}G(\eta,\eta;\boldsymbol{p})\;.$
(C.8)

Explicitly, the WFU at order $\lambda$ can be expanded as follows

$\Psi[\bar{\phi};\eta_{\rm
f}]=\mathcal{N}\,\Bigl{(}1+\log\Psi_{2}[\bar{\phi};\eta_{\rm
f}]+\log\Psi_{4}[\bar{\phi};\eta_{\rm
f}]\Bigr{)}\exp\left[\frac{1}{2!}\int\frac{\,\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}\psi_{2}(\boldsymbol{k};\eta_{\rm
f})\,\bar{\phi}(\boldsymbol{k})\bar{\phi}(-\boldsymbol{k})\right]\,,$ (C.9)

where
$\mathcal{N}^{-1}\equiv\int\mathcal{D}\bar{\phi}\,|\Psi[\bar{\phi};\eta_{\rm
f}]|^{2}$ is a normalization constant. The two-point correlator in momentum
space is then obtained using the standard Born rule of quantum mechanics:

$\braket{\phi(\eta_{\rm f},\boldsymbol{p})\phi(\eta_{\rm
f},\boldsymbol{p}^{\prime})}=\int\mathcal{D}\bar{\phi}\;\bar{\phi}(\boldsymbol{p})\bar{\phi}(\boldsymbol{p}^{\prime})\,\bigl{|}\Psi[\bar{\phi};\eta_{\rm
f}]\bigr{|}^{2}\;.$ (C.10)

[width=0.8]figures/bndry_loops

Figure 10:  The two one-loop contributions to $\braket{\phi(\eta_{\rm
f},\boldsymbol{p})\phi(\eta_{\rm f},\boldsymbol{p}^{\prime})}$. The grey
surface represents the boundary at $\eta=\eta_{\rm f}$. On the left, $(a)$, we
have the one-loop diagram on the bulk, $\log\Psi_{2}$, that contributes
directly to the WFU. On the right, $(b)$, we have the boundary diagram where
we average over two boundary fields in $\log\Psi_{4}$.

In performing this path integral at order $\lambda$ we encounter various types
of contributions. There are, for instance, bubble diagrams, that are however
cancelled once the normalization $\mathcal{N}$ is taken into account. In
principle, there could be disconnected diagrams, but at this order they do not
appear. Finally, there are two types of connected diagrams. First, from
$\delta\psi_{2}$, a diagram where the external $\bar{\phi}(\boldsymbol{p})$,
$\bar{\phi}(\boldsymbol{p}^{\prime})$ connect to the internal
$\bar{\phi}(\boldsymbol{k})\bar{\phi}(-\boldsymbol{k})$. This can be
interpreted as a loop contribution from the bulk (given the origin of
$\delta\psi_{2}$). Second, from $\psi_{4}$, we have a loop diagram where two
internal $\bar{\phi}(\boldsymbol{k}_{i})$’s are contracted among each other.
Therefore, this corresponds to a boundary loop. These two one-loop diagrams
are represented respectively in Fig. 10 (a) and (b).

One can straightforwardly check that

${\rm Re}\,\psi_{2}(\boldsymbol{k};\eta_{\rm
f})=-\frac{k^{3}}{H^{2}\bar{\varphi}_{+}(k)\bar{\varphi}_{-}(k)}\;.$ (C.11)

At this point, we are ready to compute the order $\lambda$ correction in Eq.
(C.10) (the integral can be performed, for instance, by first introducing a
generating functional and by taking its functional derivatives). We obtain

$\displaystyle\braket{\phi(\eta_{\rm f},\boldsymbol{p})\phi(\eta_{\rm
f},\boldsymbol{p}^{\prime})}^{\prime}_{(1)}=\frac{1}{2\left({\rm
Re\,}\psi_{2}(\boldsymbol{p};\eta_{\rm f})\right)^{2}}\left[{\rm
Re}\,\delta\psi_{2}(\boldsymbol{p};\eta_{\rm
f})-\frac{1}{4}\int\frac{\,\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{{\rm
Re\,}\psi_{4}(\boldsymbol{p},\boldsymbol{p},\boldsymbol{k},\boldsymbol{k};\eta_{\rm
f})}{{\rm Re\,\psi_{2}(\boldsymbol{k};\eta_{\rm f})}}\right]$
$\displaystyle=\frac{\lambda}{4\left({\rm
Re\,}\psi_{2}(\boldsymbol{p};\eta_{\rm f})\right)^{2}}\,{\rm
Re}\,\int_{-\infty}^{\eta_{\rm
f}}\frac{\,\mathrm{d}\eta}{(H\eta)^{4}}K_{\boldsymbol{p}}(\eta)^{2}\int\frac{\,\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}\left[G(\eta,\eta;\boldsymbol{k})+\frac{iK_{\boldsymbol{k}}^{2}(\eta)}{2\,{\rm
Re\,}\psi_{2}(\boldsymbol{k};\eta_{\rm f})}\right]\;.$ (C.12)

In the correlator, the subscript $(1)$ stands for the first correction in
$\lambda$, while the prime means we remove
$(2\pi)^{3}\delta^{(3)}(\boldsymbol{p}+\boldsymbol{p}^{\prime})$. In going to
the second line we used Eqs. (C.6) and (C.8). Notice that, although the two
contributions of the first line come with different symmetry factors, when
expressed in terms of bulk objects this difference cancels. Indeed, as it is
clear from Fig. 10, both diagrams have the same combinatorics from a bulk
perspective. In the second line on the RHS of (C.12), the two terms in the
square bracket are given by

$\displaystyle G(\eta,\eta;\boldsymbol{k})$
$\displaystyle=-\frac{iH^{2}}{2k^{3}}\left[\varphi_{+}(k,\eta)\varphi_{-}(k,\eta)-\frac{\bar{\varphi}_{-}(k)}{\bar{\varphi}_{+}(k)}\varphi_{+}(k,\eta)\varphi_{+}(k,\eta)\right]\;,$
(C.13) $\displaystyle\frac{iK_{\boldsymbol{k}}^{2}(\eta)}{2\,{\rm
Re\,}\psi_{2}(\boldsymbol{k};\eta_{\rm f})}$
$\displaystyle=-\frac{iH^{2}}{2k^{3}}\bar{\varphi}_{+}(k)\bar{\varphi}_{-}(k)\left(\frac{\varphi_{+}(k,\eta)}{\bar{\varphi}_{+}(k,\eta)}\right)^{2}\;.$
(C.14)

By combining these two pieces, the term that imposes the boundary conditions
at $\eta_{\rm f}$ in $G(\eta,\eta,\boldsymbol{k})$ (the last in Eq. (C.13)) is
cancelled by the boundary loop (C.14). Therefore, we arrive at the overall
one-loop contribution to the correlator

$\braket{\phi(\eta_{\rm f},\boldsymbol{p})\phi(\eta_{\rm
f},\boldsymbol{p}^{\prime})}^{\prime}_{(1)}=\frac{\lambda H^{2}}{8\left({\rm
Re\,}\psi_{2}(\boldsymbol{p};\eta_{\rm f})\right)^{2}}\,{\rm
Im}\,\int_{-\infty}^{\eta_{\rm
f}}\frac{\,\mathrm{d}\eta}{(H\eta)^{4}}K_{\boldsymbol{p}}(\eta)^{2}\int\frac{\,\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{1}{k^{3}}\varphi_{+}(k,\eta)\varphi_{-}(k,\eta)\;,$
(C.15)

where indeed we recognize the propagator of the in-in relations inside the
loop. Notice also that the WFU loop, $\delta\psi_{2}$ in Eq. (C.8), is
convergent in the IR, while the loop for the correlator (C.15) is divergent
[16]. Moreover, given that the loops are typically divergent, in order to have
finite correlators after renormalization we arrive at the conclusion that the
WFU is necessarily a divergent object on itself. The divergences from the
additional path integral in Eq. (C.10) will then make the final result UV
finite. One can also check that Eq. (C.15) agrees with the in-in computation
[36]. The cancellation between bulk and boundary terms is completely
independent on the spacetime, the theory and the evaluation time $\eta_{\rm
f}$, as it is based on the properties of the in-in formalism. Therefore, the
conclusion also applies to the model of resonant non-Gaussianities.

##  Appendix D Estimate of one-loop Witten diagrams at
$\mathcal{O}(\tilde{b})$

In this appendix we want to show that the one-loop Witten diagrams with $n$
external legs at first order in $\tilde{b}$ (see Fig. 1) scale as
$\sqrt{\alpha}\,\tilde{b}\,\bar{\zeta}^{2}(\alpha^{2}\bar{\zeta})^{n-2}(\alpha^{2}P_{\zeta})/P_{\zeta}$.
In comparison with the tree-level diagrams with the same number of external
legs, these loop corrections are suppressed by $\alpha^{2}P_{\zeta}\lesssim
1$. Note that this condition is weaker than the one of perturbative unitarity,
$\alpha^{4}P_{\zeta}\lesssim 1$, see the discussion around Eq. (3.25).

For illustrative purposes, let us focus on the one-loop Witten diagram with
two external legs (second graph in the second line of Fig. 1). The ingredients
for computing such a diagram are two bulk-to-boundary propagators, one bulk-
to-bulk propagator and the four-point vertex which is given by the quartic
interaction:

$\displaystyle\mathcal{L}_{4}=\frac{\tilde{b}\alpha^{2}\zeta^{2}}{4\eta^{2}P_{\zeta}}\bigg{[}\zeta^{\prime
2}-(\partial_{i}\zeta)^{2}\bigg{]}\cos\left(\alpha\log(\eta/\eta_{\star})-\tilde{\delta}\right)\;.$
(D.1)

The coupling above can be straightforwardly obtained from expanding the action
(3.13) for small $\zeta$ up to quartic order. 242424In general, the non-linear
self-coupling of $\zeta$’s at first order in $\tilde{b}$ reads
$\displaystyle\mathcal{L}_{n+2}=-\frac{\tilde{b}\alpha^{n}\zeta^{n}}{2\eta^{2}P_{\zeta}n!}\bigg{[}\zeta^{\prime
2}-(\partial_{i}\zeta)^{2}\bigg{]}\left\\{\begin{matrix}(-1)^{n/2}\cos\left(\alpha\log(\eta/\eta_{\star})-\tilde{\delta}\right)\;,&n+2\in{\rm
even}\\\
(-1)^{(n+1)/2}\sin\left(\alpha\log(\eta/\eta_{\star})-\tilde{\delta}\right)\;,&n+2\in{\rm
odd}\end{matrix}\right.\;.$ (D.2)

Derivatives in the action can act on internal or external legs. Let us focus
first on the case where the two derivatives act on the two bulk-to-boundary
propagators: this is the leading contribution (highest power of $\alpha$) as
we will verify later. Using the vertex (D.1), the one-loop coefficient of the
WFU with two external legs is given by 252525We define the wavefunction
coefficients as in Eq. (C.9).

$\displaystyle\psi^{\rm
1-loop}_{2}(\boldsymbol{k}_{1},\boldsymbol{k}_{2})\supset\frac{\tilde{b}\alpha^{2}}{4P_{\zeta}}\int^{\eta_{\rm
f}}_{-\infty}\frac{\,\mathrm{d}\eta}{\eta^{2}}K^{\prime}_{\boldsymbol{k}_{1}}(\eta)K^{\prime}_{\boldsymbol{k}_{2}}(\eta)\cos(\alpha\log(\eta/\eta_{\star})-\tilde{\delta})\int\frac{\,\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}G(\eta,\eta;\boldsymbol{k})\;,$
(D.3)

where the subscript denotes the number of external legs, $\boldsymbol{k}_{1}$
and $\boldsymbol{k}_{2}$ are the external momenta, and $\boldsymbol{k}$ is the
internal momentum. Note that in the expression (D.3) we considered two time
derivatives acting on $K_{\boldsymbol{k}}(\eta)$. One gets the same scaling in
$\alpha$ also for spatial derivatives: this is straightforward to verify using
the saddle-point approximation ($\alpha\gg 1$) to the $\eta$-integral in the
early-time limit. 262626Subleading corrections in the limit $|k\eta|\gg 1$
contain less power of $\eta$: on the saddle point of the $\eta$ integral,
their contributions have fewer power of $\alpha$ compared to the leading
contributions. Using the formula (C.13) and the fact that
$K^{\prime}_{\boldsymbol{k}}(\eta)=\eta k^{2}e^{ik\eta}/\bar{\varphi}_{+}(k)$,
Eq. (D.3) becomes

$\displaystyle\psi^{\rm 1-loop}_{2}(\boldsymbol{k}_{1},\boldsymbol{k}_{2})$
$\displaystyle\supset-\frac{i\tilde{b}\alpha^{2}H^{2}}{8}\frac{k_{1}^{2}k_{2}^{2}}{\bar{\varphi}_{+}(k_{1})\bar{\varphi}_{+}(k_{2})}\int^{\eta_{\rm
f}}_{-\infty}\,\mathrm{d}\eta\,e^{ik_{\rm
t}\eta}\cos\big{(}\alpha\log(\eta/\eta_{\star})-\tilde{\delta}\big{)}\times$
$\displaystyle\hskip
12.80365pt\times\int\frac{\,\mathrm{d}^{3}\boldsymbol{k}}{(2\pi)^{3}}\frac{1}{k^{3}}\left[\varphi_{+}(k,\eta)\varphi_{-}(k,\eta)-\frac{\bar{\varphi}_{-}(k)}{\bar{\varphi}_{+}(k)}\varphi_{+}(k,\eta)\varphi_{+}(k,\eta)\right]\;,$
(D.4)

where $k_{\rm t}\equiv k_{1}+k_{2}$. Setting $\eta_{\rm f}=0$ (272727Since one
gets a non-zero result for $\eta_{\rm f}=0$, it is safe to neglect subleading
terms. As we will discuss, it is useful to keep $\eta_{\rm f}\neq 0$ to study
the Minkowski limit of the result.) we have

$\displaystyle\psi^{\rm 1-loop}_{2}(\boldsymbol{k}_{1},\boldsymbol{k}_{2})$
$\displaystyle\supset-\frac{i\tilde{b}\alpha^{2}H^{2}}{8}k_{1}^{2}k_{2}^{2}\int_{-\infty}^{0}\,\mathrm{d}\eta~{}~{}e^{ik_{\rm
t}\eta}\cos\big{(}\alpha\log(\eta/\eta_{\star})-\tilde{\delta}\big{)}\times$
$\displaystyle\hskip 12.80365pt\times\frac{4\pi}{(2\pi)^{3}}\int_{0}^{\Lambda
a(\eta)}\frac{\,\mathrm{d}k}{k}\left[(1+k^{2}\eta^{2})-(1-ik\eta)^{2}e^{2ik\eta}\right]\;,$
(D.5)

where $\Lambda$ is a fixed physical cutoff. Notice that this physical cutoff
$\Lambda$ appears together with the scale factor at time $\eta$, cutting off
the integral in comoving momentum $k$ (see a more detailed discussion below
Eq. (26) of [36]). Then, we change the variable using $k_{\rm p}\equiv
k/a(\eta)=-\eta Hk$, where $k_{\rm p}$ denotes a physical momentum. Thus, the
integral above becomes

$\displaystyle\psi^{\rm 1-loop}_{2}(\boldsymbol{k}_{1},\boldsymbol{k}_{2})$
$\displaystyle\supset-\frac{i\tilde{b}\alpha^{2}H^{2}}{16\pi^{2}}k_{1}^{2}k_{2}^{2}\int_{-\infty}^{0}\,\mathrm{d}\eta~{}e^{ik_{\rm
t}\eta}\cos\big{(}\alpha\log(\eta/\eta_{\star})-\tilde{\delta}\big{)}\times$
$\displaystyle\hskip
12.80365pt\times\int_{0}^{\Lambda}\frac{\,\mathrm{d}k_{\rm p}}{k_{\rm
p}}\left[\bigg{(}1+\frac{k_{\rm p}^{2}}{H^{2}}\bigg{)}-\bigg{(}1+\frac{ik_{\rm
p}}{H}\bigg{)}^{2}e^{-2ik_{\rm p}/H}\right]\;.$ (D.6)

We now see that the integral over $k_{\rm p}$ in the second line does not
depend on $\eta$ and in fact it potentially leads to quadratic and logarithmic
divergences. More precisely, we obtain

$\displaystyle\log\Psi_{2}^{\rm 1-loop}$
$\displaystyle\supset\sqrt{\alpha}~{}\tilde{b}\,\bar{\zeta}^{2}\alpha^{2}\frac{k_{1}^{2}k_{2}^{2}}{k_{\rm
t}}\,e^{-i[\alpha\log(X_{\star})+\tilde{\delta}]}\bigg{[}\log\bigg{(}\frac{\Lambda}{H}\bigg{)}+\frac{\Lambda^{2}}{2H^{2}}+\,{\rm
finite\,terms}\bigg{]}\;,$ (D.7)

where we have used $X_{\star}=-k_{\rm t}\eta_{\star}$, and we have evaluated
the $\eta$ integral on the perturbative saddle point, i.e.
$\eta_{s}=-\alpha/k_{\rm t}$. Note that the expression above is an estimate,
assuming that $\bar{\zeta}(\boldsymbol{k})\simeq\bar{\zeta}$ and it is peaked
at some momentum. From Eq. (D.7), we see that the quadratic divergence, as
expected, can be removed by adding a local counter-term to the WFU, resulting
in, for example, renormalization of the mass term. On the other hand, the
logarithmic divergence of the form $\log(\Lambda/H)$, as argued in [36],
cannot be removed by a local counter-term since for the modes inside the
horizon such a divergence becomes non-local, $\log((a(t)\Lambda)/k)$. With
this reasoning, this logarithmic divergence is physical and we can read off
the dependence on $\alpha$ for $\psi^{\rm 1-loop}_{2}$.

Alternatively, one could draw the same conclusion using the dimensional
regularization (dim. reg.) to compute $\psi^{\rm 1-loop}_{2}$. In dim. reg.
one usually performs the $k$ integral in $d=3+\varepsilon$ dimensions. We use
$d$ to denote the number of spatial dimensions. Additionally, in
$(d+1)$-dimensions the free mode function becomes the Hankel function
$(-H\eta)^{d/2}H^{(1)}_{d/2}(-k\eta)$. Since the integral involving the Hankel
function is very complicated, one can then use a trick, proposed in [37], to
obtain a simple form of the mode function in $d=3+\varepsilon$ dimensions. One
considers a massive scalar field in $(d+1)$ dimensional de Sitter space, and
analytically continues both the number of spatial dimensions and the mass of
the field in such a way that the index of the Hankel function remains $3/2$.
By doing so, the mode function in $d=3+\varepsilon$ dimensions takes a simple
form, $(-H\eta)^{\varepsilon/2}(1-ik\eta)e^{ik\eta}$. We see that this mode
function is different from the one in $d=3$ dimensions by the overall
normalization factor $(-H\eta)^{\varepsilon/2}$. Using such a simple form of
the mode function it is then straightforward to perform the integral in $k$
analytically in the limit $\eta_{\rm f}=0$. At this point, one can see that
there are several terms appearing after the integration: a finite term, a term
going as $1/\varepsilon$ and a term that contains $\log(H/\mu)$ where $\mu$ is
a renormalization scale. Note that the finite term and the term with
$1/\varepsilon$ can be altogether removed in some renormalization scheme.
After that, similar to the calculation in (D.7), one can apply the saddle-
point approximation to perform the $\eta$ integral. Finally, we obtain the
same $\alpha$ dependence of $\psi^{\rm 1-loop}_{2}$ together with the
logarithmic divergence of the form $\log(H/\mu)$ which, as explained earlier,
cannot be removed by adding local counter-terms. 282828 We can convince
ourselves that $\log(H/\mu)$ originates from a non-analytic term in $k$ by
performing the calculation in dim. reg. in the early-times limit $|\eta_{\rm
f}|\gg 1$, where all the modes are inside the horizon and the mode functions
are approximated by the Minkowski ones (by time-translational invariance in
this limit, we can set the final time $t_{\rm f}=0$). The time and momentum
integrals in $\psi_{2}^{1-{\rm loop}}$ in this limit take the following form
$\sum_{\sigma}\int_{-\infty}^{0}\,\mathrm{d}t\,e^{-it(k_{\rm
tp}-\sigma\omega)}\int\frac{\,\mathrm{d}^{d}\boldsymbol{k}_{\rm
p}}{(2\pi)^{d}}\frac{\mu^{\delta}}{k_{\rm p}}\left[1-e^{-2ik_{\rm
p}t}\right]\propto\frac{k_{\rm
tp}^{2}+\omega^{2}}{4\delta}+\frac{1}{8}\left[(k_{\rm
tp}^{2}+\omega^{2})\log\frac{4\mu^{2}}{k_{\rm tp}^{2}-\omega^{2}}+k_{\rm
tp}\omega\log\frac{k_{\rm tp}-\omega}{k_{\rm tp}+\omega}\right]\;,$ (D.8)
where in the first step we expanded the cosine in exponentials ($\sigma=\pm
1$) and in the second step we performed the summation over $\sigma$. The
subscript $\rm p$ stands for physical momentum. We first performed the time
integral and then the one over $k_{\rm p}$. Notice that the boundary term
$e^{-2ik_{\rm p}t}$ mixes the $t$ and $k_{\rm p}$ dependences and does not
vanish in dim. reg. . Note the non-analytic dependence of the finite terms in
$\delta$.

The result (D.7) indicates that this one-loop correction is suppressed
compared to the tree-level one by $\alpha^{2}P_{\zeta}\ll 1$. This confirms
the estimates of Sec. 3.2. In addition, following the same method we used, it
is straightforward to generalise the scaling of $\log\Psi_{2}^{\rm 1-loop}$ to
$\log\Psi_{n}^{\rm 1-loop}$ where $n$ is the number of external legs:

$\displaystyle\log\Psi^{\rm
1-loop}_{n}\sim\tilde{b}\sqrt{\alpha}\,\frac{\bar{\zeta}^{2}}{P_{\zeta}}\,(\alpha^{2}\bar{\zeta})^{n-2}(\alpha^{2}P_{\zeta})\;.$
(D.9)

The suppression with respect to tree-level is always $\alpha^{2}P_{\zeta}\ll
1$.

Let us discuss what happens when derivatives do not act on the external legs:
contributions such as
$K_{\boldsymbol{k}_{1}}(\eta)K^{\prime}_{\boldsymbol{k}_{1}}(\eta)G^{\prime}(\eta,\eta;\boldsymbol{k})$
and
$K_{\boldsymbol{k}_{1}}(\eta)K_{\boldsymbol{k}_{1}}(\eta)G^{\prime\prime}(\eta,\eta;\boldsymbol{k})$
to the one-loop coefficient $\psi_{2}^{\rm 1-loop}$. Following the same
procedure as above, one can straightforwardly show that these contributions
are suppressed by $\alpha$, compared to Eq. (D.7). More explicitly, we have

$\displaystyle\log\Psi_{2}^{\rm
1-loop}(KK^{\prime}G^{\prime})\sim\frac{\tilde{b}\sqrt{\alpha}}{\alpha}\,\frac{\bar{\zeta}^{2}}{P_{\zeta}}\,(\alpha^{2}P_{\zeta})\;,\quad\log\Psi_{2}^{\rm
1-loop}(KKG^{\prime\prime})\sim\frac{\tilde{b}\sqrt{\alpha}}{\alpha^{2}}\,\frac{\bar{\zeta}^{2}}{P_{\zeta}}\,(\alpha^{2}P_{\zeta})\;.$
(D.10)

The reason for this $\alpha$ suppression is the following. In this case, there
is at least one time derivative acting on the bulk-to-bulk propagator,
implying that it will generate a factor of internal momentum $k$. Then,
changing the variables to the physical momentum ($k_{\rm p}=-\eta Hk$) leads
to an additional factor of $1/\eta$, so that the integrand in $\eta$ contains
fewer powers of $\eta$ compared with the
$K^{\prime}_{\boldsymbol{k}_{1}}(\eta)K^{\prime}_{\boldsymbol{k}_{1}}(\eta)G(\eta,\eta;\boldsymbol{k})$
contribution. Therefore, evaluating the $\eta$ integral on the saddle point
($k_{\rm t}\eta_{s}\sim-\alpha$) we find that in comparison with Eq. (D.7)
these contributions are suppressed in $\alpha$ as shown in Eq. (D.10).

Applying the same technique to the higher-loop corrections, one can deduce
that the $\ell$-loop wavefunction coefficient at $\mathcal{O}(\tilde{b})$ with
$n$ external legs scales as

$\displaystyle\log\Psi^{\rm\ell-Loop}_{n}\sim\tilde{b}\sqrt{\alpha}\,\frac{\
\bar{\zeta}^{2}}{P_{\zeta}}\,(\alpha^{2}\bar{\zeta})^{n-2}\,(\alpha^{2}P_{\zeta})^{\ell}\;.$
(D.11)

##  Appendix E Numerical methods

In this appendix we explain the numerical methods used to compute the on-shell
action in Sec. 5.2. The first step towards numerically integrating the action
is to obtain an accurate value for the free solution $\zeta(\tau,r)$ with the
late-time boundary condition. We then provide two methods to obtain such a
solution. The first method makes use of the analytical expression for the
Fourier transform $\zeta(\tau,k)$, i.e. the Euclidean version of Eq. (3.14).
Then, the real-space solution is obtained by performing an inverse-Fourier
transform at each time step, which is implemented numerically via a fast
Fourier transform (FFT). Since our system possesses spherical symmetry, the
solution then only depends on $k$ and the inverse-Fourier transform simply
becomes one dimensional. After that, we plug the real-space solution into the
action and evaluate the integral numerically, although there are some
technical points to be careful about, as we will explain below. The second
method, on the other hand, relies on solving the linear PDE for
$\zeta(\tau,r)$ with prescribed boundary conditions. Then, we can
straightforwardly compute the action on such a numerical solution. This method
in fact is similar to what was implemented in [1] (in this reference the
equation for $\zeta(\tau,r)$ was however non-linear). We will refer to these
two methods as FFT and PDE methods, respectively.

Before entering into the details of the numerical integration, let us comment
on the Gaussian profile at late times, see Eq. (5.18). In fact, in this
specific case the solution $\zeta(\tau,r)$ can be obtained analytically in
terms of exponential-integral functions, using the inverse-Fourier transform
of the (Euclidean-rotated) $\zeta(\tau,k)$ in (3.14), giving

$\zeta(\tau,r)/\bar{\zeta}=(1-4\tau^{2})\,{\rm
Re}\,W(z)+\frac{2\tau}{r}(\tau^{2}-r^{2})\,{\rm
Im}\,W(z)-\frac{2\tau}{\sqrt{\pi}}\;,$ (E.1)

where $W(z)\equiv e^{-z^{2}}(1-i\,{\rm erfi}(z))$, ${\rm erfi}(z)\equiv-i{\rm
erf}(iz)$ is the imaginary error function and $z\equiv r+i\tau$. It should be
noted that dealing with the above expression is not always straightforward.
Indeed, an accurate evaluation of Eq. (E.1) requires the implementation of
arbitrary-precision numerics (otherwise large numerical errors appear when
evaluating $W(z)$ for complex argument in the regions $|\tau|\gg 1$, $r\gg
1$), resulting in very long evaluation times. Actually, we employ this
analytical expression (E.1) only for negative $\bar{\zeta}$ since high
accuracy/precision is needed, as discussed in Sec. 5.2. However, this
procedure of finding an analytical expression for $\zeta(\tau,r)$ cannot be
generally applied to an arbitrary late-time configuration.

Here we explain in detail how to perform the numerical integration over $\tau$
and $r$. First, note that $\tau\in(-\infty,0)$ and $r\in(0,+\infty)$. In order
to obtain a better accuracy/precision, we divide our integral (5.1) into three
pieces,

$\displaystyle\Delta S_{\rm E,1}=\Delta S_{1,{\rm early}}+\Delta
S_{1,\rm{grid}}+\Delta S_{1,\rm{late}}\;,$ (E.2)

where $\Delta S_{1,{\rm early}}$ is the integration over
$\tau\in(-\infty,\tau_{\rm min})$, $\Delta S_{1,\rm{grid}}$ the one over
$\tau\in(\tau_{\rm min},\tau_{\rm max})$, and $\Delta S_{1,\rm{late}}$ the one
over $\tau\in(\tau_{\rm max},0)$. The reason for this separation, as we will
see below, is essentially the fact that the integral over $\tau$ in both
$\Delta S_{1,{\rm early}}$ and $\Delta S_{1,{\rm late}}$ can be performed
analytically, which indeed improves the matching with the saddle-point
approximation. Therefore, we are left with only the numerical integration in
$r$ for $\Delta S_{1,{\rm early}}$ and $\Delta S_{1,{\rm late}}$. Of course,
we still have to do the numerical integrations over $\tau$ and $r$ in $\Delta
S_{1,\rm{grid}}$, 292929In fact, when performing the numerical integration it
is more useful to change coordinates $\\{\tau,r\\}$ to
$\\{\tilde{t},\tilde{r}\\}$, defined as  $\tilde{t}\equiv-\log(-\tau
H)\;,\quad\tilde{r}\equiv\log(Hr)\;.$ (E.3) Notice that both $\tilde{t}$ and
$\tilde{r}$ run from $-\infty$ to $+\infty$. We see that the oscillation of
the on-shell action is periodic in $\tilde{t}$ with frequency $\alpha$,
whereas its frequency increases as a function of $\tau$ at late times.
Therefore, this suggests that performing the numerical integral over the
variable $\tilde{t}$ is more robust (the sampling can be done on a uniform
grid), compared to the integral over $\tau$. but this is less subtle since
both $\tau_{\rm min}$ and $\tau_{\rm min}$ are finite. Below, we provide a
detailed analysis of both $\Delta S_{1,{\rm early}}$ and $\Delta S_{1,{\rm
late}}$.

Let us analyse the late-time contribution $\Delta S_{1,{\rm late}}$. In this
regime, expanding the solution $\zeta(\tau,r)$ for small $|\tau|$, see Eq.
(3.16), gives rise to the late-time limit action up to the leading
contributions,

$\displaystyle\Delta S_{1,{\rm
late}}\simeq\frac{2\pi}{P_{\zeta}}\int_{0}^{+\infty}\,\mathrm{d}r\,r^{2}\int_{\tau_{\rm
max}}^{0}\,\mathrm{d}\tau$
$\displaystyle\bigg{[}\left(\nabla^{2}\bar{\zeta}(r)+\partial_{r}\bar{\zeta}(r)\,\partial_{r}\nabla^{2}\bar{\zeta}(r)\right)\cos
q(\tau,r)$
$\displaystyle\qquad\qquad\qquad-\frac{\alpha}{2}\nabla^{2}\bar{\zeta}(r)\,(\partial_{r}\bar{\zeta}(r))^{2}\sin
q(\tau,r)\bigg{]}\;,$ (E.4)

where we have used the fact that the late-time profile is spherically
symmetric and we have defined
$q(\tau,r)\equiv\alpha\left(\log(\tau/\eta_{\star})+\bar{\zeta}(r)-i\pi/2\right)-\tilde{\delta}$.
It is important to note that the limit $\tau\to 0^{-}$ is regular by
construction. Performing the integral over $\tau$ in (E.4) analytically, we
therefore obtain

$\displaystyle\Delta S_{1,{\rm late}}\simeq$
$\displaystyle-\frac{2\pi\tau_{\rm
min}}{(1+\alpha^{2})P_{\zeta}}\int_{0}^{+\infty}\,\mathrm{d}r\,r^{2}\bigg{[}\left(\nabla^{2}\bar{\zeta}(r)+\partial_{r}\bar{\zeta}(r)\,\partial_{r}\nabla^{2}\bar{\zeta}(r)\right)(\cos
q_{\rm max}(r)+\alpha\sin q_{\rm max}(r))$ $\displaystyle\hskip
142.26378pt-\frac{\alpha}{2}\nabla^{2}\bar{\zeta}(r)\,(\partial_{r}\bar{\zeta}(r))^{2}(\sin
q_{\rm max}(r)-\alpha\cos q_{\rm max}(r))\bigg{]}\;,$ (E.5)

where $q_{\rm max/min}(r)\equiv q(\tau_{\rm max/min},r)$. Notice that in the
formula shown above we have kept the terms up to $\mathcal{O}(\tau_{\rm
max})$. Actually, in our numerical implementation we include up to order
$\mathcal{O}(\tau_{\rm max}^{12})$ and when choosing $\tau_{\rm max}$ we check
that additional corrections are negligible. The integration over $r$ in (E.5)
is then performed numerically.

Finally, let us consider the early-time contribution $\Delta S_{1,{\rm
early}}$. Following the same procedure as for the late-time limit action, we
expand the solution $\zeta(\tau,r)$ for large $|\tau|$ to obtain

$\displaystyle\Delta
S_{1,\rm{early}}\simeq-\frac{2\pi}{P_{\zeta}}\int_{0}^{+\infty}\,\mathrm{d}r\,r^{2}(\partial_{r}\bar{\zeta}(r))^{2}\int_{-\infty}^{\tau_{\rm
min}}\frac{\,\mathrm{d}\tau}{\tau^{2}}\cos(q(\tau,r))\;.$ (E.6)

Notice that in the expression above the fact that $\zeta(\tau,r)$ quickly
decays for large $|\tau|$ implies that the dominant contribution in the action
at early times comes from the counter term in Eq. (3.21). Evaluating the
integral over $\tau$ analytically, we thus obtain

$\displaystyle\Delta S_{1,\rm{early}}=\frac{2\pi}{\tau_{\rm
min}(1+\alpha^{2})P_{\zeta}}\int_{0}^{+\infty}\,\mathrm{d}r\,r^{2}(\partial_{r}\bar{\zeta}(r))^{2}\left(\cos
q_{\rm min}(r)-\alpha\sin q_{\rm min}(r)\right)\;.$ (E.7)

As before, we are left with the radial integration which can be done
numerically. Notice that this contribution decays very slowly at early times,
as $\sim 1/\tau_{\rm min}$. This suggests that including $\Delta
S_{1,\rm{early}}$ allows us to choose moderately large values for $\tau_{\rm
min}$, without considering an extremely large grid.

## References

  * [1] M. Celoria, P. Creminelli, G. Tambalo, and V. Yingcharoenrat, “Beyond perturbation theory in inflation,” JCAP 06 (2021) 051, 2103.09244.
  * [2] P. Creminelli, S. Dubovsky, A. Nicolis, L. Senatore, and M. Zaldarriaga, “The Phase Transition to Slow-roll Eternal Inflation,” JHEP 09 (2008) 036, 0802.1067.
  * [3] S. Dubovsky, L. Senatore, and G. Villadoro, “The Volume of the Universe after Inflation and de Sitter Entropy,” JHEP 04 (2009) 118, 0812.2246.
  * [4] S. Dubovsky, L. Senatore, and G. Villadoro, “Universality of the Volume Bound in Slow-Roll Eternal Inflation,” JHEP 05 (2012) 035, 1111.1725.
  * [5] 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, 1804.07315.
  * [6] 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, 1806.05202.
  * [7] G. Panagopoulos and E. Silverstein, “Primordial Black Holes from non-Gaussian tails,” 1906.02827.
  * [8] G. Panagopoulos and E. Silverstein, “Multipoint correlators in multifield cosmology,” 2003.05883.
  * [9] G. A. Palma and S. Sypsas, “Non-Gaussian statistics of de Sitter spectators: A perturbative derivation of stochastic dynamics,” 2309.16474.
  * [10] X. Chen, R. Easther, and E. A. Lim, “Generation and Characterization of Large Non-Gaussianities in Single Field Inflation,” JCAP 04 (2008) 010, 0801.3295.
  * [11] R. Flauger, L. McAllister, E. Pajer, A. Westphal, and G. Xu, “Oscillations in the CMB from Axion Monodromy Inflation,” JCAP 06 (2010) 009, 0907.2916.
  * [12] R. Flauger and E. Pajer, “Resonant Non-Gaussianity,” JCAP 01 (2011) 017, 1002.0833.
  * [13] S. R. Behbahani, A. Dymarsky, M. Mirbabayi, and L. Senatore, “(Small) Resonant non-Gaussianities: Signatures of a Discrete Shift Symmetry in the Effective Field Theory of Inflation,” JCAP 12 (2012) 036, 1111.3373.
  * [14] C. Duaso Pueyo and E. Pajer, “A Cosmological Bootstrap for Resonant Non-Gaussianity,” 2311.01395.
  * [15] L. McAllister, E. Silverstein, and A. Westphal, “Gravity Waves and Linear Inflation from Axion Monodromy,” Phys. Rev. D 82 (2010) 046003, 0808.0706.
  * [16] V. Gorbenko and L. Senatore, “$\lambda\phi^{4}$ in dS,” 1911.00022.
  * [17] L. Pinol, S. Renaux-Petel, and Y. Tada, “A manifestly covariant theory of multifield stochastic inflation in phase space: solving the discretisation ambiguity in stochastic inflation,” JCAP 04 (2021) 048, 2008.07497.
  * [18] V. Vennin, Stochastic inflation and primordial black holes.  PhD thesis, U. Paris-Saclay, 6, 2020.  2009.08715.
  * [19] S. Céspedes, A.-C. Davis, and D.-G. Wang, “On the IR Divergences in de Sitter Space: loops, resummation and the semi-classical wavefunction,” 2311.17990.
  * [20] C. Cheung, P. Creminelli, A. L. Fitzpatrick, J. Kaplan, and L. Senatore, “The Effective Field Theory of Inflation,” JHEP 0803 (2008) 014, 0709.0293.
  * [21] E. Pajer, G. L. Pimentel, and J. V. S. Van Wijck, “The Conformal Limit of Inflation in the Era of CMB Polarimetry,” JCAP 06 (2017) 009, 1609.06993.
  * [22] L. Leblond and E. Pajer, “Resonant Trispectrum and a Dozen More Primordial N-point functions,” JCAP 01 (2011) 035, 1010.4565.
  * [23] J. M. Maldacena, “Non-Gaussian features of primordial fluctuations in single field inflationary models,” JHEP 0305 (2003) 013, astro-ph/0210603.
  * [24] M. H. G. Lee, C. McCulloch, and E. Pajer, “Leading loops in cosmological correlators,” JHEP 11 (2023) 038, 2305.11228.
  * [25] A. Hook and R. Rattazzi, “Softening the UV without new particles,” Phys. Rev. D 108 (2023), no. 11 115019, 2306.12489.
  * [26] G. Tambalo, M. Zumalacárregui, L. Dai, and M. H.-Y. Cheung, “Lensing of gravitational waves: Efficient wave-optics methods and validation with symmetric lenses,” Phys. Rev. D 108 (2023), no. 4 043527, 2210.05658.
  * [27] M. Serone, G. Spada, and G. Villadoro, “The Power of Perturbation Theory,” JHEP 05 (2017) 056, 1702.04148.
  * [28] A. Falkowski and R. Rattazzi, “Which EFT,” JHEP 10 (2019) 255, 1902.05936.
  * [29] S. Chang and M. A. Luty, “The Higgs Trilinear Coupling and the Scale of New Physics,” JHEP 03 (2020) 140, 1902.05556.
  * [30] N. S. M. de Santi et. al., “Field-level simulation-based inference with galaxy catalogs: the impact of systematic effects,” 2310.15234.
  * [31] F. Arroja and T. Tanaka, “A note on the role of the boundary terms for the non-Gaussianity in general k-inflation,” JCAP 1105 (2011) 005, 1103.1102.
  * [32] C. Burrage, R. H. Ribeiro, and D. Seery, “Large slow-roll corrections to the bispectrum of noncanonical inflation,” JCAP 1107 (2011) 032, 1103.4126.
  * [33] G. Rigopoulos, “Gauge invariance and non-Gaussianity in Inflation,” Phys. Rev. D84 (2011) 021301, 1104.0292.
  * [34] S. Garcia-Saenz, L. Pinol, and S. Renaux-Petel, “Revisiting non-Gaussianity in multifield inflation with curved field space,” JHEP 01 (2020) 073, 1907.10403.
  * [35] D. Anninos, T. Anous, D. Z. Freedman, and G. Konstantinidis, “Late-time Structure of the Bunch-Davies De Sitter Wavefunction,” JCAP 11 (2015) 048, 1406.5490.
  * [36] L. Senatore and M. Zaldarriaga, “On Loops in Inflation,” JHEP 12 (2010) 008, 0912.2734.
  * [37] S. Melville and E. Pajer, “Cosmological Cutting Rules,” JHEP 05 (2021) 249, 2103.09832.