Primordial Non-Gaussianity in the Cosmic Microwave Background

In the last few decades, advances in observational cosmology have given us a standard model of cosmology. We know the content of the universe to within a few percent. With more ambitious experiments on the way, we hope to move beyond the knowledge of what the universe is made of, to why the universe is the way it is. In this review paper we focus on primordial non-Gaussianity as a probe of the physics of the dynamics of the universe at the very earliest moments. We discuss 1) theoretical predictions from inflationary models and their observational consequences in the cosmic microwave background (CMB) anisotropies; 2) CMB--based estimators for constraining primordial non-Gaussianity with an emphasis on bispectrum templates; 3) current constraints on non-Gaussianity and what we can hope to achieve in the near future; and 4) non-primordial sources of non-Gaussianities in the CMB such as bispectrum due to second order effects, three way cross-correlation between primary-lensing-secondary CMB, and possible instrumental effects.

bispectrum based estimators to constrain primordial non-Gaussianity (f NL ). In Section V we discuss the current constraints on f NL by CMB bispectrum and what we can hope to achieve in near future. We also discuss non-primordial sources of non-Gaussianity which contaminate primordial bispectrum signal. In section VI we discuss other methods for constraining f NL besides CMB bispectrum. Finally in Section VII we summarize with concluding remarks.

II. INTRODUCTION: THE EARLY UNIVERSE
One of the most promising paradigms of the early universe is inflation [1,2,4], which apart from solving the flatness, homogeneity and isotropy problem, also gives a mechanism for producing the seed perturbations for structure formation, and other testable predictions 1 (for a recent review of inflationary cosmology see [43]). During inflation, the universe goes through an exponentially expanding phase. From the Friedman equation, the condition for the accelerated expansion is ρ + 3p < 0. (1) For both matter and radiation this condition is not satisfied. But it turns out that for a scalar field, the above condition can be achieved. For a spatially homogeneous scalar field, φ, moving in a potential, V (φ), the energy density is given by and the pressure is given by Hence the condition for accelerated expansion of the universe dominated with scalar field φ is Physically this condition corresponds to situations where kinetic energy of the field is much smaller than its potential energy. This condition is referred to slowly-rolling of the scalar field. During such slow-roll, the Hubble parameter, H(t) = d ln a/dt, is nearly constant in time, and the expansion scale factor, a(t) is given by This exponential expansion drives the observable universe spatially flat, homogeneous and isotropic. A toy model is shown in Fig. 1. In the slow-roll phase, φ rolls down on V (φ) slowly, satisfying Eq. (4) and hence driving the universe to expand exponentially. Near the minima of the potential, φ oscillates rapidly and inflation ends. After inflation ends, interactions of φ with other particles lead φ to decay with a decay rate of Γ φ , producing particles and radiation. This is called a reheating phase of the universe, as φ converts its energy density into heat by the particle production.
Not only inflation solves the flatness, homogeneity and isotropy problem, it also gives a mechanism for generating seed perturbations. During inflation the quantum fluctuation in the field φ are exponentially stretched due to the rapid expansion phase. The proper wavelength of the fluctuations are stretched out of the Hubble-horizon scale to that time, H −1 . Once outside the horizon, the characteristic r.m.s. amplitude of these fluctuations is |φ| rms ∼ H/(2π). These fluctuations do not change in time while outside the horizon. After inflation, and reheating, the standard hot-big scenario starts. As the universe decelerates, at some point the fluctuations re-enter the Hubble horizon, seeding matter and radiation fluctuations in the universe. Figure 2 summarizes the evolution of characteristic length scales. 1 Although inflation is the most popular theory for the early universe, other mechanisms, for example, ekpyrotic models [40] and cyclic models [41,42] have been proposed for generating nearly scale invariant Gaussian perturbations, while retaining homogeneity and flatness. In the cyclic universe, there is no beginning of time, and our expansion of the universe is one out of the infinite number of such cycles. Each cycle consists of the following phases: (1) A hot big bang phase, during which a structure formation takes place. (2) An accelerated expansion phase which dilutes the matter and radiation energy density. Since observations suggest that our universe is going through an accelerated expansion phase, in the cyclic model interpretation, we are presently going through this phase. (3) A decelerating phase, which makes the universe flat, and generates nearly Gaussian and scale invariant density perturbations. (4) A big crunch/bang transition phase during which matter and radiation is created. Although the mechanism is different, the outcome of phase (3) of the cyclic model is in some sense analogous to a slow-roll expansion phase of inflation; and phase (4) will correspond with the reheating phase in the inflationary scenario. As we will discuss in the next section these two scenarios can be distinguished by their different predictions about the gravitational waves, and non-Gaussianity. Cyclic models predict negligible contribution of gravitational waves while inflationary models can produce large gravitational wave contribution, which can be detected by next generation experiments. Second, cyclic models produce much larger non-Gaussianity (of local type) in comparison to the standard slow-roll inflationary scenario.
FIG. 1: A toy scenario for the dynamics of the scalar field during inflation. During the flat part of potential, universe expand exponentially. When field reaches near the minima of the potential, the field oscillates and the radiation is generated.

Primordial Perturbations
We use linearly perturbed conformal Friedmann Lematre Robertson Walker (FLRW) metric of the form, where all the metric perturbations, A, B, H L , and H T , are ≪ 1, and functions of conformal time τ . The spatial coordinate dependence of the perturbations is described by the scalar harmonic eigenfunctions, Q, Q i , and Q ij , that satisfy δ ij Q ,ij = −k 2 Q, Q i = −k −1 Q ,i , and Q ij = k −2 Q ,ij + 1 3 δ ij Q. Note that Q ij is traceless: δ ij Q ij = 0. Lets consider two new perturbation variables [8,44], and ζ ≡ − aḢ which are Gauge invariant. Here R ≡ H L + 1 3 H T , is perturbations in the intrinsic spatial curvature. While u reduces to δφ in the spatially flat gauge (R ≡ 0), or to −(φ/aH)R in the comoving gauge (δφ ≡ 0), its value is invariant under any gauge transformation. Similarly ζ, which reduces to R in the comoving gauge, and to −(aH/φ)δφ in the spatially flat gauge, is also gauge invariant. The perturbation variable ζ helps the perturbation analysis not only because of being gauge invariant, but also because it is conserved on super-horizon scales throughout the cosmic evolution.
The quantum fluctuations generate the gauge-invariant perturbation, u, that reduces to either δφ or (φ/aH)R depending on which gauge we use, either the spatially flat gauge or the comoving gauge. Hence, δφ flat and (φ/aH)R com are equivalent to each other at linear order. The benefit of using u is that it relates these two variables unambiguously, simplifying the transformation between δφ flat and R com .
The solution for ζ is valid throughout the cosmic history regardless of whether a scalar field, radiation, or matter dominates the universe; thus, once created and leaving the Hubble horizon during inflation, ζ remains constant in time throughout the subsequent cosmic evolution until reentering the horizon. The amplitude of ζ is fixed by the quantum-fluctuation amplitude in u This is the spectrum of ζ on super-horizon scales.
FIG. 2: Evolution of comoving horizon and generation of perturbations in the inflationary universe. Figure from Ref. [43].

From Primordial Perturbations to CMB Anisotropies
The metric perturbations perturb CMB, producing the CMB anisotropy on the sky. Among the metric perturbation variables, the curvature perturbations play a central role in producing the CMB anisotropy.
As we have shown in the previous subsection, the gauge-invariant perturbation, ζ, does not change in time on super-horizon scales throughout the cosmic evolution regardless of whether a scalar field, radiation, or matter dominates the universe. The intrinsic spatial curvature perturbation, R, however, does change when equation of state of the universe, w ≡ p/ρ, changes. Since ζ remains constant, it is useful to write the evolution of R in terms of ζ and w; however, R is not gauge invariant itself, but ζ is gauge invariant, so that the relation between R and ζ may look misleading. In 1980, Bardeen [45] introduced another gauge-invariant variable, Φ (or Φ H in the original notation), which reduces to R in the zero-shear gauge, or the Newtonian gauge, in which B ≡ 0 ≡ H T . Φ is given by Here, the terms in the parenthesis represent the shear, or the anisotropic expansion rate, of the τ = constant hypersurfaces. While Φ represents the curvature perturbations in the zero-shear gauge, it also represents the shear in the spatially flat gauge in which R ≡ 0. Using Φ, we may write ζ as where the terms in the parenthesis represent the gauge-invariant fluid velocity. We use Φ in rest of the paper because it gives the closest analogy to the Newtonian potential, which we have some intuition of. Φ reduces to R in the zero-shear gauge (or the Newtonian gauge) in which the metric (Eq.(6)) becomes just like the Newtonian limit of the general relativity.
The gauge-invariant velocity term, v − k −1Ḣ T , differentiates ζ from Φ. Since this velocity term depends on the equation of state of the universe, w = p/ρ, the velocity and Φ change as w changes, while ζ is independent of w. The evolution of Φ on super-horizon scales in cosmological linear perturbation theory gives the following [46], for adiabatic fluctuations, and hence Φ = 2 3 ζ in the radiation era (w = 1/3), and Φ = 3 5 ζ in the matter era (w = 0). Φ then perturbs CMB through the so-called (static) Sachs-Wolfe effect [47] ∆T At the decoupling epoch, the universe has already been in the matter era in which w = 0, so that we observe adiabatic temperature fluctuations of ∆T /T = − 1 3 Φ = − 1 5 ζ, and the CMB fluctuation spectrum of the Sachs-Wolfe effect, ∆ 2 SW (k), is By projecting the 3-dimension CMB fluctuation spectrum, ∆ 2 SW (k), on the sky, we obtain the angular power spectrum 2 , C l [48], where τ 0 and τ dec denote the conformal time at the present epoch and at the decoupling epoch, respectively, and n s ≡ 1 + d ln ∆ 2 (k)/d ln k is a spectral index which is conventionally used in the literature. On small angular scales (ℓ > 10), the Sachs-Wolfe approximation breaks down, and the acoustic physics in the photon-baryon fluid system modifies the primordial radiation spectrum [49]. To calculate the anisotropies at all the scales, one has to solve the Boltzmann photon transfer equation together with the Einstein equations. These equations can be solved numerically with the Boltzmann code such as CMBFAST [50]. The CMB power spectrum then can be written as Here g Tℓ (k) is called the radiation transfer function, and it contains all the physics which modifies the primordial power spectrum ∆ Φ to generate CMB power spectrum C ℓ . For the adiabatic initial conditions, in the Sachs-Wolfe limit, is called the dimensionless power spectrum. If Φ were exactly Gaussian, all the statistical properties of Φ would be encoded in the two-point function or in C ℓ in the spherical harmonic space. Since Φ is directly related to ζ through Eq. (12), all the information of ζ is also in-coded in C ℓ . Although ζ which is related to a Gaussian variable, u, through ζ = −(aH/φ)u, in the linear order ζ also obeys Gaussian statistics; however the non-linear relation between ζ and u makes ζ (and hence Φ and CMB anisotropies) slightly non-Gaussian. The non-linear relation between ζ and Φ is not the only source of non-Gaussianity in the CMB anisotropies. For example, at the second order, the relationship between Φ and ∆T /T is also non-linear.

Probes of the Cosmological Initial Conditions
The main predictions of a canonical inflation model are: • spatial flatness of the observable universe, • homogeneity and isotropy on large angular scales of the observable universe, • seed scalar and tensor perturbation with primordial density perturbations being (a) nearly scale invariant, (b) nearly adiabatic, and (c) very close to Gaussian.
At the time of writing, these predictions are consistent with all current observations. This represents a major success for the inflationary paradigm. On the other hand, the inflationary paradigm can be realized by a large 'zoo' 3 of models. In addition, 2 For the scale invariant (n = 1) case, C SW somewhat surprisingly, there exist scenarios where the Universe first contracts and then expands (such as the ekpyrotic/cyclic model), which (up to theoretical uncertainties regarding the precise mechanics of the bounce) also reproduce Universes with the properties described above. What we would like to do is to find observables that allows us to distinguish between members of the inflationary zoo. The exciting fact is that upcoming experiments will have the sensitivity to achieve this goal. Tilt and Running: Inflationary models very generically predict a slight deviation from completely flat spectrum. If we write the primordial power spectrum as ∆ Φ (k) = A(k 0 ) k k0 ns−1 , then n s = 1 correspond to flat spectrum and the quantity |n s − 1| is called a tilt, which characterizes the deviation from scale invariant spectrum. Although the deviations from the scale invariance are predicted to be small, the exact amount of deviation depends on the details of the inflationary model. For example in most slow roll models |n − 1| is of order 1/N e , where N e ∼ 60 is a number of e-folds to the end of inflation. Ghost inflation, however, predicts negligible tilt. Hence characterizing the tilt of the scalar spectral index is a useful probe of the early universe. Currently the most stringent constraints on tilt come from the WMAP 5-year data, n s = 0.960 +0.014 −0.013 [51], which already disfavors inflationary models with 'blue spectral index' (n s > 1). The 1σ error on n s will reduce to ∆n s = 0.0036 for upcoming Planck satellite and to ∆n s = 0.0016 for futuristic CMBPol like satellite [52].
Apart from the tilt in the primordial power spectrum, inflationary models also predict n s to be slightly scale dependent. This scale dependence is referred to as 'running' of the spectral index n s , and is defined as dn s /d ln k. The constraints on the running from the WMAP 5-year data are −0.090 < dn s /d ln k < 0.0019 [51]. The 1σ error will reduce to ∆(dn s /d ln k) = 0.0052 for upcoming Planck satellite and to ∆(dn s /d ln k) = 0.0036 for a fourth-generation satellite such as CMBPol [52].
Primordial Gravitational Waves: Inflation also generates tensor perturbations (gravitational waves), which although small compared to scalar component, are still detectable, in principle. So far primordial gravitational waves have not been detected. There are upper limits on their amplitude; see Ref. [53] for a current observational bounds on the level for primordial gravitational waves. Detection of these tensor perturbations or primordial gravitational waves is considered a 'smoking gun' for the inflationary scenario. In contrast to inflation, ekpyrotic (cyclic) models predict an amount of gravitational waves that is much smaller than polarized foreground emission would allow us to see even for an ideal CMB experiment. Primordial scalar perturbations create only E-modes of the CMB 4 , while primordial tensor perturbations generate both parity even E-modes and parity odd B-modes polarization [56][57][58]. The detection of primordial tensor B-modes in the CMB would confirm the existence of tensor perturbations in the early universe. This primordial B-mode signal is directly related to the Hubble parameter H during inflation, and thus a detection would establish the energy scale at which inflation happened. Various observational efforts are underway to detect such B-mode signal of the CMB [59]. Search for primordial B-modes is challenging. Apart from foreground subtraction challenges, and the challenge of reaching the instrumental sensitivity to detect primordial B-modes, there are several non-primordial sources such as weak lensing of CMB by the large scale structure [60,61], rotation of the CMB polarization [62][63][64][65], and instrumental systematics that generate B-modes which contaminate the inflationary signal [66][67][68]. The amplitude of gravitational waves is parametrized as the ratio of the amplitude of tensor and scalar perturbations, r. The limit from WMAP 5-year data is r < 0.22 (2σ) [51].
Isocurvature Modes: Inflationary models with a single scalar field predict primordial perturbations to be adiabatic. Hence detection of isocurvature density perturbations is a "smoking gun" for multi-field models of inflation. A large number of inflationary models with multiple scalar fields predict some amount of isocurvature modes [69][70][71][72][73][74][75][76][77][78][79][80][81]. For example, curvaton models predict the primordial perturbations to be a mixture of adiabatic and isocurvature perturbations. Isocurvature initial conditions specify perturbations in the energy densities of two (or more) species that add up to zero. It does not perturb the spatial curvature of comoving slice (i.e R is zero, hence the name isocurvature). In general, there can be four types of isocurvature modes, namely: baryon isocurvature modes, CDM isocurvature modes, neutrino density isocurvature modes and, neutrino velocity isocurvature modes. These perturbations imprint distinct signatures in the CMB temperature and E-polarization anisotropies [82]. The contribution of isocurvature modes is model dependent, and different models predict different amounts of it. There exists an upper limit on the allowed isocurvature modes using CMB temperature anisotropies [83,84], a characterization (or detection of any) of isocurvature modes has a potential of discriminating between early Universe models.
Primordial Non-Gaussianity: Canonical inflationary models predict primordials perturbations to be very close to Gaussian [5][6][7][8][9], and any non-Gaussianity predicted by the canonical inflation models is very small [14,15]. However models with nonlinearity [10,13,85], interacting scalar fields [12,86], and deviation from ground state [87,88] can generate large non-Gaussian perturbations. The amplitude of the non-Gaussian contribution to the perturbation is often referred to as f NL even if the nature of the non-Gaussianities can be quite different. Different models of inflation predict different amounts of f NL , starting from very close to zero for almost Gaussian perturbations, to f NL ≈ 100 for large non-Gaussian perturbations. For example, the canonical inflation models with slow roll inflation, where only a couple of derivatives of potential are responsible for inflationary dynamics, predict f NL ∼ 0.05 [15]. In models where higher order derivatives of the potential are important, the value of f NL FIG. 3: Shapes of Non-Gaussianity. The shape function F (k1, k2, k3) forms a triangle in Fourier space. The triangles are parametrized by rescaled Fourier modes, x2 = k2/k1 and x3 = k3/k1. Figure from Ref. [43] varies from f NL ∼ 0.1 where higher order derivatives are suppressed by a low UV cutoff [89] to f NL ∼ 100 based on Dirac-Born-Infeld effective action. Ghost inflation, where during inflation, the background has a constant rate of change as opposed to the constant background in conventional inflation, is also capable of giving f NL ∼ 100 [90]. The additional field models generating inhomogeneities in non-thermal species [91] can generate f NL ∼ 5 [92]; while curvaton models, where isocurvature perturbations in second field during the inflation generate adiabatic perturbations after the inflation, can have f NL ∼ 10 [93].
In the following we will see that non-Gaussianity, far from being merely a test of standard inflation, may reveal detailed information about the state and physics of the very early Universe, if it is present at the level suggested by the theoretical arguments above.

III. PRIMORDIAL NON-GAUSSIANITY
Large primordial non-Gaussianity can be generated if any of the following condition is violated [94] • Single Field. Only one scalar field is responsible for driving the inflation and the quantum fluctuations in the same field is responsible for generating the seed classical perturbations.
• Canonical Kinetic Energy. The kinetic energy of the field is such that the perturbations travel at the speed of light.
• Slow Roll. During inflation phase the field evolves much slowly than the Hubble time during inflation.
• Initial Vacuum State. The quantum field was in the Bunch-Davies vacuum state before the quantum fluctuation were generated.
To characterize the non-Gaussianity one has to consider the higher order moments beyond two-point function, which contains all the information for Gaussian perturbations. The 3-point function which is zero for Gaussian perturbations contains the information about non-Gaussianity. The 3-point correlation function of Bardeen's curvature perturbations, Φ(k), can be simplified using the translational symmetry to give where F (k 1 , k 2 , k 3 ) tells the shape of the bispectrum in momentum space while the amplitude of non-Gaussianity is captured dimensionless non-linearity parameter f NL . The shape function F (k 1 , k 2 , k 3 ) correlates fluctuations with three wave-vectors and form a triangle in Fourier space. Depending on the physical mechanism responsible for the bispectrum, the shape of the 3-point function, F (k 1 , k 2 , k 3 ) can be broadly classified into three classes [95]. The local, "squeezed," non-Gaussianity where F (k 1 , k 2 , k 3 ) is large for the configurations in which k 1 ≪ k 2 ≈ k 3 . Most of the studied inflationary and Ekpyrotic models produce non-Gaussianity of local shape (eg. [91,93,[96][97][98][99][100][101][102][103][104][105][106][107][108][109][110][111][112][113]). Second, the non-local, "equilateral," non-Gaussianity where F (k 1 , k 2 , k 3 ) is large for the configuration when k 1 ≈ k 2 ≈ k 3 . Finally the folded [114,115] Figure 3 shows these three shapes.
Non-Gaussianity of Local Type: The local form of non-Gaussianity may be parametrized in real space as 5 [13,116,117]: where ζ L (r) is the linear Gaussian part of the perturbations, and f NL characterizes the amplitude of primordial non-Gaussianity. Different inflationary models predict different amounts of f NL , starting from O(1) to f NL ∼ 100, beyond which values have been excluded by the Cosmic Microwave Background (CMB) bispectrum of WMAP temperature data. The bispectrum in this model can be written as where ∆ Φ is the amplitude of the primordial power spectrum. The local form arises from a non-linear relation between inflaton and curvature perturbations [10,11,13], curvaton models [93], or the New Ekpyrotic models [118,119]. Models with fluctuations in the reheating efficiency [9,10] and multi-field inflationary models [17] also generate non-Gaussianity of local type.
Being local in real space, non-Gaussianity of local type describes correlations among Fourier modes of very different k. In the limit in which one of the modes becomes of very long wavelength [120], k 3 → 0, (i.e. the other two k's become equal and opposite), ζ k3 freezes out much before k 1 and k 2 and behaves as a background for their evolution. In this limit F local is proportional to the power spectrum of the short and long wavelength modes As an example, for canonical single field slow-roll inflationary models, the three point function is given by [15] F slow-roll ( k 1 , k 2 , k3) Higher Deriv.  [95] where ǫ and η are the usual slow-roll parameters and are assumed to be much smaller than unity. Taking the limit k 3 → 0 gives the local form as in Eq. (20). To show this point, Figure. 4 compares the non-Gaussianity shape for local type and for slow-roll model. Although in this limit, slow-roll models do predict no-Gaussianity of local type but as evident from Eq. (21), the bispectrum of inflaton perturbations yields a non-trivial scale dependence of f NL [12,15]. However in the slow roll limit η, ǫ << 1 and hence the amplitude is too small to detect.
Non-Gaussianity of equilateral type: While vast number of inflationary models predict non-Gaussianity of local type, this model, for instance, fails completely when non-Gaussianity is localized in a specific range in k space, the case that is predicted from inflation models with higher derivative kinetic terms [90,115,[121][122][123][124]. In these models the correlation is among modes with comparable wavelengths which go out of the horizon nearly at the same time. The shape function for the equilateral shape can be written as [25] The models of this kind have large F (k 1 , k 2 , k 3 ) for the configurations where k 1 ≈ k 2 ≈ k 3 . The equilateral form arises from non-canonical kinetic terms such as the Dirac-Born-Infeld (DBI) action [122], the ghost condensation [90], or any other single-field models in which the scalar field acquires a low speed of sound [115,124].
As an example, models with higher derivative operators in the usual inflation scenario and a model of inflation based on the Dirac-Born-Infeld (DBI) action produce a bispectrum of the form The above model uses 1 8Λ 4 (∇φ) 2 (∇φ) 2 as a leading order operator. DBI inflation, which can produce large non-Gaussianity, f NL ∼ 100, also has F (k 1 , k 2 , k 3 ) of a similar form.
Ghost inflation, where an inflationary de Sitter phase is obtained with a ghost condensate, produces a bispectrum of the following form [90] where α and β are free parameters of order unity, and Ghost inflation also produces large non-Gaussianity, f NL ∼ 100. Figure III shows the shape of non-Gaussianity of equilateral type by showing F (k 1 , k 2 , k 3 ) for ghost inflation and for a model with a higher derivative term. Folded Shape: So far the 3-point functions were calculated assuming the regular Bunch-Davis vacuum state, giving rise to either local or equilateral type non-Gaussianity. However if the bispectrum is calculated by dropping the assumption of Bunch-Davis initial state give rise to bispectrum shape which peaks for the folded shape, k 1 ≈ 2k 2 ≈ 2k 3 , with shape function given as [114,115,125] where β kj are the Bogoliubov coefficients which encode information about the initial conditions, η 0 is the initial conformal time

IV. THE COSMIC MICROWAVE BACKGROUND BISPECTRUM
Since the discovery of CMB by Penzias and Wilson in 1965 [126] and the first detection of CMB temperature anisotropies on large scales by the COBE DMR [127], the space satellite WMAP and over a dozens of balloon and ground based experiments have characterized the CMB temperature anisotropies to a high accuracy and over a wide range of angular scales. The space satellite Planck which launched in 2009 will soon characterize the temperature anisotropies to even higher accuracy up to angular scales of ℓ max ≈ 2500. The CMB power spectrum is obtained by reducing all the information of N pix (∼ 10 6 for WMAP and ∼ 10 7 for Planck). Such reduction is justified to obtain a fiducial model, given the non-Gaussianities are expected to be small. With high quality data on the way, the field of non-Gaussianity is taking off. CMB bispectrum contains information which is not present in the power-spectrum and as we say in the previous section, is a unique probe of the early universe.
The harmonic coefficients of the CMB anisotropy a lm = T −1 d 2n ∆T (n)Y * ℓm can be related to the primordial fluctuation Φ as where Φ(k) is the primordial curvature perturbations, for a comoving wavevector k, g p ℓ (r) is the radiation transfer function where the index p refers to either temperature (T ) or E-polarization (E) of the CMB. A beam function b ℓ and the harmonic coefficient of noise n ℓm are instrumental effects. Eq. (27) is written for a flat background, but can easily be generalized.
Any non-Gaussianity present in the primordial perturbations Φ(k) gets transferred to the observed CMB via Eq. (27). The most common way to look for non-Gaussianity in the CMB is to study the bispectrum, the three-point function of temperature and polarization anisotropies in harmonic space. The CMB angular bispectrum is defined as and the angular-averaged bispectrum is where the matrix is the Wigner 3J symbol imposing selection rules which makes bispectrum zero unless Using Eq. (27) the bispectrum can be written as where Φ(k 1 )Φ(k 2 )Φ(k 3 ) is the primordial curvature three-point function as defined in Eq. (17).
To forecast constraints on non-Gaussianity using CMB data, we will perform a Fisher matrix analysis. The Fisher matrix for the parameters p a can be written as [21,34,117] The indices a and b run over all the parameters bispectrum depends on, we will assume all the cosmological parameters except f NL to be known. Indices ijk and pqr run over all the eight possible ordered combinations of temperature and polarization given by T T T , T T E, T ET , ET T , T EE, ET E, EET and EEE; the combinatorial factor ∆ ℓ1ℓ2ℓ3 equals 1 when all ℓ's are different, 6 when ℓ 1 = ℓ 2 = ℓ 3 , and 2 otherwise. The covariance matrix Cov is obtained in terms of C T T ℓ , C EE ℓ , and C T E ℓ (see [20,34]) by applying Wick's theorem.
For non-Gaussianity of the local type, for which the functional form F (k 1 , k 2 , k 3 ) is given by Eq. (19), we have where the functions α and β are given by In the expression above we use the dimensionless power spectrum amplitude ∆ Φ , which is defined by where n s is the tilt of the primordial power spectrum. One can compute the transfer functions g T ℓ (k) and g E ℓ (k) using publicly available codes such as CMBfast [50] and CAMB [128] In a similar way, from Eq. (22), one can derive the following expressions for the bispectrum derivatives in the equilateral case, where the functions δ, and γ are given by Recently a new bispectrum template shape, an orthogonal shape has been introduced [129] which characterizes the size of the signal (f ortho N L ) which peaks both for equilateral and flat-triangle configurations. The shape of non-Gaussianities associated with f ortho NL is orthogonal to the one associated to f equil N L . The bispectrum for orthogonal shape can be written as [129] ∂B ijk ℓ1ℓ2ℓ3 ∂f ortho

A. Estimator
An unbiased bispectrum-based minimum variance estimator for the nonlinearity parameter in the limit of full sky and homogeneous noise can be written as [17,21,25] where B ℓ1ℓ2ℓ3 is angle averaged theoretical CMB bispectrum for the model in consideration. The normalization N can be calculated to require the estimator to be unbiased, f NL = f NL . If the bispectrum B ℓ1ℓ2ℓ3 is calculated for f NL = 1 then the normalization takes the following form The estimator for non-Gaussianity, Eq. (39), can be simplified using Eq. (27) to yield where F (k 1 , k 2 , k 3 ) is a shape of 3-point function as defined in Eq. (17). Given the shape F (k 1 , k 2 , k 3 ), one is interested in, it is conceptually straightforward to constrain the non-linearity parameter from the CMB data. Unfortunately the computation time for the estimate scales as N

5/2
pix , which is computationally challenging as even for the WMAP data the number of pixels is of order N pix ∼ 10 6 . The scaling can be understood by noting that the each spherical harmonic transform scales as N 3/2 pix and the estimator requires ℓ 2 (∝ N pix ) number of spherical harmonic transforms.
The computational cost decreases if the shape can be factorized as with which the estimator simplifies tô and computational cost now scales as N 3/2 pix . For Planck (N pix ∼ 5 × 10 7 ) this translates into a speed-up by factors of millions, reducing the required computing time from thousands of years to just hours and thus making f NL estimation feasible for future surveys. The speed of the estimator now allows sufficient number of Monte Carlo simulations to characterize its statistical properties in the presence of real world issues such as instrumental effects, partial sky coverage, and foreground contamination. Using the Monte Carlo simulations it has been shown that estimator is indeed optimal, where optimality is defined by saturation of the Cramer Rao bound, if noise is homogeneous. Note that even for the non-factorizable shapes, by using the flat sky approximation and interpolating between the modes, one can estimating f NL in a computationally efficient way [130].
The extension of the estimator of f NL from the temperature data [17] to include both the temperature and polarization data of the CMB is discussed in Yadav et al. [20,21,34,35]. Summarizing briefly, we construct a cubic statistic as a combination of (appropriately filtered) temperature and polarization maps which is specifically sensitive to the primordial perturbations. This is done by reconstructing a map of primordial perturbations, and using that to define a fast estimator. We also show that this fast estimator is equivalent to the optimal estimator by demonstrating that the inverse of the covariance matrix for the optimal estimator [34] is the same as the product of inverses we get in the fast estimator. The estimator still takes only N 3/2 pix operations in comparison to the full bispectrum calculation which takes N 5/2 pix operations. For a given shape the estimator for non-linearity parameter can be written asf NL =Ŝ shape N shape , where for the equilateral, local and orthogonal shapes, the S shape can be written aŝ S local = 1 f sky r 2 dr d 2n B(n, r)B(n, r)A(n, r) S orthogonal = 9 f sky r 2 dr d 2n [B(n, r)B(n, r)A(n, r) + 8 9 D(n, r) 3 − 2B(n, r)C(n, r)D(n, r)] with B(n, r) ≡ ip lm and f sky is a fraction of sky. Index i and p can either be T or E.
Indices i, j, k, p, q and r can either be T or E. Here, ∆ ℓ1ℓ2ℓ3 is 1 when ℓ 1 = ℓ 2 = ℓ 3 , 6 when ℓ 1 = ℓ 2 = ℓ 3 , and 2 otherwise, B pqr,prim ℓ1ℓ2ℓ3 is the theoretical bispectrum for f NL = 1 [20]. It has been shown that the above estimators defined in Eq. (46) are minimum variance amongst bispectrum-based estimators for full sky coverage and homogeneous noise [20]. To be able to deal with the realistic data, the estimator has to be able to deal with the inhomogeneous noise and foreground masks. The estimator can be generalized to deal with partial sky coverage as well as inhomogeneous noise by adding a linear term toŜ prim :Ŝ prim →Ŝ prim +Ŝ linear prim . For the temperature only case, this has been done in [25]. Following the same argument, we find that the linear term for the combined analysis of CMB temperature and polarization data is given bŷ where A sim (n, r) and B sim (n, r) are the A and B maps generated from Monte Carlo simulations that contain signal and noise, and .. denotes the average over the Monte Carlo simulations. The generalized estimator is given byf Note that Ŝ linear prim MC = − Ŝ prim MC , and this relation also holds for the equilateral shape. Therefore, it is straightforward to find the generalized estimator for the equilateral shape: first, find the cubic estimator of the equilateral shape,Ŝ equil. , and take the Monte Carlo average, Ŝ equil. MC . Let us suppose thatŜ equil. contains terms in the form of ABC, where A, B, and C are some filtered maps. Use the Wick's theorem to re-write the average of a cubic product as ABC MC = A MC BC MC + B MC AC MC + C MC AB MC . Finally, remove the MC average from single maps, and replace maps in the product with the simulated maps A MC BC MC + B MC AC MC + C MC AB MC → A B sim C sim MC + B A sim C sim MC + C A sim B sim MC . This operation gives the correct expression for the linear term, both for the local form and the equilateral form.
The main contribution to the linear term comes from the inhomogeneous noise and sky cut. For the temperature only case, most of the contribution to the linear term comes from the inhomogeneous noise, and the partial sky coverage does not contribute much to the linear term. This is because the sky-cut induces a monopole contribution outside the mask. In the analysis one subtracts the monopole from outside the mask before measuringŜ prim , which makes the linear contribution from the mask small [25]. For a combined analysis of the temperature and polarization maps, however, the linear term does get a significant contribution from a partial sky coverage. Subtraction of the monopole outside of the mask is of no help for polarization, as the monopole does not exist in the polarization maps by definition. (The lowest relevant multipole for polarization is l = 2.) The estimator is still computationally efficient, taking only N 3/2 pix (times the r sampling, which is of order 100) operations in comparison to the full bispectrum calculation which takes N 5/2 pix operations. Here N pix refers to the total number of pixels. For Planck, N pix ∼ 5 × 10 7 , and so the full bispectrum analysis is not feasible while our analysis is.

A. Current Status
Currently the the Wilkinson Microwave Anisotropy Probe (WMAP) satellite provides the "best" (largest number of signal dominated modes) CMB data for non-Gaussianity analysis. Over the course of WMAP operation the field of non-Gaussianity has made vast progress both in terms of theoretical predictions of non-Gaussianities from inflation and improvement in the bispectrum based estimators. At the time of WMAP's first data release in 2003 the estimator was sub-optimal in the presence of partial sky coverage and/or inhomogeneous noise. With the sub-optimal estimator one could not use the entirety of WMAP data and only the data up to l max 350 were used to obtain the constraint f local NL = 38 ± 96(2σ) [23]. These limits were around 30 times better than the previous constraints of |f NL | < 1.5 × 10 3 from the Cosmic Background Explorer (COBE) satellite [? ].
By the time of second WMAP release in 2007 the estimator was generalized by adding a linear to the KSW estimator which allows to deal with partial sky coverage and inhomogeneous noise. The idea of adding a linear term to reduce excess variance due to noise inhomogeneity was introduced in [25]. Applied to a combination of the Q, V and W channels of the WMAP 3-year data up to l max ∼ 400 this estimator had yielded the tightest constraint at the time on f NL as: −36 < f NL < 100(2σ) [26]. This estimator was further generalized to utilize both the temperature and E-polarization information in [21], where it was pointed out that the linear term had been incorrectly implemented in Eq. 30 of [25]. Using Monte-Carlo simulations it has been shown that this corrected estimator is nearly optimal and enables analysis of the entire WMAP data without suffering from a blow-up in the variance at high multipoles 6 . The first analysis using this estimator shows an evidence of non-Gaussianity of local type at around 2.8σ in the WMAP 3-year data. Independent analysis shows the evidence of non-Gaussianity at lower significance, around 2.5σ (see Table VI).
By the time of the third WMAP data release (with 5-year obsevational data) in 2008 the f NL estimation technique was improved further by implementing the covariance matrix including inhomogeneoous noise to make the estimator completely optimal [131]. Using the optimal estimator and using the entirety of WMAP 3-year data there is an evidence for non-Gaussianity of local type at around 2.5σ level f NL ≈ 58 ± 23(1σ) [131]. However with WMAP 5-year data the significance goes down from ∼ 2.5σ to ∼ 1.8σ [131]. Table VI compares the constraints obtained by different groups using WMAP 3-year and WMAP 5-year data. Fig 6 shows this comparison in more detail, showing the constraints also as a function of maximum multipole l max used in the analysis. Few comments are in place: 1) constraints on f NL from WMAP 3-year data as a function of l max show a trend where the mean value rises at around l max = 450, below which data is consistent with Gaussianity and above which there is deviation from Gaussianity at above 2σ. The result becomes roughly independent of ℓ max > 550 with evidence for non-Gaussianity at around 2.5σ level, 2) independent analysis and using different estimators (optimal and near-optimal with linear term) sees this deviation from non-Gaussianity at around 2.5σ in WMAP 3-year data, 3) significance of non-Gaussianity goes down to around 2σ with WMAP 5-year data. The drop in the mean value between WMAP 3-year and 5-year data can be attributed to statistical shift.

B. Future Prospects
Now we discuss the future prospects of using the bispectrum estimators for constraining the non-linearity parameter f NL for local and equilateral shapes. We compute the Fisher bounds for three experimental setups, (1) cosmic variance limited experiment with perfect beam (ideal experiment hereafter), (2) Planck satellite with and noise sensitivity ∆ p = 56µK-arcmin and beam FWHM σ = 7 ′ (3) a futuristic CMBPol like satellite experiment with noise sensitivity ∆ p = 1.4µK-arcmin and beam FWHM σ = 4 ′ (CMBPol hereafter). Beside f NL we fix all the other cosmological parameters to a standard fiducial model with a flat ΛCDM cosmology, with parameters described by the best fit to WMAP 5-year results [51], given by Ω b = 0.045, Ω c = 0.23, H 0 = 70.5, n s = 0.96, n t = 0.0, and τ = 0.08. We calculate the theoretical CMB transfer functions and power spectrum from publicly available code CMBFAST [50]. We also neglect any non-Gaussianity which can be generated during recombination or there after. We discuss the importance and effect of these non-primordial non-Gaussianities in next section.
The scaling of signal-to-noise as a maximum multipole l max for the local [133,134] and equilateral model [135] are In principle one could go to arbitrary high l max but in reality secondary signals will certainly overwhelm primary signal beyond l max > 3000, we restrict to the analysis to l max = 3000. In Figure 7 we show the 1σ Fisher bound as function of maximum multipole ℓ max , for local and equilateral type of non-Gaussianity. For both local and equilateral case we show the Fisher bound for the analysis using only the CMB temperature information (TTT), only the CMB polarization information (EEE), and the combined temperature and polarization analysis. Note that by having both the temperature and E-polarization information one can improve the sensitivity by combining the information. Apart from combining the T and E signal, one can also do crosschecks and diagnostics by independently analysing the data. Temperature and polarization will have different foregrounds and instrumental systematics. Smith et al. [131] The difference between the results by Ref. [28] and [131] for WMAP 3-year data using the Kp0 mask can be a result of different choices of weighting in the near-optimal estimator. The optimal estimator has a unique weighting scheme.
A CMBPol like experiment will be able to achieve the sensitivity of ∆f local NL ≃ 2(1σ) for non-Gaussianity of local type and ∆f equil. NL ≃ 13(1σ) for non-Gaussianity of equilateral type. For the local type of non-Gaussianity this amounts to an improvement of about a factor of 2 over the Planck satellite and about a factor of 12 over current best constraints. These estimates assume that foreground cleaning can be done perfectly, i.e. the effect of residual foregrounds has been neglected. Also the contribution from unresolved point sources and secondary anisotropies such as ISW-lensing and SZ-lensing has been ignored.
Running non-Gaussianity: The primordial non-Gaussian parameter f NL has been shown to be scale-dependent in several models of inflation with a variable speed of sound, such as Dirac-Born-Infeld (DBI) models. Starting from a simple ansatz for a scale-dependent amplitude of the primordial curvature bispectrum for primordial non-Gaussianity, where K ≡ (k 1 k 2 k 3 ) 1/3 and k p is a pivot point. The primordial bispectrum is therefore determined in terms of two parameters: the amplitude f NL and the new parameter n N G quantifying its running. One can generalize the Fisher matrix analysis of the bispectra of the temperature and polarization of the CMB radiation and derive the expected constraints on the parameter n N G that quantifies the running of f NL (k) for current and future CMB missions such as WMAP, Planck and CMBPol. We will consider some non-zero f NL as our fiducial value for the Fisher matrix evaluation. Clearly, in order to be able to constrain a scale-dependence of f NL , its amplitude must be large enough to produce a detection. If f NL is too small to be detected (f NL < 2 is a lowest theoretical limit even for the ideal experiment), we will obviously not be able to measure any of its features, either. In the following we will then always consider a fiducial value of f NL large enough to enable a detection. Figure 8 shows the 1 − σ joint constraints on f NL and n N G . In the event of a significant detection of the non-Gaussian component, corresponding to f NL = 50 for the local model and f NL = 100 for the equilateral model of non-Gaussianity, is able to determine n N G with a 1 − σ uncertainty of ∆n N G ≃ 0.1 and ∆n N G ≃ 0.3, respectively, for the Planck mission and a factor of two better for CMBPol. In addition to CMB one can include the information of the galaxy power spectrum, galaxy bispectrum, and cluster number counts as a probe of non-Gaussianity on small scales to further constrain the two parameters [136].

C. Contaminations
A detection of non-Gaussianity has profound implications on our understanding of the early Universe. Hence it is important to know and quantify all the possible sources of non-Gaussianities in the CMB. Here we highlight some sources of non-Gaussianities due to second-order anisotropies after last scattering surface and during recombination. The fact that Gaussian initial conditions imply Gaussianity of the CMB is only true at linear order. We will also discuss the effects of instrumental effects and uncertainties in the cosmological parameters on the bispectrum estimate.

Secondary non-Gaussianities
Current analysis of the CMB data ignore the contributions from the secondary non-Gaussianities. For WMAP resolution it may not be a bad approximation. Studies of the dominant secondary anisotropies conclude that they are negligible for the analysis of the WMAP data for l max < 800 [117,137]. However on smaller angular scales several effects starts to kick in, for example, 1) the bispectrum contribution due to unresolved point source like thermal Sunyaev-Zeldovich clusters or standard radio sources, 2) three way correlations between primary CMB, lensed CMB and secondary anisotropies. We will refer to the bispectrum generated due to this three way correlation as B secondary−κ , where some secondaries are, the integrated Sachs-Wolfe (ISW) B ISW −κ , Sunyaev-Zeldovich signal and Rees-Sciama [23,117,[137][138][139][140].  [28]. The green triangles are using the the near-optimal estimator by Smith et al. [131]. The blue square results are obtained using either the optimal estimator by Smith et al. [131]. For all the three analysis Kp0 mask was used. Bottom panel: Comparison between 5-year results (optimal estimator, raw maps) reported in Komatsu et al. [51] and results obtained using the optimal or suboptimal estimator by Smith et al. [131].
For Future experiments such as Planck and CMBPOl the joint estimation of primordial and secondary bispectrum will be required. The observed bispectrum in general would take the following form The amplitude of bispectrum due to primary-lensing-secondary cross-correlation is proportional to the product of primary CMB power-spectrum and power spectrum of cross-correlation between secondary and lensing signal. The reduced bispectrum from the residual point sources (assuming Poisson distributed) is constant i.e. b ps ℓ1ℓ2ℓ3 = constant. The value of the constant will depend on the flux limit at which the point source can be detected and on assumed flux and frequency distribution of the sources.  Depending on the shape of primordial bispectrum in consideration, some secondary bispectra are more dangerous than others. For example, ISW-lensing B ISW −κ peaks at the "local" configurations, hence is more dangerous for local primordial shape than the equilateral primordial shape. For example for the Planck satellite if the secondary bispectrum is not incorporated in the analysis, the ISW-lensing contribution will bias the estimate for the local f NL by around ∆f local NL ≈ 10 [141]. The bispectrum contribution from primary-lensing-Rees-Sciama signal also peaks at squeezed limit and contribute to effective local f local NL ≈ 10 [142]. For Planck sensitivity the point source will contamination the local non-Gaussianity by around ∆f local NL ∼ 1 [143]. A recent analysis of the full second order Boltzmann equation for photons [144] claims that second order effects add a contamination ∆f NL ∼ 5.
The generalization of the Fisher matrix given by Eq. (31) to include multiple bispectrum contribution is where the additional X and Y indices denote a component such as primordial, point-sources, ISW-lensing etc. For fixed cosmological parameters, the signal-to-noise (S/N ) i for the component i is
The dominant bispectrum due to perturbative recombination comes from perturbations in the electron density. The amplitude of perturbations of the free electron density δ e is around a factor of 5 larger than the baryon density perturbations [165]. The bispectrum generated due to δ e peaks around the "local" configuration with corresponding effective non-linearity amplitude f NL ∼ few [166][167][168]. The bispectrum contribution due to second order terms which are the products of the first order perturbations is calculated in Ref. [169]. The bispectrum contribution from these terms which also peak for the squeezed triangles is small and can be neglected in the analysis. For example the signal-to-noise is about 0.4 at l max = 2000 for a full-sky, cosmic variance limited experiment.
Another contribution to bispectrum which peaks for the equilateral configurations comes from the non-linear evolution of the second order gravitational potential. Because of this effect the minimum detectable non-Gaussianity parameter f equil. NL changes by ∆f equil.

NL
= O(10) for Planck like experiment [135]. The bispectrum peaks for the equilateral shape because the the growth of potential happens on scales smaller than the horizon size.
On large scales, in the absence of primordial non-Gaussianities and assuming matter domination (so that the early and late ISW can be neglected), it has been shown in Ref. [162] that for the squeezed limit the effective f NL generated by second order gravitational effects on the CMB is f NL = −1/6 − cos(2θ) (also see [150][151][152]). Here θ is the angle between the short and the long modes. The angle dependent contribution comes from lensing.

Effect of Cosmological parameter uncertainties
Impact of uncertainties on the cosmological parameters effect the error bar on f NL . The effect of cosmological parameters have been discussed in Ref. [25,26,28,170]. The cosmological parameters are determined using the 2-point statistics of the CMB and therefor we expect the largest effect of f NL would come from those parameters which leave the CMB power spectrum unchanged while change the bispectrum. The expectation value of the estimator changes with the change in cosmological parameters. HereB ℓ1ℓ2ℓ3 is the true CMB bispectrum. When changing the parameters the normalization N should be changed to make the estimator unbiased. In general for a set of cosmological paramerets {p i }, the error in f NL is given by [170] δf Here the average parameter valuesp i and their covariance matrix Cov(p i , p j ) can be determined using CMB-likelihood analysis tools.
If the parameters are allowed to vary in the analysis then for WMAP this increases the 1σ uncertainty in f NL by δf local NL /f NL ≈ 16% for the local shape and δf equil NL /f NL ≈ 14% for the equilateral shape. For Planck experiment the increases in 1σ uncertainity is δf local NL /f NL ≈ 5% for local shape and δf equil.

NL
/f NL ≈ 4% for the equilateral shape. Most of the contribution to the error comes from three cosmological parameters, the amplitude of scalar perturbations ∆ Φ , the tilt of the power spectrum of the saclar perturbations n S , and re-ionization optical depth τ .
For modes inside the horizon during reionization, the reionization optical depth τ appears as a multiplicative factor e −τ in front of transfer function g i ℓ . For local model one of the mode is outside so the effect on bispectrumb local ℓ1ℓ2ℓ3 = exp(−2τ )b local ℓ1ℓ2ℓ3 and for equilateral model all the modes are inside the horizon sob equil. ℓ1ℓ2ℓ3 = exp(−3τ )b equil. ℓ1ℓ2ℓ3 . This reduces to δf local NL ≃ −2f NL τ for local model and δf equil.
The effect of amplitude of perturbations can be seen by noting that the level of non-Gaussianity is given by f NL · ∆ 1/2 Φ . Hence the decrease (increase) in the amplitude of perturbations relax (tighten) the constraints on f NL . The effect of red tilt (n s < 1) can be thought of as a reduction in power on at scales shorter than first peak and enhancement of power on larger scales. The effect of blue tilt is just opposite of red tilt. For local shape the limit on f NL becomes tighter proportional to ∆ 1/2 long [26]. Note that Ref. [170] show that the effect of cosmological parameters is negligible if the parameters are allowed to vary in the analysis and then marginalize over.

Instrumental effects and distortions along the line of sight
Here we point out that any cosmological or instrumental effect that can be modelled as a line of a sight CMB distortions of the primary CMB do not generate new bispectrum contribution. Although they can modify the the primordial bispectrum. A general model of line of sight distortions of the primary CMB are described in Ref. [66,68,171] where the changes in the Stokes parameter of the CMB due to distortions along the line-of-sight can be written as The first line captures the distortions in a single perfectly known directionn. The distortions in second line capture mixing of the polarization fields in a local region of length scale σ aroundn. We Taylor expand the CMB fields Q, U, and T around the pointn and consider the leading order terms. HereQ,Ũ , andT stands for primordial (un-distorted) CMB fields. Since (Q ± iU )(n) is spin ±2 field, a(n) is a scalar field that describes modulation in the amplitude of the fields in a given directionn; ω(n) is also a scalar field that describes the rotation of the plane of polarization, (f 1 ± if 2 ) are spin ±4 fields that describe the coupling between two spin states (spin-flip), and (γ 1 ± iγ 2 )(n) are spin ±2 fields that describe leakage from the temperature to polarization (monopole leakage hereon). Distortions in the second line of Eqn. (59), (p 1 ± p 2 ), (d 1 ± d 2 ), and q are measured in the units of the length scale σ. The field (p 1 ± ip 2 )(n) is a spin ±1 field and describes the change in the photon direction; we will refer to it as a deflection field. Finally (d 1 ± d 2 )(n) and q(n) describe leakage from temperature to polarization, (d 1 ± d 2 )(n) is spin ±1 field and we will refer to it as dipole leakage; q(n) is a scalar field that we will call quadrupole leakage.
These distortions can be produced by various cosmological processes such as weak gravitational lensing of the CMB, screening effects from patchy reionization, rotation of the plane of polarization due to magnetic fields or parity violating physics and various instrumental systematics such as gain fluctuations, pixel rotation, differential gain, pointing, differential ellipticity are also captured via line of sight distortions. All these distortions modify the primordial bispectrum as where W is a window which depends on the type of distortion in consideration and tells how the primordial CMB bispectrum modes are coupled to the distortion field power spectrum C DD ℓ . The effect of the distortions on the bispectrum is to smooth out the acoustic features. These effects for the case of lensing have been shown to be small and can be neglected [172,173].
In Ref. [174] the impact of the 1/f noise and asymmetric beam on local f local NL has been found insignificant in the context of a Planck-like experiment.

VI. OTHER PROBES OF NON-GAUSSIANITY IN THE CMB
Although using the full bispectrum is the most sensitive cubic statistic other statistical methods may be sensitive to different aspects of non-Gaussianity and, more importantly, different methods have different systematic effects. Therefore it is important to study various probes. In this section we will discuss some of the methods which have been recently used or developed to test for primordial non-Gaussianities in the CMB.

Trispectrum
The four-point function in harmonic space is called trispectrum, which can be written as where T l1l2 l3l4 (L) is the angular averaged trispectrum, L is the length of a diagonal that forms triangles with l 1 and l 2 and with l 3 and l 4 , and the matrix is the Wigner 3-j symbol. The trispectrum contains unconnected part, T G , which comes from the Gaussian part of the perturbations, and the connected part T c which contains non-Gaussian signatures.
Here, the matrix is the Wigner 6-j symbol, and t l1l2 l3l4 (L) is called the reduced trispectrum, which contains all the physical information about non-Gaussianities. For non-Gaussianity of local-type for which both f NL and g NL contribute to the trispectrum, but only f NL contributes to the bispectrum. Tripectrum based estimators for measuring f NL and g NL have been developed [175][176][177][178][179]. For local template, the bispectrum nearly contains all the information on f NL [177], however if the non-Gaussianity is seen in bispectrum, trispectrum can serve as a important cross-check.
Generically for single field slow-roll models the trispectrum is small and un-observable [180] however for more general single field models whenever the equilateral bispectrum is large, the trispectrum is large as well [181][182][183]. For example for equilateral non-Gaussianity Ref. [184] study how to tune the model parameters to get large trispectrum and small bispectrum. For multi-field inflation one can construct models that predicts small f NL but large g NL , for example Ref. [185] discusses the local form from a multi-field inflation, and briefly mentioned the condition in their class of models to get the large trispectrum and small bispectrum. Joint constraints on both f NL and g NL have the potential to add to the specificity of the search for primordial non-Gaussianity. For a given model these two numbers will often be predicted in terms of a single model parameter, such as a coupling constant, see e.g. [186] for the case of ekyprotic models. Using WMAP 5-year data, the constraints on g NL using the trispectrum are −7.4 < g NL /10 5 < 8.2 at 2σ [187].

Minkowski Functionals
Minkowski Functionals (MFs) describe morphological properties (such as area, circumference, Euler characteristic) of fluctuating fields [188][189][190][191]. For a d-dimensional fluctuating field, f , the k-th Minkowski Functionals of weakly non-Gaussian fields in, V (d) k (ν) for a given threshold ν = f /σ 0 can be written as [192,193] V (d) where σ 0 ≡ f 2 1/2 is the variance of the fluctuating field, H n (ν) are the Hermite polynomials, ω k ≡ π k/2 /Γ(k/2 + 1), and finally S (i) are the "skewness parameters" defined as which characterize the skewness of fluctuating fields and their derivatives. Here σ i characterizes the variance of the fluctuating field and is given by For CMB, for which d = 2 and f = ∆T /T , the skewness parameters are [194] where B l1l2l3 is the CMB bispectrum, W l represents a smoothing kernel which depends on the experiment beam and G m1m2m3 l1,l2,l3 is the usual Gaunt function.
Since MFs can be determined as weighted sum of the bispectrum, they contain less information than the bispectrum. MFs can still be useful because they perhaps suffer from different systematics, though they are less specific to primordial non-Gaussianity since they measure a smaller number of independent bispectrum modes. Also, the bispectrum is defined in Fourier (or harmonic) space while the MFs are defined in real space. Limits on non-Gaussianity of local-type from the MFs of the WMAP 5-year temperature data are −70 < f NL < 91(2σ) [195]. The MFs from the Planck temperature data should be sensitive to f NL ∼ 20 at 1σ level [196] in contrast to bispectrum which is sensitive to f NL ∼ 5 at 1σ level. Note that polarization data further improves the sensitivity.

Wavelets
Several studies have used wavelet representations of the WMAP maps to search for a non-Gaussian signal [197][198][199][200][201][202][203]. In most of these studies, wavelets were used as a tool for blind searches of non-Gaussian anomalies in a basis with resolution in both scale and location. However, in some more recent studies, wavelets were tuned to look for non-Gaussianity of a particular type. In the context of searches for primordial non-Gaussianity of local type, wavelet based estimators for f NL have been built by extracting a signature of local non-Gaussianity that is cubic in the wavelet coefficients from simulations of non-Gaussian skies and searching for this signature in data. This ability to calibrate on a set of simulations makes the wavelet approach very flexible. While not optimal in a least-squared sense, using a wavelet representation can be thought of as a generalized cubic statistic with a different weighting scheme to the optimal bispectrum estimator. Using such estimators therefore provides a useful exploration of nearly optimal cubic estimators similar to the full bispectrum estimator. Any believable detection of non-Gaussianity should be robust to such changes in the analysis. Similarly, contaminating non-Gaussianity from astrophysical and instrumental systematics will propagate through the analysis in a different way to the bispectrum-based analysis.
There are several constraints on local f NL using wavelet based estimators. For example, using the COBE data the constraints are |f NL | < 2200(1σ) [204]. Using an estimator based on the skewness of the wavelet coefficients, Mukherjee and Wang constrain the f NL value for WMAP 1-yr data obtaining f NL = 50 ± 160(2σ) [29]. Using an extension of the previous estimator by combining wavelet coefficients at different contiguous scales, Curto et al. obtain −8 < f NL < 111(2σ) [205]. Recently, using a generalized third order estimator based on the wavelet coefficients, Curto et al. obtain −18 < f NL < 80(2σ) [206].

Needlet Bispectrum
Needlets are a family of spherical wavelets which are localized and asymptotically uncorrelated [207,208]. The needlet based statistics as been considered for testing Gaussianity and isotropy (for example see Ref. [209][210][211][212][213][214][215]. Using the bispectrum of needlet coefficient, the constraints on non-Gaussianity of local-type using WMAP 5-year data yields f NL = 73 ± 62(2σ) [216,217]. As is clear, the needlet based bispectrum is not as sensitive as the CMB bispectrum discussed in Sec. IV, however again in the event of detection the needlet based methods can be calibrated on simulations and represent a different weighting scheme for handling the sky mask and anisotropic noise. Finally, needlets and wavelets allow for the possibility to analyze spatially localized regions in the sky.

Probing non-Gaussianity using Bayesian statistics
A somewhat different approach to searching for non-Gaussianity is provided by the Bayesian approach. Here, the starting point is an explicit physical or statistical model for the data and the goal is to evaluate the posterior density of the parameters of the model and/or the relative probability of the Gaussianity and non-Gaussianity.
On large scales, in the Sachs-Wolfe regime, one can simplify the Bayesian approach by modeling directly the temperature anisotropy. Rocha et al. (2001) [218] discuss a Bayesian exploration of a model where each spherical harmonic coefficient is drawn from a non-Gaussian distribution. In this regime, the simple form of the non-Gaussian potential for the local model Eq. (18) also translates into a simple model for the temperature anisotropy. Ref. [219] develop several results for it, including an analytical expression of the evidence ratio of the Gaussian and non-Gaussian models. At the level of current data, this approximation is too restrictive, since most of the information about f NL is contained near the resolution limit of the experiment, where most of the measured perturbation modes are concentrated.
A full implementation of a physical non-Gaussian model must include the effect of Boltzmann transport. In the context of local non-Gaussianity, the model equation Eq. (18) suggests that a full Bayesian treatment may be feasible. At the time of writing, no fully Bayesian analysis for local f NL has been published. The effort has focused on developing approximations to the full Bayesian problem.
Using a perturbative analysis, Ref. [220] relate the frequentist bispectrum estimator to moments of the Bayesian posterior distribution. Ref. [221] described approximations to the full Bayesian treatment that simplify the analysis for high signal-tonoise maps and compare these to the full Bayesian treatment for a simple 1-D toy model of non-Gaussian sky where this analysis is feasible.

VII. SUMMARY
The physics of the early universe responsible for generating the seed perturbations in the CMB is not understood. Inflation which is perhaps the the most promising paradigm for generating seed perturbations allow for vast number of inflationary models that are compatible with data based on 2-point statistics like CMB power spectrum. Moreover, the alternatives to inflation such as cyclic models are also compatible with the data. Characterizing the non-Gaussianity in the primordial perturbations has emerged as probe for discriminating between different models of the early universe. Models based on slowly rolling single field produce undetectable amount of non-Gaussianity. Single field models without the slow roll can generate large (detectable with future experiments) non-Gaussianities but 1) can not produce large non-Gaussianity of local type unless inflation started with with excited vacuum state 2) if non-Gaussianity is produced it would naturally be as bispectrum while higher order such as trispectrum can be generated, it requires fine tuning. The bispectrum of the CMB is one of the most promising tool for connecting the non-Gaussianities in the cosmic microwave background and the models of inflation. Bispectrum-based estimator which saturates Cramer-Rao bound has been developed and well characterised using non-Gaussian Monte-Carlos. Other statistics although not as sensitive to non-Gaussianity as an optimally weighted bispectrum estimator, do provide independent checks and have different systematics. While Bayesian analysis has been applied in the context of non-Gaussianity analysis, this still appears to be an open area for fruitful investigation.
Given the importance of detecting primordial non-Gussianity, it is crucial to characterise any non-primordial sources of non-Gaussianities. We describe several sources of non-Gaussianities such as from second order anisotropies after last scattering surface and during recombination.
With Planck launched and taking data, we look forward to the next few years as an exciting time in the exploration of primordial non-Gaussianity in the cosmic microwave background.