#  One loop to rule them all:
Perturbativity in the presence of ultra slow-roll dynamics

Gabriele Franciolini<EMAIL_ADDRESS>Dipartimento di Fisica,
Sapienza Università di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy  INFN,
Sezione di Roma, Piazzale Aldo Moro 2, 00185, Roma, Italy  Antonio Junior
Iovino<EMAIL_ADDRESS>Dipartimento di Fisica, Sapienza
Università di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy  INFN, Sezione di
Roma, Piazzale Aldo Moro 2, 00185, Roma, Italy  Marco Taoso
<EMAIL_ADDRESS>Istituto Nazionale di Fisica Nucleare, Sezione di
Torino, via P. Giuria 1, I–10125 Torino, Italy  Alfredo Urbano
<EMAIL_ADDRESS>Dipartimento di Fisica, Sapienza Università di
Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy  INFN, Sezione di Roma,
Piazzale Aldo Moro 2, 00185, Roma, Italy

###### Abstract

We discuss the issue of perturbativity in single-field inflationary models
with a phase of ultra slow-roll (USR) tailor suited to generate an order-one
abundance of primordial black holes (PBHs). More in detail, we impose the
condition that loop corrections made up of short-wavelength modes enhanced by
the USR dynamics do not alter the tree-level power spectrum of curvature
perturbations. In our analysis, the USR phase is preceded and followed by two
stages of ordinary slow-roll (SR), and we model the resulting SR/USR/SR
dynamics using both instantaneous and smooth transitions. Focusing on scales
relevant for CMB observations, we find that it is not possible, with these
arguments, to rule out the scenario of PBH formation via USR, not even in the
limit of instantaneous transition. However, we also find that loop corrections
of short modes on the power spectrum of long modes, even though not large
enough to violate perturbativity requirements, remain appreciable and, most
importantly, are not tamed in realistic realisations of smooth SR/USR/SR
transitions. This makes perturbativity a powerful theoretical tool to
constrain USR dynamics. We extend the analysis at any scale beyond those
relevant for CMB observations. We find that loop corrections of short modes
remain within the few percent if compared to the tree-level power spectrum.
However, we also find one notable exception of phenomenological relevance: we
show that the so-called dip in the power spectrum of curvature perturbation is
an artifact of the tree-level computation.

###### Contents

  1. I Introduction
  2. II Set-up of the computation using the “in-in” formalism
    1. II.1 Conventions
    2. II.2 The (minimal) dynamics of ultra slow-roll
    3. II.3 The cubic action
    4. II.4 Beyond the cubic action
  3. III One-loop computation
    1. III.1 Loop correction with a large hierarchy of scales
    2. III.2 Loop correction at any scales
  4. IV Time integration beyond the instantaneous transition and at any scales
    1. IV.1 Loop evaluation at the CMB scales
      1. IV.1.1 The instantaneous transition
      2. IV.1.2 Dynamics during USR
      3. IV.1.3 Dynamics at the SR/USR/SR transition
    2. IV.2 Loop evaluation at any scales
  5. V Discussion and outlook
  6. A Dynamics of curvature modes, some essential results

##  I Introduction

In this work, we consider the so-called standard scenario of primordial black
hole (PBH) formation Hawking (1971); Zel’dovich and Novikov (1967); Carr and
Hawking (1974). In this scenario, the formation of PBHs during the early
Universe is an exceptional phenomenon in which extremely dense regions of
radiation energy are tightly packed to the point of gravitational collapse
Ivanov (1998); Garcia-Bellido  _et al._ (1996); Ivanov  _et al._ (1994).
General relativity and the inflationary stage preceding the radiation epoch
offer a mechanism to generate such over-densities: small-scale curvature
perturbations, stretched way beyond the horizon by the inflationary expansion,
are transferred to the radiation fluid after the end of inflation at around
the time of their horizon re-entry. In order to truly trigger a gravitational
collapse of the radiation fluid, the amplitude of such small-scale curvature
perturbations needs to be greatly enhanced by some dynamics during the
inflationary stage. This statement can be made more quantitative by
introducing the dimensionless power spectrum $\mathcal{P}(k)$ which gives the
contribution to the variance of the curvature perturbation field per bin of
$\log k$, with $k$ the comoving wavenumber in Fourier space. At scales
relevant for Cosmic Microwave Background (CMB) observations (that is,
$0.005\lesssim k\,[\textrm{Mpc}^{-1}]\lesssim 0.2$) we typically have
$\mathcal{P}(k)=O(10^{-9})$; at smaller scales ($1.5\times 10^{13}\lesssim
k\,[\textrm{Mpc}^{-1}]\lesssim 1.5\times 10^{14}$) asteroid-mass PBHs may
comprise the totality of dark matter (DM) observed in the Universe but their
formation requires $\mathcal{P}(k)=O(10^{-2})$ (for recent reviews see Sasaki
_et al._ (2018); Escrivà  _et al._ (2022); Özsoy and Tasinato (2023); Carr and
Kuhnel (2020); Green and Kavanagh (2021)). Theory-side, therefore, what makes
the formation of PBHs an exceptional phenomenon is the fact that it requires a
seven-order-of-magnitude enhancement of the small-scale power spectrum with
respect to the value observed at large scales.

In the context of single-field models of inflation, the above-mentioned
enhancement can be dynamically realized by introducing a phase of ultra slow-
roll (USR) Ballesteros  _et al._ (2020); Bhaumik and Jain (2021); Inomata  _et
al._ (2017); Iacconi  _et al._ (2022); Kawai and Kim (2021); Bhaumik and Jain
(2020); Cheong  _et al._ (2021); Inomata  _et al._ (2018); Dalianis  _et al._
(2019); Kannike  _et al._ (2017); Hertzberg and Yamada (2018); Ballesteros and
Taoso (2018); Garcia-Bellido and Ruiz Morales (2017); Karam  _et al._ (2023);
Rasanen and Tomberg (2019); Zhao  _et al._ (2023); Frolovsky and Ketov (2023);
Dimopoulos (2017); Germani and Prokopec (2017) during which the inflaton
field, after the first conventional phase of slow-roll (SR) that is needed to
fit large-scale cosmological observations, almost stops the descent along its
potential (typically because of the presence of a quasi-stationary inflection
point) before starting rolling down again in a final stage of SR dynamics that
eventually ends inflation. In this work, we shall refer to this three-stage
dynamics as SR/USR/SR.

A very legitimate question is whether the USR dynamics is consistent with
perturbativity. Technically speaking, the dimensionless power spectrum of
curvature perturbation $\mathcal{P}(k)$ is typically computed within the free
theory. However, curvature perturbations, being gravitational in nature,
feature an intricate architecture of non-linear interactions. The effect of
non-linear interactions is twofold. On the one hand, they generate, in
addition to the variance, non-zero higher-order cumulants that may leave a
peculiar non-Gaussian pattern to the statistics of the curvature field. On the
other hand, the variance itself gets corrected with respect to the value
computed within the free theory. In this paper, we focus on the second effect
and, on very general grounds, we define the perturbative criterion

$\displaystyle\mathcal{P}(k)\equiv\mathcal{P}_{\text{tree}}(k)\left[1+\Delta{\cal
P}_{\text{1-loop}}(k)\right]~{}~{}~{}\Longrightarrow~{}~{}~{}\Delta{\cal
P}_{\text{1-loop}}(k)\overset{!}{<}1\,,$ (1)

meaning that the power spectrum computed within the free theory (the “tree-
level” power spectrum in the above equation) must be larger than the
corrections $\Delta\mathcal{P}$ introduced by the presence of interactions.
Such corrections can be organized in a formal series expansion, and we will
focus in particular on the first-order term, dubbed “1-loop” in the above
equation.

The gut feeling is that, unless one is led to consider $\mathcal{P}_{\rm
tree}(k)=O(1)$, perturbativity should be under control. However, ref.
Kristiano and Yokoyama (2022) made the bold claim that, in the presence of
USR, the perturbativity condition in eq. (1) could be violated at scales $k$
relevant for CMB observations; even more strikingly, ref. Kristiano and
Yokoyama (2022) argues that USR dynamics tailored to generate a sizable
abundance of asteroid- or solar-mass PBHs are ruled-out. What makes the claim
of ref. Kristiano and Yokoyama (2022) so hard to accept is that it basically
says that loop of short modes alters the correlation of long CMB modes. This
is counter-intuitive since it clashes with the intuition that physics at
different scales should be decoupled.

Given the above, it is not surprising that ref. Kristiano and Yokoyama (2022)
sparked an intense debate mostly exclusively polarized on defending or
disproving the claim that PBH formation from single-field inflation is ruled
out Riotto (2023a); Kristiano and Yokoyama (2023); Riotto (2023b); Firouzjahi
(2023); Firouzjahi and Riotto (2023); Choudhury  _et al._ (2023a, b, c). In
this paper, we bury the hatchet and critically examine the consequences of eq.
(1) in the presence of single-field inflation with USR dynamics.

Our analysis is structured as follows. In section II, we set the ground for
the one-loop computation; in particular, we define all our conventions in
section II.1, the SR/USR/SR background dynamics in section II.2 and the
interaction Hamiltonian in section II.3. In section III, we compute the one-
loop correction to the curvature power spectrum within the setup described in
section II.1; in particular, in section III.1, we focus on the case in which
there is a large hierarchy between the momenta running in the loop and the
external ones while in section III.2 we consider the case in which the
external momenta are generic. In section IV, we discuss the implications of
the perturbative bound in eq. (1). In particular, in section IV.1, we consider
the case in which the external momenta are long CMB modes. In this section, we
critically compare our result with those of ref. Kristiano and Yokoyama
(2022), and discuss a number of crucial generalization. In section IV.2, we
extend the computation to the case in which the external momenta are short
modes. Finally, we conclude in section V.

##  II Set-up of the computation using the “in-in” formalism

###  II.1 Conventions

First, we set our conventions. We set the reduced Planck mass to one; $t$ is
the cosmic time (with $\dot{}\equiv d/dt$) and $\tau$ the conformal time (with
${}^{\prime}\equiv d/d\tau$) with $dt/d\tau=a$ being $a$ the scale factor of
the flat Friedmann-Lemaître-Robertson-Walker (FLRW) metric
$ds^{2}=dt^{2}-a^{2}(t)d\vec{x}^{2}$, with $\vec{x}$ comoving coordinates. The
Hubble rate is $H\equiv\dot{a}/a$. The $e$-fold time $N$ is defined by
$dN=Hdt$ from which we also have $dN/d\tau=aH$. The Hubble-flow parameters
$\epsilon_{i}$ (for $i\geqslant 1$) are defined by the recursive relation

$\displaystyle\epsilon_{i}\equiv\frac{\dot{\epsilon}_{i-1}}{H\epsilon_{i-1}}\,,~{}~{}~{}~{}~{}\textrm{with:}~{}~{}~{}\epsilon_{0}\equiv\frac{1}{H}\,.$
(2)

As customary, we simply indicate as $\epsilon$ the first Hubble parameter,
$\epsilon\equiv\epsilon_{1}=-\dot{H}/H^{2}$. Instead of the second Hubble
parameter $\epsilon_{2}$, sometimes it is useful to introduce the Hubble
parameter $\eta$ defined by111We remark that in ref. Kristiano and Yokoyama
(2022) the symbol $\eta$ refers to the second Hubble parameter $\epsilon_{2}$.

$\displaystyle\eta\equiv-\frac{\ddot{H}}{2H\dot{H}}=\epsilon-\frac{1}{2}\frac{d\log\epsilon}{dN}\,,~{}~{}~{}~{}~{}~{}\textrm{with:}~{}~{}~{}~{}~{}\epsilon_{2}=2\epsilon-2\eta\,.$
(3)

We consider the theory described by the action

$\displaystyle\mathcal{S}=\int
d^{4}x\sqrt{-g}\bigg{[}\frac{1}{2}R(g)+\frac{1}{2}g^{\mu\nu}(\partial_{\mu}\phi)(\partial_{\nu}\phi)-V(\phi)\bigg{]}\,.$
(4)

$R(g)$ is the scalar curvature associated with the space-time whose geometry
is described by the metric $g$ with line element
$ds^{2}=g_{\mu\nu}dx^{\mu}dx^{\nu}$. The classical background evolves in the
flat FLRW universe and the background value of the scalar field is a function
of time, $\phi(t)$. We tacitly assume that the scalar potential features an
approximate stationary inflection point so as to trigger the transition
SR/USR/SR.

We focus on scalar perturbations. We consider the perturbed metric in the
following generic form

$\displaystyle ds^{2}=N^{2}dt^{2}-h_{ij}(N^{i}dt+dx^{i})(N^{j}dt+dx^{j})\,,$
(5)

and choose the gauge in which

$\displaystyle N=1+\delta
N(\vec{x},t)\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}N^{i}=\delta^{ij}\partial_{j}B(\vec{x},t)\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}h_{ij}=a^{2}(t)e^{2\zeta(\vec{x},t)}\delta_{ij}\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}\delta\phi(\vec{x},t)=0\,.$
(6)

The field $\zeta(\vec{x},t)$ is the only independent scalar degree of freedom
since $N$ and $N^{i}$ are Lagrange multipliers subject to the momentum and
Hamiltonian constraints. It is important to stress that the variable $\zeta$
as defined in eq. (6) is constant outside the horizon (more in general,
outside the horizon and after the end of possible non-adiabatic phases) and
represents the correct non-linear generalization of the Bardeen variable
Maldacena (2003).

At the quadratic order in the fluctuations, the action is

$\displaystyle\mathcal{S}_{2}=\int
d^{4}x\,\epsilon\,a^{3}\left[\dot{\zeta}^{2}-\frac{(\partial_{k}\zeta)(\partial^{k}\zeta)}{a^{2}}\right]\,.$
(7)

Comoving curvature perturbations are quantized by introducing the free
operator

$\displaystyle\hat{\zeta}(\vec{x},\tau)=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\hat{\zeta}(\vec{k},\tau)e^{i\vec{x}\cdot\vec{k}}\,,~{}~{}~{}~{}~{}\textrm{with:}~{}~{}~{}\hat{\zeta}(\vec{k},\tau)=\zeta_{k}(\tau)a_{\vec{k}}+\zeta_{k}^{*}(\tau)a^{{\dagger}}_{-\vec{k}}\,,$
(8)

and

$[a_{\vec{k}},a_{\vec{k}^{\prime}}]=[a_{\vec{k}}^{{\dagger}},a_{\vec{k}^{\prime}}^{{\dagger}}]=0\,,~{}~{}~{}~{}~{}~{}~{}[a_{\vec{k}},a_{\vec{k}^{\prime}}^{{\dagger}}]=(2\pi)^{3}\delta^{(3)}(\vec{k}-\vec{k}^{\prime})\,,~{}~{}~{}~{}~{}~{}~{}a_{\vec{k}}|0\rangle=0\,,$
(9)

where the last condition defines the vacuum of the free theory $|0\rangle$. We
define the comoving wavenumber $k\equiv|\vec{k}|$. The scale factor in the
FLRW universe corresponds to a rescaling of the spatial coordinate;
consequently, physically sensible results should be invariant under the
rescaling Senatore and Zaldarriaga (2010)

$\displaystyle a\to\lambda
a\,,~{}~{}~{}~{}~{}~{}~{}\vec{x}\to\vec{x}/\lambda\,,~{}~{}~{}~{}~{}~{}~{}\vec{k}\to\lambda\vec{k}\,,~{}~{}~{}~{}~{}~{}k\to|\lambda|k\,,~{}~{}~{}~{}~{}~{}{\textrm{with}}~{}~{}\lambda\in\mathbb{R}\,.$
(10)

Furthermore, if we consider the conformal time $\tau$ instead of the cosmic
time $t$, we also have

$\displaystyle\tau\to\tau/\lambda\,.$ (11)

Notice that, under the above rescaling, we have $a_{\vec{k}}\to
a_{\vec{k}}/|\lambda|^{3/2}$ (from the scaling property of the three-
dimensional $\delta$ function) and, consequently,
$\zeta_{k}\to\zeta_{k}/|\lambda|^{3/2}$ so that
$\hat{\zeta}(\vec{k},\tau)\to\hat{\zeta}(\vec{k},\tau)/|\lambda|^{3}$ and
$\hat{\zeta}(\vec{x},\tau)$ invariant. In the case of free fields, we have

$\displaystyle\langle
0|\hat{\zeta}(\vec{k}_{1},\tau_{1})\hat{\zeta}(\vec{k}_{2},\tau_{2})|0\rangle=(2\pi)^{3}\delta(\vec{k}_{1}+\vec{k}_{2})\zeta_{k_{1}}(\tau_{1})\zeta_{k_{2}}^{*}(\tau_{2})\,.$
(12)

In the presence of a time-derivative, we simply have

$\displaystyle\langle
0|\hat{\zeta}^{\prime}(\vec{k}_{1},\tau_{1})\hat{\zeta}(\vec{k}_{2},\tau_{2})|0\rangle=(2\pi)^{3}\delta(\vec{k}_{1}+\vec{k}_{2})\zeta^{\prime}_{k_{1}}(\tau_{1})\zeta_{k_{2}}^{*}(\tau_{2})\,.$
(13)

Note that the time dependence occurs in the mode function, not in the
raising/lowering operator. The mode function $\zeta_{k}(\tau)$ is related to
the linear-order Mukhanov-Sasaki (M-S) equation. More in detail, if we define
$\zeta_{k}(\tau)=u_{k}(\tau)/z(\tau)$ with $z(\tau)\equiv
a(\tau)\sqrt{2\epsilon(\tau)}$, the mode $u_{k}(\tau)$ verifies the equation

$\displaystyle\frac{d^{2}u_{k}}{d\tau^{2}}+\left(k^{2}-\frac{1}{z}\frac{d^{2}z}{d\tau^{2}}\right)u_{k}=0\,.$
(14)

We are interested in the computation of the quantity

$\displaystyle\lim_{\tau\to
0^{-}}\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}P(k)e^{i\vec{k}\cdot(\vec{x}_{1}-\vec{x}_{2})}\,,$
(15)

at some late time $\tau\to 0^{-}$ at which curvature perturbations become
constant at super-horizon scales. Equivalently, we write

$\displaystyle\lim_{\tau\to
0^{-}}\langle\hat{\zeta}(\vec{x},\tau)\hat{\zeta}(\vec{x},\tau)\rangle=\int\frac{dk}{k}\underbrace{\left[\frac{k^{3}}{2\pi^{2}}P(k)\right]}_{\equiv\mathcal{P}(k)}=\int\frac{dk}{k}\mathcal{P}(k)\,,$
(16)

where $\mathcal{P}(k)$ is the a-dimensional power spectrum. At the level of
the quadratic action, we find

$\displaystyle\langle
0|\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)|0\rangle$
$\displaystyle=\int\frac{d^{3}\vec{k}_{1}}{(2\pi)^{3}}\,d^{3}\vec{k}_{2}\delta(\vec{k}_{1}+\vec{k}_{2})\,\zeta_{k_{1}}(\tau)\zeta_{k_{2}}^{*}(\tau)e^{i(\vec{x}_{1}\cdot
k_{1}+\vec{x}_{2}\cdot
k_{2})}=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}|\zeta_{k}(\tau)|^{2}e^{i\vec{k}\cdot(\vec{x}_{1}-\vec{x}_{2})}\,,$
(17)

which of course gives the familiar result

$\displaystyle\mathcal{P}(k)=\lim_{\tau\to
0^{-}}\frac{k^{3}}{2\pi^{2}}|\zeta_{k}(\tau)|^{2}\,.$ (18)

The goal is to compute corrections that arise from the presence of
interactions. This means that, in eq. (15), the vacuum expectation value
should refer to the vacuum $|\Omega\rangle$ of the interacting theory and the
dynamics of the operator $\hat{\zeta}(\vec{x},\tau)$ is described by the full
action that also includes interactions.

We compute the left-hand side of eq. (15) by means of the “in-in” formalism
(see e.g. Calzetta and Hu (1987); Jordan (1986); Adshead  _et al._ (2009a)).
Correlators are given by

$\displaystyle\langle\Omega|\hat{\mathcal{O}}(\tau)|\Omega\rangle\equiv\langle\hat{\mathcal{O}}(\tau)\rangle=\langle
0|\left\\{\bar{T}\exp\left[{i\int_{-\infty(1+i\epsilon)}^{\tau}d\tau^{\prime}\hat{H}_{\textrm{int}}(\tau^{\prime})}\right]\right\\}\hat{\mathcal{O}}_{I}(\tau)\left\\{T\exp\left[{-i\int_{-\infty(1-i\epsilon)}^{\tau}d\tau^{\prime}\hat{H}_{\textrm{int}}(\tau^{\prime})}\right]\right\\}|0\rangle\,,$
(19)

where on the right-hand side all fields appearing in
$\hat{\mathcal{O}}_{I}(\tau)$ and $\hat{H}_{\textrm{int}}(\tau^{\prime})$ are
free fields in the interaction picture. We shall indicate free fields in the
interaction picture with the additional subscript I. It should be noted that
the latter are nothing but the operators of the free theory that we quantized
in eq. (8).

$T$ and $\bar{T}$ are the time and anti-time ordering operator, respectively.
As customary, the small imaginary deformation in the integration contour
guarantees that $|\Omega\rangle\to|0\rangle$ as $\tau\to-\infty$ where
$|\Omega\rangle$ is the vacuum of the interacting theory. On the left-hand
side, the operator $\hat{\mathcal{O}}(\tau)$ is the equal-time product of
operators at different space points, precisely like in eq. (15). We expand in
the interaction Hamiltonian, so we use the Dyson series

$\displaystyle
T\exp\left[-i\int_{-\infty_{-}}^{\tau}d\tau^{\prime}\hat{H}_{\textrm{int}}(\tau^{\prime})\right]=1-i\int_{-\infty_{-}}^{\tau}d\tau^{\prime}\hat{H}_{\textrm{int}}(\tau^{\prime})+i^{2}\int_{-\infty_{-}}^{\tau}d\tau^{\prime}\int_{-\infty_{-}}^{\tau^{\prime}}d\tau^{\prime\prime}\hat{H}_{\textrm{int}}(\tau^{\prime})\hat{H}_{\textrm{int}}(\tau^{\prime\prime})+\dots\,,$
(20)

where, for simplicity, we introduce the short-hand notation
$\infty_{\pm}\equiv\infty(1\pm i\epsilon)$. Each order in
$\hat{H}_{\textrm{int}}$ is an interaction vertex, and carries both a time
integral and the space integral (enclosed in the definition of
$\hat{H}_{\textrm{int}}$) which in Fourier space enforces momentum
conservation.

It is crucial to correctly identify the interaction Hamiltonian. Before
proceeding in this direction, let us clarify our notation. We expand the
action in the form

$\displaystyle\mathcal{S}=\int
d^{3}\vec{x}dt\,\underbrace{\mathcal{L}[\zeta(\vec{x},t),\dot{\zeta}(\vec{x},t),\partial_{k}\zeta(\vec{x},t)]}_{\equiv\,\mathcal{L}[\zeta(\vec{x},t)]}=\underbrace{\int
d^{3}\vec{x}dt\,\mathcal{L}_{2}(\vec{x},t)}_{\equiv\,\mathcal{S}_{2}}+\underbrace{\int
d^{3}\vec{x}dt\,\mathcal{L}_{3}[\zeta(\vec{x},t)]}_{\equiv\,\mathcal{S}_{3}}+\underbrace{\int
d^{3}\vec{x}dt\,\mathcal{L}_{4}[\zeta(\vec{x},t)]}_{\equiv\,\mathcal{S}_{4}}+\dots\,,$
(21)

with $\mathcal{S}_{2}$ defined in eq. (7). We also define (as a function of
conformal time)

$\displaystyle H_{\textrm{int}}^{(k)}(\tau)\equiv\int
d^{3}\vec{x}\,\mathcal{H}_{k}[\zeta(\vec{x},\tau)]~{}~{}~{}~{}~{}\Longrightarrow~{}~{}~{}~{}~{}\hat{H}_{\textrm{int}}^{(k)}(\tau)\equiv\int
d^{3}\vec{x}\,\mathcal{H}_{k}[\hat{\zeta}_{I}(\vec{x},\tau)]\,.$ (22)

At the cubic order, we simply have

$\displaystyle\mathcal{H}_{3}[\zeta(\vec{x},\tau)]=-\mathcal{L}_{3}[\zeta(\vec{x},\tau)]\,.$
(23)

We shall construct the relevant cubic interaction Hamiltonian in section II.3.
At the quartic order, simply writing $\mathcal{H}_{4}=-\mathcal{L}_{4}$ does
not capture the correct result if the cubic Lagrangian features interactions
that depend on the time derivative of $\zeta$ since the latter modify the
definition of the conjugate momentum.

Using, at the operator level, the notation introduced in eq. (22), we
schematically write at the first order in the Dyson series expansion

$\displaystyle\langle\hat{\zeta}$
$\displaystyle(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{1^{\textrm{st}}}=$
$\displaystyle\langle
0|\hat{\zeta}_{I}(\vec{x}_{1},\tau)\hat{\zeta}_{I}(\vec{x}_{2},\tau)\bigg{[}-i\int_{-\infty_{-}}^{\tau}d\tau^{\prime}\hat{H}^{(4)}_{\textrm{int}}(\tau^{\prime})\bigg{]}|0\rangle+\langle
0|\bigg{[}i\int_{-\infty_{+}}^{\tau}d\tau^{\prime}\hat{H}^{(4)}_{\textrm{int}}(\tau^{\prime})\bigg{]}\hat{\zeta}_{I}(\vec{x}_{1},\tau)\hat{\zeta}_{I}(\vec{x}_{2},\tau)|0\rangle\,.$
(24)

At the first order, therefore, the first non-zero quantum correction involves
the quartic Hamiltonian. At the second order in the Dyson series expansion and
considering again terms with up to eight fields in the vacuum expectation
values, we write schematically

$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}$
$\displaystyle=\langle
0|\hat{\zeta}_{I}(\vec{x}_{1},\tau)\hat{\zeta}_{I}(\vec{x}_{2},\tau)\bigg{[}-\int_{-\infty_{-}}^{\tau}d\tau^{\prime}\int_{-\infty_{-}}^{\tau^{\prime}}d\tau^{\prime\prime}\hat{H}^{(3)}_{\textrm{int}}(\tau^{\prime})\hat{H}^{(3)}_{\textrm{int}}(\tau^{\prime\prime})\bigg{]}|0\rangle$
$\displaystyle+\langle
0|\bigg{[}-\int_{-\infty_{+}}^{\tau}d\tau^{\prime}\int_{-\infty_{+}}^{\tau^{\prime}}d\tau^{\prime\prime}\hat{H}^{(3)}_{\textrm{int}}(\tau^{\prime\prime})\hat{H}^{(3)}_{\textrm{int}}(\tau^{\prime})\bigg{]}\hat{\zeta}_{I}(\vec{x}_{1},\tau)\hat{\zeta}_{I}(\vec{x}_{2},\tau)|0\rangle$
$\displaystyle+\langle
0|\bigg{[}i\int_{-\infty_{+}}^{\tau}d\tau^{\prime}\hat{H}^{(3)}_{\textrm{int}}(\tau^{\prime})\bigg{]}\hat{\zeta}_{I}(\vec{x}_{1},\tau)\hat{\zeta}_{I}(\vec{x}_{2},\tau)\bigg{[}-i\int_{-\infty_{-}}^{\tau}d\tau^{\prime\prime}\hat{H}^{(3)}_{\textrm{int}}(\tau^{\prime\prime})\bigg{]}|0\rangle\,.$
(25)

The vacuum expectation values of interacting-picture fields can be computed
using Wick’s theorem. Schematically, eqs. (24, 25) give rise to the following
connected diagrams.

(26)

From the above classification, we see that, at the same loop order, we have
three classes of connected diagrams that, in principle, should be discussed
together. Notice that, contrary to the first two, the last diagram is not of
1-Particle-Irreducible (1PI) type since it consists of a tadpole attached to a
two-point propagator.

To proceed further, we need to specify the background dynamics, that shape the
time evolution of the Hubble parameters, and the interaction Hamiltonian, that
specifies which terms in the Dyson expansion contribute at a given
perturbative order.

###  II.2 The (minimal) dynamics of ultra slow-roll

We start with a discussion of the USR dynamics.

Figure 1: Classical (top panel) and quantum (central panel) dynamics in the
context of an explicit single-field model of inflation that exhibits the
presence of a phase of USR in between the time interval $N_{\textrm{\rm
in}}<N<N_{\textrm{\rm end}}$ (cf. ref. Ballesteros  _et al._ (2020)). In this
specific realization, we have $N_{\textrm{\rm in}}=36.3$ and $N_{\textrm{\rm
end}}=38.8$. Top panel: we plot the evolution of the background quantities
$\epsilon$, $\epsilon_{2}$ and $\epsilon_{2}^{\prime}$ (cf. eqs. (2, 3))
together with the evolution of the Hubble rate (normalized with respect to the
value $H_{*}\equiv H(N_{*})$ and scaled by a factor 10 for ease of
comparison). Central panel: we plot the solutions of the M-S equation in eq.
(33) for two different curvature modes. The mode in black (blue) exits the
Hubble horizon well before (during) the USR phase. Bottom panel: we plot the
classicality parameter $C_{k}$ for the same two modes (cf. the main text for
details). In the case of the black mode ($N_{k}\ll N_{\rm in}$) the
classicality parameter quickly vanishes after horizon crossing, and remains
negligible also during the USR phase. In the case of the blue mode
($N_{q}\approx N_{\rm end}$) the classicality parameter remains sizable during
the USR phase, signalling that this mode retains its quantum nature during
USR.

In order to make our discussion more concrete, in the following we shall refer
to fig. 1, see the caption for details. In this figure, we plot both the
classical (top panel) and quantum (bottom panel) dynamics that characterize a
realistic model of single-field inflation that features a phase of USR because
of the presence of an approximate stationary inflection point in the
inflationary potential. We refer to ref. Ballesteros  _et al._ (2020) for more
details about the model; we remark that this model, without including loop
corrections to the computation of the curvature power spectrum, is compatible
with CMB constraints and gives $\approx 100\%$ of dark matter in the form of
asteroid-mass PBHs.222The model predicts the tensor-to-scalar ratio $r\simeq
0.037$ which is still (barely) compatible, at 95% confidence level, with the
latest results released by the BICEP and Keck collaboration Ade  _et al._
(2021). In fig. 1, the relation between $e$-fold time $N$ (bottom $x$-axis)
and conformal time $\tau$ (top $x$-axis) is given by the integral

$\displaystyle\tau=-\frac{1}{k_{\star}}\int_{N}^{N_{0}}\frac{H(N_{\star})}{H(N^{\prime})}e^{N_{\star}-N^{\prime}}dN^{\prime}\,,$
(27)

where $N_{0}$ indicates the end of inflation, with $\tau$ conventionally set
to $0$ at this time, and $N_{\star}$ the instant of time at which the comoving
scale $k_{\star}=0.05$ Mpc-1 crosses the Hubble horizon. In fig. 1, we set
$N_{\star}=0$, and the model gives $N_{0}-N_{\star}\simeq 52$. We can
highlight few crucial properties of the dynamics presented above:

  * a)

We start from the classical analysis. During USR, $\epsilon_{2}(\tau)$ changes
according to the schematic

$\displaystyle\epsilon_{2}\approx 0\ \ \overset{\textrm{SR/USR at time
}\tau_{\textrm{in}}}{\Longrightarrow}\ \ |\epsilon_{2}|>3\ \
\overset{\textrm{USR/SR at time }\tau_{\textrm{end}}}{\Longrightarrow}\ \
\epsilon_{2}\approx O(1)\,,$ (28)

thus making $\epsilon_{2}^{\prime}(\tau)$ non-zero at around the two
transitions at conformal times $\tau_{\textrm{in}}$ and $\tau_{\textrm{end}}$
(equivalently, at $e$-fold times $N_{\textrm{in}}$ and $N_{\textrm{end}}$).
The evolution of $\epsilon_{2}$ and $\epsilon^{\prime}_{2}$ are shown in the
top panel of fig. 1, with dotted green and dashed orange lines, respectively.

  * b)

During USR, the Hubble parameter $\epsilon$ decreases exponentially fast (the
inflaton almost stops its classical motion). The evolution of $\epsilon$ is
shown in the top panel of fig. 1 (dot-dashed magenta line); in addition, we
also plot the time evolution of the Hubble rate $H$.

  * c)

We now consider the USR dynamics at the quantum level. It is crucial to
understand the typical behaviour of curvature modes (solid lines in the
central panel of fig. 1) and their time derivatives (dashed lines). In the
central panel of fig. 1, we plot two representative cases: the black lines
correspond to the case of a mode $\zeta_{k}$ that exits the Hubble horizon at
some time $N_{k}$ well before the USR phase (like a CMB mode) while the blue
lines correspond to a curvature mode $\zeta_{q}$ that exits the Hubble horizon
at some time $N_{q}$ during the USR phase. We notice that the derivative
$|d\zeta_{k}/dN|$ decays exponentially fast, and, soon after horizon crossing,
becomes negligible, while $|\zeta_{k}|$ settles to a constant value.
Consequently, we expect that interaction terms that involve the time
derivative of CMB modes will be strongly suppressed.

  * d)

Finally, we consider the issue of the quantum-to-classical transition. We
define the so-called classicality parameter Assassi  _et al._ (2013)
$C_{k}=|\zeta_{k}\dot{\zeta}_{k}^{*}-\zeta_{k}^{*}\dot{\zeta}_{k}|/|\zeta_{k}\dot{\zeta}_{k}|$
which goes to zero in the classical limit. In the case of conventional SR
dynamics, the classicality parameter scales according to $C_{k}\sim 2k\tau$,
and it vanishes exponentially fast right after the horizon crossing time. In
the bottom panel of fig. 1, we plot the classicality parameter for two
representative modes of our dynamics. The black modes experiences its horizon
crossing well before the USR phase ($N_{k}\ll N_{\rm in}$). Its classicality
parameter quickly vanishes and remains $\ll 1$ during the subsequent USR
phase. The blue line, on the contrary, represents the classicality parameter
for a mode that experiences its horizon crossing during the USR phase. Its
classicality parameter remains sizable during the USR phase, signalling that
this short mode can not be treated classically during USR.

With the aim of facilitating the numerical computations of the following
sections, instead of working with a numerical description of USR, we now
introduce a simple semi-analytical model Byrnes  _et al._ (2019); Taoso and
Urbano (2021); Franciolini and Urbano (2022). We define the hyperbolic tangent
parametrization

$\displaystyle\eta(N)$ $\displaystyle=\frac{1}{2}\left[-\eta_{\rm
II}+\eta_{\rm II}\tanh\left(\frac{N-N_{\rm in}}{\delta
N}\right)\right]+\frac{1}{2}\left[\eta_{\rm II}+\eta_{\rm III}+(\eta_{\rm
III}-\eta_{\rm II})\tanh\left(\frac{N-N_{\rm end}}{\delta N}\right)\right]\,,$
(29)

where the parameter $\delta N$ controls the width of the two transitions at
$N_{\rm in}$ and $N_{\rm end}$. The limit $\delta N\to 0$ reproduces the step-
function approximation. Using the definition

$\displaystyle\delta(x)=\lim_{\epsilon\to
0}\frac{1}{2\epsilon\cosh^{2}(x/\epsilon)}\,,$ (30)

we find

$\displaystyle\lim_{\delta N\to 0}\frac{d\eta}{dN}=(-\eta_{\rm II}+\eta_{\rm
III})\delta(N-N_{\rm end})+\eta_{\rm II}\delta(N-N_{\rm in})\,.$ (31)

Using $\eta\simeq-(1/2)d\log\epsilon/dN$, we find the following expression

$\displaystyle\frac{\epsilon(N)}{\epsilon_{\textrm{ref}}}=e^{-\eta_{\rm
III}(N-N_{\rm ref})}\left[\cosh\left(\frac{N-N_{\rm end}}{\delta
N}\right)\cosh\left(\frac{N-N_{\rm in}}{\delta N}\right)\right]^{-\frac{\delta
N\eta_{\rm III}}{2}}\left[\cosh\left(\frac{N_{\rm ref}-N_{\rm end}}{\delta
N}\right)\cosh\left(\frac{N_{\rm ref}-N_{\rm in}}{\delta
N}\right)\right]^{\frac{\delta N\eta_{\rm III}}{2}}$
$\displaystyle\left[\cosh\left(\frac{N-N_{\rm end}}{\delta N}\right){\rm
sech}\left(\frac{N-N_{\rm in}}{\delta N}\right)\right]^{\delta
N\left(\eta_{\rm II}-\frac{\eta_{\rm
III}}{2}\right)}\left[\cosh\left(\frac{N_{\rm ref}-N_{\rm end}}{\delta
N}\right){\rm sech}\left(\frac{N_{\rm ref}-N_{\rm in}}{\delta
N}\right)\right]^{\frac{\delta N}{2}(-2\eta_{\rm II}+\eta_{\rm III})},\,$ (32)

where $\epsilon_{\textrm{ref}}\ll 1$ is the value of $\epsilon$ at some
initial reference time $N_{\textrm{ref}}$. For future reference, we define
$\bar{\epsilon}(N)\equiv\epsilon(N)/\epsilon_{\textrm{ref}}$.

Figure 2: Schematic evolution of $\eta(N)$ in eq. (29) (left panel),
$\epsilon(N)$ in eq. (32) (central panel) and $d\epsilon_{2}/dN$ (right panel)
as function of the number of $e$-folds $N$. We explore different values of
$\delta N$ with the limit $\delta N\to 0$ that corresponds to instantaneous
transitions SR/USR at $N=N_{\textrm{\rm in}}$ and USR/SR at $N=N_{\textrm{\rm
end}}$. In the right panel, the limit $\delta N\to 0$ corresponds to $\delta$
function transitions at $N_{\textrm{\rm in}}$ and $N_{\textrm{\rm end}}$.
Furthermore, notice that the lines corresponding to $d\epsilon_{2}/dN$ and
$-2d\eta/dN$ superimpose, showing that $-2d\eta/dN$ is a perfect approximation
of $d\epsilon_{2}/dN$.

In this way we have an analytical description of the background dynamics; most
importantly, eqs. (29, 32) are almost all that we need to know to solve the
M-S equation Sasaki (1986); Mukhanov (1988)

$\displaystyle\frac{d^{2}u_{k}}{dN^{2}}$
$\displaystyle+(1-\epsilon)\frac{du_{k}}{dN}+\left[\frac{k^{2}}{(aH)^{2}}+(1+\epsilon-\eta)(\eta-2)-\frac{d}{dN}(\epsilon-\eta)\right]u_{k}=0\,.$
(33)

For consistency with eq. (32), we consider the Hubble rate as a function of
time according to

$\displaystyle
H(N)=H(N_{\textrm{ref}})\exp\left[-\int_{N_{\textrm{ref}}}^{N}\epsilon(N^{\prime})dN^{\prime}\right]\,.$
(34)

We shall use the short-hand notation $a(N_{\textrm{i}})\equiv a_{\textrm{i}}$
and $H(N_{\textrm{i}})\equiv H_{\textrm{i}}$. Consequently, we rewrite eq.
(33) in the form

$\displaystyle\frac{d^{2}u_{k}}{dN^{2}}$
$\displaystyle+(1-\epsilon)\frac{du_{k}}{dN}+\left[\bar{k}^{2}\bigg{(}\frac{H_{\textrm{in}}}{H}\bigg{)}^{2}e^{2(N_{\textrm{in}}-N)}+(1+\epsilon-\eta)(\eta-2)-\frac{d}{dN}(\epsilon-\eta)\right]u_{k}=0\,,$
(35)

with $\bar{k}\equiv{k}/({a_{\textrm{in}}H_{\textrm{in}}})$. We solve the above
equation for different $\bar{k}$ with Bunch-Davies initial conditions

$\displaystyle\sqrt{k}\,u_{k}(N)=\frac{1}{\sqrt{2}}\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\sqrt{k}\,\frac{du_{k}}{dN}(N)=-\frac{i}{\sqrt{2}}\frac{k}{a(N)H(N)}\,,$
(36)

at some arbitrary time $N\ll N_{k}$ with $k=a(N_{k})H(N_{k})$. Modes with
$\bar{k}\approx O(1)$ are modes that exit the Hubble horizon at about the
beginning of the USR phase, modes with $\bar{k}\ll 1$ are modes that exit the
Hubble horizon well before the beginning of the USR phase, modes with
$\bar{k}\gg 1$ are modes that exit the Hubble horizon well after the beginning
of the USR phase.

Figure 3: Left panel: Tree-level power spectrum in the minimal dynamics of
section II.2. The numerical values of the other parameters are
$\eta_{\textrm{\rm II}}=3.5$, $\eta_{\textrm{\rm III}}=0$ and $N_{\textrm{\rm
end}}-N_{\textrm{\rm in}}=2.5$. In our parametrization, we go beyond the
instantaneous transition approximation and we explore different values of
$\delta N$. The vertical double-arrow indicates the growth of the power
spectrum given by the naïve scaling $\Delta\mathcal{P}=(k_{\rm end}/k_{\rm
in})^{2\eta_{\rm II}}=e^{2\eta_{\rm II}\Delta N_{\rm{USR}}}$. This scaling
captures well the amplitude of the transition from the initial to the final SR
phase but it does not give a reliable estimate of the peak amplitude of the
power spectrum, which can easily be one order of magnitude larger. Right
panel: Time evolution of two representative modes with $\bar{k}=1$ and
$\bar{k}=e^{\Delta N_{\rm{USR}}}$ for $\delta N\in[0.1\div 0.5]$ (from darker
to lighter colors, respectively). The black lines represent the limit $\delta
N\to 0$.

In the left panel of fig. 3, we show the tree-level power spectrum that we
obtain by numerically solving eq. (35) and using eq. (18). Thanks to our
parametrization in eq. (29) we control the sharpness of the transition varying
$\delta N$.

In order to make contact with the analysis of ref. Kristiano and Yokoyama
(2022), we set $\eta_{\rm III}=0$. However, it should be noted that in more
realistic models we need $\eta_{\textrm{III}}\neq 0$ and negative so that the
power spectrum decreases for modes with $\bar{k}\gg 1$. This feature is
necessary if we want to connect the USR phase to a subsequent SR phase that
ends inflation. Since we are considering single-field models of inflation, in
our analysis this is a necessary requirement. Consequently, the power spectrum
at small scales – both before and after the peak – does not respect the
property of scale invariance. Before the peak, the power spectrum of the short
modes grows with a maximum slope given by
$\mathcal{P}(\bar{k})\sim\bar{k}^{4}$; after the peak, the power spectrum of
the short modes decays approximately as
$\mathcal{P}(\bar{k})\sim\bar{k}^{2\eta_{\textrm{III}}}$. After the peak,
therefore, the power spectrum becomes approximately scale invariant only if we
take $\eta_{\textrm{III}}\approx 0$; however, in such case $\epsilon$ remains
anchored to the tiny value reached during the USR phase and inflation never
ends.

In ref. Kristiano and Yokoyama (2022), the loop integration is restricted to
the interval of modes
$\bar{k}\in[\bar{k}_{\textrm{in}},\bar{k}_{\textrm{end}}]$ where
$\bar{k}_{\textrm{in}}=1$ and $\bar{k}_{\textrm{end}}=e^{\Delta
N_{\textrm{USR}}}(H_{\textrm{end}}/H_{\textrm{in}})\simeq e^{\Delta
N_{\textrm{USR}}}$ with $\Delta N_{\textrm{USR}}\equiv
N_{\textrm{end}}-N_{\textrm{in}}$. This interval of modes is limited by the
two vertical dashed lines in the left panel of fig. 3. In ref. Kristiano and
Yokoyama (2022), limiting the integration to the range
$\bar{k}\in[\bar{k}_{\textrm{in}},\bar{k}_{\textrm{end}}]$ is justified by the
fact that the power spectrum of short modes peaks in this window of modes.

For future reference, let us stress one more important point. In the left
panel of fig. 3 we indicate the growth of the power spectrum given by the
scaling $\Delta\mathcal{P}=(k_{\rm end}/k_{\rm in})^{2\eta_{\rm II}}$. This
result immediately follows from the application of the SR formula
$\mathcal{P}(k)=H^{2}/8\pi^{2}\epsilon$ if one accounts for the exponential
decay $\epsilon\sim e^{-2\eta_{\rm II}N}$ during USR and converts $N$ into $k$
by means of the horizon-crossing condition $k=aH$. Therefore, not
surprisingly, the scaling $\Delta\mathcal{P}=(k_{\rm end}/k_{\rm
in})^{2\eta_{\rm II}}$ captures well the growth of the power spectrum if one
directly jumps from the initial to the final SR phase. However, as shown in
the left panel of fig. 3, the above estimate does not accurately describe the
amplitude of the power spectrum at the position of its peak; the latter can
easily be one order of magnitude larger than what suggested by
$\Delta\mathcal{P}=(k_{\rm end}/k_{\rm in})^{2\eta_{\rm II}}$. This features
has important consequences when estimating the PBH abundance, which rather
sensitive to the spectral amplitude. We will come back to this point in the
next section.

Finally, it is possible to check numerically that neglecting the time
dependence of the Hubble rate as in eq. (34) has a negligible impact. In the
following, therefore, we shall keep $H$ constant (that is,
$H=H_{\textrm{ref}}$ does not evolve in time). Furthermore, if we take $H$
constant and in the limit $\delta N=0$, it is possible to get, for some
special values of $\eta_{\textrm{II}}$ and $\eta_{\textrm{III}}$, a complete
analytical description of the SR/USR/SR dynamics Byrnes  _et al._ (2019);
Ballesteros  _et al._ (2020).

###  II.3 The cubic action

At the cubic order in the fluctuations, the action is

$\displaystyle\mathcal{S}_{3}=\int d^{4}x\bigg{\\{}$
$\displaystyle\epsilon^{2}a^{3}\dot{\zeta}^{2}\zeta+\epsilon^{2}a\zeta(\partial_{k}\zeta)(\partial^{k}\zeta)-2\epsilon^{2}a^{3}\dot{\zeta}(\partial_{k}\zeta)\partial^{k}(\partial^{-2}\dot{\zeta})+\frac{\epsilon\dot{\epsilon_{2}}}{2}a^{3}\dot{\zeta}\zeta^{2}-\frac{a^{3}\epsilon^{3}}{2}\big{[}\dot{\zeta}^{2}\zeta-\zeta\partial_{k}\partial_{l}(\partial^{-2}\dot{\zeta})\partial^{k}\partial^{l}(\partial^{-2}\dot{\zeta})\big{]}$
$\displaystyle+\bigg{[}\frac{d}{dt}\left(\epsilon
a^{3}\dot{\zeta}\right)-\epsilon
a\partial_{k}\partial^{k}\zeta\bigg{]}\bigg{[}\frac{\epsilon_{2}}{2}\zeta^{2}+\frac{2}{H}\dot{\zeta}\zeta-\frac{1}{2a^{2}H^{2}}(\partial_{k}\zeta)(\partial^{k}\zeta)+\frac{1}{2a^{2}H^{2}}\partial^{-2}\partial_{k}\partial_{l}(\partial^{k}\zeta\partial^{l}\zeta)$
$\displaystyle\hskip
116.65646pt+\frac{\epsilon}{H}(\partial_{k}\zeta)\partial^{k}(\partial^{-2}\dot{\zeta})-\frac{\epsilon}{H}\partial^{-2}\partial_{k}\partial_{l}\partial^{k}\zeta\partial^{l}(\partial^{-2}\dot{\zeta})\bigg{]}\bigg{\\}}\,.$
(37)

As shown in ref. Maldacena (2003), it is possible to simplify the cubic action
by means of a field redefinition that introduces a non-linear shift in the
original field. Concretely, if we define

$\displaystyle\zeta\equiv\zeta_{n}+f(\zeta_{n})\,,$ (38)

with

$\displaystyle f($
$\displaystyle\zeta)\equiv\frac{1}{2}\bigg{[}\frac{\epsilon_{2}}{2}\zeta^{2}+\frac{2}{H}\dot{\zeta}\zeta-\frac{(\partial_{k}\zeta)(\partial^{k}\zeta)}{2a^{2}H^{2}}+\frac{1}{2a^{2}H^{2}}\partial^{-2}\partial_{k}\partial_{l}(\partial^{k}\zeta\partial^{l}\zeta)+\frac{\epsilon}{H}(\partial_{k}\zeta)\partial^{k}(\partial^{-2}\dot{\zeta})-\frac{\epsilon}{H}\partial^{-2}\partial_{k}\partial_{l}\partial^{k}\zeta\partial^{l}(\partial^{-2}\dot{\zeta})\bigg{]}\,,$
(39)

we find, by direct computation, that at the quadratic order the field
$\zeta_{n}$ is described by the action

$\displaystyle\mathcal{S}_{2}(\zeta_{n})=\int
d^{4}x\,\epsilon\,a^{3}\left[\dot{\zeta}_{n}^{2}-\frac{(\partial_{k}\zeta_{n})(\partial^{k}\zeta_{n})}{a^{2}}\right]\,,$
(40)

which has the same structure as the quadratic action for the original variable
$\zeta$. However, at the cubic order, we find

$\displaystyle\mathcal{S}_{3}(\zeta_{n})=\int d^{4}x\bigg{\\{}$
$\displaystyle\epsilon^{2}a^{3}\dot{\zeta}_{n}^{2}\zeta_{n}+\epsilon^{2}a\zeta_{n}(\partial_{k}\zeta_{n})(\partial^{k}\zeta_{n})-2\epsilon^{2}a^{3}\dot{\zeta}_{n}(\partial_{k}\zeta_{n})\partial^{k}(\partial^{-2}\dot{\zeta}_{n})+\frac{\epsilon\dot{\epsilon_{2}}}{2}a^{3}\dot{\zeta}_{n}\zeta_{n}^{2}$
$\displaystyle-\frac{a^{3}\epsilon^{3}}{2}\big{[}\dot{\zeta}_{n}^{2}\zeta_{n}-\zeta_{n}\partial_{k}\partial_{l}(\partial^{-2}\dot{\zeta}_{n})\partial^{k}\partial^{l}(\partial^{-2}\dot{\zeta}_{n})\big{]}\bigg{\\}}\,,$
(41)

in which, thanks to the above field redefinition, the second and third lines
in eq. (37) cancel out.

If we neglect terms with spatial derivatives and interactions suppressed by
two or more powers of the Hubble parameter $\epsilon$, we find

$\displaystyle\mathcal{S}_{3}(\zeta_{n})\ni\int
d^{4}x\,\frac{\epsilon\dot{\epsilon_{2}}}{2}a^{3}\dot{\zeta}_{n}\zeta_{n}^{2}\,.$
(42)

Notice that we do not count the coupling $\epsilon_{2}$ as a slow-roll
suppression since we are interested in the USR phase during which
$|\epsilon_{2}|>3$ and $\dot{\epsilon}_{2}\neq 0$. Eq. (42) is the only
interaction included in ref. Kristiano and Yokoyama (2022). This means that,
implicitly, ref. Kristiano and Yokoyama (2022) computes the two-point function
for the field $\zeta_{n}$. This is because, in terms of the dynamical variable
$\zeta$, there is another interaction of order $\epsilon\epsilon_{2}$ that
should be included, that is the one in the second line of eq. (37).

However, as stressed in ref.Maldacena (2003), $\zeta_{n}$ is not the right
dynamical variable to consider since it is not conserved outside the horizon.
This is a trivial consequence of eq. (38). Since $\zeta$ is conserved outside
the horizon, $\zeta_{n}$ can not be conserved simply because various
coefficients in the non-linear relation in eq. (38) are time-dependent.
Alternatively, as discussed in ref.Maldacena (2003), the above fact is also
evident from the very same structure of the interactions that appear in eq.
(41). The interaction $\epsilon\dot{\epsilon}_{2}\dot{\zeta}_{n}\zeta_{n}^{2}$
only has one time-derivative acting on the field $\zeta_{n}$; consequently, it
alters the value of $\zeta_{n}$ on super-horizon scales (if one computes the
equation of motion for $\zeta_{n}$, it is easy to see that the constant
solution is not stable). Let us make the above considerations more concrete.
Eventually, we are interested in the computation of the two-point function for
the original curvature field. Given the field redefinition in eq. (38), we
write

$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle=$
$\displaystyle\langle\big{\\{}\hat{\zeta}_{n}(\vec{x}_{1},\tau)+f[\hat{\zeta}_{n}(\vec{x}_{1},\tau)]\big{\\}}\big{\\{}\hat{\zeta}_{n}(\vec{x}_{2},\tau)+f[\hat{\zeta}_{n}(\vec{x}_{2},\tau)]\big{\\}}\rangle$
$\displaystyle=$
$\displaystyle\langle\hat{\zeta}_{n}(\vec{x}_{1},\tau)\hat{\zeta}_{n}(\vec{x}_{2},\tau)\rangle+$
(43)
$\displaystyle\langle\hat{\zeta}_{n}(\vec{x}_{1},\tau)f[\hat{\zeta}_{n}(\vec{x}_{2},\tau)]\rangle+\langle
f[\hat{\zeta}_{n}(\vec{x}_{1},\tau)]\hat{\zeta}_{n}(\vec{x}_{1},\tau)\rangle+$
(44) $\displaystyle\langle
f[\hat{\zeta}_{n}(\vec{x}_{1},\tau)]f[\hat{\zeta}_{n}(\vec{x}_{2},\tau)]\rangle\,,$
(45)

The first term, eq. (43), corresponds to the two-point function for the
shifted curvature field whose cubic action is given by eq. (41);
$\langle\hat{\zeta}_{n}(\vec{x}_{1},\tau)\hat{\zeta}_{n}(\vec{x}_{2},\tau)\rangle$
can be computed perturbatively by means of the “in-in” formalism sketched at
the end of section II.1. Eqs. (44, 45) account for the difference between
$\zeta$ and $\zeta_{n}$ at the non-linear level. Notice that the first term in
the functional form in eq. (39) does not die off in the late-time limit
$\tau\to 0^{-}$ (in which the power spectrum must be eventually evaluated) if
we consider the case in which $\epsilon_{2}\neq 0$ after the USR phase (as
expected in realistic single-field models, cf. section II.2). However, if we
limit to the case in which $\eta_{\rm III}=0$ the contribution from the field
redefinition vanishes. This limit was considered in ref. Kristiano and
Yokoyama (2022). In order to make contact with the analysis presented in ref.
Kristiano and Yokoyama (2022), we shall also adopt in the bulk of this work
the assumption $\eta_{\rm III}=0$.

Let us now come back to the schematic in eq. (26). The cubic Hamiltonian
interaction that follows from eq. (42) gives rise to the last two topologies
of connected diagrams illustrated in eq. (26). As in ref. Kristiano and
Yokoyama (2022), we will only focus on the 1PI diagram, that is, the central
diagram in eq. (26). The last diagram in eq. (26) consists of a tadpole that
is attached to a $\zeta$-propagator and affects at one-loop its two-point
correlation function. The correct way to deal with tadpoles is by changing the
background solution, cf. ref. Sloth (2007) for a discussion in the case of
ordinary SR inflation and ref. Senatore and Zaldarriaga (2010) for the case in
which there are additional spectator fields. Recently, ref. Inomata  _et al._
(2023) estimated the tadpole correction to the background evolution in the
context of a model in which there is a resonant amplification of field
fluctuations. Imposing the condition that such modification is negligible
could give rise to an additional perturbativity bound. We postpone a
comprehensive exploration of this issue in the context of realistic USR
dynamics to future work, cf. section V.

###  II.4 Beyond the cubic action

Before proceeding, we comment about quartic interactions since, as
qualitatively discussed in eq. (26), they give rise to one-loop corrections
which are of the same order if compared to those generated by cubic
interaction terms. The derivation of the fourth-order action has been
discussed in ref. Jarnhus and Sloth (2008). Based on this result, ref.
Kristiano and Yokoyama (2022, 2023) claims that the relevant quartic
interaction in the case of USR dynamics (that is, the quartic interaction
proportional to $\dot{\epsilon_{2}}$) gives a vanishing contribution when
inserted in eq. (24). On the contrary, ref. Firouzjahi (2023) does include
quartic interactions using an approach based on the effective field theory of
inflation and claims that they give a non-trivial contribution to the loop-
corrected power spectrum. However, ref. Firouzjahi (2023) does not clarify the
origin of the discrepancy with ref. Kristiano and Yokoyama (2022). Generally
speaking, we expect cubic and quartic interactions to be inextricably linked.
For instance, the quartic Hamiltonian receives a contribution that arises from
the modification of the conjugate momentum if there are cubic interactions
which depend on $\dot{\zeta}$. Similarly, cubic interactions with spatial
derivatives are paired with quartic interactions induced by a residual spatial
conformal symmetry of the perturbed metric Senatore and Zaldarriaga (2010). En
route, we notice that interactions with spatial derivatives are usually
neglected for modes that are super-horizon. However, in the spirit of the loop
computation in ref. Kristiano and Yokoyama (2022), the momenta over which the
loop is integrated cross the horizon during the USR phase, and, naïvely, their
spatial derivatives do not pay any super-horizon suppression.

In this work, as a preliminary step towards a more complete analysis and in
order to compare our results with the claim made in refs. Riotto (2023a);
Kristiano and Yokoyama (2023); Riotto (2023b); Firouzjahi (2023); Firouzjahi
and Riotto (2023); Choudhury  _et al._ (2023a, b, c), we only focus on the
cubic interaction in eq. (42). However, we stress that all the arguments
listed above motivate the need for a more comprehensive analysis. We postpone
this task to future work, cf. section V.

##  III One-loop computation

We consider in this section the cubic interaction Hamiltonian given by (we
omit the subscript I in the interaction-picture fields)

$\displaystyle\hat{H}_{\rm int}^{(3)}(\tau)=\frac{1}{2}\int
d^{3}\vec{x}\,\epsilon(\tau)\epsilon_{2}^{\prime}(\tau)a^{2}(\tau)\zeta^{\prime}(\vec{x},\tau)\zeta(\vec{x},\tau)^{2}\,.$
(46)

We consider eq. (25); this can be written in the compact form

$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)$
$\displaystyle\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}=\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}^{(1,1)}-2{\rm
Re}\left[\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}^{(0,2)}\right]$
(47)

where

$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}^{(1,1)}\equiv$
$\displaystyle\int_{-\infty(1+i\epsilon)}^{\tau}d\tau_{1}\int_{-\infty(1-i\epsilon)}^{\tau}d\tau_{2}\langle
0|\hat{H}_{\textrm{int}}(\tau_{1})\hat{\zeta}_{I}(\vec{x}_{1},\tau)\hat{\zeta}_{I}(\vec{x}_{2},\tau)\hat{H}_{\textrm{int}}(\tau_{2})|0\rangle\,,$
(48)
$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}^{(0,2)}\equiv$
$\displaystyle\int_{-\infty(1-i\epsilon)}^{\tau}d\tau_{1}\int_{-\infty(1-i\epsilon)}^{\tau_{1}}d\tau_{2}\langle
0|\hat{\zeta}_{I}(\vec{x}_{1},\tau)\hat{\zeta}_{I}(\vec{x}_{2},\tau)\hat{H}_{\textrm{int}}(\tau_{1})\hat{H}_{\textrm{int}}(\tau_{2})|0\rangle\,.$
(49)

This expansion is consistent with Eq. (16) of Ref. Senatore and Zaldarriaga
(2010). Consider the first contribution in eq. (48), one finds

$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}^{(1,1)}=\frac{1}{4}$
$\displaystyle\int_{-\infty_{+}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{-\infty_{-}}^{\tau}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int
d^{3}\vec{y}d^{3}\vec{z}$
$\displaystyle\int\left[\prod_{i=1}^{8}\frac{d^{3}\vec{k}_{i}}{(2\pi)^{3}}\right]e^{i\vec{y}\cdot(\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3})}e^{i(\vec{x}_{1}\cdot\vec{k}_{4}+\vec{x}_{2}\cdot\vec{k}_{5})}e^{i\vec{z}\cdot(\vec{k}_{6}+\vec{k}_{7}+\vec{k}_{8})}$
$\displaystyle\langle
0|\hat{\zeta}^{\prime}_{I}(\vec{k}_{1},\tau_{1})\hat{\zeta}_{I}(\vec{k}_{2},\tau_{1})\hat{\zeta}_{I}(\vec{k}_{3},\tau_{1})\hat{\zeta}_{I}(\vec{k}_{4},\tau)\hat{\zeta}_{I}(\vec{k}_{5},\tau)\hat{\zeta}^{\prime}_{I}(\vec{k}_{6},\tau_{2})\hat{\zeta}_{I}(\vec{k}_{7},\tau_{2})\hat{\zeta}_{I}(\vec{k}_{8},\tau_{2})|0\rangle\,.$
(50)

The $36$ connected Wick contractions can be expressed as

$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}^{(1,1)}=$
$\displaystyle\int_{-\infty_{+}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{-\infty_{-}}^{\tau}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\frac{d^{3}\vec{q}}{(2\pi)^{3}}\,e^{i(\vec{x}_{1}-\vec{x}_{2})\cdot(\vec{k}+\vec{q})}$
$\displaystyle\big{[}\left|\zeta_{k+q}(\tau)\right|^{2}\big{\\{}\zeta_{k}(\tau_{1})\zeta^{\prime}_{k+q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{k+q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+$
$\displaystyle\hskip
56.9055pt\zeta_{k}(\tau_{1})\zeta^{\prime}_{k+q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k+q}^{*}(\tau_{2})\big{[}\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}+$
$\displaystyle\hskip
56.9055pt\zeta_{k}^{*}(\tau_{2})\zeta^{\prime\,*}_{k+q}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{k+q}(\tau_{1})\big{[}\zeta_{k}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})+\zeta_{k}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}+$
$\displaystyle\hskip
56.9055pt\zeta_{k}^{\prime}(\tau_{1})\zeta_{k+q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k+q}^{*}(\tau_{2})\big{[}\zeta_{q}^{*}(\tau_{2})\zeta_{k}^{\prime\,*}(\tau_{2})+\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}+$
$\displaystyle\hskip
56.9055pt\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{k+q}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{k+q}(\tau_{1})\big{[}\zeta_{q}(\tau_{1})\zeta_{k}^{\prime}(\tau_{1})+\zeta_{k}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}\big{\\}}\big{]}\,.$
(51)

Consider now eq. (49). One can write it in the form

$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}^{(0,2)}=\frac{1}{4}$
$\displaystyle\int_{-\infty_{-}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{-\infty_{-}}^{\tau_{1}}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int
d^{3}\vec{y}d^{3}\vec{z}$
$\displaystyle\int\left[\prod_{i=1}^{8}\frac{d^{3}\vec{k}_{i}}{(2\pi)^{3}}\right]e^{i\vec{y}\cdot(\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3})}e^{i(\vec{x}_{1}\cdot\vec{k}_{4}+\vec{x}_{2}\cdot\vec{k}_{5})}e^{i\vec{z}\cdot(\vec{k}_{6}+\vec{k}_{7}+\vec{k}_{8})}$
$\displaystyle\langle
0|\hat{\zeta}_{I}(\vec{k}_{4},\tau)\hat{\zeta}_{I}(\vec{k}_{5},\tau)\hat{\zeta}^{\prime}_{I}(\vec{k}_{1},\tau_{1})\hat{\zeta}_{I}(\vec{k}_{2},\tau_{1})\hat{\zeta}_{I}(\vec{k}_{3},\tau_{1})\hat{\zeta}^{\prime}_{I}(\vec{k}_{6},\tau_{2})\hat{\zeta}_{I}(\vec{k}_{7},\tau_{2})\hat{\zeta}_{I}(\vec{k}_{8},\tau_{2})|0\rangle\,.$
(52)

After Wick contractions, we find

$\displaystyle\langle\hat{\zeta}(\vec{x}_{1},\tau)\hat{\zeta}(\vec{x}_{2},\tau)\rangle_{2^{\textrm{nd}}}^{(0,2)}=$
$\displaystyle\int_{-\infty_{-}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{-\infty_{-}}^{\tau_{1}}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\frac{d^{3}\vec{q}}{(2\pi)^{3}}e^{i(\vec{x}_{1}-\vec{x}_{2})\cdot(\vec{k}+\vec{q})}$
$\displaystyle\big{[}\zeta_{k+q}^{2}(\tau)\big{\\{}\zeta_{k}(\tau_{1})\zeta^{\prime\,*}_{k+q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{k+q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+$
$\displaystyle\hskip
42.67912pt\zeta_{k}(\tau_{1})\zeta^{\prime\,*}_{k+q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k+q}^{*}(\tau_{2})\big{[}\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}+$
$\displaystyle\hskip
42.67912pt\zeta_{k}^{*}(\tau_{2})\zeta^{\prime\,*}_{k+q}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{k+q}^{*}(\tau_{1})\big{[}\zeta_{k}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})+\zeta_{k}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}+$
$\displaystyle\hskip
42.67912pt\zeta_{k}^{\prime}(\tau_{1})\zeta_{k+q}^{*}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k+q}^{*}(\tau_{2})\big{[}\zeta_{q}^{*}(\tau_{2})\zeta_{k}^{\prime\,*}(\tau_{2})+\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}+$
$\displaystyle\hskip
42.67912pt\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{k+q}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{k+q}^{*}(\tau_{1})\big{[}\zeta_{q}(\tau_{1})\zeta_{k}^{\prime}(\tau_{1})+\zeta_{k}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}\big{\\}}\big{]}\,.$
(53)

At this point we shift the momentum following the prescription $k\rightarrow
k-q$ in such a way that $k$ is identified with the external momentum. The
power spectrum at one loop can be therefore written as

$\mathcal{P}(k)=\lim_{\tau\to
0^{-}}\left(\frac{k^{3}}{2\pi^{2}}\right)\left\\{\left|\zeta_{k}(\tau)\right|^{2}+\frac{1}{(4\pi)^{2}}\left[\Delta
P_{1}(k,\tau)+\Delta P_{2}(k,\tau)\right]\right\\}\,,$ (54)

with

$\displaystyle\Delta P_{1}(k,\tau)$ $\displaystyle\equiv
4\int_{-\infty_{+}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{-\infty_{-}}^{\tau}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int_{0}^{\infty}dq\,q^{2}\,d(\cos\theta)\left|\zeta_{k}(\tau)\right|^{2}$
$\displaystyle\times\big{\\{}\zeta_{k-q}(\tau_{1})\zeta^{\prime}_{k}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k-q}^{*}(\tau_{2})\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+$
$\displaystyle\hskip
17.64056pt\zeta_{k-q}(\tau_{1})\zeta^{\prime}_{k}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\big{[}\zeta_{k-q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+\zeta_{k-q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}+$
$\displaystyle\hskip
17.64056pt\zeta_{k-q}^{*}(\tau_{2})\zeta^{\prime\,*}_{k}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{k}(\tau_{1})\big{[}\zeta_{k-q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})+\zeta_{k-q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}+$
$\displaystyle\hskip
17.64056pt\zeta_{k-q}^{\prime}(\tau_{1})\zeta_{k}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\big{[}\zeta_{q}^{*}(\tau_{2})\zeta_{k-q}^{\prime\,*}(\tau_{2})+\zeta_{k-q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}+$
$\displaystyle\hskip
17.64056pt\zeta_{k-q}^{\prime\,*}(\tau_{2})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{k}(\tau_{1})\big{[}\zeta_{q}(\tau_{1})\zeta_{k-q}^{\prime}(\tau_{1})+\zeta_{k-q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}\big{\\}}\,,$
(55) $\displaystyle\Delta P_{2}(k,\tau)$ $\displaystyle\equiv-8{\rm
Re}\Big{[}\int_{-\infty_{-}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{-\infty_{-}}^{\tau_{1}}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int_{0}^{\infty}dq\,q^{2}\,d(\cos\theta)\zeta_{k}^{2}(\tau)$
$\displaystyle\times\big{\\{}\zeta_{k-q}(\tau_{1})\zeta^{\prime\,*}_{k}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k-q}^{*}(\tau_{2})\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+$
$\displaystyle\hskip
17.64056pt\zeta_{k-q}(\tau_{1})\zeta^{\prime\,*}_{k}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\big{[}\zeta_{k-q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+\zeta_{k-q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}+$
$\displaystyle\hskip
17.64056pt\zeta_{k-q}^{*}(\tau_{2})\zeta^{\prime\,*}_{k}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{k}^{*}(\tau_{1})\big{[}\zeta_{k-q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})+\zeta_{k-q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}+$
$\displaystyle\hskip
17.64056pt\zeta_{k-q}^{\prime}(\tau_{1})\zeta_{k}^{*}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\big{[}\zeta_{q}^{*}(\tau_{2})\zeta_{k-q}^{\prime\,*}(\tau_{2})+\zeta_{k-q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}+$
$\displaystyle\hskip
17.64056pt\zeta_{k-q}^{\prime\,*}(\tau_{2})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{k}^{*}(\tau_{1})\big{[}\zeta_{q}(\tau_{1})\zeta_{k-q}^{\prime}(\tau_{1})+\zeta_{k-q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}\big{\\}}\Big{]}\,.$
(56)

###  III.1 Loop correction with a large hierarchy of scales

First, we will be concerned with external momenta that describe the large CMB
scales, while the USR takes place when modes $k_{\text{\tiny USR}}\gg k$ cross
the horizon. The situation is summarized in the following schematic

(57)

in which the blue horizontal band represents the interval of modes that cross
the horizon during the USR phase, the vertical band shaded in magenta. In
other words, as we will restrict the integration over momenta $q\in[q_{\rm
in},q_{\rm end}]$ that are enhanced by the USR phase, we can assume $q\gg k$.
Consequently, as in ref. Kristiano and Yokoyama (2022), we approximate

$\displaystyle k-q=\sqrt{k^{2}+q^{2}-2kq\cos(\theta)}\approx
q\,,~{}~{}~{}~{}~{}~{}~{}\int_{-1}^{+1}d(\cos\theta)=2\,.$ (58)

With these assumptions, we can further simplify the expressions. We collect
each contribution depending on the number of time derivatives acting on the
long mode $\zeta_{k}$. In each expression, the first line indicates terms with
no derivative on the long modes, the second one those with one derivative,
while the last with two. One finds

$\displaystyle\Delta P_{1}(k,\tau)$ $\displaystyle\equiv 8\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{\tau_{\rm
in}}^{\tau}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int_{q_{\rm
in}}^{q_{\rm end}}dq\,q^{2}\,\left|\zeta_{k}(\tau)\right|^{2}$
$\displaystyle\times\big{\\{}4\zeta_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})+$
$\displaystyle\hskip
17.64056pt2\zeta^{\prime}_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+2\zeta_{k}(\tau_{1})\zeta^{\prime\,*}_{k}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+$
$\displaystyle\hskip
17.64056pt\zeta^{\prime}_{k}(\tau_{1})\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\big{\\}}\,,$
(59) $\displaystyle\Delta P_{2}(k,\tau)$ $\displaystyle\equiv-16{\rm
Re}\Big{[}\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int_{q_{\rm
in}}^{q_{\rm end}}dq\,q^{2}\,\zeta_{k}(\tau)^{2}$
$\displaystyle\times\big{\\{}4\zeta_{k}^{*}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})+$
$\displaystyle\hskip
17.64056pt2\zeta^{\prime\,*}_{k}(\tau_{1})\zeta^{*}_{k}(\tau_{2})\zeta_{q}(\tau_{1})^{2}\zeta_{q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+2\zeta^{*}_{k}(\tau_{1})\zeta^{\prime\,*}_{k}(\tau_{2})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})^{2}+$
$\displaystyle\hskip
17.64056pt\zeta^{\prime\,*}_{k}(\tau_{1})\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}(\tau_{1})^{2}\zeta_{q}^{*}(\tau_{2})^{2}\big{\\}}\Big{]}\,.$
(60)

We can combine the two contributions using the properties of symmetric
integrals for holomorphic symmetric functions
$f(\tau_{1},\tau_{2})=f(\tau_{2},\tau_{1})$ Adshead  _et al._ (2009b)

$\int_{\tau_{\rm in}}^{\tau}d\tau_{1}\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}f\left(\tau_{1},\tau_{2}\right)=\frac{1}{2}\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}\int_{\tau_{\rm
in}}^{\tau}d\tau_{2}f\left(\tau_{1},\tau_{2}\right).$ (61)

To shorten the notation, we introduce $\Delta P(k,\tau)=\Delta
P_{1}(k,\tau)+\Delta P_{2}(k,\tau)$ and collect the individual contribution
order by order in derivatives:

  * •

0th order in time derivatives of the long mode. For ease of reading, we
introduce the short-hand notation
$\epsilon(\tau)\epsilon_{2}^{\prime}(\tau)a^{2}(\tau)\equiv g(\tau)$. Consider
the sum of the two integrals

$\displaystyle\Delta P_{0{\rm th}}$ $\displaystyle(k,\tau)=32\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau}d\tau_{2}g(\tau_{2})\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\,\left|\zeta_{k}(\tau)\right|^{2}\zeta_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})$
(62) $\displaystyle-64{\rm Re}\Big{[}\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}g(\tau_{2})\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\,\zeta_{k}(\tau)^{2}\zeta_{k}^{*}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\Big{]}\,.$
(63)

We notice that, in the first integral in eq. (62), the exchange
$\tau_{1}\leftrightarrow\tau_{2}$ transforms

$\displaystyle\zeta_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\overset{\tau_{1}\leftrightarrow\tau_{2}}{\Longrightarrow}$
$\displaystyle\zeta_{k}(\tau_{2})\zeta_{k}^{*}(\tau_{1})\zeta_{q}(\tau_{2})\zeta_{q}^{\prime}(\tau_{2})\zeta_{q}^{*}(\tau_{1})\zeta_{q}^{\prime\,*}(\tau_{1})=$
$\displaystyle[\zeta_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})]^{*}\,.$
(64)

Therefore, the first integral in eq. (62) is fully symmetric under the
exchange $\tau_{1}\leftrightarrow\tau_{2}$, and we rewrite $\Delta P_{0{\rm
th}}(k,\tau)$ as

$\displaystyle\Delta P_{0{\rm th}}$ $\displaystyle(k,\tau)=$ $\displaystyle
32\int_{\tau_{\rm in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau}d\tau_{2}g(\tau_{2})\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\,\left|\zeta_{k}(\tau)\right|^{2}{\rm
Re}\big{[}\zeta_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}$
(65) $\displaystyle-64\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}g(\tau_{2})\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\,{\rm
Re}\big{[}\zeta_{k}(\tau)^{2}\zeta_{k}^{*}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}\,.$
(66)

and apply to the first integral in eq. (65) the identity in eq. (61). We
arrive at

$\displaystyle\Delta P_{0{\rm th}}(k,\tau)=64\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}g(\tau_{2})$ $\displaystyle\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\left\\{{\rm
Re}\big{[}\zeta_{k}(\tau)\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}\right.$
$\displaystyle-\left.{\rm
Re}\big{[}\zeta_{k}(\tau)^{2}\zeta_{k}^{*}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}\right\\}\,.$
(67)

We are now in the position of combining the two integrand functions.
Schematically, we define the two combinations

$\displaystyle
X\equiv\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{1})\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}Y\equiv\zeta_{k}(\tau)\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\,,$
(68)

such that the integrand in eq. (67) becomes

$\displaystyle{\rm Re}(XY)-{\rm Re}(X^{*}Y)=-2{\rm Im}(X){\rm Im}(Y)\,.$ (69)

We thus arrive at the result

Contribution with no time derivatives on the long mode $k$
$\displaystyle\Delta P_{0{\rm th}}(k,\tau)=$ $\displaystyle-128\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int_{q_{\rm
in}}^{q_{\rm end}}dq\,q^{2}$ $\displaystyle\times{\rm
Im}\left[\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{1})\right]{\rm
Im}\left[\zeta_{k}(\tau)\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\right].$
(70)

Given that we are interested in modes $k$ that are much smaller than the USR-
enhanced ones, they are super-horizon at the time of USR phase. Thus, for any
time $\tau\gtrsim\tau_{\rm in}$ of relevance for both time integrations, one
has that

${\rm Im}\left[\zeta_{k}(\tau)\zeta_{k}^{*}(\tau_{1})\right]\simeq{\rm
Im}\left[|\zeta_{k}(\tau)|^{2}\right]=0,$ (71)

which makes the above contribution negligible.

  * •

1st order in time derivatives of the long mode. Starting from the second lines
of eqs. (59,60), we now consider the sum

$\displaystyle\Delta P_{\rm 1st}(k,\tau)=$ $\displaystyle 16\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau}d\tau_{2}g(\tau_{2})\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\,\left|\zeta_{k}(\tau)\right|^{2}$
$\displaystyle~{}~{}~{}~{}\big{[}\zeta^{\prime}_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+\zeta_{k}(\tau_{1})\zeta^{\prime\,*}_{k}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\big{]}$
$\displaystyle-$ $\displaystyle 32{\rm Re}\Big{\\{}\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}g(\tau_{2})\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\zeta_{k}(\tau)^{2}$
$\displaystyle~{}~{}~{}~{}\big{[}\zeta^{\prime\,*}_{k}(\tau_{1})\zeta^{*}_{k}(\tau_{2})\zeta_{q}(\tau_{1})^{2}\zeta_{q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+\zeta^{*}_{k}(\tau_{1})\zeta^{\prime\,*}_{k}(\tau_{2})\zeta_{q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})^{2}\big{]}\Big{\\}}\,.$
(72)

Manipulations analogue to those discussed in the previous point allow one to
combine the two integrals together. We find

$\displaystyle\Delta P_{\rm 1st}(k,\tau)=-64\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}g(\tau_{2})$ $\displaystyle\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\left\\{{\rm
Im}\left[\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{1})\right]{\rm
Im}\left[\zeta_{k}(\tau)\zeta_{k}^{*}(\tau_{2})^{2}\zeta_{q}(\tau_{1})\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}^{\prime}(\tau_{1})\right]\right.$
$\displaystyle+\left.{\rm
Im}\left[\zeta_{k}^{*}(\tau)\zeta_{k}^{\prime}(\tau_{1})\right]{\rm
Im}\left[\zeta_{k}(\tau)\zeta_{q}(\tau_{1})^{2}\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\right]\right\\}.$
(73)

Again, since we are interested in modes $k$ that are much smaller than the
USR-enhanced ones, and are super-horizon at the time of USR phase, the
contribution within the curly brackets in the first line vanishes thanks to
eq. (71). This leaves us with

Contribution with one time derivatives on the long mode $k$
$\displaystyle\Delta P_{\rm 1st}(k,\tau)$
$\displaystyle\equiv-64\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int_{q_{\rm
in}}^{q_{\rm end}}dq\,q^{2}\,$ $\displaystyle\times{\rm
Im}\left[\zeta_{k}^{*}(\tau)\zeta_{k}^{\prime}(\tau_{1})\right]{\rm
Im}\left[\zeta_{k}(\tau)\zeta_{q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\right].$
(74)

  * •

2nd order in time derivatives of the long mode. Analogue manipulations give

Contribution with two time derivatives on the long mode $k$
$\displaystyle\Delta P_{\rm 2nd}(k,\tau)$
$\displaystyle\equiv-32\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}\epsilon(\tau_{1})\epsilon_{2}^{\prime}(\tau_{1})a^{2}(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}\epsilon(\tau_{2})\epsilon_{2}^{\prime}(\tau_{2})a^{2}(\tau_{2})\int_{q_{\rm
in}}^{q_{\rm end}}dq\,q^{2}\,$ $\displaystyle\times{\rm
Im}\left[\zeta_{k}^{*}(\tau)\zeta_{k}^{\prime}(\tau_{1})\right]{\rm
Im}\left[\zeta_{k}(\tau)\zeta_{q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\right].$
(75)

We stress that the only approximation employed so far is to take the external
momentum to be much smaller than the one in the loop, i.e. $k\ll q$, which is
justified in presence of a large hierarchy between the CMB and the USR scales.

###  III.2 Loop correction at any scales

It will be useful in the following to remove the assumption that the external
momentum is much smaller than the modes in the loop, i.e. the large separation
of scales $k\ll q$. Starting again from eqs. (54,55,56), we can proceed with
analogous steps as in the previous section and define

$\displaystyle X_{1}$
$\displaystyle\equiv\zeta_{k}^{*}(\tau)\zeta_{k}^{\prime}(\tau_{1})\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}Y_{1}\equiv\zeta_{k}(\tau)\zeta_{k-q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k-q}^{*}(\tau_{2})\zeta_{k}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\,,$
$\displaystyle X_{2}$
$\displaystyle\equiv\zeta_{k}^{*}(\tau)\zeta_{k}^{\prime}(\tau_{1})\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}Y_{2}\equiv\zeta_{k}(\tau)\zeta_{k-q}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\big{[}\zeta_{k-q}^{\prime\,*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})+\zeta_{k-q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}\,,$
$\displaystyle X_{3}$
$\displaystyle\equiv\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{1})\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}Y_{3}\equiv\zeta_{k}(\tau)\zeta_{k-q}^{*}(\tau_{2})\zeta^{\prime\,*}_{k}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\big{[}\zeta_{k-q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})+\zeta_{k-q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})$
$\displaystyle X_{4}$
$\displaystyle\equiv\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{1})\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}Y_{4}\equiv\zeta_{k}(\tau)\zeta_{k-q}^{\prime}(\tau_{1})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\big{[}\zeta_{q}^{*}(\tau_{2})\zeta_{k-q}^{\prime\,*}(\tau_{2})+\zeta_{k-q}^{*}(\tau_{2})\zeta_{q}^{\prime\,*}(\tau_{2})\big{]}$
$\displaystyle X_{5}$
$\displaystyle\equiv\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{1})\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}Y_{5}\equiv\zeta_{k}(\tau)\zeta_{k-q}^{\prime\,*}(\tau_{2})\zeta_{k}^{*}(\tau_{2})\zeta_{q}^{*}(\tau_{2})\big{[}\zeta_{q}(\tau_{1})\zeta_{k-q}^{\prime}(\tau_{1})+\zeta_{k-q}(\tau_{1})\zeta_{q}^{\prime}(\tau_{1})\big{]}$
(76)

in such a way that $\Delta P\equiv\Delta P_{1}+\Delta P_{2}$ can be written in
the schematic form

Generic loop correction at any scale $k$ $\displaystyle\Delta P(k,\tau)$
$\displaystyle\equiv-16\int_{\tau_{\rm
in}}^{\tau}d\tau_{1}g(\tau_{1})\int_{\tau_{\rm
in}}^{\tau_{1}}d\tau_{2}g(\tau_{2})\int_{q_{\rm in}}^{q_{\rm
end}}dq\,q^{2}\,\int_{-1}^{1}d(\cos\theta)\times\sum_{i=1}^{5}\mbox{Im}\,(X_{i})\mbox{Im}\,(Y_{i}),$
(77)

thanks to the identity in eq. (69) and where we again introduced
$\epsilon(\tau)\epsilon_{2}^{\prime}(\tau)a^{2}(\tau)\equiv g(\tau)$. This
expression is much more intricate than the one obtained in the limit of a
large hierarchy of scales between the mode $k$ and the USR loop momenta. It
will allow us to seize the loop correction to the power spectrum also at the
USR scales where the peak of the power spectrum is generated.

##  IV Time integration beyond the instantaneous transition and at any scales

###  IV.1 Loop evaluation at the CMB scales

Let us try to simplify the structure of eq. (54) in light of the
approximations introduced so far. First of all, let us write eq. (54) in the
form

$\displaystyle\mathcal{P}(k)=\frac{H^{2}}{8\pi^{2}\epsilon_{\textrm{ref}}}\left\\{1+\lim_{\tau\to
0^{-}}\frac{4\epsilon_{\textrm{ref}}k^{3}}{H^{2}(4\pi)^{2}}\Delta P_{\rm
1st}(k,\tau)+\lim_{\tau\to
0^{-}}\frac{4\epsilon_{\textrm{ref}}k^{3}}{H^{2}(4\pi)^{2}}\Delta P_{\rm
2nd}(k,\tau)\right\\}\,,$ (78)

where we used the slow-roll approximation for the first term in eq. (54) given
that $k$ is of the order of the CMB pivot scale. We focus on the leading
correction given by $\Delta P_{\rm 1st}(k,\tau)$. Using the number of
$e$-folds as the time variable, we find that it can be written in the compact
form (cf. our definition in eq. (1))

$\displaystyle\Delta\mathcal{P}_{\rm 1-loop}(k_{*})\equiv\lim_{\tau\to
0^{-}}\frac{4\epsilon_{\textrm{ref}}k^{3}}{H^{2}(4\pi)^{2}}\Delta P_{\rm
1st}(k,\tau)=$ $\displaystyle
32\left(\frac{H^{2}}{8\pi^{2}\epsilon_{\textrm{ref}}}\right)\int_{N_{\rm
in}-\Delta N}^{N_{\rm end}+\Delta N}dN_{1}\frac{d\eta}{dN}(N_{1})\int_{N_{\rm
in}-\Delta
N}^{N_{1}}dN_{2}\bar{\epsilon}(N_{2})\frac{d\eta}{dN}(N_{2})e^{3(N_{2}-N_{\rm
in})}\int\frac{d\bar{q}}{\bar{q}^{4}}{\rm
Im}\left[\bar{\zeta}_{q}(N_{1})^{2}\bar{\zeta}_{q}^{*}(N_{2})\frac{d\bar{\zeta}_{q}^{*}}{dN}(N_{2})\right]$
(79)

where we introduced the following manupulations:

  * i)

We use the approximation $\epsilon_{2}(N)\approx-2\eta(N)$. This is because in
the relevant range of $N$ over which we integrate $\epsilon\ll 1$ while
$\eta=O(1)$, cf. the right panel of fig. 2.

  * ii)

We define $\bar{q}\equiv q/a_{\textrm{in}}H$. Furthermore, we use the two
relations

$\displaystyle\frac{a(N_{1})H}{k}=e^{N_{1}-N_{k}}\,,~{}~{}~{}~{}\textrm{with:}~{}~{}~{}~{}a(N_{k})H=k\,,~{}~{}~{}~{}~{}~{}~{}\textrm{and}~{}~{}~{}~{}\frac{a(N_{2})H}{q}=\frac{e^{N_{2}-N_{\rm
in}}}{\bar{q}}\,.$ (80)

  * iii)

We introduce the short-hand notation

$\displaystyle\bar{\zeta}_{q}(N)\equiv\frac{\epsilon_{\textrm{ref}}^{1/2}q^{3/2}\zeta_{q}(N)}{H}\,.$
(81)

The virtue of this definition is that $\bar{\zeta}_{q}(N)$ is precisely the
quantity we compute numerically by solving the M-S equation, cf. the right
panel of fig. 3. Furthermore, it should be noted that the definition in eq.
(81) is automatically invariant under the rescaling in eq. (10). The same
comment applies to the definition of $\bar{q}$ and the ratios in eq. (80).
Consequently, an expression entirely written in terms of barred quantities is
automatically invariant under the rescaling in eq. (10).

  * iv)

Importantly, in the derivation of eq. (79) we use (cf. appendix A)

$\displaystyle\textrm{Im}\bigg{[}\bar{\zeta}_{k}^{*}(N)\frac{d\bar{\zeta}_{k}}{dN}(N_{1})\bigg{]}\simeq\textrm{Im}\bigg{[}\bar{\zeta}_{k}^{*}(N_{1})\frac{d\bar{\zeta}_{k}}{dN}(N_{1})\bigg{]}=-\frac{\bar{k}^{3}}{4\bar{\epsilon}(N_{1})}e^{3(N_{\textrm{in}}-N_{1})}\,,$
(82)

with $\epsilon(N)$ given by eq. (32) for generic $\delta N$. This is because
$N_{1,2}$ vary at around $N_{\rm end}$, and in this time interval modes with
comoving wavenumbers $k\approx k_{*}$ are way outside the horizon and stay
constant. For the very same reason, we also use the slow-roll approximation

$\displaystyle\bar{\zeta}_{k}(N_{1})\bar{\zeta}^{*}_{k}(N_{2})=\frac{1}{4}\,.$
(83)

  * v)

The range of integration in eq. (79) is as follows. In the case of a smooth
transition, we take $N_{1}\in[N_{\textrm{in}}-\Delta N,N_{\textrm{end}}+\Delta
N]$ and $N_{2}\in[N_{\textrm{in}}-\Delta N,N_{1}]$ where $\Delta N$ should be
large enough to complete the SR/USR/SR transition (that is, $\Delta
N\gtrsim\delta N$). In the limit of instantaneous transition, we set
$N_{1}=N_{2}=N_{\rm end}$, which corresponds to consider the dominant
contribution given by the first $\delta$ function in eq. (31). Moreover, we
include a factor $1/2$ since, with respect to the integration over $N_{2}$,
the argument of the $\delta$ function in eq. (31) picks up the upper limit of
the integration interval at $N_{\rm end}$. The integration over $q$, on the
contrary, is limited by $\bar{q}\in[1,e^{\Delta N_{\textrm{USR}}}]$.

####  IV.1.1 The instantaneous transition

We consider the instantaneous limit (dubbed $\delta N\to 0$ in the following)
of eq. (79). We find

$\displaystyle\lim_{\delta N\to 0}\Delta\mathcal{P}_{\rm
1-loop}(k_{*})=\left(\frac{H^{2}}{8\pi^{2}\epsilon_{\textrm{ref}}}\right)\eta_{\textrm{II}}^{2}\left(\frac{k_{\textrm{end}}}{k_{\textrm{in}}}\right)^{-2\eta_{\textrm{II}}+3}\lim_{\delta
N\to 0}(16)\int_{1}^{e^{\Delta
N_{\textrm{USR}}}}\frac{d\bar{q}}{\bar{q}^{4}}|\bar{\zeta}_{q}(N_{\rm
end})|^{2}{\rm Im}\left[\bar{\zeta}_{q}(N_{\rm
end})\frac{d\bar{\zeta}_{q}^{*}}{dN}(N_{\rm end})\right]\,.$ (84)

This expression can be further simplified using (cf. eq. (123))

$\displaystyle{\rm Im}\left[\bar{\zeta}_{q}(N_{\rm
end})\frac{d\bar{\zeta}_{q}^{*}}{dN}(N_{\rm
end})\right]=\frac{\bar{q}^{3}}{4}e^{(2\eta_{\textrm{II}}-3)(N_{\textrm{end}}-N_{\textrm{in}})}=\frac{\bar{q}^{3}}{4}\left(\frac{k_{\rm
end}}{k_{\rm in}}\right)^{2\eta_{\rm II}-3}\,,$ (85)

so that we write

$\displaystyle\lim_{\delta N\to 0}\Delta\mathcal{P}_{\rm
1-loop}(k_{*})=\left(\frac{H^{2}}{8\pi^{2}\epsilon_{\textrm{ref}}}\right)\,4\eta_{\textrm{II}}^{2}\,\lim_{\delta
N\to 0}\int_{1}^{e^{\Delta
N_{\textrm{USR}}}}\frac{d\bar{q}}{\bar{q}}|\bar{\zeta}_{q}(N_{\rm
end})|^{2}\,.$ (86)

We remark that this expression is valid for generic values of $\eta_{\rm II}$
during USR.

We consider now the computation of the last integral. The factor
$\bar{\zeta}_{q}$ grows exponentially during the USR phase. In the case of
sub-horizon modes, we have $\bar{\zeta}_{q}(N)\sim e^{-(1-\eta_{\rm II})N}$
while in the case of super-horizon modes we find $\bar{\zeta}_{q}(N)\sim
e^{-(3-2\eta_{\rm II})N}$ (cf. appendix A). However, the precise estimate of
the integral in eq. (86) is complicated by the fact that curvature modes
$\bar{\zeta}_{q}$ with $\bar{q}\in[1,\exp(\Delta N_{\textrm{USR}})]$ are
neither sub- nor super-horizon but they exit the horizon during the USR phase,
thus making the analytical estimate of the argument of their exponential
growth more challenging.

The situation simplifies if we consider some special values of
$\eta_{\textrm{II}}$. We consider the case $\eta_{\textrm{II}}=3$ (that is
$\epsilon_{2}=-6$, it should be noted that this is also the case studied in
ref. Kristiano and Yokoyama (2022)). In this case, everything can be computed
analytically. We find the scaling

$\displaystyle\lim_{\delta N\to 0}\int_{1}^{e^{\Delta
N_{\textrm{USR}}}}\frac{d\bar{q}}{\bar{q}}|\bar{\zeta}_{q}(N_{\rm
end})|^{2}\approx\frac{e^{6\Delta N_{\rm USR}}}{4}(1+\Delta N_{\rm
USR})=\frac{1}{4}\left(\frac{k_{\rm end}}{k_{\rm
in}}\right)^{6}\left[1+\log\left(\frac{k_{\rm end}}{k_{\rm
in}}\right)\right]\,,$ (87)

which becomes more and more accurate for larger $k_{\rm end}/k_{\rm in}$. The
final result is

Leading one-loop correction at CMB scales in the instantaneous SR/USR/SR
transition $\displaystyle\lim_{\delta N\to 0}\Delta\mathcal{P}_{\rm
1-loop}(k_{*})\approx\left(\frac{H^{2}}{8\pi^{2}\epsilon_{\textrm{ref}}}\right)\eta_{\rm
II}^{2}\left(\frac{k_{\textrm{end}}}{k_{\textrm{in}}}\right)^{6}\left[1+\log\left(\frac{k_{\rm
end}}{k_{\rm
in}}\right)\right]\,,~{}~{}~{}~{}~{}~{}~{}\textrm{with}~{}~{}\eta_{\textrm{II}}=3,\,\,\eta_{\rm
III}=0$ (88)

which perfectly agrees with the findings of ref. Kristiano and Yokoyama (2022)
in the same limit.

The above result has a number of limitations, which we address separately in
the following subsections:

  * $\circ$

Dynamics during USR, section IV.1.2. We modify the assumption $\eta_{\rm
II}=3$ and we take $\delta N\to 0$ and $\eta_{\rm III}=0$.

  * $\circ$

Dynamics at the SR/USR/SR transition, section IV.1.3. We consider $\delta
N\neq 0$, with generic $\eta_{\rm II}$ but $\eta_{\rm III}=0$. Considering a
non-zero value of $\delta N$ is very important because it corresponds to a
more realistic smooth SR/USR/SR transition, as opposed to the instantaneous
limit with $\delta N=0$.

####  IV.1.2 Dynamics during USR

We compute eq. (86) for generic values of $\eta_{\rm II}$, still keeping
$\delta N\to 0$ and $\eta_{\rm III}=0$. From the computation of the tree-level
power spectrum (cf. section II.2 and fig. 3) we define

$\displaystyle\frac{\mathcal{P}_{\rm USR}}{\mathcal{P}_{\rm
CMB}}\equiv\frac{\mathcal{P}(\bar{k}_{\rm max})}{\mathcal{P}(\bar{k}\ll
1)}\,,$ (89)

where $\bar{k}_{\rm max}$ represent the position of the max of
$\mathcal{P}(\bar{k})$ after the growth due to the USR dynamics.

Figure 4: In both panels, we consider a generic USR dynamics with varying
$\eta_{\rm II}$ ($x$-axis) and $\Delta N_{\rm USR}$ ($y$-axis). We take
$\eta_{\rm III}=0$ and the instantaneous limit $\delta N=0$. We plot in solid
red contours of constant $\lim_{\delta N\to 0}\Delta\mathcal{P}_{\rm
1-loop}(k_{*}),$ defined in eq. (79) and computed according to eq. (86) with
$H^{2}/8\pi^{2}\epsilon_{\rm ref}=2.1\times 10^{-9}$. The region hatched in
red is defined by the condition $\lim_{\delta N\to 0}\Delta\mathcal{P}_{\rm
1-loop}(k_{*})>1$. Left panel: We superimpose contours of constant
$\mathcal{P}_{\rm USR}/\mathcal{P}_{\rm CMB}$ as defined in eq. (89) (dashed
blue). Right panel: We superimpose the line defined by the condition $f_{\rm
PBH}=1$. Along this line, we get 100% of DM in the form of asteroid-mass PBHs.

We compare in the left panel of fig. 4 contours of constant $\mathcal{P}_{\rm
USR}/\mathcal{P}_{\rm CMB}$ (dashed blue) and constant $\lim_{\delta N\to
0}\Delta\mathcal{P}_{\rm 1-loop}(k_{*})$ (solid red). We take
$H^{2}/8\pi^{2}\epsilon_{\rm ref}=2.1\times 10^{-9}$. Our analysis shows that
enhancements $\mathcal{P}_{\rm USR}/\mathcal{P}_{\rm CMB}\gtrsim 10^{8}$ are
barely compatible with the perturbativity condition $\lim_{\delta N\to
0}\Delta\mathcal{P}_{\rm 1-loop}(k_{*})<1$, which roughly means “loops $<$
tree level”. The region $\lim_{\delta N\to 0}\Delta\mathcal{P}_{\rm
1-loop}(k_{*})>1$ is hatched in red in fig. 4.

We can actually do better and compare with a careful computation of the PBH
abundance. The parameters of the dynamics in section II.2 (with $\eta_{\rm
III}=0$ and $\delta N\to 0$) are chosen in such a way that the integral

$\displaystyle f_{\rm PBH}\equiv\frac{\Omega_{\rm PBH}}{\Omega_{\rm CDM}}=\int
f_{\rm PBH}(M_{\rm PBH})d\log M_{\rm PBH}\approx 1\,,$ (90)

which means that we get $\approx 100\%$ of DM in the form of PBHs. More in
detail, we tune, for each $\eta_{\rm II}$, the value of $\Delta N_{\rm USR}$
so to get $f_{\rm PBH}\approx 1$; we choose the numerical value of $k_{\rm
in}$ in such a way that the peak of the PBH mass distribution $f_{\rm
PBH}(M_{\rm PBH})$ falls within the interval $M_{\rm
PBH}/M_{\odot}\in[10^{-14},\,10^{-13}]$ in which the condition $f_{\rm
PBH}\approx 1$ is experimentally allowed, the so-called asteroid-mass PBHs
Carr and Kuhnel (2020). We compute eq. (90) using threshold statistics and
including the full non-linear relation between the curvature and the density
contrast fields (cf. Young  _et al._ (2019); De Luca  _et al._ (2019)). The
interested reader can find more details on the computation of the abundance in
ref. Franciolini and Urbano (2022) and refs. therein.333It is possible to
further improve our analysis by including the presence of primordial non-
gaussianity (e.g. Namjoo  _et al._ (2013); Chen  _et al._ (2013); Cai  _et
al._ (2018); Pi and Sasaki (2022); Passaglia  _et al._ (2019); Figueroa  _et
al._ (2021); Biagetti  _et al._ (2021); Figueroa  _et al._ (2022); Tomberg
(2023)). In the case of local non-guassianity parametrized by a positive non-
Guassian parameter $f_{\rm NL}$, as expected in the case of USR, we
generically expect a larger abundance of PBHs compared to the Gaussian case
Atal  _et al._ (2019); Ferrante  _et al._ (2023); Gow  _et al._ (2022). This
means that, in order to achieve the same abundance of PBHs, one needs a power
spectrum with a smaller peak amplitude. This argument implies that the
presence of primordial non-Gaussianity will tend to decrease the relevance of
the one-loop corrections.  In the right panel of fig. 4 we plot the line
defined by the condition $f_{\rm PBH}\approx 1$. The comparison between the
left- and right-hand side of fig. 4 shows that, in order to fulfil the
condition $f_{\rm PBH}\approx 1$, one needs $\mathcal{P}_{\rm
USR}/\mathcal{P}_{\rm CMB}=O(10^{7})$.444There is some difference between peak
theory and threshold statistics in the computation of the abundance, already
present at the Gaussian level (see, e.g., refs. Green  _et al._ (2004); Young
_et al._ (2014); De Luca  _et al._ (2019)). The approach based on peak theory
usually requires slightly smaller values of $\mathcal{P}_{\rm
USR}/\mathcal{P}_{\rm CMB}$ in order to get the same abundance of PBHs, thus
making our findings, based on threshold statistics, even stronger.

We conclude that the condition $f_{\rm PBH}\approx 1$ lies within the region
in which perturbativity is still applicable. This is in contrast with the
conclusion reached in refs. Kristiano and Yokoyama (2022, 2023); Firouzjahi
(2023) in the limit of instantaneous SR/USR/SR transition. The origin of the
difference is the more accurate calculation of the PBHs abundance performed in
our work. In particular, in previous analyses, estimates of $f_{\rm
PBH}\approx 1$ are based on requiring $\mathcal{P}_{\rm USR}\simeq 10^{-2}$,
and on the scaling $\Delta\mathcal{P}=(k_{\rm end}/k_{\rm in})^{2\eta_{\rm
II}}$ in order to capture the growth of the power spectrum at small scales.
However, as explained in section II.2, this scaling does not accurately
describe the amplitude of the power spectrum at its peak, see the left panel
of fig. 3. Ref. Motohashi and Tada (2023) computed one-loop corrections in the
limit of instantaneous SR/USR/SR transitions in scenarios with $\eta_{\rm
II}\leq 3$, finding that the perturbativity bound is relaxed for $\eta_{\rm
II}$ smaller than 3. At a qualitative level, similar results are obtained in
the left panel of fig. 4. In ref. Motohashi and Tada (2023), the
perturbativity bound is traduced in an upper limit on the power-spectrum at
the scale $k_{\rm end}.$ However, as explained above, this procedure
underestimates the maximum amplitude of the power spectrum, see again the left
panel of fig. 3.

####  IV.1.3 Dynamics at the SR/USR/SR transition

We go beyond the instantaneous transition to check if there are cancellations
that affect the order-of-magnitude of the result in eq. (88). There are indeed
compelling reasons to believe that this is the case, as advocated in refs.
Riotto (2023b); Firouzjahi (2023); Firouzjahi and Riotto (2023). The story
goes as follows (the original argument was proposed in ref. Cai  _et al._
(2018) in which the role of non-Gaussianity from non-attractor inflation
models was considered). From the Hubble parameters in eq. (2) and the
background dynamics that follow from the action in eq. (4), it is possible to
calculate the potential and its derivatives exactly. Up to the third order in
field derivatives, we find (see also ref. Leach  _et al._ (2002))

$\displaystyle V(\phi)$ $\displaystyle=H^{2}(3-\epsilon)\,,$ (91)
$\displaystyle V^{\prime}(\phi)$
$\displaystyle=\frac{H^{2}}{\sqrt{2}}\epsilon^{1/2}\left(6-2\epsilon+\epsilon_{2}\right)\,,$
(92) $\displaystyle V^{\prime\prime}(\phi)$
$\displaystyle=H^{2}\left(6\epsilon-2\epsilon^{2}-\frac{3\epsilon_{2}}{2}+\frac{5\epsilon\epsilon_{2}}{2}-\frac{\epsilon_{2}^{2}}{4}-+\frac{\epsilon_{2}\epsilon_{3}}{2}\right)\,,$
(93) $\displaystyle V^{\prime\prime\prime}(\phi)$
$\displaystyle=\frac{H^{2}}{2\sqrt{2\epsilon}}\left[-8\epsilon^{3}+6\epsilon^{2}(4+3\epsilon_{2})-\epsilon\epsilon_{2}(18+6\epsilon_{2}+7\epsilon_{3})+\epsilon_{2}\epsilon_{3}(3+\epsilon_{2}+\epsilon_{3}+\epsilon_{4})\right]$
(94)
$\displaystyle=\frac{1}{2\sqrt{2\epsilon}}\left\\{H^{2}\left[\frac{\ddot{\epsilon_{2}}}{H^{2}}+(3+\epsilon_{2})\frac{\dot{\epsilon_{2}}}{H}\right]+O(\epsilon)\right\\}\,,$
(95)

where in eq. (95) we expanded in the parameter $\epsilon$ and wrote
$\epsilon_{3,4}$ in terms of $\epsilon_{2}$. Consider the flat gauge in which
curvature perturbations are entirely encoded into field fluctuations
$\delta\phi$ by means of the relation
$\zeta=H\delta\phi/\dot{\phi}=-\delta\phi/\sqrt{2\epsilon}$. In this gauge,
the interactions come from Taylor-expanding the quadratic action in field
fluctuations and, at the cubic order, one expects

$\displaystyle\mathcal{L}_{3}\supset\frac{a^{3}}{6}V^{\prime\prime\prime}\delta\phi^{3}=-\frac{a^{3}\epsilon}{3}(\sqrt{2\epsilon}V^{\prime\prime\prime})\zeta^{3}=-\frac{a^{3}\epsilon}{6}\left\\{H^{2}\left[\frac{\ddot{\epsilon_{2}}}{H^{2}}+(3+\epsilon_{2})\frac{\dot{\epsilon_{2}}}{H}\right]+O(\epsilon)\right\\}\zeta^{3}\,.$
(96)

As shown in ref. Chen  _et al._ (2013), the above interaction agrees (modulo a
surface term) with eq. (42) if we integrate by parts

$\displaystyle\int
d^{4}x\,\frac{\epsilon\dot{\epsilon_{2}}}{2}a^{3}\dot{\zeta}\zeta^{2}~{}~{}~{}\to~{}~{}~{}-\int
d^{4}x\,\frac{1}{6}\frac{d}{dt}\left(\epsilon\dot{\epsilon_{2}}a^{3}\right)\zeta^{3}=-\int
d^{4}x\,\frac{a^{3}\epsilon}{6}H^{2}\left[\frac{\ddot{\epsilon_{2}}}{H^{2}}+(3+\epsilon_{2})\frac{\dot{\epsilon_{2}}}{H}\right]\zeta^{3}\,,$
(97)

where in the last step we used the exact identity

$\displaystyle\frac{d}{dt}\left(\epsilon\dot{\epsilon_{2}}a^{3}\right)=a^{3}\epsilon
H^{2}\left[\frac{\ddot{\epsilon_{2}}}{H^{2}}+(3+\epsilon_{2})\frac{\dot{\epsilon_{2}}}{H}\right]\,.$
(98)

The cubic interaction in eq. (97) agrees with eq. (96) up to
$\epsilon$-suppressed terms. Rewriting the interaction as in eq. (96) is quite
instructive. From eq. (98), it seems plausible that drastic variations in time
of $\epsilon_{2}$ could enhance the cubic interaction. However, eq. (96) shows
that these interactions are ultimately controlled by $V^{\prime\prime\prime}$
so that in the case with a smooth SR/USR/SR transition in which
$V^{\prime\prime\prime}$ is expected to be “small”, there must be
cancellations at work within the combination in eq. (98) so that the relevant
coupling in eq. (96) reduces to the term that is SR suppressed. This is the
main argument that was put forth in refs. Riotto (2023b); Firouzjahi (2023);
Firouzjahi and Riotto (2023).

We shall elaborate further on this point. First of all, let us clarify what
“$V^{\prime\prime\prime}$ small” means. We rewrite eq. (95) as follows (we
omit the $O(\epsilon)$ terms and, for clarity’s sake, we write explicitly the
reduced Planck mass)

$\displaystyle\frac{V^{\prime\prime\prime}}{H}=\left(\frac{H}{\bar{M}_{\rm
Pl}}\right)\frac{1}{2\sqrt{2\epsilon}}\left[\frac{\ddot{\epsilon_{2}}}{H^{2}}+(3+\epsilon_{2})\frac{\dot{\epsilon_{2}}}{H}\right]\,.$
(99)

On the left-hand side, the quantity $V^{\prime\prime\prime}/H$ has the
dimension of a coupling. Consequently, imposing the condition
$V^{\prime\prime\prime}/H<1$ corresponds to a weak coupling regime while
$V^{\prime\prime\prime}/H>1$ corresponds to a strongly coupled one. Said
differently, from the perspective of the right-hand side of eq. (99), the
condition $V^{\prime\prime\prime}/H>1$ corresponds to a situation in which the
a-dimensional factor in front of $H/\bar{M}_{\rm Pl}$ becomes so large that it
overcomes the natural suppression given by $H/\bar{M}_{\rm Pl}\ll 1$. In the
left panel of fig. 5, we compute the ratio $V^{\prime\prime\prime}/H$ for two
benchmark SR/USR/SR dynamics with different values of $\delta N$. In the case
in which $\delta N\to 0$ (sharp transition), we observe that
$V^{\prime\prime\prime}/H$ dangerously grows towards the strongly coupled
regime while in the case of a smooth transition it safely takes $O(\ll 1)$
values. As anticipated at the beginning of this section, this argument
confirms that in the case of a smooth transition we expect a reduction in the
size of the trilinear interaction controlled by the factor in eq. (98).

Figure 5: Left panel: Graph of $V^{\prime\prime\prime}/H$ as function of the
background field value $\phi$ for two representative dynamics with,
respectively, $\delta N=0.025$ (dashed black) and $\delta N=0.4$ (solid
black). Starting from the dynamics defined as in section II.2, we compute the
potential by means of the reverse engineering approach described in ref.
Franciolini and Urbano (2022). The values $V_{\star}$ and $H_{\star}$ of,
respectively, the potential and the Hubble rate at CMB scales are chosen in
such a way that both dynamics are consistent with CMB observations (namely,
$V_{\star}\simeq 3\times 10^{-9}$ and $H_{\star}\simeq 3\times 10^{-5}$ with
the reduced Planck mass set to 1). On the right (left) side of the field value
$\phi=\phi_{\rm in}$, $V^{\prime\prime\prime}/H$ is negative (positive). Right
panel: Left-side $y$-axis: time evolution of the curvature modes
$|\bar{\zeta}_{q}(N)|$ for $\bar{q}=2$ in the case $\delta N=0.025$ (dashed
black line) and $\delta N=0.4$ (solid black line). Right-side $y$-axis:
profile of $\eta$ in the case $\delta N=0.025$ (dashed red line) and $\delta
N=0.4$ (solid red line). The region shaded in red highlights the difference
between the sharp and the smooth transition in terms of $\eta$: in the case of
a sharp transition, the curvature mode has more time to grow under the effect
of the negative friction phase implied by the condition $\eta>3/2$.

With this motivation in mind, we go back to the analysis in section IV.1.1 and
we perform the following calculation.

We compute numerically the integral in eq. (79) in order to check the validity
of the scaling in eq. (88) beyond the limit of instantaneous transition. We
define the quantity

$\displaystyle\mathcal{J}_{\delta N}$ $\displaystyle(\eta_{\rm II},\Delta
N_{\rm USR})\equiv\Delta\mathcal{P}_{\rm 1-loop}(k_{*})/\mathcal{P}_{\rm
tree}(k_{*})=$ $\displaystyle 32\int_{N_{\rm in}-\Delta N}^{N_{\rm end}+\Delta
N}dN_{1}\frac{d\eta}{dN}(N_{1})\int_{N_{\rm in}-\Delta
N}^{N_{1}}dN_{2}\bar{\epsilon}(N_{2})\frac{d\eta}{dN}(N_{2})e^{3(N_{2}-N_{\rm
in})}\int_{1}^{e^{\Delta N_{\rm USR}}}\frac{d\bar{q}}{\bar{q}^{4}}{\rm
Im}\left[\bar{\zeta}_{q}(N_{1})^{2}\bar{\zeta}_{q}^{*}(N_{2})\frac{d\bar{\zeta}_{q}^{*}}{dN}(N_{2})\right]\,,$
(100)

that we can directly compare, in the case $\eta_{\rm II}=3$, with $\eta_{\rm
II}^{2}(k_{\rm end}/k_{\rm in})^{6}[1+\log(k_{\rm end}/k_{\rm in})]$ in eq.
(88) using the fact that $k_{\rm end}/k_{\rm in}=e^{\Delta N_{\rm USR}}$.
First, we set $\delta N$ to a very small number, in order to mimic the limit
$\delta N\to 0$, and evaluate $\mathcal{J}_{\delta N}(3,\Delta N_{\rm USR})$
as function of $\Delta N_{\rm USR}$.

Figure 6: Left panel: Comparison between the value of the full integral in eq.
(100) and the analytical estimate in eq. (88). To mimic the instantaneous
transition we take $\delta N=0.025$. Right panel: We plot the ratio
$\mathcal{J}_{\delta N}(3,\Delta N_{\rm USR})$ as function of $\delta N$. In
both figures we take $\eta_{\rm II}=3$.

The comparison is shown in the left panel of fig. 6. We find an excellent
agreement in particular for large $\Delta N_{\rm USR}.$ This is expected,
since the approximation in eq. (87) is more accurate for larger $k_{\rm
end}/k_{\rm in}$. Then, we set $\Delta N_{\rm USR}=3$ and compare the value of
$\mathcal{J}_{\delta N\to 0}(3,3)$ with $\mathcal{J}_{\delta N}(3,3)$ as
function of $\delta N$. We plot the ratio $\mathcal{J}_{\delta
N}(3,3)/\mathcal{J}_{\delta N\to 0}(3,3)$ in the right panel of fig. 6.

Realistic single-field models that feature the presence of a phase of USR
dynamics typically have $\delta N=0.4\div 0.5$ (cf., e.g., ref. Ballesteros
_et al._ (2020); Taoso and Urbano (2021)). This means that, according to our
result in the right panel of fig. 6, we expect that in realistic models the
size of the loop correction gets reduced by one order of magnitude with
respect to what is obtained in the limit of instantaneous SR/USR/SR
transition. This confirms the intuition presented in refs. Riotto (2023b).

It should be noted, however, as evident from our discussion in section II.2,
that in the case of smooth SR/USR/SR transition the amplitude of the power
spectrum gets reduced with respect to the $\delta N\to 0$ limit (cf. the left
panel of fig. 3). The origin of this effect becomes evident if we consider the
right panel of fig. 5. In this figure, we plot the time evolution of the
curvature mode $|\bar{\zeta}_{q}|$ with $\bar{q}=2$ in the two cases of a
sharp and smooth transition (dashed and solid lines, respectively – see
caption for details). In the case of a sharp transition, the curvature mode
experiences a longer USR phase, and its final amplitude is larger with respect
to the case of a smooth transition. As a consequence, therefore, we expect
that the smaller size of the loop correction will be, at least partially,
compensated by the fact that finite $\delta N$ also reduces the amplitude of
the tree-level power spectrum. In order to quantify this information, we
repeat the analysis done in section IV.1.2 but now for finite $\delta N$. We
plot our result in fig. 7. For definiteness, we consider the benchmark value
$\delta N=0.4$ while we keep $\eta_{\rm II}$ and $\Delta N_{\rm USR}$ generic
as in fig. 4.

Figure 7: We consider a generic USR dynamics with varying $\eta_{\rm II}$
($x$-axis) and $\Delta N_{\rm USR}$ ($y$-axis). We take $\eta_{\rm III}=0$ and
the smooth limit $\delta N=0.4$. The region hatched in red corresponds to
$\Delta\mathcal{P}_{\rm 1-loop}(k_{*})>0$. Along the line defined by the
condition $f_{\rm PBH}=1$, we get 100% of DM in the form of asteroid-mass
PBHs. The dotted blue line and the red dashed line correspond, respectively,
to the conditions $f_{\rm PBH}=1$ and $\lim_{\delta N\to
0}\Delta\mathcal{P}_{\rm 1-loop}(k_{*})>0$ as derived in the limit of
instantaneous transition.

Our numerical analysis mirrors the previous intuition. The perturbativity
bound (the region hatched in red corresponds to the condition
$\Delta\mathcal{P}_{\rm 1-loop}(k_{*})>0$) gets weaker because of the partial
cancellation illustrated in the right panel of fig. 6. However, as previously
discussed, the drawback is that taking $\delta N\neq 0$ also reduces the peak
amplitude of the power spectrum. Consequently, the condition $f_{\rm PBH}=1$
requires, for fixed $\eta_{\rm II}$, larger $\Delta N_{\rm USR}$.

As for the limit of instantaneous transition, the condition $f_{\rm PBH}=1$
does not violate the perturbativity bound since the two above-mentioned
effects nearly compensate each other. However, our analysis reveals an
interesting aspect: modelling the SR/USR/SR transition (and, in particular,
the final USR/SR one) beyond the instantaneous limit reduces the impact of the
loop correction but, at the same time, lowers the peak amplitude of the tree-
level power spectrum, which must be compensated by a larger $\Delta N_{\rm
USR}$ see fig. 8. As illustrated in fig. 7, both these effects must be
considered together in order to properly quantify the impact of loop
corrections and the consequent perturbativity bound.

Figure 8: Left: Value of $\Delta N_{\text{\rm USR}}$ required in order to have
$f_{\text{\rm PBH}}=1$ for $\eta_{\text{\rm II}}=3$. Right: Different examples
of evolution of $\eta(N)$ responsible for the USR, assuming various $\delta N$
and fixing $\eta_{\text{\rm II}}=3$. Dashed lines reports the scenario where
$\delta N$ is increased while $\Delta N_{\text{\rm USR}}$ is kept fixed to the
value imposed to have unit PBH abundance in the limit $\delta N\to 0$. Solid
lines report the result when $\Delta N_{\text{\rm USR}}$ is instead adjusted
to keep $f_{\text{\rm PBH}}=1$ fixed. We see that smoother transitions results
in longer USR phases.

This is an interesting point. Refs. Riotto (2023b); Firouzjahi (2023);
Firouzjahi and Riotto (2023) argue that if one goes beyond the limit of
instantaneous transition then the loop correction to the CMB power spectrum
becomes effectively harmless. Technically speaking, in our analysis the role
of the parameter $-6<h<0$ that in Riotto (2023b); Firouzjahi (2023);
Firouzjahi and Riotto (2023) (see also ref. Cai  _et al._ (2018)) controls the
sharpness of the transition is played by our parameter $\delta N$ (with
$h\to-6$ that corresponds to our $\delta N\to 0$ and $h\to 0$ that corresponds
to increasing values of $\delta N$).

In light of our analysis, a very important remark naturally arises: There is a
non-trivial and crucial interplay between the detail of the USR/SR transition
and the amplitude of the tree-level power spectrum that must be properly
included before drawing any conclusion about the relative size of the loop
corrections. On the one hand, it is true that a smooth USR/SR transition
reduces the size of the loop correction; on the other one, the same smoothing
also reduces the amplitude of the power spectrum so that, in order to keep
$f_{\rm PBH}$ fixed, one is forced to either increase the duration of the USR
phase or the magnitude of $\eta$ during the latter. In the end, the two
effects tend to compensate each other if one imposes the condition $f_{\rm
PBH}=1$ (cf. fig. 7).

###  IV.2 Loop evaluation at any scales

We evaluate the loop correction at a generic external momentum $k$, thus
alleviating the assumption $k\ll q$. The dominant modes contributing to the
loop integration remain the ones that cross the horizon during the USR phase
$q\in[k_{\text{in}},k_{\text{end}}]$. As done in the previous section, we are
interested in comparing the one-loop correction with the tree level power
spectrum at the end of inflation, and therefore we perform the late time limit
$\tau\to 0^{-}$. Following the notation introduced in eq. (54), we define

$\mathcal{P}(k)=\lim_{\tau\to
0^{-}}\left(\frac{k^{3}}{2\pi^{2}}\right)\left[\left|\zeta_{k}(\tau)\right|^{2}+\frac{1}{(4\pi)^{2}}\Delta
P(k,\tau)\right]\equiv\mathcal{P}_{\text{tree}}(k)\left(1+\Delta{\cal
P}_{\text{1-loop}}\right)\,,$ (101)

In order to simplify the computation we consider the instantaneous limit
$\delta N\to 0$ of eq. (77). We perform both time integrations keeping the
dominant contribution given by the first Dirac delta in eq. (31). This implies
that we evaluate the integrand function at $\tau_{1}=\tau_{2}=\tau_{\rm end}$.
Notice also that, since the second integration only gets contributions from
half of the Dirac delta domain, we additionally include a factor of $1/2$.
Finally, the jump in $\epsilon_{2}$ leaves a factor $(2\eta_{\text{II}})$ for
each time integration. Therefore, we find

$\displaystyle\Delta P(k,\tau)\equiv-32\eta_{\text{II}}^{2}$
$\displaystyle[\epsilon(\tau_{\text{end}})a^{2}(\tau_{\text{end}})]^{2}\int_{q_{\text{in}}}^{q_{\text{end}}}dq\,q^{2}\,\int_{-1}^{1}d(\cos\theta)\times\Big{\\{}$
$\displaystyle{\rm
Im}\big{[}\zeta_{k}^{*}(\tau)\zeta_{k}^{\prime}(\tau_{\text{end}})\big{]}\times\Big{[}\rm$
$\displaystyle{\rm
Im}\big{[}\zeta_{k}(\tau)\zeta_{k}^{\prime*}(\tau_{\text{end}})|\zeta_{q}(\tau_{\text{end}})|^{2}|\zeta_{k-q}(\tau_{\text{end}})|^{2}\big{]}+$
$\displaystyle{\rm
Im}\big{[}\zeta_{k}(\tau)\zeta_{k}^{*}(\tau_{\text{end}})\big{(}|\zeta_{q}(\tau_{\text{end}})|^{2}\zeta_{k-q}(\tau_{\text{end}})\zeta_{k-q}^{\prime\,*}(\tau_{\text{end}})+|\zeta_{k-q}(\tau_{\text{end}})|^{2}\zeta_{q}(\tau_{\text{end}})\zeta_{q}^{\prime\,*}(\tau_{\text{end}})\big{)}\big{]}\Big{]}\rm+$
$\displaystyle{\rm
Im}[\zeta_{k}^{*}(\tau)\zeta_{k}(\tau_{\text{end}})]\times\Big{[}\rm$
$\displaystyle{\rm
Im}\big{[}\zeta_{k}(\tau)\zeta_{k}^{\prime*}(\tau_{\text{end}})\big{(}|\zeta_{q}(\tau_{\text{end}})|^{2}\zeta_{k-q}^{*}(\tau_{\text{end}})\zeta_{k-q}^{\prime}(\tau_{\text{end}})+|\zeta_{k-q}(\tau_{\text{end}})|^{2}\zeta^{*}_{q}(\tau_{\text{end}})\zeta_{q}^{\prime}(\tau_{\text{end}})\big{)}\big{]}+$
$\displaystyle{\rm
Im}\big{[}\zeta_{k}(\tau)\zeta_{k}^{*}(\tau_{\text{end}})\big{(}|\zeta_{k-q}^{\prime}(\tau_{\text{end}})|^{2}|\zeta_{q}(\tau_{\text{end}})|^{2}+\zeta_{k-q}^{\prime}(\tau_{\text{end}})\zeta_{k-q}^{*}(\tau_{\text{end}})\zeta_{q}(\tau_{\text{end}})\zeta_{q}^{\prime\,*}(\tau_{\text{end}})\big{)}\big{]}+$
$\displaystyle{\rm
Im}\big{[}\zeta_{k}(\tau)\zeta_{k}^{*}(\tau_{\text{end}})\big{(}|\zeta_{k-q}^{\prime}(\tau_{\text{end}})|^{2}|\zeta_{q}(\tau_{\text{end}})|^{2}+\zeta_{k-q}^{\prime\,*}(\tau_{\text{end}})\zeta_{k-q}(\tau_{\text{end}})\zeta^{*}_{q}(\tau_{\text{end}})\zeta_{q}^{\prime}(\tau_{\text{end}})\big{)}\big{]}\Big{]}\Big{\\}}\rm.$
(102)

We have collected the pieces such that each line corresponds to the $i$-th
term in the sum of eq. (77) and $k-q\equiv\sqrt{k^{2}+q^{2}-2kq\cos(\theta)}$
as in the previous section.

In the left panel of fig. 9, we show the resulting 1-loop correction as a
function of the wavenumber $k$ for a representative set of parameters leading
to $f_{\rm BH}\approx 1$: $\eta_{\text{II}}=3$ and $\Delta
N_{\text{USR}}=2.2$. We find values of $\Delta{\cal P}_{\text{1-loop}}$ of the
order of few percent, barring small oscillatory features. A notable exception
is the scale where the tree level power spectrum presents a dip, see fig. 3,
$k_{\text{dip}}/k_{\text{in}}\approx\sqrt{5/4}e^{-3\Delta N_{\text{USR}}/2}$
Byrnes  _et al._ (2019). At that scale the 1-loop correction dominates,
resulting in a spike in $\Delta{\cal P}_{\text{1-loop}}$. As a consequence,
the dip is only realised if the 1-loop correction is neglected, see the right
panel of fig. 9. We also observe that in the limit of small $k\ll
k_{\text{in}}$ the result quickly converges towards the one discussed in the
previous section, as expected. Finally, it is also interesting to notice that
the correction $\Delta{\cal P}_{\text{1-loop}}$ stays almost the same at any
scale, except around $k_{\text{dip}}$. For this reason, we expect that a
generalisation of this calculation to the case for $\delta N\neq 0$ will lead
to results similar to ones presented in the previous section for $k\ll q$.

Figure 9: In both panels, we consider a USR dynamics with $\eta_{\rm II}=3$,
$\Delta N_{\rm USR}=2.2$, $\eta_{\rm III}=0$ and the instantaneous limit
$\delta N=0$. These values corresponds to a scenario producing $f_{\text{\rm
PBH}}\simeq 1$. The vertical gridlines corresponds to $k=k_{\text{\rm in}}$
and $k_{\text{\rm end}}$ in both panels. Left panel: correction to the tree
level power spectrum as a function of $k$ in the limit of $\tau\to 0^{-}$.
Right panel: tree level power spectrum (black) compared to the 1-loop
correction (red line) and their sum (blue dashed line).

At first sight, our result that loop corrections impact the tree-level power
spectrum at the percent level seems at odds with the findings of ref. Inomata
_et al._ (2023) in which it was found that the one-loop power spectrum could
dominate over the tree-level one, thus indicating the breakdown of the
perturbation theory. Upon a closer look, however, there is no contradiction.
Ref. Inomata  _et al._ (2023) considers a particular instance of background
dynamics in which curvature perturbations are resonantly amplified due to a
specific pattern of oscillatory features in the inflaton potential. In such a
model, we checked that the condition $V^{\prime\prime\prime}/H\ll 1$ (cf. eq.
(99) and the related discussion) is not verified and, therefore, it is not
unexpected to find an amplification of loop effects.

It is instructive to consider also a different limit. Since we are assuming
that the USR is followed by a second period of slow roll, characterised by a
negligible $\eta_{\text{III}}$ and a small $\epsilon$, modes in the range
$q\in[k_{\text{in}},k_{\text{end}}]$ freeze around $\tau_{\text{end}}$.
Therefore, the loop correction at $\tau_{\text{end}}$ is very close to its
limit at $\tau\to 0^{-}$, as we verified through a numerical calculation. For
this reason, we set $\tau\to\tau_{\text{end}}$ in eq. (IV.2) and drop the
factors proportional to ${\rm
Im}[\zeta_{k}^{*}(\tau_{\text{end}})\zeta_{k}(\tau_{\text{end}})]$ which
vanish identically. Next, we switch to the barred fields and momenta notation
introduced in sec. IV.1 and simplify the expression using the Wroskian
identity (82). Finally, we arrive at the expression

One-loop correction at generic scales in the instantaneous SR/USR/SR
transition $\displaystyle\lim_{\delta N\to 0}\Delta\mathcal{P}_{\rm
1-loop}(k,\tau_{\text{end}})\approx
4\eta_{\text{II}}^{2}\left(\frac{H^{2}}{8\pi^{2}\epsilon_{\text{ref}}}\right)\frac{1}{|\bar{\zeta}_{k}(N_{\text{end}})|^{2}}\int_{1}^{e^{\Delta
N_{\text{USR}}}}\frac{d\bar{q}}{\bar{q}}\int_{-1}^{+1}d(\cos\theta)\times$
$\displaystyle\left\\{\frac{\bar{k}^{3}}{(\bar{k}-\bar{q})^{3}}|\bar{\zeta}_{q}(N_{\text{end}})|^{2}|\bar{\zeta}_{k-q}(N_{\text{end}})|^{2}+|\bar{\zeta}_{k}(N_{\text{end}})|^{2}\left[|\bar{\zeta}_{q}(N_{\text{end}})|^{2}+\frac{\bar{q}^{3}}{(\bar{k}-\bar{q})^{3}}|\bar{\zeta}_{k-q}(N_{\text{end}})|^{2}\right]\right\\}$
(103)

We stress that this result is exact in the limit of sharp transition $\delta
N\to 0$ and only neglects the contribution from the integration of the step-
function in $N_{\text{in}}$, which is numerically sub-leading.

A number of important comments are in order.

  * i)

In the case in which $k$ is a long CMB mode ($k\ll q$), the first term in the
curly brackets is suppressed by the factor $(k/q)^{3}$. This is nothing but
the number of independent Hubble patches of size $q^{-1}$ in a box of radius
$k^{-1}$. Intuitively, therefore, this contribution represents the situation
in which random small-scale fluctuations lead by chance to a large-scale
fluctuation, and the suppression factor $(k/q)^{3}$ simply indicates that it
is very unlucky for short-mode to be coherent over long scales. The meaning of
this term is very clear, and the above argument is so compelling that it
forces the intuition to believe that there is no way in which CMB modes can be
affected by small-scale ones.

It is worth emphasizing that the computation of the one-loop correction to the
correlation of long-wavelength modes $k$ due to short modes $q$ running in the
loop can be thought of as solving the non-linear evolution equation for the
long mode, cf. ref. Senatore and Zaldarriaga (2010); Riotto (2023a). In the
language of this source method, the first term in the curly brackets of eq.
(103) corresponds to the so-called cut-in-the-middle diagrams, cf. ref.
Pimentel  _et al._ (2012). It also corresponds to the Poisson-suppressed term
identified in ref. Riotto (2023a) while, in the language of ref. Kristiano and
Yokoyama (2023), it corresponds to the correlation of two inhomogeneous
solutions. More in detail, the evolution of the long mode in the presence of
interactions reads

$\displaystyle\hat{\zeta}(\vec{k},N_{\textrm{f}})=$
$\displaystyle\underbrace{\hat{\zeta}(\vec{k},N_{\textrm{in}})+a^{3}(N_{\textrm{in}})\epsilon(N_{\textrm{in}})\frac{d\hat{\zeta}}{dN}(\vec{k},N_{\textrm{in}})\int_{N_{\textrm{in}}}^{N_{\textrm{f}}}dN\frac{1}{a^{3}(N)\epsilon(N)}}_{\textrm{homogeneous
solution (free evolution)}}$ (104) $\displaystyle\underbrace{-\frac{\eta_{\rm
II}}{2}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\hat{\zeta}(\vec{q},N_{\rm
end})\hat{\zeta}(-\vec{q},N_{\rm end})+\frac{\eta_{\rm
II}}{3}\int\frac{d^{3}\vec{q}}{(2\pi)^{3}}\frac{d\hat{\zeta}}{dN}(\vec{q},N_{\rm
end})\hat{\zeta}(-\vec{q},N_{\rm end})}_{\textrm{inhomogeneous solution
(interactions) }}\,,$ (105)

where $N_{\rm f}$ represents some final $e$-fold time after the end of the USR
phase. Consider the terms in the first line, eq. (104). This is the standard
result in the absence of interactions (the homogeneous solution in ref.
Kristiano and Yokoyama (2023)). Eq. (104) tells us that the long mode stays
constant unless the duration of the USR phase is so long to overcome the
smallness of the time derivative of the long mode, which decayed exponentially
fast during the phase preceding the USR. Eq. (104), therefore, gives the
standard tree-level power spectrum if one computes the correlator
$\langle\hat{\zeta}(\vec{k},N_{\textrm{f}})\hat{\zeta}(-\vec{k},N_{\textrm{f}})\rangle$.
The two terms in eq. (105) corresponds to the inhomogeneous solution in ref.
Kristiano and Yokoyama (2023), and encode the effect of the interactions in
the evolution of the long mode. The first term in the curly brackets of our
eq. (103) corresponds to the correlator of two inhomogeneous solutions. The
equivalence between the source method and the “in-in” formalism has been
discussed explicitly in ref. Kristiano and Yokoyama (2023).

  * ii)

Consider now the second term in the curly brackets of eq. (103). Notice that
this term always factorises $|\zeta_{k}|^{2}$ that cancels the denominator in
front of the integral, which is present because of our definition of
$\Delta\mathcal{P}_{\rm 1-loop}$, cf. eq. (101). In the case in which $k$ is a
long CMB mode ($k\ll q$), this term does not pay any $(k/q)^{3}$ suppression.
In the language of the source method, it corresponds to the so-called cut-in-
the-side diagrams, cf. ref. Senatore and Zaldarriaga (2010). These diagrams
represent the evolution of the long mode due to the effect that the long mode
itself has on the expectation value of quadratic operators made of short modes
Senatore and Zaldarriaga (2010) (cf. section V). In ref. Kristiano and
Yokoyama (2023), the second term in the curly brackets of eq. (103) follows
from the correlation between inhomogeneous and homogeneous solutions.

##  V Discussion and outlook

In this work, we discussed the implications of perturbativity in the context
of single-field inflationary models that feature the presence of a transient
phase of USR. More in detail, we defined the perturbativity condition

$\displaystyle\mathcal{P}(k)\equiv\mathcal{P}_{\text{tree}}(k)\left[1+\Delta{\cal
P}_{\text{1-loop}}(k)\right]~{}~{}~{}\Longrightarrow~{}~{}~{}\Delta{\cal
P}_{\text{1-loop}}(k)\overset{!}{<}1\,,$ (106)

in which the one-loop correction is integrated over the short modes that are
enhanced by the USR dynamics. We explored the consequences of eq. (106) at any
scale $k$ even though the main motivation for our analysis was the recent
claim of ref. Kristiano and Yokoyama (2022) according to which the relative
size of the loop correction at scales relevant for CMB observations (that is,
$k=O(k_{\star})$ with $k_{\star}=0.05$ Mpc-1) threatens the validity of
perturbativity at the point of ruling out the idea of PBH formation via USR
dynamics in single-field inflation.

In this section, we summarize the main results and limitations of our analysis
and we will discuss future prospects.

  * $\circ$

In the limit of instantaneous SR/USR/SR transition, we confirm the computation
of the 1-loop corrections of ref. Kristiano and Yokoyama (2022). However, we
provide a more detailed and precise discussion of the theoretical bound that
can be obtained by imposing the perturbativity condition in eq. (106) on the
power spectrum of curvature perturbations at CMB scales.

As far as this part of the analysis is concerned, the key difference with
respect to ref. Kristiano and Yokoyama (2022) is that we compare the size of
the loop correction with an accurate computation of the PBH abundance rather
that with the order-of-magnitude estimate of the enhancement of the power
spectrum, based on the SR formula, used in ref. Kristiano and Yokoyama (2022).
Our approach, therefore, includes the following effects. i) First of all, we
generalize the USR dynamics for generic values $\eta_{\rm II}\neq 3$ (cf.
section II.2 for our parametrization of the background dynamics); ii) the
enhancement of the power spectrum at scales relevant for PBH formation is
accurately computed by numerically solving the M-S equation rather than using
the SR scaling; iii) by computing the PBH abundance $f_{\rm PBH}$, we
automatically account for the fact that the correct variable that describes
PBH formation in the standard scenario of gravitational collapse is the
smoothed density contrast rather than the curvature perturbation field, and we
include in our computation the full non-linear relation between the two.

As for this part of the analysis, our findings are summarized in fig. 4. We
find that loop corrections remain of the order of few percent and therefore it
is not possible to make the bold claim that PBH formation from single-field
inflation is ruled out – not even in the limit of instantaneous SR/USR/SR
transition.

  * $\circ$

We extend the analysis of ref. Kristiano and Yokoyama (2022) by considering a
more realistic USR dynamics. In particular, we implement a smooth description
of the SR/USR/SR transition. Recently, refs. Riotto (2023b); Firouzjahi and
Riotto (2023); Firouzjahi (2023) claimed that the presence of a smooth
transition in the final USR/SR transition makes the loop correction
effectively harmless. Our analysis shows that this conclusion could be
invalidated by the fact that there is an interplay between the size of the
loop correction and the amplitude of the tree-level power spectrum that is
needed to generate a sizable abundance of PBHs. On the one hand, it is true
that a smooth USR/SR transition reduces the size of the loop correction; on
the other one, the same smoothing also reduces the amplitude of the tree-level
power spectrum so that, in order to keep $f_{\rm PBH}$ fixed, one is forced to
either increase the duration of the USR phase or the magnitude of $\eta$
during the latter. In the end, the two effects tend to compensate each other.
As for this part of the analysis, our findings are summarized in fig. 7.

  * $\circ$

We consider the one-loop correction of short modes to the tree-level power
spectrum at any scale. We find that perturbativity is always satisfied in
models that account for the condition $f_{\rm PBH}=1$.

More quantitatively, we find that the relative size of the loop correction
with respect to the tree-level value of the power spectrum does not exceed the
level of a few percent. As for this part of the analysis, our findings are
summarized in fig. 9. We point out one notable exception of phenomenological
relevance. A generic feature of the USR dynamics is that it produces a
characteristic dip in the tree-level power spectrum, as the one observed in
the left panel of fig. 3. The phenomenological consequences of such a putative
dip range from CMB $\mu$-space distortions Özsoy and Tasinato (2022) to 21-cm
signals Balaji  _et al._ (2022). Our analysis shows that the existence of the
dip is nothing but an artifact of the tree-level computation, and disappears
after including loop corrections. This is because, due to the smallness of the
tree-level power spectrum around the characteristic wavenumbers of the dip,
the non-vanishing loop correction gives the dominant contribution. This is
illustrated in the right panel of fig. 9.

At the conceptual level, it remains true that, in the presence of USR
dynamics, loop corrections of short modes may sizably affect the power
spectrum at CMB scales. This result echoes an issue of naturalness – an
infrared quantity (the amplitude of the curvature power spectrum at CMB
scales) appears to be sensitive, via loop effects, to physics that takes place
at much shorter scales (those related to PBH formation) – and clashes with the
intuition that physics at such vastly different scales should decouple.

The coupling between short and long modes gives a physical effect for the
following reason. As discussed in section IV.2, the relevant loop correction
to the power spectrum at CMB scales comes from the correlation between
homogeneous and inhomogeneous solutions. This is most easily seen within the
source method in which one considers the correlation between a freely evolving
long mode and a second long mode which evolves in the presence of
interactions, cf. eq. (105). Borrowing from ref. Pimentel  _et al._ (2012)
(see also ref. Riotto (2023a)), we write the formal solution of the non-linear
evolution equations for a long wavelength mode $\zeta_{\rm L}$ as $\zeta_{\rm
L}=\hat{O}^{-1}[S[\zeta_{\rm S},\zeta_{\rm S},\zeta_{\rm L}]]$, where $S$
represents a generic sum of operators that are quadratic in the short
wavelength mode $\zeta_{\rm S}$ and that can also depend on $\zeta_{\rm L}$ if
one considers the short modes in the background perturbed by the long mode.
More concretely, in our case such a solution is the one given by eq. (105).
The one-loop power spectrum is given by

$\displaystyle\langle\zeta_{\rm L}\zeta_{\rm
L}\rangle\sim\langle\hat{O}^{-1}[S[\zeta_{\rm S},\zeta_{\rm S},\zeta_{\rm
L}=0]]\,\hat{O}^{-1}[S[\zeta_{\rm S},\zeta_{\rm S},\zeta_{\rm
L}=0]]\rangle+\langle\hat{O}^{-1}[S[\zeta_{\rm S},\zeta_{\rm S},\zeta_{\rm
L}]]\,\zeta_{\rm L}\rangle\,.$ (107)

The first term represents the effect of the short-scale modes in their
unperturbed state (that is, with $\zeta_{\rm L}=0$) directly on the power
spectrum of the long wavelength mode. This is our first term in eq. (103). As
discussed in section IV.2, this term does not alter the long-wavelength
correlation since it is very improbable that random short-scale fluctuations
coherently add up to induce a long-wavelength correlation. The second term in
eq. (107), on the contrary, correlates a freely evolving long mode $\zeta_{\rm
L}$ with the effect that the long mode itself has on the expectation value of
quadratic operators made of short modes. Let us explain this point, which is
crucial. Consider the schematic in fig. 10.

Figure 10: Left: Expansion in time of the unperturbed universe (time passes by
along the $y$-axis); the universe expands by the same amount at every point.
Right: Expansion in time of the perturbed universe. The long mode ($\zeta_{\rm
L}$, blue) acts as a local rescaling of the scale factor, and short scales are
modulated accordingly. More specifically, if we consider the black dots we see
that they experience a different amount of expansion depending on the value of
$\zeta_{\rm L}$.

The key point is the following. In the comoving gauge, the short modes evolve
in the background that is perturbed by the long mode. In the limit in which
the long mode $\zeta_{\rm L}$ has a wavelength much longer than the horizon,
it simply acts as a rescaling of the coordinates since it enters as a local
change of the scale factor. This is schematically illustrated in fig. 10. This
figure shows intuitively that the short scales are modulated by the presence
of the long mode. The presence of the long mode acts as a rescaling of the
coordinates and we can absorb it by rescaling the short-scale momenta
$q\to\tilde{q}=e^{\zeta_{\rm L}}q$ Riotto (2023a). If the power spectrum of
the short modes is scale-invariant, this rescaling does nothing. However, if
the power spectrum of the short modes breaks scale invariance, we
schematically have in the loop integral over the short modes, expanding at the
first order

$\displaystyle\int\frac{dq}{q}\mathcal{P}(q)\overset{q\to\tilde{q}=e^{\zeta_{\rm
L}}q}{~{}~{}~{}~{}\Longrightarrow~{}~{}~{}~{}}\int\frac{d\tilde{q}}{\tilde{q}}\mathcal{P}(\tilde{q})=\int\frac{dq}{q}\mathcal{P}(e^{\zeta_{\rm
L}}q)=\int\frac{dq}{q}\left[\mathcal{P}(q)+\zeta_{\rm
L}\frac{d\mathcal{P}}{dq}q\right]=\int\frac{dq}{q}\left[\mathcal{P}(q)+\zeta_{\rm
L}\,\mathcal{P}(q)\,\frac{d\log\mathcal{P}}{d\log q}\right]\,,$ (108)

so that the presence of the long mode affects the correlation of short modes
when their power spectrum is not scale invariant. The second term in the above
equation describes precisely the effect put forth before: the presence of the
long mode alters the expectation value of quadratic operators made of short
modes, in this case the short-mode two-point function. Back to eq. (107), one
expects the one-loop correction Riotto (2023a)

$\displaystyle\Delta\mathcal{P}_{\rm
1-loop}(k)\sim\mathcal{P}(k)\int\frac{dq}{q}\,\mathcal{P}(q)\,\frac{d\log\mathcal{P}}{d\log
q}\,.$ (109)

The above discussion shows that the one-loop corrections on long modes do not
decouple when the power spectrum of the short modes is not scale invariant.
This explains why our correction vanishes in the limit $\eta_{\rm II}=0$ in
which indeed the power spectrum does become scale invariant. The breaking of
scale invariance is the hallmark of the USR dynamics and, more importantly, a
necessary feature in all models of single-field inflation that generate an
order-one abundance of PBHs (cf. the intuitive schematic in fig. 11).

Figure 11: Illustrative schematic of the correction induced on the two-point
correlator of long modes by a loop of short modes. On the right side, we plot