#### Abstract

The most direct probe of non-Gaussian initial conditions has come from bispectrum measurements of temperature fluctuations in the Cosmic Microwave Background and of the matter and galaxy distribution at large scales. Such bispectrum estimators are expected to continue to provide the best constraints on the non-Gaussian parameters in future observations. We review and compare the theoretical and observational problems, current results, and future prospects for the detection of a nonvanishing primordial component in the bispectrum of the Cosmic Microwave Background and large-scale structure, and the relation to specific predictions from different inflationary models.

#### 1. Introduction

The standard inflationary paradigm predicts a flat Universe perturbed by nearly-Gaussian and scale-invariant primordial perturbations. These predictions have been verified to a high degree of accuracy by Cosmic Microwave Background (CMB) and Large-Scale Structure (LSS) measurements, such as those provided by the Wilkinson Microwave Anisotropy Probe (WMAP) [1], the 2dF Galaxy Redshift Survey (2dFGRS) [2], and the Sloan Digital Sky Survey (SDSS) [3]. Despite this success, it has proved to be difficult to discriminate between the vast array of inflationary scenarios that have been proposed by high-energy theoretical investigations or even to rule-out alternatives to inflation. Since most of the present constraints on the Lagrangian of the inflaton field have been obtained from measurements of the two-point function, or power spectrum, of the primordial fluctuations, a natural step to extend the available information is to look at non-Gaussian signatures in higher-order correlators.

The lowest-order additional correlator to take into account is the three-point function or its counterpart in Fourier space, the* bispectrum*. Most models of inflation are characterized by specific predictions for the bispectrum of the primordial perturbations in the gravitational potential . The bispectrum of these perturbations is defined as where we have introduced the notation so that the Dirac delta function here is . Together with the assumption of statistical homogeneity and isotropy for the primordial perturbations, this implies that the bispectrum is a function of the triplet defined by the magnitude of the wavenumbers , and forming a closed triangular configuration. The current constraints that we are able to derive on the bispectrum provide additional information about the early Universe; the possible detection of a non-vanishing primordial bispectrum in future observations would represent a major discovery, especially as it is predicted to be negligible by standard inflation.

The cosmological observable most directly related to the initial curvature bispectrum is given by the bispectrum of the CMB temperature fluctuations, which provide a map of the density perturbations at the time of decoupling, the earliest information we have about the Universe. Current measurements of individual triangular configurations of the CMB bispectrum are, however, consistent with zero. Studies of the primordial bispectrum, therefore, are usually characterized by constraints on a single-*amplitude parameter*, denoted by once a specific model for is assumed. Since most models predict a curvature bispectrum obeying the hierarchical scaling , with being the curvature power spectrum, the non-Gaussian parameter roughly quantifies the ratio defining the “strength” of the primordial non-Gaussian signal. In addition, we can write where encodes the functional dependence of the primordial bispectrum on the specific triangle configurations. For brevity, the characteristic shape-dependence of a given bispectrum is often referred to simply as the * bispectrum shape* (a precise definition of the bispectrum shape function will be given in Section 2.1). Inflationary predictions for both the amplitude and the shape of are strongly model dependent. Notice that the subscript “” stands for “nonlinear”, since a common phenomenological model for the non-Gaussianity of the initial conditions can be written as a simple nonlinear transformation of a Gaussian field. Generically, of course, non-Gaussianity is associated with nonlinearities, such as nontrivial dynamics during inflation, resonant behaviour at the end of inflation (“preheating”), or nonlinear postinflationary evolution. At the very least, future CMB and LSS observations are expected to be able to eventually detect the last of the three effects mentioned above.

Perturbations in the CMB provide a particularly convenient test of the primordial density field because CMB temperature and polarization anisotropies are small enough to be studied in the *linear regime* of cosmological perturbations. Once the effects of foregrounds are properly taken into account, a non-vanishing CMB bispectrum at large scales would be a direct consequence of a non-vanishing primordial bispectrum. As we will see, while other CMB probes of primordial non-Gaussianity are available, such as tests of the topological properties of the temperature map based on Minkowski Functionals or measurements of the CMB trispectrum, the estimator for the non-Gaussian parameter has been shown to be optimal. We will focus mostly on this bispectrum estimator in the section of this paper dedicated to the CMB.

In the standard cosmological model, the large-scale structure of the Universe, that is, the distribution of matter and galaxies on large scales, is the result of the nonlinear evolution due to gravitational instability of the same initial density perturbations responsible for the CMB anisotropies. This is, perhaps, the most important prediction of the inflationary framework which provides a common origin for the CMB and large-scale structure perturbations as the result of tiny quantum fluctuations stretched over cosmological scales during a phase of accelerated expansion. The large-scale structure we observe at low redshift, however, is characterized by large voids and small regions with very large-matter density, and it is therefore a much less direct probe of the initial conditions. The *distribution of matter becomes a highly non-Gaussian field* precisely as a result of the nonlinear growth of structures, * even for Gaussian initial conditions*. This non-Gaussianity is expressed, in particular, by a non-vanishing matter bispectrum at* any* measurable scale, including the largest scales probed by current or future redshift surveys. In this context, the effect of primordial non-Gaussianity, that is, of an initial component in the curvature bispectrum, will constitute a* correction* to the galaxy bispectrum. It follows that the possibility of constraining or detecting this initial component is strictly related to our ability to *distinguish* it from other primary sources of non-Gaussianity, that is, the nonlinear gravitational evolution, and, in the case of galaxy surveys, nonlinear bias.

The study of non-Gaussian initial conditions for large-scale structure has a relatively long history, with important contributions going back to the mid eighties. The standard picture that has been developed over the years assumed that, at large scales, the effect of primordial non-Gaussianity on the galaxy distribution is simply given in terms of an additional component to the galaxy bispectrum. This is obtained, in perturbation theory, as the linearly evolved and linearly biased initial matter bispectrum, related to the curvature bispectrum by the Poisson equation. Such component becomes subdominant as the gravity-induced non-Gaussian contribution grows in time. In this framework, as one can expect, high-redshift and large-volume galaxy surveys would constitute the best probes of the initial conditions. It has been shown, in fact, that proposed and planned redshift surveys, such as those of Euclid [4], should be able to provide constraints on the primordial non-Gaussian parameters comparable to, if not better than, those expected from CMB missions such as those of Planck. What is more important, in the event of a detection by Planck, is that confirmation by large-scale structure observations will be required.

Recent results from N-body simulations with non-Gaussian initial conditions, however, have revealed a more complex picture. The effect of primordial non-Gaussianity at large scales is not limited to an additional contribution to the galaxy bispectrum, but it quite dramatically affects the galaxy bias relation itself, that is, the relation between the matter and galaxy distributions. A surprising consequence is that it induces a large correction even for the galaxy power spectrum. Such an effect has attracted considerable recent attention and, remarkably, has placed constraints on the non-Gaussian parameter from current LSS datasets which already appear to marginally improve on CMB limits. However, from a theoretical point of view, a proper understanding of the phenomenon is not fully developed yet. For example, reliable predictions for the galaxy bispectrum are not yet available. Most importantly, as for general cosmological parameter estimation, a complete likelihood analysis aimed at constraining, or detecting, primordial non-Gaussianity in large-volume redshift surveys should involve* joint* measurements of the galaxy power spectrum and bispectrum, as well as possibly higher-order correlation functions. While we are still far from a proper assessment of what such analysis would be able to achieve, current results in this direction are very encouraging.

This review is divided into four parts. In Section 2 we will first discuss initial conditions as defined in terms of the primordial curvature bispectrum and its phenomenology. We will then review the observational consequences of primordial non-Gaussianity on the CMB bispectrum, Section 3, and on the large-scale structure bispectrum as measured in redshift surveys, Section 4. In both cases we will discuss theoretical models for the observed bispectra and technical problems related to the estimation of the non-Gaussian parameters, with the differences that naturally characterize such distinct observables. We also give an example of joint analysis using both CMB and large-scale structure when we consider the possibility of constraining a strongly scale-dependent non-Gaussian parameter , emerging in some recently proposed inflationary models.

#### 2. Initial Conditions and the Primordial Bispectrum

In this section we will briefly overview the main predictions of inflationary models regarding the non-Gaussianity (NG) of the primordial curvature perturbation field. The link between NG of primordial density fluctuations and NG of CMB and LSS will be shown in following sections. In order to provide a full description of an NG random field, all correlators beyond the 2-point function are in principle necessary. However in this review we will focus on the primordial bispectrum (i.e., three-point function in Fourier space). This is not only justified by the fact that the bispectrum is the first and simplest higher-order correlator to look at, but also by the fact that most models of inflation predict vanishingly small correlators beyond the bispectrum. In Section 2.1 we will introduce the relevant quantities, their mathematical definitions, and we will provide a general overview and classification of the bispectra predicted in different inflationary scenarios (only from a purely mathematical point of view, without linking them to the Physics originating them at this stage). Finally, a useful eigenmode expansion technique for bispectra will be introduced in Section 2.2 and applied to the calculation of correlations between different bispectra in Section 3.4. In the same section we will also show which kinds of bispectra are predicted by different models of inflation.

##### 2.1. The Primordial Bispectrum and Shape Function

The starting point for this discussion is the primordial gravitational potential perturbation which was seeded by quantum fluctuations during inflation or by some other mechanism in the very early Universe (. When characterizing the fluctuations we usually work in Fourier space with the (flat space) transform defined through The primordial power spectrum of these potential fluctuations is found using an ensemble average: where we have assumed that physical processes creating the fluctuations are statistically isotropic so that only the dependence on the wavenumber remains . Recall that, for nearly scale-invariant perturbations, the fluctuation variance on the horizon scale is almost constant , implying that .

The primordial bispectrum is found from the Fourier transform of the three-point correlator as Here, the delta function enforces the triangle condition, that is, the constraint that the wavevectors in Fourier space must close to form a triangle: . Examples of such triangles are shown in Figure 1, illustrating the basic squeezed, equilateral, and flattened triangles to which we will refer later. Note that a specific triangle can be completely described by the three lengths of its sides and so, in the isotropic case, we are able to describe the bispectrum using only the wavenumbers . The triangle condition restricts the allowed wavenumber configurations to the interior of the tetrahedron illustrated in Figure 2.

**(a)**

**(b)**

**(c)**

The most studied primordial bispectrum is the* local model* in which contributions from “squeezed” triangles are dominant, that is, with, for example, (as illustrated in Figure 1(a)). This is well motivated physically as it encompasses “superhorizon” effects during inflation when a large-scale mode (say) which has exited the Hubble radius exerts a nonlinear influence on the subsequent evolution of smaller-scale modes . Although this effect is small in single-field slow-roll inflation, it can be much larger for multifield models. In a weakly coupled regime, the potential can be split into two components: the linear term , representing a Gaussian field, giving the usual perturbation results, plus a small local non-Gaussian term [5]: where is called the nonlinearity parameter. In Fourier space, the nonlinear term is then given by the convolution From this we can infer, using (4), that the only non-vanishing contributions to the bispectrum (5) take the form In the scale-invariant case the power spectrum of the primordial potential takes the form , where defines the amplitude of primordial fluctuations at the end of inflation. Accounting for permutations, the local bispectrum then becomes Although this is a rather pathological function which diverges along the edges of the tetrahedron (i.e., when any ), we can infer from it some basic properties of the bispectrum for any model which is nearly scale invariant. For example, we can observe that the bispectrum at equal has the characteristic scaling If we remove this overall scaling by multiplying (9) by the factor , then we note that on transverse slices through the tetrahedron defined by (see Figure 2 ) the bispectrum only depends on the ratios of the wavenumbers, say, and . Indeed, it can prove convenient to characterize the bispectrum in terms of the following transverse parameters [6, 7]: with the domains and . The volume element on the regular tetrahedron of allowed wavenumbers then becomes .

These considerations lead naturally to the definition of the primordial shape function [8] where is a normalization factor which is often chosen such that is unity for the equal case; that is, (we will discuss alternatives to this rather arbitrary convention later). For example, the canonical “local” model (9) has the shape Thus it is usual to describe the primordial bispectrum in terms of an overall amplitude and a transverse two-dimensional shape , which incorporates any distinctive momentum dependence. Of course, if there is a nontrivial scale dependence, then the full three-dimensional dependence of on must be retained.

There are other physically well-motivated shapes in the literature which have also been extensively studied. The simplest shape is the *constant model * which, like the local model, has a large-angle analytic solution for the CMB bispectrum [9]. The local model tends to be the benchmark against which all other models are compared and normalized, but for practical purposes the constant model is much more useful, given its regularity at both late and early times. The *equilateral shape* is another important case with [8] While being not derived directly from a physical model, it has been chosen phenomenologically as a separable ansatz for higher-derivative models [10] and DBI inflation [11]. The equilateral shape is contrasted with the local model in Figure 3.

**(a)**

**(b)**

Another important early result was the primordial bispectrum shape for single-field slow-roll inflation derived by Acquaviva et al. [12] and Maldacena [13]: where are the usual slow-roll parameters. In the second line, we have noted that this shape can be accurately represented as the superposition of local and equilateral shapes. The coefficients in (16), which include the scalar spectral index , confirm that and so standard single slow-roll inflation cannot produce an observationally significant signal. Nevertheless, it is interesting to determine which shape is dominant in (16) and to what extent other primordial shapes are independent from one another.

Whether two different primordial shapes can be distinguished observationally can be determined from the correlation between the corresponding two CMB bispectra weighted for the anticipated signal-to-noise-ratio, as in the estimator (see next section) and the Fisher matrix analysis (see Section 3.8). However, direct calculations of the CMB bispectrum can be very computationally demanding. A much simpler approach is to determine the independence of the two shape functions and from the correlation integral (see [9], and also the study by Babich et al. in [8]) where we choose the weight function to be reflecting the primary scaling of the CMB correlator. The shape correlator is then defined by Here, the integral is over the tetrahedral region shown in Figure 2 taken out to a maximum wavenumber corresponding to the experimental range for which forecasts are sought (with with being the present-day conformal time). The weight function appropriate for mimicking the large-scale structure bispectrum estimator (see Section 4.3.2) would be different with varying scaling laws introduced by the transfer functions for wavenumbers above and below , the inverse comoving horizon at equal matter-radiation. Nevertheless, the weight given in (18) provides a compromise between these scalings, and so shape correlation results should offer a useful first approximation.

Below we will survey primordial models in the literature, showing how close the shape correlator comes to a full Fisher matrix analysis. However, here we note that the local shape (13) and the equilateral shape (54) have only a modest 46% correlation. For the natural values of the slow-roll parameters we find the somewhat surprising result that is 99.7% correlated with (and it cannot be easily tuned otherwise because is not consistent with deviations from scale invariance favored observationally ). Such strong correspondences are important in defining families of related primordial shapes, thus reducing the number of different cases for which separate observational constraints must be sought.

##### 2.2. General Primordial Bispectra and Separable Mode Expansions

The three shape functions (13), (14), and (54) quoted above share the important property of* separability*; that is, they can be written in the form or as the sum of just a few such terms. As we will see, if a shape is separable, then the computational cost of evaluating the corresponding CMB bispectrum is dramatically reduced. In fact, without this property, the task of estimating whether a nonseparable bispectrum is consistent with observation appears to be intractable (for large ). Of course, the number of models which can be expressed directly in the form of (20) is very limited, despite the usefulness of approximate ansätze such as the equilateral shape (54). Indeed, approximating nonseparable shapes by educated guesses for the separable functions is neither systematic nor computationally efficient (because arbitrary nonscaling functions create numerical difficulties, as we will explain later).

Instead, we will present a separable mode expansion approach for efficient calculations with any nonseparable bispectrum, as described in detail by Fergusson et al. in [14] (and originally proposed in [6]). Our aim will be to express any shape function as an expansion in mode functions where, here, for convenience, we have represented the symmetrized products of the separable basis functions as with a one-to-one mapping ordering the products as . The important point is that must be an independent set of well-behaved basis functions which can be used to construct complete and orthogonal three-dimensional eigenfunctions on the tetrahedral region defined by (see Figure 2) The introduction of the cutoff at is motivated by both separability and the correspondence with the observational domain In the shape correlator (19), we have already seen what is essentially an inner product between two shapes on this tetrahedral region, which we can define for two functions as

Satisfactory convergence for known bispectra can be found by using simple polynomials in the expansion (21), that is, using analogues of the Legendre polynomials on the domain (23). With unit weight, the polynomials satisfying can be found by generating functions with the first three given by [14] The first few polynomials are plotted in Figure 4, where they are contrasted with the Legendre polynomials.

The three-dimensional separable basis functions in (2) reflect the six symmetries of the bispectrum through the permuted sum of the product terms. They could have been constructed directly from simpler polynomials, such as however, the polynomials have two distinct advantages. First, the 's confer partial orthogonality on the and, secondly, these remain well behaved when convolved with transfer functions.

In order to rapidly decompose an arbitrary shape function into the coefficients , it is more convenient to work in a nonseparable orthonormal basis (. These can be derived directly from through Gram-Schmidt orthogonalization, so that with being a lower triangular matrix (see [14]). Thus we can find the unique shape function decomposition

In the orthonormal frame, Parseval's theorem ensures that the autocorrelator is simply . Hence, with a simple and efficient prescription we can construct separable and complete basis functions on the tetrahedral domain (23) providing rapidly convergent expansions for any well-behaved shape function . These eigenmode expansions will prove to be of great utility in subsequent sections. Examples of this bispectral decomposition and its rapid convergence for the equilateral and DBI models are shown in Figure 5.

**(a)**

**(b)**

##### 2.3. Families of Primordial Models and Their Correlations

We will now briefly survey the main categories of primordial models in the literature and their relative independence, closely following the discussion by Fergusson and Shellard in [9].

###### 2.3.1. The Constant Model

The constant model (14) is the simplest possible primordial shape with triangles of every configuration contributing equally to the bispectrum ; it is the equipartition model. The constant model was motivated initially by its simplicity [9] leading to an analytic solution for the large-angle CMB bispectrum, as well as due to its close correlation with equilateral models. However, the shape does have a more explicit physical motivation in at least one context [15], during multifield inflation for a slowly turning trajectory (denoted as quasi-single-field inflation). For multifield inflation, it is well known that the conversion of isocurvature fluctuations into curvature fluctuations during “corner-turning” can source significant non-Gaussianity (see, e.g., [7, 16]). In the quasi-single-field case with mass isocurvature modes, a detailed investigation of the ongoing conversion into the curvature mode demonstrated that novel shapes could be generated [15], amongst them are the shapes which were very nearly constant. Generically, these model-dependent shapes belonged to a one-parameter family which interpolated nontrivially between equilateral (54) and local (13) shapes (see also [17, 18]). This is an important caveat for the present discussion, because non-Gaussian searches could uncover shapes intermediate between the categories we will discuss below.

###### 2.3.2. Equilateral Triangles—Centre-Weighted Models

Bispectra dominated by contributions from nearly equilateral triangle configurations, can be fairly easily characterized analytically and are the most amenable to CMB searches. However, equilateral non-Gaussianity requires that the amplification of nonlinear effects around the time modes exit the horizon, which is not possible in a slow-roll single-field inflation. Instead, the kinetic terms in the effective action must be modified as in the Dirac-Born-Infeld (DBI) model [11] or by explicitly adding higher-derivative terms, such as in K-inflation (see, e.g., [19]). The resulting corrections modify the sound speed and inflation is able to take place in steep potentials. For DBI inflation, this leads to non-Gaussianity being produced with a shape function of the form [10, 11] Another example of a model with nonstandard kinetic terms is ghost inflation [20] with a derivatively coupled field driving inflation and a trilinear term in the Lagrangian creating a nonzero equilateral-type shape tending towards constant.

General non-Gaussian shapes arising from modifications to single-field inflation have been extensively reviewed in [19]. Using a Lagrangian that was an arbitrary function of the field and its first derivative, they were able to identify six distinct shapes describing the possible non-Gaussian contributions. Half of these had negligible amplitude being of the order of slow-roll parameters (with two already given in (16)). Of the remaining three shapes (see [19], and also [21]), one was believed to be subdominant and the second recovered the DBI shape (27), leaving a third distinct single-field shape which is the inverse of the local shape (13): . Finally, we recall the original equilateral shape (54), noting that it was introduced not because of a fundamental physical motivation, but as a separable approximation to the DBI shape (27) [8].

Despite the apparent visual differences between these shapes (see [9]), particularly near the edges of the tetrahedral domain, the shape correlator (19) reveals at least a 95% or greater correlation of the DBI, ghost, and single shapes to the equilateral shape (54) (consistent with results in [8, 22]). Comparative results between the shape correlator are given in Table 1 (together with the corresponding CMB correlation results brought forward and showing the efficacy of these estimates). These particular centre-weighted shapes must be regarded as a single class which would be extremely differentiate observationally, without a bispectrum detection of very high significance.

Finally, we comment on the “orthogonal” shape proposed by Smith et al. in [18], together with for characterizing single-field inflation models with an approximate shift symmetry (see also [19]). This shape is approximately , which means that it is very similar to an earlier study of flattened shapes [23] which proposed an “enfolded” shape with . From the eigenmode decomposition (26) of the equilateral model shown in Figure 5, it is clear how the degree of correlation can be altered by subtracting out the important constant term. With the specific choice of constant term used in the analysis by Smith et al. [18], one gets a correlation of about 30% between both local and equilateral shapes and the orthogonal ansatz.

###### 2.3.3. Squeezed Triangles—Corner-Weighted Models

The local shape covers a wide range of models where the non-Gaussianity is produced by local interactions. These models have their peak signal in “squeezed” states where one is much smaller than the other two due to non-Gaussianity typically being produced on superhorizon scales. We have already observed that single-field slow-roll inflation (16) is dominated by the local shape [24], though is tiny [12, 13, 24]. The production of non-Gaussianity during multiple-field inflation [7, 16, 21, 25–30] shows much greater promise through conversion of isocurvature into adiabatic perturbations (see, e.g., recent works in [15, 17, 31, 32] and references therein). The magnitude of the non-Gaussianity generated is normally around , which is at the limit for Planck detection, but models can be tuned to create larger signals. Significant can be produced in curvaton models with [33–35]. Large can also be generated at the end of inflation from massless preheating or other reheating mechanisms [36–38].

We note that local non-Gaussianity can also be created in more exotic scenarios. Models based on nonlocal field theory, such as -adic inflation, can have inflation in very steep potentials. Like single-field slow-roll inflation, the predicted “nonlocal” shape function is a combination of a dominant local shape (13) and an equilateral shape (54) (see, e.g., [39–42]). The ekpyrotic model can also generate significant [43–47]. Here the density perturbations are generated by a scalar field rolling in a negative exponential potential, so nonlinear interactions are important with .

In using the shape correlator for the local model, we must introduce a small-wavenumber cutoff, taken to be ; otherwise the shape correlator becomes infinite. This logarithmic divergence does not afflict the CMB bispectrum because we do not consider contributions below the quadrupole (a threshold which is approximated by the primordial cutoff). The local shape is modestly correlated at the 40%–55% level with the equilateral shapes, mainly through the constant term in the expansion (26). As can be seen in Table 2, this somewhat underestimates the CMB correlator. Nevertheless, a NG signal of only modest significance should be able to distinguish between these independent models.

Finally, warm inflation scenarios, that is, models in which dissipative effects play a dynamical role, are also predicted to produce significant non-Gaussianity [48, 49]. Contributions are again dominated by squeezed configurations but with a different more complex shape possessing a sign flip as the corner is approached (see Figure 6). This makes the warm and local shapes essentially orthogonal with only a 33% correlation. Again, in using the shape correlator, we need to introduce the same phenomenological cutoff as for the local model, but we also note the more serious concern which is the apparent breakdown of the approximations used to calculate the warm inflation shape near the corners and edges.

**(a)**

**(b)**

###### 2.3.4. Flattened Triangles—Edge-Weighted Models

It is possible to consider inflationary vacuum states which are more general than the Bunch-Davies vacuum, such as an excited Gaussian (and Hadamard) state (see [50], and also discussions in [19, 23]). Observations of non-Gaussianity in this case might provide insight into trans-Planckian physics. The proposed shape for the bispectrum is The bispectrum contribution from early times is dominated by flattened triangles, with, for example, and for a small sound speed can be large. Unfortunately, as the divergent analytic approximation breaks down at the boundary of the allowed tetrahedron, some form of cutoff must be imposed, as shown for the smoothed shape in Figure 6 where an edge truncation has been imposed together with a Gaussian filter. The lack of compelling physical motivation and ill-defined asymptotics make predictions for this model uncertain.

###### 2.3.5. Features—Scale-Dependent Models

There are also models in which the inflation potential has a feature, providing a break from scale invariance. This can take the form of either a step [51] or a small oscillation superimposed onto the potential [52]. Analytic forms for both by these three-point functions have been presented by Chen et al. in [53] with one approximation taking the form where is the associated scale of the feature in question and is a phase factor. Results for the shape correlator for a particular feature model (with and ) are given in Table 2, showing that it is essentially independent of all of the other shapes. Clearly, scale-dependent feature models form a distinct fifth category beyond equilateral, local, warm, and flat shapes.

#### 3. Cosmic Microwave Background

Non-Gaussian initial conditions at the end of inflation (or produced by alternative models for the generation of primordial perturbations) can produce observable signatures in both the CMB and LSS. It is clear that an eventual detection of such signatures would be of great scientific interest. Cosmological measurements of primordial NG would indeed allow to constrain and discriminate between the different candidate scenarios of primordial inflation that have been briefly reviewed in the previous section in terms of their bispectrum prediction. This section is devoted to study the CMB bispectrum produced by a non-Gaussian primordial curvature perturbation field. We will start in Section 3.1 by determining how the primordial curvature bispectrum propagates to the observed bispectra of CMB temperature and polarization anisotropies. At the end of this section, we will obtain a formula expressing the CMB bispectrum as a convolution of the primordial one with suitable *radiation transfer functions*. The latter encode all of the radiative and gravitational effects producing the observed pattern of CMB anisotropies starting from a given primordial potential. This result will not come as a surprise to the reader familiar with CMB theory, since it is completely analogous to the formula relating the power spectrum of primordial curvature perturbations to the CMB angular power spectrum. We will be simply recasting the same formalism in terms of higher-order correlator. Armed with this useful relation, we will revisit the various models described in the previous section and calculate their corresponding CMB bispectra. This will be done in Sections 3.2 and 3.3 where we will employ the useful distinction between separable and nonseparable bispectra already introduced before and discuss its various implications. We will then compare the bispectra predicted in different scenarios in order to find out whether different models produce observationally distinguishable shapes (and thus whether primordial NG is a viable tool to discriminate between different theories of the Early Universe). This will be done in Section 3.4 by means of a suitably defined * shape correlator* that will tell us “how similar” two CMB bispectrum shapes are. After this preliminary work of definition and classification of various bispectra, in Section 3.5 we will finally deal with the problem of extracting the CMB bispectrum from the data in order to produce statistical estimates of the level of primordial non-Gaussianity and compare the NG measurements to theoretical predictions.

##### 3.1. The CMB Bispectrum

In this section we will study the connection between the primordial bispectrum at the end of inflation and the observed bispectrum of CMB anisotropies . Our work will be primarily concerned with the analysis of the three-point function induced by a NG primordial gravitational potential in the CMB temperature fluctuation field. Temperature anisotropies are represented using the coefficients of a spherical harmonic decomposition of the cosmic microwave sky: Analogous expansions are performed for the -mode polarization field in order to produce polarization multipoles . For simplicity and clarity, throughout most of this review we will focus on the temperature multipoles and omit the superscript for convenience of notation. However we stress here that all of the considerations we make in the following can be readily applied to polarization multipoles and related bispectra. More discussion about this subject can be found in Section 3.8.2.

The primordial potential is imprinted on the CMB multipoles by a convolution with transfer functions representing the linear perturbation evolution, through the integral The radiation transfer functions encode all of the typical effects observed in the CMB power spectrum at linear order, that is, the Sachs-Wolfe effect, Integrated Sachs-Wolfe effect, acoustic peaks, and silk damping (see, e.g., [54, 55]). An equation identical to (31) produces the E-mode polarization CMB multipoles starting from the primordial temperature fluctuation field, provided that polarization transfer functions replace temperature transfer functions in the convolution above. It is sometimes useful to rewrite (31) in position, rather than Fourier, space. In this case it is straightforward to show that (31) becomes where, starting from the primordial potential , we transform from Cartesian into polar coordinates and defined In this expression is the spherical Bessel function of order . The CMB bispectrum is the three-point correlator of , so by substituting, we obtain where in the last line we have integrated over the angular parts of the three , having inserted the exponential integral form for the delta function in the bispectrum definition (5). The last integral over the angular part of is known as the Gaunt integral, which can be expressed in terms of Wigner- symbols as (for more details on these functions and their properties, see, e.g., [56] and references therein)

Given that most theories we will consider are assumed to be isotropic, the -dependence can be factorized out of the physically relevant part of the bispectrum [57]. It is then usual to work with the angle-averaged bispectrum, Or the even more convenient reduced bispectrum which removes the geometric factors associated with the Gaunt integral From the previous two formulae we also derive the following useful relations between the full, averaged, and reduced bispectra:

The reduced bispectrum from (34) then takes the much simpler form This is the key equation in this section, since it explicitly relates the primordial bispectrum, predicted by inflationary theories, to the reduced bispectrum observed in the cosmic microwave sky. This formula is entirely analogous to the well-known relation linking the primordial curvature power spectrum and the CMB angular power spectrum that is,

Finally, it is important to note that the Gaunt integral in (41) encodes several constraints on the angle-averaged bispectrum which are no longer transparent in the reduced bispectrum . These are as follows. (1)The sum of the three multipoles must be even (to ensure parity invariance). (2)The 's satisfy the triangle condition (to enforce rotational invariance).

Analogous to the wavenumber constraint (23), the second condition tells us that the only multipole configurations giving nonzero contributions to the bispectrum are those that form a closed triangle in harmonic (-)space. For wavenumbers, the triangle condition is enforced through the -integral over the three spherical Bessel functions which evaluates to zero if the 's cannot form a triangle, whereas in multipole space it is enforced by the angular integration over the spherical harmonics in (40).

##### 3.2. Separable Primordial Shapes and CMB Bispectrum Solutions

In terms of the shape function (12), the reduced bispectrum (44) can be rewritten as The expression above can be simplified, and simple analytic solutions can sometimes be obtained for the very important class of separable shapes obeying the ansatz , as in (20). Substituting (20) into (46), we find that where we have defined the quantities

Instead of the three-dimensional integral of (46), we now have to deal with a much more tractable product of three one-dimensional integrals. Moreover, if we work at large angular scales in the Sachs-Wolfe approximation, the transfer functions become , where and represents respectively, the present-day conformal time and the conformal time at decoupling. The presence of a product of spherical Bessel functions in the integrals above can lead in some cases to simple analytic solutions.

Let us demonstrate this for the separable primordial shapes considered in Section 2. The simplest possible shape, the constant model (14) with , has a large-angle analytic solution for the reduced bispectrum [9]:

The large-angle solution (49) is an important benchmark with which to compare the shape of late-time CMB bispectra from other models (note the scaling). The more general constant solution does not have an analytic solution because the transfer functions cannot be expressed in a simple form, but it can be evaluated numerically from the expression

The numerical solution is shown in Figure 7, exhibiting a regular pattern of acoustic peaks introduced by the oscillating transfer functions.

**(a)**

**(b)**

For the local shape (13), the Sachs-Wolfe approximation also yields a large-angle analytic solution where the divergences for the squeezed triangles () in the primordial shape (13) are also reflected in It is straightforward, in principle, to calculate the full bispectrum from the separable expressions arising from (13): where the separated integrals analogous to (50) become However, we note that these highly oscillatory integrals must be evaluated numerically with considerable care.

For the equilateral shape (54) we first make its separability explicit by expanding the expression in the form: While there is no simple large-angle analytic solution known for the equilateral model, it can be evaluated from the simplified expression where are given in (53) and , in the scale-invariant case, are defined by (compare with the local case)

Before concluding this section we would like to note how all of the solutions presented here present a characteristic scaling. This is just a direct consequence of the scaling of the primordial bispectrum in the scale-invariant case, and it is a model-independent result. It is analogous to the typical scaling displayed by the angular power spectrum in correspondence to a scale-invariant primordial power spectrum .

##### 3.3. Nonseparable Bispectra Revisited

Recall the mode expansion (21) of a general nonseparable primordial shape. If we substitute this into the expression for the reduced bispectrum (46), then the separability of the expansion leads to the same efficient calculation route discussed in the previous section through [14]: where simply result from convolving the basis functions with the transfer functions The computationally costly 3D integrals have again reduced to a sum over products of 1D integrals; we note that this economy arises because the triangle condition is enforced in (56) through the product of Bessel functions, resulting in a manifestly separable form in which we can interchange orders of integration. With this mode expansion, all nonseparable theoretical CMB bispectra become efficiently calculable provided that there is a convergent expansion for the shape function.

In the same way that we decomposed an arbitrary primordial shape in Section 2.2, it is possible to construct analogous late-time separable basis functions and orthonormal modes with which to describe the CMB bispectrum [6, 9]. The tetrahedral domain defined by the triangle condition for multipole configurations is essentially identical to that for wavenumbers (23), except that only even cases contribute . However, the appropriate weight function now incorporates Wigner- symbols arising from bispectrum products: where in the second expression we have exploited the freedom to divide by a separable function and use a weight which makes the bispectrum functions more scale invariant (eliminating an factor—see below). The inner product between two functions and is altered from the primordial wavenumber integral (59) into a sum over multipoles on the tetrahedral domain; that is, But for the change in the weight (which only affects configurations near the edges of the tetrahedron), the 1D polynomials and the 3D separable product basis functions (), as well as the resulting orthonormal modes , are nearly identical to their primordial counterparts , and defined in Section 2.2.

We can now expand an arbitrary CMB bispectrum in both the separable and orthonormal mode expansions, which is achieved in the following form: where the variance term reflects the signal-to-noise weighting expected in the CMB estimator (see Section 3.5). Again, the coefficients in the expansions are determined, first, from the orthonormal inner products and, secondly, the separable are found with the transformation matrix analogous to (26). Examples of the convergence of these mode expansions for equilateral, DBI, and cosmic string CMB bispectra are given in Figure 5.

##### 3.4. CMB Bispectrum Calculations and Correlations

Prior to the systematic mode expansion approach (56) being implemented, robust hierarchical schemes were developed to calculate any nonseparable CMB bispectrum (46) directly [6, 9]. These use the transverse coordinate system given in (11) and employ adaptive methods on a triangular grid to accurately determine the oscillatory 2D -integrations, with important efficiencies also coming from the flat sky approximation, binning, and interpolation schemes. Precision to greater than 1% across the full Planck domain was established by direct comparison with analytic solutions such as (49) and (51). Examples of nonseparable (and separable) CMB bispectra found using these hierarchical coarse-graining methods are shown in Figures 7 and 8. While the CMB bispectra retain the qualitative features of the primordial shape functions , they are overlaid with the oscillatory transfer functions which give rise to a coherent pattern of acoustic peaks. These direct bispectrum calculations revealed that typical primordial models could be described by eigenmode or other expansions using only a limited number of terms.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

Motivated by the form of the CMB estimator, we can define the following correlator to determine whether or not two competing theoretical bispectra can be distinguished by an ideal experiment: where the normalization is defined as The emergence of the inner product (59) in the expression (61) means that substitution of the mode expansions (60) for the theoretical bispectra reduces the correlator to While the late-time correlator (61) is the best measure of whether two CMB bispectra are truly independent, it can be demonstrated that for the majority of models the shape correlator (19) introduced earlier is sufficient to determine independence.

On the basis of the direct calculation of the bispectrum results and the CMB correlator, we can now quantitatively check the forecasting accuracy of the primordial shape correlator proposed previously (again closely following the discussion in [9]).

###### 3.4.1. Nearly Scale-Invariant Models

For nearly scale-invariant models, the centre values for the bispectrum all have roughly the same profile but with different normalisations. As we see from Figure 8, the oscillatory properties of the transfer functions for the CMB power spectrum create a series of acoustic peaks for any combinations involving the following multipole values: . Of course, to observe the key differences between the scale-invariant models we must study the bispectrum in the plane orthogonal to the -direction, that is, the directions reflecting changes in the primordial shape functions. To plot the bispectrum (see Figures 7 and 8), we consistently divide by the large-angle CMB bispectrum solution for the constant model (14). This is analogous to multiplying the power spectrum 's by , because it serves to remove the overall scaling of the bispectrum, flattening while preserving the transverse momentum-dependence primordial shape, and the effects of the oscillating transfer functions.

The starting point is the constant model (14) which, despite its apparent simplicity, has a CMB bispectrum revealing a nontrivial and coherent pattern of acoustic peaks that we have already noted (see Figure 7). Given that the constant model has no momentum dependence, we stress that the resulting bispectrum is the three-dimensional analogue of the angular power spectrum for a scale-invariant model. The largest (primary) peak, for example, is located where all three (corresponding to the large blue region near the origin). We can interpret Figure 7, therefore, as the pure window function or beam effect of convolving any model with the radiation transfer functions while transforming from Fourier to harmonic space.

The CMB bispectrum for the equilateral model is plotted in Figure 8, showing how the centre weighting from the primordial shape is well preserved despite the convolution with the oscillating transfer functions. For the full CMB correlator (61), the DBI, ghost, and single shapes are generally even more closely correlated with the equilateral model, presumably because distinctive features are “washed out” by the transformation from Fourier to harmonic space. Comparative results between the shape correlator and the Fisher matrix analysis are given in Table 1, establishing that these models are highly correlated and difficult to set apart observationally.

The CMB bispectrum for the local model is also shown in Figure 8, demonstrating a marked contrast with the equilateral model which reflects their different primordial shapes shown in Figure 3. The dominance of the signal in the squeezed limit creates strong parallel ridges of acoustic peaks which connect up and emanate along the corner edges of the tetrahedron (see [58] for further details). The 51% CMB correlation between the local and equilateral models is underestimated by the shape correlator at 41%, presumably because of effective smoothing due to the harmonic analysis. Reflecting their distinctive primordial properties, the CMB bispectra for the flat and warm models are poorly correlated with most of the other models, though the flat shape could be susceptible to confusion with the local CMB bispectrum with which it has a larger correlation (see Table 2). It is clear that the local, equilateral, warm, and flat shapes form four distinguishable categories among the scale-invariant models.

###### 3.4.2. Scale-Dependent Models, Cosmic Strings, and Other Late-Time Phenomena

Models which have a nontrivial scaling, such as the feature models, can have starkly contrasting bispectra as illustrated in Figure 8. For example, instead of having the same pattern of acoustic peaks which characterise the scale-invariant models, the feature model can become entirely anticorrelated so that the primary peak has the opposite sign. Later, for this particular choice of in (29), for increasing the phase of the oscillations becomes positively correlated by the second and third peaks. This can lead to small correlation with the other primordial shapes, all below 45% as shown in Table 2 for this and . Clearly, these nonseparable feature models form a distinct fifth category beyond the four scale-invariant shapes noted above and, of course, there are many possible model dependencies which can lead to further subdivision.

By way of further illustration of the breadth of other possible nonseparable CMB bispectra, we present the late-time CMB bispectrum predicted analytically for cosmic strings [59] as where , , , and Here, is a model-dependent amplitude with measuring the string tension relative to the Planck scale. The cutoffs around in (51) are associated with the string correlation length at decoupling (perturbations with can only be causally seeded after last scattering). Here, the nonseparable nature and very different scaling of the string CMB bispectrum are clear from a comparison with (51). Moreover, given the late-time origin of this signal from string metric perturbations, the modulating effect of acoustic peaks from the transfer functions is absent, as is clear from Figure 8. This is just one example of late-time phenomena such as gravitational lensing, secondary anisotropies, and contaminants which are accessible to analysis using the more general CMB mode expansions (60).

##### 3.5. The Estimation of from CMB Bispectra

In light of the previous discussion, it is evident how measurements of the bispectrum from CMB experimental datasets are able to provide information about the primordial three-point function of the cosmological curvature perturbation field at the end of inflation. This in turn allows us to put significant constraints on inflationary models or on alternative models for the generation of cosmological perturbations. We will now start dealing with the problem of bispectrum estimation in the CMB as a test of primordial non-Gaussianity.

Let us assume that we have measured the three-point function of a given CMB dataset. There are now two general ways to exploit this information.(1)*Tests of the Gaussian Hypothesis*. By comparing the measured three-point function to its expected distribution obtained from Gaussian simulation we can detect whether some configurations present a significant deviation from Gaussian expectations. The issue with this approach is that it is sensitive not only to primordial non-Gaussianity, but also to any other possible source of NG, including those of noncosmological origin. Original bispectrum tests of this kind on COBE maps [60] revealed significant deviations from Gaussianity in the data. This NG signature in the three-point function seemed to be localized in harmonic space around multipoles and was object of much scrutiny (see, e.g., [61–65]). It was then finally ascertained that the detected signal was not cosmological in origin, but due to a systematic artifact [66]. Moreover, the overall statistical significance of the result disappeared in a later analysis involving the measurement of all of the bispectrum modes available in the map [67] (only a subset of all the configurations had been studied before). General tests of Gaussianity are very useful to identify unexpected effects in the data, and to monitor systematics. However, as long as we are interested in a *primordial* NG signal, it is better to follow the approach of making an ansatz for the bispectrum we expect from the theory under study and obtain a *quantitative* constraint on a given model. This approach is outlined in the point that follows.(2)* Estimation*. In this case we choose the primordial model that we want to test, characterizing it through its bispectrum shape. We then estimate the corresponding amplitude from the data. If the final estimate is consistent wih , then we conclude that no significant detection of the given shape is produced by the data, but we still determine important constraints on the allowed range of . Note that ideally we would like to do more than just constrain the overall amplitude, and reconstruct the entire shape from the data by measuring single configurations of the bispectrum. However, the expected primordial signal is too small to allow the signal from a single-bispectrum triangle to emerge over the noise. For this reason we study the cumulative signal from all of the configurations that are sensitive to .

Since in this review we are concerned with the study of the primordial bispectrum, we will take the latter approach and deal with the problem of estimation from measurements of the bispectrum in CMB maps. We will first present a cubic estimator that optimally extracts the information from the data contained in the bispectrum (Section 3.5.1). We will then address the issue of understanding whether this optimal cubic statistic extracts all of the possible information available on * in the data* or whether there is enough additional information beyond the three-point function to allow more precise measurements using non-bispectrum-based estimators of (Section 3.5.2). We will then discuss concrete numerical implementations of bispectrum estimators (Section 3.6) and review the experimental constraints on obtained from bispectrum analysis of WMAP data (Section 3.7). Using a standard Fisher matrix analysis, forecasts on the error bars are achievable for future CMB surveys (Section 3.8). Following, we will study the NG signals in the map that could contaminate the primordial NG measurement and how they are dealt with when analyzing the data. Finally we will describe algorithms for the simulation of primordial NG CMB maps that are useful for testing and validation of estimators before applying them to real data.

In the following, we assume that the reader is familiar with essential concepts in statistical estimation theory, such as the definition of a statistical estimator, the role played by maximum-likelihood estimators in statistics, the definitions of unbiasedness and optimality, and the definition and main applications of the Fisher information matrix. The reader unfamiliar with these concepts can consult the appendix of this review and references therein.

###### 3.5.1. Bispectrum Estimator of

In this section we are concerned with the statistical inference of from measurements of the bispectrum of the CMB anisotropies. We recall that we defined earlier as the amplitude of the bispectrum of the primordial potential. In principle, we can include both temperature and polarization multipoles in the analysis, in order to maximize the available data. However, for clarity we will consider only temperature multipoles in the following and omit the superscript in , for simplicity of notation. The extension to polarization is conceptually straightforward and will be discussed in a following paragraph. We will start by considering a simple cubic (The estimator is dubbed cubic due to the fact that it contains the third power of the random variable ) statistic written in the form In the previous equation represents the statistical estimate of from the data, are the multipoles of the observed CMB temperature fluctuations, are some weight functions, and is a normalization factor that has to be chosen to make the estimator * unbiased,* that is, to ensure that

We now want to find the weights that provide the best estimator (i.e., the minimum error bar estimator)* within the class of cubic statistics* written in the form of (66). It is a well-known result (see Appendix) that the best unbiased estimator of a parameter from a given dataset is the maximum-likelihood estimator. In order to answer our question we then have to write the bispectrum likelihood as a function of the parameter and maximize with respect to .

In the assumption that the bispectrum configurations are characterized by a Gaussian distribution (This is not strictly true, but it is a good approximation. The same approach applies to most cosmological observables.), maximizing the likelihood is equivalent to minimizing the following : where is the observed angular averaged bispectrum, that is, by definition and is the bispectrum variance, that is, the six-point function We will now make the assumption that we are working in the *weak non-Gaussian limit; that is, * is small and the distribution of can be approximated as Gaussian in the calculation of the variance. The implications of this approximation will be discussed in greater detail in the following sections; for the moment it will suffice to point out that the weak non-Gaussian approximation is generally a good one since most inflationary models predict to be small, and because the level of primordial non-Gaussianity is already constrained to be small by WMAP measurements [1, 68]. After restricting indices so that and , the six-point function above can be calculated using Wick's theorem, yielding [57] In the last formula is a permutation factor that takes the value of when all 's are different, when two 's are equal, and when all 's are equal. We can now substitute (71) into (69) and differentiate with respect to to get an explicit expression for the optimal cubic statistic we were looking for: where is the reduced bispectrum and is the Gaunt integral defined by equation (40); is the normalization factor mentioned at the beginning of the paragraph that guarantees the unbiasedness of the estimator.

Note that the noise and window function of the experiment are included in the and that appear in the formula above, with the following replacements: where is the window function (not to be confused with the weights ) and is the noise power spectrum (constant for uncorrelated white noise). The noise is assumed to be Gaussian, thus characterized by a vanishing three-point function. Comparing our result (72) to the initial ansatz (66), we then see that the optimal weights are In other words we are weighting the observed bispectrum by its expected signal-to-noise ratio.

We have now constructed a statistic that optimally extracts the information about from the bispectrum of the map. The question now is the following: is there additional information about in the map that is not contained in the bispectrum? This issue will be investigated in the following sections. For the impatient reader we anticipate that the answer is no: the bispectrum statistic built here is actually the minimum error bar estimator of from CMB data.

###### 3.5.2. Optimality of the Cubic Estimator

In this section we address the issue of whether the cubic statistic (72) optimally extracts* all* of the information contained in or whether other statistical estimators (e.g., four-point function, or pixel space statistics such as the Minkowski functionals, or again wavelet estimators, just to mention a few among many possible examples) are able to produce smaller error bars and are thus more efficient than the bispectrum.

In a non-Gaussian primordial CMB map, the likelihood depends on the NG parameter . We will indicate it with , where indicates a vector including all of the 's (we will assume that all other cosmological parameters are fixed and concentrate on ). It is a well-known result in parameter estimation theory that there is a *lower limit* on the error bars that can be assigned to a given parameter (in our case ). Such lower limit, also known as the * Rao-Cramer bound*, is defined in terms of the *Fisher matrix * as (We again refer the reader unfamiliar with these concepts to the brief summary provided in appendix.) We remind the reader that the Fisher matrix is defined as If we can show that the bispectrum estimator of the previous section saturates the Rao-Cramer bound for the Fisher matrix above, then we conclude that it provides the best (i.e., minimum variance) estimate of *from the data*, rather than just the best estimate * from the bispectrum* of the data. In other words, no more information about could be extracted from than the information contained in the bispectrum. The aim of this section is to show that this is actually the case.

The issue of the optimality of bispectrum estimators of was addressed in great detail by Babich in [69]. In this section we will basically review the main results of that study, referring the reader to the original paper for their complete derivation.

As we mention in appendix, there is a sufficient and necessary condition for an estimator to saturate the Rao-Cramer bound, expressed by formula (A.5). This condition, applied to our case, reads Our aim is to show that the bispectrum statistic (72) satisfies this condition. We then need to start from a computation of the full likelihood for a general primordial non-Gaussian model. Following Babich [69], we will start by limiting ourselves to the particular case of the local model.

We recall from a previous section that local NG is the only case for which an explicit expression for the primordial potential is provided. In real space, Starting from this formula, it is possible to obtain a likelihood function for , dependent on the parameter . This is done by means of an* expansion* in terms of the order parameter . The full expression for the Probability Density Function (PDF) (see [69]) can be expanded around its Gaussian expectation for and schematically written as where is the covariance matrix of the Gaussian part of the potential that is,

Formula (80) is then telling us that the logarithm of the full likelihood can be decomposed into the sum of a Gaussian likelihood , plus a NG term that depends linearly on and that this decomposition is accurate up to terms of order ; *that is, we are assuming that NG is weak*, as we did in the previous section.

After computing , one has to account for 2D projection and radiative transfer in order to obtain the required likelihood . As shown by Babich in [69], this can be achieved by expanding the PDF (80) in spherical harmonics and performing the functional integration: where is the Dirac delta function of dimension and due to the 2D projection (As noted by Babich in [69], the additional degrees of freedom do not affect the CMB anisotropies and can therefore be integrated out.) The previous formula can be derived by recalling (32) together with the well-known formula in probability theory as follows: where is again the Dirac delta function, and are random variables linked by the functional relation , and , are the PDFs of and respectively. Solving the functional integral (82) yields [69]

In the previous formula we can recognize the standard PDF, valid in the standard Gaussian case, in the first term on the r.h.s. Added to this, we find a first-order -correction proportional to the CMB angular bispectrum. Higher-order correlators are not present at order . For reasons that will become clear shortly, although we have not computed it, we have explicitly denoted the term in the expansion with . Note that, besides assuming weak NG in this formula, we are also assuming * rotational invariance* (this is evident from the fact that the covariance matrix appearing in the Gaussian piece of (84) is diagonal and equal to ). Rotational invariance is a general property of the CMB sky, but it is broken when we deal with real CMB measurement characterized by inhomogeneous noise patterns and sky cuts. We will investigate these effects in the following section. For the moment we consider the purely ideal case described by (84). Armed with the PDF expression for local NG and recalling the necessary and sufficient condition (78), we can finally determine whether the estimator (72) is optimal or not. First of all we see that

We then see from combining (85) and (72) that We now see that, in order for the necessary and sufficient condition for optimality (78) to be verified, we need the term to be exactly equal to . The second-order quantity should then be calculated explicitly in the expansion (84) in order to complete the calculation and verify whether, or under which conditions, this is true. However this turns out not to be necessary if we consider the following “regularity condition for a PDF”. (Condition (87) can be easily derived remembering that, for a given random variable with probability density , we have by definition , and substituting in the previous expression, one then finds that the regularity condition (87) holds, provided the order of integration and differentiation can be exchanged (hence the “regularity condition” qualification). ) For a general PDF of a random variable depending on a parameter we have

Since this regularity condition must be valid for each value of the parameter ( in our case), it is clear that it must hold * term-by-term*, that is, at each order, in the expansion (84). By taking the average value of equation (86), keeping in mind that the estimator is unbiased and imposing (87), we then find that the average value of must be exactly equal to . If we could then replace in (86) with its average value, then we would exactly obtain the condition for optimality and conclude that the cubic estimator (72) saturates the Rao-Cramer bound. For present CMB experiments, the terms in the expansion (84) are evaluated summing over a large number of -modes ( for WMAP, for Planck in the signal-dominated regime), or, equivalently in pixel space, averaging over a large number of pixels (10^{6} and for WMAP and Planck, resp.). For this reason we expect that the error made by replacing the -order term in (86) with its average value will be very small. In [69], an estimate of this error has been done in the approximation of neglecting radiative transfer and projection effects (i.e., working in 3D with the primordial potential, rather than with a CMB map). The conclusion was that for a number of observations the approximation above works very well. Moreover the variance of the -order term scales like . In the full radiative transfer case we expect the scaling to be unchanged, although the coefficients in front of it that led to the estimate might change. However, as noted above, the number of pixels in present-day experiments is many orders of magnitude larger than . That leads us to conclude that the approximation of replacing the average of the first-order-term in equation (86) is a very good one. We then reach the following important conclusion.

* For a rotational invariant CMB sky, in the limit of weak NG, the cubic estimator defined by formula* (72)* is the best unbiased CMB estimator of *

Let us now move to problem of generalizing the last conclusion to shapes different from local. In this case a full expression of the primordial potential is not available. The steps that lead to the conclusion that the local estimator is optimal can thus not be reproduced. However it was pointed out by Babich in [69] that, in the limit of weak NG, the full CMB NG likelihood can still be expressed in terms of its power spectrum and bispectrum by means of an Edgeworth expansion, regardless of its full expression. The Edgeworth expansion is basically a way to express a NG PDF as a series expansion around its Gaussian part [70–72]. For CMB anisotropies one finds, at the end of the calculation, that