Quasi Single Field Inflation in the non-perturbative regime

In quasi single field inflation there are massive fields that interact with the inflaton field. If these other fields are not much heavier than the Hubble constant during inflation ($H$) these interactions can lead to important consequences for the cosmological energy density perturbations. The simplest model of this type has a real scalar inflaton field that interacts with another real scalar $S$ (with mass $m$). In this model there is a mixing term of the form $\mu {\dot \pi} S$, where $\pi$ is the Goldstone fluctuation that is associated with the breaking of time translation invariance by the time evolution of the inflaton field during the inflationary era. In this paper we study this model in the region $(\mu/H )^2 +(m/H)^2>9/4$ and $m/H \sim {\cal O}(1)$ or less. For a large part of the parameter space in this region standard perturbative methods are not applicable. Using numerical and analytic methods we derive a number of new results. In addition we study how large $\mu/H$ has to be for the large $\mu/H$ effective field theory approach to be applicable.


I. INTRODUCTION
There is very strong evidence that the universe was once in a radiation dominated era followed by a matter dominated era. Today the universe is dominated by vacuum energy density and we are entering an inflationary era where the scale factor a(t) ∝ e H 0 t , with H 0 near the Hubble constant today. It is widely believed that at very early times there was another inflationary era where the energy density was dominated by false vacuum energy giving rise to a Robertson Walker scale factor with time dependence a(t) ∝ e Ht , where H is the Hubble constant during that inflationary era [1,2]. After more than about 60 e-folds, this inflationary era ends and the universe reheats to a radiation dominated (Robertson Walker) Universe. If this is the case then the horizon and flatness problems [2] can be solved and in addition there is an attractive mechanism based on quantum fluctuations for generating density perturbations with wavelengths that were once outside the horizon [3] (see Ref. [4] for a review of inflation). It has been argued that it requires tuning to enter the inflationary era [5] (see however [6]) and furthermore that there are issues with its predictability [7] (see also [8] for a recent discussion of these issues). Nevertheless, because of the simplicity of the dynamics of the inflationary universe paradigm and the ability within it to do explicit calculations of the properties of the cosmological energy density perturbations [3] and primordial gravitational waves [9], it seems worth studying particular inflationary models in some detail. The simplest inflationary model is standard slow roll inflation with only a single real scalar field, the inflaton φ(x). It is conventional to work in a gauge where fluctuations in the inflaton field about the classical slow roll solution φ 0 (t) vanish. Then using the Stückelberg arXiv:1706.09971v2 [hep-ph] 18 Dec 2017 trick the curvature fluctuations that are constant outside the horizon and become the density perturbations when they reenter the horizon (in the radiation and matter dominated eras) arise from quantum correlations in the Goldstone mode π(x) calculated during the de-Sitter inflationary era 1 . In this model non gaussianities in cosmological density correlations arise because of connected higher point correlations of π, but they are very small [11].
Larger non-gaussianities can be achieved if there are other fields with masses around or less than the inflationary Hubble constant 2 , that couple to π (see Ref. [14] for a review). In quasi-single field inflation these extra fields do not directly influence the classical evolution of the inflaton field but impact the cosmological density perturbations since they couple to the inflaton as "virtual particles" and hence affect the the correlations of π [15]. To simplify matters we will assume an approximate shift symmetry on the inflaton field, φ(x) → φ(x)+c (where c is a constant) that is only broken by the potential, V φ , for φ. Furthermore, we assume an unbroken discrete symmetry, φ(x) → −φ(x). The simplest quasi-single field model introduced by Chen and Wang [14] has a single additional (beyond the inflaton) real scalar field S. The Lagrange density in this model contains an unusual kinetic mixing of the form µπS .
Throughout this paper we treat µ as a constant independent of time. There has been a study of the case where µ changes suddenly with time, becoming large momentarily [27].
In this paper we focus on the region of parameter space where (µ/H) 2 + (m/H) 2 > 9/4 and m/H ∼ O (1) or less (recall m is the mass term for S). In this region, non-gaussianities have an interesting oscillatory behavior [19]. We use numerical non-perturbative methods similar to those developed in [26] and the effective field theory for large µ/H to study the model in this region of parameter space. We study how large µ/H must be for the effective field theory method to be quantitatively correct. In addition we derive the n S , r plot for the model with inflaton potential V φ = m 2 φ φ 2 /2 and derive the limit on µ/H and the S potential parameter V S from Planck limits on non-gaussianity.
In section II we discuss the Lagrange density of the model we use in detail. Section III reviews quantization of the free part of the Lagrange density in flat space-time. Even this theory is non-trivial because of the unusual Lorentz non-invariant kinetic mixing between the Goldstone field π and the excitations of the massive scalar S. The massless mode has an unusual energy momentum relation that, for momentum in the range m q µ, has a non-relativistic flavor, E q = q 2 /µ [24]. The other mode is heavy with a mass µ 2 + m 2 . The fact that this mode's mass does not go to zero as m → 0 is what regularizes the divergences that occur at m = 0 when one treats µ perturbatively.
Quantization of the free field theory in de-Sitter space-time is discussed in section IV. In de-Sitter space-time a mode's physical momentum q evolves with time. At early times modes have wavelengths much less than the horizon 1/H but at later times the wavelengths get red-shifted outside the horizon. The mode functions are calculated non-perturbatively by numerically solving the differential equations they satisfy in the region of parameter space, (µ/H) 2 + (m/H) 2 > 9/4 and m/H ∼ O(1) or less. Quantum fluctuations in the field S fall off rapidly for wavelengths outside the horizon and it is the quantum fluctuations in the field π that determine the curvature and density fluctuations just as in standard slow roll single field inflation. Nevertheless, these quantum fluctuations are influenced by π's couplings to S.
In section IV we analyze (in the non-perturbative regime) the curvature perturbation power spectrum in this model focusing on the transition between the perturbative regime and the regime where the effective theory applies. Section V derives the n S , r plot in this theory for the simple inflaton potential V φ = m 2 φ φ 2 /2. Non-gaussianities are discussed in Sec VI. We calculate the bispectrum in the the equilateral and squeezed configurations in the non-perturbative region numerically. In the large µ region we show that the numerical results agree with the results from the effective theory. We derive the constraints on µ/H and the S potential parameter V S from Planck limits on non-gaussianity.
In section VII we review the derivation of the effective field theory for large µ/H and the derivation of the power spectrum using it. We then compute the bispectrum in this effective field theory including a contribution from the potential for S that was not previously presented in the literature.
Our conclusions are summarized in Sec. VIII.

II. THE MODEL
The simplest quasi-single field inflation model has a real scalar inflaton field φ that interacts with another real scalar field S. We impose a φ → −φ symmetry and an approximate shift symmetry φ → φ + c, where c is a constant. The shift symmetry is only broken by the inflaton potential V φ (φ). The Lagrangian we use has the form Interactions between the inflaton φ and the massive field S first occur at dimension 5 and if we neglect operators with dimension higher than this the interaction Lagrangian is One natural choice for the mass scale Λ is the Planck mass. This higher dimensional operator would then arise from the transition from the theory of quantum gravity to a quantum field theory. In this case the non-gaussianities are very small. However, another possibility is that there is physics at a scale Λ that is large compared to the Hubble constant during inflation but well below the Planck scale. Integrating out this physics can give rise to such an operator. We work in a gauge where the inflaton field is only a function of time, φ(x) = φ 0 (t) and take the background metric to have the form, ds 2 = dt 2 − a(t) 2 dx 2 , with the scale factor a(t) = e Ht . The Goldstone boson associated with the time translation invariance breaking by the classical evolution φ 0 (t) is denoted by π(x). The curvature perturbation is proportional to this field, ζ = −Hπ. We expand S about a background classical value S(x) = S 0 + s(x) and assume that the background solution S 0 is independent of time. This assumption is consistent with the dynamical equations of evolution for the fields provided we neglect second time derivatives of φ 0 (t). With those assumptions φ 0 (t) and S 0 satisfy, The dynamics for the fluctuations π(x) and s(x) are controlled by the Lagrange density, where the free part of the Lagrange density for the fields π and s is, where m 2 = V (S 0 ). Throughout this paper we assume that the mass parameter m for the additional scalar s is of order the Hubble constant during inflation or smaller. The interaction part of the Lagrange density is It is convenient to introduce a rescaled π that has a properly normalized kinetic term, where,φ In terms of these rescaled fields the gravitational curvature perturbation becomes, (2.11) The free and interacting Lagrange densities, after introducing a redefined scaleΛ = (1 + 2S 0 /Λ)Λ, are and In eq. (2.12) we have introduced µ = 2φ 0 /Λ. (2.14) and in eq. (2.13) only explicitly kept those terms that play a role in the calculations performed in this paper. In the following sections we will drop the tilde on the Goldstone field π to simplify the notation. Moreover, we adopt sign conventions for φ and S so thatφ 0 and µ are positive. As mentioned in the introduction the purpose of this paper is to study this model in the region of parameter space where (µ 2 + m 2 ) 1/2 /H > 3/2 and m ∼ O(H) or smaller. Some of this region, i.e. where µ/H is small or very large have been previously studied. We will compare with those results to find out how small and how large µ/H has to be for the approximate methods used in those regions to be accurate.
First let's imagine that S 0 =0. This can always be arranged by tuning the linear term in the potential V S (S) to cancel the linear term in S from the 1/Λ interaction term. Then µ/H = (2φ 0 /H 2 )(H/Λ). The measured power spectrum for the curvature perturbations implies thatφ 0 /H 2 is very large so even for small H/Λ one can achieve large values for µ/H.
Next we allow a non zero S 0 but simplify the potential so it contains no terms with more than two powers of S, explicitly V S = V S S + m 2 S 2 /2. In this case µ/H can be written as, Therefore, without tuning the tadpole in V S to cancelφ 2 0 S/Λ, it is not possible to have the mass parameter m of order the Hubble constant (or smaller) and µ/H large. Nonetheless it seems worth studying this region of parameter space since there are some novel features that arise there.
Naive dimensional analysis suggests that higher dimension operators that couple derivatives of φ to a single S are smaller than the dimension 5 operator we kept provideḋ φ 0 /Λ 2 = (µ/H) 2 (H 2 /φ 0 ) < 1. The higher powers of S will be small if in addition S 0 /Λ < 1. Since the measured amplitude of the density perturbations implies that H 2 /φ 0 is quite small the ratio µ/H can be large in the region of parameter space where the operator expansion in powers of 1/Λ is justified. Indeed, comparing the calculated power spectrum at large µ/H given in (4.13) with it's measured value, the upper limit for µ/H for power counting in the 1/Λ expansion to be valid is µ/H < ∼ 300. Of course, this is just a naturalness constraint and can be violated without the model being inconsistent.

III. FREE FIELD THEORY IN FLAT SPACE-TIME
In this section we review, for pedagogical reasons, quantization in flat space-time of the free field theory with Lagrange density in eq. (2.12). The results presented here have, by in large, been noted previously in [24,25].
Dropping the tildes and setting a(t) = 1 the Lagrange density in eq. (2.12) becomes, This corresponds to normal kinetic terms for two real scalar fields but with an unusual Lorentz non-invariant kinetic mixing. The Lagrange density has the shift symmetry π → π + c for the Goldstone field π.
The classical equations of motion for the fields π and s are, Quantization proceeds by expanding the fields in modes, The annihilation operators a (1,2) (q) and creation operators a (1,2) (q) † satisfy the usual commutation relations 4 . The time dependence of the mode functions π (1,2) q (t) and s (1,2) q (t) are determined by solving the classical equations of motion and their normalization is fixed by the canonical commutation relations of the fields with their canonical momenta. A difference from the usual case where there is no Lorentz non-invariant mixing is that the canonical momentum for the field π is notπ but ratherπ + µs. Soπ andṡ don't commute at equal time but rather satisfy [π( The time dependence of the modes has the usual exponential form π . The dispersion relations for the energies is determined by solving the classical equations of motion. This yields, which is a massless mode that we label by (1) corresponding to the minus sign and a massive mode that we label by (2) corresponding to the plus sign. The mass of mode (2) is m 2 + µ 2 . Because this mode remains massive even when m = 0 there will be no divergences in our calculations in de-Sitter space. We now focus on the large mixing region of parameter space, µ q, m. As discussed in the literature [24] the dispersion relations of the two modes can be written as The (1) mode is massless but for q much larger than m the energy grows not linearly with q but rather quadratically (like a non relativistic particle). The other mode is massive with mass µ. For very small momentum, q m, the massive scalar s only contains the massive (2) mode, i.e., |s On the other hand the Goldstone field π contains equal amounts of the (1) and (2) modes.
The infrared, q → 0, behavior of the mode function s (1,2) q changes in the special case that m = 0. Then integrating-by-parts, the kinetic mixing term in Eq. (3.1) can be recast as −µπṡ, and so it is clear that there is also a shift symmetry in s. For m = 0 the scalar field s also contains equal amounts of the two modes.
Since for large µ the second mode is heavy it is appropriate for the physics at low momentum q µ to integrate it out from the theory and write an effective Lagrange density in terms of a single massless field. For the light mode a time derivative gives factors of 1/µ and (for m = 0) at very large µ the s field contains only a small amount of that massless mode. Hence eq.(3.3) imples that, Putting this into the Lagrange density in eq. (3.1) and dropping terms suppressed by powers of 1/µ (recall a time derivative on π is suppressed by 1/µ) yields the effective Lagrange density for the massless mode, which yields the dispersion relation for the massless mode given in eq. (3.7).
In the next section we perform the quantization in curved de-Sitter space-time (with Hubble constant H). Then the physics of the massless (1) mode should be similar to that in flat space-time when the momentum and energy for that mode are large compared to H i.e., q > H and E (1) q > H. In the flat space-time large µ discussion we assumed q < µ. The energy condition E (1) q > H implies q must also satisfy q > √ µH in order for our de-Sitter space-time computations to resemble the flat space-time large µ case discussed in this subsection.

IV. FREE FIELD THEORY IN DE-SITTER SPACE TIME
Introducing conformal time, τ = −e −Ht /H, and including the measure factor √ −g in the Lagrange density so that the action is equal to d 3 xdτ L we have As in flat space we expand the quantum fields in terms of creation and annihilation operators.
Introducing η = kτ we write, and The mode functions obey the classical equations of motion, and where a " " represents an η derivative.

A. Numerical results
In the mode expansion for for the fields s and π, k is the magnitude of the comoving wavevector. The physical wavevector has magnitude q = k/a = −Hη. Hence the condition that a mode have wavelength well within the de-Sitter horizon 1/H is q/H 1 which is equivalent to −η 1. At fixed k as time evolves a mode goes from physical wavelength well within the horizon to outside the horizon.
In the region well within the horizon, −η µ/H and −η 1, the differential equations (4.4) and (4.5) simplify to Here we suppressed the superscripts (1, 2) that label mode type. The leading behavior of the mode functions is and so it is convenient to represent the general solution in the region deeply inside the horizon as A and B are functions of η with |A /A|, |B /B| 1. Substituting π k and s k back into (4.4) and (4.5) and keeping only the leading order terms in η −1 we find Therefore, in this region the canonically normalized form of π k and s k can be written as where the factor H/ √ 4k 3 is determined by the canonical commutation relations. Eq. (4.11) is used to determine the initial conditions π at a value of η 0 that is large in magnitude. The differential equations in (4.4) and (4.5) can then be solved numerically and used to determine the power spectrum for the curvature perturbation in this model.
The correction to the power spectrum ∆P ζ is defined by, ∆P ζ = P ζ − P is the power spectrum of the curvature perturbation in usual slow roll single field inflation. ∆P ζ is shown in Fig. 1. In the region of µ H ∆P ζ goes like µ 2 which agrees with the perturbative calculation [15]. In the region where µ is larger than about 10H the power spectrum P ζ grows as µ 1/2 and can be approximated by, where C = 16π Γ 2 (−1/4) 2.09. (4.14) Corrections to eqs. (4.13) and (4.14) become negligible as µ → ∞. The power spectrum in the large µ limit was calculated using the large µ effective field theory in [25]. For completeness we briefly review that calculation in Sec. VII. As shown in Ref. [15], the perturbative result diverges in the limit of m → 0. From the red curve shown in Fig. 1 we can see that the curvature perturbation is well defined at m = 0. Perturbation theory can be very misleading at modest values of m and values of µ not very much larger than unity. For example for m = 0.5H and µ = 1.5H it gives a value for ∆P ζ (in the units used for Fig. 1) equal to 310 while our numerical result is 6.2.
For the curvature perturbations one calculates the power spectrum of the π field as −η → 0. However the power spectra for the fields can be calculated at any η. For µ/H > 1 the power spectrum for the s field P s (k) falls off rapidly as −η falls below unity. The numerical results of the power spectrum of the s field P s (k) in units of H 2 /2k 3 as a function of η for a few values of µ and m are shown in Fig. 2. One can see that all the curves decrease with −η and become small as −η falls below unity.
In the usual single field inflation model P π goes to unity in units of H 2 /2k 3 as −η → 0. In this model of quasi-single field inflation, as shown in Fig. 2 for the µ = 10H, m = 2H case the asymptotic value of P π is much larger than unity. This is due to the change in the dispersion relation of the π field and can be understood using the large µ effective theory. From Fig. 2 we see that the asymptotic value of P π for the case µ = 1.2H, m = 0.9H is also much larger than 1.

B. Qualitative analysis
We can understand qualitatively the shape of the mode functions analytically. In the region well outside the horizon, −η 1, eqs. (4.4) and (4.5) can be simplified to which is invariant under the transformation π k → λ 2 π k , s k → λ 2 s k , η → λη . Therefore, the general form of the solution can be written as Putting this back into the differential equations gives equations for the power α and the coefficients Q k and R k To have nontrivial solutions for Q k and R k requires There are four solutions to this equation For the region of parameter space we focus on, α ± are complex, which can have observational consequences for the non-gaussianities [19].
For large values of µ/H the infrared behavior of the mode functions π (1,2) k and s (1,2) k match directly onto the solutions in eq. (4.20). This is shown in Fig. 3 using m = 2H and µ = 10H. The α 1 = 0 mode is constant outside the horizon. The α 2 = 3 behavior vanishes outside the horizon and can be thought of as a subdominant contribution to the massless mode. The α ± solutions correspond to the mode functions for a free scalar field with mass equal to (m 2 + µ 2 ) 1/2 . They play an important role in the calculation of non-gaussianities. For m = 2H and µ = 10H the behavior of this mode is shown by the blue dot-dashed curves in Fig. 3. One can see that it oscillates logarithmically with frequency (m 2 + µ 2 ) 1/2 , and decreases with a power of 3/2 for small −η. To get the curves shown in Fig. 3 we solve the differential equations (4.4) and (4.5) with the initial conditions (4.11). The π (2) mode shown in the left panel of Fig. 3 eventually goes to a constant as −η gets smaller. Similarly, the absolute value of the s (1) mode eventually goes like (−η) 3/2 for very small −η.
In this paragraph we focus on the α 1 = 0 solution. Putting α 1 = 0 back Eq. (4.18) we find that R k = 0. Since there is no shift symmetry in the s field it should not contain the massless mode in the far infrared. We can get the leading behavior of the s k mode function outside the horizon by putting π k = Q k back into the exact differential equation (4.4). This gives the first order inhomogeneous differential equation with general solution This behavior is shown by the red dashed curves in Fig. 3.

C. The large µ/H region
In this subsection we focus on some properties of the solutions for the mode functions that only apply for very large µ/H. We find that the curvature perturbation goes to a constant when −η < (µ/H) 1/2 instead of the usual condition that it be outside the horizon, i.e., −η < 1. This is illustrated in Fig. 4 which shows the numerical results for the power spectrum of P π as a function of η.
Examining eq. (4.5), in the region −η < (µ/H) 1/2 it is clear that the last term on the left hand side is the largest. Neglecting the other terms the solution in this region satisfies which implies that π k is constant and s k is proportional to η 2 , as in eq. (4.22). In the region (µ/H) 1/2 < −η < µ/H one can show that the differential equations for the mode functions are solved approximately by The physical wavevector of a mode with comoving wavevector k is Therefore the change of the phase of these solutions within a small time period ∆η can be written as where ∆t = a∆τ has been used. This agrees with the dispersion relation in flat space given in eq. (3.7) for the massless mode. From Fig. 4, one can see that it is in this region the solution for µ H starts to deviate from the standard slow roll solution, which corresponds to µ = 0 in the model we are studying. This is because in this region the solutions in de-Sitter space should resemble those in flat space and the light mode has a flat space dispersion relation E q = q 2 /µ which is quite different from a single massless field with dispersion relation E q = q.
Putting the solution we have found back into the differential equations (4.4) and (4.5), one can see that the terms −π k + 2π k η and − s k + 2s k η (4.27) are suppressed, which means that the terms (∂ τ π) 2 and (∂ τ s) 2 (4. 28) in the Lagrange density (4.1) can be neglected. After neglecting these two terms, there are no terms in (4.1) that contain time derivatives of s. This indicates that s has become a Lagrange multiplier and can be replaced in the Lagrange density using its classical equation of motion to express it in terms of π. This amounts to summing the tree graphs that contain virtual s propogators and is the origin of the effective theory approach developed in Refs. [24] and [25] for the behavior of π in this region. We will briefly review the basic setup for this effective field theory and use it to calculate the two-and three-point functions of the curvature perturbation in the large µ limit in Sec. VII.

V. IMPACT ON OBSERVABLES
The dimensionless power spectrum is defined as [28] ∆ 2 ζ (k) = where f is a function of the µ and m. f − 1 is shown in Fig. 1 as a function of µ/H for fixed values of m. Throughout this section we neglect the impact of the time dependence ofφ 0 on the value of S 0 since, as was discussed in section II, it is suppressed by a power of Λ.
In terms of the slow-roll parameter can be written as The tilt of the power spectrum is defined as 4) and can be written as 5) where N is the number of e-folds between when the modes of interest exit the horizon and inflation ends. From Eq. (5.1) we have where the standard results of slow-roll inflation have been used [4], and η is the other slow-roll parameter defined as −φ 0 /(Hφ 0 ).μ andm are defined aŝ Up to leading order in the slow-roll parameters we have that, Therefore at leading order in slow roll parameters Another important observable is the tensor-scalar ratio. Since the gravitational wave production is only related to the structure of the de-Sitter metric, the dimensionless tensor spectrum can still be written as Then the tensor-scalar ratio can be written Here we use the model where the inflaton potential V φ = m 2 φ φ 2 /2 as an example to discuss the effect of large µ on the observables. In this simple model, we have 13) where N cmb is the number of e-folds between when CMB scale leaves the horizon and when slow roll inflation ends. The n S , r plot for this model is shown in Fig. 5. The dotted regions are for µ from 0 to 100H and m from 0 to 6H with (µ/H) 2 + (m/H) 2 > 9/4. On these curves as µ increases r decreases, so the uppermost point of the curves corresponds to standard slow roll inflation. The constraints on the m − µ parameter space for N CMB = 50 and 60 are also shown in Fig. 6 where the regions below the curves are excluded. Clearly larger values of µ improve the agreement of the model's predictions with the measured value of n S and the bound on r.  5. Impact on the scalar spectrum index n S and the tensor-to-scalar ratio r for the φ 2 inflation model with µ from 0 to 100H and m from 0 to 6H, and (µ 2 + m 2 ) 1/2 > 3H/2. The blue and red regions are for N cmb = 50 and 60 respectively. The dotted, dashed and solid curves are for m fixed to be 0, 3H/2 and 6H respectively. The gray regions are the one-sigma and two-sigma constraints from the combination of the Planck data and the BICEP2/Keck data [29].  [29], where the blue curves are for N CMB = 50 and the red curves N CMB = 60. The regions above the curves are allowed.

VI. NON-GAUSSIANITIES
In this section we calculate the dependence of the inflaton three-point function as a function of µ and m. The small µ behavior of the bispectrum was first studied in [15]. The effective field theory for large µ was used to compute the contribution from the ∂π∂πs interaction to the bispectrum [24,25]. Here we use the numerical mode functions to extend the analysis to other values of µ.
The curvature perturbation bispectrum B ζ (k 1 , k 2 , k 3 ) is defined by 1) and we can define B π (k 1 , k 2 , k 3 ) analogously. They can be computed using the in-in formalism [31] using the interaction Lagrangian in eq. (2.13).
In this section we focus mostly on the O(V S ) term (where V S ≡ V S (S 0 )) which, for V S ∼ O(H), typically dominates over the contribution from the ∂π∂πs term. We express the O(V S ) contribution to the bispectrum in terms of the mode functions discussed earlier.
Evaluating the correlator in the far future τ = 0, we find Equation (6.2) is true for all values of k i , however we are mostly interested in its behavior in the so-called equilateral and squeezed limits. In the equilateral limit, the external momenta all have equal magnitude k i ≡ k. In this case, the integral's dependence on k can be factored out of the integral by rescaling the integration variable from τ to η = kτ : 5 We can compute this integral numerically using the numeric mode functions, but there are a couple of subtleties in its evaluation that need to be addressed. The integrand in (6.3) is highly oscillatory at large τ . For m/H and µ/H values of order one or larger, the magnitude of these oscillations does not decay quickly and it becomes difficult to perform the numerical integrations by brute force. We can alleviate this problem by Wick rotating the integral, thereby transforming the rapid oscillations into exponential decay.
Before Wick rotating it is convenient to factor out the oscillatory behavior from the mode functions. The large τ limit given in eq. (4.8) suggests that we should extract the oscillatory behavior by factorizing the mode functions as π Plugging this factorization into B equil π (k) gives In the second line we used Cauchy's theorem to rotate the region of integration from the real to the imaginary axis and changed the integration variable from η to x = −iη.
The numerical solutions found previously for A and (4.5) (see [26]). After factoring out the oscillatory behavior and changing variables to x = −iη, we find that the analytically continued functions A (i) k and B (i) k obey where a prime denotes a derivative with respect to x and we have dropped the superscripts for simplicity. The solutions should asymptote at large −x to 6 A (1) These solutions and their derivatives with respect to x give the initial conditions for numerical integration of the differential equations for A k and B k . Note that A (k)/(P ζ (k)) 2 due to the ∂π∂πs and s 3 interaction terms respectively 7 . Moreover, we have superimposed a dotted line which corresponds to the prediction of the effective field theory appropriate for large µ (which will be discussed in detail in section VII). Of course, the numerical results converge to the effective field theory results in the large µ limit. However, the effective field theory is only a good approximation of these non-gaussianities for µ > ∼ 10H. This further suggests that there is a substantial portion of the parameter space in µ that is described neither by the large µ effective theory description nor the small µ perturbative description.
The Planck collaboration has derived constraints on the magnitude of the bispectrum of the curvature perturbations using various models/templates for its dependence on the wavevectors [30]. These are usually expressed in terms of the quantity f NL . Although the model we are discussing is different from the equilateral model/template used to derive the constraint f equil NL = 4 ± 43 by the Planck collaboration in Ref. [30], we use this constraint to estimate a bound on V S . Furthermore we estimate f equil NL using just the equilateral configuration where the three wavevectors have the same magnitude taking, To determine upper bounds for V S we assume that each interaction s 3 and s∂π∂π is separately constrained by f equil NL and thus ignore any possible tuning between the two terms that may make these bounds weaker. Figure 9 shows the 2σ upper bounds for a variety of s masses, as well as the upper bound predicted in the large µ effective theory. The squeezed limit of (6.2) occurs when k 1 ≈ k 2 ≡ k k 3 . In this limit, define the ratio c ≡ k 3 /k, where c 1, and introduce the notation B sq π (k, c) for B π . We again rescale the integration variable to η = kτ to find We can analyze the leading behavior of (6.9) in c by replacing s (6.11) We can compute β − , β + , and β 2 by fitting the numerical mode functions s Then, rearranging (6.10) gives We plot Im [λ + + λ − ] in figure 10. The sine term is usually smaller and so we have not displayed it in a figure. Equation (6.12) shows that the squeezed limit of the three-point function oscillates logarithmically as a function of c. This behavior is illustrated in figure  11. Note that that the dependence of Im [α + ] = m 2 /H 2 + µ 2 /H 2 − 9/4 on µ has an important effect on the oscillations. This impacts the two point function of biased objects, see for example [32].
The oscillatory terms in eq. (6.12) are enhanced by a factor of c −1/2 , but are suppressed in the large µ limit.

VII. CALCULATING NON-GAUSSIANITY IN THE EFFECTIVE THEORY
A. Brief review of the effective theory for large µ In this subsection we begin with a brief review the effective theory approach to the case when µ/H is large. In terms of π and s the Lagrange density is As discussed in Sec. III, in flat space-time with large mixing µ there is a very massive mode and a massless mode. When µ H and k/a < µ, one may integrate out the heavy mode to get an effective theory just involving π which can be used to calculate curvature perturbations. As discussed in Sec. IV, for that purpose the (∂ τ s) 2 and (∂ τ π) 2 terms in eq. (7.1) can be neglected. Since we assume m ∼ O(H) or smaller m can also be neglected in eq. (7.1). With these approximations the equation of motion for s becomes Up the second order in π, the solution for s is Putting this solution back into eq. (7.1), the quadratic and cubic terms of the effective Lagrangian of π can be written as Quantizing the free field part of this effective theory we write for the field operator, The mode function π k (η) satisfies the classical equation of motion, which can be solved analytically for the mode function π k (η). The normalization of π k (η) is determined by the canonical commutation relations. This yields, The power spectrum of the curvature perturbation is This result was originally derived in Ref. [24,25].
The plot of P ζ as a function of µ was shown in Fig. 1. The result from the effective theory is shown by the black dashed line. One can see that for µ > 10H the result from the effective theory agrees with the numerical result.

C. Non-gaussianity of squeezed configuration
For the squeezed configuration we consider k = k 1 k 2 k 3 = ck. Taking the contribution from the 1/Λ term in the interaction Lagrange density we have 1/4 (x) oscillate rapidly when x > 1. Therefore, the integral is mainly supported in the region x < 1, which means c 2 x 1. Around y = 0 we have For c 1, B 3 and B 4 go like c 2 and we have that in the squeezed limit B sq ζ ∼ c −1 . Even though this contribution is enhanced by a power of 1/c, it is still suppressed compared to what local non-gaussianity would give which is proportional to P ζ (k 1 )P ζ (k 2 )+P ζ (k 2 )P ζ (k 3 )+ P ζ (k 3 )P ζ (k 1 ) ∼ c −3 . This c −1 behavior in the squeezed limit is also seen in equilateral non-Gaussianity.
For the contribution proportional to V S we find 5/4 (y)) dy y→c 2 x (7.21) Therefore, the V S interaction also gives a c −1 contribution to B sq ζ .

VIII. CONCLUDING REMARKS
We studied a simple quasi-single field inflation model where the inflaton couples to another scalar field S. The model contains an unusual mixing term between the inflaton and the new scalar characterized by a dimensionful parameter µ. It has been extensively studied in the literature using perturbation theory in the region where the parameter µ/H is small and using an effective field theory approach in the region of large µ/H. It has also been studied using numerical methods in other regions of parameter space. When the mass parameter m of the additional scalar field is zero perturbation theory diverges.
We numerically calculated the power spectrum and the bispectrum of the curvature perturbations when µ and the mass m satisfy (µ/H) 2 + (m/H) 2 > 9/4 with m ∼ O(H) or smaller. In much of this region, perturbation theory and the effective field theory approach are not applicable. We found that typically the effective field theory approach is valid for µ/H > 10. The numerical approach is non-perturbative in µ/H and there are no divergences at m = 0. This occurs because the heavy mode has mass m 2 + µ 2 which does not vanish as m → 0.
In the case where the inflaton potential is m 2 φ φ 2 /2, we derived constraints on the parameters m and µ from n S and r for N cmb = 50 and N cmb = 60. Larger values of µ make this inflaton potential more compatible with the data.
We computed the contributions from the ∂π∂πs and the s 3 interactions to the equilateral limit of the bispectrum of the curvature perturbations numerically and compared it with the results from the effective theory. Using these results and the Planck bounds on f N L we derived upper bounds on V S and µ.
We also analyzed the squeezed limit of the bispectrum, showing that in this model it is much smaller than for local non-gaussianity. The contribution to the squeezed bispectrum proportional to V S exhibits interesting oscillatory behavior as a function of the ratio of the small momenta to the larger one. 9 We noted that the oscillation wavelength has µ dependence that is not evident in perturbation theory. This behavior could potentially be observed in future experiments.
For small µ and m, there are potentially interesting observational consequences of the behavior of the four point function on the wavevectors that characterize its shape. We will present results on this in a further publication.

ACKNOWLEDGMENTS
HA would like to thank Asimina Arvanitaki, Cliff Burgess and Yi Wang for useful comments and discussions. This work was supported by the DOE Grant DE-SC0011632. We are also grateful for the support provided by the Walter Burke Institute for Theoretical Physics.