Primordial Non-Gaussianities from Inflation Models

This is a pedagogical review on primordial non-Gaussianities from inflation models. We introduce formalisms and techniques that are used to compute such quantities. We review different mechanisms which can generate observable large non-Gaussianities during inflation, and distinctive signatures they leave on the non-Gaussian profiles. They are potentially powerful probes to the dynamics of inflation. We also provide a non-technical and qualitative summary of the main results and underlying physics.


Introduction
An ambitious goal of modern cosmology is to understand the origin of our Universe, all the way to its very beginning. To what extend this can be achieved largely depends on what type of observational data we are able to get. Thanks to many modern experiments, we are really making progress in this direction.
One of the representative experiments is the Wilkinson Microwave Anisotropy Probe (WMAP) satellite [1][2][3][4][5][6]. It measures the light coming from the last scattering surface about 13.7 billions years ago. This cosmic microwave background (CMB) is emitted at about 379,000 years after the Big Bang, when electrons and protons combine to form neutral hydrogen atoms and photons start to travel freely through the space. Our Universe was very young at that moment and the large scale fluctuations were still developing at linear level. So the CMB actually carries valuable information much earlier than itself, which can potentially tell us about the origin of the Big Bang.
There are two amazing facts about the CMB temperature map. On the one hand, it is extremely isotropic, despite the fact that the causally connected region at the time when CMB formed spans an angle of only about 0.8 degree on the sky today. On the other hand, we do observe small fluctuations, with ∆T /T ∼ 10 −5 .
The inflationary scenario [7][8][9] naturally solves the above two puzzles. It was proposed nearly thirty years ago to address some of the basic problems of the Big Bang cosmology, namely, why the universe is so homogeneous and isotropic. In this scenario, our universe was once dominated by dark energy and had gone through an accelerated expansion phase, during which a Hubble size patch was stretched by more than 60 efolds or so. Inhomogeneities and large curvature were stretched away by this inflationary epoch, making our current observable universe very homogenous and flat. In the mean while, the fields that were responsible for and participated in this inflationary phase did have quantum fluctuations. These fluctuations also got stretched and imprinted at superhorizon scales. Later they reentered the horizon and seeded the large scale structures today [10][11][12][13][14].
The inflationary scenario has several generic predictions on the properties of the density perturbations that seed the large scale structures: • They are primordial. Namely, they were laid down at superhorizon scales and entering the horizon after the Big Bang.
• They are approximately scale-invariant. This is because, during the 60 efolds, each mode experiences the similar expansion when they are stretched across the horizon.
• They are approximately Gaussian. In simplest slow-roll inflation models, the inflaton is freely propagating in the inflationary background at the leading order. This is also found to be true in most of the other models and for different inflationary mechanisms. So the tiny primordial fluctuations can be treated as nearly Gaussian.
The CMB temperature anisotropy is the ideal data that we can use to test these predictions. The obvious first step is to analyze their two-point correlation functions, i.e. the power spectrum. All the above predictions are verified to some extent [1]. The presence of the baryon acoustic oscillations proves that the density perturbations are indeed present at the superhorizon scales and reentering the horizon as the horizon expands after the Big Bang. The spectral index, n s = 0.963 ± 0.012, is very close to one, therefore, the density perturbations are nearly scale invariant. Several generic types of non-Gaussianities are constrained to be less than one thousandth of the leading Gaussian component.
But is this enough? Experimentally, the amplitude and the scale-dependence of the power spectrum consist of about 1000 numbers for WMAP. For the Planck satellite, this will be increased up to about 3000. However, we have about one million pixels in the WMAP temperature map alone, and 60 millions for Planck. So the information that we got so far is highly compressed comparing to what the data could offer in principle. This high compression is only justified if the density perturbations are Gaussian within the ultimate sensitivities of our experiments, so all the properties of the map is determined by the two-point function. Otherwise, we are expecting a lot more information in the non-Gaussian components.
Theoretically, inflation still remains as a paradigm. We do not know what kind of fields are responsible for the inflation. We do not know their Lagrangian. We also would like to distinguish inflation from other alternatives. Being our very first data on quantum gravity, we would like to extract the maximum number of information from the CMB map to understand aspects of the quantum gravity. All these motivate us to go beyond the power spectrum.
To give an analogy, in particle physics, two-point correlation functions of fields describe freely propagating particles in Minkowski spacetime. More interesting objects are their higher order correlations. Measuring these are the goals of particle colliders. Similarly, the power spectrum here describes the freely propagating particles in the inflationary background. To find out more about their interaction details and break the degeneracies among models, we need higher order correlation functions, namely non-Gaussianities. So the role non-Gaussianities play for the very early universe is similar to the role colliders play for particle physics.
With these motivations in mind, in this review, we explore various mechanisms that can generate potentially observable primordial non-Gaussianities, and are consistent with the current results of power spectrum. We will not take the approach of reviewing models one by one. Rather, we divide them into different categories, such that models in each category share the same physical aspect which leaves a unique fingerprint on primordial non-Gaussianities. On the one hand, if any such non-Gaussianity is observed, we know what we have learned concretely in terms of fundamental physics; on the other hand, explicit forms of non-Gaussianities resulted from this exploration provide important clues on how they should be searched in data. Even if the primordial density perturbations were perfectly Gaussian, to test it, we would still go through these analyses until various well-motivated non-Gaussian forms are properly constrained.

Road map
The following is the outline of the review. For readers who would like to get a quick and qualitative understanding of the main results instead of technical details, we also provide a shortcut after the outline.
In Sec. 2, we review the essential features of the inflation model and density perturbations. In Sec. 3, we review the first-principle in-in formalism and related techniques that will be used to calculate the correlation functions in time-dependent background.
In Sec. 4, we compute the scalar three-point function in the simplest slow-roll model. We list the essential assumptions that lead to the conclusion that the non-Gaussianities in this model is too small to be observed.
In Sec. 5, we review aspects of inflation model building, emphasizing various generic problems which suggest that the realistic model may not be the algebraically simplest. We also introduce some terminologies used to describe properties of non-Gaussianities.
Sec. 6, 7 and 8 contain the main results of this review. We study various mechanisms that can lead to large non-Gaussianities, and their distinctive predictions in terms of the non-Gaussian profile.
In Sec. 9, we give a qualitative summary of the main results in this review. Before conclusion, we discuss several useful relations among different non-Gaussianities.
Here is a road map for readers who wish to have a non-technical explanation and understanding of our main results. After reading the short review on the inflation model and density perturbations in Sec. 2, one may read the first and the last paragraph of Sec. 4 to get an idea of the no-go statement, and then directly proceed to read Sec. 5. After that, one may jump to Sec. 9 where the main results are summarized in non-technical terms.
The subject of the primordial non-Gaussianities is a fast-growing one. There exists many nice reviews and books in this and closely related subjects. The introductions to inflation and density perturbations can be found in many textbooks [15][16][17][18][19][20] and reviews [21][22][23][24][25]. Inflationary model buildings in particle physics, supergravity and string theory are reviewed in Ref. [26][27][28][29][30][31][32]. Comprehensive reviews on the developments of theories and observations of primordial non-Gaussianities up to mid 2004 can be found in Ref. [33,34]. Recent comprehensive reviews on theoretical and observational developments on the bispectrum detection in CMB and large scale structure has appeared in Ref. [35,36]. A recent comprehensive review on non-Gaussianities from the second order post inflationary evolution of CMB, which acts as contaminations of the primordial non-Gaussianities, has appeared in Ref. [37]. A recent review on how primordial non-Gaussianities can be generated in alternatives to inflation can be found in Ref. [38].

Inflation and density perturbations
In this section, we give a quick review on basic elements of inflation and density perturbations. We consider the simplest slow-roll inflation. The action is The first term is the Einstein-Hilbert action. The second term describes a canonical scalar field coupled to gravity through the metric g µν . This is the inflaton, which stays on the potential V (φ) and creates the vacuum energy that drives the inflation. We first study the zero-mode background evolution of the spacetime and the inflaton. The background metric is where a(t) is the scale factor and x is the comoving spatial coordinates. The equations of motion are 3)

4)
φ 0 + 3Hφ 0 + V ′ = 0 . (2.5) The first equation determines the Hubble parameter H, which is the expansion rate of the universe. The second equation is the continuity condition. The third equation describes the evolution of the inflaton. Only two of them are independent. The requirement of having at least O(60) e-fold of inflation imposes some important conditions. By definition, to have this amount of inflation, the Hubble parameter cannot change much within a Hubble time H −1 . This gives the first condition We also require that the parameter ǫ does not change much within a Hubble time, In principle, η can be close to O(1) but ǫ kept small. In such a case, ǫ grows exponentially with efolds and the inflation period tends to be shorter. More importantly, such a case will not generate a scale-invariant spectrum, as we will see shortly, thus cannot be responsible for the CMB. The two conditions (2.6) and (2.7) are called the slow-roll conditions. Using the background equations of motion, we can see that the slow-roll conditions impose restrictions on the rolling velocity of the inflaton. The first condition (2.6) implies thaṫ So the energy driving the inflation on the right hand side of (2.3) is dominated by the potential. Adding the second condition (2.7) further implies thaẗ So the first termφ 0 in (2.5) is negligible and the evolution of the zero-mode inflaton is determined by the attractor solution 3Hφ 0 + V ′ = 0 . (2.10) Using (2.10), the slow-roll conditions can also be written in a form that restricts the shape of the potential, They are related to ǫ and η by ǫ = ǫ V , η = −2η V + 4ǫ V . (2.12) So the shape of the potential has to be rather flat relative to its height. We emphasize that, although in this example several definitions of the slow-roll conditions are all equivalent, the definition (2.6) and (2.7) are more general. In other cases that we will encounter later in this review, these two conditions are still necessary to ensure a prolonged inflation and generate a scale-invariant spectrum, but the others no longer have to be satisfied. For example, the shape of potential can be steeper, or the inflationary energy can be dominated by the kinetic energy. Now let us study the perturbations. To keep things simple but main points illustrated, in this section, we will ignore the perturbations in the gravity sector and only perturb the inflaton, φ(x, t) = φ 0 (t) + δφ(x, t) . (2.13) We also ignore terms suppressed by the slow-roll parameters, which we often denote collectively as O(ǫ). For example, the mass of the inflaton is V ′′ ∼ O(ǫ)H 2 and will be ignored.
The quadratic Lagrangian for the perturbation theory is simply 14) and the equation of motion follows, where we have written it in the comoving momentum space, The solution to the differential equation (2.15), u(k, t), is called the mode function. It is not difficult to check that To quantize the perturbations according to the canonical commutation relations between δφ and its conjugate momentum δπ ≡ ∂L/∂δφ, with the commutation relations One can check that the commutation relations (2.18) and (2.21) are equivalent because of (2.17), given that the constant on the right hand side of (2.17) is specified to be i. This gives the normalization condition for the mode function. We now write down the mode function explicitly by solving (2.15), where we have used the conformal time τ defined as dt ≡ adτ . The infinite past corresponds to τ → −∞ and the infinite future τ → 0. We also used the relation τ = −1/Ha+O(ǫ). This mode function is a superposition of two linearly independent solutions with the normalization condition followed from (2.17). Consider the limit in which the mode is well within the horizon, i.e. its wavelength a/k much shorter than the Hubble length 1/H, and consider a time period much shorter than a Hubble time. In these limits, the mode effectively feels the Minkowski spacetime, and the first component in (2.22) with the positive frequency asymptotes to the vacuum mode of the Minkowski spacetime as we can see from (2.23). We choose this component as our vacuum choice, and it is usually called the Bunch-Davies state. The annihilation operator a p annihilates the corresponding Bunch-Davies vacuum, a p |0 = 0. The mode function has the following important properties. It is oscillatory within the horizon k|τ | ≫ 1. As it gets stretched out of the horizon k|τ | ≪ 1, the amplitude becomes a constant and frozen.
Physically this means that, if we look at different comoving patches of the universe that have the superhorizon size, and ignore the shorter wavelength fluctuations, they all evolve classically but with different δφ. This difference makes them arrive at φ f , the location of the end of inflation, at different times. This space-dependent time difference δt ≈ δφ/φ 0 leads to the space-dependent inflationary e-fold difference Again we ignore terms that are suppressed by the slow-roll parameters. This e-fold difference is the conserved quantity after the mode exits the horizon, and remains so until the mode re-enters the horizon sometime after the Big Bang. It is the physical quantity that we can measure, for example, by measuring the temperature anisotropy in the CMB, ζ ≈ −5∆T /T [39]. The information about the primordial inflation is then encoded in the statistical properties of this variable. So we would like to calculate the correlation functions of this quantity. Using (2.25), (2.19), (2.24) and (2.8), we get the following two-point function, where P ζ is defined as the power spectrum and in this case it is The spectrum index is defined to be where the relation d ln k = Hdt is used. If n s = 1, the spectrum is scale invariant. The current data from CMB tells us n s = 0.963 ± 0.012 [1]. So as we have mentioned, this requires a small η, which is also a value that tends to give more e-folds of inflation. If this were the end of story, all the primordial density perturbations would be determined by this two-point function and they are Gaussian. The rest of the review will be devoted to making the above procedure rigorous and to the calculations of higher order non-Gaussian correlation functions in this and various other models.

In-in formalism and correlation functions
In this section, we review the in-in formalism and the related techniques that are used to calculate the correlation functions in time-dependent background. The main procedure is summarized in the last subsection.
We are interested in the correlation functions of superhorizon primordial perturbations generated during inflation. So our goal is to calculate the expectation value of an operator Q, which is a product in terms of field perturbations δφ a and δπ a , at the end of inflation. The subscript a labels different fields. In inflation models, these fields are, for example, the fluctuations of the scalars and metric and their conjugate momenta. In the Heisenberg picture, where t is the end of inflation, |Ω is the vacuum state for this interacting theory at the far past t 0 . We start by looking at how the time-dependence in Q(t) is generated. The Hamiltonian of the system is a functional of the fields φ a (x, t) and their conjugate momenta π a (x, t) at a fixed time t. On the left hand side of (3.2) we have suppressed the variable x and index a which are integrated or summed over. The φ a (x, t) and π a (x, t) satisfy the canonical commutation relations and their evolution is generated by H following the equations of motion, We consider a time-dependent background,φ a (x, t) andπ a (x, t) which are c-numbers and commute with everything, and the perturbations, δφ a (x, t) and δπ a (x, t), (3.5) The background evolution is determined by the classical equations of motion, The commutation relations (3.3) become those for the perturbations, We expand the Hamiltonian as where we useH to denote terms of quadratic and higher orders in perturbations. Using (3.6), (3.7) and (3.8), the equations of motion (3.4) become (3.9) So the evolution of the perturbations, δφ a and δπ a , is generated byH. It is straightforward to verify that the solutions for (3.9) are with the condition at an initial time t 0 being To have a systematic scheme to do the perturbation theory, we splitH into two parts, The H 0 is the quadratic kinematic part, which in the perturbation theory will describe the leading evolution of fields. Fields whose evolution are generated by H 0 are called the interaction picture fields. We add a superscript "I" to label such fields. They satisfy The solutions are So the idea is to encode the leading kinematic evolution in terms of the interaction picture fields, and calculate the effects of the interaction through the series expansion in terms of powers of H I . To do this, we rewrite (3.1) as Using (3.11), (3.16) and (3.13), we get The solution to (3.20) and (3.21) can be written in the following way, where the operator T means that, in each term in the Taylor series expansion of the exponential, the time variables have to be time-ordered. The operatorT will be used to mean the reversed time-ordering. In summary, the expectation value (3.1) is Notice that in all the field perturbations are in the interaction picture.
The perturbation theory is also often done in terms of the Lagrangian formalism. In the following, we show that they are equivalent. In the above, we perform perturbations on the Hamiltonian, and define δπ a by perturbing π a ≡ ∂L/∂φ a , (here we use ∂ to denote the functional derivative,) The HamiltonianH is defined by (3.8). So using the definition together with the classical equations of motions (3.6) andπ a = ∂L/∂φ a , the definition (3.8) forH becomes If we perturb the Lagrangian directly, we keep the part of the Lagrangian that is quadratic and higher in perturbations δφ a and δφ a , The δπ a is defined directly as where in the second step Eq. (3.29) has been used. So these two definitions of δπ a are equivalent. The HamiltonianH is defined throughL, Again, using (3.29), we can see that the two definitions ofH are equivalent.

Mode functions and vacuum
The Hamiltonian H 0 in the above formalism is typically chosen to be the quadratic kinematic terms for field perturbations δφ a without mixing, (3.32) So they describe free fields propagating in the time-dependent background. The A, B and C are some time-dependent background fields, and they are all positive. The solutions to the equations of motion (3.14) in momentum space, u a (k, t), are called the mode functions, where k denotes the comoving momentum. They satisfy the Wronskian condition Note that we have specified the time-independent constant on the right hand side of (3.33) to be i for the same reason that we see in Sec. 2. Namely, we decompose δφ I a as where the annihilation and creation operators satisfy the following relations, These commutation relations are equivalent to (3.7) because of (3.33), but the constant needs to be i. This gives the normalization condition for the mode functions.
Being the solutions of the second order differential equation, generally the mode function is a linear superposition of two independent solutions. So we need to specify the initial condition. For inflation models, as long as the field theory applies, one can always find an early time at which the physical momentum of the mode is much larger than the Hubble parameter and study a time interval much less than a Hubble time. Under these conditions, the equations of motion approach to those in the Minkowski limit, in which the mode function is a linear superposition of two independent plane waves, one with positive frequency and another negative. The ground state in the Minkowski spacetime is the positive one. The mode function which approaches this positive frequency state in the Minkowski limit is called the Bunch-Davies state. In physical coordinates, this limit is proportional to e −ik ph t , (for k ph ≫ m), where k ph is the physical momentum. In terms of the conformal time τ ≡ dt/a(t) and the comoving momentum coordinate k ≡ k ph /a(t) which we often use, this limit is proportional to e −ikτ . We have seen an example in Sec. 2 and will see more similar examples later with different A, B and C. The corresponding vacuum |0 is the Bunch-Davies vacuum and annihilated by a a (k) defined in (3.34), a a (k)|0 = 0.
We also would like to write the vacuum of the interacting theory (3.23) in terms of the vacuum of the free theory |0 defined above. Unlike the scattering theory where the vacuum of the free theory is generally different from the vacuum of the interaction theory, the process that we are studying here do not generate any non-trivial vacuum fluctuations through interactions. This is a direct consequence of the identity (3.36) So we can replace |Ω in (3.23) with the Bunch-Davies vacuum |0 that we have specified above.
The integrand H I (t) in (3.22) is highly oscillatory in the infinite past due to the behavior of the mode function ∝ e −ikτ . Their contribution to the integral is averaged out. For the Bunch-Davies vacuum, this regulation can be achieved by introducing a small tilt to the integration contour τ 0 → −∞(1 + iǫ) or performing a Wick rotation τ → iτ . The imaginary component turns the oscillatory behavior into exponentially decay, making the integral well defined.

Contractions
When evaluating (3.23), one encounters (anti-)time-ordered integrals, of which the integrands are products of fields, such as δφ I a and δπ I a , or δφ I a and δφ I a , sandwiched between the vacua. In contrast to the Minkowski space, in the inflationary background, we do not have a simple analogous Feynman propagator which takes care of the time-ordering. Therefore we will just evaluate the integrands, but leave the complication of the time-ordering to the final integration.
To evaluate the integrand, one can shift around the orders of fields in that product, following the rules of the commutation relations. A contraction is defined to be a non-zero commutator between the following components of two fields, δφ + a , δφ − b , where δφ + a and δφ − b denote the first and second term on the right hand side of (3.34), respectively. After normal ordering, namely moving annihilation operators to the right-most and creation operators to the left-most so that they give zeros hitting the vacuum, it is not difficult to convince oneself that the only terms left are those with all fields contracted. Feynman diagrams can be used to keep track of what kind of contractions are necessary.
In the following we demonstrate this using a simple example. We consider a field δφ I and quantize it as usual, So a contraction between the two terms, δφ(k, t ′ ) on the left and δφ(p, t ′′ ) on the right, is defined to be For example, we want to compute a contribution to the four-point function δφ 4 from a tree-diagram containing two three-point interactions of the following form, These two H I 's come from expanding F −1 or F in (3.23). The corresponding Feynman diagram is Fig. 1. In Fig. 1, the two cubic vertices each represent the three-point interaction (3.39). Each line represents a contraction. The four outgoing legs connect to the four δφ(p i , t) (i = 1, 2, 3, 4) in δφ 4 . The following is a term from the perturbative series expansion of (3.23). We demonstrate in the following one set of contractions represented by the diagram in Fig. 1, Note that all terms are contracted. The result can be further evaluated using (3.38). After integration over momenta indicated in (3.39), the final momentum conservation will always manifest itself as (2π) 3 i (k i ). There are other sets of contractions represented by the same diagram for the same term. Namely, there are three ways of picking two of the three p i 's (q i 's), so we have a symmetry factor 9; also, there are 24 permutations of the four k i 's. We need to sum over all these possibilities. We also need to sum over all possible terms containing two H I 's in the perturbative series, which are not listed here, with their corresponding time-ordered integral structure.

Three forms
Now we deal with the time-ordered integrals in the series expansion. There are two ways to expand (3.23).
In the first form, we simply expand the exponential in (3.22). For example, for an even n, the nth order term is (3.40) Each term in the above summation contains two factors of multiple integrals, one from F −1 and another from F . Each multiple integral is time-ordered or anti-time-ordered, but there is no time-ordering between the two. We call this representation the factorized form.
In the second form, we rearrange the factorized form so that all the time variables are time-ordered, and all the integrands are under a unique integral. The nth order term in this form is [45] i n t t 0 (3.41) We call this representation the commutator form. Each representation has its computational advantages and disadvantages. The factorized form is most convenient to achieve the UV (t i → t 0 ) convergence. As mentioned, after we tilt or rotate the integration contour into the positive imaginary plane for the left integral, and negative imaginary plane for the right integral, all the oscillatory behavior in the UV becomes well-behaved exponential decay. However this form is not always convenient to deal with the IR (t i → t) behavior. For cases where the correlation functions have some non-trivial evolution after modes exit the horizon, as typically happens for inflation models with multiple fields, the convergence in the IR is slow. Cancellation of spurious leading contributions from different terms in the sum (3.40) can be very implicit in this representation, and could easily lead to wrong leading order results in analytical estimation or numerical evaluation.
The commutator form is most convenient to get the correct leading order behavior in the IR. The mutual cancellation between different terms are made explicit in terms of the nested commutators, before the multiple integral is performed. However, such a regrouping of integrands makes the UV convergence very implicit. Recall that the contour deformation is made to damp the oscillatory behavior in the infinite past. In the commutator form, for any individual term in the integrand, we can still generically choose a unique convergence direction in terms of contour deformation. Although the directions are different for different terms, they can be separately chosen for each of them. But now the problem is, if these different terms have to be grouped as in the nested commutator so that the IR cancellation is explicit, the two directions get mixed. Hence the explicit IR cancellation is incompatible with the explicit UV convergence in this case.
To take advantage of both forms, we introduce a cutoff t c and write the IR part of the in-in formalism in terms of the commutator form and the UV part in terms of the factorized form [46], This representation is called the mixed form. This form is particularly efficient in numerical computations when combined with the Wick-rotations in the UV.
We will not always encounter all these subtleties in every model, but there does exist such interesting examples, as we will see in Sec. 7.1.

Summary
To end this section, we summarize the procedure that we need to go through to calculate the correlation functions in the in-in formalism.
Starting with the Lagrangian L[φ(t),φ(t)], we perturb it around the homogenous solutions φ a andφ a of the classical equations of motion, φ a (x, t) =φ a (t) + δφ a (x, t) ,φ a (x, t) =φ a (t) + δφ a (x, t) . (3.43) Keep the part of the Lagrangian that is quadratic and higher in perturbations and denote it asL. Define the conjugate momentum densities as δπ a = ∂L/∂(δφ a ). We can also equivalently expand the Hamiltonian H[φ(t), π(t)] by perturbing φ a (x, t) and π a (x, t).
Work out the Hamiltonian in terms of δφ a and δπ a , and separate them into the quadratic kinematic part H 0 , which describes the free fields in the time-dependent background, and the interaction part H I . Relabel δφ a 's and δπ a 's in the Hamiltonian density as the interaction picture fields, δφ I a 's and δπ I a 's. These latter variables satisfy the equations of motion followed from the H 0 . We quantize δφ I a and δπ I a in terms of the annihilation and creation operators as in (3.34) and (3.35). The mode functions u a (k, t) are solutions of the equations of motion from H 0 , normalized according to the Wronskian conditions (3.33) and specified by an initial condition such as the Bunch-Davies vacuum. The correlation function for Q(t) is given by where Q(t) is a product in terms of δφ I a (x, t) and δπ I (x, t). If we want to work with δφ I a and δφ I a instead of δφ I a and δπ I a , we replace δπ I a withδφ I a using the relationδφ I a = ∂H 0 /∂(δπ I a ). Choose appropriate forms in Sec. 3.4 and series-expand the integrand in powers of H I to the desired orders. Perform contractions defined in Sec. 3.3 for each term in this expansion. Each term gives a non-zero contribution only when all fields are contracted. Draw Feynman diagrams that represent the correlation functions, and use them as a guidance to do contractions. Finally sum over all possible contractions and perform the time-ordered integrations.

A no-go theorem
Simplest inflation models generate negligible amount of non-Gaussianities that are well below our current experimental abilities [47,48]. By simplest, we mean This list is extracted based on Maldacena's computation of three-point functions in an explicit slow-roll model [47]. We now review this proof. The notations here follow Ref. [49,50] and will be consistently used later in this review.
The Lagrangian for the single scalar field inflation with canonical kinetic term is the following, where φ is the inflaton field, X = − 1 2 g µν ∂ µ φ∂ ν φ is the canonical kinetic term and V is the slow-roll potential. The first term is the Einstein gravity and M P = (8πG) −1/2 is the reduced Planck mass. For convenience we will set the reduced Planck mass M P = 1. The signature of the metric is (−1, 1, 1, 1).
The inflaton starts near the top of the potential and slowly rolls down. As we have reviewed in Sec. 2, to ensure that the inflation lasts for at least O(60) efolds, the potential is required to be flat so that the slow-roll parameters (2.11) are both much less than one most of the time. The energy of the universe is dominated by the potential energy, and the inflaton follows the slow-roll attractor solution (2.10). Also as discussed in Sec. 2, we will use the following more general slow-roll parameters, To study the perturbation theory, it is convenient to use the ADM formalism, in which the metric takes the form The action becomes In the ADM formalism, the variables N and N i are Lagrangian multipliers whose equations of motion are easy to solve. In single field inflation, we have only one physical scalar perturbation [51]. We choose the uniform inflaton gauge (also called the comoving gauge) in which the scalar perturbation ζ appears in the three dimensional metric h ij in the following form, and the inflaton fluctuation δφ vanishes. The a(t) is the homogeneous scale factor of the universe, so ζ is a space-dependent rescaling factor. In this review we do not consider the tensor perturbations. We plug (4.3) and (4.6) into the action (4.4) and solve the constraint equations for the Lagrangian multipliers N and N i . We then plug them back to the action and expand up to the cubic order in ζ in order to calculate the three-point functions. To do this, in the ADM formalism, it is enough to solve N and N i to the first order in ζ. This is because their third order perturbations will multiply the zeroth order constraint equation which vanishes, and their second order perturbations will multiply the first order constraint equation which again vanishes. After some lengthy algebra, we obtain the following expansions, Here ∂ −2 is the inverse Laplacian and δL/δζ| 1 is the variation of the quadratic action with respect to the perturbation ζ. We now can follow Sec. 3 and proceed to calculate the correlation functions. For simplicity, we will always neglect the superscript "I" on various interaction picture fields. We restrict to the case where the slow-roll parameters are always small and featureless. We first look at the quadratic action. In this case, we can analytically solve the equation of motion followed from (4.7) in terms of the Fourier mode of ζ, and get the mode function where τ ≡ dt/a ≈ −(aH) −1 is the conformal time. The normalization is determined by the Wronskian condition (3.33). We have chosen the Bunch-Davies vacuum by imposing the condition that the mode function approaches the vacuum state of the Minkowski spacetime in the short wavelength limit k/a ≫ 1/H, (4.14) The dynamical behavior of ζ that has been explained around Eq. (2.24) and (2.25) is made precise here. In particular, ζ is exactly massless without dropping any O(ǫ) suppressed terms. In addition, from (4.6), we can see that, for superhorizon modes, the only effect of ζ is to provide a homogeneous spatial rescaling. And ζ is the only scalar perturbation. So the fact that ζ is frozen after horizon exit will not be changed by higher order terms.
If we choose the spatially flat gauge, we make ζ disappear and the scalar in this perturbation theory becomes the perturbation of φ. The relation between ζ and δφ in (2.25) (with O(ǫ) corrections) is thus a gauge transformation through a space-dependent time-shift.
We quantize the field as We can easily compute the two-point function at the tree level, where Since ζ remains constant after it exits the horizon, the H and ǫ are both evaluated near the horizon exit. We next look at the cubic action. For single field models, H I,3 = −L 3 . Keeping in mind that χ is proportional to ǫ, one can see that the first line of (4.8) is proportional to ǫ 2 . For the featureless potential,η = O(ǫ 2 ), where ǫ collectively denotes either ǫ or η. So the second line of (4.8) is proportional to ǫ 3 , and negligible. The third line can be absorbed by a field redefinition ζ → ζ n + f (ζ n ). The only term in f (ζ n ) that will contribute to the correlation function is written out explicitly in (4.11). All the others involve derivatives of ζ so vanish outside the horizon. Thus this redefinition eventually introduces an extra term According to (3.44), we expand the exponential to the first order in H I,3 to get the leading result, To estimate the order of magnitude of the bispectrum, we only need to keep track of the factors of H and ǫ. For example, from the first term in (4.8), we have dtH 3 (t) ⊃ − dx 3 dτ a 2 ǫ 2 ζζ ′2 , where we used the conformal time τ and the prime denotes the derivative to τ . Using a ∝ H −1 , ζ ∝ ζ ′ ∝ H/ √ ǫ, we see that this three-point vertex contributes Together with the three external legs ζ 3 and the definition P ζ ∝ H 2 /ǫ, we get Similar results can be obtained for the other two terms in the first line of (4.8). As we will define more carefully later, the size of the three-point function is conventionally characterized by the number f N L , which is defined as ζ 3 ∼ f N L P 2 ζ . So the contribution from the first line of (4.8) is f N L = O(ǫ). The extra term due to the redefinition (4.18) This completes the order-of-magnitude estimate. To get the full non-Gaussian profile, we need to compute the integrals explicitly and get where K = k 1 + k 2 + k 3 and the permutations stand for those among k 1 , k 2 and k 3 . The slow-roll parameters are of order O(0.01), so f N L ∼ O(0.01) for these models. Even if we start with Gaussian primordial perturbations, non-linear effects in CMB evolution will generate f N L ∼ O(1) [37], and a similar number for large scale structures due to the non-linear gravitational evolution or the galaxy bias [35]. It seems unlikely that we can disentangle all these contaminations and detect such small primordial non-Gaussianities in the near future.
5 Beyond the no-go

Inflation model building
The following are two examples of slow-roll potentials in the simplest inflation models that we studied in Sec. 4, The first type (5.1) belongs to the small field inflation models. The slow-roll conditions (2.11) require the potential to be flat enough relative to its height, i.e. the mass of the inflaton should satisfy m ≪ H. The second type (5.2) belongs to the large field inflation models. The potential also needs to be flat relative to its height, but here one achieves this by making the field range φ very large, typically φ ≫ M P . The other conditions that we listed in Sec. 4 should also be satisfied by these models. These are the classic examples, which exhibit algebraic simplicities and illustrate many essential features of inflation. However, when it comes to the more realistic model building in a UV complete setup, such as in supergravity and string theory, situations get much more complicated. For example, it is natural that we encounter multiple light and heavy fields, and the potentials for them form a complex landscape. These multiple fields live in an internal space, whose structure can be very sophisticated. In string theory, some of them manifest themselves as extra dimensions and can have intricate geometry and warping. All these elements have to coexist with the inflationary background that introduces profound back-reactions.
Even with varieties of model building ingredients, it has been proven to be very subtle to construct an explicit and self-consistent inflation model. Indeed various problems have been noticed over the years in the course of the inflation model building. For example, • The η-problem for slow-roll inflation [52]. As we have seen, in order to have slow-roll inflation [8,9], the mass of the inflaton field has to be light enough, m ≪ H, to maintain a flat potential. However, in the inflationary background, the natural mass of a light particle is of order H. This can be seen in many ways, and in some ideal situations they are equivalent to each other. For example, one way to see this is to consider the coupling between the Ricci scalar and the inflaton, ∼ Rφ 2 . In the inflationary background R ∼ H 2 . Unless we have good reasons to set the coefficient of this term to be much less than one, it will give inflaton a mass of order H, spoiling the inflation. Another way to see this is to note that the effective potential in supergravity takes the form V = V 0 exp(K/M 2 P ) × other terms. Here schematically K ∼ φ 2 + · · · is the Kahler potential and its dependence on φ is normalized as such to give the canonical kinetic term for φ. So the first term in the expansion of V is of order V 0 φ 2 /M 2 P ∼ H 2 φ 2 and model independent. Therefore, either symmetry needs to be imposed or other tuning contributions introduced to solve this η-problem.
• The h-problem for DBI inflation [53]. DBI inflation [54] is invented to generate inflation by a different mechanism. It makes use of the warped space in the internal field space [55,56]. These warped space impose speed-limits for the scalar field, so even if the potential is steep, the inflaton is not allowed to roll down the potential very quickly. A canonical example of warped space is where r is the extra-dimension (or internal space), h(r) = r/R is the warp factor and R is the length scale of the warped space. The position of a 3+1 dimensional brane in the r-coordinate is the inflaton. So the inflaton velocity is limited by the speed-limit in the r-direction, h 2 . In order to provide a speed-limit that is small enough for inflation, the warp factor has to be small enough, h ≪ HR. However one of the Einstein equations with the metric (5.3) takes the following form, where the second term on the left hand side is due to the back-reaction of the inflationary spacetime. It is easy to see that the naive h = r/R should be modified for h < HR, precisely where the inflation is supposed to happen. Without contributions from other source terms, such a deformed geometry closes up too quickly and leads to an unacceptable probe-brane back-reaction if we demand the inflaton still follow the speed-limit. Therefore, either symmetry, or tuning using other source terms from the right hand side of (5.4), is necessary to solve this h-problem. The η-problem and h-problem are closely related in an AdS/CFT setup.
• The field range bound [57,58]. Large field inflation models require the field range to be much larger than M P . In supergravity and string theory, starting from a tendimensional theory with 10-dim Planck mass M 10 , the 4-dim Planck mass M P is obtained by integrating out the six extra-dimensions, where we use L and V (6) ∼ L 6 to denote the size and volume of the extra-dimensions, respectively. The field range ∆φ often appears as the distance in the extra dimensions, ∆φ ∼ ∆L · M 2 (10) , with the factor M 2 (10) being the proportional coefficient. Clearly, ∆L L. If the field range manifests itself within a warped throat with a length scale R, we still require R < L, and so ∆L L. Together with M P = M 4 (10) L 3 , we get We further note that the microscopic length scale L has to be much larger than the 10-dim Planck length M −1 (10) for the field theory to make sense. So M (10) L ≫ 1, and the field range ∆φ in these models is generically sub-Planckian. For example, for a warped throat with charge N, (M (10) L) 2 We have ignored a detailed numerical coefficient appearing on the right hand side of (5.7), which is model dependent. For example, considering the volume V (6) to be the sum of the throat and a generic bulk volume, it is O(0.01) [57]; considering an extreme case where the throat does not attach to a bulk, it is O(1) [58]. Notice that, due to the dependence of M P on the volume V (6) , increasing the volume only makes the bound tighter.
• The variation of potential [59]. Even in cases where there is no fundamental restriction on the excursion of fields, one encounters problems constructing the large field inflationary potential. Large field potentials that arise from a fundamental theory take the following general from, where m fund represents typical scales in the theory. For field theory descriptions to hold, such scales are much less than M P . For example, m fund can be the higher dimensional Planck mass, string mass, or their warped scales. The λ n 's are dimensionless couplings of order O(1). Unless some symmetries are present to forbid an infinite number of terms in (5.8), or a high degree of fine-tuning is assumed, the shape of potential (5.8) varies over a scale of order m fund ≪ M P . This variation is too dramatic for the potential to be a successful large field slow-roll potential.
None of the arguments in the above list is meant to show that the specific type of inflation is impossible. In fact, these have been the driving forces for the ingenuity and creativity in the field of inflation model building. This list is used to demonstrate some typical examples of complexities in reality. Often times, solving one problem will be companied by other structures that make the model step beyond the simplest one. So we may want to keep an open mind that the algebraic simplicity may not mean the simplicity in Nature.
Following is a partial list of possibilities that allow us to go beyond the no-go theorem in Sec. 4.
• Instead of single field inflation, we can consider quasi-single field or multifield inflation models (Sec. 7 & 8).
• Instead of canonical kinetic terms, there are models where the higher derivative kinetic terms dominate the dynamics (Sec. 6.1).
• Instead of following the attractor solution such as the slow-roll precisely, features can be present in the potentials or internal space, that temporarily break the attractor solution, or cause small but persistent perturbations on the background evolution (Sec. 6.2 & 6.3).
• Instead of staying in the Bunch-Davies vacuum, other excitations can exist due to, for example, boundary conditions or low scales of new physics (Sec. 6.4).
• Although strong constraints, from experimental results and theoretical consistencies, exist on non-Einstein gravities, early universe may provide an opportunity for their appearance. We use this category to include a variety of possibilities, such as modified gravities, non-commutativity, non-locality and models beyond field theories.
There are also strong motivations from data analyses for us to search and study different large non-Gaussianities. The signal-to-noise ratio in the CMB data is not large enough for us to detect primordial non-Gaussianities model-independently. A well established method is to start with a theoretical non-Gaussian ansatz, and construct optimal estimators that compare theory and data by taking into accounts all momenta configurations. This then gives constraints on the parameters characterizing the theoretical ansatz. Therefore, the following two important possibilities exist. First, the primordial non-Gaussianities exist in data could be missed if we did not start with a right theoretical ansatz. Second, even if a non-Gaussian signal were detected with one ansatz, it does not mean that we have found the right one. So different well-motivated non-Gaussian templates are needed for clues on how corresponding data analyses should be formed. From a different perspective, even if the primordial density perturbations were Gaussian, we would still do the similar amount of work and reach the conclusion after various well-motivated non-Gaussian forms are properly constrained.

Shape and running of bispectra
In this review, we will be mainly interested in the three-point correlation functions of the scalar primordial perturbation ζ. They are also called the bispectra. In this subsection, we introduce some simple terminologies that we often encounter in studies of bispectra.
The three-point function is a function of three momenta, k 1 , k 2 and k 3 , which form a triangle due to the translational invariance. Assuming also the rotational invariance, we are left with three variables, which are their amplitudes, k 1 , k 2 and k 3 , satisfying the triangle inequalities. The information is encoded in a function S(k 1 , k 2 , k 3 ) that we define as below, whereP ζ is the fiducial power spectrum, and we fix it to be a constantP ζ ≡ P ζ (k wmap ) = 6.1 × 10 −9 , where k wmap = 0.027Mpc −1 . We have chosen the above definition so that it can be uniformly applied to different types of bispectra that we will encounter in this review. In literature, different notations have been used. The differences are simple and harmless. For example, different functions such as A = k 1 k 2 k 3 S or F = S/(k 1 k 2 k 3 ) 2 are sometimes defined. We choose S since it is dimensionless and, for scale-invariant bispectra, it is invariant under a rescaling of all momenta. This quantity is the combination that is used to plot the profiles of bispectra in literature any way, despite of different conventions. Also, the precise power spectrum P ζ instead ofP ζ is often used in the definition (5.9). Here we absorb the momentum dependence of P ζ in S. This is because the three-point function is an independent statistic relative to the two-point. In cases where both the power spectrum and bispectrum have strong scale dependence, it is not convenient if they are defined in an entangled way. Under different circumstances, different properties of S are emphasized. The conventions involved may not always be precisely consistent with each other, since they are chosen to best describe the case at hand. Following are some typical examples.
The dependence of S on k 1 , k 2 and k 3 is usually split into two kinds.
One is called the shape of the bispectrum. This refers to the dependence of S on the momenta ratio k 2 /k 1 and k 3 /k 1 , while fixing the overall momentum scale K = k 1 + k 2 + k 3 . Several special momentum configurations are shown in Fig. 2.
Another is called the running of the bispectrum. This refers to the dependence of S on the overall momentum scale K = k 1 + k 2 + k 3 , while fixing the ratio k 2 /k 1 and k 3 /k 1 .
For bispectra that are approximately scale-variant, the shape is a more important property [60,50]. We will encounter such cases in Sec. 6.1, 7.1 & 8.1. The amplitude, also called the size, of the bispectra is often denoted as f N L by matching In this case, f N L is approximately a constant but can also have a mild running, i.e. a weak dependence on the overall momentum K [61,62]. An index n N G − 1 ≡ d ln f N L /d ln k is introduced to describe this scale dependence. The power spectrum also has a mild running, In this review, when we give explicit forms of S in the approximately scale-invariant cases, for simplicity we mostly ignore these mild scale dependence and concentrate on shapes. Shapes of bispectra have been given names according to the overall dependence of S on momenta. For example, for the equilateral bispectrum, S peaks at the equilateral triangle limit and vanishes as ∼ k 3 /k 1 in the squeezed triangle limit (k 3 ≪ k 1 = k 2 ). The local bispectrum peaks at the squeezed triangle limit in the form ∼ (k 3 /k 1 ) −1 , such as the two shape components in (4.22). To visualize the shapes, we often draw 3D plots S(1, x 2 , x 3 ), where x 2 and x 3 vary from 0 to 1 and satisfy the triangle inequality There are also cases where the running becomes the most important property, while the shape is relatively less important [64,65]. In such cases, the bispectra are mostly functions of K. So f N L defined in (5.10) has strong scale dependence. Instead, one can choose a constant f N L to describe the overall running amplitude. We will encounter such cases in Sec. 6.2 & 6.3. In these cases, the shape plot S(1, x 2 , x 3 ) may look nontrivial but this is because it does not fix K.
The above dissection will become less clean for cases where both properties become important.
One thing is clear. The f N L , that is always used to quantify the level of non-Gaussianities, is only sensible with an extra label that specifies, at least qualitatively, the profile of the momentum dependence, such as shapes and runnings.
It is useful to quantify the correlations between different non-Gaussian profiles, because as we mentioned in data analyses an ansatz can pick up signals that are not completely orthogonal to it. In real data analyses this is performed in the CMB l-space. To have a simple but qualitative analogue in the k-space, we define the inner product of the two shapes as and normalize it to get the shape correlator [60,63] Following Ref. [63], we choose the weight function to be so that the k-scaling is close to the l-scaling in the data analyses estimator. Later in this review, when we use this correlator to estimate the correlations between shapes, we take the ratio between the smallest and largest k to be 2/800, close to that in WMAP. A more precise correlator should be computed in the l-space in the same way that the estimator is constructed. We refer to Ref. [35] for more details.
In typical data analyses [66][67][68][69][70], the estimator involves a triple integral of the bispectrum over the three momenta k i . To have practical computational costs, it is necessary that this integral can be factorized into a multiplication of three integrals, each involves only an individual k i . This requires the bispectrum to be of the form , or a sum of such forms. Such a form is called the factorizable form or separable form. The factor K −n may be tolerable since it can be written as (1/Γ(n)) ∞ 0 t n−1 e −Kt . If the analytical result is too complicated, to make contact with experiments we will try to construct simple factorizable ansatz or template to capture the main features of the original one. New methods that are applicable to non-factorizable bispectrum forms and are more model-independent are under active development [71].

Single field inflation
In this section, we relax several restrictions of the no-go theorem on single field inflation models and study how large non-Gaussianities can arise. We present the formalisms and compute the three-point functions. We emphasize how different physical processes during inflation are imprinted as distinctive signatures in non-Gaussianities. Obviously, any mechanism that works for single field inflation can be generalized to multi-field inflation models.

Equilateral shape: higher derivative kinetic terms
In this subsection, we study large non-Gaussianities generated by non-canonical kinetic terms in general single field inflation models, following Ref. [50].
Consider the following action for the general single field inflation [72], Comparing to (4.1), we have replace the canonical form X − V with an arbitrary function of X ≡ − 1 2 g µν ∂ µ φ∂ ν φ and φ. This is the most general Lorentz-invariant Lagrangian as a function of φ and its first derivative. It is useful to define several quantities that characterize the differential properties of P with respect to X [72,49], where c s is called the sound speed and the subscript "X" denotes the derivative with respect to X. The third derivative is enough since we will only study the three-point function here. It is a non-trivial question which forms of P will give rise to inflation. The modelindependent approach we take here is to list the conditions that an inflation model has to satisfy, no matter which mechanism is responsible for it. Namely, we generalize the slow-roll parameters in (4.2) to the following slow-variation parameters and require them to be small most of the time during the inflation. The smallness of these parameters guarantees the Hubble constant H, the parameter ǫ and the sound speed c s to vary slowly in terms of the Hubble time. Similar to the arguments given in the case of slow-roll inflation in Sec. 2, these are necessary to ensure a prolonged inflation as well as an approximately scale-invariant power spectrum that we observed in the CMB. Following the same procedure that is outlined in Sec. 4, we get the quadratic and cubic action for the scalar perturbation ζ [47,49,50]. The quadratic part is If the slow-variation parameters are always small and featureless, we can analytically solve the equation of motion followed from (6.6) and get the following mode function, Notice the appearance of c s comparing to (4.13). The two-point function is with the power spectrum where the variables are evaluated at the horizon crossing of the corresponding k-mode.
To calculate the bispectrum, we look at the cubic action. In the following, we list three terms that are most interesting for this subsection, The full terms can be found in Eq. (4.26)-(4.28) in Ref. [50]. The order of magnitude contribution from these three terms can be estimated similarly as we did in (4.20), but now we not only keep factors of H and ǫ, but also factors of c s . Take the first term as an example, we write it in terms of the conformal time, Multiplying the three external legs ζ 3 , and using the definition and 14) The other two terms are similar. A detailed calculation reveals where we have decomposed the shape of the three-point function into six parts. The first two come from the leading order terms that we listed in (6.10), In terms of f N L their sizes are The next four terms come from the subleading terms that we did not list explicitly in (6.10) as well as the subleading contributions from the first two terms. Their orders of magnitude are The detailed profiles can be found in Ref. [50].  The full results we obtained can be used in different regimes.
• If we look at the limit, c s ≪ 1 or λ/Σ ≫ 1, the leading order results give two shape components, S λ and S c . This result can also be obtained using a simple method of considering only the fluctuations in scalar field while neglecting those in gravity [74,75]. Intuitively, this is because the higher derivative terms are responsible for the generation of large non-Gaussianities, and the gravity contribution is expected to be small as we saw in Sec. 4. Therefore one expands P (X, φ) using φ(x, t) = φ 0 (t) + δφ(x, t). The derivatives of P with respect to φ are ignored because the inflation and scale invariance imposes an approximate shift symmetry on P in terms of inflaton φ. We then get two terms in the cubic Lagrangian density This gives two leading bispectra the same as (6.16) and (6.17). The approach that we present here gives a rigorous justification to such an method. The subleading order component S o may be observable as well. At this limit where the higher derivative terms of the inflaton field are dominant, the Lagrangian of the above effective field theory are generalized [76] to include, for example, the ghost inflation [77] whose Lagrangian cannot be written in a form of P (X, φ). Another two slightly different equilateral shapes arise. However it is worth to mention that, generally in single field models and Einstein gravity, going beyond P (X, φ) requires adding either terms that explicitly break the Lorentz symmetry, or terms with higher time derivatives on φ which cannot be eliminated by partial integration, such as ( φ) 2 . Different treatment of such terms and discussions on their effects can be found in [78][79][80].
• If we take the opposite, slow-roll limit, c s → 1 and λ/Σ → 0, we recover the two shape components S ǫ and S η that we got in Sec. 4, with unobservable size f N L ∼ O(ǫ).
• We can also look at the intermediate parameter space. In slow-roll inflation models, one can also add higher derivative terms [73,49]. But in order not to spoil the slow-roll mechanism, the effect of these terms can only be subdominant. This corresponds to c s ≈ 1 and λ/Σ < O(1). Using the full results, we can see that the size of the non-Gaussianity is f N L < O(1). Therefore it is important to emphasize that, for the class of models we consider here, non-slow-roll inflationary mechanisms, such as the example that will be given below, are necessary to generate observable large non-Gaussianities.
• The other terms that we did not list in (6.10) (see Ref. [50]) and their canonical limit (4.8) are also useful. These terms are exact for arbitrary values of ǫ, η and s, so the usage of the action is beyond the category of models that we focus on in this subsection. As we will see in Sec. 6.2 and Sec. 6.3, it can be applied to the cases of sharp or periodic features where these parameters do not always remain small. In the rest of this subsection, we focus on the first case. In Fig. 3 and 4, we draw the shapes of S λ and S c . The two shapes are similar. They both peak at the equilateral limit, and behave as S ∼ k 3 /k 1 in the squeezed limit k 3 ≪ k 1 = k 2 . We call these shapes the equilateral shapes. There are some small differences between S λ and S c , e.g. around the folded triangle limit k 2 + k 3 = k 1 . A factorizable shape ansatz for the equilateral shape that is often used in data analyses is the following [81]: 22) and is shown in Fig. 5. As we can see, it represents the most important features of Fig. 3 & 4.
The shape of S o is more complicated, but we expect they have the similar shapes as the equilateral one because their squeezed limits behave the same [50]. The three other shapes S ǫ , S η and S s are all close to the local shapes as their squeezed limit scale as k 1 /k 3 for The scale dependence in P ζ , c s and λ/Σ will introduce mild running for the three-point function. We usually regard only the contributions from c s and λ/Σ as the running of the non-Gaussianity.
The underlying physics of the equilateral shape can be readily understood in terms of their generation mechanism. In single field inflation, the long wavelength mode that exits the horizon are frozen and can have little interaction with modes within the horizon. The large interaction only occurs among modes that are crossing the horizon at about the same time. These modes then have similar wavelengths. This is why the shape of the non-Gaussianity peaks at the equilateral limit in momentum space.
This physical origin also suggests the caveat that, as long as there are large interactions involving modes with similar wavelengths, an equilateral-like shape may arise. For example, such cases can happen in multifield models where there are particle creation [82,83]. (See however [84]).
• An example: Dirac-Born-Infeld (DBI) inflation. An explicit example of the above general results is the DBI inflation [54,74,[85][86][87][88][89][90][91][92]. These inflation models describe a 3+1 dimensional brane moving in warped extra dimensions. The location of the brane is a scalar field in 4d effective field theory, and it is the inflaton. The warped extra dimensions provide a non-trivial internal field space for the inflaton. In terms of the 4d effective field theory, the action is The non-trivial part is the kinetic term involving the square-root. It can be understood as a generalization of the following two familiar situations. It is a higher dimensional generalization of the action of a relativistic point particle where the speed of light f −1/2 varies with x. It is also a relativistic generalization of the usual canonical kinetic term in the non-relativistic limit |f (φ)g µν ∂ µ φ∂ ν φ| ≪ 1, Because the speed-limit of the inflaton f −1/2 can vary in the internal space, if it can be made small enough near the top of potential where the inflaton is about to roll down, the warped space restricts the rolling velocity even if the potential is too steep for slow-roll inflation to happen. So the inflaton rolls ultra-relativistically, but with very small velocity, and this generates the DBI inflation. The physical consequence is now easy to obtain using the general results in this subsection. In our notation the Lagrangian is The sound speed is which is the inverse of the Lorentz boost factor γ, so c s ≪ 1. The component (6.18) vanishes identically, and we have a large bispectrum of shape A c with size (6.19). DBI inflation is still driven by the potential energy. The general single field inflation models also include the k-inflation [93], where the inflation is driven by the inflaton kinetic energy. Model construction of single field k-inflation can be found in Ref. [93][94][95]. The bispectra for such models are computed in Ref. [50,94].
Multifield generalization have been studied in Ref. [96][97][98][99][100][101][102], where this type of kinetic terms are generalized to multiple fields. The three-point functions involving these different fields have the same or similar shapes.

Sinusoidal running: sharp feature
Although various slow-variation parameters in (6.5) have to be small most of the time during inflation, they can become temporarily large. Such cases can happen if there are sharp features in inflaton potentials or internal field space, so the behavior of inflatons temporarily deviates from the attractor solution, and then relaxes back within several Hubble time, or stay longer but with small deviation amplitudes. Motivations for such models include the following. It may be possible explanations for features in power spectrum [106][107][108][109], and if so the associated non-Gaussianity is a cross-check. And there are brane inflation models that are very sensitive to sharp features present in the potential or in the internal space [110].
As an example, we study a sharp feature in the slow-roll potential. The fact that a sharp feature in potential can enhance non-Gaussianities has long been anticipated and qualitative estimates have been made by different methods [103][104][105]. The precise method of analyzing the size, running and shape of such non-Gaussianities [64,65] is made possible with the developments of the formalisms that we reviewed in Sec. 3.1, 4 and 6.1. This will be the subject of this subsection.
We start by studying the behavior of the slow-roll parameters. We use a small step in potential as an example and will ignore numerical coefficients. We use c ∼ ∆V /V to denote the relative height of the step, and d the width of the step. In the attractor solution, the inflaton velocity is given byφ ∼ V ′ /H ∼ √ ǫV . As it falls down the step, the potential energy cV gets converted to the kinetic energy, so we havė The amplitude of density perturbations is given by P ζ ∼ H 4 /φ 2 , so such a sharp feature causes glitches in the power spectrum. It will leave a dip with relative size ∆P ζ /P ζ ∼ 1 + c/ǫ − 1 sinceφ increases first, followed by oscillations caused by a non-attractor component of the mode function before it settles down again in the attractor solution. To fit the CMB data,φ cannot change much. As we can see, the sensitivity of the power spectrum to the step size c is proportional to ǫ, and we need c/ǫ 1. Reducing the width d of sharp feature increases the amplitude of the glitches, but this is only for a large d over which the inflaton spends more than one e-fold to cross. Further reducing d will not change the amplitude of the glitches since (6.28) is saturated; but the sharpness will determine how deep within the horizon the modes are affected.
So ǫ does not change much, ∆ǫ ∼ ∆(φ 2 )/H 2 ∼ c. But it changes within a very short period, ∆t ∼ ∆φ/φ ∼ d/ V (c + ǫ). So η can be very large, It is also clear that the feature is associated with a characteristic physical scale and generates a scale-dependent power spectrum and higher order correlation functions. With these qualitative behavior in mind, we now study the three-point function. An important fact of the formalisms in Sec. 4 and 6 is that the expansion is exact in terms of the slow-variation parameters. So it is valid even if these parameters are not always small, as long as the expansion in ζ ∼ O(10 −5 ) is perturbative.
In all terms in the cubic expansion (4.8), ζ appears at most with one time derivative; the field redefinition gives a term that is proportional to η at the end of the inflation; and the other terms are all suppressed by powers of ǫ, which remains small even in the presence of a sharp feature. So the most important term is in which the coupling is proportional toη. The correlation function ζ(k 1 )ζ(k 2 )ζ(k 3 ) is dominated by Precise evaluation of this expression has to be done numerically. But it is not difficult to see the generic properties of bispectra associated with a sharp feature. For long wavelength modes that already crossed the horizon at the time of the sharp feature, k i τ ≪ 1, the mode function is already frozen and the integration (6.31) gives vanishing contribution. For short wavelength modes that are still well within the horizon, the modes are not affected if their momenta are larger than the inverse of the time scale characterizing the sharpness of changes in slow-roll parameters. The modes most affected are those which are near the horizon crossing. These modes are all oscillatory, ∼ e −ik i τ . As we have studied, η ′ is temporarily boosted, so it can be roughly approximated as several hat-functions that satisfy dτ η ′ = 0. Examples of such behavior are shown in Fig. 6. If we simply approximate the hat-functions by several delta-functions, η ′ ∝ δ(τ − τ * ), the integration (6.31) will give something like where k * ≡ 1/τ * is the scale corresponding to the location of feature, φ 0 is a phase and Comparing with the effect on the power spectrum, one can keep the size of glitches in the power spectrum small while make f N L large, for example, by fixing c/ǫ and decreasing d. This ansatz describes the most important running behavior of this bispectrum. Notice that the oscillatory frequency in the k-space is of order 1/k * , which is the scale of the feature. A rescale in k * can be compensated by a rescale in all k i . Also notice that the oscillatory frequency, 3/k * , along the k 1 = k 2 = k 3 ≡ k direction is 3/2 of that in the power spectrum 1 , 2/k * . In practice, (6.32) is a crude ansatz that needs to be refined. First of all, we have only considered the modes that have not exited the horizon. For those did, as we mentioned, their correlation function is as small as usual. The ansatz needs to be cut off for the long wavelength modes K/k * ≪ 1. A more detailed analysis [65] reveals, using the hat functions as an approximation of the slow-roll parameter behavior, that the bispectrum falls off as K 2 for these long wavelength modes. Secondly, the fact that in (6.32) all short wavelength modes are equally affected is due to the sharp-change approximation. Smoother functions will only affect a finite range of modes within the horizon. So the amplitude of the ansatz should decay and how fast depends on the sharpness of feature. To take into account both effects, empirically, we can multiply (6.32) with an envelop function where n and m are parameters chosen to fit the numerical results. For example, n = 2, m = 5 for Fig. 7. Lastly, in the very squeezed limit, k 3 ≪ k 1 |K/k * |, S can no longer be approximated as a function of K only and starts to have a non-trivial shape [65]. Here we concentrate on the signature running behavior. A numerical result with is shown in Fig. 7. A subtlety encountered in the numerical integration is how to handle the oscillatory behavior at τ → −∞. One can do a tilt into the imaginary plane, −τ → −∞(1 + iǫ), as prescribed in the analytical procedure in Sec. 3.2; or more efficiently, perform integration by part to increase the convergence of the integrand at the τ → −∞ end. One may also try the method of Wick rotation, but this will first require solving the background equations of motion in the Wick-rotated space, since here we do not have the analytical expression for the mode function. Sharp features can also appear elsewhere instead of potentials, for example, in the internal warped space for DBI inflation [110]. The qualitative running behavior in bispectrum is similar, and overall large non-Gaussianities become a superposition of the approximate scaleinvariant equilateral shape and the sinusoidal running.
Non-attractor initial conditions can be included as a case of sharp features, except that we only observe the relaxation part.

Resonant running: periodic features
In this subsection, we consider a different type of features. These features may or may not be sharp, but the most important property is their periodicity. Such features will induce an oscillatory component to the background evolution, in particular, to the couplings in the interaction terms. We denote this oscillatory frequency as ω. We know that each mode oscillates when it starts the life well within the horizon. This frequency keeps on decreasing as the mode gets stretched by the inflation, until it reaches H when the mode becomes frozen. So the mode scans through all frequencies that is larger than H, up to some very high cutoff scale such as M P . Therefore, as long as ω > H, (6.36) the oscillatory frequency of the modes in the integral will hit ω at some point during the inflation. This cause a resonance between the couplings and modes, hence a constructive contribution to the correlation function [65]. Without the resonance, as we encountered previously, the highly oscillatory modes simply average out within the horizon. In contrast to the previous mechanisms, here the non-Gaussianities are generated when modes are subhorizon. We now study the properties of such a non-Gaussianity, following Ref. [65].
To estimate the integral, we use the unperturbed mode function. Similar to the sharp feature case, we get In this case, we are interested in the region |Kτ | ≫ 1 in order to have resonance. So the last term dominates as long as the momentum triangle is not too squeezed so one of the k i 's becomes < 1/τ at the resonance point. The oscillatory coupling is dominantly contributed by η ′ . The integral is proportional to dτ τ sin(ωt) exp(iKτ ) . (6.38) This integral can be done analytically using the relation t ≈ −H −1 ln(−Hτ ). But its most important properties can be understood as follows in terms of the physical picture that we described.
First, let us look at its oscillatory running in K-space. The phase of the background repeats itself after ∆N e = 2πH/ω e-fold, during which the wave-number K changes by −K∆N e . So the running of the non-Gaussianity in K-space is also oscillatory with the period given by ∆K = K∆N e = 2πKH/ω . (6.39) Note that this period is changing with K in a specific way that we will see more clearly in a moment. Next, let us look at the size of the non-Gaussianity. Each K-mode briefly resonates with the oscillatory coupling when its frequency sweeps through the resonance frequency ω. Once its frequency differs from ω by ∆ω, the integration in the 3pt starts to cancel if is performed over ∆t 1 ∼ π/∆ω. In the meanwhile it takes ∆t 2 ∼ ∆ω/(ωH) to stretch the mode and change its frequency from ω to ω − ∆ω. Equating ∆t 1 and ∆t 2 gives the time period over which the resonance occurs for this mode, ∆t ∼ π ωH . (6.40) This corresponds to the number of oscillation periods that we need to integral over to estimate the resonance contribution. Note that one period in the integral (6.38) for K = ω contributes πτ * /K, where τ * is evaluated at the resonant point. Multiplying the total number of the resonant periods (6.41), using the definition (5.9) andP ζ ≈ H 2 /(8π 2 ǫ), we see that the amplitude of S(k 1 , k 2 , k 3 ) is Slow-roll parameters acquire small oscillatory components, and here η A denotes the amplitude of such an oscillation. Other pre-factors of k i are cancelled according the definition of S and the S turns out to be a function of K only. In the last step of (6.42), we have listed the accurate numerical number, which differs from the estimate by a factor of √ 2. Summarizing both the running behavior and the amplitude, we get the following ansatz for the bispectrum S res ansatz = f res N L sin (C ln(K/k * )) , (6.43) where C = 2πK/∆K = ω/H (6.44) and k * gives a phase. The argument C ln K in (6.43) appears because of (6.39). This gives a scale dependent oscillatory frequency in the K-space. In fact, this kind of dependence makes the density perturbations in the resonance model semi-scale-invariant. We call it periodic-scale-invariant -they are invariant under a discrete subgroup of rescaling. Namely, the ansatz (6.43) is invariant if we rescale all k i by n∆K/K = 2πnH/ω e-fold, where n is an integer. Other rescaling causes a phase shift. This property is a direct consequence of the symmetry of the Lagrangian. It is periodic, so invariant under a discrete shift of the inflaton field. This periodic-scale-invariance should also be respected by the full bispectrum results, as well as other correlation functions.
As mentioned, we have derived this ansatz from the last term in (6.37). Other terms will become important in the squeezed limit. The full integration (6.37) has been worked out in [111], and the leading order results are S res = f res N L sin (C ln(K/k * )) + where O(1/C 2 ) terms are neglected because we need 1/C = H/ω ≪ 1 for large resonance. The numerical coefficient in (6.42) turns out to be √ π/(8 √ 2). As we can see, the extra terms satisfy the symmetry we mentioned and indeed give large corrections in the very squeezed limit, e.g. k 3 < k 1 H/ω. These terms also ensure a consistency condition that we will study in Sec. 9.2. An example is plotted in Fig. 8. The spike at the very squeezed limit is due to the second term in (6.45). Overall, we see that the leading shape of this bispectrum is quite trivial, being almost a function of K only, until it gets to the very squeezed limit. The most distinctive feature of this type of non-Gaussianities is the running behavior captured in (6.43). Unfortunately, this ansatz is not factorizable if the K-range is too large.
More arbitrary scale-dependence can be introduced if the features are applied over a finite range, or with varying periodicity and amplitude.
As a useful comparison, the resonant running here and sinusoidal running that we studied in the last subsection are clearly distinguishable from each other observationally. The resonant running oscillates with periods that are always much smaller than the local scale, ∆K ≪ K; the frequency has a specific scale-dependence, ∆K/K = const.; and the frequency in the power spectrum (∼ sin(C ln k) in k-space) is exactly the same as that in the bispectrum (∼ sin(C ln K) in K-space). In contrast, the bispectrum of the sinusoidal run-ning oscillates with a fixed period that approximately equals to the scale at the location of the sharp feature, ∆K ∼ k * ; the frequency is scale-independent; and the power spectrum (∼ sin(2k/k * ) in k-space) has twice an oscillatory frequency of the bispectrum (∼ sin(K/k * ) in K-space).
As an illustration, we look at an example, In this example, the inflaton is rolling over the small but periodic ripple laid on the potential. This induces an oscillatory component in the slow-roll parameters with an amplitudeη A ≈ √ 6cmφ/Λ 2 and a frequency ω ≈φ/Λ ≈ 2m/( √ 6Λ). So we have and C ≈ 2/(φΛ) .

(6.48)
A numerical example is shown in Fig. 9. As we can see, the ansatz (6.43) gives a very accurate fit to the actual running behavior. The mode function and power spectrum are the superposition of the usual unperturbed solution and a small oscillatory component [65,112,113]. We can choose parameters so that the size of the ripples on the power spectrum is small, but bispectrum is made large. This is because the non-Gaussianities rise more quickly if we increase the frequency, while the mode function has difficulty responding efficiently when the external source oscillates too fast. This mechanism may be realized in terms of brane inflation [110] where the periodic feature comes from duality cascade in warped throat [114], or the monodromy inflation [115,116] where the periodic feature comes from instanton effects [112,113].

Folded shape: a non-standard vacuum
In this subsection, we study the effect of non-standard vacuum on the primordial non-Gaussianities. We consider a different wave-function from the Bunch-Davies vacuum when modes are well within the horizon. To start, let us first discuss several motivations for this case.
• A non-Bunch-Davies vacuum can actually occur much more simply than it might sound like. Any deviation from the attractor solution of the inflaton generically generates a component of non-Bunch-Davies vacuum. This is because a general mode function is a superposition of two components, c 1 (k)u(k, t) + c 2 (k)u * (k, t), and in attractor solution we choose one of the component asymptotic to the Bunch-Davies vacuum.
A disturbance will generically introduce a mixture with the other component. In this sense, we have already encountered such a case when we studied the effect of a sharp feature in Sec. 6.2. Indeed, after the inflaton crosses the sharp feature, the oscillatory behavior in the power spectrum is precisely due to the superposition of the second non-Bunch-Davies component for some finite k-range. For an infinitely sharp change, such a disturbance with a small amplitude extends to all k that have not exited the horizon at the time of sharp feature. An analytical illustration can be found in Sec. 5.3 of Ref. [110]. The location of the sharp feature can become superhorizon at the present time, but its influence has extended to much smaller scales and becomes observable. The resonance case in Sec. 6.3 is another type of example. An analytical illustration can be found in Sec. 3.3 of Ref. [112]. For non-Gaussianities studied in Sec. 6.2 & 6.3, we only concentrated on the effects caused by the features in slow-roll parameters. In analytical analyses, we approximated the mode function by the Bunch-Davies component and ignored the disturbance. The study of this subsection can be regarded as the complementary analyses on the effect of a different mode component.
• In inflationary background, modes can be quantized in terms of time-dependent creation and annihilation operators, a k (t) and a † −k (t). The Bunch-Davies vacuum is defined as the vacuum annihilated by a k (t) as t → −∞. If a different adiabatic vacuum is defined which is annihilated by a k (t 0 ) at a finite t 0 , we introduce a non-Bunch-Davies component. For example, see [117,118]. The origin and magnitude of such a component have been debated and studied by many papers, often under the name of the "trans-Planckian effect". See Ref. [119,120] for summary and references.
• There are inflation models where the scale of new physics can be very low. In particular, in warped space it is proportional to the exponentially small warp factor. In some DBI inflation models [86,92], After these discussions, let us now focus on a specific simple problem [50]. We modify the wave-function of the Bunch-Davies vacuum by a small second component and examine its consequence for the three-point function calculated in Sec. 6.1. We consider the general single field inflation with a small sound speed c s or a large λ/Σ [50,121].
So the mode function is a specific time or specific energy scale. To see a common feature without addressing these model-dependent issues, we look at the simple limit where the τ in (6.49) can go all the way to −∞. The computation of the correlation function is essentially the same as in Sec. 6.1. The leading order correction to the bispectra is obtained by replacing one of the three u k (τ ) in the integrand by its C − component. So it simply replaces one of the k i 's in the shapes with −k i . For example, the correction to S λ is (6.50) The shape ofS λ is shown in Fig. 10. The most important feature of this shape is the enhancement at the folded triangle limit, e.g. k 1 + k 2 − k 3 = 0. The detailed form of enhancement is model dependent. For example, it is different for another shape S c . The divergence in this folded limit occurs due to our simple limit of taking τ to −∞. Imposing some kind of cutoff at the lower limit of τ will eliminate this divergence, although as mentioned the detailed modification will be highly model dependent. For example, a simple constant cutoff τ c will introduce a factor of 1 + ( 1 2 x 2 c − 1) cos x c − x c sin x c for each of the three terms in (6.50), where x c ≡ (k 1 + k 2 − k 3 )c s τ c or its cyclic. Very close to the folded limit,Kc s |τ c | ≪ 1 (K = k 1 + k 2 − k 3 or its cyclic), this regulates away the divergence; away from the folded limit,Kc s |τ c | ≫ 1, these extra factors are unity on average but with oscillations. These oscillation can be either physical, or regarded to be zero if x c is within a regulation scale which exists since the non-Bunch-Davies component is present for a finite time in the past. The case for slow-roll inflation is qualitatively similar, and more examples of the bispectra shapes and the observational prospects are discussed in [122,123]. In this case, the proportional parameter for the bispectra amplitude is no longer enhanced by 1/c 2 s or λ/Σ, but < O(1).
In order to facilitate the data analyses, a simple ansatz has been proposed in Ref. [123], S fold ansatz,1 = 6 k 2 1 k 2 k 3 + 2 perm. − 6 k 1 k 2 + 5 perm. + 18 , (6.51) which represents certain important features of this kind of bispectra. It has a smooth rising behavior in the folded limit. This ansatz is plotted in Fig. 11(a). Since the real shape has a model dependent cutoff, it remains open questions how sensitive this is to data analyses and how well the ansatz (6.51) represents it. We can also write down an ansatz which is more directly motivated from the example (6.50) and the comments after that equation, where the cutoff scale k c = 1/(c s τ c ) is a parameter. For k 1 + k 2 − k 3 ≫ k c and cyclic, we have neglected the oscillatory part and only taken the average. In this ansatz we can change the powers in the numerator and denominator to model model-dependent variations. The relation (k c + k 1 + k 2 − k 3 ) −n = (Γ(n)) −1 ∞ 0 dt t n−1 e −(kc+k 1 +k 2 −k 3 )t may be used to factorize the ansatz. This ansatz is plotted in Fig. 11(b).
Another type of non-Bunch-Davies vacuum, namely an n-particle state built on the normal Bunch-Davies vacuum, was studied in [124,125] and the non-Gaussianities were found to be unobservable.

Quasi-single field inflation
Having considered single field inflation, we now relax the condition on the number of fields. At least during inflation, we only need to consider quantum fluctuations of light fields, since if the mass of fields are very heavy, (here the relevant scale is m ≫ H), they contribute only classically and determine the classical inflaton trajectory. Multiple light fields can arise naturally if we consider the inflation models as the consequence of a UV completed framework. However, as discussed in Sec. 5.1, due to the back-reaction from the inflationary background, the mass of light fields are naturally of order H. The potential with such a shape is too steep for slow-roll inflation.
Therefore, as a natural step beyond the single field, let us consider slow-roll models with one inflationary direction, and one or more other directions that have mass neither much heavier nor much lighter than H. We will call the quanta in the inflationary direction as the inflaton and its mode the curvature mode, and the others isocurvaton and isocurvature modes. We call these models the quasi-single field inflation models [126,46].
Note that the thematic order in this review is not chronological. The non-Gaussianities in this type of models were not computed until very recently for a couple of reasons. If the mass of particles is of order O(H) or larger, the amplitude of these fields decay exponentially in time after horizon-exit. So they would not seem to be important for super-horizon perturbations even if they couple to the curvature mode. As we will see, however, their amplitudes at or near the horizon-exit are enough to make them interesting. What really suppresses their contribution is the fast oscillation behavior present for m ≫ H. Methodologically, isocurvature-to-curvature transition for non-Gaussianities was studied restrictively in the regime of super-horizon classical evolution in multi-field space [127][128][129][130][131][132][133][134], which we shall explain in more details in the next section. However, for quasi-single field inflation models, a full quantum computation in the in-in formalism is necessary to properly include the contributions from both the horizon exit and the superhorizon evolution.

Intermediate shapes: massive isocurvatons
There are potentially different ways massive isocurvatons can be coupled to the inflaton. We currently do not have a general approach in terms of model building. So what we shall do is to first study this problem through a simple example, and then discuss the features of the results that can be regarded as generic signatures of this class of models [126,46].
We consider the case where the inflaton is turning constantly by going around (a fraction of) a circle with radius R in the angular θ direction. See Fig. 12. All the parameters, such as R and couplings, are assumed to be constant during the turning. We call this assumption the constant turn case. In the θ direction the potential is the usual slow-roll potential V sr (θ). The field in the radial direction is denoted as σ and has mass of order H, and lifted by the potential V (σ). For such a turning trajectory, it is convenient to write the action in terms Figure 12: Quasi-single field inflation with turning trajectory. The field θ and σ are in the polar coordinates. The θ is the inflationary direction with a slow-roll potential. The σ is the isocurvature direction, which typically has mass of order H.
of fields in the polar coordinates, θ and σ, instead of in the Cartesian coordinates, The potential V (σ) balances off the centrifugal force necessary for the turning and traps the field at the bottom of the effective potential, V eff (σ) = − 1 2θ 2 0 (R + σ) 2 + V (σ). We define the minimum of this effective potential to be σ = 0. We expand the effective potential as whereθ 0 is the turning angular velocity and the primes on V denote derivatives with respective to σ.
To study the perturbation theory, we perturb the fields in the spatially flat gauge, and obtain the following Hamiltonian,

4)
H I 2 = −c 2 a 3 δσ Iδ θ I , (7.5) where are all constants. Terms suppressed by O(ǫ) have been ignored in this gauge. The curvature perturbation ζ is most transparent in another gauge, the uniform inflaton gauge, where θ(x, t) = θ 0 (t) , σ(x, t) = σ 0 (t) + δσ(x, t) , (7.8) and the spatial metric is In this gauge, ζ appears in the metric as the space-dependent rescale factor and the fluctuations in the inflaton is shifted away. The relation between ζ and δθ is the gauge transformation. At the leading order this is We will calculate the correlation functions in terms of δθ and then use this relation to convert them to those of ζ. The full perturbation theory that one obtains in the uniform inflaton gauge justifies the above omission of several O(ǫ) terms in the spatially flat gauge [46]. There are several important points for this Hamiltonian. First, the kinematic Hamiltonian (7.4) describes two free fields in the inflationary background. One is massless and has the familiar mode function, Another is massive and the mode function is where ν = 9/4 − m 2 /H 2 . and is equivalent to factors of Boltzmann-like suppression ∼ e −m/H . We will consider the case 0 ≤ ν ≤ 3/2. Second, there is a sharp contrast between the V ′′′ sr for the slow-roll inflaton field and the V ′′′ for the massive field σ in the non-inflationary direction. The former has to be very small, ∼ O(ǫ 2 )P 1/2 ζ H, in order to maintain the smallness of the slow-roll parameters. (Here we use ǫ to denote collectively all slow-roll parameters, ǫ ≡ −Ḣ/H 2 , η ≡ǫ/ǫH, and ξ ≡η/ηH.) Consequently it contributes O(ǫ 2 ) to the f N L of bispectrum in slow-roll inflation, generally smaller than the O(ǫ) contributions from the other terms in the same model. However, for quasi-single field inflation, the direction orthogonal to slow-roll does not have to satisfy the slow-roll conditions, and V ′′′ is almost unconstrained. For example, in the inflationary background, it can be of order H; and similarly, V ′′′′ can be of order one, etc. This isocurvaton self-interaction (7.6) becomes the source of large non-Gaussianities.
Third, the coupling between the isocurvaton and inflaton appears as a form of a two-point vertex operator in (7.5). We treat this term as part of the interaction Hamiltonian, and it is represented by the transfer vertex in Fig. 13 (a). The strength of the coupling is determined by the turning angular velocityθ 0 in this model. This coupling is responsible for the transformation of the isocurvature perturbations, in particular their large non-Gaussianities, to the curvature perturbation.
We calculate correlation functions corresponding to the Feynman diagrams Fig. 13 in terms of the in-in formalism, which we reviewed in Sec. 3.1. As an illuminating example to illustrate the different advantages of the three forms of the in-in formalism, we recall from Sec. 3.1 that the three-point function can be written in the following forms. The original definition (3.44) and (3.40), which we refer to as the factorized form, leads to × (2π) 3 δ 3 ( i p i ) + 9 other similar terms + 5 permutations of p i . (7.14) The perturbation theory here starts from the fourth order. The reorganized commutator form (3.41) leads to In the IR (τ → 0), each of the ten terms in the factorized form diverge as τ 3−6ν for 3/2 > ν > 1/2 (0 < m < √ 2H); while in the commutator form, various subtractions off the complex conjugates and the requirement that the final result has to be real make such divergence explicitly disappear.
In the UV (τ → −∞), each factor of the multiple integral that integrates from −∞ to 0 has a definite convergent direction if we choose one of the two contour tilts, τ i → −∞(1±iǫ), accordingly. Or more efficiently, by a Wick rotation τ i → ±iz i . This would have been the case for the commutator form if we can break up the integrand into individual terms. However in order to achieve the explicit IR convergence, as we saw above, these terms have to be grouped; but then they have contradicting convergence directions.
To take advantage of both forms, we introduce a cutoff τ c , and write the IR part (τ c < τ ≤ 0) of the integrals in terms of the commutator form, and the UV part (τ < τ c ) in terms of the factorized form, in the following mixed form (3.42),  To better understand the shapes analytically, we can work out the squeezed limit (p 3 ≪ p 1 = p 2 ) of the three-point function, where s(ν) is a ν-dependent numerical number.
Recall that the squeezed limit of S for the equilateral shape goes as p 3 /p 1 , while for the local shape (p 3 /p 1 ) −1 . Both the numerical results in Fig. 14 and the analytical results in Eq. (7.17) show that here we have a one-parameter family of shapes, ∼ (p 3 /p 1 ) 1/2−ν , lie between the two. We call them the intermediate shapes.
The physical origin of such shapes can be understood as follows, and should be a generic signature for the quasi-single field inflation models. As we have seen, the large equilateral non-Gaussianity arises because the interacting modes cross the horizon around the same time. The shape of bispectrum peaks at the equilateral limit where the modes all have comparable wavelengths. As we will see in Sec. 8, the large local non-Gaussianity arises due to the classical non-linear evolution of superhorizon modes in the multifield space; so the interactions are causally disconnected and behave local in position space. This is non-local from the momentum space point of view. So the shape of bispectrum peaks at the squeezed limit. Now for quasi-single field inflation, the large non-Gaussianities come from the massive isocurvaton. Depending on the mass, these modes either decay right away after they exit the horizon (for m > √ 2H), or survive for a long time at the super-horizon scales (for m < √ 2H). In the former case, the generation and transfer of non-Gaussianities maximize for modes that are exiting the horizon around the same time, resulting in quasi-equilateral shapes; in the latter case, the generation and transfer of non-Gaussianities happen in a superhorizon fashion, resulting in quasi-local shapes. In this regard, let us look more closely at the special limit m/H → 0 (ν → 3/2).
In this massless limit, an infrared cutoff to the integrals are necessary. Otherwise the transfer will last for ever for the constant turn case. The cutoff corresponds to the ending of the turning. Let us discuss the following two cases. First, we still keep V ′′′ large. Our analyses still apply in this case. Interestingly, the shape of the bispectrum goes to that of the local form in this limit. As we will explain in Sec. 8.1, this is a generic signature of a massless isocurvaton. The infrared e-fold cutoff will introduce some running in the f int N L because different modes experience different turning e-folds. Second, we would like to make the isocurvature directions flat so this becomes a two-field slow-roll inflation models. Such models were intensively studied and it is known that the isocurvature modes can be transferred to the curvature mode by turning. However, since V ′′′ ∼ O(ǫ 3/2 )H 2 /M P is required to maintain the small slow-roll parameters, the contribution we computed here generates too small non-Gaussianity. We expect contributions from other terms are small as well. So it is much more difficult to generate large non-Gaussianities in such models, essentially because imposing the slow-roll conditions in all directions are too restrictive.
To connect with data analyses, guided by the numerical results and analytical squeezed limit, we can use the following ansatz to describe the full family of shapes, These shapes are shown in Fig. 15. Comparing with Fig. 14, we can see that they match quite well except near ν = 0, around which it is better represented by another ansatz in Ref. [46]. The size of this bispectrum is where P ζ ≈ 6.1×10 −9 andθ 0 is the turning angular velocity. The α(ν) is a positive numerical number which, depending on ν, can give an additional enhancement factor of order N f (N f is the turning e-folds). Since (θ 0 /H) 2 and V ′′′ /H are the expansion parameters in the perturbation theory, they have to be small to trust our calculation. Nonetheless this is not the model-building requirement. The fluctuations of more massive (m > O(H)) fields may become important if they play a role later in the reheating [135,136]. Such cases typically require some tunings for special conditions, so that the highly suppressed fluctuation amplitude can become important.

Multifield inflation
As we have seen in Sec. 7.1, if we take the isocurvaton mass to zero in quasi-single field inflation while keep the nonlinear self-couplings of the isocurvaton V ′′′ large, the shape of the large bispectrum in the squeezed limit approaches the local form. The local form is in fact the earliest and most well-studied example of non-Gaussianities [127,137,138,66], although it was first found to be small as we have seen in Sec. 4. As we will explain in this section, a large local form is a signature of massless isocurvatons that have large non-linear evolution in multifield space. We have arrived this shape from the in-in formalism by taking the massless limit. But if we stay in this limit, there is an easier formalism, the δN formalism [139][140][141], in which the underlying physics of the local shape becomes transparent.

Local shape: massless isocurvatons
We recall that, in single field inflation, if we use the uniform inflaton gauge where there are no fluctuations in the inflaton field, the scalar perturbation ζ enters in the scale factor as a 2 e 2ζ . For superhorizon modes, ζ is frozen. If we look at the different comoving superhorizon patches, they are causally disconnected from each other. So they evolve independently and locally in space. In such a gauge, the only difference is a space-dependent scale factor. This is also called the separate universe picture. The primordial curvature perturbations manifest themselves as the different number of expansion e-fold, δN, at different positions.
We would like to generalize this picture to the multifield case in the following δN formalism. We will resort to a simple version of δN formalism stated below, which is of course a consequence of the in-in formalism, but formulated from a simple perspective which clearly illustrates the points in this section. Otherwise, as we will explain, in the most general sense the δN formalism should be defined as the in-in formalism written in terms of specified gauges.
• We consider a set of scalars φ i during inflation. Inflaton is one of them but can be different linear combinations of φ i 's as a function of time, and the other orthogonal fields are called the isocurvatons. All the modes that we are eventually interested in should all have become superhorizon when the initial slice (specified below) is chosen. We look at different horizon-size patches and label them with the coarse-grained comoving coordinate x. In the in-in formalism, the superhorizon modes behave as the c-number time-dependent background for each comoving patch. So we evolve these patches independently and classically.
• We pick an initial spatially flat slice, on which there is no scalar fluctuations in the metric and all the fluctuations are in the scalar fields φ 0i + δφ i (x). We assume that we know the statistics of such fluctuations.
• We pick the final uniform density slices. Relative to the unperturbed and perturbed initial spatially flat slices, we have, respectively, the unperturbed and perturbed final uniform density slices. For single field inflation, these two final surfaces are the same. For multifield models, they are generally different. Such final slices have the properties that the universe has the same energy densities and field configurations everywhere on them. They can be chosen during either the inflation or the reheating. After that, every separated universe will have the same evolution. The only difference is the scale factor. This is the analogy of the uniform inflaton gauge in single field inflation. We study the cases where such slices exist.
• We evolve the unperturbed φ 0i in the initial slice classically to the unperturbed final slice, and denote the number of e-folds as N 0 (φ 0i ). This is the unperturbed e-fold number. We evolve the perturbed φ 0i + δφ i (x) in the initial slice classically to the perturbed final slice, and denote the number of e-folds as N(φ 0i + δφ i (x)). The difference between them is the curvature perturbation ζ. Here N 0 is a constant that can be shifted to make δN zero.
• We expand where the subscripts on N denote the partial derivatives with respect to φ i . For example, N ij = (∂N/∂φ i )(∂N/∂φ j ). Repeated indices are summed over. The correlation functions of ζ can then be computed as the classical averages of the products, such as, • We have assumed that the statistics of the δφ i (x) are known on the initial slice. But this is not always easy to get. So we will consider the simple cases where this statistics can be approximated as Gaussian. Otherwise, calculating such initial statistics requires using the full quantum mechanical in-in formalism.
Most generally, one identifies δN with the scalar curvature ζ in the uniform inflaton gauge; and the relation between δN and δφ(x) in the δN formalism is the gauge transformation between the uniform inflaton gauge and the spatially flat gauge. Calculating the correlation functions for ζ becomes calculating those for δφ(x) using the in-in formalism. An example is the one we have seen in Sec. 7.1.
• So far we have not used the condition that the isocurvatons are massless (m ≪ H). If they are massive, after horizon exit the modes decay. So in the superhorizon classical regime, where the δN formalism is supposed to be useful, we are back in the single field inflation. Sub and near horizon perturbations should be computed by the full in-in formalism. Therefore having massless isocurvatons opens up classical multifield space in which we can have sizable δN defined in (8.1). Now let us consider the Gaussian fluctuations δφ i . From Sec. 2, we know that for massless scalars, where H * is the Hubble parameter when the corresponding mode exits the horizon. If the scalars are not exactly massless, H will have a running dependence on k 1 caused by the decay of the amplitude. Using we get the power spectrum where 8) and the bispectrum According to the definition (5.9), the shape function is As usual, we have ignored a mild running from the power spectrum P ζ . The shape (8.10) is called the local shape, which we plot in Fig. 16. This form is already factorizable. The physics of this shape can be understood from the derivation above. As explicitly demonstrated in (8.1)-(8.4), this non-Gaussianity is generated locally in position space for superhorizon modes. After Fourier transform, it becomes non-local in momentum space. That is the reason that the shape peaks at the squeezed limit.
If the perturbation δφ(x, t) on the initial spatially flat slice cannot be approximated as Gaussian, the shapes of final bispectra can be more complicated.
In this model, we assume that during inflation there is another light field σ with the potential and m σ ≪ H. This field is called the curvaton field for reasons that will become clear shortly. The energy density of the curvaton field is negligible initially. During inflation, it fluctuates and obtains the primordial amplitude σ * = σ 0 +δσ(x), where * denotes its value at the horizon-exit and after that the amplitude is approximately frozen. These perturbations are Gaussian for the potential (8.12) with the canonical kinetic term, but can be more complicated otherwise. Here we study the simple Gaussian case. After inflation, these σmodes remain frozen until the Hubble parameter drops below m σ . Then the σ-field starts to oscillate around the bottom of the potential and behavior as matter. The Universe is still radiation-dominated. The fraction of the matter energy density quickly grows, because the matter dilutes as a −3 while radiation a −4 . The σ-field decays to radiation when it reaches its lifetime.
Another assumption of the curvaton model is that the primordial fluctuations in the inflaton field is much smaller than what is needed to achieve ζ ∼ O(10 −5 ), although their total energy density may still be the dominant one. So the primordial curvature perturbation is contributed by the fluctuations in the σ field, hence the name curvaton field.
At the initial spatially flat slice t 0 , we denote the radiation and curvaton density as ρ r0 , ρ m0 , respectively, and the scale factor as a 0 . Both components initially redshift as radiation. This lasts until the Hubble parameter reaches m σ at t 1 . We denote the scale factor at t 1 as a 1 . The Friedman equation at t 1 is After this the curvaton starts to oscillate and behave as matter. Denote the decay rate of the curvaton as Γ. We use the sudden decay approximation and assume that they decay instantaneously at the epoch H = Γ, because a process that has the time scale T falls into the Hubble expansion epoch with H = 1/T . We denote this instant as t 2 . The Friedman equation at t 2 is Because at t 2 the universe has the same Hubble parameter (hence the same energy density) everywhere, this is the final uniform density slice. After that, both components become radiation and the evolution everywhere is the same. So we want to work out the expansion efolds N from t 0 to t 2 as a function of the initial field value σ. From (8.14), we get e −4N + e −3N α = const. , (8.15) where α ≡ (a 0 /a 1 )(ρ m0 /ρ r0 ). From (8.13) we can solve for a 0 /a 1 which is independent of σ at the leading order, and from (8.12) we know that ρ m0 is proportional to σ 2 . Therefore α is proportional to σ 2 . Also note that r ρ ≡ (a 2 /a 0 )α = e N α = ρ m /ρ r | t 2 is the ratio of the energy density between the curvaton and the rest of the radiation at t 2 . Using these simple facts, we can differentiate (8.15) with respect to σ once and twice, and get In terms of the definition r ≡ 3ρ m /(4ρ r + 3ρ m )| t 2 = 3r ρ /(4 + 3r ρ ) often used in the literature, So Large local non-Gaussianity arises if r ≪ 1. Note that although (8.15) only depends on σ, this is a multifield model because the curvaton takes effects during the reheating. In some simple models in which the curvaton leads to non-adiabatic perturbations between dark matter and photons, r is tightly constrained by observations [142].
The large local form has been studied most extensively in the past. Variety of possibilities exist. They all share the common feature that non-Gaussianities are generated by some massless isocurvaton fields which acquire the superhorizon evolution during the inflation. For example, in multifield slow-roll inflation a turning trajectory [152] can transfer non-Gaussianities from other directions to the inflaton direction [127][128][129][130][131][132][133][134]. But it is found to be very difficult to make non-Gaussianities large essentially because the very restrictive slow-roll conditions in all directions. In modulated reheating [153,154] or preheating [155][156][157][158][159] scenarios, the role of isocurvatons are played by the massless fields which control the couplings during the reheating or preheating. Thus they create a large local non-Gaussianity in a similar fashion as the curvaton model.
9 Summary and discussions

Summary
In this subsection, we summarize the main results of Sec. 6 -8. Non-Gaussianities, conceptually being the expectation values of perturbations in a time-dependent background, are defined by the first-principle in-in formalism. Physically, having large primordial non-Gaussianities means that there are large non-linear interactions of some sort determined by certain dynamics during inflation. Measuring them tells us the nature of the dynamics.
• Equilateral shape and higher derivative kinetic terms.
In single field inflation, the long wavelength modes that already exited the horizon are frozen. They cannot have large interactions with short wavelength modes that are still within the horizon. When modes are well within the horizon, they oscillate and the contributions to non-Gaussianities average out. Therefore the only chance to have large interaction is when all modes have similar wavelengths and exit the horizon at about the same time. Theories with higher derivative kinetic terms provide such interaction terms. This is why the resulting bispectrum shape peaks at the equilateral limit in the momentum space. It drops to zero at the squeezed limit k 3 ≪ k 1 = k 2 as k 3 /k 1 . It happens that, when these higher derivative terms become important enough so that the inflationary mechanism is no longer slow-roll, these large non-Gaussianities become observable. The forms of the bispectra are given in Eq. (6.16) and (6.17), and the shapes are plotted in Fig. 3 and 4. The factorizable ansatz that is used to represent them in data analyses is given in Eq. (6.22) and plotted in Fig. 5.
• Sinusoidal running and sharp feature.
A sharp feature, in a potential or internal field space, introduces a sharp change in the slow-roll parameters, or the generalized slow-variation parameters. This can boost the magnitudes of time-derivatives of some parameters by several orders of magnitude while still keep the power spectrum viable. These time-derivatives act as couplings in the interaction terms. So they enhance the non-Gaussianities among the modes which are near the horizon-exit. How deep they affect the modes inside the horizon depends on how sharp the changes are.
The changes in these parameters can be roughly approximated as delta-functions in time. Correlation functions involve integrations of products of the slow-variation parameters and the mode functions. The latter contain oscillatory behavior ∼ e −iKτ , where the comoving momentum K is k 1 + k 2 + k 3 for bispectra and τ is the conformal time. The delta-function specifies a scale τ * . This is why after integration the bispectrum contains a sinusoidal factor ∼ sin(K/k * ), where k * = −1/τ * is the momentum of the mode that is near the horizon-exit at the time of the feature. So the most important property of this type of non-Gaussianity is this characteristic running. A numerical result of the running behavior is plotted in Fig. 7. An ansatz is given in Eq. (6.32) and (6.34).
• Resonant running and periodic features.
The periodic features do not have to be sharp. They introduce a small background oscillatory component in the slow-variation parameters. On the other hand, the mode functions are also oscillatory before they exit the horizon. Their frequencies are high when they lie deep inside the horizon and become lower as their wavelengths get stretched by the inflation. They are frozen after the wavelengths become comparable with the horizon size H −1 . This means that their frequencies continuously scan through the range from M P (or some other large fundamental scale) to H. Therefore as long as the background oscillatory frequency ω satisfies H < ω < M P , at some point during the evolution the small oscillatory component in the slow-variation parameters will resonant with the oscillatory behavior of the mode functions, and cause a large constructive contribution to the integration.
The periodicity of the features leads to a periodic-scale-invariance in density perturbations. Namely, they are scale invariant if we rescale all momenta by a discrete e-fold 2πnH/ω, where n is an integer. This is why the most important feature of this non-Gaussianity is a running behavior ∼ sin(C ln K + phase), where C = ω/H. This leads to the ansatz (6.43). The full expression is given in (6.45) and plotted in Fig. 8. A numerical result is plotted in Fig. 9.
• Folded shape and non-Bunch-Davies vacuum.
The usual mode function of the Bunch-Davies vacuum has the positive energy mode ∼ e −ikτ . Now we consider a non-Bunch-Davies vacuum by adding a small component of negative energy mode ∼ e ikτ . The three-point function involves an integration of the product of three mode functions with momentum k 1 , k 2 and k 3 . So the leading correction to the Bunch-Davies results is to replace one of the k i 's with −k i . The usual K = k 1 + k 2 + k 3 in e −iKτ becomesK = k 1 + k 2 − k 3 and its cyclic. This effect is most important if factors ofK appear in the denominators after the τ -integration. Hence the most important feature of this type of modification is to enhance the non-Gaussianity in the folded triangle limit. An example of these bispectra is given in Eq. (6.50) and plotted in Fig. 10. Ansatz that partially capture this feature are given in Eq. (6.51) and (6.52), and plotted in Fig. 11.
All mechanisms discussed so for single field inflation apply to multifield inflation. We now consider new effects caused by introducing more fields to inflation models. These extra fields are called isocurvatons.
Since light fields typically acquire a mass of order H, the Hubble parameter, we first consider the quasi-single field inflation models where there is one massless inflaton while the isocurvatons have mass of order H instead of massless.
Unlike multifield slow-roll inflation, where each flat direction only has small non-linear terms in order to satisfy the slow-roll conditions, massive directions are not inflationary direction and are free to have large non-linear self-interactions. These non-linear interactions can be transferred to the curvature mode through couplings and source the large non-Gaussianity.
The massive isocurvaton eventually decays after horizon exit simply because they are diluted by the expansion. How fast it decays depends on its mass. If the mass is heavier, m > √ 2H, it decays faster. So the interactions can only happen when all modes are all closer to the horizon exit. This is closer to the case of the equilateral shape that we encountered above, and results in bispectra with quasi-equilateral shapes. If the mass is lighter, m < √ 2H, it decays slower. More non-Gaussianity is generated in the super-horizon scales. This is closer to the case of the local shape that we will come to below, and results in bispectra with quasi-local shapes. Overall, at the squeezed limit k 3 ≪ k 1 = k 2 , the bispectrum shapes behave as (k 3 /k 1 ) 1/2−ν , where ν goes from 0 to 3/2 (corresponding to m from 3H/2 to 0) in the example we studied. In particular, if we take the massless limit while keeping the cubic self-interactions of isocurvaton large, we get a large bispectrum that has the same squeezed limit shape as the local one. Therefore, we have a one-parameter family of shapes that lie between the local and equilateral shape.
The numerical results of these shapes are presented in Fig. 14. A simple ansatz is given in Eq. (7.18) that represents this family of shapes quite well, and is plotted in Fig. 15.
The fluctuation amplitudes of massless scalars do not decay after the horizon exit, and therefore this opens up a multifield space for the superhorizon evolution. For superhorizon modes, we can use the separate universe picture and study the classical behavior of different patches of universe. These patches are separated by horizons and evolve independently of each other. So the evolution is local in space.
Non-Gaussianities are generated when this multifield evolution is nonlinear, and any nonlinearity arising in the separate universe picture should also be local in space. A locality in position space translates to a non-locality in momentum space. This is why the resulted local shape bispectrum peaks at the squeezed limit. The behavior is (k 3 /k 1 ) −1 for k 3 ≪ k 1 = k 2 . This bispectrum is given in Eq. (8.10) and the shape is plotted in Fig. 16.
In all cases, the power spectra are either approximately scale-invariant so indistinguishable from the simplest slow-roll models, or modified with features that can be made small enough to satisfy the current observational constraints.
It is certain that this list will grow in future works, providing more refined and diverse connections between theories and experiments.

A consistency condition
As we have seen, in single field inflation, the mode that has exited the horizon is frozen. This is characterized by a constant ζ over a horizon size patch. The physical meaning of the constant ζ is a small rescaling of the scale factor. This is the only effect that the superhorizon mode has on modes with much shorter wavelength. This fact is used by Maldacena to derive a consistency condition [47] for the three-point correlation functions in the squeezed limit for single field inflation.
• Consistency condition. In the squeezed limit k 3 ≪ k 1 = k 2 , k 3 is the superhorizon mode that exited the horizon and acts as a zero-mode modulation to the two remaining modes. The correlation ζ k 1 ζ k 2 ζ k 3 is an average of the following quantity ζ k 1 ζ k 2 ζ k 3 ζ k 3 (9.1) over different ζ k 3 . We will shift the average ζ k 3 to zero by definition. Here the two-point average ζ k 1 ζ k 2 ζ k 3 is evaluated with different local scalings determined by the shift ζ k 3 . If the two-point function is exactly scale-invariant, ζ k 1 ζ k 2 ζ k 3 is just a constant. So the 3pt vanishes because ζ k 3 = 0. The non-zero contribution comes from the breaking of the scaleinvariance. To see this, we expand the two-point average in terms of a long wavelength mode ζ k 4 near the scale ζ k 4 = 0, Multiply this with ζ k 3 and average over ζ k 3 . The first term contributes zero since this is the scale-invariant component. The second term gives To get non-zero average, k 3 + k 4 = 0 is needed. Using the relation dζ k 4 = −d ln k, we get The higher order terms in (9.2) give 1 2 where k 3 + 2k 4 = 0, k 3 + 3k 4 = 0, and so on have to be satisfied respectively for each term to get non-zero average. If we only consider the tree-level three-point function, these higher order terms can be truncated since they involve more factors of P ζ and should be related to the loop diagram contributions to ζ k 1 ζ k 2 ζ k 3 . The tree diagram is O(P 2 ζ ). To connect the averages we used here with the correlation functions that we defined in previous sections, we need to restore the phase factors. Here the two-point average ζ k 1 ζ k 2 here is performed with the special point k 1 = k 2 in the phase space. To connect this with the previous definition of ζ k 1 ζ k 2 , we need to include the phase space in the neighborhood. Namely, ζ k 1 ζ k 2 = ζ k 1 ζ k 2 here (2π) 3 δ 3 (k 1 + k 2 ), so from the definition (2.26) we have ζ k 1 ζ k 2 here = (2π) 2 P ζ (k 1 )/2k 3 1 . Similarly, ζ 3 = ζ 3 here (2π) 3 δ 3 ( k i ). With the usual definition of the spectrum index n s − 1 ≡ d ln P ζ /d ln k, from (9.4) we get the following consistency condition [47] Although originally derived for slow-roll inflation, the only assumption is the single field. So this applies to any single field inflation models and has important physical implications that we discuss shortly [198]. Note that the derivation of this relation (9.6) does not rely on the smallness of the slow-variation parameters either. For the general single field inflation models that we studied in Sec. 6.1, at tree level this has been checked with explicit results to three different orders [199,50] including the slow-roll limit [47]. For resonance models, this is checked to the leading order [111].
There are three types of interesting corrections to the condition (9.6). Firstly, as mentioned, the right hand side of (9.6) should receive corrections from loop contributions. These loop contributions are associated with higher derivatives of the twopoint function. The terms (9.5), together with (9.4), provide the corresponding consistency conditions at different orders of P ζ . Note that for such orders, the correlation functions such as ζ k 1 ζ k 2 on the right hand side should also include loop corrections.
Secondly, when we assume that the only effect of the frozen superhorizon mode on the much shorter scale is a constant background rescaling, we assume that there is no interaction when these modes are all within the horizon. 2 However, large subhorizon interaction is possible in some cases, such as in Sec. 6.3 and 6.4. Such interactions disappear below a new length scale at subhorizon, only then the above assumption becomes valid. For example in the resonance model, for H/ω < k 3 /k 1 ≪ 1, the modes k 1 , k 2 and k 3 are guaranteed to resonant with the oscillatory background at some point when all of them are within the horizon. Only if k 3 /k 1 ≪ H/ω, such resonance will not happen. This is why the consistency condition is satisfied only in the very squeezed limit (k 3 ≪ k 1 /C with C = ω/H ≫ 1). For the squeezed region k 1 /C < k 3 ≪ k 1 , the left hand side of the condition is larger than the right hand side by a factor of Ck 3 /2k 1 . For the non-Bunch-Davies vacuum case, a similar scale is determined by the UV cutoff τ c , from which the non-Bunch-Davies vacuum starts to take effect.
Thirdly, even after the long wavelength mode exits the horizon, as long as k 3 /k 1 is not infinitely small, there is still dependence of the two-point function on the derivative of this long wavelength mode, in addition to the overall constant shift. This introduces a different type of finite k 3 /k 1 corrections. They start from the second order in k 3 /k 1 , because the first order corresponds to the first spatial derivative of the long wavelength mode and should vanish due to isotropy [198]. These corrections will be amplified by the associated amplitude f non−loc N L , and give an additive correction f non−loc N L (k 3 /k 1 ) 2 to the n s − 1 on the right hand side of the condition. For a large f non−loc N L , therefore, the condition needs to be satisfied in a very squeezed limit. The equilateral bispectra (6.16) and (6.17) are this type of examples.
The consistency condition (9.6) can be straightforwardly generalized to higher order correlation functions [192,187]. We emphasize that this condition only applies to single field inflation models. For inflation models involving more than one field, as we have seen, non-Gaussianities can be transferred from the isocurvature directions which do not respect this relation.
• Physical implication. Besides providing consistency checks for analytical computations, the condition also has interesting physical implications. In the following, we discuss the scale invariant cases [198], as well as the feature cases and loop corrections, ending with some cautionary remarks.
This consistency relation implies that the tree-level bispectrum in the squeezed limit is determined by the power spectrum and spectral index. We distinguish the following two cases. For the scale-invariant case, n s − 1 is of order O(ǫ) and the right hand side of (9.6) takes the local form. Indeed, as we have seen, for the single field inflation models where the non-Gaussianities are large, they take the equilateral forms which vanish in the infinitely squeezed limit. For the non-scale-invariant case, especially the highly oscillatory case such as the resonance model, the power spectrum can be highly oscillatory and n s − 1 becomes large. This can still be consistent with observations since the large n s − 1 is also highly oscillatory and therefore may escape a detection so far. But such a running non-Gaussianity is orthogonal to the scale-invariant forms.
For the loop diagrams, in the scale-invariant case, these terms are suppressed by higher orders of slow-variation parameters from, e.g. d 2 ζ k 1 ζ k 2 /(d ln k) 2 , and higher orders of ζ from, e.g. ζ 3 k 3 ; in the non-scale-invariant example, the extra terms are still highly oscillatory. In summary, a detection of an approximately scale-invariant local non-Gaussianity in the infinitely squeezed triangle limit with f loc N L > O(ǫ) can rule out all single field inflation models.
In experiments, however, the triangle cannot be perfectly squeezed. So it is an important question how squeezed it should be to achieve the above goal. For example, in the third type of corrections we discussed previously in this subsection, we need f non−loc N L (k 3 /k 1 ) 2 to be smaller than n s − 1 for the consistency condition to hold, so that the contaminations from whatever non-local f N L is small. Assuming the primordial local form is practically detectable only if f loc N L > O(1), we at least need f non−loc . For the class of the general single field models we studied in Sec. 6, if the other forms of non-Gaussianities, such as the equilateral one, can be constrained below f non−loc N L ∼ O(10), a squeezed configuration with k 3 /k 1 < 0.1 will be enough for our purpose. However, a completely model-independent statement is much trickier, because there may be bispectra with very large amplitude but orthogonal to any known bispectra that have been constrained experimentally. Besides that, in the second type of corrections, large finite-k 3 /k 1 corrections can also arise due to subhorizon interactions. Therefore, as a cautionary remark, if we would like to rule out all single field inflation models in a rigorous model-independent fashion with a detection of scale-invariant local non-Gaussianity, we have to keep in mind the caveat that there may be single field models which only respect the consistency condition in a very squeezed region beyond the experimental reach.

Superpositions
Different shapes and runnings of non-Gaussianities can be superimposed in inflation models. For example, • Mixing shapes.
It is possible that different non-Gaussianity generation mechanisms are from different components in a model, or at different stages during inflation. So two or more different shapes can get mixed, and the final shape can be rather different. For example, in Figure 17: A mixing of the equilateral (Fig. 4) and local shape (Fig. 16). Fig. 17, we plot a mixed shape between the local and equilateral shape. Notice that this is different from the intermediate shapes, since obviously the squeezed limit is always dominated by the local form. Examples of such models are discussed in Ref. [200].
• Mixing shape and running.
The shapes can also be mixed with runnings. Same as the power spectrum, the non-Gaussianities generically have some mild scale dependence. But a more dramatic case is the superposition with a strong running, such as the sinusoidal or resonant running. For example, an inflaton passing through features frequently and turning constantly at the same time on a potential landscape can generate a bispectrum which is a superposition of the resonant running and intermediate shape, as we illustrate in Fig. 18. Clearly, these two signals are orthogonal to each other very well, and have to be picked up separately through different methods in data analyses.
If a non-Gaussianity is the linear superposition of several base components, one can generally perform a change of bases to make the new bases orthogonalized. For example, as we have seen in Sec. 6.1, the leading large bispectrum has two components, S λ and S c . The two shapes are very similar, and represented by the equilateral ansatz in data analyses. However since they do have small difference, one can subtract their similarities and get a new orthogonalized base component [173]. The orthogonalization is defined by the shape correlator such as (5.12). Using this definition, the new bases can be chosen as S 1 ≈ S λ + 0.22S c (9.7) and S 2 = S c . [Note that the S λ and S c used here do not include the prefactors (1/c 2 s − 1 − 2λ/Σ) and (1/c 2 s − 1) in (6.16) and (6.17).] Their shapes are shown in Fig. 19. Notice that S 1 is half positive and half negative. Because S 1 is not of the simplest factorizable type, the following simple ansatz has been proposed to represent S 1 in data analyses [173], We plot the shape of this ansatz in Fig. 20. The current CMB constraint on this orthogonal ansatz is −410 < f orth N L < 6 [1]. For known examples of general single field inflation, such as the DBI and k-inflation, we generically get equilateral shapes. This is also clear from their physical origin that Figure 19: Orthogonalization of two shapes in Sec. 6.1 (Fig. 3 and Fig. 4) through a change of bases, c λ S λ + c c S c = c 1 S 1 + c 2 S 2 . Figure 20: An ansatz −S orth ansatz in (9.8) for the orthogonal shape S 1 in (9.7). Note we added a minus sign in this plot. we have emphasized. The orthogonal shape relies on a delicate cancellation between the two generic shapes. In principle one can do this since the required parameter space is allowed in our effective field theory of general single field inflation in Sec. 6.1, and this may provide guidance to future model building. For example, one may finetune the parameters in the k-inflation models [93,94]. Therefore, unlike the previous cases, the direct motivation here is more oriented to data analyses. The advantage of this operation is that it makes full use of data, which impose constraints on both components. In addition, as a bonus, the ansatz for the equilateral (6.22), folded (6.51) and orthogonal (9.8) shapes are not linearly independent. As we can see, they all happen to be the equilateral ansatz shifted by a constant shape ansatz (S = const.) [63]. Constraining two orthogonal bases provide efficient constraints on all three of them.
Let us do a more data-analysis oriented exercise. We would like to construct an ansatz that is orthogonal to both local and equilateral ansatz, since both were well constrained by data. (Note that S orth ansatz in (9.8) is not quite orthogonal to the local ansatz, with a correlation ∼ −0.48). To do this we start with a trial shape S trial , and demand the new orthogonal shape S orth,2 ansatz = S trial + c 1 S loc + c 2 S eq ansatz (9.9) be orthogonal to both the local and equilateral ansatz, S orth,2 ansatz · S loc = 0 , S orth,2 ansatz · S eq ansatz = 0 . (9.10) The simplest factorizable trial shapes can be either the constant shape or the local-like shape k 1 /k 2 + 5 perm., and both give the same result. Let us use the constant shape S trial = 1 as an example. Solving the conditions (9.10) gives c 1 = −0.0953 and c 2 = −0.204. So the new orthogonal ansatz is S orth,2 ansatz = 1.19 k 2 1 k 2 k 3 + 2 perm. − 1.22 k 1 k 2 + 5 perm. + 3.44 . (9.11) The numerical details may change slightly depending on the detailed definition and computation of the inner product (5.11). The shape is plotted in Fig. 21. It is somewhat exotic but the ansatz is simple. By construction, this ansatz is much more orthogonal to the local form (with correlation ∼ 0) than the S orth ansatz currently used in Ref. [173,1]. It also happens to have reasonably large correlation (∼ 0.86) with the orthogonal shape in single field inflation (the S 1 shown in Fig. 19), similar to that (∼ −0.91) between S 1 and S orth ansatz . Obviously, other choices of trial shapes can result in more exotic orthogonal shapes.
One can perform a similar orthogonalization for the two shapes in (4.22), now they are both local to start with. More generally, if a non-Gaussianity has more base components, we can orthonormalize all of them one by one, in the sense of the Gram-Schmidt process.

Conclusion
The field of primordial non-Gaussianity is growing rapidly in recent years, with simultaneous progress from the experimental results, data analyses methods, non-linear cosmology theories, physical model buildings, computational techniques, and theoretical formalisms. The progress that we have seen so far is no doubt just a beginning.
In this review, we have studied the primordial non-Gaussianities coming from the inflation models, especially various mechanisms that can produce observable large non-Gaussianities with viable power spectra. We emphasized the fingerprints that different underlying physics leave on non-Gaussian profiles, which break the degeneracy of model building. We described the physical pictures and presented their effective Lagrangians to the extent that they can be recognized when encountered in the inflation model building in a more fundamental theory. We also derived the resulting bispectra and represented them in terms of simple ansatz to the extent that they can be useful to data analyses. With the current rapid progress, we anticipate much more future developments along these lines through refinements and discoveries in both theories and experiments.
The standard model of cosmology -the Big Bang theory with ΛCDM -is now established better than ever, with the precision data coming from the cosmic microwave background and large scale structures. New data will continue to flow from many ongoing and forthcoming experiments. Although Nature does not seem to be obligated to provide us any more information beyond the standard model, exciting possibilities exist that would help us to understand the origin of the Big Bang. These include the more detailed deviations from the scale-invariance of the power spectrum, the primordial gravitational waves that we may detect from the CMB polarization, the isocurvature perturbations between matter and radiation, and the primordial non-Gaussianities. Without these types of data, the number of theoretical models with degenerate observational consequences proliferate with time and it will be hard to understand the microscopic nature of the inflation beyond our current knowledge, as well as to distinguish inflation from other possible alternatives. As we have reviewed, primordial non-Gaussianity -the collider in the very early universe -is one of the few hopes. It is becoming a target of many modern experiments. We do not know which cards Nature is hiding from us, but we are hoping and preparing for the best.