# Perturbative region on non-Gaussian parameter space in single-field
inflation

Jason Kristiano  and Jun’ichi Yokoyama

###### Abstract

We calculate one-loop correction to the two-point functions of curvature
perturbation in single-field inflation generated by cubic self-interaction.
Incorporating the observed red-tilted spectrum of curvature perturbation, the
relevant one-loop correction takes a finite value and inversely proportional
to the spectral tilt. Requiring one-loop correction to be much smaller than
the tree-level contribution leads to an upper bound on primordial non-
Gaussianity. While observationally allowed region of non-Gaussian parameter
space is found to be entirely included by the region, where one-loop
correction is smaller than the tree-level contribution, an appreciably large
region has one-loop correction larger than 1% or even 10% of the latter. If
future observations conclude non-Gaussianity falls in such a region, then it
would be important to incorporate higher-order corrections to the spectrum in
order to achieve precise cosmology. In some extreme cases, where one-loop
correction has a comparable magnitude to the tree-level contribution, it might
indicate breakdown of the cosmological perturbation theory in the context of
single-field inflation.

##  1 Introduction

Observation of cosmic microwave background (CMB) shows that its temperature is
very isotropic with small fluctuations of order of $10^{-5}$ times the
background temperature [1, 2]. The homogeneity of CMB can be explained if our
Universe underwent a very rapid accelerated expansion before the Big Bang,
known as inflation [3, 4, 5]. At the same time, small fluctuations are
generated from quantum fluctuations of spacetime and stretched by inflation
[6, 7, 8, 9, 10]. This framework is called cosmological perturbation theory.
Observation also shows the power spectrum of fluctuations is almost scale-
invariant with a highly Gaussian distribution. Deviation from scale-invariant
power spectrum is called spectral-tilt and has been observed fairly well.
Deviation from Gaussian distribution is called primordial non-Gaussianity, and
so far only upper-bound is obtained observationally [11].

Theoretically, such quantum fluctuations are analyzed within the framework of
quantum field theory (QFT). Power spectrum corresponds to the vacuum
expectation value (VEV) of two-point functions of the fluctuations. At the
lowest order, primordial non-Gaussianity corresponds to the VEV of three-point
functions of fluctuations. This originates from cubic self-interactions of the
fluctuations. At the same time, such self-interaction can generate
backreaction to the two-point functions. In the terminology of QFT, the
backreaction is called loop correction.

We are interested in models predicting a sizable amount of non-Gaussianity
close to the observational bound [11], as such models may induce large loop
corrections of inflationary power spectrum. The simplest inflation model,
namely canonical slow-roll inflation, predicts a small non-Gaussianity [12].
Its bispectrum is proportional to the slow-roll parameters. Therefore, we have
to look for other inflation models. There are some inflation models which can
realize a large non-Gaussianity by modifying kinetic terms, such as $k$\- and
$G$-inflation [13, 14, 15, 16], ghost condensate [17, 18], and Dirac-Born-
Infeld inflation [19]. For inflation models that generate large loop
correction, we must require that the one-loop correction is small enough
compared to the tree-level two-point functions, otherwise we cannot trust the
perturbation theory because the power spectrum would be infinitely large if we
calculated higher-order loop correction.

Weinberg first realized the case that loop correction can become significantly
large [20, 21]. He calculated one-loop correction generated by coupling
between curvature perturbations and massless scalar fields. For perturbation
with comoving wavenumber $k$, he found the one-loop correction proportional to
$\log k$. Another notable calculation was performed by Senatore and
Zaldarriaga [22], where they calculated one-loop correction generated by a
cubic self-interaction of curvature perturbations. Assuming a scale-invariant
curvature perturbation, they performed the calculation by two different
regularization methods, physical cutoff and dimensional regularization. They
found that both methods lead to the same one-loop correction results
proportional to $\log H$, where $H$ is Hubble scale during inflation. After
that, full calculation of one-loop correction generated by all possible self-
interaction of curvature perturbations was done by Bartolo, et. al. [23],
where they relied on the regularization method by Senatore and Zaldarriaga.
Another way to obtain one-loop correction is by bootstrap method as
demonstrated by [24], relying on the cosmological optical theorem [25].

In our previous letter [26], we calculated one-loop correction generated by a
cubic self-interaction of curvature perturbations. By properly incorporating
the observed red tilt of the spectrum of curvature perturbation, we found that
the cosmologically relevant one-loop correction has a finite contribution
inversely proportional to the spectral-tilt. Note that Senatore and
Zaldarriaga [22] assumed the scale-invariant curvature perturbation, and
consequently they found a divergent one-loop correction. They performed
regularization and renormalization to remove the divergence and extract the
finite part of it. In our approach, we do not perform regularization and
renormalization because the relevant one-loop correction is finite. Then, we
obtained an upper bound on primordial non-Gaussianity by requiring one-loop
correction to be much smaller than the tree-level two-point functions.

In this paper, we perform the full calculation of one-loop correction
generated by cubic self-interaction of curvature perturbations from the
effective field theory of inflation [27]. In section 2, we review the standard
cosmological perturbation theory. In section 3, we calculate the one-loop
correction to the power spectrum. In section 4, we plot the perturbative
region in the non-Gaussian parameter space and compare it with observational
constraint. Finally in section 5, we conclude our paper. We show the detailed
calculation in the appendix.

##  2 Cosmological Perturbation Theory

In this section, we briefly review the standard cosmological perturbation
theory. We consider the simplest single-field inflation model that can
generate a large non-Gaussianity, namely $k$ inflation or $P(X,\phi)$
inflation. It is an inflation model whose Lagrangian is a general function of
$\phi$ and $X\equiv-\nabla_{\mu}\phi\nabla^{\mu}\phi/2$, with the action

$S=\frac{1}{2}\int\mathrm{d}^{4}x\sqrt{-g}\left[M_{\mathrm{pl}}^{2}R+2P(X,\phi)\right],$
(2.1)

where $g=\mathrm{det}(g_{\mu\nu})$, $g_{\mu\nu}$ is a metric tensor, and $R$
and $M_{\mathrm{pl}}$ are the Ricci scalar and the reduced Planck scale,
respectively. The Einstein equations in the homogeneous and isotropic
background with the metric
$\mathrm{d}s^{2}=-\mathrm{d}t^{2}+a^{2}(t)\mathrm{d}\mathbf{x}^{2}$ read

$H^{2}=\frac{1}{3}\left(2\bar{X}\bar{P}_{,X}-\bar{P}\right)~{}~{}\mathrm{and}~{}~{}\dot{H}=-\bar{X}\bar{P}_{,X},$
(2.2)

where $H=\dot{a}/a$, a dot denotes time derivative, a comma represents a
partial derivative, and $\bar{X}=\dot{\bar{\phi}}^{2}/2$ is the kinetic energy
of the background scalar field with $\bar{P}=P(\bar{X},\bar{\phi})$.

###  2.1 Second-order Action

We now incorporate small perturbation to the homogeneous and isotropic
background, expressing the scalar field as
$\phi(\mathbf{x},t)=\bar{\phi}(t)+\delta\phi(\mathbf{x},t)$. It is convenient
to use the Arnowitt-Deser-Misner (ADM) formalism where the perturbed metric is
expressed as

$ds^{2}=-N^{2}\mathrm{d}t^{2}+\gamma_{ij}(\mathrm{d}x^{i}+N^{i}\mathrm{d}t)(\mathrm{d}x^{j}+N^{j}\mathrm{d}t).$
(2.3)

Here $\gamma_{ij}$ is a three-dimensional metric on slices of constant $t$,
$N$ is a lapse function, and $N^{i}$ is a shift vector. Then, the action
becomes

$S=\frac{1}{2}\int\mathrm{d}^{4}x\sqrt{\gamma}\left[N{}^{(3)}R+2NP(X,\phi)+N^{-1}(E_{ij}E^{ij}-E^{2})\right],$
(2.4)

where $\gamma=\mathrm{det}~{}\gamma_{ij}$, ${}^{(3)}R$ is Ricci scalar of
$\gamma_{ij}$,

$E_{ij}={\textstyle\frac{1}{2}}(\dot{\gamma}_{ij}-\nabla_{i}N_{j}-\nabla_{j}N_{i}),~{}~{}\mathrm{and}~{}~{}E=E^{i}_{i}.$
(2.5)

$\nabla_{i}$ is three-dimensional covariant derivative of the metric
$\gamma_{ij}$. However, such perturbations are very general and contain gauge
freedom. We can choose to work in the comoving gauge

$\delta\phi(\mathbf{x},t)=0,~{}\gamma_{ij}(\mathbf{x},t)=a^{2}(t)[1+2\zeta(\mathbf{x},t)]\delta_{ij},$
(2.6)

where $\zeta(\mathbf{x},t)$ is the comoving curvature perturbation. Then, $N$
and $N^{i}$ are obtained by solving the Hamiltonian and momentum constraint
equations

$\displaystyle\nabla_{i}[N^{-1}(E^{i}_{j}-\delta^{i}_{j}E)]=0,$
$\displaystyle{}^{(3)}R+2P-4XP_{,X}-N^{-2}(E_{ij}E^{ij}-E^{2})=0,$ (2.7)

respectively. $N_{i}$ is decomposed into two parts,
$N_{i}=\tilde{N}_{i}+\partial_{i}\psi$ where $\partial_{i}\tilde{N}^{i}=0$.
Then, $N$ and $N^{i}$ are expanded in terms of powers of $\zeta$ as

$\displaystyle N=1+\alpha_{1}+\alpha_{2}+\dots$
$\displaystyle\tilde{N}_{i}=\tilde{N}_{i}^{(1)}+\tilde{N}_{i}^{(2)}+\dots$
$\displaystyle\psi=\psi_{1}+\psi_{2}+\dots$ (2.8)

where $\alpha_{n},\tilde{N}^{(n)},\psi_{n}\sim O(\zeta^{n})$. The first-order
solution of the constraint equations are

$\alpha_{1}=\frac{\dot{\zeta}}{H},~{}~{}\tilde{N}_{i}^{(1)}=0,~{}~{}\mathrm{and}~{}~{}\psi_{1}=-\frac{\zeta}{H}+a^{2}\frac{\epsilon}{c_{s}^{2}}\partial^{-2}\dot{\zeta},$
(2.9)

where inverse Laplacian $\partial^{-2}$ is defined as
$\partial^{-2}(\partial^{2}\zeta)=\zeta$ and $c_{s}$ is the sound speed given
by

$c_{s}^{2}=\frac{\bar{P}_{,X}}{\bar{P}_{,X}+2\bar{X}\bar{P}_{,XX}}.$ (2.10)

Here $\epsilon\equiv-\dot{H}/{H^{2}}$ is a slow-roll parameter of the Hubble
parameter $H$, which has a very weak time dependence during inflation. Later,
we will also use another slow-roll parameter $\eta=\dot{\epsilon}/\epsilon H$,
which represents the change of slow-roll parameter $\epsilon$ during
inflation. After substituting the first-order solution for $N$ and $N^{i}$ to
the action and performing cumbersome partial integration, we obtain the
following second-order action of the curvature perturbation

$S^{(2)}=M_{\mathrm{pl}}^{2}\int\mathrm{d}t~{}\mathrm{d}^{3}x~{}a^{3}\frac{\epsilon}{c_{s}^{2}}\left[\dot{\zeta}^{2}-\frac{c_{s}^{2}}{a^{2}}(\partial_{i}\zeta)^{2}\right].$
(2.11)

By introducing the Mukhanov-Sasaki (MS) variable [7, 28]

$v\equiv
M_{\mathrm{pl}}z\zeta,~{}~{}~{}~{}z\equiv\frac{a}{c_{s}}\sqrt{2\epsilon},$
(2.12)

the second-order action (2.11) becomes canonically normalized

$S^{(2)}=\frac{1}{2}\int\mathrm{d}\tau~{}\mathrm{d}^{3}x\left[(v^{\prime})^{2}-c_{s}^{2}(\partial_{i}v)^{2}+\frac{z^{\prime\prime}}{z}v^{2}\right],$
(2.13)

where $\mathrm{d}\tau=a~{}\mathrm{d}t$, $\tau$ is conformal time, and the
prime denotes derivative with respect to $\tau$. The classical equation of
motion of $v(\mathbf{x},\tau)$ is

$v^{\prime\prime}(\mathbf{x},\tau)+\left(-c_{s}^{2}\partial_{i}\partial^{i}+\frac{z^{\prime\prime}}{z}\right)v(\mathbf{x},\tau)=0.$
(2.14)

Performing Fourier transformation

$v(\mathbf{x},\tau)=\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}}v(\mathbf{p},\tau)e^{-i\mathbf{p}\cdot\mathbf{x}},$
(2.15)

the equation of motion in momentum space is

$v^{\prime\prime}(\mathbf{p},\tau)+\left(c_{s}^{2}p^{2}-\frac{z^{\prime\prime}}{z}\right)v(\mathbf{p},\tau)=0,$
(2.16)

where $p$ is defined as $p=|\mathbf{p}|$. Such equation is called MS equation.
To solve it, we need the time dependence of $z$, which up to the lowest order
in the slow-roll parameters is given by

$\frac{z^{\prime\prime}}{z}=\frac{1}{\tau^{2}}\left(2+3\epsilon+\frac{3\eta}{2}\right).$
(2.17)

Quantization is performed by promoting the MS variable to an operator

$\hat{v}(\mathbf{p},\tau)=M_{\mathrm{pl}}z\zeta(\mathbf{p},\tau)=v_{p}(\tau)\hat{a}_{\mathbf{p}}+v^{*}_{p}(\tau)\hat{a}_{-\mathbf{p}}^{\dagger},$

where $v_{p}(\tau)$ is the mode function, superscript $*$ denotes complex
conjugate operation, $\hat{a}_{-\mathbf{p}}^{\dagger}$ and
$\hat{a}_{\mathbf{p}}$ are creation and annihilation operators with the
commutation relation

$\left[\hat{a}_{\mathbf{p}},\hat{a}_{-\mathbf{q}}^{\dagger}\right]=(2\pi)^{3}\delta^{3}(\mathbf{p}+\mathbf{q}).$
(2.18)

The mode function is the solution of MS equation

$v_{p}(\tau)=\left(-\frac{\pi\tau}{4}\right)^{1/2}\exp\left[\frac{i\pi}{4}(2\nu+1)\right]H_{\nu}^{(1)}(-pc_{s}\tau),~{}~{}\nu=\frac{3}{2}+\frac{\eta}{2}+\epsilon,$
(2.19)

which corresponds to the Bunch-Davies vacuum $\ket{0}$ at early time defined
by $\hat{a}_{\mathbf{p}}\ket{0}=0$. In general, solution of the MS equation is
a linear combination of Hankel function of the first kind
$H_{\nu}^{(1)}(-pc_{s}\tau)$ and second kind $H_{\nu}^{(2)}(-pc_{s}\tau)$.
Because we are looking for solution which goes as Bunch-Davies mode function
$\exp(-ipc_{s}\tau)$ at early times $\tau\rightarrow-\infty$, the coefficient
of $H_{\nu}^{(2)}(-pc_{s}\tau)$ vanishes.

Observable quantities of cosmological perturbations are their correlation
function evaluated at the end of inflation. The two-point function of
curvature perturbation and power spectrum at a late time during inflation
$\tau_{0}$ are defined as

$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(\mathbf{q})\right\rangle=(2\pi)^{3}\delta^{3}(\mathbf{p}+\mathbf{q})\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle,$
(2.20)
$\displaystyle\Delta^{2}_{s}(p)\equiv\frac{p^{3}}{2\pi^{2}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle,$
(2.21)

and take a limit $\tau_{0}\rightarrow 0$, when all the relevant modes are in
the superhorizon regime and the observed power spectrum is evaluated. Here the
bracket denotes the VEV, $\langle\cdots\rangle=\bra{0}\cdots\ket{0}$, and
$\Delta^{2}_{s}(p)$ is the power spectrum multiplied by the phase space
density.

To obtain the late time limit of curvature perturbation, we need to know the
explicit time dependence of $z$. It is given by

$\frac{z}{z_{H}}=\left(\frac{\tau}{\tau_{H}}\right)^{-\left(1+\epsilon+\frac{\eta}{2}\right)},$
(2.22)

where $\tau_{H}$ is an arbitrary time reference and $z_{H}=z(\tau_{H})$.
Choosing $\tau_{H}$ as sound horizon crossing time $c_{s}k=aH$ so

$\displaystyle\tau_{H}$
$\displaystyle=-\left[\frac{1}{aH(1-\epsilon)}\right]_{H}=-\frac{1}{c_{s}k(1-\epsilon)_{H}},$
(2.23) $\displaystyle z_{H}$
$\displaystyle=\left(\frac{a}{c_{s}}\sqrt{2\epsilon}\right)_{H}=\left(\frac{k}{H}\sqrt{2\epsilon}\right)_{H},$
(2.24)

then $z$ as a function of $\tau$ is

$z=\left(\frac{k}{H}\sqrt{2\epsilon}\right)_{H}\left[-c_{s}k\tau(1-\epsilon)_{H}\right]^{-\left(1+\epsilon+\frac{\eta}{2}\right)}.$
(2.25)

Hence, up to the lowest order in the slow-roll parameters, the mode function
is

$\zeta_{p}(\tau)=\frac{v_{p}(\tau)}{zM_{\mathrm{pl}}}=\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}^{\frac{1}{2}}\frac{e^{-ic_{s}p\tau}}{p^{3/2}}(1+ic_{s}p\tau),$
(2.26)

and the power spectrum is

$\Delta^{2}_{s(0)}(p)=\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}c_{s}\epsilon}\right)_{H}=\Delta^{2}_{s(0)}(p_{*})\left(\frac{p}{p_{*}}\right)^{n_{s}-1}.$
(2.27)

Here subscripts $(0)$ denotes tree-level contribution and $p_{*}$ is an
arbitrary pivot momentum. The power spectrum is almost scale invariant with
the deviation parametrized by the spectral index

$n_{s}-1=\frac{\mathrm{d}\log\Delta^{2}_{s(0)}}{\mathrm{d}\log
p}=-2\epsilon-\eta$ (2.28)

Indeed, there is a weak momentum dependence because of time dependence of the
quantities in the parenthesis at the horizon crossing of each mode.

###  2.2 Third-order Action

Expanding more until third-order, in the limit of $c_{s}\ll 1$, we get [29,
30, 31]

$S_{\mathrm{int}}=\int\mathrm{d}t~{}\mathrm{d}^{3}x\left[-\frac{2\lambda}{H^{3}}a^{3}\dot{\zeta}^{3}+\frac{\epsilon
M_{\mathrm{pl}}^{2}}{Hc_{s}^{2}}a\dot{\zeta}(\partial_{i}\zeta)^{2}\right],$
(2.29)

where $\lambda=X^{2}P_{,XX}+(2/3)X^{3}P_{,XXX}$. Such third-order action can
be regarded as cubic self-interaction term of the scalar perturbation. The
interaction terms are also predicted by the effective field theory of
inflation [27]. To compare with observational constraint [11], we have to
reparametrize $\lambda$ to $\tilde{c}_{3}$ as

$\lambda=-\frac{\epsilon M_{\mathrm{pl}}^{2}H^{2}}{3c_{s}^{4}}\tilde{c}_{3}.$
(2.30)

Therefore, the third-order action becomes

$S_{\mathrm{int}}=\int\mathrm{d}\tau~{}\mathrm{d}^{3}x~{}\frac{\epsilon
M_{\mathrm{pl}}^{2}}{Hc_{s}^{2}}\left[\frac{2\tilde{c}_{3}}{3c_{s}^{2}}a^{4}\dot{\zeta}^{3}+a^{2}\dot{\zeta}(\partial_{i}\zeta)^{2}\right],$
(2.31)

where we define Lagrangian density $\mathcal{L}_{\mathrm{int}}$ as the
integrand of (2.31). Although $\mathcal{L}_{\mathrm{int}}$ contain time
derivative of the field, the Hamiltonian density is still given by
$\mathcal{H}_{\mathrm{int}}=-\mathcal{L}_{\mathrm{int}}$ as shown by [32, 23].
The interaction Hamiltonian reads

$H_{\mathrm{int}}(\tau)=-\int\mathrm{d}^{3}x~{}\frac{\epsilon
M_{\mathrm{pl}}^{2}}{Hc_{s}^{2}}\left[\frac{2\tilde{c}_{3}}{3c_{s}^{2}}a^{4}\dot{\zeta}^{3}+a^{2}\dot{\zeta}(\partial_{i}\zeta)^{2}\right].$
(2.32)

Note that the dot still denotes derivative with respect to $t$. We define
$H_{\mathrm{int}}$ as function of $\tau$ for use in the following section.

##  3 One-Loop Correction

We are interested in the statistic of fluctuations in the CMB. To calculate
higher-order correction, we use the in-in perturbation theory which gives the
expectation value of an operator $\mathcal{O}(\tau)$ as

$\langle\mathcal{O(\tau)}\rangle=\left\langle\left[\bar{\mathrm{T}}\exp\left(i\int_{-\infty}^{\tau}\mathrm{d}\tau^{\prime}H_{\mathrm{int}}(\tau^{\prime})\right)\right]\mathcal{\hat{O}}(\tau)\left[\mathrm{T}\exp\left(-i\int_{-\infty}^{\tau}\mathrm{d\tau^{\prime}}H_{\mathrm{int}}(\tau^{\prime})\right)\right]\right\rangle,$
(3.1)

where $\mathrm{T}$ and $\bar{\mathrm{T}}$ denote time and antitime ordering.
Because we are interested in the higher-order correction of two-point
correlation, the relevant operator is $\zeta(\mathbf{p})\zeta(-\mathbf{p})$
evaluated at $\tau=\tau_{0}~{}(\rightarrow 0)$. Its first-order expansion
vanishes, yielding an odd-point correlation function. Second-order expansion
of the perturbation theory reads

$\displaystyle\langle\mathcal{O(\tau)}\rangle=\langle\mathcal{O(\tau)}\rangle_{(0,2)}^{\dagger}+\langle\mathcal{O(\tau)}\rangle_{(1,1)}+\langle\mathcal{O(\tau)}\rangle_{(0,2)},$
$\displaystyle\langle\mathcal{O(\tau)}\rangle_{(1,1)}=\int_{-\infty}^{\tau}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau}\mathrm{d}\tau_{2}\left\langle
H_{\mathrm{int}}(\tau_{1})\mathcal{\hat{O}}(\tau)H_{\mathrm{int}}(\tau_{2})\right\rangle,$
(3.2)
$\displaystyle\langle\mathcal{O(\tau)}\rangle_{(0,2)}=-\int_{-\infty}^{\tau}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\left\langle\mathcal{\hat{O}}(\tau)H_{\mathrm{int}}(\tau_{1})H_{\mathrm{int}}(\tau_{2})\right\rangle.$
(3.3)

They are proportional to the product of two interaction Hamiltonian. Because
each interaction Hamiltonian has two interaction terms, namely
$\dot{\zeta}^{3}$ and $\dot{\zeta}(\partial\zeta)^{2}$, multiplication of them
leads to four terms. If we label $\dot{\zeta}^{3}$ and
$\dot{\zeta}(\partial\zeta)^{2}$ as $a$ and $b$, respectively, the four terms
are $[aa]$, $[bb]$, $[ab]$, and $[ba]$.

We start our calculation from the $(1,1)$ contribution of $[aa]$ term that is
given by

$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle_{(1,1)}^{[aa]}=$
$\displaystyle
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{4}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{0}\mathrm{d}\tau_{2}\frac{a^{4}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\prod_{a=1}^{6}\left[\frac{\mathrm{d}^{3}k_{a}}{(2\pi)^{3}}\right]\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})~{}\delta^{3}(\mathbf{k}_{4}+\mathbf{k}_{5}+\mathbf{k}_{6})$
$\displaystyle\left\langle\dot{\zeta}(\mathbf{k}_{1},\tau_{1})\dot{\zeta}(\mathbf{k}_{2},\tau_{1})\dot{\zeta}(\mathbf{k}_{3},\tau_{1})\zeta(\mathbf{p})\zeta(-\mathbf{p})\dot{\zeta}(\mathbf{k}_{4},\tau_{2})\dot{\zeta}(\mathbf{k}_{5},\tau_{2})\dot{\zeta}(\mathbf{k}_{6},\tau_{2})\right\rangle.$
(3.4)

For small loop momentum $k\ll p$, we find

$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle_{(1,1)}^{[aa]}=\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\zeta_{p}(\tau_{0})\zeta_{p}(\tau_{0})$
$\displaystyle\int_{-\infty}^{0}\mathrm{d}\tau_{1}~{}\tau_{1}^{2}\int_{-\infty}^{0}\mathrm{d}\tau_{2}~{}\tau_{2}^{2}e^{-2ic_{s}p(\tau_{1}-\tau_{2})}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\frac{kp^{2}c_{s}^{6}}{2}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}.$ (3.5)

In the limit of loop momentum $k\rightarrow 0$, the momentum integration
becomes

$\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}k\rightarrow 0.$ (3.6)

Thanks to the derivative self-interaction, there is no infrared divergence in
the one-loop correction. The other terms, $[bb]$, $[ab]$, and $[ba]$, also
generate the same loop momentum integration in the limit of small loop
momentum. The $(0,2)$ contribution has a slightly different time integration
with the same momentum integration.

Hence, the loop correction comes from other domain, namely, $k>p$. In this
domain, the divergent integral comes from the $(0,2)$ contribution suffers
from a divergence, while the $(1,1)$ contribution converges in a similar way
to [22]. The $(0,2)$ contribution of $[aa]$ term is given by

$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle_{(0,2)}^{[aa]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{4}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{4}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\prod_{a=1}^{6}\left[\frac{\mathrm{d}^{3}k_{a}}{(2\pi)^{3}}\right]\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})~{}\delta^{3}(\mathbf{k}_{4}+\mathbf{k}_{5}+\mathbf{k}_{6})$
$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\dot{\zeta}(\mathbf{k}_{1},\tau_{1})\dot{\zeta}(\mathbf{k}_{2},\tau_{1})\dot{\zeta}(\mathbf{k}_{3},\tau_{1})\dot{\zeta}(\mathbf{k}_{4},\tau_{2})\dot{\zeta}(\mathbf{k}_{5},\tau_{2})\dot{\zeta}(\mathbf{k}_{6},\tau_{2})\right\rangle.$
(3.7)

Performing Wick contraction, it becomes

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[aa]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\zeta_{p}(\tau_{0})\zeta_{p}(\tau_{0})\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{4}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{4}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle
36\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}_{k}^{*}(\tau_{2})\dot{\zeta}_{p}^{*}(\tau_{1})\dot{\zeta}_{p}^{*}(\tau_{2})\dot{\zeta}_{q}(\tau_{1})\dot{\zeta}_{q}^{*}(\tau_{2}),$
(3.8)

where $\mathbf{q}=\mathbf{k}-\mathbf{p}$. After substituting mode function
(2.26) and performing time integration, we are left with integration over loop
momentum $k$ of a complicated function. The time integration is performed by
adding a small $i\epsilon$ prescription to $\tau\rightarrow\tau(1+i\epsilon)$,
so that the integral at early time $\tau\rightarrow-\infty$ converges. Because
the integration domain is $k>p$, we can approximate that complicated function
by expanding $p\ll k$. As a result, we will get a function with some different
power of $k$. For cosmological interest, the most important term is the one
with the integrand proportional to $k^{-3}$ up to the slow-roll parameters.
Other terms with power of integrand higher than $k^{-3}$ lead to positive
power of $k$ after integration. These terms will generate polynomial
divergences and should be regularized and renormalized in the same way as QFT
in the flat spacetime. Extracting terms with integrand approximately
proportional to $k^{-3}$ yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[aa]}$
$\displaystyle\cong\frac{3}{40}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle=\frac{3}{40}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(3.9)

The detailed calculation is shown in appendix A. The $(0,2)$ contribution of
$[ab]$ term is given by

$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle_{(0,2)}^{[ab]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{4}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{2}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\prod_{a=1}^{6}\left[\frac{\mathrm{d}^{3}k_{a}}{(2\pi)^{3}}\right]\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})~{}\delta^{3}(\mathbf{k}_{4}+\mathbf{k}_{5}+\mathbf{k}_{6})(\mathbf{k}_{5}\cdot\mathbf{k}_{6})$
$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\dot{\zeta}(\mathbf{k}_{1},\tau_{1})\dot{\zeta}(\mathbf{k}_{2},\tau_{1})\dot{\zeta}(\mathbf{k}_{3},\tau_{1})\dot{\zeta}(\mathbf{k}_{4},\tau_{2})\zeta(\mathbf{k}_{5},\tau_{2})\zeta(\mathbf{k}_{6},\tau_{2})\right\rangle.$
(3.10)

Performing Wick contraction, it becomes

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ab]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\zeta_{p}(\tau_{0})\zeta_{p}(\tau_{0})\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{4}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{2}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left[-12(\mathbf{k}\cdot\mathbf{q})\dot{\zeta}_{p}^{*}(\tau_{1})\dot{\zeta}_{p}^{*}(\tau_{2})\dot{\zeta}_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\dot{\zeta}_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.\hphantom{\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}}-8(\mathbf{p}\cdot\mathbf{k})\dot{\zeta}_{q}(\tau_{1})\dot{\zeta}_{q}^{*}(\tau_{2})\dot{\zeta}_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\dot{\zeta}_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\right.$
$\displaystyle\left.\hphantom{\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}}+16(\mathbf{p}\cdot\mathbf{q})\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}_{k}^{*}(\tau_{2})\dot{\zeta}_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\dot{\zeta}_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\right].$
(3.11)

The $(0,2)$ contribution of $[ba]$ term is given by

$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle_{(0,2)}^{[ba]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{2}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{4}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\prod_{a=1}^{6}\left[\frac{\mathrm{d}^{3}k_{a}}{(2\pi)^{3}}\right]\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})~{}\delta^{3}(\mathbf{k}_{4}+\mathbf{k}_{5}+\mathbf{k}_{6})(\mathbf{k}_{2}\cdot\mathbf{k}_{3})$
$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\dot{\zeta}(\mathbf{k}_{1},\tau_{1})\zeta(\mathbf{k}_{2},\tau_{1})\zeta(\mathbf{k}_{3},\tau_{1})\dot{\zeta}(\mathbf{k}_{4},\tau_{2})\dot{\zeta}(\mathbf{k}_{5},\tau_{2})\dot{\zeta}(\mathbf{k}_{6},\tau_{2})\right\rangle.$
(3.12)

Performing Wick contraction, it becomes

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ba]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\zeta_{p}(\tau_{0})\zeta_{p}(\tau_{0})\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{2}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{4}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left[24(\mathbf{p}\cdot\mathbf{q})\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}_{k}^{*}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\dot{\zeta}_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\dot{\zeta}_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.\hphantom{\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}}+12(\mathbf{k}\cdot\mathbf{q})\dot{\zeta}_{p}^{*}(\tau_{1})\dot{\zeta}_{p}^{*}(\tau_{2})\zeta_{k}(\tau_{1})\dot{\zeta}_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\dot{\zeta}_{q}^{*}(\tau_{2})\right].$
(3.13)

Extracting terms with integrand approximately proportional to $k^{-3}$, we get
the same result for $[ab]$ and $[ba]$. The detailed calculation is shown in
appendix B. Summing them yields

$\displaystyle
2\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ab]}$
$\displaystyle\cong\frac{9}{40c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle=\frac{9}{40c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(3.14)

The $(0,2)$ contribution of $[bb]$ term is given by

$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle_{(0,2)}^{[bb]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{2}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{2}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\prod_{a=1}^{6}\left[\frac{\mathrm{d}^{3}k_{a}}{(2\pi)^{3}}\right]\delta^{3}(\mathbf{k}_{1}+\mathbf{k}_{2}+\mathbf{k}_{3})~{}\delta^{3}(\mathbf{k}_{4}+\mathbf{k}_{5}+\mathbf{k}_{6})(\mathbf{k}_{2}\cdot\mathbf{k}_{3})(\mathbf{k}_{5}\cdot\mathbf{k}_{6})$
$\displaystyle\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\dot{\zeta}(\mathbf{k}_{1},\tau_{1})\zeta(\mathbf{k}_{2},\tau_{1})\zeta(\mathbf{k}_{3},\tau_{1})\dot{\zeta}(\mathbf{k}_{4},\tau_{2})\zeta(\mathbf{k}_{5},\tau_{2})\zeta(\mathbf{k}_{6},\tau_{2})\right\rangle.$
(3.15)

Performing Wick contraction, it becomes

$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb]}=-M_{\mathrm{pl}}^{4}\zeta_{p}(\tau_{0})\zeta_{p}(\tau_{0})\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{2}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{2}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}f(p,\tau_{1},\tau_{2}),\\\
$ $\displaystyle
f(p,\tau_{1},\tau_{2})=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}$
$\displaystyle\left[8(\mathbf{p}\cdot\mathbf{q})^{2}\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}^{*}_{k}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.-8(\mathbf{p}\cdot\mathbf{k})(\mathbf{p}\cdot\mathbf{q})\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}^{*}_{q}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\right.$
$\displaystyle\left.-8(\mathbf{p}\cdot\mathbf{q})(\mathbf{k}\cdot\mathbf{q})\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}^{*}_{p}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.-8(\mathbf{p}\cdot\mathbf{q})(\mathbf{k}\cdot\mathbf{q})\dot{\zeta}^{*}_{p}(\tau_{1})\dot{\zeta}^{*}_{k}(\tau_{2})\zeta_{k}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.+4(\mathbf{k}\cdot\mathbf{q})^{2}\dot{\zeta}^{*}_{p}(\tau_{1})\dot{\zeta}^{*}_{p}(\tau_{2})\zeta_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right].$
(3.16)

Extracting terms with integrand approximately proportional to $k^{-3}$ yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb]}$
$\displaystyle\cong\frac{51}{80c_{s}^{4}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle=\frac{51}{80c_{s}^{4}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(3.17)

The detailed calculation is shown in appendix C. Summing all terms $[aa]$,
$[bb]$, $[ab]$, and $[ba]$, the total $(0,2)$ contribution to the one-loop
correction is

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}=\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\left[\frac{51}{80c_{s}^{4}}+\frac{3}{20c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{30}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right]$
$\displaystyle\times\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(3.18)

For $n_{s}<1$, which is the case in our Universe, the momentum integration
converges as [26]

$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}=\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\frac{\Delta^{2}_{s(0)}(p)}{1-n_{s}}\left[\frac{51}{80c_{s}^{4}}+\frac{3}{20c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{30}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right].$
(3.19)

Hence, adding the $(2,0)$ contribution, the total one-loop correction to the
power spectrum is

$\Delta^{2}_{s(1)}(p)=\frac{[\Delta^{2}_{s(0)}(p)]^{2}}{1-n_{s}}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right].$
(3.20)

Perturbativity requirement $\Delta^{2}_{s(1)}(p)\ll\Delta^{2}_{s(0)}(p)$ reads

$\frac{\Delta^{2}_{s(0)}(p)}{1-n_{s}}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right]\ll
1.$ (3.21)

The total power spectrum as a summation of tree-level and one-loop correction
is

$\Delta^{2}_{s}(p)=\Delta^{2}_{s(0)}(p)\left\\{1+\frac{\Delta^{2}_{s(0)}(p)}{1-n_{s}}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right]\right\\}.$
(3.22)

We can write the explicit momentum dependence of total power spectrum as

$\Delta^{2}_{s}(p)=\Delta^{2}_{s}(p_{*})\left(\frac{p}{p_{*}}\right)^{N_{s}-1},$
(3.23)

where

$\displaystyle N_{s}-1$
$\displaystyle=\frac{\mathrm{d}\log\Delta^{2}_{s}}{\mathrm{d}\log p}$
$\displaystyle=\frac{\mathrm{d}\log\Delta^{2}_{s(0)}}{\mathrm{d}\log
p}+\frac{\mathrm{d}}{\mathrm{d}\log
p}\log\left\\{1+\frac{\Delta^{2}_{s(0)}}{1-n_{s}}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right]\right\\}$
$\displaystyle\approx\frac{\mathrm{d}\log\Delta^{2}_{s(0)}}{\mathrm{d}\log
p}+\frac{\mathrm{d}}{\mathrm{d}\log
p}\left\\{\frac{\Delta^{2}_{s(0)}}{1-n_{s}}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right]\right\\}$
$\displaystyle=n_{s}-1-\Delta^{2}_{s(0)}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right].$
(3.24)

One might wonder whether we can have a smooth exact scale-invariant limit
$n_{s}=1$ of the total power spectrum (3.22). For the exactly scale-invariant
case, one-loop correction is proportional to

$\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}k^{3}}$ (3.25)

Performing dynamical dimensional regularization [33, 34, 35]

$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}k^{3}}\rightarrow\mu^{\epsilon}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}k^{3+\epsilon}}$
$\displaystyle=\frac{\mu^{\epsilon}}{2\pi^{2}}\int_{p}^{\infty}\frac{\mathrm{d}k}{k^{1+\epsilon}}$
$\displaystyle=\frac{1}{2\pi^{2}\epsilon}\left(\frac{p}{\mu}\right)^{-\epsilon}\approx\frac{1}{2\pi^{2}}\left(\frac{1}{\epsilon}-\log\frac{p}{\mu}\right),$
(3.26)

where we introduce a very small positive regulator $\epsilon$ and momentum
scale $\mu$ to make the integral dimensionless. The one-loop corrected power
spectrum is

$\Delta^{2}_{s}(p)=\Delta^{2}_{s(0)}\left\\{1+\left(\frac{1}{\epsilon}-\log\frac{p}{\mu}\right)\Delta^{2}_{s(0)}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right]\right\\},$
(3.27)

which contains divergence $1/\epsilon$. To absorb the divergence, we introduce
a renormalization factor $\tilde{\Delta}^{2}_{s}(p)=Z\Delta^{2}_{s}(p)$ with

$Z=1-\left(\frac{1}{\epsilon}-\log\frac{\tilde{p}}{\mu}\right)\Delta^{2}_{s(0)}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right],$
(3.28)

and perform renormalization:

$\tilde{\Delta}^{2}_{s}(p)=\Delta^{2}_{s(0)}\left\\{1-\Delta^{2}_{s(0)}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right]\log\frac{p}{\tilde{p}}\right\\}.$
(3.29)

This renormalization fixes the total power spectrum at scale $\tilde{p}$ as
$\tilde{\Delta}^{2}_{s}(\tilde{p})=\Delta^{2}_{s(0)}$. Observable
$\tilde{\Delta}^{2}_{s}(p)$ must not depend on arbitrary parameter
$\tilde{p}$. It can be realised if $\Delta^{2}_{s(0)}$ runs as a function of
$\tilde{p}$, which satisfies the following differential equation

$\displaystyle 0$
$\displaystyle=\frac{\partial}{\partial\log\tilde{p}}\log\tilde{\Delta}^{2}_{s}(p)$
$\displaystyle\approx\frac{\partial}{\partial\log\tilde{p}}\left\\{\log\Delta^{2}_{s(0)}-\Delta^{2}_{s(0)}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right]\log\frac{p}{\tilde{p}}\right\\}$
$\displaystyle\approx\frac{\partial\log\Delta^{2}_{s(0)}}{\partial\log\tilde{p}}+\Delta^{2}_{s(0)}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right],$
$\displaystyle\frac{\partial\log\Delta^{2}_{s(0)}}{\partial\log\tilde{p}}$
$\displaystyle=-\Delta^{2}_{s(0)}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right].$
(3.30)

By identifying $\Delta^{2}_{s(0)}=\tilde{\Delta}^{2}_{s}(\tilde{p})$ and
redefine $\tilde{p}\rightarrow p$, we get

$\tilde{\Delta}^{2}_{s}(p)=\tilde{\Delta}^{2}_{s}(p_{*})\left(\frac{p}{p_{*}}\right)^{N_{s}-1},$
(3.31)

where

$N_{s}-1=-\Delta^{2}_{s(0)}\left[\frac{51}{40c_{s}^{4}}+\frac{3}{10c_{s}^{2}}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)+\frac{1}{15}\left(\frac{\tilde{c}_{3}}{c_{s}^{2}}\right)^{2}\right].$
(3.32)

Therefore, there is a smooth limit in which we recover the momentum dependence
of renormalized power spectrum of the exact scale-invariant result, by
directly substituting $n_{s}=1$ to (3.24).

##  4 Constraint on Parameter Space

Perturbativity requirement (3.21) leads to a constraint on parameter space
$\tilde{c}_{3}$ and $c_{s}$. In Figure 1, we compare our theoretical
constraint to observational constraint [11]. We know that
$\Delta_{s(0)}^{2}\equiv A_{2}=2.1\times 10^{-9}$ at pivot scale
$k_{\ast}=0.05~{}\mathrm{Mpc}^{-1}$ and $n_{s}=0.9649\pm 0.0042$ based on
$TT$, $TE$, $EE$, low $E$, and lens of Planck 2018 [2]. The curve
$\Delta^{2}_{s(1)}=\Delta^{2}_{s(0)}$ is shown by the red curve. Outside the
red curve, perturbativity no longer works and we must not trust perturbation
theory anymore. The parameter space that is allowed by Planck 2018 observation
is enclosed by the blue curve. We can see that the allowed parameter space is
inside the red curve, so the observationally allowed parameters have one-loop
correction less than tree-level contribution and perturbativity treatment
works.

However, it is important to require strong inequality on the one-loop
correction compared to the tree-level contribution, namely
$\Delta^{2}_{s(1)}\ll\Delta^{2}_{s(0)}$. Hence, we also draw the curve
$\Delta^{2}_{s(1)}=0.1\Delta^{2}_{s(0)}$ and
$\Delta^{2}_{s(1)}=0.01\Delta^{2}_{s(0)}$, shown by the magenta and green
curve, respectively. For the region inside the green curve, we can safely
neglect one-loop correction because it is very small. For the region enclosed
by the magenta and the green curve, we might consider to add one-loop
correction, depending on the precision that we want to achieve. For
observationally allowed parameter between red and magenta curve, one-loop
correction is large enough, greater than $0.1$ times tree-level contribution.
Hence, if future observations squeeze the parameter space into the region
between red and magenta curve, then considering one-loop correction in the
theory will be important to achieve precise cosmology.

Figure 1:  Comparison between theoretical and observational constraint on
$\tilde{c}_{3}$ and $c_{s}$ parameter space. Region inside blue curve is
Planck 2018 observationally allowed parameters with $99.7\%$ confidence level
[11]. Theoretical constraint when one-loop correction $\Delta^{2}_{s(1)}$
equals to tree-level contribution $\Delta^{2}_{s(0)}$ is shown by the red
curve. Another constraint when $\Delta^{2}_{s(1)}=0.1\Delta^{2}_{s(0)}$ and
$\Delta^{2}_{s(1)}=0.01\Delta^{2}_{s(0)}$ are shown by the magenta and green
curve, respectively.

##  5 Conclusion

Perturbativity may not be taken for granted. We showed that it imposes non-
trivial constraints to the theory. Requiring one-loop correction to the two-
point functions of curvature perturbation to be much smaller than its tree-
level contribution leads to an upper bound on primordial non-Gaussianity. We
compared this theoretical constraint to observational constraint on non-
Gaussian parameter space. Although all observationally allowed parameters have
one-loop correction less than tree-level contribution, a quite large region in
the parameter space suffers from one-loop correction greater than $0.1$ times
tree-level contribution. Therefore, if future observations find non-
Gaussianity inside this region, it will be important to consider one-loop
correction in the cosmological perturbation theory.

##  Appendix A Calculation of $\dot{\zeta}^{3}\times\dot{\zeta}^{3}$ term

In this appendix, we calculate the one-loop correction from
$\dot{\zeta}^{3}\times\dot{\zeta}^{3}$ term, that we label as $[aa]$ term.
From (3.8), it is given by

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[aa]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{4}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{4}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle
36\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}_{k}^{*}(\tau_{2})\dot{\zeta}_{p}^{*}(\tau_{1})\dot{\zeta}_{p}^{*}(\tau_{2})\dot{\zeta}_{q}(\tau_{1})\dot{\zeta}_{q}^{*}(\tau_{2}).$
(A.1)

Substituting mode function (2.26) leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[aa]}=$
$\displaystyle-\frac{9c_{s}^{6}}{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}kpq$
$\displaystyle\int_{-\infty}^{0}\mathrm{d}\tau_{1}~{}\tau_{1}^{2}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}~{}\tau_{2}^{2}e^{ic_{s}(p-k-q)\tau_{1}}e^{ic_{s}(p+k+q)\tau_{2}}.$
(A.2)

Then, performing time integration yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[aa]}=$
$\displaystyle~{}\frac{9c_{s}^{6}}{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}kpq$
$\displaystyle\times\frac{3k^{2}+9kp+6kq+8p^{2}+9pq+3q^{2}}{4c_{s}^{6}p^{5}(k+p+q)^{3}}.$
(A.3)

Substituting $q^{2}=k^{2}-2\mathbf{p}\cdot\mathbf{k}+p^{2}$, expanding with
respect to the power of $p/k\ll 1$, and extracting terms that have integrand
approximately proportional to $k^{-3}$ leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[aa]}=$
$\displaystyle-\frac{9c_{s}^{6}}{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{1}{c_{s}^{6}p^{4}k^{3}}\left[\frac{21s^{4}}{64}-\frac{11p^{2}s^{2}}{32}+\frac{p^{4}}{64}\right],$
(A.4)

where $s=-\mathbf{p}\cdot\hat{\mathbf{k}}$. Next, substituting
$\hat{\mathbf{p}}\cdot\hat{\mathbf{k}}=\cos\theta$ and decomposing volume
element yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[aa]}=$
$\displaystyle-\frac{9}{8}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\int_{-1}^{1}\mathrm{d}\cos\theta\left[\frac{21}{64}\cos^{4}\theta-\frac{11}{32}\cos^{2}\theta+\frac{1}{64}\right].$
(A.5)

We will use the same definition of $s$ and $\theta$ in the next appendix.
Finally, after performing angular integration we get

$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[aa]}=\frac{3}{40}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)^{2}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}.$ (A.6)

##  Appendix B Calculation of
$\dot{\zeta}^{3}\times\dot{\zeta}(\partial\zeta)^{2}$ term

In this appendix, we calculate the one-loop correction from the
$\dot{\zeta}^{3}\times\dot{\zeta}(\partial\zeta)^{2}$ cross term, that we
label as $[ab]$ and $[ba]$. From (3.11), the $[ab]$ term is given by

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ab]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{4}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{2}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left[-12(\mathbf{k}\cdot\mathbf{q})\dot{\zeta}_{p}^{*}(\tau_{1})\dot{\zeta}_{p}^{*}(\tau_{2})\dot{\zeta}_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\dot{\zeta}_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.\hphantom{\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}}-8(\mathbf{p}\cdot\mathbf{k})\dot{\zeta}_{q}(\tau_{1})\dot{\zeta}_{q}^{*}(\tau_{2})\dot{\zeta}_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\dot{\zeta}_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\right.$
$\displaystyle\left.\hphantom{\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}}+16(\mathbf{p}\cdot\mathbf{q})\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}_{k}^{*}(\tau_{2})\dot{\zeta}_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\dot{\zeta}_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\right].$
(B.1)

Substituting mode function (2.26) leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ab]}=$
$\displaystyle-\frac{1}{16}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int_{-\infty}^{0}\mathrm{d}\tau_{1}~{}\tau_{1}^{2}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{c_{s}^{2}}{pkq}e^{ic_{s}(p-k-q)\tau_{1}}e^{ic_{s}(p+k+q)\tau_{2}}$
$\displaystyle\hphantom{\int_{-\infty}^{0}\mathrm{d}\tau_{1}~{}\tau_{1}^{2}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}}\left[-12p^{2}(\mathbf{k}\cdot\mathbf{q})(1-ic_{s}k\tau_{2})(1-ic_{s}q\tau_{2})\right.$
$\displaystyle\hphantom{\int_{-\infty}^{0}\mathrm{d}\tau_{1}~{}\tau_{1}^{2}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}}\left.-8q^{2}(\mathbf{p}\cdot\mathbf{k})(1-ic_{s}p\tau_{2})(1-ic_{s}k\tau_{2})\right.$
$\displaystyle\hphantom{\int_{-\infty}^{0}\mathrm{d}\tau_{1}~{}\tau_{1}^{2}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}}\left.+16k^{2}(\mathbf{p}\cdot\mathbf{q})(1-ic_{s}p\tau_{2})(1-ic_{s}q\tau_{2})\right].$
(B.2)

Then, performing time integration yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ab]}=$
$\displaystyle-\frac{1}{16c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{1}{pkq}$
$\displaystyle\frac{(-k)}{2p^{4}(k+p+q)^{3}}\left[s\left(q^{2}\left(98k^{2}p+54k^{3}+15kp^{2}+16p^{3}\right)\right.\right.$
$\displaystyle\left.+3q\left(71k^{2}p^{2}+54k^{3}p+12k^{4}+36kp^{3}+9p^{4}\right)-2q^{4}(9k+5p)\right.$
$\displaystyle\left.+p(k+p)\left(37k^{2}p+20k^{3}+21kp^{2}+6p^{3}\right)-3pq^{3}(18k+5p)\right)$
$\displaystyle\left.+kp\left(q^{2}\left(36k^{2}+153kp+154p^{2}\right)+3q\left(39k^{2}p+6k^{3}+84kp^{2}+53p^{3}\right)\right.\right.$
$\displaystyle\left.\left.+p(k+p)\left(9k^{2}+41kp+34p^{2}\right)+9q^{3}(2k+5p)\right)\right].$
(B.3)

Substituting $q^{2}=k^{2}-2\mathbf{p}\cdot\mathbf{k}+p^{2}$, expanding with
respect to the power of $p/k\ll 1$, and extracting terms that have integrand
approximately proportional to $k^{-3}$ leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ab]}$
$\displaystyle=-\frac{1}{16c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\hphantom{=}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\left[\frac{341s^{4}}{16k^{3}p^{4}}-\frac{357s^{2}}{16k^{3}p^{2}}+\frac{11}{8k^{3}}\right]$
$\displaystyle=-\frac{1}{32c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\hphantom{=}\int_{-1}^{1}\mathrm{d}\cos\theta\left[\frac{341}{16}\cos^{4}\theta-\frac{357}{16}\cos^{2}\theta+\frac{11}{8}\right].$
(B.4)

Finally, after performing angular integration we get

$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ab]}=\frac{9}{80c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}.$ (B.5)

Next, we calculate the $[ba]$ term. From (3.13), it is given by

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ba]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{2}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{4}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left[24(\mathbf{p}\cdot\mathbf{q})\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}_{k}^{*}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\dot{\zeta}_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\dot{\zeta}_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.\hphantom{\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}}+12(\mathbf{k}\cdot\mathbf{q})\dot{\zeta}_{p}^{*}(\tau_{1})\dot{\zeta}_{p}^{*}(\tau_{2})\zeta_{k}(\tau_{1})\dot{\zeta}_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\dot{\zeta}_{q}^{*}(\tau_{2})\right].$
(B.6)

Substituting mode function (2.26) leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ba]}=$
$\displaystyle-\frac{1}{16}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int_{-\infty}^{0}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}~{}\tau_{2}^{2}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{c_{s}^{2}}{pkq}e^{ic_{s}(p-k-q)\tau_{1}}e^{ic_{s}(p+k+q)\tau_{2}}$
$\displaystyle\hphantom{\int_{-\infty}^{0}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}~{}\tau_{2}^{2}}\left[24k^{2}(\mathbf{p}\cdot\mathbf{q})(1-ic_{s}p\tau_{1})(1+ic_{s}q\tau_{1})\right.$
$\displaystyle\hphantom{\int_{-\infty}^{0}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}~{}\tau_{2}^{2}}\left.+12p^{2}(\mathbf{k}\cdot\mathbf{q})(1+ic_{s}k\tau_{1})(1+ic_{s}q\tau_{1})\right].$
(B.7)

Then, performing time integration yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ba]}=$
$\displaystyle-\frac{1}{16c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{1}{pkq}$
$\displaystyle\frac{3k}{2p^{4}(k+p+q)^{3}}\left[s\left(k^{3}\left(-33p^{2}+30pq+36q^{2}\right)+2k^{4}(9q-5p)\right.\right.$
$\displaystyle\left.+k^{2}\left(9p^{2}q-42p^{3}+34pq^{2}+18q^{3}\right)+3kp(p+q)^{2}(p-2q)\right.$
$\displaystyle\left.+p^{2}\left(3p^{2}q-14p^{3}+8pq^{2}+3q^{3}\right)\right)$
$\displaystyle\left.+kp\left(3q^{2}\left(-4k^{2}+9kp+18p^{2}\right)+3q^{3}(7p-2k)\right.\right.$
$\displaystyle\left.\left.-p\left(2k^{2}p-3k^{3}+33kp^{2}+64p^{3}\right)+3q(k+p)^{2}(7p-2k)\right)\right].$
(B.8)

Substituting $q^{2}=k^{2}-2\mathbf{p}\cdot\mathbf{k}+p^{2}$, expanding with
respect to the power of $p/k\ll 1$, and extracting terms that have integrand
approximately proportional to $k^{-3}$ leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ba]}$
$\displaystyle=-\frac{1}{16c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\hphantom{=}\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\left[\frac{243s^{4}}{8k^{3}p^{4}}-\frac{513s^{2}}{16k^{3}p^{2}}+\frac{45}{16k^{3}}\right]$
$\displaystyle=-\frac{1}{32c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\hphantom{=}\int_{-1}^{1}\mathrm{d}\cos\theta\left[\frac{243}{8}\cos^{4}\theta-\frac{513}{16}\cos^{2}\theta+\frac{45}{16}\right].$
(B.9)

Finally, after performing angular integration we get

$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[ba]}=\frac{9}{80c_{s}^{2}}\left(\frac{2\tilde{c}_{3}}{3c_{s}^{2}}\right)\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}.$ (B.10)

##  Appendix C Calculation of
$\dot{\zeta}(\partial\zeta)^{2}\times\dot{\zeta}(\partial\zeta)^{2}$ term

In this appendix, we calculate the one-loop correction from
$\dot{\zeta}(\partial\zeta)^{2}\times\dot{\zeta}(\partial\zeta)^{2}$ term,
that we label as $[bb]$ term. From (3.16), it is given by

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{2}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{2}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}f(p,\tau_{1},\tau_{2}),$
$\displaystyle f(p,\tau_{1},\tau_{2})=\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}$
$\displaystyle\left[8(\mathbf{p}\cdot\mathbf{q})^{2}\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}^{*}_{k}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.-8(\mathbf{p}\cdot\mathbf{k})(\mathbf{p}\cdot\mathbf{q})\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}^{*}_{q}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\right.$
$\displaystyle\left.-8(\mathbf{p}\cdot\mathbf{q})(\mathbf{k}\cdot\mathbf{q})\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}^{*}_{p}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.-8(\mathbf{p}\cdot\mathbf{q})(\mathbf{k}\cdot\mathbf{q})\dot{\zeta}^{*}_{p}(\tau_{1})\dot{\zeta}^{*}_{k}(\tau_{2})\zeta_{k}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right.$
$\displaystyle\left.+4(\mathbf{k}\cdot\mathbf{q})^{2}\dot{\zeta}^{*}_{p}(\tau_{1})\dot{\zeta}^{*}_{p}(\tau_{2})\zeta_{k}(\tau_{1})\zeta_{k}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2})\right].$
(C.1)

There are five terms in (C.1), and we will calculate them one by one. The
first term is

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb1]}=$
$\displaystyle-
M_{\mathrm{pl}}^{4}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{-\infty}^{0}\mathrm{d}\tau_{1}\frac{a^{2}(\tau_{1})\epsilon(\tau_{1})}{H(\tau_{1})c_{s}^{2}}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}\frac{a^{2}(\tau_{2})\epsilon(\tau_{2})}{H(\tau_{2})c_{s}^{2}}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}(\mathbf{p}\cdot\mathbf{q})^{2}\dot{\zeta}_{k}(\tau_{1})\dot{\zeta}^{*}_{k}(\tau_{2})\zeta_{p}^{*}(\tau_{1})\zeta_{p}^{*}(\tau_{2})\zeta_{q}(\tau_{1})\zeta_{q}^{*}(\tau_{2}).$
(C.2)

Substituting mode function (2.26) leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb1]}=-\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}(p^{4}+(\mathbf{p}\cdot\mathbf{k})^{2}-2p^{2}\mathbf{p}\cdot\mathbf{k})$
$\displaystyle\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{k}{2c_{s}^{2}p^{3}q^{3}}f_{1}(p,q,k),$ (C.3)
$\displaystyle
f_{1}(p,q,k)=\int_{-\infty}^{0}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}$
$\displaystyle~{}e^{ic_{s}(p-k-q)\tau_{1}}e^{ic_{s}(p+k+q)\tau_{2}}$
$\displaystyle(1-ic_{s}p\tau_{1})(1-ic_{s}p\tau_{2})(1+ic_{s}q\tau_{1})(1-ic_{s}q\tau_{2}).$
(C.4)

Then, performing time integration yields

$\displaystyle f_{1}(p,q,k)=-\frac{1}{4c_{s}^{2}p^{3}(k+p+q)^{3}}$
$\displaystyle\left[-q^{2}\left(7k^{2}+21kp+8p^{2}\right)+3p^{2}q(5k+8p)\right.$
$\displaystyle\left.+p^{2}(k+p)(5k+8p)-7q^{3}(2k+3p)-7q^{4}\right].$ (C.5)

Substituting $q^{2}=k^{2}-2\mathbf{p}\cdot\mathbf{k}+p^{2}$, expanding with
respect to the power of $p/k\ll 1$, and extracting terms that have integrand
approximately proportional to $k^{-3}$ leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb1]}=$
$\displaystyle-\frac{1}{2c_{s}^{4}p^{6}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[p^{4}+(\mathbf{p}\cdot\mathbf{k})^{2}-2p^{2}\mathbf{p}\cdot\mathbf{k}]\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\left[\frac{35s^{2}}{16k^{3}}-\frac{245s^{3}}{64k^{4}}+\frac{441s^{4}}{164k^{5}}-\frac{27p^{2}}{16k^{3}}+\frac{355p^{2}s}{64k^{4}}+\frac{149p^{4}}{64k^{5}}-\frac{231p^{2}s^{2}}{16k^{5}}\right].$
(C.6)

Next, performing phase-volume integration yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb1]}=$
$\displaystyle-\frac{1}{2c_{s}^{4}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\frac{2\pi}{4\pi}\int_{-1}^{1}\mathrm{d}\cos\theta\left[-\frac{27}{16}+\left(\frac{35}{16}+\frac{149}{64}+2\frac{355}{64}\right)\cos^{2}\theta\right.$
$\displaystyle\left.\hphantom{\frac{2\pi}{4\pi}\int_{-1}^{1}\mathrm{d}\cos\theta}-\left(\frac{231}{16}+2\frac{245}{64}\right)\cos^{4}\theta+\frac{441}{164}\cos^{6}\theta\right],$
(C.7)
$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb1]}=-\frac{1}{4c_{s}^{4}}\frac{13}{80}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(C.8)

Doing similar calculation for the second term, substituting mode function
(2.26) leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb2]}=-\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[-(\mathbf{p}\cdot\mathbf{k})^{2}+p^{2}\mathbf{p}\cdot\mathbf{k}]$
$\displaystyle\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{1}{2c_{s}^{2}p^{3}kq}f_{2}(p,q,k),$ (C.9)
$\displaystyle
f_{2}(p,q,k)=\int_{-\infty}^{0}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}$
$\displaystyle~{}e^{ic_{s}(p-k-q)\tau_{1}}e^{ic_{s}(p+k+q)\tau_{2}}$
$\displaystyle(1-ic_{s}p\tau_{1})(1-ic_{s}p\tau_{2})(1+ic_{s}q\tau_{1})(1-ic_{s}k\tau_{2}).$
(C.10)

Then, performing time integration yields

$\displaystyle f_{2}(p,q,k)=-\frac{1}{8c_{s}^{2}p^{3}(k+p+q)^{3}}$
$\displaystyle\left[-q^{2}\left(28k^{2}+51kp+12p^{2}\right)+p(3k+p)\left(3k^{2}+13kp+16p^{2}\right)\right.$
$\displaystyle\left.+q^{3}(-(14k+9p))+q(p-2k)(k+p)(7k+13p)\right].$ (C.11)

Substituting $q^{2}=k^{2}-2\mathbf{p}\cdot\mathbf{k}+p^{2}$, expanding $p\ll
k$, and extracting terms that have integrand approximately proportional to
$k^{-3}$ leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb2]}=$
$\displaystyle-\frac{1}{2c_{s}^{4}p^{6}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[-(\mathbf{p}\cdot\mathbf{k})^{2}+p^{2}\mathbf{p}\cdot\mathbf{k}]\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\left[\frac{49s^{4}}{64k^{5}}+\frac{63p^{4}}{64k^{5}}-\frac{27s^{2}p^{2}}{8k^{5}}-\frac{35s^{3}}{64k^{4}}+\frac{127sp^{2}}{64k^{4}}\right].$
(C.12)

Next, performing phase-volume integration yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb2]}=$
$\displaystyle-\frac{1}{2c_{s}^{4}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\frac{2\pi}{4\pi}\int_{-1}^{1}\mathrm{d}\cos\theta\left[-\left(\frac{63}{64}+\frac{127}{64}\right)\cos^{2}\theta+\left(\frac{27}{8}+\frac{35}{64}\right)\cos^{4}\theta-\frac{49}{64}\cos^{6}\theta\right],$
(C.13)
$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb2]}=\frac{1}{4c_{s}^{4}}\frac{151}{240}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(C.14)

Doing similar calculation for the third term, substituting mode function
(2.26) leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb3]}=-\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[p^{2}k^{2}-(p^{2}+k^{2})\mathbf{p}\cdot\mathbf{k}+(\mathbf{p}\cdot\mathbf{k})^{2}]$
$\displaystyle\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{1}{2c_{s}^{2}pkq^{3}}f_{3}(p,q,k),$ (C.15)
$\displaystyle
f_{3}(p,q,k)=\int_{-\infty}^{0}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}$
$\displaystyle~{}e^{ic_{s}(p-k-q)\tau_{1}}e^{ic_{s}(p+k+q)\tau_{2}}$
$\displaystyle(1-ic_{s}p\tau_{1})(1-ic_{s}q\tau_{2})(1+ic_{s}q\tau_{1})(1-ic_{s}k\tau_{2}).$
(C.16)

Then, performing time integration yields

$\displaystyle f_{3}(p,q,k)=-\frac{1}{8c_{s}^{2}p^{4}(k+p+q)^{3}}$
$\displaystyle\left[q^{4}(-(9k+5p))+3q^{2}\left(-11k^{2}p-3k^{3}-11kp^{2}+p^{3}\right)\right.$
$\displaystyle\left.+6p^{2}q\left(2k^{2}+8kp+3p^{2}\right)+2p^{2}(k+p)\left(2k^{2}+8kp+3p^{2}\right)\right.$
$\displaystyle\left.-2q^{3}\left(9k^{2}+19kp+7p^{2}\right)\right]$ (C.17)

Substituting $q^{2}=k^{2}-2\mathbf{p}\cdot\mathbf{k}+p^{2}$, expanding with
respect to the power of $p/k\ll 1$, and extracting terms that have integrand
approximately proportional to $k^{-3}$ leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb3]}=$
$\displaystyle-\frac{1}{2c_{s}^{4}p^{6}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[p^{2}k^{2}-(p^{2}+k^{2})\mathbf{p}\cdot\mathbf{k}+(\mathbf{p}\cdot\mathbf{k})^{2}]\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\left[-\frac{17sp}{32k^{4}}+\frac{115s^{2}p}{128k^{5}}-\frac{203s^{3}}{128k^{6}}-\frac{131p^{3}}{128k^{5}}+\frac{439sp^{3}}{128k^{6}}\right].$
(C.18)

Next, performing phase-volume integration yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb3]}=$
$\displaystyle-\frac{1}{2c_{s}^{4}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\frac{2\pi}{4\pi}\int_{-1}^{1}\mathrm{d}\cos\theta\left[-\frac{131}{128}+\left(\frac{115}{128}-\frac{17}{32}-\frac{131}{128}+\frac{439}{128}\right)\cos^{2}\theta\right.$
$\displaystyle\left.\hphantom{\frac{2\pi}{4\pi}\int_{-1}^{1}\mathrm{d}\cos\theta}+\left(\frac{115}{128}-\frac{203}{128}\right)\cos^{4}\theta\right]$
(C.19)
$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb3]}=\frac{1}{4c_{s}^{4}}\frac{227}{480}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(C.20)

Doing similar calculation for the fourth term, substituting mode function
(2.26) leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb4]}=-\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[p^{2}k^{2}-(p^{2}+k^{2})\mathbf{p}\cdot\mathbf{k}+(\mathbf{p}\cdot\mathbf{k})^{2}]$
$\displaystyle\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{1}{2c_{s}^{2}pkq^{3}}f_{4}(p,q,k),$ (C.21)
$\displaystyle
f_{4}(p,q,k)=\int_{-\infty}^{0}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}$
$\displaystyle~{}e^{ic_{s}(p-k-q)\tau_{1}}e^{ic_{s}(p+k+q)\tau_{2}}$
$\displaystyle(1+ic_{s}k\tau_{1})(1-ic_{s}q\tau_{2})(1+ic_{s}q\tau_{1})(1-ic_{s}p\tau_{2}).$
(C.22)

Then, performing time integration yields

$\displaystyle f_{4}(p,q,k)=\frac{1}{8c_{s}^{2}p^{4}(k+p+q)^{3}}$
$\displaystyle\left[-2q^{3}\left(9k^{2}+8kp-8p^{2}\right)-3q^{2}(3k+p)\left(k^{2}+2kp-p^{2}\right)+\right.$
$\displaystyle\left.6p^{2}q\left(2k^{2}-5p^{2}\right)-2p^{2}(k+p)\left(5p^{2}-2k^{2}\right)+q^{4}(5p-9k)\right].$
(C.23)
$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb4]}=$
$\displaystyle-\frac{1}{2c_{s}^{4}p^{6}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[p^{2}k^{2}-(p^{2}+k^{2})\mathbf{p}\cdot\mathbf{k}+(\mathbf{p}\cdot\mathbf{k})^{2}]\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\left[-\frac{17sp}{32k^{4}}+\frac{115s^{2}p}{128k^{5}}-\frac{203s^{3}}{128k^{6}}-\frac{131p^{3}}{128k^{5}}+\frac{439sp^{3}}{128k^{6}}\right].$
(C.24)

The fourth term has the same contribution as the third term. Doing similar
calculation for the fifth term, substituting mode function (2.26) leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb5]}=-\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[k^{4}+(\mathbf{p}\cdot\mathbf{k})^{2}-2k^{2}\mathbf{p}\cdot\mathbf{k}]$
$\displaystyle\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}\frac{p}{4c_{s}^{2}k^{3}q^{3}}f_{5}(p,q,k),$ (C.25)
$\displaystyle
f_{5}(p,q,k)=\int_{-\infty}^{0}\mathrm{d}\tau_{1}\int_{-\infty}^{\tau_{1}}\mathrm{d}\tau_{2}$
$\displaystyle~{}e^{ic_{s}(p-k-q)\tau_{1}}e^{ic_{s}(p+k+q)\tau_{2}}$
$\displaystyle(1+ic_{s}q\tau_{1})(1-ic_{s}q\tau_{2})(1+ic_{s}k\tau_{1})(1-ic_{s}k\tau_{2}).$
(C.26)

Then, performing time integration yields

$\displaystyle f_{5}(p,q,k)=-\frac{1}{4c_{s}^{2}p^{5}(k+p+q)^{3}}$
$\displaystyle\left[q^{4}\left(3k^{2}-p^{2}\right)+3q^{3}(2k+p)\left(k^{2}+kp-p^{2}\right)\right.$
$\displaystyle\left.+q^{2}\left(4k^{2}p^{2}+9k^{3}p+3k^{4}-6kp^{3}+2p^{4}\right)\right.$
$\displaystyle\left.-3p^{2}q\left(2k^{2}p+k^{3}-4kp^{2}-2p^{3}\right)\right.$
$\displaystyle\left.+p^{2}(k+p)\left(-2k^{2}p-k^{3}+4kp^{2}+2p^{3}\right)\right].$
(C.27)

Substituting $q^{2}=k^{2}-2\mathbf{p}\cdot\mathbf{k}+p^{2}$, expanding with
respect to the power of $p/k\ll 1$, and extracting terms that have integrand
approximately proportional to $k^{-3}$ leads to

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb5]}=$
$\displaystyle-\frac{1}{4c_{s}^{4}p^{6}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}$
$\displaystyle\int\frac{\mathrm{d}^{3}k}{(2\pi)^{3}}[k^{4}+(\mathbf{p}\cdot\mathbf{k})^{2}-2k^{2}\mathbf{p}\cdot\mathbf{k}]\left(\frac{H^{2}}{4M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\left[-\frac{15s^{2}}{16k^{5}}+\frac{105s^{3}}{64k^{6}}-\frac{189s^{4}}{64k^{7}}+\frac{11p^{2}}{16k^{5}}-\frac{125sp^{2}}{64k^{6}}+\frac{77s^{2}p^{2}}{16k^{7}}-\frac{99p^{4}}{64k^{7}}\right].$
(C.28)

Next, performing phase-volume integration yields

$\displaystyle\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb5]}=$
$\displaystyle-\frac{1}{4c_{s}^{4}}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\left(\frac{H^{2}}{8\pi^{2}M_{\mathrm{pl}}^{2}\epsilon
c_{s}}\right)_{H}$
$\displaystyle\frac{2\pi}{4\pi}\int_{-1}^{1}\mathrm{d}\cos\theta\left[-\frac{99}{64}+\left(\frac{11}{16}+\frac{77}{16}-2\frac{125}{64}\right)\cos^{2}\theta\right.$
$\displaystyle\left.\hphantom{\frac{2\pi}{4\pi}\int_{-1}^{1}\mathrm{d}\cos\theta}+\left(-\frac{15}{16}+2\frac{105}{64}-\frac{189}{64}\right)\cos^{4}\theta\right],$
(C.29)
$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb5]}=\frac{1}{8c_{s}^{4}}\frac{91}{40}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(C.30)

Finally, after adding all terms we get

$\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0,2)}^{[bb]}=\frac{1}{4c_{s}^{4}}\frac{51}{20}\left\langle\\!\left\langle\zeta(\mathbf{p})\zeta(-\mathbf{p})\right\rangle\\!\right\rangle_{(0)}\int_{p}^{\infty}\mathrm{d}k~{}k^{2}\frac{1}{k^{3}}\Delta^{2}_{s(0)}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{s}-1}.$
(C.31)

## Acknowledgments

J. K. acknowledges the support from JSPS KAKENHI Grant No. 22J20289 and Global
Science Graduate Course (GSGC) program of The University of Tokyo. J. Y. is
supported by JSPS KAKENHI Grant No. 20H05639 and Grant-in-Aid for Scientific
Research on Innovative Areas 20H05248.

## References

  * [1] Planck collaboration, _Planck 2018 results. I. Overview and the cosmological legacy of Planck_ , _Astron. Astrophys._ 641 (2020) A1 [1807.06205].
  * [2] Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [3] A.A. Starobinsky, _A New Type of Isotropic Cosmological Models Without Singularity_ , _Phys. Lett. B_ 91 (1980) 99.
  * [4] K. Sato, _First Order Phase Transition of a Vacuum and Expansion of the Universe_ , _Mon. Not. Roy. Astron. Soc._ 195 (1981) 467.
  * [5] A.H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys. Rev._ D23 (1980) 347.
  * [6] A.A. Starobinsky, _Spectrum of relict gravitational radiation and the early state of the universe_ , _JETP Lett._ 30 (1979) 682.
  * [7] V.F. Mukhanov and G.V. Chibisov, _Quantum Fluctuations and a Nonsingular Universe_ , _JETP Lett._ 33 (1981) 532.
  * [8] A.A. Starobinsky, _Dynamics of Phase Transition in the New Inflationary Universe Scenario and Generation of Perturbations_ , _Phys. Lett. B_ 117 (1982) 175.
  * [9] S.W. Hawking, _The Development of Irregularities in a Single Bubble Inflationary Universe_ , _Phys. Lett. B_ 115 (1982) 295.
  * [10] A.H. Guth and S.Y. Pi, _Fluctuations in the New Inflationary Universe_ , _Phys. Rev. Lett._ 49 (1982) 1110.
  * [11] Planck collaboration, _Planck 2018 results. IX. Constraints on primordial non-Gaussianity_ , _Astron. Astrophys._ 641 (2020) A9 [1905.05697].
  * [12] J.M. Maldacena, _Non-Gaussian features of primordial fluctuations in single field inflationary models_ , _JHEP_ 05 (2003) 013 [astro-ph/0210603].
  * [13] C. Armendariz-Picon, T. Damour and V.F. Mukhanov, _k - inflation_ , _Phys. Lett. B_ 458 (1999) 209 [hep-th/9904075].
  * [14] J. Garriga and V.F. Mukhanov, _Perturbations in k-inflation_ , _Phys. Lett. B_ 458 (1999) 219 [hep-th/9904176].
  * [15] T. Kobayashi, M. Yamaguchi and J. Yokoyama, _G-inflation: Inflation driven by the Galileon field_ , _Phys. Rev. Lett._ 105 (2010) 231302 [1008.0603].
  * [16] T. Kobayashi, M. Yamaguchi and J. Yokoyama, _Generalized G-inflation: Inflation with the most general second-order field equations_ , _Prog. Theor. Phys._ 126 (2011) 511 [1105.5723].
  * [17] N. Arkani-Hamed, H.-C. Cheng, M.A. Luty and S. Mukohyama, _Ghost condensation and a consistent infrared modification of gravity_ , _JHEP_ 05 (2004) 074 [hep-th/0312099].
  * [18] N. Arkani-Hamed, P. Creminelli, S. Mukohyama and M. Zaldarriaga, _Ghost inflation_ , _JCAP_ 04 (2004) 001 [hep-th/0312100].
  * [19] M. Alishahiha, E. Silverstein and D. Tong, _DBI in the sky_ , _Phys. Rev. D_ 70 (2004) 123505 [hep-th/0404084].
  * [20] S. Weinberg, _Quantum contributions to cosmological correlations_ , _Phys. Rev. D_ 72 (2005) 043514 [hep-th/0506236].
  * [21] S. Weinberg, _Quantum contributions to cosmological correlations. II. Can these corrections become large?_ , _Phys. Rev. D_ 74 (2006) 023508 [hep-th/0605244].
  * [22] L. Senatore and M. Zaldarriaga, _On Loops in Inflation_ , _JHEP_ 12 (2010) 008 [0912.2734].
  * [23] N. Bartolo, E. Dimastrogiovanni and A. Vallinotto, _One-Loop Corrections to the Power Spectrum in General Single-Field Inflation_ , _JCAP_ 11 (2010) 003 [1006.0196].
  * [24] S. Melville and E. Pajer, _Cosmological Cutting Rules_ , _JHEP_ 05 (2021) 249 [2103.09832].
  * [25] H. Goodhew, S. Jazayeri and E. Pajer, _The Cosmological Optical Theorem_ , _JCAP_ 04 (2021) 021 [2009.02898].
  * [26] J. Kristiano and J. Yokoyama, _Why Must Primordial Non-Gaussianity Be Very Small?_ , _Phys. Rev. Lett._ 128 (2022) 061301 [2104.01953].
  * [27] C. Cheung, P. Creminelli, A.L. Fitzpatrick, J. Kaplan and L. Senatore, _The Effective Field Theory of Inflation_ , _JHEP_ 03 (2008) 014 [0709.0293].
  * [28] M. Sasaki, _Large Scale Quantum Fluctuations in the Inflationary Universe_ , _Prog. Theor. Phys._ 76 (1986) 1036.
  * [29] X. Chen, B. Hu, M.-x. Huang, G. Shiu and Y. Wang, _Large Primordial Trispectra in General Single Field Inflation_ , _JCAP_ 08 (2009) 008 [0905.3494].
  * [30] X. Chen, H. Firouzjahi, M.H. Namjoo and M. Sasaki, _A Single Field Inflation Model with Large Local Non-Gaussianity_ , _EPL_ 102 (2013) 59001 [1301.5699].
  * [31] D. Seery and J.E. Lidsey, _Primordial non-Gaussianities in single field inflation_ , _JCAP_ 06 (2005) 003 [astro-ph/0503692].
  * [32] P. Adshead, R. Easther and E.A. Lim, _Cosmology With Many Light Scalar Fields: Stochastic Inflation and Loop Corrections_ , _Phys. Rev. D_ 79 (2009) 063504 [0809.4008].
  * [33] D. Green and A. Premkumar, _Dynamical RG and Critical Phenomena in de Sitter Space_ , _JHEP_ 04 (2020) 064 [2001.05974].
  * [34] T. Cohen and D. Green, _Soft de Sitter Effective Theory_ , _JHEP_ 12 (2020) 041 [2007.03693].
  * [35] C.P. Burgess, L. Leblond, R. Holman and S. Shandera, _Super-Hubble de Sitter Fluctuations and the Dynamical RG_ , _JCAP_ 03 (2010) 033 [0912.1608].