© ESO, 2016

1. Introduction

This paper, one of a set associated with the 2015 release of data from the Planck¹ mission (Planck Collaboration I 2016), describes the constraints on primordial non-Gaussianity (NG) obtained using the cosmic microwave background (CMB) maps from the full Planck mission, including a first analysis of some of the Planck polarization data.

Primordial NG is one of the most powerful tests of inflation, and more generally of high-energy early Universe physics (for some reviews, see Bartolo et al. 2004a; Liguori et al. 2010; Chen 2010b; Komatsu 2010; Yadav & Wandelt 2010). In fact, the simplest models of inflation (characterized by a single scalar field slowly rolling along a smooth potential) predict the generation of primordial fluctuations that are almost Gaussian distributed, with a tiny deviation from Gaussianity of the order of the slow-roll parameters (Acquaviva et al. 2003; Maldacena 2003). The 2013 Planck results on primordial NG are consistent with such a prediction, being compatible with Gaussian primordial fluctuations: the standard scenario of single-field, slow-roll inflation has survived its most stringent test to date. For example, in Planck Collaboration XXIV (2014) we obtained ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{2.7}\mathrm{\pm }\mathrm{5.8}$, ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{-}\mathrm{42}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{75}$, and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}\mathrm{=}\mathrm{-}\mathrm{25}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{39}$ for the amplitudes of three of the most well-studied shapes of primordial NG. On the other hand, it is well known that any deviations from the standard picture of inflation have the potential to produce distinctive NG signatures at a detectable level in the CMB anisotropies². Therefore, as already shown in Planck Collaboration XXIV (2014; see also Planck Collaboration XXII 2014) improved NG constraints allow severe limits to be placed on various classes of inflationary models that extend the simplest paradigm, in a way that is strongly complementary to the power-spectrum constraints (i.e., scalar spectral index of curvature perturbations and tensor-to-scalar amplitude ratio).

One of the main goals of this paper is to improve NG constraints using mainly the angular bispectrum of CMB anisotropies, i.e., the harmonic transform of the 3-point angular correlation function. We also investigate higher-order NG correlators like the trispectrum. We follow the same notation as Planck Collaboration XXIV (2014). The CMB angular bispectrum is related to the primordial bispectrum $\mathrm{⟨}\mathrm{\Phi }\mathrm{\left(}{{k}}_{\mathrm{1}}\mathrm{\right)}\mathrm{\Phi }\mathrm{\left(}{{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{\Phi }\mathrm{\left(}{{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{⟩}\mathrm{=}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}{\mathit{\delta }}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}\mathrm{\right(}{{k}}_{\mathrm{1}}\mathrm{+}{{k}}_{\mathrm{2}}\mathrm{+}{{k}}_{\mathrm{3}}\mathrm{\left)}{\mathit{B}}_{\mathrm{\Phi }}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{)}\mathit{,}$(1)where the field Φ, related to the comoving curvature perturbation ζ on super-horizon scales by Φ ≡ (3 / 5)ζ, is such that in the matter era, and on super-horizon scales, it reduces to Bardeen’s gauge-invariant gravitational potential (Bardeen 1980). The bispectrum B_Φ(k₁,k₂,k₃) measures the correlation among three perturbation modes. If translational and rotational invariance are assumed, it depends only on the magnitude of the three wavevectors. In general the bispectrum can be written as $\begin{array}{ccc}{\mathit{B}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{NL}}\mathit{F}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathit{,}& & \end{array}$(2)where we have introduced the dimensionless “nonlinearity parameter” f_NL (Gangui et al. 1994; Wang & Kamionkowski 2000; Komatsu & Spergel 2001; Babich et al. 2004), measuring the NG amplitude. The bispectrum is obtained by sampling triangles in Fourier space. The dependence of the function F(k₁,k₂,k₃) on the type of triangle (i.e., the configuration) formed by the three wavevectors describes the shape (and the scale dependence) of the bispectrum (Babich et al. 2004), which encodes much physical information. Different NG shapes are linked to distinctive physical mechanisms that can generate such NG fingerprints in the early Universe.

In this paper the limits on primordial NG are mainly improved through the use of the full mission data, as well as by exploiting the polarization information.

Planck results on primordial NG also provide a reconstruction of the full CMB bispectrum through different techniques (see Sect. 6.2). This complements (and adds to) the extraction of single amplitudes f_NL for specific bispectrum shapes. Such a reconstruction can point to interesting features in the bispectrum signal that go beyond the usual standard scale-invariant shapes (such as the well known “local” and “equilateral” configurations).

As we have seen, the Planck 2013 NG paper (Planck Collaboration XXIV 2014) significantly improved constraints on the standard primordial NG models with scale-invariant local, equilateral or orthogonal shapes. The Planck NG paper also included constraints from the modal estimator on a variety of other primordial models, including DBI inflation, non-Bunch-Davies models (excited initial states), directionally-dependent vector inflation models, warm inflation, and scale-dependent feature and resonance models. All scale-invariant bispectra were strongly constrained, with the possible exception of highly flattened non-Bunch-Davies models. On the other hand, the preliminary investigation of primordial oscillatory models seemed to be more promising, in that two specific feature models appeared to produce fits of some significance. One aim of the present work is to expand the detail and scope of investigations of feature and resonant models and to examine the significance of these results with a more careful analysis of the “look elsewhere” effect, through exploring multi-parameter results using large ensembles of Gaussian simulations. Also we will thoroughly analyse or re-analyse other primordial NG signals that are theoretically well-motivated and those which have appeared in the literature since the first data release. These include primordial NG arising in the context of inflation models where vector fields play a non-negligible role or primordial NG generated in the tensor (gravitational waves) perturbations. Each of these primordial NG signals carry distinctive signatures that may have been imprinted at the inflationary epoch, thus opening up a new window into the detailed physics of inflation.

The paper is organized as follows. In Sect. 2 we briefly discuss the primordial NG models that we test in this paper. Section 3 summarizes the optimal statistical estimators used to constrain the CMB bispectrum and trispectrum from Planck temperature and polarization data. In Sect. 4 we discuss the non-primordial contributions to the CMB bispectrum and trispectrum, including foreground residuals after component separation and focusing on the f_NL bias induced by the ISW-lensing bispectrum. We also analyse the impact on primordial NG estimation from the residuals of the deglitching processing. Section 5 describes an extensive suite of tests performed on realistic simulations to validate the different estimator pipelines, and compare their performance. Using simulations, we also quantify the impact on f_NL of using a variety of component-separation techniques. In Sect. 6 we derive constraints on f_NL for the local, equilateral, and orthogonal bispectra and present a reconstruction of the CMB bispectrum. We also present a reconstruction of the primordial curvature fluctuations. In Sect. 7 we validate these results by performing a series of null tests on the data to assess the robustness of our results. Section 8 investigates scale-dependent NG models and other selected bispectrum shapes. Section 9 presents the Planck limits on the CMB trispectrum. In Sect. 10 we provide constraints on CMB local bispectrum and trispectrum from Minkowski functionals. In Sect. 11 we discuss the main implications of Planck’s constraints on primordial NG for early Universe models. We conclude in Sect. 12. Appendix A contains the derivation of the statistical estimator for the amplitudes characterizing a “direction-dependent” primordial non-Gaussianity. Appendix B contains some details about Minkowski functionals.

2. Models

In this section we briefly highlight the classes of inflationary models investigated in this paper, and describe the distinctive NG they generate. Within each class a common underlying physical process gives rise to the corresponding NG shape, illustrated by concrete realizations of inflationary models. For each class we therefore provide the explicit form of the bispectrum shapes chosen for the data analysis, emphasizing extensions with variants and distinctly new shapes beyond those already described in Planck Collaboration XXIV (2014).

2.1. General single-field models of inflation

This class of models includes inflationary models with a non-standard kinetic term (or more general higher-derivative interactions), in which the inflaton fluctuations propagate with an effective sound speed c_s which can be smaller than the speed of light. For example, models with a non-standard kinetic term are described by an inflaton Lagrangian ℒ = P(X,φ), where X = g^μν∂_μφ,∂_νφ, with at most one derivative on φ, and the sound speed is given by ${\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}\mathrm{(}\mathit{\partial P}\mathit{/}\mathit{\partial X}\mathrm{)}\mathit{/}\mathrm{(}\mathit{\partial P}\mathit{/}\mathit{\partial X}\mathrm{+}\mathrm{2}\mathit{X}\mathrm{(}{\mathit{\partial }}^{\mathrm{2}}\mathit{P}\mathit{/}\mathit{\partial }{\mathit{X}}^{\mathrm{2}}\mathrm{\left)}\mathrm{\right)}$.

The NG parameter space of this class of models is generically well described by two NG shapes – “equilateral” and “orthogonal” (Senatore et al. 2010) – since usually there are two dominant interaction terms of the inflaton field giving rise to the overall NG signal. One of these typically produces a bispectrum very close to the equilateral type with ${\mathit{f}}_{\mathrm{NL}}\mathrm{~}{\mathit{c}}_{\mathrm{s}}^{-2}$ in the limit c_s ≪ 1 (Chen et al. 2007b; Senatore et al. 2010).

The equilateral-type NG is well approximated by the template (Creminelli et al. 2006) $\begin{array}{ccc}& & {\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{equil}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\mathrm{6}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\\ & & \mathrm{\times }\{\mathrm{-}\frac{\mathrm{1}}{{\mathit{k}}_{\mathrm{1}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}{\mathit{k}}_{\mathrm{2}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{k}}_{\mathrm{2}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}{\mathit{k}}_{\mathrm{3}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}}\mathrm{-}\frac{\mathrm{1}}{{\mathit{k}}_{\mathrm{3}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}{\mathit{k}}_{\mathrm{1}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}}\mathrm{-}\frac{\mathrm{2}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}\mathrm{(}\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}\mathrm{)}\mathit{/}\mathrm{3}}}\\ & & \end{array}$(3)Here P_Φ(k) = A/k^{4 − n_s} represents Bardeen’s gravitational potential power spectrum, A² being the normalization and n_s the scalar spectral index. DBI inflationary models based on string theory (Silverstein & Tong 2004; Alishahiha et al. 2004) provide physically well-motivated examples of the P(X,φ)-model. They are characterized by an almost equilateral NG with ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{35}\mathit{/}\mathrm{108}\mathrm{)}{\mathit{c}}_{\mathrm{s}}^{-2}$ for c_s ≪ 1, which typically is ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathit{<}\mathrm{-}\mathrm{5}$.

The “orthogonal” shape template is (Senatore et al. 2010) $\begin{array}{ccc}& & {\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{ortho}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\mathrm{6}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}\\ & & \mathrm{\times }\{\mathrm{-}\frac{\mathrm{3}}{{\mathit{k}}_{\mathrm{1}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}{\mathit{k}}_{\mathrm{2}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}}\mathrm{-}\frac{\mathrm{3}}{{\mathit{k}}_{\mathrm{2}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}{\mathit{k}}_{\mathrm{3}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}}\mathrm{-}\frac{\mathrm{3}}{{\mathit{k}}_{\mathrm{3}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}{\mathit{k}}_{\mathrm{1}}^{\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}}}\mathrm{-}\frac{\mathrm{8}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}\mathrm{(}\mathrm{4}\mathrm{-}{\mathit{n}}_{\mathrm{s}}\mathrm{)}\mathit{/}\mathrm{3}}}\\ & & \end{array}$(4)Equilateral and orthogonal shapes emerge also in models characterized by more general higher-derivative interactions, such as ghost inflation (Arkani-Hamed et al. 2004), effective field theories of inflation (Cheung et al. 2008; Senatore et al. 2010; Bartolo et al. 2010a), or the so “Galileon-like” models of inflation (see, e.g., Burrage et al. 2011). The latter model is constructed starting from some specific underlying symmetry for the inflaton field, and is characterized by strongly constrained derivative interactions.

2.2. Multi-field models

This class of models is characterized by the presence of additional light scalar degrees of freedom besides the inflaton, whose fluctuations give rise, or contribute, to the final primordial curvature perturbation at the end of inflation. This includes the case of “multiple-field inflation”, where inflation is driven by more than one scalar field, as well as scenarios in which additional scalar fields remain subdominant during the inflationary expansion. From the point of view of primordial NG, the element in common to all these models is that a potentially detectable level of NG in the curvature perturbation is generated via a transfer of super-horizon non-Gaussian isocurvature perturbations in the second field (not necessarily the inflaton) to the adiabatic (curvature) density perturbations, accompanied by nonlinearities in the transfer mechanism. This process typically takes place on super-horizon scales, thus implying a local form of NG in real space. When going to Fourier space, this leads to a correlation between large and small scale modes. The bispectrum for this class of models is indeed largest on so-called “squeezed” triangles (k₁ ≪ k₂ ≃ k₃). The local bispectrum is (Falk et al. 1993; Gangui et al. 1994; Gangui & Martin 2000; Verde et al. 2000; Wang & Kamionkowski 2000; Komatsu & Spergel 2001) $\begin{array}{ccc}{\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{local}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& \mathrm{2}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}[{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}\mathrm{+}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\\ & & \mathrm{+}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}]\\ & \mathrm{=}& \end{array}$(5)There is a broad literature on examples and specific realizations of this transfer mechanism from isocurvature to adiabatic perturbations (Bartolo et al. 2002; Bernardeau & Uzan 2002; Vernizzi & Wands 2006; Rigopoulos et al. 2006, 2007; Lyth & Rodriguez 2005; Tzavara & van Tent 2011; for a review on NG from multiple-field inflation models, see Byrnes & Choi 2010). An alternative, important possibility is the curvaton model (Mollerach 1990; Linde & Mukhanov 1997; Enqvist & Sloth 2002; Lyth & Wands 2002; Moroi & Takahashi 2001). In this type of scenario, a second light scalar field, subdominant during inflation, decays after inflation, generating primordial density perturbations with a potentially high level of NG (e.g., Lyth & Wands 2002; Lyth et al. 2003; Bartolo et al. 2004c). In the (simplest) adiabatic curvaton models, the local f_NL parameter was found to be (Bartolo et al. 2004c,b) ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{(}\mathrm{5}\mathit{/}\mathrm{4}{\mathit{r}}_{\mathrm{D}}\mathrm{)}\mathrm{-}\mathrm{5}{\mathit{r}}_{\mathrm{D}}\mathit{/}\mathrm{6}\mathrm{-}\mathrm{5}\mathit{/}\mathrm{3}$, when the curvaton field has a quadratic potential (Lyth & Wands 2002; Lyth et al. 2003; Lyth & Rodriguez 2005; Malik & Lyth 2006; Sasaki et al. 2006). In the previous formula, r_D = [ 3ρ_curvaton/ (3ρ_curvaton + 4ρ_radiation) ] _D is the “curvaton decay fraction” evaluated at the epoch of the curvaton decay in the sudden decay approximation. It is then easy to see that, for low values of r_D, a high level of NG can be generated³.

2.3. Isocurvature non-Gaussianity

Isocurvature NG, which was only sketched from the purely theoretical point of view in the 2013 paper, can now be analysed thanks to the polarization information.

In most of the models mentioned above, the main focus is on the level of primordial NG in the final curvature perturbation ζ. However, in inflationary scenarios where different scalar fields play a non-negligible role, residual isocurvature perturbation modes can remain after inflation. Isocurvature modes are usually investigated by considering their contribution to the power spectrum. However, if present, they would also contribute to the bispectrum, producing in general both a pure isocurvature bispectrum and mixed bispectra because of the cross-correlation between isocurvature and adiabatic perturbations (Komatsu 2002; Bartolo et al. 2002; Komatsu et al. 2005; Kawasaki et al. 2008, 2009; Langlois et al. 2008; Hikage et al. 2009; Langlois & Lepidi 2011; Langlois & van Tent 2011, 2012; Kawakami et al. 2012). While one might expect isocurvature NG to be negligible, since both (linear) isocurvature modes and (adiabatic) NG appear to be very small, and searches for isocurvature NG using WMAP data did not lead to any detections (Hikage et al. 2013a,b), this expectation can be tested at significantly higher precision by Planck. Moreover, there exist inflation models (Langlois & Lepidi 2011) where isocurvature modes, while remaining a small fraction in the power spectrum, would dominate the bispectrum.

At the time of recombination there are in principle four possible distinct isocurvature modes (in addition to the adiabatic mode): cold dark matter (CDM); baryon; neutrino density; and neutrino velocity isocurvature modes (Bucher et al. 2000). In this paper we will only consider isocurvature NG of the local type and always limit ourselves to considering the adiabatic mode together with just one type of isocurvature mode (considering each of the four types separately). Otherwise the number of free parameters becomes so large that no meaningful limits can be derived. Moreover, we assume the same spectral index for the primordial isocurvature power spectrum and the isocurvature-adiabatic cross-power spectrum as for the adiabatic power spectrum. Under those assumptions, as shown by Langlois & van Tent (2011), we have in principle six independent f_NL parameters: the usual purely adiabatic one; a purely isocurvature one; and four correlated ones.

The primordial shape templates are a generalization of Eq. (5), see Langlois & van Tent (2011, 2012): $\begin{array}{ccc}{\mathit{B}}^{\mathit{IJK}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& \mathrm{2}{\mathit{f}}_{\mathrm{NL}}^{\mathit{I,JK}}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}\mathrm{2}{\mathit{f}}_{\mathrm{NL}}^{\mathit{J,KI}}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\\ & & \end{array}$(6)where I,J,K label the different modes (adiabatic and isocurvature). The invariance under the simultaneous exchange of two of these indices and the corresponding momenta means that ${\mathit{f}}_{\mathrm{NL}}^{\mathit{I,JK}}\mathrm{=}{\mathit{f}}_{\mathrm{NL}}^{\mathit{I,KJ}}$, hence reducing the number of independent parameters from eight to six, in the case of two modes. The different bispectra vary most importantly through the fact that different types of radiation transfer functions ${\mathit{g}}_{\mathit{\ell }}^{\mathit{I}}\mathrm{\left(}\mathit{k}\mathrm{\right)}$ are used to project the primordial template onto the CMB: the reduced bispectra are of the form $\begin{array}{ccc}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathit{I,JK}}\mathrm{=}\mathrm{6}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{\infty }}{\mathit{r}}^{\mathrm{2}}\mathrm{d}\mathit{r,}{\mathit{\alpha }}_{\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{I}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{\beta }}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{J}}\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{\beta }}_{{\mathit{\ell }}_{\mathrm{3}}\mathrm{)}}^{\mathit{K}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}& & \end{array}$(7)with $\begin{array}{ccc}{\mathit{\alpha }}_{\mathit{\ell }}^{\mathit{I}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{\equiv }& \frac{\mathrm{2}}{\mathit{\pi }}\mathrm{\int }{\mathit{k}}^{\mathrm{2}}\mathrm{d}\mathit{k,}{\mathit{j}}_{\mathit{\ell }}\mathrm{\left(}\mathit{kr}\mathrm{\right)}\mathit{,}{\mathit{g}}_{\mathit{\ell }}^{\mathit{I}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}\\ {\mathit{\beta }}_{\mathit{\ell }}^{\mathit{I}}\mathrm{\left(}\mathit{r}\mathrm{\right)}& \mathrm{\equiv }& \frac{\mathrm{2}}{\mathit{\pi }}\mathrm{\int }{\mathit{k}}^{\mathrm{2}}\mathrm{d}\mathit{k,}{\mathit{j}}_{\mathit{\ell }}\mathrm{\left(}\mathit{kr}\mathrm{\right)}\mathit{,}{\mathit{g}}_{\mathit{\ell }}^{\mathit{I}}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{,}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\left(}\mathit{k}\mathrm{\right)}\mathit{.}\end{array}$Here j_ℓ is the spherical bessel function and we use the notation (ℓ₁ℓ₂ℓ₃) ≡ [ ℓ₁ℓ₂ℓ₃ + 5perm. ] / 3 !. In addition to the isocurvature index, each transfer function carries a polarization index that we do not show here. It is important to note that, unlike the case of the purely adiabatic mode, the inclusion of polarization improves the constraints on the isocurvature NG significantly, as predicted by Langlois & van Tent (2011, 2012).

2.4. Resonance and axion monodromy models

Oscillatory models for NG are physically well-motivated. Large-field inflation faces an inherent UV completion problem because the inflaton field is required to move over large distances in field space relative to the Planck mass m_Pl. An effective shift symmetry can enforce potential flatness and this can be naturally implemented in a string theory context with axions and a periodically modulated potential, so-called “axion monodromy” models. This periodicity can generate resonances in the inflationary fluctuations with logarithmically-spaced oscillations, creating imprints in the power spectrum, the bispectrum and trispectrum (Chen et al. 2008; Flauger et al. 2010; Hannestad et al. 2010; Flauger & Pajer 2011). On the other hand, sharp features or corners in an inflationary potential can temporarily drive the inflaton away from slow-roll; these large changes in the field and derivatives can create evenly-spaced oscillations, to be discussed in the next subsection. However, in multifield models residual oscillations after corner-turning can also lead to log-spaced oscillations, just as in the resonance models (Chen 2011; Achúcarro et al. 2011; Battefeld et al. 2013; Chen et al. 2015). A preliminary search for bispectrum resonance signals was performed in the first Planck analysis (Planck Collaboration XXIV 2014) and our purpose here is to substantially increase the frequency range and number of models investigated.

Simple resonance model: periodic features in the inflationary potential can induce oscillations with frequency ω that can resonate through any interactions with the inflationary fluctuations, contributing to the bispectrum. Provided that ω>H, this mode starts inside the horizon but its frequency decreases as it is stretched by inflation, until frozen when ω ≃ H. Thus periodic features introduce a driving force which can scan across a wide range of frequencies. The simplest basic behaviour of such resonant models yields logarithmic stretching and can be described by the non-scale-invariant shape (see Chen et al. 2008; Chen 2010b) $\begin{array}{ccc}{\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{res}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{6}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{res}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\sin \left[\mathit{C}\ln \mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}\mathit{\phi }\right]\mathit{,}& & \end{array}$(10)where the constant C = 1 / ln(3k_∗), k_∗ is a wavenumber associated with the periodicity, and φ is a phase. These oscillations constructively and destructively interfere with the oscillations created by the CMB transfer functions, introducing additional nodal points in the CMB bispectrum.

Generalized resonance models: in a more general context, it is possible to have more complicated resonant shapes and envelopes. Resonant single-field models with varying sound speed c_s generate three leading-order bispectrum terms (Chen 2010a): $\begin{array}{ccc}& & {\mathit{B}}^{\mathrm{res}\mathrm{-}\mathrm{cs}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{6}{\mathit{A}}^{\mathrm{2}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\{{\mathit{f}}_{\mathrm{NL}}^{\mathrm{res}\mathrm{1}}\sin \left[\mathit{C}\ln \mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}\mathit{\phi }\right]\\ & & \mathrm{+}\mathrm{3}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{res}\mathrm{2}}\frac{{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}^{\mathrm{2}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\cos \left[\mathit{C}\ln \mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}\mathit{\phi }\right]\\ & & \mathrm{+}\mathrm{27}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{res}\mathrm{3}}\frac{{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{3}}}\sin \left[\mathit{C}\ln \mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}\mathit{\phi }\right]\}\mathit{.}\end{array}$(11)The first term on the right-hand side of Eq. (11) is the basic resonant shape given in Eq. (10), while the second and third terms have the same oscillatory behaviour, but modulated by a (mildly) flattened shape, and an equilateral shape respectively. The third term is in fact the second generic shape arising in effective field theory and correlates well with the equilateral shape in Eq. (12). The second term in Eq. (11) weakly favours flattened triangles, but there are regimes for resonant models that can generate much stronger flat shapes. If the resonance begins very deep inside the horizon, then the second (negative energy) mode can also make a significant contribution that is associated with enfolded or flat bispectra; this is similar to having an excited initial state or non-Bunch-Davies (NBD) vacuum.

With these two physical motivations in mind we also investigate classes of models with resonant oscillations modulated by both the equilateral and flattened shapes, defined by $\begin{array}{ccc}{\mathit{S}}^{\mathrm{eq}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{1}}\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{2}}\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{3}}\end{array}}{{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}}\mathit{,}{\mathit{S}}^{\mathrm{flat}}\mathrm{=}\mathrm{1}\mathrm{-}{\mathit{S}}^{\mathrm{eq}}\mathit{,}& & \end{array}$(12)where $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{k}}_{\mathrm{1}}\end{array}\mathrm{\equiv }{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{-}{\mathit{k}}_{\mathrm{1}}$ (here, for simplicity we ignore the spectral index dependence of the equilateral shape in Eq. (3)). The corresponding equilateral and flattened resonant bispectra ansätze are then $\begin{array}{ccc}{\mathit{B}}^{\mathrm{res}\mathrm{-}\mathrm{eq}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{\equiv }& {\mathit{S}}^{\mathrm{eq}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\times }{\mathit{B}}^{\mathrm{res}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\\ & & \\ {\mathit{B}}^{\mathrm{res}\mathrm{-}\mathrm{flat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{\equiv }& \end{array}$We note that typically non-Bunch-Davies bispectra can be much more sharply peaked in the flattened or squeezed limits than Eq. (14), but our purpose here is to determine if this type of resonant model is favoured by the Planck data; that is, whether Eq. (14) warrants further investigation with other flattened profiles.

2.5. Scale-dependent oscillatory feature models

Temporary violations of slow-roll inflation can occur if there are sharp features in the inflationary potential (Chen et al. 2007a), as well as changes in the sound speed c_s or sharp turns in field space in multifield inflation. The inflaton field makes temporary departures from the attractor solution, which typically have a strong scale-dependent running modulated by a sinusoidal oscillation; there are model-dependent counterparts in the power spectrum, bispectrum, and trispectrum. For example, sharper or narrower features induce a relatively larger signal in the bispectrum (see e.g., Chen 2010b). An example is the analytic envelope solutions predicted for both the power spectrum and bispectrum for the single field models with a specific inflaton feature shape (Adshead et al. 2012); a search for these was presented previously in the Planck Inflation paper (Planck Collaboration XXII 2014) and likewise no significant signal was found using the corresponding bispectrum envelopes at the available modal resolution (Planck Collaboration XXIV 2014). In this new analysis, we will emphasize the search for generic oscillatory behaviour in the data over a larger range in modal resolution, although we will also look for the shapes predicted for simple features in single field models.

Constant feature model: In the previous investigation of Planck data using a coarse parameter grid (Planck Collaboration XXIV 2014), we searched for the simplest ansatz for an oscillatory bispectrum signal (Chen et al. 2007a): $\begin{array}{ccc}{\mathit{B}}^{\mathrm{feat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{6}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{feat}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\sin \left[\mathit{\omega }\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{)}\mathrm{+}\mathit{\phi }\right]\mathit{,}& & \end{array}$(15)where φ is a phase factor and ω is a frequency associated with the specific shape of the feature in the potential that disrupts the slow-roll evolution. In the earlier analysis, we also considered a damping envelope, which slightly increased the apparent significance of the best-fit feature models, though at the cost of an additional parameter (see single-field solutions below).

Generalized feature models: here, we again search for oscillatory signals in a model-independent manner. We will modulate the bispectrum cross-sections with the physically motivated equilateral and flattened shapes, reflecting the physical contexts in which they could have been generated, as for the resonant models described above in Eq. (11). If there are potential features in a model with a varying sound speed, then we can expect there to be oscillatory contributions to the bispectrum signal with a dominant equilateral shape. Motivated by the equilateral resonance model in Eq. (11), we will search for the following equilateral feature ansatz: $\begin{array}{ccc}{\mathit{B}}^{\mathrm{feat}\mathrm{-}\mathrm{eq}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{\equiv }& {\mathit{S}}^{\mathrm{eq}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\times }{\mathit{B}}^{\mathrm{feat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\\ & \mathrm{=}& \frac{\mathrm{6}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{feat}\mathrm{-}\mathrm{eq}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\frac{\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{1}}\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{2}}\end{array}\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{3}}\end{array}}{{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}}\sin \left[\mathit{\omega }\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{)}\mathrm{+}\mathit{\phi }\right]\mathit{.}\end{array}$For extremely sharp features, it is possible to excite the inflationary fluctuations as if there were a non-Bunch Davies vacuum: the oscillatory signal becomes modulated with a flattened shape (Chen et al. 2007a). Again, motivated by the enfolded resonance model in Eq. (14), we take the following simple flattened ansatz: $\begin{array}{ccc}{\mathit{B}}^{\mathrm{feat}\mathrm{-}\mathrm{flat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\equiv }{\mathit{S}}^{\mathrm{flat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{\times }{\mathit{B}}^{\mathrm{feat}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathit{.}& & \end{array}$(18)Although the exact profile of the flattened shape can be much more highly peaked on the faces in these NBD models, this ansatz should be adequate for testing whether these models are favoured. We note that while the power spectrum is insensitive to the underlying scenario creating the features, the bispectrum shape will reveal whether features arise from varying sound speed or highly excited features in the potential.

Single field feature solutions: here we use the full analytic bispectrum solution given by Adshead et al. (2012), but the dominant leading-order behaviour takes the form $\begin{array}{ccc}{\mathit{B}}^{{\mathit{K}}^{\mathrm{2}}\cos }\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{6}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{{\mathrm{K}}^{\mathrm{2}}\cos }}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}{\mathit{K}}^{\mathrm{2}}\mathit{D}\mathrm{\left(}\mathit{\alpha \omega K}\mathrm{\right)}\cos \mathrm{\left(}\mathit{\omega K}\mathrm{\right)}\mathit{,}& & \end{array}$(19)where K = k₁ + k₂ + k₃ and D(αωK) = αω/ (Ksinh(αωK)) is an envelope function, with parameter α setting an overall cut-off for the bispectrum at large wavenumbers or multipoles. This envelope and the overall K² scaling distinguishes this realistic case from the simple separable constant feature ansatz of Eq. (15). We shall allow the envelope parameter α to vary from α = 0, with no envelope (the infinitely thin limit for a feature in the potential) through to large α, with a narrow domain for the bispectrum. Alternative analytic solutions where the bispectrum is created by a variation in the sound speed c_s are dominated by the Ksin(ωK) term, as in $\begin{array}{ccc}{\mathit{B}}^{\mathit{K}\sin }\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{6}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{K}\sin }}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\mathit{K,D}\mathrm{\left(}\mathit{\alpha \omega K}\mathrm{\right)}\sin \mathrm{\left(}\mathit{\omega K}\mathrm{\right)}\mathit{.}& & \end{array}$(20)For the simplest models there is a predicted relationship between the power spectrum and bispectrum amplitude (e.g., see also Achucarro et al. 2013 for a two-field model). We note that typically the power spectrum has larger signal-to-noise at low frequency (i.e., below ω ≃ 1000) while the bispectrum dominates at higher frequency.

2.6. Non-Gaussianity from excited initial states

It is well known that if the initial vacuum state for inflation is excited and deviates from the standard Bunch-Davies vacuum, then measurable non-Gaussianities can be produced (Chen et al. 2007b; Holman & Tolley 2008; Meerburg et al. 2009; Ashoorioon & Shiu 2011). These models generically lead to non-Gaussianity that peaks in the flattened limit, where k₁ + k₂ ≈ k₃, and also often has oscillatory behaviour. Here we constrain the same selection of templates found in the 2013 Planck analysis, namely the flat model in Eq. (12), Non-Bunch-Davies (NBD; Chen et al. 2007b), NBD1 and NBD2 models (Agullo & Parker 2011; now called “NBD1 cos” and “NBD2 cos”) and NBD3 (Chen 2010b). We also introduce three new templates, NBD sin which is motivated by (Chen 2010a) and takes the form $\begin{array}{ccc}{\mathit{B}}^{\mathrm{NBD}\mathrm{-}\sin }\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{)}& \mathrm{=}& \frac{\mathrm{2}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{NBD}\mathrm{-}\sin }}{{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}\mathrm{)}}^{\mathrm{2}}}\left({\mathrm{e}}^{\mathrm{-}\mathit{\omega }\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{1}}\end{array}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}\mathit{\omega }\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{2}}\end{array}}\mathrm{+}{\mathrm{e}}^{\mathrm{-}\mathit{\omega }\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{3}}\end{array}}\right)\\ & & \mathrm{\times }\sin \left(\mathit{\omega K}\mathrm{+}\mathit{\phi }\right)\mathit{,}\end{array}$(21)where again K = k₁ + k₂ + k₃ and $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{k}}_{\mathit{i}}\end{array}\mathrm{=}\mathit{K}\mathrm{-}\mathrm{2}{\mathit{k}}_{\mathit{i}}$. The other two templates are extensions of the NBD1 cos and NBD2 cos models found in Agullo & Parker (2011) and take the form $\begin{array}{ccc}{\mathit{B}}^{\mathrm{NBD}\mathit{i}\mathrm{-}\sin }\mathrm{(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{)}& \mathrm{=}& \frac{\mathrm{2}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{NBD}\mathit{i}\mathrm{-}\sin }}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}[{\mathit{f}}_{\mathit{i}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{;}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\\ & & \mathrm{\times }\sin \mathrm{\left(}\mathit{\omega }\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{1}}\end{array}\mathrm{\right)}\mathit{/}\begin{array}{c}\mathrm{˜}\\ {\mathit{k}}_{\mathrm{1}}\end{array}\mathrm{+}\mathrm{2}\mathrm{perm}\mathrm{.}]\mathit{,}\end{array}$(22)where ${\mathit{f}}_{\mathrm{1}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{;}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{(}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{2}$, which is dominated by squeezed configurations, and ${\mathit{f}}_{\mathrm{2}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{;}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}^{\mathrm{2}}$, which has a flattened shape.

2.7. Directional-dependence motivated by gauge fields

Some models where primordial vector fields are present during inflation predict interesting NG signatures. This is the case of a coupling of the inflaton field ϕ to the kinetic term of a gauge field A^μ, ℒ contains − I²(ϕ)F², where F_μν = ∂_μA_ν − ∂_νA_μ and the coupling I²(ϕ)F² is chosen so that scale invariant vector perturbations are produced on superhorizon scales (Barnaby et al. 2012b; Bartolo et al. 2013a). The bispectrum turns out to be the sum of two contributions: one of the local shape; and another that is also enhanced in the squeezed limit (k₁ ≪ k₂ ≃ k₃), but featuring a non-trivial dependence on the angle between the small and the large wave vectors through the parameter ${\mathit{\mu }}_{\mathrm{12}}\mathrm{=}\hat{{k}}{}_{\mathrm{1}}\mathrm{·}\hat{{k}}{}_{\mathrm{2}}$ (where $\hat{{k}}\mathrm{=}{k}\mathit{/}\mathit{k}$) as ${\mathit{\mu }}_{\mathrm{12}}^{\mathrm{2}}$. Also, primordial magnetic fields sourcing curvature perturbations can cause a dependence on both μ and μ² (Shiraishi et al. 2012).

We can parametrize these shapes as variations on the local shape (Shiraishi et al. 2013a), as $\begin{array}{ccc}{\mathit{B}}_{\mathrm{\Phi }}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{L}}{\mathit{c}}_{\mathit{L}}\mathrm{\left[}{\mathit{P}}_{\mathit{L}}\mathrm{\right(}{\mathit{\mu }}_{\mathrm{12}}\mathrm{\left)}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\right(}{\mathit{k}}_{\mathrm{1}}\mathrm{\left)}{\mathit{P}}_{\mathrm{\Phi }}\mathrm{\right(}{\mathit{k}}_{\mathrm{2}}\mathrm{)}\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{perm}\mathrm{.}\mathrm{]}\mathit{,}& & \end{array}$(23)where P_L(μ) is the Legendre polynomial with P₀ = 1,P₁ = μ, and ${\mathit{P}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}\mathrm{3}{\mathit{\mu }}^{\mathrm{2}}\mathrm{-}\mathrm{1}\mathrm{)}$. For example, for L = 1 we have the shape ${\mathit{B}}_{\mathrm{\Phi }}^{\mathit{L}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{2}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathit{L}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}}{\mathrm{(}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathrm{)}}^{\mathrm{2}}}\mathit{,}\left[\frac{{\mathit{k}}_{\mathrm{3}}^{\mathrm{2}}}{{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}}\mathrm{(}{\mathit{k}}_{\mathrm{1}}^{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}^{\mathrm{2}}\mathrm{-}{{\mathit{k}}_{\mathrm{3}}^{\mathrm{2}}}^{\mathrm{)}}\mathrm{+}\mathrm{2}\hspace{0.17em}\mathrm{perm}\mathrm{.}\right]\mathit{.}$(24)The local template corresponds to c_i = 2f_NLδ_i0. Here and in the following the nonlinearity parameters ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}}$ are related to the c_L coefficients by ${\mathit{c}}_{\mathrm{0}}\mathrm{=}\mathrm{2}{\mathit{f}}_{\mathrm{NL}}^{\mathit{L}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{0}}$, ${\mathit{c}}_{\mathrm{1}}\mathrm{=}\mathrm{-}\mathrm{4}{\mathit{f}}_{\mathrm{NL}}^{\mathit{L}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}$, and ${\mathit{c}}_{\mathrm{2}}\mathrm{=}\mathrm{-}\mathrm{16}{\mathit{f}}_{\mathrm{NL}}^{\mathit{L}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{2}}$. The L = 1,2 shapes exhibit sharp variations in the flattened limit, for example for k₁ + k₂ ≈ k₃, while in the squeezed limit, L = 1 is suppressed whereas L = 2 grows like the local bispectrum shape (i.e., the L = 0 case). The I²(ϕ)F² models predict c₂ = c₀/ 2, while primordial curvature perturbations sourced by large-scale magnetic fields generate non-vanishing c₀, c₁, and c₂. Quite interestingly, in the proposed “solid inflation” scenario (Endlich et al. 2013, 2014; see also Bartolo et al. 2013b, 2014; Sitwell & Sigurdson 2014) bispectra similar to Eq. (23) can be generated, in this case with c₂ ≫ c₀ (Endlich et al. 2013, 2014). Therefore, measurements of the c_i coefficients can be an efficient probe of some detailed aspects of the inflationary mechanism, such as the existence of primordial vector fields during inflation (or a non-trivial symmetry structure of the inflaton fields, as in solid inflation).

2.8. Non-Gaussianity from gauge-field production during axion inflation

The same shift symmetry that leads to axion (monodromy) models of inflation (Sect. 2.4) naturally allows (from an effective field theory point of view) for a coupling between a pseudoscalar axion inflaton field and a gauge field of the type $\mathrm{ℒ}\mathrm{\supset }\mathrm{-}\mathrm{(}\mathit{\alpha }\mathit{/}\mathrm{4}\mathit{f}\mathrm{)}\mathit{\phi }{\mathit{F}}^{\mathit{\mu \nu }}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{F}}_{\mathit{\mu \nu }}\end{array}$, where the parameter α is dimensionless and f is the axion decay constant ($\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{F}}^{\mathit{\mu \nu }}\end{array}\mathrm{=}{\mathit{ϵ}}^{\mathit{\mu \nu \gamma \beta }}{\mathit{F}}_{\mathit{\gamma \beta }}\mathit{/}\mathrm{2}$). This scenario has a rich and interesting phenomenology both for scalar and tensor primordial fluctuations (see, e.g., Barnaby & Peloso 2011; Sorbo 2011; Barnaby et al. 2011, 2012c; Linde et al. 2013; Meerburg & Pajer 2013; Ferreira & Sloth 2014). Gauge field quanta are produced by the background motion of the inflaton field, and these in turn source curvature perturbations through an inverse decay process of the gauge field. A bispectrum of curvature fluctuations is generated as (Barnaby et al. 2011; Meerburg & Pajer 2013)⁴$\begin{array}{ccc}{\mathit{B}}^{\mathrm{inv}\mathit{.}\mathrm{dec}}\mathrm{=}\mathrm{6}{\mathit{A}}^{\mathrm{2}}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{inv}\mathit{.}\mathrm{dec}}\frac{{\sum }_{\mathit{i}}{\mathit{k}}_{\mathit{i}}^{\mathrm{3}}}{\mathrm{􏽑}{\mathit{k}}_{\mathit{i}}^{\mathrm{3}}}\frac{{\mathit{f}}_{\mathrm{3}}\mathrm{(}{\mathit{\xi }}_{\mathrm{\ast }}\mathit{,}\frac{{\mathit{k}}_{\mathrm{2}}}{{\mathit{k}}_{\mathrm{1}}}\mathit{,}{\frac{{\mathit{k}}_{\mathrm{3}}}{{\mathit{k}}_{\mathrm{1}}}}^{\mathrm{)}}}{{\mathit{f}}_{\mathrm{3}}\mathrm{\left(}{\mathit{\xi }}_{\mathrm{\ast }}\mathit{,}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{\right)}}\mathit{,}& & \end{array}$(25)where the exact expression for the function f₃ can be found in Eq. (3.29) of Barnaby et al. (2011; see also Meerburg & Pajer 2013). Here ξ characterizes the coupling strength of the axion to the gauge field . The inverse decay bispectrum peaks for equilateral configuration, since δϕ is mostly sourced by the inverse decay (δA + δA → δϕ), when two modes of the vector fields are of comparable magnitude (the correlation with the equilateral template is 94% and with the orthogonal one is 4%). We do however constrain the exact shape in Eq. (25), without resorting to the equilateral template. Another interesting observational signature that can shed light on the role played by pseudo-scalars in the early Universe is provided by tensor NG, to which we turn next.

2.9. Parity-violating tensor non-Gaussianity motivated by pseudo-scalars

While the majority of the studies on primordial and CMB NG focus on the scalar mode, tensor-mode NG has been attracting attention as a probe of high-energy theories of gravity (e.g., Maldacena & Pimentel 2011; McFadden & Skenderis 2011; Soda et al. 2011; Shiraishi et al. 2011; Gao et al. 2011) or primordial magnetic fields (Brown & Crittenden 2005; Shiraishi et al. 2012; Shiraishi 2012)⁵.

Recently, the possibility of observable tensor bispectra has been vigorously discussed in a model where the inflaton couples to a pseudoscalar field (Barnaby et al. 2012a; Cook & Sorbo 2013; Ferreira & Sloth 2014). In this model, through the gravitational coupling to the U(1) gauge field, gravitational waves (${\mathit{h}}_{\mathit{ij}}\mathrm{\equiv }\mathit{\delta }{\mathit{g}}_{\mathit{ij}}^{\mathit{TT}}\mathit{/}{\mathit{a}}^{\mathrm{2}}\mathrm{=}{\sum }_{\mathit{s}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{\pm }}{\mathit{h}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}{\mathit{e}}_{\mathit{ij}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}$) receive NG corrections, where only one of the two spin states is enhanced. The bispectrum, generally formed as $\begin{array}{ccc}⟨\underset{\mathit{i}\mathrm{=}\mathrm{1}}{\overset{\mathrm{3}}{\mathrm{􏽙}}}{\mathit{h}}^{\mathrm{\left(}{\mathit{s}}_{\mathit{i}}\mathrm{\right)}}\mathrm{\left(}{{k}}_{\mathit{i}}\mathrm{\right)}⟩\mathrm{=}\mathrm{\left(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}{\mathit{\delta }}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}\mathrm{\right(}{{k}}_{\mathrm{1}}\mathrm{+}{{k}}_{\mathrm{2}}\mathrm{+}{{k}}_{\mathrm{3}}\mathrm{\left)}{\mathit{B}}_{\mathrm{h}}^{{\mathit{s}}_{\mathrm{1}}{\mathit{s}}_{\mathrm{2}}{\mathit{s}}_{\mathrm{3}}}\mathrm{\right(}{{k}}_{\mathrm{1}}\mathit{,}{{k}}_{\mathrm{2}}\mathit{,}{{k}}_{\mathrm{3}}\mathrm{)}\mathit{,}& & \end{array}$(26)is accordingly polarized, with ${\mathit{B}}_{\mathrm{h}}^{\mathrm{+}\mathrm{+}\mathrm{+}}\mathrm{\gg }{\mathit{B}}_{\mathrm{h}}^{\mathrm{+}\mathrm{+}\mathrm{-}}\mathit{,}{\mathit{B}}_{\mathrm{h}}^{\mathrm{+}\mathrm{-}\mathrm{-}}\mathit{,}{\mathit{B}}_{\mathrm{h}}^{\mathrm{-}\mathrm{-}}$. This NG enhancement is a sub-horizon effect and therefore ${\mathit{B}}_{\mathrm{h}}^{\mathrm{+}\mathrm{+}\mathrm{+}}$ is maximized at the equilateral limit (k₁ ≃ k₂ ≃ k₃) (Cook & Sorbo 2013).

A model-independent template of the equilateral-type polarized tensor bispectrum is given by (Shiraishi et al. 2013b, 2015) $\begin{array}{ccc}{\mathit{B}}_{\mathrm{h}}^{\mathrm{+}\mathrm{+}\mathrm{+}}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathit{,}{{k}}_{\mathrm{2}}\mathit{,}{{k}}_{\mathrm{3}}\mathrm{\right)}& \mathrm{=}& {\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}{\mathit{F}}_{\mathit{\zeta }}^{\mathrm{equil}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\\ & & \end{array}$(27)with the polarization tensor ${\mathit{e}}_{\mathit{ij}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}$ obeying ${\mathit{e}}_{\mathit{ij}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{\left(}\hat{{k}}\mathrm{\right)}{\mathit{e}}_{\mathit{ij}}^{\mathrm{\left(}{\mathit{s}}^{\mathrm{\prime }}\mathrm{\right)}}\mathrm{\left(}\hat{{k}}\mathrm{\right)}\mathrm{=}\mathrm{2}{\mathit{\delta }}_{\mathit{s,}\mathrm{-}{\mathit{s}}^{\mathrm{\prime }}}$ and ${\mathit{e}}_{\mathit{ij}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}\mathrm{\ast }}\mathrm{\left(}\hat{{k}}\mathrm{\right)}\mathrm{=}{\mathit{e}}_{\mathit{ij}}^{\mathrm{(}\mathrm{-}\mathit{s}\mathrm{)}}\mathrm{\left(}\hat{{k}}\mathrm{\right)}\mathrm{=}{\mathit{e}}_{\mathit{ij}}^{\mathrm{\left(}\mathit{s}\mathrm{\right)}}\mathrm{(}\mathrm{-}\hat{{k}}\mathrm{)}$. We here have introduced a tensor nonlinearity parameter, by normalizing with the equilateral bispectrum template of curvature perturbations (${\mathit{F}}_{\mathit{\zeta }}^{\mathrm{equil}}\mathrm{\equiv }\mathrm{(}\mathrm{5}\mathit{/}\mathrm{3}{\mathrm{)}}^{\mathrm{3}}{\mathit{F}}_{\mathrm{\Phi }}^{\mathrm{equil}}\mathrm{=}\mathrm{(}\mathrm{5}\mathit{/}\mathrm{3}{\mathrm{)}}^{\mathrm{3}}{\mathit{B}}_{\mathrm{\Phi }}^{\mathrm{equil}}\mathit{/}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}$) in the equilateral limit, yielding $\begin{array}{ccc}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}\mathrm{\equiv }\underset{{\mathit{k}}_{\mathit{i}}\mathrm{\to }\mathit{k}}{\lim }\frac{{\mathit{B}}_{\mathrm{h}}^{\mathrm{+}\mathrm{+}\mathrm{+}}\mathrm{\left(}{{k}}_{\mathrm{1}}\mathit{,}{{k}}_{\mathrm{2}}\mathit{,}{{k}}_{\mathrm{3}}\mathrm{\right)}}{{\mathit{F}}_{\mathit{\zeta }}^{\mathrm{equil}}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathit{,}{\mathit{k}}_{\mathrm{2}}\mathit{,}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}}\mathrm{·}& & \end{array}$(28)The template Eq. (27) can adequately reconstruct the tensor bispectra created in the pseudoscalar inflation models⁶ (Shiraishi et al. 2013b), and thus the amplitude ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}$ is directly connected with the model parameters, e.g., the coupling strength of the pseudoscalar field to the gauge field ξ (for details see Sect. 11).

The CMB temperature and E-mode bispectra sourced by the parity-violating tensor NG have not only the usual parity-even (ℓ₁ + ℓ₂ + ℓ₃ = even) signals but also parity-odd (ℓ₁ + ℓ₂ + ℓ₃ = odd) contributions, which cannot be sourced by known scalar bispectra (Kamionkowski & Souradeep 2011; Shiraishi et al. 2011). Moreover, their shapes are mostly distinct from the scalar templates, due the different radiation transfer functions; hence they can be measured essentially independently of the scalar NG (Shiraishi et al. 2013b). The analysis of the WMAP temperature data distributed in ℓ₁ + ℓ₂ + ℓ₃ = odd configurations leads to an observational limit ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}\mathrm{=}\mathrm{(}\mathrm{0.8}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{1.1}\mathrm{)}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{4}}$ (Shiraishi et al. 2015). This paper updates the limit, by analysing both parity-even and parity-odd signals in the Planck temperature and E-mode polarization data.

3. Statistical estimation of the CMB bispectrum for polarized maps

We now provide a brief overview of the main statistical techniques that we use to estimate the nonlinearity parameter f_NL from temperature and polarization CMB data, followed by a description of the data set that will be used in our analysis.

The CMB temperature and polarization fields are characterized using the multipoles of a spherical harmonic decomposition of the CMB maps: $\begin{array}{ccc}& & \frac{\mathrm{\Delta }\mathit{T}}{\mathit{T}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell m}}{\mathit{a}}_{\mathit{\ell m}}^{\mathit{T}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{,}\\ & & \mathit{E}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell m}}{\mathit{a}}_{\mathit{\ell m}}^{\mathit{E}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{.}\end{array}$(29)At linear order, the relation between the primordial perturbation field and the CMB multipoles is (e.g., Komatsu & Spergel 2001) $\begin{array}{ccc}{\mathit{a}}_{\mathit{\ell m}}^{\mathit{X}}\mathrm{=}\mathrm{4}\mathit{\pi }\mathrm{(}\mathrm{-}\mathit{i}{\mathrm{)}}^{\mathit{\ell }}\mathrm{\int }\frac{{\mathrm{d}}^{\mathrm{3}}\mathit{k}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}}\hspace{0.17em}\mathrm{\Phi }\mathrm{(}{k}\mathrm{\left)}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\right(}\hat{{k}}\mathrm{\left)}{\mathrm{\Delta }}_{\mathit{\ell }}^{\mathit{X}}\mathrm{\right(}\mathit{k}\mathrm{)}\mathit{,}& & \end{array}$(30)where X = { T,E } denotes either temperature or E-mode polarization, Φ is the primordial gravitational potential, and ${\mathrm{\Delta }}_{\mathit{\ell }}^{\mathit{X}}$ represents the linear CMB radiation transfer function.

The CMB angular bispectrum is the three-point correlator of the a_ℓms ${\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}\mathit{,}{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}\mathrm{\equiv }\mathrm{⟨}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}}\mathrm{⟩}\mathit{,}$(31)where X_i = { T,E }. If the CMB sky is rotationally invariant, and the bispectra we are considering have even parity (which is true for combinations of T and E), then the angular bispectrum can be factorized as $\mathrm{⟨}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}}\mathrm{⟩}\mathrm{=}{\mathrm{𝒢}}_{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}^{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\hspace{0.17em}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}\mathit{,}$(32)where ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}$ is the so-called reduced bispectrum, and ${\mathrm{𝒢}}_{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}^{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}$ is the Gaunt integral, defined as the integral over the solid angle of the product of three spherical harmonics, ${\mathrm{𝒢}}_{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}^{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{\equiv }\mathrm{\int }{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\hspace{0.17em}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\hspace{0.17em}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\hspace{0.17em}{\mathrm{d}}^{\mathrm{2}}\hat{{n}}\mathit{.}$(33)The Gaunt integral (often written in terms of Wigner 3j-symbols) enforces rotational symmetry, and restricts attention to a tetrahedral domain of multipole triplets { ℓ₁,ℓ₂,ℓ₃ }, satisfying both a triangle condition and a limit given by some maximum resolution ℓ_max (the latter being defined by the finite angular resolution of the experiment under study).

Our goal is to extract the nonlinearity parameter f_NL from the data, for different primordial shapes. To achieve this, we essentially fit a theoretical CMB bispectrum ansatz b_ℓ1ℓ2ℓ3 to the observed 3-point function. Theoretical predictions for CMB angular bispectra arising from early Universe primordial models can be obtained by applying Eq. (30) to the primordial bispectra of Sect. 2, (see e.g., Komatsu & Spergel 2001). Optimized cubic bispectrum estimators were introduced by Heavens (1998), and it has been shown that for small NG the general optimal polarized f_NL estimator can be written as (Creminelli et al. 2006) $\begin{array}{ccc}\mathit{f̂}\mathrm{NL}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{N}}\sum _{{\mathit{X}}_{\mathit{i}}\mathit{,}{\mathit{X}}_{\mathit{i}}^{\mathrm{\prime }}}\sum _{{\mathit{\ell }}_{\mathit{i}}\mathit{,}{\mathit{m}}_{\mathit{i}}}\sum _{{\mathit{\ell }}_{\mathit{i}}^{\mathrm{\prime }}\mathit{,}{\mathit{m}}_{\mathit{i}}^{\mathrm{\prime }}}{\mathrm{𝒢}}_{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}^{\hspace{0.17em}\hspace{0.17em}{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}\mathit{,}\hspace{0.17em}\mathrm{th}}\\ & & \mathrm{\times }\{[{\left({\mathrm{C}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{1}}^{\mathrm{\prime }}}^{-1}\right)}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{1}}^{\mathrm{\prime }}}^{{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}}\hspace{0.17em}{\left({\mathrm{C}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{2}}^{\mathrm{\prime }}}^{-1}\right)}^{{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{2}}^{\mathrm{\prime }}}^{{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}}\\ & & \mathrm{\times }{\left({\mathrm{C}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{3}}^{\mathrm{\prime }}}^{-1}\right)}^{{\mathit{X}}_{\mathrm{3}}{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}^{\mathrm{\prime }}{\mathit{m}}_{\mathrm{3}}^{\mathrm{\prime }}}^{{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}}]\\ & & \end{array}$(34)where N is a suitable normalization chosen to produce unit response to ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{th}}$. Note that we are implicitly defining a suitable normalization convention so that ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}{\mathit{f}}_{\mathrm{NL}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{th}}$, and ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{th}}$ is the value of the theoretical template when f_NL = 1. C^-1 is the inverse of the block matrix: $\mathrm{C}\mathrm{=}\left(\begin{array}{c}\\ {\mathrm{C}}^{\mathit{TT}}& {\mathrm{C}}^{\mathit{TE}}\\ {\mathrm{C}}^{\mathit{ET}}& {\mathrm{C}}^{\mathit{EE}}\end{array}\right)\mathit{,}$(35)and the blocks represent the full TT, TE, and EE covariance matrices, with C^ET being the transpose of C^TE. All quantities in the previous equation (i.e., CMB multipoles, bispectrum template and covariances matrices) are assumed to properly incorporate instrumental beam and noise.

As standard for these estimators, we note in square brackets (below) the presence of two contributions. One is cubic in the observed a_ℓms, and correlates the bispectrum of the data to the theoretical fitting template ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{th}}$. This is generally called the “cubic term” of the estimator. The other contribution is linear in the observed a_ℓms (“linear term”). This part corrects for mean-field contributions to the error bars, introduced by rotational invariance-breaking features, such as a mask or anisotropic/correlated instrumental noise (Creminelli et al. 2006; Yadav et al. 2008).

The inverse covariance filtering operation implied by Eq. (34) is a challenging numerical task, which has been successfully performed only recently (Smith et al. 2009; Elsner & Wandelt 2012). This step can be avoided by working in the “diagonal covariance approximation”. In this approach, the estimator is built by neglecting off-diagonal entries of the covariance matrix in the cubic term in Eq. (34), and then finding the linear term that minimizes the variance for this specific cubic statistic. Applying such a procedure yields (Yadav et al. 2007) $\begin{array}{ccc}\mathit{f̂}\mathrm{NL}& \mathrm{=}& \frac{\mathrm{1}}{\mathit{N}}\sum _{{\mathit{X}}_{\mathit{i}}\mathit{,}{\mathit{X}}_{\mathit{i}}^{\mathrm{\prime }}}\sum _{{\mathit{\ell }}_{\mathit{i}}\mathit{,}{\mathit{m}}_{\mathit{i}}}{\mathrm{𝒢}}_{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}^{\hspace{0.17em}\hspace{0.17em}{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}{\mathrm{(}{\mathrm{C}}^{-1}\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}}{\mathrm{(}{\mathrm{C}}^{-1}\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathrm{(}{\mathrm{C}}^{-1}\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}\mathit{,}\hspace{0.17em}\mathrm{th}}\\ & & \mathrm{\times }[{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}}\mathrm{-}{\mathrm{C}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}}\mathrm{-}{\mathrm{C}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}}\\ & & \end{array}$(36)where ${\mathrm{C}}_{\mathit{\ell }}^{-1}$ is the inverse of the 2 × 2 matrix ${\mathrm{C}}_{\mathit{\ell }}\mathrm{=}\left(\begin{array}{c}\\ {\mathrm{C}}_{\mathit{\ell }}^{\mathit{TT}}& {\mathrm{C}}_{\mathit{\ell }}^{\mathit{TE}}\\ {\mathrm{C}}_{\mathit{\ell }}^{\mathit{ET}}& {\mathrm{C}}_{\mathit{\ell }}^{\mathit{EE}}\end{array}\right)\mathit{.}$(37)This expression can also be written as $\mathit{f̂}\mathrm{NL}\mathrm{=}\frac{\mathrm{⟨}{\mathit{b}}^{\mathrm{th}}\mathit{,}{\mathit{b}}^{\mathrm{obs}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{b}}^{\mathrm{th}}\mathit{,}{\mathit{b}}^{\mathrm{th}}\mathrm{⟩}}\mathit{,}$(38)where the observed (reduced) bispectrum includes the linear correction term and the inner product is defined as $\begin{array}{ccc}& & \mathrm{⟨}{\mathit{b}}^{\mathit{A}}\mathit{,}{\mathit{b}}^{\mathit{B}}\mathrm{⟩}\mathrm{=}\\ & & \sum _{{\mathit{X}}_{\mathit{i}}\mathit{,}{\mathit{X}}_{\mathit{i}}^{\mathrm{\prime }}}\sum _{{\mathit{\ell }}_{\mathit{i}}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}\mathit{,}\hspace{0.17em}\mathit{A}}{\mathit{h}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{2}}\mathrm{\left(}{\mathrm{C}}^{-1}{\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}}\mathrm{\right(}{\mathrm{C}}^{-1}{\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}}\mathrm{(}{\mathrm{C}}^{-1}{\mathrm{)}}_{{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}\mathit{,}\hspace{0.17em}\mathit{B}}\end{array}$(39)with ${\mathit{h}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}\sqrt{\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{3}}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& {\mathit{\ell }}_{\mathrm{3}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{,}$(40)where the last term is the Wigner 3j symbol. The denominator in Eq. (38), ⟨ b^th,b^th ⟩ is the normalization constant N.

The price to pay for the simplification obtained in Eq. (36) is, in principle, loss of optimality. However, in practice we found in our previous temperature analysis (Planck Collaboration XXIV 2014) that error bars obtained with this simplified procedure are very close to optimal, provided the a_ℓms are pre-filtered with a simple diffusive inpainting technique (see Planck Collaboration XXIV 2014 for details). We find that this still holds true when we include polarization and pre-inpaint the T, Q, U input maps. Given its practical advantages in terms of speed and simplicity, we adopt this method in the following analysis.

A well-known, major issue with both Eqs. (34) and (36) is that their direct implementation would require evaluation of all the bispectrum configurations from the data. The computational cost of this would scale like ${\mathit{\ell }}_{\max }^{\mathrm{5}}$ and be totally prohibitive for high-resolution CMB experiments like Planck. The different bispectrum estimation techniques applied to our analysis are essentially defined by the approach adopted to circumvent this problem. The advantage of having multiple independent implementations of the optimal bispectrum estimator is twofold. First, by cross-validating and comparing outputs of different pipelines, it strongly improves the robustness of the results. Second, different methods are complementary, in the sense that they have specific capabilities which go beyond simple f_NL estimation. For example, the skew-C_ℓ method defined below facilitates the monitoring of NG foreground contamination, while the binned and modal estimators allow model-independent reconstruction of the data bispectrum, and so on. The skew-C_ℓ method enables the nature of any detected NG to be determined. Thus, the simultaneous application of all these techniques also allows us to increase the range and scope of our analysis.

In the following, we briefly outline the main features of the three optimal bispectrum estimation pipelines that are used for Planck measurements of f_NL. We will only provide a short summary here, focused on the extension to polarization data, referring the reader who is interested in more technical aspects to our previous analysis of temperature data (Planck Collaboration XXIV 2014).

3.1. KSW and skew-C_ℓ estimators

KSW and skew-C_ℓ estimators (Komatsu et al. 2005; Munshi & Heavens 2010) can be used for bispectrum templates that are written in factorizable (separable) form, i.e., as a linear combination of separate products of functions⁷ of k₁, k₂, and k₃. This allows reduction of the three-dimensional integration over the bispectrum configurations into a product of three separate one-dimensional sums over ℓ₁, ℓ₂, ℓ₃. This leads to a massive reduction in computational time (O(N_pix), where N_pix is the number of pixels in the map). The main difference between the KSW and skew-C_ℓ pipelines is that the former estimates the f_NL amplitude directly, whereas the latter initially estimates the so called “bispectrum-related power spectrum” (in short, “skew-C_ℓ”) function. Roughly speaking, the skew-C_ℓ associates, with each angular wavenumber ℓ, the contribution to the amplitude f_NL (for each given shape) extracted from all triangles with one fixed side of size ℓ. After resumming over the contributions from each ℓ-bin, the final point-like f_NL estimate is obtained exactly as KSW. Equipping the KSW estimator with a skew-C_ℓ extension can be particularly useful in the presence of (expected) spurious NG contaminants in the data. The slope of the skew-C_ℓ statistic is in fact shape-dependent and can be used to separate multiple NG components in the map.

3.2. Modal estimators

Modal estimators (Fergusson et al. 2010a, 2012) are based on decomposing the bispectrum (both from theory and from data) into a sum of uncorrelated separable templates, forming a complete basis in bispectrum space, and measuring the amplitude of each. The evaluation of the amplitude for each template can be sped up by using a KSW approach (since the templates themselves are separable by construction). All amplitudes form a vector, also referred to as the “mode spectrum”. It is then possible to measure the correlation of the observed data mode spectrum with the theoretical mode spectra for different primordial shapes, in order to obtain estimates of the primordial f_NL. Note also that the observed mode spectrum from data is theory-independent, and contains all the information from the data. Correlating the observed mode spectrum to theoretical mode vectors then allows the extraction of all the f_NL amplitudes simultaneously. This makes modal estimators naturally suited for NG analyses, both when there are a large number of competing models to analyse, or when a model has free parameters through which we wish to scan (more than 500 shapes were analysed when applying this technique to Planck data). Another advantage is that by expanding into separable basis templates, the modal estimator does not require separability of the starting theoretical shape in order to be applicable. Finally, after obtaining the data mode spectrum, it is possible to build a linear combination of the basis templates, using the measured amplitudes as coefficients, thus obtaining a model-independent full reconstruction of the bispectrum of the data. Of course the reconstructed bispectrum will be smoothed, as the estimator must use a finite number of basis templates.

For this analysis, the modal method is implemented in two ways. One of them generalizes our previous temperature modal pipeline by expanding, for each shape, the corresponding TTT, EEE, TTE and EET bispectra. We then exploit separability to build the covariance matrix of these expanded bispectra (Liguori, in prep.), and to measure f_NL efficiently using Eq. (36). This modal pipeline will be referred to throughout the paper as the “Modal 1” pipeline.

The other implementation, which we will refer to as “Modal 2”, utilizes a novel approach where the ${\mathit{a}}_{\mathit{\ell m}}^{\mathit{T}}$ and ${\mathit{a}}_{\mathit{\ell m}}^{\mathit{E}}$ are first orthogonalized to produce new uncorrelated unit variance â_ℓm coefficients, $\begin{array}{ccc}\mathit{â}\begin{array}{c}\mathit{T}\\ \mathit{\ell m}\end{array}& \mathrm{=}& \frac{{\mathit{a}}_{\mathit{\ell m}}^{\mathit{T}}}{\sqrt{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}}}\\ \mathit{â}\begin{array}{c}\mathit{E}\\ \mathit{\ell m}\end{array}& \mathrm{=}& \frac{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}{\mathit{a}}_{\mathit{\ell m}}^{\mathit{E}}\mathrm{-}{\mathit{C}}_{\mathit{\ell }}^{\mathit{TE}}{\mathit{a}}_{\mathit{\ell m}}^{\mathit{T}}}{\sqrt{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}}\sqrt{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}{\mathit{C}}_{\mathit{\ell }}^{\mathit{EE}}\mathrm{-}{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TE}}}^{\mathrm{2}}}}\mathrm{·}\end{array}$We then decompose the new bispectra as $\begin{array}{ccc}\mathit{b̂}\begin{array}{c}{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}\\ {\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}\end{array}\mathrm{=}\sum _{{\mathit{m}}_{\mathit{i}}}\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\hspace{0.17em}\mathrm{⟨}\mathit{â}\begin{array}{c}{\mathit{X}}_{\mathrm{1}}\\ {\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\end{array}\mathit{â}\begin{array}{c}{\mathit{X}}_{\mathrm{2}}\\ {\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}\end{array}\mathit{â}{\begin{array}{c}{\mathit{X}}_{\mathrm{3}}\\ {\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}\end{array}}^{\mathrm{⟩}}\mathit{,}& & \end{array}$(43)which can be constrained independently, since they are uncorrelated. In this case the estimator then takes on a particularly simple form (Fergusson 2014). This new form is mathematically equivalent to the previous modal method, but involves significantly fewer terms in the estimator. However, due to the orthogonalization procedure we cannot constrain the full EEE bispectrum without further processing, just the additional part which is orthogonal to temperature. For this reason, although the “Modal 2” T+E results incorporate all the polarization information, the EEE results alone are not presented here.

In our analysis, both modal techniques (together with all the other estimators described in this section) were used to measure f_NL for the three main shapes i.e., local, equilateral, and orthogonal. Besides this, we optimized the two pipelines for different purposes. The “Modal 1” estimator was adopted to perform a large number of robustness tests of our results, especially in relation to the local, equilateral and orthogonal measurements. The “Modal 2” pipeline was instead mostly used to study a large number of “non-standard” primordial shapes (e.g., oscillatory bispectra). For this reason, each pipeline uses a different set of basis templates. The “Modal 1” estimator starts from a polynomial basis with 600 modes, and includes nine more modes that are the contributions from last scattering to the exact radial KSW expansion of the local, equilateral and orthogonal templates. The “Modal 2” expansion uses a high-resolution basis with 2000 polynomial modes, augmented with a Sachs-Wolfe local bispectrum template, in order to improve convergence efficiency in the squeezed limit. In this way, the high resolution estimator provides the ability to scan across a wide variety of non-separable and oscillatory shapes, while the lower resolution pipeline gives efficient convergence in the f_NL measurements for the standard local, equilateral, and orthogonal shapes, offering rapid analysis for validation purposes. The “Modal 1” pipeline can also be generalized for the estimation of parity-odd bispectra, which is included in our analysis of non-standard shapes.

3.3. Binned bispectrum estimator

One can also use a binned estimator (Bucher et al. 2010, 2015), exploiting the fact that the theoretical bispectra of interest are generally smooth functions in ℓ-space. As a result, data and templates can be binned in ℓ with minimal loss of information, but with large computational gains from data compression. The data bispectrum in the binning grid is then computed and compared to the binned primordial shapes to obtain f_NL. No KSW-like approach, which requires separability and mixing of theoretical and observational bispectra in the computation, is required. Instead, the binned data bispectrum and the binned theoretical bispectrum and covariance are computed and stored completely independently, and only combined at the very last stage in a sum over the bins to obtain f_NL. This means that it is very easy to test additional shapes or different cosmologies, and the data bispectrum can also be studied on its own in a non-parametric approach. In particular the smoothed binned bispectrum approach, also used in this paper, investigates the (smoothed) binned bispectrum of the map divided by its expected standard deviation, to test if there is a significant bispectral NG of any type in the map. Another advantage of the binned bispectrum estimator is that the dependence of f_NL on ℓ can be investigated for free, simply by leaving out bins from the final sum.

In more detail, the computation for the binned bispectrum estimator is based on Eqs. (38) and (39). However, instead of using the reduced bispectrum ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}$, all expressions start from the alternative rotationally-invariant reduced bispectrum ${\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}\mathrm{=}{\mathit{h}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{2}}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}$, where h is defined in Eq. (40). The expression in Eq. (39) for the inner product remains the same when replacing b by B, except that the h² becomes h^-2. The importance of ${\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}$ is that it can be determined directly from maximally-filtered maps; ${\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\hat{{n}}\hspace{0.17em}{\mathit{M}}_{{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\mathit{M}}_{{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\mathit{M}}_{{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{,}$(44)where ${\mathit{M}}_{\mathit{\ell }}^{\mathit{X}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{m}}{\mathit{a}}_{\mathit{\ell m}}^{\mathit{X}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{.}$(45)Binning is then implemented by adding a sum over all ℓ inside a bin to the expression for the filtered map given in Eq. (45), thus obtaining the observed binned bispectrum of the map ${\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}\mathit{,}\hspace{0.17em}\mathrm{obs}}$, with bin indices i₁,i₂,i₃. The linear correction term is obtained in a similar way (and subtracted from the cubic term), but with two of the maps in Eq. (44) replaced by Gaussian simulations, and taking the average over a large number of those.

The theoretical templates are binned simply by summing the exact expression over the ℓs inside a bin⁸, and the same is true for the covariance matrix ${\mathit{h}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{2}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{1}}^{\mathrm{\prime }}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{2}}^{\mathrm{\prime }}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}{\mathit{X}}_{\mathrm{3}}^{\mathrm{\prime }}}$. The binning is optimized in such a way as to maximize the overlap, defined using the inner product of Eq. (39) between the binned and the exact template for all shapes under consideration. Finally the estimate of f_NL is computed using Eq. (38), where the inner product now contains a sum over bin indices i instead of multipoles ℓ, and the bispectra and covariance matrix are replaced by their binned versions.

3.4. Data set and simulations

In the following, we will apply our bispectrum estimation pipelines both to simulations and data, and consider a large number of shapes, either primordial or non-primordial in origin. Simulations will be used for a wide range of purposes, from comparisons of the outcomes of different estimators, to tests of instrumental systematics and foreground contamination, as well as Monte Carlo evaluation of error bars. For this reason, many different sets of simulated maps will be used, with features that will vary, depending on the specific application, and will be described case by case throughout the paper. Most of the time, however, we will use the FFP8 simulation data set described in Planck Collaboration XII (2016), or mock data sets obtained by processing FFP8 maps in various ways. These are the most realistic simulations available, modelling the CMB sky and the instrumental effects of Planck to the best of our current knowledge. They have passed through the same steps of the component separation pipelines as the real sky map and are the same maps as used for the final validation of the estimators in Sect. 5.3.

As far as actual data are concerned (Planck Collaboration I 2016; Planck Collaboration II 2016; Planck Collaboration III 2016; Planck Collaboration IV 2016; Planck Collaboration V 2016; Planck Collaboration VI 2016; Planck Collaboration VII 2016; Planck Collaboration VIII 2016) the maps analysed in this work are the Planck 2015 sky map, both in temperature and in E polarization, as cleaned with the four component separation methods SMICA, SEVEM, NILC, and Commander (Planck Collaboration IX 2016). As explained in Planck Collaboration VII (2016), the polarization map has had a high-pass filter applied to it, since the characterization of systematics and foregrounds in low-ℓ polarization is not yet satisfactory. This filter removes the scales below ℓ = 20 completely, and those between ℓ = 20 and ℓ = 40 partially. In all our analyses we use ℓ_min = 40 for polarization, in order to be independent of the details of this filter. For temperature, we use ℓ_min = 2. All the final cleaned maps are smoothed with a 5′ Gaussian beam in temperature, and a 10′ Gaussian beam in polarization.

The maps are masked to remove the brightest parts of the Galaxy as well as significant point sources. The masks used are the common masks of the Planck 2015 release in temperature and polarization, which are the union of the confidence masks for the four component separation methods⁹ (Planck Collaboration IX 2016). The sky coverages are respectively f_sky = 0.76 in temperature and f_sky = 0.74 in polarization. The stability of our results as a function of the mask is investigated in Sect. 7.2, where we show that our temperature and joint temperature plus polarization results do not change significantly when we consider a larger sky coverage.

In Sects. 7.1 and 7.3 we also compare the performance of different component separation methods, and conclude that, with respect to bispectrum estimation, the most accurate results are obtained using the SMICA data set. As already done for the 2013 release, we will thus consider SMICA as our main data set, using the other methods for important cross-checking purposes.

If we consider only temperature, current SMICA data become noise dominated at ℓ ≃ 2000, while previous nominal mission data were noise dominated at ℓ ≃ 1700. The mask used for the 2013 release was also slightly larger than the current one (f_sky = 0.73 in 2013 vs. f_sky = 0.76 in 2015). Since the f_NL signal-to-noise ratio, as quantified by the Fisher Information Matrix, scales as $\mathrm{(}\mathit{S}\mathit{/}\mathit{N}\mathrm{)}\mathrm{\propto }\mathit{\ell }\sqrt{{\mathit{f}}_{\mathrm{sky}}}$ in the signal dominated regime and saturates in the noise-dominated regime, we expect an improvement in our f_NL temperature constraints of about 20% when going from the 2013 nominal mission release to the current results. Adding polarization and accounting for all possible TTE, EET, and EEE bispectra produces further improvements. Since we are neglecting the first 40 polarization multipoles, such improvements are expected to be fairly small for shapes peaking in the squeezed limit, and more pronounced for equilateral type bispectra. A Fisher matrix approach shows that error bars are expected to improve by about 10% for the local shape and about 40% for the equilateral shape. This is in good agreement with our actual measurements, as can be seen from the results presented in Sect. 6 onwards.

3.5. Data analysis settings

Now we detail the general setup adopted for the analysis of Planck 2015 data by the four different optimal bispectrum estimation pipelines, described in previous sections.

As already explained, inpainting of the masked regions of the sky is a preliminary data filtering operation that all pipelines must perform, in order to retain optimality. We found that the inpainting method used in 2013 (Planck Collaboration XXIV 2014) for temperature maps still works well when polarization is included (note that it is the original T, Q, and U maps that should be inpainted, not the derived E map). We adopt a simple diffusive inpainting method. First the masked regions of the map are filled with the average value of the rest of the map. Then the value of each masked pixel is replaced by the average value of its (generally eight) direct neighbour pixels. The latter step is repeated a fixed number of times (2000)¹⁰. Relevant, final computations in map space (see e.g., Eq. (44)) are always done after remasking, so that the inpainted regions of the map are not used directly. The relevance of the inpainting procedure is that it reduces the effect of the sharp edges and the lack of large-scale power inside the mask leaking into the rest of the map during harmonic transforms.

For the linear correction term and to determine error bars, we use the FFP8 simulations (see Planck Collaboration XII 2016, and Sect. 3.4), filtered through the different component separation pipelines, using the same weights as used for the actual data when co-adding frequency channels. To compute all theoretical quantities (like the bispectrum templates and the ISW-lensing bias) we use the Planck 2015 best-fit cosmological parameters as our fiducial cosmology. However, results are quite insensitive to small changes in these parameters.

As pointed out in Sect. 3.4, low-ℓ multipoles are filtered out of the input polarization data set, so that all of our analyses will use ℓ_min = 40 in polarization, and ℓ_min = 2 in temperature. The choice of ℓ_max is dictated by the angular resolution of the cleaned maps, which is 5′ in temperature, and 10′ in polarization, and by the fact that the temperature data become noise-dominated at ℓ ≃ 2000, while the polarization information saturates around ℓ ≃ 1000. The KSW and binned estimators use ℓ_max = 2500 for temperature, while the modal estimators use ℓ_max = 2000. As shown explicitly in Sect. 7.4, results are completely stable between ℓ = 2000 and ℓ = 2500, so that this has no impact on f_NL. Similarly the binned estimator uses ℓ_max = 2000 for polarization, while the other estimators use ℓ_max = 1500, but again Sect. 7.4 shows that this difference is unimportant. The estimators also differ in the number of maps used to compute the linear correction term and the error bars, but generally it is of the order of 200. This difference is due to the different convergence properties of the estimators, some converging faster than others.

The binned bispectrum estimator uses 57 bins¹¹ for the analysis, which were determined by optimizing the correlation between the exact and the binned templates for the different shapes in temperature and polarization, as well as the full combined case. This is equivalent to minimizing the variance of the different f_NL parameters, where we focus in particular on the primordial shapes.

As previously explained, we use two different versions of the polarized modal pipelines, called “Modal 1” and “Modal 2” in the paper. Besides technical and conceptual implementation differences, the two modal estimators also use different sets of basis templates. The “Modal 1” pipeline uses 600 polynomial modes, plus nine “KSW radial modes”, computed at last scattering, while “Modal 2” has a basis formed by 2000 polynomial modes, augmented with a Sachs-Wolfe local bispectrum template. Due to the way polarization is implemented in the “Modal 2” pipeline, it cannot determine results for E-only. More details and explanations of the different choices are provided in Sect. 3.2.

As already stressed, the use of several independent bispectrum estimators, and several completely independent component separation methods allows a remarkable level of cross-validation of our results in order to establish their robustness. The fact that the bispectrum estimators are statistically equivalent and produce practically optimal results will be established in Sect. 5. The validation of the component separation methods is described in Planck Collaboration IX (2016) and Sect. 7.

4. Non-primordial contributions to the CMB bispectrum

Here we investigate several bispectra of non-primordial origin that are expected to be present in the data, and quantify their impact on our f_NL results. We devote particular attention to assessing potential biases that these NG signals might induce on the primordial bispectra. When forecasting such biases, we assume the data analysis settings discussed in Sects. 3.4 and 3.5.

4.1. Non-Gaussianity from the lensing-ISW bispectrum

The correlation between the gravitational lensing of the CMB anisotropies and the integrated Sachs-Wolfe (ISW) effect gives rise to a secondary CMB bispectrum – characterized by an oscillatory behaviour and peaked on squeezed configurations – that is a well-known contaminant to the primordial NG signal (Hanson & Lewis 2009; Mangilli & Verde 2009; Lewis et al. 2011; Mangilli et al. 2013). The temperature-only 2013 Planck results (Planck Collaboration XXIV 2014; Planck Collaboration XIX 2014; Planck Collaboration XVII 2014) showed evidence for the first time for the lensing-ISW CMB bispectrum and associated bias. Based on the same methodology used for the 2013 Planck data analysis (Planck Collaboration XXIV 2014), here we update the computation of the lensing-ISW bispectrum and its bias to include the full mission temperature and polarization data.

As shown by Cooray & Melchiorri (2006), the direct lensing-ISW correlation in E-polarization due to rescattering of the temperature quadrupole generated by the ISW effect is negligible. However, as explained in Lewis et al. (2011), there is an important correlation between the lensing potential and the large-scale E-polarization generated by scattering at reionization. Because the lensing potential is highly correlated with the ISW signal, this also leads to a non-zero lensing-ISW bispectrum in polarization.

To determine ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{LISW}}$, the amplitude parameter of the lensing-ISW bispectrum, one simply inserts the theoretical template for this shape into the general expression of Eq. (38). The template is given by (Hu 2000; Lewis et al. 2011) $\begin{array}{ccc}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}\mathit{,}\hspace{0.17em}\mathrm{LISW}}& \mathrm{=}& {\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{2}}\mathit{\phi }}\stackrel{}{{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{3}}}}_{\mathrm{˜}}}{\mathit{f}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{1}}}\mathrm{+}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}\mathit{\phi }}\stackrel{}{{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}}}_{\mathrm{˜}}}{\mathit{f}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{X}}_{\mathrm{1}}}\\ & & \mathrm{+}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}\mathit{\phi }}\stackrel{}{{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{2}}{\mathit{X}}_{\mathrm{3}}}}_{\mathrm{˜}}}{\mathit{f}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{2}}}\mathrm{+}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{X}}_{\mathrm{3}}\mathit{\phi }}\stackrel{}{{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{1}}{\mathit{X}}_{\mathrm{2}}}}_{\mathrm{˜}}}{\mathit{f}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{X}}_{\mathrm{2}}}\\ & & \end{array}$(46)Here ${\mathit{C}}_{\mathit{\ell }}^{\mathit{T\phi }}$ and ${\mathit{C}}_{\mathit{\ell }}^{\mathit{E\phi }}$ are the temperature/polarization-lensing potential cross power spectra, and the tilde on $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\end{array}$, $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{TE}}\end{array}$, and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{C}}_{\mathit{\ell }}^{\mathit{EE}}\end{array}$ indicates that it is the lensed TT, TE, or EE power spectrum. The functions ${\mathit{f}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathit{T,E}}$ are defined by $\begin{array}{ccc}{\mathit{f}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathit{T}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\left[{\mathit{\ell }}_{\mathrm{2}}\mathrm{(}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{+}{\mathit{\ell }}_{\mathrm{3}}\mathrm{(}{\mathit{\ell }}_{\mathrm{3}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{-}{\mathit{\ell }}_{\mathrm{1}}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\right]\mathit{,}\\ {\mathit{f}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathit{E}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{2}}\left[{\mathit{\ell }}_{\mathrm{2}}\mathrm{(}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{+}{\mathit{\ell }}_{\mathrm{3}}\mathrm{(}{\mathit{\ell }}_{\mathrm{3}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{-}{\mathit{\ell }}_{\mathrm{1}}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\right]\\ & & \mathrm{\times }\left(\begin{array}{c}\\ & & \\ & & \end{array}\right){\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)}^{-1}\mathit{,}\end{array}$(47)if ℓ₁ + ℓ₂ + ℓ₃ is even and ℓ₁,ℓ₂,ℓ₃ satisfy the triangle inequality, and zero otherwise.

In this paper our main concern with the lensing-ISW bispectrum is not so much to determine its amplitude (although that is also of great interest), but to compute its influence on the primordial shapes. The bias $\mathrm{\Delta }{\mathit{f}}_{\mathrm{NL}}^{\mathrm{P}}$ due to the lensing-ISW bispectrum on the estimation of a given primordial amplitude ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{P}}$ is given by $\mathrm{\Delta }{\mathit{f}}_{\mathrm{NL}}^{\mathrm{P}}\mathrm{=}\frac{\mathrm{⟨}{\mathit{b}}^{\mathrm{LISW}}\mathit{,}{\mathit{b}}^{\mathrm{P}}\mathrm{⟩}}{\mathrm{⟨}{\mathit{b}}^{\mathrm{P}}\mathit{,}{\mathit{b}}^{\mathrm{P}}\mathrm{⟩}}\mathit{,}$(48)where the inner product is defined in Eq. (39).

The values for the bias are given in Table 1. It should be noted that these are the results as computed exactly with Eq. (48). They can differ slightly from the ones used in e.g., Table 10, where each estimator adopts values computed using the approximations appropriate to the method. However, these differences are completely insignificant. As seen already in Planck Collaboration XXIV (2014), for T-only the bias is very significant for local and to a lesser extent for orthogonal NG. For local NG the bias is larger than the error bars on f_NL. We see that for E-only the effect is non-zero but not significant. For the full T+E case, the bias is smaller than for T-only, but large enough that it is important to take into account.

The results for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{LISW}}$ can be found in Table 2. The polarized version of the template has only been implemented in the binned bispectrum estimator. Error bars have been determined based on FFP8 simulations as usual¹². The KSW estimator implements the lensing-ISW template exactly, while the binned and modal estimators use approximations, as explained in Sect. 3. In particular for the binned estimator the correlation between the binned and exact lensing-ISW template is relatively low, since it is a difficult template to bin (unlike all the other templates considered in this paper), which is reflected in the larger error bars. Tests performed on FFP8, as well as other tests, demonstrate that the lower correlation does not lead to a bias compared to the other estimators. We will use the KSW results to draw our conclusions.

We see that temperature results from the full mission are consistent with the 2013 nominal mission (Planck Collaboration XXIV 2014). Including polarization yields results that also appear consistent and decrease the error bars. However, for now the T+E conclusions should be considered preliminary, for the reasons related to polarization data discussed in detail in Sects. 6 and 7. The error bars will also improve when measured with the other bispectrum estimators. As already seen in 2013, the values for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{LISW}}$ are slightly low compared to the expected value of 1, but not significantly so. On the other hand, the detection of the lensing-ISW bispectrum is significant, even with our conservative rescaling of the error bars. The hypothesis of having no lensing-ISW bispectrum is excluded at 2.8σ using temperature alone, and improves to 3.0σ with the current preliminary result when including polarization. As mentioned above, the latter result is likely to improve with further analysis of the Planck data. In Fig. 1 we present the results of the skew-C_ℓ analysis for lensing-ISW NG for the T map, which illustrates that the instrument and data processing are not removing this expected NG signal from the data.

4.2. Non-Gaussianity from extragalactic point sources

The auto-bispectra of extragalactic point sources are a potential contaminant to primordial NG estimates at Planck frequencies. The basic modelling and methodology of this section follows the corresponding section in Planck Collaboration XXIV (2014).

Extragalactic point sources are divided into populations of unclustered and clustered sources. The former are radio and late-type infrared galaxies (see e.g., Toffolatti et al. 1998; González-Nuevo et al. 2005), while the latter are dusty star-forming galaxies constituting the cosmic infrared background (CIB; Lagache et al. 2005). The contamination due to both types of sources in NG estimators is handled via dedicated bispectrum templates which are fitted jointly with the primordial NG templates.

The unclustered sources have a white noise distribution, and hence constant polyspectra. Their reduced angular bispectrum template is thus ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{unclust}}\mathrm{=}\mathrm{const}\mathit{.}$(49)This constant is usually noted b_PS or b_src in the literature (e.g., Komatsu & Spergel 2001). This constant template is valid in polarization as well as temperature, since the polarization angles of point sources are less clustered than the source density. However, since not all these point sources are polarized, we do not measure the same sources in temperature and in polarization. In fact, there is no detection of the bispectrum of unclustered point sources in the cleaned Planck polarization map, unlike in the temperature map, where Table 3 (binned bispectrum estimator) and Fig. 2 (skew-C_ℓs) show a clear detection.

The clustered sources (CIB) have a more complex bispectrum in temperature, reflecting the distribution of the large-scale structure and the clustering of galaxies in dark matter halos (Argüeso et al. 2003; Lacasa et al. 2012; Crawford et al. 2014). The Planck results have allowed the measurement of the CIB bispectrum at frequencies 217, 353, 545 GHz in the range ℓ ≃ 200 − 700 (Planck Collaboration XXX 2014). In this multipole range, a power law was found to fit the measurement, with an exponent consistent between frequencies. However, at lower multipoles theoretical models for the CIB power spectrum (e.g., Planck Collaboration XXX 2014) and bispectrum (Lacasa et al. 2014; Pénin et al. 2014) predict a flattening of the CIB power. We thus take the TTT CIB bispectrum template to be a broken power law, ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{CIB}}\mathrm{\propto }{\left[\frac{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\ell }}_{\mathrm{1}}\mathit{/}{\mathit{\ell }}_{\mathrm{break}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\ell }}_{\mathrm{2}}\mathit{/}{\mathit{\ell }}_{\mathrm{break}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\ell }}_{\mathrm{3}}\mathit{/}{\mathit{\ell }}_{\mathrm{break}}\mathrm{)}}{\mathrm{(}\mathrm{1}\mathrm{+}{\mathit{\ell }}_{\mathrm{0}}\mathit{/}{\mathit{\ell }}_{\mathrm{break}}{\mathrm{)}}^{\mathrm{3}}}\right]}^{\mathit{q}}\mathit{,}$(50)where the index is q = 0.85, the break is located at ℓ_break = 70, and ℓ₀ = 320 is the pivot scale for normalization. Dusty star-forming galaxies emit with a low polarization fraction, and polarization correlates only over the smallest scales, so that the CIB is negligibly polarized. We thus take vanishing templates for its polarized bispectra ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{CIB}\mathit{,TTE}}\mathrm{=}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{CIB}\mathit{,TEE}}\mathrm{=}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{CIB}\mathit{,EEE}}\mathrm{=}\mathrm{0.}$(51)Both point source templates, Eqs. (49) and (50), have been implemented in the binned bispectrum estimator described in Sect. 3. The results for these two templates applied to the Planck temperature map, cleaned with the four component separation methods, can be found in Table 3. Since the two templates are highly correlated, the results have been determined in a combined analysis. The results have also been determined jointly with the primordial local, equilateral, and orthogonal templates, and the lensing-ISW bias has been subtracted, but all of this makes a negligible difference. Contamination from unclustered sources is detected in all component-separated maps. However, A_CIB is not detected.

The order of magnitude of the bispectrum amplitudes found in Table 3 is consistent with expectations. Indeed, for radio sources at 217 GHz and with a flux cut based on the Planck ERCSC (Planck Collaboration VII 2011), Lacasa & Aghanim (2014), the forecasted b_PS is approximately 2 × 10^-28. For the CIB, the Planck 2013 measurement (Planck Collaboration XXX 2014) at 217 GHz gives A_CIB ≃ 6 × 10^-27 when translated into dimensionless units. The results reported in Table 3 are consistent at the order-of-magnitude level with these estimates, although they are lower, because we are analysing cleaned maps.

The unclustered point source and CIB templates are highly correlated, at 93%. For this reason it was not deemed a priority for the other bispectrum estimators to implement the CIB template as well. Moreover, both point source templates are negligibly correlated with the primordial NG templates and the lensing-ISW template (the maximum being the correlation between equilateral and CIB templates at 2.7%, while correlations with the unclustered point source template are well below 1%). For this reason, and despite the detection of point sources in the cleaned maps, it makes no difference for the primordial results if point sources are included in a joint analysis or completely neglected.

An additional contaminant to the cosmological CMB bispectrum arises from the correlation between the gravitational lensing of the CMB anisotropies and the CIB anisotropies. This correlation was detected in Planck Collaboration XVIII (2014) using an optimal cross-spectrum estimator. The CIB-lensing bispectrum might couple with any of the primordial shapes. However, the amplitude of the CIB bispectrum is predicted to be small in the cleaned Planck maps and it has actually not been detected (see Table 3). The CIB-lensing bispectrum signal is frequency dependent, and it is mostly dominant in the very high Planck frequencies (see e.g., Curto et al. 2015).

4.3. Non-Gaussianity from residuals of the deglitching processing

Cosmic rays interacting with the cryogenic detectors induce spikes in timelines. These high-amplitude, fast-rising signals are followed by a decay tail. We observe three families of glitches, characterized by their temporal shape. The amplitude and time constants of the decays depend on which part of the satellite is hit (Catalano et al. 2014; Planck Collaboration X 2014). These random events are Poisson-distributed in time and produce highly non-Gaussian systematics.

A method has been developed to remove them directly at the time-ordered information (TOI) level. This process is performed iteratively, and is described in detail in Planck Collaboration X (2014). The short glitches are just flagged in the data, whereas for the long ones only the fast part is flagged, and the long tail is substracted from the timeline. This procedure is not perfect, and there are residuals from the potentially biased errors in the fit, and the undetected glitches under the threshold of 3.2σ of the TOI noise rms. They could in principle produce a non-Gaussian signal in the final map. In addition, these residuals could interact with the mapmaking procedure at the destriping level, since the error on the offset determination could be non-Gaussian due to undetected glitches or a possible bias in the errors of the removal of tails. This is important, because in more than 95% of the TOI data, tails have been subtracted.

To estimate the effect of these residuals on the determination of NG, we created two sets of simulations (one including glitches and the other not) for every bolometer of the 143 GHz channel. We generated Gaussian CMB maps, and applied the full TOI processing with a realistic instrumental noise (Planck Collaboration VII 2016). In the simulations with glitches, we added glitches at the TOI level, following the properties measured in the data, and cleaned them with the procedure applied to the data. For the simulations without glitches, we have the same CMB and noise realization, but no glitches added at the TOI level.

We estimated the bias caused by glitches on the measurement of f_NL using the binned bispectrum estimator. The bias on f_NL induced by the glitch residuals g on a map T, including noise and CMB is given by , where the noise in the weighting of the estimator is determined from the simulations with glitches (as it would be for the data). Results are shown in Table 4. For most shapes, we detect no significant bias. The higher signal and high dispersion for the local shape might be due to a mis-calibration of the linear correction. In any case, for all shapes the bias due to glitches is a negligible correction to the value of f_NL, given its error bars, and we will not take it into account in the remainder of the paper.

5. Validation tests

During the work for the 2013 release, culminating in the NG results of Planck Collaboration XXIV (2014), the advantage of having multiple independent bispectrum estimator implementations was amply demonstrated. This allows for very useful cross-checking of results, both during development and for the final analysis, thus greatly improving the robustness of and confidence in the final results. For this new release we followed the same procedure, with the same three principal bispectrum estimators: KSW; binned; and modal, all of which had their pipelines updated to handle polarization data in addition to temperature.

Beyond the usefulness of cross-checking, the three estimators complement each other and have different strengths. The KSW estimator can treat separable bispectrum templates without approximation, but it is more work to add new templates and non-separable templates cannot be handled at all. The binned and modal estimators can reconstruct the full bispectrum (smoothed in different domains), while the skew-C_ℓ extension of the KSW estimator can be used to investigate the bispectrum beyond f_NL. The binned bispectrum estimator is the fastest on a single map or a set of unrelated maps, but becomes slower than the other two on a large set of realizations based on the same settings, because the linear correction term cannot be precomputed. The modal estimator can investigate a wide selection of oscillating or otherwise rapidly changing bispectrum templates that would be difficult to bin, while the binned bispectrum estimator can quickly implement and determine the f_NL of an additional template or the effect of a different cosmology if the binned bispectrum of the maps has already been computed. The binned estimator gets the dependence of f_NL on ℓ for free with its results, while the modal estimator allows for a statistical investigation of the mode coefficients.

In this section we show some of the validation tests, in particular for polarization. In Sect. 5.1 we investigate the agreement between estimators, map-by-map, on sets of successively more realistic maps. In Sect. 5.2 we show that the estimators are unbiased in the presence of a non-zero f_NL. Finally, in Sect. 5.3 we show that the estimators are essentially optimal on a set of the most realistic Planck simulations available, which are those used to compute the error bars on our final results.

5.1. Agreement between estimators on a map-by-map basis

The maps used in this subsection are realistic simulations of the CMB (at resolution N_side = 2048) but without any foregrounds. They do not contain any primordial NG, but do include ISW-lensing. Since the final FFP8 simulations were not yet available, the main goal was to make sure that the estimators agreed with each other, not only on average, but also on a map-by-map basis. For this purpose it was enough to look at only 49 maps. Establishing optimality of the estimators requires a larger number of maps, and is shown on the FFP8 simulations in Sect. 5.3.

In our first test we include the effect of the 143 GHz beam, but in other respects the simulations are ideal (no noise, and no mask). The analysis used ℓ_max = 2000 for both T and E. The results for the average over the maps for the KSW, binned, and both modal estimators, as well as for the difference between each estimator and KSW, are shown in Table 5. The shapes are assumed to be independent in this analysis, which means that the bias on the local shape due to the ISW-lensing effect is clearly visible. Results are shown for T-only, E-only, and the full combined T+E analysis. Note that the second modal implementation cannot compute results for E alone. One clearly sees that the results agree very well. It is also interesting to note that in this ideal noiseless case, one can actually determine f_NL more accurately from polarization alone than from temperature alone. This is due to the narrower transfer function in polarization, so that the primordial bispectrum is less smoothed in its projection to two-dimensional harmonic space.

The second test is identical to the first, except that we add realistic anisotropic noise realizations to the full-sky maps, based on the 143 GHz channel. The estimators now require the use of the linear correction term, and results are shown in Table 6. The agreement is still very good, although slightly worse than in the ideal case, as expected. The fact that the error bars for the T-only local case here are actually a bit smaller than in the ideal case is an artefact of the small number of maps; i.e., the error bars have not completely converged yet. On the other hand, the fact that the error bars for E-only are much larger than in the ideal case is a real effect; the Planck single-frequency polarization maps are noise-dominated.

Finally, the third test is identical to the second, except that we now also add a mask. The mask chosen is realistic, based on the union of the confidence masks provided by the SMICA, NILC, SEVEM, and Commander methods for this particular set of simulations. It contains both a Galactic and a point source part. The temperature mask leaves 79% of the sky unmasked, while the polarization mask leaves 76%. The results are shown in Table 7, while the map-by-map comparison is given in Fig. 3. From the table we see that the agreement between the different bispectrum estimators is still very good and only slightly degraded when compared to the previous case. The typical discrepancy between the bispectrum estimators, even in this most realistic case, is less than about a third of the uncertainty on f_NL. This is apparent in the map-by-map comparison of Fig. 3.

5.2. Validation of estimators in the presence of primordial non-Gaussianity

After the map-by-map comparison of the previous section, we next want to make sure that the estimators are unbiased. For this purpose we prepared a different set of 100 T and E CMB simulations, still with cosmological parameters as determined by Planck. This time ISW-lensing is not present, but there is a nonzero local f_NL = 12. To these maps we add the same beam, anisotropic noise, and mask as before. We again take ℓ_max = 2000, and the results are given in Table 8.

We see that all the estimators correctly recover the input value, both in temperature and in polarization. The results for the equilateral and orthogonal shapes are consistent with the fact that those templates have a non-zero correlation with the local shape (the table gives the results for an analysis where all shapes are assumed independent). For example, a joint analysis of the T+E binned estimator gives ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{11.5}\mathrm{\pm }\mathrm{6.4}$, ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{-}\mathrm{7.5}\mathrm{\pm }\mathrm{51}$, and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}\mathrm{=}\mathrm{-}\mathrm{0.4}\mathrm{\pm }\mathrm{29}$. Except for the first modal estimator in E-only (due to an insufficient number of maps in the linear correction term), we also find that the error bars for the bispectrum-based estimators are very close to the Fisher errors. Note that a slight increase in the error bars compared to Fisher estimates is expected for the local shape in T-only and in T+E, due to the signal being significantly different from zero there (the Fisher error bars for the local case are 5.8 for T-only, 26 for E-only, and 5.0 for T+E). Hence the estimators are effectively optimal, as will be illustrated in more detail in the next section.

5.3. Validation of estimators on realistic Planck simulations

As a final validation test, we ran our estimators on a large set of the most realistic simulations available. These are the FFP8 simulations (Planck Collaboration XII 2016) using SMICAfor foreground separation. They are the same simulations we use to determine the error bars on our final SMICA results in Sect. 6. They contain the Collaboration’s best estimates of the CMB sky and of Planck’s noise and beam effects, and have been cleaned by SMICA in the same way as the real sky map. The mask used is the same common mask defined for the real data analysis. For this test the estimators were all processed with the same settings used for the final data analysis.

Here we take 159 of these maps, and process these using all the estimators. By contrast, for the final results in Sect. 6, the convergence of the error bars of each estimator was carefully checked, using more maps if required. This explains why there are some small differences between the error bars in Sect. 6 and the ones presented here.

The results are shown in Table 9. Note that these are the results from an independent analysis, without subtracting the ISW-lensing bias. We also show the results from Minkowski functionals (for the local case only)¹³. We see that there is very good agreement between the bispectrum estimators even on this most complex and realistic set of simulations. The standard deviation of the difference between bispectrum estimators generally stays below one third of the error bar on f_NL, the only exception being the T-only equilateral result for the Modal 1 pipeline, which is still smaller than one half of the error bar. We see that the results from Minkowski functionals are consistent, but clearly suboptimal for f_NL. They are however a valuable, independent check.

The exact Fisher error bars for the nine shapes considered in the table are, in the same order as the table: 5.4, 69, 35; 31, 131, 74; 4.7, 43, 21. Taking into account the relative error in the standard deviation of $\mathrm{1}\mathit{/}\sqrt{\mathrm{2}\mathrm{(}\mathit{n}\mathrm{-}\mathrm{1}\mathrm{)}}$, which is 5.6% for 159 maps, we see that all bispectrum estimators are effectively optimal on all shapes, except for the E-only equilateral case where they appear slightly suboptimal. The small suboptimality of the Modal 2 pipeline for the local shape seen here disappears once more maps are used (see the results in Sect. 6).

In conclusion, all these validation tests show that we have very good agreement between the results from the different bispectrum estimators, not just on average, but also on a map-by-map basis. In addition we see that, despite the approximations made in the pipelines, and the simple treatment of the masked part of the maps (diffusive inpainting method and f_sky factor), the bispectrum estimators we use are all essentially optimal.

6. Results

6.1. Constraints on local, equilateral, and orthogonal f_NL

In this section we investigate the local, equilateral, and orthogonal primordial templates. These are now established as the standard shapes to study first when investigating the bispectrum (see Sect. 2 for a theoretical motivation and description of these shapes). However, they represent only the tip of the bispectral iceberg, and many more shapes are investigated in Sect. 8, while full model-independent reconstructions of the bispectrum are presented in Sect. 6.2.

For a complete description of the Planck data set and the bispectrum estimator configurations we have used, we refer the reader in particular to Sects. 3.4 and 3.5. To summarize the overall analysis methodology, we have employed four independent bispectrum estimators on the full mission Planck temperature and polarization maps obtained from the four different component separation pipelines, SMICA, SEVEM, NILC  and Commander. The bispectrum estimators are the KSW estimator with its skew-C_ℓ extension using exact separable templates (Sect. 3.1), the Binned estimator using fixed multipole bins (Sect. 3.3), and the Modal 1 and Modal 2 estimators, which both use separable eigenmode expansions (Sect. 3.2). Temperature is analysed over the multipole range ℓ_min = 2 to ℓ_max = 2000 or above and polarization is analysed from ℓ_min = 40 to ℓ_max = 1500 or above (Sect. 3.5). By employing inpainting and a linear term, all these estimators essentially achieve optimality (as shown by comparison with Fisher matrix forecasts). The linear term in Eq. (36) and the uncertainties are determined using the FFP8 simulations (Sect. 3.5), also processed through the different foreground-separation pipelines. Our thorough validation campaign for these estimators is presented in Sect. 5.

The results of the analysis of the four cleaned maps with the four estimators, for T-only, E-only, and full T+E, are shown in Table 10, which is one of the main results of this paper. Results are determined while assuming all shapes to be independent, and are shown both with and without subtraction of the ISW-lensing bias (see Sect. 4.1 for more details about ISW-lensing). This bias is most important (relative to the size of the error bars) for the local shape, but also non-negligible for the orthogonal shape. Results here have not been marginalized over the point source contributions. While Sect. 4.2 shows that there is still a significant contamination by unclustered point sources in the cleaned maps, the correlation with the primordial templates is so small that this has no impact on the results reported here (as checked explicitly).

While Table 10 is the main result of this section, in order to simplify the use of the Planck results by the wider scientific community, we also present in Table 11 the results that can be considered the final Planck 2015 results for the local, equilateral, and orthogonal shapes. As in 2013, we select the combination of the KSW estimator and the SMICA map for this. The SMICA map consistently performs well in all data validation tests, which are discussed in detail in Sect. 7. The KSW estimator, while unable to deal with non-separable templates, treats separable templates exactly, and the local, equilateral, and orthogonal template are all separable. On the other hand, the binned and modal estimators can deal with non-separable shapes and have other advantages as well (like full bispectrum reconstruction), but at the price of using approximations for the templates. However, they have all been optimized in such a way that the correlation with the exact templates for the three primordial shapes is close to perfect, so that in the end the results by the different estimators are statistically equivalent. Compared to the corresponding values in Table 10, the difference in the numbers in the last column of Table 11 is due to the fact that in the latter the equilateral and orthogonal f_NL have been determined jointly.

Focusing on the results for temperature-only and the full temperature plus polarization (T+E) results, we see that there is no evidence for a non-zero bispectrum with any of these three primordial shapes (local, equilateral, and orthogonal). After ISW-lensing subtraction, all f_NL values are consistent with 0 at 68% CL. The temperature results are all very similar to the ones from the nominal mission data published in 2013 (Planck Collaboration XXIV 2014), with very minor improvements in the error bars due to the additional temperature data. We also see that results are quite consistent when including polarization, with error bars shrinking by about 15% for local, 35% for equilateral, and 40% for orthogonal.

Table 10 displays very good agreement between the results from the different estimators, at the level expected from the validation tests in Sect. 5. We also note how the error bars, which are determined using the FFP8 simulations, are statistically indistinguishable from the optimal Fisher expectation.

Different component separation methods also show a good level of agreement when looking at temperature-only and combined temperature plus polarization results. The accuracy of this statement will be shown and quantified in detail in Sect. 7. However, in the same section, we will also show how the agreement between f_NL extracted from different cleaned maps becomes significantly degraded when considering polarization-only results¹⁴. The reasons behind this loss of internal consistency are not fully understood at present. Polarization data are, however, much noisier than temperature data, implying that the EEE bispectra have a close to negligible weight in the final combined measurement, which is dominated by the TTT and TTE configurations. In fact, as just mentioned above, the combined measurement looks perfectly self-consistent: local, equilateral and orthogonal f_NL measurements in the T+E column of Table 10 pass all our tests of robustness.

We can thus conclude that, while highly challenging from a technical point of view, the inclusion of polarization in our estimator pipelines has been a success, allowing for a significant tightening of the constraints on the three standard primordial bispectrum shapes. On the other hand, in light of the outstanding issues in E-only results, we present our results conservatively, and recommend the reader to consider all f_NL constraints that make use of polarization data throughout this paper as preliminary at the current stage. We stress again that this is a conservative choice, which is made despite the fact that no test to date shows any evidence of leakage of the issues in EEE bispectra into the T+E measurements. A detailed description of all the data validation tests, which lead to the robustness-related assessments summarized here, can be found in Sect. 7 (for readers less interested in the technical details, the main results and conclusions of all these tests are summarized in Sect. 7.6).

6.2. Bispectrum reconstruction

6.2.1. Modal bispectrum reconstruction

The starting point for modal bispectrum estimation is the robust extraction of the modal coefficients β_n from each of the full mission foreground-separated maps, that is, SMICA, SEVEM, NILC, and Commander. The β_n-coefficients are obtained for each of the temperature, polarization, and mixed bispectrum components, TTT, TTE, TEE, and EEE. Their cross-correlation between cleaning methods is an important validation of their accuracy, as we shall discuss in the next section, with excellent correspondence for temperature and some differences remaining in polarization. The modal basis number n_max = 2001 for the full mission analysis has been substantially increased offering a higher effective resolution when compared to the 2013 Planck Data Release where n_max = 601 modes were used. Several different basis functions have been used, including trigonometric functions, sinlog basis functions, and polynomials (closely related to Legendre functions), with the latter chosen because of excellent convergence in the squeezed and flattened limits.

We can reconstruct the full three-dimensional Planck bispectrum, obtained using these basis functions, to visualize its main properties and to determine robustness. A comparison between the temperature-only bispectra from the Nominal Mission and full mission at the same n_max = 601 modal resolution is shown in Fig. 4. Note the excellent agreement with all the main features replicated in the new data. In Fig. 4 in the third bispectrum, we also demonstrate the much higher bispectrum resolution achieved with the full n_max = 2001 modes. The tetrapyd shape reflects the constraints on the wavenumbers ℓ₁,ℓ₂, and ℓ₃, with the squeezed configuration appearing on the axes that lie along one ℓ_i = 0. The expected ISW-lensing bispectrum is an oscillating signal in the squeezed limit along the tetrapyd edges; it is now measured with a significance of 3.0σ (see Sect. 4.1). This ISW-lensing signal sets an interesting benchmark or threshold against which to compare the other strong features observed in the bispectrum and now defined with greater precision. The original “plus-minus” feature, with a large positive red peak around ℓ ≈ 150 followed by a larger negative peak near ℓ ≈ 250, remains though with more substructure, together with a broad negative peak in the equilateral limit around ℓ ≈ 900, which can be associated with the third acoustic peak from the transfer functions. Oscillatory models, which can connect these three peaks, achieve higher significance. The apparent signal observed in the flattened limit remains, with a distinct pattern of blue and red features on the surface of the tetrapyd.

We also include a comparison with WMAP-9 in Fig. 5, where we have restricted the reconstructions to ℓ_max = 600 for comparison with n_max = 601 modes. These plots, using identical isosurfaces, show the same bispectrum structure including the “plus-minus” feature clearly bisecting the main ℓ = 200 peak and the first oscillation of the ISW-lensing bispectrum visible along the lower tetrapyd edges. The WMAP-9 reconstruction only shows significant differences from Planck in the top right region, where the higher noise levels in WMAP-9 make its reconstruction less reliable.

All four components of the temperature and polarization bispectrum reconstruction obtained from SMICA are shown in Fig. 6. A direct comparison of the EEE polarization bispectrum for SEVEM, NILC and Commander, is shown in Fig. 7, where we note that these are orthogonalized E-mode contributions (see the Modal 2 discussion in Sect. 3). It is interesting to observe patterns of features evident in the polarization bispectra from the different foreground-cleaned maps, which, although inherently noisier, have qualitative similarities. At a quantitative level, however, the polarization bispectra modes from different methods are less correlated in polarization than in temperature, as we discuss in Sect. 7.

6.2.2. Binned bispectrum reconstruction

The (reconstructed) binned bispectrum of a given map is a natural product of the binned bispectrum estimator code (see Sect. 3.3). To test if any bin has a significant NG signal, we study the binned bispectrum divided by its expected standard deviation, a quantity for which we will use the symbol ℬ_i1i2i3. With the binning used in the estimator, the pixels are dominated by noise. We thus smooth in three dimensions with a Gaussian kernel of a certain width σ_bin. To avoid edge effects, due to the sharp boundaries of the domain of definition of the bispectrum, we renormalize the smoothed bispectrum, so that the pixel values would be normal-distributed for a Gaussian map.

In Figs. 8 and 9, we show slices of this smoothed binned signal-to-noise bispectrum ℬ_i1i2i3 with a Gaussian smoothing of σ_bin = 2, as a function of ℓ₁ and ℓ₂. Very red or very blue regions correspond to a significant NG of any type. The two figures only differ in the value chosen for the ℓ₃-bin, which is [ 518,548 ] for the first figure, and [ 1291,1345 ] for the second. We have defined two cross-bispectra here: ${\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{T}\mathrm{2}\mathit{E}}\mathrm{\equiv }{\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{TTE}}\mathrm{+}{\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{TET}}\mathrm{+}{\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{ETT}}$; and ${\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{TE}\mathrm{2}}\mathrm{\equiv }{\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{TEE}}\mathrm{+}{\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{ETE}}\mathrm{+}{\mathit{B}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{EET}}$. These two cross-bispectra are then divided by their respective standard deviations (taking into account the covariance terms) to produce the corresponding ${\mathrm{ℬ}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{T}\mathrm{2}\mathit{E}}$ and ${\mathrm{ℬ}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{TE}\mathrm{2}}$. Those three different permutations are not equal a priori due to the condition i₁ ≤ i₂ ≤ i₃, which is implemented in the code to reduce computations by a factor of six. However, part of the smoothing procedure involves adding the other five identical copies, so that in the end the plots are symmetric under interchange of ℓ₁ and ℓ₂ (and ℬ_i1i2i3 is symmetric under interchange of all its indices). The grey areas in the plots are regions where the bispectrum is not defined, either because it is outside of the triangle inequality, or because of the limitation ${\mathit{\ell }}_{\max }^{\mathit{E}}\mathrm{=}\mathrm{2000}$. Given that in both plots ℓ₃ is fixed at less than 2000, ${\mathrm{ℬ}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{TE}\mathrm{2}}$ is not defined if both ℓ₁ and ℓ₂ are larger than 2000, while ${\mathrm{ℬ}}_{{\mathit{i}}_{\mathrm{1}}{\mathit{i}}_{\mathrm{2}}{\mathit{i}}_{\mathrm{3}}}^{\mathit{EEE}}$ is undefined if either ℓ₁ or ℓ₂ (or both) are larger than 2000.

Results are shown for the four component separation methods SMICA, SEVEM, NILC, and Commander, and for TTT, T2E, TE2, and EEE. In addition we show on the second line of each figure the result for TTT with the radio (unclustered) and CIB (clustered) point source bispectra subtracted according to their jointly measured amplitudes. It is clear, in particular in the second figure, that at higher ℓ there is a very significant point source contamination in the cleaned TTT bispectra, in agreement with the results of Table 3. However, after removing it we do not see a clear signal of any other residual NG. Of course this is for the moment only a qualitative statement; more quantitative tools for studying the amount of NG in these smoothed bispectra are in development.

Looking at the polarized bispectra in the high-ℓ₃ slices, in particular TE2 and EEE, we do see some bluer and redder regions that might indicate residual NG. This agrees with statements made earlier, and discussed in greater detail in the next section, that the Planck polarized bispectrum is for the moment not as clean and well-understood as the temperature one. We also see a very good qualitative agreement between the four component separation methods in temperature, which worsens somewhat when mixing in more and more polarization; in particular SMICA and NILC give very similar results.

6.3. Primordial curvature reconstruction

In this section, we compress the information in T and E maps into maps of projected primordial curvature fluctuations, ζ. To the extent that the primary CMB temperature and polarization are Gaussian fields and the early universe ζ fluctuations are presumed Gaussian, all of the information in T,E is encoded in the mean field (Wiener-filtered map) plus fluctuations, characterized by a covariance matrix; it is just a different compression/re-expression of the original statistics. Thus, one could look for non-Gaussian deviations by directly evaluating the 3-point and higher statistics of this re-expression. Since fluctuations as well as the mean are given, a full description of the errors follows. That a Gaussian assumption is made does not mean that the Wiener-filtered map is in fact Gaussian, since the constraining data may drive it from a map consistent with Gaussianity. Hence one can search in the statistics of the Wiener-filtered map for evidence for non-Gaussianity.

The weighting of the temperature and E-mode polarization associated with the Wiener filter is an optimal one (inverse-total-covariance weighted), and estimators for non-Gaussianity constructed using this expression are nearly optimal. However, since a (fiducial) primordial power spectrum is assumed for the Wiener filter (usually a uniform n_sζ-spectrum, using the best fit Planck 2015 scalar spectral index and power amplitude A_s), this approach is best suited for perturbative non-Gaussianity of the sort we treat in this paper. In practice, the estimators used here for non-Gaussianity act directly on temperature and polarization data, rather than through the intermediary of optimal ζ-maps and their fluctuations.

The scalar fluctuations can be expressed in terms of the curvature variable ζ(x) = lna(x) on uniform total density hypersurfaces, where a is the inhomogeneous expansion factor, or, equivalently, by its wavenumber transform ζ(k). In turn, ζ(k) can be expanded in multipoles, ζ_LMκ(k). Instead of the magnitude of the wavenumber k, a mixed representation gives a multipole expansion at each comoving distance from our location, χ, ζ_LMκ(χ), M = 0,...,L. Here κ = c or s, with κ = c referring to the real part of ζ_LM and κ = s to the imaginary part. For M = 0 there is no s component, only c. The mean of the ζ field given the temperature field T and its covariance describing allowed fluctuations about that mean are $\begin{array}{ccc}& & \\ & & \mathrm{\times }\mathrm{⟨}\mathit{\delta }{\mathit{\zeta }}_{\mathit{LM\kappa }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\mathit{\delta }{\mathit{\zeta }}_{{\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}{\mathit{\kappa }}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\mathit{T}\mathrm{⟩}\mathrm{=}{\mathit{\delta }}_{\mathit{L}{\mathit{L}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{M}{\mathit{M}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{\kappa }{\mathit{\kappa }}^{\mathrm{\prime }}}\\ & & \mathrm{\times }\{{\mathit{C}}_{\mathit{L}}^{\mathit{\zeta \zeta }}\mathrm{\left(}\mathit{\chi ,}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{-}{\mathit{C}}_{\mathit{L}}^{\mathit{\zeta T}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathrm{[}{{\mathit{C}}_{\mathit{L}}^{\mathit{TT}}}^{\mathrm{]}}}^{-1}{\mathit{C}}_{\mathit{L}}^{\mathit{T\zeta }}\mathrm{\left(}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{\right)}\}\mathit{.}\end{array}$In these expressions, C^TT is the total covariance for the temperature, including both signal and noise variances. The specific forms assume the C^TT matrix is diagonal in multipole space, depending only on ${\mathit{C}}_{\mathit{L}}^{\mathit{TT}}$, but if the noise is inhomogeneous it will have off-diagonal components and the equations become matrix equations. In this section, we assume the noise is diagonal. Equation (52) shows the unsurprising result that for each LMκ only one mode is determined by T. Replacing T by E gives $\mathrm{⟨}{\mathit{\zeta }}_{\mathit{LM\kappa }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\hspace{0.17em}\mathrm{|}\hspace{0.17em}{\mathit{a}}_{\mathit{LM\kappa }}^{\mathit{E}}\mathrm{⟩}$ and ⟨ δζδζ^† | E ⟩. Although T and E are correlated, the uncorrelated part of E delivers a different mode for ζ from the one given by T. Thus when ζ is constrained by both T and E, the two modes deliver substantially more information than for T and E alone, and the fluctuations about the mean are thereby diminished. This is further helped by the acoustic oscillations of polarization being out of phase with those for T, so when T is down, E is not, and vice versa. The interplay of the two modes is quantified by signal-to-noise in Fig. 10.

When both T and E are included as constraints, a two by two matrix appears which includes the C^TE correlation in the off-diagonal as well as the C^TT and C^EE along the diagonal: $\begin{array}{ccc}& & \mathrm{⟨}{\mathit{\zeta }}_{\mathit{LM\kappa }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\hspace{0.17em}\mathrm{|}\hspace{0.17em}{\mathit{a}}_{\mathit{LM\kappa }}^{\mathit{T}}\mathit{,}{\mathit{a}}_{\mathit{LM\kappa }}^{\mathit{E}}\mathrm{⟩}\mathrm{=}\\ & & {\mathrm{[}\begin{array}{c}\\ & \end{array}\mathrm{]}}^{†}{\left(\begin{array}{c}\\ & \\ & \end{array}\right)}^{-1}\left[\begin{array}{c}\\ \\ \end{array}\right]\mathit{.}\\ & & \mathrm{⟨}\mathit{\delta }{\mathit{\zeta }}_{\mathit{LM\kappa }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\mathit{\delta }{\mathit{\zeta }}_{{\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}{\mathit{\kappa }}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{\right)}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\mathit{T,E}\mathrm{⟩}\mathrm{=}{\mathit{\delta }}_{\mathit{L}{\mathit{L}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{M}{\mathit{M}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{\kappa }{\mathit{\kappa }}^{\mathrm{\prime }}}\{{\mathit{C}}_{\mathit{L}}^{\mathit{\zeta \zeta }}\mathrm{\left(}\mathit{\chi ,}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{\right)}\\ & & \mathrm{-}{\mathrm{[}\begin{array}{c}\\ & \end{array}\mathrm{]}}^{†}{\left(\begin{array}{c}\\ & \\ & \end{array}\right)}^{-1}\left[\begin{array}{c}\\ \\ \end{array}\right]\hspace{0.17em}\begin{array}{c}\mathrm{⎫}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎬}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎪}\\ \mathrm{⎭}\end{array}\mathit{.}\end{array}$Figure 11 shows an all-sky mean field (Wiener-filtered) map of the curvature variable ζ constructed in (densely-packed) shells from the multipoles ζ_LMκ(χ). The figure actually shows ζ projected onto differential visibility, w = de^{− τ}/ dχ. The map ζ_w(θ,φ) is the spherical transform of ζ_LMκw = ^∫w(χ)dχζ_LMκ(χ). The top left map is ⟨ ζ_w | T ⟩ using the all-sky maps SMICA DX11. The middle left map is ⟨ ζ_w | E ⟩, based on just the polarization information. (The E-maps did not have high pass filtering imposed, so that an indication of the full results obtainable with E alone can be seen.) The bottom left panel shows ⟨ ζ_w | T,E ⟩. Visually the combined two-mode constraint has more detail than either of the one-mode constrained maps, enhancing the estimation of ζ_w. In each of the three cases, two sample realizations are shown which include allowed fluctuations about the mean map. These illustrate the level of uncertainty in the patterns found in the mean maps, necessary ingredients to the full statistical description. Note that the fluctuations about the mean are larger for T or E alone than for the combined T,E.

In all cases, the fluctuations become very large indeed if we allow for arbitrary projections, not conditioned by the visibility. Another approach similar to the differential visibility projection is to make a slice at χ = χ_dec; i.e., making the projection a delta-function in χ at the peak where dw/ dlnχ = 0, which defines χ_dec. This single slice case is shown in Fig. 12. The figures look quite similar in structure to the projected differential visibility plots. Although χ ~ χ_dec is where most information is, the fluctuations of a single slice are somewhat larger than for the differential visibility projection. The differential visibility includes a subdominant “reionization bump”, with a small weight focused at low L where sample variance is large; restricting the projection of ζ to be just over the reionization region results in a low L mean map that is largely swamped by the allowed fluctuations about it. The fluctuations about the mean become much larger if we project over all χ, since there are vast terrains in χ in which there is very little CMB temperature or polarization information; in particular between recombination and reionization. For that uncharted territory, a realization of the reconstructed ζ-map reverts to a realization of the fiducial ζ-power.

Figure 10 quantifies where the information resides in L − χ space. Here, the signal variance is the ensemble average of the square of the mean field, assuming the ζ are drawn from a Gaussian with fiducial covariance, $\mathrm{⟨}\mathrm{⟨}{\mathit{\zeta }}_{\mathit{LM\kappa }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\mathit{T}\mathrm{⟩}\mathrm{⟨}{\mathit{\zeta }}_{\mathit{LM\kappa }}\mathrm{(}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{)}\hspace{0.17em}\mathrm{|}\hspace{0.17em}\mathit{T}\mathrm{⟩}\mathrm{⟩}\mathrm{=}{\mathit{C}}_{\mathit{L}}^{\mathit{\zeta T}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathrm{[}{{\mathit{C}}_{\mathit{L}}^{\mathit{TT}}}^{\mathrm{]}}}^{-1}{\mathit{C}}_{\mathit{L}}^{\mathit{T\zeta }}\mathrm{\left(}{\mathit{\chi }}^{\mathrm{\prime }}\mathrm{\right)}\mathit{,}$(56)and the noise (i.e., fluctuation) variance is the covariance matrix, as given by Eq. (54) for T and the equivalent for E. Figure 10 plots the signal-to-noise at each individual χ slice (hence ignores the components off-diagonal in χ): $\mathrm{[}\mathit{S}\mathit{/}\mathit{N}{\mathrm{]}}_{\mathit{L}}^{\mathit{T}}\mathrm{(}\mathit{\chi }\mathrm{)}\mathrm{=}\frac{\mathit{\rho }}{\sqrt{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{\rho }}^{\mathrm{2}}\mathrm{)}}}\hspace{0.17em}\mathrm{;}\mathit{\rho }\mathrm{\equiv }\frac{{\mathit{C}}_{\mathit{L}}^{\mathit{\zeta T}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}}{\sqrt{{\mathit{C}}_{\mathit{L}}^{\mathit{TT}}{\mathit{C}}_{\mathit{L}}^{\mathit{\zeta \zeta }}\mathrm{\left(}\mathit{\chi ,\chi }\mathrm{\right)}}}\hspace{0.17em}$(57)for T, with an equivalent expression for E. The signal-to-noise structure in χ can be contrasted with the differential visibility plotted in the bottom panel: Fig. 10 shows that the L–χ information has a reach beyond the differential visibility structure, especially for low multipoles, because the associated waves can straddle the last scattering surface. The “reionization bump” in signal-to-noise seen in Fig. 10 is not that prominent. The Integrated Sachs Wolfe impact on the S/N is evident, but is relatively low, with the consequence that we cannot draw out high significance results from the ISW effect for ζ reconstruction (or equivalently for gravitational potential reconstruction).

The top panels in Fig. 10 are for an ideal experiment, with no noise (apart from the cosmic variance “noise” in the Wiener map fluctuations). The bottom panels of Fig. 10 include realistic Planck noise, as estimated from the FFP8 simulations: the high signal to noise in E in the top panels is noticeably diminished around L ~ 100 over what a cosmic variance limited experiment would give.

Figures 11 and 12 are Gaussian-filtered on a relatively large L ~ 200 scale (40′ FWHM). To see what happens at higher resolution, Figs. 13 and 14 zoom in on a typical 20° × 20° patch, with long waves removed using a filter W_L = sin²(L − L_c) / ΔL for L_c<L<L_c + ΔL. The specific choices are L_c = 20 and ΔL = 20, the values used for the Planck high pass for polarization maps. To allow direct comparison, we have just done the same high pass filtering for the T map. We also removed the means of the maps. The resolution is L ~ 400 (20′ FWHM). Figures 13 and 14 illustrate that the fluctuations play a larger role in this higher resolution regime. Note that the χ_dec slice has about the same fluctuation level as the differential visibility projection.

Lensing effects are not taken into account in these ζ-maps. In principle one could de-lens the temperature and polarization maps before forming the Wiener filter. In practice this would be a highly noisy operation with current Planck data. Hence the contaminating influence of lensing on these ζ-reconstructions would be treated by comparing simulated ζ-maps with and without lensing. Such corrections are expected to be a subdominant bias, as it has been modelled in the rest of the paper.

7. Validation of Planck results

In this section, we perform a battery of tests aimed at verifying the robustness of the results obtained in the previous section. Table 10 shows excellent agreement with our 2013 analysis of nominal mission data (Planck Collaboration XXIV 2014). The agreement using different component separation methods in temperature is also generally very good. Our focus here is thus on polarization bispectra. Redundancy is perhaps the most important element in our analysis, as far as robustness is concerned. We devote considerable attention to comparing the outcomes of different estimators and component separation pipelines, and assess their level of internal consistency. We also verify the stability of our results in the harmonic and pixel domains, by considering different sky cuts and multipole intervals. Given the large computational requirements of these tests, and since results from different optimal estimators agree very well, as shown in previous sections, we will variously use the KSW, binned, or modal pipeline for different tests. In doing this we will also exploit the complementarity of different decompositions, which might make some of them more suited for different tests than others (for example, the binned pipeline directly works with a harmonic space decomposition of the bispectrum, thus making it perfectly suited for tests of ℓ-dependence, the modal pipeline can perform quick model-independent tests by working on a relatively small subset of bispectrum modes, and so on).

7.1. Dependence on foreground cleaning method

7.1.1. Comparison between f_Nl measurements

In Table 10 we show f_NL results for the local, equilateral and orthogonal shapes, using four different optimal estimators, and four different foreground cleaning pipelines. The agreement between different estimators, on a given map, is within a fraction of a standard deviation, in line with theoretical expectations and simulations, as reported in Sect. 5. This level of agreement applies to all of the T+E, (TTT), and (EEE) bispectra.

The overall picture becomes more complex when comparing outputs across different foreground cleaning methods and estimators. Whereas for TTT and T+E results the agreement also seems quite good in this case (being at the level of half a sigma or better, for nearly all combinations of cleaned maps and shapes), larger discrepancies are present in the EEE bispectrum measurements. The most notable differences are found for the equilateral shape, where SMICA and NILC find values of f_NL consistent with 0 within 1σ, while SEVEM and Commander measure a roughly 2σ deviation from Gaussianity. The largest discrepancy is found for the pair Commander–NILC, using the binned pipeline (see Table 10). This estimator recovers ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{369}\mathrm{\pm }\mathrm{160}$ using the Commander E-map and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{97}\mathrm{\pm }\mathrm{141}$ using the NILC E-map. Other pipelines, and different choices of component separation methods, show slightly smaller but similar discrepancies, at a level of about 1.5σ. The same shape and estimator analysis of temperature maps shows good agreement: Commander recovers ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{-}\mathrm{36}\mathrm{\pm }\mathrm{73}$, while NILC gives ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{-}\mathrm{45}\mathrm{\pm }\mathrm{71}$. The combined T+E measurement, still for the same modal pipeline and equilateral shape, yields ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{26}\mathrm{\pm }\mathrm{50}$ for Commander and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{-}\mathrm{4}\mathrm{\pm }\mathrm{46}$ for NILC, corresponding to about half a standard deviation difference. This general trend is seen for other shapes and estimators.

Simulations were used to give insight into the expected level of disagreement. For each of the four component separation methods, we generated 100 FFP8-based Gaussian simulations with realistic beams and noise levels. These simulations start from the same initial single frequency realizations, and are processed through the four different foreground cleaning pipelines. The starting maps do not include any foreground component (the same map generation procedure is used in the Monte Carlo determination of error bars). The differences in final simulations are thus generated only by the different data filtering and coadding operations performed either in pixel, harmonic, or needlet domains by the various foreground cleaning methods, and by additional manipulations of the maps which are required for f_NL estimation, such as inpainting. Therefore, the average scattering in f_NL, measured from these realizations, provides us with a baseline assessment of the expected discrepancies between different methods when foreground residuals and other spurious sources of NG are negligible. We can then compare them with differences observed on the data to establish whether they are consistent with expectations, or are too large. The latter would raise the concern that foreground contamination, or other systematics, might be affecting the results.

Results are shown in Table 12, for EEE and T+E and two different sky coverages. The scatter between f_NL values from simulations is about 0.5σ for both T+E and EEE. This is smaller than the differences in the Planckf_NL values obtained from EEE analysis of different foreground-cleaned maps, especially for the equilateral shape. However, for the final combined T+E measurement, observed differences are in good agreement with expectations from simulations for the majority of cases. Another important point is that the consistency shown in Table 12 for T+E measurements is stable to the change of sky coverage (in polarization) from f_sky = 0.74 to f_sky = 0.64. This will be confirmed by additional tests later in this section. For the SMICA–SEVEM pair we also verified stability using an even larger mask with f_sky = 0.52.

Residual foregrounds may be responsible for at least some of the observed excess of scatter in EEE-derived f_NL between different cleaning algorithms. This is supported by the fact that several EEE results for this test change significantly for different masks, and that discrepancies are alleviated by using a larger mask, especially for equilateral shapes (see e.g., SEVEM–NILC and SMICA–SEVEM in Table 12). However, modal coefficients and their correlations are stable to a change of mask (see below), as are values of f_NL for a given component separation method (see Table 13).

Another possible contributor is a mismatch between the noise model (used to build the estimator normalization, weights and linear term), and the actual noise in the data. Polarization data are very noisy, and it is a known problem that the model assumed underestimates the true noise. This means that the error bars for EEE f_NL results, quoted in Table 10, are somewhat underestimated, which does not seem to be a problem for the final T+E results, since the weight of the EEE bispectrum in the final combined measurement is very low. This is confirmed by the results of this test. Indeed, we investigate EEE in detail because it is a useful and sensitive indicator of various systematics in the polarized maps (which could eventually leak into the TTE and TEE bispectra), rather than for its statistical weight in the final measurement. It is then fair to say that issues in the EEE bispectra, and related f_NL measurements are not yet fully understood and will require further investigation. Even though the T+E are consistent, we recommend that results that include polarization data are regarded as preliminary at this stage.

7.1.2. Comparison between reconstructed bispectra

It is important to stress that the conclusions reached at the end of the previous subsection refer to the three main bispectra in our analysis, defined by the standard scale-invariant local, equilateral, and orthogonal primordial shapes. These shapes select a specific subset of configurations in the overall bispectrum domain (essentially squeezed, equilateral, and flattened triangles). Therefore, testing consistency between methods for these shapes does not guarantee that results for the many other NG models considered in this work (such as, for example, the oscillatory bispectra of Sect. 8) will display the same level of agreement. For this reason we decided to perform a model-independent test of consistency between methods, based on comparisons between the β_n eigenmodes used for bispectrum reconstruction in Sect. 6.2. We also reconstruct the bispectrum starting from a binned ℓ-decomposition, and this will be used in Sect. 7.4 to study stability of the results in the harmonic domain. For the β_n study we consider a simple test based on measuring the correlation coefficient between modes extracted from different foreground-cleaned maps. The correlation is defined, as usual, by ${\mathit{r}}_{\mathit{ij}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{cov}{\mathrm{(}{\mathit{\beta }}_{\mathit{n}}^{\mathit{i}}\mathit{,}{{\mathit{\beta }}_{\mathit{n}}^{\mathit{j}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathrm{\left(}{\mathit{\sigma }}_{\mathit{n}}^{\mathrm{2}}{\mathrm{)}}_{\mathit{i}}\mathrm{\right(}{\mathit{\sigma }}_{\mathit{n}}^{\mathrm{2}}{\mathrm{)}}_{\mathit{j}}}\mathit{,}$(58)

and we measure it for each combination of the SMICA, SEVEM, NILC, and Commander maps, labelled by the indices i,j. Results are given in Tables 14 and 15 for the two modal pipelines and are illustrated in Fig. 15. These results show an excellent degree of correlation between different maps in temperature (especially for SMICA, SEVEM, and NILC), which is reduced when polarization is considered. In fact the correlation for polarization is not much lower than temperature for SMICA and NILC, while it reduces the correlation for the pairs SMICA–SEVEM, and NILC–SEVEM, and for Commander when paired with any other method. This is consistent with previous findings of our f_NL-based test. To test if these results are due to foreground residuals (or other effects that are not included in the simulations), we evaluate the same mode-mode correlations on the same sets of 100 realistic, foreground-free, Gaussian simulations as previously used, and processed through each of the different component separation pipelines. For this analysis we consider TTT and EEE bispectra, expanded via the low-resolution modal estimator. Our results are reported in Table 14, in the simulation column, and they clearly show that the trend in the simulations is consistent with what we see in the Planck data. In particular, EEE results show a lower degree of correlation in simulated maps, for the same pairs of methods. The observed loss of correlation in polarization does not seem to come from unresolved foregrounds or other unaccounted systematics, but rather something intrinsic to the foreground-removal algorithms. They are substantially different, as SMICA and NILC both perform the cleaning in harmonic space, at the level of E and B multipoles, whereas SEVEM is essentially a pixel-space template fitting method, performing the subtraction on Q and U maps, or inpainted before f_NL estimation. These issues will be studied in greater detail in future work, using Wiener-filtered, as well as inpainted maps for f_NL estimation. However, we have already seen that the larger scatter between modes from different foreground cleaning methods does not have a serious impact on f_NL estimation, at least for the standard local, equilateral, and orthogonal shapes. The non-standard shapes need to be analysed separately to check robustness of NG polarization results. This is the approach we will take for the various non-standard NG models.

7.2. Dependence on sky coverage

For each of the four component separation methods, we have used two different polarization masks, namely the same polarization mask as in Sect. 6, with f_sky = 0.74 (defined as the polarization “common mask” in Sect. 3.4), and an extended mask with f_sky = 0.64. The temperature mask is kept unchanged in this test, and it covers a sky fraction f_sky = 0.76 (the temperature “common mask” of Sect. 3.4). We report the variation in f_NL for the three standard shapes in Table 13, which shows insensitivity to f_sky, in agreement with earlier results on T+E. In this case, however, the EEE results also seem quite stable, supporting the view that foreground residuals are not affecting our local, equilateral, and orthogonal f_NL results, especially for the final, combined T+E measurements. Tests on FFP8 simulations including foregrounds (see Sect. 7.3) suggest that f_NL measurements obtained from the SMICA and SEVEM maps are the most accurate under the current choice of mask. As a further check of these two methods we consider a third polarization mask, with f_sky = 0.53, and repeat the combined T+Ef_NL measurement, also finding stable results. For SMICA we find ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{5.6}\mathrm{\pm }\mathrm{5.4}$, ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{65}\mathrm{\pm }\mathrm{58}$, and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}\mathrm{=}\mathrm{-}\mathrm{30}\mathrm{\pm }\mathrm{26}$, while for SEVEM we obtain ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{9.4}\mathrm{\pm }\mathrm{5.4}$, ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}\mathrm{=}\mathrm{75}\mathrm{\pm }\mathrm{59}$, ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}\mathrm{=}\mathrm{-}\mathrm{50}\mathrm{\pm }\mathrm{30}$.

We also perform model-independent checks by looking at the correlation coefficient between different sets of bispectrum modes, in a similar way to Sect. 7.1.2, but now changing the polarization mask. Results are reported in Table 16, and confirm that the data and simulations behave similarly, and that polarization modes display a lower correlation level than temperature.

7.3. Tests on simulations

We consider two realistic data simulations, one of which is Gaussian, while the other includes local NG. We start with a foreground-free realization, add foregrounds according to the Planck Sky Model, and finally process through the four component separation pipelines. By estimating f_NL in the input foreground-free simulation, for each method, and comparing to f_NL recovered from the cleaned maps (or with the input local f_NL, for the NG case), we can assess both the impact of foregrounds on our measurement before subtraction and which method gives the highest accuracy. The necessity to clean is very apparent in the middle set of columns in Table 17, where no cleaning has been performed.

SMICA and SEVEM give the best results, both in the G and NG tests. In the G test, reported in Table 17, SMICA results show an agreement between the input and the cleaned map at the level σ_fNL/ 2 for all shapes, and for all of TTT, EEE, and T+E. SEVEM displays a similar level of accuracy, except for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$, where the difference is larger, but within one standard deviation. NILC and Commander clearly perform worse for the local shape, with NILC showing a 2σ_fNL difference, and Commander even larger than that. In the NG test, reported in Table 18, SEVEM gives the most accurate results, recovering the input with σ_fNL/ 2 accuracy or better. Results for SMICA are accurate at the 1σ_fNL level for the TTT constraint, and worse (about 2σ_fNL) in EEE. However, the combined T+E measurement is again very good (σ_fNL/ 2). NILC is also performing very well in TTT and T+E, while the EEE result is more than 2σ_fNL off.

The test described here has several limitations, the main and most obvious one being that it has been performed on just two maps (due to lack of availability of a large sample of this type of simulation at this stage). Another clear issue is that some methods, in particular Commander, seem to perform much better on actual data than on these simulations. On the other hand some important trends, observed in the data in previous tests, are clearly reproduced here, like the good stability of SMICA and SEVEM, especially for the combined T+E results and, most notably, the fact that the clear degradation in the accuracy of the EEE measurement for some methods does not seem to propagate to T+E.

7.4. Dependence on multipole number

In this section we discuss another stability test of our results, namely the dependence of the results for f_NL on the maximum and minimum multipole number used in the analysis. This test is most easily performed with the binned bispectrum estimator, since it gets the dependence of f_NL on ℓ for free with its standard analysis, simply by leaving out bins in the final sum when computing f_NL (the binned equivalent of Eq. (39)).

The dependence on ℓ_max of the results for the three standard primordial shapes (local, equilateral, and orthogonal modes), is shown in Fig. 16, for T-only, E-only, and full T+E. As mentioned in Sect. 6, the KSW and binned estimators use ℓ_max = 2500 for temperature, while the modal estimators use ℓ_max = 2000. As can be seen in the figure, both the T-only and T+E results are basically unchanged between ℓ = 2000 and ℓ = 2500 for all three shapes, showing that this difference has no impact on the results (as was to be expected from the excellent agreement between estimators in Table 10). In fact, the values of f_NL for T and T+E are reasonably stable (given their error bars) down to much lower values of ℓ_max, of about 1000.

On the polarization side we can draw a similar conclusion. The binned estimator uses ℓ_max = 2000 for polarization, while the other estimators use ℓ_max = 1500, but we see that results for E remain basically unchanged between ℓ = 1500 and ℓ = 2000. Central values and error bars for E for all three shapes have clearly converged by ℓ = 1500, and are in fact reasonably stable down to much lower values of about 700.

As we noted in the 2013 analysis (Planck Collaboration XXIV 2014), when going to the much lower WMAP resolution of ℓ_max ≃ 500, we agree with the slightly high value of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ that the WMAP team reported (Bennett et al. 2013). This value is also confirmed when including polarization. One can clearly see the value of the higher resolution of Planck.

The dependence on ℓ_min is shown in Fig. 17. Here all estimators used the same values, ${\mathit{\ell }}_{\min }^{\mathit{T}}\mathrm{=}\mathrm{2}$ and ${\mathit{\ell }}_{\min }^{\mathit{E}}\mathrm{=}\mathrm{40}$. As explained in Planck Collaboration VII (2016), not all systematic and foreground uncertainties in the low-ℓ HFI polarization data have been fully characterized yet, and hence it was decided to filter out these data.

For equilateral and orthogonal shapes the values for f_NL and their error bars are quite stable as a function of ℓ_min up to about ℓ = 100 (and ℓ ≃ 300 for E-only), which is not surprising, since these templates have little weight at low ℓ. The local template, on the other hand, depends very strongly on the lowest multipoles, which is reflected in the very rapid growth of the error bars when ℓ_min increases. Looking at T-only and T+E we see a very similar pattern, with ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ being reasonably stable, although there are some jumps. The local result for E-only wanders a bit more outside of the ± 1σ region, in agreement with the other tests in this section, which also indicate that E-only is not as stable as T-only and T+E. However, that is still quite acceptable, given the small weight of E-only in the full T+E result.

One can work out quite generally that when Y is a subset of a data set X, and P_X and P_Y are the values of a parameter P determined from these two data sets, then the variance of the difference P_Y − P_X is equal to | Var(P_Y) − Var(P_X) | Hence we can determine how likely the jumps in f_NL as a function of ℓ_min are. It turns out that the jump in the T-only value of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ between ℓ_min = 40 and ℓ_min = 53 is a 2.5σ effect (using the values of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ before subtraction of the ISW-lensing bias, which also depends on ℓ). Similarly, the jump in the T-only value of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}$ between ℓ_min = 154 and ℓ_min = 211 is (by chance) also a 2.5σ effect. Given the fact that there are 57 bins, having such a jump appears to be consistent from a statistical point of view.

7.5. A directional analysis with a needlet-based modal estimator

The validation tests on simulations, described in Sect. 7.3, point to SMICA and SEVEM as the best foreground cleaning methods for f_NL estimation. Results in Table 10 show that SMICA also has slightly smaller error bars, thus making it the method of choice for our final results.

As a further check of residual foreground contamination in the SMICA map, in this section we investigate the possible directional dependence of SMICA-derived third-order statistics by means of a needlet-based modal estimator (i.e., an estimator based on the decomposition described in Sect. 3.2, and references therein, where we use cubic combinations of needlets as our basis modes). In other words, we analyse the behaviour of the needlet bispectrum (see Lan & Marinucci 2008; Rudjord et al. 2010; Donzelli et al. 2012) on separate patches of the sky, and we study the fluctuations of the corresponding residuals.

Rather than assuming a specific anisotropic model, we instead calculate the contribution to the local f_NL from different regions of the sky and look for evidence of anisotropy in the result¹⁵.

Our modal needlet estimator has been validated with respect to the procedures considered in Sect. 5, showing excellent agreement. Since in this paper we use the needlet estimator only in this section, and as a diagnostic tool, we will not explicitly report the outcome of these validation tests here, for the sake of conciseness. The advantages of using needlet-based modal estimators have been advocated in Lan & Marinucci (2008), Rudjord et al. (2010), Fergusson et al. (2010a, 2012); we refer to these papers for more discussion and details. In short, however, they can be summarized as follows:

• 1.

  it is possible to achieve a strong data compression, i.e., to investigate cubic statistics by means of a small number of modes (needlet frequencies) for many different bispectrum templates;

• 2.

  needlet transforms have good correlation properties in pixel space, which allows study of the pixel contribution to the f_NL signal for the templates under investigation, by treating different directions independently.

In our analysis, we first divide the sky into several large “regions”, with boundaries defined by the pixels of a HEALPix grid (Górski et al. 2005). at lower resolution than the starting map (which is at N_side = 2048). For the low resolution grid we consider N_side values between 2 and 8. For each pixel in the coarse grid, we then compute the local f_NL using our modal needlet estimator, and neglect contributions from external regions.

The correlation matrices between f_NL measurements in different regions were computed via Monte Carlo simulations. N_side = 4 was chosen, as providing the best tradeoff between having a large number of regions for directional analysis (i.e., the total number of pixels in the low resolution grid), and a low correlation between different regions. This is shown in Fig. 18. It is readily seen that, at N_side = 4, the correlation is largely concentrated in one or two points near the main diagonals (where it is still low, never exceeding 34%), and falls off rapidly for all other pixels. Note that the results here refer to temperature only. The EEE local polarization error bars are large, even for the full sky analysis, making this directional approach uninformative for Planck polarization data. We concentrate on the TTT bispectrum here, complementing other validation tests in this section, which are mostly focused on polarization.

Having obtained our correlation matrices, and having shown that different regions are essentially uncorrelated, we can then proceed to extract f_NL for each region in the actual SMICA map. As a test of directional-dependent contamination of the f_NL measurement, we can also compare our results, region by region, with the fluctuations expected by looking at the diagonal of our Monte Carlo correlation matrix.

The results of this analysis are shown in Fig. 19. In the left panel, we represent the directional local f_NL map, extracted with this method. In the right panel, we report the f_NL values, region by region, and compare them to expectations from simulations. The red line gives the expected standard deviation, while the blue one gives estimates on the component-separated maps with the Monte Carlo error bars. Our estimator is normalized in such a way that the sum of all these contributions would yield exactly the f_NL estimator for the full map, so these results can be viewed as a partition of the estimates along the different directions. It is readily seen that no significant fluctuation occurs, so that our results are consistent with the absence of directionally-dependent features (which could occur due to, for instance, residual foreground contamination). As an additional check, we have also investigated the possible presence of a dipole in these data, and found that our results are consistent with Gaussian isotropic simulations.

7.6. Summary of the main validation results

Throughout Sect. 7 we have shown a battery of tests aimed at evaluating the robustness of our data set, from the point of view of bispectrum estimators, focusing especially on the polarizaton part. We studied the stability of our results (local, equilateral, and orthogonal f_NL measurements, model-independent bispectrum reconstruction) under a change of sky coverage, multipole range, and choice of component separation methods. We also considered simulated data sets and studied the ability of different component separation methods to recover the input f_NL after foreground subtraction. Our main conclusions from these tests can be summarized as follows:

• TTT and T+E results are stable both in the pixel and harmonic domains, for different component separation methods. For SMICA, we also checked that TTT temperature constraints on local f_NL show no evidence of a directional variation via a needlet-based analysis.

• SMICA and SEVEM perform better than NILC and Commander at recovering the original f_NL in foreground-cleaned simulations. At the same time, SMICA allows slightly better constraints on f_NL than SEVEM, due to a lower (by a small amount) noise level.

• EEE bispectra, and related f_NL measurements, have some problems, and do not pass all the tests. Different component separation methods show a low level of consistency (especially when comparing pixel-based cleaning methods to harmonic-based cleaning methods). This disagreement is only partly alleviated by choosing a larger E-mask, so that residual foregrounds do not seem to fully explain all issues. An important caveat, already pointed out previously, is that the noise model in polarized FFP8 simulations is known to underestimate the actual noise level in the data, leading to some degree of underestimation of the error bars in our EEE results. We stress again, however, that this has little impact on the final T+E constraints, due to the high noise level and consequent low statistical weight of EEE bispectra. This was verified in detail both on data and on simulations.

In light of the above analysis, we conclude that, as far as bispectrum estimation is concerned, our best cleaned map is the one produced by SMICA, in line with our previous 2013 analysis. We also conclude that our TTT-based f_NL constraints, summarized in Tables 10 and 11, are robust. Joint T+E constraints pass all our validation tests. On the other hand, in the light of the remaining issues in the EEE bispectra and in the FFP8 polarized simulations, as we stressed at the end of Sect. 6.1, we suggest that all measurements that include polarization data in this paper should be regarded as preliminary.

8. Other non-Gaussianity shapes for f_NL

This section discusses new searches for NG beyond standard single-field inflation. The focus here is on extensions to the analysis undertaken in Planck Collaboration XXIV (2014) with new limits on isocurvature models, further oscillatory models over a broader frequency range, and parity-violating tensor NG. However, we also briefly revisit all the non-standard models constrained in Planck Collaboration XXIV (2014), including effective field theory, non-Bunch Davies, and directionally-dependent models, in particular noting the impact of the new (preliminary) polarization data on the previous constraints.

8.1. Isocurvature non-Gaussianity

Here we show the results obtained for a study of the isocurvature NG in the Planck 2015 SMICA map using the binned bispectrum estimator. As explained in Sect. 2.3, we only investigate isocurvature NG of the local type, and in addition always consider only one isocurvature mode (either cold dark matter, neutrino density, or neutrino velocity isocurvature) in addition to the adiabatic mode. Note that the baryon isocurvature mode behaves identically to the cold dark matter one, only rescaled by factors of Ω_b/ Ω_c, so there is no need to consider it separately. In that case there are six different f_NL parameters: a purely adiabatic one (a,aa), which corresponds to the result from Sect. 6), a purely isocurvature one (i,ii); and four mixed ones (see Sect. 2.3 for an explanation of the notation).

The results are given in Table 19¹⁶ and show no clear signs of any isocurvature NG. There are a few values that deviate from zero by up to about 2.5σ, but such a small deviation, in particular given the large number of results, cannot be considered a detection. We do see that many constraints are tightened considerably when including polarization, by up to the predicted factor of about six for the cold dark matter a,ii, i,ai, and i,ii modes in the joint analysis. As discussed in detail in Sect. 7, results including polarization data should be considered preliminary, and that is even more important here, since these results depend so strongly on the additional information from polarization.

In the results so far we allowed for a possible correlation between the isocurvature and adiabatic modes. However, if we assume that they are completely uncorrelated, with a zero cross power spectrum, then there are only two f_NL parameters, the a,aa and i,ii ones. In Table 20 we give the results for this uncorrelated case. The independent results are the same as in the previous table, while in the joint results one can clearly see the difference between the neutrino density mode (the bispectrum template of which has a large overlap with the adiabatic one), and the cold dark matter and neutrino velocity modes (with templates that are very different from the adiabatic one). Again there is no evidence for any isocurvature NG; the almost 2σ result for the neutrino velocity isocurvature mode in the temperature-only case does not survive the addition of polarization.

8.2. Feature models

An important and well-motivated class of scale-dependent bispectra is the feature model, characterized by linear oscillations described by Eq. (15) and its variants in Eqs. (16) and (18). In Planck Collaboration XXIV (2014) we performed an initial search for a variety of feature models using the Modal estimator. This earlier search was limited to ω< 200 by the native resolution of our implementation of the Modal estimator (using 600 modes), roughly the same range as the initial WMAP bispectrum feature model searches at lower precision (with only 50 eigenmodes (Fergusson et al. 2012). Note, in the previous Planck release we used wavenumber k_c in line with the theory literature, but here we switch to frequency ω, in line with more recent observational power spectrum searches; the two are related by ω = 2π/ 3k_c. With the improved estimator resolution (now using 2001 modes) we are able to achieve convergence over a broader range up to ω = 350. We perform a frequency scan of 350 sampling points between ω = 10 and ω = 350, i.e., 35 independent frequencies and 10 phases. We also extend the number and variety of feature and resonance models that are investigated, essentially probing the resolution domain in which we have obtained a reliable Modal bispectrum reconstruction (see Fig. 4).

Constant feature ansatz: for the constant feature shape of Eq. (15), we can extend the frequency range much further with another approach. As the bispectrum in Eq. (15) is separable, it allows the construction of a KSW estimator (Münchmeyer et al. 2014) for direct bispectrum estimation at any given frequency. The bispectrum can be written as a sum of sine and cosine components which can be estimated separately (equivalent to measuring the amplitude and phase above) and this method was used to constrain frequencies up to ω = 3000. The range where the two estimators overlap provides validation of the two methods and excellent agreement was seen (see Fig. 20).

Apart from cross-validation with two estimators, we have undertaken further tests to determine the robustness of the results to foreground and noise effects. In Fig. 21 (left panel), we show the effect on feature model results of the subtraction of the simple point source bispectrum, as well as the ISW-lensing bispectrum. This study was a major motivation for adopting the more conservative “extended” common mask, because the consistency between different component separation methods improved markedly for low frequencies, with the original common mask requiring much larger point-source subtractions (e.g., for NILC subtraction the maximum raw significance has reduced from σ = 4.0 to σ_clean = 2.2 at ω = 110). After cleaning these signals, the SMICA, SEVEM and NILC results are in good agreement, and also consistent with each other when polarization is included (while the Commander results generally have a larger variance and so are not included in the plotted averages). Fortunately, the effect of subtracting ISW-lensing and point source bispectra diminishes rapidly at higher frequencies ω> 200 and should be negligible; subtraction was only undertaken in the Modal region ℓ< 350. In Fig. 21 (right panel), we show the effect on the averaged significance of reducing the Planck domain from the usual ℓ_max = 2000 to ℓ_max = 1500 (ℓ_max = 1500 to ℓ_max = 1000 for E-modes). Despite the non-trivial change in overall signal-to-noise entailed, there is no strong evidence for an ℓ-dependent signal, as might be expected if there was substantial NG feature contamination in the noise-dominated region. Finally, we note that most peaks at low ω show some correlation between T-only and T+E, although there are notable exceptions, such as the peak at ω = 180, which is removed after inclusion of polarization (see also the phase plots in Fig. 20 before marginalization). The temperature feature peaks observed in Fig. 21 at ω ≈ 110, 150, and 180 are consistent with the peaks identified previously in Planck Collaboration XXIV (2014).

In Fig. 20, the full set of frequency results, 0 <ω< 3000, for the constant feature model in Eq. (15) is shown for both the Modal and KSW estimators. There is good agreement for regions where the estimators overlap and PS and ISW lensing bispectra subtraction is not necessary (ω> 200). Generally, there is tighter consistency between temperature-only results, than with polarization where there is additional scatter between foreground separation methods. Scanning across the full frequency range, there is no strong evidence for any large maximum that might indicate unequivocal evidence for a feature model signal. The maximum peak significance obtained with either T-only or T+E is consistent with expectations for a Gaussian model over this frequency range. In particular, for the KSW estimator using the SMICA data, the highest significance found in the range 200 <ω< 3000 is 3.2σ in T-only and 2.9σ in T+E. To gauge the likelihood of these results occurring randomly, realistic Gaussian SMICA simulations were analysed with the KSW estimator and found to typically produce a highest peak with 3.1( ± 0.3)σ over the same frequency range.

Generalized feature models: we have also deployed the Modal estimator to look at (non-separable) feature models with equilateral and flattened cross-sections, as motivated by varying sound speed scenarios and those with highly excited states, respectively. In the left panels in Fig. 22 we show results from the equilateral feature model of Eq. (16), including the frequency/phase contours before marginalization for the SMICA T-only map. Multiple peaks are apparent in the temperature signal across the Modal range up to an average maximum 3.3σ raw significance. However, from the lower panel it is clear that the polarization signal is not well correlated with the temperature in the equilateral case, generally reducing peak heights with the maximum now about 2.6σ (while eliminating the ω = 180 peak altogether). This temperature peak remains present with the ℓ_max = 1500 cutoff, where the signal is slightly higher, but the polarization in this case is less well correlated (using ${\mathit{\ell }}_{\max }^{\mathit{E}}\mathrm{=}\mathrm{1000}$). For Gaussian noise we would not expect polarization to reinforce a high temperature signal on average. It may also be that the equilateral temperature signal has some residual diffuse point-source contamination. The equilateral feature model is the most affected by the removal of point sources, so the presence of a more complex correlated PS bispectrum (not removed by the constant PS template subtraction) remains for future investigation.

Results for the flattened feature model in Eq. (18) are shown in Fig. 22 (right panels), displaying more coherence between temperature and polarization. The temperature signal with a 2.6σ peak between 50 <ω< 150 is reinforced by polarization and merges to make a broad 3σ peak around ω = 90, together with another distinct peak at ω ≈ 140. Such broad frequency peaks are not expected, because neighbouring feature models should be nearly uncorrelated over a range Δω_eft ≈ 13 (as we discuss below). As the phase plots in Fig. (22) indicate, this breadth in frequency ω may reflect the neighbouring feature models adjusting phase φ to match an underlying NG signal of a related, but different, nature. We note also that the frequency region for ω< 100 is susceptible to some degeneracy with cosmological parameters. We shall consider a “look-elsewhere” statistical analysis of these results below.

Single-field feature solutions: we have also searched for the specific analytic solutions predicted for single-field inflation models with step-like potential features, as given in Eqs. (19) and (20), with results shown in Fig. 23. The highest peaks for the K²-cosine model occur around 2.5σ with temperature-only, then rises to 2.7σ when polarization data are included, again with peaks at other distinct frequencies apparent. The K-sine model shows a similar apparent signal level, with a maximum T-only 2.7σ peak, dropping to 2.4σ with polarization. One further difficulty with a positive interpretation of these bispectrum results in this low frequency range is that stronger S/N counterparts in the power spectrum are predicted for these simple models (Adshead et al. 2011), whereas no significant correlated oscillation signals are apparent at the relevant peak frequencies Planck Collaboration XX (2016).

Feature model peak statistics: in order to determine consistency with Gaussianity for these feature model results, we can apply a number of statistical tests developed for this specific purpose (Fergusson et al. 2015) and, if warranted, also apply these jointly in combination with power spectrum results, as for the WMAP polyspectra analysis (Fergusson et al. 2015). When scanning across the (ω,φ) parameter-dependent feature models, we are searching through independent models for which Gaussian noise, by chance, can lead to a large apparent signal. We must correct for this multiplicity of tests when determining the actual significance of results for a given model – this is a quantitative correction for any model with free parameters, distinct from the a posteriori choice of models to test. The simplest approach is to determine whether the maximum peak is consistent with Gaussian expectations, which can be determined from Monte Carlo simulations. However, in Fergusson et al. (2015) it was recognized that these feature models can be accurately characterized analytically with a χ-distribution having two degrees of freedom¹⁷. Taking ℓ_max = 2000 and using this analysis, the frequency step size between models that are uncorrelated is approximately Δω_eff = 13, so we have an effective number of independent feature models N_eff ≈ 27 for the Modal frequency range (with N_eff ≈ 230 across the larger constant feature survey range). Accordingly appropriate look-elsewhere corrections have been applied to find an adjusted significance for the maximum peak signal found in all the feature model searches undertaken, which is shown in Table 21. Given that this feature model survey is over many independent frequency models and combinations of data, even the highest raw significances above 3σ (“Raw” column in Table 21) are reduced to a corrected significance below 2σ (“Single” column). Hence, there appears to be no evidence from maximum peak statistics for feature model deviations from Gaussianity.

Nevertheless, we also examine the possibility that multiple feature models are contributing to a NG signal, given the apparent emergence of several preferred frequency peaks. This integrated multi-peak statistic can also be accurately approximated analytically (Eq. (60)) using a χ-distribution; essentially we sum over all independent frequencies using the single peak significance adjusted for the “look-elsewhere” effect (see Eq. (59)). Most of the signal surveys exhibit an unusual number of broad overlapping peaks within the accessible frequency domain, so the multi-peak statistic does yield a much higher significance, with many models above 2σ after “look-elsewhere” correction. Notable cases are the temperature-only signal for the equilateral feature model which yields an average significance of 3.4σ across the foreground-cleaned maps with concordant bispectrum results (i.e., SMICA, SEVEM and NILC); however, this interesting multi-peak significance drops to only 1.6σ when the polarization data are included (assuming the reliability of E results). On the other hand, the flattened feature model has an average multi-peak significance of 1.7σ in temperature only, which rises to 3.4σ with polarization included (higher at 3.7σ if Commander data were to be included in the average). In this case, beyond the number of peaks, it is also their width that contributes, with the main signal around ω ≈ 90, much broader than Δω ≈ 13. Finally, after look-elsewhere effects are taken into account, the K²-cosine single-field solutions yield multi-peak statistics that rise with the inclusion of polarization data from 2.0σ to 2.5σ, while the K-sine falls from 3.0σ (T) to 1.9σ (T+E).

An interesting, but not entirely coherent, picture emerges from these searches for non-standard models in the new Planck temperature data, especially when combined with the additional (preliminary) polarization information. In Planck Collaboration XXIV (2014), we noted that the feature model searches provided interesting hints of NG. This more rigorous statistical analysis confirms this view, allowing for several alternative feature model explanations of the apparently high NG signal observed in the bispectrum reconstructions (see Fig. 4). However, there is no strong evidence for a single large feature model at a particular frequency; rather, the high multi-peak statistics indicate signal that is spread across several broad peaks. Given the variability between different feature models and polarization component-separation methods, we note the caveat that the integrated multi-peak statistic could be sensitive to calibration issues and foreground contamination. For this reason, we do not make strong claims for these non-standard signals at this stage, but we note that oscillatory models will continue to be investigated thoroughly over a wider frequency domain and using the more reliable polarization data available in the final Planck data release.

8.3. Resonance/axion monodromy models

Generalized resonance models: using the Modal expansion, we have embarked on a survey of the simplest resonance model (Eq. (10)), as well as the equilateral and flattened variants proposed in the literature, i.e., described by Eqs. (13) and (14), respectively. The raw significance for the resonance models for both temperature-only and temperature plus polarization data are shown in Fig. 24; these are the maximal results marginalized over the phase parameter φ. Previously, the resonance model was studied in Planck Collaboration XXIV (2014) using the Modal expansion over a narrower frequency range yielding no results above a raw significance of 1σ. In this extended analysis over a wider frequency range, the constant sin(log ) model (Eq. (10)) produces 2.2σ peaks for T-only, and 2.6σ for T+E. The equilateral resonance model (Eq. (13)) achieves a maximum 2.8σ in T-only at ω ≈ 35, rising to a more impressive average 3.2σ for T+E. For the flattened case (Eq. (14)) we have 2.5σ and 3.0σ, respectively at ω ≈ 12. Qualitatively, the results shown in Fig. 24, exhibiting broad peaks, are similar to those for feature models.

Resonance model peak statistics: to determine the statistical significance of these results given the look-elsewhere effect of scanning across the parameters (ω,φ), we have used the two peak statistics defined above in Eqs. (59) and (60) for feature models (Fergusson et al. 2015). In this case, the maximum peak statistic for the constant resonance model of 2.6σ (T+E) is readjusted to an unremarkable “look-elsewhere” single peak significance of 0.9σ. Likewise the apparently significant results above 3σ for the equilateral and flattened models now fall below 2.0σ with T+E. Using the single peak statistic alone, we would conclude that there is no strong evidence for any individual resonance model. Resonance models also generate oscillations in the power spectrum, and an analysis based on the 2015 temperature and polarization likelihood is presented in (Planck Collaboration XX 2016). A combined statistical treatment of resonance model power spectrum and bispectrum results will be reported in the future.

The multi-peak statistic in Eq. (60) integrates the resonance model signal across all frequencies to determine consistency with Gaussianity. The constant resonance model has a modest multi-peak signal but, like the feature models, the equilateral and flattened resonance shapes offer stronger hints. The multi-peak equilateral signal rose from 1.9σ (T-only) to 3.1σ (T+E) after adjusting for the “look-elsewhere” effect, while the flattened signal went from 2.4σ (T-only) to 3.2σ (T+E). These interesting results, reflecting those obtained for feature models, suggests the fit to any underlying NG signal might await alternative, but related, oscillatory models for a more compelling explanation.

High frequency resonance model estimator: we have further surveyed the simple resonance model (Eq. (10)) with a second approach. Using a model-specific expansion in terms of linear oscillation, proposed in Münchmeyer et al. (2015), it is possible to extend the frequency range of the analysis considerably. In this approach one exploits the fact that any bispectrum shape which is a function of k₁ + k₂ + k₃ can be expanded in Fourier modes of k₁ + k₂ + k₃, resulting in an effectively one-dimensional expansion, as opposed to the general Modal expansion. In the present implementation we use 800 sine and cosine modes, which cover the full frequency space of the power spectrum search, i.e., 0 < ω < 1100. The significances found with this method are presented in Fig. 25. As was the case for the feature model, we find that SMICA, SEVEM and NILC results are in good agreement in both temperature and polarization, while Commander results generally have larger variance and so are not included in the plotted averages. We find that the largest peak of the average is 3.6σ in temperature and 3.1σ in temperature and polarization combined. The results of the Modal expansion discussed in the previous paragraph are in good agreement with the high frequency resonance model estimator in the overlapping frequency range.

Due to the high computational demands of this analysis, we have only exactly assessed the look-elsewhere effect for SMICAT data, for which we find an expected maximum peak of 3.5σ ± 0.4σ in the case of Gaussian maps, to be compared to 3.7σ in the SMICAT data, demonstrating that the results are fully consistent with Gaussianity. The expected maximum peak for Gaussian maps was calculated from the Fisher matrix with the method described in Meerburg et al. (2016). The average over component separation methods as well as the T+E data is even less significant. For the high frequency estimator we have assessed the significance of multiple peaks in the following way. Define the multi-peak amplitude A_M as the sum of squares of the M highest significances σ_i in the frequency range, i.e., ${\mathit{A}}_{\mathit{M}}\mathrm{=}{{}^{\mathrm{(}}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{\mathit{M}}{{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}}^{\mathrm{)}}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, where only approximately independent peaks with Δω > 10 are considered. One can then compare A_M to its distribution in the Gaussian case and get an individual significance σ_M for each number of peaks M, where we assume M ≤ 10. The multi-peak statistic is then obtained by maximizing over M, leading to an additional look-elsewhere effect that we also accounted for. In this way we find that in the SMICAT data the raw peak maximum 3.7σ, is reduced to 0.5σ for the single-peak statistic and to 0.6σ in the multi-peak case. This large reduction in significance is due to the large number of independent frequencies as well as to the maximization over phases. One may argue that the frequency range ω< 1100 is too large, as EFT arguments for resonance non-Gaussianities (Behbahani et al. 2012) limit the frequency range to . As an example, we have therefore also calculated the look-elsewhere-corrected significances when we limit the analysis to ω< 250. In this case we find a single-peak significance of 0.6σ and a multi-peak significance of 0.9σ. Clearly the results are fully consistent with Gaussianity.

8.4. Equilateral-type models and the effective field theory of inflation

There is considerable interest in equilateral-type models because they are physically well-motivated, through e.g., varying sound speed scenarios. There are generic predictions available from the effective field theory of inflation, notably the two specific effective field theory (EFT1 and EFT2) shapes that give rise to the equilateral and orthogonal approximations. These models were previously constrained in Planck Collaboration XXIV (2014) and the reader is referred to Sect. 2 of that paper for analytic expressions for the specific shapes constrained here. In Table 23, we list the main equilateral-type models in the literature, giving constraints for T-only and T+E. All these models correlate well with the equilateral ansatz (Eq. (12)) and likewise do not show a significant signal. However, despite this correlation, it is interesting to note the variation between models, largely due to the difference between these shapes in the flattened limit. The implications of these results are discussed in Sect. 11.

8.5. Models with excited initial states (non-Bunch-Davies vacua)

Non-Bunch-Davies (NBD) or excited initial states are models which produce flattened (or squeezed) bispectrum shapes. The wide variety of NBD models that have been proposed are briefly classified and labelled in Sect. 2, following a more extensive overview in Sect. 2 of Planck Collaboration XXIV (2014) where more analytic forms and the first constraints were presented. The latest Planck constraints for these models are listed in Table 24, obtained using the Modal 2 estimator with polarization, Despite the apparent “flattened” signal seen in the Planck bispectrum reconstructions (Fig. 4), this is generally not matched well by the specific modulation induced by the acoustic peaks for these scale-invariant NBD models. Tight constraints emerge for most models. The largest signal obtained is from the NBD sinusoidal shape which gives a 1.6σT-only raw significance, rising to 2.1σ for T+E; this is hardly an impressive correspondence given the number of models surveyed and the parameter freedom used in maximizing the signal. However, an important caveat for NBD models is that the predicted shapes can be very narrow in the flattened limit, in which case solutions have been smoothed to match the current Modal resolution (though this has improved considerably since the Planck 2013 NG analysis). An improved match to the warm inflation shape means that the final constraint shown in Table 24 is more robust, with further implications discussed in Sect. 11.

8.6. Direction-dependent primordial non-Gaussianity

We impose observational limits on direction-dependent primordial NG parametrized by Eq. (23). Rather than using c₁ and c₂ we instead choose to work with the non-linearity parameters ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}\mathrm{=}\mathrm{1}}\mathrm{=}\mathrm{-}{\mathit{c}}_{\mathrm{1}}\mathit{/}\mathrm{4}$ and ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{2}}\mathrm{=}\mathrm{-}{\mathit{c}}_{\mathrm{2}}\mathit{/}\mathrm{16}$ (chosen to match a primordial bispectrum that is equal to the equilateral shape in the equilateral limit) keeping the notation from the 2013 results. We estimated the ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}}$ values from temperature data and high-pass filtered polarization data from the four foreground-cleaned CMB maps SMICA, NILC, SEVEM, and Commander, where we apply the common mask. The details of the KSW estimator and its derivation is presented in Appendix A. For temperature data, we use the common mask as adopted in Planck Collaboration XII (2014), which has more conservative foreground masking than the newly available mask. We choose the more conservative foreground masking, considering the fact that anisotropic NG is more sensitive to residual foregrounds. We set the maximum multipole to 2000 and 1000 for temperature and polarization data, respectively. Validating our analysis pipeline with realistic simulations, we find that the asymmetry of the Planck beam, coupled with the Planck scanning pattern, inflates the statistical fluctuations of the ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}}$ significantly. Noting the large angular scale of artificial anisotropy produced by the beam asymmetry, we set the minimum multipole to 101, and find that the statistical fluctuation of estimation from simulated data is close to the theoretical expectations.

These two shapes are also constrained using the Modal 2 estimator, which is not affected by the beam asymmetry and is used in the same form as elsewhere in the paper with multipoles from 2 to 2000 and 30 to 1500 being used for temperature and polarization, respectively. The present constraints are consistent with those found for T-only in Planck Collaboration XXIV (2014), but at higher resolution convergence has improved considerably, reflected in the lower variance.

We find that the ISW-lensing bispectrum and the unresolved point-sources bispectrum bias the estimation of the ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}}$, in particular in the analysis of temperature data. For our final values, we subtract both these biases from our estimation. In Table 25, we report the estimated value of ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}}$ from the foreground-cleaned CMB maps. For L = 1 the effect of the differing ℓ-ranges between the two estimators is not so significant and the results are quite consistent. For L = 2, which has significant signal in the squeezed configuration, the effect of removing small scales from the KSW estimator is more pronounced, resulting in significantly enlarged error bars. In light of this, the differences seen between the central values for the two methods is to be expected and does not indicate any inconsistencies between the two approaches. The slight differences between the results from different foreground-cleaned temperature maps are within the likely range of statistical fluctuations, estimated from realistic simulations of CMB and noise propagated through the pipelines of foreground-cleaned map making. As seen in Table 25, we find that the estimated values of ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}}$ from Planck 2015 temperature plus polarization data are consistent with zero.

8.7. Parity-violating tensor non-Gaussianity

We present observational limits on the parity-violating tensor nonlinearity parameter ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}$ from the temperature and E-mode polarization data. Unlike the usual scalar-mode templates, the CMB bispectra sourced from the tensor NG (Eq. (27)) are written in non-factorizable forms (Shiraishi et al. 2013b); hence, we use the separable Modal pipeline in our bispectrum estimations.

The parity-violating NG under examination induces non-vanishing signals not only in parity-even configurations (ℓ₁ + ℓ₂ + ℓ₃ = even) but also in the parity-odd ones (ℓ₁ + ℓ₂ + ℓ₃ = odd) in the temperature and E-mode polarization bispectra (Shiraishi et al. 2013b). The optimal estimator, including all (even + odd) bispectrum signals, is expressed by the linear combination of the parity-even and parity-odd estimators, reading (Liguori, in prep.) $\begin{array}{ccc}\mathit{f̂}\begin{array}{c}\mathrm{all}\\ \mathrm{NL}\end{array}\mathrm{=}\frac{{\mathit{N}}^{\mathrm{even}}\mathit{f̂}\begin{array}{c}\mathrm{even}\\ \mathrm{NL}\end{array}\mathrm{+}{\mathit{N}}^{\mathrm{odd}}\mathit{f̂}\begin{array}{c}\mathrm{odd}\\ \mathrm{NL}\end{array}}{{\mathit{N}}^{\mathrm{even}}\mathrm{+}{\mathit{N}}^{\mathrm{odd}}}\mathit{,}& & \end{array}$(61)where N^{even / odd} is the normalization factor (related to the Fisher matrix as N^{even / odd} = 6F^{even / odd}) defined for ℓ₁ + ℓ₂ + ℓ₃ = even / odd. The parity-even estimator $\mathit{f̂}\begin{array}{c}\mathrm{even}\\ \mathrm{NL}\end{array}$ can be computed using the original Modal methodology (Fergusson et al. 2010a, 2012; Fergusson 2014; Liguori, in prep.), while in computations of the parity-odd estimator $\mathit{f̂}\begin{array}{c}\mathrm{odd}\\ \mathrm{NL}\end{array}$, we follow the extended spin-weighted pipeline (Shiraishi et al. 2014, 2015; Liguori, in prep.).

Our ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}$ estimations (with both temperature and polarization data) are based on the resolution of ℓ_max = 500 and HEALPix N_side = 512, leading to feasible computational costs. This choice is not expected to change the results significantly, in comparison to the analysis at higher resolution, e.g., ℓ_max = 2000 and Healpix N_side = 2048, since the cosmic variance and instrumental noise are already far higher than the signals for ℓ ≳ 300(Shiraishi et al. 2013b). Only in the polarization data analysis is an effective ℓ_min also adopted, which is motivated by the high-pass filtering process for ℓ ≤ 40 in component separation.

Within the above ℓ ranges, the theoretical bispectrum templates are decomposed with the eigenbasis composed of polynomials and some special functions reconstructing the tensor-mode features, e.g., temperature enhancement due to the ISW effect (ℓ ≲ 100), and two E-mode peaks created by reionization (ℓ ≲ 10) and recombination (ℓ ≃ 100). The resulting factorized templates are more than 95% correlated with the original ones. The validity of our numerical computations has been confirmed through the map-by-map comparisons of $\mathit{f̂}\begin{array}{c}\mathrm{even}\mathit{/}\mathrm{odd}\\ \mathrm{NL}\end{array}$ at very low resolution, showing the consistency between the values from the Modal methodology and those obtained by the brute-force summations like Eq. (36). We have also checked that our parity-even estimator successfully leads to the constraints on ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$, ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}$, and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}$ at ℓ_max = 500, compatible with the results from the binned estimator.

Our limits estimated from the foreground-cleaned temperature and high-pass filtered polarization data (SMICA, SEVEM, and NILC) are summarized in Table 26. The data and MC simulations used here, including all experimental features, i.e., beam, anisotropic noise levels and partial sky mask, have been inpainted using the identical recursive process adopted in the standard f_NL estimations (see Sect. 3.5). The sky fractions of the temperature and polarization masks adopted here are, respectively, f_sky = 0.76 and f_sky = 0.74. Although the error bars and the linear terms have been computed using 160 MC simulations, the resulting error bars are close to the expected values, (f_skyF)^{− 1 / 2}.

We have confirmed the stability of the T-only constraints, and significant scatter of the E-only constraints both in the parity-even case and in the parity-odd one, when changing f_sky. Such E-mode instability has given insignificant effects on our T+E constraints in the parity-even case, as they are determined almost exclusively by TTT, like the scalar NG analyses. In contrast, our parity-odd T+E results vary a lot, due to the E-mode scatter (quantitatively speaking, only a few percent change of f_sky has shifted ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}$ by more than 1σ), because TTE and TEE contribute significantly to the signal-to-noise ratio in the parity-odd case (Shiraishi et al. 2013b). Table 27 presents the correlations of the bispectra reconstructed from the component separated maps, also indicating the robustness of the T-only constraints and the instability of the E-only results. We report only stable results in Table 26 and conclude that there is no evidence at > 2σ of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}$ in the parity-even, parity-odd or their whole domain.

The parity-odd components of the TTT and EEE bispectra extracted model-independently from the SMICA data are visually represented in Fig. 26. It is apparent from this figure that the PlanckTTT bispectrum has similar features to the WMAP one (Shiraishi et al. 2015), e.g., distinctive signals distributed around ℓ₁ ≈ ℓ₂ ≈ ℓ₃. As indicated by the roughly 70% correlation between the SMICA and WMAP bispectra (see Table 27), the Planck T-only limits in Table 26 are close to the WMAP ones (68% CL): ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}\mathit{/}{\mathrm{10}}^{\mathrm{2}}\mathrm{=}\mathrm{4}\mathrm{\pm }\mathrm{16}$ for parity-even; and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}\mathit{/}{\mathrm{10}}^{\mathrm{2}}\mathrm{=}\mathrm{80}\mathrm{\pm }\mathrm{110}$ for parity-odd (Shiraishi et al. 2015).

9. Limits on the primordial trispectrum

So far, we have considered a variety of physically motivated possibilities for the inflationary 3-point function, or bispectrum. A similar phenomenology exists for the 4-point function, or trispectrum. Our constraints on the trispectrum will use CMB temperature only; we do not use polarization in this section. We start by briefly reviewing the inflationary physics and classifying the signals we will search for.

First, some notation: a “primed” ζ-trispectrum ⟨ ζ_k1ζ_k2ζ_k3ζ_k4 ⟩ ′ denotes the connected trispectrum without its momentum-conserving delta function, i.e., $\mathrm{⟨}{\mathit{\zeta }}_{{{k}}_{\mathrm{1}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{2}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{3}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{4}}}\mathrm{⟩}\mathrm{=}\mathrm{⟨}{\mathit{\zeta }}_{{{k}}_{\mathrm{1}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{2}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{3}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{4}}}{\mathrm{⟩}}^{\mathrm{\prime }}\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{3}}{\mathit{\delta }}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}(\begin{array}{c}\sum \end{array}{{k}}_{\mathit{i}})\mathrm{+}\mathrm{disc}\mathrm{.}\mathit{,}$(62)where “+ disc.” denotes the disconnected contributions to the 4-point function.

One possible signal is the “local” trispectrum ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$, which arises if the non-Gaussian adiabatic curvature ζ is of cubic-type form (see, e.g., Okamoto & Hu 2002), i.e., $\mathit{\zeta }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{\mathit{\zeta }}_{\mathrm{G}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{+}\frac{\mathrm{9}}{\mathrm{25}}\hspace{0.17em}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\hspace{0.17em}{\mathit{\zeta }}_{\mathrm{G}}\mathrm{\left(}{x}\mathrm{\right)}3$(63)where ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ is a free parameter and ζ_G is a Gaussian field. In this model, the bispectrum is zero (since there is a ζ → − ζ symmetry) and the trispectrum is given by $\begin{array}{ccc}\mathrm{⟨}{\mathit{\zeta }}_{{{k}}_{\mathrm{1}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{2}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{3}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{4}}}{\mathrm{⟩}}^{\mathrm{\prime }}\mathrm{=}\frac{\mathrm{54}}{\mathrm{25}}\hspace{0.17em}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{[}{\mathit{P}}_{\mathit{\zeta }}\mathrm{\left(}{\mathit{k}}_{\mathrm{1}}\mathrm{\right)}{\mathit{P}}_{\mathit{\zeta }}\mathrm{\left(}{\mathit{k}}_{\mathrm{2}}\mathrm{\right)}{\mathit{P}}_{\mathit{\zeta }}\mathrm{\left(}{\mathit{k}}_{\mathrm{3}}\mathrm{\right)}\mathrm{+}\mathrm{3}\mathrm{perm}\mathrm{.}\mathrm{]}\mathit{.}& & \end{array}$(64)Analogously to the case of the local bispectrum, ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$, the observational signal-to-noise for ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ is largest in the “squeezed” limit, k₁ ≪ min(k₂,k₃,k₄), and there is a consistency relation which shows that in single-field inflation, the four-point function is always small in the squeezed limit (e.g., Senatore & Zaldarriaga 2012a). Thus ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ can only be detectably large in multi-field models. Conversely, there are multi-field models where ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ is detectable. The main obstacle here is technical naturalness, i.e., ensuring that radiative corrections do not generate an observationally larger bispectrum. This can be the case if a large bispectrum is forbidden by a Z₂ symmetry, or by supersymmetry (Senatore & Zaldarriaga 2012b).

A further category of four-point signals can be obtained by adding quartic interactions to the inflationary action. Following Smith et al. (2015), we will concentrate on the simplest possibility, by considering quartic operators consistent with the symmetries of inflation and having the lowest possible number of derivatives (Bartolo et al. 2010b; Senatore & Zaldarriaga 2011, 2012b). There are three such operators, of the schematic form , , and (∂_iσ)²(∂_jσ)². By a short calculation using the in-in formalism (Maldacena 2003), the associated four-point functions are (Smith et al. 2015; see also Huang & Shiu 2006): $\begin{array}{ccc}\mathrm{⟨}{\mathit{\zeta }}_{{{k}}_{\mathrm{1}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{2}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{3}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{4}}}{\mathrm{⟩}}^{\mathrm{\prime }}& \mathrm{=}& \frac{\mathrm{9216}}{\mathrm{25}}{\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}\mathrm{4}}{\mathit{A}}_{\mathit{\zeta }}^{\mathrm{3}}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{0}}\mathrm{d}{\mathit{\tau }}_{\mathit{E}}\hspace{0.17em}{\mathit{\tau }}_{\mathit{E}}^{\mathrm{4}}\left(\underset{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}{\overset{\mathrm{4}}{\mathrm{􏽙}}}\frac{{\mathit{e}}^{{\mathit{k}}_{\mathit{i}}{\mathit{\tau }}_{\mathit{E}}}}{{\mathit{k}}_{\mathit{i}}}\right)\\ & \mathrm{=}& \\ \mathrm{⟨}{\mathit{\zeta }}_{{{k}}_{\mathrm{1}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{2}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{3}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{4}}}{\mathrm{⟩}}^{\mathrm{\prime }}& \mathrm{=}& \mathrm{-}\frac{\mathrm{13}\hspace{0.17em}\mathrm{824}}{\mathrm{325}}{\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}\mathrm{2}\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{2}}}{\mathit{A}}_{\mathit{\zeta }}^{\mathrm{3}}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{0}}\mathrm{d}{\mathit{\tau }}_{\mathit{E}}\hspace{0.17em}{\mathit{\tau }}_{\mathit{E}}^{\mathrm{2}}\\ & & \mathrm{\times }\hspace{0.17em}[\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{k}}_{\mathrm{3}}{\mathit{\tau }}_{\mathit{E}}\mathrm{)}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{k}}_{\mathrm{4}}{\mathit{\tau }}_{\mathit{E}}\mathrm{)}}{{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}^{\mathrm{3}}{\mathit{k}}_{\mathrm{4}}^{\mathrm{3}}}\mathrm{(}{{k}}_{\mathrm{3}}\mathrm{·}{{k}}_{\mathrm{4}}\mathrm{)}{\mathit{e}}^{\sum {\mathit{k}}_{\mathit{i}}{\mathit{\tau }}_{\mathit{E}}}\\ & & \mathrm{+}\mathrm{5}\mathrm{perm}\mathrm{.}]\\ & \mathrm{=}& \mathrm{-}\frac{\mathrm{27}\hspace{0.17em}\mathrm{648}}{\mathrm{325}}{\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}\mathrm{2}\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{2}}}{\mathit{A}}_{\mathit{\zeta }}^{\mathrm{3}}\\ & & \mathrm{\times }\hspace{0.17em}[\frac{{\mathit{K}}^{\mathrm{2}}\mathrm{+}\mathrm{3}\mathrm{(}{\mathit{k}}_{\mathrm{3}}\mathrm{+}{\mathit{k}}_{\mathrm{4}}\mathrm{)}\mathit{K}\mathrm{+}\mathrm{12}{\mathit{k}}_{\mathrm{3}}{\mathit{k}}_{\mathrm{4}}}{{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}^{\mathrm{3}}{\mathit{k}}_{\mathrm{4}}^{\mathrm{3}}{\mathit{K}}^{\mathrm{5}}}\mathrm{(}{{k}}_{\mathrm{3}}\mathrm{·}{{k}}_{\mathrm{4}}\mathrm{)}\\ & & \\ \mathrm{⟨}{\mathit{\zeta }}_{{{k}}_{\mathrm{1}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{2}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{3}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{4}}}{\mathrm{⟩}}^{\mathrm{\prime }}& \mathrm{=}& \frac{\mathrm{82}\hspace{0.17em}\mathrm{944}}{\mathrm{2575}}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}{\mathit{A}}_{\mathit{\zeta }}^{\mathrm{3}}{\mathrm{\int }}_{\mathrm{-}\mathrm{\infty }}^{\mathrm{0}}\mathrm{d}{\mathit{\tau }}_{\mathit{E}}\\ & & \mathrm{\times }\hspace{0.17em}\left[\underset{\mathit{i}\hspace{0.17em}\mathrm{=}\hspace{0.17em}\mathrm{1}}{\overset{\mathrm{4}}{\mathrm{􏽙}}}\frac{\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{k}}_{\mathit{i}}{\mathit{\tau }}_{\mathit{E}}\mathrm{)}{\mathit{e}}^{{\mathit{k}}_{\mathit{i}}{\mathit{\tau }}_{\mathit{E}}}}{{\mathit{k}}_{\mathit{i}}^{\mathrm{3}}}\right]\\ & & \mathrm{\times }\hspace{0.17em}\left[\mathrm{(}{{k}}_{\mathrm{1}}\mathrm{·}{{k}}_{\mathrm{2}}\mathrm{)}\mathrm{(}{{k}}_{\mathrm{3}}\mathrm{·}{{k}}_{\mathrm{4}}\mathrm{)}\mathrm{+}\mathrm{2}\mathrm{perm}\mathrm{.}\right]\\ & \mathrm{=}& \frac{\mathrm{165}\hspace{0.17em}\mathrm{888}}{\mathrm{2575}}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}{\mathit{A}}_{\mathit{\zeta }}^{\mathrm{3}}\\ & & \mathrm{\times }\hspace{0.17em}\left(\frac{\mathrm{2}{\mathit{K}}^{\mathrm{4}}\mathrm{-}\mathrm{2}{\mathit{K}}^{\mathrm{2}}\sum {\mathit{k}}_{\mathit{i}}^{\mathrm{2}}\mathrm{+}\mathit{K}\sum {\mathit{k}}_{\mathit{i}}^{\mathrm{3}}\mathrm{+}\mathrm{12}{\mathit{k}}_{\mathrm{1}}{\mathit{k}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{3}}{\mathit{k}}_{\mathrm{4}}}{{\mathit{k}}_{\mathrm{1}}^{\mathrm{3}}{\mathit{k}}_{\mathrm{2}}^{\mathrm{3}}{\mathit{k}}_{\mathrm{3}}^{\mathrm{3}}{\mathit{k}}_{\mathrm{4}}^{\mathrm{3}}{\mathit{K}}^{\mathrm{5}}}\right)\\ & & \end{array}$Here K = ∑ k_i and we have introduced parameters ${\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}{\mathrm{4}}^{}}$, ${\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}{\mathrm{2}}^{}\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{2}}}$, and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}$ to parametrize the amplitude of each trispectrum.

The normalization of the g_NL-parameters is chosen so that $\mathrm{⟨}{\mathit{\zeta }}_{{{k}}_{\mathrm{1}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{2}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{3}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{4}}}\mathrm{⟩}\mathrm{=}\mathrm{(}\mathrm{216}\mathit{/}\mathrm{25}\mathrm{)}\hspace{0.17em}{\mathit{g}}_{\mathrm{NL}}{\mathit{A}}_{\mathit{\zeta }}^{\mathrm{3}}\mathit{/}{\mathit{k}}^{\mathrm{9}}$ for “tetrahedral” configurations, with | k_i | = k and (k_i·k_j) = − k³/ 3. This is the analogue of the commonly-used normalization for the bispectrum, where f_NL parameters are defined so that $\mathrm{⟨}{\mathit{\zeta }}_{{{k}}_{\mathrm{1}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{2}}}{\mathit{\zeta }}_{{{k}}_{\mathrm{3}}}\mathrm{⟩}\mathrm{=}\mathrm{(}\mathrm{18}\mathit{/}\mathrm{5}\mathrm{)}\hspace{0.17em}{\mathit{f}}_{\mathrm{NL}}{\mathit{A}}_{\mathit{\zeta }}^{\mathrm{2}}\mathit{/}{\mathit{k}}^{\mathrm{6}}$ for equilateral configurations with | k_i | = k.

For simplicity in Eqs. (65)–(67) we have assumed a scale-invariant initial power spectrum P_ζ(k) = A_ζ/k³. In order to analyse Planck data, we must slightly generalize this to a power-law spectrum P_ζ(k) ∝ k^{n_s − 4}. Our scheme for doing this follows Appendix C of Smith et al. (2015).

A Fisher matrix analysis shows that there is one large correlation among the three trispectra in Eqs. (65)–(67), so that to an excellent approximation we can treat only two of the trispectra as independent. To quantify this, in Smith et al. (2015) it is shown that the shape is 98.6 % correlated to a linear combination of the shapes and (∂σ)⁴. Therefore, we will only search for the parameters ${\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}{\mathrm{4}}^{}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}$.

We note that the analysis that leads to the trispectrum shapes ${\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}{\mathrm{4}}^{}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}$ is very similar to the analysis that leads to the “standard” bispectrum shapes ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equi}}$ and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}$. However, there are some minor differences as follows. In the bispectrum case, one considers the cubic operators and , but it is conventional to define observables ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equi}}\mathit{,}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}$ which are related to the operator coefficients by a linear transformation. This is done because the two cubic operators are about 90 % correlated, so it is convenient to orthogonalize. In the trispectrum case the correlation is smaller (around 60 % for Planck), and we have chosen to omit the orthogonalization step. Another reason to omit the orthogonalization step is that the trispectrum shape (∂σ)⁴ is a signature of multi-field inflation. In single field inflation, the (∂σ)⁴ trispectrum is not technically natural; radiative corrections generate cubic operators of the form or , which generate a bispectrum with larger signal-to-noise.

There are more trispectrum shapes one might consider. For example, classifying Galilean invariant quartic operators leads to higher-derivative trispectra, which are not highly correlated to the trispectra considered above (Bartolo et al. 2013c; Arroja et al. 2013). We have only considered “contact” diagrams arising from one power of a quartic operator, and it would be interesting to study “exchange” diagrams arising from two cubic operators and exchange of a light particle (e.g., Chen et al. 2009; Arroja et al. 2009; Chen & Wang 2010; Bartolo et al. 2010a; Baumann & Green 2012). We leave these as extensions for future work.

Summarizing, we will search for the following trispectrum signals: $\left\{{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}\mathrm{4}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}\right\}$(68)defined by Eqs. (64), (65), and (67) above.

9.1. Data analysis

Turning now to data analysis, we use the machinery from Smith et al. (2015). The first step is to represent each trispectrum as a small sum of factorizable terms as follows. The angular CMB trispectrum can be written either as an integral over comoving distance r (in the case of ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$) or as a double integral over (τ,r) where τ is conformal time during inflation (in the case of the or (∂σ)⁴ trispectra). We approximate the integral by a finite sum, which represents the CMB trispectrum as a sum of terms that are factorizable in a sense defined in Smith et al. (2015). A large number of sampling points are needed to obtain a good approximation to the integral, leading to a large number of terms in the factorizable representation. However, there exists an optimization algorithm, which takes as input a trispectrum that has been represented as a sum of many factorizable terms, and outputs a representation with fewer terms. The reduction can be quite dramatic, as shown in Table 28. The optimization algorithm guarantees that the output trispectrum accurately approximates the input trispectrum, in the sense that the two are nearly observationally indistinguishable.

Armed with “small” factorizable representations for each trispectrum, the next step is to run an analysis pipeline that estimates the amplitude of each trispectrum from Planck data. We use the “pure MC” pipeline from Smith et al. (2015), which compares the trispectrum of the data to the mean trispectrum of an ensemble of simulations. This pipeline requires a filtering operation $\mathit{d}\mathrm{\to }\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{a}}_{\mathit{\ell m}}\end{array}$ which processes the pixel-space CMB data d and generates a harmonic-space map $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{a}}_{\mathit{\ell m}}\end{array}$. Our filtering operation is defined by the following steps:

• 1.

  Starting from the data d, we compute (with uniform pixel weighting) a best-fit monopole and dipole outside the Galactic mask. We use the temperature “common mask”, the union of the confidence masks for the SMICA, SEVEM, NILC, and Commander component separation methods (Planck Collaboration IX 2016).

• 2.

  The mask defines a few “islands”, i.e., isolated groups of pixels that are unmasked, but contained in a larger masked region. We slightly enlarge the mask so that it removes the islands.

• 3.

  We classify the components of the masked part of the sky into “small” masked regions with ≤ 1000 pixels (at HEALPix resolution N_side = 2048), and “large” regions with > 1000 pixels. Small regions usually correspond to point sources, and large regions typically correspond to areas of diffuse galactic emission. In small regions, we inpaint the CMB by assigning the unique map that agrees with the data on boundary pixels, and whose value in each interior pixel is the average of the neighboring pixels.

• 4.

  In large regions, we do not inpaint the CMB, but rather apodize the boundary of the large region using cosine apodization with 12′ radius.

• 5.

  We apply a spherical harmonic transform to the inpainted, apodized CMB map to obtain a harmonic-space map a_ℓm with ℓ_max = 1600. We then take the final filtered map $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{a}}_{\mathit{\ell m}}\end{array}$ to be $\begin{array}{c}\mathrm{˜}\\ {\mathit{a}}_{\mathit{\ell m}}\end{array}\mathrm{=}\frac{{\mathit{a}}_{\mathit{\ell m}}}{{\mathit{b}}_{\mathit{\ell }}{\mathit{C}}_{\mathit{\ell }}\mathrm{+}{\mathit{b}}_{\mathit{\ell }}^{-1}{\mathit{N}}_{\mathit{\ell }}}$(69)where b_ℓ is the beam, C_ℓ is the fiducial CMB power spectrum, and N_ℓ is the sky-averaged noise power spectrum (without beam deconvolution). To motivate this choice of ℓ-weighting, we note that for an ideal all-sky experiment with isotropic noise, we have a_ℓm = b_ℓs_ℓm + n_ℓm where s_ℓm,n_ℓm are signal and noise realizations. In this case, Eq. (69) weights the signal as ${\mathit{s}}_{\mathit{\ell m}}\mathit{/}\mathrm{(}{\mathit{C}}_{\mathit{\ell }}\mathrm{+}{\mathit{b}}_{\mathit{\ell }}^{-2}{\mathit{N}}_{\mathit{\ell }}\mathrm{)}$, which is optimal.

In our pipeline, we apply this filter to the component-separated SMICA maps (Planck Collaboration IX 2016), obtaining a harmonic-space map $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{a}}_{\mathit{\ell m}}\end{array}$. We apply the same filter to 1000 Monte Carlo simulations to obtain an ensemble of harmonic-space maps. Our pipeline has the property that it always estimates the trispectrum of the data in excess of the trispectrum in the simulations. Since the simulations include lensing, this means that lensing bias will automatically be subtracted from our g_NL estimates.

Now that the filter, data realization, and Monte Carlo simulations have been fully specified, the details of the pipeline are described in Sect. IX.B of Smith et al. (2015). For each trispectrum, the pipeline outputs an estimate of g_NL and an estimate of the statistical error. Our basic results are: $\begin{array}{ccc}& & {\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{9.0}\mathrm{\pm }\mathrm{7.7}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{\mathrm{4}}\mathrm{;}\\ & & \\ & & {\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{0.1}\mathrm{\pm }\mathrm{3.8}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}\mathit{.}\end{array}$(70)No deviation from Gaussian statistics is seen. These results significantly improve the previous best constraints on the trispectrum from WMAP (Vielva & Sanz 2010; Smidt et al. 2010; Fergusson et al. 2010b; Hikage & Matsubara 2012; Sekiguchi & Sugiyama 2013; Regan et al. 2013; Smith et al. 2015) and large-scale structure (Desjacques & Seljak 2010; Giannantonio et al. 2014; Leistedt et al. 2014).

A constraint on ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ from Planck 2013 data was recently reported by Feng et al. (2015), who find ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{13}\hspace{0.17em}\mathrm{\pm }\hspace{0.17em}\mathrm{18}\mathrm{)}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}{\mathrm{10}}^{\mathrm{4}}$. Our central value in Eq. (70) agrees well with this result, but the statistical error is smaller by a factor of 2.3. This improvement is partly due to the lower noise levels in Planck 2015 data, and partly due to the use of a better estimator.

Each line in Eq. (70) is a “single-g_NL” constraint; i.e., the constraint on one g_NL parameter with the other g_NL-parameters held fixed. For joint constraints, one needs to know the full covariance matrix. The correlation between ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ and the other two parameters is negligble, and the ${\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}{\mathrm{4}}^{}}$-${\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}$ correlation is: $\mathrm{Corr}\left({\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}\mathrm{4}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}\right)\mathrm{=}\mathrm{0.61.}$(71)Multi-field models of inflation will generally give a linear combination of , , and (∂_iσ)²(∂_jσ)² trispectra. In this case we proceed as follows. First, if the coefficient is non-zero, we can use the near-degeneracy with a linear combination of the other two operators to absorb it into the effective values of ${\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}{\mathrm{4}}^{}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}$. A Fisher matrix analysis shows that the coefficients of this linear combination are $\begin{array}{ccc}& & \mathrm{(}{\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}\mathrm{4}}{\mathrm{)}}_{\mathrm{eff}}\mathrm{=}\mathrm{0.59}\hspace{0.17em}{\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}\mathrm{2}\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{2}}}\\ & & \end{array}$(72)It is convenient to define the two-component parameter vector: ${\mathit{g}}_{\mathit{i}}\mathrm{=}\left(\begin{array}{c}\\ {\mathit{g}}_{\mathrm{NL}}^{\mathit{\sigma ̇}\mathrm{4}}\\ {\mathit{g}}_{\mathrm{NL}}^{\mathrm{(}\mathit{\partial \sigma }{\mathrm{)}}^{\mathrm{4}}}\end{array}\right)\mathit{.}$(73)We also compute a two-by-two Fisher matrix F_ij, whose diagonal is given by ${\mathit{F}}_{\mathit{ii}}\mathrm{=}\mathrm{1}\mathit{/}{\mathit{\sigma }}_{\mathit{i}}^{\mathrm{2}}$, where σ_i is the single-g_NL statistical error in Eq. (70), and whose off-diagonal is ${\mathit{F}}_{\mathrm{12}}\mathrm{=}\mathit{r}{\mathit{F}}_{\mathrm{11}}^{\mathrm{1}\mathit{/}\mathrm{2}}{\mathit{F}}_{\mathrm{22}}^{\mathrm{1}\mathit{/}\mathrm{2}}$, where r is the correlation in Eq. (71). This procedure gives: ${\mathit{F}}_{\mathit{ij}}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{3.3}& \mathrm{9.2}\\ \mathrm{9.2}& \mathrm{68.7}\end{array}\right)\mathrm{\times }{\mathrm{10}}^{-13}\mathit{.}$(74)For a given parameter vector g_i, we can define a trispectrum-χ² by ${\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{g}\mathrm{\right)}\mathrm{=}\mathrm{[}{\mathit{F}}_{\mathit{ii}}\mathit{ĝ}\mathit{i}\mathrm{-}\mathrm{(}\mathit{Fg}{\mathrm{)}}_{\mathit{i}}\mathrm{]}\hspace{0.17em}{\mathit{F}}_{\mathit{ij}}^{-1}\hspace{0.17em}\mathrm{[}{\mathit{F}}_{\mathit{jj}}\mathit{ĝ}\mathit{j}\mathrm{-}\mathrm{\left(}\mathit{Fg}{\mathrm{)}}_{\mathit{j}}\mathrm{\right]}$(75)where ĝ_i = ( − 0.21 × 10⁶, − 0.10 × 10⁵) is the vector of best-fit single-g_NL values from Eq. (70). This definition of χ² follows from the observation that (F_iiĝ_i) is an estimator with expectation value (Fg)_i and covariance matrix Cov(F_iiĝ_i,F_jjĝ_j) = F_ij.

The inflationary implications of these trispectrum constraints are discussed in Sect. 11.5 below.

10. Minkowski functionals results

In this section, we present constraints on local NG at first and second order (${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$) obtained with Minkowski functionals (MFs) on temperature and polarization E maps. MFs (Mecke et al. 1994; Schmalzing & Buchert 1997; Schmalzing & Gorski 1998; Winitzki & Kosowsky 1998) are a measure of fields’ local morphology used to constrain their stationarity, isotropy and Gaussianity. Mostly probing general NG in a frequentist fashion in two-dimensions on CMB maps (Eriksen et al. 2004; Komatsu et al. 2005; Modest et al. 2013; Natoli et al. 2010; Curto et al. 2008) or three-dimensions on LSS data (Park et al. 2005; Wiegand et al. 2014), they have also been used to measure specific NG targets with Bayesian methods, such as ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ (Hikage et al. 2006, 2008; Ducout et al. 2013; Planck Collaboration XXIV 2014), other bispectrum and trispectrum shapes (Hikage & Matsubara 2012) and topological defects (Planck Collaboration XXV 2014). New developments have been made recently, using needlets (Fantaye et al. 2015), neural networks (Novaes et al. 2015) or allowing scale-dependent measurements (Munshi et al. 2013).

MF-based limits are well known to be suboptimal for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$, but they provide an independent cross-check of bispectrum and trispectrum-based estimators. They are complementary to optimal estimators: they are weighted integrals of the polyspectra and are sensitive to any source of NG, including foregrounds and secondaries. While MFs are not always able to distinguish between these different sources and systematics, they allow upper bounds to be put on them.

The most recent constraints (1σ) from MFs on ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ have been obtained with WMAP (Hikage & Matsubara 2012) and Planck (Planck Collaboration XXIV 2014): ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{4.2}\mathrm{\pm }\mathrm{20.5}\mathrm{;} {\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{(}\mathrm{1.9}\mathrm{\pm }\mathrm{6.4}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{\mathrm{5}}\mathit{.}$(76)

10.1. Method and definition of MFs

For a smoothed two-dimensional field δ of zero mean and of variance ${\mathit{\sigma }}_{\mathrm{0}}^{\mathrm{2}}$, defined on the sphere, we consider an excursion set of height ν = δ/σ₀, i.e., the set of points where the field exceeds the threshold ν. We use four functionals denoted by V_k(ν)(k = 0,1,2,3). The first three correspond to MFs: V₀ is the fractional Area of the regions above the threshold, V₁ is the Perimeter of these regions and V₂ is the Genus, defined as the total number of connected components of the excursion above the threshold minus the total number of connected components under the threshold. The fourth, V₃, is the Number of clusters (also referred to as N_clusters). This is the number of connected regions above the threshold for positive thresholds and below the threshold for negative thresholds. Precise definitions and formulae for the quantities V_k as well as their expectation values for Gaussian fields are summarized in Appendix B.

We calculate the four normalized¹⁸ functionals v_k(ν) on n_th = 26 thresholds ν, between ν_min = − 3.5 and ν_max = + 3.5.

For this analysis, we used the same temperature and polarization E data, simulations and masks described in Sect. 3.4 for consistency with the bispectrum estimators. In addition, the maps are filtered to optimize constraints on local NG (Ducout et al. 2013), the filters used being similar to Wiener filters for T and E (W_M), and for the first (W_D1) and second (W_D2) derivatives of these fields (Fig. 27): $\begin{array}{ccc}{\mathit{W}}_{\mathrm{D}\mathrm{1}}& \mathrm{\propto }& \sqrt{\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}\hspace{0.17em}{\mathit{W}}_{\mathrm{M}}\hspace{0.17em}\mathrm{;}\\ {\mathit{W}}_{\mathrm{D}\mathrm{2}}& \mathrm{\propto }& \end{array}$For the temperature map, known point sources in the mask are inpainted.

We define the vector y as any combination $\mathit{y}\mathrm{=}\mathrm{\left\{}{\mathit{v}}_{\mathit{k}}^{\mathit{A,W}}\mathrm{\right\}}$ with k = { 0,3 }, A = { T,E }, W = { W_M,W_D1,W_D2 }, ŷ being the vector measured on the data.

From these measurements, we then use a Bayesian method to jointly estimate ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$, $\begin{array}{ccc}\mathit{P}\mathrm{\left(}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}\hspace{0.17em}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\right|}\hspace{0.17em}\mathit{ŷ}\mathrm{)}\mathrm{=}\frac{\mathit{P}\mathrm{\left(}\mathit{ŷ}\hspace{0.17em}\mathrm{\right|}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\left)}\mathit{P}\mathrm{\right(}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{)}}{\mathrm{\int }\mathit{P}\mathrm{\left(}\mathit{ŷ}\hspace{0.17em}\mathrm{\right|}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\left)}\mathit{P}\mathrm{\right(}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{)}\mathrm{d}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{d}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}}\mathrm{·}& & \end{array}$(79)We take a uniform prior for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ in the range − 400 to 400, and for ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ in the range − 4 × 10⁶ to + 4 × 10⁶, while the evidence is just considered as a normalization.

Assuming MFs are multi-variate Gaussian-distributed we obtain the posterior distribution for $\mathrm{\left(}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\right)}$ with a χ² test $\mathit{P}\mathrm{(}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}\hspace{0.17em}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{|}\hspace{0.17em}\mathit{ŷ}\mathrm{)}\mathrm{\propto }\exp \left[\mathrm{-}\frac{{\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{ŷ}\mathit{,}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\right)}}{\mathrm{2}}\right]\mathit{,}$(80)with $\begin{array}{ccc}& & {\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{ŷ}\mathit{,}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\right)}\mathrm{\equiv }\\ & & \end{array}$(81)For this test, we use two types of simulations to first construct a model including primordial NG $\mathit{y̅}{}_{\mathrm{sim}\mathrm{1}}\mathrm{\left(}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\right)}$ and secondly a covariance matrix ${\mathit{C}}_{\mathrm{sim}\mathrm{2}}\mathrm{\equiv }\mathrm{⟨}\left({\mathit{y}}_{\mathrm{sim}\mathrm{2}}\mathrm{-}\mathit{y̅}\mathrm{sim}\mathrm{2}\right){\left({\mathit{y}}_{\mathrm{sim}\mathrm{2}}\mathrm{-}\mathit{y̅}\mathrm{sim}\mathrm{2}\right)}^{\mathrm{T}}\mathrm{⟩}\mathit{,}$(82)with averaged over the simulations. We now describe the details of these simulations.

• Simulations 1: Non-Gaussian model For the firsttype of simulation, we included all possible sources of NG,assuming that the total and individual levels of NG are smallenough that MFs are linear with respect to those NG levels(Ducout et al. 2013). The threekinds of NG we included are foreground residuals (Galacticresiduals with scalable amplitude α, as well as radio sources, CIB anisotropies, secondaries (SZ, lensing and ISW-lensing, but not SZ-lensing) and primordial NG (${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$, ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$): $\begin{array}{ccc}{\mathrm{sim}\mathrm{1}}^{\mathit{i}}& \mathrm{=}& {\mathrm{map}}_{\mathrm{lensed}}^{\mathit{i}}\mathrm{\left(}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathit{,}{\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\right)}\\ & & \mathrm{+}\mathrm{ma}{\mathrm{p}}^{\mathrm{fg}}\mathrm{\left(}\mathrm{radio}\hspace{0.17em}\mathrm{sources}\mathit{,}\mathrm{CIB}\mathit{,}\mathrm{SZ}\mathrm{\right)}\\ & & \mathrm{+}\mathit{\alpha }\mathrm{\times }\mathrm{ma}{\mathrm{p}}^{\mathrm{fg}}\mathrm{\left(}\mathrm{Galactic}\hspace{0.17em}\mathrm{residuals}\mathrm{\right)}\mathit{.}\end{array}$(83)We tried to reproduce all instrumental effects, with realistic effective beams (isotropic window functions), noise from FFP8 simulations (Planck Collaboration XII 2016), filtered with component separation weights. We checked the accuracy of these simulations by comparing them to FFP8 simulations, using no foreground and no primordial NG. The astrophysical models are provided by the Planck Sky Model (PSM, Delabrouille et al. 2013), while the primordial NG simulations are computed as in Elsner & Wandelt (2009). The lensing uses LensPix¹⁹. The power spectrum used for these NG simulations and the lensing is the best-fit power spectrum from Planck 2013+ACT/SPT+BAO (Planck Collaboration XXII 2014). We created n₁ = 200 simulations i, using n₁ maps for the primordial NG, while we had only one astrophysical foreground simulation.

• Simulations 2: FFP8 (Planck Collaboration XII 2016) MC simulations Since NG is weak, the covariance matrix C is computed with n₂ = 10⁴ simulations, including no primordial NG and no foregrounds. These simulations reproduce realistic instrumental effects (anisotropy of beams in particular), realistic noise and component separation filtering. The only NG still present in these simulations are lensing and the ISW-lensing correlation.

Validation of the estimator

Part of the validation for the MFs estimator is described in Sect. 5.3 for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$, to compare the results to bispectrum estimators on realistic simulations (FFP8 MC, second item above). In addition we present in Table 29 the results obtained on the same realistic simulations for ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$, and on simulations containing primordial NG (first item above), with ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{10}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}{\mathrm{10}}^{\mathrm{5}}$, with 200 simulations used in each case. These tests have been performed using the SMICA method with lensing bias removed.

10.2. Results

Results for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ estimation with MFs on the four component separated maps in temperature and polarization are presented in Table 30. The results for polarization E-only maps are not quoted, since these results were not sufficiently stable (cf. Sect. 7.6). No deviation from Gaussianity is detected. T+E analysis generally finds higher values for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$, but remains consistent with Gaussianity.

The posteriors for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ from SMICA are shown in Fig. 28. One interesting point is that the estimates of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ and ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ are almost uncorrelated (r< 0.1); this can be inferred when we consider the parity of MF deviations from Gaussianity, which is different for the two parameters (Matsubara 2010).

Foreground and secondary biases are removed from these estimates, since the NG model directly includes them. However, an estimation of their contribution in the map is reported in Table 31.

Foreground and secondary biases

Foreground residuals are generally negligible, in particular in the T analysis. This is different from the Planck 2013 results, where the residuals were more important; this can be explained by the beam correction applied to these previous estimates which exaggerated signals from small scales.

Lensing has a significant signature in MF estimation of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$, but is even stronger in ${\mathit{g}}_{\mathrm{NL}}^{\mathrm{local}}$ (the four-point correlation signature) and could be detected (and not treated just as a bias) with this estimator. The Wiener filters enhance the scales where lensing is dominant.

11. Implications for early Universe physics

The NG constraints obtained in this paper show consistency of Planck data with Gaussian primordial fluctuations, thus confirming the results obtained in the 2013 release (Planck Collaboration XXIV 2014) and improving them through the inclusion of CMB polarization data. The standard single-field slow-roll models of inflation have therefore been confirmed as a viable scenario for inflation, passing one of their most stringent tests, based on lack of measurable deviations from Gaussianity. The constraints obtained on local, equilateral, and orthogonal NG, after accounting for various contaminants, strongly limit different mechanisms proposed as alternatives to the standard inflationary models to explain the seeds of cosmological perturbations. Measurements on deviations from Gaussianity for other primordial bispectral shapes help to shed light on more subtle effects about the detailed physics of inflation.

As in Planck Collaboration XXIV (2014), in the following we derive limits on parameters of the models from the NG constraints in the following way (unless explicitly stated otherwise): we construct a posterior based on the assumption that the sampling distribution is Gaussian (as supported by Gaussian simulations); the likelihood is approximated by the sampling distribution, but centred on the NG estimate (see Elsner & Wandelt 2009); we employ uniform or Jeffreys’ priors, over intervals of the parameters values that are physically meaningful, or as otherwise stated; and in the cases when two or more parameters are involved, we marginalize the posterior to provide one-dimensional constraints on the parameter considered.

11.1. General single-field models of inflation

DBI models: DBI models of inflation (Silverstein & Tong 2004; Alishahiha et al. 2004), characterized by a non-standard kinetic term of the inflaton field, predict a non-linearity parameter ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{DBI}}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{35}\mathit{/}\mathrm{108}\mathrm{)}\mathrm{(}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{-}\mathrm{1}\mathrm{)}$, where c_s is the sound speed of the inflaton perturbations (Silverstein & Tong 2004; Alishahiha et al. 2004; Chen et al. 2007b). The corresponding bispectrum shape is very close to the equilateral shape. Nonetheless we have constrained the exact theoretical (non-separable) shape (see Eq. (7) of Planck Collaboration XXIV 2014). The constraint we obtain ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{DBI}}\mathrm{=}\mathrm{2.6}\mathrm{\pm }\mathrm{61.6}$ from temperature data (${\mathit{f}}_{\mathrm{NL}}^{\mathrm{DBI}}\mathrm{=}\mathrm{15.6}\mathrm{\pm }\mathrm{37.3}$ from temperature and polarization) at 68% CL (with ISW-lensing and point sources subtracted, see Table 23) implies ${\mathit{c}}_{\mathrm{s}}^{\mathrm{DBI}}\mathrm{\ge }\mathrm{0.069} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathrm{-}\mathrm{only}\mathrm{)}\mathit{,}$(84)and ${\mathit{c}}_{\mathrm{s}}^{\mathrm{DBI}}\mathrm{\ge }\mathrm{0.087} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathit{+}\mathit{E}\mathrm{)}\mathit{.}$(85)In Planck Collaboration XXIV (2014) we constrained the so-called infrared (IR) DBI models (Chen 2005b,a), which arise in string frameworks. We focused on a minimal setup, considering a regime where stringy effects are negligible and predictions for primordial perturbations are built within standard field theory. In the companion paper Planck Collaboration XX (2016) we present an analysis of a more general class of IR DBI models which accounts for stringy signatures (see Bean et al. 2008) by combining Planck power spectrum and bispectrum constraints.

Implications for effective field theory of inflation: in this subsection we use the effective field theory approach to inflation in order to translate the contraints on ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}$ and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}$ into limits on the parameters of the Lagrangian of general single-field models of inflation (of the type P(X,ϕ) models). In particular we derive the most conservative bound on the sound speed of the inflaton perturbations for this class of models.

The effective field theory approach (Cheung et al. 2008; Weinberg 2008) provides an efficient way to constrain inflationary perturbations for various classes of models that incorporate deviations from the standard single-field slow-roll scenario. In this approach the Lagrangian of the system is expanded into the (lowest dimension) operators obeying the underlying symmetries. We consider general single-field models described by the following action $\begin{array}{ccc}\mathit{S}& \mathrm{=}& \int {\mathrm{d}}^{\mathrm{4}}\mathit{x}\sqrt{\mathrm{-}\mathit{g}}[\mathrm{-}\frac{{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}\mathit{Ḣ}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\left(\mathit{\pi ̇}\mathrm{2}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\frac{\mathrm{(}{\mathit{\partial }}_{\mathit{i}}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}}\right)\\ & & \mathrm{-}{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}\mathit{Ḣ}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{)}\mathit{\pi ̇}\frac{\mathrm{(}{\mathit{\partial }}_{\mathit{i}}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}{{\mathit{a}}^{\mathrm{2}}}\mathrm{+}\left({\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}\mathit{Ḣ}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{)}\mathrm{-}\frac{\mathrm{4}}{\mathrm{3}}{\mathit{M}}_{\mathrm{3}}^{\mathrm{4}}\right)\mathit{\pi ̇}\mathrm{3}]\mathit{,}\end{array}$(86)where the curvature perturbation is related to the scalar field π as ζ = − Hπ. The inflaton interaction terms and generate two kind of bispectra with amplitudes, respectively, ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{EFT}\mathrm{1}}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{85}\mathit{/}\mathrm{324}\mathrm{)}\mathrm{(}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{-}\mathrm{1}\mathrm{)}$ and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{EFT}\mathrm{2}}\mathrm{=}\mathrm{-}\mathrm{(}\mathrm{10}\mathit{/}\mathrm{243}\mathrm{)}\mathrm{(}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{-}\mathrm{1}\mathrm{)}{}^{\mathrm{[}}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathrm{+}\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}^{\mathrm{]}}$, where M₃ is the amplitude of the operator (Senatore et al. 2010; see also Chen et al. 2007b; Chen 2010b). These two bispectra both peak for equilateral triangles in Fourier space. Nevertheless, they are sufficiently different, and the total NG signal turns out to be a linear combination of the two, leading also to an orthogonal shape. We can put constraints on c_s and the dimensionless parameter ${\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{c}\end{array}}_{\mathrm{3}}\mathrm{(}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{=}\mathrm{2}{\mathit{M}}_{\mathrm{3}}^{\mathrm{4}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}\mathit{Ḣ}{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}\mathrm{\right)}$ (Senatore et al. 2010). Notice that DBI inflationary models corresponds to having $\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathrm{=}\mathrm{3}\mathrm{(}\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\mathit{/}\mathrm{2}$, while c_s = 1 and M₃ = 0 (or ${\begin{array}{c}{}_{\mathrm{˜}}\\ \mathit{c}\end{array}}_{\mathrm{3}}\mathrm{(}{\mathit{c}}_{\mathrm{s}}^{-2}\mathrm{-}\mathrm{1}\mathrm{)}\mathrm{=}\mathrm{0}\mathrm{)}$ represent the non-interacting (vanishing NG) case.

The mean values of the estimators for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}$ and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}$ are expressed in terms of c_s and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{c}}_{\mathrm{3}}\end{array}$ by $\begin{array}{ccc}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}& \mathrm{=}& \frac{\mathrm{1}\mathrm{-}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{[}\mathrm{-}\mathrm{0.275}\mathrm{-}\mathrm{0.0780}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{-}\mathrm{(}\mathrm{2}\mathit{/}\mathrm{3}\mathrm{)}\mathrm{\times }\mathrm{0.780}{\begin{array}{c}\mathrm{˜}\\ {\mathit{c}}_{\mathrm{3}}\end{array}}^{\mathrm{]}}\\ {\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}& \mathrm{=}& \end{array}$(87)Here the coefficients come from the Fisher matrix between the equilateral and orthogonal templates and the theoretical bispectra predicted by the two operators and . We use a χ² statistic given by ${\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathit{,}{\mathit{c}}_{\mathrm{s}}\mathrm{\right)}\mathrm{=}{{v}}^{\mathrm{T}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathit{,}{\mathit{c}}_{\mathrm{s}}\mathrm{\right)}{\mathit{C}}^{-1}{v}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathit{,}{\mathit{c}}_{\mathrm{s}}\mathrm{\right)}$, where ${\mathit{v}}^{\mathit{i}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathit{,}{\mathit{c}}_{\mathrm{s}}\mathrm{\right)}\mathrm{=}{\mathit{f}}^{\mathit{i}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathit{,}{\mathit{c}}_{\mathrm{s}}\mathrm{\right)}\mathrm{-}{\mathit{f}}_{\mathit{P}}^{\mathit{i}}$ (i = {equilateral, orthogonal}), ${\mathit{f}}_{\mathit{P}}^{\mathit{i}}$ being the joint estimates of equilateral and orthogonal f_NL (see Table 11), C the covariance matrix of the joint estimators, and ${\mathit{f}}^{\mathit{i}}\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathit{,}{\mathit{c}}_{\mathrm{s}}\mathrm{\right)}$ is provided by Eq. (87). As an example in Fig. 29 we show the 68%,95%, and 99.7% confidence regions for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}$ and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}$ obtained from the T + E constraints, requiring χ² ≤ 2.28,5.99, and 11.62 respectively (corresponding to a χ² variable with two degrees of freedom). In Fig. 30 we show the corresponding confidence regions in the $\mathrm{\left(}\begin{array}{c}{}_{\mathrm{˜}}\\ {\mathit{c}}_{\mathrm{3}}\end{array}\mathit{,}{\mathit{c}}_{\mathrm{s}}\mathrm{\right)}$ parameter space. Marginalizing over $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{c}}_{\mathrm{3}}\end{array}$ we find ${\mathit{c}}_{\mathrm{s}}\mathrm{\ge }\mathrm{0.020} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathrm{-}\mathrm{only}\mathrm{)}\mathit{,}$(88)and ${\mathit{c}}_{\mathrm{s}}\mathrm{\ge }\mathrm{0.024} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathit{+}\mathit{E}\mathrm{)}\mathit{.}$(89)The constraints improve by a few percent in T-only and by up to 25% by including polarization, in comparison with those of Planck Collaboration XXIV (2014).

Galileon models of inflation: Galileon models of inflation (Burrage et al. 2011; Kobayashi et al. 2010; Mizuno & Koyama 2010; Ohashi & Tsujikawa 2012) are well motivated models based on the so called “Galilean symmetry” (Nicolis et al. 2009). They are characterized by stability properties that are quite well understood (ghost-free, and stable against quantum corrections) and can arise naturally within fundamental physics models (de Rham & Gabadadze 2010b,a). Moreover they are an interesting example of models where gravity is modified on large scales and we focus on them also as a typical example of a more general class of modified gravity theories that are ghost-free (the so called Horndeski theories, Horndeski 1974; see also Planck Collaboration XIII 2016 where these models are discussed in the context of dark energy/modified gravity scenarios for the late time evolution of the universe). The predictions for the primordial perturbations are very rich. Bispectra can be generated with the same shapes as the “EFT1” and “EFT2” bispectra (see also discussion in Creminelli et al. 2011), however the amplitude(s) scale with the fluctuation sound speed as ${\mathit{c}}_{\mathrm{s}}^{-4}$, differently from the general single-field models of inflation considered in the above subsection. They can be written as (at the lowest-order in slow-roll parameters) $\begin{array}{ccc}{\mathit{f}}_{\mathrm{NL}}^{\mathrm{EFT}\mathrm{1}}& \mathrm{=}& \frac{\mathrm{17}}{\mathrm{972}}\left(\mathrm{-}\frac{\mathrm{5}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}}\mathrm{+}\frac{\mathrm{30}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{-}\frac{\mathrm{40}}{{\mathit{c}}_{\mathrm{s}}\mathit{c̅}\mathrm{s}}\mathrm{+}\mathrm{15}\right)\\ {\mathit{f}}_{\mathrm{NL}}^{\mathrm{EFT}\mathrm{2}}& \mathrm{=}& \frac{\mathrm{1}}{\mathrm{243}}(\frac{\mathrm{5}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{4}}}\mathrm{+}\frac{\mathrm{30}\mathit{/}\mathit{A}\mathrm{-}\mathrm{55}}{{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}}\mathrm{+}\frac{\mathrm{40}}{{\mathit{c}}_{\mathrm{s}}\mathit{c̅}\mathrm{s}}\mathrm{-}\mathrm{320}\frac{{\mathit{c}}_{\mathrm{s}}}{\mathit{c̅}\mathrm{s}}\mathrm{-}\frac{\mathrm{30}}{\mathit{A}}\mathrm{+}\mathrm{275}\\ & & \mathrm{-}\mathrm{225}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{+}\mathrm{280}\frac{{\mathit{c}}_{\mathrm{s}}^{\mathrm{3}}}{\mathit{c̅}\mathrm{s}})\mathrm{·}\end{array}$Here A, , and c_s are dimensionless parameters of the models. In particular c_s is the sound speed of the Galileon scalar field, while is a parameter that appears to break the standard consistency relation for the tensor-to-scalar perturbation ratio (, n_T being the tensor spectral index)²⁰. Accordingly to Eq. (87) a linear combination of these two bispectra generate equilateral and orthogonal bispectra templates²¹. From the Planck constraints on ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{equil}}$ and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{ortho}}$ (see Table 11),

we derive constraints on these model parameters following the procedure described at the beginning of this section. We choose log-constant priors in the ranges 10^-4 ≤ A ≤ 10⁴, and , together with a uniform prior 10^-4<c_s< 1. These priors have been choosen essentially on the basis of perturbative regime validity of the theory and to allow for a quite wide range of parameter values. In Fig. 31, as an example, probability contours are shown in the parameter space from the T + E constraints, after marginalizing over the parameter A. Marginalizing over both A and we find ${\mathit{c}}_{\mathrm{s}}^{\mathrm{Galileon}}\mathrm{\ge }\mathrm{0.21} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathrm{-}\mathrm{only}\mathrm{)}\mathit{,}$(92)and ${\mathit{c}}_{\mathrm{s}}^{\mathrm{Galileon}}\mathrm{\ge }\mathrm{0.23} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathit{+}\mathit{E}\mathrm{)}\mathit{.}$(93)Notice that interestingly enough the parameter can be negative in principle, corresponding to a blue spectral tilt of inflationary gravitational waves (without any kind of instability). We therefore also explore this branch, with a log-constant prior (for ), , and the same priors for the other parameters as above. Figure 32 shows the probability contours in the plane, after marginalizing over the parameter A, for the n_T< 0 branch. Marginalizing over both A and gives ${\mathit{c}}_{\mathrm{s}}^{\mathrm{Galileon}}\mathrm{\ge }\mathrm{0.19} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathrm{-}\mathrm{only}\mathrm{)}\mathit{,}$(94)and ${\mathit{c}}_{\mathrm{s}}^{\mathrm{Galileon}}\mathrm{\ge }\mathrm{0.21} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathit{+}\mathit{E}\mathrm{)}\mathit{.}$(95)A combined analysis of the Planck bispectrum and power spectrum constraints on the Galileon models is presented in the companion Planck paper on inflation (Planck Collaboration XX 2016).

11.2. Multi-field models

Curvaton models: the simplest adiabatic curvaton models predict local NG with an amplitude (Bartolo et al. 2004c,b) ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\frac{\mathrm{5}}{\mathrm{4}{\mathit{r}}_{\mathrm{D}}}\mathrm{-}\frac{\mathrm{5}{\mathit{r}}_{\mathrm{D}}}{\mathrm{6}}\mathrm{-}\frac{\mathrm{5}}{\mathrm{3}}\mathit{,}$(96)for a quadratic potential of the curvaton field (Lyth & Wands 2002; Lyth et al. 2003; Lyth & Rodriguez 2005; Malik & Lyth 2006; Sasaki et al. 2006), where r_D = [ 3ρ_curvaton/ (3ρ_curvaton + 4ρ_radiation) ] _D is the “curvaton decay fraction” at the time of the curvaton decay in the sudden decay approximation. Assuming a uniform prior, 0 <r_D< 1, our constraint ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{2.5}\mathrm{\pm }\mathrm{5.7}$ at 68% CL (see Table 11) yields ${\mathit{r}}_{\mathrm{D}}\mathrm{\ge }\mathrm{0.16} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathrm{-}\mathrm{only}\mathrm{)}\mathit{,}$(97)while accounting for temperature and polarization data (${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{0.8}\mathrm{\pm }\mathrm{5.0}$ at 68% CL) gives ${\mathit{r}}_{\mathrm{D}}\mathrm{\ge }\mathrm{0.19} \mathrm{95}\mathrm{\%CL}\mathrm{(}\mathit{T}\mathit{+}\mathit{E}\mathrm{)}\mathit{,}$(98)improving the previous Planck bound which was previously r_D ≥ 0.15 (95% CL; Planck Collaboration XXIV 2014). In Planck Collaboration XX (2016), assuming there is some relic isocurvature fluctuations in the curvaton field, a limit on r_D is derived from the bounds on isocurvature fluctuations. In this restricted case, the limit r_D> 0.98 (95% CL) is derived, which is consistent with the constraint given here.

Notice that the above expression of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ Eq. (96) is valid under the assumption that there is no significant decay of the inflaton into curvaton particles. In general one should account for such a possibility. For example, if the classical curvaton field survives and starts to dominate, then the curvaton particles produced during reheating (which have the same equation of state as the classical curvaton field) are expected to survive and dominate over other species at the epoch of their decay. The classical curvaton field and the curvaton particles decay at the same time, inevitably producing adiabatic perturbations (for a detailed discussion see Linde & Mukhanov 2006). A general formula for ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$, accounting for the possibility that the inflaton field decays into curvaton particles, is provided in Sasaki et al. (2006): ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{(}\mathrm{1}\mathrm{+}{\mathrm{\Delta }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{)}\frac{\mathrm{5}}{\mathrm{4}{\mathit{r}}_{\mathrm{D}}}\mathrm{-}\frac{\mathrm{5}{\mathit{r}}_{\mathrm{D}}}{\mathrm{6}}\mathrm{-}\frac{\mathrm{5}}{\mathrm{3}}\mathit{,}$(99)where ${\mathrm{\Delta }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{=}{\mathit{\rho }}_{\mathrm{curv}\mathit{.}\hspace{0.17em}\hspace{0.17em}\mathrm{particles}}\mathit{/}{\mathit{\rho }}_{\mathrm{curv}\mathit{.}\mathrm{field}}$ measures the ratio of the energy density of curvaton particles to the energy density of the classical curvaton field (Linde & Mukhanov 2006; Sasaki et al. 2006) and ρ_curvaton in the expression for r_D is given by ρ_curvaton = ρ_{curv.particles} + ρ_curv.field. Using the uniform prior 0 <r_D< 1 and $\mathrm{0}\mathit{<}{\mathrm{\Delta }}_{\mathrm{s}}^{\mathrm{2}}\mathit{<}{\mathrm{10}}^{\mathrm{2}}$ our measurements of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ constrain ${\mathrm{\Delta }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\le }\mathrm{8.5}$ at 95% CL (T) and ${\mathrm{\Delta }}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\le }\mathrm{6.9}$ (T+E).

11.3. Non-standard inflation models

Directional dependence motivated by gauge fields: in Table 25 we constrained directionally-dependent bispectra (Eq. (23)). This kind of NG is generated by inflationary models characterized by the presence of gauge fields. An actual realization of this type of scenario can be obtained with a coupling between the inflaton and the gauge field(s), via the kinetic term of the field(s), i.e., ℒ = − I²(φ)F²/ 4. In this formula, F² represents the strength of the gauge field, while I(φ) is a function of the inflaton field with an appropriate time dependence (see, e.g., Ratra 1992). In this type of scenario, vector fields can be generated during inflation, and this in turn determines the excitation of L = 0 and 2 modes in the bispectrum, with ${\mathit{f}}_{\mathrm{NL}}^{\mathit{L}}\mathrm{=}{\mathit{X}}_{\mathit{L}}\mathrm{\left(}\mathrm{\right|}{\mathit{g}}_{\mathrm{\ast }}\mathrm{|}\mathit{/}\mathrm{0.1}\mathrm{)}\hspace{0.17em}\mathrm{(}{\mathit{N}}_{{\mathit{k}}_{\mathrm{3}}}\mathit{/}\mathrm{60}\mathrm{)}$, where X_{L = 0} = (80 / 3) and X_{L = 2} = − (10 / 6) (Barnaby et al. 2012b; Bartolo et al. 2013a; Shiraishi et al. 2013a). The parameter g_∗, appearing in the equations above, represents the amplitude of a quadrupolar anisotropy in the power spectrum (see, e.g., Ackerman et al. 2007), while N defines the number of e-folds, before the end of inflation, when the relevant scales exit the horizon. It is thus clear that these models predict both a degree of statistical anisotropy in the power spectrum, and a potentially non-negligible bispectrum, as well as a direct relation between the two.

Starting from our SMICA constraints from T(T+E) in Table 25, marginalizing over a uniform prior 50 ≤ N ≤ 70, and assuming uniform priors on − 1 ≤ g_∗ ≤ 1, we obtain the limits − 0.050 <g_∗< 0.050 (− 0.040 <g_∗< 0.040), and − 0.31 <g_∗< 0.31 (− 0.29 <g_∗< 0.29), from the L = 0, L = 2 modes, respectively (95% CL) (considering g_∗ as scale-independent). We note that these constraints refer to all models in which curvature perturbations are sourced by an I²(φ)F² term (see references in Shiraishi et al. 2013a). The constraints we obtain are consistent with the tighter (model-independent) limits obtained in Planck Collaboration XX (2016) for the case of a scale-independent g_∗, from an analysis of quadrupolar anisotropies in the CMB power spectrum.

NG from gauge field production during axion inflation: we have also constrained the inverse decay NG of Eq. (25) arising typically in models where the inflaton field is a pseudoscalar axion that couples to a gauge field. Using the modal estimator we get the following constraints (removing ISW-lensing bias): ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{inv}\mathit{.}\mathrm{dec}}\mathrm{=}\mathrm{17}\mathrm{\pm }\mathrm{43} \mathrm{68}\mathrm{\%CL}\mathit{,}$(100)for temperature only; and ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{inv}\mathit{.}\mathrm{dec}}\mathrm{=}\mathrm{23}\mathrm{\pm }\mathrm{26} \mathrm{68}\mathrm{\%CL}\mathit{,}$(101)from temperature+polarization (see Table 23). The NG amplitude is given by ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{inv}\mathit{.}\mathrm{dec}}\mathrm{=}{\mathit{f}}_{\mathrm{3}}\mathrm{\left(}{\mathit{\xi }}_{\mathrm{\ast }}\mathit{,}\mathrm{1}\mathit{,}\mathrm{1}\mathrm{\right)}{\mathrm{𝒫}}_{\mathrm{\ast }}^{\mathrm{3}}{\mathit{e}}^{\mathrm{6}\mathit{\pi }{\mathit{\xi }}_{\mathrm{\ast }}}\mathit{/}{\mathrm{𝒫}}_{\mathit{\zeta }}^{\mathrm{2}}\mathrm{\left(}{\mathit{k}}_{\mathrm{\ast }}\mathrm{\right)}$, where is the power spectrum of vacuum-mode curvature perturbations (i.e., the power spectrum predicted without the coupling to gauge fields), ${\mathrm{𝒫}}_{\mathit{\zeta }}^{\mathrm{2}}\mathrm{\left(}{\mathit{k}}_{\mathrm{\ast }}\mathrm{\right)}$ is the dimensionless scalar power spectrum of curvature perturbations (a star denoting evaluation at the pivot scale). The NG parameter is exponentially sensitive to the strength of the coupling between the axion and the gauge field. From Eq. (101) we limit the strength of the coupling to ξ ≤ 2.5   (95% CL). The details together with constraints on the axion decay constant can be found in Planck Collaboration XX (2016), where an overview of various observational limits on axion (monodromy) models of inflation is presented. This limit is in agreement with the one that can be derived from tensor non-Gaussianities (see below).

Tensor NG and pseudoscalars: in inflationary scenarios associated with a pseudoscalar coupling to a U(1) gauge field, the tensor bispectrum generated via the gravitational interaction with the gauge field is expressed by Eq. (27), and the amplitude ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}$ depends on the following: the coupling strength of the pseudoscalar field to the gauge field (ξ); a slow-roll parameter for the inflaton (ϵ); and the power spectrum of vacuum-mode curvature perturbations (). The expression is ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}\mathrm{\approx }\mathrm{6.4}\mathrm{\times }{\mathrm{10}}^{\mathrm{11}}{\mathrm{𝒫}}^{\mathrm{3}}{\mathit{ϵ}}^{\mathrm{3}}{\mathit{e}}^{\mathrm{6}\mathit{\pi \xi }}\mathit{/}{\mathit{\xi }}^{\mathrm{9}}$ (Cook & Sorbo 2013; Shiraishi et al. 2013b). The constraints on ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{tens}}$ presented in Table 26 can then be used to constrain the model parameters. Clearly there are strong degeneracies, but if we marginalize over a uniform prior and set ϵ = 0.01, then assuming a uniform prior 0.1 ≤ ξ ≤ 7, from the SMICA (T-only or T+E) limit, we obtain ξ< 3.3 (95% CL).

Warm inflation: we update the constraints on warm inflation models in the strongly dissipative regime, when dissipative effects are relevant. In this regime ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{warm}}\mathrm{=}\mathrm{-}\mathrm{15}\ln {}^{\mathrm{(}}\mathrm{1}\mathrm{+}{\mathit{r}}_{\mathrm{d}}\mathit{/}{\mathrm{14}}^{\mathrm{)}}\mathrm{-}\mathrm{5}\mathit{/}\mathrm{2}$ (Moss & Xiong 2007) with a large dissipation parameter r_d = Γ / (3H) (where Γ is a friction term for the inflaton evolution describing the energy transfer from the inflaton field to radiation).. The limit from the 2013 Planck release is log ₁₀r_d ≤ 2.6 (95% CL) (Planck Collaboration XXIV 2014). Assuming a constant prior 0 ≤ log ₁₀r_d ≤ 4, the new SMICA constraint ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{warmS}}\mathrm{=}\mathrm{-}\mathrm{23}\mathrm{\pm }\mathrm{36}$ at 68% CL from T (${\mathit{f}}_{\mathrm{NL}}^{\mathrm{warmS}}\mathrm{=}\mathrm{-}\mathrm{14}\mathrm{\pm }\mathrm{23}$ from T+E), see Table 24, yields a limit on the dissipation parameter of log ₁₀r_d ≤ 3.3 (log ₁₀r_d ≤ 2.5) at 95% CL, with the T+E constraints (in brackets), slightly improving the 2013Planck limits. Values of r_d ≳ 2.5 (strongly-dissipative regime) are still allowed; however, the Planck constraint puts the model in a regime where there might be an overproduction of gravitinos (see Hall & Peiris 2008 and references therein). Unlike the strong dissipative regime, in the intermediate and weak dissipative regimes (r_D ≤ 1) the NG level strongly depends on the microscopic parameters (T/H and r_D), giving rise to a new additional bispectrum shape (for details see Bastero-Gil et al. 2014).

11.4. Alternatives to inflation

Ekpyrotic/cyclic models have been proposed as an alternative to inflation (for a review, see Lehners 2010). Local NG is generated from the conversion of “intrinsic” non-Gaussian entropy perturbation modes into curvature fluctuations. Models based on a conversion taking place during the ekpyrotic phase (the so called “ekpyrotic conversion mechanism”) are already ruled out (Koyama et al. 2007; Planck Collaboration XXIV 2014). Ekpyrotic models where “kinetic conversion” occurs after the ekpyrotic phase predict a local bispectrum with ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{=}\mathrm{(}\mathrm{3}\mathit{/}\mathrm{2}\mathrm{)}\hspace{0.17em}{\mathit{\kappa }}_{\mathrm{3}}\sqrt{\mathit{ϵ}}\mathrm{\pm }\mathrm{5}$, where the sign depends on the details of the conversion process (Lehners & Steinhardt 2008, 2013; Lehners 2010), where ϵ ≃ 50 or greater are typical values. If we take ϵ ≃ 100 and a uniform prior on − 5 <κ₃< 5 the constraints on ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ from T-only (see Table 11), yield − 0.91 <κ₃< 0.58 and − 0.25 <κ₃< 1.2 at 95% CL, for the plus and minus sign in ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ respectively. From the T+E ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ constraints (Table 11) we obtain − 0.94 <κ₃< 0.38 and − 0.27 <κ₃< 1.0 at 95% CL, for the plus and minus sign in ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}$ respectively. If we consider ϵ ≃ 50 we derive the following limits: − 1.3 <κ₃< 0.81 and − 0.35 <κ₃< 1.8 at 95% CL from T-only; − 1.3 <κ₃< 0.53 and − 0.38 <κ₃< 1.5 at 95% CL from T+E. Another variant of the ekpyrotic models has been investigated in Qiu et al. (2013), Li (2013), Fertig et al. (2014), where the intrinsic NG is zero and NG is generated only by non-linearities in the conversion mechanism, reaching a value of ${\mathit{f}}_{\mathrm{NL}}^{\mathrm{local}}\mathrm{\simeq }\mathrm{\pm }\mathrm{5}$.

11.5. Inflationary interpretation of CMB trispectrum results

We briefly interpret the trispectrum constraints in an inflationary context. First, consider the case of single field inflation. The action for the Goldstone boson π is highly constrained by residual diffeomorphism invariance (see, e.g., Smith et al. 2015). To lowest order in the derivative expansion, the most general action is: $\begin{array}{ccc}{\mathit{S}}_{\mathit{\pi }}& \mathrm{=}& \int {\mathrm{d}}^{\mathrm{4}}\mathit{x}\sqrt{\mathrm{-}\mathit{g}}\hspace{0.17em}\{\mathrm{-}{\mathit{M}}_{\mathrm{Pl}}^{\mathrm{2}}\mathit{Ḣ}{\mathrm{(}{\mathit{\partial }}_{}}^{}\end{array}$