{'fraction_non_alphanumeric': 0.054709555992813755, 'fraction_numerical': 0.027696980066729403, 'mean_word_length': 4.124616396317404, 'pattern_counts': {'":': 0, '<': 17, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 18, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hunting for Primordial Non-Gaussianity in the Cosmic Microwave Background

Since the first limit on the (local) primordial non-Gaussianity parameter, fNL, was obtained from COBE data in 2002, observations of the CMB have been playing a central role in constraining the amplitudes of various forms of non-Gaussianity in primordial fluctuations. The current 68% limit from the 7-year WMAP data is fNL=32+/-21, and the Planck satellite is expected to reduce the uncertainty by a factor of four in a few years from now. If fNL>>1 is found by Planck with high statistical significance, all single-field models of inflation would be ruled out. Moreover, if the Planck satellite finds fNL=30, then it would be able to test a broad class of multi-field models using the four-point function (trispectrum) test of tauNL>=(6fNL/5)^2. In this article, we review the methods (optimal estimator), results (WMAP 7-year), and challenges (secondary anisotropy, second-order effect, and foreground) of measuring primordial non-Gaussianity from the CMB data, present a science case for the trispectrum, and conclude with future prospects.


Introduction
The physics of the very early, primordial universe is best probed by measurements of statistical properties of primordial fluctuations. The primordial fluctuations are the seeds for the temperature and polarization anisotropies of the CMB and the large-scale structure of the universe that we observe today. Therefore, both the CMB and the large-scale structure are excellent probes of the primordial fluctuations. In this article, we shall focus on the CMB. See the article by V. Desjacques and U. Seljak in this volume for the corresponding review on the large-scale structure as a probe of the primordial fluctuations.
This article reviews a recent progress on our using the CMB as a probe of a particular statistical aspect of primordial fluctuations called "non-Gaussianity." Reviews on this subject were written in 2001 [1] and 2004 [2]. The former review would be most useful for those who are new to this subject.
In this article, we focus on the new discoveries that have been made since 2004. Particularly notable ones include: 1. It has been proven that all inflation models (not just simple ones [3,4]) based upon a single scalar field would be ruled out regardless of the details of models [5], if the primordial non-Gaussianity parameter called f NL (more precisely, the "local type" f NL as described later) is found to be much greater than unity.
3. The most serious contamination of the local-type f NL due to the secondary CMB anisotropy, the coupling between the Integrated Sachs-Wolfe (ISW) effect and the weak gravitational lensing, has been identified [8,[12][13][14][15][16]. However, note that the astrophysical contamination such as the Galactic foreground emission and radio point sources may still be the most serious contaminant of f NL . These effects would pose a serious analysis challenge to measuring f NL from the Planck data.
4. The importance of distinguishing different triangle configurations of the three-point function of the CMB was realized [17,18] and has been fully appreciated. It has been shown by many researchers that different configurations probe distinctly different aspects of the physics of the primordial universe. The list of possibilities is long, and a terribly incomplete list of references on recent work (since ∼ 2004) is: [19] on a general analysis of various shapes; [20][21][22][23][24][25] on the local shape (k 3 ≪ k 1 ≈ k 2 ); [17,[26][27][28] on the equilateral shape (k 1 ≈ k 2 ≈ k 3 ); [29,30] on the flattend (or folded) shape (k 1 ≈ 2k 2 ≈ 2k 3 ); [31,32] on the orthogonal shape (which is nearly orthogonal to both local and equilateral shapes); [33][34][35][36][37] on combinations of different shapes; and [38,39] on oscillating bispectra. Also see references therein. 5. The connected four-point function of primordial fluctuations has been shown to be an equally powerful probe of the physics of the primordial universe. In particular, a combination of the three-and four-point functions may allow us to further distinguish different scenarios. Many papers have been written on this subject over the last few years: [40][41][42][43] on single-field models; [35, on multi-field models; and [70,71] on isocurvature perturbations. CMB data are expected to provide useful limits on the parameters of the "local-form trispectrum," τ NL and g NL [72,73]. Preliminary limits on these parameters have been obtained from the WMAP data by [74,75].
The number of researchers working on primordial non-Gaussianity has increased dramatically: Science White Paper on non-Gaussianity submitted to Decadal Survey Astro2010 was co-signed by 61 scientists [76].

What do we mean by "Gaussianity"?
What do we mean by "Gaussian fluctuations"? Let us consider the distribution of temperature anisotropy of the CMB that we observe on the sky, ∆T (n). The temperature anisotropy is Gaussian when its probability density function (PDF) is given by where ∆T i ≡ ∆T (n), ξ ij ≡ ∆T i ∆T j is the covariance matrix (or the two-point correlation function) of the temperature anisotropy, |ξ| is the determinant of the covariance matrix, and N pix is the number of pixels on the sky. We often work in harmonic space by expanding ∆T using spherical harmonics: ∆T (n) = lm a lm Y lm (n). The PDF for a lm is given by where C lm,l ′ m ′ ≡ a * lm a l ′ m ′ , and N harm is the number of l and m. When a lm is statistically homogeneous and isotropic (which is not always the case because of, e.g., non-uniform noise), one finds C lm,l ′ m ′ = C l δ ll ′ δ mm ′ , and thus the PDF simplifies to Here, C l is the angular power spectrum. The latest determination of C l of the CMB temperature anisotropy is shown in Figure 1.
The important property of a Gaussian distribution is that the PDF is fully specified by the covariance matrix. In other words, the covariance matrix contains all the information on statistical properties of Gaussian fluctuations. When the PDF is given by equation 3, the power spectrum, C l , contains all the information on a lm . This is not true for non-Gaussian fluctuations, for which one needs information on higher-order correlation functions.
Let us close this subsection by noting that a non-zero deviation of the covariance matrix from the diagonal form, ∆C lm,l ′ m ′ = C lm,l ′ m ′ − C l δ ll ′ δ mm ′ , does not imply non-Gaussianity: the PDF can be a Gaussian with a non-diagonal covariance matrix as given in equation 2. A non-zero ∆C lm,l ′ m ′ may arise in cosmological models that violate statistical isotropy. Such models may yield anisotropic Gaussian fluctuations; thus, one must distinguish between non-Gaussianity and a violation of statistical isotropy.

What do we mean by "non-Gaussianity"?
What do we mean by "non-Gaussian fluctuations"? Any deviation from a Gaussian distribution (such as equation 1 or 2) is called non-Gaussianity. When fluctuations in the CMB are non-Gaussian, one cannot generally write down its PDF, unless one considers certain models (e.g., inflation). Nevertheless, when non-Gaussianity is weak, one may expand the PDF around a Gaussian distribution [80] and obtain Figure 1. The angular power spectrum of the CMB temperature anisotropy, C l , measured from the WMAP 7-year data [77], along with the temperature power spectra from the ACBAR [78] and QUaD [79] experiments. The solid line shows the best-fitting 6-parameter flat ΛCDM model to the WMAP data alone. The angular power spectrum contains all the information on fluctuations in the CMB, if fluctuations are Gaussian. If fluctuations are non-Gaussian, one must use the higher-order correlation functions (such as three-and four-point functions) to fully exploit the cosmological information contained in the CMB. This figure is adopted from [11].
Here, the expansion is truncated at the three-point function (bispectrum) of a lm , and thus we have assumed that the connected four-point and higher-order correlation functions are negligible compared to the power spectrum and bispectrum. (This condition is not always satisfied.) By evaluating the above derivatives, one obtains ‡ This formula is useful, as it tells us how to estimate the angular bispectrum, a l 1 m 1 a l 2 m 2 a l 3 m 3 , optimally from a given data by maximizing this PDF. In practice, we usually parametrize the bispectrum using a few parameters (e.g., f NL ), and estimate those parameters from the data by maximizing the PDF with respect to the parameters. ‡ Babich [81] derived this formula for C lm,l ′ m ′ = C l δ ll ′ δ mm ′ .
In the limit that the contribution of the connected four-point function (trispectrum) to the PDF is negligible compared to those of the power spectrum and bispectrum, equation 5 contains all the information on non-Gaussian fluctuations characterized by the covariance matrix, C l 1 m 1 ,l 2 m 2 = a * l 1 m 1 a l 2 m 2 , and the angular bispectrum, a l 1 m 1 a l 2 m 2 a l 3 m 3 . This approach can be extended straightforwardly to the trispectrum if necessary.
3. Extracting f NL from the CMB data

General formula
We have not defined what we mean by "f NL ." For the moment, let us loosely define it as the amplitude of a certain shape of the angular bispectrum: where the function, b (i) l 1 l 2 l 3 , is called the "reduced angular bispectrum" [82] and defines the shape of the angular bispectrum for a given model denoted by an index i (which may refer to, e.g., "local," "equilateral," "orthogonal," etc), and G m 1 m 2 m 3 The physical role of the Gaunt integral is to assure that (l 1 , m 1 ), (l 2 , m 2 ), and (l 3 , m 3 ) form a triangle. In the small-angle limit, the Gaunt integral becomes a 2-d delta function: G m 1 m 2 m 3 l 1 l 2 l 3 → (2π) 2 δ D (l 1 + l 2 + l 3 ) [83]. Given this parametrization, one can then maximize the PDF given in equation 5 with respect to f One can use this formula to determine multiple amplitudes of angular bispectra simultaneously.
Here, S i are given by the data as where 1/6 is included such that F ij in equation 8 becomes the Fisher matrix of f (i) NL . In other words, the covariance matrix of f (i) NL is given by the inverse of F ij , i.e., The 68% uncertainty in f where Here, the summation over m can be done using the Fast Fourier Transform (FFT), as Y lm (θ, φ) ∝ e imφ . This technique is used by the HEALPix package [84], and thus one may use HEALPix to do this summation. To compute d ll ′ (n), one may use Monte Carlo simulations. Namely, as C lm,l ′ m ′ = a * lm a l ′ m ′ , we have the exact relation between d ll ′ and e l : d ll ′ (n) = e l (n)e l ′ (n) . One can evaluate the ensemble average using the Monte Carlo simulation of the CMB and the instrumental noise. Let us denote this operation by d ll ′ (n) = e l (n)e l ′ (n) MC . The final formula for S i is which is valid for general forms of b (i) l 1 l 2 l 3 . Note that the integral, d 2n , must be done over the full sky, even in the presence of the mask: the information on the mask is included in the calculation of the Fisher matrix, F ij . The only assumptions that we have made so far are: (1) each angular bispectrum component has only one free parameter, i.e., the amplitude, and (2) non-Gaussianity (if any) is weak, and the PDF of a lm is given by equation 5.
Finally, the explicit form of the Fisher matrix is given by where f sky is the fraction of the sky outside of the mask. When the covariance matrix is diagonal, the expression simplifies to where Equation 16 may also be written as where ∆ l 1 l 2 l 3 = 1, 2, and 6 when all of l i 's are different, two of l ′ i are the same, and all of l i 's are the same, respectively.

Poisson bispectrum
As a warm up, let us consider the simplest example: point sources randomly distributed over the sky. As mentioned already, this is a contamination of the primordial non-Gaussianity parameters, and thus an accurate measurement of this component is quite important, especially for Planck as well as for the South Pole Telescope (SPT) and Atacama Cosmology Telescope (ACT), which are working at high frequencies (ν > 100 GHz) where star-forming galaxies dominate a lm at l > 1000.
For the Poisson distribution, the reduced bispectrum is independent of multipoles, b src l 1 l 2 l 3 = 1, in the absence of window functions, and is given by b src in the presence of window functions. Here, w l is an experimental window function (a product of the beam transfer function and the pixel window function). Let us then use b src (instead of f NL because this component has nothing to do with primordial fluctuations) to denote the amplitude of the Poisson bispectrum.
From the data, we measure S src given by where a map E(n) is defined by [85] where e l (n) is given by equation (12). The E(n) map is a Wiener-filtered map of point sources randomly distributed on the sky.
The observed angular bispectrum is related to the 3-dimensional bispectrum of primordial curvature perturbations, In the linear order, the primordial curvature perturbation is related to Bardeen's curvature perturbation [103] in the matter-dominated era, Φ, by ζ = 5 3 Φ [104]. The CMB temperature anisotropy in the Sachs-Wolfe limit [105] is given by ∆T /T = − 1 3 Φ = − 1 5 ζ. We write the bispectrum of Φ as There is a useful way of visualizing the shape dependence of the bispectrum. We can study the structure of the bispectrum by plotting the magnitude of F (k 1 , k 2 , k 3 )(k 2 /k 1 ) 2 (k 3 /k 1 ) 2 as a function of k 2 /k 1 and k 3 /k 1 for a given k 1 , with a condition that k 1 ≥ k 2 ≥ k 3 is satisfied. In order to classify various shapes of the triangles, let us use the following names: Figure. 2 for the visual representations of these triangles.
We shall explore 3 different shapes of the primordial bispectrum: "local," "equilateral," and "orthogonal." They are defined as follows: 1. Local form. The local form bispectrum is given by [82,106,107] where P Φ = A/k 4−ns is the power spectrum of Φ with a normalization factor A. This form is called the local form, as this bispectrum can arise from the curvature perturbation in the form of Φ = Φ L + f local NL Φ 2 L , where both sides are evaluated at the same location in space (Φ L is a linear Gaussian fluctuation). § The local form, F local (k 1 , k 2 , k 3 )(k 2 /k 1 ) 2 (k 3 /k 1 ) 2 , peaks at the so-called "squeezed" triangle § However, Φ = Φ L + f local NL Φ 2 L is not the only way to produce this type of bispectrum. One can also produce this form from multi-scalar field inflation models where scalar field fluctuations are nearly as a function of k 2 /k 1 and k 3 /k 1 for a given k 1 , with a condition that k 3 ≤ k 2 ≤ k 1 is satisfied. As the primordial bispectra shown here are (nearly) scale invariant, the shapes look similar regardless of the values of k 1 . The amplitude is normalized such that it is unity at the point where for which k 3 ≪ k 2 ≈ k 1 [18]. See the top-left panel of Figure 3. In this limit, we obtain How large is f local NL from inflation? The earlier calculations showed that f local NL from single-field slow-roll inflation would be of order the slow-roll parameter, ǫ ∼ 10 −2 [106,114,115]. More recently, Maldacena [3] and Acquaviva et al. [4] found that the coefficient of P Φ (k 1 )P Φ (k 3 ) from the simplest single-field slow-roll inflation with the canonical kinetic term in the squeezed limit is given by Comparing this result with the form predicted by the f local NL model, one obtains f local NL = (5/12)(1 − n s ), which gives f local NL = 0.015 for n s = 0.963.
2. Equilateral form. The equilateral form bispectrum is given by [7] F equil (k 1 , k 2 , k 3 ) This function approximates the bispectrum forms that arise from a class of inflation models in which scalar fields have non-canonical kinetic terms. One example is the so-called Dirac-Born-Infeld (DBI) inflation [26,116], which gives f equil NL ∝ −1/c 2 s in the limit of c s ≪ 1, where c s is the effective sound speed at which scalar field fluctuations propagate relative to the speed of light. There are various other models that can produce f equil NL [19,27,28,117,118]. The equilateral form, Figure 3. The local and equilateral forms are nearly orthogonal to each other, which means that both can be measured nearly independently.
3. Orthogonal form. The orthogonal form, which is constructed such that it is nearly orthogonal to both the local and equilateral forms, is given by [31] F orthog (k 1 , k 2 , k 3 ) This form approximates the forms that arise from a linear combination of higher-derivative scalar-field interaction terms, each of which yields forms similar to the equilateral shape. Senatore, Smith and Zaldarriaga [31] found that, using the "effective field theory of inflation" approach [117], a certain linear combination of similarly equilateral shapes can yield a distinct shape which is orthogonal to both the local and equilateral forms. The orthogonal form, Note that these are not the most general forms one can write down, and there are other forms which would probe different aspects of the physics of inflation [19,29,33,37,69,119].
Of these forms, the local form bispectrum has special significance. Creminelli and Zaldarriaga [5] showed that not only models with the canonical kinetic term, but all single-inflation models predict the bispectrum in the squeezed limit given by equation 25, regardless of the form of potential, kinetic term, slow-roll, or initial vacuum state. Also see [19,28,117]. This means that a convincing detection of f local NL would rule out all single-field inflation models.

Optimal estimator for f local NL
Given the form of Φ, one can calculate the harmonic coefficients of temperature and E-mode polarization anisotropies as where g T l (k) and g P l (k) are the radiation transfer functions of the temperature and polarization anisotropies, respectively, which can be calculated by solving the linearized Boltzmann equations. One may use the publicly-available Boltzmann codes such as CMBFAST [120] or CAMB [121] for computing the radiation transfer functions. From now on, we shall focus on the temperature anisotropy, largely for simplicity. (See [122,123] for the treatment of polarization in the angular bispectrum.) The limits on f NL expected from Planck are dominated by the temperature information, and thus the polarization information is not expected to yield competitive limits over the next, say, > 5 years.
For the local-form bispectrum given in equation 23, the reduced bispectrum (equation 6) is given by [82] b local where Using this form in equation 14, one finds S local as which can be measured from the data. Here, maps A(n, r) and B(n, r) are defined by [6] A(n, r) ≡ where e l (n) is given by equation (12).
A CMBFAST-based code for computing g T l (k) and g P l (k) is available at http://gyudon.as.utexas.edu/∼komatsu/CRL. A recent version of CAMB has an option to calculate these functions (http://camb.info).

Optimal estimator for f equil NL
For the equilateral-form bispectrum given in equation 26, the reduced bispectrum is given by [7] b equil where γ l (r) = 2 π k 2 dk P which can be measured from the data. Here, maps C(n, r) and D(n, r) are defined by [7] C(n, r) ≡ l w l γ l (r)e l (n) = lm w l γ l (r)(C −1 a) lm Y lm (n), 3.6. Optimal estimator for f orthog NL For the equilateral-form bispectrum given in equation 27, the reduced bispectrum is given by [31] b orthog l 1 l 2 l 3 = 3b equil l 1 l 2 l 3 − 12 r 2 dr δ l 1 (r)δ l 2 (r)δ l 3 (r)w l 1 w l 2 w l 3 .
Using this form in equation 14, one finds S orthog as which can be measured from the data.

General formula for the lensing-secondary coupling
Given the special importance of the local-form bispectrum, we must understand what other (non-primordial) effects might also produce the local form, potentially preventing us from measuring f local NL .
The local-form bispectrum is generated when the power spectrum of shortwavelength fluctuations is modulated by long-wavelength fluctuations; thus, a mechanism that couples small scales to large scales can potentially generate the localform bispectrum.
The weak gravitational lensing provides one such mechanism. The local-form bispectrum may then be generated when long-and short-wavelength fluctuations are coupled by the lensing. To see how this might happen, let us write the observed temperature anisotropy in terms of the original (unlensed) contribution from the last scattering surface at z = 1090, ∆T P (where "P " stands for "primary"), the lensing potential, φ, and the secondary anisotropy generated between z = 1090 and z = 0, ∆T S (where "S" stands for "secondary"): where and r * is the comoving distance out to z = 1090, and Φ is Bardeen's curvature perturbation, which is related to the usual Newtonian gravitational potential by Φ = −Φ Newton . Transforming this into harmonic space and computing the reduced bispectrum, one obtains [12] b lens−S l 1 l 2 l 3 = l 1 (l 1 + 1) − l 2 (l 2 + 1) + l 3 (l 3 + 1) 2 where C P l is the power spectrum of the CMB from the decoupling epoch only (i.e., no ISW), and C φS l ≡ φ * lm a S lm is the lensing-secondary cross-correlation power spectrum. From this result, one finds that a non-zero bispectrum is generated when the secondary anisotropy traces the large-scale structure (i.e., Φ). Various secondary effects have been studied in the literature: the Sunyaev-Zel'dovich effect [12], cosmic reionization [124], point sources [125], and ISW [12]. It has been shown that the last one, the ISW-lensing coupling, is the most dominant contamination of f local NL [14]. Using equation 46 in equation 14, one finds S lens−S as which can be measured from the data. Here, the "linear terms" contain 9 terms with MC , such as −P (n) [∂ 2 E](n)Q(n) MC , etc. The map E(n) is given by equation 21, and the other maps are defined by The maps with ∂ 2 are given by ∂ 2 P = − l l(l + 1)w l C P l e l (n), etc. The map P (n) is a Wiener-filtered map of the primary temperature anisotropy from z = 1090.

Lensing-ISW coupling
A change in the curvature perturbation yields a secondary temperature anisotropy via the ISW effect [105]: where r is the comoving distance and r * is the comoving distance out to z = 1090. Here, note again Φ = −Φ Newton . The cross-power spectrum of φ and the ISW effect is then given by where P ΦΦ ′ (k, r) is the cross-power spectrum of Φ and Φ ′ ≡ ∂Φ/∂r, which can be calculated from the power spectrum of Φ, P Φ (k, r), as P ΦΦ ′ (k, r) = 1 2 [∂P Φ (k, r)/∂r] [13,126]. Here, P Φ (k, r) is not the primordial power spectrum, but it includes the linear transfer function, T (k), and the growth factor of Φ, g(r): Using this, one finds P ΦΦ ′ (k, r) = (g ′ /g)P Φ (k, r). Note that g(r) is normalized such that g(r) = 1 during the matter-dominated era. With this result, it is easy to see why the lensing-ISW coupling yields the squeezed configuration: on very large scales, where T (k) → 1, C φ,ISW l ∝ 1/l 3 . On smaller scales, T (k) declines with k, and thus C φ,ISW l falls faster than 1/l 3 . The lensing coupling includes l(l + 1)C φ,ISW l , which falls faster than 1/l, i.e., the largest power comes from the smallest l.
A recent estimate by Hanson et al. [15] showed that the lensing-ISW coupling, if not included in the parameter estimation, would bias f local NL by ∆f local NL = 9.3. The expected bias for WMAP is ∆f local NL = 2.7 [11]. One can remove this bias by including the lensing-ISW coupling (or any other lensing-secondary couplings) using the optimal estimator given by equation 48.

General discussion
So far, we have assumed that one can use equation 28: to convert the primordial curvature perturbation to the temperature anisotropy. (Here, the subscript "p" stands for "primordial," by which we mean Φ p = 3 5 ζ without the linear transfer function.) However, this equation is valid only for linear theory. As any non-linear effects can produce non-Gaussianity, one has to study the impacts of various non-linear effects on the observed non-Gaussianity. The origin of the linear radiation transfer function is the linearized Boltzmann equation: where η is the conformal time, µ ≡k ·n, ∆ (1) ≡ 4[∆T (1) (k, µ, η)/T ] is the perturbation in the photon energy density, and S (1) is the linear source function, which depends on the metric perturbations as well as on the density, velocity, pressure, and stress perturbations of matter and radiation in the universe and the photon polarization. The second-order Boltzmann equation is then similarly written as where ∆ (2) ≡ 8[∆T (2) (k,n, η)/T ] + 12[∆T (1) (k, µ, η)/T ] 2 , and S (2) is the second-order source function. Note that the azimuthal symmetry is lost at the second order, and thus the perturbations depend on the directions of k andn independently. In this case, the second-order a lm is given by [127] a (2) where F l ′ m ′ lm is the second-order radiation transfer function, whose form is determined by the second-order source function, S (2) , in the Boltzmann equation.
The shape of the second-order bispectrum, a l 3 m 3 , is determined by the shape of the second-order radiation transfer function. If the second-order radiation transfer function vanishes in the squeezed limit, i.e., F l ′ m ′ lm (k ′ , k ′′ , k) → 0 for k → 0, then the CMB bispectrum would not peak at the squeezed configuration, and thus the resulting f local NL would be small. The second-order source function is quite complicated [128][129][130][131][132][133][134][135][136][137][138], but it can be divided into two parts ¶: (1) The terms given by the products of the first-order perturbations, such as [Φ (1) ] 2 .
The intrinsically second-order terms are sourced by products of the first-order perturbations, and thus it is created by the late-time evolution of cosmological perturbations, whereas the terms in (1) are set by the initial conditions. ¶ This decomposition is not gauge invariant, and thus which terms belong to (1) or (2) depends on the gauge that one chooses. Therefore, one must specify the gauge when making such a decomposition. Our discussion in this section is based on the gauge choice made by Bartolo, Matarrese and Riotto [128,129] and Pitrou, Uzan and Bernardeau [136,137], which reduces to the Newtonian gauge at the linear order. This seems a convenient gauge, as the products of the first-order terms only give |f local NL | < 1 [127].
The contamination of f local NL due to the terms in (1) is small, |f local NL | < 1 [127]. Recently, Pitrou, Bernardeau and Uzan [137] reported a surprising result that the terms in (2) would give f local NL ∼ 5 for the Planck data (l max = 2000). Why surprising? As the intrinsically second-order terms arise as a consequence of the late-time evolution of the cosmological perturbations, they are generated by the causal mechanism, i.e., gravity and hydrodynamics. It is difficult for the causal mechanism to generate the bispectrum in the squeezed configuration, as it requires very long wavelength perturbations to be coupled to short wavelength ones.

Newtonian calculation
As an example, let us consider the well-known second-order solution for Φ (2) in the sub-horizon limit, i.e., k ≫ aH, which is equivalent to taking the non-relativistic (Newtonian) limit. Here, the second-order Bardeen curvature perturbation is defined by Φ = Φ (1) + 1 2 Φ (2) . The explicit solution is [136] where the linear perturbation, Φ (1) , on the right hand side is constant during the matter dominated era, and the symmetrized function, F (s) 2 , is defined as Note that F (s) 2 is related to the function G given in equation (8.9) of [129] as The function F (s) 2 (k 1 , k 2 ) vanishes in the squeezed limit, k 1 = −k 2 , and thus the CMB bispectrum generated from Φ (2) in the Newtonian limit is not given by the local form. To see this, let us calculate and P Φ (k) = AT 2 (k)/k 4−ns . The shape dependence of F 2nd (k 1 , k 2 , k 3 )(k 2 /k 1 ) 2 (k 3 /k 1 ) 2 is shown in Figure 4 for various values of k 1 (because F 2nd is not scale invariant). The CMB data are sensitive to k 1 < k max ∼ 0.2 h Mpc −1 (l max /2000). We find that the bispectrum peaks at the equilateral configuration on large scales (k 1 < 10 −2 h Mpc −1 ), and it peaks along the elongated configurations on small scales (k 1 ∼ 0.1 h Mpc −1 ). It peaks at the squeezed configuration on a very small scale (k 1 ∼ 1 h Mpc −1 ), but these scales are not accessible by the CMB due to the Silk damping. Note that the most squeezed + Note that our Φ is (−1) times Φ used in equation 72 of [136].
configuration shown in this Figure has k 1 = k 2 = 100k 3 . The dominant shape changes with scales, as the linear transfer function, T (k), declines with k, with the small-scale limit given by T (k) ∝ ln k/k 2 . From these results, we expect the second-order effect in the Newtonian limit to yield only a small contamination of f local NL . The dominant contribution to the second-order temperature anisotropy in the subhorizon limit is given by the second-order Sachs-Wolfe effect [136]: * where R * ≡ 3ρ b /(4ρ γ ) is the baryon-photon ratio at the decoupling epoch. (Here, a factor of 1/2 comes from our way of defining the second-order temperature anisotropy, ∆T = ∆T (1) + ∆T (2) and the second-order curvature perturbation, Φ = Φ (1) + 1 2 Φ (2) . This definition follows from Ref. [127].) The corresponding second-order a lm is where the linear primordial perturbation, Φ (1) p , is related to Φ (1) as Φ (1) (k) = Φ (1) p (k)T (k). Comparing this with equation 56, we identify the term inside the square bracket as the second-order radiation transfer function, F l ′ m ′ lm .♯ With this result, one can calculate the reduced bispectrum of the Newtonian secondorder effect, b 2nd l 1 l 2 l 3 . The resulting f local NL is always less than unity regardless of the angular scales (D. Nitta 2010, private communication; also see [140,141]). The calculation was done for l ≤ 2000.
How can we reconcile this result with those found in [137]? The calculations given above (equations 57 and 61) are valid only in the sub-horizon limits, and thus they are not suitable for calculating the contribution to the squeezed-limit bispectrum, which can correlate super-and sub-horizon fluctuations. Therefore, the difference between these results seems to imply: 1. The dominant contamination of f local NL comes from the general relativistic (post Newtonian) evolution of Φ (2) that is not captured by the above Newtonian calculation (equation 57).
2. The full second-order radiation transfer function beyond the sub-horizon approximation (equation 61) gives the dominant contribution to f local NL . * The CMB bispectrum from the second-order ISW effect was considered in [139,140].
♯ Incidentally, in the notation of [127] (see their equation 2.29), S . Shapes of the second-order bispectrum due to the second-order curvature perturbations in the Newtonian limit given in equation 60. Each panel shows the normalized amplitude of F 2nd (k 1 , k 2 , k 3 )(k 2 /k 1 ) 2 (k 3 /k 1 ) 2 as a function of k 2 /k 1 and k 3 /k 1 for a given k 1 , with a condition that k 3 ≤ k 2 ≤ k 1 is satisfied. The amplitude is normalized such that it is unity at the point where The CMB data are sensitive to k 1 < k max ∼ 0.2 h Mpc −1 (l max /2000), where the second-order bispectrum peaks at the equilateral configuration on large scales, and peaks along the elongated configurations on a smaller scale (k 1 ∼ 0.1 h Mpc −1 ). On a very small scale (k 1 ∼ 1 h Mpc −1 ), it peaks at the squeezed configuration. Note that the most squeezed configuration shown here has k 1 = k 2 = 100k 3 .
Perhaps both contributions are important. This is yet to be confirmed; however, if this is true, one should be able to construct a simple template for the second-order bispectrum, and use it to remove the contamination by including its amplitude, f 2nd NL , in the fit.
6. Four-point function: local-form trispectrum test of multi-field models Widely used notation for the "local-form trispectrum" is where k ij ≡ |k i + k j |. (In this section, we use Φ for the primordial perturbation, i.e., Φ = Φ p = 3 5 ζ.) When the curvature perturbation is given by the simplest local form, L , one finds the above trispectrum with τ N L = (6f local NL /5) 2 [44]. However, in general τ NL is different from (6f local NL /5) 2 .
To see this, let us consider a broad class of multi-field models in which the primordial curvature perturbation, ζ, is given in terms of the field derivatives of the number of efolds, N = ln a, and the perturbation in the I-th scalar field, δφ I , as where δφ I δφ J = 0 for I = J. This expansion is known as the "δN formalism" [22,88,114,142,143]. In this case, f local NL and τ NL are given by [22,46,144] where N ,I ≡ ∂N/∂φ I and N ,IJ ≡ ∂ 2 N/∂φ I ∂φ J . Suyama and Yamaguchi [49] showed that the Cauchy-Schwarz inequality implies that the following inequality, is satisfied. To derive this result, use the Cauchy-Schwarz inequality: with The equality, τ NL = (6f local NL /5) 2 is satisfied for the simplest local-form model, ζ = N 1 δφ 1 + (N 11 /2)δφ 2 1 . Note, however, that one finds a different relation between f local NL and τ NL when the Cauchy-Schwarz inequality becomes trivial, i.e., 0 = 0. For example, when ζ is given by [44] and δφ 1 δφ 2 = 0, one finds τ NL ∼ 10 3 (f local NL ) 4/3 [50]. The Cauchy-Schwarz inequality becomes 0 = 0 because a I = 0 for all I, and thus both the bispectrum and the trispectrum come from the second term in equation 71: the bispectrum is given by the 6-point function of δφ 2 , and the trispectrum is given by the 8-point function of δφ 2 .
In this case, whether τ NL ≥ (6f local NL /5) 2 is satisfied depends on the value of f local NL . For this particular example, the current limit of f local NL < 74 implies that τ NL ≥ (6f local NL /5) 2 is still satisfied parametrically.
The inequality is valid also for the "quasi-single field inflation" of [145]; see equation (7.6) of [69]. (Note that τ SE N L and τ CI N L in [69] correspond to τ NL and g NL in this paper, respectively.) ) and the point source bispectrum amplitude, b src (in units of 10 −5 µK 3 sr 2 ), from the WMAP 7-year temperature maps. This table is adopted from [11].
Therefore, if observations indicate τ NL < (6f local NL /5) 2 , a broad class of multi-field models satisfying the above conditions would be ruled out. This property makes the trispectrum a powerful probe of the physics of multi-field models. Note that the simplest local-form limit, τ NL = (6f local NL /5) 2 , has no special significance, as this is just one of many possibilities of multi-field models. (All single-field models predict a non-detectable level of the primordial f local NL , and thus the observational test using the above relation between f local NL and τ NL has no relevance to single-field models.) The expected 95% uncertainties in τ NL from the 7-year WMAP data (l max ∼ 500) and the 1-year Planck data (l max ∼ 1500) are 5000 and 560, respectively [73]. If the Planck finds f local NL ∼ 30, then it would be able to test if the measured τ NL would satisfy equation 67. This provides an excellent science case for the trispectrum that would be measured by Planck.
The expected uncertainties in g NL have not been calculated yet, although we expect them to be much greater than those for τ NL , as g NL is the coefficient of the cubic-order term (i.e., g NL is much more difficult to constrain than τ NL ) [146].
The local-form trispectrum is not the only possibility. Various other inflation models predict distinctly different quadrilateral shape dependence. For some general analyses of shapes, see [35,42,58,69].
Finally, while we do not discuss the large-scale structure of the universe in this article, the most promising probe of the local-form trispectrum seems to be the bispectrum of galaxies. See [86,147,148] for details.

Bispectrum
(Most of this subsection is adopted from Section 6.2 of [11].) In 2002, the first limit on f local NL was obtained from the COBE 4-year data [149] by [150], using the angular bispectrum. The limit was improved by an order of magnitude when the WMAP first year data were used to constrain f local NL [151]. Since then the limits have improved steadily as WMAP collects more years of data and the bispectrum method for estimating f local NL has improved [6,7,10,85,[152][153][154][155][156]. Using the optimal estimators described in Section 3, we have constrained the primordial non-Gaussianity parameters as well as the point-source bispectrum using the WMAP 7-year data. The 7-year data and results are described in Refs. [11,77,[157][158][159][160].
We use the V-and W-band maps at the HEALPix resolution N side = 1024. As the optimal estimator weights the data optimally at all multipoles, we no longer need to choose the maximum multipole used in the analysis, i.e., we use all the data. We use both the raw maps (before cleaning foreground) and foreground-reduced (clean) maps to quantify the foreground contamination of f NL parameters. For all cases, we find the best limits on f NL parameters by combining the V-and W-band maps, and marginalizing over the synchrotron, free-free, and dust foreground templates [158]. As for the mask, we always use the KQ75y7 mask [158].
In Table 7.1, we summarize our results: 1. Local form results. The 7-year best estimate of f local NL is f local NL = 32 ± 21 (68% CL). The 95% limit is −10 < f local NL < 74. When the raw maps are used, we find f local NL = 59±21 (68% CL). When the clean maps are used, but foreground templates are not marginalized over, we find f local NL = 42 ± 21 (68% CL). These results (in particular the clean-map versus the foreground marginalized) indicate that the foreground emission makes a difference at the level of ∆f local NL ∼ 10. † † We find that the V+W result is lower than the V-band or W-band results. This is possible, as the V+W result contains contributions from the cross-correlations of V and W such as VVW and VWW . As for the point-source bispectrum, we do not detect b src in V, W, or V+W. In [85], we estimated that the residual sources could bias f local NL by a small positive amount, and applied corrections using Monte Carlo simulations. In this paper, we do not attempt to make such corrections, but we note that sources could give ∆f local NL ∼ 2 (note that the simulations used by [85] likely overestimated the effect of sources by a factor of two). As the estimator has changed from that used by [85], extrapolating the previous results is not trivial. Source corrections to f equil NL and f orthog NL could be larger [85], but we have not estimated the magnitude of the effect for the 7-year data.
As we described in Section 4, among various sources of secondary non-Gaussianities which might contaminate measurements of primordial non-Gaussianity (in particular f local NL ), a coupling between the ISW effect and the weak gravitational lensing is the most dominant source of confusion for f local NL . Calabrese et al. [162] used the skewness power spectrum method of [163] to search for this term in the WMAP 5-year data and found a null result.

Trispectrum
The optimal estimators for the trispectrum have not been implemented, largely because they are computationally demanding. While the first measurements of the angular bispectrum were made from the the COBE 4-year data [149] by [1,164] in 2001, limits on the physical parameters have not been obtained from the direct trispectrum analysis.

Other statistical methods
While the optimal estimators for the f NL parameters (with the minimum variance) must be constructed from the PDF as in Section 3, there are various other ways of constraining non-Gaussianity. While these other methods are usually sub-optimal, they serve as useful diagnosis tools of the results obtained from the direct bispectrum and trispectrum methods. In some cases, they are easier to implement than the optimal estimators. A major progress in the topological Gaussianity test using the Minkowski functionals [166][167][168][169][170] since 2004 is the derivation and implementation of the analytical formula for the Minkowski functionals of the CMB [171,172]. This method has been applied to the WMAP data [173,174] as well as to the BOOMERanG data [175]. The Planck data are expected to reach the 68% limit of ∆f local NL = 20 [171], which is worse than the limit from the optimal method, ∆f local NL = 5 [82]. An advantage of the Minkowski functionals is that the measurements of the Minkowski functionals do not depend on the models, and thus the computational cost is the same for all models. This allows one to obtain limits on various models, for which the optimal estimators are difficult to implement. For example, a limit on the primordial non-Gaussianity in the isocurvature perturbation is currently available only from the Minkowski functionals [174].
Instead of expanding the temperature anisotropy into spherical harmonics, one may choose to expand it using a different basis. One popular basis used in the CMB community is the so-called Spherical Mexican Hat Wavelet (SMHW). See [176,177] for reviews on this method. A major progress in this method is the realization that the 3-point function of the wavelet coefficients made of large and small smoothing scales is nearly an optimal estimator for the local-form bispectrum: when only the adjacent scales are included, Curto et al. [178] found −8 < f local NL < 111 (95% CL) from the WMAP 5year data. When all the scales (including large-small scale combinations) are included in the analysis, the limit improved significantly to −18 < f local NL < 80 (95% CL) [179], which is similar to the optimal limit from the 5-year data, −4 < f local NL < 80 (95% CL) [10]. An advantage of the SMHW is that it retains information on the spatial distribution of the signal. This property can be used to measure f local NL as a function of positions on the sky [179]. In addition, the analytical formula for the SMHW as a function of f local NL has been derived and implemented (A. Curto 2009, private communication). See [180,181] for earlier limits on f local NL from the SMHW. Another form of spherical wavelets that has been used to constrain f local NL is the spherical needlets [182]. The limits on f local NL are reported in [161,183,184]. This method also allows one to look for a spatial variation in f local NL , and the results are reported in [185,186]. For the other types of wavelets considered in the literature, see [187,188] and references therein.
Many other statistical methods have been proposed and used for constraining f local NL in the literature. An incomplete list of references is: [189,190] on the real-space 3-point function; [191] on the integrated bispectrum; [192] on the 2-1 cumulant correlator; [181] on the local curvature; and [74,193] on the N-point PDF. Also see references therein. While we have not listed the statistical methods that have not been used to constrain the primordial non-Gaussianity parameters yet, there are many other methods proposed in a general context in the literature.

Conclusion
Since the last review articles on signatures of primordial non-Gaussianity in the CMB were written in 2001 [1] and 2004 [2], a lot of progress has been made in this field. The current standard lore may be summarized as follows: 1. Shape and physics. Different aspects of the physics of the primordial universe appear in different shapes of three-and four-point functions.
2. Importance of local shape. Of these shapes, the local shapes have special significance: a significant detection of the local-form bispectrum (with f local NL ≫ 1) would rule out all single-field inflation models, and the local-form trispectrum can be used to rule out a broad class of multi-field models (if not all multi-field models) by testing τ NL ≥ (6f local NL /5) 2 .
3. Optimal estimators. The optimal estimators of the bispectrum and trispectrum can be derived systematically from the expansion of the PDF. The optimal bispectrum estimator has been implemented. 4. Secondary. The most serious contamination of f local NL is due to the lensing-ISW coupling, which can be removed by using the template given in Section 4.

5.
Foreground. The Galactic foreground contamination is minimal for f equil NL and f orthog NL , but it can be as large as f local NL ∼ 10 for the local-form bispectrum. This must be carefully studied and eliminated in the Planck data analysis. The random (Poisson) point-source contamination can be removed by using the template given in Section 3.2. Some outstanding issues for the "CMB and primordial non-Gaussianity" include: 1. Second order. (In Newtonian gauge) the products of the first-order terms and the intrinsically second-order terms in the sub-horizon limit do not contaminate the local-form bispectrum very much (∆f local NL < 1). However, would the post-Newtonian effect give ∆f local NL ∼ 5, as found by [137]? If so, we need to construct a template for this effect.
2. More foreground. How can we model the non-Poisson (clustered) point source bispectrum? How about the foreground (and secondary) contamination of the primordial trispectrum?
3. Trispectrum estimators. How can we implement the optimal trispectrum estimators for both local and non-local shapes?
These issues would become important when the Planck data are analyzed in search of primordial non-Gaussianity. The Planck is expected to reduce the uncertainty in f local NL by a factor of four compared to the current limit, f local NL = 32 ± 21 (68% CL). If the Planck detected f local NL ∼ 30, then the trispectrum would provide an important test of multi-field models. In particular, if f local NL ≫ 1 and τ NL < (6f local NL /5) 2 are found, then all single-field models and many (if not all) multi-field models would be ruled out, and thus the standard paradigm of inflation as the origin of fluctuations would face a serious challenge.
However, do not despair even if the Planck did not detect the primordial bispectrum or trispectrum -while the CMB may end its leading role as a probe of primordial non-Gaussianity (unless the next-generation, comprehensive CMB satellite which can measure both the temperature and polarization to the cosmic-variance-limited precision is funded [194]), the large-scale structure of the universe would eventually take over and substantially reduce the uncertainties in the local-form parameters such as f local NL , τ NL , and g NL (see Desjacques and Seljak's article in this volume). presentations, D. Jeong for making Figure 3 and 4, D. Nitta for calculating f local NL from the second-order Φ in the Newtonian limit, and T. Takahashi for valuable comments on Section 6. This work is supported in part by an NSF grant PHY-0758153.