EDP Sciences
Planck 2015 results
Free Access
Issue
A&A
Volume 594, October 2016
Planck 2015 results
Article Number A17
Number of page(s) 66
Section Cosmology (including clusters of galaxies)
DOI https://doi.org/10.1051/0004-6361/201525836
Published online 20 September 2016

© ESO, 2016

1. Introduction

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

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

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

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

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

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

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

2. Models

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

2.1. General single-field models of inflation

This class of models includes inflationary models with a non-standard kinetic term (or more general higher-derivative interactions), in which the inflaton fluctuations propagate with an effective sound speed cs which can be smaller than the speed of light. For example, models with a non-standard kinetic term are described by an inflaton Lagrangian ℒ = P(X,φ), where X = gμνμφ,νφ, with at most one derivative on φ, and the sound speed is given by .

The NG parameter space of this class of models is generically well described by two NG shapes – “equilateral” and “orthogonal” (Senatore et al. 2010) – since usually there are two dominant interaction terms of the inflaton field giving rise to the overall NG signal. One of these typically produces a bispectrum very close to the equilateral type with in the limit cs ≪ 1 (Chen et al. 2007b; Senatore et al. 2010).

The equilateral-type NG is well approximated by the template (Creminelli et al. 2006) (3)Here PΦ(k) = A/k4 − ns represents Bardeen’s gravitational potential power spectrum, A2 being the normalization and ns the scalar spectral index. DBI inflationary models based on string theory (Silverstein & Tong 2004; Alishahiha et al. 2004) provide physically well-motivated examples of the P(X,φ)-model. They are characterized by an almost equilateral NG with for cs ≪ 1, which typically is .

The “orthogonal” shape template is (Senatore et al. 2010) (4)Equilateral and orthogonal shapes emerge also in models characterized by more general higher-derivative interactions, such as ghost inflation (Arkani-Hamed et al. 2004), effective field theories of inflation (Cheung et al. 2008; Senatore et al. 2010; Bartolo et al. 2010a), or the so “Galileon-like” models of inflation (see, e.g., Burrage et al. 2011). The latter model is constructed starting from some specific underlying symmetry for the inflaton field, and is characterized by strongly constrained derivative interactions.

2.2. Multi-field models

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

2.3. Isocurvature non-Gaussianity

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

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

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

The primordial shape templates are a generalization of Eq. (5), see Langlois & van Tent (2011, 2012): (6)where I,J,K label the different modes (adiabatic and isocurvature). The invariance under the simultaneous exchange of two of these indices and the corresponding momenta means that , hence reducing the number of independent parameters from eight to six, in the case of two modes. The different bispectra vary most importantly through the fact that different types of radiation transfer functions are used to project the primordial template onto the CMB: the reduced bispectra are of the form (7)with Here j is the spherical bessel function and we use the notation (123) ≡ [ 123 + 5perm. ] / 3 !. In addition to the isocurvature index, each transfer function carries a polarization index that we do not show here. It is important to note that, unlike the case of the purely adiabatic mode, the inclusion of polarization improves the constraints on the isocurvature NG significantly, as predicted by Langlois & van Tent (2011, 2012).

2.4. Resonance and axion monodromy models

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

Simple resonance model: periodic features in the inflationary potential can induce oscillations with frequency ω that can resonate through any interactions with the inflationary fluctuations, contributing to the bispectrum. Provided that ω>H, this mode starts inside the horizon but its frequency decreases as it is stretched by inflation, until frozen when ωH. Thus periodic features introduce a driving force which can scan across a wide range of frequencies. The simplest basic behaviour of such resonant models yields logarithmic stretching and can be described by the non-scale-invariant shape (see Chen et al. 2008; Chen 2010b) (10)where the constant C = 1 / ln(3k), k is a wavenumber associated with the periodicity, and φ is a phase. These oscillations constructively and destructively interfere with the oscillations created by the CMB transfer functions, introducing additional nodal points in the CMB bispectrum.

Generalized resonance models: in a more general context, it is possible to have more complicated resonant shapes and envelopes. Resonant single-field models with varying sound speed cs generate three leading-order bispectrum terms (Chen 2010a): (11)The first term on the right-hand side of Eq. (11) is the basic resonant shape given in Eq. (10), while the second and third terms have the same oscillatory behaviour, but modulated by a (mildly) flattened shape, and an equilateral shape respectively. The third term is in fact the second generic shape arising in effective field theory and correlates well with the equilateral shape in Eq. (12). The second term in Eq. (11) weakly favours flattened triangles, but there are regimes for resonant models that can generate much stronger flat shapes. If the resonance begins very deep inside the horizon, then the second (negative energy) mode can also make a significant contribution that is associated with enfolded or flat bispectra; this is similar to having an excited initial state or non-Bunch-Davies (NBD) vacuum.

With these two physical motivations in mind we also investigate classes of models with resonant oscillations modulated by both the equilateral and flattened shapes, defined by (12)where (here, for simplicity we ignore the spectral index dependence of the equilateral shape in Eq. (3)). The corresponding equilateral and flattened resonant bispectra ansätze are then We note that typically non-Bunch-Davies bispectra can be much more sharply peaked in the flattened or squeezed limits than Eq. (14), but our purpose here is to determine if this type of resonant model is favoured by the Planck data; that is, whether Eq. (14) warrants further investigation with other flattened profiles.

2.5. Scale-dependent oscillatory feature models

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

Constant feature model: In the previous investigation of Planck data using a coarse parameter grid (Planck Collaboration XXIV 2014), we searched for the simplest ansatz for an oscillatory bispectrum signal (Chen et al. 2007a): (15)where φ is a phase factor and ω is a frequency associated with the specific shape of the feature in the potential that disrupts the slow-roll evolution. In the earlier analysis, we also considered a damping envelope, which slightly increased the apparent significance of the best-fit feature models, though at the cost of an additional parameter (see single-field solutions below).

Generalized feature models: here, we again search for oscillatory signals in a model-independent manner. We will modulate the bispectrum cross-sections with the physically motivated equilateral and flattened shapes, reflecting the physical contexts in which they could have been generated, as for the resonant models described above in Eq. (11). If there are potential features in a model with a varying sound speed, then we can expect there to be oscillatory contributions to the bispectrum signal with a dominant equilateral shape. Motivated by the equilateral resonance model in Eq. (11), we will search for the following equilateral feature ansatz: For extremely sharp features, it is possible to excite the inflationary fluctuations as if there were a non-Bunch Davies vacuum: the oscillatory signal becomes modulated with a flattened shape (Chen et al. 2007a). Again, motivated by the enfolded resonance model in Eq. (14), we take the following simple flattened ansatz: (18)Although the exact profile of the flattened shape can be much more highly peaked on the faces in these NBD models, this ansatz should be adequate for testing whether these models are favoured. We note that while the power spectrum is insensitive to the underlying scenario creating the features, the bispectrum shape will reveal whether features arise from varying sound speed or highly excited features in the potential.

Single field feature solutions: here we use the full analytic bispectrum solution given by Adshead et al. (2012), but the dominant leading-order behaviour takes the form (19)where K = k1 + k2 + k3 and D(αωK) = αω/ (Ksinh(αωK)) is an envelope function, with parameter α setting an overall cut-off for the bispectrum at large wavenumbers or multipoles. This envelope and the overall K2 scaling distinguishes this realistic case from the simple separable constant feature ansatz of Eq. (15). We shall allow the envelope parameter α to vary from α = 0, with no envelope (the infinitely thin limit for a feature in the potential) through to large α, with a narrow domain for the bispectrum. Alternative analytic solutions where the bispectrum is created by a variation in the sound speed cs are dominated by the Ksin(ωK) term, as in (20)For the simplest models there is a predicted relationship between the power spectrum and bispectrum amplitude (e.g., see also Achucarro et al. 2013 for a two-field model). We note that typically the power spectrum has larger signal-to-noise at low frequency (i.e., below ω ≃ 1000) while the bispectrum dominates at higher frequency.

2.6. Non-Gaussianity from excited initial states

It is well known that if the initial vacuum state for inflation is excited and deviates from the standard Bunch-Davies vacuum, then measurable non-Gaussianities can be produced (Chen et al. 2007b; Holman & Tolley 2008; Meerburg et al. 2009; Ashoorioon & Shiu 2011). These models generically lead to non-Gaussianity that peaks in the flattened limit, where k1 + k2k3, and also often has oscillatory behaviour. Here we constrain the same selection of templates found in the 2013 Planck analysis, namely the flat model in Eq. (12), Non-Bunch-Davies (NBD; Chen et al. 2007b), NBD1 and NBD2 models (Agullo & Parker 2011; now called “NBD1 cos” and “NBD2 cos”) and NBD3 (Chen 2010b). We also introduce three new templates, NBD sin which is motivated by (Chen 2010a) and takes the form (21)where again K = k1 + k2 + k3 and . The other two templates are extensions of the NBD1 cos and NBD2 cos models found in Agullo & Parker (2011) and take the form (22)where , which is dominated by squeezed configurations, and , which has a flattened shape.

2.7. Directional-dependence motivated by gauge fields

Some models where primordial vector fields are present during inflation predict interesting NG signatures. This is the case of a coupling of the inflaton field ϕ to the kinetic term of a gauge field Aμ, contains I2(ϕ)F2, where Fμν = μAννAμ and the coupling I2(ϕ)F2 is chosen so that scale invariant vector perturbations are produced on superhorizon scales (Barnaby et al. 2012b; Bartolo et al. 2013a). The bispectrum turns out to be the sum of two contributions: one of the local shape; and another that is also enhanced in the squeezed limit (k1k2k3), but featuring a non-trivial dependence on the angle between the small and the large wave vectors through the parameter (where ) as . Also, primordial magnetic fields sourcing curvature perturbations can cause a dependence on both μ and μ2 (Shiraishi et al. 2012).

We can parametrize these shapes as variations on the local shape (Shiraishi et al. 2013a), as (23)where PL(μ) is the Legendre polynomial with P0 = 1,P1 = μ, and . For example, for L = 1 we have the shape (24)The local template corresponds to ci = 2fNLδi0. Here and in the following the nonlinearity parameters are related to the cL coefficients by , , and . The L = 1,2 shapes exhibit sharp variations in the flattened limit, for example for k1 + k2k3, while in the squeezed limit, L = 1 is suppressed whereas L = 2 grows like the local bispectrum shape (i.e., the L = 0 case). The I2(ϕ)F2 models predict c2 = c0/ 2, while primordial curvature perturbations sourced by large-scale magnetic fields generate non-vanishing c0, c1, and c2. Quite interestingly, in the proposed “solid inflation” scenario (Endlich et al. 2013, 2014; see also Bartolo et al. 2013b, 2014; Sitwell & Sigurdson 2014) bispectra similar to Eq. (23) can be generated, in this case with c2c0 (Endlich et al. 2013, 2014). Therefore, measurements of the ci coefficients can be an efficient probe of some detailed aspects of the inflationary mechanism, such as the existence of primordial vector fields during inflation (or a non-trivial symmetry structure of the inflaton fields, as in solid inflation).

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

The same shift symmetry that leads to axion (monodromy) models of inflation (Sect. 2.4) naturally allows (from an effective field theory point of view) for a coupling between a pseudoscalar axion inflaton field and a gauge field of the type , where the parameter α is dimensionless and f is the axion decay constant (). This scenario has a rich and interesting phenomenology both for scalar and tensor primordial fluctuations (see, e.g., Barnaby & Peloso 2011; Sorbo 2011; Barnaby et al. 2011, 2012c; Linde et al. 2013; Meerburg & Pajer 2013; Ferreira & Sloth 2014). Gauge field quanta are produced by the background motion of the inflaton field, and these in turn source curvature perturbations through an inverse decay process of the gauge field. A bispectrum of curvature fluctuations is generated as (Barnaby et al. 2011; Meerburg & Pajer 2013)4(25)where the exact expression for the function f3 can be found in Eq. (3.29) of Barnaby et al. (2011; see also Meerburg & Pajer 2013). Here ξ characterizes the coupling strength of the axion to the gauge field . The inverse decay bispectrum peaks for equilateral configuration, since δϕ is mostly sourced by the inverse decay (δA + δAδϕ), when two modes of the vector fields are of comparable magnitude (the correlation with the equilateral template is 94% and with the orthogonal one is 4%). We do however constrain the exact shape in Eq. (25), without resorting to the equilateral template. Another interesting observational signature that can shed light on the role played by pseudo-scalars in the early Universe is provided by tensor NG, to which we turn next.

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

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

Recently, the possibility of observable tensor bispectra has been vigorously discussed in a model where the inflaton couples to a pseudoscalar field (Barnaby et al. 2012a; Cook & Sorbo 2013; Ferreira & Sloth 2014). In this model, through the gravitational coupling to the U(1) gauge field, gravitational waves () receive NG corrections, where only one of the two spin states is enhanced. The bispectrum, generally formed as (26)is accordingly polarized, with . This NG enhancement is a sub-horizon effect and therefore is maximized at the equilateral limit (k1k2k3) (Cook & Sorbo 2013).

A model-independent template of the equilateral-type polarized tensor bispectrum is given by (Shiraishi et al. 2013b, 2015) (27)with the polarization tensor obeying and . We here have introduced a tensor nonlinearity parameter, by normalizing with the equilateral bispectrum template of curvature perturbations () in the equilateral limit, yielding (28)The template Eq. (27) can adequately reconstruct the tensor bispectra created in the pseudoscalar inflation models6 (Shiraishi et al. 2013b), and thus the amplitude is directly connected with the model parameters, e.g., the coupling strength of the pseudoscalar field to the gauge field ξ (for details see Sect. 11).

The CMB temperature and E-mode bispectra sourced by the parity-violating tensor NG have not only the usual parity-even (1 + 2 + 3 = even) signals but also parity-odd (1 + 2 + 3 = odd) contributions, which cannot be sourced by known scalar bispectra (Kamionkowski & Souradeep 2011; Shiraishi et al. 2011). Moreover, their shapes are mostly distinct from the scalar templates, due the different radiation transfer functions; hence they can be measured essentially independently of the scalar NG (Shiraishi et al. 2013b). The analysis of the WMAP temperature data distributed in 1 + 2 + 3 = odd configurations leads to an observational limit (Shiraishi et al. 2015). This paper updates the limit, by analysing both parity-even and parity-odd signals in the Planck temperature and E-mode polarization data.

3. Statistical estimation of the CMB bispectrum for polarized maps

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

The CMB temperature and polarization fields are characterized using the multipoles of a spherical harmonic decomposition of the CMB maps: (29)At linear order, the relation between the primordial perturbation field and the CMB multipoles is (e.g., Komatsu & Spergel 2001) (30)where X = { T,E } denotes either temperature or E-mode polarization, Φ is the primordial gravitational potential, and represents the linear CMB radiation transfer function.

The CMB angular bispectrum is the three-point correlator of the aℓms (31)where Xi = { T,E }. If the CMB sky is rotationally invariant, and the bispectra we are considering have even parity (which is true for combinations of T and E), then the angular bispectrum can be factorized as (32)where is the so-called reduced bispectrum, and is the Gaunt integral, defined as the integral over the solid angle of the product of three spherical harmonics, (33)The Gaunt integral (often written in terms of Wigner 3j-symbols) enforces rotational symmetry, and restricts attention to a tetrahedral domain of multipole triplets { 1,2,3 }, satisfying both a triangle condition and a limit given by some maximum resolution max (the latter being defined by the finite angular resolution of the experiment under study).

Our goal is to extract the nonlinearity parameter fNL from the data, for different primordial shapes. To achieve this, we essentially fit a theoretical CMB bispectrum ansatz b123 to the observed 3-point function. Theoretical predictions for CMB angular bispectra arising from early Universe primordial models can be obtained by applying Eq. (30) to the primordial bispectra of Sect. 2, (see e.g., Komatsu & Spergel 2001). Optimized cubic bispectrum estimators were introduced by Heavens (1998), and it has been shown that for small NG the general optimal polarized fNL estimator can be written as (Creminelli et al. 2006) (34)where N is a suitable normalization chosen to produce unit response to . Note that we are implicitly defining a suitable normalization convention so that , and is the value of the theoretical template when fNL = 1. C-1 is the inverse of the block matrix: (35)and the blocks represent the full TT, TE, and EE covariance matrices, with CET being the transpose of CTE. All quantities in the previous equation (i.e., CMB multipoles, bispectrum template and covariances matrices) are assumed to properly incorporate instrumental beam and noise.

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

The inverse covariance filtering operation implied by Eq. (34) is a challenging numerical task, which has been successfully performed only recently (Smith et al. 2009; Elsner & Wandelt 2012). This step can be avoided by working in the “diagonal covariance approximation”. In this approach, the estimator is built by neglecting off-diagonal entries of the covariance matrix in the cubic term in Eq. (34), and then finding the linear term that minimizes the variance for this specific cubic statistic. Applying such a procedure yields (Yadav et al. 2007) (36)where is the inverse of the 2 × 2 matrix (37)This expression can also be written as (38)where the observed (reduced) bispectrum includes the linear correction term and the inner product is defined as (39)with (40)where the last term is the Wigner 3j symbol. The denominator in Eq. (38), bth,bth is the normalization constant N.

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

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

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

3.1. KSW and skew-C estimators

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

3.2. Modal estimators

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

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

The other implementation, which we will refer to as “Modal 2”, utilizes a novel approach where the and are first orthogonalized to produce new uncorrelated unit variance âℓm coefficients, We then decompose the new bispectra as (43)which can be constrained independently, since they are uncorrelated. In this case the estimator then takes on a particularly simple form (Fergusson 2014). This new form is mathematically equivalent to the previous modal method, but involves significantly fewer terms in the estimator. However, due to the orthogonalization procedure we cannot constrain the full EEE bispectrum without further processing, just the additional part which is orthogonal to temperature. For this reason, although the “Modal 2” T+E results incorporate all the polarization information, the EEE results alone are not presented here.

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

3.3. Binned bispectrum estimator

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

In more detail, the computation for the binned bispectrum estimator is based on Eqs. (38) and (39). However, instead of using the reduced bispectrum , all expressions start from the alternative rotationally-invariant reduced bispectrum , where h is defined in Eq. (40). The expression in Eq. (39) for the inner product remains the same when replacing b by B, except that the h2 becomes h-2. The importance of is that it can be determined directly from maximally-filtered maps; (44)where (45)Binning is then implemented by adding a sum over all inside a bin to the expression for the filtered map given in Eq. (45), thus obtaining the observed binned bispectrum of the map , with bin indices i1,i2,i3. The linear correction term is obtained in a similar way (and subtracted from the cubic term), but with two of the maps in Eq. (44) replaced by Gaussian simulations, and taking the average over a large number of those.

The theoretical templates are binned simply by summing the exact expression over the s inside a bin8, and the same is true for the covariance matrix . The binning is optimized in such a way as to maximize the overlap, defined using the inner product of Eq. (39) between the binned and the exact template for all shapes under consideration. Finally the estimate of fNL is computed using Eq. (38), where the inner product now contains a sum over bin indices i instead of multipoles , and the bispectra and covariance matrix are replaced by their binned versions.

3.4. Data set and simulations

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

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

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

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

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

3.5. Data analysis settings

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

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

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

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

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

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

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

4. Non-primordial contributions to the CMB bispectrum

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

4.1. Non-Gaussianity from the lensing-ISW bispectrum

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

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

To determine , the amplitude parameter of the lensing-ISW bispectrum, one simply inserts the theoretical template for this shape into the general expression of Eq. (38). The template is given by (Hu 2000; Lewis et al. 2011) (46)Here and are the temperature/polarization-lensing potential cross power spectra, and the tilde on , , and indicates that it is the lensed TT, TE, or EE power spectrum. The functions are defined by (47)if 1 + 2 + 3 is even and 1,2,3 satisfy the triangle inequality, and zero otherwise.

In this paper our main concern with the lensing-ISW bispectrum is not so much to determine its amplitude (although that is also of great interest), but to compute its influence on the primordial shapes. The bias due to the lensing-ISW bispectrum on the estimation of a given primordial amplitude is given by (48)where the inner product is defined in Eq. (39).

Table 1

Bias in the three primordial fNL parameters due to the lensing-ISW signal for the four component separation methods.

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

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

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

thumbnail Fig. 1

Skew-C spectrum for the lensing-ISW effect (red line with data points), from the temperature map. The blue curve is the theoretically-expected spectrum. Note that the points beyond = 1500 are significantly correlated.

Open with DEXTER

Table 2

Results for the amplitude of the lensing-ISW bispectrum from the SMICA, SEVEM, NILC, and Commander foreground-cleaned maps, for different bispectrum estimators.

4.2. Non-Gaussianity from extragalactic point sources

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

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

The unclustered sources have a white noise distribution, and hence constant polyspectra. Their reduced angular bispectrum template is thus (49)This constant is usually noted bPS or bsrc in the literature (e.g., Komatsu & Spergel 2001). This constant template is valid in polarization as well as temperature, since the polarization angles of point sources are less clustered than the source density. However, since not all these point sources are polarized, we do not measure the same sources in temperature and in polarization. In fact, there is no detection of the bispectrum of unclustered point sources in the cleaned Planck polarization map, unlike in the temperature map, where Table 3 (binned bispectrum estimator) and Fig. 2 (skew-Cs) show a clear detection.

thumbnail Fig. 2

Skew-C spectrum for unclustered point sources (red line with data points), from the temperature map. The blue curve is the theoretical spectrum, given the amplitude determined with the KSW estimator.

Open with DEXTER

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

Table 3

Joint estimates of the bispectrum amplitudes of unclustered and clustered point sources in the cleaned Planck temperature map, determined with the binned bispectrum estimator.

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

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

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

4.3. Non-Gaussianity from residuals of the deglitching processing

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

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

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

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

Table 4

Results on the impact of cosmic ray residuals on the estimation of fNL at 143GHz, determined using the binned bispectrum estimator.

Table 5

Results from the different fNL estimators for the set of CMB simulations described in Sect. 5.1 in the ideal case without noise or mask.

5. Validation tests

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

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

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

Table 6

As Table 5, but with noise and no mask.

Table 7

As Table 5, but with noise and a mask.

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

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

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

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

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

thumbnail Fig. 3

Map-by-map comparison of the results from the different estimators for local (left), equilateral (centre), and orthogonal (right) fNL (taking the shapes to be independent), for the third set of simulations described in Sect. 5.1, including both noise and a mask. Results are shown for T-only (top), E-only (centre), and the full combined T+E case (bottom). The legend for the estimators can be found in the top right figure. The horizontal solid line is the average value of all maps for KSW, and the dashed and dotted horizontal lines correspond to ± 1σ and ± 2σ deviations, respectively.

Open with DEXTER

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

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

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

Table 8

Results from the different estimators for fNL for the set of CMB simulations with described in Sect. 5.2.

5.3. Validation of estimators on realistic Planck simulations

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

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

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

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

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

Table 9

Results from the different estimators for fNL for the set of SMICA simulations based on FFP8 described in Sect. 5.3.

Table 10

Results for the fNL parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW, binned and modal estimators from the SMICA, NILC, SEVEM, and Commander foreground-cleaned maps.

6. Results

6.1. Constraints on local, equilateral, and orthogonal fNL

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

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

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

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

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

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

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

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

6.2. Bispectrum reconstruction

6.2.1. Modal bispectrum reconstruction

thumbnail Fig. 4

Modal bispectrum reconstruction for Planck 2013 (top left) and 2015 (top right) temperature-only data, both using the SMICA maps. Here, we restrict the 2015 resolution to the same as 2013, using similar polynomials with nmax = 601. The two bispectra are very close to being in complete agreement in the signal-dominated regime shown up to max = 1500. In the lower panel, we show the Planck 2015 temperature bispectrum at high resolution using the full nmax = 2001 polynomial modes. Large-scale features in the top panels become subdivided but the main 2013 signals remain, notably a stronger measurement of the ISW-lensing signal (the regular oscillations in the squeezed limit).

Open with DEXTER

thumbnail Fig. 5

Modal reconstruction for the WMAP-9 bispectrum (left) and the Planck SMICA 2015 T-only bispectrum (right) plotted for the domain ≤ 450, using identical isosurface levels. Here, we employed the full 2001 eigenmodes for both the Planck analysis at max = 2000 and for WMAP-9 analysis at max = 600, but for comparison purposes we have only used the first 600 eigenmodes in order to obtain a comparable resolution. The main features in the WMAP-9 bispectrum have counterparts in the Planck version, revealing an oscillatory pattern in the central region, as well as features on the tetrapyd surface. The WMAP-9 bispectrum has a much larger noise signal beyond = 350 than the more sensitive Planck experiment, leading to apparent residuals in this region.

Open with DEXTER

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

Table 11

Results for the fNL parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW estimator from the SMICA foreground-cleaned map.

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

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

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

thumbnail Fig. 6

CMB temperature and polarization bispectrum reconstructions for PlanckSMICA maps using the full set of polynomial modes with nmax = 2001 and with signal-to-noise weighting. The top bispectra are the symmetric pure temperature TTT (left) plotted with ≤ 1500 and E-mode polarization EEE (right) shown for 30 ≤ ≤ 1100. Below are the mixed temperature/polarization bispectra with TTE on the left (with E multipoles in the z-direction) and TEE on the right (with T multipoles in the z-direction). All S/N thresholds are the same.

Open with DEXTER

thumbnail Fig. 7

Comparison of CMB polarization bispectrum EEE reconstructions for PlanckNILC, SEVEM, and Commander foreground-separated maps with signal-to-noise weighting. Note that these results are not as internally consistent between the four methods, also comparing SMICA shown in Fig. 6, which is closest to NILC. We will compare the underlying modal coefficients below to demonstrate these differences quantitatively.

Open with DEXTER

thumbnail Fig. 8

Smoothed binned signal-to-noise bispectrum for the Planck 2015 cleaned sky map, as determined with the binned estimator, as a function of 1 and 2 for a fixed 3-bin [518, 548]. From left to right results are shown for the four component separation methods SMICA, SEVEM, NILC, and Commander. From top to bottom are shown: TTT, TTT cleaned from radio and CIB point sources; T2E, TE2; and EEE. The colour range is in signal-to-noise from − 4 to + 4. The light grey regions are where the bispectrum is not defined, either because it is outside the triangle inequality or because of the cut .

Open with DEXTER

thumbnail Fig. 9

Similar to Fig. 8, but with 3 ∈ [ 1291,1345 ].

Open with DEXTER

6.2.2. Binned bispectrum reconstruction

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

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

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

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

6.3. Primordial curvature reconstruction

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

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

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

thumbnail Fig. 10

Signal-to-noise for T and E as a function of multipole L and comoving distance χ for an ideal cosmic-variance limited experiment (upper panels) and for the Planck experiment, with noise determined from FFP8 simulations (lower panels). The complementary nature of the information provided by T and E is evident. The differential visibility is shown for comparison.

Open with DEXTER

thumbnail Fig. 11

Mean field and fluctuations for differential visibility projected ζ, as described in the text. The filter used in these maps is a 40′ FWHM Gaussian. The grey regions are masked out.

Open with DEXTER

When both T and E are included as constraints, a two by two matrix appears which includes the CTE correlation in the off-diagonal as well as the CTT and CEE along the diagonal: Figure 11 shows an all-sky mean field (Wiener-filtered) map of the curvature variable ζ constructed in (densely-packed) shells from the multipoles ζLMκ(χ). The figure actually shows ζ projected onto differential visibility, w = deτ/ dχ. The map ζw(θ,φ) is the spherical transform of ζLMκw = w(χ)dχζLMκ(χ). The top left map is ζw | T using the all-sky maps SMICA DX11. The middle left map is ζw | E, based on just the polarization information. (The E-maps did not have high pass filtering imposed, so that an indication of the full results obtainable with E alone can be seen.) The bottom left panel shows ζw | T,E. Visually the combined two-mode constraint has more detail than either of the one-mode constrained maps, enhancing the estimation of ζw. In each of the three cases, two sample realizations are shown which include allowed fluctuations about the mean map. These illustrate the level of uncertainty in the patterns found in the mean maps, necessary ingredients to the full statistical description. Note that the fluctuations about the mean are larger for T or E alone than for the combined T,E.

thumbnail Fig. 12

Mean field and fluctuations for a χdec slice of ζ, as described in the text. The filter used in these maps is a 40′ FWHM Gaussian. The grey regions are masked out.

Open with DEXTER

thumbnail Fig. 13

Mean field and fluctuations for differential visibility projected ζ, for a 20 deg by 20 deg patch, as described in the text. The filter used in these maps is a 20′ FWHM Gaussian. Sources in grey are masked out.

Open with DEXTER

thumbnail Fig. 14

Mean field and fluctuations for a χdec slice of ζ, for a 20 deg by 20 deg patch, as described in the text. The filter used in these maps is a 20′ FWHM Gaussian. Sources in grey are masked out.

Open with DEXTER

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

Figure 10 quantifies where the information resides in Lχ space. Here, the signal variance is the ensemble average of the square of the mean field, assuming the ζ are drawn from a Gaussian with fiducial covariance, (56)and the noise (i.e., fluctuation) variance is the covariance matrix, as given by Eq. (54) for T and the equivalent for E. Figure 10 plots the signal-to-noise at each individual χ slice (hence ignores the components off-diagonal in χ): (57)for T, with an equivalent expression for E. The signal-to-noise structure in χ can be contrasted with the differential visibility plotted in the bottom panel: Fig. 10 shows that the Lχ information has a reach beyond the differential visibility structure, especially for low multipoles, because the associated waves can straddle the last scattering surface. The “reionization bump” in signal-to-noise seen in Fig. 10 is not that prominent. The Integrated Sachs Wolfe impact on the S/N is evident, but is relatively low, with the consequence that we cannot draw out high significance results from the ISW effect for ζ reconstruction (or equivalently for gravitational potential reconstruction).

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

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

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

7. Validation of Planck results

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

7.1. Dependence on foreground cleaning method

7.1.1. Comparison between fNl measurements

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

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

Table 12

Comparison between local, equilateral, and orthogonal fNL results, obtained using the four different component separation pipelines.

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

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

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

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

7.1.2. Comparison between reconstructed bispectra

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

Table 13

For each of the four foreground cleaned maps, we compute fNL for the local, equilateral, and orthogonal modes using two different polarization masks, one with fsky = 0.74 and the other with fsky = 0.64, while for temperature we use a single mask with fsky = 0.76.

Table 14

Correlation coefficients between pairs of bispectrum modes, extracted using two different component-separated maps.

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

thumbnail Fig. 15

Scatter plots showing correlations between bispectrum modes extracted from the different Planck foreground-cleaned maps, for all possible pairs of component separation methods. Top: TTT bispectrum modes. Bottom: EEE bispectrum modes. While temperature shows a strong correlation, the loss of correlation in polarization between the different methods is evident in these plots, as discussed in the text and quantified in Tables 14 and 15. Results here and in Table 15 have been obtained using the high resolution modal pipeline (2001 modes), while results in Table 14 have been obtained with the low resolution modal pipeline. By construction, the high resolution pipeline is measuring not the full EEE bispectrum of the map, but the component of EEE that is orthogonal to TTT. For this reason, r2 measured by the two pipelines for EEE will not be identical. With this caveat in mind, the agreement between the two modal approaches is very good.

Open with DEXTER

Table 15

r2 statistic, Eq. (58), showing the degree of correlation between measured bispectrum β coefficients for the component separation methods shown in Fig. 15 (upper three rows for temperature, lower for polarization).

7.2. Dependence on sky coverage

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

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

Table 16

Correlation coefficients between pairs of EEE bispectrum modes, extracted using two different masks for each of the four component-separated maps.

Table 17

Comparison of component separation methods, using Gaussian FFP8 simulations.

7.3. Tests on simulations

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

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

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

thumbnail Fig. 16

Evolution of the fNL parameters (solid blue line with data points) and their uncertainties (dashed lines) for the three primordial bispectrum templates as a function of the maximum multipole number max used in the analysis. From left to right the figures show, respectively local, equilateral, and orthogonal, while the different rows from top to bottom show results for T-only, E-only, and full T+E. To indicate more clearly the evolution of the uncertainties, they are also plotted around the final value of fNL (solid green lines without data points, showing the ± 1σ range around the dashed green lines). The results are for SMICA, and assume all shapes to be independent, and that the ISW-lensing bias has been subtracted. They have been determined with the binned bispectrum estimator.

Open with DEXTER

thumbnail Fig. 17

Evolution of the fNL parameters (solid blue line with data points) and their uncertainties (dashed lines) for the three primordial bispectrum templates as a function of the minimum multipole number min used in the analysis. From left to right the figures show respectively, local, equilateral, and orthogonal, while the different rows from top to bottom show results for T-only, E-only, and full T+E. To indicate more clearly the evolution of the uncertainties, they are also plotted around the final value of fNL (solid green lines without data points, with ± 1σ range around the dashed green line). The results are for SMICA, and assume all shapes to be independent, as well as that the ISW-lensing bias has been subtracted. They have been determined with the binned bispectrum estimator.

Open with DEXTER

Table 18

Same test as in Table 17, but with a NG map as input, with .

7.4. Dependence on multipole number

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

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

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

As we noted in the 2013 analysis (Planck Collaboration XXIV 2014), when going to the much lower WMAP resolution of max ≃ 500, we agree with the slightly high value of that the WMAP team reported (Bennett et al. 2013). This value is also confirmed when including polarization. One can clearly see the value of the higher resolution of Planck.

The dependence on min is shown in Fig. 17. Here all estimators used the same values, and . As explained in Planck Collaboration VII (2016), not all systematic and foreground uncertainties in the low- HFI polarization data have been fully characterized yet, and hence it was decided to filter out these data.

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

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

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

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

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

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

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

  • 1.

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

  • 2.

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

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

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

thumbnail Fig. 18

Temperature only, local shape, pixel correlation matrix from Monte Carlo analysis, at Nside = 4. The most correlated pixels are those closest to the main diagonal; however, these values are always lower than 34% in the chosen case.

Open with DEXTER

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

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

thumbnail Fig. 19

Temperature only, local fNL directional contributions from SMICA. As explained in the text, summing over all the pixel values would give the full sky fNL needlet estimator result. The left panel displays the directional fNL map. On the right, the blue points represent the fNL contibution for each direction (i.e., for each pixel in the directional map), with Monte Carlo error bars. The red line is the average from simulations, which is consistent with zero.

Open with DEXTER

7.6. Summary of the main validation results

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

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

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

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

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

8. Other non-Gaussianity shapes for fNL

Table 19

Results for local isocurvature NG, determined from the SMICAPlanck 2015 map with the binned bispectrum estimator.

Table 20

Similar to Table 19, except that we now assume that the adiabatic and isocurvature mode are completely uncorrelated.

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

8.1. Isocurvature non-Gaussianity

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

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

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

8.2. Feature models

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

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

thumbnail Fig. 20

Constant feature model results for both T-only and T+E data across a wide frequency range. The upper four panels show the feature signal in the Modal range 0 <ω< 350. The two upper left panels show contours of the raw significance σ obtained from the SMICA map as a function of the frequency ω, for T-only and T+E, respectively. The upper right panels show the maximum signal after marginalizing over phase φ for both T-only and T+E for all foreground separation models. The third and fourth panels show the maximum feature signal in both T-only and T+E across the frequency range 0 <ω< 1000, plotting both Modal results (dashed lines) and KSW results (solid lines for 200 <ω< 1000); these show good agreement in the overlap. The lower two panels give the maximum KSW results for T-only and T+E in the range 1000 <ω< 3000).

Open with DEXTER

thumbnail Fig. 21

Constant feature ansatz validation for the Modal estimator, showing the effect of ISW-lensing and point source subtraction for < 300 (left panel) and the impact of a lower maximum cutoff on the average signal (right panel), i.e., lowering max = 2000 to 1500 (T) and max = 1500 to 1000 (E). All Modal 2 results in Sect. 8 have used the extended common mask, except the validation analysis at different resolutions (right panel) which for consistency employs the common mask.

Open with DEXTER

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

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

thumbnail Fig. 22

Generalized feature models analysed at max = 2000 (E-modes max = 1500) for the different Planck foreground separation methods, SMICA (blue), SEVEM (red), NILC (green), and Commander (yellow), together with the SSN (SMICA - SEVEM - NILC ) average (black). The left three panels apply to the equilateral feature models, showing, respectively, in the top panel the full feature survey significance at each frequency and phase (temperature only), the maximum significance at each frequency for temperature only (middle), and with polarization (lower). The right three panels apply to the flattened feature models, showing the maximum significance at each frequency for temperature only (top right) and with polarization (middle right), along with significance at each frequency and phase for temperature and polarization (right lower).

Open with DEXTER

thumbnail Fig. 23

Single field feature model significance with a K2cosωK scaling dependence (Eq. (19); left panels, T-only upper and T + Elower) or with a KsinωK scaling (Eq. (20); right panels). To find the maximum signal, these results have been marginalized over the α-dependent envelope function ranging from α → 0 (no envelope) to the maximum cutoff allowed by the Modal resolution αω = 90.

Open with DEXTER

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

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

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

Table 21

Peak statistics for the different feature models showing the Raw peak maximum significance (for the given Modal survey domain), the corrected significance of this Single maximum peak after accounting for the parameter survey size (the “look-elsewhere” effect) and the Multi-peak statistic which integrates across the adjusted significance of all peaks to determine consistency with Gaussianity.

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

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

thumbnail Fig. 24

Generalized resonance models analysed at max = 2000 (E-modes max = 1500) for the different Planck foreground separation methods, SMICA (blue), SEVEM (red), NILC (green), and Commander (yellow), together with the SSN (SMICASEVEMNILC ) average (black). The upper panels apply to the constant resonance model (Eq. (10)), with T-only (left) and T+E (right), the middle panels give results for the equilateral resonance model (Eq. (13)), and the lower panels for the flattened resonance model (Eq. (14)). Both the equilateral and flattened resonance models produce broad peaks, which are reinforced with polarization (middle and bottom right panels).

Open with DEXTER

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

Table 22

Peak statistics for the resonance models showing the maximum Raw peak significance, the Single peak significance after accounting for the parameter survey “look-elsewhere” effect, and the Multi-peak statistic integrating across all peaks (also accounting for the “look-elsewhere” correction).

thumbnail Fig. 25

Standard resonance model results for both T-only and T+E across a wide frequency range using the high frequency resonance model estimator. The first and second panels show the signal in both T-only (first panel) and T+E (second panel) across the frequency range 0 <ω < 500. The lower two panels give the results for T-only (third panel) and T+E (fourth panel) in the range 500 <ω < 1100.

Open with DEXTER

8.3. Resonance/axion monodromy models

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

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

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

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

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

Table 23

Constraints on models with equilateral-type NG covering the shapes predicted by the effective field theory of inflation, together with constraints on specific non-canonical inflation models, such as DBI inflation.

Table 24

Constraints on models with excited initial states (non-Bunch-Davies models), as well as warm inflation.

Table 25

Direction-dependent NG results for both the L = 1 and L = 2 models.

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

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

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

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

8.6. Direction-dependent primordial non-Gaussianity

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

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

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

8.7. Parity-violating tensor non-Gaussianity

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

The parity-violating NG under examination induces non-vanishing signals not only in parity-even configurations (1 + 2 + 3 = even) but also in the parity-odd ones (1 + 2 + 3 = odd) in the temperature and E-mode polarization bispectra (Shiraishi et al. 2013b). The optimal estimator, including all (even + odd) bispectrum signals, is expressed by the linear combination of the parity-even and parity-odd estimators, reading (Liguori, in prep.) (61)where Neven / odd is the normalization factor (related to the Fisher matrix as Neven / odd = 6Feven / odd) defined for 1 + 2 + 3 = even / odd. The parity-even estimator can be computed using the original Modal methodology (Fergusson et al. 2010a, 2012; Fergusson 2014; Liguori, in prep.), while in computations of the parity-odd estimator , we follow the extended spin-weighted pipeline (Shiraishi et al. 2014, 2015; Liguori, in prep.).

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

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

Table 26

Results for the tensor nonlinearity parameter , estimated from the SMICA, SEVEM, and NILC temperature and high-pass filtered polarization maps.

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

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

Table 27

Correlation coefficients between pairs of bispectrum modes, extracted from two of the Planck component separated maps and the WMAP foreground-cleaned map at max = 500 resolution.

thumbnail Fig. 26

Parity-odd signals (1 + 2 + 3 = odd) of the TTT (top) and EEE (bottom) bispectra (i ≤ 500) recovered from the SMICA maps by means of the Modal decomposition with 101 simple polynomial-based eigenmodes, not including any special functions fitting the CMB tensor-mode features. In the panel for EEE, only the signals larger than = 40 are shown. The TTT and EEE bispectra shown here are rescaled with a constant Sachs-Wolfe weighting and signal-to-noise weighting, respectively.

Open with DEXTER

The parity-odd components of the TTT and EEE bispectra extracted model-independently from the SMICA data are visually represented in Fig. 26. It is apparent from this figure that the PlanckTTT bispectrum has similar features to the WMAP one (Shiraishi et al. 2015), e.g., distinctive signals distributed around 123. As indicated by the roughly 70% correlation between the SMICA and WMAP bispectra (see Table 27), the Planck T-only limits in Table 26 are close to the WMAP ones (68% CL): for parity-even; and for parity-odd (Shiraishi et al. 2015).

9. Limits on the primordial trispectrum

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

First, some notation: a “primed” ζ-trispectrum ζk1ζk2ζk3ζk4 ⟩ ′ denotes the connected trispectrum without its momentum-conserving delta function, i.e., (62)where “+ disc.” denotes the disconnected contributions to the 4-point function.

One possible signal is the “local” trispectrum , which arises if the non-Gaussian adiabatic curvature ζ is of cubic-type form (see, e.g., Okamoto & Hu 2002), i.e., (63)where is a free parameter and ζG is a Gaussian field. In this model, the bispectrum is zero (since there is a ζ → − ζ symmetry) and the trispectrum is given by (64)Analogously to the case of the local bispectrum, , the observational signal-to-noise for is largest in the “squeezed” limit, k1 ≪ min(k2,k3,k4), and there is a consistency relation which shows that in single-field inflation, the four-point function is always small in the squeezed limit (e.g., Senatore & Zaldarriaga 2012a). Thus can only be detectably large in multi-field models. Conversely, there are multi-field models where is detectable. The main obstacle here is technical naturalness, i.e., ensuring that radiative corrections do not generate an observationally larger bispectrum. This can be the case if a large bispectrum is forbidden by a Z2 symmetry, or by supersymmetry (Senatore & Zaldarriaga 2012b).

A further category of four-point signals can be obtained by adding quartic interactions to the inflationary action. Following Smith et al. (2015), we will concentrate on the simplest possibility, by considering quartic operators consistent with the symmetries of inflation and having the lowest possible number of derivatives (Bartolo et al. 2010b; Senatore & Zaldarriaga 2011, 2012b). There are three such operators, of the schematic form , , and (iσ)2(jσ)2. By a short calculation using the in-in formalism (Maldacena 2003), the associated four-point functions are (Smith et al. 2015; see also Huang & Shiu 2006): Here K = ∑ ki and we have introduced parameters , , and to parametrize the amplitude of each trispectrum.

The normalization of the gNL-parameters is chosen so that for “tetrahedral” configurations, with | ki | = k and (ki·kj) = − k3/ 3. This is the analogue of the commonly-used normalization for the bispectrum, where fNL parameters are defined so that for equilateral configurations with | ki | = k.

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

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

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

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

Summarizing, we will search for the following trispectrum signals: (68)defined by Eqs. (64), (65), and (67) above.

9.1. Data analysis

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

Table 28

Number of factorizable terms Nin needed to represent each trispectrum by direct sampling of the integral, and number of terms Nout obtained after running the optimization algorithm from Sect. VII of Smith et al. (2015).

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

  • 1.

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

  • 2.

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

  • 3.

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

  • 4.

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

  • 5.

    We apply a spherical harmonic transform to the inpainted, apodized CMB map to obtain a harmonic-space map aℓm with max = 1600. We then take the final filtered map to be (69)where b is the beam, C is the fiducial CMB power spectrum, and N is the sky-averaged noise power spectrum (without beam deconvolution). To motivate this choice of -weighting, we note that for an ideal all-sky experiment with isotropic noise, we have aℓm = bsℓm + nℓm where sℓm,nℓm are signal and noise realizations. In this case, Eq. (69) weights the signal as , which is optimal.

In our pipeline, we apply this filter to the component-separated SMICA maps (Planck Collaboration IX 2016), obtaining a harmonic-space map . We apply the same filter to 1000 Monte Carlo simulations to obtain an ensemble of harmonic-space maps. Our pipeline has the property that it always estimates the trispectrum of the data in excess of the trispectrum in the simulations. Since the simulations include lensing, this means that lensing bias will automatically be subtracted from our gNL estimates.

Now that the filter, data realization, and Monte Carlo simulations have been fully specified, the details of the pipeline are described in Sect. IX.B of Smith et al. (2015). For each trispectrum, the pipeline outputs an estimate of gNL and an estimate of the statistical error. Our basic results are: (70)No deviation from Gaussian statistics is seen. These results significantly improve the previous best constraints on the trispectrum from WMAP (Vielva & Sanz 2010; Smidt et al. 2010; Fergusson et al. 2010b; Hikage & Matsubara 2012; Sekiguchi & Sugiyama 2013; Regan et al. 2013; Smith et al. 2015) and large-scale structure (Desjacques & Seljak 2010; Giannantonio et al. 2014; Leistedt et al. 2014).

A constraint on from Planck 2013 data was recently reported by Feng et al. (2015), who find . Our central value in Eq. (70) agrees well with this result, but the statistical error is smaller by a factor of 2.3. This improvement is partly due to the lower noise levels in Planck 2015 data, and partly due to the use of a better estimator.

Each line in Eq. (70) is a “single-gNL” constraint; i.e., the constraint on one gNL parameter with the other gNL-parameters held fixed. For joint constraints, one needs to know the full covariance matrix. The correlation between and the other two parameters is negligble, and the - correlation is: (71)Multi-field models of inflation will generally give a linear combination of , , and (iσ)2(jσ)2 trispectra. In this case we proceed as follows. First, if the coefficient is non-zero, we can use the near-degeneracy with a linear combination of the other two operators to absorb it into the effective values of and . A Fisher matrix analysis shows that the coefficients of this linear combination are (72)It is convenient to define the two-component parameter vector: (73)We also compute a two-by-two Fisher matrix Fij, whose diagonal is given by , where σi is the single-gNL statistical error in Eq. (70), and whose off-diagonal is , where r is the correlation in Eq. (71). This procedure gives: (74)For a given parameter vector gi, we can define a trispectrum-χ2 by (75)where ĝi = ( − 0.21 × 106, − 0.10 × 105) is the vector of best-fit single-gNL values from Eq. (70). This definition of χ2 follows from the observation that (Fiiĝi) is an estimator with expectation value (Fg)i and covariance matrix Cov(Fiiĝi,Fjjĝj) = Fij.

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

10. Minkowski functionals results

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

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

The most recent constraints (1σ) from MFs on and have been obtained with WMAP (Hikage & Matsubara 2012) and Planck (Planck Collaboration XXIV 2014): (76)

10.1. Method and definition of MFs

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

We calculate the four normalized18 functionals vk(ν) on nth = 26 thresholds ν, between νmin = − 3.5 and νmax = + 3.5.

thumbnail Fig. 27

Filters used to optimize constraints on local NG, in harmonic space. The temperature filter WM is a smoothed version of the true Wiener filter obtained with realistic models, while the E-mode WM filter is adapted from the temperature one, with a cutoff value at ≃ 800. The formulae for the derivative filters are given in Eq. (78).

Open with DEXTER

For this analysis, we used the same temperature and polarization E data, simulations and masks described in Sect. 3.4 for consistency with the bispectrum estimators. In addition, the maps are filtered to optimize constraints on local NG (Ducout et al. 2013), the filters used being similar to Wiener filters for T and E (WM), and for the first (WD1) and second (WD2) derivatives of these fields (Fig. 27): For the temperature map, known point sources in the mask are inpainted.

We define the vector y as any combination with k = { 0,3 }, A = { T,E }, W = { WM,WD1,WD2 }, ŷ being the vector measured on the data.

From these measurements, we then use a Bayesian method to jointly estimate and , (79)We take a uniform prior for in the range − 400 to 400, and for in the range − 4 × 106 to + 4 × 106, while the evidence is just considered as a normalization.

Assuming MFs are multi-variate Gaussian-distributed we obtain the posterior distribution for with a χ2 test (80)with (81)For this test, we use two types of simulations to first construct a model including primordial NG and secondly a covariance matrix (82)with averaged over the simulations. We now describe the details of these simulations.

  • Simulations 1: Non-Gaussian model For the firsttype of simulation, we included all possible sources of NG,assuming that the total and individual levels of NG are smallenough that MFs are linear with respect to those NG levels(Ducout et al. 2013). The threekinds of NG we included are foreground residuals (Galacticresiduals with scalable amplitude α, as well as radio sources, CIB anisotropies, secondaries (SZ, lensing and ISW-lensing, but not SZ-lensing) and primordial NG (, ): (83)We tried to reproduce all instrumental effects, with realistic effective beams (isotropic window functions), noise from FFP8 simulations (Planck Collaboration XII 2016), filtered with component separation weights. We checked the accuracy of these simulations by comparing them to FFP8 simulations, using no foreground and no primordial NG. The astrophysical models are provided by the Planck Sky Model (PSM, Delabrouille et al. 2013), while the primordial NG simulations are computed as in Elsner & Wandelt (2009). The lensing uses LensPix19. The power spectrum used for these NG simulations and the lensing is the best-fit power spectrum from Planck 2013+ACT/SPT+BAO (Planck Collaboration XXII 2014). We created n1 = 200 simulations i, using n1 maps for the primordial NG, while we had only one astrophysical foreground simulation.

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

Validation of the estimator

Part of the validation for the MFs estimator is described in Sect. 5.3 for , to compare the results to bispectrum estimators on realistic simulations (FFP8 MC, second item above). In addition we present in Table 29 the results obtained on the same realistic simulations for , and on simulations containing primordial NG (first item above), with and , with 200 simulations used in each case. These tests have been performed using the SMICA method with lensing bias removed.

Table 29

Results for local NG parameters at first and second order, and , obtained with Minkowski functionals on SMICA simulations in temperature and polarization.

10.2. Results

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

The posteriors for and from SMICA are shown in Fig. 28. One interesting point is that the estimates of and are almost uncorrelated (r< 0.1); this can be inferred when we consider the parity of MF deviations from Gaussianity, which is different for the two parameters (Matsubara 2010).

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

Table 30

Results for local NG parameters at first and second order, and , obtained with Minkowski functionals on all four component separated maps in temperature and polarization.

thumbnail Fig. 28

Joint constraint on and obtained with MFs. The contour lines represent 1, 2 and 3σ limits of a 2D-Gaussian distribution. Constraints were obtained with SMICA temperature and polarization E maps.

Open with DEXTER

Foreground and secondary biases

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

Lensing has a significant signature in MF estimation of , but is even stronger in (the four-point correlation signature) and could be detected (and not treated just as a bias) with this estimator. The Wiener filters enhance the scales where lensing is dominant.

Table 31

Biases for local NG parameters at first and second order , obtained with Minkowski Functionals on SMICA in temperature and polarization.

11. Implications for early Universe physics

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

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

11.1. General single-field models of inflation

DBI models: DBI models of inflation (Silverstein & Tong 2004; Alishahiha et al. 2004), characterized by a non-standard kinetic term of the inflaton field, predict a non-linearity parameter , where cs is the sound speed of the inflaton perturbations (Silverstein & Tong 2004; Alishahiha et al. 2004; Chen et al. 2007b). The corresponding bispectrum shape is very close to the equilateral shape. Nonetheless we have constrained the exact theoretical (non-separable) shape (see Eq. (7) of Planck Collaboration XXIV 2014). The constraint we obtain from temperature data ( from temperature and polarization) at 68% CL (with ISW-lensing and point sources subtracted, see Table 23) implies (84)and (85)In Planck Collaboration XXIV (2014) we constrained the so-called infrared (IR) DBI models (Chen 2005b,a), which arise in string frameworks. We focused on a minimal setup, considering a regime where stringy effects are negligible and predictions for primordial perturbations are built within standard field theory. In the companion paper Planck Collaboration XX (2016) we present an analysis of a more general class of IR DBI models which accounts for stringy signatures (see Bean et al. 2008) by combining Planck power spectrum and bispectrum constraints.

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

The effective field theory approach (Cheung et al. 2008; Weinberg 2008) provides an efficient way to constrain inflationary perturbations for various classes of models that incorporate deviations from the standard single-field slow-roll scenario. In this approach the Lagrangian of the system is expanded into the (lowest dimension) operators obeying the underlying symmetries. We consider general single-field models described by the following action (86)where the curvature perturbation is related to the scalar field π as ζ = − . The inflaton interaction terms and generate two kind of bispectra with amplitudes, respectively, and , where M3 is the amplitude of the operator (Senatore et al. 2010; see also Chen et al. 2007b; Chen 2010b). These two bispectra both peak for equilateral triangles in Fourier space. Nevertheless, they are sufficiently different, and the total NG signal turns out to be a linear combination of the two, leading also to an orthogonal shape. We can put constraints on cs and the dimensionless parameter (Senatore et al. 2010). Notice that DBI inflationary models corresponds to having , while cs = 1 and M3 = 0 (or represent the non-interacting (vanishing NG) case.

thumbnail Fig. 29

68%, 95%, and 99.7% confidence regions in the parameter space , defined by thresholding χ2 as described in the text.

Open with DEXTER

The mean values of the estimators for and are expressed in terms of cs and by (87)Here the coefficients come from the Fisher matrix between the equilateral and orthogonal templates and the theoretical bispectra predicted by the two operators and . We use a χ2 statistic given by , where (i = {equilateral, orthogonal}), being the joint estimates of equilateral and orthogonal fNL (see Table 11), C the covariance matrix of the joint estimators, and is provided by Eq. (87). As an example in Fig. 29 we show the 68%,95%, and 99.7% confidence regions for and obtained from the T + E constraints, requiring χ2 ≤ 2.28,5.99, and 11.62 respectively (corresponding to a χ2 variable with two degrees of freedom). In Fig. 30 we show the corresponding confidence regions in the parameter space. Marginalizing over we find (88)and (89)The constraints improve by a few percent in T-only and by up to 25% by including polarization, in comparison with those of Planck Collaboration XXIV (2014).

thumbnail Fig. 30

68%, 95%, and 99.7% confidence regions in the single-field inflation parameter space , obtained from Fig. 29 via the change of variables in Eq. (87).

Open with DEXTER

Galileon models of inflation: Galileon models of inflation (Burrage et al. 2011; Kobayashi et al. 2010; Mizuno & Koyama 2010; Ohashi & Tsujikawa 2012) are well motivated models based on the so called “Galilean symmetry” (Nicolis et al. 2009). They are characterized by stability properties that are quite well understood (ghost-free, and stable against quantum corrections) and can arise naturally within fundamental physics models (de Rham & Gabadadze 2010b,a). Moreover they are an interesting example of models where gravity is modified on large scales and we focus on them also as a typical example of a more general class of modified gravity theories that are ghost-free (the so called Horndeski theories, Horndeski 1974; see also Planck Collaboration XIII 2016 where these models are discussed in the context of dark energy/modified gravity scenarios for the late time evolution of the universe). The predictions for the primordial perturbations are very rich. Bispectra can be generated with the same shapes as the “EFT1” and “EFT2” bispectra (see also discussion in Creminelli et al. 2011), however the amplitude(s) scale with the fluctuation sound speed as , differently from the general single-field models of inflation considered in the above subsection. They can be written as (at the lowest-order in slow-roll parameters) Here A, , and cs are dimensionless parameters of the models. In particular cs is the sound speed of the Galileon scalar field, while is a parameter that appears to break the standard consistency relation for the tensor-to-scalar perturbation ratio (, nT being the tensor spectral index)20. Accordingly to Eq. (87) a linear combination of these two bispectra generate equilateral and orthogonal bispectra templates21. From the Planck constraints on and (see Table 11),

thumbnail Fig. 31

68%, 95%, and 99.7% probability contours in the Galileon models for cs and parameters for the branch (tensor spectral index nT< 0).

Open with DEXTER

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

thumbnail Fig. 32

68%, 95%, and 99.7% probability contours in the Galileon models for cs and parameters for the branch (blue tensor spectral index nT> 0).

Open with DEXTER

11.2. Multi-field models

Curvaton models: the simplest adiabatic curvaton models predict local NG with an amplitude (Bartolo et al. 2004c,b) (96)for a quadratic potential of the curvaton field (Lyth & Wands 2002; Lyth et al. 2003; Lyth & Rodriguez 2005; Malik & Lyth 2006; Sasaki et al. 2006), where rD = [ 3ρcurvaton/ (3ρcurvaton + 4ρradiation) ] D is the “curvaton decay fraction” at the time of the curvaton decay in the sudden decay approximation. Assuming a uniform prior, 0 <rD< 1, our constraint at 68% CL (see Table 11) yields (97)while accounting for temperature and polarization data ( at 68% CL) gives (98)improving the previous Planck bound which was previously rD ≥ 0.15 (95% CL; Planck Collaboration XXIV 2014). In Planck Collaboration XX (2016), assuming there is some relic isocurvature fluctuations in the curvaton field, a limit on rD is derived from the bounds on isocurvature fluctuations. In this restricted case, the limit rD> 0.98 (95% CL) is derived, which is consistent with the constraint given here.

thumbnail Fig. 33

68% and 95% confidence regions in the plane, with the Lorentz invariant model in Eq. (107) shown as the dashed line.

Open with DEXTER

Notice that the above expression of Eq. (96) is valid under the assumption that there is no significant decay of the inflaton into curvaton particles. In general one should account for such a possibility. For example, if the classical curvaton field survives and starts to dominate, then the curvaton particles produced during reheating (which have the same equation of state as the classical curvaton field) are expected to survive and dominate over other species at the epoch of their decay. The classical curvaton field and the curvaton particles decay at the same time, inevitably producing adiabatic perturbations (for a detailed discussion see Linde & Mukhanov 2006). A general formula for , accounting for the possibility that the inflaton field decays into curvaton particles, is provided in Sasaki et al. (2006): (99)where measures the ratio of the energy density of curvaton particles to the energy density of the classical curvaton field (Linde & Mukhanov 2006; Sasaki et al. 2006) and ρcurvaton in the expression for rD is given by ρcurvaton = ρcurv.particles + ρcurv.field. Using the uniform prior 0 <rD< 1 and our measurements of constrain at 95% CL (T) and (T+E).

11.3. Non-standard inflation models

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

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

NG from gauge field production during axion inflation: we have also constrained the inverse decay NG of Eq. (25) arising typically in models where the inflaton field is a pseudoscalar axion that couples to a gauge field. Using the modal estimator we get the following constraints (removing ISW-lensing bias): (100)for temperature only; and (101)from temperature+polarization (see Table 23). The NG amplitude is given by , where is the power spectrum of vacuum-mode curvature perturbations (i.e., the power spectrum predicted without the coupling to gauge fields), is the dimensionless scalar power spectrum of curvature perturbations (a star denoting evaluation at the pivot scale). The NG parameter is exponentially sensitive to the strength of the coupling between the axion and the gauge field. From Eq. (101) we limit the strength of the coupling to ξ ≤ 2.5   (95% CL). The details together with constraints on the axion decay constant can be found in Planck Collaboration XX (2016), where an overview of various observational limits on axion (monodromy) models of inflation is presented. This limit is in agreement with the one that can be derived from tensor non-Gaussianities (see below).

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

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

11.4. Alternatives to inflation

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

11.5. Inflationary interpretation of CMB trispectrum results

We briefly interpret the trispectrum constraints in an inflationary context. First, consider the case of single field inflation. The action for the Goldstone boson π is highly constrained by residual diffeomorphism invariance (see, e.g., Smith et al. 2015). To lowest order in the derivative expansion, the most general action is: (102)where the parameter M4 is related to the trispectrum by (103)The constraint in Eq. (70) translates to the following parameter constraint in single field inflation: (104)This constraint is a factor 1.8 better than WMAP.

Turning now to multifield inflation, we consider an action of the form (105)where σ is a spectator field that acquires quantum fluctuations with power spectrum Pσ(k) = H2/ (2k3) and converts to adiabatic curvature via ζ = (2Aζ)1 / 2H-1σ. The trispectrum in this model is (106)so we can constrain its parameters by thresholding the χ2 defined in Eq. (75). For example, if we consider the Lorentz invariant model (107)so that the parameters Λi of the more general action in Eq. (105) are given by , then by thresholding at Δχ2 = 4 (as appropriate for one degree of freedom), we obtain the following constraint on the parameter Λ: (108)Constraints in other parameter spaces can also be obtained by thresholding the χ2 defined in Eq. (75). For example we show 68% and 95% confidence regions in the -plane, obtained by thresholding at χ2 = 2.28 and χ2 = 5.99 as appropriate for a χ2 random variable with two degrees of freedom.

DBI trispectrum: the trispectrum constraint on the shape in Eq. (70) can also be used to obtain a lower bound on the DBI model sound speed. This is because, in the small sound speed limit (Chen et al. 2009; Arroja et al. 2009), the dominant contribution to the contact interaction trispectrum (Huang & Shiu 2006) has this shape. The corresponding non-linearity parameter is . We follow the same procedure as described at the beginning of this Sect. 11 and, assuming a uniform prior in the range 0 ≤ cs ≤ 1 / 5, we can derive a constraint on cs as (109)This constraint is consistent with the ones derived from the bispectrum measurements (see Eqs. (84) and (85)) and it is only a factor of about three worse. Notice, however, that in this case we are ignoring the scalar exchange contribution, which is of the same order in cs.

Curvaton trispectrum: for the simplest curvaton scenario, the trispectrum non-linearity parameter prediction is (Sasaki et al. 2006) (110)Following the procedure described at the beginning of Sect. 11, we use the observational constraint obtained in Sect. 9 (Eq. (70)), and the same prior (0 <rD< 1) as in 11.2, to obtain a lower bound on the curvaton decay fraction as (111)This limit is consistent with the previous ones derived using the bispectrum measurements and it is a factor of about 3 to 4 worse.

12. Conclusions

In this paper we have presented the constraints on primordial NG using the full Planck mission data. The results have improved compared to the Planck 2013 release (Planck Collaboration XXIV 2014) as a consequence of including data from the full mission and taking advantage of Planck’s polarization capability – the first time that maps of the CMB polarization anisotropies have been used to constrain primordial NG.

Using temperature data alone, the constraints on the local, equilateral, and orthogonal bispectrum templates are , , and . Moving from the nominal Planck 2013 data to the full mission data yielded modest improvements of up to 15% (for the orthogonal shape). After the inclusion of full mission polarization data, our final constraints become , , and , which represents a substantial step forward relative to Planck 2013, with error bars shrinking by 14% (local), 43% (equilateral), and 46% (orthogonal). These improved limits on the standard shapes enhance our understanding of different inflationary models that can potentially lead to subtle effects beyond the simplest models of inflation.

The reason that the polarization data provide such complementary constraints on primordial curvature perturbations is due to the phase shift of the CMB polarization transfer functions compared to the temperature transfer functions. So, despite the relatively much lower signal-to-noise in the polarization maps, their inclusion leads to appreciable improvements on limits on NG parameters. Nevertheless, the full characterization of the noise properties in the polarized maps is still ongoing. In spite of the extensive testing and cross-checks validating the combined temperature and polarization results, we therefore conservatively recommend that all results that include polarization information, not just the polarization-only results, be taken as preliminary at this stage.

The complementary nature of the polarization information also represents an important cross-check on the analysis. The Planck results based on polarization alone are statistically consistent with the results based on temperature alone, with a precision comparable to that achievable in an optimal analysis of the WMAP 3-yr temperature maps (Spergel et al. 2007; Creminelli et al. 2007; Yadav & Wandelt 2008).

Our analysis was subject to an extensive validation exercise. In addition to extensive simulation tests, including for the first time a detailed test of the impact of time-domain de-glitching, our results are supported by tests for robustness under change of estimator implementations (KSW, binned bispectrum, and two modal estimators), and variations in sky coverage as well as upper and lower multipole cutoffs. We also tested for possible directional dependence using a needlet estimator. These tests form the basis of our selection of SMICA as the main foreground cleaning method for our headline results.

The Planck 2015 analysis presented here provides constraints on a greatly extended range of template families. These extensions include a tenfold increase in the range of frequencies covered in feature models, giving rise to linearly oscillating bispectra, generalized shapes for oscillating models including for logarithmic oscillations, tests for deviations from the Bunch-Davies vacuum, models of equilateral type in the context of the effective field theory of inflation, and direction-dependent primordial NG. Beyond purely scalar mode templates we also tested for parity-violating tensor NG.

Using the full mission data with polarization, we have investigated the hints of NG reported in the Planck 2013 analysis of oscillatory features. While no individual feature or resonance model rises above our detection threshold of 3σ (after inclusion of the look-elsewhere effect), the results of integrated (multi-peak) statistical tests indicate that continued investigation of oscillatory and non-scaling models is warranted.

In addition to searches for specific NG templates, we present model-independent reconstructions of the temperature and polarization bispectra using the modal and binned bispectrum approaches. These full mission reconstructions can achieve twice the resolution of the Planck 2013 results, demonstrating excellent consistency in temperature, and good agreement with the WMAP9 reconstruction in regions where this earlier data set is signal-dominated.

The inclusion of polarization information leads to significantly improved constraints on NG in primordial isocurvature perturbations, providing complementary information to 2-point function results for models where the NG in isocurvature components is more easily detectable than its contribution to the power spectrum.

A significant addition to this year’s analysis is the inclusion of detailed trispectrum results due to cubic NG. The local trispectrum is constrained by Planck temperature data to be and the other two shapes were also found to be consistent with Gaussianity. Both 3-point and 4-point constraints are consistent with the improved (though still suboptimal) constraints from Minkowski functionals, a very different estimation framework. This concordance adds confidence in our results.

We have discussed the implications of our results on the physics of the early Universe, showing that the n-point functions for n> 2 provide a significant window onto the primordial Universe beyond the power spectrum, constraining general-single field, multifield, and non-standard inflation models, as well as alternatives to inflation. Using bispectrum and trispectrum limits we updated results on the parameter space of the inflationary models (and alternatives) already tested in 2013, and constrained the parameter space of other well-motivated inflationary models (e.g., Galileon-like models of inflation, and models where axion/pseudoscalar fields are present during inflation).

The global picture that emerges is one of consistency with the premises of the ΛCDM cosmology, namely that the structure we observe today is the consequence of passive evolution of adiabatic, Gaussian, primordial seed perturbations. Nevertheless, NG at some level is expected in all inflationary models, and hence we should strive to find means to reduce errors on fNL still further.


1

Planck (http://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).

2

We refer the reader to Planck Collaboration XXIV (2014) and references therein for a detailed summary of the models and underlying physical mechanisms generating various types of primordial NG.

3

NG perturbations can arise also at the end of inflation, e.g., from nonlinearities during the (p)reheating phase (e.g., Enqvist et al. 2005; Chambers & Rajantie 2008; Barnaby & Cline 2006; see also Bond et al. 2009) or from fluctuations in the inflaton decay rate or interactions, as found in modulated (p)reheating and modulated hybrid inflation (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).

4

For simplicity we assume a scale-invariant bispectrum.

5

See Planck Collaboration XIX (2016) for the Planck constraints on magnetically-induced NG.

6

The form of the tensor bispectrum is the same whether the inflaton field is identified with the pseudoscalar field or not.

7

We note that the local, equilateral, and orthogonal templates of Sect. 2 are separable. In fact, while the theoretical local NG models are manifestly separable, the equilateral and orthogonal templates of Eqs. (3) and (4) are factorizable approximations of the original non-separable shapes, that were derived exactly with the purpose of allowing the application of this type of estimator (Creminelli et al. 2006; Senatore et al. 2010).

8

We note that the enormous computational gain of the binned bispectrum estimator comes from the binned determination of the observed bispectrum; determining the theoretical bispectrum templates is fast, even when done exactly.

9

We note that the Planck collaboration produced two slightly different sets of union masks (see Planck Collaboration IX 2016 for details). We choose to adopt the more conservative set in this paper, as we found that the agreement between different component separation methods significantly increases with these masks when we measure fNL of shapes that peak in the squeezed limit (while the differences are very small in other cases).

10

For bispectrum purposes we found no difference between the results when performing the procedure without a buffer (the so-called “Gauss-Seidel” method, where amongst the neighbours will be pixels both at the current and at the previous iteration) and with a buffer (the so-called “Jacobi” method, where all neighbour pixels will be at the previous iteration), except that the former converges faster. We found 2000 iterations to work well in the “Gauss-Seidel” case.

11

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]).

12

The average value of the lensing-ISW amplitude determined from the FFP8 simulations is around 0.85 of the expected value. This value is very consistent across bispectrum estimators and component separation methods, which provides a useful consistency test in its own right. Except for this effect, all other tests on the temperature FFP8 maps show them to be very robust and to behave as expected, for example in the determination of the lensing-ISW bias on the local shape. We took this effect into account by increasing all error bars in the table by the appropriate factor (i.e., dividing them by ≃ 0.85).

13

Since the Minkowski-functional pipeline automatically subtracts the ISW-lensing bias, the theoretical value for the bias as computed from the Fisher matrix has been added to its results, to make a direct comparison possible.

14

The E-only fNL agreement is still at a reasonable 1σfNL level in most cases. However this is larger than expectations from simulations, as described in Sect. 7.

15

Even though we focus here on directional contributions to the local shape (which was typically found as one of the most sensitive to residual foreground contamination), this type of directional analysis can be done in a model-independent way, by looking separately at different needlet modes; we leave this for future work.

16

Compared to definitions in the literature based on ζ and S (see e.g., Langlois & van Tent 2012), here we adopt definitions based on Φadi = 3ζ/ 5 and Φiso = S/ 5, in order to make the link with the standard adiabatic result more direct. Conversion factors to obtain results based on ζ and S are 6/5, 2/5, 2/15, 18/5, 6/5, and 2/5, for the six modes, respectively.

17

For the feature model with parameters (ω,φ), the adjusted significance S for the raw significance σ after accounting for the “look-elsewhere” effect is given by (Fergusson et al. 2015) (59) where Fχ,2 is the cumulative distribution function of the χ-distribution of degree two and Neff is the effective number of independent feature models. This can also be used to investigate whether feature models are contributing at several frequencies. This multi-peak statistic integrates over all peak signals using the corrected significance S, i.e., (60) where Δω is the sampling step-size and Δωeff is the effective correlation scale between independent models given by Δωeff = (ωmaxωmin) / (Neff − 1). Essentially this sums up all significant peaks that have a non-zero adjusted significance, after accounting for the look-elsewhere effect in Eq. (59).

18

Raw Minkowski functionals Vk depend on the Gaussian part of fields through a normalization factor Ak, a function only of the shape of the power spectrum. We therefore normalize functionals vk = Vk/Ak to focus on NG, see Appendix B.

19

http://cosmologist.info/lenspix

20

For the explicit expressions of these parameters in terms of the coefficients of the Galileon Lagrangian see Planck Collaboration XX (2016).

21

We note that we are neglecting corrections (where means also ); see Burrage et al. (2011) and Ribeiro & Seery (2011). These corrections will have a different shape associated with them and they are not necessarily small when compared with some of the terms displayed, e.g., the terms .

22

A component is a connected subset of the excursion.

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); 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/planck-collaboration. Some of the results in this paper have been derived using the HEALPix package. 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 IN2P3 Computer Center (http://cc.in2p3.fr) for providing a significant amount of the computing resources and services needed for the analysis with the binned bispectrum estimator. This research used resources of the National Energy Research Scientific Computing Center, 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. We also acknowledge the IAP magique3 computer facilities. Some computations were performed on the GPC 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

Appendix A: Derivation of an estimator for cL

As parameterized by Eq. (23), we express a primordial bispectrum of direction-dependence: (A.1)where is a Legendre Polynomial of order L. It can be shown that such a primordial bispectrum leads to a CMB bispectrum (A.2)where p denotes either temperature or E-mode polarization, is the reduced CMB bispectrum, and the term with large parentheses denotes the Wigner-3j symbol. The reduced CMB bispectrum is given by (Shiraishi et al. 2013a) (A.3)where δL,l denotes the Kronecker delta function, { ... } the Wigner 6j symbol, “perm.” means permutations, w is a beam window function, and Here, is the radiation transfer function for temperature or E-mode polarization, and j(x) is a spherical Bessel function. The h symbol is (A.6)By maximizing the likelihood with respect to cL, we obtain the KSW estimator for cL: (A.7)where MC denotes that the average is over Monte Carlo simulations and is a normalization constant. The normalization constant is given by (A.8)Using Eq. (A.2), we find that (A.9)where with being defined in Eq. (A.5). In the derivation above, we used the identities and Applying Eqs. (A.9) to (A.7), we find where Since , is not a real function, but a complex function. We estimate efficiently by computing the following with HEALPix (Górski et al. 2005): Here

Appendix B: Definition of Minkowski Functionals and theoretical expectations

For a field f(x) of zero average and variance defined on the two-dimensional sphere S2, an overdense excursion set is defined as (B.1)The boundary of the excursion is (B.2)Then the three Minkowski functionals on the sphere are where and dl are respectively elements of solid angles (surface) and of angle (distance), κ is the geodesic curvature. Note that the Genus can be also expressed as the number of components22 in the excursion minus the number of holes in the excursion.

The fourth functional V3(ν), is defined, for ν> 0, as the number of components in the excursion. Symmetrically, for ν< 0, it is the number of underdense components (or the number of components in the excursion { x ∈ S2 | f(x) <νσ0 }).

In the Gaussian limit, the functionals can be expressed the following way (see, e.g., Matsubara 2010; Vanmarcke 1983): (B.6)with and (B.9)The amplitude Ak depends only on the shape of the power spectrum C: where ωkπk/ 2/ Γ(k/ 2 + 1), which gives ω0 = 1, ω1 = 2, ω2 = π and σ0 and σ1 are respectively the rms of the field and its first derivatives.

All Tables

Table 1

Bias in the three primordial fNL parameters due to the lensing-ISW signal for the four component separation methods.

Table 2

Results for the amplitude of the lensing-ISW bispectrum from the SMICA, SEVEM, NILC, and Commander foreground-cleaned maps, for different bispectrum estimators.

Table 3

Joint estimates of the bispectrum amplitudes of unclustered and clustered point sources in the cleaned Planck temperature map, determined with the binned bispectrum estimator.

Table 4

Results on the impact of cosmic ray residuals on the estimation of fNL at 143GHz, determined using the binned bispectrum estimator.

Table 5

Results from the different fNL estimators for the set of CMB simulations described in Sect. 5.1 in the ideal case without noise or mask.

Table 6

As Table 5, but with noise and no mask.

Table 7

As Table 5, but with noise and a mask.

Table 8

Results from the different estimators for fNL for the set of CMB simulations with described in Sect. 5.2.

Table 9

Results from the different estimators for fNL for the set of SMICA simulations based on FFP8 described in Sect. 5.3.

Table 10

Results for the fNL parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW, binned and modal estimators from the SMICA, NILC, SEVEM, and Commander foreground-cleaned maps.

Table 11

Results for the fNL parameters of the primordial local, equilateral, and orthogonal shapes, determined by the KSW estimator from the SMICA foreground-cleaned map.

Table 12

Comparison between local, equilateral, and orthogonal fNL results, obtained using the four different component separation pipelines.

Table 13

For each of the four foreground cleaned maps, we compute fNL for the local, equilateral, and orthogonal modes using two different polarization masks, one with fsky = 0.74 and the other with fsky = 0.64, while for temperature we use a single mask with fsky = 0.76.

Table 14

Correlation coefficients between pairs of bispectrum modes, extracted using two different component-separated maps.

Table 15

r2 statistic, Eq. (58), showing the degree of correlation between measured bispectrum β coefficients for the component separation methods shown in Fig. 15 (upper three rows for temperature, lower for polarization).

Table 16

Correlation coefficients between pairs of EEE bispectrum modes, extracted using two different masks for each of the four component-separated maps.

Table 17

Comparison of component separation methods, using Gaussian FFP8 simulations.

Table 18

Same test as in Table 17, but with a NG map as input, with .

Table 19

Results for local isocurvature NG, determined from the SMICAPlanck 2015 map with the binned bispectrum estimator.

Table 20

Similar to Table 19, except that we now assume that the adiabatic and isocurvature mode are completely uncorrelated.

Table 21

Peak statistics for the different feature models showing the Raw peak maximum significance (for the given Modal survey domain), the corrected significance of this Single maximum peak after accounting for the parameter survey size (the “look-elsewhere” effect) and the Multi-peak statistic which integrates across the adjusted significance of all peaks to determine consistency with Gaussianity.

Table 22

Peak statistics for the resonance models showing the maximum Raw peak significance, the Single peak significance after accounting for the parameter survey “look-elsewhere” effect, and the Multi-peak statistic integrating across all peaks (also accounting for the “look-elsewhere” correction).

Table 23

Constraints on models with equilateral-type NG covering the shapes predicted by the effective field theory of inflation, together with constraints on specific non-canonical inflation models, such as DBI inflation.

Table 24

Constraints on models with excited initial states (non-Bunch-Davies models), as well as warm inflation.

Table 25

Direction-dependent NG results for both the L = 1 and L = 2 models.

Table 26

Results for the tensor nonlinearity parameter , estimated from the SMICA, SEVEM, and NILC temperature and high-pass filtered polarization maps.

Table 27

Correlation coefficients between pairs of bispectrum modes, extracted from two of the Planck component separated maps and the WMAP foreground-cleaned map at max = 500 resolution.

Table 28

Number of factorizable terms Nin needed to represent each trispectrum by direct sampling of the integral, and number of terms Nout obtained after running the optimization algorithm from Sect. VII of Smith et al. (2015).

Table 29

Results for local NG parameters at first and second order, and , obtained with Minkowski functionals on SMICA simulations in temperature and polarization.

Table 30

Results for local NG parameters at first and second order, and , obtained with Minkowski functionals on all four component separated maps in temperature and polarization.

Table 31

Biases for local NG parameters at first and second order , obtained with Minkowski Functionals on SMICA in temperature and polarization.

All Figures

thumbnail Fig. 1

Skew-C spectrum for the lensing-ISW effect (red line with data points), from the temperature map. The blue curve is the theoretically-expected spectrum. Note that the points beyond = 1500 are significantly correlated.

Open with DEXTER
In the text
thumbnail Fig. 2

Skew-C spectrum for unclustered point sources (red line with data points), from the temperature map. The blue curve is the theoretical spectrum, given the amplitude determined with the KSW estimator.

Open with DEXTER
In the text
thumbnail Fig. 3

Map-by-map comparison of the results from the different estimators for local (left), equilateral (centre), and orthogonal (right) fNL (taking the shapes to be independent), for the third set of simulations described in Sect. 5.1, including both noise and a mask. Results are shown for T-only (top), E-only (centre), and the full combined T+E case (bottom). The legend for the estimators can be found in the top right figure. The horizontal solid line is the average value of all maps for KSW, and the dashed and dotted horizontal lines correspond to ± 1σ and ± 2σ deviations, respectively.

Open with DEXTER
In the text
thumbnail Fig. 4

Modal bispectrum reconstruction for Planck 2013 (top left) and 2015 (top right) temperature-only data, both using the SMICA maps. Here, we restrict the 2015 resolution to the same as 2013, using similar polynomials with nmax = 601. The two bispectra are very close to being in complete agreement in the signal-dominated regime shown up to max = 1500. In the lower panel, we show the Planck 2015 temperature bispectrum at high resolution using the full nmax = 2001 polynomial modes. Large-scale features in the top panels become subdivided but the main 2013 signals remain, notably a stronger measurement of the ISW-lensing signal (the regular oscillations in the squeezed limit).

Open with DEXTER
In the text
thumbnail Fig. 5

Modal reconstruction for the WMAP-9 bispectrum (left) and the Planck SMICA 2015 T-only bispectrum (right) plotted for the domain ≤ 450, using identical isosurface levels. Here, we employed the full 2001 eigenmodes for both the Planck analysis at max = 2000 and for WMAP-9 analysis at max = 600, but for comparison purposes we have only used the first 600 eigenmodes in order to obtain a comparable resolution. The main features in the WMAP-9 bispectrum have counterparts in the Planck version, revealing an oscillatory pattern in the central region, as well as features on the tetrapyd surface. The WMAP-9 bispectrum has a much larger noise signal beyond = 350 than the more sensitive Planck experiment, leading to apparent residuals in this region.

Open with DEXTER
In the text
thumbnail Fig. 6

CMB temperature and polarization bispectrum reconstructions for PlanckSMICA maps using the full set of polynomial modes with nmax = 2001 and with signal-to-noise weighting. The top bispectra are the symmetric pure temperature TTT (left) plotted with ≤ 1500 and E-mode polarization EEE (right) shown for 30 ≤ ≤ 1100. Below are the mixed temperature/polarization bispectra with TTE on the left (with E multipoles in the z-direction) and TEE on the right (with T multipoles in the z-direction). All S/N thresholds are the same.

Open with DEXTER
In the text
thumbnail Fig. 7

Comparison of CMB polarization bispectrum EEE reconstructions for PlanckNILC, SEVEM, and Commander foreground-separated maps with signal-to-noise weighting. Note that these results are not as internally consistent between the four methods, also comparing SMICA shown in Fig. 6, which is closest to NILC. We will compare the underlying modal coefficients below to demonstrate these differences quantitatively.

Open with DEXTER
In the text
thumbnail Fig. 8

Smoothed binned signal-to-noise bispectrum for the Planck 2015 cleaned sky map, as determined with the binned estimator, as a function of 1 and 2 for a fixed 3-bin [518, 548]. From left to right results are shown for the four component separation methods SMICA, SEVEM, NILC, and Commander. From top to bottom are shown: TTT, TTT cleaned from radio and CIB point sources; T2E, TE2; and EEE. The colour range is in signal-to-noise from − 4 to + 4. The light grey regions are where the bispectrum is not defined, either because it is outside the triangle inequality or because of the cut .

Open with DEXTER
In the text
thumbnail Fig. 9

Similar to Fig. 8, but with 3 ∈ [ 1291,1345 ].

Open with DEXTER
In the text
thumbnail Fig. 10

Signal-to-noise for T and E as a function of multipole L and comoving distance χ for an ideal cosmic-variance limited experiment (upper panels) and for the Planck experiment, with noise determined from FFP8 simulations (lower panels). The complementary nature of the information provided by T and E is evident. The differential visibility is shown for comparison.

Open with DEXTER
In the text
thumbnail Fig. 11

Mean field and fluctuations for differential visibility projected ζ, as described in the text. The filter used in these maps is a 40′ FWHM Gaussian. The grey regions are masked out.

Open with DEXTER
In the text
thumbnail Fig. 12

Mean field and fluctuations for a χdec slice of ζ, as described in the text. The filter used in these maps is a 40′ FWHM Gaussian. The grey regions are masked out.

Open with DEXTER
In the text
thumbnail Fig. 13

Mean field and fluctuations for differential visibility projected ζ, for a 20 deg by 20 deg patch, as described in the text. The filter used in these maps is a 20′ FWHM Gaussian. Sources in grey are masked out.

Open with DEXTER
In the text
thumbnail Fig. 14

Mean field and fluctuations for a χdec slice of ζ, for a 20 deg by 20 deg patch, as described in the text. The filter used in these maps is a 20′ FWHM Gaussian. Sources in grey are masked out.

Open with DEXTER
In the text
thumbnail Fig. 15

Scatter plots showing correlations between bispectrum modes extracted from the different Planck foreground-cleaned maps, for all possible pairs of component separation methods. Top: TTT bispectrum modes. Bottom: EEE bispectrum modes. While temperature shows a strong correlation, the loss of correlation in polarization between the different methods is evident in these plots, as discussed in the text and quantified in Tables 14 and 15. Results here and in Table 15 have been obtained using the high resolution modal pipeline (2001 modes), while results in Table 14 have been obtained with the low resolution modal pipeline. By construction, the high resolution pipeline is measuring not the full EEE bispectrum of the map, but the component of EEE that is orthogonal to TTT. For this reason, r2 measured by the two pipelines for EEE will not be identical. With this caveat in mind, the agreement between the two modal approaches is very good.

Open with DEXTER
In the text
thumbnail Fig. 16

Evolution of the fNL parameters (solid blue line with data points) and their uncertainties (dashed lines) for the three primordial bispectrum templates as a function of the maximum multipole number max used in the analysis. From left to right the figures show, respectively local, equilateral, and orthogonal, while the different rows from top to bottom show results for T-only, E-only, and full T+E. To indicate more clearly the evolution of the uncertainties, they are also plotted around the final value of fNL (solid green lines without data points, showing the ± 1σ range around the dashed green lines). The results are for SMICA, and assume all shapes to be independent, and that the ISW-lensing bias has been subtracted. They have been determined with the binned bispectrum estimator.

Open with DEXTER
In the text
thumbnail Fig. 17

Evolution of the fNL parameters (solid blue line with data points) and their uncertainties (dashed lines) for the three primordial bispectrum templates as a function of the minimum multipole number min used in the analysis. From left to right the figures show respectively, local, equilateral, and orthogonal, while the different rows from top to bottom show results for T-only, E-only, and full T+E. To indicate more clearly the evolution of the uncertainties, they are also plotted around the final value of fNL (solid green lines without data points, with ± 1σ range around the dashed green line). The results are for SMICA, and assume all shapes to be independent, as well as that the ISW-lensing bias has been subtracted. They have been determined with the binned bispectrum estimator.

Open with DEXTER
In the text
thumbnail Fig. 18

Temperature only, local shape, pixel correlation matrix from Monte Carlo analysis, at Nside = 4. The most correlated pixels are those closest to the main diagonal; however, these values are always lower than 34% in the chosen case.

Open with DEXTER
In the text
thumbnail Fig. 19

Temperature only, local fNL directional contributions from SMICA. As explained in the text, summing over all the pixel values would give the full sky fNL needlet estimator result. The left panel displays the directional fNL map. On the right, the blue points represent the fNL contibution for each direction (i.e., for each pixel in the directional map), with Monte Carlo error bars. The red line is the average from simulations, which is consistent with zero.

Open with DEXTER
In the text
thumbnail Fig. 20

Constant feature model results for both T-only and T+E data across a wide frequency range. The upper four panels show the feature signal in the Modal range 0 <ω< 350. The two upper left panels show contours of the raw significance σ obtained from the SMICA map as a function of the frequency ω, for T-only and T+E, respectively. The upper right panels show the maximum signal after marginalizing over phase φ for both T-only and T+E for all foreground separation models. The third and fourth panels show the maximum feature signal in both T-only and T+E across the frequency range 0 <ω< 1000, plotting both Modal results (dashed lines) and KSW results (solid lines for 200 <ω< 1000); these show good agreement in the overlap. The lower two panels give the maximum KSW results for T-only and T+E in the range 1000 <ω< 3000).

Open with DEXTER
In the text
thumbnail Fig. 21

Constant feature ansatz validation for the Modal estimator, showing the effect of ISW-lensing and point source subtraction for < 300 (left panel) and the impact of a lower maximum cutoff on the average signal (right panel), i.e., lowering max = 2000 to 1500 (T) and max = 1500 to 1000 (E). All Modal 2 results in Sect. 8 have used the extended common mask, except the validation analysis at different resolutions (right panel) which for consistency employs the common mask.

Open with DEXTER
In the text
thumbnail Fig. 22

Generalized feature models analysed at max = 2000 (E-modes max = 1500) for the different Planck foreground separation methods, SMICA (blue), SEVEM (red), NILC (green), and Commander (yellow), together with the SSN (SMICA - SEVEM - NILC ) average (black). The left three panels apply to the equilateral feature models, showing, respectively, in the top panel the full feature survey significance at each frequency and phase (temperature only), the maximum significance at each frequency for temperature only (middle), and with polarization (lower). The right three panels apply to the flattened feature models, showing the maximum significance at each frequency for temperature only (top right) and with polarization (middle right), along with significance at each frequency and phase for temperature and polarization (right lower).

Open with DEXTER
In the text
thumbnail Fig. 23

Single field feature model significance with a K2cosωK scaling dependence (Eq. (19); left panels, T-only upper and T + Elower) or with a KsinωK scaling (Eq. (20); right panels). To find the maximum signal, these results have been marginalized over the α-dependent envelope function ranging from α → 0 (no envelope) to the maximum cutoff allowed by the Modal resolution αω = 90.

Open with DEXTER
In the text
thumbnail Fig. 24

Generalized resonance models analysed at max = 2000 (E-modes max = 1500) for the different Planck foreground separation methods, SMICA (blue), SEVEM (red), NILC (green), and Commander (yellow), together with the SSN (SMICASEVEMNILC ) average (black). The upper panels apply to the constant resonance model (Eq. (10)), with T-only (left) and T+E (right), the middle panels give results for the equilateral resonance model (Eq. (13)), and the lower panels for the flattened resonance model (Eq. (14)). Both the equilateral and flattened resonance models produce broad peaks, which are reinforced with polarization (middle and bottom right panels).

Open with DEXTER
In the text
thumbnail Fig. 25

Standard resonance model results for both T-only and T+E across a wide frequency range using the high frequency resonance model estimator. The first and second panels show the signal in both T-only (first panel) and T+E (second panel) across the frequency range 0 <ω < 500. The lower two panels give the results for T-only (third panel) and T+E (fourth panel) in the range 500 <ω < 1100.

Open with DEXTER
In the text
thumbnail Fig. 26

Parity-odd signals (1 + 2 + 3 = odd) of the TTT (top) and EEE (bottom) bispectra (i ≤ 500) recovered from the SMICA maps by means of the Modal decomposition with 101 simple polynomial-based eigenmodes, not including any special functions fitting the CMB tensor-mode features. In the panel for EEE, only the signals larger than = 40 are shown. The TTT and EEE bispectra shown here are rescaled with a constant Sachs-Wolfe weighting and signal-to-noise weighting, respectively.

Open with DEXTER
In the text
thumbnail Fig. 27

Filters used to optimize constraints on local NG, in harmonic space. The temperature filter WM is a smoothed version of the true Wiener filter obtained with realistic models, while the E-mode WM filter is adapted from the temperature one, with a cutoff value at ≃ 800. The formulae for the derivative filters are given in Eq. (78).

Open with DEXTER
In the text
thumbnail Fig. 28

Joint constraint on and obtained with MFs. The contour lines represent 1, 2 and 3σ limits of a 2D-Gaussian distribution. Constraints were obtained with SMICA temperature and polarization E maps.

Open with DEXTER
In the text
thumbnail Fig. 29

68%, 95%, and 99.7% confidence regions in the parameter space , defined by thresholding χ2 as described in the text.

Open with DEXTER
In the text
thumbnail Fig. 30

68%, 95%, and 99.7% confidence regions in the single-field inflation parameter space , obtained from Fig. 29 via the change of variables in Eq. (87).

Open with DEXTER
In the text
thumbnail Fig. 31

68%, 95%, and 99.7% probability contours in the Galileon models for cs and parameters for the branch (tensor spectral index nT< 0).

Open with DEXTER
In the text
thumbnail Fig. 32

68%, 95%, and 99.7% probability contours in the Galileon models for cs and parameters for the branch (blue tensor spectral index nT> 0).

Open with DEXTER
In the text
thumbnail Fig. 33

68% and 95% confidence regions in the plane, with the Lorentz invariant model in Eq. (107) shown as the dashed line.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.