© ESO 2020

1. Introduction

This paper, one of a set associated with the 2018 release (also known as “PR3”) of data from the Planck¹ mission (Planck Collaboration I 2020), resents the data analysis and constraints on primordial non-Gaussianity (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 higher-order 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, non-Gaussian models of inflation. These models span from the irreducible minimal amount of primordial NG, predicted by standard single-field models of slow-roll 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 well-known 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 re-analyse 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 model-independent 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 non-primordial 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 scale-dependent NG models and other selected bispectrum shapes, including NG in the tensor (primordial gravitational wave) sector; in this section reconstructions and model-independent 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 single-field models of inflation

The parameter space of single-field 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

(1)

while the orthogonal NG is described by

(2)

Here the potential Φ is defined in relation to the comoving curvature perturbation ζ by Φ ≡ (3/5)ζ on superhorizon scales (thus corresponding to Bardeen’s gauge-invariant gravitational potential (Bardeen 1980) during matter domination on superhorizon scales). P_Φ(k) = A/k^{4 − n_s} 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 non-standard kinetic term or more general higher-derivative interactions. In the first case the inflaton Lagrangian is , where X = g^μν∂_μϕ ∂_νϕ, with at most one derivative on ϕ (Chen et al. 2007b). Different higher-derivative interactions of the inflaton field characterize, ghost inflation (Arkani-Hamed 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 self-interactions. The Dirac-Born-Infeld (DBI) models of inflation (Silverstein & Tong 2004; Alishahiha et al. 2004) are a string-theory-motivated 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 single-field inflation. We discuss them using one example in Sect. 8. Nevertheless, we shall also explicitly search for such EFT shapes, analysing their exact non-separable predicted shapes, B^EFT1 and B^EFT2, along with those of DBI, B^DBI, and ghost inflation, B^ghost (Arkani-Hamed et al. 2004).

2.2. Multi-field models

The bispectrum for multi-field models is typically of the local type²

(3)

This usually arises when more scalar fields drive inflation and give rise to the primordial curvature perturbation (“multiple-field 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 super-horizon scales to the curvature perturbations. Non-Gaussianities if present are transferred too. This, along with non-linearities 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 small-scale modes, peaking in the “squeezed” configurations k₁ ≪ k₂ ≈ k₃. This is a consequence of the transfer mechanism taking place on superhorizon scales and thus generating a localized point-by-point 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 curvaton-like 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 non-linear 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., Bastero-Gil et al. 2014).

2.3. Isocurvature non-Gaussianity

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 cross-correlation between isocurvature and adiabatic perturbations (Komatsu 2002; Bartolo et al. 2002; Komatsu et al. 2005; Kawasaki et al. 2008, 2009; Langlois et al. 2008; Hikage et al. 2009, 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 cold-dark-matter (CDM) density, baryon-density, neutrino-density, and neutrino-velocity 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 multi-field 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 adiabatic-isocurvature cross-power 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):

(4)

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 non-Gaussianity

We briefly describe inflationary models that predict a mildly scale-dependent 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 (single-field slow-roll, with canonical action and initial conditions) predicts that both the amplitude and scale dependence are of the order of the slow-roll 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 well-motivated models make different predictions and these can be used to confirm or exclude such models. Measuring the running of the non-Gaussianity 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. Local-type scale-dependent bispectrum

We start by describing models with a local-type mildly scale-dependent 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 slow-roll approximation and the δN formalism, Byrnes et al. (2010b) found a quite general expression for the power spectrum of the primordial potential perturbation:

(5)

where the indexes a, b run over the different scalar fields.

The non-linearity parameter then reads (Byrnes et al. 2010b)

(6)

where the last line is the general result valid for any number of slow-roll 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

(7)

where n_multi, a and n_f, ab are parameters of the models that are proportional to the slow-roll 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 non-Gaussianity 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)

(8)

In this case

(9)

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 self-interactions can uniquely be probed by the running of f_NL.

The second class of models are two-field 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)

(10)

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 scale-dependent bispectra, even in that parameter region. Byrnes et al. (2010b) also computed the running of the trispectra amplitudes τ_NL and g_NL. For general single-source models they showed that n_τNL = 2n_NG, analogous to the well-known consistency relation , providing a useful consistency check (see, e.g., Smith et al. 2011).

2.4.2. Equilateral type scale-dependent bispectrum

General single-field models that can produce large bispectra having a significant correlation with the equilateral template also predict a mild running non-Gaussianity. A typical example is DBI-inflation, 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

(11)

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³. 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 single-source model of Eq. (8) and the local two-field 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 well-motivated inflationary models, including those with an imposed shift symmetry or if there are sharp features in the inflationary potential. Our first Planck temperature-only non-Gaussian 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 large-field 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 logarithmically-spaced 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)

(12)

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₁ + k₂ + k₃ constant slices, with either equilateral or flattened shapes, respectively (see e.g., Chen 2010a), which we take to have the form

(13)

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

(14)

This analysis does not exhaust resonant models associated with a non-Bunch-Davies 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 corner-turning 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. Scale-dependent 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)

(15)

where ω is a frequency determined by the specific feature properties and ϕ is a phase. The wavenumbers k_i are in units of Mpc⁻¹. A more accurate analytic bispectrum solution has been found that includes a damping envelope taking the form (Adshead et al. 2012)

(16)

where K = k₁ + k₂ + k₃ and the envelope function is given by D(αωK) = αω/(K sinh(αωK)). Here, the model-dependent 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

(17)

In order to encompass the widest range of physically-motivated feature models, we modulated the predicted signal (Eq. (15)) with equilateral and flattened shapes, as defined in Eq. (13), they are

(18)

(19)

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 non-Bunch Davies bispectra with a flattened or enfolded shape (Chen et al. 2007a).

2.6. Non-Gaussianity from excited initial states

Inflationary perturbations generated by “excited” initial states (non-Bunch-Davies vacuum states) generically create non-Gaussianities 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₁ + k₂ ≈ k₃ (and cyclic permutations). In the present analysis, we investigate all the non-Bunch-Davies (NBD) models discussed in the previous Planck non-Gaussian 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 power-law k-inflation, with excitations generated at time τ_c, yielding an oscillation period (and cut-off) k_c ≈ (τ_cc_s)⁻¹; the two leading-order shapes for excited canonical single-field inflation labelled B^NBD1-cos and B^NBD2-cos (Agullo & Parker 2011); and a non-oscillatory 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^NBD-sin (Chen 2010a), and similar extensions for the single-field excited models (Agullo & Parker 2011), labelled B^NBD1-sin and B^NBD2-sin, with the former dominated by oscillatory squeezed configurations.

2.7. Directional-dependent NG

The standard local bispectrum in the squeezed limit (k₁ ≪ k₂ ≈ k₃) has an amplitude that is the same for all different angles between the large-scale mode with wavevector k₁ and the small-scale modes parameterized by wavevector k_s = (k₂ − k₃)/2. More generally, we can consider “anisotropic” bispectra, in the sense of an angular dependence on the orientation of the large-scale and small-scale 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 non-zero odd L moments (since they no longer vanish when defined using the non-symmetrized small-scale wavevector). We then parameterize variations of local NG using (Shiraishi et al. 2013a):

(20)

where P_L(μ) is the Legendre polynomial with P₀ = 1,  P₁ = μ, and . For instance, in the L = 1 case the shape is given by

(21)

Bispectra of the directionally dependent class in general peak in the squeezed limit (k₁ ≪ k₂ ≈ k₃), but they feature a non-trivial dependence on the parameter . The local NG template corresponds to c_i = 2f_NLδ_i0. The non-linearity 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₁ + k₂ ≈ k₃, 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 long-lived 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² of a U(1) gauge field A^μ, via the interaction term I²(φ)F², where F_μν = ∂_μA_ν − ∂_νA_μ and the coupling I²(φ)F² can allow for scale invariant vector fluctuations to be generated on superhorizon scales (Barnaby et al. 2012a; Bartolo et al. 2013a)⁴. Primordial magnetic fields sourcing curvature perturbations can also generate a dependence on both μ and μ² (Shiraishi 2012). The I²(φ)F² models predict c₂ = c₀/2, while models where the primordial curvature perturbations are sourced by large-scale magnetic fields produce c₀, c₁, and c₂. The so-called “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₂ ≫ c₀ (Endlich et al. 2013, 2014). Inflationary models that break rotational invariance and parity also generate this kind of NG with the specific prediction c₀ : c₁ : c₂ = 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 non-trivial structure underlying the inflationary model (as in solid inflation).

Recently much attention has been focused on the possibility of testing the presence of higher-spin particles via their imprints on higher-order 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; Arkani-Hamed & Maldacena 2015; Baumann et al. 2018; Arkani-Hamed et al. 2018). In the case of long-lived superhorizon higher-spin (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 non-trivial dependence on the mass and spin of the particles, as well as on wavenumber (Arkani-Hamed & Maldacena 2015; Baumann et al. 2018; Moradinezhad Dizgah et al. 2018).

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

Together with the scalar mode, the NG signature in the tensor-mode sector also provides a powerful probe of inflation. In single-field 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 tensor-mode sector. In these cases a non-vanishing bispectrum of primordial gravitational waves, , arises via the non-linear interaction between the axion and the gauge field. Its magnitude varies depending on the shape of the axion-gauge coupling, and, in the best-case 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 tensor-mode production is a subhorizon event). Its size is therefore quantified by the so-called tensor non-linearity parameter,

(22)

with .

In this paper we constrain by measuring the CMB temperature and E-mode bispectra computed from (for the exact shape of see PCNG15). By virtue of their parity-violating nature, the induced CMB bispectra have non-vanishing signal for not only the even but also the odd ℓ₁ + ℓ₂ + ℓ₃ 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 foreground-cleaned temperature and high-pass filtered E-mode data, where the E-mode information for ℓ <  40 was entirely discarded in order to avoid foreground contamination. This paper updates those limits with additional data, including large-scale E-mode information.

We stress that there are also other theoretical models that yield interesting NG signals in the tensor-tensor-tensor bispectrum, as well as mixed correlations such as tensor-tensor-scalar and tensor-scalar-scalar bispectra. Such NG signals would be testable by using the B-mode polarization, while the analysis of T-mode and E-mode bispectra provides quite limited information, even with the Planck resolution because of large contamination due to the scalar-scalar-scalar 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 scalar-mode bias is negligible. This is one reason why we are now focusing on the axion-model f_NL template (Eq. (22)) and measuring it with the PlanckT- and E-mode data.

3. Estimators and data analysis procedures

3.1. Bispectrum estimators

We give here a short description of the data-analysis 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

(23)

where defines the “reduced bispectrum,” and is the Gaunt integral, tahis is the integral over solid angle of the product of three spherical harmonics,

(24)

The Gaunt integral (which can be expressed as a product of Wigner 3j-symbols) 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, {ℓ₁, ℓ₂, ℓ₃}.

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 3-point 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)

(25)

where the normalization N is fixed by requiring unit response to when f_NL = 1. C⁻¹ is the inverse of the block matrix:

(26)

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 mean-field 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 inverse-covariance 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)

(27)

Here, represents the inverse of the following 2 × 2 matrix:

(28)

As already described in PCNG13, we find that this simplification, while avoiding the covariance-inversion operation, still leads to uncertainties that are very close to optimal, provided that the multipoles are pre-filtered using a simple diffusive inpainting method. As in previous analyses, we stick to this approach here.

A brute-force 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 cross-validation 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, model-independent reconstruction of the bispectrum in different decomposition domains, precise characterization of spurious bispectrum components, monitoring direction-dependent NG signals, and so on.

3.1.1. KSW and skew-C_ℓ estimators

Komatsu-Spergel-Wandelt (KSW) and skew-C_ℓ 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 three-dimensional summation over ℓ₁, ℓ₂, ℓ₃ into a product of three separate one-dimensional sums over each multipole.

The skew-C_ℓ pipeline differs from KSW essentially in that, before collapsing the estimate into the f_NL parameter, it initially determines the so called “bispectrum-related power spectrum” (in short, “skew-C_ℓ”) function (see Munshi & Heavens 2010) for details). The slope of this function is shape-dependent, which makes the skew-C_ℓ extension very useful to separate and monitor multiple and spurious NG components in the map.

3.1.2. Running of primordial non-Gaussianity

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 non-Gaussianity, encoded in the spectral index of the non-linear 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:

(29)

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

(30)

where k_sum = k₁ + k₂ + k₃.

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):

(31)

Making one of these substitutions immediately yields the scale-dependent 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 scale-dependent 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 non-flat 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)

(32)

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

(33)

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 scale-dependent bispectra, the Fisher matrix is

(34)

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 Wigner-3j 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

(35)

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

(36)

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 scale-dependent 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 non-separable bispectra are expanded with arbitrary precision into separable basis modes. Therefore the treatment of non-separable 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 model-independent 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 KSW-like 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 non-parametric 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. High-frequency feature and resonant-model 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 high-frequency 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 component-separation 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 mean-field correction terms for our estimators, as well as for performing data-validation 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 component-separation 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 non-masked 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 component-separation 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 best-fit cosmology and making use of the CAMB computer code⁵ (Lewis et al. 2000) to compute radiation transfer functions and theoretical power spectra. The HEALPix computer code⁶ (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 high-pass 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 non-Gaussianity 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 non-standard 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 higher-resolution basis, including 2000 polynomial modes and a Sachs-Wolfe local template, to improve efficiency in the squeezed limit.

For oscillating non-Gaussianities 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 high-frequency range of shapes with the low frequency range used for cross-valiation with Modal 2. Both of these estimators are equivalent to those used in PCNG15.

4. Non-primordial contributions to the CMB bispectrum

In this section we investigate those non-primordial 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 non-primordial contaminants (that we do not detect in the cleaned maps) is part of the validation work in Sect. 6.

4.1. Non-Gaussianity 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 Sachs-Wolfe (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 small-scale 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 high-precision experiment like Planck the effect is large enough that it must be taken into account. The temperature-only 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 T-only and for the full T+E results. The template of the lensing bispectrum⁷ is given by (Hu 2000; Lewis et al. 2011),

(37)

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/polarization-lensing potential cross-power spectra. The functions are defined by

(38)

if the sum ℓ₁ + ℓ₂ + ℓ₃ is even and ℓ₁, ℓ₂, ℓ₃ 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 T-only 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 well-suited 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.

It was pointed out in Hill (2018) that the coupling between the ISW effect and the thermal Sunyaev-Zeldovich (tSZ) effect can produce a significant ISW-tSZ-tSZ temperature bispectrum, which peaks for squeezed modes, and can therefore contaminate especially our local and lensing results. The semi-analytic approach in Hill (2018) makes predictions for single frequency channels, and cannot be directly applied to the multi-frequency component-separated data we are using. However, it is interesting that such an approach shows an anti-correlation at all frequencies between the ISW-tSZ-tSZ contamination and the lensing bispectrum shape. This could be a possible explanation for the slightly low T-only value of observed in the final data, across all component-separated 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 “no-SZ” 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.

The hypothesis that ISW-tSZ-tSZ residuals are contributing to the temperature-only 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 SZ-removed (hereafter “no-SZ”) SMICA measurements of local f_NL, however, do not support this hypothesis, since they do not display any large shift, whereas a residual ISW-tSZ-tSZ bispectrum should correlate with all shapes that peak in the squeezed limit. Still, one cannot exclude the possibility that multi-frequency component-separation affects the two cases differently. Both our local no-SZ 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 ISW-tSZ-tSZ contamination on the temperature-only 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 frequency-by-frequency 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 component-separated 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 T-only data and for T+E the bias is very significant for local and to a lesser extent for orthogonal NG. For T-only 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 E-only the effect is completely negligible.

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 four-point 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 non-Gaussianity constraints beyond Planck. Furthermore, unlike the signal contamination discussed above, this will likely be equally important for polarization. Because of the shape of the lensing-induced 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 flat-sky 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 signal-induced lensing contributions discussed above.

4.2. Non-Gaussianity 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 late-type infrared galaxies (see e.g. Toffolatti et al. 1998; González-Nuevo et al. 2005), while the latter are primarily dusty star-forming 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)

(39)

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 E-only 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)

(40)

where the index is q = 0.85, the break is located at ℓ_break = 70, and ℓ₀ = 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 component-separation 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 component-separated 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 point-source 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.

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 component-separation 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 T-only, E-only, 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 best-fit Λ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 T-only 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.

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 non-separable shapes or reconstruct the full bispectrum, can treat exactly the local, equilateral, and orthogonal templates that are analysed here.

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, full-sky 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 lensing-bias 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 component-separated 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 T-only 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 already-noted 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. detector-subset) 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 model-independent 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 T-only 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 component-separation 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 E-only 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 T-only: 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 E-only analysis has not translated into a corresponding improvement in the T+E analysis. The answer is relatively simple: the EEE-bispectra 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 T-only to the full T+E analysis than the errors for the latter.

In conclusion, our current results show no evidence for non-Gaussianity 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 non-Gaussian measurements are concerned.

5.2. Further bispectrum shapes

5.2.1. Isocurvature non-Gaussianity

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 CDM-density, neutrino-density, and neutrino-velocity 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⁸. The second is that certain parameters depend more on the high-ℓ adiabatic modes (which are well-determined), 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).

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 T-only 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 neutrino-velocity modes improve by a factor of about 2 when going from the T-only to the full T+E analysis (e.g. neutrino velocity i,ii changes from 270 ± 230 to −51 ± 110), while the improvements for the neutrino-density 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 neutrino-density 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.

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 T-only 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 E-only 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 T-only errors. For the neutrino-velocity 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 cross-power 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 T-only 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 neutrino-velocity 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 E-only data). Again there is no evidence for any isocurvature NG: we do not consider the almost 2σ result for the neutrino-density 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 non-Gaussianity

In this section we present our analysis of the scale-dependent 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 one-field local model (Eq. (8)), the two-field local model (Eq. (10)), and the geometric-mean equilateral model, where is parametrized as in Eq. (31). Here we present results for both the SMICA and Commander temperature maps. Since the E-mode 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 scale-dependent 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.

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 KSW-type estimator for the broad class of resonance-type 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 cross-sectional 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 component-separation 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 cross-sections approximately 1.7 times larger and those for the local cross section 0.5 times smaller. The results are consistent across component-separation 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 “look-elsewhere” 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 look-elsewhere-corrected result being for the equilateral cross-section 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.

5.2.4. Scale-dependent 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 cross-sectional 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 cross-sections 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. Look-elsewhere-adjusted results are presented in Table 10. Here all results are consistent with Gaussianity for all models in this class.

5.2.5. High-frequency feature and resonant-model estimator

We used the high-frequency 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 component-separations 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)σ⁹.

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 non-Gaussianity in the high-frequency feature and resonance-model 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. Equilateral-type models and the effective field theory of inflation

Many physically well-motivated inflationary models produce non-Gaussianity 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.

5.2.7. Models with excited initial states (non-Bunch-Davies 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 look-elsewhere 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.

5.2.8. Direction-dependent primordial non-Gaussianity

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 component-separation methods so cannot be considered robust. We can then conclude that our results are consistent with Gaussianity.

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

In this section, we report constraints on the tensor non-linearity parameter (Eq. (22)) obtained from temperature and E-mode polarization maps. As in our 2015 analysis, we examine even and odd ℓ₁ + ℓ₂ + ℓ₃ multipole domains, employing the original parity-even Modal estimator, (Fergusson et al. 2010, 2012; Shiraishi et al. 2019), and its parity-odd version, (Shiraishi et al. 2014, 2015, 2019), respectively. The constraints obtained from both domains are also combined by computing

(41)

where F^even/odd is the Fisher matrix from ℓ₁ + ℓ₂ + ℓ₃ = 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 E-mode 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 component-separated maps are summarized in Table 14. We confirm there that the sizes of errors in the parity-odd T+E analysis reduce to almost the same level as the parity-even counterparts. This is because the low-ℓ signal of TTE, which dominates the signal-to-noise ratio from ℓ₁ + ℓ₂ + ℓ₃ = odd (Shiraishi et al. 2013b), is now taken into account. Although the errors for the parity-even case are as large as the 2015 ones, owing to the sensitivity improvement of the parity-odd part, the whole-domain constraints become more stringent.

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 foreground-removed 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 cross-correlations 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 signal-dominated regime.

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 non-Gaussianity beyond that of the explicit models tested can be found in the cleaned CMB maps. Since we want to detect specific bin-triplets standing out from the noise, we work with the linear-term-corrected signal-to-noise-ratio bispectrum. Except for the high-ℓ bins, where a point-source signal is present, the bin-triplets 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 bin-triplets that have non-Gaussian noise (due to the fact that they contain very few valid ℓ-triplets). We also have to take into account edge effects from the non-trivial 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 signal-to-noise bispectrum ℬ_i1i2i3, with a Gaussian smoothing of σ_bin  =  2, are shown in Figs. 11 and 12. These are slices for the 20th and 40th ℓ₃-bin as a function of ℓ₁ and ℓ₂. For the cross-bispectra 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 component-separation methods SMICA, SEVEM, NILC, and Commander, and for TTT, T2E, TE2, and EEE. We also show the TTT bispectra after we remove the best joint-fit unclustered and clustered point-source contribution (see Table 4). In the top row of Fig. 12, we can clearly see this significant point-source 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 joint-fit 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 component-separation 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.

In order to quantify the possible residual non-Gaussianity in these smoothed bispectra, we focus on the minimum and maximum of our bin-triplets. 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 non-trivial 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 signal-to-noise bispectrum, and compute the p-values 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⁶ Monte Carlo simulations turn out to be sufficient and so we do not need to use the semi-analytical Ansatz introduced in Bucher et al. (2016).

In Table 15, we report for the smoothing lengths σ_bin = 1, 2, and 3, the two-tailed p-value¹⁰ 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.

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 well-verified agreement on both data and simulations.

6.1. Dependence on foreground-cleaning method

6.1.1. Comparison between f_NL measurements

For this test, we consider 160 FFP10-based simulations with realistic beams and noise levels. A starting set of single-frequency Gaussian maps is processed through the different component-separation 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 component-separation 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 component-separation 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.

This leads us to a first interesting observation, namely that polarization-only 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 polarization-based 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 T-only results, even though limited to specific methods and the orthogonal shape only. The most striking example is a large difference between the orthogonal T-only f_NL measurements obtained with SMICA and Commander (Δf_NL/σ_Δ ≈ 10). Smaller but non-negligible are also the orthogonal T-only discrepancies for the NILC – Commander and SMICA – SEVEM pairs (both with Δf_NL/σ_Δ ≈ 4). These results have been cross-checked 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 T-only discrepancy. This is good news, since the main component-separation 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 mode-by-mode or bin-by-bin correlation analysis over the full bispectrum domain (see next section) show no other signs of anomalies for any component-separation 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 detector-set 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 component-separation 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 component-separation 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 component-separation pipelines in a more general, model-independent fashion. For this purpose we calculate the coefficient of determination (see Allen 1997), denoted by R², from the mode or bin amplitudes extracted from different foreground-cleaned 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.

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 component-separation 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 simulation-based 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 non-Gaussianity 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 percent-level 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 signal-dominated regime we consider for our analysis.

As a first test, we generate “extra-noise” multipoles, drawing them from a zero-mean 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 extra-noise mock data set. However, we calibrate our linear term, estimator normalization, and weights using the original noise maps. We check the effect of the mis-calibration 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 map-by-map f_NL differences, while the latter shows the map-by-map 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.

In the test we have just described we generate uncorrelated extra-noise multipoles (with a non-white 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 extra-noise 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 “extra-noise" 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 non-flat 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 hit-count 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 linear-term 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 non-primordial contributions to the Planck bispectrum, namely Galactic dust and the Sunyaev-Zeldovich 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 (second-order) non-linearities from gravity in general relativity and by non-linearities 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). Non-primordial contributions that we do detect (lensing and extragalactic point sources) were discussed in Sect. 4.

6.3.1. Non-Gaussianity of the thermal dust emission

In principle, the foreground-cleaned 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 component-separation 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 143-GHz 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.

The analysis of the T-only maps produced by the four component-separation 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)¹¹. Given the small size of the errors, this non-detection 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%).

6.3.2. Impact of the tSZ effect

The SMICA component-separation method also produces a foreground-cleaned temperature map, “SMICA no-SZ,” where, in addition to the usual foregrounds, contamination by the thermal Sunyaev-Zeldovich (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 T-only 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.

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 foreground-corrected 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 T-only 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 component-separation 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 component-separation methods regarding the T-only 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.

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 E-only 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.

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.

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 E-only 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 T-only equilateral). We see that the T-only results have stabilized by ℓ = 2000 and the E-only 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 component-separation methods and with all bispectrum estimators, separately considering T-only, E-only, and T+E results. The agreement for E-only results has significantly improved with respect to the previous release.

• The possible exception is provided by the orthogonal T-only 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 component-separation methods, using both the Modal and the Binned approaches. Polarization-only 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 component-separated CMB maps.

• Contamination from the thermal SZ effect on the standard primordial and lensing T-only 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 polarization-only bispectra with respect to the previous release. Whereas in 2015 we cautioned the reader to take polarization-based f_NL estimates as preliminary, we can now state that our T+E-based 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 four-point 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 four-point 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 local-type trispectrum arises if the initial adiabatic curvature ζ is given by the following non-Gaussian model:

(42)

where ζ_G is a Gaussian field. The trispectrum in the local model is given by

(43)

In this equation and throughout this section, a “primed” four-point function denotes the four-point function without its momentum-conserving delta function, that is,

(44)

where “+ disc.” denotes disconnected contributions to the 4-point function. It can be shown that the local-type trispectrum in Eq. (43) is always negligibly small in single-field 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σ)²(∂_jσ)² 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:

(45)

(46)

Here K = ∑_ik_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²/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 scale-invariant power spectrum P_ζ(k) = A_ζ/k³. For the Planck analysis, we slightly modify these trispectra to a power-law spectrum P_ζ(k)∝k^{n_s − 4}, as described in Appendix C of Smith et al. (2015).

To estimate each g_NL parameter from Planck data, we use the “pure-MC” 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 experiment-specific linear operation whose input consists of one or more pixel-space maps, and whose output is a single harmonic-space map, d_ℓm. In the Planck trispectrum analysis, we define the filter as follows. First, we take the single-frequency pixel-space maps, and combine them to obtain a single component-separated pixel-space map, using one of the component-separation algorithms, SMICA, SEVEM, NILC, Commander, or SMICA no-SZ (a variant of the SMICA algorithm that guarantees zero response to Compton-y sources, at the expense of slightly higher noise; see Sect. 6.3.2). Second, we subtract the best-fit 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 pixel-space map out to l_max = 1600, obtaining a harmonic-space map d_ℓm. Finally, we define the filtered map by applying the multiplicative factor:

(47)

This sequence of steps defines a linear operation, whose input is a set of single-frequency pixel-space maps, and whose output is a filtered harmonic-space 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 Planck-specific.

The results of the analysis, for all three trispectrum shapes and five different component-separation algorithms, are presented in Table 24. We do not find evidence for a nonzero primordial trispectrum.

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

(48)

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σ)²(∂_jσ)². In addition, there is a third operator that arises at the same order in the EFT expansion. For completeness, its trispectrum is given by

(49)

However, a Fisher matrix analysis shows that this trispectrum is nearly 100% correlated with the and (∂_iσ)²(∂_jσ)² trispectra. If the parameter is non-zero, then we can absorb it into the “effective” values of the parameters and as

(50)

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 two-component parameter vector

(51)

and let denote the two-component vector of single-g_NL estimates from the SMICA maps (Table 24):

(52)

We also define a two-by-two Fisher matrix F_ij, whose diagonal is given by , where σ_i is the single-g_NL statistical error in Table 24, and whose off-diagonal is , where r is the correlation in Eq. (48). This gives:

(53)

Now, given a set of “theory” g_NL values, represented by a two-vector g_i, we compare to the Planck data by computing the following quantity for the trispectrum:

(54)

In a model-building context where the g_NL quantities g_i depend on model parameters, confidence regions on model parameters can be obtained by appropriately thresholding χ². We give some examples in Sect. 8.

8. Implications for early-Universe 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 single-field models of slow-roll 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 one-dimensional limits on the parameter under consideration.

8.1. General single-field models of inflation

DBI models. Dirac-Born-Infeld (DBI) models of inflation (Silverstein & Tong 2004; Alishahiha et al. 2004) arise from high-energy string-theory constructions and generate a non-linearity 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 non-standard kinetic term of the inflaton field. Notice that we have constrained the exact theoretical (non-separable) 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

(55)

and

(56)

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 single-field slow-roll inflation. Typically the equilateral and orthogonal templates represent an accurate basis to describe the full parameter space of EFT single-field models of inflation, and therefore we used the constraints on and .

As a concrete example, let us consider the Lagrangian of general single-field models of inflation (of the form P(X, φ) models, where X = g^μν∂_μϕ ∂_νϕ) written with the EFT approach:

(57)

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₃ 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

(58)

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 non-interacting (vanishing NG) case corresponds to c_s = 1 and M₃ = 0 (or .

We then proceed as in the two previous analyses (PCNG13; PCNG15). We employ a χ² 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.28, 5.99, and 11.62, respectively (corresponding to a χ² variable with two degrees of freedom). In Fig. 20 we show the corresponding confidence regions in the parameter space. Marginalizing over we find

(59)

and

(60)

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. Multi-field 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 non-linearity parameter (Bartolo et al. 2004c,d)

(61)

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 single-field and multi-field scenarios. Our constraint at 68% CL (see Table 6) implies

(62)

The constraint at 68% CL obtained from temperature and polarization data yields the constraint

(63)

These limits indicate that in such scenarios the curvaton field has a non-negligible energy density when it decays. Meaningful improvements are achieved with respect to previous bounds in PCNG15, namely an almost 20% improvement from T-only data and a 10% improvement when including E-mode 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:

(64)

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

(65)

and

(66)

which does not exclude a contribution of curvaton particles comparable to the one from the classical curvaton field.

8.3. Non-standard inflation models

Directional-dependent NG. Table 13 shows the constraints on directionally-dependent 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 scale-invariant gauge fields are present during inflation. As summarized in Sect. 2 they are also produced from partially massless higher-spin 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² of a gauge field via a term ℒ  =   − I²(ϕ)F²/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 super-horizon vector field perturbations switches on the L = 0 and L = 2 modes in the bispectrum, with non-linearity 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 e-folds (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 vacuum-mode curvature power spectrum, ϵ is a slow-roll 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⁻⁹ <  𝒫 <  3.0 × 10⁻⁹; 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 non-linear parameter is related to the tensor-to-scalar 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⁻⁷ and 2.7 × 10⁻⁹ (95% CL) for r = 10⁻² and 10⁻³, 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)¹². Assuming a constant prior 0 ≤ log₁₀r_d ≤ 4, the SMICA constraints (at 68% CL) from T and from T+E (see Table 12), yield log₁₀r_d ≤ 3.6 and log₁₀r_d ≤ 3.5, respectively, at 95 % CL. The results show that strong-dissipative 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” non-Gaussianity 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 ¹³. 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 <  κ₃ <  5 the constraints on from T only (see Table 6), implies −1.1 <  κ₃ <  0.36 and −0.43 <  κ₃ <  1.0 at 95% CL, for the plus and minus sign in , respectively. The T+E constraints on (Table 6) yield −1.05 <  κ₃ <  0.27 and −0.38 <  κ₃ <  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 <  κ₃ <  0.51 and −0.62 <  κ₃ <  1.5 at 95% CL from T only; and −1.5 <  κ₃ <  0.39 and −0.54 <  κ₃ <  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 single-field inflationary models, using the effective action for the Goldstone boson π (see, e.g., Smith et al. 2015):

(67)

The mass scale M₄ is related to our previously-defined g_NL parameters by:

(68)

Therefore, using the limit from Table 24, we get the following constraint on single-field models:

(69)

Next consider the case of multifield inflation. Here, we consider an action of the more general form:

(70)

where σ is a light field that acquires quantum fluctuations with power spectrum P_σ(k) = H²/(2k³). We assume that σ converts to adiabatic curvature ζ, i.e. ζ = (2A_ζ)^1/2H⁻¹σ. The model parameters Λ_i are related to our previously-defined g_NL parameters by:

(71)

and can be constrained by thresholding the χ²-statistic defined in Eq. (54). For example, to constrain the parameter Λ in the Lorentz invariant model

(72)

we set obtain g_NL-values using Eq. (71), and use Eq. (54) to obtain χ² as a function of Λ. We then threshold at Δχ² = 4 (as appropriate for one degree of freedom), to obtain the following constraint:

(73)

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 small-sound-speed 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:

(74)

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 local-type trispectrum with an amplitude given by (Sasaki et al. 2006)

(75)

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

(76)

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 full-mission CMB temperature and E-mode 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 foreground-cleaning methods (SMICA, SEVEM, NILC, and Commander), studies of robustness under changes in sky coverage and multipole range, and an analysis of the impact of noise-related 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 polarization-only, 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 scale-dependent feature and resonance bispectra, running f_NL models, isocurvature primordial NG, and parity-breaking models. We set tight constraints on all these scenarios, but do not detect any significant signals.

On the other hand, the non-primordial 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 best-fit cosmological parameters, further indicating the absence of significant foreground contamination or spurious systematic effects. We also explicitly checked for the presence of various non-primordial 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 point-source 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 field-theory approach to inflation.

In addition to estimates of bispectrum and trispectrum amplitudes, we produce model-independent reconstructions and analyses of the Planck CMB bispectrum. Finally, we use our measurements to obtain constraints on early-Universe scenarios that can generate primordial NG. We consider, for example, general single-field models of inflation, curvaton models, models with axion fields producing parity-violating tensor bispectra, and inflationary scenarios generating directionally-dependent 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 polarization-based CMB bispectrum constraints and the first detection of the (non-primordial) bispectrum induced by correlations between CMB lensing and secondary anisotropies. Our stringent tests of many types of non-Gaussianity are fully consistent with expectations from the standard single-field slow-roll 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 local-type NG would enable us to either confirm or rule out a large class of multi-field models. A similar target for equilateral, orthogonal, and scale-dependent shapes would allow us to distinguish standard slow-roll from more complex single-field scenarios, such as those characterized by higher-derivative kinetic terms or slow-roll-breaking 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.

---

Acknowledgments

The Planck Collaboration acknowledges the support of: ESA; CNES and CNRS/INSU-IN2P3-INP (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. DE-AC02-05CH11231. 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

All Tables

All Figures