IFT UAM-CSIC 24-60

# One-loop power spectrum in ultra slow-roll inflation and
implications for primordial black hole dark matter

Guillermo Ballesteros1,2<EMAIL_ADDRESS>Jesús Gambín Egea2
<EMAIL_ADDRESS>1Departamento de Física Teórica, Universidad Autónoma de
Madrid (UAM), Campus de Cantoblanco, 28049 Madrid, Spain  2Instituto de Física
Teórica UAM-CSIC, Campus de Cantoblanco, 28049 Madrid, Spain

###### Abstract

We apply the in-in formalism to address the question of whether the size of
the one-loop spectrum of curvature fluctuations in ultra-slow-roll inflation
models designed for producing a large population of primordial black holes
implies a breakdown of perturbation theory. We consider a simplified piece-
wise description of inflation, in which the ultra-slow-roll phase is preceded
and followed by slow-roll phases linked by transitional periods. We work in
the $\delta\phi$-gauge, including all relevant cubic and quartic interactions
and the necessary counterterms to renormalize the ultraviolet divergences,
regularized by a cutoff. The ratio of the one-loop to the tree-level
contributions to the spectrum of curvature perturbations is controlled by the
duration of the ultra-slow-roll phase and of the transitions. Our results
indicate that perturbation theory does not necessarily break in well-known
models proposed to account for all the dark matter in the form of primordial
black holes.

###### Contents

  1. I Introduction
  2. II Model of USR inflation
  3. III In-In formalism, regularization and counterterms
  4. IV Structure of $\mathcal{P}_{\zeta}$
  5. V Discussion
  6. A In-in formalism in more detail
  7. B Action for fluctuations in the $\delta\phi$-gauge
    1. B.1 Expansion of the potential
    2. B.2 Expansions of $N$ and $N_{i}$
    3. B.3 Quadratic, cubic and quartic actions for fluctuations
  8. C Classical dynamics of super-Hubble modes
  9. D Change of gauge
    1. D.1 Scalar field transformation
    2. D.2 Metric transformation
    3. D.3 Two-point correlation function of $\zeta$
  10. E Tadpoles
  11. F Asymptotic limits
    1. F.1 Limit $k\to 0$
    2. F.2 Limit $k\to\infty$
  12. G Dimensional analysis and a quick estimate of $\mathcal{P}_{\zeta}$
    1. G.1 The limit $\delta N\to 0$
  13. H Bispectrum and consistency condition

##  I Introduction

Primordial black holes (PBHs) in the range that goes from $10^{-16}$ to
$10^{-12}$ Solar masses might account for all the dark matter (DM) of the
Universe Carr  _et al._ (2010); Niikura  _et al._ (2019); Katz  _et al._
(2018); Montero-Camacho  _et al._ (2019); Ballesteros  _et al._ (2020a); Carr
_et al._ (2021); Green and Kavanagh (2021). Although the existence of PBHs has
not been demonstrated, that possibility alone makes them worth studying.
Indeed, they have become serious contenders for the solution to the dark
matter problem, on the same footing as particle candidates such as e.g. WIMPS
and warmer candidates, the QCD axion and axion-like particles. Moreover,
whether they are discovered or ruled out in the future, PBH physics can help
us learn about the very early Universe.

A popular mechanism to form PBHs consists in the collapse of Hubble-sized
density fluctuations in the radiation epoch Carr and Hawking (1974). These
collapsing regions, with densities above a certain threshold Nakama  _et al._
(2014); Musco (2019), could have originated from specific dynamics during a
preceding inflationary period. The most studied model of that kind posits a
phase of so-called ultra slow-roll (USR) Tsamis and Woodard (2004); Kinney
(2005) during inflation. In terms of a single scalar field –the inflaton,
$\phi$– driving inflation as it descends a potential $V(\phi)$, this
corresponds to $\phi$ undergoing a significant deceleration, such that
$\ddot{\phi}\simeq-3H\dot{\phi}$, where dots indicate derivatives with respect
to cosmic time, $t$, and $H$ is the Hubble function measuring the
(accelerated) expansion of the Universe. This dynamics implies a rapid growth
of curvature fluctuations, which in turn leads to large density perturbations
during the subsequent radiation epoch.

A simple estimate assuming Gaussian primordial curvature fluctuations
indicates that their spectrum $\mathcal{P_{\zeta}}$ (which we will define
later on) must be $\sim 10^{-2}$ at comoving wavenumbers of the order of
$10^{12}$ Mpc-1 – $10^{14}$ Mpc-1 if all the dark matter is made of PBHs
formed from that mechanism in the aforementioned mass range. This spectral
value is much larger than the inferred one from the cosmic microwave
background (CMB), which is $\mathcal{P_{\zeta}}\simeq 2\cdot 10^{-9}$ Fixsen
_et al._ (1996). Since $\mathcal{P_{\zeta}}$ scales as the square of the
curvature fluctuation, $\zeta$, values of $\mathcal{P_{\zeta}}$ of order
$10^{-2}$ almost beg the question of the validity of perturbation theory. A
paper by J. Kristiano and J. Yokoyama from November 2022 Kristiano and
Yokoyama (2022) studied this question considering the one-loop correction to
$\mathcal{P_{\zeta}}$ in USR inflation, using the in-in formalism Schwinger
(1961); Keldysh (1964); Weinberg (2005). They concluded that such large values
of $\zeta$ imply the breakdown of perturbation theory at CMB scales, which,
according to their work, would severely threaten USR as a possibility for PBH
formation. Other papers on the same issue have been appearing since then, with
contradicting conclusions that have kept the question open, see Riotto (2023);
Firouzjahi (2023a); Franciolini  _et al._ (2023); Tada  _et al._ (2023);
Iacconi  _et al._ (2023) for a non-exhaustive but representative list of
references.

We weigh in this debate improving the treatment of the problem in several
respects. As in most earlier works on the topic, we use the in-in formalism to
compute $\mathcal{P}_{\zeta}$ at one-loop and compare the result to its tree-
level counterpart. However, we include the complete set of relevant
interactions (cubic, quartic and boundary terms) for fluctuations, several of
which where ignored in previous analyses.111Although there are no boundary
terms in the set of interactions we consider, these interactions are related
to boundary terms in the $\zeta$-gauge, on which the problem has been studied
in previous works. We manage to do this by working on the $\delta\phi$-gauge,
choosing spacetime coordinates for which $\zeta=0$ and the physical content of
the primordial fluctuations is contained in the perturbation of the inflaton,
$\delta\phi$. In addition, we implement a regularization procedure of
ultraviolet divergences at the level of the action for fluctuations. We use a
cutoff to regularize the ultraviolet divergences, similarly to what has been
done in previous literature on the topic, see e.g. Kristiano and Yokoyama
(2022); Franciolini  _et al._ (2023). However, we also absorb the divergences
by introducing adequate counterterms, whereas in previous works the cutoff was
given a numerical value motivated by physical intuition (which made the result
depend on its value). Although we do not implement a full renormalization
procedure, we are able to analyze the validity of perturbation theory,
providing an answer to the issue raised in Kristiano and Yokoyama (2022).

We find that whether perturbation theory breaks down depends on the duration
of the transition between slow-roll (SR) inflation and USR inflation. Our
results indicate that for $\mathcal{P_{\zeta}}\sim 10^{-2}$ and well-motivated
inflationary models considered in the literature Ballesteros and Taoso (2018);
Ballesteros  _et al._ (2020b), cosmological perturbation theory is valid, in
the sense that the one-loop spectrum is significantly smaller than the tree-
level one.

##  II Model of USR inflation

We consider a two-parameter, piece-wise, description of the inflationary
dynamics leading to a large tree-level $\mathcal{P_{\zeta}}$. We assume a
phase of SR followed by USR and then, again, SR. The transitions between these
phases are characterized by a single parameter, which controls the duration of
both transitions. The other parameter of the model is the duration of the USR
phase. Imposing a value for $\mathcal{P_{\zeta}}$ establishes a one-to-one
correspondence between both parameters. The dynamics of the SR phases is
controlled by the small quantity $0<\epsilon=-\dot{H}/H^{2}\ll 1$, whose
actual value is largely irrelevant in what follows. The quantity
$\eta=\dot{\epsilon}/(H\epsilon)$ is, by definition, exactly equal to $-6$ in
the USR phase and, will be assumed to vanish during both SR phases. The
transitions between phases are described in terms of a third quantity, $\nu$,
as we explain next.

Assuming a single inflaton $\phi$, with a canonical kinetic term and a
potential $V(\phi)$, working in conformal time, $\tau$, and neglecting terms
suppressed by powers of $\epsilon$, the action for fluctuations in the
$\delta\phi$-gauge is

$\displaystyle
S={\int}\differential\tau\,\differential^{3}\mathbf{x}\,\left[\frac{a^{2}}{2}\left(\left(\partial_{\tau}\delta\phi\right)^{2}-\left(\partial_{i}\delta\phi\right)^{2}\right)-a^{4}\sum_{n\geq
2}\frac{V_{n}\,\delta\phi\,^{n}}{n!}\right]$ (1)

and $V_{n}=\differential^{n}V/\differential\phi^{n}$. The interactions that
arise from the metric fluctuations are suppressed (see Appendix B), and only
the interactions coming from the potential survive at lowest order in
$\epsilon$:

$\displaystyle a^{2}V_{2}=-(aH)^{2}(\nu^{2}-9/4)\,,\quad
a^{2}V_{3}=-\dfrac{aH(\nu^{2})^{\prime}}{\sqrt{2\epsilon}M_{P}}\,,\quad
a^{2}V_{4}=-\dfrac{1}{2\epsilon
M_{P}^{2}}\left((\nu^{2})^{\prime\prime}-aH(\nu^{2})^{\prime}\left(1+\dfrac{\eta}{2}\right)\right)\,,$
(2)

where primes denote derivatives with respect to conformal time,

$\displaystyle\nu^{2}\equiv\frac{9}{4}+\frac{1}{2}\left(3\,\eta+\frac{\eta^{2}}{2}+\frac{\eta^{\prime}}{aH}\right)\in\mathbb{R}$
(3)

and $M_{P}^{2}=1/(8\pi\,G)$ is the reduced Planck mass squared. In both SR
($\epsilon,\,|\eta|\ll 1$) and USR, $\nu=3/2$, see also Wands (1999); Leach
and Liddle (2001); Leach  _et al._ (2001). We impose that $\nu^{2}$ is piece-
wise constant. The function $\nu^{2}$ in the transition from SR to USR and,
also in the subsequent transition from USR to SR, is then set by their (equal)
duration. This in turn determines $\eta$ and $\epsilon$ completely, which are
found integrating their respective definitions. Figure 1 shows in an example
$\nu^{2}$, $\epsilon$ and $\eta$ as functions of the number of e-folds of
inflation ($N=\int Hdt$) elapsed from the beginning of the first transition.
In terms of this variable, the duration of the USR phase is denoted $\Delta N$
and that of the transitions is $\delta N$. For $\delta N=0$, we recover the
model used in Kristiano and Yokoyama (2022). For the (well-motivated) known
potentials that lead to transient USR compatible with $\mathcal{P_{\zeta}}\sim
10^{-2}$ Ballesteros and Taoso (2018); Ballesteros  _et al._ (2020b), the
function $\nu^{2}$ is indeed approximately constant during the transitions,
which satisfy $\delta N\lesssim 1<\Delta N$.

Since $\nu$ is discontinuous at the beginning and at the end of each
transition, the self-interactions of $\delta\phi$ (proportional to $V_{3,4}$)
are Dirac deltas centered on those instants. In the limit $\delta N\to 0$,
$\nu^{2}$ satisfies $\absolutevalue{\nu^{2}}\rightarrow 3/\delta N$ during the
transitions; i.e. the interactions diverge in the limit of instantaneous
transitions between SR and USR. It is therefore important to consider smooth
transitions. This effect is not as transparent in the $\zeta$-gauge because
the dependence on $\nu^{2}$ itself does not arise so naturally.

Although we use the $\delta\phi$-gauge, we are interested in
$\mathcal{P}_{\zeta}$, defined, at any order in perturbation theory, through
the two-point correlation function:

$\left\langle\zeta(x)\zeta(y)\right\rangle=\int\dfrac{\differential^{3}\mathbf{k}}{4\pi
k^{3}}e^{i\mathbf{k}(\mathbf{x}-\mathbf{y})}\mathcal{P}_{\zeta}(\tau,k)\,.$
(4)

For modes satisfying $k\equiv\absolutevalue{\mathbf{k}}\ll aH$ in the last SR
phase ($\eta=0$) we have
$\left\langle\delta\phi(x)\delta\phi(y)\right\rangle=2\epsilon
M_{P}^{2}\left\langle\zeta(x)\zeta(y)\right\rangle$ (see Appendix D).222In a
general model where in the last phase of inflation $\eta\neq 0$, this
expression has corrections proportional to the value of $\eta$ in that phase.

Figure 1:  Evolution of $\nu^{2}$ (black dotted), $\eta$ (black continuous)
and $\epsilon$ (orange continuous, right axis) as functions of $N$ for $\delta
N=0.25$, $\Delta N=2$ and $\epsilon(N<0)=10^{-3}$. The transitions (of
duration $\delta N$) between SR ($\eta=0$) and USR ($\eta=-6$) are shown as
shaded bands.

##  III In-In formalism, regularization and counterterms

In the in-in formalism Weinberg (2005); Chen (2010), at second order in the
interaction Hamiltonian, $H_{I}$, the vacuum expectation value of an operator
$Q(\tau)$ can be obtained as (see Appendix A)

$\displaystyle\expectationvalue{Q(\tau)}$
$\displaystyle=\bra{0}Q_{I}(\tau)\ket{0}+2\Im{\int_{-\infty}^{\tau}\differential\tau^{\prime}\bra{0}Q_{I}(\tau)H_{I}(\tau^{\prime}_{-})\ket{0}}$
$\displaystyle\quad+2\Re{\int^{\tau}_{-\infty}\differential\tau^{\prime}\int^{\tau}_{\tau^{\prime}}\differential\tau^{\prime\prime}\bra{0}{\left(H_{I}(\tau^{\prime\prime}_{+})Q_{I}(\tau)-Q_{I}(\tau)H_{I}(\tau^{\prime\prime}_{-})\right)}H_{I}(\tau^{\prime}_{-})\ket{0}}\,.$
(5)

On the right hand side of this expression, the fields are in the interaction
picture (I), i.e. they are described by the free Hamiltonian and satisfy
canonical commutation rules. $\ket{0}$ is the vacuum in the interaction
picture. We use $\tau_{\pm}\equiv\tau(1\pm i\omega)$ (with $\omega>0$) to
guarantee the correct projection onto the interaction vacuum in
$\tau\to-\infty$.

In a model with cubic and quartic interactions, the topology of the diagrams
generated from equation (5) for the two-point function of $\delta\phi$ are
shown in Figure 2.333Diagrams that are suppressed by powers of $M_{P}$ are not
shown, e.g. the diagrams induced at second order by the quartic interaction.
It can be checked that we do not have to include the bubble diagrams to the
actual computation we perform. The calculation at one-loop of the two-point
function gives an ultraviolet (UV) divergence that needs to be renormalized,
see Senatore and Zaldarriaga (2010); Weinberg (2008). For the regularization
of the divergences, we use a UV cutoff $\Lambda_{\rm UV}$. Whereas the time
integrals can be done explicitly, for the spatial part we work in Fourier
space. The cutoff goes into the momentum loops as usual, but its effect on the
time integrals needs to be worked out. In practice, we remove from the
integral a domain where the time intervals are smaller than those allowed by
$\Lambda_{\rm UV}$ Senatore and Zaldarriaga (2010), as follows:

$\int^{\tau}_{-\infty}\differential\tau^{\prime}\int^{\tau}_{\tau^{\prime}}\differential\tau^{\prime\prime}\to\int^{\tau+{1/(a(\tau)\Lambda_{\rm
UV})}}_{-\infty}\differential\tau^{\prime}\int^{\tau}_{\tau^{\prime}+1/(a(\tau^{\prime})\Lambda_{\rm
UV})}\differential\tau^{\prime\prime}\,.$ (6)

Once the divergences have been regularized, we introduce counterterms to
absorb the dependence of the loops on the regulator $\Lambda_{\rm UV}$
Weinberg (2011). To obtain the counterterms we can think in terms of the
action for the inflaton, $\phi$. General covariance implies that the wave
function renormalization, $\phi\to Z_{\phi}\,\phi_{R}$ (where the subscript R
stands for renormalized) requires the quantity $Z_{\phi}\equiv
1+\delta_{\phi}$ to be constant.444This renormalization also changes the
background action. However, by introducing a renormalization of the vacuum
expectation value (VEV) of $\phi$, $\phi_{0}(t)$, one can impose that the
background dynamics remains unchanged, equivalent to what happens with the
Higgs VEV (see e.g. Hollik (1990)). Similar arguments may in principle be used
for the counterterms coming from the potential. However, in the model under
consideration, $V(\phi)$ is given by the background evolution described in the
previous section. Therefore, just as $V_{n}$ are functions of time, by
renormalizing $V_{n}\to V_{R,n}+\delta_{V_{n}}$, the counterterms
$\delta_{V_{n}}$ are functions of time. The counterterms that will affect the
renormalization of the two-point correlation will be $\delta_{\phi}$ and
$\delta_{V_{2}}\equiv\delta_{V}$. Including these counterterms in the action
(1) we extract the complete interaction Hamiltonian in the interaction
picture:555For convenience, we will omit the subscript R, always keeping in
mind that we are working with renormalized quantities.

$\displaystyle
H_{I}(\tau)=\int\differential^{3}\mathbf{x}\,\bigg{[}a^{4}\left(\dfrac{V_{3}\delta\phi^{3}}{3!}+\dfrac{V_{4}\delta\phi^{4}}{4!}\right)+(a^{2}\delta_{\phi}\,\delta\phi\,\delta\phi^{\prime})^{\prime}+2a^{2}\,\delta\phi\,\left(a^{2}\tilde{\delta}_{V}-\delta_{\phi}\nabla^{2}\right)\delta\phi\bigg{]}\,,$
(7)

where $\tilde{\delta}_{V}=\delta_{\phi}V_{2}+\delta_{V}/4$. This Hamiltonian
includes all possible interaction terms and counterterms that affect
$\mathcal{P}_{\zeta}$ at one-loop at leading order in $\epsilon$.

We quantize $\delta\phi$ –in the interaction picture– with creation and
annihilation operators satisfying standard commutation relations,

$\delta\phi(x)=\int\dfrac{\differential^{3}\mathbf{k}}{(2\pi)^{3/2}}e^{i\mathbf{k}\mathbf{x}}\left(\delta\phi_{k}(\tau)a_{\mathbf{k}}+\delta\phi_{k}^{*}(\tau)a^{\dagger}_{-\mathbf{k}}\right)\,,$
(8)

where ∗ indicates complex conjugation. The modes $\delta\phi_{k}$ obey
$\delta\phi^{\prime\prime}_{k}+2aH\delta\phi_{k}^{\prime}+(k^{2}+a^{2}V_{2})\delta\phi_{k}=0$
with Bunch-Davies initial conditions
$\delta\phi_{k}(\tau\to-\infty)=e^{-ik\tau}/\sqrt{2ka^{2}}$, which guarantee
canonical commutation rules for $\delta\phi$ Bunch and Davies (1978). Their
evolution in time is

$\delta\phi_{k}(\tau)=(-k\tau)^{3/2}\left(\alpha_{k}J_{\nu}(-k\tau)+\beta_{k}Y_{\nu}(-k\tau)\right)\,,$
(9)

where $J_{\nu}$ and $Y_{\nu}$ are the first and second kind Bessel functions,
respectively. We stress that during the first transition $\nu^{2}$ can be
negative and therefore $\nu$ can be imaginary. The coefficients $\alpha_{k}$
and $\beta_{k}$ are obtained in each phase imposing that $\delta\phi_{k}$ and
$\delta\phi_{k}^{\prime}$ are continuous across the boundaries.

Figure 2:  Diagrams arising from eq. (5) when calculating
$\expectationvalue{\delta\phi^{2}}$. We distinguish the tree-level (first
row), the one-loop (first diagrams of rows second and third) and the
counterterms (second diagram of the second row). In addition, we have a term
proportional to the tadpole, bubble diagrams and a term composed of two
tadpoles. None of the latter diagrams (in shaded boxes) contribute to
$\mathcal{P}_{\zeta}$.

##  IV Structure of $\mathcal{P}_{\zeta}$

Due to momentum conservation

$\expectationvalue{\delta\phi(\tau,\mathbf{x})\delta\phi(\tau,\mathbf{y})}\equiv\int\dfrac{\differential^{3}\mathbf{k}}{(2\pi)^{3}}e^{i\mathbf{k}(\mathbf{x}-\mathbf{y})}D(\tau,k)\,,$
(10)

at all loop orders, for some function $D(\tau,k)=D(\tau,k)_{\rm
tl}+D(\tau,k)_{\rm 1l}+D(\tau,k)_{\rm ct}\,,$ where we have split the
contributions at tree-level (tl), one-loop (1l) and those coming from the
counterterms (ct) (see Figure 2), so that
$\mathcal{P}_{\zeta}(\tau,k)=k^{3}D(\tau,k)/(4\pi^{2}M_{P}^{2}\,\epsilon)$.
The tree-level contribution is just

$\mathcal{P}_{\zeta}^{\rm
tl}(\tau,k)=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\absolutevalue{\delta\phi_{k}(\tau)}^{2}\,.$
(11)

The one-loop can be separated into terms coming from cubic and quartic
interactions. The quartic part is

$\displaystyle\mathcal{P}_{\zeta}^{\rm
1l,V_{4}}(\tau,k)=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\imaginary\left\\{\delta\phi^{2}_{k}(\tau)\int^{\tau}_{-\infty}\differential\tau^{\prime}a^{4}(\tau^{\prime})V_{4}(\tau^{\prime})\delta\phi^{*2}_{k}(\tau^{\prime})\int\dfrac{\differential^{3}\mathbf{p}}{(2\pi)^{3}}\absolutevalue{\delta\phi_{p}(\tau^{\prime})}^{2}\right\\}\,,$
(12)

Since the transitions are instantaneous, the prescription
$\tau^{\prime}_{-}=\tau^{\prime}(1-i\omega)$ has no effect on
$\mathcal{P}_{\zeta}^{\rm 1l,V_{4}}$ and we omit it. The loop integral has no
dependence on the external momentum $k$, so $\mathcal{P}_{\zeta}^{\rm
1l,V_{4}}$ can be fully absorbed by the counterterms, as we shall see below.
The cubic part is

$\displaystyle\mathcal{P}_{\zeta}^{\rm 1l,V_{3}}(\tau,k)$
$\displaystyle=\dfrac{k^{3}}{\pi^{2}M_{P}^{2}\epsilon(\tau)}\imaginary\\{\int_{-\infty}^{\tau+\frac{1}{a(\tau)\Lambda}}\differential\tau^{\prime}a^{4}(\tau^{\prime}_{-})\,V_{3}(\tau^{\prime}_{-})\,\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau^{\prime}_{-})\int_{\tau^{\prime}+\frac{1}{a(\tau^{\prime})\Lambda}}^{\tau}\differential\tau^{\prime\prime}a^{4}(\tau^{\prime\prime})V_{3}(\tau^{\prime\prime})\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau^{\prime\prime})}$
$\displaystyle\quad\quad\quad\quad\times\int\dfrac{\differential^{3}\mathbf{p}}{(2\pi)^{3}}\delta\phi_{p}(\tau^{\prime\prime})\,\delta\phi_{p}^{*}(\tau^{\prime}_{-})\,\delta\phi_{|\mathbf{k}-\mathbf{p}|}(\tau^{\prime\prime})\,\delta\phi_{|\mathbf{k}-\mathbf{p}|}^{*}(\tau^{\prime}_{-})\large\bigg{\\}}\,.$
(13)

In this case, the $i\omega$ prescription is important for the convergence of
the loop integral, and therefore we must keep it.666In principle, we should
also keep this prescription in $\tau^{\prime\prime}$. However, it has no
effect and convergence is guaranteed simply by keeping $\tau^{\prime}_{-}$.
Since the interactions $V_{3}$ are Dirac deltas at the edges of the
transitions, only a few points contribute to the integrals of (13). Let us
determine which ones are relevant. We can denote by $\tau_{1}$ the conformal
time corresponding to $N=0$ in Figure 1. Then, we denote as $\tau_{2}$ and
$\tau_{3}$ the conformal times at which the USR phase starts and ends,
respectively. The beginning of the final SR phase is identified with
$\tau_{4}$. Considering the time integrals in (13), points falling along the
diagonal ($\tau^{\prime}=\tau^{\prime\prime}$) in the plane
$\\{\tau^{\prime},\tau^{\prime\prime}\\}$ are removed by the regulator (see
eq. (6)).777If we used a regulator that does not exclude the points on the
diagonal (red line in Figure 3), such as dimensional regularization, the
result of the finite part of $\mathcal{P}_{\zeta}$ at one-loop would be
different. This difference would be compensated by the finite part of the
counterterms. See Appendix G, where the equivalence between this type of
regulators and a cutoff is shown in the limit $\delta N\to 0$. The pairs
$(\tau_{i},\tau_{j})$ that contribute to the one-loop spectrum are those that
satisfy $\tau^{\prime}=\tau_{i}<\tau^{\prime\prime}=\tau_{j}$, see Figure 3.

Figure 3:  The orange shaded region represents the integration domain of the
(cosmic) time integrals that arises calculating $\mathcal{P_{\zeta}}$ at
second-order in $H_{I}$. The red diagonal line represents
$t^{\prime}=t^{\prime\prime}$. The (thinner) black line parallel to it
restricts the integration regime due to the cutoff, see eq. (6). The
separation between the two lines is given by $\Lambda^{-1}_{\rm UV}$. The
green lines mark the times of the transitions between different phases. Their
intersections within the shaded area correspond to the points (black dots)
contributing to the time integrals in our model. At each point, the
interaction is a Dirac delta in time, which simplifies the integration
greatly.

Taking this into account,

$\displaystyle\mathcal{P}_{\zeta}^{\rm
1l,V_{3}}(\tau,k)=\dfrac{k^{3}}{2\pi^{2}M_{P}^{2}\epsilon(\tau)}\dfrac{H^{2}}{M_{P}^{2}}\sum_{i=1}^{4}\sum_{j=i+1}^{4}\dfrac{a^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{i}}\dfrac{a^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{j}}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{j})}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{i})\,l(k;j,i)}\,,$
(14)

where $\Delta\nu^{2}(\tau_{i})$ is the change in $\nu^{2}$ at $\tau_{i}$ and
$l(k;j,i)$ comes from the loop integral in the second line of (13). We can
write $l(k;j,i)$ conveniently with a change of variables Espinosa  _et al._
(2018),888We define $p=k(s-d)/2$ and $|\mathbf{k}-\mathbf{p}|=k(s+d)/2$ and we
use the fact that eq. (13) is symmetric under the exchange of $\mathbf{p}$ and
$\mathbf{k}-\mathbf{p}$ to integrate over half of the integration domain.

$\displaystyle l(k;j,i)=\dfrac{k^{3}}{16\pi^{2}}\int_{1+\Lambda_{\rm
IR}/k}^{\infty}\differential s\int_{0}^{1}\differential
d\,(s^{2}-d^{2})\,\delta\phi_{k(s-d)/2}(\tau_{j})\delta\phi^{*}_{k(s-d)/2}(\tau_{i-})\delta\phi_{k(s+d)/2}(\tau_{j})\delta\phi^{*}_{k(s+d)/2}(\tau_{i-})\,.$
(15)

We note that $l(k;i,j)$ has no UV divergence, thanks to the $\tau_{i-}$
prescription. Without this prescription, $l$ would be oscillatory in the UV
(see Appendix F). It is also crucial that the contributions from the diagonal
of Figure 3 are excluded by the regulator, as these would introduce a
divergent contribution that cannot be absorbed by the counterterms because of
the different momentum dependence of eq. (14) and the contribution from the
counterterms (see eq. (16) below).999If the points on the diagonal are
included, $l(k;i,i)$ would be real and linearly divergent in the UV. The
structure of the divergence would be $\mathcal{P}_{\zeta}^{\rm 1l,V_{3}}\sim
k^{3}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{i})}^{2}$, while the
structure of the counterterms would be $\mathcal{P}_{\zeta}^{\rm ct}\sim
k^{3}\Im{\delta\phi_{k}^{2}(\tau)\delta\phi_{k}^{*2}(\tau_{i})}$. Previous
works (see e.g. Kristiano and Yokoyama (2022)) obtained UV divergences from a
cubic interaction in the $\zeta$-gauge, whereas we obtain a finite result
thanks to the regulator we use (and the appropriate use of the $i\omega$
prescription). However, $l$ shows an infrared (IR) divergent part ($l^{\rm
IR}$), which we regularize with an IR cutoff $\Lambda_{\rm IR}$. This type of
IR divergence arises from eq. (5) in perturbation theory for massless free
fields. In this work we are only concerned with UV divergences, and assume
that IR divergences are either unphysical and disappear when calculating
physical observables Gerstenlauer  _et al._ (2011); Giddings and Sloth (2011);
Senatore and Zaldarriaga (2013) or can be addressed beyond perturbation theory
Gorbenko and Senatore (2019).101010In the latter case, we will assume that the
finite contribution of the IR effect will not be larger than the finite part
of the rest of the calculation we make. Therefore, in practice, we redefine
$l$ as follows: $l\to l-l^{\rm IR}$. Finally, the contribution of the
counterterms (ct) will be

$\displaystyle\mathcal{P}_{\zeta}^{\rm
ct}(\tau,k)=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\left(2\,\delta_{\phi}\,\absolutevalue{\delta\phi_{k}(\tau)}^{2}+8\Im{\delta\phi_{k}^{2}(\tau)\int_{-\infty}^{\tau}\differential\tau^{\prime}a^{2}(\delta_{\phi}k^{2}+\tilde{\delta}_{V})\delta\phi_{k}^{*2}\evaluated{}_{\tau^{\prime}}}\right)\,.$
(16)

Making the replacement

$\tilde{\delta}_{V}\to\tilde{\delta}_{V}-\dfrac{1}{8}a^{2}V_{4}\int\dfrac{\differential^{3}\mathbf{p}}{(2\pi)^{3}}\absolutevalue{\delta\phi_{p}}^{2}$
(17)

the (divergent and finite) contributions of $\mathcal{P}_{\zeta}^{\rm
1l,V_{4}}$ can be completely absorbed by the counterterms –as it happens in
$\lambda\phi^{4}$ in Minkowski–. We stress that we can make this redefinition
thanks to the arbitrary time dependence of $\tilde{\delta}_{V}$. The situation
is different for $\mathcal{P}_{\zeta}^{\rm 1l,V_{3}}$, where the counterterms
cannot absorb its finite part for all $k$. Therefore, the only relevant
contribution at one-loop –and the only one we include– comes from
$\mathcal{P}_{\zeta}^{\rm 1l,V_{3}}\equiv\mathcal{P}_{\zeta}^{\rm 1l}$. Since
the only UV divergence appearing in the one-loop calculation is reabsorbed
with the counterterms, the two-point correlation is, after this procedure,
completely finite.

We note that in Figure 2 there are diagrams that are not relevant in our
calculation at one-loop. These diagrams fall in one of the following three
categories: 1) bubble diagrams (whose contribution vanishes), 2) diagrams that
do not contribute to $\mathcal{P}_{\zeta}$,111111We have only one diagram in
this category, the disconnected one formed by two tadpoles (last diagram of
the last row of Figure 2). Although this diagram does affect the two-point
correlation, its momentum structure is different from that of the power
spectrum, so it does not modify $\mathcal{P}_{\zeta}$. and 3) diagrams
proportional to the tadpole (which vanishes imposing
$\expectationvalue{\delta\phi}=0$ at loop level, using counterterms, as shown
in Appendix E).

Figure 4:  The dotted curves represent lines of constant $\mathcal{P}_{\zeta}$
(with values $1$, $10^{-2}$, $10^{-4}$ and $10^{-6}$ from top to bottom),
computed at tree-level. In the shaded region, the one-loop contribution to
$\mathcal{P}_{\zeta}$ is larger than its tree-level counterpart at the scales
near the maximum of the latter. The three black dots correspond to the three
examples shown in Figure 5. The region where $\mathcal{P}_{\zeta}^{\rm
1l}>\mathcal{P}_{\zeta}^{\rm tl}$ is approximately scale invariant, as can be
seen in Figure 8.

##  V Discussion

A full renormalization procedure would require imposing a set of conditions on
$\mathcal{P}_{\zeta}$ allowing us to extract the finite part of the
counterterms that make the theoretical prediction of the total power spectrum
coincide with the one inferred from observations, which is currently
unconstrained at small enough scales. Although we have not performed a full
renormalization, which would have required working with a tractable functional
form for the potential, we can use our results to draw conclusions about the
validity of perturbation theory. The validity of eq. (5), which is the basis
of our analysis, assumes $\mathcal{P}_{\zeta}^{\rm
1l}\ll\mathcal{P}_{\zeta}^{\rm tl}$. For consistency, eq. (5) also requires
$\mathcal{P}_{\zeta}^{\rm ct}\sim\mathcal{P}_{\zeta}^{\rm
1l}\ll\mathcal{P}_{\zeta}^{\rm tl}$. The contributions to
$\mathcal{P}_{\zeta}$ coming from the counterterms must be of the same order
as those coming from the loops because otherwise the divergences cannot be
absorbed. Taking the condition $\mathcal{P}_{\zeta}^{\rm
1l}\sim\mathcal{P}_{\zeta}^{\rm tl}$ as a proxy for perturbation theory
breakdown, we find that whether perturbation theory breaks or not depends on
the values chosen for the two parameters, $\delta N$ and $\Delta N$, of the
model. The contribution to $\mathcal{P}_{\zeta}$ at one-loop depends strongly
on the width $\delta N$ of the transitions between SR and USR. Moreover, we
find that when $\delta N$ is made arbitrarily small, the one-loop contribution
diverges (see Appendix G), invalidating any prediction based on perturbation
theory through the in-in formalism. In models that can account for all DM in
the form of PBHs –which requires $\mathcal{P}_{\zeta}\sim\order{0.01}$ at the
peak in the Gaussian approximation, see e.g. Ballesteros and Taoso (2018)–,
$\delta N$ needs to be $\gtrsim 0.1$ for perturbation theory to hold, as shown
in Figures 4 and 5. In realistic USR inflationary models Ballesteros and Taoso
(2018); Ballesteros  _et al._ (2020b) this parameter is $\delta N\sim
0.4-0.5$. For such values of $\delta N$, perturbation theory does indeed
appear to hold in the model we have considered for our analysis. Even though
inflation does not end in our model and $\mathcal{P}_{\zeta}$ does not feature
a peak but rather a plateau-like feature for $k\to\infty$, we think that our
results are evidence that models in which all of the DM is comprised of PBHs
originating from USR inflation are not necessarily hampered by perturbation
theory breaking, contrary to what was argued in Kristiano and Yokoyama (2022).
In future work our analysis can be extended to include
$|\eta|\sim\mathcal{O}(1)$ in the last phase of Figure 1, as well as using
dimensional regularization to deal with UV divergences in
$\mathcal{P_{\zeta}^{\rm 1l}}$. Another possible direction consists in
exploring what happens assuming a given functional form for the potential
$V(\phi)$, such as the one in Ballesteros  _et al._ (2020b).

We note that the Gaussian approximation, and hence the benchmark
$\mathcal{P}_{\zeta}\sim\order{0.01}$, may not be enough to compute the PBH
abundance for realistic USR models accurately. We expect that the effect of
small non-Gaussian corrections around the maximum of the distribution of
$\zeta$ will not change our results significantly Taoso and Urbano (2021).
However, there exist indications that the tail of the distribution in USR
models may fall exponentially, differing substantially from a Gaussian one,
see e.g. Ezquiaga  _et al._ (2020); Figueroa  _et al._ (2021). If this were
the case, a smaller maximum $\mathcal{P}_{\zeta}$ would be sufficient to
account for all DM, which would make attaining the conditions needed for
perturbation theory to hold even easier.

The debate has been intense since Kristiano and Yokoyama (2022) appeared. An
early reaction already suggested the need of smoothing the transition between
the SR and USR phases, see Riotto (2023). More detailed calculations
implementing smooth transitions gave contradictory results, in some cases
supporting the thesis of Kristiano and Yokoyama (2022) (see e.g. Firouzjahi
(2023b)) and in others indicating that the tension could be relaxed, see e.g.
Franciolini  _et al._ (2023). Before our work, previous analyses have
primarily focused on the specific cubic interaction in the $\zeta$-gauge that
was first considered in Kristiano and Yokoyama (2022) ($S\propto\int
d^{4}x\,\epsilon\,\eta^{\prime}\,a^{2}\,\zeta^{\prime}\,\zeta^{2}$). It turns
out that the quartic interaction induced by this cubic term in the action for
$\zeta$, which has been ignored in previous works, is the origin of the
divergence in the limit of instantaneous transitions ($\delta N\to 0$) in this
gauge (as shown in Appendix G). However, the quartic action Firouzjahi (2023b)
and some boundary terms Fumagalli (2023); Tada  _et al._ (2023); Kristiano and
Yokoyama (2023); Firouzjahi (2023a) were also partially considered before our
work.121212In the cubic action in the $\zeta$-gauge Maldacena (2003), there
are several boundary terms arising from integrating by parts, see also Wang
(2014). It is necessary to include all of them in $H_{i}$ to compute
correlation functions correctly. However, the boundary terms without time
derivatives of the fields cannot contribute to any correlator, as it is argued
e.g. in Burrage  _et al._ (2011) using path integral arguments. See Weinberg
(2005) for more details on path integral methods in the in-in formalism.
Although several boundary terms disappear from the cubic action for $\zeta$
after the field redefinition of Maldacena (2003), care must be taken to always
include all relevant contributions. It is worth mentioning the work presented
in Tada  _et al._ (2023), where it was claimed that Maldacena’s consistency
relation Maldacena (2003) and the $i\omega$ prescription guarantee that the
one-loop power spectrum has to be suppressed with respect to the tree-level
one. We have however verified that the consistency relation is satisfied in
our analysis (see Appendix H).

In previous works, the loop integrals (mainly involving the cubic interaction
in the $\zeta$-gauge mentioned above) have been cut between two values of the
momentum (in particular, between the wavenumbers becoming super-Hubble at the
beginning and at the end of USR). That procedure makes the resulting
computation of the one-loop contribution to $\mathcal{P}_{\zeta}$ dependent on
those cutoffs, casting a doubt on any conclusions extracted from it. In our
analysis, we have dealt with this problem regularizing and absorbing the
relevant divergences with appropriate counterterms. We have found that, after
regularization, the cubic interactions in the $\delta\phi$-gauge give a finite
one-loop contribution to $\mathcal{P}_{\zeta}$. However, the quartic
interactions (on the same gauge) required appropriate counterterms to absorb
divergent contributions.

Figure 5:  Tree-level (red dashed) and absolute value of renormalized one-loop
(black continuous) contributions to $\mathcal{P}_{\zeta}(k)$ for three
different choices of $\delta N$ and $\Delta N$. From left to right: $\\{\delta
N=0.1\,,\Delta N=2.07\\}$, $\\{\delta N=0.2\,,\Delta N=1.99\\}$, $\\{\delta
N=0.5\,,\Delta N=1.77\\}$. The vertical green bands indicate the ranges of
momenta becoming super-Hubble during the equally-shaded transitions of Figure
1. The value $k=k_{0}$ corresponds to the Fourier mode becoming super-Hubble
at the beginning of the first transition (between SR and USR).

Acknowledgments.

This work has been funded by the following grants: 1) Contrato de Atracción de
Talento (Modalidad 1) de la Comunidad de Madrid (Spain), 2017-T1/TIC-5520 and
2021-5A/TIC-20957, 2) PID2021-124704NB-I00 funded by
MCIN/AEI/10.13039/501100011033 and by ERDF A way of making Europe, 3)
CNS2022-135613 MICIU/AEI/10.13039/501100011033 and by the European Union
NextGenerationEU/PRTR, 4) the IFT Centro de Excelencia Severo Ochoa Grant
CEX2020-001007-S, funded by MCIN/AEI/10.13039/501100011033. JGE has been
supported by a PhD contract contrato predoctoral para formación de doctores
PRE2021-100714 associated to the aforementioned Severo Ochoa program,
CEX2020-001007-S-21-3. We thank Gabriele Franciolini and Thomas Konstadin for
comments and Marco Taoso for questions. We thank Jaime Fernández Tejedor for
comments and discussions. We thank José Ramón Espinosa for discussions.

##  Appendix A In-in formalism in more detail

Let $Q(t)$ be a generic Hermitian operator whose vacuum expectation value
(VEV) for which we want to calculate

$\expectationvalue{Q(t)}\equiv\sideset{{}_{H}}{}{\mathop{\bra{\Omega}}}Q_{H}(t)\ket{\Omega}_{H},$
(18)

where $\ket{\Omega}_{H}$ is the ground state, i.e. the vacuum of the system.
This is customarily done in the so-called in-in formalism; see Schwinger
(1961); Keldysh (1964) and also Weinberg (2005); Chen (2010); Wang (2014) for
later descriptions. We work in the Heisenberg (H) picture, in which
$\ket{\Omega}_{H}$ does not evolve in time, and instead all the time evolution
falls on the operators, such as $Q_{H}(t)$. The time evolution of the
operators in this picture is given by the equation

$Q_{H}(t)=U^{-1}(t,t_{0})Q_{S}\,U(t,t_{0})\,,$ (19)

where $U^{-1}(t,t_{0})U(t,t_{0})=1$ . Here we are introducing an auxiliary
time $t_{0}$ in which all objects in the Heisenberg and Schrödinger (S)
picture coincide. $U(t,t_{0})$ is the time evolution operator that satisfies

$\partial_{t}U(t,t_{0})=-iH(t)U(t,t_{0})\,,\quad U(t_{0},t_{0})=1\,.$ (20)

The Hamiltonian $H(t)$ can be separated into a free part and an interaction
part, $H=H_{0}+H_{\rm i}$. Thus, we can define the operators in the
interaction picture (I)

$\displaystyle
Q_{I}(t)=U_{0}^{-1}(t,t_{0})Q_{S}\,U_{0}(t,t_{0})\,,\quad\text{where}\quad\partial_{t}U_{0}(t,t_{0})=-iH_{0}(t)U_{0}(t,t_{0}),\;U_{0}(t_{0},t_{0})=1\,.$
(21)

The operators in the interaction picture follow the dynamics governed by the
free (quadratic) Hamiltonian, without interactions. That is, one can
understand the operators in the interaction picture as operators in the
Heisenberg picture, but with the Hamiltonian of the system being simply the
free Hamiltonian. Therefore, the fields in the interaction picture can be
expressed in terms of the creation and annihilation operators –which will
allow us to use Wick’s theorem when computing correlators–. We can now write
the operators in the Heisenberg picture in terms of the operators in the
interaction picture,

$\displaystyle Q_{H}(t)=F^{-1}(t,t_{0})Q_{I}(t)F(t,t_{0}),\quad
F(t,t_{0})\equiv U_{0}^{-1}(t,t_{0})U(t,t_{0})\,.$ (22)

Differentiating, we obtain the equation that is satisfied by the operator
$F(t,t_{0})$:

$\displaystyle\partial_{t}F(t,t_{0})=-iH_{I}(t)F(t,t_{0})\,,$ (23)

where the Hermitian operator

$\displaystyle H_{I}(t)\equiv U_{0}^{-1}(t,t_{0})H_{\rm i}(t)U_{0}(t,t_{0})$
(24)

represents the interaction Hamiltonian in the interaction picture, and then
the fields that compose $H_{I}$ are in the interaction picture. The formal
solution of eq. (23) and its inverse are, respectively:

$\displaystyle F(t,t_{0})=T\exp\left(-i\int_{t_{0}}^{t}\differential
t^{\prime}H_{I}(t^{\prime})\right)\,,\quad
F^{-1}(t,t_{0})=\bar{T}\exp\left(i\int_{t_{0}}^{t}\differential
t^{\prime}H_{I}(t^{\prime})\right)\,,$ (25)

where $T$ and $\bar{T}$ indicate time and anti-time ordering. The VEV of
$Q(t)$ is then

$\expectationvalue{Q(t)}=\sideset{{}_{H}}{}{\mathop{\bra{\Omega}}}F^{-1}(t,t_{0})Q_{I}(t)F(t,t_{0})\ket{\Omega}_{H}\,.$
(26)

We note that all the elements of the latter expression are in the interaction
picture except for $\sideset{{}_{H}}{}{\mathop{\bra{\Omega}}}$ and
$\ket{\Omega}_{H}$. Therefore, we need to rewrite this expression using the
vacuum in the interaction picture $\ket{0}_{I}$ (i.e. the vacuum of the
Hamiltonian without interactions). Like $\ket{\Omega}_{H}$, $\ket{0}_{I}$ does
not evolve over time. Let us consider

$\displaystyle\exp\left(-i\int_{t_{0}}^{t}\differential
t^{\prime}H(t^{\prime})\right)\ket{0}_{I}=\exp\left(-i\int_{t_{0}}^{t}\differential
t^{\prime}H(t^{\prime})\right)\sum_{n}\ket{n}_{H}\sideset{{}_{H}}{{}_{I}}{\mathop{\innerproduct{n}{0}}}$
$\displaystyle=\exp\left(-i\int_{t_{0}}^{t}\differential
t^{\prime}E_{\Omega}(t^{\prime})\right)\left(\sideset{{}_{H}}{{}_{I}}{\mathop{\innerproduct{\Omega}{0}}}\ket{\Omega}_{H}+\sum_{n\neq\Omega}\exp\left(-i\int_{t_{0}}^{t}\differential
t^{\prime}\left(E_{n}(t^{\prime})-E_{\Omega}(t^{\prime})\right)\right)\>\sideset{{}_{H}}{{}_{I}}{\mathop{\innerproduct{n}{0}}}\ket{n}_{H}\right)\,.$
(27)

Since $E_{n}(t)>E_{\Omega}(t)$ for all $n$ and $t$, if we make the replacement
$t_{0}\to-\infty^{-}\equiv-\infty(1-i\omega)$ with $\omega>0$, all excited
states are suppressed and we get

$\exp\left(-i\int_{-\infty^{-}}^{t}\differential
t^{\prime}H(t^{\prime})\right)\ket{0}_{I}=\sideset{{}_{H}}{{}_{I}}{\mathop{\innerproduct{\Omega}{0}}}\exp\left(-i\int_{-\infty^{-}}^{t}\differential
t^{\prime}H(t^{\prime})\right)\ket{\Omega}_{H}.$ (28)

Applying time ordering,

$U(t,-\infty^{-})\ket{\Omega}_{H}=\dfrac{U(t,-\infty^{-})\ket{0}_{I}}{\sideset{{}_{H}}{{}_{I}}{\mathop{\innerproduct{\Omega}{0}}}}\,,$
(29)

and we obtain the following useful relation between the vacuum in the
interaction and Heisenberg picture:131313We interpret the prescription
$i\omega$ as a way to turn off the interactions in the limit $t\to-\infty$,
well motivated by the Bunch-Davies initial conditions.

$F(t,-\infty^{-})\ket{\Omega}_{H}=\dfrac{F(t,-\infty^{-})\ket{0}_{I}}{\sideset{{}_{H}}{{}_{I}}{\mathop{\innerproduct{\Omega}{0}}}}\equiv\mathcal{N}\,F(t,-\infty^{-})\ket{0}_{I}\quad\textrm{and}\quad\sideset{{}_{H}}{}{\mathop{\bra{\Omega}}}F^{-1}(t,-\infty^{+})=\mathcal{N}^{*}\sideset{{}_{I}}{}{\mathop{\bra{0}}}F^{-1}(t,-\infty^{+})\,,$
(30)

where $\infty^{+}\equiv\infty(1+i\omega)$. This allows us to write all the
elements of the VEV in the interaction picture:

$\expectationvalue{Q(t)}=\absolutevalue{\mathcal{N}}^{2}\sideset{{}_{I}}{}{\mathop{\bra{0}}}F^{-1}(t,-\infty^{+})Q_{I}(t)F(t,-\infty^{-})\ket{0}_{I}\,,$
(31)

where $F^{-1}(t,-\infty^{+})$ can be read from (25). The last piece we need to
understand is what is $\absolutevalue{\mathcal{N}}^{2}$. Assuming that
$\ket{\Omega}_{H}$ is properly normalized, i.e.
$\sideset{{}_{H}}{{}_{H}}{\mathop{\innerproduct{\Omega}{\Omega}}}=1$, if we
calculate the VEV of the identity operator ($Q=1$), we get that

$\frac{1}{\absolutevalue{\mathcal{N}}^{2}}=\sideset{{}_{I}}{}{\mathop{\bra{0}}}F^{-1}(t,-\infty^{+})F(t,-\infty^{-})\ket{0}_{I}\,.$
(32)

Making a perturbative expansion, it becomes clear that
$\absolutevalue{\mathcal{N}}^{-2}$ is the sum over the bubble-type diagrams,
defined as the set of diagrams with at least one part disconnected from any
external leg. Thus, the final formula to be used for the calculation of the
VEV of any operator $Q(t)$ is (see also Senatore (2017))

$\expectationvalue{Q(t)}=\dfrac{\sideset{{}_{I}}{}{\mathop{\bra{0}}}F^{-1}(t,-\infty^{+})Q_{I}(t)F(t,-\infty^{-})\ket{0}_{I}}{\sideset{{}_{I}}{}{\mathop{\bra{0}}}F^{-1}(t,-\infty^{+})F(t,-\infty^{-})\ket{0}_{I}}=\sideset{{}_{I}}{}{\mathop{\bra{0}}}F^{-1}(t,-\infty^{+})Q_{I}(t)F(t,-\infty^{-})\ket{0}_{I}\evaluated{}_{\textrm{no
bubbles}}\,.$ (33)

For brevity, we write $\ket{0}_{I}$ as $\ket{0}$ in what follows (and in the
main body of this work) and, even if we will not write it explicitly, bubble
diagrams will be assumed to be removed.

The main expression we use in this paper, eq. (5), is a second order expansion
in the interaction Hamiltonian of the in-in formalism for the calculation of a
VEV. To derive it we proceed by first expanding $F(t,-\infty^{-})$ from eq.
(25):

$\displaystyle F(t,-\infty^{-})=1-i\int_{-\infty^{-}}^{t}\differential
t^{\prime}H_{I}(t^{\prime})-\dfrac{1}{2}T\int_{-\infty^{-}}^{t}\differential
t^{\prime}\int_{-\infty^{-}}^{t}\differential
t^{\prime\prime}H_{I}(t^{\prime})H_{I}(t^{\prime\prime})\,.$ (34)

Now, we make a redefinition of time such that the prescription $i\omega$
passes from the lower integration limit to the time of the Hamiltonian,
$t\rightarrow t_{-}=t(1-i\omega)$, such that

$\displaystyle F(t,-\infty^{-})=1-i\int_{-\infty}^{t}\differential
t^{\prime}H_{I}(t^{\prime}_{-})-\dfrac{1}{2}T\int_{-\infty}^{t}\differential
t^{\prime}\int_{-\infty}^{t}\differential
t^{\prime\prime}H_{I}(t^{\prime}_{-})H_{I}(t^{\prime\prime}_{-})\,.$ (35)

In principle, the upper integration limit should also be modified, but this
change has no effect when, at the end of the calculation, we take $\omega\to
0$. The only effect of the prescription $i\omega$ is observed in the limit
$t^{\prime},t^{\prime\prime}\to-\infty$, introducing a damping in the time
integrals, which is interpreted as a projection onto the interaction vacuum.
Now we use the definition of time ordering to rewrite

$\dfrac{1}{2}T\int_{-\infty}^{t}\differential
t^{\prime}\int_{-\infty}^{t}\differential
t^{\prime\prime}H_{I}(t^{\prime}_{-})H_{I}(t^{\prime\prime}_{-})=\int_{-\infty}^{t}\differential
t^{\prime}\int_{-\infty}^{t^{\prime}}\differential
t^{\prime\prime}H_{I}(t^{\prime}_{-})H_{I}(t^{\prime\prime}_{-})=\int_{-\infty}^{t}\differential
t^{\prime}\int_{t^{\prime}}^{t}\differential
t^{\prime\prime}H_{I}(t^{\prime\prime}_{-})H_{I}(t^{\prime}_{-})\,.$ (36)

Therefore, $F$ and its inverse are

$\displaystyle F(t,-\infty^{-})$
$\displaystyle=1-i\int_{-\infty}^{t}\differential
t^{\prime}H_{I}(t^{\prime}_{-})-\int_{-\infty}^{t}\differential
t^{\prime}\int_{t^{\prime}}^{t}\differential
t^{\prime\prime}H_{I}(t^{\prime\prime}_{-})H_{I}(t^{\prime}_{-})\,,$ (37)
$\displaystyle F^{-1}(t,-\infty^{+})$
$\displaystyle=1+i\int_{-\infty}^{t}\differential
t^{\prime}\left(H_{I}(t^{\prime}_{-})\right)^{\dagger}-\int_{-\infty}^{t}\differential
t^{\prime}\int_{t^{\prime}}^{t}\differential
t^{\prime\prime}\left(H_{I}(t^{\prime\prime}_{-})H_{I}(t^{\prime}_{-})\right)^{\dagger}\,,$
(38)

so that the VEV of $Q$ can be written as

$\displaystyle\expectationvalue{Q(t)}$
$\displaystyle=\bra{0}Q_{I}(t)\ket{0}+2\Im{\int_{-\infty}^{t}\differential
t^{\prime}\bra{0}Q_{I}(t)H_{I}(t^{\prime}_{-})\ket{0}}+\int_{-\infty}^{t}\differential
t^{\prime}\int_{-\infty}^{t}\differential
t^{\prime\prime}\bra{0}H_{I}(t^{\prime\prime}_{+})Q_{I}(t)H_{I}(t^{\prime}_{-})\ket{0}$
$\displaystyle\quad-2\Re{\int_{-\infty}^{t}\differential
t^{\prime}\int_{t^{\prime}}^{t}\differential
t^{\prime\prime}\bra{0}Q_{I}(t)H_{I}(t^{\prime\prime}_{-})H_{I}(t^{\prime}_{-})\ket{0}}\,.$
(39)

Using that

$\bra{0}H_{I}(t^{\prime}_{+})Q_{I}(t)H_{I}(t^{\prime\prime}_{-})\ket{0}=\bra{0}H_{I}(t^{\prime\prime}_{+})Q_{I}(t)H_{I}(t^{\prime}_{-})\ket{0}^{*}$
(40)

we can rewrite the last term of the first line of (39) as

$\int_{-\infty}^{t}\differential t^{\prime}\int_{-\infty}^{t}\differential
t^{\prime\prime}\bra{0}H_{I}(t^{\prime\prime}_{+})Q_{I}(t)H_{I}(t^{\prime}_{-})\ket{0}=2\Re{\int_{-\infty}^{t}\differential
t^{\prime}\int_{t^{\prime}}^{t}\differential
t^{\prime\prime}\bra{0}H_{I}(t^{\prime\prime}_{+})Q_{I}(t)H_{I}(t^{\prime}_{-})\ket{0}}\,,$
(41)

so that the final expression is

$\displaystyle\expectationvalue{Q(t)}$
$\displaystyle=\bra{0}Q_{I}(t)\ket{0}+2\Im{\int_{-\infty}^{t}\differential
t^{\prime}\bra{0}Q_{I}(t)H_{I}(t^{\prime}_{-})\ket{0}}$
$\displaystyle\quad+2\Re{\int_{-\infty}^{t}\differential
t^{\prime}\int_{t^{\prime}}^{t}\differential
t^{\prime\prime}\bra{0}\left(H_{I}(t^{\prime\prime}_{+})Q_{I}(t)-Q_{I}(t)H_{I}(t^{\prime\prime}_{-})\right)H_{I}(t^{\prime}_{-})\ket{0}}\,.$
(42)

We emphasize that this expression does not include the bubble diagrams and we
stress the importance of the prescription $i\omega$ (represented by $t_{\pm}$)
for the convergence of the time integrals.

##  Appendix B Action for fluctuations in the $\delta\phi$-gauge

In the ADM formalism Arnowitt  _et al._ (2008); Maldacena (2003); Wang (2014)
the spacetime metric is written as follows:

$\displaystyle ds^{2}=-N^{2}\differential
t^{2}+\gamma_{ij}\left(N^{i}\differential t+\differential
x^{i}\right)\left(N^{j}\differential t+\differential x^{j}\right)\,,$ (43)

where the lapse $N$ and the shift $N^{i}$ are non-dynamical quantities acting
as Lagrange multipliers in the action, and are therefore determined in terms
of the dynamical quantities of the system. We use the scalar-vector-tensor
decomposition of the metric

$\displaystyle\gamma_{ij}=a^{2}\left[e^{\Gamma}\right]_{ij}\,,\quad\Gamma_{ij}=2\zeta\delta_{ij}+\partial_{ij}E+\partial_{(i}E_{j)}+h_{ij}\,,\quad
N=1+\alpha\,,\quad N_{i}=\partial_{i}\beta+\beta_{i}\,,$ (44)

where
$\partial^{i}E_{i}=\partial^{i}\beta_{i}=\partial^{i}h_{ij}=\tensor{h}{{}^{i}_{i}}=0$.
We also decompose the inflaton field $\phi$ into a homogeneous background
$\phi_{0}$ and fluctuations:

$\displaystyle\phi(t,x)=\phi_{0}(t)+\delta\phi(t,x)\,.$ (45)

Doing an infinitesimal coordinate (gauge) transformation we can eliminate two
scalar and two vector degrees of freedom, leaving only one scalar and two
tensor physical degrees of freedom. Although we are primarily interested in
scalar modes, for completeness we will also include tensor modes.141414By
including the tensor modes, we can obtain the action describing the
gravitational waves induced at second order by scalar fluctuations. One
possible choice of gauge is: $\delta\phi=E=0$, which is known as the
$\zeta$-gauge. The variable $\zeta$, defined via eq. (44), eventually becomes
constant for $k\ll H$ Maldacena (2003); Weinberg (2003), which allows us to
relate it to observable quantities in a simple way, see e.g. Baumann (2022).
However, the one-loop calculation of $\mathcal{P}_{\zeta}$ in the model we are
interested in is simplified if we work in the $\delta\phi$-gauge, where
$\zeta=E=0$. This is because in the inflationary model we consider, the scalar
interactions coming from the Ricci scalar $R$ in the action

$S=\dfrac{1}{2}\int\sqrt{-g}\>\differential^{4}x\bigg{\\{}M_{P}^{2}R-g^{\mu\nu}\partial_{\mu}\phi\partial_{\nu}\phi-2V(\phi)\bigg{\\}}\,,$
(46)

are suppressed with respect to the interactions coming from the potential
$V(\phi)$. The $\delta\phi$-gauge allows us to study simultaneously cubic and
quartic interactions, avoiding writing interactions in the form of boundary
terms that may become relevant in the loop calculation Arroja and Tanaka
(2011); Fumagalli (2023).151515In ref. Arroja and Tanaka (2011), the role of
the boundary terms in the tree-level calculation of the bispectrum of $\zeta$
is analyzed, concluding that there is only one relevant boundary term and that
it is cancelled by the field redefinition introduced by Maldacena in Maldacena
(2003). However, in the calculation of other observable quantities, and
specifically at loop level, there is a larger set of relevant boundary terms
that needs to be included in the calculations. This difficulty is circumvented
using the $\delta\phi$-gauge.

Let us see which are the most relevant interactions of the action (46) in the
$\delta\phi$-gauge under the assumption that $\epsilon\equiv-\dot{H}/H^{2}\ll
1$ and $\eta\equiv\dot{\epsilon}/(H\epsilon)\sim 1$ at most. We write the
components of the metric and its inverse as

$\displaystyle g_{\mu\nu}=\left(\begin{array}[]{cc}-N^{2}+N^{k}N_{k}&N_{i}\\\
N_{j}&\gamma_{ij}\end{array}\right),\quad
g^{\mu\nu}=\left(\begin{array}[]{cc}-1/N^{2}&N^{i}/N^{2}\\\
N^{j}/N^{2}&\gamma^{ij}-N^{i}N^{j}/N^{2}\end{array}\right).$ (51)

The Ricci scalar decomposes according to

$\displaystyle
R=R^{(3)}-\dfrac{1}{N^{2}}(K^{2}-K^{ij}K_{ij})+\textrm{{boundary
term}}\,,\quad
K_{ij}=\dfrac{1}{2}\left(\dot{\gamma}_{ij}-\nabla_{i}N_{j}-\nabla_{j}N_{i}\right)\,,\quad
K=\tensor{K}{{}^{i}_{i}}\,.$ (52)

We use that the spatial indices (which we have been denoting with Latin
letters) are risen and lowered with $\gamma_{ij}$, so that
$N^{i}=\gamma^{ij}N_{j}$, where $\gamma^{ij}$ is the inverse of $\gamma_{ij}$.
In our notation, $R^{(3)}$ is the Ricci scalar for $\gamma_{ij}$, $\nabla_{i}$
is the covariant derivative associated to $\gamma_{ij}$ and $N^{-1}K_{ij}$ is
the extrinsic curvature. In addition, we now have $\sqrt{-g}=N\sqrt{\gamma}$.
Therefore, the action (46) is rewritten as

$\displaystyle
S=\int\differential^{4}xN\sqrt{\gamma}\left\\{\dfrac{M_{P}^{2}}{2}\left(R^{(3)}-\dfrac{1}{N^{2}}\left(K^{2}-K^{ij}K_{ij}\right)\right)+\dfrac{1}{2N^{2}}(\dot{\phi}-N^{i}\partial_{i}\phi)^{2}-\dfrac{1}{2}\gamma^{ij}\partial_{i}\phi\partial_{j}\phi-V(\phi)\right\\},$
(53)

where we have eliminated the boundary terms so that the variational principle
is well defined Gibbons and Hawking (1977); York (1972); Dyer and
Hinterbichler (2009). This action is valid at any order in perturbation
theory. As anticipated, the lapse $N$ and the shift $N_{i}$ have no time
derivatives, so these fields behave as Lagrange multipliers. Varying the
action with respect to them, we obtain that

$\displaystyle
M_{P}^{2}\left(\partial_{j}\dfrac{K}{N}-\gamma_{jk}\nabla_{i}\dfrac{K^{ik}}{N}\right)+\dfrac{1}{N}(\dot{\phi}-N^{i}\partial_{i}\phi)\partial_{j}\phi$
$\displaystyle=0\,,$ (54)
$\displaystyle\dfrac{M_{P}^{2}}{2}\left(R^{(3)}+\dfrac{1}{N^{2}}\left(K^{2}-K^{ij}K_{ij}\right)\right)-\dfrac{1}{2N^{2}}(\dot{\phi}-N^{i}\partial_{i}\phi)^{2}-\dfrac{1}{2}\gamma^{ij}\partial_{i}\phi\partial_{j}\phi-V(\phi)$
$\displaystyle=0\,.$ (55)

where

$\displaystyle\nabla_{i}\dfrac{K^{ik}}{N}\equiv\partial_{i}\dfrac{K^{ik}}{N}+\dfrac{K^{ik}}{N}\Gamma^{j}_{ij}+\dfrac{K^{ij}}{N}\Gamma^{k}_{ij}\,.$
(56)

By solving these algebraic equations, we can reintroduce $N$ and $N_{i}$ into
the action, eliminating the Lagrange multipliers. We are interested in the
action in the $\delta\phi$-gauge up to order four in powers of $\delta\phi$
and $h_{ij}$. In order to get there, we start by studying next the properties
of the potential, the lapse and the shift.

###  B.1 Expansion of the potential

The interaction terms coming from the potential are obtained by making the
expansion

$V(\phi)=V(\phi_{0}+\delta\phi)=\sum_{n=0}^{\infty}V_{n}(\phi_{0})\dfrac{\delta\phi^{n}}{n!}\;\;\textrm{where}\;\;V_{n}(\phi_{0})=\dfrac{\partial^{n}V(\phi_{0})}{\partial\phi_{0}^{n}}=\dfrac{1}{\dot{\phi_{0}}}\dfrac{\differential
V_{n-1}(\phi_{0})}{\differential t}\,.$ (57)

Using that $\dot{\phi_{0}}^{2}=2H^{2}M_{P}^{2}\epsilon$ and
$V_{0}=H^{2}M_{P}^{2}(3-\epsilon)$ and defining $\epsilon=-\dot{H}/H^{2}$,
$\eta=\dot{\epsilon}/(H\epsilon)=\epsilon_{2}$ and, in general,
$\epsilon_{n}=\dot{\epsilon}_{n-1}/(H\epsilon_{n-1})$ for $n\geq 2$, the
derivatives of the potential are

$\displaystyle V_{1}$
$\displaystyle=\frac{H^{2}M_{P}\sqrt{\epsilon}}{\sqrt{2}}(-6+2\epsilon-\eta)\,,$
(58) $\displaystyle V_{2}$
$\displaystyle=-\frac{1}{4}H^{2}\left(8\epsilon^{2}-2\epsilon(12+5\eta)+\eta(6+\eta+2\epsilon_{3})\right)\simeq-H^{2}\left(\nu^{2}-\dfrac{9}{4}\right)\,,$
(59) $\displaystyle V_{3}$
$\displaystyle=\frac{H^{2}}{2\sqrt{2}\sqrt{\epsilon}M_{P}}\left(8\epsilon^{3}-6\epsilon^{2}(4+3\eta)+\epsilon\eta(18+6\eta+7\epsilon_{3})-\eta\epsilon_{3}(3+\epsilon_{4}+\eta+\epsilon_{3})\right)\simeq-\dfrac{H}{\sqrt{2\epsilon}M_{P}}\dot{\left(\nu^{2}\right)}\,,$
(60) $\displaystyle V_{4}$ $\displaystyle=-\frac{H^{2}}{8\epsilon
M_{P}^{2}}\left(32\epsilon^{4}-16\epsilon^{3}(6+7\eta)+2\epsilon^{2}\eta(72+39\eta+32\epsilon_{3})-\epsilon\eta\left(6\eta^{2}+\eta(18+35\epsilon_{3})+6\epsilon_{3}(8+3\epsilon_{4}+3\epsilon_{3})\right)\right.$
$\displaystyle\left.\quad+\eta\epsilon_{3}\left(-\eta^{2}+\eta(-3+\epsilon_{4}+3\epsilon_{3})+2\epsilon_{4}(3+\epsilon_{4}+\epsilon_{5})+6\epsilon_{3}(1+\epsilon_{4})+2\epsilon_{3}^{2}\right)\right)\simeq-\dfrac{1}{2\epsilon
M_{P}^{2}}\left(\ddot{\left(\nu^{2}\right)}-\dfrac{H\eta}{2}\dot{\left(\nu^{2}\right)}\right)\,,$
(61)

where

$\displaystyle\nu^{2}\equiv\frac{9}{4}+\frac{1}{2}\left(3\,\eta+\frac{\eta^{2}}{2}+\frac{\eta^{\prime}}{a\,H}\right)\in\mathbb{R}\,,$
(62)

and we recall that primes denote derivatives with respect to conformal time
and dots to cosmic time. In the approximations leading to the above
expressions, we have applied that $\epsilon\ll 1$ and $\epsilon_{n}\sim 1$ for
$n\geq 2$. By changing to conformal time we recover the expressions given in
eq. (2). We note that $V_{n}\sim\epsilon^{1-n/2}$ for $n\geq 1$. This implies
that in the limit $\epsilon\ll 1$ the interactions coming from the potential
dominate over those coming from the purely gravitational part of the action,
since these are not accompanied by this enhancement, as we shall see below.

###  B.2 Expansions of $N$ and $N_{i}$

Since we want to calculate the action up to fourth order, we must obtain the
fluctuations of $N$ and $N_{i}$ up to second order. This is because when
calculating the action at order $n$, the fluctuations at order $n$ and $n-1$
of both the lapse and the shift cancel out Wang (2014). We define

$N=1+\alpha=1+\alpha_{1}+\alpha_{2}+\dots\,,\quad
N_{i}=\partial_{i}\beta+\beta_{i}=\partial_{i}\beta_{1}+\partial_{i}\beta_{2}+\beta_{1,i}+\beta_{2,i}+\dots,$
(63)

where the subscript n (e.g. in $\alpha_{n}$) refers to the perturbative order.
In this way, we seek to obtain $\alpha$, $\beta$ and $\beta_{i}$, for which we
only have to use the eqs. (54) and (55) and expand them perturbatively.
Starting with the first order, we obtain

$\alpha_{1}=\sqrt{\dfrac{\epsilon}{2}}\dfrac{\delta\phi}{M_{P}}\,,\quad\dfrac{1}{a^{2}}\partial^{2}\beta_{1}=\dfrac{\sqrt{\epsilon}}{2\sqrt{2}M_{P}}\left(H\eta\delta\phi-2\dot{\delta\phi}\right)\,,\quad\beta_{1,i}=0.$
(64)

We stress that $\alpha_{1}$ and $\beta_{1}$ are order $\sqrt{\epsilon}$. We
now study the lapse and the shift at second order,

$\displaystyle 0=\,$ $\displaystyle
2H^{2}M_{P}^{2}(\epsilon-3)\alpha_{2}-\dfrac{2HM_{P}^{2}}{a^{2}}\partial^{2}\beta_{2}-\dfrac{2-\epsilon}{4}\dot{\delta\phi}^{2}-\dfrac{H\epsilon\eta}{4}\delta\phi\dot{\delta\phi}$
$\displaystyle+\dfrac{1}{2a^{2}}\left(2\sqrt{2\epsilon}HM_{P}\partial^{i}\beta_{1}-\partial^{i}\delta\phi\right)\partial_{i}\delta\phi+\dfrac{H^{2}}{16}\left(-8\epsilon^{2}+\epsilon\left(24-4\eta+\eta^{2}\right)+8\left(\nu^{2}-\dfrac{9}{4}\right)\right)\delta\phi^{2}$
$\displaystyle-\dfrac{M_{P}^{2}}{8}\left(\dot{h}^{ij}\dot{h}_{ij}+\dfrac{1}{a^{2}}\partial^{k}h^{ij}\partial_{k}h_{ij}\right)+\dfrac{M_{P}^{2}}{2a^{2}}\left(4Hh^{ij}+\dot{h}^{ij}-\dfrac{1}{a^{2}}\partial^{ij}\beta_{1}\right)\partial_{ij}\beta_{1}\,,$
(65) $\displaystyle 0=\,$
$\displaystyle\dfrac{M_{P}^{2}}{2a^{2}}\partial^{2}\beta_{2,j}-2HM_{P}^{2}\partial_{j}\alpha_{2}+\dfrac{M_{P}^{2}}{4}\left(\dot{h}^{ik}\left(\partial_{j}h_{ik}-\partial_{k}h_{ij}\right)+h^{ik}\partial_{k}\dot{h}_{ij}\right)-\dfrac{M_{P}^{2}}{2a^{2}}\partial^{2}h_{ij}\partial^{i}\beta_{1}+\dfrac{\sqrt{\epsilon}M_{P}}{2\sqrt{2}}\dot{h}_{ij}\partial^{i}\delta\phi$
$\displaystyle-\dfrac{\sqrt{\epsilon}M_{P}}{\sqrt{2}a^{2}}\partial^{i}\delta\phi\partial_{ij}\beta_{1}+\dfrac{1}{4}\left(H\epsilon\left(4+\eta\right)\delta\phi-2\left(\epsilon-2\right)\dot{\delta\phi}\right)\partial_{j}\delta\phi\,.$
(66)

Even though the expressions for the lapse and the shift fluctuations at second
order are cumbersome, since our goal is to obtain the action in the limit
$\epsilon\ll 1$ and $\eta\sim\mathcal{O}(1)$ (at most), it is enough to note
that
$\alpha_{2}\sim\beta_{2}\sim\beta_{2,i}\sim\delta\phi^{2}+\sqrt{\epsilon}\,h\,\varphi+h^{2}$.

###  B.3 Quadratic, cubic and quartic actions for fluctuations

To obtain the quadratic, cubic and quartic actions, we have to introduce the
elements studied above in the action (53) and make a perturbative expansion in
powers of the fields. In addition, we will systematically apply integration by
parts with spatial derivatives, but not with time derivatives, due to known
subtleties regarding boundary terms Arroja and Tanaka (2011) (which affect the
n-point correlation functions). The quadratic action is then

$\displaystyle S_{2}$
$\displaystyle=\int\differential^{4}x\,\bigg{\\{}\dfrac{M_{P}^{2}a^{3}}{8}\left(\dot{h}^{ij}\dot{h}_{ij}-\dfrac{1}{a^{2}}\partial^{k}h^{ij}\partial_{k}h_{ij}\right)+\dfrac{a^{3}}{2}\left(\dot{\delta\phi}^{2}-\dfrac{1}{a^{2}}\partial^{i}\delta\phi\partial_{i}\delta\phi+\dfrac{2}{a^{3}}\partial_{t}\left(Ha^{3}\epsilon\right)\delta\phi^{2}-V_{2}\delta\phi^{2}\right)$
$\displaystyle\quad\quad\quad\quad-\dfrac{1}{2}\partial_{t}\left(Ha^{3}\epsilon\delta\phi^{2}\right)\bigg{\\}}$
(67)
$\displaystyle\simeq\int\differential^{4}x\,\bigg{\\{}\dfrac{M_{P}^{2}a^{3}}{8}\left(\dot{h}^{ij}\dot{h}_{ij}-\dfrac{1}{a^{2}}\partial^{k}h^{ij}\partial_{k}h_{ij}\right)+\dfrac{a^{3}}{2}\left(\dot{\delta\phi}^{2}-\dfrac{1}{a^{2}}\partial^{i}\delta\phi\partial_{i}\delta\phi-
V_{2}\delta\phi^{2}\right)\bigg{\\}}\,.$ (68)

In the second line we have taken the limit $\epsilon\ll 1$ to keep only the
terms we are interested in. As we have commented above, the self-interactions
of $\delta\phi$ coming from the purely gravitational part of the action are
suppressed with respect to those coming from the potential. Therefore, at this
level, the action for the scalar field is simply the action of a free scalar
field in a FLRW universe.

The complete cubic action (without any approximations) is

$\displaystyle S_{3}$
$\displaystyle=\int\differential^{4}x\,\Bigg{\\{}\dfrac{M_{P}^{2}a}{24}h^{ij}\left(4h^{kl}\partial_{lj}h_{ik}-3h^{kl}\partial_{kl}h_{ij}+2\partial_{k}h_{jl}\partial^{l}\tensor{h}{{}_{i}^{k}}\right)-\dfrac{M_{P}a^{3}\sqrt{\epsilon}}{8\sqrt{2}}\delta\phi\left(\dot{h}^{ij}\dot{h}_{ij}+\dfrac{1}{a^{2}}\partial^{k}h^{ij}\partial_{k}h_{ij}\right)$
$\displaystyle\quad\quad\quad\quad-\dfrac{1}{4}M_{P}^{2}a\dot{h}^{jk}\partial_{i}h_{jk}\partial^{i}\beta_{1}+\dfrac{M_{P}a\sqrt{\epsilon}}{2\sqrt{2}}\dot{h}^{ij}\delta\phi\partial_{ij}\beta_{1}+\dfrac{M_{P}^{2}}{2a}\partial_{i}h_{jk}\partial^{i}\beta_{1}\partial^{jk}\beta_{1}+\dfrac{a}{2}h^{ij}\partial_{i}\delta\phi\partial_{j}\delta\phi$
$\displaystyle\quad\quad\quad\quad+\dfrac{H^{2}a^{3}\epsilon^{3/2}\left(24-8\epsilon+4\eta+\eta^{2}\right)}{16\sqrt{2}M_{P}}\delta\phi^{3}-\dfrac{a^{3}\sqrt{\epsilon}V_{2}}{2\sqrt{2}M_{P}}\delta\phi^{3}-\dfrac{a^{3}}{6}V_{3}\delta\phi^{3}-\dfrac{Ha^{3}\epsilon^{3/2}(\eta-2)}{4\sqrt{2}M_{P}}\delta\phi^{2}\dot{\delta\phi}$
$\displaystyle\quad\quad\quad\quad+\dfrac{a^{3}\sqrt{\epsilon}(\epsilon-2)}{4\sqrt{2}M_{P}}\delta\phi\dot{\delta\phi}^{2}-a\dot{\delta\phi}\partial_{i}\delta\phi\partial^{i}\beta_{1}-\dfrac{M_{P}\sqrt{\epsilon}}{2\sqrt{2}a}\delta\phi\partial_{ij}\beta_{1}\partial^{ij}\beta_{1}-\dfrac{a\sqrt{\epsilon}}{2\sqrt{2}M_{P}}\delta\phi\partial_{i}\delta\phi\partial^{i}\delta\phi\Bigg{\\}}\,.$
(69)

Considering that $\beta_{1}\sim\sqrt{\epsilon}\,\delta\phi$ and taking the
limit $\epsilon\ll 1$, we obtain

$\displaystyle S_{3}$
$\displaystyle\simeq\int\differential^{4}x\,\bigg{\\{}\dfrac{M_{P}^{2}a}{24}h^{ij}\left(4h^{kl}\partial_{lj}h_{ik}-3h^{kl}\partial_{kl}h_{ij}+2\partial_{k}h_{jl}\partial^{l}\tensor{h}{{}_{i}^{k}}\right)-\dfrac{M_{P}a^{3}\sqrt{\epsilon}}{8\sqrt{2}}\left(\dot{h}^{ij}\dot{h}_{ij}+\dfrac{1}{a^{2}}\partial^{k}h^{ij}\partial_{k}h_{ij}\right)\delta\phi$
$\displaystyle\quad\quad\quad\quad-\dfrac{1}{4}M_{P}^{2}a\dot{h}^{jk}\partial_{i}h_{jk}\partial^{i}\beta_{1}+\dfrac{a}{2}h^{ij}\partial_{i}\delta\phi\partial_{j}\delta\phi-\dfrac{a^{3}}{6}V_{3}\delta\phi^{3}\bigg{\\}}\,.$
(70)

Again, since $V_{3}\sim\epsilon^{-1/2}$, the term $\sim V_{3}\delta\phi^{3}$
dominates over any other self-interaction of the scalar field.

We now turn to study the quartic action. Due to the large number of terms that
appear at this order, we will keep only those that are of our interest. For
the calculation of the scalar power spectrum, we need the terms that scale as
$\sim\delta\phi^{4}$. We will also keep the terms $\sim h^{2}\delta\phi^{2}$,
which contribute to the tensor power spectrum induced by scalar modes. We
obtain:

$\displaystyle S_{4}^{\delta\phi^{4}}$
$\displaystyle+S_{4}^{h^{2}\delta\phi^{2}}\simeq\int\differential^{4}x\,\Bigg{\\{}-\dfrac{a}{4}\tensor{h}{{}_{i}^{k}}\tensor{h}{{}_{j}{}_{k}}\partial^{i}\delta\phi\partial^{j}\delta\phi-\dfrac{a^{3}V_{4}}{24}\delta\phi^{4}$
$\displaystyle-a^{3}\alpha_{2}\left(\dfrac{M_{P}^{2}}{8}\left(\dot{h}^{ij}\dot{h}_{ij}+\dfrac{1}{a^{2}}\partial^{k}h^{ij}\partial_{k}h_{ij}\right)+\dfrac{1}{2}\left(\dot{\delta\phi}^{2}+\dfrac{1}{a^{2}}\partial^{i}\delta\phi\partial_{i}\delta\phi+V_{2}\delta\phi^{2}\right)+3H^{2}M_{P}^{2}\alpha_{2}+\dfrac{2HM_{P}^{2}}{a^{2}}\partial^{2}\beta_{2}\right)$
$\displaystyle-a\left(\dfrac{M_{P}^{2}}{4}\dot{h}^{jk}\partial^{i}h_{jk}+\dot{\delta\phi}\partial^{i}\delta\phi\right)\left(\partial_{i}\beta_{2}+\beta_{2,i}\right)+\dfrac{M_{P}^{2}a}{2}\tensor{h}{{}^{k}{}^{[i}}\tensor{\dot{h}}{{}^{j]}_{k}}\partial_{i}\beta_{2,j}-\dfrac{M_{P}^{2}}{4a}\beta_{2,i}\partial^{2}\beta_{2}^{i}\Bigg{\\}}\,.$
(71)

The second and third lines are simplified using the equations of motion of
$\alpha_{2}$, $\beta_{2}$ and $\beta_{2,i}$ at lower order in $\epsilon$,

$\displaystyle 0$
$\displaystyle=6H^{2}M_{P}^{2}\alpha_{2}+\dfrac{2HM_{P}^{2}}{a^{2}}\partial^{2}\beta_{2}+\dfrac{M_{P}^{2}}{8}\left(\dot{h}^{ij}\dot{h}_{ij}+\dfrac{1}{a^{2}}\partial^{k}h^{ij}\partial_{k}h_{ij}\right)+\dfrac{1}{2}\left(\dot{\delta\phi}^{2}+\dfrac{1}{a^{2}}\partial^{i}\delta\phi\partial_{i}\delta\phi+V_{2}\delta\phi^{2}\right)\,,$
(72) $\displaystyle 0$
$\displaystyle=\dfrac{M_{P}^{2}}{2a^{2}}\partial^{2}\beta_{2,j}-2HM_{P}^{2}\partial_{j}\alpha_{2}+\dfrac{M_{P}^{2}}{4}\left(\dot{h}^{ik}\left(\partial_{j}h_{ik}-\partial_{k}h_{ij}\right)+h^{ik}\partial_{k}\dot{h}_{ij}\right)+\dot{\delta\phi}\partial_{j}\delta\phi\,,$
(73)

obtaining

$\displaystyle S_{4}^{\delta\phi^{4}}+S_{4}^{h^{2}\delta\phi^{2}}$
$\displaystyle\simeq\int\differential^{4}x\,\Bigg{\\{}-\dfrac{a}{4}\tensor{h}{{}_{i}^{k}}\tensor{h}{{}_{j}{}_{k}}\partial^{i}\delta\phi\partial^{j}\delta\phi-\dfrac{a^{3}V_{4}}{24}\delta\phi^{4}+3a^{3}H^{2}M_{P}^{2}\alpha_{2}^{2}\Bigg{\\}}$
$\displaystyle\simeq\int\differential^{4}x\,\Bigg{\\{}-\dfrac{a}{4}\tensor{h}{{}_{i}^{k}}\tensor{h}{{}_{j}{}_{k}}\partial^{i}\delta\phi\partial^{j}\delta\phi-\dfrac{a^{3}V_{4}}{24}\delta\phi^{4}+\dfrac{3a^{3}}{8}\partial^{-2}\partial^{i}\left(\dot{h}^{jk}\partial_{i}h_{jk}\right)\partial^{-2}\partial^{l}\left(\dot{\delta\phi}\partial_{l}\delta\phi\right)\Bigg{\\}}\,.$
(74)

We conclude that the set of important self-interactions of $\delta\phi$ up to
fourth order are those coming from the potential. This conclusion can be
generalized to any perturbative order. The reason is that the interactions
coming from the potential at order $n$ depend on $V_{n}$, while the
interactions coming from the metric depend on $N_{(n-2)}$ and $N_{(n-2),i}$ at
most, and therefore on $V_{(n-2)}$. As we have seen above
$V_{(n-2)}/V_{n}\sim\epsilon$. It is important to stress that this only takes
into account that $\epsilon\ll 1$, assuming that the rest of SR parameters are
$\sim 1$ –at least in some phase of inflation that we are interested in–.
Hence, this result does not apply to standard SR but does apply to USR.

##  Appendix C Classical dynamics of super-Hubble modes

We study now the classical dynamics of the fields $\delta\phi$ and $h_{ij}$.
First, let us define a variable with respect to which all calculations will be
much simpler,

$\varphi\equiv-H\dfrac{\delta\phi}{\dot{\phi}_{0}}\,.$ (75)

We work in the super-Hubble limit, i.e. $k\ll aH$ (and not in the $\epsilon\ll
1$ limit, which we will not use again), in which we neglect gradients. We are
also interested in taking a late time limit, which is formally the limit
$a\to\infty$. Moreover, we will assume that $1/(a^{3}\epsilon)\to 0$. Under
these simplifications, and working with the variable $\varphi$, the quadratic
action simplifies to

$\displaystyle
S_{2}=\int\differential^{4}x\Bigg{\\{}\dfrac{M_{P}^{2}a^{3}}{8}\dot{h}^{ij}\dot{h}_{ij}+M_{P}^{2}a^{3}\epsilon\dot{\varphi}^{2}+\partial_{t}\left(HM_{P}^{2}a^{3}\epsilon\left(-\epsilon+\dfrac{\eta}{2}\right)\varphi^{2}\right)\Bigg{\\}}\,.$
(76)

The equations of motion are

$\partial_{t}\left(a^{3}\dot{h}_{ij}\right)=\order{2}\,,\quad\partial_{t}\left(a^{3}\epsilon\dot{\varphi}\right)=\order{2},$
(77)

where we write $\order{2}$ because we have not included the effects of the
cubic action. That is, these equations of motion have corrections that are
quadratic in the fields. Neglecting them, for $k\ll aH$ and at late times,
both $\varphi$ and $h_{ij}$ are constant, i.e. $\dot{\varphi}=\order{2}$ and
$\dot{h}_{ij}=\order{2}$.

We will now consider what those quadratic corrections are, for which we have
to study the cubic action under the same approximations. When neglecting
gradients, special care must be taken with terms of the form
$\sim\partial^{2}\beta_{1}$, since $\beta_{1}$ is defined with the inverse of
the Laplacian, see eq. (64). Terms with two time derivatives contribute to the
equations of motion as a third order correction, so we can neglect
them.161616For instance, for $\mathcal{L}\supset
f(t)\varphi\dot{\varphi}^{2}$, we obtain
$\sim\partial_{t}(f(t)\varphi\dot{\varphi})+\order{4}\sim\order{3}$. With
these considerations, the cubic action simplifies to

$\displaystyle
S_{3}=\int\differential^{4}x\,\frac{HM_{P}^{2}a^{3}\epsilon}{6}\bigg{\\{}$
$\displaystyle\left(2\epsilon^{3}-3\epsilon^{2}(2+3\eta)+\epsilon\,\eta\,(9+6\eta+4\epsilon_{3})-\eta\,\epsilon_{3}\,(3+\eta+\epsilon_{4}+\epsilon_{3})\right)H\,\varphi^{3}$
$\displaystyle\quad\quad\quad\quad\quad-$ $\displaystyle
3\epsilon\,(2\epsilon-3\eta)\,\varphi^{2}\,\dot{\varphi}\bigg{\\}}\,.$ (78)

Taking into account simultaneously $S_{2}$ and $S_{3}$, the updated equations
of motion are

$\partial_{t}\left(a^{3}\dot{h}_{ij}\right)=\order{3},\;\partial_{t}\left(a^{3}\epsilon\partial_{t}\left(\varphi+\dfrac{\eta}{4}\varphi^{2}\right)\right)=\order{3}\,,$
(79)

where we have used the fact that $\varphi^{2}$ enters and leaves the time
derivative as if it were a constant, since the corrections left by this
operation are of cubic order. Now, neglecting gradients and at late times we
have that

$\displaystyle\partial_{t}\left(\varphi+\eta\,\varphi^{2}/4\right)=\order{3}\quad\text{and
}\quad\dot{h}_{ij}=\order{3}\,.$ (80)

As we will see in Appendix D , the variable $\zeta$ of Maldacena is
$\zeta=\varphi+\eta/4\varphi^{2}$ at second order in fluctuations. That is,
obtaining that $\varphi+\eta\,\varphi^{2}/4$ is constant for super-Hubble
modes is consistent with the result that $\dot{\zeta}=0$ at all orders in the
limit $k\ll aH$ Maldacena (2003); Weinberg (2003).

Finally, we analyze the lapse and the shift in this limit in terms of
$\varphi$. At first order,

$\alpha_{1}=-\epsilon\varphi\,,\quad\dfrac{1}{a^{2}}\partial^{2}\beta_{1}=\epsilon\dot{\varphi}\,,\quad\beta_{1,i}=0\,.$
(81)

By using the equation of motion of $\varphi$, we obtain that
$\partial^{2}\beta_{1}/a^{2}=-H\,\epsilon\,\eta\,\epsilon_{3}\,\varphi^{2}/4=\order{2}$.
At second order, neglecting gradients and, again, using the equations of
motion of $\varphi$ and $h$ (i.e. imposing that
$\dot{\varphi}\sim\dot{h}\sim\order{2}$), we obtain that

$\alpha_{2}=\dfrac{\epsilon(2\epsilon+\eta)}{4}\varphi^{2}\,,\quad\dfrac{1}{a^{2}}\partial^{2}\beta_{2}=\dfrac{H\,\epsilon\,\eta\,\epsilon_{3}}{4}\varphi^{2}\,,\quad\beta_{2,i}=0\,.$
(82)

##  Appendix D Change of gauge

Starting from the decomposition of the metric and the scalar field of the eqs.
(44) and (45), we can choose the coordinate system so that we are in the
$\zeta$-gauge ($\delta\phi=E=E_{i}=0$) or the $\delta\phi$-gauge
($\zeta=E=E_{i}=0$). The question we are going to solve is what is the
relation between the variables in one gauge and the other at third order, see
also e.g. Maldacena (2003); Malik and Wands (2009); Christopherson and Malik
(2009). In particular, starting from the $\delta\phi$-gauge, we are going to
make a coordinate transformation that leads to the $\zeta$-gauge. At the end
of this appendix we will explain why we must go up to third order to be
consistent with the calculation we do of the power spectrum of $\zeta$.

In what follows we use a tilde to distinguish among the two gauges.
Concretely, the transformation we are looking for is such that

$\delta\phi\to\tilde{\delta\phi}=0\,,\quad\zeta=0\to\tilde{\zeta}\,,\quad
E=0\to\tilde{E}=0\,,\quad E_{i}=0\to\tilde{E}_{i}=0\,,\quad
h_{ij}\to\tilde{h}_{ij}\,.$ (83)

Both the lapse and the shift will also transform under the change of gauge,
but they will still be algebraic variables and therefore we do not need to
worry about their transformation properties. The fields in the new gauge will
be a series expansion of the fields in the old gauge, e.g.
$\tilde{\zeta}\equiv\sum_{n}\tilde{\zeta}_{n}$, where the subscript n
indicates the order of the expansion. The coordinate change that sends one
gauge to the other is

$x^{\mu}\to\tilde{x}^{\mu}=x^{\mu}-\xi^{\mu}(x)\,,\quad\textrm{where}\quad\xi^{\mu}(x)\equiv\left(\xi^{0}(x),\xi^{i}(x)=\partial^{i}\xi(x)+\mathbf{\xi}^{i}(x)\right)={\sum_{n}\xi^{\mu}_{n}}\,.$
(84)

We also need the inverse transformation $\tilde{x}^{\mu}\to
x^{\mu}=\tilde{x}^{\mu}+\tilde{\xi}^{\mu}(\tilde{x})$, where
$\tilde{\xi}^{\mu}(\tilde{x})$ is related to $\xi^{\mu}(x)$ by the expression

$\tilde{\xi}^{\mu}(\tilde{x})=\xi^{\mu}(x)=\xi^{\mu}(\tilde{x}+\tilde{\xi}(\tilde{x}))=\xi^{\mu}(\tilde{x})+\left(\xi^{\rho}(\tilde{x})+\xi^{\sigma}(\tilde{x})\partial_{\sigma}\xi^{\rho}(\tilde{x})\right)\partial_{\rho}\xi^{\mu}(\tilde{x})+\dfrac{1}{2}\xi^{\rho}(\tilde{x})\xi^{\sigma}(\tilde{x})\partial_{\rho\sigma}\xi^{\mu}(\tilde{x})+\order{4}\,.$
(85)

###  D.1 Scalar field transformation

Being the inflaton field $\phi$ a scalar, it transforms as
$\phi(x)\to\tilde{\phi}(\tilde{x})=\phi(x)$, which implies that
$\tilde{\phi}_{0}(\tilde{t})=\phi_{0}(t)+\delta\phi(x)$. From here, we obtain

$\displaystyle\delta\phi(x)=\tilde{\phi}_{0}(t)-\phi_{0}(t)-\xi^{0}(x)\dot{\tilde{\phi}}_{0}(t)+\dfrac{1}{2}\left(\xi^{0}(x)\right)^{2}\ddot{\tilde{\phi}}_{0}(t)-\dfrac{1}{3!}\left(\xi^{0}(x)\right)^{3}\dddot{\tilde{\phi}}_{0}(t)\,.$
(86)

At zeroth order, $\phi_{0}(t)=\tilde{\phi}_{0}(t)$. To study higher orders, we
have to expand $\xi^{0}(x)={\sum_{n}\xi^{0}_{n}(x)}$ in the expression above:

$\delta\phi=-\xi^{0}_{1}\dot{\phi}_{0}+\left(-\xi^{0}_{2}\dot{\phi}_{0}+\dfrac{1}{2}\left(\xi^{0}_{1}\right)^{2}\ddot{\phi}_{0}\right)+\left(-\xi^{0}_{3}\dot{\phi}_{0}+\xi^{0}_{1}\xi^{0}_{2}\ddot{\phi}_{0}-\dfrac{1}{3!}\left(\xi^{0}_{1}\right)^{3}\dddot{\phi}_{0}\right)$
(87)

so that, up to order three we obtain

$\displaystyle\xi^{0}_{1}(x)=-\dfrac{\delta\phi(x)}{\dot{\phi}_{0}(t)}=\dfrac{\varphi(x)}{H}\,,\quad\xi^{0}_{2}(x)=\dfrac{1}{2}\left(\xi^{0}_{1}(x)\right)^{2}\dfrac{\ddot{\phi}_{0}}{\dot{\phi}_{0}}\,,\quad\xi^{0}_{3}(x)=\dfrac{\left(\xi^{0}_{1}(x)\right)^{3}}{2}\left(\left(\dfrac{\ddot{\phi}_{0}}{\dot{\phi}_{0}}\right)^{2}-\dfrac{1}{3}\dfrac{\dddot{\phi}_{0}}{\dot{\phi}_{0}}\right)\,,$
(88)

where the variable $\varphi$ was defined in (75) and, in terms of slow-roll
parameters:

$\displaystyle\dot{\phi}_{0}=HM_{P}\sqrt{2\epsilon}\,,\quad\dfrac{\ddot{\phi}_{0}}{\dot{\phi}_{0}}=H\left(-\epsilon+\dfrac{\eta}{2}\right)\,,\quad\dfrac{\dddot{\phi}_{0}}{\dot{\phi}_{0}}=H^{2}\left(2\epsilon^{2}+\dfrac{\eta^{2}}{4}-\dfrac{5}{2}\epsilon\eta+\dfrac{\dot{\eta}}{2H}\right)\,.$
(89)

We note that the scalar field transformation only constrains $\xi^{0}$ but
leaves $\xi^{i}$ completely free, which allows us to impose the rest of the
conditions of the gauge transformation.

###  D.2 Metric transformation

The metric transforms as

$\displaystyle g_{\mu\nu}(x)\to\tilde{g}_{\mu\nu}(\tilde{x})=\dfrac{\partial
x^{\rho}}{\partial\tilde{x}^{\mu}}\dfrac{\partial
x^{\sigma}}{\partial\tilde{x}^{\nu}}g_{\rho\sigma}(x)=\left(\delta_{\mu}^{\rho}+\partial_{\mu}\tilde{\xi}^{\rho}(\tilde{x})\right)\left(\delta_{\nu}^{\sigma}+\partial_{\nu}\tilde{\xi}^{\sigma}(\tilde{x})\right)g_{\rho\sigma}(\tilde{x}+\tilde{\xi}(\tilde{x}))\,,$
(90)

where
$\partial_{\mu}\tilde{\xi}^{\rho}(\tilde{x})=\partial\tilde{\xi}^{\rho}(\tilde{x})/\partial\tilde{x}^{\mu}$
and
$g_{\rho\sigma}(\tilde{x}+\tilde{\xi}(\tilde{x}))=\sum_{n}\frac{1}{n!}\tilde{\xi}^{\delta_{1}}(\tilde{x})\dots\tilde{\xi}^{\delta_{n}}(\tilde{x})\partial_{\delta_{1}\dots\delta_{n}}g_{\rho\sigma}(\tilde{x})$
. Since we have written all the fields at the same point $\tilde{x}$, for
simplicity we will no longer write the spatial dependence. By expanding eq.
(90), we will be able to impose the rest of the conditions of our gauge
transformation, obtaining $\tilde{\xi}^{i}$, $\tilde{\zeta}$ .and
$\tilde{h}_{ij}$.

Let us consider an object $M_{ij}=M_{ji}$ with two spatial indices, which we
can decompose as follows:

$\displaystyle
M_{ij}=\delta_{ij}A+2\partial_{ij}B+2\partial_{(i}C_{j)}+D_{ij}\,,$ (91)

where $\partial^{i}C_{i}=0$ and $\delta^{ij}D_{ij}=\partial^{i}D_{ij}=0$. To
extract each of these variables from $M_{ij}$, we define ad-hoc projectors

$\displaystyle\hat{P}^{kl}\,M_{kl}$
$\displaystyle\equiv\dfrac{1}{2}\left(\delta^{kl}-\partial^{-2}\partial^{kl}\right)M_{kl}=A\,,$
(92) $\displaystyle\hat{Q}^{kl}\,M_{kl}$
$\displaystyle\equiv\dfrac{1}{4}\partial^{-2}\left(3\partial^{-2}\partial^{kl}-\delta^{kl}\right)M_{kl}=B\,,$
(93) $\displaystyle\hat{W}_{i}^{kl}\,M_{kl}$
$\displaystyle\equiv\partial^{-2}\left(\tensor{\delta}{{}_{i}^{l}}\partial^{k}-\partial^{-2}\tensor{\partial}{{}_{i}^{k}{}^{l}}\right)\,M_{kl}=C_{i}\,,$
(94) $\displaystyle\hat{F}_{ij}^{kl}\,M_{kl}$
$\displaystyle\equiv\left(\tensor{\delta}{{}_{i}^{k}}\tensor{\delta}{{}_{j}^{l}}-\dfrac{1}{2}\delta_{ij}\delta^{kl}+\dfrac{1}{2}\partial^{-2}\left(\delta_{ij}\partial^{kl}-4\tensor{\delta}{{}_{\left(i\right.}^{l}}\tensor{\partial}{{}_{\left.j\right)}^{k}}+\partial_{ij}\delta^{kl}+\partial^{-2}\tensor{\partial}{{}_{i}{}_{j}^{k}{}^{l}}\right)\right)M_{kl}=D_{ij}\,.$
(95)

These projectors will allow us to analyze the gauge transformations of the
fields. Before expanding eq. (90) in fluctuations, let us write down the
metric fluctuations we will need, in both gauges. Starting from
$g_{ij}=\gamma_{ij}$, we have

$g_{0,ij}=a^{2}\delta_{ij}\,,\quad g_{1,ij}=a^{2}h_{ij}\,,\quad
g_{2,ij}=\dfrac{1}{2}a^{2}\tensor{h}{{}_{i}^{k}}\tensor{h}{{}_{k}{}_{j}}\,,\quad
g_{3,ij}=\dfrac{1}{6}a^{2}\tensor{h}{{}_{i}^{k}}\tensor{h}{{}_{k}^{l}}\tensor{h}{{}_{l}{}_{j}}\,,$
(96)

where, again, the subscript n denotes the perturbative order of $g_{n,ij}$ and
we use a comma just to separate it from the spacetime indices. The expressions
of $\tilde{g}_{ij}$ are:

$\displaystyle\tilde{g}_{0,ij}$ $\displaystyle=a^{2}\delta_{ij}\,,$ (97)
$\displaystyle\tilde{g}_{1,ij}$
$\displaystyle=a^{2}\left(2\tilde{\zeta}_{1}\delta_{ij}+\tilde{h}_{1,ij}\right)\,,$
(98) $\displaystyle\tilde{g}_{2,ij}$
$\displaystyle=a^{2}\left(2\left(\tilde{\zeta}_{2}+\tilde{\zeta}_{1}^{2}\right)\delta_{ij}+2\tilde{\zeta}_{1}\tilde{h}_{1,ij}+\tilde{h}_{2,ij}+\dfrac{1}{2}\tensor{\tilde{h}}{{}_{1,i}^{k}}\tilde{h}_{1,kj}\right),$
(99) $\displaystyle\tilde{g}_{3,ij}$
$\displaystyle=a^{2}\bigg{\\{}2\left(\tilde{\zeta}_{3}+2\tilde{\zeta}_{2}\tilde{\zeta}_{1}+\dfrac{2}{3}\tilde{\zeta}_{1}^{3}\right)\delta_{ij}+\left(2\tilde{h}_{2,ij}+\tensor{\tilde{h}}{{}_{1,i}^{k}}\tilde{h}_{1,kj}\right)\tilde{\zeta}_{1}$
$\displaystyle\quad\quad\quad+2\left(\tilde{\zeta}_{2}+\tilde{\zeta}_{1}^{2}\right)\tilde{h}_{1,ij}+\tilde{h}_{3,ij}+\tensor{\tilde{h}}{{}_{2,\left(i\right.}^{k}}\tilde{h}_{1,k\left.j\right)}+\dfrac{1}{6}\tensor{\tilde{h}}{{}_{1,i}^{k}}\tensor{\tilde{h}}{{}_{1,k}^{l}}\tilde{h}_{1,lj}\bigg{\\}}\,.$
(100)

We stress that it is here that we have made use of the gauge freedom we had
left ($\xi_{i}$) to impose that $\tilde{E}=\tilde{E}_{i}=0$. In addition, we
will also need to know $g_{00}=-N^{2}+\gamma^{ij}N_{i}N_{j}$ and
$g_{0i}=N_{i}$,

$\displaystyle g_{0,00}=-1\,,\quad g_{1,00}=-2\alpha_{1}\,,$ (101)
$\displaystyle g_{0,0i}=0\,,\quad g_{1,0i}=\partial_{i}\beta_{1}\,,\quad
g_{2,0i}=\partial_{i}\beta_{2}+\beta_{2,i}\,,$ (102)

where $\alpha$ and $\beta$ are introduced in eq. (63). Now we can expand eq.
(90), so that at zeroth order $\tilde{g}_{0,\mu\nu}=g_{0,\mu\nu}$. At first
order,

$\displaystyle\tilde{g}_{1,\mu\nu}=g_{1,\mu\nu}+\xi_{1}^{\delta}\partial_{\delta}g_{0,\mu\nu}+2\partial_{\left(\mu\right.}\xi^{\rho}_{1}\,g_{0,\left.\nu\right)\rho}\,.$
(103)

Since we are only interested in obtaining $\tilde{\zeta}$ and $\tilde{h}$, we
can focus on the component $ij$,

$\displaystyle\tilde{g}_{1,ij}=g_{1,ij}+2a^{2}\varphi\delta_{ij}+2a^{2}\partial_{(i}\xi_{1,j)}\,.$
(104)

To obtain $\xi^{i}$, $\tilde{\zeta}$ and $\tilde{h}_{ij}$, we use the
projectors introduced in eqs. (92)–(95),

$\displaystyle\tilde{\zeta}_{1}=\varphi\,,\quad\tilde{h}_{1,ij}=h_{ij}\,,\quad\xi^{i}_{1}=0\,.$
(105)

At second order,

$\displaystyle\tilde{g}_{2,ij}=g_{2,ij}+2a^{2}\left(\left(1+\dfrac{\eta}{4}\right)\varphi^{2}+\dfrac{\varphi\dot{\varphi}}{H}\right)\delta_{ij}+2a^{2}\partial_{(i}\xi_{2,j)}+2a^{2}\varphi
h_{ij}-\chi_{2,ij}\,,$ (106)
$\displaystyle\chi_{2,ij}\equiv\partial_{i}\xi^{0}_{1}\partial_{j}\xi^{0}_{1}-2\partial_{\left(i\right.}\beta_{1}\partial_{\left.j\right)}\xi^{0}_{1}-\dfrac{a^{2}}{H}\varphi\dot{h}_{ij}\,.$
(107)

Projecting, we obtain

$\displaystyle\tilde{\zeta}_{2}=\dfrac{\eta}{4}\varphi^{2}+\dfrac{\varphi\dot{\varphi}}{H}-\dfrac{1}{2a^{2}}\hat{P}^{kl}\,\chi_{2,kl}\,,\quad\tilde{h}_{2,ij}=-\dfrac{1}{a^{2}}\hat{F}^{kl}_{ij}\,\chi_{2,kl}\,,\quad\xi_{2}=\dfrac{1}{a^{2}}\hat{Q}^{kl}\chi_{2,kl}\quad\textrm{and}\quad\mathbf{\xi}_{2,i}=\dfrac{1}{a^{2}}\hat{W}^{kl}_{i}\,\chi_{2,kl}\,.$
(108)

At third order, the number of terms involved in the gauge transformation
becomes very large. Therefore, we use the super-Hubble limit, where
$\dot{\varphi}=\dot{h}=\order{2}$ and gradients are neglected. Then, the
transformation at third order is

$\displaystyle\tilde{g}_{3,ij}=g_{3,ij}+2a^{2}\left(\dfrac{2}{3}+\dfrac{\eta}{2}+\dfrac{\eta^{2}}{12}+\dfrac{\eta\epsilon_{3}}{6}\right)\varphi^{3}\delta_{ij}+2a^{2}\left(1+\dfrac{\eta}{4}\right)\varphi^{2}h_{ij}+a^{2}\varphi\tensor{h}{{}_{i}^{k}}\tensor{h}{{}_{k}{}_{j}}+2a^{2}\partial_{(i}\xi_{3,j)}-\chi_{3,ij}\,,$
(109)

where $\chi_{3,ij}=-2\partial_{(i}\beta_{2}\partial_{j)}\xi^{0}_{1}$, which is
the only non-zero term for $k\ll aH$, because $\beta_{2}$ depends on the
inverse of the Laplacian (eq. (82)). Again, we project to obtain that

$\displaystyle\tilde{\zeta}_{3}=\dfrac{\eta^{2}+2\eta\epsilon_{3}}{12}\varphi^{3}-\dfrac{1}{2a^{2}}\hat{P}^{kl}\,\chi_{3,kl}\,,\quad\tilde{h}_{3,ij}=-\dfrac{1}{a^{2}}\hat{F}^{kl}_{ij}\,\chi_{3,kl}\,,\quad\xi_{3}=\dfrac{1}{a^{2}}\hat{Q}^{kl}\,\chi_{3,kl}\quad\textrm{and}\quad\mathbf{\xi}_{3,i}=\dfrac{1}{a^{2}}\hat{W}^{kl}_{i}\,\chi_{3,kl}\,.$
(110)

Although it seems that $\chi_{2,ij}$ and $\chi_{3,ij}$ have an effect on the
gauge transformation, their contribution vanishes for super-Hubble modes, as
we will demonstrate. Let us consider $\tilde{\zeta}$,

$\displaystyle\tilde{\zeta}=\tilde{\zeta}_{1}+\tilde{\zeta}_{2}+\tilde{\zeta}_{3}=\varphi+\dfrac{\eta}{4}\varphi^{2}+\dfrac{\varphi\dot{\varphi}}{H}+\dfrac{\eta^{2}+2\eta\epsilon_{3}}{12}\varphi^{3}-\dfrac{1}{2a^{2}}\hat{P}^{kl}\left(\chi_{2,kl}+\chi_{3,kl}\right)\,.$
(111)

Because we are working at third order, and because
$\varphi\dot{\varphi}=\order{3}$, we cannot neglect this term. In the above
expression, we get the combination

$\displaystyle\chi_{2,kl}+\chi_{3,kl}=-2\partial_{(i}\xi^{0}_{1}\partial_{j)}\left(\beta_{1}+\beta_{2}\right)-\dfrac{a^{2}}{H}\varphi\dot{h}_{ij}\,.$
(112)

Since $\dot{h}_{ij}=\order{3}$, we can neglect the term
$\sim\varphi\,\dot{h}_{ij}$. On the other hand,

$\displaystyle\beta_{1}+\beta_{2}=a^{2}\epsilon\partial^{-2}\left(\dot{\varphi}+\dfrac{H\eta\epsilon_{3}}{4}\varphi^{2}\right)=a^{2}\epsilon\partial^{-2}\partial_{t}\left(\varphi+\dfrac{\eta}{4}\varphi^{2}\right)+\order{3}=\order{3}\,,$
(113)

where we have used eq. (80). We conclude that $\chi_{2}+\chi_{3}=\order{4}$
and therefore we can neglect its effect at third order,

$\displaystyle\tilde{\zeta}=\varphi+\dfrac{\eta}{4}\varphi^{2}+\dfrac{\eta^{2}-\eta\,\epsilon_{3}}{12}\varphi^{3}\,.$
(114)

The same applies to $\tilde{h}_{ij}$ and $\xi_{i}$, so that finally

$\displaystyle\tilde{h}_{ij}$
$\displaystyle=\tilde{h}_{1.ij}+\tilde{h}_{2,ij}+\tilde{h}_{3,ij}=h_{ij},$
(115) $\displaystyle\xi_{i}$
$\displaystyle=\xi_{1,i}+\xi_{2,i}+\xi_{3,i}=0\,.$ (116)

An important conclusion is that at third order $h_{ij}$ is invariant under a
change of gauge from the $\zeta$-gauge to the $\delta\phi$-gauge. Although we
do not prove it, this result suggests it can be generalized at higher orders
in perturbation theory.

Let us assume that eq. (116) can be extended to all orders in perturbation
theory. We stress that this is not a given. As we have just seen eq. (116)
holds thanks to a subtle cancellation involving the inverse of the Laplacian,
see eqs. (112) and (113). Therefore in the transformation of the metric (90)
we are going to neglect $\partial_{i}\tilde{\xi}^{\rho}$:

$\tilde{g}_{ij}(\tilde{x})=g_{ij}(\tilde{x}+\tilde{\xi}(\tilde{x}))=\sum_{n}\dfrac{1}{n!}\tilde{\xi}^{\delta_{1}}(\tilde{x})\dots\tilde{\xi}^{\delta_{n}}(\tilde{x})\partial_{\delta_{1}\dots\delta_{n}}g_{ij}(\tilde{x})\simeq\sum_{n}\dfrac{1}{n!}\left(\tilde{\xi}^{0}(\tilde{x})\right)^{n}\dfrac{\partial^{n}}{\partial
t^{n}}g_{ij}(\tilde{x})=g_{ij}(\tilde{t}+\tilde{\xi}^{0}(\tilde{x}),\tilde{\mathbf{x}})\,,$
(117)

where we have also neglected the gradients of the metric,
$\partial_{k}g_{ij}$. Then, from
$\tilde{\gamma}_{ij}(\tilde{x})=\gamma_{ij}(\tilde{t}+\tilde{\xi}^{0}(\tilde{x}),\tilde{\mathbf{x}})$
we get

$2\left(\log
a(\tilde{t})+\tilde{\zeta}(\tilde{x})\right)\delta_{ij}+\tilde{h}_{ij}(\tilde{x})=2\log
a(\tilde{t}+\tilde{\xi}^{0}(\tilde{x}))\delta_{ij}+h_{ij}(\tilde{t}+\tilde{\xi}^{0}(\tilde{x}),\tilde{\mathbf{x}})\,.$
(118)

Taking into account that $h_{ij}$ is constant for $k\ll H$
($h_{ij}(\tilde{t}+\tilde{\xi}^{0}(\tilde{x}),\tilde{\mathbf{x}})\simeq
h_{ij}(\tilde{t},\tilde{\mathbf{x}})$), we have finally that171717See e.g.
Sugiyama  _et al._ (2013) for the connection with the $\delta N$-formalism for
large scale modes of $\zeta$.

$\tilde{\zeta}(\tilde{x})=\log\dfrac{a(t)}{a(\tilde{t})}\,,\quad\tilde{h}_{ij}(\tilde{x})=h_{ij}(\tilde{x})\,.$
(119)

We stress that the decomposition (44) has been instrumental in the derivation
of this result. A different decomposition would mix $h_{ij}$ with lower
helicity variables.

###  D.3 Two-point correlation function of $\zeta$

At late times, the relation between the variables $\zeta$ and $\delta\phi$ in
their respective gauges is

$\zeta=-\dfrac{\delta\phi}{\sqrt{2\epsilon}M_{P}}+\dfrac{\eta}{4}\left(\dfrac{\delta\phi}{\sqrt{2\epsilon}M_{P}}\right)^{2}-\dfrac{\eta^{2}-\eta\epsilon_{3}}{12}\left(\dfrac{\delta\phi}{\sqrt{2\epsilon}M_{P}}\right)^{3}+\order{\dfrac{\delta\phi}{M_{P}}}^{4}\,.$
(120)

This relation is exact in the sense that it does not make any SR
approximation. We compute $\mathcal{P}_{\zeta}$ at one-loop using correlations
of $\delta\phi$. Knowing the expression above, we can explain why it is
necessary to make the expansion up to third order. The two-point correlation
of $\zeta$ is, schematically,

Figure 6:  Complete set of diagrams that affect the calculation at one-loop of
$\expectationvalue{\zeta^{2}}$ in terms of correlations of $\delta\phi$.
Bubble diagrams and terms proportional to tadpoles have been omitted.

$\displaystyle\expectationvalue{\zeta^{2}}\sim\dfrac{1}{M_{P}^{2}}\expectationvalue{\delta\phi^{2}}+\dfrac{1}{M_{P}^{3}}\expectationvalue{\delta\phi^{3}}+\dfrac{1}{M_{P}^{4}}\expectationvalue{\delta\phi^{4}}+\order{\dfrac{H}{M_{P}}}^{5}\,,$
(121)

where we introduce $H$ as it is the only energy scale of inflation. Assuming
the set of interactions we work with in this paper, the contributions to the
two-point correlation of $\zeta$ at order $(H/M_{P})^{4}$ are shown in Figure
6. Since $V_{3}\sim H(H/M_{P})$ and $V_{4}\sim(H/M_{P})^{2}$, the one-loop
diagrams that contribute to $\langle\delta\phi^{2}\rangle$ are of the same
order in $(H/M_{P})$ as the contribution to $\langle\zeta^{2}\rangle$ coming
from $\langle\delta\phi^{3}\rangle$ and $\langle\delta\phi^{4}\rangle$.

In the model we consider, $\eta=0$ at late times and therefore only the first
term of the expansion (120) survives. However, we stress that (120) is valid
for any inflationary model as long as at late times $\dot{\varphi}\propto
1/(a^{3}\epsilon)\to 0$.

##  Appendix E Tadpoles

In this appendix we discuss why we do not need to include tadpoles in our
calculations. Since we work with a theory with cubic interactions, we
inevitably get a non-zero one-loop contribution to the correlator
$\expectationvalue{\delta\phi(x)}$, as shown in Figure 7. In principle, this
type of interaction implies the need of considering contributions to the power
spectrum beyond those studied in the main body of this work. However, by using
counterterms we can avoid having to include these effects, imposing that
$\expectationvalue{\delta\phi(x)}=0$.

The interaction Hamiltonian we are going to use is

$H_{I}(\tau)=\int\differential^{3}\mathbf{x}\left[a^{4}\dfrac{V_{3}\delta\phi^{3}}{3!}+\delta_{\rm
tp}\delta\phi\right]\,.$ (122)

The coupling of the counterterm, $\delta_{\rm tp}$, will be responsible for
the renormalization of the tadpole. Since this counterterm comes from
renormalizing the potential (in particular, from making $V_{1}\to
V_{1}+\delta_{V_{1}}$), $\delta_{\rm tp}$ will have an arbitrary time
dependence, as well as $\delta_{V}\equiv\delta_{V_{2}}$. The reason why we do
not include the quartic interaction here is that it contributes to the next
order in loops. Using eq. (5), we see that if we impose

$\displaystyle\expectationvalue{\delta\phi(x)}$
$\displaystyle=2\Im{\int_{-\infty}^{\tau}\differential\tau^{\prime}\bra{0}\delta\phi(x)\int\differential^{3}\mathbf{x}^{\prime}\left[a^{4}(\tau^{\prime}_{-})\dfrac{V_{3}(\tau^{\prime}_{-})\delta\phi^{3}(x^{\prime}_{-})}{3!}+\delta_{\rm
tp}(\tau^{\prime}_{-})\delta\phi(x^{\prime}_{-})\right]\ket{0}}$
$\displaystyle=\int_{-\infty}^{\tau}\differential\tau^{\prime}\Im{\delta\phi_{0}(\tau)\delta\phi_{0}^{*}(\tau^{\prime}_{-})\left[a^{4}(\tau^{\prime}_{-})V_{3}(\tau^{\prime}_{-})\int\dfrac{\differential^{3}\mathbf{p}}{(2\pi)^{3}}\absolutevalue{\delta\phi_{p}(\tau^{\prime}_{-})}^{2}+2\delta_{\rm
tp}(\tau^{\prime}_{-})\right]}=0\,,\quad\forall\,\tau\,,$ (123)

in order for this equation to be satisfied, we need the counterterm to be

$2\delta_{\rm
tp}(\tau)=-a^{4}(\tau)V_{3}(\tau)\int\dfrac{\differential^{3}\mathbf{p}}{(2\pi)^{3}}\absolutevalue{\delta\phi_{p}(\tau)}^{2}\,.$
(124)

Considering this functional form of $\delta_{\rm tp}$, it is easy to see that
the tadpole contributions to the two-point correlation (diagrams in Figure 7)
through the eq. (5) vanish.

Figure 7:  Diagrams that contribute to the tadpole
$\expectationvalue{\delta\phi(x)}$ (left) or are proportional to it and
contribute to two-point correlation
$\expectationvalue{\delta\phi(x)\delta\phi(y)}$ (right).

##  Appendix F Asymptotic limits

In this appendix we analyze in detail some aspects of the calculation of
$\mathcal{P}_{\zeta}$ at one-loop. First, we consider the limit $k\to 0$,
concluding that all the contributions to $\mathcal{P}_{\zeta}$ are scale
invariant in this limit, but different from zero in general. For completeness,
we also discuss briefly the limit $k\to\infty$. Finally, we study the
parameter dependence of $\mathcal{P}_{\zeta}$ and use this dependence to
estimate its magnitude.

###  F.1 Limit $k\to 0$

Looking at the equation of motion for $\delta\phi_{k}$,

$\delta\phi_{k}^{\prime\prime}+2Ha\delta\phi_{k}^{\prime}+(k^{2}+a^{2}V_{2})\delta\phi_{k}=0\,,$
(125)

and the Bunch-Davies initial conditions $\delta\phi_{k}(\tau\to-\infty)\to
e^{-ik\tau}/\sqrt{2ka^{2}}$, we note that $\delta\phi_{k}$ admits a series
expansion around $k=0$

$\delta\phi_{k}(\tau)=\dfrac{1}{k^{3/2}}\sum_{n=0}^{\infty}k^{n}\delta\phi_{n}(\tau)\,.$
(126)

Each term $\delta\phi_{n}$ can be obtained by solving the equation of motion
and imposing the initial (and boundary) conditions order by order in $k$. In
particular, the solution of the modes in the first SR phase is

$\delta\phi_{k}(\tau<\tau_{1})=\dfrac{iHe^{-ik\tau}}{\sqrt{2k^{3}}}(1+ik\tau)\,{.}$
(127)

Then, it is straightforward to obtain the evolution of $\delta\phi_{n}$ along
inflation. Since $\delta\phi_{k}/\sqrt{\epsilon}$ is constant in the limit
$k\to 0$,

$\delta\phi_{0}(\tau)=\sqrt{\dfrac{\epsilon(\tau)}{\epsilon_{SR}}}\dfrac{iH}{\sqrt{2}}\,,$
(128)

where $\epsilon_{SR}\equiv\epsilon(\tau<\tau_{1})$ is constant. We can now
analyze the limit $k\to 0$ of the loop integral $l(k;i,j)$, which is IR
divergent. Starting from the general expression, eq. (15),

$l(k;i,j)=\dfrac{k^{3}}{16\pi^{2}}\int_{1+\Lambda_{\rm
IR}/k}^{\infty}\differential s\int_{0}^{1}\differential
d\>(s^{2}-d^{2})\delta\phi_{k(s-d)/2}(\tau_{i})\delta\phi^{*}_{k(s-d)/2}(\tau_{j-})\delta\phi_{k(s+d)/2}(\tau_{i})\delta\phi^{*}_{k(s+d)/2}(\tau_{j-})\,,$
(129)

we observe that the IR divergence comes from the region $s\sim 1$ and $d\sim
s$, given that

$\delta\phi_{k(s-d)/2}(\tau_{i})=\dfrac{1}{(k(s-d)/2)^{3/2}}\left(\delta\phi_{0}(\tau_{i})+\order{k^{2}(s-d)^{2}}\right)\,.$
(130)

The prescription $\tau_{j-}$ has an effect only on the UV, so we can omit it
for the calculation of the IR divergence. Then, the IR divergent part of
$l(k;i,j)$ is

$\displaystyle l^{\rm IR}(k;i,j)$
$\displaystyle=\dfrac{\delta\phi_{0}(\tau_{i})\delta\phi_{0}^{*}(\tau_{j})}{\pi^{2}}\int_{1+\Lambda_{\rm
IR}/k}\differential s\int_{0}^{1}\differential
d\>\dfrac{s}{(s-d)^{2}}\left(\delta\phi_{ks}(\tau_{i})\delta\phi_{ks}^{*}(\tau_{j})+\order{s-d}^{2}\right)$
$\displaystyle=\dfrac{\delta\phi_{0}(\tau_{i})\delta\phi_{0}^{*}(\tau_{j})}{\pi^{2}}\int_{1+\Lambda_{\rm
IR}/k}\differential
s\dfrac{1}{s-1}\left(\delta\phi_{ks}(\tau_{i})\delta\phi_{ks}^{*}(\tau_{j})+\textrm{finite}\right)$
$\displaystyle=-\dfrac{\log\Lambda_{\rm
IR}/k}{\pi^{2}}\delta\phi_{0}(\tau_{i})\delta\phi_{0}^{*}(\tau_{j})\delta\phi_{k}(\tau_{i})\delta\phi_{k}^{*}(\tau_{j})+\textrm{finite}\,.$
(131)

The IR divergence is $\sim\log k/\Lambda_{\rm IR}$, a standard result in this
type of calculations. As already mentioned, it is not the aim of our work to
deal with IR divergences, so we eliminate them by hand, by redefining
$l(k;i,j)$ subtracting from it its IR divergent part, $l^{\rm IR}(k;i,j)$.
With this redefinition there are no IR divergences (by construction), but
neither UV divergences due to the $\tau_{j-}$ prescription.

Now, we compute the limit $k\to 0$,

$\displaystyle\lim_{k\to 0}l(k;i,j)$
$\displaystyle=\dfrac{4}{\pi^{2}k^{3}}\delta\phi_{0}^{2}(\tau_{i})\delta\phi_{0}^{*2}(\tau_{j})\int_{1+\Lambda_{\rm
IR}/k}^{\infty}\differential s\int_{0}^{1}\dfrac{\differential
d}{\left(s^{2}-d^{2}\right)^{2}}+\dfrac{\log\Lambda_{\rm
IR}/k}{\pi^{2}k^{3}}\delta\phi_{0}^{2}(\tau_{i})\delta\phi_{0}^{*2}(\tau_{j})$
$\displaystyle=\dfrac{\log
2-1}{\pi^{2}k^{3}}\delta\phi_{0}^{2}(\tau_{i})\delta\phi_{0}^{*2}(\tau_{j})=\dfrac{\epsilon(\tau_{i})\epsilon(\tau_{j})H^{4}}{4\pi^{2}\epsilon_{SR}^{2}k^{3}}\left(\log
2-1\right)+\order{1/k^{2}}\,,$ (132)

which is real and its dependence on $k$ is $\Re{l(k)}\sim k^{-3}$. The first
imaginary contribution is $\Im{l(k)}\sim k^{0}$. However, calculating it is
complicated because it receives a contribution from both the IR and the UV, so
we cannot do an expansion around $k=0$ as we have done to calculate the
dominant real part. Nevertheless, the imaginary contribution plays an
important role in the IR limit ($k\to 0$) of the power spectrum at one-loop.

We are now in a position to analyze each of the contributions to the power
spectrum in this limit. Starting with the tree-level,

$\mathcal{P}^{\rm
tl}_{\zeta}(\tau,k)=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\absolutevalue{\delta\phi_{k}(\tau)}^{2}\to\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\dfrac{1}{k^{3}}\absolutevalue{\delta\phi_{0}(\tau)}^{2}=\dfrac{H^{2}}{8\pi^{2}M_{P}^{2}\epsilon_{SR}}+\order{k^{2}}\,,$
(133)

we recover the usual scale invariant spectrum. The contribution of the
counterterms is

$\displaystyle\mathcal{P}^{\rm ct}_{\zeta}(\tau,k)$
$\displaystyle=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\left(2\,\delta_{\phi}\,\absolutevalue{\delta\phi_{k}(\tau)}^{2}+8\Im{\delta\phi_{k}^{2}(\tau)\int_{-\infty}^{\tau}\differential\tau^{\prime}a^{2}(\delta_{\phi}k^{2}+\tilde{\delta}_{V})\delta\phi_{k}^{*2}\evaluated{}_{\tau^{\prime}}}\right)$
$\displaystyle\to 2\,\delta_{\phi}\,\mathcal{P}^{\rm
tl}_{\zeta}(\tau,k)+\dfrac{2k^{3}}{\pi^{2}M_{P}^{2}\epsilon(\tau)}\Im{\delta\phi_{k}^{2}(\tau)\int_{-\infty}^{\tau}\differential\tau^{\prime}a^{2}\tilde{\delta}_{V}\delta\phi_{k}^{*2}\evaluated{}_{\tau^{\prime}}}\,.$
(134)

We have to analyze the term
$\Im{\delta\phi_{k}^{2}(\tau)\delta\phi_{k}^{2*}(\tau^{\prime})}=2\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau^{\prime})}\Re{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau^{\prime})}$.
Making an expansion around $k\to 0$, we get

$\displaystyle\dfrac{1}{\sqrt{\epsilon(\tau)\epsilon(\tau^{\prime})}}\Re{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau^{\prime})}$
$\displaystyle\to\dfrac{H^{2}}{2k^{3}\epsilon_{SR}}+\order{1/k}\,,$ (135)
$\displaystyle\dfrac{1}{\sqrt{\epsilon(\tau)\epsilon(\tau^{\prime})}}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau^{\prime})}$
$\displaystyle\to\dfrac{1-\left(a(\tau^{\prime})/a(\tau)\right)^{3+\eta}}{2(3+\eta)a^{3}(\tau^{\prime})H\epsilon(\tau^{\prime})}+\order{k^{2}}\,,$
(136)

where we have made a set of simplifications, taking $\eta$ as a constant and
assuming that $\tau$ and $\tau^{\prime}$ belong in the same phase. We note
that the imaginary part is $\sim k^{0}$, because it satisfies the boundary
conditions $\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau)}=0$ and
$\Im{\delta\phi^{\prime}_{k}(\tau)\delta\phi_{k}^{*}(\tau)}=-1/(2a^{2}(\tau))$.
We also note that this imaginary part has two solutions: one constant and one
decaying (with the well-known exception occurring when $\eta<-3$). When the
simplifications we have made are no longer valid, the solution for the real
part remains the same, but imaginary part is more complicated. However, its
dependence with the momentum $k$ in the IR limit will still be
$\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau^{\prime})}\to\order{k^{0}}$.
Keeping this in mind, we observe that the leading contribution of the
counterterms to the power spectrum is also scale invariant, analogous to the
tree-level one. Finally, we analyze the one-loop contribution:

$\displaystyle\mathcal{P}^{\rm
1l}_{\zeta}(\tau,k)=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\dfrac{2H^{2}}{M_{P}^{2}}\sum_{i=1}^{4}\sum_{j=i+1}^{4}\dfrac{a^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{i}}\dfrac{a^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{j}}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{j})}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{i})\,l(k;j,i)}\,{.}$
(137)

As with the counterterms, we note that
$\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{j})}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{i})\,l(k;j,i)}\to\order{1/k^{3}}$,
indicating that the one-loop contribution is scale invariant as well in this
limit. Consequently, we deduce that all contributions to the total power
spectrum (tree-level, one-loop and counterterms) are scale invariant in the
limit $k\to 0$.

###  F.2 Limit $k\to\infty$

For completeness, we will analyze the behavior in the UV limit $k\to\infty$ of
the one-loop contribution to the power spectrum. In this limit, the modes are
sub-Hubble, they do not feel the effects of the background and then they are
still in Bunch-Davies, $\delta\phi_{k\to\infty}(\tau)\to
e^{-ik\tau}/\sqrt{2ka^{2}}$. Then, the loop integral is

$\displaystyle\lim_{k\to\infty}l(k;i,j)$
$\displaystyle=\dfrac{k}{16\pi^{2}a_{i}^{2}a_{j}^{2}}\int_{1+\Lambda_{\rm
IR}/k}^{\infty}\differential s\int_{0}^{1}\differential
d\>e^{-iks(\tau_{i}-\tau_{j-})}+\dfrac{\log\Lambda_{\rm
IR}/k}{\pi^{2}}\dfrac{\sqrt{\epsilon_{i}\epsilon_{j}}H^{2}e^{-ik(\tau_{i}-\tau_{j})}}{4\,\epsilon_{SR}\,k\,a_{i}\,a_{j}}$
$\displaystyle=\dfrac{ie^{-ik(\tau_{i}-\tau_{j})}}{16\pi^{2}a_{i}^{2}a_{j}^{2}(\tau_{j}-\tau_{i})}+\order{1/k}\,.$
(138)

Here we see the importance of the prescription $\tau_{j-}$ in the convergence
of the integral. If it were missing, we would obtain a result that oscillates
in the UV. We also note that if we had not subtracted the IR divergence, it
would not play any role since it is suppressed by a factor $\sim\order{1/k}$.
Because in the UV limit $l(k;i,j)\sim k^{0}$, the one-loop contribution to the
power spectrum in this limit is

$\displaystyle\mathcal{P}^{\rm 1l}_{\zeta}(\tau,k)$
$\displaystyle=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\dfrac{2H^{2}}{M_{P}^{2}}\sum_{i=1}^{4}\sum_{j=i+1}^{4}\dfrac{a^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{i}}\dfrac{a^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{j}}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{j})}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{i})\,l(k;j,i)}$
$\displaystyle\sim
k^{2}\delta\phi^{2}_{k}(\tau)\sim\dfrac{1}{k}\mathcal{P}^{\rm
tl}_{\zeta}(\tau,k)\,.$ (139)

That is, the contribution of the one-loop is suppressed with respect to that
of the tree-level.

##  Appendix G Dimensional analysis and a quick estimate of
$\mathcal{P}_{\zeta}$

Let us define a dimensionless comoving wave number $\kappa\equiv k/k_{0}$. In
terms of $k_{0}$, $k_{0}\tau=-\dfrac{\tau}{\tau_{0}}=-e^{N_{0}-N}$, where
$k_{0}\tau_{0}\equiv-1$. We use the freedom we have to set $N_{0}=0$.
Examining the equations of motion for $\delta\phi_{k}(\tau)$, we observe that
its solution will always be a function of $x\equiv k\tau=-\kappa e^{-N}$.
Therefore, $\delta\phi_{k}(\tau)=\mathcal{C}\tilde{\delta\phi}_{\kappa}(x)$,
where $\tilde{\delta\phi}_{\kappa}(x)$ is a dimensionless function and
$\mathcal{C}$ is a constant that may depend on $k$ and captures all the
dimensions of $\delta\phi_{k}(\tau)$. To determine $\mathcal{C}$, we analyze
the initial conditions of $\delta\phi_{k}(\tau)$,
$\delta\phi_{k}(\tau\to-\infty)\propto H\tau e^{-ix}/\sqrt{k}$, concluding
that

$\delta\phi_{k}(\tau)\equiv\dfrac{H}{k_{0}^{3/2}}\tilde{\delta\phi}_{\kappa}(x)\,.$
(140)

Now, the loop integral is $l(k;i,j)\sim k^{3}\int\differential s\differential
d\>\delta\phi^{4}\equiv\dfrac{H^{4}}{k_{0}^{3}}\tilde{l}(\kappa;i,j)$, where
again $\tilde{l}(\kappa;i,j)$ is dimensionless. We will also need
$a(\tau)=-\dfrac{1}{H\tau}=\dfrac{k_{0}\tau_{0}}{H\tau}=\dfrac{k_{0}}{H}\tilde{a}(x/\kappa)$.
With all that, we can move on to analyze the contribution to the power
spectrum of the tree-level

$\mathcal{P}^{\rm
tl}_{\zeta}(\tau,k)=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\absolutevalue{\delta\phi_{k}(\tau)}^{2}=\dfrac{H^{2}}{M_{P}^{2}}\dfrac{\kappa^{3}}{4\pi^{2}\epsilon(\tau)}\absolutevalue{\tilde{\delta\phi}_{\kappa}(x)}^{2}\,,$
(141)

and the one-loop

$\displaystyle\mathcal{P}^{\rm 1l}_{\zeta}(\tau,k)$
$\displaystyle=\dfrac{k^{3}}{4\pi^{2}M_{P}^{2}\epsilon(\tau)}\dfrac{2H^{2}}{M_{P}^{2}}\sum_{i=1}^{4}\sum_{j=i+1}^{4}\dfrac{a^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{i}}\dfrac{a^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{j}}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{j})}\Im{\delta\phi_{k}(\tau)\delta\phi_{k}^{*}(\tau_{i})\,l(k;j,i)}$
$\displaystyle=\dfrac{H^{4}}{M_{P}^{4}}\dfrac{\kappa^{3}}{2\pi^{2}\epsilon(\tau)}\sum_{i=1}^{4}\sum_{j=i+1}^{4}\dfrac{\tilde{a}^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{i}}\dfrac{\tilde{a}^{3}\Delta\nu^{2}}{\sqrt{\epsilon}}\evaluated{}_{\tau_{j}}\Im{\tilde{\delta\phi}_{\kappa}(x)\tilde{\delta\phi}^{*}_{\kappa}(x_{j})}\Im{\tilde{\delta\phi}_{\kappa}(x)\tilde{\delta\phi}^{*}_{\kappa}(x_{i})\,\tilde{l}(\kappa;j,i)}\,.$
(142)

We obtain, as it should be, that both contributions to the power spectrum are
dimensionless. The factor that controls (naively) the perturbative expansion
is $H^{2}/M_{P}^{2}$. Moreover, we find that the dependence on $k_{0}$
disappears. This was expected since $k_{0}$ is only a reference value without
physical meaning.

Finally, we will make an estimate of the value of $\mathcal{P}^{\rm
1l}_{\zeta}/\mathcal{P}^{\rm tl}_{\zeta}$ in a general scenario, i.e. for any
value of $\delta N$ and $\Delta N$ (always being $\delta N\ll\Delta N$). From
the previous analysis $\mathcal{P}^{\rm 1l}_{\zeta}/\mathcal{P}^{\rm
tl}_{\zeta}\sim H^{2}/M_{P}^{2}\sim\mathcal{P}^{\rm tl}_{\zeta}$. Moreover,
since $\mathcal{P}^{\rm 1l}_{\zeta}$ includes two factors of
$\Delta\nu^{2}\sim 1/\delta N$, we expect

$\dfrac{\mathcal{P}^{\rm 1l}_{\zeta}(k)}{\mathcal{P}^{\rm
tl}_{\zeta}(k)}\propto\dfrac{1}{\delta N^{2}}\mathcal{P}^{\rm
tl}_{\zeta}(k)\,.$ (143)

However, it is challenging to determine the $k$-dependence because it requires
an analytical estimate of the loop integral $l(k;i,j)$. Hence, to estimate the
magnitude of the ratio $\mathcal{P}^{\rm 1l}_{\zeta}/\mathcal{P}^{\rm
tl}_{\zeta}$ we resort to numerical calculations. In Figure 8 we illustrate
its variation with different values of $\delta N$ and $\mathcal{P}^{\rm
tl}_{\zeta}(k_{\rm peak})$. We see that

$\dfrac{\mathcal{P}^{\rm 1l}_{\zeta}(k)}{\mathcal{P}^{\rm
tl}_{\zeta}(k)}=\order{1}\dfrac{\mathcal{P}^{\rm tl}_{\zeta}(k_{\rm
peak})}{\delta N^{2}}$ (144)

is a good estimate for all different values of $\delta N$ and $\Delta N$ and
all the scales $k$.

Figure 8:  Each graphic shows, for nine different choices of $\delta N$ and
$\mathcal{P}^{\rm tl}_{\zeta}(k_{\rm peak})$: At the top, tree-level (red
dashed) and absolute value of renormalized one-loop (black continuous); at the
bottom, absolute value of renormalized one-loop divided by tree-level (black
continuous) and $\mathcal{P}^{\rm tl}_{\zeta}(k_{\rm peak})/\delta N^{2}$
(black dashed). The rows, from left to right, represent $\delta
N=\\{0.1,\,0.2,\,0.5\\}$, while the columns, from top to bottom, represent
$\mathcal{P}^{\rm tl}_{\zeta}(k_{\rm
peak})=\\{10^{-2},\,10^{-4},\,10^{-6}\\}$.

###  G.1 The limit $\delta N\to 0$

As we have seen, by analysing the dimensions of $\mathcal{P}^{\rm 1l}_{\zeta}$
and through a numerical estimation, we obtain that $\mathcal{P}^{\rm
1l}_{\zeta}\propto 1/\delta N^{2}$. In the limit $\delta N\to 0$,
$\mathcal{P}^{\rm 1l}_{\zeta}$ is therefore divergent, a difference with
respect to previous results (see e.g. Kristiano and Yokoyama (2022);
Franciolini  _et al._ (2023)). It is natural to ask what is the origin of this
difference.

In order to compare with the literature, we abandon in this Appendix the
cutoff regularization that we have been using so far. This means that we now
have to include the points along the symmetrical (red) diagonal of Figure 3,
which up until now we have left out due to the regulator. Doing so,
$\mathcal{P}^{\rm 1l}_{\zeta}$ becomes UV divergent, and its therefore
necessary to introduce another regulator that does not eliminate the diagonal
contributions, e.g. dimensional regularization or Pauli-Villars. We will left
the specific nature of the regulator unspecified and simply assume that we can
indeed regularize the integrals in the UV. We will work directly in the limit
$\delta N\to 0$, and for simplicity we will focus exclusively on the last
transition from USR to SR, i.e. for times
$\tau\in\left[\tau_{3}\,,\tau_{4}\right]$. As we shall see, this is enough to
understand the behavior of $\mathcal{P}^{\rm 1l}_{\zeta}$ in the limit $\delta
N\to 0$. With these considerations, we are going to calculate the two-point
correlation of $\zeta$. As a consistency check, we will first use the
$\delta\phi$-gauge and then the $\zeta$-gauge. We will find that both gauges
lead to the same result.

Starting with the calculation via the $\delta\phi$-gauge, we have

$\displaystyle\expectationvalue{\zeta(x)\zeta(y)}^{\rm 1l}_{\delta\phi}$
$\displaystyle=\dfrac{1}{2\epsilon(\tau)M_{p}^{2}}\Bigg{(}2\Im{\int_{-\infty}^{\tau}\differential\tau^{\prime}\expectationvalue{\delta\phi(x)\delta\phi(y)H_{I}(\tau^{\prime}_{-})}}$
$\displaystyle+2\Re{\int_{-\infty}^{\tau}\differential\tau^{\prime}\int_{\tau^{\prime}}^{\tau}\differential\tau^{\prime\prime}\expectationvalue{\left(H_{I}(\tau^{\prime\prime}_{+})\delta\phi(x)\delta\phi(y)-\delta\phi(x)\delta\phi(y)H_{I}(\tau^{\prime\prime}_{-})\right)H_{I}(\tau^{\prime}_{-})}}\Bigg{)}\,,$
(145)

where

$H_{I}(\tau)=\int\differential^{3}\mathbf{x}\,a^{4}\left(\dfrac{V_{3}}{3!}\delta\phi^{3}+\dfrac{V_{4}}{4!}\delta\phi^{4}\right)\,.$
(146)

Since in the limit $\delta N\to 0$ we have
$\Delta\nu^{2}(\tau_{3})=-\Delta\nu^{2}(\tau_{4})=3/\delta N$, we see that the
time dependence of $V_{4}$ relevant for the calculation is of the form
$V_{4}\sim\frac{1}{\delta N}\delta(\tau-\tau_{i})$.181818In fact, $V_{4}$ also
has contributions $\sim\delta(\tau-\tau_{i})^{\prime}$. However, because the
contribution from $V_{4}$ to $\mathcal{P}^{1l}_{\zeta}$ is integrated over all
time, we can integrate by parts to move the time derivative from the Dirac
delta to the fields. Therefore, in practice the argument by which we neglect
the quartic contribution applies. The contribution to the power spectrum of
$\delta\phi$ coming from $V_{4}$ (and computed in the $\delta\phi$-gauge)
turns out to scale as

$\mathcal{P}^{\rm 1l,V_{4}}_{\delta\phi}\sim\dfrac{1}{\delta
N}\left(f(\tau,\tau_{3})-f(\tau,\tau_{4})\right)\,,$ (147)

where $f$ is a continuous function whose explicit form we do not need to
write. In the limit $\tau_{4}\to\tau_{3}$, there is a cancellation in the
numerator which makes $\mathcal{P}^{\rm
1l,V_{4}}_{\delta\phi}\sim\order{\delta N^{0}}$ and so does not diverge as
$\delta N\to 0$. As we shall see, the contribution coming from $V_{3}$ does
instead diverge in that limit, and therefore $\mathcal{P}^{\rm
1l,V_{4}}_{\delta\phi}$ is a subdominant contribution that we can neglect.
Thus, we have

$\expectationvalue{\zeta(x)\zeta(y)}^{\rm
1l}_{\delta\phi}=\dfrac{1}{\epsilon(\tau)M_{p}^{2}}\Re{\int_{-\infty}^{\tau}\differential
x^{\prime}\dfrac{a^{4}V_{3}}{3!}\int_{\tau^{\prime}}^{\tau}\differential
x^{\prime\prime}\dfrac{a^{4}V_{3}}{3!}\expectationvalue{\left(\delta\phi^{3}(x^{\prime\prime}_{+})\delta\phi(x)\delta\phi(y)-\delta\phi(x)\delta\phi(y)\delta\phi^{3}(x^{\prime\prime}_{-})\right)\delta\phi^{3}(x^{\prime}_{-})}}\,.$
(148)

This correlation induces the topology of diagrams depicted in the third row of
Figure 2. Including all the diagrams that affect the power spectrum (the first
and the second topologies),191919As shown in Appendix E, by imposing that
$\expectationvalue{\zeta(x)}=0$ the effect of the diagram proportional to the
tadpole (second of the third row of Fig. 2) vanishes. However, we note that
here we have not yet imposed any conditions on the counterterms, and therefore
all contributions to $\expectationvalue{\zeta(x)\zeta(y)}$ must be included.
Moreover, our aim is to compare the result with that obtained in another
gauge, so the only way to do this consistently is to blindly include all
diagrams, without the need to talk about counterterms. we arrive at

$\displaystyle\expectationvalue{\zeta(x)\zeta(y)}^{\rm
1l}_{\delta\phi}=\int\dfrac{\differential^{3}\mathbf{k}}{(2\pi)^{3}}e^{i\mathbf{k}(\mathbf{x}-\mathbf{y})}$
$\displaystyle\int\dfrac{\differential^{3}\mathbf{p}}{(2\pi)^{3}}\dfrac{36M_{P}^{4}H^{2}}{\delta
N^{2}}\real\\{a_{3}^{6}\epsilon_{3}^{2}\zeta_{k}(\tau)\zeta_{k}^{*}(\tau_{3})2i\Im{\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{3})}\absolutevalue{\zeta_{p}(\tau_{3})}^{2}\absolutevalue{\zeta_{\absolutevalue{\mathbf{k}-\mathbf{p}}}(\tau_{3})}^{2}+(3\to
4)$
$\displaystyle-2a_{3}^{3}\epsilon_{3}a_{4}^{3}\epsilon_{4}\zeta_{k}(\tau)\zeta_{k}^{*}(\tau_{3})2i\Im{\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{4})}\zeta_{p}(\tau_{4})\zeta_{p}^{*}(\tau_{3-})\zeta_{\absolutevalue{\mathbf{k}-\mathbf{p}}}(\tau_{4})\zeta_{\absolutevalue{\mathbf{k}-\mathbf{p}}}^{*}(\tau_{3-})$
$\displaystyle-2a_{3}^{3}\epsilon_{3}a_{4}^{3}\epsilon_{4}i\Im{\zeta_{k}^{2}(\tau_{4})\zeta_{k}^{*2}(\tau)}\absolutevalue{\zeta_{p}(\tau_{3})}^{2}\lim_{l\to
0}\zeta_{l}(\tau_{4})\zeta_{l}^{*}(\tau_{3})\Big{\\}}\,,$ (149)

where $a_{i}=a(\tau_{i})$ and $\epsilon_{i}=\epsilon(\tau_{i})$. The first two
lines correspond to the first diagram, while the third line corresponds to the
second one. Because the fields are in the interaction picture, i.e. they are
determined by the free action, we have used the linear relation between the
classical fields:
$\delta\phi_{k}(\tau)=-\sqrt{2\epsilon(\tau)}M_{p}\zeta_{k}(\tau)$. We have
also omitted the prescription $i\omega$ wherever it plays no role. Of course,
eq. (G.1) is divergent in the UV, and we stress that we would need to
regularize it in order to obtain a finite result. We will now take the limit
$\delta N\to 0$. To do so, we note that
$a^{3}_{4}\epsilon_{4}=a^{3}_{3}\epsilon_{3}+\order{\delta N^{2}}$ and
$\tau_{4}-\tau_{3}=\delta N/(a_{3}H)$. Moreover, since $\zeta$ is smooth in
this limit, we can expand $\zeta_{k}(\tau_{4})$ in series around $\tau_{3}$.
Then, we arrive at

$\displaystyle\expectationvalue{\zeta(x)\zeta(y)}^{\rm
1l}_{\delta\phi}=\int\dfrac{\differential^{3}\mathbf{k}}{(2\pi)^{3}}e^{i\mathbf{k}(\mathbf{x}-\mathbf{y})}\int\dfrac{\differential^{3}\mathbf{p}}{(2\pi)^{3}}\dfrac{54M_{P}^{2}H}{\delta
N}a_{3}^{3}\epsilon_{3}\Im{\zeta_{k}^{2}(\tau)\zeta_{k}^{*2}(\tau_{3})}\absolutevalue{\zeta_{p}(\tau_{3})}^{2}\,.$
(150)

where we have used
$\Im{\zeta_{k}(\tau)\zeta_{k}^{*^{\prime}}(\tau)}=1/(4\epsilon(\tau)M_{P}^{2}a^{2}(\tau))$,
which is a general relation derived from the quantization condition. The two-
point correlation obtained is valid for the whole range of external comoving
momentum $k$ in the limit $\delta N\to 0$. We conclude that, using a regulator
that does not forbid the points along the (symmetric) diagonal of Figure 3,
the naive estimate of the divergence of the two-point correlation in the limit
$\delta N\to 0$ relaxes from $1/\delta N^{2}$ to $1/\delta N$. This, however,
still contradicts the results in earlier works, see e.g. Kristiano and
Yokoyama (2022); Franciolini  _et al._ (2023) (obtained in the $\zeta$-gauge)
where no such divergence was identified. We will now try to address this
apparent inconsistency by calculating the two-point correlation directly in
the $\zeta$-gauge. As we will see, this calculation agrees with eq. (150).

In the $\zeta$-gauge, the quadratic action together with the dominant term of
the cubic action in the inflationary model under consideration is Maldacena
(2003)

$S=\int\differential\tau\differential^{3}\mathbf{x}\,M_{P}^{2}a^{2}\epsilon\left(\left(\zeta^{\prime}\right)^{2}-\left(\partial\zeta\right)^{2}+\dfrac{\eta^{\prime}}{2}\zeta^{\prime}\zeta^{2}\right)\,.$
(151)

We are leaving aside boundary terms which, of course, may have a relevant
contribution to the correlation functions and, as a rule, always have to be
included Arroja and Tanaka (2011); Burrage  _et al._ (2011); Braglia and Pinol
(2024). The reason why we do not include these boundary terms is because none
of them has a coupling that presents a divergence in the limit $\delta N\to
0$. In other words, they are subdominant terms in this limit. In order to
compute the two-point correlation using the in-in formalism, we need to
calculate the interaction Hamiltonian in the interaction picture. The
canonical conjugate momentum of $\zeta$ is

$\pi_{\zeta}=\dfrac{\partial\mathcal{L}}{\partial\zeta^{\prime}}=2M_{P}^{2}a^{2}\epsilon\left(\zeta^{\prime}+\dfrac{\eta^{\prime}}{4}\zeta^{2}\right)\,,\quad\zeta^{\prime}=\dfrac{\pi_{\zeta}}{2M_{P}^{2}a^{2}\epsilon}-\dfrac{\eta^{\prime}}{4}\zeta^{2}\,,$
(152)

and therefore the Hamiltonian is

$\mathcal{H}=\pi_{\zeta}\zeta^{\prime}-\mathcal{L}=M_{P}^{2}a^{2}\epsilon\left(\dfrac{\pi_{\zeta}^{2}}{4M_{P}^{4}a^{4}\epsilon^{2}}+\left(\partial\zeta\right)^{2}-\dfrac{\eta^{\prime}}{4M_{P}^{2}a^{2}\epsilon}\pi_{\zeta}\zeta^{2}+\dfrac{\left(\eta^{\prime}\right)^{2}}{16}\zeta^{4}\right)\,.$
(153)

We identify the interaction Hamiltonian as

$\mathcal{H}_{\rm
i}=M_{P}^{2}a^{2}\epsilon\left(-\dfrac{\eta^{\prime}}{4M_{P}^{2}a^{2}\epsilon}\pi_{\zeta}\zeta^{2}+\dfrac{\left(\eta^{\prime}\right)^{2}}{16}\zeta^{4}\right)\,.$
(154)

Now, we take into account an important detail (see e.g. Chen  _et al._ (2006)
as well): in the in-in formalism we work in the interaction picture, so the
conjugate momentum in this picture only sees the free action, i.e.
$\pi_{\zeta}^{I}=2M_{P}^{2}a^{2}\epsilon\left(\zeta^{I}\right)^{\prime}$.
Thus, the interaction Hamiltonian in the interaction picture is

$H_{I}(\tau)=\int\differential^{3}\mathbf{x}\,M_{P}^{2}a^{2}\epsilon\left(-\dfrac{\eta^{\prime}}{2}\zeta^{\prime}\zeta^{2}+\dfrac{\left(\eta^{\prime}\right)^{2}}{16}\zeta^{4}\right)\,.$
(155)

We stress that $H_{I}$, which has a cubic and a quartic interaction, is
induced exclusively by the cubic term in the original action. The quartic
interaction is usually neglected in the literature and, as we shall see, it is
precisely the one that gives rise to the divergence in the limit $\delta N\to
0$. Knowing $H_{I}$, we can calculate the two-point correlation

$\displaystyle\expectationvalue{\zeta(x)\zeta(y)}^{\rm 1l}_{\zeta}$
$\displaystyle=2\Im{\int_{-\infty}^{\tau}\differential\tau^{\prime}\expectationvalue{\zeta(x)\zeta(y)H_{I}(\tau^{\prime}_{-})}}$
$\displaystyle+2\Re{\int_{-\infty}^{\tau}\differential\tau^{\prime}\int_{\tau^{\prime}}^{\tau}\differential\tau^{\prime\prime}\expectationvalue{\left(H_{I}(\tau^{\prime\prime}_{+})\zeta(x)\zeta(y)-\zeta(x)\zeta(y)H_{I}(\tau^{\prime\prime}_{-})\right)H_{I}(\tau^{\prime}_{-})}}\,.$
(156)

In the limit $\delta N\to 0$, we have that
$\eta^{\prime}(\tau)=-6a_{3}H/\delta N$ for
$\tau\in\left[\tau_{3}\,,\tau_{4}\right]$. This coupling (squared) introduces
a factor $1/\delta N^{2}$, but the time integrals (over an integrand that is
smooth in this limit) cover an integration region with area $\sim\delta
N^{2}$. Therefore, the (purely) cubic contribution does not diverge, as
obtained in previous works; see in particular Franciolini  _et al._ (2023).
However, the coupling of the (induced) quartic contribution also introduces a
factor $1/\delta N^{2}$ and, since it appears inside just one time integral
(from the first line of (G.1), which gives a factor $\sim\delta N$), this
leads to the $1/\delta N$ divergence that we also find in the
$\delta\phi$-gauge. We make the calculation explicitly to show that the
results in the $\zeta$\- and $\delta\phi$-gauges are consistent.

$\displaystyle\expectationvalue{\zeta(x)\zeta(y)}^{\rm 1l}_{\zeta}$
$\displaystyle=2\Im{\int_{-\infty}^{\tau}\differential^{4}x^{\prime}\,M_{P}^{2}a^{2}\epsilon\dfrac{\left(\eta^{\prime}\right)^{2}}{16}\expectationvalue{\zeta(x)\zeta(y)\zeta^{4}(x^{\prime}_{-})}}$
$\displaystyle=\int\dfrac{\differential^{3}\mathbf{k}}{(2\pi)^{3}}e^{i\mathbf{k}(\mathbf{x}-\mathbf{y})}\int\dfrac{\differential^{3}\mathbf{p}}{(2\pi)^{3}}\dfrac{54M_{P}^{2}H}{\delta
N}a_{3}^{3}\epsilon_{3}\Im{\zeta_{k}^{2}(\tau)\zeta_{k}^{*2}(\tau_{3})}\absolutevalue{\zeta_{p}(\tau_{3})}^{2}\,.$
(157)

The equations (150) and (G.1) coincide. It appears that the $\delta
N\rightarrow 0$ divergence may have been missed in previous works due to an
incorrect derivation of $H_{I}$.

To conclude, we discuss the discrepancy between this result and the estimate
we obtained earlier regularizing with a cutoff, which was a scaling
$\mathcal{P}_{\zeta}^{\rm 1l}\sim 1/\delta N^{2}$. As we have seen, the reason
why the divergence relaxes to $1/\delta N$ after removing the cutoff is a
cancellation involving the contributions coming from the interactions located
along the (symmetric) diagonal of Figure 3. If instead we keep the cutoff –see
eq. (6)–, this cancellation does not occur and the divergence remains as
$1/\delta N^{2}$. We stress that in any case, in the limit $\delta N\to 0$ we
obtain that $\mathcal{P}_{\zeta}^{\rm 1l}$ grows unbounded; i.e. perturbation
theory is broken in that limit. Therefore, we cannot trust the estimate of
$\mathcal{P}_{\zeta}^{\rm 1l}$ obtained through any regulator.

##  Appendix H Bispectrum and consistency condition

The cubic interaction we consider in this work (see eq. (7)) implies a tree-
level contribution to the three-point function of $\zeta$ that we can
calculate using the in-in formalism (eq. (5)):

$\displaystyle\expectationvalue{\zeta(x)\zeta(y)\zeta(z)}$
$\displaystyle=-\dfrac{1}{(\sqrt{2\epsilon(\tau)}M_{P})^{3}}\expectationvalue{\delta\phi(x)\delta\phi(y)\delta\phi(z)}$
(158)
$\displaystyle=-\dfrac{1}{(\sqrt{2\epsilon(\tau)}M_{P})^{3}}2\Im{\int_{-\infty}^{\tau}\differential^{4}x^{\prime}\dfrac{a^{4}(\tau^{\prime})V_{3}(\tau^{\prime})}{3!}\bra{0}\delta\phi(x)\delta\phi(y)\delta\phi(z)\delta\phi^{3}(x^{\prime})\ket{0}}\,,$
(159)

where, as usual,
$x^{\prime}=\left(\tau^{\prime},\,\mathbf{x}^{\prime}\right)$. The
prescription $i\omega$ plays no role due to the instantaneous interactions we
consider, and so we have omitted it in the last expression. We also stress
that the linear relation we have used between $\zeta$ and $\delta\phi$ (eq.
(120)) is only valid at late times. We obtain

$\expectationvalue{\zeta(x)\zeta(y)\zeta(z)}=\int\dfrac{\differential^{3}\mathbf{p_{1}}}{(2\pi)^{3}}\dfrac{\differential^{3}\mathbf{p_{2}}}{(2\pi)^{3}}\dfrac{\differential^{3}\mathbf{p_{3}}}{(2\pi)^{3}}e^{i(\mathbf{p_{1}}\mathbf{x}+\mathbf{p_{2}}\mathbf{y}+\mathbf{p_{3}}\mathbf{z})}(2\pi)^{3}\delta(\mathbf{p_{1}}+\mathbf{p_{2}}+\mathbf{p_{3}})B_{\zeta}(\tau;\mathbf{p_{1}},\mathbf{p_{2}},\mathbf{p_{3}})$
(160)

where, at tree-level, the bispectrum is

$B^{\rm
tl}_{\zeta}(\tau;\mathbf{p_{1}},\mathbf{p_{2}},\mathbf{p_{3}})=\sum_{i=1}^{4}\dfrac{Ha^{3}_{i}\Delta\nu^{2}_{i}}{2M_{P}^{4}\left(\epsilon(\tau)^{3}\epsilon_{i}\right)^{1/2}}\Im{\delta\phi_{p_{1}}(\tau)\delta\phi^{*}_{p_{1}}(\tau_{i})\delta\phi_{p_{2}}(\tau)\delta\phi^{*}_{p_{2}}(\tau_{i})\delta\phi_{p_{3}}(\tau)\delta\phi^{*}_{p_{3}}(\tau_{i})}\,.$
(161)

In the squeezed limit, $p_{1}\equiv k\ll p_{2}\simeq p_{3}\equiv p$,

$B^{\rm
tl}_{\zeta}(\tau;\mathbf{k},\mathbf{p},-\mathbf{p})=\sum_{i=1}^{4}\dfrac{Ha^{3}_{i}\Delta\nu^{2}_{i}}{2M_{P}^{4}\epsilon^{2}(\tau)}\Im{\left(\delta\phi_{p}(\tau)\delta\phi^{*}_{p}(\tau_{i})\right)^{2}}\absolutevalue{\delta\phi_{k}(\tau)}^{2}\,,$
(162)

where we have used that
$\delta\phi_{k}(\tau_{i})=\sqrt{\epsilon_{i}/\epsilon_{k}(\tau)}\delta\phi_{k}(\tau)\left(1+\order{k^{-2}}\right)$
in the limit $\absolutevalue{k\tau}\ll 1$. Defining

$P^{\rm tl}_{\zeta}(\tau,k)\equiv\dfrac{2\pi^{2}}{k^{3}}\mathcal{P}^{\rm
tl}_{\zeta}(\tau,k)=\dfrac{\absolutevalue{\delta\phi_{k}(\tau)}^{2}}{2\epsilon(\tau)M_{P}^{2}}\,,$
(163)

we finally arrive at the expression

$B^{\rm tl}_{\zeta}(\tau;\mathbf{k},\mathbf{p},-\mathbf{p})=2\,H\,P^{\rm
tl}_{\zeta}(\tau,p)\,P^{\rm
tl}_{\zeta}(\tau,k)\sum_{i=1}^{4}a^{3}_{i}\Delta\nu^{2}_{i}\Im{\dfrac{\left(\delta\phi_{p}(\tau)\delta\phi^{*}_{p}(\tau_{i})\right)^{2}}{\absolutevalue{\delta\phi_{p}(\tau)}^{2}}}\,.$
(164)

Maldacena’s consistency relation Maldacena (2003) imposes a constraint on the
bispectrum in this squeezed limit,202020See Bravo  _et al._ (2018) for a
generalization when $\zeta$ is not constant for super-Hubble modes.

$B_{\zeta}(\tau;\mathbf{k},\mathbf{p},-\mathbf{p})=-\dfrac{\differential\log\mathcal{P}_{\zeta}(\tau,p)}{\differential\log
p}P_{\zeta}(\tau,p)P_{\zeta}(\tau,k)\,.$ (165)

Restricting the calculation to tree-level, the RHS of eq. (165) only depends
on the quadratic action through the dynamics of the free fields. The LHS,
computed in the in-in formalism and explicitly written in eq. (164), depends
on the coupling of the cubic action. Therefore, this is a non-trivial check of
the relation between the quadratic and cubic actions, which holds in our
model, as shown in Figure 9.

Figure 9:  The ratio ${B^{\rm
tl}_{\zeta}(\tau;\mathbf{k},\mathbf{p},-\mathbf{p})}/{(P^{\rm
tl}_{\zeta}(\tau,p)P^{\rm tl}_{\zeta}(\tau,k))}$ in the squeezed limit (black
continuous) and
$-\differential\log\mathcal{P}_{\zeta}(\tau,p)/{\differential\log p}$ (red
dashed) for the three models $\\{\delta N,\Delta N\\}$ of Figure 5. Both
quantities coincide for all values of $p/p_{0}$. The consistency relation is
satisfied and the bispectrum is small for large $p/p_{0}$.

As in Appendix G, it is important to consider the limit $\delta N\to 0$, since
the model with instantaneous transitions originally considered in Kristiano
and Yokoyama (2022) is recovered. It was found in Kristiano and Yokoyama
(2023); Tada  _et al._ (2023), working on the $\zeta$-gauge, that in an
analogous model defined in terms of $\eta$ the consistency relation is also
satisfied and, moreover, the bispectrum is finite. The limit of instantaneous
transitions must be treated with care to compute the bispectrum of $\zeta$
using the $\delta\phi$-gauge, since there is a subtle cancellation that
removes the naive divergence (as $\delta N\rightarrow 0$) that might be
expected to be introduced by the transitions. Being specific, in this limit
$\Delta\nu^{2}(\tau_{1})=\Delta\nu^{2}(\tau_{4})=-\Delta\nu^{2}(\tau_{2})=-\Delta\nu^{2}(\tau_{3})=-3/\delta
N+\mathcal{O}(1)$, and eq. (161) simplifies to

$\displaystyle B_{\zeta}^{\rm
tl}(\tau;\mathbf{p}_{1},\mathbf{p}_{2},\mathbf{p}_{3})$
$\displaystyle=-\dfrac{12HM_{P}^{2}}{\delta
N}\imaginary\\{\dfrac{\delta\phi_{p_{1}}(\tau)}{\sqrt{2\epsilon(\tau)}M_{P}}\dfrac{\delta\phi_{p_{2}}(\tau)}{\sqrt{2\epsilon(\tau)}M_{P}}\dfrac{\delta\phi_{p_{3}}(\tau)}{\sqrt{2\epsilon(\tau)}M_{P}}$
(166)
$\displaystyle\times\left(a_{1}^{3}\epsilon_{1}\dfrac{\delta\phi^{*}_{p_{1}}(\tau_{1})}{\sqrt{2\epsilon_{1}}M_{P}}\dfrac{\delta\phi^{*}_{p_{2}}(\tau_{1})}{\sqrt{2\epsilon_{1}}M_{P}}\dfrac{\delta\phi^{*}_{p_{3}}(\tau_{1})}{\sqrt{2\epsilon_{1}}M_{P}}+(\tau_{1}\rightarrow\tau_{4})-(\tau_{1}\rightarrow\tau_{2})-(\tau_{1}\rightarrow\tau_{3})\right)\Bigg{\\}}\,.$

In order to facilitate a comparison with the earlier literature mentioned
above, we define now $\tau_{1}\equiv\tau_{s}$ and $\tau_{3}\equiv\tau_{e}$ so
that $\tau_{2}-\tau_{1}=\delta N/(a_{s}H)+\mathcal{O}(\delta N^{2})$ and
$\tau_{4}-\tau_{3}=\delta N/(a_{e}H)+\mathcal{O}(\delta N)^{2}$. Care must be
taken when taking the limit $\delta N\to 0$, since the time derivative of
$\delta\phi_{k}(\tau)$ is not continuous. However, we stress that
$\zeta\propto\delta\phi/\sqrt{\epsilon}$ has a smooth derivative, and hence
$\zeta_{k}(\tau_{2})=\zeta_{k}(\tau_{s})+(\tau_{2}-\tau_{1})\zeta^{\prime}_{k}(\tau_{s})$
and
$\zeta_{k}(\tau_{4})=\zeta_{k}(\tau_{e})+(\tau_{4}-\tau_{3})\zeta^{\prime}_{k}(\tau_{e})$.
Again, $\epsilon$ must be treated with care, since its derivative is not
continuous in this limit, and we obtain
$a^{3}_{2}\epsilon_{2}=a^{3}_{s}\epsilon_{s}(1+\order{\delta N^{2}})$ and
$a^{3}_{4}\epsilon_{4}=a^{3}_{e}\epsilon_{e}(1+\order{\delta N^{2}})$. In this