# Non-Gaussianities in the Extended EFT of Inflation

Amjad Ashoorioon  Ghazal Geshnizjani  Hyung J. Kim

###### Abstract

In earlier works, we studied the validity of Extended Effective Field Theory
of Inflation (EEFToI) in the regime where initial conditions are set with
dispersion relations $\omega^{2}\propto k^{6}$. We had also evaluated and
examined the power spectrum for some interesting corners of the parameter
space. In this paper, we compute the bispectrum in the EEFToI, take a closer
look at the strong coupling constraints and calculate the size of the non-
Gaussianities in those regions of parameter space. We also investigate the
shape of triangles that contribute to the enhancement of non-Gaussianities in
this regime. We find that there are allowed parts of parameter spaces where
EEFToI description with initial conditions set with $\omega^{2}\propto k^{6}$
is sensible and interesting.

##  1 Introduction

The inflationary cosmology is one of the most widely accepted paradigms of the
early universe cosmology. An early phase of quasi-de Sitter expansion in which
the universe accelerated, was developed in order to address some existing
issues of the Standard Big Bang (SBB) Cosmology, such as the flatness and the
horizon problem[1, 2]. Soon it was also recognised that it can also provide a
framework for generating the seeds of large scale structures (LSS) from sub-
horizon vacuum quantum fluctuations [3]. However, with all its successes it is
important to note Inflation is not just one specific model manifested from a
fundamental theory of particle physics and gravity. In fact, as we do not yet
have a well understood theory of quantum gravity, every inflationary model
should be viewed as an Effective Field Theory in the semi-classical gravity
regime. Therefore, it is important to make sure the effective field theory
description is self consistent in its regime of its validity [4, 5].
Meanwhile, as a healthy scientific exercise we should also explore if there
are other sensible alternative scenarios [6, 7] to inflationary paradigm that
can describe the initial conditions of the universe. This is specially
relevant in the light of some of the less desirable outcomes of the inflation
such as eternally inflationary spacetimes[1, 2] or being geodesically past
incomplete [8]. The latter implies inflation does not address the singularity
problem [9].

As for probing Inflation through a consistent effective field theory
framework, which is the topic of this paper, in 2006 Creminelli et al. [10]
introduced a systematic way to write the effective field theory for the
cosmological perturbations of single field models without explicitly
specifying the details of the action of the matter field. In this formalism,
coined as the Effective Field Theory of Inflation (EFToI), the action for
curvature perturbations is constructed by requiring the original theory to be
four dimensional diffeomorphism invariant with only one dynamical scalar
degree of freedom (inflaton field). The action is constructed in the unitary
gauge, i.e. the preferred spatial slicing of the inflaton, where it is taken
to be the clock. In this gauge, the time diffeomorphism is spontaneously
broken and is no longer a linearly realized symmetry. Due to the reduced
symmetry in this gauge, there can be more terms in the action as all the 3D
spatial diffeomorphism invariant terms can also be included. This action is
further simplified by expanding it perturbatively around a quasi-de Sitter
Friedman-Robertson-Walker (FRW) background and in the decoupling limit, where
mixing of the scalar mode with gravity can be ignored. The detailed derivation
of the action will be reviewed later in the next section. One of the strengths
of this approach is that it can unite and reproduce all the other single field
inflationary models that had previously been encountered such as the canonical
single field model, ghost inflation, and DBI inflation, by adjusting the
parameters under one framework. In particular, it allows exploring some of the
non-standard features of single field models, such as the deviations of the
speed of propagation from speed of light or modified dispersion relations in
light of the newly introduced operators in unitary gauge.

Having noticed the potential of describing higher order dispersion relations
within EFToI, the very early works of EFToI [11, 10, 12], explored their
implications for inflation up to the order of $\omega^{2}\propto k^{4}$.
Nevertheless, terms leading to $\omega^{2}\propto k^{n}$ with $n\geq 6$ were
not studied, due to the speculations of strong coupling at low energy limit.
However, as pointed out in our papers [13, 14, 15] using a power counting
argument for all the the relevant operators, EFT does not necessarily face any
strong coupling problem in the low energy as it naturally renders back to the
$k^{4}$ and $k^{2}$ regimes, when the mode gets stretched due to the universe
expansion. Therefore, we introduced the notion of Extended Effective Field
Theory of Inflation (EEFToI) where operators that could produce sixth order
$k$ terms or higher orders can also be included in the action. Our focus in
those works was primarily on the implications of such models for the two point
function and we left out the rigorous calculation of the three point function
for later. In particular, we examined the parameter space of the modified
dispersion relations of EEFToI, where adiabatic conditions were mildly
violated but still under control, and computed the enhancement to the
amplitude of the curvature perturbations compared to the cases initiating from
the Bunch Davies initial conditions. This analysis allowed us to find some
regions in the parameter space that can lead to one to two orders of magnitude
enhancement of the scalar power spectrum which in return imply lower energy
scales for inflation itself 111This framework can also be further pushed into
the regions of the parameter space where the power spectrum gets modulated by
many more orders of magnitude and $\Delta{{}_{\zeta}}\sim 10^{-1}$ at the
relevant scales, to explain the origin of enhancement in the scalar power
spectrum needed to explain the formation of primordial black holes (PBHs)
[16]. There it is argued that the effective field theory criterion of
avoidance of the strong coupling could be satisfied too, as the modes for such
dispersion relations are at most dominated by the quartic part of the
dispersion relation at horizon crossing. However, in this work, we mainly
focused on the regions of the parameter space that was discussed in [14] where
the enhancement is at most about an order of magnitude from the usual
Lorentzian dispersion relation with the Bunch-Davies vacuum.. Having done that
study, we now return to the analysis of the three point function and the
strong coupling constraints, to obtain a more clear picture of the tuning
required for the coefficients of the higher dimensional operators in order to
have a sensible EFT description with initial conditions set via
$\omega^{2}\propto k^{6}$ at early times. In fact since the observation of the
fluctuations in the cosmic microwave background (CMB) and the measurements of
$f_{{}_{\rm NL}}$ parameters impose a much tighter constraints on interaction
terms in the action, we use them to further scrutinise the effects of the
modifications to dispersion relations 222For study of non-Gaussianity in the
EFToI refer to [12, 11]. Furthermore, if the primordial non-Gaussianities are
eventually observed, calculating the non-Gaussianity signatures of the theory
allows us to distinguish between different scenarios of inflation.

The outline of the paper is as follows. In section (2), we will briefly review
EFToI and its extension to EEFToI. In section (3), we will examine the
interaction terms and calculate the bispectrum. Section (4) summarises our
numerical result for non-Gausianities and constraints on parameter space. We
will make our concluding remarks in section (5).

##  2 Review of Effective Field Theory and Extended Effective Field Theory of
Inflation

In the context of inflationary cosmology, the basic idea of the single field
EFToI is that even though the action is invariant under all diffeomorphisms
(temporal and spatial), the cosmological solution has a preferred time
slicing, in which the fluctuations of the inflaton field vanish. In this
gauge, which is known as the unitary gauge, the surviving symmetries of the
action are only the spatial diffeomorphisms. In the unitary gauge, the
dynamical scalar degree of freedom which was sourced by the scalar field, is
manifested in the longitudinal component of the metric through the Goldstone
boson corresponding to the broken time symmetry. As expected this Goldstone
boson transforms non-linearly under the time diffeomorphism and one can
restore the time diffeomorphism back in the action, using the Stueckelberg
trick. This property provides a framework to write down the most general
theory of fluctuations around a time-evolving background where the time diffs
are non-linearly realized.

To start, one writes the most general action in unitary gauge using all the
possible operators that satisfy spatial diffeomorphisms. These operators, in
addition to the 4D diff invariant terms such as Reimann tensor and its
covariant derivatives, can also include terms that are pure functions of time,
extrinsic curvature terms $K_{\mu\nu}$ and $g^{00}$ operators. It has been
shown in details that using this formalism, the Lagrangian of single scalar
field inflationary models around flat Friedmann-Lemaître-Robertson-Walker
(FLRW) metric can be written as [11, 14, 15, 15]

$\mathcal{L}=M_{{Pl}}^{2}\left[\frac{1}{2}\,R+\dot{H}\,g^{00}-\left(3\,H^{2}+\dot{H}\right)\right]+\sum_{m\geq
2}\,\mathcal{L}_{m}(g^{00}+1,\delta K_{\mu\nu},\delta
R_{\mu\nu\rho\sigma};\nabla_{\mu};t)\ ,$ (2.1)

where $\mathcal{L}_{m}$ represents the spatially invariant terms of
perturbative order $m$ in arguments $g^{00}+1,~{}\delta K_{\mu\nu}$ and
$\delta R_{\mu\nu\rho\sigma}$. This action can in principle include many terms
each describing perturbations in a different inflationary model. However, in
the original derivation of EFToI [11], due to the reasons that will be
explained soon, only the terms of mass dimensions up to 2 were considered in
the sum. Later, in the EEFToI [14, 15], the terms up to $[{\rm mass}]^{4}$
were also examined and one could even argue more terms can be included too.
Once a selection of terms are written down, then the Stückelberg trick is used
to make the Goldstone boson manifest. For that, first the time transformation
$t\rightarrow t+\xi^{0}({x^{\mu}})$ is taken and then time diffeomorphisms was
restored by substituting $\xi^{0}(x^{\mu})$ with a field $-\pi(x^{\mu})$. Now
asserting that under time diff $\pi\rightarrow\pi-\xi^{0}({x^{\mu}})$, action
becomes invariant under all diffs at all orders. Using this procedure, one can
retrieve the Lagrangian density of the canonical single field slow-roll
inflation from the EFToI by taking the terms of perturbative order of $m=0$.
That entails, expanding $m=0$ terms, perturbatively in $\pi$ and taking the
decoupling limit where $\dot{H}\rightarrow 0$ while $M_{{}_{\rm
Pl}}^{2}\dot{H}$ is kept fixed, leading to the following Lagrangian density
for $\pi$ 333where we have chosen $(-,+,+,+)$ convention,

$\mathcal{L}_{{}_{\rm slow-roll}}=-M_{{}_{\rm
Pl}}^{2}\,\dot{H}\left(\dot{\pi}^{2}-{\frac{(\partial_{i}\pi)^{2}}{a^{2}}}\right)\,.$
(2.2)

One can substitute $\pi$ with the more familiar conserved quantity
$\zeta=-H\pi$ to obtain the standard slow-roll inflationary action for
$\zeta$.

Next, including more terms with $m\geq 2$ and operators up to mass order two
($d=2$) a larger class of inflationary scenarios were identified and studied
[11],

$\displaystyle\mathcal{L}_{{}_{\rm EFToI}}$ $\displaystyle=$ $\displaystyle
M_{{}_{\rm
Pl}}^{2}\left[\frac{1}{2}\,R+\dot{H}\,g^{00}-\left(3\,H^{2}+\dot{H}\right)\right]$
(2.3)
$\displaystyle+\frac{M_{2}^{4}}{2!}\,(g^{00}+1)^{2}+\frac{M_{3}^{4}}{3!}\,(g^{00}+1)^{3}+\frac{\bar{M}_{1}^{3}}{2}(g^{00}+1)\delta
K^{\mu}_{\ \mu}-\frac{\bar{M}_{2}^{2}}{2}\,(\delta K^{\mu}_{\ \mu})^{2}$
$\displaystyle-\frac{\bar{M}_{3}^{2}}{2}\,\delta K^{\mu}_{\ \nu}\,\delta
K^{\nu}_{\ \mu}\ .$

Here, the upper indices are exponents of $M^{i}$ and $\bar{M}^{i}$ and are
adjusted to reflect the mass dimension of the coefficients. Note that the
Lagrangian density above is not an exhaustive list of all the possible
operators with mass dimensions $d\leq 2$. Since $(g^{00}+1)$ is dimension
zero, multiples of the terms listed above and powers of $(g^{00}+1)$ will not
change their mass dimension. Since $g^{00}+1\sim 2\dot{\pi}$, these new terms
will not contribute to powerspectrum but may change bispectrum or higher order
correlation functions. There is also another mass dimension squared operator,
$\bar{M}^{2}\nabla^{\mu}g^{00}\nabla_{\mu}g^{00}$ which after the
implementation of the Stueckelberg technique yields Ostragradski ghosts at
second order for the Goldstone boson, $\pi$. That operator in the unitary
gauge action leads to

$\bar{M}^{2}\left(-\ddot{\pi}^{2}+\frac{k^{2}}{a^{2}}\dot{\pi}^{2}\right)\,.$
(2.4)

in the Lagrangian for the Goldstone boson. As we will see, such ghosts are
generated in the Extended EFT of inflation too which limits the size and value
of their corresponding coupling constants 444The second term,
$\frac{k^{2}}{a^{2}}\dot{\pi}^{2}$ leads to time-dependent modification of the
speed of sound, making it vanish at infinite past, when the wavelength of the
modes are small. This will lead to reducing the cutoff below the Hubble
parameter at some point in the evolution of the modes and lead to strong
coupling, please see [15] for a similar situation that arises in the context
of Extended EFT of inflation..

The Lagrangian of the EFToI (2.3) in the decoupling limit, where
$\dot{H}\rightarrow 0$ and $M_{{}_{\rm Pl}}^{2}\dot{H}$ is kept fixed, results
into the following quadratic and cubic expressions for $\pi$ [11, 17]:

$\displaystyle\mathcal{L}_{{}_{\rm EFToI}}^{(\pi)}$ $\displaystyle=$
$\displaystyle-M_{{}_{\rm
Pl}}^{2}\,\dot{H}\left(\dot{\pi}^{2}-{\frac{(\partial_{i}\pi)^{2}}{a^{2}}}\right)+2M_{2}^{4}\left(\dot{\pi}^{2}+\dot{\pi}^{3}-\dot{\pi}\frac{(\partial_{i}\pi)^{2}}{2a^{2}}\right)-\frac{4}{3}M_{3}^{4}\dot{\pi}^{3}-\bar{M}_{1}^{3}\bigg{(}3H\dot{\pi}^{2}$
$\displaystyle-$ $\displaystyle
H\frac{(\partial_{i}\pi)^{2}}{2a^{2}}-\frac{(\partial_{j}^{2}\pi)(\partial_{i}\pi)^{2}}{2a^{4}}\bigg{)}-\frac{\bar{M}_{2}^{2}}{2}\bigg{(}9H^{2}\dot{\pi}^{2}-3H^{2}\frac{(\partial_{i}\pi)^{2}}{a^{2}}+\frac{(\partial_{i}^{2}\pi)^{2}}{a^{4}}+\frac{H(\partial_{i}^{2}\pi)(\partial_{j}\pi)^{2}}{a^{4}}$
$\displaystyle+$
$\displaystyle\frac{2\dot{\pi}\partial_{i}^{2}\partial_{j}\pi\partial_{j}\pi}{a^{4}}\bigg{)}-\frac{\bar{M}_{3}^{2}}{2}\bigg{(}3H^{2}\dot{\pi}^{2}-H^{2}\frac{(\partial_{i}\pi)^{2}}{a^{2}}+\frac{(\partial_{i}^{2}\pi)^{2}}{a^{4}}+\frac{2H(\partial_{i}^{2}\pi)(\partial_{j}\pi)^{2}}{a^{4}}$
$\displaystyle+$
$\displaystyle\frac{2\dot{\pi}\partial_{i}^{2}\partial_{j}\pi\partial_{j}\pi}{a^{4}}\bigg{)}+\dots\
.$

First, note that any term that contributes to the coefficient of
$\dot{\pi}^{2}$ in the approximate plain wave regime, leads to deviations from
$c_{{}_{S}}=1$. However, any terms that results in spatial derivatives for
$\pi$, can potentially modify the dispersion relation. Some modifications
could even lead to non-hyperbolic equations of motion for $\pi$ that do not
have propagating solutions or are not even well posed. As we see, the first
term in above expression, $\frac{M_{2}^{4}}{2!}\,(g^{00}+1)^{2}$, only
contributes to modifying $c_{{}_{S}}$ [11],

$c_{{}_{S}}^{-2}=1-\frac{2M_{2}^{4}}{M_{\rm Pl}^{2}\dot{H}}\,.$ (2.6)

This shows that including only this term, the superluminal propagation can be
avoided for $M_{2}^{4}>0$, since during inflation $\dot{H}<0$ 555There are
some arguments that superluminal propagation could lead to the prohibition of
UV completion [18].. Through a quick substitution for $M_{2}$, we see the
cubic term in $\pi$ from this same operator, appears with a factor
proportional to $(1-\frac{1}{c_{{}_{S}}^{2}})$. This agrees with the generic
feature that in the single field inflationary models, smaller values of the
speed of sound result in the enhancement of the non-gaussianities. This in
return, through observational bounds, further constraints the value of
$M_{2}$.

Next, the term proportional to $\bar{M}_{1}$ in (2.3) is odd under time
reversal and if we assume theory has time reversal symmetry we could set it to
zero. Nevertheless, it also preserves the form of dispersion relation as
$\omega^{2}=c_{{}_{S}}^{2}k^{2}$ but contributes to changes in speed of sound.
In the exact de-Sitter limit, $\dot{H}=0$, with $M_{2}\neq 0$ and
$\bar{M}_{1}\sim-M_{2}\sim\bar{M}$, the speed of sound turns out to be quite
small, $c_{{}_{S}}^{2}=\frac{H}{4\bar{M}}\ll 1$, but if $M_{2}=0,~{}M_{3}<0$,
then $c_{{}_{S}}^{2}=\frac{1}{6}$.

Finally, the last two terms would not only add corrections to the sound speed
for $\pi$ field, but also induce corrections proportional to $k^{4}$ to the
dispersion relation. The effect of such quartic corrections to the dispersion
relation on the power spectrum and bispectrum has been worked out both
analytically and numerically [19, 17].

While in the original work of [11], it was already noticed that including
other possible operators in (2.1) in the action, could lead to dispersion
relations beyond the quartic form, it was argued that those possibilities
would lead to the violation of the EFT criterion. The rationale of the
argument is that, assuming $\omega^{2}\propto k^{2n}$ with $n\geq 3$, the
power of the energy scaling of $\pi$ would be $-1/2+3/2n$ which means the
scaling dimension of the cubic operator $\dot{\pi}(\partial_{i}\pi)^{2}$ would
be $(7-3n)/(2n)$. For $n\geq 3$, i.e. beyond quartic form, this scaling
dimension becomes negative which leads to this operator becoming relevant at
low energies and hence the effective field theory regime breaks down. However,
as we argued in [14] this statement would be true if $\omega^{2}\propto
k^{2n}$ dispersion relation remained valid up to horizon crossing and there
were no other intermediate EFT regime that dominated the dispersion relation
at low energies. We pointed out first, interestingly there are always finite
number of relevant operators for every given $n$. Second, the operators
$(\partial_{i}^{2}\pi)^{2}$ and $(\partial_{i}\pi)^{2}$ that change the
dispersion relation are among relevant operators with negative scaling but
larger in magnitude powers. Therefore, they can make the theory transition
into a different regime of dispersion relation before it becomes strongly
coupled. For example, in the case of $\omega^{2}\sim k^{6}$, the scaling power
of $\dot{\pi}(\partial_{i}\pi)^{2}$ is $-1/3$, but the scaling power of
$(\partial_{i}^{2}\pi)^{2}$ and $(\partial_{i}\pi)^{2}$ are respectively
$-2/3$ and $-4/3$. This implies that the dispersion relation would change from
$\omega^{2}\sim k^{6}$ to $\omega^{2}\sim k^{4}$ or $\omega^{2}\sim k^{2}$ at
some energy scale $\Lambda_{\rm dis}$. Hence, if the coefficients of the
higher dimensional operators are such that different cutoffs of each EFT
regime follow a hierarchy of successive scales, theory remains weakly coupled
throughout its evolution down to low energies[14]. This in principle can allow
for an interesting extension of the EFT regime to $\omega^{2}\propto k^{6}$.

In the case of EFT of inflation, the extrinsic curvature squared terms in the
action not only induced corrections proportional to $k^{4}$, but they induce
time dependent corrections to $k^{2}$. Therefore, pure $k^{4}$ dispersion
relation can come about only through fine-tuning the coefficients of the
operators in the unitary gauge action. Likewise, addition of higher mass
dimension operators in the unitary gauge action, not only produces corrections
beyond the quartic order to the dispersion relation, they will at the same
time induce lower order corrections to the dispersion relation which will
become important at the low energy scales. Hence, it seems that the theory
benefits from some self-healing properties.

Motivated by this reasoning, in the EEFToI papers [13, 14], we extended the
EFToI formalism to study the consequences of including terms up to order $d=4$
that can further change the dispersion relation in addition to the terms that
had previously been studied in EFToI:

$\mathcal{L}_{{}_{\rm EEFToI}}=\mathcal{L}_{{}_{\rm
EFToI}}+\sum_{d=3}^{4}\mathcal{L}_{2,d}$ (2.7)

The additional terms up to order $d=4$ that were included in EEFToI are given
by

$\displaystyle\sum_{d=3}^{4}\mathcal{L}_{2,d}$ $\displaystyle=$
$\displaystyle\frac{\bar{M}_{4}}{2}\nabla^{\mu}g^{00}\nabla^{\nu}\delta
K_{\mu\nu}-\frac{\delta_{1}}{2}\,(\nabla_{\mu}\delta
K^{\nu\gamma})(\nabla^{\mu}\delta
K_{\nu\gamma})-\frac{\delta_{2}}{2}\,(\nabla_{\mu}\delta K^{\nu}_{\ \nu})^{2}$
(2.8) $\displaystyle-\frac{\delta_{3}}{2}\,(\nabla_{\mu}\delta K^{\mu}_{\
\nu})(\nabla_{\gamma}\delta
K^{\gamma\nu})-\frac{\delta_{4}}{2}\,\nabla^{\mu}\delta
K_{\nu\mu}\nabla^{\nu}\delta K_{\sigma}^{\sigma}.$

Here, $\delta_{i}$ are dimensionless coefficients.

The new terms included in the Extended EFToI (2.8) for mass dimension 3 and 4,
lead to the following corrections to the quadratic action in Foureir space:

$\displaystyle\mathcal{L}_{2,3}^{(\pi_{k})}=\frac{\bar{M}_{4}}{2}\left(\frac{k^{4}H\pi_{k}^{2}}{a^{4}}+\frac{k^{2}H^{3}\pi_{k}^{2}}{a^{2}}-9H^{3}\dot{\pi}_{k}^{2}\right)\,,$
(2.9)

and

$\displaystyle\mathcal{L}_{2,4}^{(\pi_{k})}$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\delta_{1}\left(\frac{k^{6}\pi_{k}^{2}}{a^{6}}-\frac{3H^{2}k^{4}\pi_{k}^{2}}{a^{4}}-\frac{k^{4}\dot{\pi}_{k}^{2}}{a^{4}}+\frac{4H^{4}k^{2}\pi_{k}^{2}}{a^{2}}-6H^{4}\dot{\pi}_{k}^{2}-3H^{2}\ddot{\pi_{k}}^{2}\right)$
$\displaystyle-\frac{1}{2}\delta_{2}\left(\frac{k^{6}\pi_{k}^{2}}{a^{6}}+\frac{H^{2}k^{4}\pi_{k}^{2}}{a^{4}}-\frac{k^{4}\dot{\pi}_{k}^{2}}{a^{4}}+\frac{6H^{4}k^{2}\pi_{k}^{2}}{a^{2}}-9H^{2}\ddot{\pi_{k}}^{2}\right)$
$\displaystyle-\frac{1}{2}\delta_{3}\left(\frac{k^{6}\pi_{k}^{2}}{a^{6}}+\frac{3H^{2}k^{4}\pi_{k}^{2}}{a^{4}}+\frac{H^{2}k^{2}{\dot{\pi}_{k}}^{2}}{a^{2}}-9H^{4}\dot{\pi}_{k}^{2}\right)$
$\displaystyle-\frac{1}{2}\delta_{4}\left(\frac{k^{6}\pi_{k}^{2}}{a^{6}}+\frac{H^{2}k^{4}\pi_{k}^{2}}{2a^{4}}+\frac{9H^{4}k^{2}\pi_{k}^{2}}{2a^{2}}+\frac{3H^{2}k^{2}\dot{\pi}_{k}^{2}}{a^{2}}+\frac{27}{2}H^{4}\dot{\pi}_{k}^{2}\right)\,.$

As we see $\delta_{1}$ and $\delta_{2}$ expressions contain $(\ddot{\pi})^{2}$
or Ostrogradsky ghosts. In the context of an effective field theory, ghosts
are not necessarily a pathology if they are small corrections and under
control. This is because ghosts maybe an artifact of the truncation and the
effective theory is only meant to capture the relevant features of the full
theory [20, 21]. However, in this work we will not study their consequences
and for simplicity focus on the ghost free possibilities by setting the
coefficients $\delta_{1}=\delta_{2}=0$ 666Setting $\delta_{1}+3\delta_{2}=0$,
one could in principle avoid the $\ddot{\pi}^{2}$ in the action for the
Goldstone boson, $\pi$. The term proportional to $\delta_{1}$,
$(\nabla_{\mu}\delta K^{\nu\gamma})(\nabla^{\mu}\delta K_{\nu\gamma})$, gives
rise to ghosts for the tensor perturbations too [14]. Imposing lack of ghosts
in both the scalar and tensor sectors of the theory one has to impose
$\delta_{1}=\delta_{2}=0$ 777One can show that if $\nabla^{\mu}$ operators are
contracted with the induced metric, $h^{\mu\nu}$, there will be no ghosts in
the scalar and tensor part of the theory. However, here we keep to the
original framework of the Effective Field Theory of inflation, in which
operators built out of $\nabla_{\mu}$ are allowed as they are invariant under
spatial diffeomorphisms[14]..

Also we like to point out that similar to EFToI, (2.9) is not the exhaustive
list of all $d=3$ and $d=4$ operators. To begin with, we are only including
terms that are quadratics in $\delta g^{00}$ and $\delta{K}$, as only those
terms lead to corrections for the dispersion relation. Next, in the derivation
of the above expression terms that could be absorbed into other ones through
integration by part are not included. For example, the terms
$(\nabla_{\mu}\nabla_{\nu}\delta g^{00})\delta K^{\mu\nu}$, and $\delta
g^{00}(\nabla_{\mu}\nabla_{\nu}\delta K^{\mu\nu})$ are absorbed into the
$\bar{M}_{4}$ term. Finally, there are more terms such as
$(\nabla_{\mu}\nabla^{\mu}\delta g^{00})^{2}$,
$(\nabla_{\mu}\nabla^{\mu}\delta g^{00})\text{Tr}(\delta K)$ that have been
discarded as they also result in $(\ddot{\pi})^{2}$ in the quadratic
Lagrangian for $\pi$ 888Similarly one can avoid the ghosts in these operators
by contracting the $\nabla_{\mu}$ operators with the induced metric, rather
then the four dimensional metric.. There is also the possibility that in a
particular corner of the parameter space, the coefficients of these terms are
non-zero, yet they conspire in a way that all the $(\ddot{\pi})^{2}$ terms
cancel out each other and the theory is ghost free. In such a case, assuming
there is any interesting physics, contributions to terms such as
$k^{2}\dot{\pi}_{k}^{2}$ and $k^{4}\dot{\pi}_{k}^{2}$ and the potential strong
coupling regimes need to be studied carefully [15].

Going forward we also only explore models where $\delta_{3}=-3\delta_{4}$. The
reason for this additional assumption is that we are more interested in
investigating models where the effects of $k^{6}\pi_{k}^{2}$ operators is
dominant and setting $\delta_{3}=-3\delta_{4}$ will automatically cancel off
$k^{2}\dot{\pi}_{k}^{2}$ terms. Deviating from this condition, one has to
carefully take into account the competing impact of $k^{2}\dot{\pi}_{k}^{2}$
term which can introduce a different cut off and dispersion relation in UV.
However, it can still produce a well-defined finite power spectrum under
certain criteria [15].

Under the assumptions mentioned above and after we reorganize all the terms
contributing to the Lagrangian, it will have the following form in Fourier
space:

$\int dt\,d^{3}k\sqrt{-g}\,\mathcal{L}_{{}_{\rm EEFToI}}=\int
dt\,d^{3}k\,a^{3}\left(\frac{1}{2}A_{1}\dot{\pi}_{k}^{2}+\frac{1}{2}\left[C_{1}\frac{k^{6}}{a^{6}}+D_{1}\frac{k^{4}}{a^{4}}+F_{1}\frac{k^{2}}{a^{2}}\right]\pi_{k}^{2}\right),$
(2.11)

where the the coefficients $A_{1},C_{1},D_{1}$ and $F_{1}$ are defined as

$\displaystyle A_{1}=-2M_{{}_{\rm
Pl}}^{2}\dot{H}+4M_{2}^{4}-6\bar{M}_{1}^{3}H-9H^{2}\bar{M}_{2}^{2}-3H^{2}\bar{M}_{3}^{2}+\frac{27}{2}H^{4}\delta_{3}-9H^{3}\bar{M}_{4}\,,$
$\displaystyle C_{1}=\frac{2}{3}\delta_{3}\,,$ $\displaystyle
D_{1}=\bar{M}_{2}^{2}+\bar{M}_{3}^{2}+\frac{17}{18}H^{2}\delta_{3}-\bar{M}_{4}H\,,$
$\displaystyle F_{1}=-2M_{{}_{\rm
Pl}}^{2}\dot{H}-\bar{M}_{1}^{3}H-3H^{2}\bar{M}_{2}^{2}-\bar{M}_{3}^{2}H^{2}+\frac{1}{2}H^{4}\delta_{3}-\bar{M}_{4}H^{3}\,.$
(2.12)

Next, expanding the Canonical variable $a\sqrt{A_{1}}\pi_{k}$ in terms of the
corresponding mode functions $u_{k}$, and computing the equation of motion for
$u_{k}$ yields,

$\displaystyle
u_{k}^{\prime\prime}+\left(\frac{C_{1}}{A_{1}}\frac{k^{6}}{a^{4}}+\frac{D_{1}}{A_{1}}\frac{k^{4}}{a^{2}}+\frac{F_{1}}{A_{1}}k^{2}-\frac{a^{\prime\prime}}{a}\right)u_{k}=0\,,$
(2.13)

where the primes denote derivatives with respect to the conformal time $\tau$
such that $f^{\prime}=\frac{d}{d\tau}f=a\frac{d}{dt}f=a\dot{f}$. Note that in
this equation for the positive slow varying values of $F_{1}/A_{1}$, in the
regime where $\frac{F_{1}}{A_{1}}k^{2}$ is the dominant term ($\omega\propto
k$) the speed of sound for scalar perturbations is well defined and real, with

$c_{{}_{S}}^{2}=\frac{F_{1}}{A_{1}}.$ (2.14)

In the de Sitter limit where $a=-1/(H\tau)$ and defining $y=c_{{}_{S}}k\tau$
and $\tilde{u}(y)=\sqrt{c_{{}_{S}}k}\,u_{k}(\tau)$ to preserve the
normalization for Wronskian condition, we can rewrite equation(2.13) as

$\frac{d^{2}\tilde{u}(y)}{dy^{2}}+\left(1+\alpha_{0}y^{2}+\beta_{0}y^{4}-\frac{2}{y^{2}}\right)\tilde{u}(y)=0\,,$
(2.15)

where =

$\displaystyle\alpha_{0}=\frac{D_{1}A_{1}H^{2}}{F_{1}^{2}},$
$\displaystyle\text{ \hskip 14.22636ptand \hskip
14.22636pt}\beta_{0}=\frac{C_{1}A_{1}^{2}H^{4}}{F_{1}^{3}}.$ (2.16)

We can solve the mode equation (2.15) numerically, taking the positive
frequency WKB mode in the infinite past,

$\tilde{u}(y\rightarrow-\infty)\simeq\frac{1}{2}\left(-\frac{\pi
y}{3}\right)^{1/2}H^{(1)}_{\frac{1}{6}}\left(-\sqrt{\beta_{0}}\frac{y^{3}}{3}\right).$
(2.17)

Noting that curvature perturbations $\zeta_{k}$ are related to $\tilde{u}$ via
$\zeta_{k}(\tau)=-\frac{H}{a\sqrt{A_{1}c_{{}_{S}}k}}\tilde{u}(y)$, the effect
of sixth order polynomial dispersion relation of EEFToI on the scalar power
spectrum can be obtained. This result can be represented as a modulation
factor, $\gamma_{{}_{S}}$, on the power spectrum expected from the Bunch-
Davies vacuum $P_{\zeta}^{B.D.}$,

$P_{{}_{\zeta}}=\gamma_{{}_{S}}P_{{}_{\zeta}}^{B.D.}\,.$ (2.18)

Since the dimensionless power spectrum is related to $P_{{}_{\zeta}}^{B.D}$ as

$\displaystyle\Delta^{B.D.}_{\zeta}\equiv\frac{k^{3}}{2\pi^{2}}P_{\zeta}^{B.D.}=\frac{1}{4\pi^{2}}\frac{H^{4}}{A_{1}c_{{}_{S}}},$
(2.19)

we also get

$\displaystyle\Delta_{{}_{\zeta}}=\frac{\gamma_{{}_{S}}}{4\pi^{2}}\frac{H^{4}}{A_{1}c_{{}_{S}}}.$
(2.20)

As can be expected, value of $\gamma_{{}_{S}}$ depends on values of
$\alpha_{0}$ and $\beta_{0}$. Under the assumption that the dispersion
relation does not become tachyonic before horizon crossing (i.e.
$1+\alpha_{0}y^{2}+\beta_{0}y^{4}>0$), it was shown in [13] that the greatest
enhancements to the Bunch-Davies power spectrum is obtained when
$\alpha_{0}<0$. In particular, with the assumption that the mode never becomes
tachyonic while it is inside the horizon, we obtained the maximum enhancement
for $\alpha_{0}=-0.18$ and $\beta_{0}=\alpha_{0}^{2}/4$. For $\alpha_{0}<0$
and $\beta_{0}\leq\alpha_{0}^{2}/4$, the dispersion relation can become
tachyonic, allowing for larger enhancements of the power spectrum [14, 16].
For $\alpha_{0},~{}\beta_{0}\geq 0$, the scalar power spectrum get suppressed
with respect to its Bunch-Davies counterpart, $\gamma_{{}_{S}}\lesssim 1$.
This observation along with the fact that for pure quartic positive correction
to the dispersion relation, $\gamma_{{}_{S}}<1$ too [19], drives us toward
this conclusion that all positive correction to the dispersion relation has
the effect of diminishing the power spectrum, For small values of
$0\leq\alpha_{0},~{}\beta_{0}\ll 1$, $\gamma_{{}_{S}}\approx 1$. This agrees
with our intuition that if WKB approximation is not violated and the mode
changes adiabatically from one regime to another, we do not expect modes to
get very excited and they should remain close to the adiabatic vacuum. Figure
(1) illustrates some of these effects for different values of $\alpha_{0}$ and
$\beta_{0}$.

Now that we have reviewed the main findings of EEFToI at the leve of power
spectrum, we can move on to the main goal of this paper, which is computing
the third order terms of EEFToI action, bispectrum and finally the non-
Gaussianity signatures or additional constraints for EEFToI models.

(a) $\gamma_{{}_{S}}$ as a function of $\beta_{0}/\alpha_{0}^{2}$ when
$\alpha_{0}<0$ and $\beta_{0}>0$ and $0.01\leq|\alpha_{0}|\leq 2$ and
$\frac{\alpha^{2}}{4}\leq\beta_{0}$.

(b) Plot of $\gamma_{{}_{S}}$ as a function of $\beta_{0}$ when $\alpha_{0}$
and $\beta_{0}$ are both positive and $1\leq\alpha_{0}\leq 2$ with
$0\leq\beta_{0}\leq 2$.

Figure 1: Plot of the amplification of power spectrum compared to the Bunch
Davies power spectrum $P_{s}^{B.S}$ with varying values of $\alpha_{0}$ and
$\beta_{0}$.

##  3 Calculation of the Bispectrum and Estimating the non-Gaussianity

In the previous section, we reviewed the impact of some interesting regions of
parameter space in the EEFToI on the power spectrum, where the perturbations
start with modified dispersion relation $\omega^{2}\propto k^{6}$. The
computation of bispectrum and thus the estimation of non-Gaussianity for
EEFToI can help us explore more distinctive signatures of these models or
further constrain the parameter space of the EEFToI theory. In Fourier space,
the bispectrum is related to the vacuum expectation value of the three point
function for $\zeta_{\textbf{k}}$ as

$\langle\zeta_{\textbf{k}_{1}}(t)\zeta_{\textbf{k}_{2}}(t)\zeta_{\textbf{k}_{3}}(t)\rangle=(2\pi)^{3}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})B(k_{1},k_{2},k_{3})\,,$
(3.1)

and the three point function in the interaction picture can be computed as

$\langle\zeta_{\textbf{k}_{1}}(t)\zeta_{\textbf{k}_{2}}(t)\zeta_{\textbf{k}_{3}}(t)\rangle=-i\int^{t}_{t_{in}}dt^{\prime}\;\langle[\zeta_{\textbf{k}_{1}}(t)\zeta_{\textbf{k}_{2}}(t)\zeta_{\textbf{k}_{3}}(t),H_{int}(t^{\prime})]\rangle\,,$
(3.2)

where in these calculations, $\pi_{\textbf{k}}$ can be substituted for
$\zeta_{\textbf{k}}$ using the relation
$\zeta_{\textbf{k}}=-H\pi_{\textbf{k}}$ and $t_{in}$ is the initial time. The
interaction Hamiltonian is given by $H_{int}=-L_{int}$ which includes terms
beyond quadratic order in the action[11].

Before we do a robust computation of the three point functions, let us perform
a simple analysis to see how the interaction terms in the action act under the
energy scaling when the dispersion relation is governed by $\omega^{2}\propto
k^{6}$. The kinetic term in the action (2.11) in real space has the following
form for $\pi$

$S=\int dt\;d^{3}x\;[\frac{1}{2}A_{1}\dot{\pi}^{2}].$ (3.3)

Rescaling energy by a factor of $\Sigma$, $E\rightarrow\Sigma E$, means that
$t\rightarrow\Sigma^{-1}t$, which given the dispersion relation
$\omega^{2}\propto k^{6}$, implies $x\rightarrow\Sigma^{-1/3}x$. Assuming
$A_{1}$ variation in time is negligible (due to generalized slow-roll
conditions) the dimensional scaling implies that for the kinetic term to
remain invariant, $\pi$ must remain scale-invariant as well,
$\pi\rightarrow\pi$.

We now check the scaling of a generic term in action, $\int
dt~{}d^{3}x\;\mathcal{M}_{m,s,p}$ where $m$ denotes the number of time
derivatives, $s$ number of spatial derivatives and $p$ the overall powers of
$\pi$. Since $x\rightarrow\Sigma^{-1/3}x$, the spatial derivative goes as
$\partial_{i}\rightarrow\Sigma^{1/3}\partial_{i}$ so the contribution to
action scales as

$\mathcal{M}_{m,s,p}\;dt\;dx^{3}\rightarrow\Sigma^{-2+m+s/3}\mathcal{M}_{m,s,p}\;dt\;d^{3}x.$
(3.4)

When $m+s/3<2$, this term diverges in the low energy limit of
$\Sigma\rightarrow 0$ and becomes relevant. Note the powers of $\pi$ cannot be
more than the number of derivatives, as the time diffeomorphism invariance of
the action does not allow Goldstone boson to be produced without derivatives.
This gives us the relations $p\leq s+m$. Also, $p=1$ corresponds to linear
terms, which cancel out for on-shell solutions. All possible relevant
operators thus satisfy the following constraints,

$\displaystyle 2\leq p\leq s+m\,,$ (3.5) $\displaystyle m+s/3<2\;.$ (3.6)

Furthermore, the spatial diffeomorphism invariance enforces $s$ to be even.
Taking into account all of these constraints, the only possible combinations
of $(m,s,p)$ that remain are $(0,2,2)$, $(1,2,2)$, $(0,4,2)$, $(0,4,3)$,
$(1,2,3)$ and $(0,4,4)$. This implies the number of relevant operators is
finite and excludes the existence of any relevant operators beyond quartic
terms as well999This exercise can be done for any $\omega^{2}\propto k^{2n}$
and there are always finite number of relevant operators.. Also note that
$\mathcal{M}_{1,2,2}\propto\partial\dot{\pi}\partial_{i}\pi=\frac{1}{2}\frac{\partial}{\partial
t}~{}(\partial_{i}\pi)^{2}$ term becomes
$\mathcal{M}_{0,2,2}\propto(\partial_{i}\pi)^{2}$ after integrating by parts,
which can be absorbed into the canonical kinetic term. The scaling power of
$\mathcal{M}_{0,2,2}\propto(\partial_{i}\pi)^{2}$ and
$\mathcal{M}_{0,4,2}\propto(\partial_{i}^{2}\pi)^{2}$ which modify the
dispersion relation, as we had pointed out earlier are respectively $-4/3$ and
$-2/3$. The scaling power of
$\mathcal{M}_{1,2,3}\propto\dot{\pi}(\partial_{i}\pi)^{2}$ is $-1/3$ while
both $\mathcal{M}_{0,4,2}\propto\partial_{j}^{2}\pi(\partial_{i}\pi)^{2}$ and
$\mathcal{M}_{0,4,4}\propto(\partial_{i}\pi)^{4}$ scale as $-2/3$ as well. The
term $\mathcal{M}_{0,4,4}\propto(\partial_{i}\pi)^{4}$ contributes to
trispectrum which we will not discuss in this paper101010See [22] for
calculation of trispectrum in EFToI, with at most quartic dispersion
relation.. Therefore, to calculate the strongest contributions to bispectrum
we only need to consider operators $(0,4,3)$ and $(1,2,3)$, which are

$(\partial_{i}\pi)^{2}\dot{\pi}$ (3.7)
$\partial^{2}_{j}\pi(\partial_{i}\pi)^{2}.$ (3.8)

Note that both of these terms are also produced in (LABEL:Lpi) from the terms
included in the EFToI as well:

$\frac{1}{2}M_{2}^{4}(g^{00}+1)^{2}\xrightarrow{\\{\mathcal{O}(\pi^{3})\\}}-2M_{2}^{4}\frac{\dot{\pi}(\partial_{i}\pi)^{2}}{a^{2}}$
(3.9)

and

$\displaystyle-\frac{1}{2}\bar{M}_{1}^{3}(g^{00}+1)\delta
K^{\mu}_{\mu}-\frac{\bar{M}_{2}^{2}}{2}\,(\delta K^{\mu}_{\
\mu})^{2}-\frac{\bar{M}_{3}^{2}}{2}\,\delta K^{\mu}_{\ \nu}\,\delta K^{\nu}_{\
\mu}$
$\displaystyle\xrightarrow{\\{\mathcal{O}(\pi^{3})\\}}-\frac{1}{2}\left(\bar{M}_{1}^{3}+\bar{M}_{2}^{2}H+2\bar{M}_{3}^{2}H\right)\frac{\partial_{i}^{2}\pi(\partial_{j}\pi)^{2}}{a^{4}}.$
(3.10)

These cubic terms get additional contributions in the extended theory. The
detailed calculation can be found in the appendix (A) and the corrections
after setting $\delta_{3}=-3\delta_{4}$ are the following

$\displaystyle-\frac{\delta_{3}}{2}\nabla^{\mu}\delta$ $\displaystyle
K_{\nu\mu}\nabla^{\nu}\delta{K^{\sigma}}_{\sigma}-\frac{\delta_{4}}{2}\nabla_{\mu}\delta{K^{\mu}}_{\nu}\nabla_{\gamma}\delta
K^{\gamma\nu}$
$\displaystyle\xrightarrow{\\{\mathcal{O}(\pi^{3})\\}}\left(\frac{\delta_{3}H^{3}}{3}\right)\frac{\partial_{i}^{2}\pi(\partial_{j}\pi)^{2}}{a^{4}}-\delta_{3}H^{4}\frac{\dot{\pi}(\partial_{i}\pi)^{2}}{a^{2}}$
(3.11)

and

$\frac{\bar{M}_{4}}{2}\nabla^{\mu}\delta g^{00}\nabla^{\nu}\delta
K_{\mu\nu}\xrightarrow{\\{\mathcal{O}(\pi^{3})\\}}9\bar{M}_{4}H^{3}\frac{\dot{\pi}(\partial\pi)^{2}}{a^{2}}+\frac{1}{2}\bar{M}_{4}H^{2}\frac{\partial_{i}^{2}\pi(\partial_{j}\pi)^{2}}{a^{4}}$
(3.12)

In the limit where the parameters of EFToI and EEFToI satisfy

$-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4}=0,$ (3.13)

and

$\frac{1}{3}\delta_{3}H^{3}-\frac{1}{2}\left(\bar{M}_{1}^{3}+\bar{M}_{2}^{2}H+2\bar{M}_{3}^{2}H-\bar{M}_{4}H^{2}\right)=0,$
(3.14)

these terms will not lead to any strong coupling or non-Gaussianities.
However, aside from these particular corners of parameter space, we show below
that we do not need to worry about the impact of these terms. We start by
making a crude estimation of the contribution of these terms by directly
comparing them to quadratic terms as they appear in action. Around freezing we
have $\omega\sim H$ and $k\sim H/c_{{}_{S}}$ and approximating
$\dot{\pi}\sim\omega\pi$, the ratio of contributions of the cubic term
$(\partial_{i}\pi)^{2}\dot{\pi}$ to the quadratic term can be estimated as

$\displaystyle\frac{\mathcal{L}_{\dot{\pi}(\partial\pi)^{2}}}{\mathcal{L}_{2}}$
$\displaystyle\sim$
$\displaystyle(-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4})\frac{\dot{\pi}(\partial\pi)^{2}}{\frac{A_{1}}{2}{\dot{\pi}}^{2}}\xrightarrow{\pi_{k}}2(-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4})\frac{\omega
k^{2}\pi_{k}^{3}}{A_{1}\omega^{2}\pi_{k}^{2}}$ (3.15) $\displaystyle\sim$
$\displaystyle
2(-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4})\frac{H\pi(\frac{H^{2}}{c_{{}_{S}}^{2}}\pi^{2})}{A_{1}H^{2}\pi^{2}}\sim
2(-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4})\frac{H}{A_{1}c_{{}_{S}}^{2}}\pi,$
$\displaystyle\sim$ $\displaystyle
2\frac{(-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4})}{A_{1}c_{{}_{S}}^{2}}\zeta\sim
8\pi^{2}\left(-\delta_{3}+9\frac{\bar{M}_{4}}{H}-2\frac{M_{2}^{4}}{H^{4}}\right)\frac{\Delta_{{}_{\zeta}}}{c_{{}_{S}}\gamma_{{}_{S}}}\zeta$

and similarly for $\partial^{2}_{j}\pi(\partial_{i}\pi)^{2}$, we get

$\displaystyle\frac{\mathcal{L}_{(\partial\pi)^{2}\partial^{2}\pi}}{\mathcal{L}_{2}}$
$\displaystyle\sim$
$\displaystyle\left(\frac{2}{3}\delta_{3}H^{3}+\bar{M}_{4}H^{2}-\left(\bar{M}_{1}^{3}+\bar{M}_{2}^{2}H+2\bar{M}_{3}^{2}H\right)\right)\frac{\left(\frac{H}{c_{{}_{S}}}\pi\right)^{2}\left(\frac{H}{c_{{}_{S}}}\right)^{2}\pi}{A_{1}H^{2}\pi^{2}}$
$\displaystyle\sim$
$\displaystyle\left(\frac{2}{3}\delta_{3}H^{3}+\bar{M}_{4}H^{2}-\bar{M}_{1}^{3}-\bar{M}_{2}^{2}H-2\bar{M}_{3}^{2}H\right)\frac{H^{2}}{A_{1}c_{{}_{S}}^{4}}\pi$
$\displaystyle\sim$
$\displaystyle\left(\frac{2}{3}\delta_{3}H^{3}+\bar{M}_{4}H^{2}-\bar{M}_{1}^{3}-\bar{M}_{2}^{2}H-2\bar{M}_{3}^{2}H\right)\frac{H}{A_{1}c_{{}_{S}}^{4}}\zeta$
$\displaystyle\sim$ $\displaystyle
4\pi^{2}\left(\frac{2}{3}\delta_{3}-\frac{\bar{M}_{1}^{3}}{H^{3}}-\frac{\bar{M}_{2}^{2}}{H^{2}}-2\frac{\bar{M}_{3}^{2}}{H^{2}}+\frac{\bar{M}_{4}}{H}\right)\frac{\Delta_{{}_{\zeta}}}{c_{{}_{S}}^{3}\gamma_{{}_{S}}}\zeta$

The perturbation theory breaks down when the above ratios become of order one
but given that amplitude of curvature perturbation is of order $10^{-5}$
implying $\zeta\Delta_{{}_{\zeta}}\sim 10^{-15}$, then it is reasonable to
expect that there is a large viable window of parameter space where strong
coupling is avoided. However, as we discuss later CMB observations set a much
stronger constraints on this ratios due to the limits they set on the non-
Gaussianity.

To see the restrictions on the size of non-Gaussianity from the CMB
observations, we go back to the more rigorous calculation of bispectrum and
the the contribution of the cubic terms through the interaction Hamiltonian
and equations (3.1) and (3.2). Note that the Lagrangian and the action are
written in terms of $\pi$, and we are interested in the comoving curvature
perturbation, $\zeta=-H\pi$, and its bispectrum $<\zeta^{3}>$. Therefore, we
transform the terms $H_{int}(\pi)\rightarrow H_{int}(\zeta)$ and also for
convenience switch to conformal time $\tau$ instead of $t$. As we pointed out
earlier, the most dangerous contributions from $H_{int}$ or $L_{int}$ are
coming from terms $\dot{\pi}(\partial_{i}\pi)^{2}$ and
$\partial_{i}^{2}\pi(\partial_{j}\pi)^{2}$. In terms of $\zeta$, these terms
become

$\displaystyle\dot{\pi}(\partial\pi)^{2}\rightarrow-\frac{1}{aH^{3}}\zeta^{\prime}(\partial\zeta)^{2}$
(3.17)
$\displaystyle\partial^{2}\pi(\partial\pi)^{2}\rightarrow-\frac{1}{H^{3}}\partial^{2}\zeta(\partial\zeta)^{2}.$
(3.18)

Next, similar to calculating power spectrum, $\zeta_{k}$ is written as a
quantum field operator in terms of the creation and annihilation operators
$a_{\textbf{k}}^{\dagger}$, $a_{-\textbf{k}}$ satisfying the canonical
commutation relations
$[a_{\textbf{k}},a_{\textbf{k}}^{\dagger}]=(2\pi)^{3}\delta^{3}(\textbf{k}-\textbf{k}^{\prime})$
and the mode functions $f_{k}(\tau)$:

$\zeta(\textbf{k},\tau)=a_{-\textbf{k}}f_{k}(\tau)+a_{\textbf{k}}^{\dagger}f^{*}_{k}(\tau).$
(3.19)

In practice, we carried out the calculations numerically in terms of the modes
functions
$\tilde{u}(c_{{}_{S}}k\tau)=-\frac{\sqrt{c_{{}_{S}}kA_{1}}a}{H}f_{k}(\tau)$
which we had computed previously from quadratic action and performed the
integration in terms of the dimensionless variable $x_{i}=c_{{}_{S}}k_{i}\tau$
111111See appendix (B) for more details.. Under these changes of variables,
contributions from the dominant cubic interaction terms to the three point
function were found to be as follows,

1\. contributions from the $\dot{\pi}(\partial\pi)^{2}$ term to
$\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>$
is:

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>_{\dot{\pi}(\partial\pi)^{2}}=$
$\displaystyle(2\pi)^{3}\delta^{3}(\textbf{k}_{3}+\textbf{k}_{2}+\textbf{k}_{1})\frac{H^{12}}{c_{{}_{S}}^{8}A_{1}^{3}}\left[-\delta_{3}+9\frac{\bar{M}_{4}}{H}-2\left(\frac{M_{2}}{H}\right)^{4}\right]$
$\displaystyle\left[i\frac{(1-\lambda^{2}-\theta^{2})x_{f}^{3}}{k_{1}^{6}\lambda\theta}\left\\{\int_{-\infty}^{x_{f}}dx_{1}^{\prime}\;(x_{1}^{\prime})^{2}\right.\right.\tilde{u}(x_{f}){\tilde{u}^{*}}_{,x}(x_{1}^{\prime})$
$\displaystyle\tilde{u}(\lambda\;x_{f}){\tilde{u}^{*}}(\lambda\;x_{1}^{\prime})\tilde{u}(\theta\;x_{f}){\tilde{u}^{*}}(\theta\;x_{1}^{\prime})$
$\displaystyle+\int_{-\infty}^{x_{f}}dx_{1}^{\prime}\;x_{1}^{\prime}\;\tilde{u}(x_{f}){\tilde{u}^{*}}(x_{1}^{\prime})\tilde{u}(\lambda\;x_{f}){\tilde{u}^{*}}(\lambda\;x_{1}^{\prime})$
$\displaystyle\tilde{u}(\theta\;x_{f}){\tilde{u}^{*}}(\theta\;x_{1}^{\prime})+{\rm
c.c.}\bigg{\\}}\bigg{]}+{\rm cycl.}$ (3.20)

where we are defining parameters $\lambda\equiv\frac{k_{2}}{k_{1}}$ and
$\theta\equiv\frac{k_{3}}{k_{1}}$ that determine the shape of the triangles
formed by $\textbf{k}_{1},~{}\textbf{k}_{2}$ and $\textbf{k}_{3}$.

2\. contributions from the $\partial^{2}\pi(\partial\pi)^{2}$ term to
$\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>$
is:

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>_{\partial^{2}\pi(\partial\pi)^{2}}=$
$\displaystyle(2\pi)^{3}\delta^{3}(\textbf{k}_{3}+\textbf{k}_{2}+\textbf{k}_{1})$
$\displaystyle\frac{H^{12}}{c_{{}_{S}}^{10}A_{1}^{3}}\left[\frac{1}{3}\delta_{3}-\frac{1}{2}\left(\left(\frac{\bar{M}_{1}}{H}\right)^{3}+\left(\frac{\bar{M}_{2}}{H}\right)^{2}+2\left(\frac{\bar{M}_{3}}{H}\right)^{2}-\frac{\bar{M}_{4}}{H}\right)\right]$
$\displaystyle\left[i\frac{(1-\lambda^{2}-\theta^{2})x_{f}^{3}}{k_{1}^{6}\lambda\theta}\int_{-\infty}^{x_{f}}dx_{1}^{\prime}(x_{1}^{\prime})^{3}\;\tilde{u}(x_{f}){\tilde{u}^{*}}(x_{1}^{\prime})\tilde{u}(\lambda
x_{f})\right.$ $\displaystyle{\tilde{u}^{*}}(\lambda
x_{1}^{\prime})\tilde{u}(\theta x_{f}){\tilde{u}^{*}}(\theta
x_{1}^{\prime})+{\rm c.c.}\bigg{]}+{\rm cycl.}$ (3.21)

The above derivations allows to numerically evaluate bispectrum for different
shapes of triangles at different scales and for given values of parameters. To
simplify connecting these expressions to our numerical calculations, we break
down their contribution to bispectrum from each interaction, $I$, as

$(4\pi^{2}\Delta^{B.D.}_{\zeta})^{2}~{}\tilde{q}_{I}~{}\left[\frac{1}{k_{1}^{6}}G_{I}(x_{f},\lambda,\theta)+{\rm
sym}\right]\longrightarrow B_{I}(k_{1},k_{2},k_{3})$ (3.22)

where
$(4\pi^{2}\Delta^{B.D.}_{\zeta})^{2}=\frac{H^{8}}{A_{1}^{2}c_{{}_{S}}^{2}}$
from (2.19) and

$\displaystyle\tilde{q}_{\partial^{2}\pi(\partial\pi)^{2}}$
$\displaystyle\equiv\frac{H^{4}}{c_{{}_{S}}^{6}A_{1}}\left[-\delta_{3}+9\frac{\bar{M}_{4}}{H}-2\left(\frac{M_{2}}{H}\right)^{4}\right]\,,$
(3.23) $\displaystyle\tilde{q}_{\dot{\pi}(\partial\pi)^{2}}$
$\displaystyle\equiv\frac{H^{4}}{c_{{}_{S}}^{6}A_{1}}\left[\frac{1}{3}\delta_{3}-\frac{1}{2}\left(\left(\frac{\bar{M}_{1}}{H}\right)^{3}+\left(\frac{\bar{M}_{2}}{H}\right)^{2}+2\left(\frac{\bar{M}_{3}}{H}\right)^{2}-\frac{\bar{M}_{4}}{H}\right)\right]\,.$
(3.24)

We can now numerically calculate contributions to bispectrum for both
operators at various triangle shapes. Figure (2) displays one of our results
for $B_{I}(k_{1},k_{2},k_{3})k_{1}^{2}k_{2}^{2}k_{3}^{2}$, evaluated at
$k_{1}=1$ and normalized over $B_{I}(1,1,1)$, for various triangle shapes.

(a) The shape of bispectrum for operator $\dot{\pi}(\partial\pi)^{2}$

(b) The shape of bispectrum for operator $\dot{\pi}(\partial\pi)^{2}$

Figure 2: Calculation of the quantity
$B_{I}(k_{1},k_{2},k_{3})k_{1}^{2}k_{2}^{2}k_{3}^{2}$, evaluated at $k_{1}=1$
and normalized over $B_{I}(1,1,1)$ for various triangle shapes. This result
corresponds to value of $\alpha_{0}=-0.18$ and $\beta_{0}=\alpha_{0}^{2}/4$.

As we see the bispectrum seems to be enhanced for triangle configurations
other than just folded [23, 24] or squeezed type [25, 26]. This generically
agrees to what one expects from non-Bunch-Davies vacuum initial conditions
although each manifestation leads to distinct bispectrum enhancement for
different shapes and makes it possible to distinguish them from each other
121212The calculation of bispectrum due to these operators for EFToI models
with $\omega^{2}\propto k^{4}$ have also been shown to have non-analytical
forms and result in very different shapes from the folded and local kind[27,
17].. Now to make the connection with observations, we need to turn
observational constraints on bispectrum, into bounds on parameters of our
models by substituting for different mass and coupling parameters in
$\tilde{q}_{i}$ and other variables determining $B_{I}$. We will further
explore these bounds and discuss our numerical calculations in the next
section.

##  4 Numerical Results and Constraints on the Models

We know provide more details on our numerical calculations and how we compared
them with observtions in order to interpret the constraints on our model.

To quantify the strength of non-gaussianity due to terms discussed in last
section and connect them to observations, it is conventional to express the
result in terms of the quantities known as $f_{{}_{\rm NL}}$ in the
literature.

Historically, $f_{{}_{\rm NL}}$ was first defined to characterise the local
non-Gaussianity in CMB data [28]. In that case, one assumes the non-Gaussian
corrections for curvature perturbation can be expressed locally as

$\zeta=\zeta_{{}_{\rm L}}+\frac{3}{5}f_{{}_{\rm NL}}\zeta_{{}_{\rm L}},$ (4.1)

where $\zeta_{\rm L}$ is the gaussian scalar perturbation. Now under this
assumption as well as near scale invariance (i.e. $n_{{}_{\rm S}}-1\ll 1$),
one can obtain the following relationship between the three point function in
k space and the dimensionless power spectrum of $\zeta$ 131313See appendix (C)
for details.,

$\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>=(2\pi)^{3}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})[(2\pi)^{4}\Delta_{{}_{\zeta}}^{2}]~{}f_{{}_{\rm
NL}}^{{}^{\rm
local}}\times\left\\{\frac{3}{10}~{}\bigg{(}\frac{k_{1}^{3}+k_{2}^{3}+k_{3}^{3}}{k_{1}^{3}k_{2}^{3}k_{3}^{3}}\bigg{)}\right\\}.$
(4.2)

Therefore, assuming a local form bispectrum and comparing (3.1) and (4.2), one
can construct a dimensionless quantity directly from bispectrum itself
consistent with the definition in (4.1),

$f_{{}_{\rm NL}}^{{}^{\rm local}}=\left(\frac{10}{3}\frac{\prod
k_{i}^{3}}{\sum
k_{i}^{3}}\right)\times\frac{1}{(2\pi)^{4}\Delta_{{}_{\zeta}}^{2}}B(\textbf{k}_{1},\textbf{k}_{2},\textbf{k}_{3}).$
(4.3)

In practice, even if the theory produces exact local form bispectrum, the data
is going to be contaminated with noise and other secondary effects so direct
application of above formula to data would lead to extracting different values
of $f_{{}_{\rm NL}}^{{}^{\rm local}}$ for different triangles in momentum
spaces,

$f_{{}_{\rm NL}}^{{}^{\rm
local}}(\textbf{k}_{1},\textbf{k}_{2},\textbf{k}_{3})=\left(\frac{10}{3}\frac{\prod
k_{i}^{3}}{\sum
k_{i}^{3}}\right)\times\frac{1}{(2\pi)^{4}\Delta_{{}_{\zeta}}^{2}}B(\textbf{k}_{1},\textbf{k}_{2},\textbf{k}_{3}).$
(4.4)

Therefore, finding the observed value of $f_{{}_{\rm NL}}^{{}^{\rm local}}$ is
a more subtle task of optimisation of the fitting between the template in
(4.2) and all the data points and it involves statistical schemes beyond the
scope of this work. It can easily be seen that the particular momentum
dependence of (4.2) diverges in squeezes limits where one of the modes
$k_{i}\to 0$, so this template is suited for fitting models where
contributions from those corners are enhanced.

However, in general the three point function for different interactions and in
different models does not necessarily lead to k-dependence presented in (4.2)
and it can have crests and troughs at different locations. Therefore, over the
years there have been different templates developed to best extract possible
non-gaussianities in data141414Readers are referred to [29] for details on
different templates fitting for different models. in order to test different
types of models. The next famous example are non-canonical single scalar field
inflationary models [30, 31] that usually peak at the limits of equilateral
shape triangles [32] and can be approximated by

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>^{{}^{\rm
equil.}}$ $\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})[(2\pi)^{4}\Delta_{{}_{\zeta}}^{2}]~{}f_{{}_{\rm
NL}}^{{}^{\rm equil.}}\times\bigg{(}-\frac{9}{10}\bigg{)}$ (4.5)
$\displaystyle\left\\{\left[\frac{1}{k_{1}^{3}k_{2}^{3}}+\text{cycl.}\right]+\frac{2}{k_{1}^{2}k_{2}^{2}k_{3}^{2}}-\left[\frac{1}{k_{1}k_{2}^{2}k_{3}^{3}}+(5~{}\text{perm.})\right]\right\\}.$

For agnostic bispectra interpretation one can perform a more sophisticated
mathematical approach of projecting the observed bispectra on different known
and independent templates by defining an inner product between bispectrum and
these templates [27, 33] and then use joint estimators to simultaneously fit
for several types of non-gaussinities. Of course, given that these functions
do not form a complete basis for all the possible models they may not
necessary capture all aspects of non-gaussianity. Non-Bunch-Davies bispectra
have been found to get enhanced for both the flattened and squeezed shapes
[32, 23, 25, 24, 26]. Planck collaboration [29] has used variety of these
templates as model estimators in order to test non-Bunch-Davies (NBD) vacuum
initial conditions. For example, one of the simple flatten triangles shaped
enhancements that can rise due to NBD models [32, 23] is

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>^{{}^{\rm
NBD}}$ $\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})[(2\pi)^{4}\Delta_{{}_{\zeta}}^{2}]~{}f_{{}_{\rm
NL}}^{{}^{\rm flat}}$ (4.6)
$\displaystyle\times\left\\{\frac{3}{10}~{}\bigg{(}\frac{1}{\prod
k_{i}}\bigg{)}\left[\frac{1}{(k_{2}+k_{3}-k_{1})^{3}}+\text{cycl.}\right]\right\\}+\dots\,,$

which diverges as $k_{1}\to k_{2}+k_{3}$. Another template that has been used
for testing EFToI models is the orthogonal template defined as

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>^{{}^{\rm
ortho.}}$ $\displaystyle=$
$\displaystyle(2\pi)^{3}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})[(2\pi)^{4}\Delta_{{}_{\zeta}}^{2}]~{}f_{{}_{\rm
NL}}^{{}^{\rm ortho}}\times\bigg{(}-\frac{9}{10}\bigg{)}$ (4.7)
$\displaystyle\left\\{\left[\frac{3}{k_{1}^{3}k_{2}^{3}}+\text{cycl.}\right]+\frac{8}{k_{1}^{2}k_{2}^{2}k_{3}^{2}}-\left[\frac{3}{k_{1}k_{2}^{2}k_{3}^{3}}+(5\,\text{perm.})\right]\right\\}.$

For our purposes here, since we want to make only order of magnitude
estimations and since $f_{{}_{\rm NL}}$ constraints from observations are only
reported for these known template, we will make a crude estimation by making a
linear fit to three templates of local, equilateral and orthogonal templates
and then compare our result to observational constraints reported for these
shapes. However, as figure (2) indicates the shape of bispectrum can not be
described through just linear combination of these templates and there is some
residual non-gausianities left. A more robust way of constraining our model
which is beyond the scope of this work, would be direct comparison of our
numerical template with the data. However, here after substituting for
$B_{I}(\textbf{k}_{1},\textbf{k}_{2},\textbf{k}_{3})$ from (3.22) and that
$\Delta_{{}_{\zeta}}=\gamma_{{}_{S}}\Delta^{B.D}_{\zeta}$, we approximate
different contributions to $f_{{}_{NL}}^{{}^{\rm
local}},~{}f_{{}_{NL}}^{{}^{\rm equil.}}$ and $f_{{}_{NL}}^{{}^{\rm ortho.}}$
in terms of $G_{I}$ and $\tilde{q}_{I}$ as following linear combination

$\frac{\tilde{q}_{I}~{}}{\gamma_{{}_{S}}^{2}}\left[\frac{1}{k_{1}^{6}}G_{I}(x_{f},\lambda,\theta)+\text{cycl}\right]=f_{{}_{\rm
NL}}^{{}^{(I)\rm local}}~{}T^{{}^{\rm local}}+f_{{}_{\rm NL}}^{{}^{(I)\rm
equil.}}~{}T^{{}^{\rm equil.}}+f_{{}_{NL}}^{{}^{(I)\rm ortho.}}~{}T^{{}^{\rm
ortho.}}.$ (4.8)

The template functions $T^{(\dots)}(k_{1},k_{2},k_{3})$ are the templates
mentioned in (4.2), (4.5) and (4.7):

$\displaystyle T^{{}^{\rm local}}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle\frac{3}{10}~{}\left\\{\frac{k_{1}^{3}+k_{2}^{3}+k_{3}^{3}}{k_{1}^{3}k_{2}^{3}k_{3}^{3}}\right\\}$
$\displaystyle T^{{}^{\rm equil.}}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle-\frac{9}{10}\left\\{\left[\frac{1}{k_{1}^{3}k_{2}^{3}}+\text{cycl.}\right]+\frac{2}{k_{1}^{2}k_{2}^{2}k_{3}^{2}}-\left[\frac{1}{k_{1}k_{2}^{2}k_{3}^{3}}+(5\,\text{perm.})\right]\right\\}$
(4.9) $\displaystyle T^{{}^{\rm ortho.}}(k_{1},k_{2},k_{3})$ $\displaystyle=$
$\displaystyle-\frac{9}{10}\left\\{\left[\frac{3}{k_{1}^{3}k_{2}^{3}}+\text{cycl.}\right]+\frac{8}{k_{1}^{2}k_{2}^{2}k_{3}^{2}}-\left[\frac{3}{k_{1}k_{2}^{2}k_{3}^{3}}+(5\,\text{perm.})\right]\right\\}.$
(4.10)

Given that $k_{2}/k_{1}=\lambda$ and $k_{3}/k_{1}=\theta$ the problem reduces
to finding the best fit for the coefficients $\frac{f_{{}_{\rm
NL}}^{(I)\dots}\gamma_{{}_{S}}^{2}}{\tilde{q}_{I}~{}}$ in this relation

$\displaystyle\left[\frac{1}{k_{1}^{6}}G_{I}(x_{f},\lambda,\theta)+\text{cycl.}\right]_{k_{1}=1}\approx$
$\displaystyle\frac{f_{{}_{\rm NL}}^{{}^{(I)\rm
local}}\gamma_{{}_{S}}^{2}}{\tilde{q}_{I}~{}}T^{{}^{\rm
local}}(1,\lambda,\theta)+\frac{f_{{}_{\rm NL}}^{{}^{(I)\rm
equil.}}\gamma_{{}_{S}}^{2}}{\tilde{q}_{I}~{}}T^{{}^{\rm
equil}}(1,\lambda,\theta)$ (4.11) $\displaystyle+\frac{f_{{}_{\rm
NL}}^{{}^{(I)\rm ortho.}}\gamma_{{}_{S}}^{2}}{\tilde{q}_{I}~{}}T^{{}^{\rm
ortho.}}(1,\lambda,\theta).$

Next, we provide a brief summary of our numerical scheme for readers. Note
that as it was described in section (2), we had computed power spectrum which
itself involved numerical computation of the canonical mode function
$\tilde{u}(y)$ where $y=c_{{}_{S}}k\tau$. This computation was carried out for
different values of $\alpha_{0}$ and $\beta_{0}$ that themselves are functions
of the parameters in the EEFToI action. We had found empirically that in the
regime that the dispersion relation does not become tachyonic inside the
horizon, the amplification to the power spectrum $\gamma_{{}_{S}}$ was the
largest (see Figure (1)) in the following region

$-0.5\leq\alpha_{0}\leq-0.1\text{\hskip 14.22636ptand\hskip
14.22636pt}0\leq\frac{\alpha_{0}^{2}}{4}\leq\beta_{0}\leq\frac{\alpha_{0}^{2}}{3}.$
(4.12)

In fact, within that region the enhancement was found to be the greatest
($\gamma_{{}_{S}}\geq 14.6$) when setting $\alpha_{0}=-0.18$ and
$\beta_{0}=\alpha_{0}^{2}/4$. From physics point of view, a more interesting
applications of the EEFToI and EFToI is when they result in larger values of
$\gamma_{{}_{S}}$ or significant deviations from the Bunch Davis initial
conditions. Therefore, it is also more interesting to study the impact of
observational constraints on bispectrum in this particular region. To do that,
numerically computed values of $\tilde{u}(y)$ for different values of
$\alpha_{0}$ and $\beta_{0}$ were substituted in the equations (3) and (3) and
then numerically integrated to obtain the contributions to bispectrum.

In order to set the initial conditions at $x\to-\infty$, we chose the value of
$x_{initial}$ to be where $\beta_{0}x^{4}$ is far greater than
$\alpha_{0}x^{2}$ ($x_{initial}\sim-50|\alpha_{0}/\beta_{0}|^{1/2}$). This
guarantees that the initial condition is set in the regime where
$\omega^{2}\propto k^{6}$ dominates the dispersion relation. For the upper
limit of the integration, we choose $x_{f}=-0.00001$. This corresponds to
setting different values of $\tau_{f}$ for different $k_{i}$, but ensures all
the modes that were considered had crossed the freezing radius at that point.
We found the result of the computation was not very sensitive to changing
$x_{initial}$ further into the past or $x_{f}$ into the present. The triangle
shapes that were explored corresponded to values of $0.2\leq\lambda\leq 1$ and
$0.2\leq\theta\leq 1$. The lower bounds for $\lambda$ and $\theta$ should in
principle either be set by ratio of the observed modes in CMB ($\sim 10^{3}$)
or theoretically by validity of the effective field theory regime. Note, the
extreme squeezed limit $\theta\to 0$ or $\lambda\to 0$ can push the initial
conditions for two of the modes into ultraviolet scales or the other into
infrared. Therefore, one needs to ensure all of the mode functions such as
$\tilde{u}(\lambda x)$ also initiate in the same regime. In practice, setting
$0.2$ as the lower bound for both $\lambda$ and $\theta$ was due to the
limitation of our numerical method. That is, we chose the values of $\lambda$
and $\theta$ such that their ratios varied within $0.2$ and $1$ consistent
with our criteria but not as far as the observational and theoretical
threshold would allow it in the squeezed limit. However, our numerical
solution did not indicate any divergence in the asymptotic behaviour of
bispectrum.

Figure (3) represents the shape of the best fit linear model for
$\dot{\pi}(\partial\pi)^{2}$ to $T^{{}^{\rm local}}$, $T^{{}^{equil.}}$ and
$T^{{}^{ortho.}}$ with coefficients $\frac{f_{{}_{\rm
NL}}^{{}^{(I)\dots}}\gamma_{{}_{S}}^{2}}{\tilde{q}_{I}}$ given in table (1)
and the residual bispectrum after subtracting the best fit. The mode functions
in this case correspond to the values of $\alpha_{0}=-0.18$ and
$\beta_{0}=\alpha_{0}^{2}/4$. Similarly, figure (4) displays the best linear
fit model to the same templates for $\partial^{2}\pi(\partial\pi)^{2}$
operator and the residual bispectrum, corresponding to the same values of
$\alpha_{0}$ and $\beta_{0}$ where obtained values of $\frac{f_{{}_{\rm
NL}}^{{}^{(I)\dots}}\gamma_{{}_{S}}^{2}}{\tilde{q}_{I}~{}}$ are given in table
(2).

Figure 3: Left plot represents the shape of the best fit linear model for $\dot{\pi}(\partial\pi)^{2}$ to $T^{{}^{\rm local}}$, $T^{{}^{\rm equil.}}$ and $T^{{}^{\rm ortho.}}$ with coefficients given in table (1) and right plot is the residual bispectrum after subtracting the best fit. The mode functions correspond to values of $\alpha_{0}=-0.18$ and $\beta_{0}=\alpha_{0}^{2}/4$. $\dot{\pi}(\partial\pi)^{2}$ |  Estimate |  Standard Error
---|---|---
$f_{{}_{\rm NL}}^{{}^{\rm(I)local}}\gamma_{{}_{S}}^{2}/\tilde{q}_{I}$ |  -65 |  2.9
$f_{{}_{\rm NL}}^{{}^{\rm(I)equil.}}\gamma_{{}_{S}}^{2}/\tilde{q}_{I}$ |  881 |  16
$f_{{}_{\rm NL}}^{{}^{\rm(I)ortho}}\gamma_{{}_{S}}^{2}/\tilde{q}_{I}$ |  -339 |  7
Table 1: The coefficients of the best fit linear model for
$\dot{\pi}(\partial\pi)^{2}$ to $T^{{}^{\rm local}}$, $T^{{}^{\rm equil.}}$
and $T^{{}^{ortho.}}$ corresponding to values of $\alpha_{0}=-0.18$ and
$\beta_{0}=\alpha_{0}^{2}/4$.

Figure 4: Left plot represents the shape of best fit linear model for $\partial^{2}\pi(\partial\pi)^{2}$ to $T^{{}^{\rm local}}$, $T^{{}^{\rm equil.}}$ and $T^{{}^{\rm ortho.}}$ with coefficients given in table (2) and right plot is the residual bispectrum after subtracting the best fit. The mode functions correspond to values of $\alpha_{0}=-0.18$ and $\beta_{0}=\alpha_{0}^{2}/4$. $\partial^{2}\pi(\partial\pi)^{2}$ |  Estimate |  Standard Error
---|---|---
$f_{{}_{\rm NL}}^{{}^{(I)\rm local}}\gamma_{{}_{S}}^{2}/\tilde{q}_{I}$ |  -282 |  29
$f_{{}_{\rm NL}}^{{}^{(I)\rm equil.}}\gamma_{{}_{S}}^{2}/\tilde{q}_{I}$ |  3610 |  344
$f_{{}_{\rm NL}}^{{}^{(I){\rm ortho.}}}\gamma_{{}_{S}}^{2}/\tilde{q}_{I}$ |  -1033 |  99
Table 2: The coefficients of the best fit linear model for $\partial^{2}\pi(\partial\pi)^{2}$ to $T^{{}^{\rm local}}$, $T^{{}^{\rm equil.}}$ and $T^{{}^{\rm ortho.}}$ corresponding to values of $\alpha_{0}=-0.18$ and $\beta_{0}=\alpha_{0}^{2}/4$. Constraints |  Note
---|---
$\alpha_{0}<0$ |  Regions where the power spectrum is enhanced. Note that for $\alpha_{0}>0$ the power spectrum is suppressed (not interesting)
$\beta_{0}>0$ |  To avoid tachyonic instabilities (note $\delta_{3}>0$)
$\Delta_{{}_{\zeta}}\sim 10^{-9}$ |  The constraint from power spectrum of the CMB observation ($1\lesssim\gamma_{{}_{S}}\lesssim 14$).
$H>\frac{\sqrt{A_{1}}}{M_{{}_{\rm Pl}}}$ |  Decoupling limit
$0\leq c_{{}_{T}}^{2}\leq 1$ |  Speed of sound for the tensor modes (Note: $\bar{M}_{3}^{2}<0$)
$0\leq c_{{}_{S}}^{2}=\frac{F_{1}}{A_{1}}\leq 1$ |  To avoid superluminality in the regime $\omega^{2}\sim k^{2}$
$A_{1}\geq 0$ |  To avoid ghosts at the leading order
$\dot{H}<0$ |  Inflation satisfies the null energy condition
$\tilde{q}_{i}\leq 10^{-1}$ |  Constraints from $f_{{}_{\rm NL}}$ measurements
Table 3: List of the constraints applied to the EEFToI models

We now try to explore the allowed regions of the parameters for EEFToI models
or rather impose constraints on them, based on the observed constraints on the
size of the non-Gaussianity. The most recent CMB observation by Planck
collaboration [29] the following bounds on the size of the non-Gaussianities
for the local, equilateral and orthogonal templates: $f^{{}^{\rm
local}}_{{}_{\rm NL}}\simeq-0.9\pm 5.1$, $f^{{}^{\rm equil.}}_{{}_{\rm
NL}}\simeq-26\pm 47$ and $f^{{}^{\rm ortho.}}_{{}_{\rm NL}}\simeq-38\pm 24$.
Now to make the connection to our model, we make a comparison between these
observational constraints for different shapes and the best linear fit values
of $f_{{}_{\rm NL}}^{{}^{\rm(I)}}\gamma_{{}_{S}}^{2}/\tilde{q}_{I}$ to infer
bounds on value of $\tilde{q}_{I}$. We then translate this into constraints on
the parameters of our models in the action by substituting for different mass
and coupling parameters in $\gamma_{{}_{S}}^{2}$ and $\tilde{q}_{i}$. We
require for each shape estimated $f_{{}_{\rm NL}}^{{}^{(I)}}$ be less than
reported value, which leads to

$\left|\frac{f_{{}_{\rm
NL}}^{{}^{(}I)}\gamma_{{}_{S}}^{2}}{\tilde{q}_{I}}\right|^{{}^{\rm
estimate}}<\frac{|f_{{}_{\rm NL}}^{{}^{\rm
reported}}|\gamma_{{}_{S}}^{2}}{|\tilde{q}_{I}|}\implies|\tilde{q}_{I}|<\frac{|f_{{}_{\rm
NL}}^{{}^{\rm reported}}|}{\left|\frac{f_{{}_{\rm
NL}}^{(I)}\gamma_{{}_{S}}^{2}}{\tilde{q}_{I}}\right|^{{}^{\rm
estimate}}}\gamma_{{}_{S}}^{2}\lesssim 10^{-2}-10^{-3}\gamma_{{}_{S}}^{2}$
(4.13)

Figure 5: Parameters search of EEFToI model under the constraints listed in
table (3). We searched for allowed values of
$M_{2},\bar{M}_{1},\bar{M}_{2},\bar{M}_{3}^{2},\bar{M}_{4},\delta_{3},H,\dot{H}$.
Furthermore, we set an additional cutoff at $M_{{}_{\rm Pl}}$, except for
$\delta_{3}$ which is dimensionless. Although the points in the plots are
scattered, this is due to the method carried out for searching through
parameters. Instead of searching the entire grid, which would take enormous
computation time, we took a random search method. Note that the histogram does
not represent the probability of each parameter. It just means that the
algorithm searched in those respective regions more frequently.

Therefore, we can crudely set a bound on
$\tilde{q}_{I}\lesssim~{}10^{-2}-10^{-3}\ \gamma_{{}_{S}}^{2}$. In this case
substituting larger than unity values such as $\gamma_{{}_{S}}\sim 14$ results
into a not very tight bound of $|\tilde{q}_{I}|\lesssim 10^{-1}$. Now to
address what this constraint mean in terms of the original couplings and
coefficients in EEFToI action, and whether there are corners of the parameter
space that lead to scenarios we have explored, we run a search algorithm in
the parameter space under the list of all the criteria we have set for the
model. Table (3) shows the list of all these conditions and constraints. The
primary parameters that were included in EEFToI were
$M_{2},\bar{M}_{1},\bar{M}_{2},\bar{M}_{3},\bar{M}_{4},\delta_{3},H,M_{{}_{\rm
Pl}}^{2}\dot{H}$.

Figure (5) displays the result of a search algorithm we performed to find
valid values of parameters of EEFToI imposing conditions (3). We also have
imposed an additional cutoff on all of the mass dimension parameters to be
below the Planck scale. That is to avoid break down of semi-classical gravity
approximation. Figure (5), visualizes the result of the search algorithm
indicating the values of parameters that were found to satisfy the above
constraints. Therefore, we can infer that there are allowed regions in the
mass and coupling parameter space, arising from the sixth order polynomial
dispersion relations [13] and leading to a sensible effective field theory
without violating the current observations constraints on non-gaussianity.

##  5 Conclusion

In earlier works, we had examined the power spectrum for some interesting
corners of the parameter space in the EEFToI models with the dispersion
relation $\omega^{2}\propto k^{6}$ when the wavelength of the mode is very
small [13, 14]. The dispersion relation then evolves into quartic and
quadratic before the modes exit the horizon. In this work, we studied the
bispectrum for these models in the regime where the amplification of the power
spectrum due to the modified dispersion relation is more significant. One
implications of taking into account such possibilities in EEFToI is that they
can lead to lower values of $H$ during inflation (2.20) and still be
consistent with the observed value of amplitude for scalar perturbations. This
widens the possibilities of viable low energy inflationary models.

Our power counting arguments in section (3) indicated that the most pronounced
and possibly dangerous contribution of the cubic interaction terms are coming
from $\dot{\pi}(\partial_{i}\pi)^{2}$ and
$\partial_{i}^{2}\pi(\partial_{j}\pi)^{2}$, we computed their impact. We found
the triangle shapes corresponding to peaks and cresses in the bispectrum
(Figure 2) are quite distinct from other inflationary models. In particular,
we made a linear fit of the shapes of the bispectra for these two operators,
to the three known and widely used templates of local, equilateral and
orthogonal shapes. Interestingly, the major contributions to the bispectra for
these two cubic interaction operators was projected on the equilateral
configuration, whereas in the other EFT of inflation with non-Bunch-Davies
vacuum, it is often the orthogonal template which is dominant. This
calculation also allowed us to inspect the strong coupling constraints and the
restrictions due to the size of non-Gaussianity from the CMB observations. We
translated the reported observational constraints from Planck Collaboration
[29] in values of $f_{{}_{\rm NL}}$ in terms of the parameters of the model.
This lead to an order of magnitude estimation for the bounds on variables
$\tilde{q}_{I}$ which include the masses and couplings in the model defined in
(3.23) as $\tilde{q}_{I}\lesssim 10^{-1}$. Finally, we used a search algorithm
and found that there are allowed regions of parameter space where EEFToI
description with $\omega^{2}\propto k^{6}$ is sensible and interesting. For
future, it is intriguing to perform a direct comparison of our numerical
template for bispectrum with the data itself and investigate its signatures.
Another aspect that we leave for future work is analysing the trispectrum and
specifically the impact and signatures of the term $(\partial_{i}\pi)^{4}$ in
the Extended EFT of inflation. Furthermore, we like to add that the terms we
have included in this analysis are not the exhaustive list of all the possible
terms consistent with spatial diffeomorphisms in unitary gauge upto mass
dimension 4. Our analysis has been more focused on terms leading to
interesting $\omega^{2}\propto k^{6}$ regimes and enhanced signatures of
purely scalar bispectrum for such regimes. However, there are other possible
regimes and also contractions of $g^{00}+1,\delta K_{\mu\nu},\delta
R_{\mu\nu\rho\sigma},\nabla_{\mu}$ that through their quadratic or interaction
contributing terms in the action, could potentially lead to interesting and
distinctive observational signatures.

##  6 Acknowledgment

A.A. has received funding/support from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska -Curie grant
agreement No 860881-HIDDeN. G.G. and H.J.K. research is supported by the
Discovery Grant from Natural Science and Engineering Research Council of
Canada (NSERC) and G.G. is supported partly by Perimeter Institute (PI) as
well. Research at PI is supported by the Government of Canada through the
Department of Innovation, Science and Economic Development Canada and by the
Province of Ontario through the Ministry of Research, Innovation and Science.

##  Appendix A The Interaction Lagrangian for the EEFToI

In the section we summarize all the cubic terms in the action $\int
d^{4}x\sqrt{g}\sum_{d=3}^{4}\mathcal{L}_{2,d}$ generated through the new
EEFToI terms (2.8) after setting $\delta_{1}=\delta_{2}=0$.

The contributions from the operator proportional to $\bar{M}_{4}$ to the cubic
interactions in the action, taking into account $\sqrt{g}$, goes as:

$a^{3}\nabla^{\mu}\delta g^{00}\nabla^{\nu}\delta
K_{\mu\nu}\longrightarrow-\frac{3H(\partial_{i}\pi)^{2}\partial_{j}^{2}\dot{\pi}}{a}+\frac{H^{2}(\partial_{i}\pi{})^{2}\partial_{j}^{2}\pi}{a}+18aH^{3}\dot{\pi}(\partial_{i}\pi{})^{2}-14aH^{2}\dot{\pi}\partial_{i}\pi\partial_{i}\dot{\pi}-8aH\dot{\pi}\ddot{\pi}\partial_{i}^{2}\pi+18a^{3}H^{2}\dot{\pi}^{2}\ddot{\pi}+\frac{\partial_{j}\pi\partial_{i}\partial_{j}\pi\partial_{i}\partial_{k}^{2}\pi}{a^{3}}+\frac{\partial_{i}\pi\partial_{i}\partial_{j}\pi\partial_{k}^{2}\partial_{j}\pi}{a^{3}}+\frac{2\ddot{\pi}\partial_{i}\partial_{j}\pi{}^{2}}{a}+\frac{\partial_{j}\pi\partial_{i}\dot{\pi}\partial_{i}\partial_{j}\dot{\pi}}{a}+\frac{5\partial_{i}\dot{\pi}\partial_{i}\partial_{j}\pi\partial_{j}\dot{\pi}}{a}+\frac{\partial_{i}\pi\partial_{i}\dot{\pi}\partial_{j}\partial_{j}\dot{\pi}}{a}$
(A.1)

The contributions from the operator proportional to the $\delta_{3}$ to the
cubic interaction terms in the action goes as:

${a^{3}}\nabla^{\mu}\delta
K_{\nu\mu}\nabla^{\nu}\delta{K^{\sigma}}_{\sigma}\longrightarrow
3a\dot{\pi}\left(\partial_{i}\pi\right){}^{2}H^{4}+18a^{3}\dot{\pi}^{2}\ddot{\pi}H^{3}-21a\dot{\pi}\partial_{i}\pi\partial_{i}\dot{\pi}H^{3}+6a\dot{\pi}^{2}\partial_{i}^{2}\pi
H^{3}-\frac{\left(\partial_{i}\pi\right){}^{2}\partial_{j}^{2}\pi
H^{3}}{a}+3a\dot{\pi}\partial_{i}\pi\partial_{i}\ddot{\pi}H^{2}+\frac{4\dot{\pi}\partial_{i}\pi\partial_{i}\pi\partial_{j}^{2}\pi
H^{2}}{a}-6a\dot{\pi}\ddot{\pi}\partial_{i}^{2}\pi
H^{2}-3a\dot{\pi}^{2}\partial_{i}^{2}\dot{\pi}H^{2}+\frac{2\partial_{i}\dot{\pi}\partial_{i}\partial_{j}\pi\partial_{j}\pi
H^{2}}{a}+\frac{9\partial_{i}^{2}\pi\partial_{j}\pi\partial_{j}\dot{\pi}H^{2}}{a}-\frac{4\dot{\pi}\partial_{i}^{2}\pi\partial_{j}^{2}\pi
H^{2}}{a}+18a^{3}\dot{\pi}^{3}\dot{H}H^{2}-\frac{3}{2}a\dot{\pi}\left(\partial_{i}\pi\right){}^{2}\dot{H}H^{2}+\frac{3\ddot{\pi}\left(\partial_{i}\partial_{j}\pi\right){}^{2}H}{a}+\frac{3\partial_{i}^{2}\pi\left(\partial_{j}\dot{\pi}\right){}^{2}H}{2a}+\frac{5\dot{\pi}\partial_{i}\dot{\pi}\partial_{i}\partial_{j}^{2}\pi
H}{a}+\frac{5\partial_{i}\dot{\pi}\partial_{i}\left(\partial_{j}\dot{\pi}\right)\partial_{j}\pi
H}{2a}+\frac{17\partial_{i}\dot{\pi}\partial_{i}\left(\partial_{j}\pi\right)\partial_{j}\dot{\pi}H}{2a}+\frac{3\partial_{i}^{2}\dot{\pi}\partial_{j}\pi\partial_{j}\dot{\pi}H}{2a}-\frac{4\partial_{i}^{2}\pi\partial_{j}\pi\partial_{j}\partial_{k}^{2}\pi
H}{a^{3}}-\frac{\partial_{i}\pi\partial_{i}\ddot{\pi}\partial_{j}^{2}\pi
H}{a}-\frac{\ddot{\pi}\partial_{i}^{2}\pi\partial_{j}^{2}\pi
H}{a}+\frac{2\dot{\pi}\partial_{i}^{2}\pi\partial_{j}^{2}\dot{\pi}H}{a}+\frac{4\partial_{i}^{2}\pi\left(\partial_{j}\partial_{k}\pi\right)^{2}\pi
H}{a^{3}}+\frac{2\left(\partial_{i}\partial_{j}\pi\right){}^{2}\partial_{k}^{2}\pi
H}{a^{3}}-\frac{2\partial_{i}\pi\partial_{i}\partial_{j}^{2}\pi\partial_{k}^{2}\pi
H}{a^{3}}-27a\dot{\pi}^{2}\partial_{i}^{2}\pi\dot{H}H+\frac{3\left(\partial_{i}\pi\right){}^{2}\partial_{j}^{2}\pi\dot{H}H}{2a}\\\
-\frac{4\dot{\pi}\left(\partial_{i}\partial_{j}^{2}\pi\right)^{2}}{a^{3}}-\frac{7\partial_{i}\partial_{j}\pi\partial_{i}\partial_{k}^{2}\pi\partial_{j}\dot{\pi}}{2a^{3}}-\frac{3\partial_{i}\pi\partial_{i}\partial_{k}\dot{\pi}\partial_{j}^{2}\partial_{k}\pi}{2a^{3}}-\frac{\partial_{i}^{2}\dot{\pi}\partial_{j}^{2}\partial_{k}\pi\partial_{k}\pi}{a^{3}}-\frac{2\partial_{i}^{2}\dot{\pi}\partial_{j}\partial_{k}\pi\partial_{k}\partial_{j}\pi}{a^{3}}-\frac{3\partial_{i}^{2}\partial_{j}\pi\partial_{j}\dot{\pi}\partial_{k}^{2}\pi}{2a^{3}}-\frac{\left(\partial_{i}\partial_{j}\pi\right){}^{2}\partial_{k}^{2}\dot{\pi}}{a^{3}}-\frac{\partial_{i}^{2}\partial_{j}\pi\partial_{j}\pi\partial_{k}^{2}\dot{\pi}}{2a^{3}}-\frac{2\dot{\pi}\left(\partial_{i}^{2}\partial_{j}\pi\right)^{2}}{a^{3}}+\frac{15\dot{\pi}(\partial_{j}\partial_{i}\pi)^{2}\dot{H}}{a}$

The contribution of the operator proportional to the $\delta_{4}$ in
Lagrangian to the cubic interactions is given by:

${a^{3}}\nabla_{\mu}\delta{K^{\mu}}_{\nu}\nabla_{\gamma}\delta
K^{\gamma\nu}\longrightarrow\frac{2H\partial_{j}^{2}\pi\left(\partial_{i}\partial_{k}\pi\right){}^{2}}{a^{3}}-\frac{2H\partial_{i}\pi\partial_{j}^{2}\pi\partial_{i}\partial_{k}^{2}\pi}{a^{3}}-\frac{6\dot{\pi}H^{2}\left(\partial_{i}\partial_{j}\pi\right){}^{2}}{a}+\frac{4H^{2}\partial_{i}\pi\partial_{j}^{2}\pi\partial_{i}\dot{\pi}}{a}-\frac{2\dot{\pi}H^{2}\partial_{i}^{2}\pi\partial_{j}^{2}\pi}{a}+\frac{4\dot{\pi}H\partial_{i}\partial_{j}^{2}\pi\partial_{i}\dot{\pi}}{a}+\frac{H\partial_{j}\pi\partial_{i}\dot{\pi}\partial_{i}\partial_{j}\dot{\pi}}{a}+\frac{5H\partial_{i}\partial_{j}\pi\partial_{i}\dot{\pi}\partial_{j}\dot{\pi}}{a}+\frac{H\partial_{j}^{2}\pi\left(\partial_{i}\dot{\pi}\right){}^{2}}{a}+\frac{H\partial_{i}\pi\partial_{i}\dot{\pi}\partial_{j}^{2}\dot{\pi}}{a}-\frac{H^{3}\left(\partial_{i}\pi\right){}^{2}\partial_{j}^{2}\pi}{a}+3\dot{\pi}aH^{4}\left(\partial_{i}\pi\right){}^{2}-6\dot{\pi}aH^{3}\partial_{i}\pi\partial_{i}\dot{\pi}+12\dot{\pi}^{2}aH^{3}\partial_{i}^{2}\pi-2\dot{\pi}aH^{2}\left(\partial_{i}\dot{\pi}\right){}^{2}-\frac{2\dot{\pi}(\partial_{i}\partial_{j}^{2}\pi)^{2}}{a^{3}}-\frac{\partial_{j}\pi\partial_{i}\partial_{k}^{2}\pi\partial_{i}\partial_{j}\dot{\pi}}{a^{3}}-\frac{5\partial_{i}\partial_{j}\pi\partial_{i}\partial_{k}^{2}\pi\partial_{j}\dot{\pi}}{a^{3}}-\frac{\partial_{j}^{2}\pi\partial_{i}\partial_{k}^{2}\pi\partial_{i}\dot{\pi}}{a^{3}}-\frac{\partial_{j}\pi\partial_{i}^{2}\partial_{j}\pi\partial_{k}^{2}\dot{\pi}}{a^{3}}+18\dot{\pi}^{3}a^{3}H^{4}$
(A.2)

##  Appendix B More Details about the Numerical Calculations of Bispectrum and
the Variable Transformations

To calculate different contributions to
$\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>$,
we first substituted $\zeta$ for $\pi$ using $\pi\rightarrow\zeta=-H\pi$ and
then transformed from cosmic time to conformal time $dt\rightarrow a\;d\tau$.
For the $\dot{\pi}(\partial\pi)^{2}$ term, this results in factoring factors
of $H$ and $a$ into the coefficients in the following way,

$C_{\dot{\pi}(\partial\pi)^{2}}\dot{\pi}(\partial\pi)^{2}\xrightarrow{}-\frac{1}{H^{3}}C_{\dot{\pi}(\partial\pi)^{2}}\dot{\zeta}(\partial\zeta)^{2}\xrightarrow{}-\frac{1}{aH^{3}}C_{\dot{\pi}(\partial\pi)^{2}}\zeta^{\prime}(\partial\zeta)^{2}.$
(B.1)

Therefore, we took track of these factors by defining

$C_{\zeta^{\prime}(\partial\zeta)^{2}}=-\frac{1}{aH^{3}}C_{\dot{\pi}(\partial\pi)^{2}}=-\frac{1}{H^{3}}[-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4}].$
(B.2)

After evaluating all the commutation relations between the operators in the
interaction terms and the external vertices, we found the contribution to the
bispectrum to be

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>_{\dot{\pi}(\partial\pi)^{2}}=-i\int_{-\infty}^{\tau_{f}}a(\tau^{\prime})d\tau^{\prime}C_{\zeta^{\prime}(\partial\zeta)^{2}}$
$\displaystyle\;2(2\pi)^{3}\delta^{3}(\textbf{k}_{3}+\textbf{k}_{2}+\textbf{k}_{1})\times$
$\displaystyle
f_{k_{1}}(\tau_{f}){f_{k_{1}}^{*}}^{\prime}(\tau^{\prime})f_{k_{2}}(\tau_{f}){f_{k_{2}}^{*}}(\tau^{\prime})$
$\displaystyle
f_{k_{3}}(\tau_{f}){f_{k_{3}}^{*}}(\tau^{\prime}){\textbf{k}}_{2}\cdot{\textbf{k}}_{3}+\text{\text{cycl.}}+\text{c.c.},$
(B.3)

where $f$ as we had defined in (3.19) represents the mode function for
$\zeta(\textbf{k})$. Next, we substitute for $f$ and $f^{\prime}$ in terms of
the canonical variable $u_{k}(\tau)$ which satisfies the Wronskian condition
and had been computed numerically from quadratic action. The substitution
relations are given by

$\displaystyle\sqrt{A_{1}}f$ $\displaystyle=-\frac{H}{a}u$
$\displaystyle\sqrt{A_{1}}f^{\prime}$
$\displaystyle=-\frac{H}{a}\frac{du}{d\tau}+H^{2}u$

to obtain

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>_{\dot{\pi}(\partial\pi)^{2}}$
$\displaystyle=$
$\displaystyle-i\int_{-\infty}^{\tau_{f}}a(\tau^{\prime})d\tau^{\prime}\left(\frac{H^{6}}{a(\tau_{f})^{3}a(\tau^{\prime})^{3}A_{1}^{3}}\right)\;C_{\zeta^{\prime}(\partial\zeta)^{2}}\;2(2\pi)^{3}\delta^{3}(\textbf{k}_{3}+\textbf{k}_{2}+\textbf{k}_{1})$
$\displaystyle
u_{k_{1}}(\tau_{f}){u_{k_{1}}^{*}}^{\prime}(\tau^{\prime})u_{k_{2}}(\tau_{f}){u_{k_{2}}^{*}}(\tau^{\prime})u_{k_{3}}(\tau_{f}){u_{k_{3}}^{*}}(\tau^{\prime}){\textbf{k}}_{2}\cdot{\textbf{k}}_{3}+\text{cycl.}+\text{c.c.}$
$\displaystyle+i\int_{-\infty}^{\tau_{f}}a(\tau^{\prime})d\tau^{\prime}\left(\frac{H^{7}}{a(\tau_{f})^{3}a(\tau^{\prime})^{2}A_{1}^{3}}\right)C_{\zeta^{\prime}(\partial\zeta)^{2}}\;2(2\pi)^{3}\delta^{3}(\textbf{k}_{3}+\textbf{k}_{2}+\textbf{k}_{1})$
$\displaystyle
u_{k_{1}}(\tau_{f}){u_{k_{1}}^{*}}(\tau^{\prime})u_{k_{2}}(\tau_{f}){u_{k_{2}}^{*}}(\tau^{\prime})u_{k_{3}}(\tau_{f}){u_{k_{3}}^{*}}(\tau^{\prime}){\textbf{k}}_{2}\cdot{\textbf{k}}_{3}+\text{cycl.}+\text{c.c.}$
(B.4)

Next, changing the integration variable to the dimensionless variable
$x_{i}=c_{s}k_{i}\tau$, and substituting for $a=-1/H\tau=-c_{s}k_{i}/Hx_{i}$
and defining $\tilde{u}(x)=\sqrt{c_{s}k}\,u_{k}(\tau)$151515This is due to the
symmetry of our ODE and normalization under this change of variable.  as well
as $x_{f}=c_{s}k_{1}\tau_{f}$, we arrive at.

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>_{\dot{\pi}(\partial\pi)^{2}}=$
$\displaystyle
i2(2\pi)^{3}\delta^{3}(\textbf{k}_{3}+\textbf{k}_{2}+\textbf{k}_{1})\frac{1}{\sqrt{c_{s}^{6}k_{1}^{2}k_{2}^{2}k_{3}^{2}}}\frac{H^{11}x_{f}^{3}}{(c_{s}k_{1})^{5}A_{1}^{3}}\left[\int_{-\infty}^{x_{f}}(x_{1}^{\prime})^{2}dx_{1}^{\prime}\;\right.$
$\displaystyle
C_{\zeta^{\prime}(\partial\zeta)^{2}}\;\tilde{u}(x_{f}){\tilde{u}^{*}}_{,x}(x_{1}^{\prime})\tilde{u}(\frac{k_{2}}{k_{1}}x_{f}){\tilde{u}^{*}}(\frac{k_{2}}{k_{1}}x_{3}^{\prime})\tilde{u}(\frac{k_{3}}{k_{1}}x_{f}){\tilde{u}^{*}}(\frac{k_{2}}{k_{1}}x_{1}^{\prime})$
$\displaystyle{\textbf{k}}_{2}\cdot{\textbf{k}}_{3}+\int_{-\infty}^{x_{f}}(x_{1}^{\prime})dx_{1}^{\prime}\;C_{\zeta^{\prime}(\partial\zeta)^{2}}\tilde{u}(x_{f}){\tilde{u}^{*}}(x_{1}^{\prime})\tilde{u}(\frac{k_{2}}{k_{1}}x_{f}){\tilde{u}^{*}}(\frac{k_{2}}{k_{1}}x_{1}^{\prime})$
$\displaystyle\tilde{u}(\frac{k_{3}}{k_{1}}x_{f}){\tilde{u}^{*}}(\frac{k_{3}}{k_{1}}x_{1}^{\prime}){\textbf{k}}_{2}\cdot{\textbf{k}}_{3}\bigg{]}+\text{cycl.}+\text{c.c.}$
(B.5)

Applying the Cosine Law,

$k_{1}^{2}=k_{2}^{2}+k_{3}^{2}-2k_{3}k_{2}\cos(\vartheta)\Rightarrow\textbf{k}_{2}\cdot\textbf{k}_{3}=k_{3}k_{2}\cos(\vartheta)=-\frac{k_{1}^{2}-k_{3}^{2}-k_{2}^{2}}{2},$
(B.6)

and substituting for $C_{\zeta^{\prime}(\partial\zeta)^{2}}$ lead to

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>_{\dot{\pi}(\partial\pi)^{2}}=$
$\displaystyle
i(2\pi)^{3}\delta^{3}(\textbf{k}_{3}+\textbf{k}_{2}+\textbf{k}_{1})\frac{(k_{1}^{2}-k_{3}^{2}-k_{2}^{2})}{\sqrt{c_{s}^{6}k_{1}^{2}k_{2}^{2}k_{3}^{2}}}\frac{H^{8}x_{f}^{3}}{(c_{s}k_{1})^{5}A_{1}^{3}}$
$\displaystyle[-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4}]\left[\int_{-\infty}^{x_{f}}(x_{1}^{\prime})^{2}dx_{1}^{\prime}\;\right.\tilde{u}(x_{f}){\tilde{u}^{*}}_{,x}(x_{1}^{\prime})$
$\displaystyle\tilde{u}(\frac{k_{2}}{k_{1}}x_{f}){\tilde{u}^{*}}(\frac{k_{2}}{k_{1}}x_{3}^{\prime})\tilde{u}(\frac{k_{3}}{k_{1}}x_{f}){\tilde{u}^{*}}(\frac{k_{3}}{k_{1}}x_{1}^{\prime})$
$\displaystyle+\int_{-\infty}^{x_{f}}x_{1}^{\prime}dx_{1}^{\prime}\;\tilde{u}(x_{f}){\tilde{u}^{*}}(x_{1}^{\prime})\tilde{u}(\frac{k_{2}}{k_{1}}x_{f}){\tilde{u}^{*}}(\frac{k_{2}}{k_{1}}x_{1}^{\prime})$
$\displaystyle\tilde{u}(\frac{k_{3}}{k_{1}}x_{f}){\tilde{u}^{*}}(\frac{k_{3}}{k_{1}}x_{1}^{\prime})\bigg{]}+\text{cycl.}+\text{c.c.}$
(B.7)

Now by defining parameters $\lambda\equiv\frac{k_{2}}{k_{1}}$ and
$\theta\equiv\frac{k_{3}}{k_{1}}$ which determine the shape of the triangles
formed by $\textbf{k}_{1},~{}\textbf{k}_{2}$ and $\textbf{k}_{3}$ we can
rewrite the above expression as

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>_{\dot{\pi}(\partial\pi)^{2}}=$
$\displaystyle(2\pi)^{3}\delta^{3}(\textbf{k}_{3}+\textbf{k}_{2}+\textbf{k}_{1})\frac{H^{8}}{c_{s}^{8}A_{1}^{3}}[-\delta_{3}H^{4}+9\bar{M}_{4}H^{3}-2M_{2}^{4}]$
$\displaystyle\left[i\frac{(1-\lambda^{2}-\theta^{2})x_{f}^{3}}{k_{1}^{6}\lambda\theta}\left\\{\int_{-\infty}^{x_{f}}(x_{1}^{\prime})^{2}dx_{1}^{\prime}\;\right.\right.\tilde{u}(x_{f}){\tilde{u}^{*}}_{,x}(x_{1}^{\prime})$
$\displaystyle\tilde{u}(\lambda\;x_{f}){\tilde{u}^{*}}(\lambda\;x_{1}^{\prime})\tilde{u}(\theta\;x_{f}){\tilde{u}^{*}}(\theta\;x_{1}^{\prime})$
$\displaystyle+\int_{-\infty}^{x_{f}}(x_{1}^{\prime})dx_{1}^{\prime}\;\tilde{u}(x_{f}){\tilde{u}^{*}}(x_{1}^{\prime})\tilde{u}(\lambda\;x_{f}){\tilde{u}^{*}}(\lambda\;x_{1}^{\prime})$
$\displaystyle\tilde{u}(\theta\;x_{f}){\tilde{u}^{*}}(\theta\;x_{1}^{\prime})+\text{c.c.}\bigg{\\}}\bigg{]}+\text{cycl.}$
(B.8)

Similar procedure was performed for evaluating the contributions from the
$\partial^{2}\pi(\partial\pi)^{2}$ term.

##  Appendix C Calculations of Bispectrum for the Local-type Non-Gaussianity

For the local-type NG, we characterise the size of the non-Gaussianty by
parameter $f_{{}_{\rm NL}}$ as

$\zeta(x)=\zeta_{{}_{\rm L}}(x)+\frac{3}{5}f_{{}_{\rm NL}}\zeta_{{}_{\rm
L}}^{2}(x),$ (C.1)

with

$\left<\zeta_{{}_{\rm L}}(x)\right>=0.$ (C.2)

The Fourier transform of $\zeta_{{}_{\rm L}}$ is then obtained by

$\displaystyle\zeta(\textbf{k})$ $\displaystyle=$ $\displaystyle\int
d^{3}x\;e^{-ikx}\zeta(x)=\int d^{3}x\;e^{-ikx}\left(\zeta_{{}_{\rm
L}}(x)+\frac{3}{5}f_{{}_{\rm NL}}\zeta_{{}_{\rm L}}^{2}(x)\right)$ (C.3)
$\displaystyle=$ $\displaystyle\int
d^{3}x\;e^{-ikx}\left(\frac{1}{(2\pi)^{3}}\int
d^{3}\tilde{k}\;e^{i\tilde{k}x}\zeta_{{}_{\rm
L}}(\tilde{\textbf{k}})+\frac{3}{5}f_{{}_{\rm NL}}\frac{1}{(2\pi)^{6}}\int
d^{3}\tilde{k}\;e^{i\tilde{k}x}\zeta_{{}_{\rm L}}(\tilde{\textbf{k}})\int
d^{3}k^{\prime}\;e^{i{k^{\prime}}x}\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right)$ $\displaystyle=$ $\displaystyle\left(\int
d^{3}\tilde{k}\;\delta^{3}(\textbf{k}-\tilde{\textbf{k}})\zeta_{{}_{\rm
L}}(\tilde{\textbf{k}})+\frac{3}{5}f_{{}_{\rm NL}}\frac{1}{(2\pi)^{3}}\int
d^{3}\tilde{k}d^{3}k^{\prime}\;\delta^{3}(\textbf{k}-\textbf{k}^{\prime}-\tilde{\textbf{k}})\zeta_{{}_{\rm
L}}(\tilde{\textbf{k}})\;\zeta_{{}_{\rm L}}(\textbf{k}^{\prime})\right)$
$\displaystyle=$ $\displaystyle\zeta_{{}_{\rm
L}}(\textbf{k})+\frac{3}{5}f_{{}_{\rm NL}}\frac{1}{(2\pi)^{3}}\int
d^{3}k^{\prime}\;\zeta_{{}_{\rm
L}}(\textbf{k}-\textbf{k}^{\prime})\;\zeta_{{}_{\rm L}}(\textbf{k}^{\prime})$

Substituting this for $\zeta(\textbf{k})$ in three-point function leads to

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>=$
$\displaystyle\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\zeta_{{}_{\rm L}}(\textbf{k}_{3})\right>$
$\displaystyle+\frac{3}{5}f_{{}_{\rm
NL}}\frac{1}{(2\pi)^{3}}\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\int d^{3}k^{\prime}\;\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\;\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>+\dots$ (C.4)

Using Isserlis’ theorem, the first term becomes

$\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\zeta_{{}_{\rm
L}}(\textbf{k}_{3})\right>=\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{1})\right>\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\zeta_{{}_{\rm L}}(\textbf{k}_{3})\right>=\dots=0.$ (C.5)

The second term can also be calculated using Wick’s theorem

$\displaystyle\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\int d^{3}k^{\prime}\;\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\;\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>=$ $\displaystyle\int
d^{3}k^{\prime}\;\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\;\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>$ $\displaystyle=$ $\displaystyle\int
d^{3}k^{\prime}\;\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\right>\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\;\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>$ $\displaystyle+\int
d^{3}k^{\prime}\;\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\right>$ $\displaystyle+\int
d^{3}k^{\prime}\;\left<\zeta_{{}_{\rm L}}(\textbf{k}_{2})\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\right>$ (C.6)

Next, taking into account that

$\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\right>=(2\pi)^{3}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2})P_{{}_{\zeta}}(k_{1})$
(C.7)

and

$\left<\zeta_{{}_{\rm L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\;\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>=(2\pi)^{3}\delta^{3}(\textbf{k}_{3})P_{{}_{\zeta}}(k^{\prime}).$
(C.8)

the first term reduced to

$\int d^{3}k^{\prime}\;\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\right>\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\;\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>=\int
d^{3}k^{\prime}\;(2\pi)^{6}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2})\delta^{3}(\textbf{k}_{3})P_{{}_{\zeta}}(k_{1})P_{{}_{\zeta}}(k^{\prime})$
(C.9)

which does not contribute to three point function for $k_{3}\neq 0$. The
second and the third terms are the same as

$\displaystyle\int d^{3}k^{\prime}\;\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\right>=$ $\displaystyle\int
d^{3}k^{\prime}\;(2\pi)^{6}\delta^{3}(\textbf{k}_{1}+\textbf{k}^{\prime})P(k_{1})\delta^{3}(\textbf{k}_{2}+\textbf{k}_{3}-\textbf{k}^{\prime})P_{{}_{\zeta}}(k_{2})$
$\displaystyle=$
$\displaystyle(2\pi)^{6}\delta^{3}(\textbf{k}_{2}+\textbf{k}_{3}+\textbf{k}_{1})P_{{}_{\zeta}}(k_{1})P_{{}_{\zeta}}(k_{2}),$
$\displaystyle\int d^{3}k^{\prime}\;\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>\left<\zeta_{{}_{\rm
L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\right>=$ $\displaystyle\int
d^{3}k^{\prime}\;(2\pi)^{6}\delta^{3}(\textbf{k}_{2}+\textbf{k}^{\prime})P(k_{2})\delta^{3}(\textbf{k}_{1}+\textbf{k}_{3}-\textbf{k}^{\prime})P(k_{1})$
$\displaystyle=$
$\displaystyle(2\pi)^{6}\delta^{3}(\textbf{k}_{2}+\textbf{k}_{3}+\textbf{k}_{1})P_{{}_{\zeta}}(k_{1})P_{{}_{\zeta}}(k_{2}).$

Therefore, adding everything we get

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>$
$\displaystyle=+\frac{3}{5}f_{{}_{\rm
NL}}\frac{1}{(2\pi)^{3}}\left<\zeta_{{}_{\rm L}}(\textbf{k}_{1})\zeta_{{}_{\rm
L}}(\textbf{k}_{2})\int d^{3}k^{\prime}\;\zeta_{{}_{\rm
L}}(\textbf{k}_{3}-\textbf{k}^{\prime})\;\zeta_{{}_{\rm
L}}(\textbf{k}^{\prime})\right>+{\rm sym.}$
$\displaystyle=+\frac{6}{5}(2\pi)^{3}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})f_{{}_{\rm
NL}}\bigg{(}P_{{}_{\zeta}}(k_{1})P_{{}_{\zeta}}(k_{2})+P_{{}_{\zeta}}(k_{3})P_{{}_{\zeta}}(k_{2})+P_{{}_{\zeta}}(k_{1})P_{{}_{\zeta}}(k_{3})\bigg{)}.$
(C.10)

Substituting for dimensionless power spectrum,
$\Delta_{{}_{\zeta}}(k)=P_{{}_{\zeta}}(k)\frac{k^{3}}{2\pi^{2}}$ we can
express above expression as:

$\displaystyle\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>=\frac{3}{10}(2\pi)^{7}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})f_{{}_{\rm
NL}}\bigg{(}\frac{\Delta_{{}_{\zeta}}(k_{1})\Delta_{{}_{\zeta}}(k_{2})}{k_{1}^{3}k_{2}^{3}}+\frac{\Delta_{{}_{\zeta}}(k_{3})\Delta_{{}_{\zeta}}(k_{2})}{k_{3}^{3}k_{2}^{3}}+\frac{\Delta_{{}_{\zeta}}(k_{1})\Delta_{{}_{\zeta}}(k_{3})}{k_{3}^{3}k_{1}^{3}}\bigg{)}$

If scale dependence of $\Delta_{{}_{\zeta}}(k)$ is negligible, $|n_{s}-1|\ll
1$ then we get:

$\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>=\frac{3}{10}(2\pi)^{7}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})f_{{}_{\rm
NL}}~{}[\Delta_{{}_{\zeta}}(k_{1})]^{2}\bigg{(}\frac{k_{1}^{3}+k_{2}^{3}+k_{3}^{3}}{k_{1}^{3}k_{2}^{3}k_{3}^{3}}\bigg{)}\,.$

After reordering

$\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>=(2\pi)^{7}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})\frac{[\Delta_{{}_{\zeta}}(k_{1})]^{2}}{k_{1}^{3}k_{2}^{3}k_{3}^{3}}\times\left\\{+\frac{3}{10}f_{{}_{\rm
NL}}~{}\bigg{(}k_{1}^{3}+k_{2}^{3}+k_{3}^{3}\bigg{)}\right\\}$

Therefore, in the case of (C.1) when collecting different contributions to
bispectrum as

$\left<\zeta(\textbf{k}_{1})\zeta(\textbf{k}_{2})\zeta(\textbf{k}_{3})\right>=(2\pi)^{3}\delta^{3}(\textbf{k}_{1}+\textbf{k}_{2}+\textbf{k}_{3})B(k_{1},k_{2},k_{3})$
(C.11)

we get

$f_{{}_{\rm
NL}}=\frac{1}{(2\pi)^{4}\Delta_{{}_{\zeta}}^{2}}B(k_{1},k_{2},k_{3})\left\\{\frac{10}{3}\frac{\prod
k_{i}^{3}}{\sum k_{i}^{3}}\right\\}.$ (C.12)

## References

  * [1] A. D. Linde, ETERNAL CHAOTIC INFLATION, Mod. Phys. Lett. A 1 (1986) 81.
  * [2] A. H. Guth, Eternal inflation and its implications, Journal of Physics A: Mathematical and Theoretical 40 (Jun, 2007) 6811–6826.
  * [3] V. F. Mukhanov and G. V. Chibisov, Quantum Fluctuations and a Nonsingular Universe, JETP Lett. 33 (1981) 532–535.
  * [4] I. Babic, C. P. Burgess, and G. Geshnizjani, Keeping an eye on DBI: power-counting for small-cs cosmology, JCAP 05 (2020) 023, [1910.05277].
  * [5] C. P. Burgess, Intro to Effective Field Theories and Inflation, 1711.10592.
  * [6] J. L. Kim and G. Geshnizjani, Spectrum of Cuscuton Bounce, JCAP 03 (2021) 104, [2010.06645].
  * [7] S. S. Boruah, H. J. Kim, M. Rouben, and G. Geshnizjani, Cuscuton bounce, JCAP 08 (2018) 031, [1802.06818].
  * [8] A. Borde and A. Vilenkin, Eternal inflation and the initial singularity, Phys. Rev. Lett. 72 (1994) 3305–3309, [gr-qc/9312022].
  * [9] A. Borde, A. H. Guth, and A. Vilenkin, Inflationary space-times are incompletein past directions, Phys. Rev. Lett. 90 (2003) 151301, [gr-qc/0110012].
  * [10] P. Creminelli, M. A. Luty, A. Nicolis, and L. Senatore, Starting the Universe: Stable Violation of the Null Energy Condition and Non-standard Cosmologies, JHEP 12 (2006) 080, [hep-th/0606090].
  * [11] C. Cheung, P. Creminelli, A. Fitzpatrick, J. Kaplan, and L. Senatore, The Effective Field Theory of Inflation, JHEP 03 (2008) 014, [0709.0293].
  * [12] N. Bartolo, M. Fasiello, S. Matarrese, and A. Riotto, Large non-gaussianities in the effective field theory approach to single-field inflation: the bispectrum, Journal of Cosmology and Astroparticle Physics 2010 (Aug, 2010) 008–008.
  * [13] A. Ashoorioon, R. Casadio, G. Geshnizjani, and H. J. Kim, Getting Super-Excited with Modified Dispersion Relations, JCAP 09 (2017) 008, [1702.06101].
  * [14] A. Ashoorioon, R. Casadio, M. Cicoli, G. Geshnizjani, and H. J. Kim, Extended Effective Field Theory of Inflation, JHEP 02 (2018) 172, [1802.03040].
  * [15] A. Ashoorioon, Non-Unitary Evolution in the General Extended EFT of Inflation \& Excited Initial States, JHEP 12 (2018) 012, [1807.06511].
  * [16] A. Ashoorioon, A. Rostami, and J. T. Firouzjaee, EFT compatible PBHs: effective spawning of the seeds for primordial black holes during inflation, JHEP 07 (2021) 087, [1912.13326].
  * [17] N. Bartolo, M. Fasiello, S. Matarrese, and A. Riotto, Large non-Gaussianities in the Effective Field Theory Approach to Single-Field Inflation: the Bispectrum, JCAP 08 (2010) 008, [1004.0893].
  * [18] A. Adams, N. Arkani-Hamed, S. Dubovsky, A. Nicolis, and R. Rattazzi, Causality, analyticity and an IR obstruction to UV completion, JHEP 10 (2006) 014, [hep-th/0602178].
  * [19] A. Ashoorioon, D. Chialva, and U. Danielsson, Effects of Nonlinear Dispersion Relations on Non-Gaussianities, JCAP 06 (2011) 034, [1104.2338].
  * [20] C. P. Burgess and M. Williams, Who You Gonna Call? Runaway Ghosts, Higher Derivatives and Time-Dependence in EFTs, JHEP 08 (2014) 074, [1404.2236].
  * [21] J. Garriga and A. Vilenkin, Living with ghosts in Lorentz invariant theories, JCAP 01 (2013) 036, [1202.1239].
  * [22] N. Bartolo, M. Fasiello, S. Matarrese, and A. Riotto, Large non-Gaussianities in the Effective Field Theory Approach to Single-Field Inflation: the Trispectrum, JCAP 09 (2010) 035, [1006.5411].
  * [23] R. Holman and A. J. Tolley, Enhanced Non-Gaussianity from Excited Initial States, JCAP 05 (2008) 001, [0710.1302].
  * [24] A. Ashoorioon and G. Shiu, A Note on Calm Excited States of Inflation, JCAP 03 (2011) 025, [1012.3392].
  * [25] I. Agullo and L. Parker, Non-gaussianities and the Stimulated creation of quanta in the inflationary universe, Phys. Rev. D 83 (2011) 063526, [1010.5766].
  * [26] A. Ashoorioon, K. Dimopoulos, M. M. Sheikh-Jabbari, and G. Shiu, Reconciliation of High Energy Scale Models of Inflation with Planck, JCAP 02 (2014) 025, [1306.4914].
  * [27] L. Senatore, K. M. Smith, and M. Zaldarriaga, Non-Gaussianities in Single Field Inflation and their Optimal Limits from the WMAP 5-year Data, JCAP 01 (2010) 028, [0905.3746].
  * [28] E. Komatsu and D. N. Spergel, Acoustic signatures in the primary microwave background bispectrum, Physical Review D 63 (Feb, 2001).
  * [29] Y. A. et al., Planck 2018 results. ix. constraints on primordial non-gaussianity, 2019.
  * [30] C. Armendariz-Picon, T. Damour, and V. F. Mukhanov, k - inflation, Phys. Lett. B 458 (1999) 209–218, [hep-th/9904075].
  * [31] M. Alishahiha, E. Silverstein, and D. Tong, DBI in the sky, Phys. Rev. D 70 (2004) 123505, [hep-th/0404084].
  * [32] X. Chen, M.-x. Huang, S. Kachru, and G. Shiu, Observational signatures and non-Gaussianities of general single field inflation, JCAP 01 (2007) 002, [hep-th/0605045].
  * [33] P. D. Meerburg, J. P. van der Schaar, and P. S. Corasaniti, Signatures of Initial State Modifications on Bispectrum Statistics, JCAP 05 (2009) 018, [0901.4044].