Non-Gaussianities in the Extended EFT of Inflation

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.


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 [10][11][12], explored their implications for inflation up to the order of ω 2 ∝ k 4 . Nevertheless, terms leading to ω 2 ∝ k n with n ≥ 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 1 . 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 ω 2 ∝ k 6 at early times. In fact since the observation of the fluctuations in the cosmic microwave background (CMB) and the measurements of f 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 2 . Furthermore, if the primordial non-Gaussianities are eventually observed, calculating the non-Gaussianity 1 This 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 ∆ ζ ∼ 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.
2 For study of non-Gaussianity in the EFToI refer to [11,12] 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).

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 µν 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] where L m represents the spatially invariant terms of perturbative order m in arguments g 00 + 1, δK µν and δR µνρσ . 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 [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 → t + ξ 0 (x µ ) is taken and then time diffeomorphisms was restored by substituting ξ 0 (x µ ) with a field −π(x µ ). Now asserting that under time diff π → π − ξ 0 (x µ ), 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 π and taking the decoupling limit whereḢ → 0 while M 2 PlḢ is kept fixed, leading to the following Lagrangian density for π 3 , One can substitute π with the more familiar conserved quantity ζ = −Hπ to obtain the standard slow-roll inflationary action for ζ. Next, including more terms with m ≥ 2 and operators up to mass order two (d = 2) a larger class of inflationary scenarios were identified and studied [11], Here, the upper indices are exponents of M i andM 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 ≤ 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 ∼ 2π, 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,M 2 ∇ µ g 00 ∇ µ g 00 which after the implementation of the Stueckelberg technique yields Ostragradski ghosts at second order for the Goldstone boson, π. That operator in the unitary gauge action leads tō 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 4 . The Lagrangian of the EFToI (2.3) in the decoupling limit, whereḢ → 0 and M 2 PlḢ is kept fixed, results into the following quadratic and cubic expressions for π [11,17]: First, note that any term that contributes to the coefficient ofπ 2 in the approximate plain wave regime, leads to deviations from c S = 1. However, any terms that results in spatial derivatives for π, can potentially modify the dispersion relation. Some modifications could 4 The second term, k 2 a 2π 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. even lead to non-hyperbolic equations of motion for π that do not have propagating solutions or are not even well posed. As we see, the first term in above expression, M 4 2 2! (g 00 + 1) 2 , only contributes to modifying c S [11], This shows that including only this term, the superluminal propagation can be avoided for M 4 2 > 0, since during inflationḢ < 0 5 . Through a quick substitution for M 2 , we see the cubic term in π from this same operator, appears with a factor proportional to (1 − 1 ). 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 toM 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 ω 2 = c 2 S k 2 but contributes to changes in speed of sound. In the exact de-Sitter limit,Ḣ = 0, with M 2 = 0 andM 1 ∼ −M 2 ∼M , the speed of sound turns out to be quite small, c 2 S = H

4M
1, but if M 2 = 0, M 3 < 0, then c 2 S = 1 6 . Finally, the last two terms would not only add corrections to the sound speed for π 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 [17,19].
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 ω 2 ∝ k 2n with n ≥ 3, the power of the energy scaling of π would be −1/2 + 3/2n which means the scaling dimension of the cubic operatoṙ π(∂ i π) 2 would be (7 − 3n)/(2n). For n ≥ 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 ω 2 ∝ 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 (∂ 2 i π) 2 and (∂ i π) 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 ω 2 ∼ k 6 , the scaling power ofπ(∂ i π) 2 is −1/3, but the scaling power of (∂ 2 i π) 2 and (∂ i π) 2 are respectively −2/3 and −4/3. This implies that the dispersion relation would change from ω 2 ∼ k 6 to ω 2 ∼ k 4 or ω 2 ∼ k 2 at some energy scale Λ 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 ω 2 ∝ 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: The additional terms up to order d = 4 that were included in EEFToI are given by Here, δ 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: (2.10) As we see δ 1 and δ 2 expressions contain (π) 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 δ 1 = δ 2 = 0 6 7 .
6 Setting δ1 +3δ2 = 0, one could in principle avoid theπ 2 in the action for the Goldstone boson, π. The term proportional to δ1, (∇µδK νγ )(∇ µ δKνγ), 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 δ1 = δ2 = 0 7 One can show that if ∇ µ operators are contracted with the induced metric, h µν , there will be no ghosts 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 δg 00 and δ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 (∇ µ ∇ ν δg 00 )δK µν , and δg 00 (∇ µ ∇ ν δK µν ) are absorbed into theM 4 term. Finally, there are more terms such as (∇ µ ∇ µ δg 00 ) 2 , (∇ µ ∇ µ δg 00 )Tr(δK) that have been discarded as they also result in (π) 2 in the quadratic Lagrangian for π 8 . 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 (π) 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π2 k and k 4π2 k and the potential strong coupling regimes need to be studied carefully [15].
Going forward we also only explore models where δ 3 = −3δ 4 . The reason for this additional assumption is that we are more interested in investigating models where the effects of k 6 π 2 k operators is dominant and setting δ 3 = −3δ 4 will automatically cancel off k 2π2 k terms. Deviating from this condition, one has to carefully take into account the competing impact of k 2π2 k 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: where the the coefficients A 1 , C 1 , D 1 and F 1 are defined as Next, expanding the Canonical variable a √ A 1 π k in terms of the corresponding mode functions u k , and computing the equation of motion for u k yields, where the primes denote derivatives with respect to the conformal time τ such that f = d dτ f = a d dt f = aḟ . Note that in this equation for the positive slow varying values of F 1 /A 1 , in the regime where F 1 A 1 k 2 is the dominant term (ω ∝ k) the speed of sound for scalar 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 ∇µ are allowed as they are invariant under spatial diffeomorphisms [14]. 8 Similarly one can avoid the ghosts in these operators by contracting the ∇µ operators with the induced metric, rather then the four dimensional metric.
perturbations is well defined and real, with (2.14) In the de Sitter limit where a = −1/(Hτ ) and defining y = c S kτ andũ(y) = c S k u k (τ ) to preserve the normalization for Wronskian condition, we can rewrite equation(2.13) as (2.16) We can solve the mode equation (2.15) numerically, taking the positive frequency WKB mode in the infinite past,ũ (2.18) Since the dimensionless power spectrum is related to we also get As can be expected, value of γ S depends on values of α 0 and β 0 . Under the assumption that the dispersion relation does not become tachyonic before horizon crossing (i.e. 1 + α 0 y 2 + β 0 y 4 > 0), it was shown in [13] that the greatest enhancements to the Bunch-Davies power spectrum is obtained when α 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 α 0 = −0.18 and β 0 = α 2 0 /4. For α 0 < 0 and β 0 ≤ α 2 0 /4, the dispersion relation can become tachyonic, allowing for larger enhancements of the power spectrum [14,16]. For α 0 , β 0 ≥ 0, the scalar power spectrum get suppressed with respect to its Bunch-Davies counterpart, γ S 1. This observation along with the fact that for pure quartic positive correction to the dispersion relation, γ 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 ≤ α 0 , β 0 1, γ S ≈ 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 α 0 and β 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.

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 ω 2 ∝ 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 ζ k as and the three point function in the interaction picture can be computed as where in these calculations, π k can be substituted for ζ k using the relation ζ k = −Hπ 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 ω 2 ∝ k 6 . The kinetic term in the action (2.11) in real space has the following form for π Rescaling energy by a factor of Σ, E → ΣE, means that t → Σ −1 t, which given the dispersion relation ω 2 ∝ k 6 , implies x → Σ −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, π must remain scale-invariant as well, π → π. We now check the scaling of a generic term in action, dt d 3 x M m,s,p where m denotes the number of time derivatives, s number of spatial derivatives and p the overall powers of π. Since x → Σ −1/3 x, the spatial derivative goes as ∂ i → Σ 1/3 ∂ i so the contribution to action scales as M m,s,p dt dx 3 → Σ −2+m+s/3 M m,s,p dt d 3 x. When m + s/3 < 2, this term diverges in the low energy limit of Σ → 0 and becomes relevant. Note the powers of π 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 ≤ 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,

5)
m + s/3 < 2 . is finite and excludes the existence of any relevant operators beyond quartic terms as well 9 . Also note that M 1,2,2 ∝ ∂π∂ i π = 1 2 ∂ ∂t (∂ i π) 2 term becomes M 0,2,2 ∝ (∂ i π) 2 after integrating by parts, which can be absorbed into the canonical kinetic term. The scaling power of M 0,2,2 ∝ (∂ i π) 2 and M 0,4,2 ∝ (∂ 2 i π) 2 which modify the dispersion relation, as we had pointed out earlier are respectively −4/3 and −2/3. The scaling power of M 1,2,3 ∝π(∂ i π) 2 is −1/3 while both M 0,4,2 ∝ ∂ 2 j π(∂ i π) 2 and M 0,4,4 ∝ (∂ i π) 4 scale as −2/3 as well. The term M 0,4,4 ∝ (∂ i π) 4 contributes to trispectrum which we will not discuss in this paper 10 . Therefore, to calculate the strongest contributions to bispectrum we only need to consider operators (0, 4, 3) and (1, 2, 3), which are Note that both of these terms are also produced in (2.5) from the terms included in the EFToI as well: (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 δ 3 = −3δ 4 are the following In the limit where the parameters of EFToI and EEFToI satisfy and 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 ω ∼ H and k ∼ H/c S and approximatingπ ∼ ωπ, the ratio of contributions of the cubic term (∂ i π) 2π to the quadratic term can be estimated as and similarly for ∂ 2 j π(∂ i π) 2 , we get 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 ζ∆ ζ ∼ 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 π, and we are interested in the comoving curvature perturbation, ζ = −Hπ, and its bispectrum < ζ 3 >. Therefore, we transform the terms H int (π) → H int (ζ) and also for convenience switch to conformal time τ instead of t. As we pointed out earlier, the most dangerous contributions from H int or L int are coming from termsπ(∂ i π) 2 and ∂ 2 i π(∂ j π) 2 . In terms of ζ, these terms becomė Next, similar to calculating power spectrum, ζ k is written as a quantum field operator in terms of the creation and annihilation operators a † k , a −k satisfying the canonical commutation relations [a k , a † k ] = (2π) 3 δ 3 (k − k ) and the mode functions f k (τ ): In practice, we carried out the calculations numerically in terms of the modes functions u(c S kτ ) = − √ c S kA 1 a H f k (τ ) which we had computed previously from quadratic action and performed the integration in terms of the dimensionless variable x i = c S k i τ 11 . 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π(∂π) 2 term to ζ(k 1 )ζ(k 2 )ζ(k 3 ) is: where we are defining parameters λ ≡ k 2 k 1 and θ ≡ k 3 k 1 that determine the shape of the triangles formed by k 1 , k 2 and k 3 .
2. contributions from the ∂ 2 π(∂π) 2 term to ζ(k 1 )ζ(k 2 )ζ(k 3 ) is: 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 11 See appendix (B) for more details.
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 2 1 k 2 2 k 2 3 , evaluated at k 1 = 1 and normalized over B I (1, 1, 1), for various triangle shapes. As we see the bispectrum (a) The shape of bispectrum for operatorπ(∂π) 2 (b) The shape of bispectrum for operatorπ(∂π) 2 Figure 2: Calculation of the quantity B I (k 1 , k 2 , k 3 )k 2 1 k 2 2 k 2 3 , evaluated at k 1 = 1 and normalized over B I (1, 1, 1) for various triangle shapes. This result corresponds to value of α 0 = −0.18 and β 0 = α 2 0 /4. 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 12 . 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 inq i and other variables determining B I . We will further explore these bounds and discuss our numerical calculations in the next section.

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 NL in the literature.
Historically, f 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 where ζ L is the gaussian scalar perturbation. Now under this assumption as well as near scale invariance (i.e. n S − 1 1), one can obtain the following relationship between the three point function in k space and the dimensionless power spectrum of ζ 13 , 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), 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 local NL for different triangles in momentum spaces, Therefore, finding the observed value of f local NL 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 → 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 data 14 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 ) . (4.5) 13 See appendix (C) for details. 14 Readers are referred to [29] for details on different templates fitting for different models.
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 [23][24][25][26]32]. 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 [23,32] is which diverges as k 1 → k 2 + k 3 . Another template that has been used for testing EFToI models is the orthogonal template defined as For our purposes here, since we want to make only order of magnitude estimations and since f 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 (k 1 , k 2 , k 3 ) from (3.22) and that ∆ ζ = γ S ∆ B.
Given that k 2 /k 1 = λ and k 3 /k 1 = θ the problem reduces to finding the best fit for the coefficients f (I)... (1, λ, θ). (4.11) 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ũ(y) where y = c S kτ . This computation was carried out for different values of α 0 and β 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 γ S was the largest (see Figure (1)) in the following region In fact, within that region the enhancement was found to be the greatest (γ S ≥ 14.6) when setting α 0 = −0.18 and β 0 = α 2 0 /4. From physics point of view, a more interesting applications of the EEFToI and EFToI is when they result in larger values of γ 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ũ(y) for different values of α 0 and β 0 were substituted in the equations (3.20) and (3.21) and then numerically integrated to obtain the contributions to bispectrum.
In order to set the initial conditions at x → −∞, we chose the value of x initial to be where β 0 x 4 is far greater than α 0 x 2 (x initial ∼ −50|α 0 /β 0 | 1/2 ). This guarantees that the initial condition is set in the regime where ω 2 ∝ 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 τ 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 ≤ λ ≤ 1 and 0.2 ≤ θ ≤ 1. The lower bounds for λ and θ should in principle either be set by ratio of the observed modes in CMB (∼ 10 3 ) or theoretically by validity of the effective field theory regime. Note, the extreme squeezed limit θ → 0 or λ → 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ũ(λx) also initiate in the same regime. In practice, setting 0.2 as the lower bound for both λ and θ was due to the limitation of our numerical method. That is, we chose the values of λ and θ 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.   corresponding to values of α 0 = −0.18 and β 0 = α 2 0 /4. 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   NL −38 ± 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 (I) NL γ 2 S /q I to infer bounds on value ofq 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 γ 2 S andq i . We require for each shape estimated f (I) NL be less than reported value, which leads to Therefore, we can crudely set a bound onq I 10 −2 − 10 −3 γ 2 S . In this case substituting larger than unity values such as γ S ∼ 14 results into a not very tight bound of |q I | 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 ( Regions where the power spectrum is enhanced. Note that for α 0 > 0 the power spectrum is suppressed (not interesting) β 0 > 0 To avoid tachyonic instabilities (note δ 3 > 0) ∆ ζ ∼ 10 −9 The constraint from power spectrum of the CMB observation Speed of sound for the tensor modes (Note:M 2 To avoid ghosts at the leading ordeṙ H < 0 Inflation satisfies the null energy conditioñ q i ≤ 10 −1 Constraints from f NL measurements Table 3: List of the constraints applied to the EEFToI models Furthermore, we set an additional cutoff at M Pl , except for δ 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. 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.

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 ω 2 ∝ 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π(∂ i π) 2 and ∂ 2 i π(∂ j π) 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 NL in terms of the parameters of the model. This lead to an order of magnitude estimation for the bounds on variables q I which include the masses and couplings in the model defined in (3.23) asq I 10 −1 . Finally, we used a search algorithm and found that there are allowed regions of parameter space where EEFToI description with ω 2 ∝ 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 (∂ i π) 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 ω 2 ∝ 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, δK µν , δR µνρσ , ∇ µ that through their quadratic or interaction contributing terms in the action, could potentially lead to interesting and distinctive observational signatures.

A The Interaction Lagrangian for the EEFToI
In the section we summarize all the cubic terms in the action d 4 x √ g 4 d=3 L 2,d generated through the new EEFToI terms (2.8) after setting δ 1 = δ 2 = 0.
The contributions from the operator proportional toM 4 to the cubic interactions in the action, taking into account √ g, goes as: The contributions from the operator proportional to the δ 3 to the cubic interaction terms in the action goes as: The contribution of the operator proportional to the δ 4 in Lagrangian to the cubic interactions is given by:

B More Details about the Numerical Calculations of Bispectrum and the Variable Transformations
To calculate different contributions to ζ(k 1 )ζ(k 2 )ζ(k 3 ) , we first substituted ζ for π using π → ζ = −Hπ and then transformed from cosmic time to conformal time dt → a dτ . For thė π(∂π) 2 term, this results in factoring factors of H and a into the coefficients in the following way, (B.1) Therefore, we took track of these factors by defining 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 where f as we had defined in (3.19) represents the mode function for ζ(k). Next, we substitute for f and f in terms of the canonical variable u k (τ ) which satisfies the Wronskian condition and had been computed numerically from quadratic action. The substitution relations are given by Now by defining parameters λ ≡ k 2 k 1 and θ ≡ k 3 k 1 which determine the shape of the triangles formed by k 1 , k 2 and k 3 we can rewrite the above expression as ζ(k 1 )ζ(k 2 )ζ(k 3 ) π(∂π) 2 =(2π) 3 δ 3 (k 3 + k 2 + k 1 ) u(θ x f )ũ * (θ x 1 ) + c.c. + cycl. (B.8) Similar procedure was performed for evaluating the contributions from the ∂ 2 π(∂π) 2 term.