Issue 
A&A
Volume 571, November 2014
Planck 2013 results



Article Number  A24  
Number of page(s)  58  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201321554  
Published online  29 October 2014 
Planck 2013 results. XXIV. Constraints on primordial nonGaussianity
^{1} APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de Paris, Sorbonne Paris Cité, 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France
^{2} Aalto University Metsähovi Radio Observatory, Metsähovintie 114, 02540 Kylmälä, Finland
^{3} African Institute for Mathematical Sciences, 68 Melrose Road, 7945 Muizenberg, Cape Town, South Africa
^{4} Agenzia Spaziale Italiana Science Data Center, via del Politecnico snc, 00133 Roma, Italy
^{5} Agenzia Spaziale Italiana, viale Liegi 26, 00998 Roma, Italy
^{6} Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{7} Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, 4000 Durban, South Africa
^{8} CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{9} CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{10} California Institute of Technology, Pasadena, California, USA
^{11} Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
^{12} Centro de Estudios de Física del Cosmos de Aragón (CEFCA), Plaza San Juan, 1, planta 2, 44001 Teruel, Spain
^{13} Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{14} Consejo Superior de Investigaciones Científicas (CSIC), 2800 Madrid, Spain
^{15} DSM/Irfu/SPP, CEASaclay, 91191 GifsurYvette Cedex, France
^{16} DTU Space, National Space Institute, Technical University of Denmark, Elektrovej 327, 2800 Kgs. Lyngby, Denmark
^{17} Département de Physique Théorique, Université de Genève, 24 Quai E. Ansermet, 1211 Genève 4, Switzerland
^{18} Departamento de Física Fundamental, Facultad de Ciencias, Universidad de Salamanca, 37008 Salamanca, Spain
^{19} Departamento de Física, Universidad de Oviedo, Avda. Calvo Sotelo s/n, 33007 Oviedo, Spain
^{20} Department of Astronomy and Astrophysics, University of Toronto, 50 Saint George Street, Toronto, Ontario, Canada
^{21} Department of Astrophysics/IMAPP, Radboud University Nijmegen, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{22} Department of Electrical Engineering and Computer Sciences, University of California, Berkeley, California, USA
^{23} Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{24} Department of Physics and Astronomy, Dana and David Dornsife College of Letter, Arts and Sciences, University of Southern California, Los Angeles CA 90089, USA
^{25} Department of Physics and Astronomy, University College London, London WC1E 6BT, UK
^{26} Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{27} Department of Physics, Florida State University, Keen Physics Building, 77 Chieftan Way, Tallahassee, Florida, USA
^{28} Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, 00014 Helsinki, Finland
^{29} Department of Physics, Princeton University, Princeton, New Jersey, USA
^{30} Department of Physics, University of California, Berkeley, California, USA
^{31} Department of Physics, University of California, One Shields Avenue, Davis, California, USA
^{32} Department of Physics, University of California, Santa Barbara, California, USA
^{33} Department of Physics, University of Illinois at UrbanaChampaign, 1110 West Green Street, Urbana, Illinois, USA
^{34} Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, via Marzolo 8, 35131 Padova, Italy
^{35} Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, via Saragat 1, 44122 Ferrara, Italy
^{36} Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, 00185 Roma, Italy
^{37} Dipartimento di Fisica, Università degli Studi di Milano, via Celoria, 16, 20133 Milano, Italy
^{38} Dipartimento di Fisica, Università degli Studi di Trieste, via A. Valerio 2, 34127 Trieste, Italy
^{39} Dipartimento di Fisica, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, 00133 Roma, Italy
^{40} Dipartimento di Matematica, Università di Roma Tor Vergata, via della Ricerca Scientifica 1, 00133 Roma, Italy
^{41} Discovery Center, Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark
^{42} Dpto. Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{43} European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo s/n, Urbanización Villafranca del Castillo, 28691 Villanueva de la Cañada, Madrid, Spain
^{44} European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{45} Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, 00014 Helsinki, Finland
^{46} INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, 35122 Padova, Italy
^{47} INAF – Osservatorio Astronomico di Roma, via di Frascati 33, 00040 Monte Porzio Catone, Italy
^{48} INAF – Osservatorio Astronomico di Trieste, via G.B. Tiepolo 11, 34143 Trieste, Italy
^{49} INAF Istituto di Radioastronomia, via P. Gobetti 101, 40129 Bologna, Italy
^{50} INAF/IASF Bologna, via Gobetti 101, 40129 Bologna, Italy
^{51} INAF/IASF Milano, via E. Bassini 15, 20133 Milano, Italy
^{52} INFN, Sezione di Bologna, via Irnerio 46, 40126, Bologna, Italy
^{53} INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy
^{54} IPAG: Institut de Planétologie et d’Astrophysique de Grenoble, Université Joseph Fourier, Grenoble 1/CNRSINSU, UMR 5274, 38041 Grenoble, France
^{55} IUCAA, Post Bag 4, Ganeshkhind, Pune University Campus, 411 007 Pune, India
^{56} Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London, SW7 2AZ, UK
^{57} Infrared Processing and Analysis Center, California Institute of Technology, Pasadena CA 91125, USA
^{58} Institut Néel, CNRS, Université Joseph Fourier Grenoble I, 25 rue des Martyrs, 38042 Grenoble, France
^{59} Institut Universitaire de France, 103 bd SaintMichel, 75005 Paris, France
^{60} Institut d’Astrophysique Spatiale, CNRS (UMR 8617) Université ParisSud 11, Bâtiment 121, 91405 Orsay, France
^{61} Institut d’Astrophysique de Paris, CNRS (UMR 7095), 98bis Boulevard Arago, 75014 Paris, France
^{62} Institute for Space Sciences, 077125 BucharestMagurale, Romania
^{63} Institute of Astronomy and Astrophysics, Academia Sinica, 106 Taipei, Taiwan
^{64} Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{65} Institute of Theoretical Astrophysics, University of Oslo, Blindern, 0315 Oslo, Norway
^{66} Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, 38200 La Laguna, Tenerife, Spain
^{67} Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, 39005 Santander, Spain
^{68} Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, California, USA
^{69} Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester, M13 9PL, UK
^{70} Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge, CB3 0HA, UK
^{71} LAL, Université ParisSud, CNRS/IN2P3, 91383 Orsay, France
^{72} LERMA, CNRS, Observatoire de Paris, 61 Avenue de l’Observatoire, 75014 Paris, France
^{73} Laboratoire AIM, IRFU/Service d’Astrophysique – CEA/DSM – CNRS – Université Paris Diderot, Bât. 709, CEASaclay, 91191 GifsurYvette Cedex, France
^{74} Laboratoire Traitement et Communication de l’Information, CNRS (UMR 5141) and Télécom ParisTech, 46 rue Barrault, 75634 Paris Cedex 13, France
^{75} Laboratoire de Physique Subatomique et de Cosmologie, Université Joseph Fourier Grenoble I, CNRS/IN2P3, Institut National Polytechnique de Grenoble, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{76} Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{77} Lawrence Berkeley National Laboratory, Berkeley, California, USA
^{78} MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{79} McGill Physics, Ernest Rutherford Physics Building, McGill University, 3600 rue University, Montréal, QC, H3A 2T8, Canada
^{80} MilliLab, VTT Technical Research Centre of Finland, Tietotie 3, 02044 Espoo, Finland
^{81} National University of Ireland, Department of Experimental Physics, Maynooth, Co. Kildare, Ireland
^{82} Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark
^{83} Observational Cosmology, Mail Stop 36717, California Institute of Technology, Pasadena CA, 91125, USA
^{84} Optical Science Laboratory, University College London, Gower Street, London, UK
^{85} SBITPLPPC, EPFL, 1015 Lausanne, Switzerland
^{86} SISSA, Astrophysics Sector, via Bonomea 265, 34136 Trieste, Italy
^{87} School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff, CF24 3AA, UK
^{88} School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{89} Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, 117997 Moscow, Russia
^{90} Space Sciences Laboratory, University of California, Berkeley, California, USA
^{91} Special Astrophysical Observatory, Russian Academy of Sciences, Nizhnij Arkhyz, Zelenchukskiy region, 369167, KarachaiCherkessian Republic, Russia
^{92} Stanford University, Dept of Physics, Varian Physics Bldg, 382 via Pueblo Mall, Stanford, California, USA
^{93} SubDepartment of Astrophysics, University of Oxford, Keble Road, Oxford OX1 3RH, UK
^{94} Theory Division, PHTH, CERN, 1211 Geneva 23, Switzerland
^{95} UPMC Univ. Paris 06, UMR7095, 98bis Boulevard Arago, 75014 Paris, France
^{96} Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{97} University of Granada, Departamento de Física Teórica y del Cosmos, Facultad de Ciencias, 18071 Granada, Spain
^{98} Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
Received: 22 March 2013
Accepted: 16 December 2013
The Planck nominal mission cosmic microwave background (CMB) maps yield unprecedented constraints on primordial nonGaussianity (NG). Using three optimal bispectrum estimators, separable templatefitting (KSW), binned, and modal, we obtain consistent values for the primordial local, equilateral, and orthogonal bispectrum amplitudes, quoting as our final result f_{NL}^{local} = 2.7 ± 5.8, f_{NL}^{equil} = 42 ± 75, and f_{NL}^{orth} = 25 ± 39 (68% CL statistical). NonGaussianity is detected in the data; using skewC_{ℓ} statistics we find a nonzero bispectrum from residual point sources, and the integratedSachsWolfelensing bispectrum at a level expected in the ΛCDM scenario. The results are based on comprehensive crossvalidation of these estimators on Gaussian and nonGaussian simulations, are stable across component separation techniques, pass an extensive suite of tests, and are confirmed by skewC_{ℓ}, wavelet bispectrum and Minkowski functional estimators. Beyond estimates of individual shape amplitudes, we present modelindependent, threedimensional reconstructions of the Planck CMB bispectrum and thus derive constraints on earlyUniverse scenarios that generate primordial NG, including general singlefield models of inflation, excited initial states (nonBunchDavies vacua), and directionallydependent vector models. We provide an initial survey of scaledependent feature and resonance models. These results bound both general singlefield and multifield model parameter ranges, such as the speed of sound, c_{s} ≥ 0.02 (95% CL), in an effective field theory parametrization, and the curvaton decay fraction r_{D} ≥ 0.15 (95% CL). The Planck data significantly limit the viable parameter space of the ekpyrotic/cyclic scenarios. The amplitude of the fourpoint function in the local model τ_{NL}< 2800 (95% CL). Taken together, these constraints represent the highest precision tests to date of physical mechanisms for the origin of cosmic structure.
Key words: cosmic background radiation / cosmology: observations / cosmology: theory / early Universe / inflation / methods: data analysis
© ESO, 2014
1. Introduction
This paper, one of a set associated with the 2013 release of data from the Planck^{1} mission (Planck Collaboration I 2014), describes the constraints on primordial nonGaussianity (NG) obtained using the cosmic microwave background (CMB) maps derived from the data acquired by Planck during its nominal operations period, i.e., between 12 August 2009 and 27 November 2010.
Primordial NG is one of the most informative fingerprints of the origin of structure in the Universe, probing physics at extremely high energy scales inaccessible to laboratory experiments. Possible departures from a purely Gaussian distribution of the CMB anisotropies provide powerful observational access to this extreme physics (Allen et al. 1987; Salopek & Bond 1990; Falk et al. 1993; Gangui et al. 1994; Verde et al. 2000; Gangui & Martin 2000b; Wang & Kamionkowski 2000; Komatsu & Spergel 2001; Acquaviva et al. 2003; Maldacena 2003; Babich et al. 2004; for recent reviews Bartolo et al. 2004a; Liguori et al. 2010; Chen 2010b; Komatsu 2010; Yadav & Wandelt 2010). A robust detection of primordial NG – or a strong constraint on it – discriminates among competing mechanisms for the generation of the cosmological perturbations in the early Universe. Different inflationary models, firmly rooted in modern theoretical particle physics, predict different amplitudes, shapes, and scale dependence of NG. As a result, primordial NG is complementary to the scalarspectral index of curvature perturbations and the tensortoscalar amplitude ratio, distinguishing between inflationary models that are degenerate on the basis of their power spectra alone. Even in the simplest models of inflation, consisting of a single slowlyrolling scalar field, a small (but calculable) level of NG is predicted (Acquaviva et al. 2003; Maldacena 2003); this is undetectable in presentquality CMB and largescale structure measurements. However, as demonstrated by a large body of work in recent years, extending this simplest paradigm will generically lead to detectable levels of NG in CMB anisotropies. Critically, a robust detection of primordial NG would rule out all canonical singlefield slowroll models of inflation, pointing to physics beyond the simplest “textbook” picture of inflation. Conversely, significant improvements in the constraints on primordial NG strongly limit extensions to the simplest paradigm, thus providing powerful clues to the physical mechanism that generated cosmic structure.
If the primordial fluctuations are Gaussiandistributed, then they are completely characterised by their twopoint correlation function, or equivalently, their power spectrum. If they are nonGaussian, there is additional statistical information in the higherorder correlation functions, which is not captured by the twopoint correlation function. In particular, the 3point correlation function, or its Fourier counterpart, the bispectrum, is important because it is the lowestorder statistic that can distinguish between Gaussian and nonGaussian perturbations. One of the main goals of this paper is to constrain the amplitude and shape of primordial NG using the angular bispectrum of the CMB anisotropies. The CMB angular bispectrum is related to the primordial bispectrum defined by (1)Here we define the potential Φ in terms of the comoving curvature perturbation ζ on superhorizon scales by Φ ≡ (3/5)ζ. In matter domination, on superhorizon scales, Φ is equivalent to Bardeen’s gaugeinvariant gravitational potential (Bardeen 1980), and we adopt this notation for historical consistency. The bispectrum B_{Φ}(k_{1},k_{2},k_{3}) measures the correlation among three perturbation modes. Assuming translational and rotational invariance, it depends only on the magnitudes of the three wavevectors. In general the bispectrum can be written as (2)Here, f_{NL} is the socalled “nonlinearity parameter” (Gangui et al. 1994; Wang & Kamionkowski 2000; Komatsu & Spergel 2001; Babich et al. 2004), a dimensionless parameter measuring the amplitude of NG. The bispectrum is measured by sampling triangles in Fourier space. The dependence of the function F(k_{1},k_{2},k_{3}) on the type of triangle (i.e., the configuration) formed by the three wavevectors describes the shape of the bispectrum (Babich et al. 2004), which encodes much physical information. It can also encode the scale dependence, i.e., the running, of the bispectrum (Chen 2005c)^{2}. Different NG shapes are linked to distinctive physical mechanisms that can generate such nonGaussian fingerprints in the early Universe. For example, the socalled “local” NG (Gangui et al. 1994; Verde et al. 2000; Wang & Kamionkowski 2000; Komatsu & Spergel 2001) is characterized by a signal that is maximal for “squeezed” triangles with k_{1} ≪ k_{2} ≃ k_{3} (or permutations; Maldacena 2003) which occurs, in general, when the primordial NG is generated on superhorizon scales. Conversely, “equilateral” NG (Babich et al. 2004) peaks for equilateral configurations k_{1} ≈ k_{2} ≈ k_{3}, due to correlations between fluctuation modes that are of comparable wavelengths, which can occur if the three perturbation modes mostly interact when they cross the horizon approximately at the same time. Other relevant shapes include the socalled “folded” (or flattened) NG (Chen et al. 2007b), which is due to correlations between perturbation modes that are enhanced for k_{1} + k_{2} ≈ k_{3}, or the “orthogonal” NG (Senatore et al. 2010) that generates a signal with a positive peak at the equilateral configuration and a negative peak at the folded configuration.
We now sketch how nonGaussian information in the initial conditions is transferred to observable quantities (in this instance, the CMB anisotropies) in the context of inflation. Primordial perturbations in the inflaton field(s) φ(x,t) = φ_{0}(t) + δφ(x,t) (where δφ denotes quantum fluctuations about the background value φ_{0}(t)) can be characterized by the comoving curvature perturbation ζ, since this is conserved on superhorizon scales for adiabatic perturbations. The inflaton fluctuations δφ (in the flat gauge) induce a curvature perturbation^{3} at linear order; however, nonlinearities induce corrections to this relation. The primordial NG in the curvature perturbation ζ is intrinsically nonlinear, so that its contribution to the CMB anisotropies is transferred linearly at leading order. In particular, at the linear level, the curvature perturbation ζ is related to Bardeen’s gravitational potential Φ during the matterdominated epoch by Φ = (3/5)ζ and ΔT/T ~ gζ, where g is the linear radiation transfer function; thus, any primordial NG will be transferred to the CMB even at linear order. For example, in the largeangular scale limit, the linear SachsWolfe effect reads ΔT/T = −Φ/3 = −ζ/ 5. Further, any other field excited during the inflationary phase which develops quantum fluctuations contributing to the primordial curvature perturbation – whether or not it is driving inflation – can leave its nonGaussian imprint in the CMB anisotropies.
Thus the bispectrum of Eq. (1) measures the fundamental (self) interactions of the scalar field(s) involved in the inflationary phase and/or generating the primordial curvature perturbation, as well as measuring nonlinear processes occurring during or after inflation. It therefore brings insights into the fundamental physics behind inflation, possibly allowing for the first time a reconstruction of the inflationary Lagrangian itself. For example, in a large class of inflationary models which involve additional light field(s) different from the inflaton, the superhorizon evolution of the fluctuations in the additional field(s) and their transfer to the adiabatic curvature perturbations can generate a large primordial NG of the local type. This is the case for curvatontype models (Linde & Mukhanov 1997; Lyth & Wands 2002; Lyth et al. 2003) where the latetime decay of a scalar field, belonging to the noninflationary sector of the theory, induces curvature perturbations; models where the curvature perturbation is generated by the local fluctuations of the inflaton’s coupling to matter during the reheating phase (Kofman 2003; Dvali et al. 2004a); and multifield models of inflation (see, e.g., Bartolo et al. 2002; Bernardeau & Uzan 2002; Vernizzi & Wands 2006; Rigopoulos et al. 2006, 2007; Lyth & Rodriguez 2005; Byrnes & Choi 2010). Since the nonlinear processes take place on superhorizon scales, the form of NG is local in real space and thus, in Fourier space, the bispectrum correlates large and small Fourier modes. “Equilateral” NG (Babich et al. 2004) is a generic feature of singlefield models with a noncanonical kinetic term, which can also generate the “orthogonal” type of NG (Senatore et al. 2010). In general, these models are characterized by higherderivative interactions of the inflaton field. The correlation between the fluctuation modes is suppressed when one of the modes is on superhorizon scales, because the derivative terms are redshifted away, so that the correlation is maximal for three modes of comparable wavelengths that cross the horizon at the same time. An example of “folded” NG is the one generated in a class of singlefield models with nonBunchDavies vacuum (Chen et al. 2007b; Holman & Tolley 2008). Indeed, these and other types of primordial NG can also be produced in other models, and we refer to Sect. 2 for more details. All these models can easily yield primordial NG with an amplitude much bigger than the one predicted in the standard models of singlefield slowroll inflation, for which the NG amplitude turns out to be proportional to the usual slowroll parameters (Acquaviva et al. 2003; Maldacena 2003).
Given that a robust detection of primordial NG would represent a breakthrough in the understanding of the physics governing the Universe during its very first stages, it is crucial that all sources of contamination are sufficiently understood to firmly control their effects. In particular, any nonlinearity in the postinflationary Universe can introduce NG into perturbations that were initially Gaussian. Therefore, one must ensure that a primordial origin is not ascribed to a nonprimordial contaminant; however, estimators of (primordial) NG from CMB data will also typically be sensitive to such contaminating signals. Potential nonprimordial sources of NG can be classified into four broad categories: instrumental systematic effects (see e.g., Donzelli et al. 2009); residual foregrounds and point sources; secondary CMB anisotropies, such as the SunyaevZeldovich (SZ) effect (Zeldovich & Sunyaev 1969), gravitational lensing (see Lewis & Challinor 2006, for a review), the integrated SachsWolfe (ISW) effect (Sachs & Wolfe 1967) or the ReesSciama effect (Rees & Sciama 1968); and effects arising from nonlinear (secondorder) perturbations in the Boltzmann equations (due to the nonlinear nature of General Relativity and the nonlinear dynamics of the photonbaryon fluid at recombination). Among the secondary anisotropies, the crosscorrelation of the ISW/ReesSciama and lensing (Goldberg & Spergel 1999) produces the dominant contamination to the (local) primordial NG. The impact is mainly on the local type of primordial NG, because the ISWlensing correlation couples the largescale gravitational potential fluctuations sourcing the ISW effect with the smallscale lensing effects of the CMB, thus producing a bispectrum which peaks on the squeezed configurations, as for the local shape. Detailed analyses have shown that the ISWlensing bispectrum can introduce a bias to local primordial NG, while the bias to equilateral primordial NG is negligible (see Serra & Cooray 2008; Smith & Zaldarriaga 2011; Hanson et al. 2009b; Lewis et al. 2011, 2012; Mangilli & Verde 2009; Junk & Komatsu 2012; Mangilli et al. 2013). In our analysis we have carefully accounted for this effect (we report the values of the ISWlensing bias in Sect. 5.2, and demonstrate the detection of the effect with skewC_{ℓ}s), as well as validating our results through an extensive suite of simulations and null tests in order to quantify the effects of systematic effects and diffuse and pointsource foregrounds. Finally, a consistent treatment of weak NG in the CMB must account for additional contributions that arise at the nonlinear (secondorder) level both in the gravitational perturbations after inflation ends, and for the evolution of the CMB anisotropies at secondorder in perturbation theory at large and small angular scales. It has been shown that these secondorder CMB effects yield negligible contamination to primordial NG for Planckquality data (Bartolo et al. 2004b, 2005, 2010c, 2012; Creminelli & Zaldarriaga 2004b; Boubekeur et al. 2009; Nitta et al. 2009; Senatore et al. 2009; Khatri & Wandelt 2009; Bartolo & Riotto 2009; Khatri & Wandelt 2010; Bartolo et al. 2010c; Creminelli et al. 2011b; Bartolo et al. 2012; Huang & Vernizzi 2013; Su et al. 2012; Pettinari et al. 2013).
Previous constraints on various shapes of primordial NG come from the WMAP9 data (Bennett et al. 2013). For the local shape they find (68% CL). For equilateraltype NG, they obtain (68% CL), while for the orthogonal shape (68% CL). Other analyses employing different estimators give compatible constraints. Limits on other shapes, such as e.g. flattened and feature models, have also been obtained (Fergusson et al. 2012).
Before concluding this section let us point out the connection between the analyses presented here and in the companion paper Planck Collaboration XXIII (2014) on the statistical and isotropy properties of the CMB. Statistical anisotropy and NG are essentially two alternative descriptions of the same phenomenon on the sky (Ferreira & Magueijo 1997). Specifically any Gaussian but statistically anisotropic model becomes, after averaging over the possible (a priori unknown) orientations of the anisotropy, a statistically isotropic nonGaussian model. For example local NG can be generated by largescale field fluctuations that couple to the smallscale power. For the given fixed realization of largescale modes that we see, the smallscale anisotropies look anisotropic on the sky, and it is equally valid to describe this as a Gaussian anisotropic model (assuming the largescale modes are Gaussian). In this paper we mostly focus on the nonGaussian interpretation of various physically motivated models, although it is useful to bear both perspectives in mind, in particular when considering what forms of nonprimordial signal might cause contamination. Planck Collaboration XXIII (2014) consider a broad class of more general phenomenological forms of anisotropy, which are complementary to the analysis presented here.
This paper is organized as follows. In Sect. 2, we present models generating primordial NG that have been tested in this paper. Section 3 summarizes the statistical estimators used to constrain the CMB bispectrum from Planck data and the methods for the reconstruction of the CMB bispectrum. Section 4 summarizes the statistical estimator used to constrain the CMB trispectrum. In Sect. 5, we discuss the nonprimordial contributions to the CMB bispectrum and trispectrum, including foreground residuals after component separation and focusing on the f_{NL} bias induced by the ISWlensing bispectrum. Section 6 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 componentseparation techniques. Section 7 contains our main results: we present constraints on f_{NL} for the local, equilateral, and orthogonal bispectra, and a selected set of other bispectrum shapes; we show a reconstruction of the CMB bispectrum, and give limits on the CMB trispectrum. In Sect. 8 we validate these results by performing a series of null tests on the data to assess the robustness of our results. We also evaluate the impact of the Planck data processing on the primordial NG signal. In Sect. 9, we discuss the main implications of Planck’s constraints on primordial NG for early Universe models. We conclude in Sect. 10. The realistic Planck simulations used in various steps of the analysis and validation tests are described in Appendix A. Appendix B contains a derivation of the expected scatter between f_{NL} results on the same map from different estimators used in the validation tests of Sect. 6, while Appendix C presents a comparison of constraints on some selected nonstandard bispectrum shapes using different foregroundcleaned maps.
2. Inflationary models for primordial nonGaussianity
There is a simple reason why standard singlefield models of slowroll inflation predict a tiny level of NG, of the order of the usual slowroll parameters ^{4}: in order to achieve an accelerated period of expansion, the inflaton potential must be very flat, thus suppressing the inflaton (self)interactions and any sources of nonlinearity, and leaving only its weak gravitational interactions as the main source of NG. This fact leads to a clear distinction between the simplest models of inflation, and scenarios where a significant amplitude of NG can be generated (e.g., Komatsu 2010), as follows. The simplest inflationary models are based on a set of minimal conditions: (i) a single weaklycoupled neutral single scalar field (the inflaton, which drives inflation and generates the curvature perturbations); (ii) with a canonical kinetic term; (iii) slowly rolling down its (featureless) potential; (iv) initially lying in a BunchDavies (ground) vacuum state. In the last few years, an important theoretical realization has taken place: a detectable amplitude of NG with specific triangular configurations (corresponding broadly to wellmotivated classes of physical models) can be generated if any one of the above conditions is violated (Bartolo et al. 2004a; Liguori et al. 2010; Chen 2010b; Komatsu 2010; Yadav & Wandelt 2010):

“local” NG, where the signal peaks in “squeezed” triangles (k_{1} ≪ k_{2} ≃ k_{3})(e.g., multifield models of inflation);

“equilateral” NG, peaking for k_{1} ≈ k_{2}≈ k_{3}. Examples of this class include singlefield models with noncanonical kinetic term (Chen et al. 2007b), such as kinflation (ArmendarizPicon et al. 1999; Chen et al. 2007b) or DiracBornInfield (DBI) inflation (Silverstein & Tong 2004; Alishahiha et al. 2004), models characterized by more general higherderivative interactions of the inflaton field, such as ghost inflation (ArkaniHamed et al. 2004), and models arising from effective field theories (Cheung et al. 2008);

“folded” (or flattened) NG. Examples of this class include: singlefield models with nonBunchDavies vacuum (Chen et al. 2007b; Holman & Tolley 2008) and models with general higherderivative interactions (Senatore et al. 2010; Bartolo et al. 2010a);

“orthogonal” NG which is generated, e.g., in singlefield models of inflation with a noncanonical kinetic term (Senatore et al. 2010), or with general higherderivative interactions.
All these models naturally predict values of  f_{NL}  ≫ 1. A detection of such a signal would rule out the simplest models of singlefield inflation, which, obeying all the conditions above, are characterized by weak gravitational interactions with  f_{NL}  ≪ 1.
The above scheme provides a general classification of inflationary models in terms of the corresponding NG shapes, which we adopt for the data analysis presented in this paper:

1.
“general” singlefield inflationary models (tested using theequilateral, orthogonal and folded shapes);

2.
multifield models of inflation (tested using the local shape).
In each class, there exist specific realizations of inflationary models which are characterized by the same underlying physical mechanism, generating a specific NG shape. We will investigate these classes of inflationary models by constraining the corresponding NG content, focusing on amplitudes and shapes. We also perform a survey of nonstandard models giving rise to alternative specific shapes of NG. Different NG shapes are observationally distinguishable if their crosscorrelation is sufficiently low; almost all of the shapes analysed in this paper are highly orthogonal to each other (e.g., Babich et al. 2004; Fergusson & Shellard 2007).
There are exceptional cases which evade this classification: for example, some exotic nonlocal singlefield theories of inflation produce local NG (Barnaby & Cline 2008), while some multifield models can produce equilateral NG, e.g., if some particle production mechanism is present (examples include trapped inflation Green et al. 2009, and some models of axion inflation Barnaby & Peloso 2011; Barnaby et al. 2011, 2012b). Another example arises in a class of multifield models where the second scalar field is not light, but has a mass m ≈ H, of the order of the Hubble rate during inflation. Then NG with an intermediate shape, interpolating between local and equilateral, can be produced – “quasisingle field” models of inflation (Chen & Wang 2010a,b) – for which the NG shape is similar to the socalled constant NG of Fergusson & Shellard (2007). Furthermore, there is the possibility of a superposition of shapes (and/or running of NG), generated if different mechanisms sourcing NG act simultaneously during the inflationary evolution. For example, in multifield DBI inflation, equilateral NG is generated at horizon crossing from the higherderivative interactions of the scalar fields, and it adds to the local NG arising from the superhorizon nonlinear evolution (e.g., Langlois et al. 2008a,b; Arroja et al. 2008; RenauxPetel 2009).
In the following subsections, we discuss each of these possibilities in turn. The reader already familiar with this background material may skip to Sect. 3.
2.1. General singlefield models of inflation
Typically in models with a nonstandard kinetic term (or more general higherderivative interactions), inflaton perturbations propagate with an effective sound speed c_{s} which can be smaller than the speed of light, and this results in a contribution to the NG amplitude in the limit c_{s} ≪ 1. For example, models with a nonstandard kinetic term are described by an inflaton Lagrangian ℒ = P(X,φ), where X = g^{μν}∂_{μ}φ∂_{ν}φ, with at most one derivative on φ, and the sound speed is .
In this case, two interaction terms give the dominant contribution to primordial NG, one of the type and the other of the type , which arise from expanding the P(X,φ) Lagrangian. Each of these two interaction terms generates a bispectrum with a shape similar to the equilateral type, with the second inflaton interaction yielding a nonlinearity parameter , independent of the amplitude of the other bispectrum. Equilateral NG is usually generated by derivative interactions of the inflaton field; derivative terms are suppressed when one perturbation mode is frozen on superhorizon scales during inflation, and the other two are still crossing the horizon, so that the correlation between the three perturbation modes will be suppressed, while it is maximal when all the three modes cross the horizon at the same time, which happens for k_{1} ≈ k_{2} ≈ k_{3}.
The equilateral type NG is well approximated by the template (Creminelli et al. 2006) (3)
where P_{Φ}(k) = A/k^{4 − ns} is the power spectrum of Bardeen’s gravitational potential with normalization A^{2} and scalar spectral index n_{s}. For example, the models introduced in the string theory framework based on the DBI action (Silverstein & Tong 2004; Alishahiha et al. 2004) can be described within the P(X,φ)class, and they give rise to an equilateral NG with an overall amplitude for c_{s} ≪ 1, which turns out typically to be ^{5}.
The equilateral shape emerges also in models characterized by more general higherderivative interactions, such as ghost inflation (ArkaniHamed et al. 2004) or models within effective field theories of inflation (Cheung et al. 2008; Senatore et al. 2010; Bartolo et al. 2010a).
Taken individually, each higherderivative interaction of the inflaton field generically gives rise to a bispectrum with a shape which is similar – but not identical to – the equilateral form (an example is provided by the two interaction terms discussed above for an inflaton with a nonstandard kinetic term). Therefore it has been shown, using an effective field theory approach to inflationary perturbations, that it is possible to build a combination of the corresponding similar equilateral shapes to generate a bispectrum that is orthogonal to the equilateral one, the socalled “orthogonal” shape. This can be approximated by the template (Senatore et al. 2010) (4)
The orthogonal bispectrum can also arise as the predominant shape in some inflationary realizations of Galileon inflation (RenauxPetel et al. 2011).
Nonseparable singlefield bispectrum shapes: while most singlefield inflation bispectra can be wellcharacterized by the equilateral and orthogonal shapes, we note that these are separable ansätze which only approximate the contributions from two leading order terms in the cubic Lagrangian. In an effective field theory approach these correspond to two shapes which can be associated directly with the inflaton field interactions and (in the language of the effective field theory of inflation the inflaton scalar degree of freedom π is related to the comoving curvature perturbation as ζ = − Hπ). They are, respectively (Senatore et al. 2010; see also Chen et al. 2007b^{6}) These shapes differ from equilateral in the flattened or collinear limit. DBI inflation gives a closely related shape of particular interest phenomenologically (Alishahiha et al. 2004), (7)For brevity, we have given the scaleinvariant form of the shape functions, without the mild power spectrum running. There are also subleading order terms which give rise to additional nonseparable shapes, but these are expected to be much smaller without special finetuning.
2.2. Multifield models
This class of models generally includes an additional light scalar field (or more fields) during inflation, which can be different from the inflaton, and whose fluctuations contribute to the final primordial curvature perturbation of the gravitational potential. It could be the case of inflation driven by several scalar fields – “multiplefield inflation” – or the one where the inflaton drives the accelerated expansion, while other scalar fields remain subdominant during inflation. This encompasses, for instance, a large class of multifield models which leads to nonGaussian isocurvature perturbations (for earlier works, see e.g., Linde & Mukhanov 1997; Peebles 1997; Bucher & Zhu 1997). More importantly, such models can also lead to crosscorrelated and nonGaussian adiabatic and isocurvature modes, where NG is first generated by large nonlinearities in some scalar (possibly noninflatonic) sector of the theory, and then efficiently transferred to the inflaton adiabatic sector(s) through the crosscorrelation of adiabatic and isocurvature perturbations^{7} (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 multiplefield inflation models, see, Byrnes & Choi 2010). Another interesting possibility is the curvaton model (Mollerach 1990; Enqvist & Sloth 2002; Lyth & Wands 2002; Moroi & Takahashi 2001), where a second light scalar field, subdominant during inflation, decays after inflation ends, producing the primordial density perturbations which can be characterized by a high NG level (e.g., Lyth & Wands 2002; Lyth et al. 2003; Bartolo et al. 2004d). NG in the curvature perturbation can be generated at the end of inflation, e.g., due to the nonlinear dynamics of (p)reheating (e.g., Enqvist et al. 2005; Chambers & Rajantie 2008; Barnaby & Cline 2006; see also Bond et al. 2009) or, as in modulated (p)reheating and modulated hybrid inflation, due to local fluctuations in the decay rate/interactions of the inflaton field (Kofman 2003; Dvali et al. 2004a,b; Bernardeau et al. 2004; Zaldarriaga 2004; Lyth 2005; Salem 2005; Lyth & Riotto 2006; Kolb et al. 2006; Cicoli et al. 2012). The common feature of all these models is that a large NG in the curvature perturbation can be produced via both a transfer of superhorizon nonGaussian isocurvature perturbations in the second field (not necessarily the inflaton) to the adiabatic density perturbations, and via additional nonlinearities in the transfer mechanism. Since, typically, this process occurs on superhorizon scales, the form of NG is local in real space. Being local in real space, the bispectrum correlates large and small scale Fourier modes. The local bispectrum is given by (Falk et al. 1993; Gangui et al. 1994; Verde et al. 2000; Wang &Kamionkowski 2000; Komatsu & Spergel 2001) (8)Most of the signaltonoise ratio in fact peaks in the squeezed configurations (k_{1} ≪ k_{2} ≃ k_{3}) (9)The typical example of a curvature perturbation that generates the bispectrum of Eq. (8) is the standard local form for the gravitational potential (Hodges et al. 1990; Kofman et al. 1991; Salopek & Bond 1990; Gangui et al. 1994; Verde et al. 2000; Wang & Kamionkowski 2000; Komatsu & Spergel 2001) (10)where Φ_{L}(x) is the linear Gaussian gravitational potential and is the amplitude of a quadratic nonlinear correction (though this is not the only possibility: e.g., the gravitational potential produced in multiplefield inflation models generally cannot be reduced to the Eq. (10)). For example, in the (simplest) adiabatic curvaton models, the NG amplitude turns out to be (Bartolo et al. 2004d,c) , 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” evaluated at the epoch of the curvaton decay in the sudden decay approximation. Therefore, for r_{D} ≪ 1, a high level of NG is imprinted.
There exists a clear distinction between multifield and singlefield models of inflation that can be probed via a consistency condition (Maldacena 2003; Creminelli & Zaldarriaga 2004a; Chen et al. 2007b; Chen 2010b): in the squeezed limit, singlefield models predict a bispectrum (11)and thus in the squeezed limit, in a modelindependent sense (i.e., not only for standard singlefield models). This means that a significant detection of local NG (in the squeezed limit) would rule out a very large class of singlefield models of inflation (not just the simplest ones). Although based on very general conditions, the consistency condition of Eq. (11) can be violated in some wellmotivated inflationary settings (we refer the reader to Chen 2010b; Chen et al. 2013 and references therein for more details).
Quasisingle field inflation: quasisingle field inflation has an extra field (or fields) with mass m close to the Hubble parameter H during inflation; these models evolve quiescently, producing a calculable nonGaussian signature (Chen & Wang 2010b). The resulting oneparameter bispectrum smoothly interpolates between local and equilateral models, though in a nontrivial manner: (12)where ν = (9/4 − m^{2}/H^{2})^{1/2} and N_{ν} is the Neumann function of order ν. Quasisingle field models can also produce an essentially “constant” bispectrum defined by . The constant model is the simplest possible nonzero primordial shape, with all its latetime CMB structure simply reflecting the behaviour of the transfer functions.
Alternatives to inflation: local NG can also be generated in some alternative scenarios to inflation, for instance in cyclic/ekpyrotic models (for a review, see Lehners 2010), due to the same basic curvaton mechanism described above. In this case, typical values of the nonlinearity parameter can easily reach .
2.3. Nonstandard models giving rise to alternative specific forms of NG
NonBunchDavies vacuum and higherderivative interactions: another interesting bispectrum shape is the folded one, which peaks in flattened configurations. To facilitate data analyses, the flat shape has been usually parametrized by the template (Meerburg et al. 2009) (13)The initial quantum state of the inflaton is usually specified by requiring that, at asymptotically early times and short distances, its fluctuations behave as in flat space. Deviations from this standard “BunchDavies” vacuum can result in interesting features in the bispectrum. Models with an initial nonBunchDavies vacuum state (Chen et al. 2007b; Holman & Tolley 2008; Meerburg et al. 2009; Ashoorioon & Shiu 2011) can generate sizeable NG similar to this type. NG highly correlated with such a template can be produced in singlefield models of inflation from higherderivative interactions (Bartolo et al. 2010a), and in models where a “Galilean” symmetry is imposed (Creminelli et al. 2011a). In both cases, cubic inflaton interactions with two derivatives of the inflaton field arise. Singlefield inflation models with a small sound speed, studied in Senatore et al. (2010), can generate the flat shape, as a result of a linear combination of the orthogonal and equilateral shapes. In fact, from a simple parametrization point of view, the flat shape can be always written as F_{flat}(k_{1},k_{2},k_{3}) = [F_{equil}(k_{1},k_{2},k_{3}) − F_{ortho}(k_{1},k_{2},k_{3})] /2 (Senatore et al. 2010). Despite this, we provide constraints also on the amplitude of the flat bispectrum shape of Eq. (13).
For models with excited (i.e., nonBunchDavies) initial states, the resulting NG shapes are modeldependent, but they are usually characterized by the importance of flattened or collinear triangles, with k_{3} ≈ k_{1} + k_{2} along the edges of the tetrapyd. We will denote the original flattened bispectrum shape, given in Eqs. (6.2) and (6.3) of Chen et al. (2007b), by ; it is generically much more flattened than the “flat” model of Eq. (13). Although this shape was derived specifically for powerlaw kinflation, it encapsulates several different shapes, with amplitudes which can vary between different phenomenological models. These shapes are also typically oscillatory, being regularized by a cutoff scale k_{c} giving the oscillation period; this cutoff k_{c} ≈ (c_{s}τ_{c})^{1} is determined by the (finite) time τ_{c} in the past when the nonBunchDavies component was initially excited. For excited canonical singlefield inflation, the two leading order shapes can be described (Agullo & Parker 2011) by the ansatz (14)where is dominated by squeezed configurations, has a flattened shape, and i = 1,2. Note that for all oscillatory shapes, the relevant bispectrum equation defines the normalisation of f_{NL}. The flattened signal is most easily enhanced in the limit of small sound speed c_{s}, for which a regularized ansatz is given by (Chen 2010b) (15)
Scaledependent feature and resonant models: oscillating bispectra can be generated from violation of a smooth slowroll evolution (“feature” or “resonant” NG). These models have the distinctive property of a strong running NG, which breaks approximate scaleinvariance. A sharp feature in the inflaton potential forces the inflaton field away from the attractor solution, and causes oscillations as it relaxes back; these oscillations can appear in the bispectrum (Wang & Kamionkowski 2000; Chen et al. 2007a, 2008), as well as the power spectrum and other correlators. An analytic form for the oscillatory bispectrum for these feature models is (Chen et al. 2008) (16)where φ is a phase factor and k_{c} is a scale associated with the feature, which is linked in turn to an effective multipole periodicity ℓ_{c} of the CMB bispectrum. Typically, these oscillations will decay with an envelope of the form exp [ − (k_{1} + k_{2} + k_{3}) /mk_{c}] for a modeldependent parameter m.
Closely related “resonant” bispectra can be created by periodic features superimposed on a smooth inflation potential (Chen et al. 2008; Flauger & Pajer 2011); these induce small periodic features in the background evolution, with which the quantum inflaton fluctuations can resonate while still inside the horizon. Resonant models are particularly relevant in the context of axion inflation models (e.g., Flauger et al. 2010; Flauger & Pajer 2011; Barnaby et al. 2012b). These mechanisms also create oscillatory behaviour in the bispectrum, but with a more constant amplitude and a wavelength that becomes logarithmically stretched. Here, the resonant oscillations for most models can be represented in the form (17)where the constant C = 1/ln(3k_{c}) and φ is a phase.
Finally, we note that periodic features in the inflationary potential can excite the vacuum state, as well as perturbing the background inflation trajectory (Chen 2010a). Such models offer the intriguing possibility of combining the flattened nonBunchDavies shape with periodic oscillations: (18)
This ansatz represents the dominant folded resonant contribution in inflationary models with noncanonical kinetic terms, which competes with resonant (Eq. (17)) and equilateral (Eq. (3)) contributions; however, for slowroll singlefield inflation, there are additional terms.
Directional dependence motivated by gauge fields: additional variations of the bispectrum shape have been proposed for models with vector fields, which can have an additional directional dependence through the parameter where . For example, primordial magnetic fields sourcing curvature perturbations can cause a dependence on both μ and μ^{2} (Shiraishi et al. 2012), and a coupling between the inflaton φ and the gauge field strength F^{2} can yield a μ^{2} dependence (Barnaby et al. 2012a; Bartolo et al. 2013). We can parameterize these shapes as variations on the local shape, following Shiraishi et al. (2013), as (19)where P_{L}(μ) is the Legendre polynomial with P_{0} = 1,P_{1} = μ and . For example, for L = 1 we have the shape (20)Also the recently introduced “solid inflation” model (Endlich et al. 2013) generates bispectra similar to Eq. (19). Here and in the following the nonlinearity parameters are related to the c_{L} coefficients by , , and . The L = 1,2 shapes exhibit sharp variations in the flattened limit for e.g., k_{1} + k_{2} ≈ k_{3}, while in the squeezed limit, L = 1 is suppressed whereas L = 2 grows like the local bispectrum shape (i.e., the L = 0 case). Whether or not the underlying gauge field models prove robust, this directional dependence on the wave vectors is a generic feature which yields distinct bispectrum families, deserving closer study.
Warm inflation: in warm inflation (Berera 1995), where dissipative effects are important, a nonGaussian signal can be generated (e.g., Moss & Xiong 2007) that peaks in the squeezed limit – but with a more complex shape than the local one – and exhibiting a low crosscorrelation with the other shapes (see references in Liguori et al. 2010).
2.4. Higherorder nonGaussianity: the trispectrum
The connected fourpoint functions of CMB anisotropies (or the harmonic counterpart, the socalled trispectrum) can also provide crucial information about the mechanism that gave rise to the primordial curvature perturbations (Okamoto & Hu 2002). The primordial trispectrum is usually characterised by two amplitudes τ_{NL} and g_{NL}: τ_{NL} is most often related to type contributions, while g_{NL} is the amplitude of intrinsic cubic nonlinearities in the primordial gravitational potential (corresponding, in terms of field interactions, to a scalarexchange and to a contact interaction term, respectively). They correspond to “soft” limits of the full fourpoint function, with respectively the diagonal and one side of the general wavevector trapezoid being much smaller than the others. In the CMB maps they appear respectively approximately as a spatial variation in amplitude of the smallscale fluctuations, and a spatial variations in the value of f_{NL} correlated with the largescale temperature. In addition to possible primordial signals that are the focus of this paper there is also expected to be a large lensing trispectrum (of very different shape), discussed in detail in Planck Collaboration XVII (2014).
The simplest local trispectrum is given by
where k_{ij} ≡  k_{i} + k_{j} . Previous constraints on τ_{NL} and g_{NL} have been derived, e.g., by Smidt et al. (2010) who obtained −7.4 × 10^{5}<g_{NL}< 8.2 × 10^{5} and −0.6 × 10^{4}<τ_{NL}< 3.3 × 10^{4} (at 95% CL) analysing WMAP5 data; for the same datasets Fergusson et al. (2010b) obtained −5.4 × 10^{5}<g_{NL}< 8.6 × 10^{5} (68% CL). This kind of trispectrum typically arises in multifield inflationary models where large NG arise from the conversion of isocurvature perturbations on superhorizon scales. If the curvature perturbation is the standard local form, in real space one has . In this case, ; however, in general the trispectrum amplitude can be larger.
The trispectrum is a complementary observable to the CMB bispectrum as it can further distinguish different inflationary scenarios. This is because the same interactions that lead to the bispectrum might be responsible also for a large trispectrum, so that the different NG parameters can be related to each other in a welldefined way within specific models. If there is a nonzero squeezedshape bispectrum there must necessarily be a trispectrum, with (Suyama & Yamaguchi 2008; Suyama et al. 2010; Sugiyama et al. 2011; Sugiyama 2012; Lewis 2011; Smith et al. 2011; Assassi et al. 2012; Kehagias & Riotto 2012). In the simplest inflationary scenarios the prediction would be , but larger values would indicate more complicated dynamics. Several inflationary scenarios have been found in which the bispectrum is suppressed, thus leaving the trispectrum as the largest higherorder correlator in the data. A detection of a large trispectrum and a negligible bispectrum would be a smoking gun for these models. This is the case, for example, of certain curvaton and multifield field models of inflation (Byrnes et al. 2006; Sasaki et al. 2006; Byrnes & Choi 2010), which for particular parameter choice can produce a significant τ_{NL} and g_{NL} and small f_{NL}. Large trispectra are also possible in singlefield models of inflation with higherderivative interactions (see, e.g., Chen et al. 2009; Arroja et al. 2009; Senatore & Zaldarriaga 2011; Bartolo et al. 2010b), but these would be suppressed in the squeezed limit since they are generated by derivative interactions at horizoncrossing, and hence only project weakly onto the local shapes. These equilateral trispectra arise can be welldescribed by some template forms (Fergusson et al. 2010b). Naturally, higherorder correlations could also be considered, but are not directly studied in this paper.
3. Statistical estimation of the CMB bispectrum
In this section, we review the statistical techniques that we use to estimate the nonlinearity parameter f_{NL}. We begin by fixing some notation and describing the CMB angular bispectrum in Sect. 3.1. We then introduce in Sect. 3.2 the optimal f_{NL} bispectrum estimator. From Sect. 3.2.1 onwards we describe in detail the different implementations of the optimal estimator that were developed and applied to Planck data.
3.1. The CMB angular bispectrum
Temperature anisotropies are represented using the a_{ℓm} coefficients of a spherical harmonic decomposition of the CMB map, (22)we write C_{ℓ} = ⟨  a_{ℓm}  ^{2} ⟩ for the angular power spectrum and Ĉ_{ℓ} = (2ℓ + 1)^{1} ∑ _{m}  a_{ℓm}  ^{2} for the corresponding (ideal) estimator; hats “” denote estimated quantities. The CMB angular bispectrum is the threepoint correlator of the a_{ℓm}: (23)If the CMB sky is rotationally invariant, the angular bispectrum can be factorized as follows: (24)where b_{ℓ1ℓ2ℓ3} is the so called reduced bispectrum, and is the Gaunt integral, defined as: where h_{ℓ1ℓ2ℓ3} is a geometrical factor, (27)The Wigner3j symbol in parentheses enforces rotational symmetry, and allows us to restrict attention to a tetrahedral domain of multipole triplets { ℓ_{1},ℓ_{2},ℓ_{3} }, 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). This threedimensional domain of allowed multipoles, sometimes referred to in the following as a “tetrapyd”, is illustrated in Fig. 1 and it is explicitly defined by
(28)Here, is the isotropic subset of the full space of bispectra, denoted by .
One can also define an alternative rotationallyinvariant reduced bispectrum B_{ℓ1ℓ2ℓ3} in the following way: (29)Note that this B_{ℓ1ℓ2ℓ3} is equal to h_{ℓ1ℓ2ℓ3} times the angleaveraged bispectrum as defined in the literature. From Eqs. (24) and (25), and the fact that the sum over all m_{i} of the Wigner3j symbol squared is equal to 1, it is easy to see that B_{ℓ1ℓ2ℓ3} is related to the reduced bispectrum by (30)The interest in this bispectrum B_{ℓ1ℓ2ℓ3} is that it can be estimated directly from maximallyfiltered maps of the data: (31)where the filtered maps are defined as: (32)This can be seen by replacing the in Eq. (29) by its estimate a_{ℓ1m1}a_{ℓ2m2}a_{ℓ3m3} and then using Eq. (25) to rewrite the Wigner symbol in terms of a Gaunt integral, which in its turn is expressed as an integral over the product of three spherical harmonics.
Fig. 1
Permitted observational domain of Eq. (28) for the CMB bispectrum b_{ℓ1ℓ2ℓ3}. Allowed multipole values (ℓ_{1},ℓ_{2},ℓ_{3}) lie inside the shaded “tetrapyd” region (tetrahedron+pyramid), satisfying both the triangle condition and the experimental resolution ℓ <L ≡ ℓ_{max}. 
3.2. CMB bispectrum estimators
The full bispectrum for a highresolution map cannot be evaluated explicitly because of the sheer number of operations involved, , as well as the fact that the signal will be too weak to measure in individual multipoles with any significance. Instead, we essentially use a leastsquares fit to compare the bispectrum of the observed CMB multipoles with a particular theoretical bispectrum b_{ℓ1ℓ2ℓ3}. We then extract an overall “amplitude parameter” f_{NL} for that specific template, after defining a suitable normalization convention so that we can write , where is defined as the value of the theoretical bispectrum ansatz for f_{NL} = 1.
Optimal 3point estimators (introduced by Heavens 1998; see also Gangui & Martin 2000a), are those which saturate the CramérRao bound. Taking into account the fact that instrument noise and masking can break rotational invariance, it has been shown that the general optimal f_{NL} estimator can be written as (Babich 2005; Creminelli et al. 2006; Senatore et al. 2010; Verde et al. 2013): (33)where C^{1} is the inverse of the covariance matrix C_{ℓ1m1,ℓ2m2} ≡ ⟨ a_{ℓ1m1}a_{ℓ2m2} ⟩ and N is a suitable normalization chosen to produce unit response to .
In the expression of the optimal estimator above we note the presence of two contributions, one (hereafter defined the “cubic term” of the estimator) is cubic in the observed a_{ℓm} and correlates the bispectrum of the data to the theoretical fitting template , while the other is linear in the observed a_{ℓm} (hereafter, the “linear term”), which is zero on average. In the rotationallyinvariant case the linear term is proportional to the monopole in the map, which has been set to zero, so in this case the estimator simply reduces to the cubic term. However, when rotational invariance is broken by realistic experimental features such as a Galactic mask or an anisotropic noise distribution, the linear term has an important effect on the estimator variance. In this case, the coupling between different ℓ would in fact produce a spurious increase in the error bars (coupling of Fourier modes due to statistical anisotropy can be “misinterpreted” by the estimator as NG). The linear term correlates the observed a_{ℓm} to the power spectrum anisotropies and removes this effect, thus restoring optimality (Creminelli et al. 2006; Yadav et al. 2008, 2007).
The actual problem with Eq. (33) is that its direct implementation to get an optimal f_{NL} estimator would require measurement of all the bispectrum configurations from the data. As already mentioned at the beginning of this section, the computational cost of this would scale like and be totally prohibitive for highresolution CMB experiments. Even taking into account the constraints imposed by isotropy, the number of multipole triples { ℓ_{1},ℓ_{2},ℓ_{3} } is of the order of 10^{9} at Planck resolution, and the number of different observed bispectrum configurations is of the order of 10^{15}. For each of them, costly numerical evaluation of the Wigner symbol is also required. This is completely out of reach of existing supercomputers. It is then necessary to find numerical solutions that circumvent this problem and in the following subsections we will show how the different estimators used for the f_{NL}Planck data analysis address this challenge. Before entering into a more accurate description of these different methods, we would like however to stress again that they are all going to be different implementations of the optimal f_{NL} estimator defined by Eq. (33); therefore they are conceptually equivalent and expected to produce f_{NL} results that are in very tight agreement. This will later on allow for stringent validation tests based on comparing different pipelines. On the other hand, it will soon become clear that the different approaches that we are going to discuss also open up a range of additional applications beyond simple f_{NL} estimation for standard bispectra. Such applications include, for example, full bispectrum reconstruction (in a suitably smoothed domain), tests of directional dependence of f_{NL}, and other ways to reduce the amount of data, going beyond simple singlenumber f_{NL} estimation. So different methods will also provide a vast range of complementary information.
Another important preliminary point, to notice before discussing different techniques, is that none of the estimators in the following sections implement exactly Eq. (33), but a slightly modified version of it. In Eq. (33) the CMB multipoles always appear weighted by the inverse of the full covariance matrix. Inverse covariance filtering of CMB data at the high angular resolutions achieved by experiments like WMAP and Planck is another very challenging numerical issue, which was fully addressed only recently (Smith et al. 2009; Komatsu et al. 2011; Elsner & Wandelt 2013). For our analyses we developed two independent inversecovariance filtering pipelines. The former is based on an extension to Planck resolution of the algorithm used for WMAP analysis (Smith et al. 2009; Komatsu et al. 2011); the latter is based on the algorithm described in Elsner & Wandelt (2013). However, detailed comparisons interestingly showed that our estimators perform equally well (i.e., they saturate the CramérRao bound) if we approximate the covariance matrix as diagonal in the filtering procedure and we apply a simple diffusive inpainting procedure to the masked areas of the input CMB maps. A more detailed description of our inpainting and Wiener filtering algorithms can be found in Sect. 3.3.
In the diagonal covariance approximation, the minimum variance estimator is obtained by making the replacement (C^{1}a)_{ℓm} → a_{ℓm}/C_{ℓ} in the cubic term and then including the linear term that minimizes the variance for this class of cubic estimator (Creminelli et al. 2006). This procedure leads to the following expression: (34)where the tilde denotes the modification of C_{ℓ} and b_{ℓ1ℓ2ℓ3} to incorporate instrument beam and noise effects, and indicates that the multipoles are obtained from a map that was masked and preprocessed through the inpainting procedure detailed in Sect. 3.3. This means that (35)where b_{ℓ} denotes the experimental beam, and N_{ℓ} is the noise power spectrum. For simplicity of notation, in the following we will drop the tilde and always assume that beam, noise and inpainting effects are properly included.
Using Eqs. (29) and (30) we can rewrite Eq. (34) in terms of the bispectrum B_{ℓ1ℓ2ℓ3}: (36)In the above expression, B^{th} is the theoretical template for B (with f_{NL} = 1) and B^{obs} denotes the observed bispectrum (the cubic term), extracted from the (inpainted) data using Eq. (31). B^{lin} is the linear correction, also computed using Eq. (31) by replacing two of the filtered temperature maps by simulated Gaussian ones and averaging over a large number of them (three permutations). The variance V in the inversevariance weights is given by (remember that these should be viewed as being the quantities with tildes, having beam and noise effects included) with g_{ℓ1ℓ2ℓ3} a permutation factor (g_{ℓ1ℓ2ℓ3} = 6 when all ℓ are equal, g_{ℓ1ℓ2ℓ3} = 2 when two ℓ are equal, and g_{ℓ1ℓ2ℓ3} = 1 otherwise). Both Eqs. (34) and (36) will be used in the following. Equation (34) will provide the starting point for the KSW, skewC_{ℓ} and modal estimators, while Eq. (36) will be the basis for the binned and wavelets estimators.
Next, we will describe in detail the different methods, and show how they address the numerical challenge posed by the necessity to evaluate a huge number of bispectrum configurations. To summarize loosely: the KSW estimator, the skewC_{ℓ} approach and the separable modal methodology achieve massive reductions in computational costs by exploiting structural properties of b^{th}, e.g., separability. On the other hand, the binned bispectrum and the wavelet approaches achieve computational gains by data compression of B^{obs}.
3.2.1. The KSW estimator
To understand the rationale behind the KSW estimator (Komatsu et al. 2005, 2003; Senatore et al. 2010; Creminelli et al. 2006; Yadav et al. 2008, 2007; Yadav & Wandelt 2008; Smith & Zaldarriaga 2011), assume that the theoretical reduced bispectrum can be exactly decomposed into a separable structure, e.g., there exist some sequences of functions α(ℓ,r),β(ℓ,r) such that we can approximate b_{ℓ1ℓ2ℓ3} as (37)where r is a radial coordinate. This assumption is fulfilled in particular by the local shape (Komatsu & Spergel 2001; Babich et al. 2004), with α(ℓ,r) and β(ℓ,r) involving integrals of products of spherical Bessel functions and CMB radiation transfer functions. Let us consider the optimal estimator of Eq. (33) and neglect for the moment the linear part. Exploiting Eq. (37) and the factorizability property of the Gaunt integral (Eq. (25)), the cubic term of the estimator can be written as: (38)where (39)and (40)From the formulae above we see that the overall triple integral over all the configurations ℓ_{1},ℓ_{2},ℓ_{3} has been factorized into a product of three separate sums over different ℓ. This produces a massive reduction in computational time, as the problem now scales like instead of the original . Moreover, the bispectrum can be evaluated in terms of a cubic statistic in pixel space from Eq. (38), and the functions , are obtained from the observed a_{ℓm} by means of Fast Harmonic Transforms.
It is easy to see that the linear term can be factorized in analogous fashion. Again considering the local shape type of decomposition of Eq. (37), it is possible to find: (41)where ⟨·⟩_{MC} denotes a Monte Carlo (MC) average over simulations accurately reproducing the properties of the actual data set (basically we are taking a MC approach to estimate the product between the theoretical bispectrum and the a_{ℓm} covariance matrix appearing in the linear term expression).
The estimator can be finally expressed as a function of S_{cub} and S_{lin}: (42)Whenever it can be applied, the KSW approach makes the problem of f_{NL} estimation computationally feasible, even at the high angular resolution of the Planck satellite. One important caveat is that factorizability of the shape, which is the starting point of the method, is not a general property of theoretical bispectrum templates. Strictly speaking, only the local shape is manifestly separable. However, a large class of inflationary models can be extremely well approximated by separable equilateral and orthogonal templates (Babich et al. 2004; Creminelli et al. 2006; Senatore et al. 2010). The specific expressions of cubic and linear terms are of course templatedependent, but as long as the template itself is separable their structure is analogous to the example shown in this section, i.e., they can be written as pixel space integrals of cubic products of suitablyfiltered CMB maps (involving MC approximations of the a_{ℓm} covariance for the linear term). For a complete and compact summary of KSW implementations for local, equilateral and orthogonal bispectra see Komatsu et al. (2009, Appendix).
3.2.2. The skewC_{ℓ} extension
The skewC_{ℓ} statistics were introduced by Munshi & Heavens (2010) to address an issue with estimators such as KSW which reduce the map to an estimator of f_{NL} for a given type of NG. This level of data compression, to a single number, has the disadvantage that it does not allow verification that a NG signal is of the type which has been estimated. KSW on its own cannot tell if a measurement of f_{NL} of given type is actually caused by NG of that type, or by contamination from some other source or sources. The skewC_{ℓ} statistics perform a less radical data compression than KSW (to a function of ℓ), and thus retain enough information to distinguish different NG signals. The desire to find a statistic which is able to fulfil this rôle, but which is still optimal, drives one to a statistic which is closely related to KSW, and indeed reduces to it when the scaledependent information is not used. A further advantage of the skewC_{ℓ} is that it allows joint estimation of the level of many types of NG simultaneously. This requires a large number of simulations for accurate estimation of their covariance matrix, and they are not used in this role in this paper. However, they do play an important part in identifying which sources of NG are clearly detected in the data, and which are not.
We define the skewC_{ℓ} statistics by extending from KSW, as follows: from Eq. (38), the numerator ℰ can be rewritten as (43)where (44)and (45)The skewC_{ℓ} approach allows for the full implementation of the KSW procedure, when the sum in Eq. (43) is fully evaluated; furthermore, it allows for extra degrees of flexibility, e.g., by restricting the sum to subsets of the multipole space, which may highlight specific features of the NG signal. Each form of NG considered has its own α,β, hence its own set of skewC_{ℓ}, denoted , and we have chosen to illustrate here with the local form, but as with KSW the method can be extended to other separable shapes, and some skewC_{ℓ} do not involve integrals, such as the ISWlensing skew statistic. Note that in this paper we do not fit the S_{ℓ} directly, but instead we estimate the NG using KSW, and then verify (or not) the nature of the NG by comparing the skewC_{ℓ} with the theoretical expectation. No further free parameters are introduced at this stage. This procedure allows investigation of KSW detections of NG of a given type, assessing whether or not they are actually due to NG of that type.
3.2.3. Separable modal methodology
Primordial bispectra need not be manifestly separable (like the local bispectrum), or be easily approximated by separable ad hoc templates (equilateral and orthogonal), so the direct KSW approach above cannot be applied generically (nor to latetime bispectra). However, we can employ a highlyefficient generalization by considering a complete basis of separable modes describing any latetime bispectrum (see Fergusson & Shellard 2007; Fergusson et al. 2010a), as applied to WMAP7 data for a wide variety of separable and nonseparable bispectrum models (Fergusson et al. 2012). See also Planck Collaboration XXIII (2014) and Planck Collaboration XXV (2014). We can achieve this by expanding the signaltonoiseweighted bispectrum as (46)where the (nonorthogonal) separable modes Q_{n} are defined by (47)It is more efficient to label the basis as Q_{n}, with the subscript n representing an ordering of the { i,j,k } products (e.g., by distance i^{2} + j^{2} + k^{2}). The are any complete basis functions up to a given resolution of interest and they can be augmented with other special functions adapted to target particular bispectra. The modal coefficients in the bispectrum of Eq. (46) are given by the inner product of the weighted bispectrum with each mode (48)where the modal transformation matrix is (49)In the following, the specific basis functions we employ include either weighted Legendrelike polynomials or trigonometric functions. These are combined with the SachsWolfe local shape and the separable ISW shape in order to obtain high correlations to all known bispectrum shapes (usually in excess of 99%).
Substituting the separable mode expansion of Eq. (46) into the estimator and exploiting the separability of the Gaunt integral (Eq. (25)), yields (50)where the represent versions of the original CMB map filtered with the basis function q_{p} (and the weights (), that is, (51)The maps incorporate the same mask and a realistic model of the inhomogeneous instrument noise; a large ensemble of these maps, calculated from Gaussian simulations, is used in the averaged linear term in the estimator of Eq. (50), allowing for the subtraction of these important effects. Defining the integral over these convolved product maps as cubic and linear terms respectively, the estimator reduces to a simple sum over the mode coefficients (54)where . The estimator sum in Eq. (54) is now straightforward to evaluate because of separability, since it has been reduced to a product of three sums over the observational maps (Eq. (50)), followed by a single 2D integral over all directions (Eq. (52)). The number of operations in evaluating the estimator sum is only .
For the purposes of testing a wide range inflationary models, we can also define a set of primordial basis functions for wavenumbers satisfying the triangle condition (again we will order the { i,j,k } with n). This provides a separable expansion for an arbitrary primordial shape function S(k_{1},k_{2},k_{3}) = B(k_{1},k_{2},k_{3})/(k_{1}k_{2}k_{3})^{2}, that is, (55)Using the same transfer functions as in the KSW integral (38), we can efficiently project forward each separable primordial mode to a corresponding latetime solution (essentially the projected CMB bispectrum for ). By finding the inner product between these projected modes and the CMB basis functions Q_{n}(ℓ_{1},ℓ_{2},ℓ_{3}), we can obtain the transformation matrix (Fergusson et al. 2010a,b) (56)which projects the primordial expansion coefficients to latetime: (57)When Γ_{np} is calculated once we can efficiently convert any given primordial bispectrum B(k_{1},k_{2},k_{3}) into its latetime CMB bispectrum counterpart using Eq. (46). We can use this to extend the KSW methodology and to search for the much broader range of primordial models beyond local, equilateral and orthogonal, having validated on these standard shapes.
3.2.4. Binned bispectrum
In the binned bispectrum approach (Bucher et al. 2010), the computational gains are achieved by data compression of the observed . This is quite feasible, because like the power spectrum the bispectrum is a rather smooth function, with features on the scale of the acoustic peaks. In this way one obtains an enormous reduction of the computational resources needed at the cost of only a tiny increase in variance (typically 1%).
More precisely, the following statistic is considered, (58)where the Δ_{i} are intervals (bins) of multipole values [ℓ_{i},ℓ_{i + 1} − 1], for i = 0,...,(N_{bins} − 1), with ℓ_{0} = ℓ_{min} and ℓ_{Nbins} = ℓ_{max} + 1, and the other bin boundaries chosen in such a way as to minimize the variance of . The binned bispectrum is then obtained by using T_{i} instead of T_{ℓ} in the expression for the sample bispectrum of Eq. (31), to obtain: (59)The linear term B^{lin} is binned in an analogous way, and the theoretical bispectrum template B^{th} and variance V are also binned by summing them over the values of ℓ inside the bin. Finally f_{NL} is determined using the binned version of Eq. (36), i.e., by replacing all quantities by their binned equivalent and replacing the sum over ℓ by a sum over bin indices i. An important point is that the binned bispectrum estimator does not mix the observed bispectrum and the theoretical template weights until the very last step of the computation of (the final sum over bin indices). Thus, the (binned) bispectrum of the map is also a direct output of the code. Moreover, one can easily study the ℓdependence of the results by omitting bins from this final sum.
The full binned bispectrum allows one to explore the bispectral properties of maps independent of a theoretical model. Binned bispectra have been used to compare component separation maps and singlefrequency maps dominated by foregrounds, as presented in Sects. 3.4.2 and 7. In its simplest implementation, which is used in this paper, the binned estimator uses tophat filters in harmonic space, which makes the Gaussian noise independent between different bins; however, slightly overlapping bins could be used to provide better localization properties in pixel space. In this sense the binned estimator is related to the wavelet estimators, which we discuss below.
3.2.5. Wavelet f_{NL} estimator
Wavelet methods are wellestablished in the CMB literature and have been applied to virtually all areas of the data analysis pipeline, including mapmaking and component separation, point source detection, search for anomalies and anisotropies, crosscorrelation with large scale structure and the ISW effect, and many others (see for instance Antoine & Vandergheynst 1998; MartínezGonzález et al. 2002; Cayon et al. 2003; McEwen et al. 2007a,b; Pietrobon et al. 2006; Starck et al. 2006; Cruz et al. 2007; Faÿ et al. 2008; Feeney et al. 2011; Geller & Mayeli 2009a,b; Starck et al. 2010; Scodeller et al. 2011; FernándezCobos et al. 2012). These methods have the advantage of possessing localization properties both in real and harmonic space, allowing the effects of masked regions and anisotropic noise to be dealt with efficiently.
In terms of the current discussion, wavelets can be viewed as a way to compress the sample bispectrum vector by a careful binning scheme in the harmonic domain. See also Planck Collaboration XXIII (2014). In particular, the wavelet bispectrum can be rewritten as (60)where (61)Here b is the position on the sky at which the wavelet coefficient is evaluated and σ is is the dispersion of the wavelet coefficient map w(R;b) for the scale R. The filters ω_{ℓ}(R) can be seen as the coefficients of the expansion into spherical harmonic of the Spherical Mexican Hat Wavelet (SMHW) filter (62)Here is a normalizing constant and y = 2tan(θ/ 2) represents the distance between x and n, evaluated on the stereographic projection on the tangent plane at n; θ is the corresponding angular distance, evaluated on the spherical surface.
The implementation of the linearterm correction can proceed in analogy with the earlier cases. However, note that, in view of the realspace localization properties of the wavelet filters, the linear term here is smaller than for KSW and related approaches, although not negligible. Moreover, it can be wellapproximated by a termbyterm samplemean subtraction for the wavelet coefficients, which allows for a further reduction of computational costs. Further details can be found in Curto et al. (2011a,b, 2012); Regan et al. (2013; see also Lan & Marinucci 2008; Rudjord et al. 2009; Pietrobon et al. 2009, 2010; Donzelli et al. 2012, for related needletbased procedures).
3.3. Wiener filtering
As discussed in Sect. 3.2, the f_{NL} bispectrum estimator requires, in principle, inverse covariance filtering of the data to achieve optimality (equivalent to Wiener filtering up to a trivial multiplication by the inverse of the signal power spectrum).
We have used the iterative method of Elsner & Wandelt (2013) for Wiener filtering simulations and data. The algorithm makes use of a messenger field, introduced to mediate between the pixel space and harmonic space representation, where noise and signal properties can be specified most directly. Since the Wiener filter is the maximum a posteriori solution, we monitor the χ^{2} of the current solution as a convergence diagnostics. We terminate the algorithm as soon as the improvement in the posterior between two consecutive steps has dropped below a threshold of Δχ^{2} ≤ 10^{4}σ_{χ2}, where σ_{χ2} is the standard deviation of χ^{2}distributed variables with a number of degrees of freedom given by the number of observed pixels. Results of this method have been crossvalidated using an independent conjugate gradient inversion algorithm with multigrid preconditioning, originally developed for the analysis of WMAP data in Smith et al. (2009). Applying this estimator to simulations preprocessed using the above mentioned algorithms yielded optimal error bars as expected.
However, we found that optimal error bars could also be achieved for all shapes using a much simpler diffusive inpainting prefiltering procedure that can be described as follows: all masked areas of the sky (both Galactic and point sources) are filled in with an iterative scheme. Each pixel in the mask is filled with the average of all surrounding pixels, and this is repeated 2000 times over all masked pixels (we checked on simulations that convergence of all f_{NL} estimates was achieved with 2000 iterations). Note that the effect of this inpainting procedure, especially visible for the Galactic mask, is effectively to apodize the mask, reproducing smallscale structure near the edges and only largescale modes in the interior. This helps to prevent propagating any sharpedge effects or lack of largescale power in the interior of the mask to the unmasked regions during harmonic transforms.
Any bias and/or excess variance arising from the inpainting procedure were assessed through MC validation (see Sect. 6) and found to be negligible. Since the inpainting procedure is particularly simple to implement, easily allows inclusion of realistic correlatednoise models in the simulations, and retains optimality, we chose inpainting as our data filtering procedure for the f_{NL} analysis.
3.4. CMB bispectrum reconstruction
3.4.1. Modal bispectrum reconstruction
Modal and related estimators can be used to reconstruct the full bispectrum from the modal coefficients β_{ijk} obtained from map filtering with the separable basis functions Q_{ijk}(ℓ_{1},ℓ_{2},ℓ_{3}) (Eq. (52)) (Fergusson et al. 2010a). It is easy to show that the expectation value for β_{ijk} (or equivalently β_{n}) for an ensemble of nonGaussian maps generated with a CMB bispectrum shape given by α_{n} (Eq. (48)) (with amplitude f_{NL}) is (63)where γ_{np} is defined in Eq. (49). Hence, we expect the best fit coefficients for a particular α_{n} realization to be given by the β_{n}s themselves (for a sufficiently large signal). Assuming that we can extract the β_{n} coefficients accurately from a particular experiment, we can directly reconstruct the CMB bispectrum using the expansion of Eq. (47), that is, where, for convenience, we have also defined orthonormalized basis functions R_{n}(ℓ_{1},ℓ_{2},ℓ_{3}) with coefficients and such that ⟨ R_{n},R_{p} ⟩ = δ_{np}. This method has been validated for simulated maps, showing the accurate recovery of CMB bispectra, and it has been applied to the WMAP7 data to reconstruct the full 3D CMB bispectrum (Fergusson et al. 2012).
To quantify whether or not there is a modelindependent deviation from Gaussianity, we can consider the total integrated bispectrum. By summing over all multipoles, we can define a total integrated nonlinearity parameter which, with the orthonormal modal decomposition of Eq. (65) becomes (Fergusson et al. 2012) (66)where N_{loc} is the normalization for the local f_{NL} = 1 model. The expectation value contains more than just the threepoint correlator, with contributions from products of the twopoint correlators and even higherorder contributions, (67)Here is the full 6point function over the unconnected Gaussian part, i.e., the product of 3C_{ℓ}. So, for a Gaussian input this recovers an average of 1 per mode. In the nonGaussian case the leadingorder contributions after the Gaussian part are the bispectrum squared and the product of the power spectrum and the trispectrum, which enter at the same order (for additional explanations see Fergusson et al. 2012).
3.4.2. Smoothed binned bispectrum reconstruction
As explained in Sect. 3.2.4, the full binned bispectrum of the maps under study is one of the products of the binned estimator code.
Given the relatively fine binning (about 50 bins up to ℓ = 2000 or 2500), most of the measurement in any single bintriplet is noise. If combinations of maps are chosen so that the CMB primordial signal dominates, most of this noise is Gaussian, reflecting the fact that even when ⟨ xyz ⟩ = 0, for a particular statistical realization, xyz is almost certainly nonzero. If our goal is to test whether there is any statistically significant signal, then it makes sense to normalize by defining a new field ℬ_{i1i2i3}, which is the binned bispectrum divided by its expected standard deviation (computed in the standard way assuming Gaussian statistics). The distribution of ℬ_{i1i2i3} in any one bintriplet is very nearly Gaussian as a result of the central limit theorem, and there is almost no correlation between different bins.
We could study the significance of the extreme values at this fine resolution, but it also makes sense to smooth in order to detect features coherent over a wider range of ℓ. If the threedimensional domain over which the bispectrum is defined (consisting of those triplets in the range [ℓ_{min},ℓ_{max}] satisfying the triangle inequality and parity condition) did not have boundaries, the smoothing would be more straightforward. We smooth with a Gaussian kernel of varying width in Δ(lnℓ), normalized so that the smoothed function in each pixel would be normaldistributed if the input map were Gaussian. In this way, based on the extreme values, it is possible to decide on whether there is a statistically significant NG signal in the map in a “blind” (nonparametric) way.
4. Statistical estimation of the CMB trispectrum
4.1. The squeezeddiagonal trispectrum: τ_{NL}
The 4point function, or equivalently the trispectrum, of the CMB, can also place interesting constraints on inflationary physics. There are several physically interesting “shapes” of the trispectrum (e.g., Huang & Shiu 2006; Byrnes et al. 2006; Fergusson et al. 2010b; Izumi et al. 2012), in analogy to the bispectrum case. In the simplest nonGaussian models, the CMB bispectrum has larger signaltonoise than the trispectrum, but there are examples of technically natural models in which the trispectrum has larger signaltonoise (e.g., Senatore & Zaldarriaga 2011; Baumann & Green 2012; see also Bartolo et al. 2010b). This can happen in models in which the field modulating the fluctuation amplitude is only weakly correlated to the observed largescale curvature perturbation.
Analysis of the trispectrum is more challenging than that of the bispectrum, due to the increased range of systematic effects and secondary signals which can contribute. For example, gravitational lensing of the CMB generates a manysigma contribution to the trispectrum, though it has a distinctive anisotropic shape that differs from primordial NG modulated by scalar fields. As an instrumental example, any mismatch between the true covariance of the observed CMB plus noise and the covariance which is assumed in the analysis (due, for example, to mischaracterisation of the pointing, beams, or noise properties) will generally lead to biases in the estimated trispectrum. Due to these challenges, we have deferred a full analysis of the primordial trispectrum to a future paper, and here focus on the simplest squeezed shape that can provide useful constraints on primordial models, τ_{NL}.
τ_{NL}is most easily understood as measuring the largescale modulation of smallscale power. The constraints on f_{NL} show that such a modulation must be small if correlated with the temperature. However it is possible for multifield inflation models to produce squeezedshape modulations which are uncorrelated with the largescale curvature perturbations. Such models can be constrained by the trispectrum, conventionally parameterized by τ_{NL} in the squeezeddiagonal shape.
For example, consider the case where a smallscale Gaussian curvature perturbation ζ_{0} is modulated by another field φ so that the primordial perturbation is given by (68)where φ(x) is a largescale modulating field (with amplitude ≪1). The largescale modes of φ can be measured from the modulation they induce in the smallscale ζ power spectrum. If φ has a nearly scaleinvariant spectrum, the nearlywhite cosmic variance noise on the reconstruction dominates on smallscales, so only the very largest modes can be reconstructed (Kogo & Komatsu 2006). A reconstruction of φ is going to be limited to only very largescale variations, in which case the scale of the variation is very large compared to the width of the lastscattering surface; i.e., in any particular direction a largescale modulating field will modulate all smallscale perturbations through the lastscattering surface by approximately the same amount. This approximation is good at the percent level, and can readily be related to the full trispectrum estimator, as discussed in more detail in Pearson et al. (2012, Sect. IV; see also Okamoto & Hu 2002; Munshi et al. 2011; Smidt et al. 2010).
A largescale power modulation therefore translates directly into a largescale modulation of the smallscale CMB temperature: (69)where T_{g} are the usual smallscale Gaussian CMB temperature anisotropies and r_{∗} is the radial distance to the lastscattering surface. We can quantify the trispectrum as a function of modulation scale by using the power spectrum of the modulation, (70)As is conventional, we normalize relative to , the power spectrum of the primordial curvature perturbation at the location of the recombination surface. The field f is directly observable, but is not, since the curvature perturbation can only be constrained very indirectly on very large scales. We shall therefore give constraints on f, which is directly constrained by Planck, but also on τ_{NL} for comparison with the inflation literature. Note that τ_{NL} ~ 500 corresponds to an modulation.
A general quadratic estimator methodology for reconstructing f was developed in Hanson & Lewis (2009), which we broadly follow here. The structure is essentially identical to that for lensing reconstruction (Planck Collaboration XVII 2014), where here instead of reconstructing a lensing potential (or deflection angle), we are reconstructing a scalar modulation field. The quadratic maximum likelihood estimator for the largescale modulation field f (assuming it is small) is given by (71)where is a quadratic function of the filtered data that can be calculated quickly in real space: (72)Here is an inversevariance filtering sky map (which accounts for sky cuts and inhomogeneous noise), and is the lensed C_{ℓ}. The “mean field” can be estimated from simulations, along with the Fisher normalization ℱ that is given by the covariance of . The i,j indices are included here, since we shall be using different sky maps with independent noise to avoid noise biases at the level of modulation field reconstruction. For low L and high ℓ_{max} the reconstruction noise is very nearly constant (white, because each small patch of sky gives a nearlyuncorrelated but noisy estimate of the smallscale power), and the reconstruction is very local.
In practice the inversevariance filtering is imperfect, the noise cannot be modelled exactly, and the normalizing Fisher matrix ℱ_{ℓmℓ′m′} evaluated from simulations would be inaccurate. Instead we focus on directly, which is approximately an inversevariance weighed reconstruction of the modulation, and is manifestly very local in real space (and hence zero in the cut part of the sky). Since the reconstruction noise, which also approximately determines the normalization, is nearly white (constant in L), in real space has an expectation nearly proportional to the underlying modulation outside the mask.
We then define an estimator of the modulation power spectrum (73)where (74)is a noise bias for zero signal estimated from simulations. The normalization A_{L} is the analytic ideal fullsky normalization which is very close to a constant, and k_{L} is a calibration factor determined from withsignal simulations. On small scales k_{L} ∝ f_{sky}^{1}, but has some scale dependence: it increases towards k_{L} ∝ f_{sky}^{2} at L = 0 at very low L. We shall sometimes plot , corresponding to the uncalibrated reconstruction of the power modulation, which is very local in real space.
For each value of the modulation scale L, Eq. (73) defines a separate estimator for τ_{NL}(75)We can combine estimators from all ℓ by constructing (76)where and cov_{LL′} is the covariance of from simulations with τ_{NL} = 0. On the fullsky the estimators from each L would be independent, but the mask introduces significant coupling between the very low multipoles and this form of the estimator allows us to account for this. In the fullsky uncorrelated approximation, with a nearly scaleinvariant primordial spectrum and using the whiteness of the reconstruction noise, the estimator for τ_{NL} simplifies to (Pearson et al. 2012) (77)where N ≡ L_{min}^{2} − (L_{max} + 1)^{2}. This result does not require many simulations to estimate the covariance accurately for inversion, and is typically expected to give very similar results with an error bar that is less than 10% larger. We calculate both as a crosscheck, but report results for because it is more robust, and in our simulation results actually has slightly lower tails (though larger variance). Mean fields and the bias are estimated in all cases from 1000 zeroτ_{NL} simulations, and the mask used retains about 70% of the sky.
If there is a nearly scaleinvariant signal, so as expected in most multifield inflation models, the contributions fall rapidly ∝1 /L^{3}, as expected when measuring a scaleinvariant signal that has large white reconstruction noise. The signal is therefore on very large scales, with typically half the Fisher signal in the dipole modulation and 95% of the signal at L ≤ 4, justifying the squeezed approximations used. We use L_{max} = 10 for the estimators, which includes almost all of the signaltonoise but avoids excessive contamination with the “blue” spectrum of lensing contributions. However, due to the small number of modes involved, the posterior distributions of τ_{NL} can have quite broad tails corresponding to the finite probability that all the largestscale modulation modes just happen to be near zero. To improve constraints on large values it can help to include a larger range of L, and we consider L up to L_{max} = 50, which is about the limit of where the approximations are valid.
5. Nonprimordial contributions to the CMB bispectrum and trispectrum
In this subsection we present the steps followed to account for and remove the main nonprimordial contributions to CMB NG.
5.1. Foreground subtraction
Foreground emission signals in the microwave bands have a strong nonGaussian signature. Therefore any residual emission in the CMB data can give a spurious apparent primordial NG detection. In our analysis we considered Planck CMB foregroundcleaned maps created using several independent techniques, as described in Planck Collaboration XII (2014): explicit parametrization and fitting of foregrounds in real space (CommanderRuler, CR) (Eriksen et al. 2006, 2008); Spectral Matching of foregrounds implementing Independent Component Analysis (SMICA) (Delabrouille et al. 2003; Cardoso et al. 2008); Internal Linear Combination (Needlet Internal Linear Combination, NILC) (Delabrouille et al. 2009); and Internal Template Fitting (SEVEM) (FernándezCobos et al. 2012). These and other techniques underwent a prelaunch testing phase (Leach et al. 2008). Each method provides a Planck CMB foregroundcleaned map with a confidence mask, which defines the trusted cleaned region of the sky; an estimate of the noise in the output CMB map obtained from halfring difference maps; and an estimate of the beam transfer function of the processed map. The resolution reaches 5 arcminutes. In addition a union of all the confidence masks, denoted as U73, is provided. Channels from both the Low Frequency Instrument (LFI, Planck Collaboration II 2014) and the High Frequency Instrument (HFI, Planck Collaboration VI 2014) of Planck are used to achieve each of the reconstructed CMB templates. The validation of CMB reconstruction through component separation is based on the inspection of several observables, as explained in detail in Planck Collaboration XII (2014): the twopoint correlation function and derived cosmological parameters; indicators of NG including the f_{NL} results presented in the present paper; and crosscorrelation with known foreground templates. Based on various figuresofmerit, the foreground cleaning techniques performed comparably well (Planck Collaboration XII 2014).
In order to test foreground residuals a battery of simulations is required. In the simulations the foreground emission was modelled with the prelaunch version of the Planck Sky Model (PSM), based on observations of the emission from our own Galaxy and known extraGalactic sources, largely in the radio and infrared bands. The PSM is described in Delabrouille et al. (2013) and includes models of CMB (including a dipole), diffuse Galactic emissions (synchrotron, freefree, thermal dust, Anomalous Microwave Emission and CO molecular lines), emission from compact objects (thermal SZ effect, kinetic SZ effect, radio sources, infrared sources, correlated farinfrared background and ultracompact H ii regions). The sky model includes total intensity as well as polarization, which was not used in this paper.
The PSM has been used to create the sixth round of full focal plane (FFP6) simulations, a set of simulations for the 2013 data release based on detailed models of the sky and instrument (e.g., noise properties, beams, satellite pointing and mapmaking process), consisting of both Gaussian and nonGaussian CMB realizations. A description of the FFP6 simulations can be found in Appendix A (see Planck Collaboration 2013, for details). The FFP6 set has been used to test and validate the component separation algorithms employed in Planck and to establish uncertainties on the outputs (Planck Collaboration XII 2014). We will also use FFP6 simulations in the next sections in order to validate our analysis (see Sects. 6–8).
5.2. The integrated SachsWolfelensing bispectrum
One of the most relevant mechanisms that can generate NG from secondary CMB anisotropies is the coupling between weak lensing and the ISW (Sachs & Wolfe 1967) effect. This is in fact the leading contribution to the CMB secondary bispectrum with a blackbody frequency dependence (Goldberg & Spergel 1999; Verde & Spergel 2002; Giovi et al. 2005).
Weak lensing of the CMB is caused by gradients in the matter gravitational potential that distorts the CMB photon geodesics. The ISW on the other hand arise because of timevarying gravitational potentials due to the linear and nonlinear growth of structure in the evolving Universe. Both the lensing and the ISW effect are then related to the matter gravitational potential and thus are correlated phenomena. This gives rise to a nonvanishing 3point correlation function. Furthermore, lensing is related to nonlinear processes which are therefore nonGaussian. A detailed description of the signal, which accounts also for the contribution from the earlyISW effect, can be found in Lewis (2012).
The ISWlensing bispectrum takes the form: (78)where P, L, and ISW indicate primordial, lensing and ISW contributions respectively. This becomes (79)where is the Gaunt integral and is the reduced bispectrum given by (80)Here is the lensed CMB power spectrum and is the ISWlensing crosspower spectrum (Lewis 2012; Goldberg & Spergel 1999; Verde & Spergel 2002; Cooray & Hu 2000) that expresses the statistical expectation of the correlation between the lensing and the ISW effect.
As shown in Hanson et al. (2009b), Mangilli & Verde (2009), and Lewis et al. (2011), the ISWlensing bispectrum can introduce a contamination in the constraints on primordial local NG from the CMB bispectrum. Both bispectra are maximal for squeezed or nearly squeezed configurations. The bias on a primordial f_{NL} (e.g., local) due to the presence of the ISWlensing cross correlation signal is (81)with (82)where B^{ISW − L} and B^{P} refer respectively to the ISWlensing and the primordial bispectrum, and V is defined below Eq. (36).
The bias in the three primordial f_{NL} parameters due to the ISWlensing signal for the four componentseparation methods.
Results for the amplitude of the ISWlensing bispectrum from the SMICA, NILC, SEVEM, and CR foregroundcleaned maps, for the KSW, binned, and modal (polynomial) estimators; error bars are 68% CL.
The bias in the estimation of the three primordial f_{NL} from Planck is given in Table 1. As one can see, taking into account the f_{NL} statistical error bars shown, e.g., in Table 8, the local shape is most affected by this bias (at the level of more than 1σ_{local}), followed by the orthogonal shape (at the level of about 0.5σ_{ortho}), while the equilateral shape is hardly affected. In this paper we have taken into account the bias reported in Table 1 by subtracting it from the measured f_{NL}^{8}.
The results for the amplitude of the ISWlensing bispectrum from the different foregroundcleaned maps are given in Table 2. It should be noted that the binned and modal estimators are less correlated to the exact template for the ISWlensing shape than they are for the primordial shapes, hence their larger error bars compared to KSW (which uses the exact template by construction Mangilli et al. 2013). The conclusion is that we detect the ISWlensing bispectrum at a value consistent with the fiducial value of 1, at a significance level of 2.6σ (taking the SMICAKSW value as reference). For details about comparisons between different estimators and analysis of the data regarding primordial shapes we refer the reader to Sects. 6 and 7.
Fig. 2
Binned skewC_{ℓ} statistics from the SMICA map for a) ISWlensing and b) Poisson point sources. Theoretical curves are not fitted to the data shown, but are plotted with the amplitude (the only free parameter) determined from the KSW technique. The Poisson pointsource foreground is clearly detected, and the ISWlensing skewspectrum is evident and consistent with the overall 2.6σ detection. b_{ps} is the Poisson pointsource amplitude in dimensionless units of 10^{29}, and is the ISWlensing amplitude in units of that expected from the Planck bestfit cosmology. Error bars come from covariance estimates from 1000 simulated maps, and the points are mildly correlated. 
We show for the SMICA map in the top figure of Fig. 2 the measured skewC_{ℓ} spectrum (see Sect. 3.2.2) for optimal detection of the ISWlensing bispectrum, along with the bestfitting estimates of f_{NL} from the KSW method for different values of ℓ. It should be noted that the skewC_{ℓ} spectrum is not a fit to the KSW data points; its shape is fully fixed by the template under consideration, with only the overall amplitude as a free parameter. Hence the agreement between the curve and the points is good evidence that KSW is really detecting the ISWlensing effect and not some other source of NG. Note that point sources, at the level determined by their own skewspectrum, do not contribute significantly to the ISWlensing statistic). See Planck Collaboration XVII (2014) and Planck Collaboration XIX (2014) for further information about the detection by Planck of the ISWlensing signal.
5.3. Pointsources bispectrum
ExtraGalactic point sources at Planck frequencies are divided into two broad categories: radio sources with synchrotron and/or freefree emission; and infrared galaxies with thermal emission from dust heated by young stars. Radio sources are dominant at central CMB frequencies up to 143 GHz, and can be considered unclustered (Toffolatti et al. 1998; GonzálezNuevo et al. 2005). Hence their bispectrum is constant and is related to their number counts as (83)with S the flux density, dn/dS the number counts per steradian, S_{c} the flux cut and k_{ν} the conversion factor from flux to relative temperature elevation, depending on the frequency and instrumental bandpass.
Infrared galaxies become important at higher frequencies, 217 GHz and above, and are highly clustered in dark matter halos, which enhances their bispectrum on large angular scales (Lacasa et al. 2012; Curto et al. 2013). However, in the Planck context it was shown by Lacasa & Aghanim (2014) that the IR bispectrum is more than 90% correlated with the Poissonian template of the radio sources. So a joint estimation of f_{NL} with a Poissonian bispectrum template will essentially account for the IR signal, and provide quasiidentical values compared to an analysis accounting for the IR bispectrum template. Indeed, in our final optimal bispectrum constraints for primordial shapes, we will account for the potential contamination from point sources by jointly fitting primordial and Poisson templates to the data.
Our final measured pointsource bispectrum amplitudes from the data are reported in Table 3. The amplitude is expressed in dimensionless units, i.e., it has been divided by the appropriate power of the monopole temperature T_{0}, and has been multiplied by 10^{29}. As shown in Sect. 8.1, the Poisson template is the only one that still evolves significantly between ℓ = 2000 and ℓ = 2500. This explains the differences between the values of the KSW and binned (that use ℓ_{max} = 2500) and the modal (that uses ℓ_{max} = 2000) estimators. It has been shown that for the same value of ℓ_{max} all three estimators agree very well.
We finally conclude from Table 3 that we detect the pointsource bispectrum with high significance in the SMICA, NILC, and SEVEM cleaned maps, while it is absent from the CR cleaned map. The measured skewC_{ℓ} spectrum of the SMICA map in the bottom figure of Fig. 2 gives further evidence that the NG from foreground point sources is convincingly detected. The only degree of freedom in this plot is the amplitude, which is not set by a direct fit to the skewC_{ℓ}, but rather is estimated by KSW. As a result, the good agreement with the shape of this skewC_{ℓ} spectrum is powerful evidence that there is NG from point sources. However, this still turns out to be a negligible contaminant for primordial f_{NL} studies, due to the very low correlation between the Poisson bispectrum and the primordial shapes.
Results for the amplitude of the point source (Poisson) bispectrum (in dimensionless units of 10^{29}) from the SMICA, NILC, SEVEM, and CR foregroundcleaned maps, for the KSW, binned, and modal (polynomial) estimators; error bars are 68% CL.
5.4. Nonprimordial contributions to the trispectrum
The main noninstrumental source of nonprimordial signal is the kinematic modulation dipole due to the peculiar velocity of the earth, v, whose magnitude is (Challinor & van Leeuwen 2002; Kosowsky & Kahniashvili 2011; Amendola et al. 2011). If data are used to constrain τ_{NL} using the dipole modulation (which shrinks the Fisher error by a factor of two relative to starting at L = 2), the dipoleinduced signal must be subtracted, since its modulation reconstruction has signaltonoise larger than unity at Planck resolution. Confirmation that this signal is detected with the expected magnitude and direction is a good test of our methodology. The dipole signal seen by Planck is studied in detail in Planck Collaboration XXVII (2014), so we only summarize the key points here.
The local Doppler effect modulates the observed CMB temperature by at leading order, so that . The spectrum in each direction remains a blackbody, but the relative response in the intensity at the observed Planck frequencies is however frequency dependent. The effective thermodynamic fractional temperature anisotropy ΔΘ at each frequency for zero peculiar velocity is defined by (84)With peculiar velocity the temperature T depends on the second order term , so expanding the Planck function to second order then gives a change in the effective smallscale temperature anisotropy from both first and second order terms in the expansion of I_{ν}: (85)where x ≡ hν/k_{B}T_{0} (and we neglect small secondorder nonmodulation terms). Thus the anisotropies in the Planck maps have a dipolar modulation given by , where for the frequency bands we use b_{143} ≈ 2, b_{217} ≈ 3, and β ≡  v  = 1.23 × 10^{3} in the direction of CMB dipole. In addition our peculiar velocity induces kinetic aberration, which looks at leading order exactly like a dipole lensing convergence and only projects weakly into the power anisotropy estimator. For constraining τ_{NL} both of the expected kinematic signals can be included in the simulations, and hence subtracted in the mean field of the modulation reconstruction.
Secondary effects are dominated by the significant and very blue lensing signal. However unlike for the f_{NL} bispectrum lensing only overlaps with τ_{NL} at a small fraction of the error bar as long as only low modulation multipoles L ≲ 10 are used, where the τ_{NL} signal peaks (Pearson et al. 2012). We include lensing in the simulations, so lensing is straightforwardly accounted for in our analysis by its inclusion in the noise bias (Eq. (74)) and mean field.
A variety of instrumental effects can also give a spurious modulation signal if not modelled accurately. In particular the mean field due to anisotropic noise is very large (Hanson et al. 2009a). On the ultralarge scales of interest for τ_{NL}, our understanding of the noise is not adequate to calculate accurately and subtract this large signal. Instead, as for the power spectrum estimation, we use crossmap estimators that have no noise mean field on average. Both the noise and most other instrumental effects such as gain variations are expected to produce a signal with approximate symmetry about the ecliptic plane. Our modulation reconstruction methodology is especially useful here, since we can easily inspect the orientation of any signal found; for example a naive treatment of the noise not using crossmaps would give a large apparent quadrupolar modulation signal aligned with the ecliptic, corresponding to percentlevel misestimation of the noise mean field from inaccurate noise simulation.
Beam asymmetries are included in the simulation, as described in Planck Collaboration XVII (2014), but their effect is very small, since the modulation we are reconstructing is isotropic.
Since we are reconstructing a modulation of smallscale power, the estimator is totally insensitive to smooth largescale foregrounds. However largescale variation in smallscale foreground power can mimic a trispectrum modulation. We project out 857 GHz as a dust template in our inversevariance filtering procedure, as described in Appendix A of Planck Collaboration XVII (2014), but do not include any other foreground model in the trispectrum analysis. Any unmodelled foreground power variation would increase the τ_{NL} signal, so our modelling is sufficient to place a robust upper limit.
6. Validation tests
The f_{NL} results quoted in this paper have all been crossvalidated using multiple bispectrumbased estimators from different groups. Having multiple estimators was extremely useful for the entire analysis, for two main reasons. First, it allowed great improvement in the robustness of the final results. In the early stages of the work the comparison between different independent techniques helped to resolve bugs and other technical issues in the various computer codes, while during the later stages it was very useful to understand the data and find the optimal way of extracting information about the various bispectrum templates. Secondly, besides these crosschecking purposes, different estimators provide also interesting complementary information, going beyond simple f_{NL} estimation. For example, the binned and modal estimators provide a reconstruction of the full bispectrum of the data (smoothed in different domains), the skewC_{ℓ} estimator allows monitoring of the contribution to f_{NL} from different sources of NG, the wavelets reconstruction allows f_{NL} directionality tests, and so on.
In this section we are concerned with the first point above, that is, the use of multiple bispectrumbased pipelines as a way to improve the robustness of the results. For this purpose, a large amount of work was dedicated to the development and analysis of various test maps, in order to validate the estimators. This means not only checking that the various estimators recover the input f_{NL} within the expected errors, but also that the results agree on a mapbymap basis.
The section is split into two parts. Section 6.1 shows results on a set of initially fullsky, noiseless, Gaussian CMB simulations, to which we add, in several steps, realistic complications, including primordial NG, anisotropic coloured noise, and a mask, showing the impact on the results at each step. In Sect. 6.2 we show our results on a set of simulations that mimic the real data as closely as possible (except for the presence of foreground residuals, which will be studied in Sect. 8.4): no primordial NG, but NG due to the ISWlensing effect; simulated instrumental effects and realistic noise; and simulations passed through the component separation pipelines. In fact these are the FFP6 simulations (see Appendix A) that are used to determine the error bars for the final Planck results.
We present here only a small subset of the large number of validation tests that were performed. For example, we also had a number of “blind f_{NL} challenges”, in which the different groups received a simulated data set with an unknown value of input f_{NL} for a given shape and they had to report their estimated values. In addition different noise models were tested (white vs. coloured and isotropic vs. anisotropic), leading to the conclusion that it is important to make the noise in the estimator calibration as realistic as possible (coloured and anisotropic). We also tested different Galactic and point source masks, with and without inpainting, concluding that it is best to fill in both the point sources and the Galactic mask, using a sufficient number of iterations in our diffusive procedure to entirely fill in the point source gaps, while at the same time only effectively apodizing the Galactic mask (no smallscale structure in its interior). There were also various tests on FFP simulations of Planck data with Gaussian or nonGaussian CMB and all foregrounds, provided by the PSM (see Appendix A.4). These simulations were tested both before and after they had passed through the component separation pipelines. In all comparison tests the results were consistent with input f_{NL} values and differences between estimators were consistent with theoretical expectations.
6.1. Validation of estimators in the presence of primordial nonGaussianity
The aim of the first set of validation tests is threefold. First, we want to study the level of agreement from different estimators in ideal conditions (i.e., fullsky noiseless data). The expected scatter between measurements is, in this case, entirely due to the slightly imperfect correlation between weights of estimators that adopt different schemes to approximate the primordial shape templates. For this case the scatter can be computed analytically (see Appendix B for details). We can then verify that our results in ideal conditions match theoretical expectations. This is done in Sect. 6.1.1. Second, we want to make sure that the estimators are unbiased and correctly recover f_{NL} in input for local, equilateral, and orthogonal shapes. This is done in Sects. 6.1.2 and 6.1.3, where a superposition of local, equilateral and orthogonal bispectra is included in the simulations and the three f_{NL} values are estimated both independently and jointly. Finally we want to understand how much the agreement between pipelines in ideal conditions is degraded when we include a realistic correlated noise component and a sky cut, thus requiring the introduction of a linear term in the estimators in order to account for offdiagonal covariance terms introduced by the breaking of rotational invariance. Since we want to study the impact of adding noise and masking separately, we will first work on a set of fullsky maps with noise in Sect. 6.1.2, and then add a mask in Sect. 6.1.3.
The tests were applied to the KSW, binned and modal estimators. These are the bispectrum pipelines used to analyse Planck data in Sect. 7. Our goal for this set of tests is not so much to attain the tightest possible agreement between methods, as it is to address the points summarized in the above paragraph. For this reason the estimator implementations used in this specific Section were slightly less accurate but faster to compute than those adopted for the final data analysis of Sect. 7. The primary difference with respect to the main analysis is that a smaller number of simulations was used to calibrate the linear term (80–100 in these tests, as against 200 or more for the full analysis). For the modal estimator we also use a faster expansion with a smaller number of modes: 300 here versus 600 in the high accuracy version of the pipeline^{9} used in Sect. 7. Even with many fewer modes, the modal estimator is still quite accurate: the correlation coefficient for the modal expansion of the local template is 0.95, while for the equilateral and orthogonal shapes it is 0.98.
6.1.1. Ideal Gaussian simulations
As a basis for the other tests we start with the ideal case, a set of 96 simulations of a fullsky Gaussian CMB, with a Gaussian beam with FWHM 5 arcmin and without any noise, cut off at ℓ_{max} = 2000 in our analyses. The independent Fisher matrix error bars in that case are 4.2 for local NG, 56 for equilateral, and 28 for orthogonal.
Note that this test does not make sense for all estimators, and hence results are not included for all of them. For example, for the binned estimator the optimal binning depends on the noise. While this dependence is not very strong, the difference between no noise and Planck noise is sufficiently large that a completely different binning would have to be used just for this test, going against the purpose of this section to validate the estimators as used for the data analysis^{10}.
The purpose here is mostly aimed at checking consistency with the following formula (derived in Appendix B) for the expected scatter (standard deviation) between f_{NL} results of the same map from an exact and an approximate estimator: (86)Here Δ_{th} is the standard deviation of the exact estimator and r is the correlation coefficient that gives the correlation of the approximate bispectrum template with the exact one, defined as (87)where the label “th” denotes the initial bispectrum shape to fit to the data, and “exp” is the approximate expanded one. Note that this formula has been obtained under the simplifying assumptions of Gaussianity, fullsky coverage and homogeneous noise. For applications dealing with more realistic cases we might expect the scatter to become larger, while remaining qualitatively consistent.
The results averaged over the whole set of maps are given in Table 4 for the KSW and modal estimators individually, as well as for their difference. The plane wave modal expansion implemented here achieves about 98% correlation with the separable shapes used by KSW. According to the formula above we then expect a standard deviation of mapbymap differences of order 0.2Δ_{fNL} for a given shape, where Δ_{fNL} is the corresponding f_{NL} error bar. Looking at the lefthand side of Table 4, we see that the error bars are 4 for local NG, 50 for equilateral, and 30 for orthogonal. So we predict a standard deviation of mapbymap differences of 0.8, 10 and 6 for local, equilateral, and orthogonal NG, respectively. As one can see from the “ModalKSW” column, the measurements are in excellent agreement with the theoretical expectation.
6.1.2. NonGaussian simulations with realistic noise
A set of 96 fullsky nonGaussian CMB simulations was created according to the process described by Fergusson et al. (2010a), with local , equilateral , and orthogonal . The effect of a 5 arcmin beam was added, as well as realistic coloured and anisotropic noise according to the specifications of the SMICA cleaned map. The independent Fisher matrix error bars in that case are 5.3 for local, 63 for equilateral, and 33 for orthogonal NG, while the joint ones are respectively 6.0, 64, and 37.
The results averaged over the whole set are given in Table 5 for the various estimators individually, as well as for the differences with respect to KSW. Compared to the previous case we now deviate from the exact theoretical expectation for two reasons: we include a realistic correlated noise component; and we have NG in the maps. The presence of NG in the input maps will lower the agreement between estimators with respect to the Gaussian case if the correlation between weights is not exactly 100%. This is even more true in this specific case, where NG of three different kinds is present in the input maps and also crosscorrelation terms between different expanded shapes are involved (and propagated over in the joint analysis). Moreover, when noise is included the specific modal expansion used for this test is 95% correlated to the separable KSW local shape (so there is a 3% reduction of the correlation compared to the ideal case for the modal local shape); we thus expect a further degradation of the level of agreement for this specific case. Finally, in order to correct for noise effects, a linear term has to be added to the estimators. Since the linear term is obtained by MC averaging over just 80 or 96 simulations in this test (depending on the estimator), MC errors are also adding to the measured differences. Of course the MC error can be reduced by increasing the number of simulations in the linear term sample. We do this for the analysis of the real data and in Sect. 6.2, but it was computationally too expensive for this set of preliminary validation tests, so we decided here to just account for it in the final interpretation of the results.
As a consequence of the above, we can no longer expect the mapbymap f_{NL} differences to follow perfectly the theoretical expectation, obtained in the previous section in idealized conditions (fullsky, no noise, and Gaussianity). With these caveats in mind, the agreement between different pipelines remains very good, being about 0.3σ in most cases and about 0.5σ for the modalKSW difference in the local case, which can be easily explained by the fact that this is the set of weights with the lowest correlation (95%, as mentioned above). All estimators are unbiased and recover the correct input values.
6.1.3. Impact of the mask
To the simulations of Sect. 6.1.2 we now apply the Planck union mask – denoted U73 – masking both the Galaxy and the brightest point sources and leaving 73% of the sky unmasked (Planck Collaboration XII 2014). This is the same mask used to analyse Planck data in Sect. 7. The independent Fisher matrix error bars in that case (taking into account the f_{sky} correction) are 6.2 for local NG, 74 for equilateral, and 39 for orthogonal, while the joint ones are respectively 7.1, 76, and 44.
All masked areas of the sky (both Galactic and point sources) are filled in with a simple iterative method. In this simple inpainting method each pixel in the mask is filled with the average of all eight surrounding pixels, and this is repeated 2000 times over all masked pixels. The fillingin helps to avoid propagating the effect of a sharp edge and the lack of largescale power inside the mask to the unmasked regions during harmonic transforms. This inpainting method is the one that was used to produce all NG results in this paper for methods that need it (KSW, binned and modal).
The results averaged over the whole set of simulations are given in Table 6 for the various estimators individually, as well as for the differences with respect to KSW. The mapbymap results are shown in Fig. 3.
This is the most realistic case we consider in this set of tests. Besides noise, we also include a sky cut and our usual mask inpainting procedure. All the caveats mentioned for the previous case are still valid, and possibly emphasized by the inclusion of mask and inpainting. In the light of this, the agreement is still very good, worsening a bit with respect to the “fullsky + noise” case only for the local measurement, where the mask is indeed expected to have the biggest impact. In the joint analysis all estimators recover the correct input values for the local and orthogonal cases, but all estimators find a value for equilateral NG that is somewhat too low. It is unclear whether this is an effect of masking and inpainting on the equilateral measurement or just a statistical fluctuation for this set of simulations. In any case, this potential bias is small compared to the statistical uncertainty, so that it would not have a significant impact on the final results.
To summarize the results of this Sect. 6.1, we performed an extensive set of validation tests between different f_{NL} estimators using strongly, but not perfectly, correlated primordial NG templates in their weights. The test consisted in comparing the f_{NL} measured by the different estimators for different sets of simulations, on a mapbymap basis. We started from ideal conditions: fullsky Gaussian noiseless maps. In this case we computed a theoretical formula providing the expected standard deviation of the f_{NL} differences, as a function of the correlations between the input NG templates in the different estimators. Our results match this formula very well. In the other two simulation sets we added realistic features (noise, mask and inpainting) and we included a linear combination of local, equilateral and orthogonal NG. First of all we verified that all the pipelines correctly recover the three f_{NL} input values, hence they are unbiased. Moreover, we observed that adding such features produces an expected slight degradation of the level of agreement between different pipelines, that nevertheless remains very good: about 0.3–0.4σ for equilateral and orthogonal NG, and about 0.5–0.6σ for local NG, which is the shape most affected by mask and noise contamination.
6.2. Validation of estimators on realistic Planck simulations
In the tests of the previous subsection we checked the bias of the estimators and studied their level of agreement, given the correlation between their weights, in the presence of noise and a sky cut. To speed up the computation while still retaining enough accuracy for the purposes of that analysis, we used a relatively small number of maps for linear term calibrations (80–100) and used a smaller number of modes than usual in the modal estimator. In the present subsection we instead try to simulate as accurately as possible real data analysis conditions. Our goal is to obtain an accurate MCbased expectation of the scatter between different f_{NL} measurements when the pipelines are run on actual Planck maps.
Fig. 3
Mapbymap comparison of the results from the different estimators for local (top), equilateral (centre), and orthogonal (bottom) f_{NL} for the set of masked nonGaussian simulations described in Sect. 6.1.3, assuming the shapes to be independent. The horizontal solid line is the average value of all maps for KSW, and the dashed and dotted horizontal lines correspond to 1σ and 2σ deviations, respectively. 
Fig. 4
Mapbymap comparison of the results from the different estimators for local (top), equilateral (centre), and orthogonal (bottom) f_{NL} for 99 maps from a set of realistic lensed simulations passed through the SMICA pipeline, described in Sect. 6.2, assuming the shapes to be independent. The horizontal solid line is the average value of the maps for KSW, and the dashed and dotted horizontal lines correspond to 1σ and 2σ deviations, respectively. 
To this aim we use FFP6 simulation maps described in Appendix A. The original FFP6 maps were lensed using the Lenspix algorithm, and processed through the SMICA component separation pipeline. They were then multiplied by the Galactic and point source mask U73 as in the actual f_{NL} analysis, and inpainted as usual. Since our final results show full consistency with Gaussianity for local, equilateral and orthogonal shapes, we do not include any primordial f_{NL} in these maps. We note that although the simulations were passed through SMICA in order to provide a realistic filtering of the data, they did not include any foreground components. The impact of foreground residuals will be studied separately in Sect. 8.4.
The configuration of all bispectrum pipelines was the same as used for the final data analysis, which implies a correlation of 99% or better between the weights of the KSW, binned and modal estimators. Linear terms were calibrated using 200 simulations, after verifying that this number allows accurate convergence for all the shapes. For this test we also included the wavelet bispectrum pipeline described in Sect. 3. Although this last estimator turns out to be about 30% suboptimal and, in its current implementation, less correlated with the primordial templates than the other estimators, it does provide an additional interesting crosscheck of our results by introducing another decomposition basis. We thus used it to analyse SMICA data in Sect. 7, while the other three pipelines were used on all maps.
A comparison of the measured f_{NL} mapbymap for all shapes and estimators is shown in Fig. 4. As an overall figure of merit of the level of agreement achieved by different pipelines we take as usual the standard deviation of the mapbymap f_{NL} differences, σ_{δfNL}. Table 7 shows that the final agreement between the three optimal pipelines (KSW, binned, and modal) is close to saturating the ideal bound in Eq. (86) determined by the imperfect correlation of the weights, i.e., it varies from about once to twice σ_{δfNL} ≃ 0.15 Δ_{fNL} for an r = 0.99 correlation. This is very consistent with the level of agreement that we find between estimators for the final results from the data, providing a good indication that no spurious NG features are present in the actual data set when compared to our simulations. It should be noted that we found a similarly good level of agreement between estimators for the nonprimordial shapes of point sources and ISWlensing, although we chose not to present those results here in order to focus on the primordial shapes. Finally, regarding the wavelet pipeline, the lower weight correlation and suboptimal error bars produce an expected larger scatter when compared to the other estimators. Nonetheless, the level of agreement is still of order 1σ, which is quite acceptable for consistency checks of the optimal results. Again, this MC expectation agrees with what we see in our results on the real data.
7. Results
For our analysis of Planck data we considered foregroundcleaned maps obtained with the four component separation methods SMICA, NILC, SEVEM, and CR. For each map, f_{NL} amplitudes for the local, equilateral, and orthogonal primordial shapes have been measured using three (four for SMICA) bispectrum estimators described in Sect. 3. The results can be found in Sect. 7.1. These estimators, as explained earlier, basically use an expansion of the theoretical bispectrum templates in different domains, and truncate the expansion when a high level of correlation with the primordial templates is achieved. These accurate decompositions, which are highly correlated with each other, are then matched to the data in order to extract f_{NL}. The different expansions are all different implementations of the maximumlikelihood estimator given in Eq. (33). So the final estimates are all expected to be optimal, and measure f_{NL} from nearly identical fitting templates. As discussed and tested in detail on simulations in Sect. 6, central f_{NL} values from different methods are expected to be consistent with each other within about 0.3σ_{fNL}. It is then clear that comparing outputs from both different estimators and different component separation methods, as we do, allows for stringent internal consistency checks and improved robustness of the final f_{NL} results.
In addition, the binned and modal techniques produce shapeindependent full bispectrum reconstructions in their own different domains. These reconstructions, discussed in Sect. 7.2, complement the standard f_{NL} measurements in an important way, since they allow detection of possible NG features in the threepoint function of the data that do not correlate significantly with the standard primordial shapes. This advantage is shared by the skewC_{ℓ} method, also applied to the data. A detection of such features would either produce a warning that some residual spurious NG effects are still present in the data or provide an interesting hint of “nonstandard” primordial NG that is not captured by the local, equilateral and orthogonal shapes. Additional constraints for a broad range of specific models are provided in Sect. 7.3 (see also Sect. 2.3). In Sect. 7.4 we present the constraints on local NG obtained with Minkowski Functionals. Finally, in Sect. 7.5 we present our CMB trispectrum results.
7.1. Constraints on local, equilateral and orthogonal f_{NL}
Fig. 5
Binned skewC_{ℓ} statistics from the SMICA map for a) local; b) equilateral; and c) orthogonal. Theoretical curves are not fitted to the data shown, but are plotted with the amplitude (the only free parameter) determined from the KSW technique. There is no evidence for detection of primordial NG. Error bars are derived from the covariance of estimates from 1000 simulations. There are mild correlations between data points in all figures, but very strong correlations at high ℓ in the local case, reaching correlation coefficients r> 0.99 for ℓ > 1750. 
Our goal here is to investigate the standard separable local, equilateral and orthogonal templates used e.g., in previous WMAP analyses (see e.g., Bennett et al. 2013). When using the modal, binned, or wavelet estimator, these theoretical templates are expanded approximately (albeit very accurately) using the relevant basis functions or bins. On the other hand, the KSW estimator by construction works with the exact templates and, for this reason, it is chosen as the baseline to provide the final f_{NL} results for the standard shapes (local, equilateral, orthogonal), see Table 8. However, both the binned and modal estimators achieve optimal performance and an extremely high correlation for the standard templates (~99%), so they are statistically equivalent to KSW, as demonstrated in the previous section. This means that we can achieve a remarkable level of crossvalidation for our Planck NG results. We will be able to present consistent constraints for the local, equilateral and orthogonal models for all four Planck foregroundcleaned maps, using three independent optimal estimators (refer to Table 9). Regarding component separation methods, we adopt the SMICA map as the default for the final KSW results given its preferred status among foregroundseparation techniques in Planck Collaboration XII (2014). The other component separation maps will be used for important crossvalidation of our results and to evaluate potential sensitivity to foreground residuals.
All the results presented in this section were obtained using the union mask U73, which leaves 73% of the sky unmasked. The mask is the union of the confidence masks of the four different component separation methods, where each confidence mask defines the region where the corresponding CMB cleaning is trusted (see Planck Collaboration XII 2014). As will be shown in Sect. 8.2, results are robust to changes that make the mask larger, but choosing a significantly smaller mask would leave some NG foreground contamination. For the linear term CMB and noise calibration, and error bar determination, we used sets of realistic FFP6 maps that include all steps of data processing, and have realistic noise and beam properties (see Appendix A). The simulations were also lensed using the Lenspix algorithm and filtered through the component separation pipelines.
In Table 8 we show results for the combination of the KSW estimator and the SMICA map, at a resolution of ℓ_{max} = 2500. We present both “independent” singleshape results and “ISWlensing subtracted” ones. The former are obtained by directly fitting primordial templates to the data. For the latter, two additional operations have been performed. In the first place, as the name indicates, they have been corrected by subtracting the bias due to the correlation of the primordial bispectra to the latetime ISWlensing contribution (Mangilli & Verde 2009; Junk & Komatsu 2012; Hanson et al. 2009b, see Sect. 5.2). In addition, a joint fit of the primordial shape with the (Poissonian) pointsource bispectrum amplitude extracted from the data has been performed on the results marked “ISWlensing subtracted”^{11}. Since the ISWlensing bispectrum is peaked on squeezed configurations, its impact is well known to be largest for the local shape. The ISWlensing bias is also important for orthogonal measurements (there is a correlation coefficient r ~ − 0.5 between the local and orthogonal CMB templates), while it is very small in the equilateral limit. The values of the ISWlensing bias we subtract, summarized in Table 1, are calculated assuming the Planck bestfit cosmological model as our fiducial model. The same fiducial parameters were of course consistently used to compute the theoretical bispectrum templates and the estimator normalization. Regarding the point source contamination, we detect a Poissonian bispectrum at high significance in the SMICA map, see Sect. 5.3. However, marginalizing over point sources has negligible impact on the final primordial f_{NL} results, because the Poisson bispectrum template has very small correlations with all the other shapes.
In light of the discussion at the beginning of this section, we take the numbers from the KSW SMICA analysis in Table 8 as the final local, equilateral and orthogonal f_{NL} constraints for the current Planck data release. These results clearly show that no evidence of NG of the local, equilateral or orthogonal type is found in the data. After ISWlensing subtraction, all f_{NL} for the three primordial shapes are consistent with 0 at 68% CL. Note that these numbers have been crosschecked using two completely independent KSW pipelines, one of which is an extension to Planck resolution of the pipeline used for the WMAP analysis (Bennett et al. 2013).
Results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW estimator from the SMICA foregroundcleaned map. Both independent singleshape results and results marginalized over the point source bispectrum and with the ISWlensing bias subtracted are reported; error bars are 68% CL.
Results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW, binned and modal estimators from the SMICA, NILC, SEVEM, and CR foregroundcleaned maps.
Unlike other methods, the KSW technique is not designed to provide a reconstruction of the full bispectrum of the data. However, the related skewC_{ℓ} statistic described in Sect. 3.2.2 allows, for each given shape, visualization and study of the contribution to the measured f_{NL} from separate ℓbins. This is a useful tool to study potential spurious NG contamination in the data. We show for the SMICA map in Fig. 5 the measured skewC_{ℓ} spectrum for optimal detection of primordial local, equilateral and orthogonal NG, along with the bestfitting estimates of f_{NL} from the KSW method. Contrary to the case of the point source and ISWlensing foregrounds (see Sect. 5), the skewC_{ℓ} statistics do not show convincing evidence for detection of the primordial shapes, suggesting that these primordial effects are not significant sources of the NG in the map. Again, point sources contribute very little to this statistic; ISWlensing contributes, but only a small fraction of the amplitude.
As mentioned before, our analysis went beyond the simple application of the KSW estimator to the SMICA map. All f_{NL} pipelines developed for Planck analysis were applied to all componentseparated maps by SMICA, NILC, SEVEM, and CR. We found from simulations in the previous sections that the KSW, binned, and modal pipelines saturate the CramérRao bound, while the wavelet estimator in its current implementation provides slightly suboptimal results. Wavelets remain however a useful crosscheck of the other methods, also given some technical complementarities, e.g., they are the only approach that does not require inpainting, as explained in Sect. 3. Hence we include wavelet results, but only for SMICA. The f_{NL} results for the optimal KSW, binned and modal bispectrum estimators, for the four component separation methods, are summarized in Table 9, one of the main products of our analysis of Planck data. The wavelet bispectrum analysis of SMICA is reported in Table 10. In the analysis, the KSW and binned bispectrum estimators considered multipoles up to ℓ_{max} = 2500, while the modal estimator went to ℓ_{max} = 2000. As shown in Sect. 8.1 and Table 16, error bars and central values for the three primordial shapes have converged at ℓ_{max} = 2000, so the final primordial f_{NL} estimates from the three pipelines are directly comparable^{12}.
The binned bispectrum estimator used 51 bins, which were determined by optimizing the expected variance of the different f_{NL} parameters, focusing in particular on the primordial shapes.^{13} The modal estimator employed a polynomial basis (n_{max} = 600) previously described in Fergusson et al. (2010a), but augmented with a local shape mode (approximating the SW largeangle local solution) to improve convergence in the squeezed limit. The above choices for the binned and modal methods produce a very high correlation (generally 99% or better) of the expanded/binned templates with the exact ones used by the KSW estimator. The wavelet estimator is based on thirdorder statistics generated by the different possible combinations of the wavelet coefficient maps of the SMHW evaluated at certain angular scales. See for example Antoine & Vandergheynst (1998) and MartínezGonzález et al. (2002) for detailed information about this wavelet. We considered a set of 15 scales logarithmically spaced between 1.3 and 956.3 arcmin and we also included the unconvolved map. The wavelet map w(R_{i};b) (Eq. (61)) for each angular scale R_{i} has an associated extended mask generated from the mask U73 following the procedure described and extensively used in Cutro et al. (2009b,a, 2011a,b, 2012), Donzelli et al. (2012); Regan et al. (2013). The wavelet coefficient maps are later combined into the thirdorder moments q_{ijk} (Eq. (60)), for a total 816 different statistics, and these statistics are used to constrain f_{NL} through a χ^{2} test.
Results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the suboptimal wavelet estimator from the SMICA foregroundcleaned map.
The high level of agreement between results from the KSW, binned and modal f_{NL} estimators, and from all the component separation pipelines, is representative of the robustness of our results with respect to residual foreground contamination, and is fully consistent with our preliminary MC analysis shown in Sect. 6. The scatter with wavelets is a bit larger, but this was expected due to the suboptimality of the wavelet estimator and is also in agreement with our MC expectations from the tests. Therefore wavelets do provide another successful crosscheck.
7.2. Bispectrum reconstruction
As previously explained (see Sect. 3), in addition to looking in specific bispectrumspace directions and extracting the single number f_{NL} for given shapes, the binned and modal pipelines have the capability to generate a smoothed (i.e., either coarsegrained in ℓspace, or truncated at a given expansion eigenmode) reconstruction of the full bispectrum of the data. See also Planck Collaboration XXIII (2014).
7.2.1. Modal bispectrum reconstruction
Fig. 6
Full 3D CMB bispectrum recovered from the Planck foregroundcleaned maps, including SMICA (left), NILC (centre) and SEVEM (right), using the hybrid Fourier mode coefficients illustrated in Fig. 9, These are plotted in threedimensions with multipole coordinates { ℓ_{1},ℓ_{2},ℓ_{3} } on the tetrahedral domain shown in Fig. 1 out to ℓ_{max} = 2000. Several density contours are plotted with red positive and blue negative. The bispectra extracted from the different foregroundseparated maps appear to be almost indistinguishable. 
Fig. 7
Planck CMB bispectrum detail in the signaldominated regime showing a comparison between full 3D reconstruction using hybrid Fourier modes (left) and hybrid polynomials (right). Note the consistency of the main bispectrum properties which include an apparently “oscillatory” central feature for lowℓ together with a flattened signal beyond to ℓ ≲ 1400. Note also the periodic CMB ISWlensing signal in the squeezed limit along the edges of the tetrapyd. 
The modal pipeline was applied to the Planck temperature maps for the foregroundseparation techniques SMICA, NILC, and SEVEM (Fergusson et al. 2010a). For this analysis we used two alternative sets of hybrid basis functions in order to crosscheck results and identify particular signals. First, we employed trigonometric functions (n_{max} = 300) augmented with the SW local mode, together with the three separable modes contributing to the CMB ISWlensing signal. Secondly, we employed the same polynomial basis (n_{max} = 600) with local SW mode as was used for f_{NL} estimation.
The modal coefficients extracted from the Planck SMICA, NILC, and SEVEM maps are shown in Fig. 9. Here we have used the hybrid Fourier modes with local and ISWlensing modes. These amplitudes show remarkable consistency between the different maps, demonstrating that the alternative foreground separation techniques do not appear to be introducing spurious NG. Note that here the coefficients are for the orthonormalized modes R_{n} (Eq. (65)) and they have a roughly constant variance, so anomalously large modes can be easily identified. It is evident, for example, that among the low modes there are large signals, which include the ISWlensing signal and point source contributions.
Using the modal expansion of Eq. (46) with Eq. (65), we have reconstructed the full 3D Planck bispectrum. This is illustrated in Fig. 6, where we show “tetrapyd” comparisons between different foreground cleaned maps. The tetrapyd (see Fig. 1) is the region defined by the multipoles that obey the triangle condition, with ℓ ≤ ℓ_{max}. The 3D plots show the reduced bispectrum of the map, divided by a SachsWolfe CMB bispectrum solution for a constant primordial shape, S(k_{1},k_{2},k_{3}) = 1. This constant primordial bispectrum template normalizaton is carried out in order to remove an ~ℓ^{4} scaling from the starting bispectrum (it is analogous to multiplication of the power spectrum by ℓ(ℓ + 1)). To facilitate the interpretation of 3D bispectrum figures, note that squeezed configurations lie on the edges of the tetrapyd, flattened on the faces and equilateral in the interior, with b_{ℓℓℓ} on the diagonal. The colour levels are equally spaced with red denoting positive values, and blue denoting negative. Given the correspondence of the coefficients for SMICA, NILC, and SEVEM, the reconstructed 3D signals also appear remarkably consistent, showing similar contours out to ℓ ≲ 1500. At large multipoles ℓ approaching ℓ_{max} = 2000, there is increased randomness in the reconstruction due to the rise in experimental noise and some evidence for a residual point source contribution.
There are some striking features evident in the 3D bispectrum reconstruction which appear in both Fourier and polynomial representations, as shown in more detail in Fig. 7. There is an apparent oscillation at low ℓ ≲ 500 already seen in WMAP7 (Fergusson et al. 2012). Beyond out to ℓ ~ 1200 there are further distinct features (mostly “flattened” on the walls of the tetrapyd), and an oscillating ISWlensing contribution can be discerned in the squeezed limit. For comparison, the bispectrum reconstruction for one of the lensed Gaussian simulations described in Sect. 6.2 is shown in Fig. 8. ISWlensing oscillating signatures in the squeezed limit are visible also in this case. When comparing results from a single simulation to the data, it is however always useful to keep in mind that the observed tetrapyd pattern is quite realization dependent, since full bispectrum reconstructions are noisy. Whatever its origin, Gaussian or otherwise, Fig. 7 reveals the CMB bispectrum of our Universe as observed by Planck.
Fig. 8
CMB bispectrum from a Gaussian simulation including gravitational lensing. This reconstruction uses the same hybrid polynomials as for the Planck bispectrum in Fig. 7 (right), with which it can be compared. Note that an indication of a ISWlensing bispectrum signal can be seen along the edges of the tetrahedron. 
Fig. 9
Modal bispectrum coefficients for the mode expansion (Eq. (65)) obtained from Planck foregroundcleaned maps using hybrid Fourier modes. The different component separation methods, SMICA, NILC and SEVEM exhibit remarkable agreement. The variance from 200 simulated noise maps was nearly constant for each of the 300 modes, with the average ±1σ variation shown in red. 
Fig. 10
Total integrated bispectrum defined in Eq. (66) as a cumulative sum over orthonormal modal coefficients (upper panel) and over multipoles up to a given ℓ (lower panel). Above, the relative quantity is plotted, where is the mean obtained from 200 CMB Gaussian maps with the standard deviation shown as the red line. Below the square of the bispectrum is integrated over the tetrapyd out to ℓ and its significance plotted relative to the Gaussian standard deviation (1σ red line). A hybrid polynomial basis n_{max} = 600 is employed in the signaldominated region ℓ ≤ 1500. 
The cumulative sum over the squared orthonormal coefficients from Eq. (66) for the Planck data is illustrated in Fig. 10 (upper panel). The Planck bispectrum contribution can be directly compared with Gaussian expectations averaged from 200 lensed Gaussian maps with simulated residual foregrounds. It is interesting to note that the integrated bispectrum signal fairly consistently exceeds the Gaussian mean by around 2σ over much of the domain. This includes the ISW and PS contributions for which subtraction only has a modest effect. Also shown (lower panel) is the corresponding cumulative quantity as a function of multipole ℓ, for which features have visible counterparts at comparable ℓ in Fig. 7. Despite the high bispectrum signal, this χ^{2}test over the orthonormal mode coefficients is cumulatively consistent with Gaussianity for each of the three component separation methods considered. It is however important to stress that this is a model independent integrated measurement of NG. Fitting the data with specific bispectrum templates can of course enhance the signaltonoise ratio for given models, especially in light of the 1 to 2σ excess with respect to the Gaussian expectation, shown in the lower panel. This measurement is thus not in disagreement with detection of a residual point source (Poisson) bispectrum in the same maps (see Table 3), or with the results shown in our feature models survey of Sect. 7.3.3. We also note some differences in the high mode region between SEVEM and the other two methods, with the SEVEM results being closer to the Gaussian expectation. This is consistent with tests on simulations and with SEVEM measuring a slightly lower ISWlensing amplitude than the other methods (see Table 2). On the other hand the discrepancies are well within statistical bounds, as it can be seen by comparing them to the Gaussian standard deviation from simulations (red line). We thus conclude that, even accounting for these high modes deviations, the different component separation methods display a good level of internal consistency. It is also important to notice that this already good level of agreement becomes even stronger in “nonblind” f_{NL} measurements, when primordial bispectrum templates are fit to the data, and specific theoretically relevant regions of the tetrapyd are selected. This can be clearly seen from Table 9.
Fig. 11
Smoothed observed bispectrum as determined with the binned estimator divided by its expected standard deviation, as a function of ℓ_{1} and ℓ_{2}, with ℓ_{3} in the bin [610, 654]. From left to right on the top row are shown: SMICA, NILC, and SEVEM; and on the bottom row: CR and the raw 143 GHz channel. For comparison purposes the last figure on the bottom row shows the same quantity for one of the lensed Gaussian simulations described in Sect. 6.2. 
7.2.2. Binned bispectrum reconstruction
As explained in Sect. 3.4.2, it is interesting to study the smoothed observed bispectrum divided by its expected standard deviation, since this will indicate if there is a significant deviation from Gaussianity for certain regions of ℓspace. This quantity is shown in Figs. 11 and 12 as a function of ℓ_{1} and ℓ_{2}, for two different values (or rather, bins) of ℓ_{3}: the intermediate value [610, 654] in Fig. 11 and the high value [1330, 1374] in Fig. 12. Each figure shows the results for the SMICA, NILC, SEVEM, and CR cleaned maps as well as for the raw 143 GHz channel map. For comparison, the result for one of the lensed Gaussian simulations described in Sect. 6.2 is also shown. The bispectra were obtained with the binned bispectrum estimator and smoothed with a Gaussian kernel as explained in Sect. 3.4.2. Very blue or red regions indicate significant NG, regions that are less red or blue just represent expected fluctuations of a Gaussian distribution.
From Fig. 11 at an intermediate value of ℓ_{3} we can conclude that there is a very good agreement between SMICA, NILC, and SEVEM for all values of ℓ_{1} and ℓ_{2}, and with CR up to about ℓ_{1},ℓ_{2} ~ 1500. In fact, up to 1500 there is also a good agreement with the raw 143 GHz channel. We also see no significant nonGaussian features in this figure (except maybe in the CR and raw maps at ℓ_{1},ℓ_{2}> 2000). The lensed Gaussian simulation in the last panel looks quite Gaussian as well (as it should), but very different from the others. This is not surprising, as all we are seeing here are the small random fluctuations that exist in any Gaussian realization, and the simulation represents of course another realization than the real data. (Note that the ISWlensing NG, while present in the simulation, is not really visible in the particular slices shown here due to the linear scale of the axes and the fact that the ISWlensing NG peaks in the squeezed configuration.)
Figure 12 at a high value of ℓ_{3}, on the contrary, shows significant nonGaussian features in the raw map, but much less NG in the cleaned maps. In particular one can see the pointsource bispectral signal at highℓ approximately at equilateral configurations in the data. This is absent in the the lensed Gaussian simulation, which has no point sources. There is still an excellent agreement between SMICA, NILC, and SEVEM. The CR map shows less NG than the other three cleaned maps, which is consistent with the absence of a detection of the Poisson point source bispectrum for CR, see Table 3.
7.3. Constraints on specific targeted shapes
We have deployed the modal estimator to investigate a wide range of the inflationary models described in Sect. 2. This is the same validated estimator for which the standard f_{NL} results have been reported in the Sect. 7, but it is augmented with the primordial modal decomposition and projection described in Sect. 3.2.3. The resulting modalprojected local, equilateral and orthogonal shapes are ~99% correlated with those found using direct integration of Eq. (38) (as for the analysis above). Modal correlations for the other models investigated were determined for both the primordial shapes and the latetime projected decompositions and were all above 90%, unless stated otherwise. This primordial modal estimator pipeline has been applied already extensively to the WMAP7 data (Fergusson et al. 2012).
7.3.1. Nonseparable singlefield bispectrum shape results
Having characterised singlefield inflation bispectra using combinations of the separable equilateral and orthogonal ansätze, we note that the actual leadingorder nonseparable contributions (Eqs. (6), (7)) exhibit significant differences in the collinear (flattened) limit. For this reason we provide constraints on DBI inflation (Eq. (7)) and the two effective field theory shapes (Eqs. (5), (6)), as well as the ghost inflation bispectrum, which is an exemplar of higherorder derivative theories (specifically Eq. (3.8) in ArkaniHamed et al. 2004). Using the primordial modal estimator, with the SMICA foregroundcleaned data, we find: (88)where we have normalized with the usual primordial f_{NL} convention which is shapedependent (i.e., the central value of the shape function is taken such that S(k,k,k) = 1). In parentheses we also give a reweighted constraint for easier comparison with the equilateral constraint from the same modal estimator, i.e., we have rescaled using the Fisher variance for the closelyrelated equilateral shape. Given the strong crosscorrelation (above 95%) between all these models, the equilateral family results of (88) reveal larger differences around σ/3 than might be expected (and somewhat larger than observed previously in the WMAP data – Fergusson et al. 2012). The reason for this variation between the equilateral shapes in Planck appears to be the additional signal observed in the flattened limit in the bispectrum reconstruction beyond the WMAP signaldominated range (see Fig. 6). There is also a contribution from the small correlation difference between equilateral models from primordial modal and KSW methods. The results for these models for all the SMICA, NILC and SEVEM foregroundseparated maps are given in Appendix C (Table C.3).
Constraints on flattened or collinear bispectrum models (and related models) using the SMICA foregroundcleaned Planck map.
Planck bispectrum estimation results for feature models compared to the SMICA foregroundcleaned maps.
7.3.2. NonBunchDavies vacuum results
We have investigated the nonseparable shapes arising from excited initial states (nonBunchDavies vacuum models) which usually peak in the flattened or collinear limit. In particular, we have searched for the four nonseparable bispectra described in Eqs. (14) and (15), as well as the original flattened shape (Eqs. (6.2)–(3)) in Chen et al. 2007b). This entails choosing suitable cutoffs k_{c} to ensure that the signal is strongly flattened (i.e., distinct from flat in Eq. (13)), while also accurately represented by the modal expansion at both early and late times (Eqs. (55), (56)). For , we adopted the same edge truncation and mild Gaussian filter described in Fergusson et al. (2012), while for and , which are described by Eq. (14), we chose k_{c} = 0.001, and in Eq. (15) we take k_{c} = 0.01. The shape correlations for most nonBunchDavies vacua were good (above 90%), except for the strongly squeezed model with oscillations of Eq. (14) which was relatively poor (60%). Together with the orthogonal (Eq. (4)), flat (Eq. (13)) and vector (Eq. (19)) shapes, these nonBunchDavies models explore a broad range of flattened models, with a variety of different widths for picking out signals around the faces of the tetrapyd (see Fig. 1).
The f_{NL} results obtained for the nonBunchDavies models from the different foregroundcleaned map bispectra were consistent and the constraints from SMICA (for brevity) are given in Table 11. More comprehensive results from SMICA, NILC and SEVEM can be found in Table C.3 in Appendix C. Both and (Eq. (14)) produced raw results above 2σ, in part picking out the flattened signal observed in the bispectrum reconstruction in Fig. 6. However, these flattened squeezed signals are also correlated with CMB ISWlensing and so, after subtracting the predicted ISW bias (as well as the measured point source signal), most NBD f_{NL} results were reduced to 1σ or less (see “Clean f_{NL}” column in Table 11). The exception was the most flattened model which remained higher , i.e., with signals at 2.0σ, 1.8σ and 2.1σ for SMICA, NILC and SEVEM respectively.
We emphasise that this has to be considered just as preliminary study of flattened NG in the Planck data using four exemplar models. In order to reach a complete statistical assessment of constraints regarding flattened models in forthcoming analyses, we will have to undertake a systematic search for bestfit Planck NBD models using the parameter freedom available.
7.3.3. Scaledependent feature and resonant model results
We have investigated whether the Planck bispectrum reconstructions include oscillations expected in feature or resonant models (Eqs. (16), (17)). Although poorly correlated with scaleinvariant shapes, the feature and resonant models have (at least) two free parameters – the period k_{c} and the phase φ – forming a model space which must be scanned to determine if there is any significant correlation (in the absence of any physical motivation for restricting attention to specific periodicities). We have undertaken an initial survey of these models with the wavelength range defined by the native resolution of the present modal estimator (hybrid local polynomials with 600 modes), similar to the feature model search in WMAP data in Fergusson et al. (2012). For feature models of Eq. (16) we can obtain high correlations (above 95%) for the predicted CMB bispectrum if we take k_{c}> 0.01, that is, for an effective multipole periodicity ℓ_{c}> 140 feature models are accurately represented.
The results of a first survey of feature models in the Planck data is shown in Table 12 for 0.01 ≤ k_{c} ≤ 0.1 and phases φ = 0,π/ 4,π/ 2,3π/ 4 (for φ ≥ π we will identify a correlation with the opposite sign). Again, there was good consistency between the different foregroundseparation methods SMICA, NILC and SEVEM showing that the results are robust to potential residual foreground contamination in the data. For brevity we only give SMICA results here, while providing measurements from other component separation methods in Appendix C. Feature signals are typically largely uncorrelated with the ISWlensing or point sources, but nevertheless we subtract these signals and give results for the cleaned f_{NL}. The Table 12 results show that there is a parameter region around 0.01 ≤ k_{c} ≤ 0.025 for which signals well in excess of 2σ are possible (we undertook a broader search with 0.01 ≤ k_{c} ≤ 0.1 but found only a low signal beyond k> 0.3). It appears that some feature models are able to match the lowℓ “plusminus” and other features in the Planck bispectrum reconstruction (see Fig. 6). The best fit model has k_{c} = 0.0185 (ℓ_{c} ≈ 260) and phase φ = 0 with a signal −3σ. As a further validation step of our results, we also reanalysed the models with >2.5σ significance using a different modal decomposition, namely an oscillating Fourier basis (n_{max} = 300) augmented with a local SW mode (the same used for the reconstruction plots in Sect. 7). The results from this basis are shown in Appendix C and they are fully consistent with the polynomial measurements presented here. The previous bestfit WMAP feature model, k_{c} = 0.014 (ℓ_{c} ≈ 200) and phase φ = 3π/ 4, attained a 2.15σ signal with ℓ < 500 (Fergusson et al. 2012), but it only remains at this level for Planck.
We note however that the apparently high statistical significance of these results is much lower if we consider this to be a blind survey of feature models, because we are seeking several uncorrelated models simultaneously. Following what we did for our study of impact of foregrounds in Sect. 8, we considered a set of 200 realistic lensed FFP6 simulations, processed through the SMICA pipeline, and including realistic foreground residuals. If we use this accurate MC sample to search for the same grid of 52 feature models as in Table 12, we find a typical maximum signal of 2.23( ± 0.56)σ. Searching across all feature models (see below) studied here yields an expected maximum 2.37( ± 0.53)σ (whereas the survey for all 511 models from all paradigms investigated yielded 2.55( ± 0.52)σ). This means that our bestfit model from data has a statistical significance below 1.5σ above the maximum signal expectation from simulations, so we conclude that we have no significant detection of feature models from Planck data.
Fig. 13
CMB bispectrum shown for the bestfit feature model with an envelope with parameters k = 0.01875, phase φ = 0 and Δk = 0.045 (see Table 13). Compare with the Planck bispectrum reconstruction, Fig. 7. 
Feature models typically have a damping envelope representing the decay of the oscillations as the inflaton returns to its background slowroll evolution. Indeed, the feature envelope is a characteristic of the primordial mechanism producing the fluctuations, decaying as k increases for inflation while rising for contracting models like the ekpyrotic case (Chen 2011). We have made an initial survey to determine whether a decaying envelope improves the significance of any feature models. The envelope employed was a Gaussian centred at k_{c} = 0.045 with a falloff Δk = 0.015,0.02,0.025,0.03,0.035,0.04,0.045 and results for specific parameters are given in Table 13. The best fit model remains k = 0.01875 (ℓ_{c} = 265) with phase φ = 0 and the significance rises to 3.23σ, together with a second model k = 0.02 (ℓ_{c} = 285) φ = π/ 4. However the caveats about blind survey statistics previously noted also do not allow a claim of any detection in this case. A plot of the bestfit feature model with a decay envelope is shown in Fig. 13, for which the main features should be compared with those in Fig. 7. NonGaussian bispectrum signals from feature models typically produce counterparts in the power spectrum as will be described in Sect. 9. An improved statistical interpretation of the results presented in this section will be possible when this additional investigation will be completed.
We have also undertaken a survey of resonant models and the nonBunchDavies resonant models (or enfolded resonance models). With the modal estimator, we can achieve high accuracy for the predicted bispectrum for k_{c}> 0.001 (note that this has a different logarithmic dependence to feature models and a varying effective ℓ_{c}). For the resonance model shape of Eq. (18), we have not undertaken an extensive survey, except selecting a likely range for a high signal with periodicity comparable to the feature model, that is, with 0.25 <k_{c}< 0.5 and phases φ = 0,π/ 4,π/ 2,3π/ 4,π. However, no significant signal was found (all below 1σ), as can be verified in Table C.1 in Appendix C. For the enfolded resonance model shape of Eq. (18) , we have undertaken a preliminary search in the range 4 <k_{c}< 12 with the same phases. Again, no significant signal emerges from the Planck data, as shown in Table C.2 in Appendix C.
7.3.4. Directional dependence motivated by gauge fields
We have investigated whether there is significant NG from bispectrum shapes with nontrivial directional dependence (Eq. (19)), which are motivated by inflationary models with vector fields. Using the primordial modal estimator we obtained a good correlation with the L = 1 flattenedtype model, but the squeezed L = 2 model produced a relatively poor correlation of only 60%, given the complexity of the dominant squeezed limit. Preliminary constraints on these models are given in the Table 11, showing no evidence of a significant signal.
7.3.5. Warm inflation
Warm inflation produces a related shape with a sign change in the squeezed limit. This also had a poor correlation, until smoothing (WarmS) was applied as described in Fergusson et al. (2012). The resulting bispectrum shows no evidence for significant correlation with Planck data (SMICA), (89)The full list of constraints for SMICA, NILC and SEVEM models can be found for warm inflation and vector models in Table C.3 in Appendix C.
Feature model results with an envelope decay function of width Δk.
7.3.6. Quasisinglefield inflation
Finally, quasisinglefield inflation has been analysed constraining the bispectrum shape of Q (Eq. (12)), that depends on two parameters, ν and . In order to constrain this model we have calculated modal coefficients for 0 ≤ ν ≤ 1.5 in steps of 0.01 (so 151 models in total). These were then applied to the data and the one with the greatest significance was selected. Results are shown in Fig. 26. The maximum signal occurred at ν = 1.5, (0.31σ). To obtain error curves we performed a full likelihood using 2 billion simulations following the method described in Sefusatti et al. (2012). Such a large number of simulations was possible as they were generated from the modal βcovariance matrix which is calculated once from the 200 Planck realistic CMB simulations, rather than repeatedly from the CMB simulations themselves. The procedure is to take the 151 × 151 correlation matrix for the models (this is just the normalized dot product of the modal coefficients). This is then diagonalised using PCA, after which only the first 5 eigenvalues are kept as the remaining eigenvalues are <10^{10}. The βcovariance matrix is projected into the same subbasis where it is also diagonalised via PCA into 5 orthonormal modes, with the two leading modes closely correlated with local and equilateral. The procedure by which to produce a simulation is to generate five Gaussian random numbers and add the mean values obtained from the Planck data, rotating them to the subbasis where we determine the ν with the greatest significance. The result is then projected back to the original space to determine the related f_{NL}. The two billion results from this MC analysis are then converted into confidence curves plotted in Fig. 26. The curve shows that there is no preferred value for ν with all values allowed at 3σ. This reflects the results obtained from data previously, where we found the least preferred value of ν = 0.86 had only a marginally lower significance of 0.28σ (Sefusatti et al. 2012). Of course, these conclusions are directly related to the null results for both local and equilateral templates.
7.4. Constraints on local nonGaussianity with Minkowski Functionals
In this subsection, we present constraints on local NG obtained with Minkowski Functionals (MFs). MFs describe the morphological properties of the CMB field and can be used as generic estimators of NG (Komatsu et al. 2003; Eriksen et al. 2004; De Troia et al. 2007; Hikage et al. 2008; Curto et al. 2008; Natoli et al. 2010; Hikage & Matsubara 2012; Modest et al. 2013). As they are sensitive to every order of NG, they can be used to constrain different bispectrum and trispectrum shapes (Hikage et al. 2006, 2008; Hikage & Matsubara 2012). They are therefore complementary to, and a useful validation of, optimal estimators. Their precise definition and analytic formulations are presented in Planck Collaboration XXIII (2014). The MF technique is also used in the companion paper Planck Collaboration XXV (2014).
We review here the properties of MFs, as a complementary tool to polyspectrum based estimators.
Fig. 14
The Wiener filter W_{M} used to constrain with MFs. 
First, they are defined in real space, which makes MFs robust to masking effects and no linear term is needed to take into account the anisotropy of the data model. Second, as MFs are sensitive to every nonGaussian feature in the maps, they can be a useful probe of every potential bias in the bispectrum measurement, in particular the different astrophysical contaminations (foregrounds and secondaries).
There is a limitation to MF studies: they can be expressed in terms of weighted sums of the bispectrum (and trispectrum) in harmonic space (Matsubara 2010), hence the angledependence of the bispectrum is partially lost. This makes MFs suboptimal in two ways: increasing error bars for constraints on specific shapes and reducing the distinguishability of different bispectrum shapes. This lack of specificity can introduce biases, as MFs will partially confuse primordial and nonprimordial sources of NG and can introduce degeneracies between different primordial shapes. Constraints on orthogonal and equilateral shapes are quite degenerate with MFs, we therefore chose here to focus on the local bispectrum shape. We also leave trispectrum analyses for future studies.
An attractive feature of MFs is their linearity for weak NG (f_{NL}) and weak signals (such as point sources, and Galactic residuals after masking and component separation) (Ducout et al. 2013). This property can be used to estimate different known nonprimordial contributions.
7.4.1. Method
We constrain using the optimized procedure described in Ducout et al. (2013). To obtain constraints on , we apply a specific Wiener filter on the map (W_{M}), shown in Fig. 14. We do not use here the filter designed to enhance the information from the gradients of the map (), because this filter is very sensitive to smallscale effects and may be biased by foreground residuals.
Validation tests with MFs: results for obtained using the filter W_{M}, for ℓ_{max} = 2000 and N_{side} = 2048.
We use maps at HEALPix resolution N_{side} = 1024 (Górski et al. 2005) and ℓ_{max} = 2000. Our results are based on the four normalized^{14} functionals v_{k}(k = 0,3) (respectively Area, Perimeter, Genus and N_{cluster}), computed on n_{th} = 26 thresholds ν, between ν_{min} = − 3.5 and ν_{max} = + 3.5 in units of the standard deviation of the map.
We combine all functionals into one vector y (of size n = 104). We then analyse this vector in a Bayesian way to obtain a posterior for the , and hence an estimate of this parameter. The principle is to compare the vector measured on the data ŷ to the ones measured on nonGaussian simulations with the same systematic effects (realistic noise, effective beam) and data processing (Wiener filtering) as the data, . Modelling the MFs as multivariate Gaussians we obtain the posterior distribution for with a χ^{2} test: (90)with (91)Since NG is weak, the covariance matrix C is computed with 10^{4} Gaussian simulations, again reproducing effective beam, realistic noise and filtering of the data. The dependence of the MFs on , , is obtained as an average of ŷ measured on 100 simulations. The simulations used here are based on the WMAP7 bestfit power spectrum (Komatsu et al. 2011), using the procedure described in Elsner & Wandelt (2009).
7.4.2. Validation tests
We report here validation of the MFs estimator on in thoroughly realistic Planck simulations. This validation subsection is analogous to Sect. 6.1 concerning bispectrumbased estimators. The same tests (ideal Gaussian maps, fullsky nonGaussian maps with noise and nonGaussian maps with noise and mask) are performed, but different nonGaussian simulations are used. NonGaussian CMB simulations as processed in Fergusson et al. (2010a) only guarantee the correctness of the 3point correlations. Since the MFs are sensitive to higherorder npoint functions, they were validated with physical simulations (Elsner & Wandelt 2009).
The first test consists of 100 simulations of a fullsky Gaussian CMB, with a Gaussian beam with FWHM of 5 arcmin and without any noise, cut off at ℓ_{max} = 2000, with N_{side} = 2048. Here validation tests were made at N_{side} = 2048, but results (estimate and error bars) remain the same at N_{side} = 1024 as we keep the same ℓ_{max}. The second test includes nonGaussian simulations with and realistic coloured and anisotropic noise, processed through the Planck simulation pipeline and the componentseparation method SMICA. Finally, in the third test we add the union mask U73 to the previous simulations, masking both the Galaxy and the brightest point sources, and leaving 73% of the sky unmasked. Only the inpainting of the smallest holes in the point sources part of the mask was performed. For these three tests, the results are presented in Table 14. We give here the average estimate and error bar obtained on the 100 simulations, when we use the four functionals. The results show that the MF estimator is unbiased, robust, and a competitive alternative to bispectrumbased estimators. Moreover a mapbymap comparison of the results obtained on with KSW and MFs estimators showed a fair agreement between the two methods.
7.4.3. Results
For our analysis we considered a foregroundcleaned map obtained with the component separation method SMICA at N_{side} = 1024 and ℓ_{max} = 2000. As for the previous results in this section, we used the union mask U73, which leaves f_{sky} = 73% of the sky after masking Galaxy and point sources. To take into account some instrumental effects (asymmetry of beams, component separation processing) and known nonGaussian contributions (lensing), we used realistic unlensed and lensed simulations (10^{3}) of Planck data (FFP6 simulations, see Appendix A). First, MFs were applied to the unlensed simulations and the resulting curves served to calibrate the estimator, as the Gaussian part of the NG curves ^{15}. This estimate is referred to as the “raw map”. Secondly, MFs were applied to the lensed simulations, and the same procedure was applied, the result being referred as “lensingsubtracted”. We summarize the procedure in the following equation: (92)Here we assume that the MF respond linearly to lensing at first order and that primordial NG and lensing contributions are therefore additive.
Estimates of obtained with MFs on Planck data.
Fig. 15
MFs curves for SMICA at N_{side} = 1024 and ℓ_{max} = 2000, for the four functionals v_{k}: a) area, b) perimeter, c) genus, and d) N_{cluster}. The curves are the difference of each normalized MF, measured from the data, to the average from Gaussian Planck realistic simulations (not lensed). The difference curves are normalized by the maximum of the Gaussian curve. To compare the curves to the presence of primordial NG, the average difference curves for nonGaussian simulations with is also represented (100 simulations). 
Additionally, we tried to characterize other nonprimordial contributions that one could expect in masked SMICAcleaned maps covering 73% of the sky. To this end, we used simulations of extragalactic foregrounds and secondary anisotropies: uncorrelated (Poissonian) point sources; clustered CIB; and SZ clusters. These component simulations reproduce accurately the whole Planck data processing pipeline (beam asymmetry, component separation method). Using the linearity of MFs (Ducout et al. 2013), we could introduce these effects as a simple additive bias on the curves following (93)where is the average bias measured on 100 simulations. Note that the SZ simulation does not take into account the SZlensing correlation, which is expected to be negligible given the error bars.
Results are summarized in Table 15 and MFs curves are shown in Fig. 15, without including the lensing subtraction (“raw curves”). Considering the larger error bars of MF estimators, the constraints obtained are consistent with those from the bispectrumbased estimators, even without subtracting the expected nonprimordial contributions. Moreover, results are quite robust to Galactic residuals: constraints obtained with other component separation methods (NILC and SEVEM), with different sky coverage, differ from the SMICA results presented here by less than .
Fig. 16
The two upper maps show the modulation reconstruction mean field at L ≤ 100, which is essentially a map of the expected total smallscale power on the masked map as a function of position (assuming there is no primordial power modulation). The top mean field map from the 143 GHz auto estimator has a large signal from both the cut (which can be calculated accurately), and from the noise anisotropy (aligned roughly with the ecliptic, which cannot be estimated very accurately from simulations). The lower mean field is the 143 × 217 GHz crossestimator map, and does not have the contribution from the noise anisotropy (note the colour scale is adjusted). The lower plot shows the corresponding mean field power spectra compared to the reconstruction noise (connected part of the trispectrum); the reconstruction noise is much smaller than both the detector noise and mask contributions to the mean field. Since any τ_{NL} signal is all on large scales we do not reconstruct modes above L_{max} = 100. 
Fig. 17
Power spectrum of the power modulation reconstructed from 143 × 217 GHz maps. Shading shows the 68%, 95% and 99% CL intervals from simulations with no modulation or kinematic signal. The dashed lines are when the mean field simulations include no kinematic effects, showing a clear detection of a modulation dipole. The blue points show the expected kinematic modulation dipole signal from simulations, along with 1σ error bars (only first four points shown for clarity). The solid line subtracts the dipolar kinematic signal in the mean fields from simulations including the expected signal, and represents out best estimate of the nonkinematic signal (note this is not just a subtraction of the power spectra since the mean field takes out the fixed dipole anisotropy in real space before calculating the remaining modulation power). The dotted line shows the expected signal for τ_{NL} = 1000. 
Fig. 18
Distribution of estimators from Gaussian simulations (L_{max} = 10) compared to data estimates (vertical lines). The distribution is rather skewed because the main contributions are from L ≲ 4 where the modulation power spectra have skewed χ^{2} distributions with low degrees of freedom. The red line shows the predicted distribution for a weighted sum of estimators assuming the reconstructed modulation modes are Gaussian with 2L + 1 modes measured per L, which fits the full simulations well. The black vertical lines show the data estimates from L_{max} = 10, and should be compared to the green which have L_{max} = 2 and hence are insensitive to the anomalous octopole signal. The dashed lines are τ_{NL,1}, the slightly more optimal variant of the estimator. 
7.5. CMB trispectrum results
As shown in Fig. 16 the modulation reconstruction mean field has two large contributions, one from the mask and one from anisotropic noise, reflecting the fact that they both look like a large spatiallyvarying modulation of the fluctuation power. The noise we use to estimate the mean field is taken from FFP6 simulations, adjusted with an additional 10 μK arcmin white noise component to match the power spectrum in the observed maps. However this is still only an approximate description of the instrumental noise present in the data. The mean field from nonindependent maps (e.g., 143 × 143 GHz maps) shows a large noise anisotropy that is primarily quadrupolar before masking, and any mismatch between the simulated noise and reality would lead to a large error in the mean field subtraction. By instead using the modulation estimator for 143 × 217 GHz maps errors due to misestimation of the noise are avoided, and the mean field is then dominated by the shape of the Galactic cut, which is well known, and a smaller uncertainty from assumed simulation power spectrum and beam errors (see Fig. 16). For this reason for our main result we work with modulation reconstructions generated from 143 × 217 GHz maps with independent noise, which removes the leading error due to noise mean field misestimation.
Figure 17 shows the reconstructed modulation power from 143 × 217 GHz maps that we use for our analysis. We show two results: one where we do not include the expected kinematic dipole signal in the mean field subtraction (see Sect. 5.4), and one were we do so that the reconstruction should then be dominated by the primordial and any unmodelled systematic effects. In the first case the 143 × 217 result gives a clear first detection of the dipolar kinematic modulation signal of roughly the expected magnitude (see Planck Collaboration XXVII (2014) for a more detailed discussion of this signal). Including the expected kinematic signal in the simulations (and hence the mean field) removes this signal, giving a cosmological modulation reconstruction that is broadly consistent with no modulation (statistical isotropy) except for the anomalous very significant signal in the modulation octopole.
Note that only the twopoint reconstruction is free from noise bias, the fourpoint is still sensitive to noise modelling at the level of the subtraction of the reconstruction noise power spectrum (Eq. (74)). However as shown in the Fig. 17, is not that much larger than the reconstruction scatter at low multipoles where the τ_{NL} signal peaks, so the sensitivity to noise misestimation is much less than in the mean field subtraction (where the largescale noise anisotropy gives a largescale mean field in the autoestimators orders of magnitude larger, Fig. 16).
The τ_{NL} estimator from the 143 × 217 GHz modulation reconstruction gives , compared to a null hypothesis distribution at 95% CL (σ_{τNL} ≈ 335). Our quoted error bar is assuming zero signal so that there is no signal cosmic variance contribution, and the bulk of the apparent signal is coming from the high octopole seen in Fig. 17. The alternative estimator gives a slightly different weighting to the octopole, giving with an expected nullhypothesis σ_{τNL} ≈ 332. The surprisingly large difference between the estimators can be explained as due to the large octopole signal, which has . However the shape of the total signal would not be expected from a genuine τ_{NL} signal, since as shown in Fig. 17 on average the latter is expected to fall off approximately proportional to 1 /L^{2} (i.e., a large primordial τ_{NL} would be expected in most realisations to give large dipole and quadrupole signals that we do not see). If we estimate τ_{NL} using L_{max} = 2 we obtain with only a slightly larger nullhypothesis error σ_{τNL} = 349, where in this case .
Note that the distribution of is quite skewed because the signal is dominated by the verylow multipole modulation power spectra which have skewed χ^{2}like distributions due to the large cosmic variance there (see Fig. 18; Hanson & Lewis 2009; Smith & Kamionkowski 2012). The reconstruction noise acts like nearlyuncorrelated Gaussian white noise, so each of the comes from a sum of squares of ~2L + 1 modulation reconstruction modes; the shape of the distribution is consistent with what would come from calculating a weighted sum of χ^{2}distributed random variables. If we assume that any primordial modulation modes giving rise to a physical τ_{NL} signal are also Gaussian, for any given physical τ_{NL} the estimators would also have χ^{2} distributions. This allows us to evaluate the posterior distribution of τ_{NL} given the observed , in exactly the same way that one can do for the CMB temperature power spectrum. For each L the posterior distribution on the full sky would have an inversegamma distribution. We follow Hamimeche & Lewis (2008) by generalizing this to a cutsky approximation for a range of multipoles: (94)
where M_{LL′} is the covariance of the estimators calculated from null hypothesis simulations, is the estimator reconstruction noise, (95)and (96)For uncorrelated χ^{2}distributed estimators this distribution reduces to the exact distribution, a product of inversegamma distributions. This approximation to the posterior can be used to evaluate the probability of any scale dependent τ_{NL}(L), and does not rely on compression into a single estimator; it can therefore fully account for the observed distribution of modulation power between L. Here we focus on the main case of interest where τ_{NL}(L) is nearlyscaleinvariant so that for all L we have τ_{NL}(L) = τ_{NL}.
Fig. 19
Approximate posterior distributions for a range of L_{max}. The distributions have broad tails to high values because of the small number of largescale modulation modes that are measured, and hence large cosmic variance. For L_{max} ≥ 3 the distributions are pulled away from zero by the significant octopole modulation signal observed, and are gradually move back towards zero as we include more modulation modes that are inconsistent with large τ_{NL} values. As shown in Fig. 20 the octopole has significant frequency dependence and is therefore unlikely to be physical. 
Fig. 20
Comparison of the unnormalized modulation power with various combinations of frequencies. The middle plot shows the results used for our main τ_{NL} result, since it removes the significant largelyquadrupolar signal from anisotropic noise misestimation seen in the two other plots. The noticeable difference in the odd octopole signal between channels indicates that the residual signal in 143 × 217 GHz is unlikely to be physical, but we cannot currently identify its origin. 
The resulting posterior distributions of τ_{NL} are shown in Fig. 19 for a range of L_{max}. These are strongly skewed, in the same way as the posterior from the low quadrupole in the CMB temperature data. The high octopole is pulling the distributions up to higher τ_{NL} values, but increasing L_{max} can reduce the highL tail because very large τ_{NL} values are inconsistent with the low modulation power seen at L ≠ 3. With L_{max} = 2 the posterior peaks near zero, but the distribution is then very broad because there are only about 8 modes, which therefore have large cosmic variance. For L_{max} = 50 we find τ_{NL}< 2800 at 95% CL, which we take as our upper limit.
Figure 20 shows the modulation reconstructions for the 143 GHz and 217 GHz maps separately compared to the cross estimator. The picture is complicated here by the large signals from noise misestimation seen in the 143 × 143 and 217 × 217 estimators, however the fact that the octopole in 143 × 143 is lower than in the crossestimator indicates that the octopole signal is very unlikely to be mainly physical. Our measured τ_{NL} limit in practice represents a strong upper limit on the level of primordial τ_{NL} that could be present, since unmodelled varying smallscale foreground or nonconstant gain/calibration would also only serve to increase the measured estimate compared to primordial on average. The octopole signal does vary slightly with Galactic mask, though at present we cannot clearly isolate its origin. If more extensive analysis (for example using crossmap estimators at the same frequency) can identify a nonphysical origin and remove it, the quoted upper limit on τ_{NL} would become significantly tighter. For a more extensive discussion of possible foreground and systematic effect issues that can affect 4point estimators see Planck Collaboration XXVII (2014) and Planck Collaboration XVII (2014).
Fig. 21
Dependence of the dipole modulation amplitude as a function of ℓ_{max}. Upper panel: the amplitude of the dipole of the reconstructed power modulation from 143 × 217 GHz maps with the kinematic dipole subtracted; the dashed line shows that the observed signal decreases with ℓ_{max} in the same way as the rms value expected from simulations. Lower panel: corresponding confidence values for the observed dipole of the modulation power spectrum assuming it followed a chisquared distribution; this shows the anomalouslooking results for ℓ_{max} ~ 60 and ℓ_{max} ~ 600 consistent with Planck Collaboration XXVII (2014), but that on smaller scales the observed signal is consistent with isotropy as expected from the scaleinvariant constraint using ℓ_{max} = 2000 shown in Fig. 17. 
We have focussed on nearly scale invariant modulations here. As discussed in Planck Collaboration XXIII (2014) there are some potentially interesting “anomalies” in the distribution of power in the Planck data, especially the hemispherical power asymmetry at ℓ_{max} ≤ 600. If these power asymmetries are physical rather than statistical flukes, they must be strongly scaledependent, and would correspond to a scaledependent τ_{NL}like trispectrum. As we have shown here the dipole power asymmetry does not persist to small scales, except for the signal aligned with the CMB dipole expected from kinematic effects: from our analysis using ℓ_{max} = 2000 we find a primordial consistent with Gaussian simulations, corresponding to a dipole modulation amplitude  f  ~ 0.2%. To be consistent with a ~7% modulation at ℓ_{max} ~ 60 and ~1% modulation at ℓ_{max} 600, as seen in the lowerℓ anisotropies of both WMAP and Planck the primordial trispectrum would have to be strongly scaledependent on smallscale modes, so that larger smallscale modes are modulated more than the smallest ones (rather than just τ_{NL} = τ_{NL}(L)). Figure 21 shows how the allowed dipole modulation amplitude drops as ℓ_{max} increases (similar to Fig. 3 of Hanson & Lewis (2009) for WMAP but now extending to the higher ℓ_{max} available from Planck). This result is consistent with Planck Collaboration XXIII (2014), where different estimators and analysis cuts also find no evidence of primordial dipolelike power asymmetry on small scales, but confirm the largescale distribution of power seen with WMAP and also a marginallysignificant smalleramplitude (f ~ 1%) dipolelike modulation at ℓ_{max} ~ 600.
8. Validation of Planck results
Fig. 22
Evolution of the f_{NL} parameters (solid blue line with data points) and their uncertainties (dashed lines) for the five bispectrum templates as a function of the maximum multipole number ℓ_{max} used in the analysis. From left to right and top to bottom the figures show respectively local, equilateral, orthogonal, diffuse point sources (all four with the ISWlensing bias subtracted), ISWlensing and local again (the last two without subtracting the ISWlensing bias). To show better the evolution of the uncertainties, they are also plotted around the final value of f_{NL} (solid green lines without data points). The results are for SMICA, assume all shapes to be independent, and have been determined with the binned bispectrum estimator. 
Results for f_{NL} (assumed independent, without any correction for the ISWlensing bias) of the SMICA cleaned map using different values of ℓ_{max}, for the KSW and binned estimators.
Here we perform a set of tests to check the robustness and stability of our f_{NL} measurements. As these are validation tests of Planck results, and not internal comparisons of bispectrum pipelines (already shown to be in tight agreement in Sects. 6 and 7) we will not employ all the bispectrum estimators on each test. In general we choose to use two estimators on each test, in order to have a crosscheck of the outcomes without excessive redundancy.
8.1. Dependence on maximum multipole number
The dependence on the maximum multipole number ℓ_{max} of the SMICA results (assuming independent shapes) is shown in Fig. 22 (for the binned estimator) and Table 16 (for both the KSW and binned estimators). Testing the ℓ_{max} dependence is easiest for the binned estimator, where one can simply omit the highest bins in the final sum when computing f_{NL}. It is clear that we have reached convergence both for the values of f_{NL} and for their error bars at ℓ_{max} = 2500, with the possible exception of the error bars of the diffuse pointsource bispectrum. The diffuse pointsource bispectrum template is dominated by equilateral configurations at high ℓ. Moreover, for all the shapes except point sources, results at ℓ_{max} = 2000 are very close to those at ℓ_{max} = 2500, taking into account the size of the error bars.
It is very interesting to see that at ℓ_{max} ~ 500 we find a local f_{NL} result in very good agreement with the WMAP9 value of 39.8 ± 19.9 (Bennett et al. 2013) (or 37.2 ± 19.9 after ISWlensing bias subtraction). At these low ℓ_{max} values we also find negative values for orthogonal f_{NL}, although not as large or significant as the WMAP9 value (which is −245 ± 100). One can clearly see the importance of the higher resolution of Planck both for the values of the different f_{NL} parameters and for their error bars.
It is also clear that the higher resolution of Planck is absolutely crucial for the ISWlensing bispectrum; this is simply undetectable at WMAP resolution. On the other hand, the high sensitivity of Planck measurements also exposes us to a larger number of potentially spurious effects. For example we see that the bispectrum of point sources is also detected at high significance by Planck at ℓ_{max} ≥ 2000, while remaining undetectable at lower resolutions. The presence of this bipectrum in the data could in principle contaminate our primordial f_{NL} measurements. For this reason, the presence of a large point source signal has been accounted for in previous sections by always including the Poisson bispectrum in a joint fit with primordial shapes. Fortunately, it turns out that the very low correlation between the primordial templates and the Poisson one makes the latter a negligible contaminant for f_{NL}, even when the residual point source amplitude is large.
8.2. Dependence on mask and consistency between frequency channels
Results for f_{NL} (assumed independent) for the raw frequency maps at 70, 100, 143, and 217 GHz with a very large mask (f_{sky} = 0.32) compared to the SMICA result with the union mask U73 (f_{sky} = 0.73), as determined by the binned (with ℓ_{max} = 2500) and modal (with ℓ_{max} = 2000) estimators.
To test the dependence on the mask, we have analysed the SMICA maps applying four different masks. Firstly the union mask U73 used for the final results in Sect. 7, which leaves 73% of the sky unmasked. Secondly we used the confidence mask CSSMICA89 of the SMICA technique, which leaves 89% of the sky. Next, a bigger mask constructed by multiplying the union mask U73 with the Planck Galactic mask CG60, leading to a mask that leaves 56% of the sky. And finally a very large mask, leaving only 32% of the sky, which is the union of the mask CL31 – used for power spectrum estimation on the raw frequency maps – with the union mask U73 (for mask details see Planck Collaboration XII 2014 for U73, CSSMICA89, and CG60; Planck Collaboration XV 2014 for CL31). The results of this analysis are presented in Table 17 for two different estimators: binned and modal. The f_{NL} are assumed independent here. In order to correctly interpret our results and conclusions, an important point to note is that binned results have been obtained choosing ℓ_{max} = 2500, while modal results correspond to ℓ_{max} = 2000. Primordial shape and ISWlensing results and error bars saturate at ℓ_{max} = 2000 (see Sect. 8.1), so the results from the two estimators are directly comparable in this case. The Poisson (point sources) bispectrum is however dominated by highℓ equilateral configurations and the signal for this specific template still changes from ℓ = 2000 to ℓ = 2500. The differences in central values and uncertainties between the two estimators for the Poisson shape are fully consistent with the different ℓ_{max} values. Direct comparisons on data and simulations between these two estimators and the KSW estimator showed that Poisson bispectrum results match each other very well when the same ℓ_{max} is used.
Results from the modal pipeline have uncertainties determined from MC simulations, while the results from the binned pipeline (in Table 17 and the next only) are given with Fisher error bars. It is very interesting to see that even with the large f_{sky} = 0.32 mask, the simple inpainting technique still allows us to saturate the (Gaussian) CramérRao bound, except for the ISWlensing shape where we have a significant detection of NG in a squeezed configuration (so that an error estimate assuming Gaussianity is not good enough). Finally we note that only for the tests in this and in the next paragraph we adopted a faster but slightly less accurate version of the modal estimator than the one used to obtain the final f_{NL} constraints in Sect. 7. In this faster implementation we use fewer modes in order to increase computational speed, and consequently we get a slight degrading of the level of correlation of our expanded templates with the initial primordial shapes. Note that the changes are small: we go from 99% correlation for local, equilateral and orthogonal shapes in the most accurate (and slower) implementation to 98% correlation for equilateral and orthogonal snapes, and 95% correlation for local shape in the faster version. This of course still allows for very stringent validation tests for all the primordial shapes, and produces results very close to the highaccuracy pipeline, while at the same time increasing overall speed by almost a factor 2. Both versions of the modal pipeline were separately validated on simulations (see Sect. 6).
Besides confirming again the good level of agreement between the two estimators already discussed in Sects. 6 and 7, the main conclusion we draw from this analysis is that our measurements for all shapes are robust to changes in sky coverage, taking into account the error bars and significance levels, at least starting from a certain minimal mask. The f_{sky} = 0.89 mask is probably a bit too small, likely leaving foreground contamination around the edges of the mask, though even for this mask the results are consistent within 1σ, except for point sources (which might suggest the presence of residual Galactic point source contamination for the small mask). The results from the f_{sky} = 0.73 and f_{sky} = 0.56 masks are highly consistent. This conclusion does not really change when going down to f_{sky} = 0.32, although uncertainties of course start increasing significantly for this large mask.
We also investigate if there is consistency between frequency channels when the largest mask with f_{sky} = 0.32 is used, and if these results agree with the SMICA results obtained with the common mask. The results (assuming independent f_{NL}) are given both for the binned and the modal estimator in Table 18. As in the previous table, the full modal pipeline (faster but slightly less accurate version) has been run here, obtaining both central values from data and MC error bars from simulations, while the binned pipeline (which is slower in determining full error bars than the modal pipeline) is used to crosscheck the modal measurements and has error bars given by simple Fisher matrix estimates. As one can see here and as was also checked explicitly in many other cases, the error bars from different estimators are perfectly consistent with each other and saturate the CramérRao bound (except in the case of a significant nonGaussian ISWlensing detection).
A detailed analysis of Table 18 might actually suggest that the agreement between the two estimators employed for this test, although still clearly good, is slightly degraded when compared to their performance on clean maps from different component separation pipelines. If we compare e.g., SMICA results in Table 17 to raw data results in Table 18, we see that in the former case the discrepancy between the two estimators is at most of order σ_{fNL}/ 3, and smaller in most cases. In the latter case, however, we notice several measurements displaying differences of order σ_{fNL}/ 2 between the two pipelines, and the value of at 70 GHz being 1σ away. We explain these larger differences as follows. For SMICA runs we calibrated the estimator linear terms using FFP6 simulations, accurately reproducing noise properties and correlations (see Appendix A). On the other hand, for the present tests on raw frequency channels we adopted a simple noise model, based on generating uncorrelated noise multipoles with a power spectrum as extracted from the halfring null map, and remodulating the noise in pixel space according to the hitcount distribution. This approximation is expected to degrade the accuracy of the linear term calibration, and thus to produce a slightly lower agreement of different pipelines for shapes where the linear correction is most important. Those are the shapes that take significant contributions from squeezed triangles: local and ISWlensing, and to a smaller but nonnegligible extent orthogonal, i.e., exactly the shapes for which we find slightly larger differences.
We conclude that no significant fluctuations are observed when comparing measurements from different frequency channels (between themselves or with the clean and coadded SMICA map) or from different estimators on a given channel for the primordial shapes. The same is true for the ISWlensing shape, although it should be noted that in particular the 70 GHz channel (like WMAP) does not have sufficient resolution to measure either the lensing or point source contributions. The uncertainties of the point source contribution vary significantly between frequency channels, although results remain consistent between channels given the error bars (when all multipoles up to ℓ_{max} = 2500 are taken into account, as performed by the binned estimator). This is due to the fact that this shape is dominated by highℓ equilateral configurations, the signaltonoise of which depends crucially on the beam and noise characteristics, which vary from channel to channel. In the SMICA map point sources are partially removed by foreground cleaning, explaining the significantly lower value than for 217 GHz. As explained before, differences between the binned and modal estimators regarding point sources are due to the different values of ℓ_{max} (2500 for binned and 2000 for modal), which particularly affects the 217 GHz channel and the SMICA cleaned map.
8.3. Null tests
Results for f_{NL} (assumed independent) of the SMICA halfring null maps, determined by the KSW, binned and modal estimators.
Results for f_{NL} (assumed independent) of several null maps determined by the binned estimator.
To make sure there is no hidden NG in the instrumental noise, we performed a set of tests on null maps. These are noiseonly maps obtained from differences between maps with the same sky signal. In the first place we constructed halfring null maps, i.e., maps constructed by taking the difference between the first and second halves of each pointing period, divided by 2. Secondly, we constructed a survey difference map (Survey 2 minus Survey 1 divided by 2). A “survey” is defined as half a year of data, roughly the time needed to scan the full sky once; the nominal period of Planck data described by these papers contains two full surveys. Finally we constructed the detector set difference map (“detset 1” minus “detset 2” divided by 2). The four polarized detectors at each frequency from 100 to 353 GHz are split into two detector sets per frequency, in such a way that each set can measure all polarizations and the detectors in a set are aligned in the focal plane (see Planck Collaboration VI 2014; and Planck Collaboration XII 2014 for details on the null maps analysed in this section).
All these maps are analysed using the union mask U73 used for the final data results. However, in the case of the survey and detector set difference maps this mask needs to be increased by the unseen pixels. That effect only concerns a few additional pixels for the detector set null map, but is particularly important for the survey difference map, since a survey only approximately covers the full sky. The final f_{sky} of the mask used for the survey difference map is 64%.
The test consisted of extracting f_{NL} from the null maps described above, using only the cubic part of the bispectrum estimators (i.e., no linear term correction), and keeping the same weights as for the full “signal + noise” analysis. This means that the weights were not optimized for noiseonly maps, as our aim was not to study the bispectrum of the noise per se but rather to check whether the noise alone produces a threepoint function detectable by our estimators when they are run in the same configuration as for the actual CMB data analysis. For a similar reason it would have been pointless to introduce a linear term in this test. The purpose of the linear correction is in fact that of decreasing the error bars by accounting for offdiagonal covariance terms introduced by sky cuts and noise correlations when optimal weights are used, which is not the case here.
Our f_{NL} error bars for this test are obtained by running the estimators’ cubic part on Gaussian noise simulations including realistic correlation properties. In the light of the above paragraph it is clear that such uncertainties have nothing to do with the actual uncertainties from CMB data, and cannot be compared to them.
Since SMICA was the main componentseparation method for our final analysis of data, we present in Table 19 the results of our SMICA halfring study using the KSW, binned and modal estimators, i.e., all the three main pipelines used in this paper. For the cleaned maps we do not have survey or detector set difference maps. Those are, however, available for single frequency channels. Thus we also studied all three types of null maps for the raw 143 GHz channel in Table 20, using the binned estimator. In both tables all f_{NL} shapes are assumed to be independent. The binned estimator is best suited for these specific tests as its cubic part is less sensitive to masking compared to other pipelines, especially modal. Therefore in this “cubic only” test, the binned results provide the most stringent constraints in terms of final error bars.
As one can see Planck passes these null tests without any problems: all values found for f_{NL} in these null maps are completely negligible compared to the final measured results on the data maps, and consistent with zero within the error bars.
8.4. Impact of foreground residuals
In Sect. 7 we applied our bispectrum estimators to Planck data filtered through four different component separation methods: SMICA, NILC, SEVEM and CR (for a detailed description of component separation techniques used for Planck see Planck Collaboration XII 2014). The resulting set of f_{NL} measurements shows very good internal consistency both between different estimators (as expected from our MC validation tests of bispectrum pipelines described in Sect. 6) and between different foregroundcleaned maps. This already makes it clear that foreground residuals in the data are very well under control, and their impact on the final f_{NL} results is only at the level of a small fraction of the measured error bars. In this section we further investigate this issue, and validate our previous findings on data by running extensive tests in which we compare simulated data sets with and without foreground residuals from two different component separation pipelines, SMICA and NILC. The goal is to provide a MCbased assessment of the expected f_{NL} systematic error from residual foreground contamination.
For each component separation pipeline, we consider two sets of simulations. One set includes realistic Planck noise and beam, is masked and inpainted in the same way as we do for real data, and is processed through SMICA and NILC but it does not contain any foreground component. The other set is obtained by adding to the first one a number of diffuse foreground residuals: thermal and spinning dust components; freefree and synchrotron emission; kinetic and thermal SZ; CO lines and correlated CIB. These residuals have been evaluated by applying the component separation pipelines to accurate synthetic Planck datasets including foreground emission according to the PSM (Delabrouille et al. 2013), and are of course dependent on the cleaning method adopted. The simple procedure of adding foreground residuals to the initially clean simulations is made possible because we consider only linear component separation methods for our analysis. Linearity is in general an important requirement for foreground cleaning algorithms aiming at producing maps suitable for NG analyses. All maps in both samples are lensed using the LensPix algorithm. We analyse both sets using different bispectrum estimators (modal, KSW, binned) for crossvalidation purposes.
The presence of residual foreground components in the data can have two main effects on the measured f_{NL}. The first is to introduce a bias in the f_{NL} measurements due to the correlation between the foreground and primordial 3point function for a given shape. The second is to increase the error bars while leaving the bispectrum estimator asymptotically unbiased. This is a consequence of accidental correlations between primordial CMB anisotropies and foreground emission. Of course these “CMBCMBforegrounds” bispectrum terms average to zero but they do not cancel in the bispectrum variance 6point function. An interesting point to note is that a large foreground threepoint function will tend to produce a negative bias in the local bispectrum measurements. That is because foreground emission produces a positive skewness of the CMB temperature distribution (“excess of photons”), and a positive skewness is in turn related to a negative ^{16}. A large negative is thus a signature of significant foreground contamination in the map. This is indeed what we observe if we consider raw frequency maps with a small Galactic cut, which is why our frequencybyfrequency analysis in Sect. 8.2 was performed using a very large mask. For more details regarding effects of foreground contamination on primordial NG measurements in the context of the WMAP analysis see Yadav & Wandelt (2008); and Senatore et al. (2010).
Summary of our f_{NL} analysis of foreground residuals.
In our test we built maps contaminated with foreground residuals by simply adding residual components to the clean maps. That means that the difference for the ith realization in the simulated sample exactly quantifies the change in f_{NL} due to foregrounds on that realization. In order to assess the level of residual foreground contamination on primordial and ISWlensing shapes, first of all we consider the difference between average values and standard deviations of f_{NL} measured from the two map samples for each shape. As shown in Table 21 we find that neither the average nor the standard deviation shows a significant change between the two datasets. That means that foreground residuals are clearly subdominant, as they do not bias the estimator for any shape and they do not increase the variance through spurious correlations with the CMB primordial signal.
We also consider the difference on a mapbymap basis and compute its standard deviation. This is used as an estimate of the expected bias on a single realization due to the presence of residuals. As already expected from the negligible change in the standard deviations of the two samples, the variance of the mapbymap differences is also very small: Table 21 again shows that it is at most about σ_{fNL}/ 6 for any given shape, where σ_{fNL} is the f_{NL} standard deviation for that shape. As an example, in Fig. 23 we show the measured values of for the first 99 maps in both the SMICA and NILC samples, comparing results with and without residuals. It is evident also from this plot that the change in central value due to including residuals is very small. The very good agreement between the two component separation pipelines is also worth notice.
From the study shown here and from the comparison between different component separation methods in Sect. 7, we can thus conclude that the combination of foregroundcleaned maps and f_{sky} = 0.73 sky coverage we adopt in this work provide a very robust choice for f_{NL} studies.
Fig. 23
The measured for the first 99 maps in the lensed FFP6 simulation sample used for the foreground studies presented in Sect. 8.4. We show measurements from SMICA and NILC processed maps both with and without foreground residuals. The horizontal solid line is the average value of the SMICA clean maps, and the dashed and dotted horizontal lines correspond to 1σ and 2σ deviations, respectively. The plot clearly shows the very small impact of including residuals, and the very good consistency between the two component separation pipelines. 
8.5. Impact of HFI timeordered information processing on NG constraints
Our validation tests are designed to determine whether instrumental or data processing systematics could lead to a spurious detection of primordial NG. Our analysis protocol passes these tests which allows us to rule out this concern with confidence. These validation tests do not exclude the possibility that the extensive processing of HFI timeordered information (TOI) could somehow remove nonGaussian signals present on the sky and thus negatively impact Planck’s ability to detect them. If this were true it would weaken the bounds on primordial NG reported here. HFI TOI processing includes glitch fitting, gain drift and correction (ADC nonlinearity), 4K cooler line removal, telegraph noise step correction, and TOD filtering to correct for the bolometer time constant (Planck Collaboration VI 2014; Planck Collaboration VII 2014; Planck Collaboration X 2014).
The following facts constitute strong evidence that HFI TOI processing does not remove nonGaussian signals:

1.
Planck’s 2.6σ ISWlensing bispectrum measurement isconsistent with expectations from the LCDM model, and has theright skewC_{ℓ} shape. Like NG of local type, this is a squeezedbispectrum template.

2.
When frequency maps are combined to maximize the CIB signal, Planck finds a 25σ detection of the nearly squeezed CIBlensing bispectrum, consistent with the CIB 2point correlations inferred by Planck (Planck Collaboration XXX 2014).

3.
Planck detects residual point sources in the SMICA maps at 5σ, seen as a bispectrum of equilateral shape at high ℓ consistent with the expected skewspectrum shape (Fig. 2) on the angular scales where residual infrared sources and far infrared background are expected to be the dominant contaminants of the power spectrum according to the foreground residuals in FFP6 simulations of SMICA maps (Fig. E.3 in Planck Collaboration XII 2014).

4.
The trispectrum signal imprinted by Planck’s motion with respect to the CMB rest frame is seen at a sensible level and in a plausible direction (Planck Collaboration XXVII 2014).

5.
Planck detects the trispectrum signal due to lensing by large scale structure with high significance and in accordance with LCDM predictions based on the Planck parameter likelihood.

6.
Figure 22 shows that we reproduce the WMAP9 local NG results on those larger angular scales where the WMAP9 data are cosmic variance limited.
We conclude that all the forms of NG that should have been seen by Planck have been seen, in quantitative agreement with the expected level across the entire range of angular scales probed. While it cannot be excluded with absolute certainty that some unknown aspect of HFI TOI processing could have affected some unknown form of NG, the presence of these nonGaussian features in the resulting maps (in addition to signals such as the extracted map of ydistortion, and galactic and extragalactic foregrounds), gives us confidence in the force of the bounds on primordial NG described in this paper.
9. Implications for early Universe physics
The NG analyses performed in this paper show that Planck data are consistent with Gaussian primordial fluctuations. The standard models of singlefield slowroll inflation have therefore survived the most stringent tests of Gaussianity performed to date. On the other hand, the NG constraints obtained on different primordial bispectrum shapes (e.g., local, equilateral and orthogonal), after properly accounting for various contaminants, severely limit various classes of mechanisms for the generation of the primordial perturbations proposed as alternatives to the standard models of inflation.
In the following subsections, unless explicitly stated otherwise, we translate limits on NG to limits on parameters of the theories by constructing a posterior assuming the following: the sampling distribution is Gaussian (which is supported by Gaussian simulations); the likelihood is approximated by the sampling distribution but centred on the NG estimate (see Elsner & Wandelt 2009); uniform or Jeffreys’ prior (where stated), over ranges which are physically meaningful, or as otherwise stated. Where two parameters are involved, the posterior is marginalized to give onedimensional constraints on the parameter of interest.
9.1. General singlefield models of inflation
DBI models: the constraints on provide constraints on the sound speed with which the inflaton fluctuations can propagate during inflation. For example, for DBI models of inflation (Silverstein & Tong 2004; Alishahiha et al. 2004), where the inflaton field features a nonstandard kinetic term, the predicted value of the nonlinearity parameter is (Silverstein & Tong 2004; Alishahiha et al. 2004; Chen et al. 2007b). Although very similar to the equilateral shape, we have obtained constraints directly on the theoretical (nonseparable) predicted shape (Eq. 7)). The constraint at 68% CL (see Eq. (88)) implies (97)The DBI class contains two possibilities based on string constructions. In ultraviolet (UV) DBI models, the inflaton field moves under a quadratic potential from the UV side of a warped background to the infrared side. It is known that this case is already at odds with observations, if theoretical internal consistency of the model and constraints on power spectra and primordial NG are taken into account (Baumann & McAllister 2007; Lidsey & Huston 2007; Bean et al. 2007; Peiris et al. 2007). Our results strongly limit the relativistic régime of these models even without applying the theoretical consistency constraints.
It is therefore interesting to look at infrared (IR) DBI models (Chen 2005b,a) where the inflaton field moves from the IR to the UV side, and the inflaton potential is , with a wide range 0.1 <β< 10^{9} allowed in principle. In previous work (Bean et al. 2008) a 95% CL limit of β< 3.7 was obtained using WMAP. In a minimal version of IR DBI, where stringy effects are neglected and the usual field theory computation of the primordial curvature perturbation holds, one finds (Chen 2005c; Chen et al. 2007b) c_{s} ≃ (βN/ 3)^{1}, n_{s} − 1 = − 4 /N, where N is the number of efolds; further, primordial NG of the equilateral type is generated with an amplitude . For this model, the range N ≥ 60 is compatible with Planck’s 3σ limits on n_{s} (Planck Collaboration XXII 2014). Marginalizing over 60 ≤ N ≤ 90, we find (98)dramatically restricting the allowed parameter space of this model.
Powerlaw kinflation: these models (ArmendarizPicon et al. 1999; Chen et al. 2007b) predict , where n_{s} − 1 = − 3γ, . Assuming a prior of 0 <γ< 2/3, our constraint at 68% CL (see Table 9) leads to the limit γ ≥ 0.05 at 95% CL. On the other hand, Planck’s constraint on n_{s} − 1 yields the limit 0.01 ≤ γ ≤ 0.02 (Planck Collaboration XXII 2014). These conflicting limits severely constrain this class of models.
Flat Models and higher derivative interactions: flat NG can characterize inflationary models which arise from independent interaction terms different from the and discussed in Sect. 2 (see also Sect. 9.2). An example is given by models of inflation based on a Galilean symmetry (Creminelli et al. 2011a), where one of the inflaton cubic interaction terms allowed by the Galilean symmetry, , contributes to the flat bispectrum with an amplitude (Creminelli et al. 2011a). Here, π is the relevant inflaton scalar degree of freedom, ϵ the usual slowroll parameter and a the scale factor and H the Hubble parameter during inflation. Our constraint at 68% CL (see Table 11) leads to at 68% CL. These interaction terms are similar to those arising in general inflaton field models that include extrinsic curvature terms, e.g., parameterized in the Effective Field Theory approach as (Bartolo et al. 2010a), which contribute to a flat bispectrum with an amplitude . In this case, we obtain at 68% CL.
9.2. Implications for effective field theory of inflation
The effective field theory approach to inflation (Weinberg 2008; Cheung et al. 2008) provides a general way to scan the NG parameter space of inflationary perturbations. For example, one can expand the Lagrangian of the dynamically relevant degrees of freedom into the dominant operators satisfying some underlying symmetries. We will focus on general singlefield models parametrized by the following operators (up to cubic order) (99)where π is the scalar degree of freedom (ζ = − Hπ). The measurements on and can be used to constrain the magnitude of the inflaton interaction terms and which give respectively and (Senatore et al. 2010; see also Chen et al. 2007b; Chen 2010b). These two operators give rise to shapes that peak in the equilateral configuration that are, nevertheless, slightly different, so that the total NG signal will be a linear combination of the two, possibly leading also to an orthogonal shape. There are two relevant NG parameters, c_{s}, the sound speed of the the inflaton fluctuations, and M_{3} which characterizes the amplitude of the other operator . Following Senatore et al. (2010) we will focus on the dimensionless parameter . For example, DBI inflationary models corresponds to , while the noninteracting model (vanishing NG) correspond to c_{s} = 1 and M_{3} = 0 (or ).
Fig. 24
68%, 95%, and 99.7% confidence regions in the parameter space , defined by thresholding χ^{2} as described in the text. 
Fig. 25
68%, 95%, and 99.7% confidence regions in the singlefield inflation parameter space , obtained from Fig. 24 via the change of variables in Eq. (100). 
The mean values of the estimators for equilateral and orthogonal NG amplitudes are given in terms of c_{s} and by (100)where , and the coefficients are computed from the Fisher correlation matrix between the equilateral and orthogonal template bispectra and the theoretical bispectra arising from the two operators and . Given our constraints on and , and the covariance matrix C of the joint estimators, we can define a χ^{2} statistic given by , where the vector v is given by . , where i = {equilateral, orthogonal}, are the joint estimates of the equilateral and orthogonal f_{NL} measured by Planck and is given by Eq. (100). Figure 24 shows the 68%,95%, and 99.7% confidence regions for and , obtained by requiring χ^{2} ≤ 2.28,5.99, and 11.62 respectively, as appropriate for a χ^{2} variable with two degrees of freedom. The corresponding confidence regions in the parameter space are shown in Fig. 25. After marginalizing over we find the following conservative bound on the inflaton soundspeed (101)
Note that we have also looked explicitly for the nonseparable shapes in Sect. 7.3.1, in particular the two effective field theory shapes and the DBI inflation shape (see Eqs. (5)–(7)).
9.3. Multifield models
Curvaton models: Planck NG constraints have interesting implications for the simplest adiabatic curvaton models. They predict (Bartolo et al. 2004d,c) (102)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” evaluated at the epoch of the curvaton decay in the sudden decay approximation. Assuming a prior 0 <r_{D}< 1, given our constraint at 68% CL, we obtain (103)In Planck Collaboration XXII (2014) a limit on r_{D} is derived from the constraints on isocurvature perturbations under the assumption that there is some residual isocurvature fluctuations in the curvaton field. For this restricted case, they find r_{D}> 0.98 (95% CL), compatible with the constraint obtained here.
Quasisingle field inflation: it is beyond the scope of this paper to perform a general multifield analysis employing the local NG constraints. However, we have performed a detailed analysis for the quasisingle field models (see Eq. (12)). Quasisingle field (QSF) inflation models (Chen & Wang 2010b,a; Baumann & Green 2012) are a natural consequence of inflation modelbuilding in string theory and supergravity (see Sect. 2.2). In addition to the inflaton field, these models have extra fields with masses of order the Hubble parameter, which are stabilized by supersymmetry. A distinctive observational signature of these massive fields is a oneparameter family of large NG whose squeezed limits interpolate between the local and the equilateral shape. Therefore, by measuring the precise momentumdependence of the squeezed configurations in the NG, in principle, we are directly measuring the parameters of the theory naturally determined by the fundamental principle of supersymmetry. These models produce a bispectrum (Eq. (12)) depending on two parameters ν, , with a shape that interpolates between the local shape, where ν = 1.5 and the equilateral shape, where ν = 0.
Results are shown in Fig. 26 (see Sect. 7.3.6 for details of the analyses). The best fit value corresponds to ν = 1.5, f_{NL} = 4.79 which would imply, within this scenario, that the extra field different from the inflaton has a mass m ≪ H. However, the figure shows that there is no preferred value for ν with all values allowed at 3σ.
Fig. 26
68%, 95%, and 99.7% confidence intervals for ν and for quasisingle field inflation. The best fit value of ν = 1.5, is marked with an X. The contours were calculated using MC methods by creating 2 × 10^{9} simulations using the β covariance matrix around this best fit model. 
9.4. Nonstandard inflation models
Constraints on excited initial states: results from Sect. 7.3 constrain a variety of models with flattened bispectra, often in combination with a nontrivial squeezed limit. The most notable examples are bispectra produced in excited initial state models (nonBunchDavies vacua), which can be generated by strong disturbances away from background slowroll evolution or additional transPlanckian physics (Chen et al. 2007b; Holman & Tolley 2008; Meerburg et al. 2009; Agullo & Parker 2011). The constraints we have obtained are summarized in Table 11, and cover four representative cases (see Eqs. (14), (15)) in the literature. We find no strong evidence for these flattened bispectra in the Planck data after subtraction of the ISWlensing signal, with which all these models have some correlation. This is consistent with an earlier constraint on the NBD model from WMAP (Fergusson et al. 2012). However, this investigation is limited by the present resolution of the polynomial modal estimator (n_{max} = 600), so more strongly flattened models are not excluded by this analysis.
Directional dependence motivated by fields: directionallydependent bispectra (Eq. (19)), motivated by inflation with gauge fields, have also been constrained (see Table 11). For example, models with a kinetic term of the vector field(s) ℒ = − I^{2}(φ)F^{2}/ 4, where F^{2} is the strength of the gauge field, and I(φ) is a function of the inflaton field which, with an appropriate time dependence (see, e.g., Ratra 1992), can generate vector fields during inflation. For these models the L = 0,2 modes in the bispectrum are excited with , where X_{L = 0} = (80/3) and X_{L = 2} = − (10/6), respectively (Barnaby et al. 2012a; Bartolo et al. 2013; Shiraishi et al. 2013). Here g_{∗} is the amplitude of a quadrupolar anisotropy in the power spectrum (see, e.g., Ackerman et al. 2007) and N is the number of efolds before the end of inflation when the relevant scales exit the horizon. These modes therefore relate the bispectrum amplitude to the level of statistical anisotropy in the power spectrum. Marginalizing over a prior 50 ≤ N ≤ 70 assuming uniform priors on g_{∗}, our constraints in Table 11 lead to the limits −0.05 <g_{∗}< 0.05 and −0.36 <g_{∗}< 0.36 from the L = 0, L = 2 modes respectively (95% CL). Note, however, that in the current analysis only a modest correlation was possible with the shape corresponding to the L = 2 mode. These results apply to all models where curvature perturbations are sourced by a I^{2}(φ)F^{2} term (see references in Shiraishi et al. 2013).
Feature models: nonscaleinvariant oscillatory bispectrum shapes can be generated by sharp or periodic features in the inflaton potential, with particular recent interest on axion monodromy models motivated by string theory (see Sect. 2.3). We have undertaken a survey of simple feature models (Eq. (16)) which have a periodicity accessible to the polynomial modal estimator (wavenumbers K = k_{1} + k_{2} + k_{3} ≳ 0.01), a twoparameter space spanned by K and a phase φ. There are interesting best fit models for the Planck CMB bispectrum around K = 0.01875, φ = 0 (that is, with a largeℓ bispectrum periodicity around Δℓ = 260), with results shown in Table 12. We note important caveats on the statistics of parameter surveys like this in Sect. 7.3.3; given the large numbers of feature models studied, we have to apply a higher threshold for statistical significance as shown for a survey of 200 Gaussian simulations. This feature survey takes forward earlier results for the WMAP data (Fergusson et al. 2012), with the apparent fit reflecting the signal observed in the Planck CMB reconstruction (see Fig. 7). Most attention on feature models has been motivated by the simplest singlefield case for which there are correlated signals predicted in the bispectrum and power spectrum (see e.g., Chen et al. 2007a; Adshead et al. 2012). In this case, regions with small k ~ 0.001 are favoured for producing an observable bispectrum^{17}, given existing WMAP power spectrum constraints on these models. Periodicities Δℓ ≲ 20 are anticipated (see Adshead et al. 2012) which are not accessible to the present modal bispectrum estimator analysis, but which are discussed in Planck Collaboration XXII (2014). Conversely, the Planck feature model survey using the power spectrum (Planck Collaboration XXII 2014) does cover bispectrum scales and parameters investigated in this paper. An analysis of the Bayesian posterior probability (Planck Collaboration XXII 2014) does not appear to provide evidence favoring parameters associated with the current bestfit bispectrum model. More detailed analysis using the specific bispectrum envelope for the singlefield feature solution is being undertaken.
Resonance models: we have also investigated resonance models of Eq. (17) such as axion monodromy and enfolded resonance models of Eq. (18), in which nonBunchDavies vacua are excited by sharp features. This limited analysis focuses on periodicities associated with the bestfit feature model and the results are described in Tables C.1 and C.2. No significant signal was found in this domain for either of these two models. However, we note that the logarithmic dependence of the bispectrum creates challenges at low k, as we must ensure important features do not fall below the modal resolution. This restricts the present survey range, which will be extended in future. Again, we note that most attention on these models has focused on higher ℓperiodicities than those accessible to the present modal estimator (see e.g., Flauger & Pajer 2011; Peiris et al. 2013), for the same reason as the feature models. A resonance model survey using the Planck power spectrum has been undertaken and the results can be found in Planck Collaboration XXII (2014); however, this also currently excludes high frequencies that have received attention in the literature.
Warm inflation: this model, where dissipative effects are important, predicts (Moss & Xiong 2007) where the dissipation parameter r_{d} = Γ/(3H) must be large for strong dissipation. The limit from WMAP is r_{d} ≤ 2.8 × 10^{4} (Moss & Xiong 2007). Assuming a prior 0 ≤ log _{10}r_{d} ≤ 4, our constraint at 68% CL (see Sect. 7.3.5) yields a limit on the dissipation parameter of log _{10}r_{d} ≤ 2.6 (95% CL), improving the previous limit by nearly two orders of magnitude. The stronglydissipative regime with r_{d} ≳ 2.5 still remains viable; however, the Planck constraint puts the model in a regime which can lead to an overproduction of gravitinos (see Hall & Peiris 2008, and references within).
9.5. Alternatives to inflation
Perhaps the most striking example is given by the ekpyrotic/cyclic models (for a review, see Lehners 2010) proposed as alternative to inflationary models. Typically they predict a local NG . In particular, the socalled “ekpyrotic conversion” mechanism (in which isocurvature perturbations are converted into curvature perturbations during the ekpyrotic phase) yields , where c_{1} is a parameter in the potential, requiring 10 ≳ c_{1} ≳ 20 for compatibility with power spectrum constraints. This case was ~4σ discrepant with WMAP data, and with Planck it is decisively ruled out given our bounds at 68% CL (see Table 9) yielding c_{1} ≤ 4.2 at 95% CL. The predictions for the local bispectrum from other ekpyrotic models (based on the so called “kinetic conversion” taking place after the ekpyrotic phase) yield , where values −1 <κ_{3}< 5 and ϵ ~ 100 are natural (Lehners 2010). Taking ϵ ~ 100, in this case we obtain −0.9 <κ_{3}< 0.6 and −0.2 <κ_{3}< 1.3 at 95% CL, for the plus and minus sign in respectively, significantly restricting the viable parameter space of this model.
9.6. Implications of CMB trispectrum results
The nondetection of localtype f_{NL} by Planck raises the immediate question of whether there might still be a large trispectrum. In this first analysis we have focused on the local shape τ_{NL}. Most inflation models predict (Byrnes et al. 2006; Elliston et al. 2012), and hence given our tight f_{NL} limits, would be predicted to be very small. This is easily consistent with our conservative upper limit τ_{NL}< 2800, and also with the significantly smaller signals found in the modulation dipole and quadrupole. Our upper limit also restricts the freedom of curvatonlike models with quadratic terms that are nearly uncorrelated with the curvature perturbation, which could in principle have f_{NL} ~ 0 but a large trispectrum (Byrnes & Choi 2010).
10. Conclusions
In this paper we have derived constraints on primordial NG, using the CMB maps derived from the Planck nominal mission data. Using three optimal bispectrum estimators – KSW, binned, and modal – we obtained consistent values for the primordial local, equilateral, and orthogonal bispectrum amplitudes, quoting as our final result , , and (68 % CL statistical). Hence there is no evidence for primordial NG of one of these shapes. We did, however, measure the ISWlensing bispectrum expected in the ΛCDM scenario, as well as a contribution from diffuse point sources, and these contributions are clearly seen in the form of the associated skewC_{ℓ}. Indeed, the detection of ISWlensing and Poisson point source skewC_{ℓ} with the right shapes and at expected levels is good evidence that the data processing is not destroying NG in the maps, and gives confidence that the nondetections of primordial NG are robust. These results have been confirmed by measurements using the wavelet bispectrum, and Minkowski functional estimators, and demonstrated to be stable for the four different component separation techniques SMICA, NILC, SEVEM, and CR, showing their robustness against foreground residuals. They have also passed an extensive suite of tests studying the dependence on the maximum multipole number and the mask, consistency checks between frequency channels, and several null tests. In addition, we have summarized in this paper an extensive validation campaign for the three optimal bispectrum estimators on Gaussian and nonGaussian simulations.
Extending our analysis beyond estimates of individual shape amplitudes, we presented modelindependent, threedimensional reconstructions of the Planck CMB bispectrum using the modal and binned bispectrum estimators. These results were also shown to be fully consistent between the different component separation techniques even for the full bispectrum, and contained no significant NG signals of a type not captured by our standard templates at low multipole values. At high multipoles, some indications of unidentified NG signals were found, as also evidenced by the results from the skewC_{ℓ} estimator. Further study will be required to ascertain whether these indications are due to foreground residuals, beams, data processing, or a more interesting signal.
Using the modal decomposition, we have presented constraints on key primordial NG scenarios, including general singlefield models with nonseparable shapes, excited initial states (nonBunchDavies vacua), and directionallydependent vector models. We have also undertaken an initial survey of scaledependent feature and resonance models.
Moving beyond threepoint correlations, we have obtained limits from the Planck data on the amplitude of the local fourpoint function. Our trispectrum reconstruction yielded a signal consistent with zero except for an anomalously large octopole. Frequency dependence indicated that this was unlikely to be primordial, but allowing the signal to be primordial we placed an upper limit τ_{NL}< 2800 (95% CL).
We have discussed the impact of these results on the inflationary model space, and derived bispectrum constraints on a selection of specific inflationary mechanisms, including both general singlefield inflationary models (extensions to the standard singlefield slowroll case) and multifield models. Our results led to a lower bound on the speed of sound, c_{s} ≥ 0.02 (95% CL), in an effective field theory approach to inflationary models which includes models with nonstandard kinetic terms. Strong constraints on other scenarios such as IR DBI, kinflation, inflationary models involving gauge fields, and warm inflation have been obtained. Within the class of multifield models, our measurements limit the curvaton decay fraction to r_{D} ≥ 0.15 (95% CL). Ekpyrotic/cyclic scenarios were shown to be under pressure from the Planck data: we robustly ruled out the socalled “ekpyrotic conversion” mechanism, and found that the parameter space of the “kinetic conversion” mechanism is significantly limited.
With these results, the paradigm of standard singlefield slowroll inflation has survived its most stringent tests todate.
Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states (in particular the lead countries France and Italy), with contributions from NASA (USA) and telescope reflectors provided by a collaboration between ESA and a scientific consortium led and funded by Denmark.
Specifically, one can define the shape of the bispectrum as the dependence of F(k_{1},k_{2},k_{3})(k_{1}k_{2}k_{3})^{2} on the ratios of momenta, e.g., (k_{2}/k_{1}) and (k_{3}/k_{1}), once the overall scale of the triangle K = k_{1} + k_{2} + k_{3} is fixed. The scale dependence of the bispectrum can be characterized by the dependence of F(k_{1},k_{2},k_{3})(k_{1}k_{2}k_{3})^{2} on the overall scale K, once the ratios (k_{2}/k_{1}) and (k_{3}/k_{1}) are fixed (see, e.g., Chen 2010b).
For the curvature perturbation, we follow the notation and sign conventions of Komatsu et al. (2011). ζ is also sometimes denoted ℛ (see e.g., Lidsey et al. 1997; Lyth & Riotto 1999, and references therein), while the comoving curvature perturbation ℛ as defined, e.g., in Malik & Wands (2009) is such that ℛ = − ζ.
This has been shown in the pioneering research which demonstrated that perturbations produced in singlefield models of slowroll inflation are characterized by a lowamplitude NG (Salopek & Bond 1990; Falk et al. 1993; Gangui et al. 1994). Later Acquaviva et al. (2003) and Maldacena (2003) obtained a complete quantitative prediction for the nonlinearity parameter in singlefield slowroll inflation models, also showing that the predicted NG is characterized by a shape dependence which is more complex than suggested by previous results expressed in terms of the simple parameterization (Gangui et al. 1994; Verde et al. 2000; Wang & Kamionkowski 2000; Komatsu & Spergel 2001), where Φ_{L} is the linear gravitational potential.
An effectively singlefield model with a nonstandard kinetic term and a reduced sound speed for the adiabatic perturbation modes might also arise in coupled multifield systems, where the heavy fields are integrated out: see discussions in, e.g., Tolley & Wyman (2010), Achúcarro et al. (2011), Shiu & Xu (2011).
Notice that the two shapes (5) and (6) correspond to a linear combination of the two shapes found in Chen et al. (2007b).
This may happen, for instance, if the inflaton field is coupled to the other scalar degrees of freedom, as expected on particle physics grounds. These scalar degrees of freedom may have large selfinteractions, so that their quantum fluctuations are intrinsically nonGaussian, because, unlike the inflaton case, the selfinteraction strength in such an extra scalar sector does not suffer from the usual slowroll conditions.
See Kim et al. (2013) for other debiasing techniques.
While most of the modal results in this paper come from the most accurate 600 modes pipeline, a few computationally intensive data validation tests of Sect. 8 also use the fast 300 modes version; therefore the results in this section also provide a direct validation of the fast modal pipeline.
While the binning with 48 bins and ℓ_{max} = 2000 used in the validation tests of Sects. 6.1.2 and 6.1.3 is also slightly different from the binning used for the data analysis with 51 bins and ℓ_{max} = 2500, these differences are very small and the binnings have very similar correlation coefficients of 0.99 or more for local and equilateral shapes, and 0.95 for orthogonal.
The lower ℓ_{max} for the modal pipeline is also a conservative choice in view of the large survey of “nonstandard” models that will be presented in Sect. 7.3
The boundary values of the bins are: 2, 4, 10, 18, 27, 39, 55, 75, 99, 130, 170, 224, 264, 321, 335, 390, 420, 450, 518, 560, 615, 644, 670, 700, 742, 800, 850, 909, 950, 979, 1005, 1050, 1110, 1150, 1200, 1230, 1260, 1303, 1346, 1400, 1460, 1510, 1550, 1610, 1665, 1725, 1795, 1871, 1955, 2091, 2240, and 2500 (i.e., the first bin is [2, 3], the second [4, 9], etc., while the last one is [2240, 2500]).
Raw MFs V_{k} depend on the Gaussian part of fields through a normalization factor A_{k} that is a function only of the power spectrum shape. We therefore normalize functionals v_{k} = V_{k}/A_{k} to focus on NG; see Planck Collaboration XXIII (2014) and references therein.
While not rigorous, this argument captures the leading effect since Galactic foregrounds predominantly contaminate large scales. In principle, positively skewed, small scale foreground residuals (ℓ > 60), or the negatively skewed SZ effect, can contribute positive bias. Our simulations with foreground residuals demonstrate that these contributions are subdominant.
Planck Collaboration XVI (2014) confirms an anomaly in the power spectrum at 20 ≲ ℓ ≲ 40 first noted in WMAP, which leads to an improvement in likelihood when fitted with a feature in the inflationary potential (Peiris et al. 2003). Unfortunately, the bestfit parameters in this case do not lead to an observable bispectrum (Adshead et al. 2011).
Acknowledgments
The development of Planck has been supported by: ESA; CNES and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, JA and RES (Spain); Tekes, AoF and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. The primary platform for analysis with the modal and KSW estimators was the COSMOS supercomputer, part of the DiRAC HPC Facility supported by STFC and the UK Large Facilities Capital Fund. We gratefully acknowledge IN2P3 Computer Center (http://cc.in2p3.fr) for providing a significant amount of the computing resources and services needed for the analysis with the binned bispectrum estimator. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US Department of Energy under Contract No. DEAC0205CH11231. We acknowledge the use of resources from the Norwegian national super computing facilities NOTUR. We also acknowledge the IAP magique3 computer facilities as well as the curvaton computer at the Department of Physics of the University of Illinois (Urbana, IL). We further acknowledge the computer resources and technical assistance provided by the Spanish Supercomputing Network nodes at Universidad de Cantabria and Universidad Politécnica de Madrid as well as the support provided by the Advanced Computing and eScience team at IFCA.
References
 Achúcarro, A., Gong, J.O., Hardeman, S., Palma, G. A., & Patil, S. P. 2011, J. Cosmol. Astropart. Phys., 1, 30 [Google Scholar]
 Ackerman, L., Carroll, S. M., & Wise, M. B. 2007, Phys. Rev. D, 75, 083502 [NASA ADS] [CrossRef] [Google Scholar]
 Acquaviva, V., Bartolo, N., Matarrese, S., & Riotto, A. 2003, Nucl. Phys. B, 667, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Adshead, P., Hu, W., Dvorkin, C., & Peiris, H. V. 2011, Phys. Rev. D, 84, 3519 [NASA ADS] [CrossRef] [Google Scholar]
 Adshead, P., Dvorkin, C., Hu, W., & Lim, E. A. 2012, Phys. Rev. D, 85, 3531 [NASA ADS] [Google Scholar]
 Agullo, I., & Parker, L. 2011, Gen. Relativ. Gravit., 43, 2541 [Google Scholar]
 Alishahiha, M., Silverstein, E., & Tong, D. 2004, Phys. Rev. D, 70, 3505 [Google Scholar]
 Allen, T. J., Grinstein, B., & Wise, M. B. 1987, Phys. Lett. B, 197, 66 [NASA ADS] [CrossRef] [Google Scholar]
 Amendola, L., Catena, R., Masina, I., et al. 2011, J. Cosmol. Astropart. Phys., 7, 27 [Google Scholar]
 Antoine, J.P., & Vandergheynst, P. 1998, J. Math. Phys., 39, 3987 [NASA ADS] [CrossRef] [Google Scholar]
 ArkaniHamed, N., Creminelli, P., Mukohyama, S., & Zaldarriaga, M. 2004, J. Cosmol. Astropart. Phys., 4, 1 [Google Scholar]
 ArmendarizPicon, C., Damour, T., & Mukhanov, V. 1999, Phys. Lett. B, 458, 209 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Arroja, F., Mizuno, S., & Koyama, K. 2008, J. Cosmol. Astropart. Phys., 8, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Arroja, F., Mizuno, S., Koyama, K., & Tanaka, T. 2009, Phys. Rev. D, 80, 3527 [NASA ADS] [Google Scholar]
 Ashoorioon, A., & Shiu, G. 2011, J. Cosmol. Astropart. Phys., 3, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Assassi, V., Baumann, D., & Green, D. 2012, J. Cosmol. Astropart. Phys., 11, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Babich, D. 2005, Phys. Rev. D, 72, 3003 [NASA ADS] [CrossRef] [Google Scholar]
 Babich, D., Creminelli, P., & Zaldarriaga, M. 2004, J. Cosmol. Astropart. Phys., 8, 9 [Google Scholar]
 Bardeen, J. M. 1980, Phys. Rev. D, 22, 1882 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Barnaby, N., & Cline, J. M. 2006, Phys. Rev. D, 73, 106012 [NASA ADS] [CrossRef] [Google Scholar]
 Barnaby, N., & Cline, J. M. 2008, J. Cosmol. Astropart. Phys., 6, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Barnaby, N., & Peloso, M. 2011, Phys. Rev. Lett., 106, 1301 [NASA ADS] [CrossRef] [Google Scholar]
 Barnaby, N., Namba, R., & Peloso, M. 2011, J. Cosmol. Astropart. Phys., 4, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Barnaby, N., Namba, R., & Peloso, M. 2012a, Phys. Rev. D, 85, 3523 [NASA ADS] [Google Scholar]
 Barnaby, N., Pajer, E., & Peloso, M. 2012b, Phys. Rev. D, 85, 3525 [NASA ADS] [Google Scholar]
 Bartolo, N., & Riotto, A. 2009, J. Cosmol. Astropart. Phys., 3, 17 [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2002, Phys. Rev. D, 65, 103505 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Komatsu, E., Matarrese, S., & Riotto, A. 2004a, Phys. Rep., 402, 103 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2004b, J. Cosmol. Astropart. Phys., 1, 3 [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2004c, Phys. Rev. Lett., 93, 231301 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2004d, Phys. Rev. D, 69, 043503 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2005, J. Cosmol. Astropart. Phys., 8, 10 [NASA ADS] [Google Scholar]
 Bartolo, N., Fasiello, M., Matarrese, S., & Riotto, A. 2010a, J. Cosmol. Astropart. Phys., 8, 8 [Google Scholar]
 Bartolo, N., Fasiello, M., Matarrese, S., & Riotto, A. 2010b, J. Cosmol. Astropart. Phys., 9, 35 [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2010c, Adv. Astron., 2010 [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2012, J. Cosmol. Astropart. Phys., 2, 17 [Google Scholar]
 Bartolo, N., Matarrese, S., Peloso, M., & Ricciardone, A. 2013, Phys. Rev. D, 87, 023504 [NASA ADS] [CrossRef] [Google Scholar]
 Baumann, D., & Green, D. 2012, Phys. Rev. D, 85, 3520 [Google Scholar]
 Baumann, D., & McAllister, L. 2007, Phys. Rev. D, 75, 3508 [Google Scholar]
 Bean, R., Shandera, S. E., Tye, S.H. H., & Xu, J. 2007, J. Cosmol. Astropart. Phys., 5, 4 [Google Scholar]
 Bean, R., Chen, X., Peiris, H., & Xu, J. 2008, Phys. Rev. D, 77, 023527 [NASA ADS] [CrossRef] [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [Google Scholar]
 Berera, A. 1995, Phys. Rev. Lett., 75, 3218 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Bernardeau, F., & Uzan, J.P. 2002, Phys. Rev. D, 66, 103506 [NASA ADS] [CrossRef] [Google Scholar]
 Bernardeau, F., Kofman, L., & Uzan, J.P. 2004, Phys. Rev. D, 70, 083004 [NASA ADS] [CrossRef] [Google Scholar]
 Bond, J. R., Frolov, A. V., Huang, Z., & Kofman, L. 2009, Phys. Rev. Lett., 103, 071301 [NASA ADS] [CrossRef] [Google Scholar]
 Boubekeur, L., Creminelli, P., D’Amico, G., Noreña, J., & Vernizzi, F. 2009, J. Cosmol. Astropart. Phys., 8, 29 [Google Scholar]
 Bucher, M., & Zhu, Y. 1997, Phys. Rev. D, 55, 7415 [NASA ADS] [CrossRef] [Google Scholar]
 Bucher, M., van Tent, B., & Carvalho, C. S. 2010, MNRAS, 407, 2193 [NASA ADS] [CrossRef] [Google Scholar]
 Byrnes, C. T., & Choi, K.Y. 2010, Adv. Astron., 2010 [arXiv:1002.3110] [Google Scholar]
 Byrnes, C. T., Sasaki, M., & Wands, D. 2006, Phys. Rev. D, 74, 123519 [NASA ADS] [CrossRef] [Google Scholar]
 Cardoso, J., Martin, M., Delabrouille, J., Betoule, M., & Patanchon, G. 2008, IEEE Select. Topics Signal Process., 2, 735 [NASA ADS] [CrossRef] [Google Scholar]
 Cayon, L., MartínezGonzález, E., Argueso, F., Banday, A. J., & Gorski, K. M. 2003, MNRAS, 339, 1189 [NASA ADS] [CrossRef] [Google Scholar]
 Challinor, A., & van Leeuwen, F. 2002, Phys. Rev. D, 65, 3001 [Google Scholar]
 Chambers, A., & Rajantie, A. 2008, Phys. Rev. Lett., 100, 1302 [Google Scholar]
 Chen, X. 2005a, J. High Energy Phys., 8, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X. 2005b, Phys. Rev. D, 71, 3506 [Google Scholar]
 Chen, X. 2005c, Phys. Rev. D, 72, 3518 [NASA ADS] [Google Scholar]
 Chen, X. 2010a, J. Cosmol. Astropart. Phys., 12, 3 [Google Scholar]
 Chen, X. 2010b, Adv. Astron. [arXiv:1002.1416] [Google Scholar]
 Chen, X. 2011, Phys. Lett. B, 706, 111 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., & Wang, Y. 2010a, Phys. Rev. D, 81, 063511 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., & Wang, Y. 2010b, J. Cosmol. Astropart. Phys., 4, 27 [Google Scholar]
 Chen, X., Easther, R., & Lim, E. A. 2007a, J. Cosmol. Astropart. Phys., 6, 23 [Google Scholar]
 Chen, X., Huang, M.X., Kachru, S., & Shiu, G. 2007b, J. Cosmol. Astropart. Phys., 1, 2 [Google Scholar]
 Chen, X., Easther, R., & Lim, E. A. 2008, J. Cosmol. Astropart. Phys., 4, 10 [Google Scholar]
 Chen, X., Hu, B., Huang, M.X., Shiu, G., & Wang, Y. 2009, J. Cosmol. Astropart. Phys., 8, 8 [NASA ADS] [Google Scholar]
 Chen, X., Firouzjahi, H., Namjoo, M. H., & Sasaki, M. 2013, Europhys. Lett., 102, 59001 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Cheung, C., Fitzpatrick, A. L., Kaplan, J., Senatore, L., & Creminelli, P. 2008, J. High Ener. Phys., 3, 14 [Google Scholar]
 Cicoli, M., Tasinato, G., Zavala, I., Burgess, C. P., & Quevedo, F. 2012, J. Cosmol. Astropart. Phys., 5, 39 [Google Scholar]
 Cooray, A., & Hu, W. 2000, ApJ, 534, 533 [NASA ADS] [CrossRef] [Google Scholar]
 Creminelli, P., & Zaldarriaga, M. 2004a, J. Cosmol. Astropart. Phys., 10, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Creminelli, P., & Zaldarriaga, M. 2004b, Phys. Rev. D, 70, 083532 [NASA ADS] [CrossRef] [Google Scholar]
 Creminelli, P., Nicolis, A., Senatore, L., Tegmark, M., & Zaldarriaga, M. 2006, J. Cosmol. Astropart. Phys., 5, 4 [Google Scholar]
 Creminelli, P., D’Amico, G., Musso, M., Noreña, J., & Trincherini, E. 2011a, J. Cosmol. Astropart. Phys., 2, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Creminelli, P., Pitrou, C., & Vernizzi, F. 2011b, J. Cosmol. Astropart. Phys., 11, 25 [Google Scholar]
 Cruz, M., Cayon, L., MartínezGonzález, E., Vielva, P., & Jin, J. 2007, ApJ, 655, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., MaciasPerez, J. F., MartínezGonzález, E., et al. 2008, A&A, 486, 383 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Curto, A., MartínezGonzález, E., & Barreiro, R. B. 2009a, ApJ, 706, 399 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., MartínezGonzález, E., Mukherjee, P., et al. 2009b, MNRAS, 393, 615 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., MartínezGonzález, E., & Barreiro, R. B. 2011a, MNRAS, 412, 1038 [NASA ADS] [Google Scholar]
 Curto, A., MartínezGonzález, E., Barreiro, R. B., & Hobson, M. P. 2011b, MNRAS, 417, 488 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., MartínezGonzález, E., & Barreiro, R. B. 2012, MNRAS, 426, 1361 [NASA ADS] [CrossRef] [Google Scholar]
 Curto, A., Tucci, M., GonzálezNuevo, J., et al. 2013, MNRAS, 432, 728 [NASA ADS] [CrossRef] [Google Scholar]
 De Troia, G., Ade, P. A. R., Bock, J. J., et al. 2007, ApJ, 670, L73 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Cardoso, J.F., & Patanchon, G. 2003, MNRAS, 346, 1089 [NASA ADS] [CrossRef] [Google Scholar]
 Delabrouille, J., Cardoso, J., Le Jeune, M., et al. 2009, A&A, 493, 835 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Delabrouille, J., Betoule, M., Melin, J.B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Donzelli, S., Hansen, F. K., Liguori, M., & Maino, D. 2009, ApJ, 706, 1226 [NASA ADS] [CrossRef] [Google Scholar]
 Donzelli, S., Hansen, F. K., Liguori, M., Marinucci, D., & Matarrese, S. 2012, ApJ, 755, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Ducout, A., Bouchet, F. R., Colombi, S., Pogosyan, D., & Prunet, S. 2013, MNRAS, 429, 2104 [NASA ADS] [CrossRef] [Google Scholar]
 Dvali, G., Gruzinov, A., & Zaldarriaga, M. 2004a, Phys. Rev. D, 69, 083505 [NASA ADS] [CrossRef] [Google Scholar]
 Dvali, G., Gruzinov, A., & Zaldarriaga, M. 2004b, Phys. Rev. D, 69, 023505 [NASA ADS] [CrossRef] [Google Scholar]
 Elliston, J., Alabidi, L., Huston, I., Mulryne, D., & Tavakol, R. 2012, J. Cosmol. Astropart. Phys., 9, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Elsner, F., & Wandelt, B. D. 2009, ApJS, 184, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Elsner, F., & Wandelt, B. D. 2013, A&A, 549, A111 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Endlich, S., Nicolis, A., & Wang, J. 2013, J. Cosmol. Astropart. Phys., 10, 11 [Google Scholar]
 Enqvist, K., & Sloth, M. S. 2002, Nucl. Phys. B, 626, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Enqvist, K., Jokinen, A., Mazumdar, A., Multamäki, T., & Väihkönen, A. 2005, Phys. Rev. Lett., 94, 161301 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Novikov, D. I., Lilje, P. B., Banday, A. J., & Górski, K. M. 2004, ApJ, 612, 64 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Dickinson, C., Lawrence, C. R., et al. 2006, ApJ, 641, 665 [NASA ADS] [CrossRef] [Google Scholar]
 Eriksen, H. K., Jewell, J. B., Dickinson, C., et al. 2008, ApJ, 676, 10 [NASA ADS] [CrossRef] [Google Scholar]
 Falk, T., Rangarajan, R., & Srednicki, M. 1993, ApJ, 403, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Faÿ, G., Guilloux, F., Betoule, M., et al. 2008, Phys. Rev. D, 78, 083013 [NASA ADS] [CrossRef] [Google Scholar]
 Feeney, S. M., Johnson, M. C., Mortlock, D. J., & Peiris, H. V. 2011, Phys. Rev. D, 84, 043507 [NASA ADS] [CrossRef] [Google Scholar]
 Fergusson, J. R., & Shellard, E. P. S. 2007, Phys. Rev. D, 76, 083523 [NASA ADS] [CrossRef] [Google Scholar]
 Fergusson, J. R., Liguori, M., & Shellard, E. P. S. 2010a, Phys. Rev. D, 82, 023502 [NASA ADS] [CrossRef] [Google Scholar]
 Fergusson, J. R., Regan, D. M., & Shellard, E. P. S. 2010b [arXiv:1012.6039] [Google Scholar]
 Fergusson, J. R., Liguori, M., & Shellard, E. P. S. 2012, J. Cosmol. Astropart. Phys., 12, 32 [Google Scholar]
 FernándezCobos, R., Vielva, P., Barreiro, R. B., & MartínezGonzález, E. 2012, MNRAS, 420, 2162 [NASA ADS] [CrossRef] [Google Scholar]
 Ferreira, P. G., & Magueijo, J. 1997, Phys. Rev. D, 56, 4578 [NASA ADS] [CrossRef] [Google Scholar]
 Flauger, R., & Pajer, E. 2011, J. Cosmol. Astropart. Phys., 1, 17 [Google Scholar]
 Flauger, R., McAllister, L., Pajer, E., Westphal, A., & Xu, G. 2010, J. Cosmol. Astropart. Phys., 6, 9 [Google Scholar]
 Gangui, A., & Martin, J. 2000a, Phys. Rev. D, 62, 3004 [CrossRef] [Google Scholar]
 Gangui, A., & Martin, J. 2000b, MNRAS, 313, 323 [NASA ADS] [CrossRef] [Google Scholar]
 Gangui, A., Lucchin, F., Matarrese, S., & Mollerach, S. 1994, ApJ, 430, 447 [NASA ADS] [CrossRef] [Google Scholar]
 Geller, D., & Mayeli, A. 2009a, Math. Z., 262, 895 [CrossRef] [Google Scholar]
 Geller, D., & Mayeli, A. 2009b, Math. Z., 263, 235 [CrossRef] [Google Scholar]
 Giovi, F., Baccigalupi, C., & Perrotta, F. 2005, Phys. Rev. D, 71, 3009 [NASA ADS] [CrossRef] [Google Scholar]
 Goldberg, D. M., & Spergel, D. N. 1999, Phys. Rev. D, 59, 3002 [Google Scholar]
 GonzálezNuevo, J., Toffolatti, L., & Argüeso, F. 2005, ApJ, 621, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Green, D., Horn, B., Senatore, L., & Silverstein, E. 2009, Phys. Rev. D, 80, 3533 [CrossRef] [Google Scholar]
 Hall, L. M. H., & Peiris, H. V. 2008, J. Cosmol. Astropart. Phys., 1, 27 [Google Scholar]
 Hamimeche, S., & Lewis, A. 2008, Phys. Rev. D, 77, 3013 [Google Scholar]
 Hanson, D., & Lewis, A. 2009, Phys. Rev. D, 80, 3004 [Google Scholar]
 Hanson, D., Rocha, G., & Górski, K. 2009a, MNRAS, 400, 2169 [NASA ADS] [CrossRef] [Google Scholar]
 Hanson, D., Smith, K. M., Challinor, A., & Liguori, M. 2009b, Phys. Rev. D, 80, 083004 [NASA ADS] [CrossRef] [Google Scholar]
 Heavens, A. F. 1998, MNRAS, 299, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Hikage, C., & Matsubara, T. 2012, MNRAS, 425, 2187 [NASA ADS] [CrossRef] [Google Scholar]
 Hikage, C., Komatsu, E., & Matsubara, T. 2006, ApJ, 653, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Hikage, C., Matsubara, T., Coles, P., et al. 2008, MNRAS, 389, 1439 [NASA ADS] [CrossRef] [Google Scholar]
 Hodges, H. M., Blumenthal, G. R., Kofman, L. A., & Primack, J. R. 1990, Nucl. Phys. B, 335, 197 [NASA ADS] [CrossRef] [Google Scholar]
 Holman, R., & Tolley, A. J. 2008, J. Cosmol. Astropart. Phys., 5, 1 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, M. X., & Shiu, G. 2006, Phys. Rev. D, 74, 1301 [NASA ADS] [Google Scholar]
 Huang, Z., & Vernizzi, F. 2013, Phys. Rev. Lett., 110, 1303 [Google Scholar]
 Izumi, K., Mizuno, S., & Koyama, K. 2012, Phys. Rev. D, 85, 3521 [NASA ADS] [CrossRef] [Google Scholar]
 Junk, V., & Komatsu, E. 2012, Phys. Rev. D, 85, 123524 [NASA ADS] [CrossRef] [Google Scholar]
 Kehagias, A., & Riotto, A. 2012, Nucl. Phys. B, 864, 492 [NASA ADS] [CrossRef] [Google Scholar]
 Keihänen, E., Keskitalo, R., KurkiSuonio, H., Poutanen, T., & Sirviö, A. 2010, A&A, 510, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Khatri, R., & Wandelt, B. D. 2009, Phys. Rev. D, 79, 3501 [Google Scholar]
 Khatri, R., & Wandelt, B. D. 2010, Phys. Rev. D, 81, 3518 [Google Scholar]
 Kim, J., Rotti, A., & Komatsu, E. 2013, J. Cosmol. Astropart. Phys., 4, 21 [NASA ADS] [CrossRef] [Google Scholar]
 Kofman, L. 2003 [arXiv:astroph/0303614] [Google Scholar]
 Kofman, L., Blumenthal, G. R., Hodges, H., & Primack, J. R. 1991, in Largescale Structures and Peculiar Motions in the Universe, eds. D. W. Latham, & L. A. N. da Costa, ASP Conf. Ser., 15, 339 [Google Scholar]
 Kogo, N., & Komatsu, E. 2006, Phys. Rev. D, 73, 3007 [CrossRef] [Google Scholar]
 Kolb, E. W., Riotto, A., & Vallinotto, A. 2006, Phys. Rev. D, 73, 3522 [Google Scholar]
 Komatsu, E. 2010, Class. Quantum Grav., 27, 124010 [Google Scholar]
 Komatsu, E., & Spergel, D. N. 2001, Phys. Rev. D, 63, 3002 [Google Scholar]
 Komatsu, E., Kogut, A., Nolta, M. R., et al. 2003, ApJS, 148, 119 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., Spergel, D. N., & Wandelt, B. D. 2005, ApJ, 634, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., Dunkley, J., Nolta, M. R., et al. 2009, ApJS, 180, 330 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., Smith, K. M., Dunkley, J., et al. 2011, ApJS, 192, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Kosowsky, A., & Kahniashvili, T. 2011, Phys. Rev. Lett., 106, 191301 [NASA ADS] [CrossRef] [Google Scholar]
 Lacasa, F., & Aghanim, N. 2014, A&A, 569, A51 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lacasa, F., Aghanim, N., Kunz, M., & Frommert, M. 2012, MNRAS, 421, 1982 [NASA ADS] [CrossRef] [Google Scholar]
 Lan, X., & Marinucci, D. 2008, Electr. J. Stat., 2, 332 [NASA ADS] [CrossRef] [Google Scholar]
 Langlois, D., RenauxPetel, S., Steer, D. A., & Tanaka, T. 2008a, Phys. Rev. Lett., 101, 061301 [NASA ADS] [CrossRef] [Google Scholar]
 Langlois, D., RenauxPetel, S., Steer, D. A., & Tanaka, T. 2008b, Phys. Rev. D, 78, 063523 [NASA ADS] [CrossRef] [Google Scholar]
 Leach, S. M., Cardoso, J., Baccigalupi, C., et al. 2008, A&A, 491, 597 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Lehners, J.L. 2010, Adv. Astron., 2010, 67 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A. 2011, J. Cosmol. Astropart. Phys., 10, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A. 2012, J. Cosmol. Astropart. Phys., 6, 23 [Google Scholar]
 Lewis, A., & Challinor, A. 2006, Phys. Rep., 429, 1 [Google Scholar]
 Lewis, A., Challinor, A., & Hanson, D. 2011, J. Cosmol. Astropart. Phys., 3, 18 [Google Scholar]
 Lidsey, J. E., & Huston, I. 2007, J. Cosmol. Astropart. Phys., 7, 2 [Google Scholar]
 Lidsey, J. E., Liddle, A. R., Kolb, E. W., et al. 1997, Rev. Mod. Phys., 69, 373 [NASA ADS] [CrossRef] [Google Scholar]
 Liguori, M., Sefusatti, E., Fergusson, J. R., & Shellard, E. P. S. 2010, Adv. Astron., 2010, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Linde, A., & Mukhanov, V. 1997, Phys. Rev. D, 56, 535 [NASA ADS] [CrossRef] [Google Scholar]
 Lyth, D. H. 2005, J. Cosmol. Astropart. Phys., 11, 6 [Google Scholar]
 Lyth, D. H., & Riotto, A. A. 1999, Phys. Rep., 314, 1 [Google Scholar]
 Lyth, D. H., & Riotto, A. 2006, Phys. Rev. Lett., 97, 121301 [NASA ADS] [CrossRef] [Google Scholar]
 Lyth, D. H., & Rodriguez, Y. 2005, Phys. Rev. Lett., 95, 121302 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Lyth, D. H., & Wands, D. 2002, Phys. Lett. B, 524, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Lyth, D. H., Ungarelli, C., & Wands, D. 2003, Phys. Rev. D, 67, 023503 [NASA ADS] [CrossRef] [Google Scholar]
 Maldacena, J. 2003, J. High Ener. Phys., 5, 13 [Google Scholar]
 Malik, K. A., & Lyth, D. H. 2006, J. Cosmol. Astropart. Phys., 9, 8 [Google Scholar]
 Malik, K. A., & Wands, D. 2009, Phys. Rep., 475, 1 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Mangilli, A., & Verde, L. 2009, Phys. Rev. D, 80, 3007 [Google Scholar]
 Mangilli, A., Wandelt, B., Elsner, F., & Liguori, M. 2013, A&A, 555, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 MartínezGonzález, E., Gallegos, J. E., Argüeso, F., Cayon, L., & Sanz, J. L. 2002, MNRAS, 336, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Matsubara, T. 2010, Phys. Rev. D, 81, 3505 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Vielva, P., Hobson, M. P., MartínezGonzález, E., & Lasenby, A. N. 2007a, MNRAS, 376, 1211 [NASA ADS] [CrossRef] [Google Scholar]
 McEwen, J. D., Vielva, P., Wiaux, Y., et al. 2007b, J. Fourier Anal. Appl., 13, 495 [Google Scholar]
 Meerburg, P. D., van der Schaar, J. P., & Stefano Corasaniti, P. 2009, J. Cosmol. Astropart. Phys., 5, 18 [Google Scholar]
 Mitra, S., Rocha, G., Górski, K. M., et al. 2011, ApJS, 193, 5 [NASA ADS] [CrossRef] [Google Scholar]
 Modest, H. I., Räth, C., Banday, A. J., et al. 2013, MNRAS, 428, 551 [NASA ADS] [CrossRef] [Google Scholar]
 Mollerach, S. 1990, Phys. Rev. D, 42, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Moroi, T., & Takahashi, T. 2001, Phys. Lett. B, 522, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Moss, I. G., & Xiong, C. 2007, J. Cosmol. Astropart. Phys., 4, 7 [Google Scholar]
 Munshi, D., & Heavens, A. 2010, MNRAS, 401, 2406 [NASA ADS] [CrossRef] [Google Scholar]
 Munshi, D., Heavens, A., Cooray, A., et al. 2011, MNRAS, 412, 1993 [NASA ADS] [CrossRef] [Google Scholar]
 Natoli, P., de Troia, G., Hikage, C., et al. 2010, MNRAS, 408, 1658 [NASA ADS] [CrossRef] [Google Scholar]
 Nitta, D., Komatsu, E., Bartolo, N., Matarrese, S., & Riotto, A. 2009, J. Cosmol. Astropart. Phys., 5, 14 [Google Scholar]
 Okamoto, T., & Hu, W. 2002, Phys. Rev. D, 66, 3008 [Google Scholar]
 Pearson, R., Lewis, A., & Regan, D. 2012, J. Cosmol. Astropart. Phys., 3, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Peebles, P. J. E. 1997, ApJ, 483, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H. V., Komatsu, E., Verde, L., et al. 2003, ApJS, 148, 213 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H., Baumann, D., Friedman, B., & Cooray, A. 2007, Phys. Rev. D, 76, 103517 [NASA ADS] [CrossRef] [Google Scholar]
 Peiris, H. V., Easther, R., & Flauger, R. 2013, J. Cosmol. Astropart. Phys., 9, 18 [Google Scholar]
 Pettinari, G. W., Fidler, C., Crittenden, R., Koyama, K., & Wands, D. 2013, J. Cosmol. Astropart. Phys., 4, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Pietrobon, D., Balbi, A., & Marinucci, D. 2006, Phys. Rev. D, 74, 043524 [NASA ADS] [CrossRef] [Google Scholar]
 Pietrobon, D., Cabella, P., Balbi, A., de Gasperis, G., & Vittorio, N. 2009, MNRAS, 396, 1682 [NASA ADS] [CrossRef] [Google Scholar]
 Pietrobon, D., Cabella, P., Balbi, A., et al. 2010, MNRAS, 402, L34 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration 2013, The Explanatory Supplement to the Planck 2013 results, http://www.sciops.esa.int/wikiSI/planckpla/index.php?title=Main_Page (ESA) [Google Scholar]
 Planck Collaboration I. 2014, A&A, 571, A1 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2014, A&A, 571, A2 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2014, A&A, 571, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2014, A&A, 571, A4 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2014, A&A, 571, A5 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2014, A&A, 571, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2014, A&A, 571, A7 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2014, A&A, 571, A8 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2014, A&A, 571, A9 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2014, A&A, 571, A10 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2014, A&A, 571, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2014, A&A, 571, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2014, A&A, 571, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIV. 2014, A&A, 571, A14 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XV. 2014, A&A, 571, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVI. 2014, A&A, 571, A16 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2014, A&A, 571, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2014, A&A, 571, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2014, A&A, 571, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2014, A&A, 571, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2014, A&A, 571, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2014, A&A, 571, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIII. 2014, A&A, 571, A23 [NASA ADS] [CrossRef] [EDP Sciences] [MathSciNet] [PubMed] [Google Scholar]
 Planck Collaboration XXIV. 2014, A&A, 571, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXV. 2014, A&A, 571, A25 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVI. 2014, A&A, 571, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2014, A&A, 571, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVIII. 2014, A&A, 571, A28 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIX. 2014, A&A, 571, A29 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXXI. 2014, A&A, 571, A31 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Ratra, B. 1992, ApJ, 391, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Rees, M. J., & Sciama, D. W. 1968, Nature, 217, 511 [NASA ADS] [CrossRef] [Google Scholar]
 Regan, D., Mukherjee, P., & Seery, D. 2013, Phys. Rev. D, 88, 043512 [NASA ADS] [CrossRef] [Google Scholar]
 Reinecke, M., Dolag, K., Hell, R., Bartelmann, M., & Enßlin, T. A. 2006, A&A, 445, 373 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 RenauxPetel, S. 2009, J. Cosmol. Astropart. Phys., 10, 12 [NASA ADS] [Google Scholar]
 RenauxPetel, S., Mizuno, S., & Koyama, K. 2011, J. Cosmol. Astropart. Phys., 11, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Rigopoulos, G. I., Shellard, E. P. S., & van Tent, B. J. W. 2006, Phys. Rev. D, 73, 083522 [NASA ADS] [CrossRef] [Google Scholar]
 Rigopoulos, G. I., Shellard, E. P. S., & van Tent, B. J. W. 2007, Phys. Rev. D, 76, 083512 [NASA ADS] [CrossRef] [Google Scholar]
 Rudjord, O., Hansen, F. K., Lan, X., et al. 2009, ApJ, 701, 369 [NASA ADS] [CrossRef] [Google Scholar]
 Sachs, R. K., & Wolfe, A. M. 1967, ApJ, 147, 73 [NASA ADS] [CrossRef] [Google Scholar]
 Salem, M. P. 2005, Phys. Rev. D, 72, 123516 [NASA ADS] [CrossRef] [Google Scholar]
 Salopek, D. S., & Bond, J. R. 1990, Phys. Rev. D, 42, 3936 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Sasaki, M., Väliviita, J., & Wands, D. 2006, Phys. Rev. D, 74, 3003 [Google Scholar]
 Scodeller, S., Rudjord, O., Hansen, F. K., et al. 2011, ApJ, 733, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Sefusatti, E., Fergusson, J. R., Chen, X., & Shellard, E. P. S. 2012, J. Cosmol. Astropart. Phys., 8, 33 [NASA ADS] [CrossRef] [Google Scholar]
 Senatore, L., & Zaldarriaga, M. 2011, J. Cosmol. Astropart. Phys., 1, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Senatore, L., Tassev, S., & Zaldarriaga, M. 2009, J. Cosmol. Astropart. Phys., 9, 38 [Google Scholar]
 Senatore, L., Smith, K. M., & Zaldarriaga, M. 2010, J. Cosmol. Astropart. Phys., 1, 28 [Google Scholar]
 Serra, P., & Cooray, A. 2008, Phys. Rev. D, 77, 107305 [NASA ADS] [CrossRef] [Google Scholar]
 Shiraishi, M., Nitta, D., Yokoyama, S., & Ichiki, K. 2012, J. Cosmol. Astropart. Phys., 3, 41 [NASA ADS] [CrossRef] [Google Scholar]
 Shiraishi, M., Komatsu, E., Peloso, M., & Barnaby, N. 2013, J. Cosmol. Astropart. Phys., 5, 2 [Google Scholar]
 Shiu, G., & Xu, J. 2011, Phys. Rev. D, 84, 3509 [CrossRef] [Google Scholar]
 Silverstein, E., & Tong, D. 2004, Phys. Rev. D, 70, 3505 [Google Scholar]
 Smidt, J., Amblard, A., Byrnes, C. T., et al. 2010, Phys. Rev. D, 81, 3007 [NASA ADS] [Google Scholar]
 Smith, K. M., & Zaldarriaga, M. 2011, MNRAS, 417, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. M., Senatore, L., & Zaldarriaga, M. 2009, J. Cosmol. Astropart. Phys., 9, 6 [Google Scholar]
 Smith, K. M., Loverde, M., & Zaldarriaga, M. 2011, Phys. Rev. Lett., 107, 191301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Smith, T. L., & Kamionkowski, M. 2012, Phys. Rev. D, 86, 063009 [NASA ADS] [CrossRef] [Google Scholar]
 Starck, J.L., Moudden, Y., Abrial, P., & Nguyen, M. 2006, A&A, 446, 1191 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Starck, J.L., Murtagh, F., & Fadili, J. M. 2010, Sparse image and signal processing, wavelets, curvelets, morphological diversity (Cambridge: Cambridge University Press) [Google Scholar]
 Su, S.C., Lim, E. A., & Shellard, E. P. S. 2012 [arXiv:1212.6968] [Google Scholar]
 Sugiyama, N. S. 2012, J. Cosmol. Astropart. Phys., 5, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Sugiyama, N. S., Komatsu, E., & Futamase, T. 2011, Phys. Rev. Lett., 106, 1301 [CrossRef] [PubMed] [Google Scholar]
 Suyama, T., & Yamaguchi, M. 2008, Phys. Rev. D, 77, 023505 [NASA ADS] [CrossRef] [Google Scholar]
 Suyama, T., Takahashi, T., Yamaguchi, M., & Yokoyama, S. 2010, J. Cosmol. Astropart. Phys., 12, 30 [NASA ADS] [CrossRef] [Google Scholar]
 Toffolatti, L., Argueso Gomez, F., de Zotti, G., et al. 1998, MNRAS, 297, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Tolley, A. J., & Wyman, M. 2010, Phys. Rev. D, 81, 3502 [NASA ADS] [CrossRef] [Google Scholar]
 Tzavara, E., & van Tent, B. 2011, J. Cosmol. Astropart. Phys., 6, 26 [Google Scholar]
 Verde, L., & Spergel, D. N. 2002, Phys. Rev. D, 65, 3007 [Google Scholar]
 Verde, L., Wang, L., Heavens, A. F., & Kamionkowski, M. 2000, MNRAS, 313, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Verde, L., Jimenez, R., AlvarezGaume, L., Heavens, A. F., & Matarrese, S. 2013, J. Cosmol. Astropart. Phys., 6, 23 [Google Scholar]
 Vernizzi, F., & Wands, D. 2006, J. Cosmol. Astropart. Phys., 5, 19 [Google Scholar]
 Wang, L., & Kamionkowski, M. 2000, Phys. Rev. D, 61, 063504 [NASA ADS] [CrossRef] [Google Scholar]
 Weinberg, S. 2008, Phys. Rev. D, 77, 123541 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Yadav, A. P. S., & Wandelt, B. D. 2008, Phys. Rev. Lett., 100, 181301 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, A. P. S., & Wandelt, B. D. 2010, Adv. Astron., 2010, 71 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, A. P. S., Komatsu, E., & Wandelt, B. D. 2007, ApJ, 664, 680 [NASA ADS] [CrossRef] [Google Scholar]
 Yadav, A. P. S., Komatsu, E., Wandelt, B. D., et al. 2008, ApJ, 678, 578 [NASA ADS] [CrossRef] [Google Scholar]
 Zaldarriaga, M. 2004, Phys. Rev. D, 69, 043508 [NASA ADS] [CrossRef] [Google Scholar]
 Zeldovich, Y. B., & Sunyaev, R. A. 1969, Ap&SS, 4, 301 [NASA ADS] [CrossRef] [Google Scholar]
Appendix A: Full focal plane simulations
The purpose of the full focal plane (FFP) simulations is to provide a complete, coherent realization of the full Planck mission (both HFI and LFI) using the best available estimates of its characteristics, including the satellite pointing, the individual detector beams, band passes and noise properties, and the various data flags. Each FFP data set consists of three parts: a fiducial observation of an input sky (including our best estimates of all of the foreground components and a realization of the CMB sky corresponding to a chosen cosmology), together with separate MC realization sets of the CMB and noise. The maps made from these simulated data are used both to validate and verify the analysis algorithms and their implementations and to quantify uncertainties in and remove biases from the analysis of the real data. The FFP simulations have been used for validating mapmaking, component separation (for both CMB and diffuse foreground recovery), power spectra estimation, cosmological parameter estimation and measurements of primordial NG and lensing. The complete specification of any FFP simulation consists of the definition of the inputs (mission characteristics and sky), the outputs (timeordered data and maps), and the processes used to generate the latter from the former. For this data release we have used the sixth FFP simulation suite, henceforth referred to as FFP6.
Appendix A.1: Mission characteristics
The goal of FFP6 is to replicate the 2013 Planck data release (DR1) as closely as possible. The satellite pointing incorporates the PTCOR6 wobble corrections used by HFI (Planck Collaboration VI 2014), and the focal plane geometry of each instrument is provided by the corresponding DPC in the form of a reduced instrument model (RIMO) (Planck Collaboration II 2014, and Planck Collaboration VI 2014). The overall data flags are a combination of the pointing and detector flags. For both instruments, satellite repointing maneuvers and planet crossings are flagged. For pairs of detectors in the same horn, each sample flagged in one detector is flagged in the other member of the pair.
In addition to the geometric detector offsets, the RIMOs for both instruments provide a representation of the detector bandpasses. For LFI, the RIMO also provides a parameterized noise model (white noise level, knee frequency and power law index) for each detector (Planck Collaboration II 2014). LFI beams were measured from Jupiter scans, fitted using Gaussian circular and elliptical approximations, converted into a polarized beam model, and smeared to account for the satellite motion (Planck Collaboration IV 2014). For HFI, the noise power spectral densities (PSDs) are averaged and smoothed versions of the raw ringbyring spectra determined by the HFI pipeline. This processing produced the per detector mean noise spectrum over the duration of the nominal mission, calibrated to convert from preprocessed time ordered information (TOI) units to thermodynamic K (Planck Collaboration VI 2014). HFI beam maps were synthesized from the elliptical beam parameters and GaussHermite coefficients from the Mars observation (Planck Collaboration VII 2014).
Appendix A.2: Input sky
The FFP6 input sky is generated using the Planck Sky Model (PSM; Delabrouille et al. 2013) and includes the CMB together with the current best estimate foreground templates for the cosmic infrared background, CO line emission, freefree, synchrotron, spinning dust, thermal dust, kinetic and thermal SZ, and radio and infrared point sources. For simulation purposes, point sources are considered strong if their flux is more than 100 mJy at any of 30, 100, 300 or 1000 GHz; strong point sources are provided in catalogs, whilst weak sources are mapped. All maps are provided in thermodynamic K at HEALPix resolution 2048 (Górski et al. 2005) and smoothed with a 4 arcmin beam, while catalogues are in Jy. The CMB is generated as a lensed sum of a scalar, tensor and nonGaussian realization with particular values of the tensor to scalar ratio (r) and NG parameter (f_{NL}). For every detector, each foreground component and catalogue of strong point sources also takes the appropriate bandpassed from the RIMO and applies it.
Appendix A.3: FFP6 outputs
The FFP6 simulation consists of maps of the fiducial sky together with sets of MC realizations of the CMB and of the noise, each of which is produced using a distinct processing pipeline; full details on these can be found in the Explanatory Supplement (Planck Collaboration 2013).
The fiducial sky maps are made from explicit timeordered data sets which are generated using the Level S software suite (Reinecke et al. 2006) for every detector and for every component; this results in 770 timeordered datasets occupying 17TB of disk space. For every component we then make maps of the data at each frequency using the Madam/TOAST mapmaker (Keihänen et al. 2010) for all the detectors, all of the detector quadruplets and all of the unpolarized detector singlets, and using both the full data and the two halfring data subsets. In addition we make maps of the sum of all of the components’ time streams with simulated noise added. HFI data are mapped at HEALPix resolution 2048 and LFI data at 1024; this results in 2370 maps occupying 0.4TB of disk space. We also produced hit maps, binned maps and the 3 × 3 block diagonal whitenoise covariance matrices and their inverses.
The CMB MC provides a set of 1000 temperatureonly maps of the full data combining all of the detectors at each frequency. These are made directly in the map domain by calculating the effective beam at each frequency using the FEBeCoP software (Mitra et al. 2011) and applying it to each of 1000 realizations of the CMB sky drawn from the bestfit WMAP spectrum and generated using the HEALPix synfast tool (Górski et al. 2005); this results in 9000 maps occupying 1.2TB of disk space.
The MC noise maps are generated using Madam/TOAST, which has the ability to simulate the noise time stream on the fly (bypassing the otherwise cripplingly expensive writing and reading of timeordered data objects). For every detector combination at each frequency, and for both the full data and halfring data subsets, we generate 1000 noise realizations; this results in 222 000 maps occupying 33TB of disk space.
The bulk of the FFP simulations were performed at the US Department of Energy’s National Energy Research Scientific Computing Center (NERSC); the LFI noise MC for the various detector quadruplets were produced at the Finnish Supercomputing Center (CSC); in total the generation of FFP6 required some 15 million CPUhours and 50TB of disk space. All of the temperature data in the FFP6 simulation set are available for public use at NERSC.
Appendix A.4: Validation and exploitation
Specific tests for assessing our confidence in the claims were conducted for all the main results in the Planck 2013 release, namely CMB reconstruction through foreground cleaning, likelihood to cosmology, lensing and primordial NG; these are discussed in detail in the corresponding papers (Planck Collaboration XII 2014; Planck Collaboration XV 2014; Planck Collaboration XVI 2014; Planck Collaboration XVII 2014; Planck Collaboration XXIII 2014 and this work). The validation was conducted with the necessary accuracy for measuring secondorder cosmological effects like in the case of primordial NG as well as lensing. In the former case, blind tests on an FFP set with a nonzero and detectable value in the CMB component were performed and the correct value (within the error bars) was detected after foreground cleaning. In the case of lensing, null detection tests were conducted successfully after the process of foreground cleaning on maps that did not contain this signal as input.
Residual systematics which were not taken into account belong to the lowℓ and highℓ regimes for LFI and HFI, respectively. Although the latter were not translated directly into spurious NG effects, their effect was quantified. For LFI, sidelobe corrections have been quantified to a level which is comparable or below 0.05% in the beam transfer function on the range of ℓ between 200 and 1500, being increasingly subdominant with respect to noise at high ℓs. A number of HFI systematic effects and their treatment in the data processing are not included in the simulations. On large angular scales, the most important effects are the Analogue Digital Correction (ADC) nonlinearities (which principally manifest themselves as an apparent gain drift), Zodiacal light emission, and far sidelobe pickup of the Galaxy. On small angular scales the most important effects are cosmic ray hits (flagging of those from real data was taken into account into the FFP6 processing) and the 4HeJT (4K) cooler lines. The processing steps used to remove these effects and assessments of the residuals are described in Planck Collaboration VI (2014); Planck Collaboration VIII (2014); Planck Collaboration X (2014); Planck Collaboration XIV (2014).
Planck will be updating its suite of simulations for the forthcoming release including polarization. The instrumental characteristics described above will be updated, and particular care is being taken regarding the ability to simulate and control the main systematic effects affecting the polarised signal.
Appendix B: Expected level of agreement from bispectrum estimators with correlated weights
The estimator crossvalidation work presented in Sect. 6.1 was based on comparing results from different estimators using sets (typically 50 to 100 simulations in size) of both Gaussian and nonGaussian simulations. We started from idealized maps (e.g., fullsky, noiseless, Gaussian simulations) and then went on to include an increasing number of realistic features at each additional validation step. This allowed better testing and characterization of the response of different pipelines and bispectrum decompositions to various potential spurious effects in the data. As a preliminary step, we derived a general formula describing the expected level of agreement between estimators with different but strongly correlated weights, with the simplifying assumption of fullsky measurements and homogeneous noise. This theoretical expectation, summarized by Eq. (86), was then used as a benchmark against which to assess the quality of the results. The aim of this appendix is to describe in detail how we obtained Eq. (86).
Let us consider the idealized case of fullsky, noiseless CMB measurements (note that the following conclusions also work for homogeneous noise, because the pure cubic f_{NL} estimators without linear term corrections are still optimal) . Under these assumptions, the f_{NL} estimator for a given CMB shape B_{ℓ1ℓ2ℓ3} can be written simply as (see Sect. 3 for details): (B.1)where is the angleaveraged bispectrum template for a given theoretical shape, a_{ℓ1m1}a_{ℓ2m2}a_{ℓ3m3} is a bispectrum estimate constructed from the data, and F is the normalization of the estimator, provided by the Fisher matrix number (B.2)As explained in Sect. 3, the different f_{NL} estimation techniques implemented in this work can be seen as separate implementations of the optimal estimator of Eq. (B.1). Each implementation is based on expanding the angular bispectrum as a sum of basis templates defined in different domains: for example in our analyses we built templates out of products of wavelets at different scales, cubic combinations of 1dimensional polynomials and plane waves (what we call “modal estimator” in the main text), and ℓbinning of the bispectrum (what we called “binned estimator” in the main text). Our initial theoretical templates in this work are the local, equilateral and orthogonal separable cases used in the KSW and skewC_{ℓ} estimators. In this sense the KSW/skewC_{ℓ} estimator provides an “exact” estimate of f_{NL} for this choice of shapes, while the other pipelines provide an approximate result that approaches KSW measurements as the expansions get more accurate. The differences between results from different pipelines became smaller and smaller as the approximate expanded templates converge to the starting one (e.g., by increasing the number of ℓbins or the number of wavelets and polynomial triplets). The degree of convergence can be measured through the correlation coefficient r between the initial bispectrum and its reconstructed expanded version. The correlation coefficient is, as usual, defined as (B.3)where the label “th” denotes the initial bispectrum shape to fit to the data, and “exp” is the approximate expanded one. The correlator between shapes naturally defines a scalar product in bispectrum space, and in the following the operation of correlating two shapes will be denoted by the symbol ⟨ , ⟩, so that the definition above would read .
Whichever basis and separation scheme we have chosen, let us call ℛ_{n}(ℓ_{1},ℓ_{2},ℓ_{3}) the nth template in the adopted bispectrum expansion, and α_{n} the coefficients of the expansion, so that we can write: (B.4)From now on we will also call ℛ_{n}(ℓ_{1},ℓ_{2},ℓ_{3}) the “modes” of our expansion, N_{exp} is the number of modes we are using to approximate the starting template in order to obtain a correlation coefficient r. We will also assume that the modes form an orthonormal basis, that is: (B.5)where is the Kronecker delta symbol. The orthogonality assumption does not imply loss of generality since any starting set of modes can always be rotated and orthogonalized. We now consider an expansion with a number of coefficients N_{th}>N_{exp} that perfectly reproduces the initial bispectrum (i.e., r = 1), and is characterized by the same modes and coefficients as the previous one up to N_{exp}. So we can write (B.6)We now build two optimal estimators of f_{NL} for the shape B_{th}: an “exact” estimator and an “approximate” one. In the exact estimator we fit the actual B_{th} shape to the data and obtain the estimate , while in the approximate estimator we fit the expanded shape B_{exp} to obtain . We want to understand by how much the “exact” and “approximate” f_{NL} measurements are expected to differ, as a function of the correlation coefficient r between the weights B_{th} and B_{exp} that appear in the two estimators.
For each mode template ℛ_{n}(ℓ_{1},ℓ_{2},ℓ_{3}) we can build an optimal estimator (following Eq. (B.1)) in order to fit the mode to the data and get a corresponding amplitude estimate β_{n}. In other words, the observed bispectrum can then be reconstructed as in Eq. (B.4), but with coefficients β_{n} instead of α_{n}. Due to the orthonormality of the ℛ_{n}, the β coefficients have unit variance. It is then possible to show (Fergusson et al. 2010a) that the f_{NL} estimate for a given B_{th} with expansion parameters α_{n} is given by (B.7)In the light of all the above, the exact and approximate estimators are: Thanks to the orthonormality properties of the modes, we can easily relate the Fisher matrix normalization F to the expansion coefficients α:
(B.10)In analogous fashion we can derive an expression for the correlation coefficient r between the two estimators we are comparing. It is easy to check that the following holds: (B.11)and using the equation just derived above we can also write r^{2} = F_{exp}/F_{th}, i.e., the square of the correlation coefficients between the estimators equals the ratio of the normalizations.
Summary of results from the modal estimator survey of primordial models for the main nonstandard bispectrum shapes.
Crossvalidation of best fit feature model results for the SMICA, NILC and SEVEM foregroundcleaned maps.
Comparison of f_{NL} results for the hybrid polynomial and Fourier modes for a variety of nonseparable and feature models.
Armed with this preliminary notation we can now calculate the expected scatter between the exact and approximate f_{NL} measurement when we apply the two estimators to the same set of data. In order to quantify it we will consider the variable , and calculate its standard deviation. We find (B.12)where we made use of Eq. (B.11). It is easy to show that orthonormality of the ℛ templates implies no correlation of the amplitudes β, i.e., . A straightforward calculation then yields: (B.13)This, together with Eqs. (B.10, B.11) finally gives: (B.14)where Δ_{th} is the standard deviation of the exact estimator. Equation (B.14) is the starting point of our validation tests. It provides an estimate of the expected scatter between f_{NL} estimators with correlated weights, as a function of the f_{NL} error bars and of the correlation coefficient r. Note that this formula has been obtained under the simplifying assumptions of Gaussianity, fullsky coverage and homogeneous noise. For applications dealing with more realistic cases we might expect the scatter to become larger than this expectation, while remaining qualitatively consistent.
In our tests we started from sets of simulated maps and compared f_{NL} results from different pipelines and shapes on a mapbymap basis. In our validation tests the correlation levels between templates in different expansions were varying between r ~ 0.95 and r ~ 0.99, depending on the estimators and the shapes under study. Using the formula above, this corresponds to an expected scatter 0.15 Δ ≤ σ ≤ 0.3 Δ, where Δ is the f_{NL} standard deviation for the shape under study.
In Sect. 6.1 we presented several applications to simulated data sets, showing that the recovered results are fully consistent with these expectations, and thus the different pipelines are fully consistent with each other.
Appendix C: Targeted bispectrum constraints
This Appendix provides further tables of results for primordial models, extending those given in Sect. 7.3, notably for resonance models, while it also gives additional validation checks for the modal bispectrum estimator, beyond the extensive tests described already in the paper for the standard bispectrum shapes. In particular, using each of the SMICA NILC and SEVEM maps, we will quote results for the main paradigms for nonstandard bispectrum models, including comparisons for the feature model results. Remarkably consistent results are again obtained using the different foregroundseparation methodologies and using different modal basis functions.
Appendix C.1: Additional results for resonance models and enfolded resonance models
As described in Sect. 7.3, we have surveyed resonance models (Eq. (17)) in the region of most interest for feature models, that is, with comparable periodicities to those with described in Table 12 near the bestfit feature model. The results from this initial survey for the SMICA componentseparated map produced no significant signal, with the results Table C.1. We note that while the feature and resonance models proved similar for the WMAP analysis (Fergusson et al. 2012), wavelengths are much more differentiated for Planck and so it can be difficult to resolve the shortest wavelengths at low ℓ. This is the key limitation on surveys with the present estimator and will be circumvented in future, by using specific separable templates to improve the overall resolution. Just as local and ISW templates can be incorporated into the analysis, so can targeted feature modes. Note that consistent results were obtained using the NILC and SEVEM maps, though we only show results here for feature models with an envelope (see discussion below with Table C.4).
The enfolded resonant (or nonBunchDavies feature) models (Eq. (18)) were also surveyed for these periodicities and also yielded no significant signal; see Table C.2. These are interesting shapes which hold out the prospect of exhibiting both oscillatory and flattened features observed in the Planck bispectrum reconstruction, see Fig. 7. Due to similar resolution restrictions, only relatively large multipole periodicities (ℓ > 100) have been surveyed in the present work, again searching around the periodicities exhibited for feature models.
Appendix C.2: Comparison of f_{NL} results from SMICA, NILC and SEVEM foregroundcleaned maps
Tables C.3 and C.4 compare bispectrum results extracted from Planck maps created using the three different componentseparation techniques, SMICA, NILC and SEVEM. In addition to the models discussed in Sect. 7.3, we have also included the constant model which is defined by . The abbreviations EFT denotes the effective field theory single field shapes, NBD the nonBunchDavies (excited initial state) models, vector models are gauge field inflation with directional dependence, along with DBI, Ghost and Warm inflation, also described previously. The expression “cleaned” refers to removal of the predicted ISWlensing signal and the measured point source signal, unless stated otherwise.
We note that there is good consistency between the different foregroundseparation techniques for all models, whether equilateral, flattened, or squeezed as shown in Table C.3. For the nonscaling case, differences for the bestfit feature models were below 1/3σ confirming the interesting results discussed in Sect. 3.2.3, see Table C.4.
Appendix C.3: Comparison of f_{NL} results from ISW Fourier basis with hybrid polynomials
As described in Sect. 3.2.3, the modal bispectrum estimator can flexibly incorporate any suitable basis functions with which to expand the bispectrum and separably filter CMB maps (Fergusson et al. 2010a). For the Planck analysis, we have evolved two sets of basis functions to fulfil three basic criteria: first, to provide a complete basis for the bispectrum up to a given ℓresolution, secondly, to represent accurately primordial models of interest and, thirdly, to incorporate the CMB ISWlensing signal, which with diffuse point sources provides a significant secondary signal which must be subtracted. The first basis functions are n_{max} = 600 polynomials (closely related to Legendre polynomials) which are supplemented with the SachsWolfe local mode in order to represent more accurately the squeezed limit; enhanced orthogonality is preconditioned by choosing these from a larger set of polynomials. The second basis functions are n_{max} = 300 Fourier modes, augmented with the same SW local mode, together with the the separable ISWlensing modes. Both modal expansions proved useful, providing important validation and crosschecks, however, the twofold resolution improvement from the polynomials meant that most quantitative results employed this basis. This improved resolution was particularly important in probing flattened models on the edge of the tetrapyd.
In Fig. 7, we show a direct comparison between the modal reconstruction of the 3D bispectrum using the polynomial and Fourier mode expansions. The basic features of the two reconstructions are in good agreement, confirming a central feature which changes sign at low ℓ and a flattened signal beyond as discussed previously in Sect. 7.2. Notably the polynomial basis, with double the resolution, preserves the largescale features observed in the Fourier basis.
In Table C.5, results from both basis expansions are shown for a variety of nonseparable models. These demonstrate consistent results where the Fourier basis had sufficient resolution, as indicated by the ratio of the variance reflecting the Fisher ratio (i.e., above 90% correlation as indicated by the results in Appendix B). Note that we independently determine the estimator correlation between the exact solution and primordial decomposition and then at late times with the CMB decomposition; we use a polynomial basis as the overall benchmark here. This analysis also includes several feature models (phase φ = 0) showing good agreement from the Fourier basis while the Fisher estimates remain accurate. Again, the hybrid Fourier basis degrades in accuracy towards k = 0.02 as it reaches its resolution limit, when the variance disparity rises towards 10%. With n_{max} = 600 modes and ℓ_{max} = 2000, the polynomial basis retains a good correlation for all primordial feature models for k> 0.01. The accuracy and robustness of the feature model results have been verified using ℓ_{max} = 1500 for the polynomial expansion, for example, obtaining 3.1σ with the SMICA map for the bestfit model (K = 0.01875,φ = 0).
All Tables
The bias in the three primordial f_{NL} parameters due to the ISWlensing signal for the four componentseparation methods.
Results for the amplitude of the ISWlensing bispectrum from the SMICA, NILC, SEVEM, and CR foregroundcleaned maps, for the KSW, binned, and modal (polynomial) estimators; error bars are 68% CL.
Results for the amplitude of the point source (Poisson) bispectrum (in dimensionless units of 10^{29}) from the SMICA, NILC, SEVEM, and CR foregroundcleaned maps, for the KSW, binned, and modal (polynomial) estimators; error bars are 68% CL.
Results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW estimator from the SMICA foregroundcleaned map. Both independent singleshape results and results marginalized over the point source bispectrum and with the ISWlensing bias subtracted are reported; error bars are 68% CL.
Results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW, binned and modal estimators from the SMICA, NILC, SEVEM, and CR foregroundcleaned maps.
Results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the suboptimal wavelet estimator from the SMICA foregroundcleaned map.
Constraints on flattened or collinear bispectrum models (and related models) using the SMICA foregroundcleaned Planck map.
Planck bispectrum estimation results for feature models compared to the SMICA foregroundcleaned maps.
Validation tests with MFs: results for obtained using the filter W_{M}, for ℓ_{max} = 2000 and N_{side} = 2048.
Results for f_{NL} (assumed independent, without any correction for the ISWlensing bias) of the SMICA cleaned map using different values of ℓ_{max}, for the KSW and binned estimators.
Results for f_{NL} (assumed independent) for the raw frequency maps at 70, 100, 143, and 217 GHz with a very large mask (f_{sky} = 0.32) compared to the SMICA result with the union mask U73 (f_{sky} = 0.73), as determined by the binned (with ℓ_{max} = 2500) and modal (with ℓ_{max} = 2000) estimators.
Results for f_{NL} (assumed independent) of the SMICA halfring null maps, determined by the KSW, binned and modal estimators.
Results for f_{NL} (assumed independent) of several null maps determined by the binned estimator.
Summary of results from the modal estimator survey of primordial models for the main nonstandard bispectrum shapes.
Crossvalidation of best fit feature model results for the SMICA, NILC and SEVEM foregroundcleaned maps.
Comparison of f_{NL} results for the hybrid polynomial and Fourier modes for a variety of nonseparable and feature models.
All Figures
Fig. 1
Permitted observational domain of Eq. (28) for the CMB bispectrum b_{ℓ1ℓ2ℓ3}. Allowed multipole values (ℓ_{1},ℓ_{2},ℓ_{3}) lie inside the shaded “tetrapyd” region (tetrahedron+pyramid), satisfying both the triangle condition and the experimental resolution ℓ <L ≡ ℓ_{max}. 

In the text 
Fig. 2
Binned skewC_{ℓ} statistics from the SMICA map for a) ISWlensing and b) Poisson point sources. Theoretical curves are not fitted to the data shown, but are plotted with the amplitude (the only free parameter) determined from the KSW technique. The Poisson pointsource foreground is clearly detected, and the ISWlensing skewspectrum is evident and consistent with the overall 2.6σ detection. b_{ps} is the Poisson pointsource amplitude in dimensionless units of 10^{29}, and is the ISWlensing amplitude in units of that expected from the Planck bestfit cosmology. Error bars come from covariance estimates from 1000 simulated maps, and the points are mildly correlated. 

In the text 
Fig. 3
Mapbymap comparison of the results from the different estimators for local (top), equilateral (centre), and orthogonal (bottom) f_{NL} for the set of masked nonGaussian simulations described in Sect. 6.1.3, assuming the shapes to be independent. The horizontal solid line is the average value of all maps for KSW, and the dashed and dotted horizontal lines correspond to 1σ and 2σ deviations, respectively. 

In the text 
Fig. 4
Mapbymap comparison of the results from the different estimators for local (top), equilateral (centre), and orthogonal (bottom) f_{NL} for 99 maps from a set of realistic lensed simulations passed through the SMICA pipeline, described in Sect. 6.2, assuming the shapes to be independent. The horizontal solid line is the average value of the maps for KSW, and the dashed and dotted horizontal lines correspond to 1σ and 2σ deviations, respectively. 

In the text 
Fig. 5
Binned skewC_{ℓ} statistics from the SMICA map for a) local; b) equilateral; and c) orthogonal. Theoretical curves are not fitted to the data shown, but are plotted with the amplitude (the only free parameter) determined from the KSW technique. There is no evidence for detection of primordial NG. Error bars are derived from the covariance of estimates from 1000 simulations. There are mild correlations between data points in all figures, but very strong correlations at high ℓ in the local case, reaching correlation coefficients r> 0.99 for ℓ > 1750. 

In the text 
Fig. 6
Full 3D CMB bispectrum recovered from the Planck foregroundcleaned maps, including SMICA (left), NILC (centre) and SEVEM (right), using the hybrid Fourier mode coefficients illustrated in Fig. 9, These are plotted in threedimensions with multipole coordinates { ℓ_{1},ℓ_{2},ℓ_{3} } on the tetrahedral domain shown in Fig. 1 out to ℓ_{max} = 2000. Several density contours are plotted with red positive and blue negative. The bispectra extracted from the different foregroundseparated maps appear to be almost indistinguishable. 

In the text 
Fig. 7
Planck CMB bispectrum detail in the signaldominated regime showing a comparison between full 3D reconstruction using hybrid Fourier modes (left) and hybrid polynomials (right). Note the consistency of the main bispectrum properties which include an apparently “oscillatory” central feature for lowℓ together with a flattened signal beyond to ℓ ≲ 1400. Note also the periodic CMB ISWlensing signal in the squeezed limit along the edges of the tetrapyd. 

In the text 
Fig. 8
CMB bispectrum from a Gaussian simulation including gravitational lensing. This reconstruction uses the same hybrid polynomials as for the Planck bispectrum in Fig. 7 (right), with which it can be compared. Note that an indication of a ISWlensing bispectrum signal can be seen along the edges of the tetrahedron. 

In the text 
Fig. 9
Modal bispectrum coefficients for the mode expansion (Eq. (65)) obtained from Planck foregroundcleaned maps using hybrid Fourier modes. The different component separation methods, SMICA, NILC and SEVEM exhibit remarkable agreement. The variance from 200 simulated noise maps was nearly constant for each of the 300 modes, with the average ±1σ variation shown in red. 

In the text 
Fig. 10
Total integrated bispectrum defined in Eq. (66) as a cumulative sum over orthonormal modal coefficients (upper panel) and over multipoles up to a given ℓ (lower panel). Above, the relative quantity is plotted, where is the mean obtained from 200 CMB Gaussian maps with the standard deviation shown as the red line. Below the square of the bispectrum is integrated over the tetrapyd out to ℓ and its significance plotted relative to the Gaussian standard deviation (1σ red line). A hybrid polynomial basis n_{max} = 600 is employed in the signaldominated region ℓ ≤ 1500. 

In the text 
Fig. 11
Smoothed observed bispectrum as determined with the binned estimator divided by its expected standard deviation, as a function of ℓ_{1} and ℓ_{2}, with ℓ_{3} in the bin [610, 654]. From left to right on the top row are shown: SMICA, NILC, and SEVEM; and on the bottom row: CR and the raw 143 GHz channel. For comparison purposes the last figure on the bottom row shows the same quantity for one of the lensed Gaussian simulations described in Sect. 6.2. 

In the text 
Fig. 12
Similar to Fig. 11, but with ℓ_{3} in the bin [1330, 1374]. 

In the text 
Fig. 13
CMB bispectrum shown for the bestfit feature model with an envelope with parameters k = 0.01875, phase φ = 0 and Δk = 0.045 (see Table 13). Compare with the Planck bispectrum reconstruction, Fig. 7. 

In the text 
Fig. 14
The Wiener filter W_{M} used to constrain with MFs. 

In the text 
Fig. 15
MFs curves for SMICA at N_{side} = 1024 and ℓ_{max} = 2000, for the four functionals v_{k}: a) area, b) perimeter, c) genus, and d) N_{cluster}. The curves are the difference of each normalized MF, measured from the data, to the average from Gaussian Planck realistic simulations (not lensed). The difference curves are normalized by the maximum of the Gaussian curve. To compare the curves to the presence of primordial NG, the average difference curves for nonGaussian simulations with is also represented (100 simulations). 

In the text 
Fig. 16
The two upper maps show the modulation reconstruction mean field at L ≤ 100, which is essentially a map of the expected total smallscale power on the masked map as a function of position (assuming there is no primordial power modulation). The top mean field map from the 143 GHz auto estimator has a large signal from both the cut (which can be calculated accurately), and from the noise anisotropy (aligned roughly with the ecliptic, which cannot be estimated very accurately from simulations). The lower mean field is the 143 × 217 GHz crossestimator map, and does not have the contribution from the noise anisotropy (note the colour scale is adjusted). The lower plot shows the corresponding mean field power spectra compared to the reconstruction noise (connected part of the trispectrum); the reconstruction noise is much smaller than both the detector noise and mask contributions to the mean field. Since any τ_{NL} signal is all on large scales we do not reconstruct modes above L_{max} = 100. 

In the text 
Fig. 17
Power spectrum of the power modulation reconstructed from 143 × 217 GHz maps. Shading shows the 68%, 95% and 99% CL intervals from simulations with no modulation or kinematic signal. The dashed lines are when the mean field simulations include no kinematic effects, showing a clear detection of a modulation dipole. The blue points show the expected kinematic modulation dipole signal from simulations, along with 1σ error bars (only first four points shown for clarity). The solid line subtracts the dipolar kinematic signal in the mean fields from simulations including the expected signal, and represents out best estimate of the nonkinematic signal (note this is not just a subtraction of the power spectra since the mean field takes out the fixed dipole anisotropy in real space before calculating the remaining modulation power). The dotted line shows the expected signal for τ_{NL} = 1000. 

In the text 
Fig. 18
Distribution of estimators from Gaussian simulations (L_{max} = 10) compared to data estimates (vertical lines). The distribution is rather skewed because the main contributions are from L ≲ 4 where the modulation power spectra have skewed χ^{2} distributions with low degrees of freedom. The red line shows the predicted distribution for a weighted sum of estimators assuming the reconstructed modulation modes are Gaussian with 2L + 1 modes measured per L, which fits the full simulations well. The black vertical lines show the data estimates from L_{max} = 10, and should be compared to the green which have L_{max} = 2 and hence are insensitive to the anomalous octopole signal. The dashed lines are τ_{NL,1}, the slightly more optimal variant of the estimator. 

In the text 
Fig. 19
Approximate posterior distributions for a range of L_{max}. The distributions have broad tails to high values because of the small number of largescale modulation modes that are measured, and hence large cosmic variance. For L_{max} ≥ 3 the distributions are pulled away from zero by the significant octopole modulation signal observed, and are gradually move back towards zero as we include more modulation modes that are inconsistent with large τ_{NL} values. As shown in Fig. 20 the octopole has significant frequency dependence and is therefore unlikely to be physical. 

In the text 