Issue 
A&A
Volume 641, September 2020
Planck 2018 results



Article Number  A9  
Number of page(s)  47  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201935891  
Published online  11 September 2020 
Planck 2018 results
IX. Constraints on primordial nonGaussianity
^{1}
AIM, CEA, CNRS, Université ParisSaclay, Université ParisDiderot, Sorbonne Paris Cité, 91191 GifsurYvette, France
^{2}
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
^{3}
African Institute for Mathematical Sciences, 68 Melrose Road, Muizenberg, Cape Town, South Africa
^{4}
Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{5}
Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{6}
CITA, University of Toronto, 60 St. George St., Toronto, ON M5S 3H8, Canada
^{7}
CNRS, IRAP, 9 Av. colonel Roche, BP 44346, 31028 Toulouse cedex 4, France
^{8}
California Institute of Technology, Pasadena, CA, USA
^{9}
Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
^{10}
Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, CA, USA
^{11}
Département de Physique Théorique, Université de Genève, 24 quai E. Ansermet, 1211 Genève 4, Switzerland
^{12}
Département de Physique, École normale supérieure, PSL Research University, CNRS, 24 rue Lhomond, 75005 Paris, France
^{13}
Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{14}
Departamento de Física, Universidad de Oviedo, C/ Federico García Lorca, 18, Oviedo, Spain
^{15}
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 GL Nijmegen, The Netherlands
^{16}
Department of General Education, National Institute of Technology, Kagawa College, 355 Chokushicho, Takamatsu, Kagawa 7618058, Japan
^{17}
Department of Mathematics, University of Stellenbosch, Stellenbosch 7602, South Africa
^{18}
Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, British Columbia, Canada
^{19}
Department of Physics & Astronomy, University of the Western Cape, Cape Town 7535, South Africa
^{20}
Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{21}
Department of Physics, Gustaf Hällströmin katu 2a, University of Helsinki, Helsinki, Finland
^{22}
Department of Physics, Princeton University, Princeton, NJ, USA
^{23}
Department of Physics, University of California, Santa Barbara, CA, USA
^{24}
Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy
^{25}
Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{26}
Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, Roma, Italy
^{27}
Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria, 16, Milano, Italy
^{28}
Dipartimento di Fisica, Università degli Studi di Trieste, Via A. Valerio 2, Trieste, Italy
^{29}
Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy
^{30}
European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, Madrid, Spain
^{31}
European Space Agency, ESTEC, Keplerlaan 1, 2201 AZ Noordwijk, The Netherlands
^{32}
Gran Sasso Science Institute, INFN, viale F. Crispi 7, 67100 L’Aquila, Italy
^{33}
HEP Division, Argonne National Laboratory, Lemont, IL 60439, USA
^{34}
Haverford College Astronomy Department, 370 Lancaster Avenue, Haverford, PA, USA
^{35}
Helsinki Institute of Physics, Gustaf Hällströmin katu 2, University of Helsinki, Helsinki, Finland
^{36}
INAF – OAS Bologna, Istituto Nazionale di Astrofisica – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Area della Ricerca del CNR, Via Gobetti 101, 40129 Bologna, Italy
^{37}
INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, Padova, Italy
^{38}
INAF – Osservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, Trieste, Italy
^{39}
INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy
^{40}
INAF/IASF Milano, Via E. Bassini 15, Milano, Italy
^{41}
INFN – CNAF, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{42}
INFN, Sezione di Bologna, viale Berti Pichat 6/2, 40127 Bologna, Italy
^{43}
INFN, Sezione di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{44}
INFN, Sezione di Milano, Via Celoria 16, Milano, Italy
^{45}
INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy
^{46}
INFN, Sezione di Roma 2, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy
^{47}
Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK
^{48}
Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay cedex, France
^{49}
Institut d’Astrophysique de Paris, CNRS (UMR7095), 98bis boulevard Arago, 75014 Paris, France
^{50}
Institute Lorentz, Leiden University, PO Box 9506, 2300 RA Leiden, The Netherlands
^{51}
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{52}
Institute of Theoretical Astrophysics, University of Oslo, Blindern, Oslo, Norway
^{53}
Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, Tenerife, Spain
^{54}
Instituto de Astrofísica e Ciências do Espaço, Faculdade de Ciências da Universidade de Lisboa, Campo Grande, 1749016 Lisboa, Portugal
^{55}
Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, Santander, Spain
^{56}
Istituto Nazionale di Fisica Nucleare, Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy
^{57}
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, USA
^{58}
Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester M13 9PL, UK
^{59}
Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{60}
Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), UTIAS, The University of Tokyo, Chiba 277 8583, Japan
^{61}
Laboratoire d’Océanographie Physique et Spatiale (LOPS), Univ. Brest, CNRS, Ifremer, IRD, Brest, France
^{62}
Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{63}
Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{64}
Low Temperature Laboratory, Department of Applied Physics, Aalto University, Espoo 00076, Finland
^{65}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{66}
Mullard Space Science Laboratory, University College London, Surrey RH5 6NT, UK
^{67}
NAOCUKZN Computational Astrophysics Centre (NUCAC), University of KwaZuluNatal, Durban 4000, South Africa
^{68}
National Centre for Nuclear Research, ul. L. Pasteura 7, 02093 Warsaw, Poland
^{69}
Perimeter Institute for Theoretical Physics, Waterloo, ON N2L 2Y5, Canada
^{70}
Purple Mountain Observatory, No. 8 Yuan Hua Road, 210034 Nanjing, PR China
^{71}
SISSA, Astrophysics Sector, Via Bonomea 265, 34136 Trieste, Italy
^{72}
San Diego Supercomputer Center, University of California, San Diego, 9500 Gilman Drive, La Jolla CA, 92093, USA
^{73}
School of Chemistry and Physics, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{74}
School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK
^{75}
School of Physics and Astronomy, Sun Yatsen University, 2 Daxue Rd, Tangjia, Zhuhai, PR China
^{76}
School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{77}
School of Physics, Indian Institute of Science Education and Research Thiruvananthapuram, Maruthamala PO, Vithura, Thiruvananthapuram, 695551 Kerala, India
^{78}
School of Physics, The University of New South Wales, Sydney, NSW 2052, Australia
^{79}
Simon Fraser University, Department of Physics, 8888 University Drive, Burnaby, BC, Canada
^{80}
Sorbonne Université, Observatoire de Paris, Université PSL, École normale supérieure, CNRS, LERMA, 75005 Paris, France
^{81}
Sorbonne Université, UMR7095, Institut d’Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
^{82}
Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, Moscow 117997, Russia
^{83}
Space Science Data Center – Agenzia Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy
^{84}
Space Sciences Laboratory, University of California, Berkeley, CA, USA
^{85}
The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{86}
Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex 4, France
^{87}
Van Swinderen Institute for Particle Physics and Gravity, University of Groningen, Nijenborgh 4, 9747 AG Groningen, The Netherlands
^{88}
Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
Received:
15
May
2019
Accepted:
5
November
2019
We analyse the Planck fullmission cosmic microwave background (CMB) temperature and Emode polarization maps to obtain constraints on primordial nonGaussianity (NG). We compare estimates obtained from separable templatefitting, binned, and optimal modal bispectrum estimators, finding consistent values for the local, equilateral, and orthogonal bispectrum amplitudes. Our combined temperature and polarization analysis produces the following final results: f_{NL}^{local} = −0.9 ± 5.1; f_{NL}^{equil} = −26 ± 47; and f_{NL}^{ortho} = −38 ± 24 (68% CL, statistical). These results include lowmultipole (4 ≤ ℓ < 40) polarization data that are not included in our previous analysis. The results also pass an extensive battery of tests (with additional tests regarding foreground residuals compared to 2015), and they are stable with respect to our 2015 measurements (with small fluctuations, at the level of a fraction of a standard deviation, which is consistent with changes in data processing). Polarizationonly bispectra display a significant improvement in robustness; they can now be used independently to set primordial NG constraints with a sensitivity comparable to WMAP temperaturebased results and they give excellent agreement. In addition to the analysis of the standard local, equilateral, and orthogonal bispectrum shapes, we consider a large number of additional cases, such as scaledependent feature and resonance bispectra, isocurvature primordial NG, and paritybreaking models, where we also place tight constraints but do not detect any signal. The nonprimordial lensing bispectrum is, however, detected with an improved significance compared to 2015, excluding the null hypothesis at 3.5σ. Beyond estimates of individual shape amplitudes, we also present modelindependent reconstructions and analyses of the Planck CMB bispectrum. Our final constraint on the local primordial trispectrum shape is g_{NL}^{local} = (−5.8 ± 6.5) × 10^{4} (68% CL, statistical), while constraints for other trispectrum shapes are also determined. Exploiting the tight limits on various bispectrum and trispectrum shapes, we constrain the parameter space of different earlyUniverse scenarios that generate primordial NG, including general singlefield models of inflation, multifield models (e.g. curvaton models), models of inflation with axion fields producing parityviolation bispectra in the tensor sector, and inflationary models involving vectorlike fields with directionallydependent bispectra. Our results provide a highprecision test for structureformation scenarios, showing complete agreement with the basic picture of the ΛCDM cosmology regarding the statistics of the initial conditions, with cosmic structures arising from adiabatic, passive, Gaussian, and primordial seed perturbations.
Key words: cosmic background radiation / cosmology: observations / cosmology: theory / early Universe / inflation / methods: data analysis
© ESO 2020
1. Introduction
This paper, one of a set associated with the 2018 release (also known as “PR3”) of data from the Planck^{1} mission (Planck Collaboration I 2020), resents the data analysis and constraints on primordial nonGaussianity (NG), which are obtained using the Legacy Planck cosmic microwave background (CMB) maps. It also includes some implications for inflationary models driven by the 2018 NG constraints. This paper updates the earlier study based on the temperature data from the nominal Planck operations period, including the first 14 months of observations (Planck Collaboration XXIV 2014, hereafter PCNG13), and a later study that used temperature data and a first set of polarization maps from the full Planck mission–29 and 52 months of observations for the High Frequency Instrument (HFI) and the Low Frequency Instrument (LFI), respectively (Planck Collaboration XVII 2016, hereafter PCNG15). The analysis described in this paper sets the most stringent constraints to date on primordial NG, which are near what is ultimately possible from using only CMB temperature data. The results of this paper are mainly based on the measurements of the CMB angular bispectrum, which are complemented with the next higherorder NG correlation function, that is the trispectrum. For notations and conventions relating to (primordial) bispectra and trispectra we refer the reader to the two previous Planck papers on primordial NG (PCNG13; PCNG15). This paper also complements the precise characterization of inflationary models (Planck Collaboration X 2020) and cosmological parameters (Planck Collaboration VI 2020), with specific statistical estimators that go beyond the constraints on primordial power spectra. It also complements the statistical and isotropy tests on CMB anisotropies of Planck Collaboration VII (2020), focusing on the interpretation of specific, well motivated, nonGaussian models of inflation. These models span from the irreducible minimal amount of primordial NG, predicted by standard singlefield models of slowroll inflation, to various classes of inflationary models that constitute the prototypes of extensions of the standard inflationary picture and of physically motivated mechanisms that are able to generate a higher level of primordial NG measurable in the CMB anisotropies. This work establishes the most robust constraints on some of the most wellknown and studied types of primordial NG, namely the local, equilateral, and orthogonal shapes. Moreover, this 2018 analysis includes a better characterization of the constraints coming from CMB polarization data. Besides focusing on these major goals, we reanalyse a variety of other NG signals and also investigate some new aspects of primordial NG. For example, we perform for the first time an analysis of the running of NG using Planck data in the context of some well defined inflationary models. Additionally, we constrained primordial NG predicted by theoretical scenarios on which much attention has been focused recently, such as, bispectrum NG generated in the tensor (gravitational wave) sector. For a detailed analysis of oscillatory features that combines power spectrum and bispectrum constraints see Planck Collaboration X (2020). As in the last data release (“PR2”), besides extracting the constraints on NG amplitudes for specific shapes, we also provide a modelindependent reconstruction of the CMB angular bispectrum by using various methods. Such a reconstruction can help to pin down interesting features in the CMB bispectrum signal beyond those captured by existing parameterizations.
The paper is organized as follows. In Sect. 2 we recall the main primordial NG models tested in this paper. Section 3 briefly describes the bispectrum estimators that we use, as well as details about the data set and our analysis procedures. In Sect. 4 we discuss detectable nonprimordial contributions to the CMB bispectrum, namely those arising from lensing and point sources. In Sect. 5 we constrain f_{NL} for the local, equilateral, and orthogonal bispectra. We also report the results for scaledependent NG models and other selected bispectrum shapes, including NG in the tensor (primordial gravitational wave) sector; in this section reconstructions and modelindependent analyses of the CMB bispectrum are also provided. In Sect. 6 these results are validated through a series of null tests on the data, with the goal of assessing the robustness of the results. This includes, in particular, a first analysis of Galactic dust and thermal SZ residuals. Planck CMB trispectrum limits are obtained and discussed in Sect. 7. In Sect. 8 we derive the main implications of Planck’s constraints on primordial NG for some specific early Universe models. We conclude in Sect. 9.
2. Models
Primordial NG comes in with a variety of shapes, corresponding to well motivated classes of inflationary model. For each class, a common physical mechanism is responsible for the generation of the corresponding type of primordial NG. Below we briefly summarize the main types of primordial NG that are constrained in this paper, providing the precise shapes that are used for data analysis. For more details about specific realizations of inflationary models within each class, see the previous two Planck papers on primordial NG (PCNG13; PCNG15) and reviews (e.g. Bartolo et al. 2004a; Liguori et al. 2010; Chen 2010b; Komatsu 2010; Yadav & Wandelt 2010). We only give a more expanded description of those shapes of primordial NG analysed here for the first time with Planck data (e.g. running of primordial NG).
2.1. General singlefield models of inflation
The parameter space of singlefield models is well described by the so called equilateral and orthogonal templates (Creminelli et al. 2006; Chen et al. 2007b; Senatore et al. 2010). The equilateral shape is
while the orthogonal NG is described by
Here the potential Φ is defined in relation to the comoving curvature perturbation ζ by Φ ≡ (3/5)ζ on superhorizon scales (thus corresponding to Bardeen’s gaugeinvariant gravitational potential (Bardeen 1980) during matter domination on superhorizon scales). P_{Φ}(k) = A/k^{4 − ns} is the Bardeen gravitational potential power spectrum, with normalization A and scalar spectral index n_{s}. A typical example of this class is provided by models of inflation where there is a single scalar field driving inflation and generating the primordial perturbations, characterized by a nonstandard kinetic term or more general higherderivative interactions. In the first case the inflaton Lagrangian is , where X = g^{μν}∂_{μ}ϕ ∂_{ν}ϕ, with at most one derivative on ϕ (Chen et al. 2007b). Different higherderivative interactions of the inflaton field characterize, ghost inflation (ArkaniHamed et al. 2004) or models of inflation based on Galileon symmetry (e.g., Burrage et al. 2011). The two amplitudes and usually depend on the sound speed c_{s} at which the inflaton field fluctuations propagate and on a second independent amplitude measuring the inflaton selfinteractions. The DiracBornInfeld (DBI) models of inflation (Silverstein & Tong 2004; Alishahiha et al. 2004) are a stringtheorymotivated example of the P(X, ϕ) models, predicting an almost equilateral type NG with for c_{s} ≪ 1. More generally, the effective field theory (EFT) approach to inflationary perturbations (Cheung et al. 2008; Senatore et al. 2010; Bartolo et al. 2010a) yields NG shapes that can be mapped into the equilateral and orthogonal template basis. The EFT approach allows us to draw generic conclusions about singlefield inflation. We discuss them using one example in Sect. 8. Nevertheless, we shall also explicitly search for such EFT shapes, analysing their exact nonseparable predicted shapes, B^{EFT1} and B^{EFT2}, along with those of DBI, B^{DBI}, and ghost inflation, B^{ghost} (ArkaniHamed et al. 2004).
2.2. Multifield models
The bispectrum for multifield models is typically of the local type^{2}
This usually arises when more scalar fields drive inflation and give rise to the primordial curvature perturbation (“multiplefield inflation”), or when extra light scalar fields, different from the inflaton field driving inflation, determine (or contribute to) the final curvature perturbation. In these models initial isocurvature perturbations are transferred on superhorizon scales to the curvature perturbations. NonGaussianities if present are transferred too. This, along with nonlinearities in the transfer mechanism itself, is a potential source of significant NG (Bartolo et al. 2002; Bernardeau & Uzan 2002; Vernizzi & Wands 2006; Rigopoulos et al. 2006, 2007; Lyth & Rodriguez 2005; Tzavara & van Tent 2011; Jung & van Tent 2017). The bispectrum of Eq. (3) mainly correlates large with smallscale modes, peaking in the “squeezed” configurations k_{1} ≪ k_{2} ≈ k_{3}. This is a consequence of the transfer mechanism taking place on superhorizon scales and thus generating a localized pointbypoint primordial NG in real space. The curvaton model (Mollerach 1990; Linde & Mukhanov 1997; Enqvist & Sloth 2002; Lyth & Wands 2002; Moroi & Takahashi 2001) is a clear example where local NG is generated in this way (e.g., Lyth & Wands 2002; Lyth et al. 2003; Bartolo et al. 2004c). In the minimal adiabatic curvaton scenario (Bartolo et al. 2004c,d), in the case when the curvaton field potential is purely quadratic (Lyth & Wands 2002; Lyth et al. 2003; Lyth & Rodriguez 2005; Malik & Lyth 2006; Sasaki et al. 2006). Here r_{D} = [3ρ_{curv}/(3ρ_{curv} + 4ρ_{rad})]_{D} represents the “curvaton decay fraction” at the epoch of the curvaton decay, employing the sudden decay approximation. Significant NG can be produced (Bartolo et al. 2004c,d) for low values of r_{D}; a different modelling of the curvaton scenario has been discussed by Linde & Mukhanov (2006) and Sasaki et al. (2006). We update the limits on both models in Sect. 8, using the local NG constraints. More general models with a curvatonlike spectator field have also been intensively investigated recently (see, e.g., Torrado et al. 2018). Notice that through a similar mechanism to the curvaton mechanism, local bispectra can be generated from nonlinear dynamics during the preheating and reheating phases (Enqvist et al. 2005; Chambers & Rajantie 2008; Barnaby & Cline 2006; Bond et al. 2009) or due to fluctuations in the decay rate or interactions of the inflaton field, as realized in modulated (p)reheating and modulated hybrid inflationary models (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). We also explore whether there is any evidence for dissipative effects during warm inflation, with a signal which changes sign in the squeezed limit (see e.g., BasteroGil et al. 2014).
2.3. Isocurvature nonGaussianity
In most of the models mentioned in this section the focus is on primordial NG in the adiabatic curvature perturbation ζ. However, in inflationary scenarios with multiple scalar fields, isocurvature perturbation modes can be produced as well. If they survive until recombination, these will then contribute not only to the power spectrum, but also to the bispectrum, producing in general both a pure isocurvature bispectrum and mixed bispectra because of the crosscorrelation between isocurvature and adiabatic perturbations (Komatsu 2002; Bartolo et al. 2002; Komatsu et al. 2005; Kawasaki et al. 2008, 2009; Langlois et al. 2008; Hikage et al. 2009, 2013a,b; Langlois & Lepidi 2011; Langlois & van Tent 2011; Kawakami et al. 2012; Langlois & van Tent 2012).
In the context of the ΛCDM cosmology, there are at the time of recombination four possible distinct isocurvature modes (in addition to the adiabatic mode), namely the colddarkmatter (CDM) density, baryondensity, neutrinodensity, and neutrinovelocity isocurvature modes (Bucher et al. 2000). However, the baryon isocurvature mode behaves identically to the CDM isocurvature mode, once rescaled by factors of Ω_{b}/Ω_{c}, so we only consider the other three isocurvature modes in this paper. Moreover, we only investigate isocurvature NG of the local type, since this is the most relevant case in multifield inflation models, which we require in order to produce isocurvature modes. We also limit ourselves to studying just one type of isocurvature mode (considering each of the three types separately) together with the adiabatic mode, to avoid the number of free parameters becoming so large that no meaningful limits can be derived. Finally, for simplicity we assume the same spectral index for the primordial isocurvature power spectrum and the adiabaticisocurvature crosspower spectrum as for the adiabatic power spectrum, again to reduce the number of free parameters. As shown by Langlois & van Tent (2011), under these assumptions we have in general six independent f_{NL} parameters: the usual purely adiabatic one; a purely isocurvature one; and four correlated ones.
The primordial isocurvature bispectrum templates are a generalization of the local shape in Eq. (3):
where I, J, K label the different adiabatic and isocurvature modes. The invariance under the simultaneous exchange of two of these indices and the corresponding momenta means that , which reduces the number of independent parameters from eight to six in the case of two modes (and explains the notation with the comma). The different CMB bispectrum templates derived from these primordial shapes vary most importantly in the different types of radiation transfer functions that they contain. For more details, see in particular Langlois & van Tent (2012).
An important final remark is that, unlike the case of the purely adiabatic mode, polarization improves the constraints on the isocurvature NG significantly, up to a factor of about 6 as predicted by Langlois & van Tent (2011, 2012) and confirmed by the 2015 Planck analysis (PCNG15). The reason for this is that while the isocurvature temperature power spectrum (to which the local bispectrum is proportional) becomes very quickly negligible compared to the adiabatic one as ℓ increases (already around ℓ ≈ 50 for CDM), the isocurvature polarization power spectrum remains comparable to the adiabatic one to much smaller scales (up to ℓ ≈ 200 for CDM). Hence there are many more polarization modes than temperature modes that are relevant for determining these isocurvature f_{NL} parameters. For more details, again see Langlois & van Tent (2012).
2.4. Running nonGaussianity
We briefly describe inflationary models that predict a mildly scaledependent bispectrum, which is also known in the literature as the running of the bispectrum (see e.g., Chen 2005; Liguori et al. 2006; Kumar et al. 2010; Byrnes et al. 2010a,b; Shandera et al. 2011). In inflationary models this running is as natural as the running of the power spectrum, that is the spectral index n_{s}. Other models with strong scale dependence, for example, oscillatory models, are discussed in Sect. 2.5. Further possibilities for strong scale dependence exist (see e.g., Khoury & Piazza 2009; Riotto & Sloth 2011), but we do not consider these in our study. The simplest model (singlefield slowroll, with canonical action and initial conditions) predicts that both the amplitude and scale dependence are of the order of the slowroll parameters (this is true except in some very particular models, see, e.g., Chen et al. 2013), they are small and currently not observable. However, other elaborate but theoretically wellmotivated models make different predictions and these can be used to confirm or exclude such models. Measuring the running of the nonGaussianity parameters with scale is important because this running carries information about, for instance, the number of inflationary fields and their interactions. This information may not be accessible with the power spectrum alone. The first constraints on the running of a local model were obtained with WMAP7 data in Becker & Huterer (2012). Forecasts of what would be feasible with future data were performed by for instance LoVerde et al. (2008), Sefusatti et al. (2009), Becker et al. (2011, 2012), and Giannantonio et al. (2012).
2.4.1. Localtype scaledependent bispectrum
We start by describing models with a localtype mildly scaledependent bispectrum. Assuming that there are multiple scalar fields during inflation with canonical kinetic terms, that their correlators are Gaussian at horizon crossing, and using the slowroll approximation and the δN formalism, Byrnes et al. (2010b) found a quite general expression for the power spectrum of the primordial potential perturbation:
where the indexes a, b run over the different scalar fields.
The nonlinearity parameter then reads (Byrnes et al. 2010b)
where the last line is the general result valid for any number of slowroll fields. The functions f_{cd} (as well as the functions 𝒫_{ab}) can be parameterized as power laws. In the general case, f_{NL} can also be written as
where n_{multi, a} and n_{f, ab} are parameters of the models that are proportional to the slowroll parameters. It is clear that in the general case there are too many parameters to be constrained. Instead we consider two simpler cases, which are among the three models of running nonGaussianity that are analysed in Sect. 5.2.2.
Firstly, when the curvature perturbation originates from only one of the scalar fields (e.g. as in the simplest curvaton scenario) the bispectrum simplifies to (Byrnes et al. 2010b)
In this case
where n_{NG} is the running parameter which is sensitive to the third derivative of the potential. If the field producing the perturbations is not the inflaton field but an isocurvature field subdominant during inflation, then neither the spectral index measurement nor the running of the spectral index are sensitive to the third derivative. Therefore, those selfinteractions can uniquely be probed by the running of f_{NL}.
The second class of models are twofield models where both fields contribute to the generation of the perturbations but the running of the bispectrum is still given by one parameter only (by choosing some other parameters appropriately) as (Byrnes et al. 2010b)
Comparing the two templates (8) and (10) one sees that there are multiple ways to generalize (with one extra parameter) the constant local f_{NL} model, even with the same values for f_{NL} and n_{NG}. If one is able to distinguish observationally between these two shapes then one could find out whether the running originated from single or multiple field effects for example.
Byrnes et al. (2010b) further assumed that n_{fNL} ln (k_{max}/k_{min}) ≪ 1. In our case, ln(k_{max}/k_{min}) ≲ 8 and n_{NG} can be at most of order 0.1. If the observational constraints on n_{NG} using the previous theoretical templates turn out to be weaker, then one cannot use those constraints to limit the fundamental parameters of the models because the templates are being used in a region where they are not applicable. However, from a phenomenological point of view, we wish to argue that the previous templates are still interesting cases of scaledependent bispectra, even in that parameter region. Byrnes et al. (2010b) also computed the running of the trispectra amplitudes τ_{NL} and g_{NL}. For general singlesource models they showed that n_{τNL} = 2n_{NG}, analogous to the wellknown consistency relation , providing a useful consistency check (see, e.g., Smith et al. 2011).
2.4.2. Equilateral type scaledependent bispectrum
General singlefield models that can produce large bispectra having a significant correlation with the equilateral template also predict a mild running nonGaussianity. A typical example is DBIinflation, as studied, for example by Chen (2005) and Chen et al. (2007b), with a generalization within the effective field theory of inflation in Bartolo et al. (2010d). Typically in these models a running NG arises of the form
where n_{NG} is the running parameter and k_{piv} is a pivot scale needed to constrain the amplitude. For example, in the case where the main contribution comes from a small sound speed of the inflaton field, n_{NG} = −2s, where , and therefore running NG allows us to constrain the time dependence of the sound speed^{3}. The equilateral NG with a running of the type given in Eq. (11) is a third type of running NG analysed in Sect. 5.2.2 (together with the local singlesource model of Eq. (8) and the local twofield model of Eq. (10)). We refer the reader to Sect. 3.1.2 for the details on the methodology adopted to analyse these models.
2.5. Oscillatory bispectrum models
Oscillatory power spectrum and bispectrum signals are possible in a variety of wellmotivated inflationary models, including those with an imposed shift symmetry or if there are sharp features in the inflationary potential. Our first Planck temperatureonly nonGaussian analysis (PCNNG13) included a search for the simplest resonance and feature models, while the second Planck temperature with polarization analysis (PCNG15) substantially expanded the frequency range investigated, while also encompassing a much wider class of oscillatory models. These phenomenological bispectrum shapes had free parameters designed to capture the main properties of the key extant oscillatory models, thus surveying for any oscillatory signals present in the data at high significance. Our primary purpose here is to use the revised Planck 2018 data set to determine the robustness of our second analysis, so we only briefly introduce the models studied, referring the reader to our previous work (PCNG15) for more detailed information.
2.5.1. Resonance and axion monodromy
Motivated by the UV completion problem facing largefield inflation, effective shift symmetries can be used to preserve the flat potentials required by inflation, with a prime example being the periodically modulated potential of axion monodromy models. This periodic symmetry can cause resonances during inflation, imprinting logarithmicallyspaced oscillations in the power spectrum, bispectrum and beyond (Chen et al. 2008; Flauger et al. 2010; Hannestad et al. 2010; Flauger & Pajer 2011). For the bispectrum, to a good approximation, these models yield the simple oscillatory shape (see e.g., Chen 2010b)
where the constant ω is an effective frequency associated with the underlying periodicity of the model and ϕ is a phase. The units for the wavenumbers, k_{i}, are arbitrary as any specific choice can be absorbed into the phase which is marginalized over in our results. There are more general resonance models that naturally combine properties of inflation inspired by fundamental theory, notably a varying sound speed c_{s} or an excited initial state. These tend to modulate the oscillatory signal on K = k_{1} + k_{2} + k_{3} constant slices, with either equilateral or flattened shapes, respectively (see e.g., Chen 2010a), which we take to have the form
where . Note that S^{eq} correlates closely with the equilateral shape in (1) and S^{flat} with the orthogonal shape in (2), since the correction from the spectral index n_{s} is small. The resulting generalized resonance shapes for which we search are then
This analysis does not exhaust resonant models associated with a nonBunchDavies initial state, which can have a more sharply flattened shape (“enfolded” models), but it should help identify this tendency if present in the data. In addition, the discussion in Flauger et al. (2017) showed that the resonant frequency can “drift” slowly over time with a correction term to the frequency being proposed, but again we leave that for future analysis. Finally, we note that there are multifield models in which sharp cornerturning can result in residual oscillations with logarithmic spacing, thus mimicking resonance models (Achúcarro et al. 2011; Battefeld et al. 2013). However, these oscillations are more strongly damped and can be searched for by modulating the resonant shape (Eq. (12)) with a suitable envelope, as discussed for feature models below.
2.5.2. Scaledependent oscillatory features
Sharp features in the inflationary potential can generate oscillatory signatures (Chen et al. 2007a), as can rapid variations in the sound speed c_{s} or fast turns in a multifield potential. Narrow features in the potential induce a corresponding signal in the power spectrum, bispectrum, and trispectrum; to a first approximation, the oscillatory bispectrum has a simple sinusoidal behaviour given by (Chen et al. 2007a)
where ω is a frequency determined by the specific feature properties and ϕ is a phase. The wavenumbers k_{i} are in units of Mpc^{−1}. A more accurate analytic bispectrum solution has been found that includes a damping envelope taking the form (Adshead et al. 2012)
where K = k_{1} + k_{2} + k_{3} and the envelope function is given by D(αωK) = αω/(K sinh(αωK)). Here, the modeldependent parameter α determines the large wavenumber cutoff, with α = 0 for no envelope in the limit of an extremely narrow feature. Oscillatory signals generated instead by a rapidly varying sound speed c_{s} take the form
In order to encompass the widest range of physicallymotivated feature models, we modulated the predicted signal (Eq. (15)) with equilateral and flattened shapes, as defined in Eq. (13), they are
Like our survey of resonance models, this allows the feature signal to have arisen in inflationary models with (slowly) varying sound speeds or with excited initial states. In the latter case, it is known that very narrow features can mimic nonBunch Davies bispectra with a flattened or enfolded shape (Chen et al. 2007a).
2.6. NonGaussianity from excited initial states
Inflationary perturbations generated by “excited” initial states (nonBunchDavies vacuum states) generically create nonGaussianities with a distinct enfolded shape (see e.g., Chen et al. 2007b; Holman & Tolley 2008; Meerburg et al. 2009), that is, where the bispectrum signal is dominated by flattened configurations with k_{1} + k_{2} ≈ k_{3} (and cyclic permutations). In the present analysis, we investigate all the nonBunchDavies (NBD) models discussed in the previous Planck nonGaussian papers, where explicit equations can be found for the bispectrum shape functions. In the original analysis (PCNNG13), we described: the vanilla flattened shape model B^{flat} ∝ S^{flat} already given in Eq. (13); a more realistic flattened model B^{NBD} (Chen et al. 2007b) from powerlaw kinflation, with excitations generated at time τ_{c}, yielding an oscillation period (and cutoff) k_{c} ≈ (τ_{c}c_{s})^{−1}; the two leadingorder shapes for excited canonical singlefield inflation labelled B^{NBD1cos} and B^{NBD2cos} (Agullo & Parker 2011); and a nonoscillatory sharply flattened model B^{NBD3}, with large enhancements from a small sound speed c_{s} (Chen 2010b). In the second (temperature plus polarization) analysis (PCNG15), we also studied additional NBD shapes including a sinusoidal version of the original NBD bispectrum B^{NBDsin} (Chen 2010a), and similar extensions for the singlefield excited models (Agullo & Parker 2011), labelled B^{NBD1sin} and B^{NBD2sin}, with the former dominated by oscillatory squeezed configurations.
2.7. Directionaldependent NG
The standard local bispectrum in the squeezed limit (k_{1} ≪ k_{2} ≈ k_{3}) has an amplitude that is the same for all different angles between the largescale mode with wavevector k_{1} and the smallscale modes parameterized by wavevector k_{s} = (k_{2} − k_{3})/2. More generally, we can consider “anisotropic” bispectra, in the sense of an angular dependence on the orientation of the largescale and smallscale modes, where in the squeezed limit the bispectrum depends on all even powers of (where , and all odd powers vanish by symmetry even out of the squeezed limit). Expanding the squeezed bispectrum into Legendre polynomials with even multipoles L, the L > 0 shapes then be used to cleanly isolate new physical effects. In the literature it is more common to expand in the angle , which makes some aspects of the analysis simpler while introducing nonzero odd L moments (since they no longer vanish when defined using the nonsymmetrized smallscale wavevector). We then parameterize variations of local NG using (Shiraishi et al. 2013a):
where P_{L}(μ) is the Legendre polynomial with P_{0} = 1, P_{1} = μ, and . For instance, in the L = 1 case the shape is given by
Bispectra of the directionally dependent class in general peak in the squeezed limit (k_{1} ≪ k_{2} ≈ k_{3}), but they feature a nontrivial dependence on the parameter . The local NG template corresponds to c_{i} = 2f_{NL}δ_{i0}. The nonlinearity parameters are related to the c_{L} coefficients by , , and . The L = 1 and 2 shapes are characterized by sharp variations in the flattened limit, for example, for k_{1} + k_{2} ≈ k_{3}, while in the squeezed limit, L = 1 is suppressed, unlike L = 2, which grows like the local bispectrum shape (i.e. the L = 0 case).
Bispectra of the type in Eq. (20) can arise in different inflationary models, for example, models where anisotropic sources contribute to the curvature perturbation. Bispectra of this type are indeed a general and unavoidable outcome of models that sustain longlived superhorizon gauge vector fields during inflation (Bartolo et al. 2013a). A typical example is the case of the inflaton field φ coupled to the kinetic term F^{2} of a U(1) gauge field A^{μ}, via the interaction term I^{2}(φ)F^{2}, where F_{μν} = ∂_{μ}A_{ν} − ∂_{ν}A_{μ} and the coupling I^{2}(φ)F^{2} can allow for scale invariant vector fluctuations to be generated on superhorizon scales (Barnaby et al. 2012a; Bartolo et al. 2013a)^{4}. Primordial magnetic fields sourcing curvature perturbations can also generate a dependence on both μ and μ^{2} (Shiraishi 2012). The I^{2}(φ)F^{2} models predict c_{2} = c_{0}/2, while models where the primordial curvature perturbations are sourced by largescale magnetic fields produce c_{0}, c_{1}, and c_{2}. The socalled “solid inflation” models (Endlich et al. 2013; see also Bartolo et al. 2013b, 2014; Endlich et al. 2014; Sitwell & Sigurdson 2014) also predict bispectra of the form Eq. (20). In this case c_{2} ≫ c_{0} (Endlich et al. 2013, 2014). Inflationary models that break rotational invariance and parity also generate this kind of NG with the specific prediction c_{0} : c_{1} : c_{2} = 2 : −3 : 1 (Bartolo et al. 2015). Therefore, measurements of the c_{i} coefficients can test for the existence of primordial vector fields during inflation, fundamental symmetries, or nontrivial structure underlying the inflationary model (as in solid inflation).
Recently much attention has been focused on the possibility of testing the presence of higherspin particles via their imprints on higherorder inflationary correlators. Measuring primordial NG can allow us to pin down masses and spins of the particle content present during inflation, making inflation a powerful cosmological collider (Baumann & Green 2012; Chen 2010b; Chen & Wang 2010; Noumi et al. 2013; ArkaniHamed & Maldacena 2015; Baumann et al. 2018; ArkaniHamed et al. 2018). In the case of longlived superhorizon higherspin (effectively massless or partially massless higher spin fields) bispectra like in Eq. (20) are generated, where even coefficients up to c_{n = 2s} are excited, s being the spin of the field (Franciolini et al. 2018). A structure similar to Eq. (20) arises in the case of massive spin particles, where the coefficients c_{i} have a specific nontrivial dependence on the mass and spin of the particles, as well as on wavenumber (ArkaniHamed & Maldacena 2015; Baumann et al. 2018; Moradinezhad Dizgah et al. 2018).
2.8. Parityviolating tensor nonGaussianity motivated by pseudoscalars
Together with the scalar mode, the NG signature in the tensormode sector also provides a powerful probe of inflation. In singlefield inflation with Einstein gravity the tensor NG is highly suppressed. However, the tensor NG is enhanced for several modifications of Einstein gravity (e.g., Maldacena & Pimentel 2011; Gao et al. 2011; Akita & Kobayashi 2016; Domènech et al. 2017; Bartolo & Orlando 2017; Naskar & Pal 2018; Anninos et al. 2019; Ozsoy et al. 2019) or the addition of some extra source fields (in this context see, e.g., Cook & Sorbo 2012, 2013; Barnaby et al. 2012b; Senatore et al. 2014; Mirbabayi et al. 2015; Namba et al. 2016; Agrawal et al. 2018).
In some inflationary scenarios involving the axion field, there are chances to realize the characteristic NG signal in the tensormode sector. In these cases a nonvanishing bispectrum of primordial gravitational waves, , arises via the nonlinear interaction between the axion and the gauge field. Its magnitude varies depending on the shape of the axiongauge coupling, and, in the bestcase scenario, the tensor mode can be comparable in size to or dominate the scalar mode (Cook & Sorbo 2013; Namba et al. 2016; Agrawal et al. 2018).
The induced tensor bispectrum is polarized as (because the source gauge field is maximally chiral), and peaked at around the equilateral limit (because the tensormode production is a subhorizon event). Its size is therefore quantified by the socalled tensor nonlinearity parameter,
with .
In this paper we constrain by measuring the CMB temperature and Emode bispectra computed from (for the exact shape of see PCNG15). By virtue of their parityviolating nature, the induced CMB bispectra have nonvanishing signal for not only the even but also the odd ℓ_{1} + ℓ_{2} + ℓ_{3} triplets (Shiraishi et al. 2013b). Both are investigated in our analysis, yielding more unbiased and accurate results. The Planck 2015 paper (PCNG15) found a best limit of (68% CL), from the foregroundcleaned temperature and highpass filtered Emode data, where the Emode information for ℓ < 40 was entirely discarded in order to avoid foreground contamination. This paper updates those limits with additional data, including largescale Emode information.
We stress that there are also other theoretical models that yield interesting NG signals in the tensortensortensor bispectrum, as well as mixed correlations such as tensortensorscalar and tensorscalarscalar bispectra. Such NG signals would be testable by using the Bmode polarization, while the analysis of Tmode and Emode bispectra provides quite limited information, even with the Planck resolution because of large contamination due to the scalarscalarscalar bispectrum (Meerburg et al. 2016a; Domènech et al. 2017; Bartolo et al. 2019). In contrast, for models involving the axion and the gauge field, the tensor contribution dominates the scalar one at large scales and hence such a scalarmode bias is negligible. This is one reason why we are now focusing on the axionmodel f_{NL} template (Eq. (22)) and measuring it with the PlanckT and Emode data.
3. Estimators and data analysis procedures
3.1. Bispectrum estimators
We give here a short description of the dataanalysis procedures used in this paper. For additional details, we refer the reader to the primordial NG analysis associated with previous Planck releases (PCNG13; PCNG15) and to references provided below.
For a rotationally invariant CMB sky and even parity bispectra (as is the case for combinations of T and E), the angular bispectrum can be written as
where defines the “reduced bispectrum,” and is the Gaunt integral, tahis is the integral over solid angle of the product of three spherical harmonics,
The Gaunt integral (which can be expressed as a product of Wigner 3jsymbols) enforces rotational symmetry. It satisfies both a triangle inequality and a limit given by some maximum experimental resolution ℓ_{max}. This defines a tetrahedral domain of allowed bispectrum triplets, {ℓ_{1}, ℓ_{2}, ℓ_{3}}.
In order to estimate the f_{NL} value for a given primordial shape, we need to compute a theoretical prediction of the corresponding CMB bispectrum ansatz and fit it to the observed 3point function (see e.g., Komatsu & Spergel 2001).
Optimal cubic bispectrum estimators were first discussed in Heavens (1998). It was then shown that, in the limit of small NG, the optimal polarized f_{NL} estimator is described by (Creminelli et al. 2006)
where the normalization N is fixed by requiring unit response to when f_{NL} = 1. C^{−1} is the inverse of the block matrix:
The blocks represent the full TT, TE, and EE covariance matrices, with C^{ET} being the transpose of C^{TE}. CMB a_{ℓm} coefficients, bispectrum templates, and covariance matrices in the previous relation are assumed to include instrumental beam and noise.
As shown in the formula above, these estimators are always characterized by the presence of two distinct contributions. One is cubic in the observed multipoles, and computes the correlation between the observed bispectrum and the theoretical template . This is generally called the “cubic term” of the estimator. The other is instead linear in the observed multipoles. Its role is that of correcting for meanfield contributions to the uncertainties, generated by the breaking of rotational invariance, due to the presence of a mask or to anisotropic/correlated instrumental noise (Creminelli et al. 2006; Yadav et al. 2008).
Performing the inversecovariance filtering operation implied by Eq. (25) is numerically very demanding (Smith et al. 2009; Elsner & Wandelt 2012). An alternative, simplified approach, is that of working in the “diagonal covariance approximation,” yielding (Yadav et al. 2007)
Here, represents the inverse of the following 2 × 2 matrix:
As already described in PCNG13, we find that this simplification, while avoiding the covarianceinversion operation, still leads to uncertainties that are very close to optimal, provided that the multipoles are prefiltered using a simple diffusive inpainting method. As in previous analyses, we stick to this approach here.
A bruteforce implementation of Eq. (27) would require the evaluation of all the possible bispectrum configurations in our data set. This is completely unfeasible, as it would scale as . The three different bispectrum estimation pipelines employed in this analysis are characterized by the different approaches used to address this issue.
Before describing these methods in more detail in the following sections, we would like to stress here, the importance of having these multiple approaches. The obvious advantage is that this redundancy enables a stringent crossvalidation of our results. There is, however, much more than that, as different methods allow a broad range of applications, beyond f_{NL} estimation, such as, for example, modelindependent reconstruction of the bispectrum in different decomposition domains, precise characterization of spurious bispectrum components, monitoring directiondependent NG signals, and so on.
3.1.1. KSW and skewC_{ℓ} estimators
KomatsuSpergelWandelt (KSW) and skewC_{ℓ} estimators (Komatsu et al. 2005; Munshi & Heavens 2010) can be applied to bispectrum templates that can be factorized, they can be written or well approximated as a linear combination of separate products of functions. This is the case for the standard local, equilateral, and orthogonal shapes, which cover a large range of theoretically motivated scenarios. The idea is that factorization leads to a massive reduction in computational time, via reduction of the threedimensional summation over ℓ_{1}, ℓ_{2}, ℓ_{3} into a product of three separate onedimensional sums over each multipole.
The skewC_{ℓ} pipeline differs from KSW essentially in that, before collapsing the estimate into the f_{NL} parameter, it initially determines the so called “bispectrumrelated power spectrum” (in short, “skewC_{ℓ}”) function (see Munshi & Heavens 2010) for details). The slope of this function is shapedependent, which makes the skewC_{ℓ} extension very useful to separate and monitor multiple and spurious NG components in the map.
3.1.2. Running of primordial nonGaussianity
In the previous 2015 analysis, the KSW pipeline was used only to constrain the separable local, equilateral, orthogonal, and lensing templates. In the current analysis we extend its scope by adding the capability to constrain running of nonGaussianity, encoded in the spectral index of the nonlinear amplitude f_{NL}, denoted n_{NG}.
In our analysis we consider both the two local running templates, described by Eqs. (8) and (10) in Sect. 2.4.1, and the general parametrizazion for equilateral running of Sect. 2.4.2, which reads:
where n_{NG} is the running parameter and k_{piv} is a pivot scale needed to constrain the amplitude. Contrary to the two local running shapes this expression is not explicitly separable. To make it suitable for the KSW estimator (e.g. to preserve the factorizability over k_{i}), we can use a Schwinger parametrization and rearrange it as
where k_{sum} = k_{1} + k_{2} + k_{3}.
Alternatively, but not equivalently, factorizability can be preserved by replacing the arithmetic mean of the three wavenumbers with the geometric mean (Sefusatti et al. 2009):
Making one of these substitutions immediately yields the scaledependent version of any bispectrum shape. Analysis in Oppizzi et al. (2018) has shown strong correlation between the two templates, where the former behaves better numerically and is the template of choice for the running in this analysis.
A generalization of the local model, taking into account the scale dependence of f_{NL}, can be found in Byrnes et al. (2010b), as summarized in Sect. 2.4.1.
Unlike f_{NL}, the running parameter n_{NG} cannot be estimated via direct template fitting. The optimal estimation procedure, developed in Becker & Huterer (2012) and extended to all the scaledependent shapes treated here in Oppizzi et al. (2018), is based instead on the reconstruction of the likelihood function, with respect n_{NG}. The method exploits the KSW estimator to obtain estimates of for different values of the running, using explicitly separable bispectrum templates. With these values in hand, the running parameter probability density function (PDF) is computed from its analytical expression.
The computation of the marginalized likelihood depends on the choice of the prior distributions; in Becker & Huterer (2012) and Oppizzi et al. (2018) a flat prior on was assumed. This prior depends on the choice of the arbitrary pivotal scale k_{piv}, since a flat prior on defined at a certain scale, corresponds to a nonflat prior for another scale. The common solution is to select the pivot scale that minimizes the correlation between the parameters. This is in general a good choice, and would work properly in the case of a significant detection of a bispectrum signal. In the absence of a clear detection, however, it is worth noting some caveats. Since the range of scales available is obviously finite, a fit performed at a certain pivot scale tend to favour particular values of n_{NG}. Therefore, there is not a perfectly ‘fair’ scale for the fit. As a consequence, statistical artefacts can affect the estimated constraints in the case of low significance of the measured central value. To prevent this issue, we resort to two additional approaches that make the final n_{NG} PDF pivot independent: the implementation of a parametrization invariant Jeffreys prior; and frequentist likelihood profiling. Assuming that the bispectrum configurations follow a Gaussian distribution, the likelihood can be written as (see Becker & Huterer 2012, for a derivation)
where is the value of the NG amplitude recovered from the KSW estimator for a fixed n_{NG} value of the running, and N is the KSW normalization factor. Integrating this expression with respect to we obtain the marginalized likelihood. Assuming a constant prior we obtain
The Jeffreys prior is defined as the square root of the determinant of the Fisher information matrix ℐ(f_{NL}, n_{NG}). In the case of separable scaledependent bispectra, the Fisher matrix is
where θ_{α} and θ_{β} correspond to or n_{NG} (depending on the value of the index), b_{ℓ1ℓ2ℓ3} is the reduced bispectrum, and the matrix is a Wigner3j symbol.
We search an expression for the posterior distribution marginalized over . Assuming the Jeffreys prior for both parameters and integrating over f_{NL}, we obtain the marginalized posterior
The implementation of this expression in the estimator is straightforward; the only additional step is the numerical computation of the Fisher matrix determinant for each value of n_{NG} considered. The derived expression is independent of the pivot scale.
Alternatively, in the frequentist approach, instead of marginalizing over , the likelihood is sampled along its maximum for every n_{NG} value. For fixed n_{NG}, the maximum likelihood is given exactly by the KSW estimator . From Eq. (32), we see that for this condition the first exponential is set to 1 (since at the maximum), and the profile likelihood reduces to
Notice that this expression also does not depend on the pivot scale. We additionally used this expression to perform a likelihood ratio test between our scaledependent models and the standard local and equilateral shapes.
3.1.3. Modal estimators
Modal estimators (Fergusson et al. 2010, 2012) are based on constructing complete, orthogonal bases of separable bispectrum templates (“bispectrum modes”) and finding their amplitudes by fitting them to the data. This procedure can be made fast, due to the separability of the modes, via a KSW type of approach. The vector of estimated mode amplitudes is referred to as the “mode spectrum”. This mode spectrum is theory independent and it contains all the information that needs to be extracted from the data. It is also possible to obtain theoretical mode spectra, by expanding primordial shapes in the same modal basis used to analyse the data. This allows us to measure f_{NL} for any given primordial bispectrum template, by correlating the theoretical mode vectors, which can be quickly computed for any shape, with the data mode spectrum. This feature makes modal techniques ideal for analyses of a large number of competing models. Also important is that nonseparable bispectra are expanded with arbitrary precision into separable basis modes. Therefore the treatment of nonseparable shapes is always numerically efficient in the modal approach. Finally, the data mode spectrum can be used, in combination with measured mode amplitudes, to build linear combinations of basis templates, which provide a modelindependent reconstruction of the full data bispectrum. This reconstruction is of course smoothed in practice, since we use a finite number of modes. The modal bispectrum presented here follows the same approach as in 2015. In particular we use two modal pipelines, “Modal 1” and “Modal 2”, characterized both by a different approach to the decomposition of polarized bispectra and by a different choice of basis, as detailed in PCNG13, PCNG15, and at the end of Sect. 3.2.2.
3.1.4. Binned bispectrum estimator
The “Binned” bispectrum estimator (Bucher et al. 2010, 2016) is based on the exact optimal f_{NL} estimator, in combination with the observation that many bispectra of interest are relatively smooth functions in ℓ space. This means that data and templates can be binned in ℓ space with minimal loss of information, but with large computational gains. As a consequence, no KSWlike approach is required, and the theoretical templates and observational bispectra are computed and stored completely independently, and only combined at the very last stage in a sum over the bins to obtain f_{NL}. This has several advantages: the method is fast; it is easy to test additional shapes without having to rerun the maps; the bispectrum of a map can be studied on its own in a nonparametric approach (a binned reconstruction of the full data bispectrum is provided, which can additionally be smoothed); and the dependence of f_{NL} on ℓ can be investigated for free, simply by leaving out bins from the final sum. All of these advantages are used to good effect in this paper.
The Binned bispectrum estimator was described in more detail in the papers associated with the 2013 and 2015 Planck releases, and full details can be found in Bucher et al. (2016). The one major change made to the Binned estimator code compared to the 2015 release concerns the computation of the linear correction term, required to make the estimator optimal in the case that rotational invariance is broken, as it is in the Planck analysis because of the mask and anisotropic noise. The version of the code used in 2015, while fast to compute the linear correction for a single map, scaled poorly with the number of maps, as the product of the data map with all the Gaussian maps squared had to be recomputed for each data map. Hence computing real errors, which requires analysing a large set of realistic simulations, was slow. The new code can precompute the average of the Gaussian maps squared, and then quickly apply it to all the data maps. For the full Planck analysis, with errors based on 300 simulations, one gains an order of magnitude in computing time (see Bucher et al. 2016, for more details).
3.1.5. Highfrequency feature and resonantmodel estimator
In addition to the modal expansion described in Sect. 3.1.3, we have considered an extended frequency range of the constant feature model and the constant resonance model with specialized highfrequency estimators, as in PCNG15. The constant feature bispectrum in Eq. (15) is separable and thus allows for the construction of a KSW estimator (Münchmeyer et al. 2014) for direct bispectrum estimation at any given frequency. For the constant resonance model in Eq. (12), we use the method of Münchmeyer et al. (2015), which expands the logarithmic oscillations in terms of separable linear oscillatory functions. We used 800 sine and cosine modes for this expansion.
3.2. Data set and analysis procedures
3.2.1. Data set and simulations
For our temperature and polarization data analyses we use the Planck 2018 CMB maps, as constructed with the four componentseparation methods, SMICA, SEVEM, NILC, and Commander (Planck Collaboration IV 2020). We also make much use of simulated maps, for several different purposes, from computing errors to evaluating the linear meanfield correction terms for our estimators, as well as for performing datavalidation checks. Where not otherwise specified we used the FFP10 simulation data set described in Planck Collaboration II (2020), Planck Collaboration III (2020), and Planck Collaboration IV (2020), which are the most realistic Planck simulations currently available. The maps we consider have been processed through the same four componentseparation pipelines. The same weights used by the different pipelines on actual data have been adopted to combine different simulated frequency channels.
Simulations and data are masked using the common masks of the Planck 2018 release in temperature and polarization; see Planck Collaboration IV (2020) for a description of how these masks have been produced. The sky coverage fractions are, f_{sky} = 0.779 in temperature and f_{sky} = 0.781 in polarization.
3.2.2. Data analysis details
Now we describe the setup adopted for the analysis of Planck 2018 data by the four different f_{NL} estimators described earlier in this section.
In order to smooth mask edges and retain optimality, as explained earlier, we inpaint the mask via a simple diffusive inpainting method (Bucher et al. 2016). First, we fill the masked regions with the average value of the nonmasked part of the map. Then we replace each masked pixel with the average value of its neighbours and iterate this 2000 times. This is exactly the same procedure as adopted in 2013 and 2015.
Linear correction terms and f_{NL} errors are obtained from the FFP10 simulations, processed through the four componentseparation pipelines. To this end, all pipelines use all the available 300 FFP10 noise realizations, except both modal estimators, which use only 160 maps (in order to speed up the computation). The good convergence of the modal pipelines with 160 maps was thoroughly tested in previous releases. There, we showed with a large number of tests on realistic simulations that the level of agreement between all our bispectrum estimators was perfectly consistent with theoretical expectations. Accurate tests have found some mismatch between the noise levels in the data and that of the FFP10 simulations (Planck Collaboration III 2020; Planck Collaboration IV 2020). This is roughly at the 3% level in the noise power spectrum at ℓ ≈ 2000, in temperature, and at the percent level, in polarization. We find (see Sect. 6.2 for details) that this mismatch does not play a significant role in our analysis, and can safely be ignored.
All theoretical quantities (e.g. bispectrum templates and lensing bias) are computed assuming the Planck 2018 bestfit cosmology and making use of the CAMB computer code^{5} (Lewis et al. 2000) to compute radiation transfer functions and theoretical power spectra. The HEALPix computer code^{6} (Górski et al. 2005) is used to perform spherical harmonic transforms.
As far as temperature is concerned, we maintain the same multipole ranges as in the 2013 and 2015 analyses, which is 2 ≤ ℓ ≤ 2000 for the KSW and modal estimators and 2 ≤ ℓ ≤ 2500 for the Binned estimator. The different choice of ℓ_{max} does not produce any significant effect on the results, since the 2000 < ℓ ≤ 2500 range is noise dominated and the measured value of f_{NL} remains very stable in that range, as confirmed by validation tests discussed in Sect. 6. The angular resolution (beam FWHM) of both the cleaned temperature and polarization maps is 5 arcmin.
The main novelty of the current analysis is the use of the lowℓ polarization multipoles that were not exploited in 2015 (ℓ < 40 polarization multipoles were removed by means of a highpass filter). More precisely, KSW and modal estimators work in the polarization multipole range 4 ≤ ℓ ≤ 1500, while the Binned estimator considers 4 ≤ ℓ ≤ 2000. For the same reasons as above, this different choice of ℓ_{max} does not have any impact on the results. The impact of these extra polarization modes on the final results is discussed at the end of Sect. 5.1.
The choice of using ℓ_{min} = 4 for polarization, thus removing the first two polarization multipoles, is instead dictated by the presence of some anomalous results in tests on simulations. When ℓ = 2 and ℓ = 3 are included, we observe some small bias arising in the local f_{NL} measurement extracted from FFP10 maps, together with a spurious increase of the uncertainties. We also notice larger discrepancies between the different estimation pipelines than expected from either theoretical arguments or previous validation tests on simulations. This can be ascribed to the presence of some small level of nonGaussianity in the polarization noise at very low ℓ. We stress that the choice of cutting the first two polarization multipoles does not present any particular issue, since it is performed a priori, before looking at the data (as opposed to the simulations), and generates an essentially negligible loss of information.
In addition, the Binned bispectrum estimator removes from the analysis all bispectrum TEE configurations (i.e. those involving one temperature mode and two polarization modes) with the temperature mode in the bin [2, 3]. This is again motivated by optimality considerations: with these modes included the computed errors are much larger than the optimal Fisher errors, while after removing them the errors are effectively optimal. As errors are computed from simulations, this is again an a priori choice, made before looking at the data. It is not clear why only the Binned bispectrum estimator requires this additional removal, but of course the estimators are all quite different, with different sensitivities, which is exactly one of the strengths of having multiple estimators for our analyses.
The Binned bispectrum estimator uses a binning that is identical to the one in 2015, with 57 bins. The boundary values of the bins are 2, 4, 10, 18, 30, 40, 53, 71, 99, 126, 154, 211, 243, 281, 309, 343, 378, 420, 445, 476, 518, 549, 591, 619, 659, 700, 742, 771, 800, 849, 899, 931, 966, 1001, 1035, 1092, 1150, 1184, 1230, 1257, 1291, 1346, 1400, 1460, 1501, 1520, 1540, 1575, 1610, 1665, 1725, 1795, 1846, 1897, 2001, 2091, 2240, and 2500 (i.e. the first bin is [2, 3], the second [4, 9], etc., while the last one is [2240, 2500]). This binning was determined in 2015 by minimizing the increase in the theoretical variance for the primordial shapes due to the binning.
As in our 2015 analysis, we use two different polarized modal estimators. The Modal 1 pipeline expands separately the TTT, EEE, TTE, and EET bispectra (Shiraishi et al. 2019). It then writes the estimator normalization in separable expanded form and estimates f_{NL} via a direct implementation of Eq. (27). The Modal 2 pipeline uses a different approach (see Fergusson 2014, for details). It first orthogonalizes T and E multipoles to produce new, uncorrelated, and coefficients. It then builds uncorrelated bispectra out of these coefficients, which are constrained independently, simplifying the form and reducing the number of terms in the estimator. However, the rotation procedure does not allow a direct estimation of the EEE bispectrum. Direct EEE reconstruction is generally useful for validation purposes and can be performed with the Modal 1 estimator.
As in 2015, Modal 1 is used to study in detail the local, equilateral, and orthogonal shapes, as well as to perform a large number of validation and robustness tests. Modal 2 is mostly dedicated to a thorough study of nonstandard shapes having a large parameter space (like oscillatory bispectra). The two pipelines are equipped with modal bases optimized for their respective purposes. Modal 1 uses 600 polynomial modes, augmented with radial modes extracted from the KSW expansion of the local, equilateral, and orthogonal templates, in order to speed up convergence for these shapes. The Modal 2 expansion uses a higherresolution basis, including 2000 polynomial modes and a SachsWolfe local template, to improve efficiency in the squeezed limit.
For oscillating nonGaussianities the Modal 2 estimator loses resolution for models with high frequencies. In order to constrain this region of parameter space we use two specialized estimators (Münchmeyer et al. 2014, 2015) that specifically target the highfrequency range of shapes with the low frequency range used for crossvaliation with Modal 2. Both of these estimators are equivalent to those used in PCNG15.
4. Nonprimordial contributions to the CMB bispectrum
In this section we investigate those nonprimordial contributions to the CMB bispectrum that we can detect in the cleaned maps, namely lensing and extragalactic point sources. These then potentially have to be taken into account when determining the constraints on the various primordial NG shapes in Sect. 5. On the other hand, the study of other nonprimordial contaminants (that we do not detect in the cleaned maps) is part of the validation work in Sect. 6.
4.1. NonGaussianity from the lensing bispectrum
CMB lensing generates a significant CMB bispectrum (Goldberg & Spergel 1999; Hanson & Lewis 2009; Mangilli & Verde 2009; Lewis et al. 2011; Mangilli et al. 2013). In temperature, this is due to correlations between the lensing potential and the integrated SachsWolfe (ISW; Sachs & Wolfe 1967) contribution to the CMB anisotropies. In polarization, the dominant contribution instead comes from correlations between the lensing potential and E modes generated by scattering at reionization. Both the ISW and reionization contributions affect large scales, while lensing is a smallscale effect, so that the resulting bispectra peak on squeezed configurations. Therefore, they can significantly contaminate primordial NG measurements, especially for the local shape. It has been known for a while that for a highprecision experiment like Planck the effect is large enough that it must be taken into account. The temperatureonly 2013 Planck results (PCNG13, Planck Collaboration XIX 2014, Planck Collaboration XVII 2014) showed the first detection of the lensing CMB temperature bispectrum and the associated bias. This was later confirmed in the 2015 Planck results (PCNG15, Planck Collaboration XXI 2016) both for Tonly and for the full T+E results. The template of the lensing bispectrum^{7} is given by (Hu 2000; Lewis et al. 2011),
where the X_{i} are either T or E. The tilde on indicates that it is the lensed power spectrum, while and are the temperature/polarizationlensing potential crosspower spectra. The functions are defined by
if the sum ℓ_{1} + ℓ_{2} + ℓ_{3} is even and ℓ_{1}, ℓ_{2}, ℓ_{3} satisfies the triangle inequality, and zero otherwise. Unlike for all other templates, the amplitude parameter is not unknown, but should be exactly equal to 1 in the context of the assumed ΛCDM cosmology.
The results for can be found in Table 1. Error bars have been determined based on FFP10 simulations. As we have seen in the previous releases, the results for Tonly are on the low side. SMICA and NILC have remained stable compared to 2015 and are marginally consistent with the expected value at the 1σ level. SEVEM and Commander, on the other hand, have both decreased compared to 2015 and are now further than 1σ away from unity. However, when polarization is added all results increase and become mostly consistent with unity at the 1σ level. Using the SMICA map (which we often focus on in the rest of the paper, see Sect. 6 for discussion) and the Modal 1 estimator (because it is one of the two estimators for which the lensing template has been implemented in both T and E and the Binned estimator is slightly less wellsuited for this particular shape, since it is a difficult template to bin), we conclude that we have a significant detection of the lensing bispectrum; the hypothesis of having no lensing bispectrum is excluded at 3.5σ using the full temperature and polarization data.
Results for the amplitude of the lensing bispectrum from the SMICA, SEVEM, NILC, and Commander foregroundcleaned CMB maps, for different bispectrum estimators.
It was pointed out in Hill (2018) that the coupling between the ISW effect and the thermal SunyaevZeldovich (tSZ) effect can produce a significant ISWtSZtSZ temperature bispectrum, which peaks for squeezed modes, and can therefore contaminate especially our local and lensing results. The semianalytic approach in Hill (2018) makes predictions for single frequency channels, and cannot be directly applied to the multifrequency componentseparated data we are using. However, it is interesting that such an approach shows an anticorrelation at all frequencies between the ISWtSZtSZ contamination and the lensing bispectrum shape. This could be a possible explanation for the slightly low Tonly value of observed in the final data, across all componentseparated maps. We therefore investigate this issue further, by measuring the amplitude of the lensing bispectrum using SMICA maps in which the tSZ signal has been subtracted in addition to the usual components. The results for this “noSZ” map are given in Table 2. It is clear that the results have increased and are now closer to unity. Using these values, the hypothesis of having no lensing bispectrum is excluded at 3.8σ using the full temperature and polarization data, and the Modal 1 estimator.
Results for the amplitude of the lensing bispectrum from the SMICA noSZ temperature map, combined with the standard SMICA polarization map, for the Binned and Modal 1 bispectrum estimators.
The hypothesis that ISWtSZtSZ residuals are contributing to the temperatureonly lensing bispectrum amplitude result is reinforced by the fact that our measurements from T+E are systematically closer to 1 (polarization does not correlate with ISW and helps in debiasing the result). The SZremoved (hereafter “noSZ”) SMICA measurements of local f_{NL}, however, do not support this hypothesis, since they do not display any large shift, whereas a residual ISWtSZtSZ bispectrum should correlate with all shapes that peak in the squeezed limit. Still, one cannot exclude the possibility that multifrequency componentseparation affects the two cases differently. Both our local noSZ results and a further discussion of this effect are presented in Sect. 6.3.2, where the impact on other primordial shapes is also evaluated and comparisons with simulations are carried out. The final conclusion is that the evidence for ISWtSZtSZ contamination on the temperatureonly lensing bispectrum measurements is not very strong, but the possibility cannot be ruled out.
We also note here that another potentially important source of contamination for the local and lensing shapes is given by the coupling between lensing and the CIB. However, the frequencybyfrequency analysis in Hill (2018) shows that in this case the expected bias is positive at all frequencies. The systematically low values of observed in temperature seem therefore to indicate that the CIB bispectrum contamination does not leak into the final componentseparated maps, at least not at a level that is significant for our analysis. This is further reinforced by the fact that we do not detect any CIB signal directly in the cleaned maps (see Sect. 4.2).
In this paper our main concern with the lensing bispectrum is its influence on the primordial shapes. The bias due to the lensing bispectrum on the estimation of the f_{NL} parameter of another shape S is given by the inner product of the lensing bispectrum (Eq. (37)) with the bispectrum of that shape S, divided by the inner product of the bispectrum S with itself (see PCNG15, for more details, or e.g. Bucher et al. 2016 for a derivation). The values for the bias, as computed from theory, are given in Table 3. Note that the bias values that can be read off from, for example, Table 5 in Sect. 5 can differ slightly from these, because each estimator uses values computed using the approximations appropriate to the estimator. However, those differences are insignificant compared to the uncertainties. As seen already in the previous releases, for Tonly data and for T+E the bias is very significant for local and to a lesser extent for orthogonal NG. For Tonly local NG the bias is even larger than the error bars on f_{NL}. Hence it is quite important to take this bias into account. On the other hand, for Eonly the effect is completely negligible.
Bias in the three primordial f_{NL} parameters due to the lensing signal for the four componentseparation methods.
Lastly, we would like to point out that lensing can also contribute to the covariance (Babich & Zaldarriaga 2004). To lowest order, as was done in this analysis, the spectra in Eq. (27) are replaced with the lensed spectra. However, it was shown by Babich & Zaldarriaga (2004), and later confirmed by Lewis et al. (2011), that the contribution from the connected fourpoint function induced by lensing can become an important contribution to the covariance. Although the analysis here has not shown the effect of lensing to be large (since tests were performed and analysis was done on lensed simulations with no obvious degradation over the Gaussian Fisher errors), it is expected that this will become an important challenge in CMB nonGaussianity constraints beyond Planck. Furthermore, unlike the signal contamination discussed above, this will likely be equally important for polarization. Because of the shape of the lensinginduced covariance, which tends to be largest for squeezed configurations, the local shape will be most affected. The estimates performed by Babich & Zaldarriaga (2004) were done in the flatsky limit and did not include polarization; in addition, computations were truncated at linear order in the potential. Although the latter seems justified, it was later shown that for the lensing power spectrum covariance the linear contribution is actually subdominant to the quadratic contribution (Peloton et al. 2017). For future CMB analysis it will be important to address these open questions. One obvious solution to the extra covariance would be to delens the maps (Green et al. 2017) before applying the estimators. Interestingly, this would automatically remove some of the signalinduced lensing contributions discussed above.
4.2. NonGaussianity from extragalactic point sources
As seen in the previous releases, extragalactic point sources are a contaminant present in the bispectrum as measured by Planck. They are divided into populations of unclustered and clustered sources. The former are radio and latetype infrared galaxies (see e.g. Toffolatti et al. 1998; GonzálezNuevo et al. 2005), while the latter are primarily dusty starforming galaxies constituting the cosmic infrared background (CIB; Lagache et al. 2005). For both types of point sources analytic (heuristic) bispectrum templates have been determined, which can be fitted jointly with the primordial NG templates to deal with the contamination.
The templates used here are the same as those used in PCNG15 (see that paper for more information and references). The reduced angular bispectrum template of the unclustered sources is (Komatsu & Spergel 2001)
This template is valid in polarization as well as temperature. However, we do not detect all the same point sources in polarization as we detected in temperature, so that a full T+E analysis does not make sense for this template. In fact, there is no detection of unclustered point sources in the cleaned Planck polarization map, so that we do not include the Eonly values in the table either. The reduced angular bispectrum template for the clustered sources, for example the CIB, is (Lacasa et al. 2014; Pénin et al. 2014)
where the index is q = 0.85, the break is located at ℓ_{break} = 70, and ℓ_{0} = 320 is the pivot scale for normalization. This template is valid only for temperature; the CIB is negligibly polarized.
The results for both extragalactic point source templates, as determined by the Binned bispectrum estimator applied to the Planck temperature map cleaned with the four componentseparation methods, can be found in Table 4. Because the two templates are highly correlated (93%), the results have been determined through a joint analysis. Contamination from unclustered sources is detected in all componentseparated maps, although at different levels, with SEVEM having the largest contamination. The CIB bispectrum, on the other hand, is not detected in a joint analysis. Both pointsource templates are negligibly correlated with the primordial NG templates and the lensing template (all well below 1% for the unclustered point sources). For this reason, and despite the detection of unclustered point sources in the cleaned maps, it makes no difference for the primordial results in the next sections if point sources are included in a joint analysis or completely neglected.
Joint estimates of the bispectrum amplitudes of unclustered and clustered point sources in the cleaned Planck temperature maps, determined with the Binned bispectrum estimator.
5. Results
5.1. Constraints on local, equilateral, and orthogonal f_{NL}
We now describe our analysis of the standard local, equilateral, and orthogonal shapes. We employ the four bispectrum estimators described in Sects. 3.1.1, 3.1.3, and 3.1.4 on the temperature and polarization maps generated by the SMICA, SEVEM, NILC, and Commander componentseparation pipelines. Further details about our data analysis setup were provided in Sect. 3.2.2. As explained there, the main novelty, compared to the 2015 release, is the use of polarization multipoles in the range 4 ≤ ℓ < 40, which were previously excluded.
Our final results are summarized in Table 5, while data validation tests will be presented in Sect. 6. As in 2015, we show final f_{NL} estimates for Tonly, Eonly, and for the full T+E data set, with and without subtraction of the lensing bias. When we subtract the lensing bias, we assume a theoretical prediction for the lensing bispectrum amplitude based on the Planck bestfit ΛCDM cosmological parameters (see Sect. 4.1). We note that propagating uncertainties in the parameters has a negligible effect on the predicted lensing bispectrum, and the validity of the ΛCDM assumption is of course consistent with all Planck measurements. An alternative to the direct subtraction of the predicted lensing bias of the primordial shapes would be performing a full joint bispectrum analysis (accounting for the measured lensing bispectrum amplitude and propagating the uncertainty into the final primordial NG error budget). While the latter approach is in principle more conservative, we opt for the former both for simplicity and for consistency with previous analyses, as it turns out that the difference between the two methods has no significant impact on our results. If, as an example, we perform a Tonly joint bispectrum analysis using the SMICA map and the Binned estimator, we obtain and . The joint T+E analysis produces instead and . Hence the uncertainties obtained from the joint analysis are practically stable with respect to those shown in Table 5, while the small shift in the local central value (due to the low value of the measured ) would not change our conclusions in any way. More details about the lensing contribution are provided in Sect. 4, where we also discuss the negligible impact of point source contamination on primordial bispectra.
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, SEVEM, NILC, and Commander foregroundcleaned maps.
Table 5 constitutes the most important result of this section. As done in 2013 and 2015, we select the KSW estimator and the SMICA map to provide the final Planck results for the local, equilateral, and orthogonal bispectra; these results are summarized in Table 6. The motivations for this choice are as in the past: SMICA performs well in all validation tests and shows excellent stability across different data releases; the KSW estimator, while not able to deal with nonseparable shapes or reconstruct the full bispectrum, can treat exactly the local, equilateral, and orthogonal templates that are analysed here.
Results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW estimator from the SMICA foregroundcleaned map.
As for the two previous data releases, we note that the agreement between the different estimators – for all the maps and all the shapes considered, both in temperature and polarization, as well as in the full T+E results – is well in line with both our theoretical expectations and Monte Carlo studies, where we found an average measured f_{NL} scatter at the level of ≲σ_{fNL}/3 between different pipelines (this was discussed at length in PCNG13 and PCNG15). The observed scatter between two estimators (for a given realization) can be larger than what the very small reported differences in their variance might suggest. In an ideal, noiseless, fullsky experiment, the observed scatter is only due to differences in the estimator weights, coming from the use of different bispectrum expansions or binning schemes. In such an ideal case, if two different bispectrum templates have a correlation coefficient r, we showed in Appendix B of PCNG13 that the standard deviation of the expected scatter δf_{NL} is given by . This leads to differences between f_{NL} results that are a sizeable fraction of the estimator standard deviation, even for highly correlated weights. To be more explicit, a 95% correlation between different input templates leads, using the formula above, to a σ_{fNL}/3 average scatter in f_{NL} estimates; on the other hand, it would produce only a 5% difference in the final uncertainties.
Considering a more concrete example, let us focus on the current SMICA KSW T+E results, which we quote as our final, recommended bounds. Let us consider, for example, the difference between the KSW and Binned estimator for the local shape. This is relatively large, at approximately 0.3σ_{fNL} after lensingbias subtraction, compared to an error bar difference of 2%. Let us now assume a correlation r ≈ 0.98 between the weights, consistent both with what we see in simulations and with such a small difference in the errors. If we substitute this into our formula, we obtain an expected scatter of σ_{δfNL} = 0.2σ_{fNL}. The observed scatter is then 1.5σ away from this average and therefore fairly consistent with it. The chosen example displays a relatively large difference, compared to all the other combinations that can be built in Table 5; note also that the formula we are using represents an ideal case, and our validation tests on simulations in realistic conditions (e.g. masking and anisotropic noise) show, as expected, that the actual scatter between two pipelines is generally a bit larger than this ideal expectation (PCNG13; PCNG15).
If, instead of comparing central values, we look at uncertainties, we see as expected that all pipelines produce nearly optimal constraints. For the local shape, we see that the Modal 1 pipeline produces 6% smaller errors than for example KSW; however, this is within the expected Monte Carlo error and it seems to be just an effect of the selected simulation sample used to compute σ_{fNL}. One should also consider that Modal 1 uses 160 FFP10 maps to extract the standard deviation, versus 300 maps for the other pipelines. This explanation is confirmed by our many validation tests on different sets of simulations.
A high level of internal consistency is also displayed between f_{NL} estimates obtained from different componentseparated maps, as well as in the comparison of current results with those from previous releases. One small exception is provided by the orthogonal f_{NL} estimate obtained from Commander Tonly data, for which we notice both a larger fluctuation with respect to 2013 and 2015 results and a larger discrepancy with other foreground cleaned maps, in particular with the SMICA one. However, we do not find this worrisome for a number of reasons. First of all, the SMICA – Commander difference is still at the level of the 1σ orthogonal f_{NL} uncertainty and all methods, including Commander, show full consistency with . Therefore, this discrepancy does not pose any problem for the theoretical interpretation of the result. Moreover, this fluctuation completely goes away when accounting for polarization data, the reliability of which has become significantly higher with respect to our previous analysis (see Sect. 6 for details). Finally, the discrepancy is for a very specific shape and it is entirely driven by the alreadynoted fluctuation in the Commander orthogonal result with respect to 2013 and 2015. The other methods remain stable, in particular SMICA, which we take as the map of choice for our final results. The observed fluctuation in orthogonal f_{NL} from Commander can likely be explained by the unavailability of “detset” (i.e. detectorsubset) maps for this release, which constituted in the past a useful input for improving the accuracy of the Commander map. It is, however, important to stress that Commander itself shows excellent agreement with other methods, when measuring f_{NL} for all other shapes and also when correlating the bispectrum modes and bins in a modelindependent fashion, both in temperature and polarization (again see Sect. 6 for a complete discussion of these tests).
Comparing the uncertainties in Table 5 to those in the corresponding table in the 2015 analysis paper (PCNG15), and focussing on the ones for the local shape, since those are most sensitive to lowℓ modes, we see the following: for Tonly data the errors are approximately equal on average, slightly better in some cases, and slightly worse in others. A possible explanation for the slightly larger errors could be the fact that the realism of the simulations (both in temperature and polarization) has improved from FFP8 used in 2015 to FFP10 used here. For our default SMICA componentseparation method this leads to a slightly higher noise power spectrum of a few percent in the regions where noise is relevant. So the errors in 2015 might actually have been slightly underestimated. However, the differences are small enough that they could just be random fluctuations, especially given that not all estimators show the same effect. For Eonly data, on the other hand, we see a clear improvement of the errors for all estimators. That is as expected, since we are now including all the additional polarization modes with 4 ≤ ℓ < 40 in the analysis, and the local shape is quite sensitive to these lowℓ modes. Finally, for the full T+E analysis, we see similar results as for Tonly: all errors have remained the same, to within fluctuations of around a few percent, at most. So one might wonder why the improvement in the Eonly analysis has not translated into a corresponding improvement in the T+E analysis. The answer is relatively simple: the EEEbispectra only have a very small contribution to the final T+E analysis, as shown in Fig. 1. This figure also shows that the equilateral and orthogonal shapes are more sensitive to polarization modes than the local shape, which explains why the errors for the former improve more when going from Tonly to the full T+E analysis than the errors for the latter.
Fig. 1. Weights of each polarization configuration going into the total value of f_{NL} for, from left to right, local, equilateral, and orthogonal shapes. Since we impose ℓ_{1} ≤ ℓ_{2} ≤ ℓ_{3}, there is a difference between, TEE for example (smallest ℓ is temperature) and EET (largest ℓ is temperature). 
In conclusion, our current results show no evidence for nonGaussianity of the local, equilateral, or orthogonal type and are in very good agreement with the previous 2013 and 2015 analyses. We also show in Sect. 6 that the overall robustness and internal consistency of the polarization data set has significantly improved, as far as primordial nonGaussian measurements are concerned.
5.2. Further bispectrum shapes
5.2.1. Isocurvature nonGaussianity
In this section we present a study of the isocurvature NG in the Planck 2018 SMICA map using the Binned bispectrum estimator. This analysis is complementary to the one based on the power spectrum presented in Planck Collaboration X (2020). The underlying modelling approach was discussed in Sect. 2.3, and as explained there, we only investigate isocurvature NG of the local type, and in addition always consider the adiabatic mode together with only one isocurvature mode, there we consider separately CDMdensity, neutrinodensity, and neutrinovelocity isocurvature. In that case there are in general six different f_{NL} parameters: the purely adiabatic one (a,aa) from Sect. 5.1; a purely isocurvature one (i,ii); and four mixed ones.
The results can be found in Table 7, both for an independent analysis of the six parameters (i.e. assuming that only one of the six parameters is present) and for a fully joint analysis (i.e. assuming that all six parameters are present, which is clearly the correct thing to do in the correlated framework described above). The reason for the very differently sized errors is a combination of two effects. The first is a simple normalization issue, due to switching from the more natural ζ and S variables (commonly used in the inflation literature) to Φ_{adi} = −3ζ/5 and Φ_{iso} = −S/5, which are more commonly used in the CMB literature^{8}. The second is that certain parameters depend more on the highℓ adiabatic modes (which are welldetermined), while others are dominated by the much suppressed, and hence unconstrained, highℓ isocurvature modes. For example, for the joint CDM T+E case, when compensating for the normalization factor, one would find the error for a,aa and a,ai to be around 10, and the other four errors to be around 200 (see Langlois & van Tent 2012, for further discussion of these effects).
Results for f_{NL} for local isocurvature NG determined from the SMICA Planck 2018 map with the Binned bispectrum estimator.
As in our analysis of the 2015 Planck data (PCNG15), we see no clear sign of any isocurvature NG. There are a few values that deviate from zero by up to about 2.5σ, but such a small deviation cannot be considered a detection, given the large number of tests and the fact that the deviations are not consistent between Tonly and T+E data. For example, looking at the 300 Gaussian simulations that were used to determine the linear correction and the uncertainties, we find that 84 of them have at least one > 2.5σ result in the two CDM columns of Table 7, while for neutrino density and neutrino velocity the numbers are 62 and 80, respectively.
We see that many constraints are tightened considerably when including polarization, by up to the predicted factor of about 5–6 for the CDM a,ii, i,ai, and i,ii modes in the joint analysis (e.g. −1300 ± 1000 for T only, decreasing to 190 ± 180 for T+E in the CDM i,ai case). Focussing now on the independent results, where it is easier to understand their behaviour (as things are not mixed together), we see that the uncertainties of some of the CDM and neutrinovelocity modes improve by a factor of about 2 when going from the Tonly to the full T+E analysis (e.g. neutrino velocity i,ii changes from 270 ± 230 to −51 ± 110), while the improvements for the neutrinodensity modes are much smaller, of the order of what we see for the pure adiabatic mode. This can be explained if we look at the contribution of the various polarization modes to the total f_{NL} parameters, as indicated in Fig. 2 for the a,ii modes (to give an example). While the diagram for the neutrinodensity mode very much resembles the one for the pure adiabatic local case in Fig. 1, those for CDM and neutrino velocity are quite different and have a much larger contribution from polarization modes (but quite different between the two cases). Hence isocurvature NG also provides a good test of the quality of the polarization maps, and a clear motivation for extending our analyses to include polarization in addition to temperature.
Fig. 2. Weights of each polarization configuration going into the total value of the a,ii mixed f_{NL} parameter for, from left to right, CDM, neutrinodensity, and neutrinovelocity isocurvature, in addition to the adiabatic mode. Since we impose ℓ_{1} ≤ ℓ_{2} ≤ ℓ_{3}, there is a difference between, TEE for example (where the smallest ℓ is temperature) and EET (where the largest ℓ is temperature). 
Comparing the results in Table 7 to those of the corresponding table in our 2015 analysis (PCNG15) we see in the first place that the Tonly results are mostly very stable, generally shifting by much less than 1σ (e.g. joint neutrino velocity i,ii is −970 ± 1400 in 2015 versus −1000 ± 1400 now). Secondly, for the Eonly results we see that all the errors have significantly decreased (e.g. taking again joint neutrino velocity i,ii, we had 2200 ± 1600 in 2015 versus −380 ± 1000 now), in line with the fact that we are now using the additional 4 ≤ ℓ < 40 E modes. We also generally see larger changes in the central values, with several shifts of 1σ or larger; nevertheless, the results are consistent with zero. Because we have added new polarization data to the analysis compared to 2015, and the quality of the polarization data has improved overall, these larger shifts are not unexpected. Finally, looking at the full T+E analysis, we see that the central values shift a bit more than for T only, but remain consistent with zero. For CDM and neutrino density the errors are marginally larger than in 2015, similarly to what we see with the Binned estimator for the local adiabatic shape (see Sect. 5.1). As shown in Fig. 2, both of these depend very little on lowℓ polarization, and so are likely mostly driven by the marginal increase in the Tonly errors. For the neutrinovelocity mode, on the other hand, some of the errors do decrease (e.g. 480 ± 430 in 2015 versus 38 ± 300 now for joint neutrino velocity i,ii). As seen in Fig. 2, this mode depends more strongly on the ETT combination, which does involve lowℓ polarization.
In the results so far we looked at the most general case, having a possible correlation between the isocurvature and adiabatic modes. However, if we assume that the adiabatic and the isocurvature modes have a crosspower spectrum of zero and are completely uncorrelated, then there are only two free f_{NL} parameters, the a,aa and the i,ii ones. In Table 8 we present the results for this uncorrelated case. The independent results are the same as in the previous table and have been repeated for convenience. The significant increase in the Tonly and T+E errors for the neutrino density case when going from the independent to the joint analysis clearly illustrates the fact that its bispectrum template has a large overlap with the adiabatic one, something that also explains the similarity of Figs. 1 and 2 for that case. The CDM and neutrinovelocity modes, on the other hand, have templates that are very different from the adiabatic one and their errors hardly increase (except for neutrino velocity for Eonly data). Again there is no evidence for any isocurvature NG: we do not consider the almost 2σ result for the neutrinodensity isocurvature mode in the T+E joint analysis to be significant, although it will be interesting to keep an eye on this in future CMB experiments with even better polarization measurements.
5.2.2. Running nonGaussianity
In this section we present our analysis of the scaledependent bispectrum shapes described in Sect. 2.4; we obtain these results following the pipeline described in Sect. 3.1.2. In Figs. 3–5 we show the results, respectively, for the onefield local model (Eq. (8)), the twofield local model (Eq. (10)), and the geometricmean equilateral model, where is parametrized as in Eq. (31). Here we present results for both the SMICA and Commander temperature maps. Since the Emode polarization maps are not expected to significantly improve the constraints on n_{NG}, while requiring a significant growth in the computational cost of the estimation pipeline, we do not include polarization in this particular analysis. We show the PDF inferred from all three of the methodologies described. Each point is derived from a KSW estimate of the amplitude for the corresponding scaledependent template with the given running. The KSW pipeline and the map processing steps are the same as applied to the estimation of the local, equilateral, and orthogonal shapes. We include in this analysis the multipole range from 2 to 2000 and the results are corrected for the lensing bias. All curves in Figs. 3–5 are normalized to integrate to one. We consider possible values of the running in the interval n_{NG} = [ − 10, 10]. This interval is two orders of magnitude wider than the theoretical expectation of the models, which are valid in the regime of mild scale dependence, having n_{NG} ≈ 0.1.
Fig. 3. PDF of the running parameter n_{NG} for the onefield local model. Top: SMICA map. Bottom: Commander map. Blue squares give the marginalized posterior assuming a constant prior, green circles are the posterior assuming a Jeffreys prior, and red triangles are the profiled Likelihood. 
Fig. 4. PDF of the running parameter n_{NG} for the twofield local model. Top: SMICA map. Bottom: Commander map. Blue squares give the marginalized posterior assuming a constant prior, green circles are the posterior assuming a Jeffreys prior, and red triangles are the profiled Likelihood. 
Fig. 5. PDF of the running parameter n_{NG} for the geometric mean equilateral parametrization. Top: SMICA map. Bottom: Commander map. Blue squares give the marginalized posterior assuming a constant prior, green circles are the posterior assuming a Jeffreys prior, and red triangles are the profiled Likelihood. 
The effects of the choice of prior are very obvious: if we adopt a constant prior (blue squares) we can always identify a peak in the distribution and define proper constraints; however, this is not the case for the other priors. If we implement an uninformative prior (green circles), the shape of the distribution becomes complex, showing multiple peaks or even diverging on the boundaries, making it impossible to define constraints. A similar behaviour appears in the profiled likelihood approach (red triangles; see Eq. (36)). We used the likelihood to also perform a likelihood ratio test between its maximum value and the value for n_{NG} = 0. Notice that in the case of zero running, these models reduce to the usual local or equilateral shapes. We assume an acceptance threshold of α = 0.01 As expected, we do not find any evidence in favour of scale dependent models.
5.2.3. Resonance and axion monodromy
Now we present results from the Modal 2 and from an adapted KSWtype estimator for the broad class of resonancetype models. These models can be characterized by a bispectrum template having logarithmic oscillations with a scale, as introduced in Sect. 2.5.1. For the Modal 2 analysis we examine templates with frequency in the range 0 < ω ≤ 50 and with a range of possible crosssectional behaviour covering constant, flat, local, and equilateral type templates to span a broader range of possible models. The raw results have been maximised over phase and are presented in Fig. 6 for the four componentseparation methods. Since the errors both grow with frequency, due to increasing suppression by the both transfer functions and projection effects, and at each frequency vary by up to 15% with phase, we only plot the raw significance. Approximate 1 sigma error bars for the constant resonance model can be obtained from the formula σ^{T} ≈ 3.5ω + 51 and σ^{T + E} ≈ 1.7ω + 32 with errors for the equilateral and flattened crosssections approximately 1.7 times larger and those for the local cross section 0.5 times smaller. The results are consistent across componentseparation methods and are comparable with those previously presented in PCNG15, but with slight reductions in significance. Since we are surveying a large range of frequencies, we must correct our results for the “lookelsewhere” effect to asses their true significance. This is done using the optimal methods first proposed in Fergusson et al. (2015), which assess the true significance of both single peaks and also clusters of multiple peaks (which may indicate a model in this class but with a waveform that only partially correlates with the templates used). The results are presented in Table 9, with the largest lookelsewherecorrected result being for the equilateral crosssection with around 2σ for a single peak. This is not on its own significant, but may warrant further investigation with more realistic exact templates near this frequency.
Fig. 6. Generalized resonancemodel significance surveyed over the Modal 2 frequency range with the uppermost panels for the constant resonance model (Eq. 12), showing Tonly (left) and T+E (right) results, then below this the equilateral resonance model (Eq. (13)), followed by the flattened model (Eq. (14)) with the bottom panels, representing the local resonance model with a squeezed envelope. These models have been investigated using Planck temperature data to ℓ_{max} = 2000 and polarization data to ℓ_{max} = 1500, with the Planck componentseparation methods SMICA (blue), SEVEM (orange), NILC (green), and Commander (red). These results are broadly consistent with those found previously (PCNG15), with some broad peaks of moderate significance emerging in both the equilateral and flattened resonance models, somewhat enhanced by including polarization data (upper middle and lower middle right panels). 
Peak statistics, as defined in Fergusson et al. (2015), for the resonance models, showing the maximum “Raw” peak significance, the “Single” peak significance after accounting for the parameter survey lookelsewhere effect, and the “Multi”peak statistic integrating across all peaks (also accounting for the lookelsewhere correction).
5.2.4. Scaledependent oscillatory features
In this section we show the broad range of oscillatory models sourced through features in the inflationary potential or sound speed, described in Sect. 2.5.2, which can be described by a sinusoid multiplied by either a crosssectional template or a scaling function. The results are plotted after maximization over phase in Figs. 7 and over envelope widths in 8. Since the errors both grow with frequency, due to increasing suppression by the both transfer functions and projection effects, and at each frequency vary by up to 20% with phase, we only plot the raw significance. Approximate 1 sigma error bars for the feature model can be obtained from the formula σ^{T} ≈ 0.41ω + 60 and σ^{T + E} ≈ 0.22ω + 33 with errors for the equilateral and flattened crosssections approximately 1.7 times larger. However, when we vary the envelope the error bars change by several orders of magnitude so it is not possible to provide approximate formula in this case. Lookelsewhereadjusted results are presented in Table 10. Here all results are consistent with Gaussianity for all models in this class.
Fig. 7. Generalized feature model significance surveyed over the Modal 2 frequency range 0 < ω < 350, after marginalizing over the phase ϕ. Top panels: results for the constant feature model (Eq. (15)) with T only (left) and T+E (right), middle panels: for the equilateral feature model (Eq. (18)), and bottom panels: for the flattened feature model (Eq. (19)). The same conventions apply as in Fig. 6. These feature model results generally have lower significance than obtained previously (PCNG15), with polarization data not tending to reinforce apparent peaks found using temperature data only. 
Fig. 8. Top: significance of singlefield models for the potential feature case with a K^{2} cos ωK scaling dependence (Eq. (16)). Bottom: significance for rapidly varying sound speed with a K sin ωK scaling (Eq. (17)). Left panels: for T only and right panels: for T + E. These results have been marginalized over the envelope parameter α (determined by feature width and height) from α = 0 to αω = 90. There appears to be no evidence for these very specific signatures in this frequency range when polarization data are included. 
Peak statistics, as defined in Fergusson et al. (2015), for the different feature models, showing the Raw peak maximum significance (for the given Modal 2 survey domain), the corrected significance of this Single maximum peak after accounting for the parameter survey size (the lookelsewhere effect) and the Multipeak statistic which integrates across the adjusted significance of all peaks to determine consistency with Gaussianity.
5.2.5. Highfrequency feature and resonantmodel estimator
We used the highfrequency linear oscillation estimator described in Sect. 3.1.5 to probe the constant feature model with frequencies up to ω = 3000. In the overlapping frequency range we also confirmed that the results are compatible with those of the Modal pipeline. The results are shown in the upper two panels of Fig. 9. Due to the computational demands of this estimator we have only analysed SMICA maps, since it was already shown in PCNG15 that all componentseparations methods agree well for this estimator. No statistically significant peak is found, and the distribution of peaks is consistent with the 2015 analysis. The highest peak is 2.9σ in T only and 3.5σ in T + E, compared to the Gaussian expectation of 3.1(±0.3)σ^{9}.
Fig. 9. Highfrequency estimator results for feature and resonance models. Upper two panels: significance for the constant feature model (Eq. (15)) surveyed over the frequency range 0 < ω < 3000 after marginalizing over the phase ϕ, for T and T + E SMICA maps. Bottom two panels: significance for the constant resonance model (Eq. (12)) surveyed over the frequency range 0 < ω < 1000 after marginalizing over the phase ϕ, for T and T + E SMICA maps. As for the Modal expansion case, the highfrequency results generally have lower significance than obtained previously (PCNG15), with polarization data not tending to reinforce apparent peaks found using temperature data only. 
For the constant resonance model, we used the high frequency estimator described in Sect. 3.1.5 to cover the frequency range 0 < ω < 1000. The results are shown in Fig. 7 (lower two panels) for SMICA data. The highest peak is 3.1σ for T only and 3.0σ for T + E, compared to the Gaussian expectation of 3.4σ ± 0.4σ. The Gaussian expectation was obtained from the covariance matrix of the estimators using the method of Meerburg et al. (2016b). In summary we do not find evidence for nonGaussianity in the highfrequency feature and resonancemodel analysis.
We used the results obtained here to perform a joint power and bispectrum analysis, which was presented in the accompanying paper Planck Collaboration X (2020) in a search for correlated features. No evidence was found for such correlated features in either the power spectrum or the bispectrum.
5.2.6. Equilateraltype models and the effective field theory of inflation
Many physically wellmotivated inflationary models produce nonGaussianity of the equilateral type. While the equilateral template accurately approximates models in this class, it is nevertheless interesting to constrain the exact templates for two reasons. Firstly, while the equilateral template correlates well with the true models, they are only normalized in the equilateral limit, which neglects the rest of configuration space and so is only approximate. This is demonstrated by the fact that the uncertainties for the true templates differ by up to a factor of 2 compared with the equilateral limit. Secondly, even small deviations in correlation can result in much larger shifts in measurements, for example models that are 95% correlated should produce measurements that are on average 1/3σ apart. Here we constrain the exact templates for all inflationary bispectra that fit into the broad equilateral class; for details on these exact templates we refer the reader to the paper PCNG15. The results for all models in this class are presented in Table 11 and are all below 1σ, as expected due to their correlation with the equilateral template, but with significant spread due to the small differences in correlation between templates.
Constraints on models with equilateraltype NG covering the shapes predicted by the effective field theory of inflation, together with constraints on specific noncanonical inflation models, such as DBI inflation.
5.2.7. Models with excited initial states (nonBunchDavies vacua)
Inflationary models that modify the initial vacuum for the inflaton generally produce shapes in the flattened class, meaning that they peak in triangle configurations with zero area. There is a wide variety of models of this type proposed in the literature, which are described in detail in Sect. 2.6. For simplicity we also include here the results for the warm inflation template, described in Sect. 2.2, since it exhibits similar behaviour once projected, despite the mechanism behind it being quite different. The results for this class are presented in Table 12. The most significant measurement is for the NBD sin template, which produces results around 2σ. However, it is important to note that this result involves a marginalization over a frequency type parameter so there is a lookelsewhere effect that has not yet been taken into account. Because of this we expect the true significance will be lower, so we can safely claim that our results are consistent with Gaussianity, while noting that it may be interesting to revisit oscillatory NBD templates, like the NBD sin model, using future data sets.
Constraints on models with excited initial states (nonBunchDavies models), as well as warm inflation.
5.2.8. Directiondependent primordial nonGaussianity
Here we present results for inflationary models where gauge fields induce a direction dependance to a local bispectrum, as described in Sect. 2.7. Results for the L = 1 and L = 2 templates are presented in Table 13 from the Modal 2 pipeline. Due to the complicated behaviour in the squeezed limit, the convergence for these models is lower than for some others, but the Modal correlation remains above 93% so a perfect estimator may see shifts of order 0.5σ. The largest result is 1.9σ for L = 1 and SMICA, but this is not seen consistently across all componentseparation methods so cannot be considered robust. We can then conclude that our results are consistent with Gaussianity.
Directiondependent NG results for both the L = 1 and L = 2 models from the Modal 2 pipeline.
5.2.9. Parityviolating tensor nonGaussianity motivated by pseudoscalars
In this section, we report constraints on the tensor nonlinearity parameter (Eq. (22)) obtained from temperature and Emode polarization maps. As in our 2015 analysis, we examine even and odd ℓ_{1} + ℓ_{2} + ℓ_{3} multipole domains, employing the original parityeven Modal estimator, (Fergusson et al. 2010, 2012; Shiraishi et al. 2019), and its parityodd version, (Shiraishi et al. 2014, 2015, 2019), respectively. The constraints obtained from both domains are also combined by computing
where F^{even/odd} is the Fisher matrix from ℓ_{1} + ℓ_{2} + ℓ_{3} = even/odd. This analysis is performed for under the multipole ranges 2 ≤ ℓ ≤ 500 in temperature and 4 ≤ ℓ ≤ 500 in polarization. Here, the use of the first 40 multipoles of the Emode polarization data, which were disregarded in the 2015 analysis, boosts contributions of TTE, TEE, and EEE to constraining . Although the data, simulations (used for the computation of the linear terms and error bars) and analysis details (e.g. masks, beams, and noise distributions) are somewhat different, other settings for the estimation are basically the same as in the 2015 analysis (PCNG15).
The results from four different componentseparated maps are summarized in Table 14. We confirm there that the sizes of errors in the parityodd T+E analysis reduce to almost the same level as the parityeven counterparts. This is because the lowℓ signal of TTE, which dominates the signaltonoise ratio from ℓ_{1} + ℓ_{2} + ℓ_{3} = odd (Shiraishi et al. 2013b), is now taken into account. Although the errors for the parityeven case are as large as the 2015 ones, owing to the sensitivity improvement of the parityodd part, the wholedomain constraints become more stringent.
Results for the tensor nonlinearity parameter obtained from the SMICA, SEVEM, NILC, and Commander temperature and polarization maps.
Regardless of some updates, we find no > 2σ signal, which is consistent with the conclusion of the 2015 analysis. This indicates no parity violation in the primordial Universe and accordingly gives constraints on some axion inflationary models (see Sect. 8.3).
5.3. Bispectrum reconstruction
5.3.1. Modal bispectrum reconstruction
The Modal bispectrum estimator filters the Planck foregroundremoved CMB maps, for example, SMICA, SEVEM, NILC, and Commander, using n_{max} = 2001 polynomial modes to obtain model coefficients β_{n}. This procedure is undertaken to obtain all auto and crosscorrelations between the temperature and polarization components TTT, TTE, TEE, and EEE. We can then use the β_{n} coefficients with the polynomial modes to obtain a full 3D reconstruction of the Planck temperature and polarization bispectra and this is shown in Fig. 10. These bispectra are in close agreement with those published in the paper analysing the Planck 2015 results (PCNG15) when comparison is made in the signaldominated regime.
Fig. 10. Signaltonoiseweighted temperature and polarization bispectra obtained from Planck SMICA maps using the Modal reconstruction with n_{max} = 2001 polynomial modes. The bispecra are (clockwise from top left) the temperature TTT for ℓ ≤ 1500, the Emode polarization EEE for ℓ ≤ 1100, the mixed temperature and polarization TEE (with T multipoles in the zdirection), and lastly TTE (with E multipoles in the zdirection). The S/N thresholds are the same between all plots. 
5.3.2. Binned bispectrum reconstruction
The Binned bispectrum estimator can also be used to study the bispectrum itself, in addition to determining the f_{NL} amplitudes of specific templates, as in the previous sections. In particular we can investigate if any nonGaussianity beyond that of the explicit models tested can be found in the cleaned CMB maps. Since we want to detect specific bintriplets standing out from the noise, we work with the lineartermcorrected signaltonoiseratio bispectrum. Except for the highℓ bins, where a pointsource signal is present, the bintriplets are noise dominated. Instead of rebinning, we can use a Gaussian kernel to smooth the bispectrum, so that structure localized in harmonic space stands out from the noise. In this process, we mask out a few bintriplets that have nonGaussian noise (due to the fact that they contain very few valid ℓtriplets). We also have to take into account edge effects from the nontrivial domain of definition of the bispectrum. The method has been described in PCNG15 and more extensively in Bucher et al. (2016).
Slices of the smoothed binned signaltonoise bispectrum ℬ_{i1i2i3}, with a Gaussian smoothing of σ_{bin} = 2, are shown in Figs. 11 and 12. These are slices for the 20th and 40th ℓ_{3}bin as a function of ℓ_{1} and ℓ_{2}. For the crossbispectra mixing T and E modes, we defined ; and , with corresponding variances Var(B^{T2E}) = Var(TTE) + Var(TET) + Var(ETT) + 2 Cov(TTE, TET) + 2 Cov(TTE, ETT) + 2 Cov(TET, ETT), and similarly for Var(B^{TE2}), where we have omitted the bin indices for clarity. The red and blue regions correspond to a significant NG, whereas grey areas in Figs. 11 and 12 show regions where the bispectrum is not defined. Results are shown for the four componentseparation methods SMICA, SEVEM, NILC, and Commander, and for TTT, T2E, TE2, and EEE. We also show the TTT bispectra after we remove the best jointfit unclustered and clustered pointsource contribution (see Table 4). In the top row of Fig. 12, we can clearly see this significant pointsource signal at high ℓ. The T2E and TE2 bispectra do not have any obvious signals standing out, but we see some stronger NG in the EEE combination. Removing the best jointfit contribution from all shapes (Tables 4 and 5) does not reduce this region (a case we do not show). The main difference with the previous release is the better qualitative agreement between the four componentseparation methods in polarization, as was already the case for temperature. Commander and SEVEM now show similar structures as the NILC and SMICA bispectra, which have remained quite stable.
Fig. 11. Smoothed binned signaltonoise bispectra ℬ for the Planck 2018 cleaned sky maps. We show slices as a function of ℓ_{1} and ℓ_{2} for a fixed ℓ_{3}bin [518, 548]. From left to right results are shown for the four componentseparation methods SMICA, SEVEM, NILC, and Commander. From top to bottom we show: TTT; TTT cleaned of clustered and unclustered point sources; T2E; TE2; and EEE. The colour range shows signaltonoise from −4 to +4. The light grey regions are where the bispectrum is not defined, either because it lies outside the triangle inequality or because of the cut . 
In order to quantify the possible residual nonGaussianity in these smoothed bispectra, we focus on the minimum and maximum of our bintriplets. In the case of statistically independent Gaussian numbers, we can calculate the probability distribution of the extreme value statistics. However, once we introduce correlations due to the smoothing with nontrivial boundary conditions, we do not have an analytical formula. We can instead generate Monte Carlo simulations of Gaussian random numbers with the same boundary conditions, apply the smoothing as for the data signaltonoise bispectrum, and compute the pvalues of the observed extremum statistics as the fraction of simulations having a more extreme extremum than our data. This requires many simulations to study the very unlikely events. In the current analysis, 10^{6} Monte Carlo simulations turn out to be sufficient and so we do not need to use the semianalytical Ansatz introduced in Bucher et al. (2016).
In Table 15, we report for the smoothing lengths σ_{bin} = 1, 2, and 3, the twotailed pvalue^{10} of the maximum and of the minimum, defined as p = 2 min[Prob(X_{MC} ≤ X_{data}),Prob(X_{MC} ≥ X_{data})], where X is either the minimum or the maximum of the distribution. As expected, we detect a highly significant departure from Gaussian statistics in the maxima for TTT when we do not correct for the contribution from point sources, and the signal stands out more when increasing the smoothing kernel size. When looking at the SEVEM data for σ_{bin} = 2 and 3, we find no simulation with a higher maximum, but for all the other cases our analysis should be robust. Most bispectra seem to be compatible with a simple Gaussian distribution, except for the EEE bispectrum in the region shown in Fig. 12, for multipole triplets around [900, 1300, 1800]. We also see some significantly high maxima for the TTT bispectrum of SEVEM and Commander (even after correcting for point sources), located around [800, 1100, 2000] (shown in Fig. 13). The origin of these signals clustered in multipole space is not understood.
Fig. 13. Similar to Fig. 11, but for the ℓ_{3}bin [771,799], and only for TTT cleaned of clustered and unclustered point sources. 
Twotailed pvalues of the maxima and the minima of the smoothed bispectra.
6. Validation of Planck results
Two important potential sources of systematic effects in f_{NL} estimation are foreground residuals (which can also be related to the choice of mask) and an imperfect modelling of instrumental noise, either of which can lead to miscalibration of the linear correction term in the estimator. In this section, we perform a battery of tests aimed at testing the impact of these systematics. We typically chose only one of the KSW, Binned, or Modal pipelines for each of the tests described below. This is possible because of their excellent and wellverified agreement on both data and simulations.
6.1. Dependence on foregroundcleaning method
6.1.1. Comparison between f_{NL} measurements
For this test, we consider 160 FFP10based simulations with realistic beams and noise levels. A starting set of singlefrequency Gaussian maps is processed through the different componentseparation pipelines. Each pipeline combines the various frequency maps by adopting exactly the same coadding and filtering approach as done for the data. At the end, we have four sets of 160 frequency coadded realizations (the SMICA set, the SEVEM set, the NILC set and the Commander set). For each realization, in any of the four sets, we measure f_{NL} for the local, equilateral, and orthogonal shapes, using the Modal 1 pipeline. Then, for each shape and pair of methods, we build the difference between the two results and call it Δf_{NL}. To give an example, let us consider SMICA, SEVEM, and the local shape. For each of the 160 realizations, we measure from the SMICA map (), repeat the operation with the SEVEM map (), and build the difference . By repeating this for each possible pair of methods, we obtain six sets of 160 differences Δf_{NL} for a given shape.
We then extract the standard deviation of each set, which we call σ_{ΔfNL} (not to be confused with σ_{fNL}; σ_{ΔfNL} is the standard deviation of the scatter between two cleaning methods, whereas σ_{fNL} is simply the f_{NL} uncertainty for a given method). The simulations we use do not include any foreground component; therefore, the extracted σ_{ΔfNL} is only due to different frequency weighting and filtering schemes, and to further operations during NG estimation, such as inpainting. We can use these quantities to verify the consistency between the f_{NL} differences obtained from data and from simulations. A large Δf_{NL} scatter observed in the data, for a given pair of cleaned maps, compared to the corresponding expectation σ_{ΔfNL}, would signal potential residual foreground contamination in at least one of the two maps.
Our results for this test are summarized in Table 16. For each pair of componentseparation methods, and the three standard shapes, we show the measured scatter Δf_{NL} and the ratio Δf_{NL}/σ_{Δ}. In order to assess the statistical significance of the result, we need of course to take into account that a multiplicity correction is necessary, since we are considering six pairs of methods. In Table 17 we report the fraction of starting simulations for which, after the componentseparation processing, Δf_{NL}/σ_{Δ} is larger than 1, 2 or 3, for any pair of methods. We also report the largest value of Δf_{NL}/σ_{Δ}, measured for all combinations and simulations. We see from this table that measured values of Δf_{NL}/σ_{Δ} from the data (Table 16) are not particularly unusual up to Δf_{NL}/σ_{Δ} ≈ 3.
Comparison between local, equilateral, and orthogonal f_{NL} results, obtained using the four different componentseparation pipelines and the Modal 1 bispectrum estimator.
Fraction of simulations for which the differences between pairs of componentseparation methods are above various levels of σ.
This leads us to a first interesting observation, namely that polarizationonly f_{NL} results show no significant discrepancy among different cleaned maps. This is a large improvement with respect to our 2015 analysis, where in a similar test we found large differences in EEE bispectra for several combinations. Such discrepancies, together with other anomalies, led us in the previous release to warn the reader that all polarizationbased f_{NL} measurements were less robust than TTT estimates, and had to be considered as preliminary. This is no longer the case, since both this and other tests (see next section) show that the polarization data are now fully reliable for primordial NG studies. Achieving such reliability was indeed one of the main goals for our analysis in this data release.
Despite consistency in polarization, somewhat surprisingly we now find some relevant discrepancies in Tonly results, even though limited to specific methods and the orthogonal shape only. The most striking example is a large difference between the orthogonal Tonly f_{NL} measurements obtained with SMICA and Commander (Δf_{NL}/σ_{Δ} ≈ 10). Smaller but nonnegligible are also the orthogonal Tonly discrepancies for the NILC – Commander and SMICA – SEVEM pairs (both with Δf_{NL}/σ_{Δ} ≈ 4). These results have been crosschecked using the Binned pipeline, finding agreement between estimators: the Binned estimator finds Δf_{NL}/σ_{Δ} of approximately 8, 4, and 4, for these three cases, respectively.
As anticipated in Sect. 5.1, a closer inspection shows that these differences are less worrisome that they might appear at first glance. First of all, a comparison with our 2013 and 2015 results shows that SMICA and NILC orthogonal TTT measurements have remained very stable, while SEVEM and especially Commander display significant changes in this data release. Considering that all pipelines agreed very well and displayed robustness to a large number of validation tests in temperature in both previous data releases, we conclude that the latter two methods can be identified as the sources of the current orthogonal Tonly discrepancy. This is good news, since the main componentseparation method that we have focused on for f_{NL} analysis (including in PCNG13 and PCNG15) is SMICA.
It is also important to stress that all significant discrepancies are specifically confined to the orthogonal TTT case. The other shapes and also a modebymode or binbybin correlation analysis over the full bispectrum domain (see next section) show no other signs of anomalies for any componentseparation method. Even considering the largest discrepancy, arising from the SMICA – Commander pair, this still amounts to just a 1σ deviation in . This is due to the fact that the level of agreement displayed by different cleaning methods on simulations typically amounts to a small fraction of the f_{NL} errors. The implications for inflationary constraints of shifting by 1σ are essentially negligible. As mentioned earlier, the changes in orthogonal f_{NL} coming from Commander can be explained through the unavailability of detectorset maps for the current release.
The main conclusion to be drawn from the validation tests described in this section is that the reliability of our final T+Ef_{NL} results has significantly increased with respect to the previous release, thanks to a clear improvement in the robustness of the polarization data. The issues we find in the temperature data are confined to the orthogonal shape and to specific componentseparation pipelines, not affecting the final SMICA measurements used for inflationary constraints.
Of course, any considerations in this section that might lead us to a preference for specific componentseparation methods, apply only to primordial NG analysis, and not to other cosmological or astrophysical analyses. There is no generally preferred cleaning method for all applications, and separate assessments should be conducted case by case.
6.1.2. Comparison between reconstructed bispectra
The local, equilateral, and orthogonal directions already cover a significant part of the entire bispectrum domain, and deserve special attention, since they are crucial for inflationary constraints. It is nevertheless useful to also check the agreement between componentseparation pipelines in a more general, modelindependent fashion. For this purpose we calculate the coefficient of determination (see Allen 1997), denoted by R^{2}, from the mode or bin amplitudes extracted from different foregroundcleaned maps. The coefficient of determination is a standard statistical measure of what proportion of the variance of one variable can be predicted from another. For example in our case a score of 0.9 for SMICA – SEVEM would mean that 90% of the mode or bin amplitudes measured in SMICA could be explained by the mode or bin amplitudes measured in SEVEM and 10% of the amplitudes would be unexplained.
Scatter plots of modes are shown in Fig. 14, while those for bins are shown in Fig. 15. For TTT the lowest Modal coefficient of determination is 0.91 (for SEVEM – NILC), while for the Binned bispectrum values all TTT coefficients of determination are larger than 0.99. For EEE the lowest Modal value is 0.71, again for SEVEM – NILC, while the corresponding Binned value is 0.78. The lowest Binned value overall is 0.71, for SMICA – SEVEM ETT (not shown in the figure). It is also possible to see the impact of excluding the lowest [2, 3] bin from the analysis for E; if it were included that value would drop from 0.71 to 0.59.
Fig. 14. Scatter plots of the 2001 Modal 2 coefficients for each combination of componentseparation methods, labelled with the R^{2} coefficient of determination. The figures on the left are for the temperature modes and those on the right for pure polarization modes (the Modal 2 pipeline reconstructs the component of the EEE bispectrum that is orthogonal to TTT, so this is not exactly the same as the EEE bispectrum of the other estimators). 
Fig. 15. Scatter plots of the bispectrum values in each bintriplet (13 020 for TTT, all of which have been multiplied by 10^{16}, and 11 221 for EEE, all of which have been multiplied by 10^{20}) for all combinations of componentseparation methods, with the R^{2} coefficient of determination. The figures on the left are for TTT and those on the right for EEE. 
These results show a very good level of agreement between all methods. Again, we see a large improvement in the polarization maps, by comparing with a similar test performed in PCNG15. Also interesting is the fact that no anomalies show up in temperature data for any specific mode or bin. The issues discussed in the previous section, which affected only specific componentseparation pipelines, seem also completely confined to the combination of modes and bins that selects the orthogonal shape region of the bispectrum domain, in such a way as to be invisible in this more general analysis.
6.2. Testing noise mismatch
Accurate tests of the FFP10 maps have found some level of mismatch between the noise model adopted in the simulations and the actual noise levels in the data (Planck Collaboration III 2020; Planck Collaboration IV 2020). In temperature this is roughly at the 3% level in the noise power spectrum at ℓ ≈ 2000. Percent level differences are also seen in polarization. This issue raises some concern for our estimators, since we use simulations to calibrate them. We thus decided to perform some simulationbased tests to check to what extent this noise mismatch may affect our Monte Carlo errors (as long as the noise is Gaussian, noise mismatches of this kind cannot bias the estimators). For all the work in this section we consider SMICA maps only.
Uncertainties in nonGaussianity parameters might be affected by two effects. One is the suboptimality of the estimator weights in the cubic term, if the power spectra extracted from simulations do not match the data. The other is an imperfect Monte Carlo calibration of the linear correction term, leading to the inability to fully correct for anisotropic and correlated noise features. Given that we are considering percentlevel corrections, we do not expect the former of the two effects to be of particular significance; this is also confirmed by past analyses, in which the cosmological parameters have been updated several times, leading to percent changes in the fiducial power spectrum, without any appreciable difference in the f_{NL} results. The effect on the linear term is harder to assess, and to some extent it depends on the specific noise correlation properties in the data, which are possibly not fully captured in the simulations. In general, it is reasonable to expect a power spectrum mismatch of a few percent to have only a small impact, unless significant spatial correlations between large and small scales are present in the data and not captured by the simulation noise model. The linear term correction is generally dominated by the mask, in the mostly signaldominated regime we consider for our analysis.
As a first test, we generate “extranoise” multipoles, drawing them from a zeromean Gaussian distribution having as power spectrum the difference between the simulation and data noise power spectra. We then add these realizations to the original noise simulations and apply our estimators to this extranoise mock data set. However, we calibrate our linear term, estimator normalization, and weights using the original noise maps. We check the effect of the miscalibration through a comparison between the f_{NL} measurements obtained in this way and those extracted from the original realizations, without extra noise included. As a figure of merit, we consider the local f_{NL} results, since they are generally most sensitive to noise features. At the end of this analysis, we verify that the impact of the noise mismatch is very small: the uncertainties change by a negligible amount with respect to the exactly calibrated case, while the scatter between f_{NL} measured with and without extra noise in the maps is zero on average, as expected, and has a standard deviation σ_{ΔfNL} ≈ σ_{fNL}/10.
Our results are shown in detail in the first row of Table 18 and in Fig. 16. The former reports the f_{NL} error bars and the standard deviation of mapbymap f_{NL} differences, while the latter shows the mapbymap f_{NL} scatter for 50 realizations. These very small deviations were largely expected, given the small change in the amplitude of the noise power spectrum we are considering.
Fig. 16. Effects on f_{NL} scatter, simulationbysimulation, of different levels of noise mismatch, as described in Sect. 6.2. Top: EEEonly results. Bottom: combined T+E results. The magenta line represents the case in which we generate an extranoise component in harmonic space, both in temperature and polarization, with variance equal to the difference between the noise spectra in the data and the simulations. All other lines represent cases in which we generate extranoise maps in pixel space, with the rms per pixel extracted from simulations and rescaled by different factors, as specified in the legend; the label “TQU” further specifies the pixelspace test in which extra noise is added to both temperature and polarization data, unlike in the other three pixelspace cases, in which we include only a polarization extranoise component. 
In the test we have just described we generate uncorrelated extranoise multipoles (with a nonwhite spectrum). On the other hand, we know that the estimator is most sensitive to couplings between large and small scales, which require a properly calibrated linear term correction. Therefore, we decided to also test the impact of a possible linear term miscalibration of this type, by generating and studying an extranoise component directly in pixel space. In this case, ℓspace correlations between large and small scales arise due the spatially anisotropic distribution of the noise. We proceed as follows. After extracting the noise rms of the SMICA FFP10 polarization simulations, we rescale it by a fixed factor A, and we use this rescaled rms map to generate new Gaussian “extranoise" realizations in pixel space, which we add to the original noise maps. We consider different cases. Firstly, we take a rescaling factor A_{TQU} = 0.2 for both temperature and polarization noise maps. We then perform a more detailed study of the effect on polarization maps only. For this purpose, we leave the temperature noise unchanged and rescale the polarization rms noise by factors ranging from A = 0.1 up to A = 0.3.
We see again from the summary of results reported in Table 18 that the f_{NL} error change is always very small. The same can be said of the standard deviation of the f_{NL} scatter between realizations with and without extra noise, which reaches at most a value σ_{ΔfNL} ≈ σ_{fNL}/3, for a large A = 0.3. These results provide a good indication that a noise mismatch between simulations and data is not a concern for primordial NG estimation, unless the mismatch itself is well above the estimated percent level in the noise power spectrum and produces large correlations between small and large scales. It is also worth mentioning that, in the very early stages of our primordial NG analysis of Planck data, when accurate simulations were not yet available, we calibrated our estimators using very simple noise models (e.g. we initially generated noise in harmonic space, with a nonflat power spectrum consistent with the data, but neglecting any correlations between different scales; we then went to pixel space and modulated the noise map with the hitcount map to anisotropize it). Despite the simplicity of this approach, we were able to verify later, using FFP simulations, that such simple models already produced an accurate linearterm correction. This further reinforces our confidence in the robustness of f_{NL} estimators to imperfect modelling of the noise.
6.3. Effects of foregrounds
Here we look at two nonprimordial contributions to the Planck bispectrum, namely Galactic dust and the SunyaevZeldovich effect, which could potentially be present in the cleaned maps, but which we do not detect in the end. Another potential contaminant for both primordial and lensing bispectrum results is the “intrinsic bispectrum,” induced in the CMB by weak (secondorder) nonlinearities from gravity in general relativity and by nonlinearities in the recombination physics. This would set the minimal level of CMB NG present even for Gaussian initial conditions of the primordial curvature perturbation. However, this is of no particular concern to us, since for the Planck data set its expected impact is very small both in temperature (Bartolo et al. 2004b, 2005, 2010c, 2012; Creminelli & Zaldarriaga 2004; Boubekeur et al. 2009; Nitta et al. 2009; Bartolo & Riotto 2009; Senatore et al. 2009; Khatri & Wandelt 2009, 2010; Creminelli et al. 2011; Huang & Vernizzi 2013; Su et al. 2012; Pettinari et al. 2014) and polarization (Lewis 2012; Pettinari et al. 2014). Nonprimordial contributions that we do detect (lensing and extragalactic point sources) were discussed in Sect. 4.
6.3.1. NonGaussianity of the thermal dust emission
In principle, the foregroundcleaned maps should not contain any noticeable NG of Galactic origin. However, in raw observations the strongest contamination to the primordial NG is due to the thermal dust emission, which induces a large negative bias in the measurements of . Therefore, it is important to verify that this contamination has been removed entirely through the different componentseparation methods.
There is no analytical template for the dust bispectrum, unlike the extragalactic templates discussed in Sect. 4.2. A simple method using instead numerical templates for the different Galactic foregrounds has been described in Jung et al. (2018) (see also Coulton & Spergel 2019). First, using the Binned bispectrum estimator, the dust bispectral shape is computed from the thermal dust emission map at 143 GHz (the dominant frequency channel in the cleaned CMB maps) produced by the Commander technique. This is actually a map determined at higher frequencies, where the dust dominates, and then rescaled to 143 GHz. Since there was no improvement for the temperature map of this foreground in the latest release, we use the Planck 2015 temperature map here. Then, this numerical dust bispectrum can be used as a theoretical template in the analysis of any other data map with the Binned bispectrum estimator. The only condition is that the same mask, beam, and binning are used for both the determination of the dust bispectrum and the analysis itself.
An illustration of the large bias induced by the presence of dust in the map is given in Table 19. It shows the f_{NL} parameters of the primordial local shape and the thermal dust bispectrum for both an independent and a joint analysis of the raw 2018 143GHz Planck temperature map. As a reminder, in an independent analysis we assume that only one of the templates is present in the data (and we repeat the analysis for each individual template), while in a fully joint analysis we assume that all templates are present, so that we have to take their correlations into account. In the independent case, there is a strong detection of local NG (at about the 5σ level), while the dust is observed at the expected level (within the 1σ interval centred on the expected value: ). The joint analysis shows that indeed the large negative of the independent case is entirely due to the presence of dust in the map, while the dust is still detected at the expected amount.
Independent and joint estimates (see the main text for a definition) of the f_{NL} parameters of the primordial local shape and the thermal dust bispectrum in the raw Planck 143 GHz temperature map, determined using the Binned bispectrum estimator.
The analysis of the Tonly maps produced by the four componentseparation methods is the main result of this subsection. Table 20 gives the values of f_{NL} for several primordial NG shapes (local, equilateral, and orthogonal), the amplitudes of some extragalactic foreground bispectra (unclustered point sources and the CIB, as defined in Sect. 4.2) and the f_{NL} of the thermal dust emission bispectrum, after subtracting the lensing bias, in both an independent and a fully joint analysis. There is no significant detection of dust in any of the four maps (the worst case being SMICA with slightly more than a 1σ deviation)^{11}. Given the small size of the errors, this nondetection means that there is at most a few percent of dust contamination in the cleaned maps (outside the mask). However, the errors of the local shape increase significantly in the joint analysis because the dust and the local shapes are quite correlated (more than 60%).
Independent and joint estimates of the f_{NL} parameters of the indicated shapes, including in particular the dust template, for the cleaned maps produced by the four componentseparation methods, as determined with the Binned bispectrum estimator.
6.3.2. Impact of the tSZ effect
The SMICA componentseparation method also produces a foregroundcleaned temperature map, “SMICA noSZ,” where, in addition to the usual foregrounds, contamination by the thermal SunyaevZeldovich (tSZ) effect has been subtracted as well (see Planck Collaboration IV 2020). This allows us to test if the tSZ contamination has any significant impact on our primordial results (we already saw in Sect. 4.1 that it does seem to have an impact on the determination of the lensing NG). This is an important test, since there have been recent claims in the literature (Hill 2018) that it might have an effect.
The results of the analysis can be found in Table 21. Because this effect is only important in temperature, we restrict ourselves to a Tonly analysis. Results have been determined with the KSW estimator, with the Binned estimator (this time using only 150 maps for the linear correction and the error bars instead of 300, which is enough for the purposes of this test), and with the Modal 1 estimator. The mask used is the same as for the main analysis. The table also contains the difference with the result determined from the normal SMICA temperature map (without tSZ removal) and the uncertainties on this difference.
Impact of the thermal SunyaevZeldovich (tSZ) effect on f_{NL} estimators for temperature data. First three columns: results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the indicated estimators from the SMICA foregroundcleaned temperature map with additional removal of the tSZ contamination. Last three columns: difference with the result from the normal SMICA map (see Table 5), where the error bar is the standard deviation of the differences in f_{NL} for the Gaussian FFP10 simulations when processed through the two different SMICA foregroundcleaning procedures. Results have been determined using an independent singleshape analysis and are reported with subtraction of the lensing bias; error bars are 68% CL.
We see a shift of about 1σ_{ΔfNL} (hence insignificant) in the local shape result, which together with orthogonal is the shape predicted to be most contaminated by the tSZ effect (see Hill 2018). We actually see a larger shift in the equilateral result, which is supposed to be almost unaffected by this effect, while the orthogonal shift is the largest, at more than 2σ_{ΔfNL} for all estimators. However, such a marginal effect in the orthogonal shape without a corresponding effect in the local shape leads us to conclude that we do not detect any significant impact of contamination of the usual foregroundcorrected maps by the tSZ effect on our primordial NG results. In other words, while some tSZ contamination, peaking in the squeezed limit, is expected to be present in the standard temperature maps, this is too small to be clearly disentangled from the statistical fluctuations in the f_{NL} results; in other words the tSZ contamination is not bigger than effects due to the different processing of the data when tSZ is included in the foreground components for the SMICA analysis. Furthermore, all shifts discussed here are much smaller than the uncertainties on the f_{NL} values themselves, due to the fact that the f_{NL} scatter between different cleaned maps (σ_{ΔfNL}) is much smaller than the f_{NL} error (σ_{fNL}).
Finally, it should also be pointed out that, using the same criteria, we cannot strictly call the shift in discussed in Sect. 4.1 significant either. The observation that all TTT measurements are below the expected value and that adding polarization systematically shifts them up, does point to some tSZ contamination in Tonly results for the lensing bispectrum. However, the analysis discussed in this section finds again only a 2σ_{ΔfNL} effect in this case; this is slightly larger or smaller than the significance of the orthogonal f_{NL} shift, depending on the estimator. Again, intrinsic statistical uncertainties make it hard to detect this systematic effect.
6.4. Dependence on sky coverage
The temperature and polarization mask we are using have been determined to be the optimal masks for use on the maps produced by the componentseparation pipelines, according to criteria explained in Planck Collaboration IV (2020), and hence are used for all Planck analyses on those maps. However, the choice of mask can have an impact on the results for f_{NL}. In the first place the sky fraction of the mask has a direct effect on the size of the uncertainties. Potentially more important, however, is the effect the mask might have on the amount of foreground residuals. Hence we judge it important to investigate the impact of the choice of mask on our results by comparing results for several different masks. All tests are performed on the SMICA maps (with one exception, detailed below) using the Binned bispectrum estimator, using 150 maps for the linear correction and the errors.
As a first test, performed on the temperature map only, we take the union of the mask used in this paper (f_{sky} = 0.78) with the mask we used in our 2015 analysis (f_{sky} = 0.76), leading to a mask that leaves a fraction f_{sky} = 0.72 of the sky uncovered. The results for the standard primordial shapes are given in Table 22, as well as their differences with the results using the standard mask (see Table 5). The errors on the differences have been determined from the scatter among 150 simulations when analysed with the two different masks. This particular test is performed both on the SMICA and the Commander maps, since its initial purpose was a further check on the discrepancy between those two componentseparation methods regarding the Tonly orthogonal result, discussed in detail in Sect. 6.1.1. We see, however, that the orthogonal result is very stable for both methods, while the local result shifts by about 1σ_{ΔfNL} and the equilateral result moves by about 2σ_{ΔfNL}, which corresponds to about (2/3)σ_{fNL}. Checking the 150 Gaussian simulations used for determining the errors and the linear correction, we see that there are 20 that have at least one fluctuation larger than 2σ_{ΔfNL} in at least one of the three shapes, which corresponds to a 13.3% probability. This is large enough that we consider the shift consistent with a statistical fluctuation. In any case, all results remain consistent with zero. It is interesting to note that this result for the equilateral shape is much closer to the one determined for the 2015 Planck data.
Tests of dependence on choice of mask.
As a second test, to check the impact of the size of the polarization mask, we return to the standard temperature mask, but this time we change the polarization mask. It is altered as follows: each hole in the common mask is grown by a region 20 pixels in width. This reduces the sky fraction from 0.78 to 0.73. Results for this test can be found in Table 23, including the differences with the results determined with the standard mask. We see that many results shift around somewhat, but nothing appears very significant (all less than 1.5σ_{ΔfNL}). The largest is again a (2/3)σ_{fNL} shift for the Eonly equilateral case, moving it closer to zero. This test was also performed with the KSW and Modal 1 estimators, giving consistent (but not identical) results. In particular, while some values shift a bit more, the shift in the equilateral shape is smaller for both the KSW and Modal 1 estimators. However, all estimators agree on the signs of the shifts.
Tests of growing the polarization mask size.
Our third and final test of the effects of mask choice is very similar to the previous one, except that this time the polarization mask is enlarged by an additional 40 (instead of 20) pixels around every hole. This further reduces f_{sky} to 0.66. Results can also be found in Table 23, and we find similar conclusions as for the previous test.
To summarize, while we see some effects on our f_{NL} results when considering different masks, none of these appear to be significant (further reinforced by the fact that the shifts vary somewhat between estimators), nor do they change the conclusions of our paper.
6.5. Dependence on multipole number
As in previous releases we also test the dependence of the results for f_{NL} on the choice of ℓ_{min} and ℓ_{max} used in the analysis. This test is most easily performed using the Binned estimator. Results are shown in Fig. 17 for the dependence on ℓ_{min} and in Fig. 18 for the dependence on ℓ_{max}.
Fig. 17. Evolution of the f_{NL} parameters (solid blue line with data points) and their uncertainties (dotted blue lines) for the three primordial bispectrum templates as a function of the minimum multipole number ℓ_{min} used in the analysis. From left to right the panels show local, equilateral, and orthogonal shape results, while the different rows from top to bottom show results for T only, E only, and full T+E data. To indicate more clearly the evolution of the uncertainties, they are also plotted around the final value of f_{NL} (solid green lines without data points, around the horizontal dashed green line). The results here have been determined with the Binned bispectrum estimator for the SMICA map, assume all shapes to be independent, and the lensing bias has been subtracted. 
Fig. 18. Same as Fig. 17, but this time as a function of the maximum multipole number ℓ_{max} used in the analysis. 
Considering first Fig. 17, for ℓ_{min}, the plots look very similar to the ones in the paper investigating the 2015 Planck data (PCNG15), with increased stability for the Eonly local results. As explained in Sect. 3.2.2, it was decided to use ℓ_{min} = 4 for the polarization maps, since there was an issue with bias and increased variance when the two lowest multipoles for E were included. However, this still allows us to use many more lowℓ polarization modes than for the 2015 data, when it was necessary to adopt ℓ_{min} = 40 for polarization.
Turning to Fig. 18, for ℓ_{max}, we also notice good agreement with the analysis of the 2015 data, and slightly more stable results (e.g. for Tonly equilateral). We see that the Tonly results have stabilized by ℓ = 2000 and the Eonly results by ℓ = 1500, so that there is no problem with the KSW and Modal estimators using these lower values for ℓ_{max}. As before we also confirm the “WMAP excess” for the local shape at ℓ ≈ 500 (Bennett et al. 2013), even more clearly than in PCNG15. However, this does not appear to be a major outlier when considering all values of ℓ_{max} and all possible shapes.
6.6. Summary of validation tests
Throughout this section we have discussed a set of tests aimed at evaluating the robustness of our results. For convenience, we summarize here our main findings.

We find good consistency for f_{NL} local, equilateral, and orthogonal measurements, between all componentseparation methods and with all bispectrum estimators, separately considering Tonly, Eonly, and T+E results. The agreement for Eonly results has significantly improved with respect to the previous release.

The possible exception is provided by the orthogonal Tonly estimates. In this case, a comparison with simulations shows significant differences in specific cases, namely SMICA – Commander and, to a lesser extent, SMICA – SEVEM. However, these differences are coming from fluctuations in f_{NL} for Commander and SEVEM with respect to the two previous Planck releases, when all temperature maps were always in excellent agreement. Our main SMICA results are, on the other hand, completely stable. Moreover, the discrepancy becomes much less significant when adding polarization, it is limited to one specific shape, and it is in any case at a level that does not alter the interpretation of the results. Therefore, in the end, we do not consider this issue to be problematic.

Nevertheless, in light of this test, we are led to a slight preference for SMICA and NILC as methods of choice for primordial NG analysis. Given that SMICA also gave a slightly better performance for NG analysis in the two previous releases, and that results extracted from SMICA maps have been quite stable over time, we maintain SMICA as our final choice, as already justified in detail in the papers analysing the 2013 and 2015 Planck data (PCNG13; PCNG15).

We find very good consistency between fully reconstructed bispectra, for different componentseparation methods, using both the Modal and the Binned approaches. Polarizationonly bispectra again show a large improvement compared to the previous analysis of Planck data.

The observed noise mismatch between the data and the FFP10 simulations does not seem to impact our results. Our cubic statistics cannot be biased by this mismatch, and tests on simulations show that the effect on f_{NL} errors is negligible.

Results are stable to changes in sky coverage and different cuts in the multipole domain. Restricting the analysis to the multipole range probed by WMAP shows good agreement between WMAP and Planck.

We find no sign of any residual Galactic thermal dust contamination in the Planck componentseparated CMB maps.

Contamination from the thermal SZ effect on the standard primordial and lensing Tonly results is not at a significant level, compared to statistical errors.
Overall, the results display a high level of internal consistency, and are notably characterized by a large improvement in the quality of polarizationonly bispectra with respect to the previous release. Whereas in 2015 we cautioned the reader to take polarizationbased f_{NL} estimates as preliminary, we can now state that our T+Ebased constraints are fully robust. This is one of the main conclusions of this paper.
7. Limits on the primordial trispectrum
We now present constraints on three shapes for the primordial fourpoint function or trispectrum, denoted , , and , and describe them below. The details of the analysis are mostly unchanged from the paper analysing the 2015 Planck data (see Sect. 9 of PCNG15). Nevertheless, we briefly review the fourpoint analysis here; for more details, see PCNG15, or Smith et al. (2015), which contains technical details of the pipeline.
First, we describe the three different signals of interest. The localtype trispectrum arises if the initial adiabatic curvature ζ is given by the following nonGaussian model:
where ζ_{G} is a Gaussian field. The trispectrum in the local model is given by
In this equation and throughout this section, a “primed” fourpoint function denotes the fourpoint function without its momentumconserving delta function, that is,
where “+ disc.” denotes disconnected contributions to the 4point function. It can be shown that the localtype trispectrum in Eq. (43) is always negligibly small in singlefield inflation (Senatore & Zaldarriaga 2012a); however, it can be large in multifield models of inflation in which a large bispectrum is forbidden by symmetry (Senatore & Zaldarriaga 2012b).
The next two shapes , are generated by the operators and (∂_{i}σ)^{2}(∂_{j}σ)^{2} in the effective field theory (EFT) of inflation (Bartolo et al. 2010b; Senatore & Zaldarriaga 2012b; Smith et al. 2015). For data analysis purposes, they can be defined by the following trispectra:
Here K = ∑_{i}k_{i}, and numerical prefactors have been chosen so that the trispectra have the same normalization as the local shape in Eq. (43) when restricted to tetrahedral wavenumber configurations with k_{i}=k and (k_{i} ⋅ k_{j}) = − k^{2}/3. We mention in advance that there is another shape that arises at the same order in the EFT expansion, to be discussed at the end of this section.
In Eqs. (45) and (46), we have written each trispectrum in two algebraically equivalent ways, an integral representation and an “integrated” form. The integral representation arises naturally when evaluating the Feynman diagram for the EFT operator. It also turns out to be useful for data analysis, since the resulting “factorizable” representation for the trispectrum leads to an efficient algorithm for evaluating the CMB trispectrum estimator (Planck Collaboration XVII 2016; Smith et al. 2015).
For simplicity in Eqs. (45) and (46), we have assumed a scaleinvariant power spectrum P_{ζ}(k) = A_{ζ}/k^{3}. For the Planck analysis, we slightly modify these trispectra to a powerlaw spectrum P_{ζ}(k)∝k^{ns − 4}, as described in Appendix C of Smith et al. (2015).
To estimate each g_{NL} parameter from Planck data, we use the “pureMC” trispectrum estimation pipeline from section IX.B of Smith et al. (2015). In this pipeline, the data are specified as a filtered harmonic space map , and its covariance is characterized via a set of 1000 filtered signal + noise simulations.
The “filter” is an experimentspecific linear operation whose input consists of one or more pixelspace maps, and whose output is a single harmonicspace map, d_{ℓm}. In the Planck trispectrum analysis, we define the filter as follows. First, we take the singlefrequency pixelspace maps, and combine them to obtain a single componentseparated pixelspace map, using one of the componentseparation algorithms, SMICA, SEVEM, NILC, Commander, or SMICA noSZ (a variant of the SMICA algorithm that guarantees zero response to Comptony sources, at the expense of slightly higher noise; see Sect. 6.3.2). Second, we subtract the bestfit monopole and dipole, inpaint masked point sources, and apodize the Galactic plane boundary. The details of these steps are unchanged from the 2015 Planck data analysis, and are described in Sect. 9.1 of PCNG15. Third, we take the spherical transform of the pixelspace map out to l_{max} = 1600, obtaining a harmonicspace map d_{ℓm}. Finally, we define the filtered map by applying the multiplicative factor:
This sequence of steps defines a linear operation, whose input is a set of singlefrequency pixelspace maps, and whose output is a filtered harmonicspace map . This filtering operation is used as a building block in the trispectrum pipeline described in Smith et al. (2015), and is the only part of the pipeline that is Planckspecific.
The results of the analysis, for all three trispectrum shapes and five different componentseparation algorithms, are presented in Table 24. We do not find evidence for a nonzero primordial trispectrum.
Planck 2018 constraints on the trispectrum parameters , , and from different componentseparated maps.
Each entry in Table 24 is a constraint on a single g_{NL}parameter with the others held fixed. We next consider joint constraints involving multiple g_{NL}parameters. In this case, we need to know the covariance matrix between g_{NL}parameters. We find that
and the correlation between and the other two g_{NL}parameters is negligible.
Multifield models of inflation generally predict a linear combination of the quartic operators and (∂_{i}σ)^{2}(∂_{j}σ)^{2}. In addition, there is a third operator that arises at the same order in the EFT expansion. For completeness, its trispectrum is given by
However, a Fisher matrix analysis shows that this trispectrum is nearly 100% correlated with the and (∂_{i}σ)^{2}(∂_{j}σ)^{2} trispectra. If the parameter is nonzero, then we can absorb it into the “effective” values of the parameters and as
Therefore, to study joint constraints involving multiple g_{NL} parameters, it suffices to consider the two parameters, and , with correlation coefficient given in Eq. (48).
We define the twocomponent parameter vector
and let denote the twocomponent vector of singleg_{NL} estimates from the SMICA maps (Table 24):
We also define a twobytwo Fisher matrix F_{ij}, whose diagonal is given by , where σ_{i} is the singleg_{NL} statistical error in Table 24, and whose offdiagonal is , where r is the correlation in Eq. (48). This gives:
Now, given a set of “theory” g_{NL} values, represented by a twovector g_{i}, we compare to the Planck data by computing the following quantity for the trispectrum:
In a modelbuilding context where the g_{NL} quantities g_{i} depend on model parameters, confidence regions on model parameters can be obtained by appropriately thresholding χ^{2}. We give some examples in Sect. 8.
8. Implications for earlyUniverse physics
We now want to convert constraints on primordial NG into constraints on parameters of various models of inflation. This allows us to highlight the constraining power of NG measurements, as an additional complementary observable beyond the CMB power spectra. In particular NG constraints can severely limit the parameter space of models that are alternatives to the standard singlefield models of slowroll inflation, since they typically feature a higher level of NG.
Unless stated otherwise, we follow the same procedures adopted in PCNG13 and PCNG15. A posterior of the model parameters is built based on the following steps: we start from the assumption that 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); we use uniform or Jeffreys’ priors, over intervals of the model parameter space that are physically meaningful (or as otherwise stated); and in some cases where two or more parameters are involved, we marginalize the posterior to provide onedimensional limits on the parameter under consideration.
8.1. General singlefield models of inflation
DBI models. DiracBornInfeld (DBI) models of inflation (Silverstein & Tong 2004; Alishahiha et al. 2004) arise from highenergy stringtheory constructions and generate a nonlinearity parameter , where c_{s} is the sound speed of the inflaton perturbations (Silverstein & Tong 2004; Alishahiha et al. 2004; Chen et al. 2007b). The enhancement of the NG amplitude due to a possible sound speed c_{s} < 1 arises from a nonstandard kinetic term of the inflaton field. Notice that we have constrained the exact theoretical (nonseparable) shape (see Eq. (7) of PCNG13), even though it is very similar to the equilateral type. Our constraint from temperature data ( from temperature and polarization) at 68% CL (with lensing and point sources subtracted, see Table 11) implies
and
Implications for the effective field theory of inflation. Now we can update CMB limits on the speed of sound c_{s} at which inflaton fluctuations propagated in the very early Universe. A very general constraint on this inflationary parameter can be obtained by employing the EFT approach to inflation (Cheung et al. 2008; Weinberg 2008, and see Sect. 7). This approach allows us to obtain predictions for the parameter space of primordial NG through a general characterization of the inflaton field interactions. The Lagrangian of the system is expanded into the dominant operators that respect some underlying symmetries. The procedure thus determines a unifying scheme for classes of models featuring deviations from singlefield slowroll inflation. Typically the equilateral and orthogonal templates represent an accurate basis to describe the full parameter space of EFT singlefield models of inflation, and therefore we used the constraints on and .
As a concrete example, let us consider the Lagrangian of general singlefield models of inflation (of the form P(X, φ) models, where X = g^{μν}∂_{μ}ϕ ∂_{ν}ϕ) written with the EFT approach:
The scalar perturbation π generates the curvature perturbation ζ = −Hπ. In this case there are two relevant inflaton interactions, and , producing two specific bispectra with amplitudes and , respectively.
Here M_{3} is the amplitude of the operator (see Senatore et al. 2010; Chen et al. 2007b; Chen 2010b), with the dimensionless parameter (Senatore et al. 2010). The two EFT shapes can be projected onto the equilateral and orthogonal shapes, with the mean values of the estimators for and expressed in terms of c_{s} and as
where the coefficients come from the Fisher matrix between the theoretical bispectra predicted by the two operators and and the equilateral and orthogonal templates. Notice that DBI models reduce to the condition , while the noninteracting (vanishing NG) case corresponds to c_{s} = 1 and M_{3} = 0 (or .
We then proceed as in the two previous analyses (PCNG13; PCNG15). We employ a χ^{2} statistic computed as , with (i = {equilateral, orthogonal}), where are the joint estimates of equilateral and orthogonal f_{NL} values (see Table 6), while are provided by Eq. (58) and C is the covariance matrix of the joint estimators. Figure 19 shows the 68%, 95%, and 99.7% confidence regions for and , as derived from from the T + E constraints, with the requirement χ^{2} ≤ 2.28, 5.99, and 11.62, respectively (corresponding to a χ^{2} variable with two degrees of freedom). In Fig. 20 we show the corresponding confidence regions in the parameter space. Marginalizing over we find
Fig. 19. 68%, 95%, and 99.7% confidence regions in the parameter space , defined by thresholding χ^{2}, as described in the text. 
Fig. 20. 68%, 95%, and 99.7% confidence regions in the singlefield inflation parameter space , obtained from Fig. 19 via the change of variables in Eq. (58). 
and
There is a slight improvement in comparison with the constraints obtained in PCNG15 coming from the T + E data. Notice that, contrary to the situation in 2015, the final c_{s} constraint here is essentially determined using temperature data, with no significant improvement when adding the polarization data set. There are two reasons for this. Firstly, the errors on and display less improvement in going from T to T + E, in comparison to the previous release (see Sect. 9 for further comments on this point). Secondly, the shift of the central value of (and also of ) when passing from T to T + E is somewhat larger than in 2015. Such a shift is now towards more negative values and has the net effect of counterbalancing the improvement of the error bars on and when adding polarization, thus leaving the constraints on c_{s} unchanged with respect to the one from temperature data only.
8.2. Multifield models
Constraints on primordial NG of the local type lead to strong implications for models of inflation where scalar fields (different from the inflaton) are dynamically important for the generation of the primordial curvature perturbation. In the following we test two scenarios for curvaton models.
Basic curvaton models. The simplest adiabatic curvaton models predict primordial NG of the local shape with a nonlinearity parameter (Bartolo et al. 2004c,d)
in the case where the curvaton field has a quadratic potential (Lyth & Wands 2002; Lyth et al. 2003; Lyth & Rodriguez 2005; Malik & Lyth 2006; Sasaki et al. 2006). Here the parameter r_{D} = [3ρ_{curvaton}/(3ρ_{curvaton} + 4ρ_{radiation})]_{D} is the “curvaton decay fraction” at the time of the curvaton decay in the sudden decay approximation. We assume a uniform prior, 0 < r_{D} < 1. It is worth recalling that these models predict a lower bound for the level of NG, of the order of unity (corresponding precisely to ), which is considered as a typical threshold to distinguish between standard singlefield and multifield scenarios. Our constraint at 68% CL (see Table 6) implies
The constraint at 68% CL obtained from temperature and polarization data yields the constraint
These limits indicate that in such scenarios the curvaton field has a nonnegligible energy density when it decays. Meaningful improvements are achieved with respect to previous bounds in PCNG15, namely an almost 20% improvement from Tonly data and a 10% improvement when including Emode polarization.
Decay into curvaton particles. We reach similar improvements on the parameters of the second scenario of the curvaton models we consider. In this case one accounts for the possibility that the inflaton field can decay into curvaton particles (Linde & Mukhanov 2006), a possibility that is neglected in the above expression (Eq. (61)) for . It might be the case that the classical curvaton field survives and begins to dominate. In this case also the curvaton particles produced during reheating are expected to survive and dominate over other species at the epoch of their decay (since thay have the same equation of state as the classical curvaton field). Primordial adiabatic perturbations are generated given that the classical curvaton field and the curvaton particles decay at the same time (see Linde & Mukhanov 2006 for a detailed discussion). To interpret f_{NL} in this scenario we employ the general formula for derived in Sasaki et al. (2006), which takes into account the possibility that the inflaton field decays into curvaton particles:
Here the parameter is the ratio of the energy density of curvaton particles to the energy density of the classical curvaton field (Linde & Mukhanov 2006; Sasaki et al. 2006), while now ρ_{curv} in the expression for r_{D} must be replaced by the sum of the densities of the curvaton particles and curvaton field. As in PCNG15 we use uniform priors 0 < r_{D} < 1 and . Our limits on constrain
and
which does not exclude a contribution of curvaton particles comparable to the one from the classical curvaton field.
8.3. Nonstandard inflation models
Directionaldependent NG. Table 13 shows the constraints on directionallydependent bispectra (Eq. (20)). This kind of NG is predicted by several different inflationary models. For example, it is a robust and (almost unavoidable) outcome of models of inflation where scaleinvariant gauge fields are present during inflation. As summarized in Sect. 2 they are also produced from partially massless higherspin particles (Franciolini et al. 2018) or from models of solid inflation (Endlich et al. 2013, 2014; Shiraishi et al. 2013a), as well as in models of inflation that break both rotational and parity invariance (Bartolo et al. 2015). To compare with the constraints obtained in the analysis of the 2015 Planck data (PCNG15), we reconsider the specific model where the inflaton is coupled to the kinetic term F^{2} of a gauge field via a term ℒ = − I^{2}(ϕ)F^{2}/4, where I(ϕ) is a function that depends on the inflaton field, having an appropriate time evolution during inflation (see, e.g. Ratra 1992). Specifically in these models the production of superhorizon vector field perturbations switches on the L = 0 and L = 2 modes in the bispectrum, with nonlinearity parameters , with X_{L = 0} = (80/3) and X_{L = 2} = − (10/6), respectively (Barnaby et al. 2012a; Bartolo et al. 2013a; Shiraishi et al. 2013a). In these expressions g_{*} is a parameter that measures the amplitude of a quadrupolar anisotropy in the power spectrum (see, e.g. Ackerman et al. 2007), while N is the number of efolds (from the end of inflation) at which the relevant scales cross outside the Hubble scale. It is therefore interesting to set some limits on the parameter g_{*} exploiting the constraints from primordial NG of this type. Using the SMICA constraints from T (or T+E) in Table 13, marginalizing over a uniform prior 50 ≤ N ≤ 70, and assuming uniform priors on −1 ≤ g_{*} ≤ 1, we obtain the 95% bounds −0.041 < g_{*} < 0.041 (−0.036 < g_{*} < 0.036), and −0.35 < g_{*} < 0.35 (−0.30 < g_{*} < 0.30), from the L = 0 and L = 2 modes, respectively (considering g_{*} to be scale independent).
Tensor NG and pseudoscalars. Using the SMICA T+E result (68% CL), we here place constraints on two specific inflation models, including either a U(1)axion coupling or an SU(2)axion one. The former U(1) model results in , where 𝒫 is the vacuummode curvature power spectrum, ϵ is a slowroll parameter of the inflaton field, and ξ expresses the strength of the U(1)axion coupling (Cook & Sorbo 2013; Shiraishi et al. 2013b). We then fix ϵ to be 0.01 and marginalize 𝒫 with the prior, 1.5 × 10^{−9} < 𝒫 < 3.0 × 10^{−9}; assuming a prior, 0.1 < ξ < 7.0, the upper bound on ξ is derived as ξ < 3.3 (95% CL). In the latter SU(2) model, under one specific condition, the tensor nonlinear parameter is related to the tensortoscalar ratio r and the energy density fraction of the gauge field Ω_{A} as (Agrawal et al. 2018). The lower bound on Ω_{A} is estimated under a prior, 0 < Ω_{A} < 1. We find Ω_{A} > 2.3 × 10^{−7} and 2.7 × 10^{−9} (95% CL) for r = 10^{−2} and 10^{−3}, respectively.
Warm inflation. As in previous analyses (PCNG13; PCNG15), we adopt the expression (Moss & Xiong 2007). This is valid when dissipative effects are strong, this is for values r_{d} ≳ 2.5 of the dissipation parameter r_{d} = Γ/(3H) (measuring the effectiveness of the energy transfer from the inflaton field to radiation)^{12}. Assuming a constant prior 0 ≤ log_{10}r_{d} ≤ 4, the SMICA constraints (at 68% CL) from T and from T+E (see Table 12), yield log_{10}r_{d} ≤ 3.6 and log_{10}r_{d} ≤ 3.5, respectively, at 95 % CL. The results show that strongdissipative effects in warm inflation models remain allowed. This is a regime where gravitino overproduction problems can be evaded (for a discussion see Hall & Peiris 2008).
8.4. Alternatives to inflation
As an example we update the constraints on some ekpyrotic/cyclic models (e.g. see Lehners 2010 for a review). Typically local NG is produced through a conversion of “intrinsic” nonGaussianity in the entropy fluctuations into the curvature perturbation. This conversion can proceed in different ways. The “ekpyrotic conversion” models, for which the conversion acts during the ekpyrotic phase, have already been ruled out (Koyama et al. 2007, PCNG13). On the other hand, in the “kinetic conversion” models the conversion takes place after the ekpyrotic phase and a local bispectrum is generated with an amplitude ^{13}. The sign depends on the details of the conversion process (Lehners & Steinhardt 2008, 2013; Lehners 2010), and typical values of the parameter ϵ are ϵ ≈ 50 or greater. Assuming ϵ ≈ 100 and using a uniform prior on −5 < κ_{3} < 5 the constraints on from T only (see Table 6), implies −1.1 < κ_{3} < 0.36 and −0.43 < κ_{3} < 1.0 at 95% CL, for the plus and minus sign in , respectively. The T+E constraints on (Table 6) yield −1.05 < κ_{3} < 0.27 and −0.38 < κ_{3} < 0.94 at 95% CL, for the plus and minus sign, respectively. If we take ϵ ≈ 50 as an example, we obtain the following limits: −1.6 < κ_{3} < 0.51 and −0.62 < κ_{3} < 1.5 at 95% CL from T only; and −1.5 < κ_{3} < 0.39 and −0.54 < κ_{3} < 1.3 at 95% CL from T+E constraints.
8.5. Inflationary interpretation of CMB trispectrum results
We briefly analyse inflationary implications of the Planck trispectrum constraints, using the SMICA limits on g_{NL} parameters from Table 24.
First, we consider singlefield inflationary models, using the effective action for the Goldstone boson π (see, e.g., Smith et al. 2015):
The mass scale M_{4} is related to our previouslydefined g_{NL} parameters by:
Therefore, using the limit from Table 24, we get the following constraint on singlefield models:
Next consider the case of multifield inflation. Here, we consider an action of the more general form:
where σ is a light field that acquires quantum fluctuations with power spectrum P_{σ}(k) = H^{2}/(2k^{3}). We assume that σ converts to adiabatic curvature ζ, i.e. ζ = (2A_{ζ})^{1/2}H^{−1}σ. The model parameters Λ_{i} are related to our previouslydefined g_{NL} parameters by:
and can be constrained by thresholding the χ^{2}statistic defined in Eq. (54). For example, to constrain the parameter Λ in the Lorentz invariant model
we set obtain g_{NL}values using Eq. (71), and use Eq. (54) to obtain χ^{2} as a function of Λ. We then threshold at Δχ^{2} = 4 (as appropriate for one degree of freedom), to obtain the following constraint:
DBI trispectrum. We use the trispectrum constraints on the shape in Table 24 to determine a lower bound on the sound speed of the inflaton field in DBI models. In fact in these models the dominant contribution in the smallsoundspeed limit (Chen et al. 2009; Arroja et al. 2009) to the contact interaction trispectrum (Huang & Shiu 2006) produces such a shape, with an amplitude . We employ the same procedure described at the beginning of this section and, assuming a uniform prior in the range 0 ≤ c_{s} ≤ 1/5, we derive the following constraint on c_{s}:
This constraint is independent from and consistent with the bounds of Eqs. (55) and (56) obtained from the bispectrum measurements. Notice that in the trispectrum case we are ignoring the scalar exchange contribution, which turns out to be of the same order in c_{s}.
Curvaton trispectrum. A generic prediction of the simplest adiabatic curvaton scenario is also a localtype trispectrum with an amplitude given by (Sasaki et al. 2006)
Following the procedure described at the beginning of this section, we use the observational constraint obtained in Sect. 7 (see Table 24), and the same uniform prior (0 < r_{D} < 1) as in Sect. 8.2, to obtain a lower bound on the curvaton decay fraction
This limit is consistent with the previous ones derived using the bispectrum measurements and it is about a factor of 4 weaker.
9. Conclusions
In this paper, we have presented constraints on primordial NG, using the Planck fullmission CMB temperature and Emode polarization maps. Compared to the Planck 2015 release we now include the lowℓ (4 ≤ ℓ < 40) polarization multipole range.
Our analysis produces the following final results (68% CL, statistical): ; ; and . These results are overall stable with respect to our constraints from the 2015 Planck data. They show no real improvement in errors, despite the additional polarization modes. This is due to a combination of two factors. Firstly, the local shape, which is most sensitive to lowℓ modes and where one would naively expect an improvement, is actually less sensitive to polarization than the equilateral and orthogonal shapes. This means that in the end none of the three shapes are very sensitive to lowℓ polarization modes. Secondly, the temperature and polarization simulations used to determine the errors have a more realistic but slightly higher noise level than in the previous release.
On the other hand, the quality of polarization data shows a clear improvement with respect to our previous analysis. This is confirmed by a large battery of tests on our data set, including comparisons between different estimator implementations (KSW, Binned, and two Modal estimators) and foregroundcleaning methods (SMICA, SEVEM, NILC, and Commander), studies of robustness under changes in sky coverage and multipole range, and an analysis of the impact of noiserelated systematics. While in our previous release we cautioned the reader to take polarization bispectra and related constraints as preliminary, in light of these tests we now consider our results based on the combined temperature and polarization data set to be fully reliable. This also implies that polarizationonly, EEE bispectra can now be used for independent tests, leading to primordial NG constraints at a sensitivity level comparable to that of WMAP from temperature bispectra, and yielding statistical agreement.
As in the previous analyses, we go beyond the local, equilateral, and orthogonal f_{NL} constraints by considering a large number of additional cases, such as scaledependent feature and resonance bispectra, running f_{NL} models, isocurvature primordial NG, and paritybreaking models. We set tight constraints on all these scenarios, but do not detect any significant signals.
On the other hand, the nonprimordial lensing bispectrum is now detected with an improved significance compared to 2015, excluding the null hypothesis at 3.5σ. The amplitude of the signal is consistent with the expectation from the Planck bestfit cosmological parameters, further indicating the absence of significant foreground contamination or spurious systematic effects. We also explicitly checked for the presence of various nonprimordial contaminants, like unclustered extragalactic point sources, CIB, Galactic thermal dust, and the thermal SZ effect, but apart from the first, none of these were detected. The small amount of remaining pointsource signal in the cleaned maps has no impact on our other constraints because of its negligible correlations.
We update our trispectrum constraints, now finding (68% CL, statistical), while also constraining additional shapes, generated by different operators in an effective fieldtheory approach to inflation.
In addition to estimates of bispectrum and trispectrum amplitudes, we produce modelindependent reconstructions and analyses of the Planck CMB bispectrum. Finally, we use our measurements to obtain constraints on earlyUniverse scenarios that can generate primordial NG. We consider, for example, general singlefield models of inflation, curvaton models, models with axion fields producing parityviolating tensor bispectra, and inflationary scenarios generating directionallydependent bispectra (such as those involving vector fields).
In our data analysis efforts, which started with the 2013 release, we achieved a number of crucial scientific goals. In particular we reached an unprecedented level of sensitivity in the determination of the bispectrum and trispectrum amplitude parameters (f_{NL}, g_{NL}) and significantly extended the standard local, equilateral, and orthogonal analysis, encompassing a large number of additional shapes motivated by a variety of inflationary models. Moreover, we produced the first polarizationbased CMB bispectrum constraints and the first detection of the (nonprimordial) bispectrum induced by correlations between CMB lensing and secondary anisotropies. Our stringent tests of many types of nonGaussianity are fully consistent with expectations from the standard singlefield slowroll paradigm and provide strong constraints on alternative scenarios. Nevertheless, the current level of sensitivity does not allow us to rule out or confirm most alternative scenarios. It is natural at this stage to ask ourselves what should be the f_{NL} sensitivity goal for future cosmological experiments. A number of studies has identified f_{NL} ∼ 1 as a target. Achieving such sensitivity for localtype NG would enable us to either confirm or rule out a large class of multifield models. A similar target for equilateral, orthogonal, and scaledependent shapes would allow us to distinguish standard slowroll from more complex singlefield scenarios, such as those characterized by higherderivative kinetic terms or slowrollbreaking features in the inflaton potential (see e.g., Alvarez et al. 2014; Finelli et al. 2018, and references therein). With this aim in mind, the challenge for future cosmological observations will be therefore that of reducing the f_{NL} errors from this paper by at least one order of magnitude.
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 and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).
See, e.g. Byrnes & Choi (2010) for a review on this type of model in the context of primordial NG. Early papers discussing primordial local bispectra given by Eq. (3) include Falk et al. (1993), Gangui et al. (1994), Gangui & Martin (2000), Verde et al. (2000), Wang & Kamionkowski (2000), and Komatsu & Spergel (2001).
This would help in further breaking (via primordial NG) some degeneracies among the parameters determining the curvature power spectrum in these modes. For a discussion and an analysis of this type, see Planck Collaboration XXII (2014) and Planck Collaboration XX (2016).
Notice that indeed these models generate bispectra (and power spectra) that break statistical isotropy and, after an angle average, the bispectrum takes the above expression (Eq. (20)).
As a reminder, let us stress that the expression ‘lensing bispectrum’ in this paper always refers to the 3point function generated by correlations between the lensing potential and ISW or reionization contributions, as explained in the main text. We are therefore not referring here to NG lensing signatures arising from the deflection potential alone, such as those considered in the context of CMB lensing reconstruction (and producing a leading trispectrum contribution).
While a onetailed pvalue quantifies how likely it is for the maxima (a similar description holds for the minima) of the Monte Carlo simulations to be higher than the maximum of the data, the twotailed pvalue also considers how likely it is for it to be lower than that of the data. This means that a maximum that is too low yields a low pvalue too.
While no dust is detected in the cleaned maps using the template determined from the dust map, it should be pointed out that there is no guarantee that the dust residuals (or negative dust residuals in the case of an oversubtraction), after passing through the componentseparation pipelines, have exactly the same form as the original dust bispectrum. However, it seems reasonable to assume that the resulting shape would still be highly correlated with the original dust template, so that this remains a meaningful test.
The intermediate and weak dissipative regimes (r_{d} ≤ 1) predict an NG amplitude with a strong dependence on the microscopic parameters (T/H and r_{d}), giving rise to a different additional bispectrum shape (see BasteroGil et al. 2014).
There might also be the case where the intrinsic NG is vanishing and primordial NG is generated only by nonlinearities in the conversion process, reaching an amplitude (Qiu et al. 2013; Li 2013; Fertig et al. 2014).
Acknowledgments
The Planck Collaboration acknowledges the support of: ESA; CNES and CNRS/INSUIN2P3INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, 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 ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck. Some of the results in this paper have been derived using the HEALPix package. We also acknowledge the use of CAMB. Part of this work was undertaken on the STFC COSMOS@DiRAC HPC Facility at the University of Cambridge, funded by UK BIS NEI grants. We gratefully acknowledge the IN2P3 Computer Center (http://cc.in2p3.fr) for providing a significant amount of the computing resources and services needed for the analysis. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DEAC0205CH11231. Some computations were performed on the GPC and Niagara cluster at the SciNet HPC Consortium; SciNet is funded by the Canada Foundation for Innovation under the auspices of Compute Canada, the Government of Ontario, and the University of Toronto.
References
 Achúcarro, A., Gong, J.O., Hardeman, S., Palma, G. A., & Patil, S. P. 2011, JCAP, 1, 30 [Google Scholar]
 Ackerman, L., Carroll, S. M., & Wise, M. B. 2007, Phys. Rev. D, 75, 083502 [NASA ADS] [CrossRef] [Google Scholar]
 Adshead, P., Dvorkin, C., Hu, W., & Lim, E. A. 2012, Phys. Rev. D, 85, 023531 [NASA ADS] [CrossRef] [Google Scholar]
 Agrawal, A., Fujita, T., & Komatsu, E. 2018, Phys. Rev. D, 97, 103526 [CrossRef] [Google Scholar]
 Agullo, I., & Parker, L. 2011, Gen. Rel. Grav., 43, 2541 [NASA ADS] [CrossRef] [Google Scholar]
 Akita, Y., & Kobayashi, T. 2016, Phys. Rev. D, 93, 043519 [CrossRef] [Google Scholar]
 Alishahiha, M., Silverstein, E., & Tong, D. 2004, Phys. Rev. D, 70, 123505 [NASA ADS] [CrossRef] [Google Scholar]
 Allen, M. 1997, The Coefficient of Determination in Multiple Regression (Boston, MA: Springer, US), 91 [Google Scholar]
 Alvarez, M., Baldauf, T., Bond, J. R., et al. 2014, ArXiv eprints [arXiv:1412.4671] [Google Scholar]
 Anninos, D., De Luca, V., Franciolini, G., Kehagias, A., & Riotto, A. 2019, JCAP, 2019, 045 [CrossRef] [Google Scholar]
 ArkaniHamed, N., & Maldacena, J. 2015, ArXiv eprints [arXiv:1503.08043] [Google Scholar]
 ArkaniHamed, N., Creminelli, P., Mukohyama, S., & Zaldarriaga, M. 2004, JCAP, 4, 1 [Google Scholar]
 ArkaniHamed, N., Baumann, D., Lee, H., & Pimentel, G. L. 2018, ArXiv eprints [arXiv:1811.00024] [Google Scholar]
 Arroja, F., Mizuno, S., Koyama, K., & Tanaka, T. 2009, Phys. Rev. D, 80, 043527 [NASA ADS] [CrossRef] [Google Scholar]
 Babich, D., & Zaldarriaga, M. 2004, Phys. Rev. D, 70, 083005 [NASA ADS] [CrossRef] [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., Namba, R., & Peloso, M. 2012a, Phys. Rev. D, 85, 123523 [NASA ADS] [CrossRef] [Google Scholar]
 Barnaby, N., Moxon, J., Namba, R., et al. 2012b, Phys. Rev. D, 86, 103508 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., & Orlando, G. 2017, JCAP, 2017, 034 [CrossRef] [Google Scholar]
 Bartolo, N., & Riotto, A. 2009, JCAP, 3, 17 [NASA ADS] [CrossRef] [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, JCAP, 1, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2004c, Phys. Rev. D, 69, 043503 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2004d, Phys. Rev. Lett., 93, 231301 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2005, JCAP, 8, 10 [Google Scholar]
 Bartolo, N., Fasiello, M., Matarrese, S., & Riotto, A. 2010a, JCAP, 8, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Fasiello, M., Matarrese, S., & Riotto, A. 2010b, JCAP, 9, 35 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2010c, Adv. Astron., 2010, 157079 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Fasiello, M., Matarrese, S., & Riotto, A. 2010d, JCAP, 12, 026 [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., & Riotto, A. 2012, JCAP, 2, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., Peloso, M., & Ricciardone, A. 2013a, Phys. Rev. D, 87, 023504 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., Peloso, M., & Ricciardone, A. 2013b, JCAP, 8, 22 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Peloso, M., Ricciardone, A., & Unal, C. 2014, JCAP, 11, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Bartolo, N., Matarrese, S., Peloso, M., & Shiraishi, M. 2015, JCAP, 7, 039 [CrossRef] [Google Scholar]
 Bartolo, N., Orlando, G., & Shiraishi, M. 2019, JCAP, 2019, 050 [CrossRef] [Google Scholar]
 BasteroGil, M., Berera, A., Moss, I. G., & Ramos, R. O. 2014, JCAP, 12, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Battefeld, T., Niemeyer, J. C., & Vlaykov, D. 2013, JCAP, 5, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Baumann, D., & Green, D. 2012, Phys. Rev. D, 85, 103520 [NASA ADS] [CrossRef] [Google Scholar]
 Baumann, D., Goon, G., Lee, H., & Pimentel, G. L. 2018, J. High Energy Phys., 4, 140 [CrossRef] [Google Scholar]
 Becker, A., & Huterer, D. 2012, Phys. Rev. Lett., 109, 121302 [CrossRef] [Google Scholar]
 Becker, A., Huterer, D., & Kadota, K. 2011, JCAP, 1, 006 [CrossRef] [Google Scholar]
 Becker, A., Huterer, D., & Kadota, K. 2012, JCAP, 12, 034 [CrossRef] [Google Scholar]
 Bennett, C. L., Larson, D., Weiland, J. L., et al. 2013, ApJS, 208, 20 [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, JCAP, 8, 29 [NASA ADS] [CrossRef] [Google Scholar]
 Bucher, M., Moodley, K., & Turok, N. 2000, Phys. Rev. D, 62, 083508 [NASA ADS] [CrossRef] [Google Scholar]
 Bucher, M., van Tent, B., & Carvalho, C. S. 2010, MNRAS, 407, 2193 [NASA ADS] [CrossRef] [Google Scholar]
 Bucher, M., Racine, B., & van Tent, B. 2016, JCAP, 5, 055 [CrossRef] [Google Scholar]
 Burrage, C., de Rham, C., Seery, D., & Tolley, A. J. 2011, JCAP, 1, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Byrnes, C. T., & Choi, K.Y. 2010, Adv. Astron., 2010 [NASA ADS] [CrossRef] [Google Scholar]
 Byrnes, C. T., Nurmi, S., Tasinato, G., & Wands, D. 2010a, JCAP, 2, 034 [CrossRef] [Google Scholar]
 Byrnes, C. T., Gerstenlauer, M., Nurmi, S., Tasinato, G., & Wands, D. 2010b, JCAP, 10, 004 [CrossRef] [Google Scholar]
 Chambers, A., & Rajantie, A. 2008, Phys. Rev. Lett., 100, 041302 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X. 2005, Phys. Rev. D, 72, 123518 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X. 2010a, JCAP, 12, 3 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X. 2010b, Adv. Astron., 2010, 638979 [CrossRef] [Google Scholar]
 Chen, X., & Wang, Y. 2010, JCAP, 4, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., Easther, R., & Lim, E. A. 2007a, JCAP, 6, 23 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., Huang, M.X., Kachru, S., & Shiu, G. 2007b, JCAP, 1, 002 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., Easther, R., & Lim, E. A. 2008, JCAP, 0804, 010 [NASA ADS] [CrossRef] [Google Scholar]
 Chen, X., Hu, B., Huang, M.X., Shiu, G., & Wang, Y. 2009, JCAP, 8, 8 [Google Scholar]
 Chen, X., Firouzjahi, H., Namjoo, M. H., & Sasaki, M. 2013, EPL (Europhys. Lett.), 102, 59001 [CrossRef] [Google Scholar]
 Cheung, C., Fitzpatrick, A. L., Kaplan, J., Senatore, L., & Creminelli, P. 2008, J. High Energy Phys., 3, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Cicoli, M., Tasinato, G., Zavala, I., Burgess, C. P., & Quevedo, F. 2012, JCAP, 5, 39 [NASA ADS] [CrossRef] [Google Scholar]
 Cook, J. L., & Sorbo, L. 2012, Phys. Rev. D, 85, 023534 [CrossRef] [Google Scholar]
 Cook, J. L., & Sorbo, L. 2013, JCAP, 11, 47 [NASA ADS] [CrossRef] [Google Scholar]
 Coulton, W. R., & Spergel, D. N. 2019, JCAP, 10, 56 [CrossRef] [Google Scholar]
 Creminelli, P., & Zaldarriaga, M. 2004, Phys. Rev. D, 70, 083532 [NASA ADS] [CrossRef] [Google Scholar]
 Creminelli, P., Nicolis, A., Senatore, L., Tegmark, M., & Zaldarriaga, M. 2006, JCAP, 5, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Creminelli, P., Pitrou, C., & Vernizzi, F. 2011, JCAP, 11, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Domènech, G., Hiramatsu, T., Lin, C., et al. 2017, JCAP, 2017, 034 [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]
 Elsner, F., & Wandelt, B. D. 2009, ApJS, 184, 264 [NASA ADS] [CrossRef] [Google Scholar]
 Elsner, F., & Wandelt, B. D. 2012, ArXiv eprints [arXiv:1211.0585] [Google Scholar]
 Enqvist, K., & Sloth, M. S. 2002, Nucl. Phys. B, 626, 395 [NASA ADS] [CrossRef] [Google Scholar]
 Endlich, S., Nicolis, A., & Wang, J. 2013, JCAP, 10, 11 [NASA ADS] [CrossRef] [Google Scholar]
 Endlich, S., Horn, B., Nicolis, A., & Wang, J. 2014, Phys. Rev. D, 90, 063506 [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]
 Falk, T., Rangarajan, R., & Srednicki, M. 1993, ApJ, 403, L1 [NASA ADS] [CrossRef] [Google Scholar]
 Fergusson, J. 2014, Phys. Rev. D, 90, 043533 [NASA ADS] [CrossRef] [Google Scholar]
 Fergusson, J. R., Liguori, M., & Shellard, E. P. S. 2010, Phys. Rev. D, 82, 023502 [NASA ADS] [CrossRef] [Google Scholar]
 Fergusson, J. R., Liguori, M., & Shellard, E. P. S. 2012, JCAP, 12, 32 [NASA ADS] [CrossRef] [Google Scholar]
 Fergusson, J., Gruetjen, H., Shellard, E., & Liguori, M. 2015, Phys. Rev. D, 91, 023502 [NASA ADS] [CrossRef] [Google Scholar]
 Fertig, A., Lehners, J.L., & Mallwitz, E. 2014, Phys. Rev. D, 89, 103537 [NASA ADS] [CrossRef] [Google Scholar]
 Finelli, F., Bucher, M., Achúcarro, A., et al. 2018, J. Cosmol. AstroPart. Phys., 2018, 016 [NASA ADS] [CrossRef] [Google Scholar]
 Flauger, R., & Pajer, E. 2011, JCAP, 1, 17 [NASA ADS] [CrossRef] [Google Scholar]
 Flauger, R., McAllister, L., Pajer, E., Westphal, A., & Xu, G. 2010, JCAP, 6, 9 [NASA ADS] [CrossRef] [Google Scholar]
 Flauger, R., McAllister, L., Silverstein, E., & Westphal, A. 2017, JCAP, 1710, 055 [CrossRef] [Google Scholar]
 Franciolini, G., Kehagias, A., Riotto, A., & Shiraishi, M. 2018, Phys. Rev. D, 98, 043533 [CrossRef] [Google Scholar]
 Gangui, A., & Martin, J. 2000, 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]
 Gao, X., Kobayashi, T., Yamaguchi, M., & Yokoyama, J. 2011, Phys. Rev. Lett., 107, 211301 [NASA ADS] [CrossRef] [Google Scholar]
 Giannantonio, T., Porciani, C., Carron, J., Amara, A., & Pillepich, A. 2012, MNRAS, 422, 2854 [CrossRef] [Google Scholar]
 Goldberg, D. M., & Spergel, D. N. 1999, Phys. Rev. D, 59, 103002 [NASA ADS] [CrossRef] [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., Meyers, J., & van Engelen, A. 2017, JCAP, 1712, 005 [CrossRef] [Google Scholar]
 Hall, L. M. H., & Peiris, H. V. 2008, JCAP, 1, 27 [NASA ADS] [CrossRef] [Google Scholar]
 Hannestad, S., Haugbolle, T., Jarnhus, P. R., & Sloth, M. S. 2010, JCAP, 1006, 001 [Google Scholar]
 Hanson, D., & Lewis, A. 2009, Phys. Rev. D, 80, 063004 [NASA ADS] [CrossRef] [MathSciNet] [Google Scholar]
 Heavens, A. F. 1998, MNRAS, 299, 805 [NASA ADS] [CrossRef] [Google Scholar]
 Hikage, C., Koyama, K., Matsubara, T., Takahashi, T., & Yamaguchi, M. 2009, MNRAS, 398, 2188 [NASA ADS] [CrossRef] [Google Scholar]
 Hikage, C., Kawasaki, M., Sekiguchi, T., & Takahashi, T. 2013a, JCAP, 7, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Hikage, C., Kawasaki, M., Sekiguchi, T., & Takahashi, T. 2013b, JCAP, 3, 20 [NASA ADS] [CrossRef] [Google Scholar]
 Hill, J. C. 2018, Phys. Rev. D, 98, 083542 [CrossRef] [Google Scholar]
 Holman, R., & Tolley, A. J. 2008, JCAP, 0805, 001 [Google Scholar]
 Hu, W. 2000, Phys. Rev. D, 62, 043007 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, M.X., & Shiu, G. 2006, Phys. Rev. D, 74, 121301 [NASA ADS] [CrossRef] [Google Scholar]
 Huang, Z., & Vernizzi, F. 2013, Phys. Rev. Lett., 110, 101303 [CrossRef] [Google Scholar]
 Jung, G., & van Tent, B. 2017, JCAP, 1705, 019 [CrossRef] [Google Scholar]
 Jung, G., Racine, B., & van Tent, B. 2018, JCAP, 1811, 047 [CrossRef] [Google Scholar]
 Kawakami, E., Kawasaki, M., Miyamoto, K., Nakayama, K., & Sekiguchi, T. 2012, JCAP, 7, 37 [NASA ADS] [CrossRef] [Google Scholar]
 Kawasaki, M., Nakayama, K., Sekiguchi, T., Suyama, T., & Takahashi, F. 2008, JCAP, 11, 19 [NASA ADS] [CrossRef] [Google Scholar]
 Kawasaki, M., Nakayama, K., Sekiguchi, T., Suyama, T., & Takahashi, F. 2009, JCAP, 1, 42 [NASA ADS] [CrossRef] [Google Scholar]
 Khatri, R., & Wandelt, B. D. 2009, Phys. Rev. D, 79, 023501 [NASA ADS] [CrossRef] [Google Scholar]
 Khatri, R., & Wandelt, B. D. 2010, Phys. Rev. D, 81, 103518 [CrossRef] [Google Scholar]
 Khoury, J., & Piazza, F. 2009, JCAP, 7, 026 [CrossRef] [Google Scholar]
 Kofman, L. 2003, ArXiv eprints [arXiv:astroph/0303614] [Google Scholar]
 Kolb, E. W., Riotto, A., & Vallinotto, A. 2006, Phys. Rev. D, 73, 023522 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E. 2002, ArXiv eprints [arXiv:astroph/0206039] [Google Scholar]
 Komatsu, E. 2010, Class. Quant. Grav., 27, 124010 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., & Spergel, D. N. 2001, Phys. Rev. D, 63, 063002 [NASA ADS] [CrossRef] [Google Scholar]
 Komatsu, E., Spergel, D. N., & Wandelt, B. D. 2005, ApJ, 634, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Koyama, K., Mizuno, S., Vernizzi, F., & Wands, D. 2007, JCAP, 0711, 024 [NASA ADS] [CrossRef] [Google Scholar]
 Kumar, J., Leblond, L., & Rajaraman, A. 2010, JCAP, 4, 024 [CrossRef] [Google Scholar]
 Lacasa, F., Pénin, A., & Aghanim, N. 2014, MNRAS, 439, 123 [NASA ADS] [CrossRef] [Google Scholar]
 Lagache, G., Puget, J.L., & Dole, H. 2005, ARA&A, 43, 727 [NASA ADS] [CrossRef] [Google Scholar]
 Langlois, D., & Lepidi, A. 2011, JCAP, 1, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Langlois, D., & van Tent, B. 2011, Class. Quant. Grav., 28, 222001 [NASA ADS] [CrossRef] [Google Scholar]
 Langlois, D., & van Tent, B. 2012, JCAP, 7, 40 [NASA ADS] [CrossRef] [Google Scholar]
 Langlois, D., Vernizzi, F., & Wands, D. 2008, JCAP, 12, 4 [NASA ADS] [CrossRef] [Google Scholar]
 Lehners, J. L. 2010, Adv. Astron., 2010 [NASA ADS] [CrossRef] [Google Scholar]
 Lehners, J.L., & Steinhardt, P. J. 2008, Phys. Rev. D, 77, 063533 [NASA ADS] [CrossRef] [Google Scholar]
 Lehners, J.L., & Steinhardt, P. J. 2013, Phys. Rev. D, 87, 123533 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A. 2012, JCAP, 6, 023 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A., Challinor, A., & Hanson, D. 2011, JCAP, 3, 18 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473 [Google Scholar]
 Li, M. 2013, Phys. Lett. B, 724, 192 [NASA ADS] [CrossRef] [Google Scholar]
 Liguori, M., Hansen, F. K., Komatsu, E., Matarrese, S., & Riotto, A. 2006, Phys. Rev. D, 73, 043505 [CrossRef] [Google Scholar]
 Liguori, M., Sefusatti, E., Fergusson, J. R., & Shellard, E. P. S. 2010, Adv. Astron., 2010 [NASA ADS] [CrossRef] [Google Scholar]
 Linde, A., & Mukhanov, V. 1997, Phys. Rev. D, 56, 535 [NASA ADS] [CrossRef] [Google Scholar]
 Linde, A., & Mukhanov, V. 2006, JCAP, 4, 9 [NASA ADS] [CrossRef] [Google Scholar]
 LoVerde, M., Miller, A., Shandera, S., & Verde, L. 2008, JCAP, 4, 014 [NASA ADS] [CrossRef] [Google Scholar]
 Lyth, D. H. 2005, JCAP, 11, 6 [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., & Riotto, A. 2006, Phys. Rev. Lett., 97, 121301 [NASA ADS] [CrossRef] [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. M., & Pimentel, G. L. 2011, J. High Energy Phys., 9, 45 [NASA ADS] [CrossRef] [Google Scholar]
 Malik, K. A., & Lyth, D. H. 2006, JCAP, 9, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Mangilli, A., & Verde, L. 2009, Phys. Rev. D, 80, 123007 [NASA ADS] [CrossRef] [Google Scholar]
 Mangilli, A., Wandelt, B., Elsner, F., & Liguori, M. 2013, A&A, 555, A82 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Meerburg, P. D., van der Schaar, J. P., & Stefano Corasaniti, P. 2009, JCAP, 5, 18 [Google Scholar]
 Meerburg, P. D., Meyers, J., van Engelen, A., & AliHaïmoud, Y. 2016a, Phys. Rev. D, 93, 123511 [CrossRef] [Google Scholar]
 Meerburg, P. D., Münchmeyer, M., & Wandelt, B. 2016b, Phys. Rev. D, 93, 043536 [NASA ADS] [CrossRef] [Google Scholar]
 Mirbabayi, M., Senatore, L., Silverstein, E., & Zaldarriaga, M. 2015, Phys. Rev. D, 91, 063518 [CrossRef] [Google Scholar]
 Mollerach, S. 1990, Phys. Rev. D, 42, 313 [NASA ADS] [CrossRef] [Google Scholar]
 Moradinezhad Dizgah, A., Lee, H., Muñoz, J. B., & Dvorkin, C. 2018, JCAP, 5, 013 [CrossRef] [Google Scholar]
 Moroi, T., & Takahashi, T. 2001, Phys. Lett. B, 522, 215 [NASA ADS] [CrossRef] [Google Scholar]
 Moss, I. G., & Xiong, C. 2007, JCAP, 4, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Münchmeyer, M., Bouchet, F., Jackson, M., & Wandelt, B. 2014, A&A, 570, A94 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Münchmeyer, M., Meerburg, P. D., & Wandelt, B. D. 2015, Phys. Rev. D, 91, 043534 [NASA ADS] [CrossRef] [Google Scholar]
 Munshi, D., & Heavens, A. 2010, MNRAS, 401, 2406 [NASA ADS] [CrossRef] [Google Scholar]
 Namba, R., Peloso, M., Shiraishi, M., Sorbo, L., & Unal, C. 2016, JCAP, 1, 041 [CrossRef] [Google Scholar]
 Naskar, A., & Pal, S. 2018, Phys. Rev. D, 98, 083520 [CrossRef] [Google Scholar]
 Nitta, D., Komatsu, E., Bartolo, N., Matarrese, S., & Riotto, A. 2009, JCAP, 5, 14 [NASA ADS] [CrossRef] [Google Scholar]
 Noumi, T., Yamaguchi, M., & Yokoyama, D. 2013, J. High Energy Phys., 6, 51 [CrossRef] [Google Scholar]
 Oppizzi, F., Liguori, M., Renzi, A., Arroja, F., & Bartolo, N. 2018, JCAP, 5, 045 [CrossRef] [Google Scholar]
 Ozsoy, O., Mylova, M., Parameswaran, S., et al. 2019, JCAP, 9, 036 [CrossRef] [Google Scholar]
 Peloton, J., Schmittfull, M., Lewis, A., Carron, J., & Zahn, O. 2017, Phys. Rev. D, 95, 043508 [CrossRef] [Google Scholar]
 Pénin, A., Lacasa, F., & Aghanim, N. 2014, MNRAS, 439, 143 [NASA ADS] [CrossRef] [Google Scholar]
 Pettinari, G. W., Fidler, C., Crittenden, R., et al. 2014, Phys. Rev. D, 90, 103010 [CrossRef] [Google Scholar]
 Planck Collaboration XVII. 2014, A&A, 571, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIX. 2014, A&A, 571, A19 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXII. 2014, A&A, 571, A22 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXIV. 2014, A&A, 571, A24 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2016, A&A, 594, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XX. 2016, A&A, 594, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXI. 2016, A&A, 594, A21 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2020, A&A, 641, A1 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2020, A&A, 641, A2 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2020, A&A, 641, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2020, A&A, 641, A4 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2020, A&A, 641, A5 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2020, A&A, 641, A7 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2020, A&A, 641, A8 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2020, A&A, 641, A9 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2020, A&A, 641, A10 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2020, A&A, 641, A11 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2020, A&A, 641, A12 [CrossRef] [EDP Sciences] [Google Scholar]
 Qiu, T., Gao, X., & Saridakis, E. N. 2013, Phys. Rev. D, 88, 043525 [NASA ADS] [CrossRef] [Google Scholar]
 Ratra, B. 1992, ApJ, 391, L1 [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]
 Riotto, A., & Sloth, M. S. 2011, Phys. Rev. D, 83, 041301 [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]
 Sasaki, M., Valiviita, J., & Wands, D. 2006, Phys. Rev. D, 74, 103003 [NASA ADS] [CrossRef] [Google Scholar]
 Sefusatti, E., Liguori, M., Yadav, A. P. S., Jackson, M. G., & Pajer, E. 2009, JCAP, 12, 022 [CrossRef] [Google Scholar]
 Senatore, L., & Zaldarriaga, M. 2012a, JCAP, 1208, 001 [Google Scholar]
 Senatore, L., & Zaldarriaga, M. 2012b, JHEP, 1204, 024 [NASA ADS] [CrossRef] [Google Scholar]
 Senatore, L., Tassev, S., & Zaldarriaga, M. 2009, JCAP, 9, 38 [NASA ADS] [CrossRef] [Google Scholar]
 Senatore, L., Smith, K. M., & Zaldarriaga, M. 2010, JCAP, 1, 28 [NASA ADS] [CrossRef] [Google Scholar]
 Senatore, L., Silverstein, E., & Zaldarriaga, M. 2014, JCAP, 2014, 016 [CrossRef] [Google Scholar]
 Shandera, S., Dalal, N., & Huterer, D. 2011, JCAP, 3, 017 [CrossRef] [Google Scholar]
 Shiraishi, M. 2012, JCAP, 6, 15 [NASA ADS] [CrossRef] [Google Scholar]
 Shiraishi, M., Komatsu, E., Peloso, M., & Barnaby, N. 2013a, JCAP, 5, 2 [NASA ADS] [CrossRef] [Google Scholar]
 Shiraishi, M., Ricciardone, A., & Saga, S. 2013b, JCAP, 11, 51 [NASA ADS] [CrossRef] [Google Scholar]
 Shiraishi, M., Liguori, M., & Fergusson, J. R. 2014, JCAP, 5, 8 [NASA ADS] [CrossRef] [Google Scholar]
 Shiraishi, M., Liguori, M., & Fergusson, J. R. 2015, JCAP, 1, 7 [NASA ADS] [CrossRef] [Google Scholar]
 Shiraishi, M., Liguori, M., Fergusson, J. R., & Shellard, E. P. S. 2019, JCAP, 06, 046 [CrossRef] [Google Scholar]
 Silverstein, E., & Tong, D. 2004, Phys. Rev. D, 70, 103505 [NASA ADS] [CrossRef] [Google Scholar]
 Sitwell, M., & Sigurdson, K. 2014, Phys. Rev. D, 89, 123509 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. M., Senatore, L., & Zaldarriaga, M. 2009, JCAP, 9, 6 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. M., Loverde, M., & Zaldarriaga, M. 2011, Phys. Rev. Lett., 107, 191301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Smith, K. M., Senatore, L., & Zaldarriaga, M. 2015, ArXiv eprints [arXiv:1502.00635] [Google Scholar]
 Su, S.C., Lim, E. A., & Shellard, E. P. S. 2012, ArXiv eprints [arXiv:1212.6968] [Google Scholar]
 Toffolatti, L., Argueso Gomez, F., de Zotti, G., et al. 1998, MNRAS, 297, 117 [NASA ADS] [CrossRef] [Google Scholar]
 Torrado, J., Byrnes, C. T., Hardwick, R. J., Vennin, V., & Wands, D. 2018, Phys. Rev. D, 98, 063525 [CrossRef] [Google Scholar]
 Tzavara, E., & van Tent, B. 2011, JCAP, 6, 26 [NASA ADS] [CrossRef] [Google Scholar]
 Verde, L., Wang, L., Heavens, A. F., & Kamionkowski, M. 2000, MNRAS, 313, 141 [NASA ADS] [CrossRef] [Google Scholar]
 Vernizzi, F., & Wands, D. 2006, JCAP, 5, 19 [NASA ADS] [CrossRef] [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., 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]
 Yadav, A. P. S., & Wandelt, B. D. 2010, Adv. Astron., 2010 [NASA ADS] [CrossRef] [Google Scholar]
 Zaldarriaga, M. 2004, Phys. Rev. D, 69, 043508 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Results for the amplitude of the lensing bispectrum from the SMICA, SEVEM, NILC, and Commander foregroundcleaned CMB maps, for different bispectrum estimators.
Results for the amplitude of the lensing bispectrum from the SMICA noSZ temperature map, combined with the standard SMICA polarization map, for the Binned and Modal 1 bispectrum estimators.
Bias in the three primordial f_{NL} parameters due to the lensing signal for the four componentseparation methods.
Joint estimates of the bispectrum amplitudes of unclustered and clustered point sources in the cleaned Planck temperature maps, determined with the Binned bispectrum estimator.
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, SEVEM, NILC, and Commander foregroundcleaned maps.
Results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW estimator from the SMICA foregroundcleaned map.
Results for f_{NL} for local isocurvature NG determined from the SMICA Planck 2018 map with the Binned bispectrum estimator.
Peak statistics, as defined in Fergusson et al. (2015), for the resonance models, showing the maximum “Raw” peak significance, the “Single” peak significance after accounting for the parameter survey lookelsewhere effect, and the “Multi”peak statistic integrating across all peaks (also accounting for the lookelsewhere correction).
Peak statistics, as defined in Fergusson et al. (2015), for the different feature models, showing the Raw peak maximum significance (for the given Modal 2 survey domain), the corrected significance of this Single maximum peak after accounting for the parameter survey size (the lookelsewhere effect) and the Multipeak statistic which integrates across the adjusted significance of all peaks to determine consistency with Gaussianity.
Constraints on models with equilateraltype NG covering the shapes predicted by the effective field theory of inflation, together with constraints on specific noncanonical inflation models, such as DBI inflation.
Constraints on models with excited initial states (nonBunchDavies models), as well as warm inflation.
Directiondependent NG results for both the L = 1 and L = 2 models from the Modal 2 pipeline.
Results for the tensor nonlinearity parameter obtained from the SMICA, SEVEM, NILC, and Commander temperature and polarization maps.
Comparison between local, equilateral, and orthogonal f_{NL} results, obtained using the four different componentseparation pipelines and the Modal 1 bispectrum estimator.
Fraction of simulations for which the differences between pairs of componentseparation methods are above various levels of σ.
Independent and joint estimates (see the main text for a definition) of the f_{NL} parameters of the primordial local shape and the thermal dust bispectrum in the raw Planck 143 GHz temperature map, determined using the Binned bispectrum estimator.
Independent and joint estimates of the f_{NL} parameters of the indicated shapes, including in particular the dust template, for the cleaned maps produced by the four componentseparation methods, as determined with the Binned bispectrum estimator.
Impact of the thermal SunyaevZeldovich (tSZ) effect on f_{NL} estimators for temperature data. First three columns: results for the f_{NL} parameters of the primordial local, equilateral, and orthogonal shapes, determined by the indicated estimators from the SMICA foregroundcleaned temperature map with additional removal of the tSZ contamination. Last three columns: difference with the result from the normal SMICA map (see Table 5), where the error bar is the standard deviation of the differences in f_{NL} for the Gaussian FFP10 simulations when processed through the two different SMICA foregroundcleaning procedures. Results have been determined using an independent singleshape analysis and are reported with subtraction of the lensing bias; error bars are 68% CL.
Planck 2018 constraints on the trispectrum parameters , , and from different componentseparated maps.
All Figures
Fig. 1. Weights of each polarization configuration going into the total value of f_{NL} for, from left to right, local, equilateral, and orthogonal shapes. Since we impose ℓ_{1} ≤ ℓ_{2} ≤ ℓ_{3}, there is a difference between, TEE for example (smallest ℓ is temperature) and EET (largest ℓ is temperature). 

In the text 
Fig. 2. Weights of each polarization configuration going into the total value of the a,ii mixed f_{NL} parameter for, from left to right, CDM, neutrinodensity, and neutrinovelocity isocurvature, in addition to the adiabatic mode. Since we impose ℓ_{1} ≤ ℓ_{2} ≤ ℓ_{3}, there is a difference between, TEE for example (where the smallest ℓ is temperature) and EET (where the largest ℓ is temperature). 

In the text 
Fig. 3. PDF of the running parameter n_{NG} for the onefield local model. Top: SMICA map. Bottom: Commander map. Blue squares give the marginalized posterior assuming a constant prior, green circles are the posterior assuming a Jeffreys prior, and red triangles are the profiled Likelihood. 

In the text 
Fig. 4. PDF of the running parameter n_{NG} for the twofield local model. Top: SMICA map. Bottom: Commander map. Blue squares give the marginalized posterior assuming a constant prior, green circles are the posterior assuming a Jeffreys prior, and red triangles are the profiled Likelihood. 

In the text 
Fig. 5. PDF of the running parameter n_{NG} for the geometric mean equilateral parametrization. Top: SMICA map. Bottom: Commander map. Blue squares give the marginalized posterior assuming a constant prior, green circles are the posterior assuming a Jeffreys prior, and red triangles are the profiled Likelihood. 

In the text 
Fig. 6. Generalized resonancemodel significance surveyed over the Modal 2 frequency range with the uppermost panels for the constant resonance model (Eq. 12), showing Tonly (left) and T+E (right) results, then below this the equilateral resonance model (Eq. (13)), followed by the flattened model (Eq. (14)) with the bottom panels, representing the local resonance model with a squeezed envelope. These models have been investigated using Planck temperature data to ℓ_{max} = 2000 and polarization data to ℓ_{max} = 1500, with the Planck componentseparation methods SMICA (blue), SEVEM (orange), NILC (green), and Commander (red). These results are broadly consistent with those found previously (PCNG15), with some broad peaks of moderate significance emerging in both the equilateral and flattened resonance models, somewhat enhanced by including polarization data (upper middle and lower middle right panels). 

In the text 
Fig. 7. Generalized feature model significance surveyed over the Modal 2 frequency range 0 < ω < 350, after marginalizing over the phase ϕ. Top panels: results for the constant feature model (Eq. (15)) with T only (left) and T+E (right), middle panels: for the equilateral feature model (Eq. (18)), and bottom panels: for the flattened feature model (Eq. (19)). The same conventions apply as in Fig. 6. These feature model results generally have lower significance than obtained previously (PCNG15), with polarization data not tending to reinforce apparent peaks found using temperature data only. 

In the text 
Fig. 8. Top: significance of singlefield models for the potential feature case with a K^{2} cos ωK scaling dependence (Eq. (16)). Bottom: significance for rapidly varying sound speed with a K sin ωK scaling (Eq. (17)). Left panels: for T only and right panels: for T + E. These results have been marginalized over the envelope parameter α (determined by feature width and height) from α = 0 to αω = 90. There appears to be no evidence for these very specific signatures in this frequency range when polarization data are included. 

In the text 
Fig. 9. Highfrequency estimator results for feature and resonance models. Upper two panels: significance for the constant feature model (Eq. (15)) surveyed over the frequency range 0 < ω < 3000 after marginalizing over the phase ϕ, for T and T + E SMICA maps. Bottom two panels: significance for the constant resonance model (Eq. (12)) surveyed over the frequency range 0 < ω < 1000 after marginalizing over the phase ϕ, for T and T + E SMICA maps. As for the Modal expansion case, the highfrequency results generally have lower significance than obtained previously (PCNG15), with polarization data not tending to reinforce apparent peaks found using temperature data only. 

In the text 
Fig. 10. Signaltonoiseweighted temperature and polarization bispectra obtained from Planck SMICA maps using the Modal reconstruction with n_{max} = 2001 polynomial modes. The bispecra are (clockwise from top left) the temperature TTT for ℓ ≤ 1500, the Emode polarization EEE for ℓ ≤ 1100, the mixed temperature and polarization TEE (with T multipoles in the zdirection), and lastly TTE (with E multipoles in the zdirection). The S/N thresholds are the same between all plots. 

In the text 
Fig. 11. Smoothed binned signaltonoise bispectra ℬ for the Planck 2018 cleaned sky maps. We show slices as a function of ℓ_{1} and ℓ_{2} for a fixed ℓ_{3}bin [518, 548]. From left to right results are shown for the four componentseparation methods SMICA, SEVEM, NILC, and Commander. From top to bottom we show: TTT; TTT cleaned of clustered and unclustered point sources; T2E; TE2; and EEE. The colour range shows signaltonoise from −4 to +4. The light grey regions are where the bispectrum is not defined, either because it lies outside the triangle inequality or because of the cut . 

In the text 
Fig. 12. Similar to Fig. 11, but for the ℓ_{3}bin [1291, 1345]. 

In the text 
Fig. 13. Similar to Fig. 11, but for the ℓ_{3}bin [771,799], and only for TTT cleaned of clustered and unclustered point sources. 

In the text 
Fig. 14. Scatter plots of the 2001 Modal 2 coefficients for each combination of componentseparation methods, labelled with the R^{2} coefficient of determination. The figures on the left are for the temperature modes and those on the right for pure polarization modes (the Modal 2 pipeline reconstructs the component of the EEE bispectrum that is orthogonal to TTT, so this is not exactly the same as the EEE bispectrum of the other estimators). 

In the text 
Fig. 15. Scatter plots of the bispectrum values in each bintriplet (13 020 for TTT, all of which have been multiplied by 10^{16}, and 11 221 for EEE, all of which have been multiplied by 10^{20}) for all combinations of componentseparation methods, with the R^{2} coefficient of determination. The figures on the left are for TTT and those on the right for EEE. 

In the text 
Fig. 16. Effects on f_{NL} scatter, simulationbysimulation, of different levels of noise mismatch, as described in Sect. 6.2. Top: EEEonly results. Bottom: combined T+E results. The magenta line represents the case in which we generate an extranoise component in harmonic space, both in temperature and polarization, with variance equal to the difference between the noise spectra in the data and the simulations. All other lines represent cases in which we generate extranoise maps in pixel space, with the rms per pixel extracted from simulations and rescaled by different factors, as specified in the legend; the label “TQU” further specifies the pixelspace test in which extra noise is added to both temperature and polarization data, unlike in the other three pixelspace cases, in which we include only a polarization extranoise component. 

In the text 
Fig. 17. Evolution of the f_{NL} parameters (solid blue line with data points) and their uncertainties (dotted blue lines) for the three primordial bispectrum templates as a function of the minimum multipole number ℓ_{min} used in the analysis. From left to right the panels show local, equilateral, and orthogonal shape results, while the different rows from top to bottom show results for T only, E only, and full T+E data. To indicate more clearly the evolution of the uncertainties, they are also plotted around the final value of f_{NL} (solid green lines without data points, around the horizontal dashed green line). The results here have been determined with the Binned bispectrum estimator for the SMICA map, assume all shapes to be independent, and the lensing bias has been subtracted. 

In the text 
Fig. 18. Same as Fig. 17, but this time as a function of the maximum multipole number ℓ_{max} used in the analysis. 

In the text 
Fig. 19. 68%, 95%, and 99.7% confidence regions in the parameter space , defined by thresholding χ^{2}, as described in the text. 

In the text 
Fig. 20. 68%, 95%, and 99.7% confidence regions in the singlefield inflation parameter space , obtained from Fig. 19 via the change of variables in Eq. (58). 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.