The bispectrum of polarized galactic foregrounds

TL;DR

This paper analyzes the bispectrum of galactic dust and synchrotron emissions across multiple frequencies, revealing strong non-Gaussian features that are crucial for accurate CMB studies and foreground mitigation.

Contribution

It provides the first comprehensive measurement of the galactic foreground bispectrum, including polarization and parity-odd components, across a wide frequency range.

Findings

Strong temperature bispectra in dust and synchrotron emissions, especially in the squeezed limit.

Detection of significant polarized bispectra in galactic dust, peaking in the squeezed configuration.

No residual foreground bispectra found in Planck component separated maps after masking bright sources.

Abstract

Understanding the properties of the galactic emission at millimetre wavelengths is important for studies of the cosmic microwave background (CMB). In this work we explore the bispectrum, the harmonic equivalent of the three point function, from galactic dust and synchrotron emission. We investigate these effects across a broad range of frequencies using the synchrotron dominated S-band Polarization All Sky Survey (SPASS) maps at 2.3 GHz, the Planck satellite maps (30-857 GHz) and dust dominated Infrared Astronomical Satellite (IRAS) maps at 3 THz. We measure bispectra of total intensity fields, T, as well as the gradient, E, and curl modes, B, of the polarization field. We find that the synchrotron and galactic dust have strong temperature bispectra with significant contributions in the squeezed limit, which probes the correlations between two small scale modes with a large scale mode. Additionally, we find that the dust also has strong polarised bispectra that also peak in the squeezed configuration. We explore parity odd bispectra, such as BTT bispectra, and find strong parity odd bispectra for the galactic dust notably in BTT, BTE and BEE configurations. After masking bright sources, we find no evidence for polarised synchrotron bispectra and no evidence for cross bispectra between the dust and synchrotron emission. The strong foreground bispectra discussed here need to be carefully controlled to avoid biasing measurements of primordial non-Gaussianity. Finally we use these bispectra tools to test for residual foregrounds in the component separated Planck maps and find no evidence of residual foregrounds. These tools will be useful for characterizing residual foregrounds in component separated maps, particularly for experiments with less frequency coverage than the Planck satellite.

Equations

$$I_{\mathrm{dust}}\propto {\nu }^{{\beta }_{\mathrm{dust}}}{B}_{\nu }(T_{\mathrm{dust}}),$$

$$I_{\mathrm{sync}}\propto {\nu }^{{\beta }_{\mathrm{sync}}}$$

$$\Delta T(\mathbf{n})=\sum _{\ell ,m}{a}_{\ell ,m}{\mathrm{Y}}_{\ell ,m}(\mathbf{n}).$$

$$(Q\pm iU)(\mathbf{n})=\sum _{\ell ,m}{a}_{\pm 2,\ell ,m}{}_{\pm 2}{Y}_{\ell ,m}(\mathbf{n}),$$

$$a_{E,\ell ,m}=-(a_{2,\ell ,m}+a_{-2,\ell ,m})/2,$$

$$a_{B,\ell ,m}=i(a_{2,\ell ,m}-a_{-2,\ell ,m})/2.$$

$$⟨a_{{\ell }_1,m_1}^X{a}_{{\ell }_2,m_2}^Y{a}_{{\ell }_3,m_3}^Z⟩=B^{X,Y,Z}({\ell }_1,m_1,{\ell }_2,m_2,{\ell }_3,m_3).$$

$$⟨a_{{\ell }_1,m_1}^X{a}_{{\ell }_2,m_2}^Y{a}_{{\ell }_3,m_3}^Z⟩=\sqrt{\frac{(2{\ell }_1+1)(2{\ell }_2+1)(2{\ell }_3+1)}{4\pi }}\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ m_1& m_2& m_3\end{array}\right)\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ 0& 0& 0\end{array}\right)b_{{\ell }_1,{\ell }_2,{\ell }_3}^{X,Y,Z},$$

$$|{\ell }_a-{\ell }_b|\le {\ell }_c\le {\ell }_a+{\ell }_b.$$

$${\mathcal{G}}_{{\ell }_1,{\ell }_2,{\ell }_3}^{m_1,m_2,m_3}$$

$$=\sqrt{\frac{(2{\ell }_1+1)(2{\ell }_2+1)(2{\ell }_3+1)}{4\pi }}\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ m_1& m_2& m_3\end{array}\right)\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ 0& 0& 0\end{array}\right).$$

$${\hat{b}}_{{\ell }_1,{\ell }_2,{\ell }_3}^{X,Y,Z}=\frac{1}{N_{{\ell }_1,{\ell }_2,{\ell }_3}}\sum _m\int {\mathrm{d}}^2\mathbf{n}{\mathrm{Y}}_{{\ell }_1,m_1}(\mathbf{n}){\mathrm{Y}}_{{\ell }_2,m_2}(\mathbf{n}){\mathrm{Y}}_{{\ell }_3,m_3}(\mathbf{n})a_{{\ell }_1,m_1}^X{a}_{{\ell }_2,m_2}^Y{a}_{{\ell }_3,m_3}^Z,$$

$$W_i^X(\mathbf{n})=\sum _{{\ell }_i<\ell \le {\ell }_{i+1}}\sum _m{Y}_{\ell ,m}(\mathbf{n})a_{\ell ,m}^X.$$

$${\hat{b}}_{i,j,k}^{X,Y,Z}=\frac{1}{N_{i,j,k}^{\prime }}{W}_i^X(\mathbf{n})W_j^Y(\mathbf{n})W_k^Z(\mathbf{n}).$$

$$N_{i,j,k}^{\prime }=\sum _{{\ell }_i<{\ell }_1\le {\ell }_{i+1}}\sum _{{\ell }_j<{\ell }_2\le {\ell }_{j+1}}\sum _{{\ell }_k<{\ell }_3\le {\ell }_{k+1}}\frac{(2{\ell }_1+1)(2{\ell }_2+1)(2{\ell }_3+1)}{4\pi }\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ 0& 0& 0\end{array}\right)\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ 0& 0& 0\end{array}\right).$$

$$V_{i,j,k}^{X,Y,Z,X^{\prime },Y^{\prime },Z^{\prime }}=\frac{1}{N_{i,j,k}^{\prime }{N}_{i,j,k}^{\prime }}\sum _{{\ell }_i<{\ell }_1\le {\ell }_{i+1}}\sum _{{\ell }_j<{\ell }_2\le {\ell }_{j+1}}\sum _{{\ell }_k<{\ell }_3\le {\ell }_{k+1}}$$

$$\times C_{{\ell }_1}^{X,X^{\prime }}{C}_{{\ell }_2}^{Y,Y^{\prime }}{C}_{{\ell }_3}^{Z,Z^{\prime }}{g}_{{\ell }_1,{\ell }_2,{\ell }_3},$$

$$h_{{\ell }_1,{\ell }_2,{\ell }_3}^{m_1,m_2,m_3}$$

$$=\int {\mathrm{d}}^2\mathbf{n}\left({}_{-2}{Y}_{{\ell }_1,m_1}(\mathbf{n}){}_1{Y}_{{\ell }_2,m_2}(\mathbf{n}){}_1{Y}_{{\ell }_3,m_3}(\mathbf{n})-{}_2{Y}_{{\ell }_1,m_1}(\mathbf{n}){}_{-1}{Y}_{{\ell }_2,m_2}(\mathbf{n}){}_{-1}{Y}_{{\ell }_3,m_3}(\mathbf{n})\right).$$

$$⟨a_{{\ell }_1,m_1}^X{a}_{{\ell }_2,m_2}^Y{a}_{{\ell }_3,m_3}^Z⟩=\frac{1}{6}\left(h_{{\ell }_1,{\ell }_2,{\ell }_3}^{m_1,m_2,m_3}+h_{{\ell }_3,{\ell }_1,{\ell }_2}^{m_3,m_1,m_2}+h_{{\ell }_2,{\ell }_3{\ell }_1}^{m_2,m_3,m_1}\right)b_{{\ell }_1,{\ell }_2,{\ell }_3}^{odd,X,Y,Z}.$$

$${\hat{b}}_{{\ell }_1,{\ell }_2,{\ell }_3}^{odd,X,Y,Z}=\frac{1}{6N_{{\ell }_1,{\ell }_2,{\ell }_3}^o}\left(h_{{\ell }_1,{\ell }_2,{\ell }_3}^{m_1,m_2,m_3}+h_{{\ell }_3,{\ell }_1,{\ell }_2}^{m_3,m_1,m_2}+h_{{\ell }_2,{\ell }_3,{\ell }_1}^{m_2,m_3,m_1}\right)a_{{\ell }_1,m_1}^X{a}_{{\ell }_2,m_2}^Y{a}_{{\ell }_3,m_3}^Z.$$

$${}_s{W}_i^X=\sum _{{\ell }_i<\ell <{\ell }_{i+1}}\sum _m{}_s{Y}_{\ell ,m}(\mathbf{n})a_{\ell ,m}^X.$$

$${\hat{b}}_{i,j,k}^{odd,X,Y,Z}=\frac{1}{6N_{i,j,k}^{o^{\prime }}}\int {\mathrm{d}}^2\mathbf{n}({}_{-2}{W}_i^X{}_1{W}_j^Y{}_1{W}_k^Z-{}_2{W}_i^X{}_{-1}{W}_j^Y{}_{-1}{W}_k^Z+{}_{-2}{W}_k^X{}_1{W}_i^Y{}_1{W}_j^Z-{}_2{W}_k^X{}_{-1}{W}_i^Y{}_{-1}{W}_j^Z+$$

$${}_{-2}{W}_j^X{}_1{W}_k^Y{}_1{W}_i^Z-{}_2{W}_j^X{}_{-1}{W}_k^Y{}_{-1}{W}_i^Z).$$

$$N_{i,j,k}^{\prime o}=\frac{1}{6}\sum _{{\ell }_i<{\ell }_1\le {\ell }_{i+1}}\sum _{{\ell }_j<{\ell }_2\le {\ell }_{j+1}}\sum _{{\ell }_k<{\ell }_3\le {\ell }_{k+1}}\sum _{m_i}\left(h_{{\ell }_1,{\ell }_2,{\ell }_3}^{m_1,m_2,m_3}+h_{{\ell }_3,{\ell }_1,{\ell }_2}^{m_3,m_1,m_2}+h_{{\ell }_2,{\ell }_3{\ell }_1}^{m_2,m_3,m_1}\right)\left(h_{{\ell }_1,{\ell }_2,{\ell }_3}^{m_1,m_2,m_3}+h_{{\ell }_3,{\ell }_1,{\ell }_2}^{m_3,m_1,m_2}+h_{{\ell }_2,{\ell }_3{\ell }_1}^{m_2,m_3,m_1}\right),$$

$${V^o}_{{\ell }_1,{\ell }_2,{\ell }_3^{\prime }}^{X,Y,Z,X^{\prime },Y^{\prime },Z^{\prime }}=\frac{1}{N_{i,j,k}^{\prime o}{N}_{i,j,k}^{\prime o}}\sum _{{\ell }_i<{\ell }_1\le {\ell }_{i+1}}\sum _{{\ell }_j<{\ell }_2\le {\ell }_{j+1}}\sum _{{\ell }_k<{\ell }_3\le {\ell }_{k+1}}\sum _{m_i}$$

$$\times C_{{\ell }_1}^{X,X^{\prime }}{C}_{{\ell }_2}^{Y,Y^{\prime }}{C}_{{\ell }_3}^{Z,Z^{\prime }}{g}_{{\ell }_1,{\ell }_2,{\ell }_3}.$$

$${\hat{b}}_{i,j,k}^{X,Y,Z}=\frac{1}{N_{i,j,k}^{\prime }}\int {\mathrm{d}}^2\mathbf{n}{W}_i^X(\mathbf{n})W_j^Y(\mathbf{n})W_k^Z(\mathbf{n})-W_i^X(\mathbf{n})⟨{W^{(G)}}_j^Y(\mathbf{n}){W^{(G)}}_k^Z(\mathbf{n})⟩$$

$$-W_j^Y(\mathbf{n})⟨{W^{(G)}}_i^X(\mathbf{n}){W^{(G)}}_k^Z(\mathbf{n})⟩-W_k^Z(\mathbf{n})⟨{W^{(G)}}_i^X(\mathbf{n}){W^{(G)}}_j^Y(\mathbf{n})⟩,$$

$${\hat{b}}_{i,I}^{2D:X,(Y,Z)}=\frac{1}{N}\sum _{2{\ell }_I^2<{\ell }_j^2+{\ell }_k^2\le 2{\ell }_{I+1}^2}{\hat{b}}_{i,j,k}^{X,Y,Z}$$

$${{\chi }^2}_{i,I}^{2D:X,(Y,Z)}=\frac{1}{N}\sum _{2{\ell }_I^2<{\ell }_j^2+{\ell }_k^2\le 2{\ell }_{I+1}^2}{\hat{b}}_{i,j,k}^{X,Y,Z}(V^{X,Y,Z,X,Y,Z}{)}_{i,j,k}^{-1}{\hat{b}}_{i,j,k}^{X,Y,Z},$$

$${\hat{b}}_I^{1D:X,Y,Z}=\frac{1}{N}\sum _{3{\ell }_I^2<{\ell }_i^2+{\ell }_j^2+{\ell }_k^2\le 3{\ell }_{I+1}^2}{\hat{b}}_{i,j,k}^{X,Y,Z},$$

$${{\chi }^2}_I^{1D:X,Y,Z}=\frac{1}{N}\sum _{3{\ell }_I^2<{\ell }_i^2+{\ell }_j^2+{\ell }_k^2\le 3{\ell }_{I+1}}{\hat{b}}_{i,j,k}^{X,Y,Z}(V^{X,Y,Z,X,Y,Z}{)}_{i,j,k}^{-1}{\hat{b}}_{i,j,k}^{X,Y,Z}.$$

$$⟨a^X({\mathbf{\ell }}_1)a^Y({\mathbf{\ell }}_2)a^Z({\mathbf{\ell }}_3)={\delta }^{(2)}({\mathbf{\ell }}_1+{\mathbf{\ell }}_2+{\mathbf{\ell }}_3)b({\ell }_1,{\ell }_2,{\ell }_3),$$

Full text

1 Introduction

One of the main goals of current and up-coming cosmic microwave background (CMB) surveys is to detect the imprint of primordial tensor modes. Primordial tensor modes are theorised to have been generated in the early universe during inflation (Grishchuk, 1975; Starobinskij, 1979; Rubakov et al., 1982) and these tensor modes then leave an imprint on the CMB at the surface of last scattering (Fabbri & Pollock, 1983; Abbott & Wise, 1984). Their contribution to the temperature anisotropies (T mode) and curl-free component of the polarised CMB (E mode) has found to be masked by the dominant scalar modes (Spergel et al., 2007). Instead cosmologist seek to measure the tensor modes by measuring the curl component of the polarised CMB (B mode) (Seljak, 1997; Seljak & Zaldarriaga, 1997). Scalar modes do not produce B mode polarisation and so the detection of the B modes from the surface of last scattering would be strong evidence for primordial tensor modes (Polnarev, 1985; Zaldarriaga & Seljak, 1997; Kamionkowski et al., 1997) .

Whilst the level of the primordial B mode signal is unknown, the bounds from current experiments mean that it is unlikely that primordial signals will be the dominant sky signal at any frequency (Array et al., 2018). Currently there are four known sources of signal which could mask the primordial signal. These foregrounds are: polarised dust emission, polarised synchrotron emission, gravitational lensing and patchy tau. First polarized dust emission arises from asymmetric dust grains that are aligned with the galactic magnetic field. The dust emission can be described by a modified black body

[TABLE]

where $B_{\nu}(T)$ is the Planck function, $T_{\mathrm{dust}}$ is the dust temperature (and the Planck best fit value is $T_{\mathrm{dust}}=19.6K$), and $\beta_{\mathrm{dust}}$ is the dust spectral index, which is found to be $\beta_{\mathrm{dust}}=1.48\pm 0.01$ for temperature and $\beta_{\mathrm{dust}}=1.53\pm 0.02$ for polarization (Planck Collaboration et al., 2018). Dust emision is the dominant foreground at frequencies above $\sim 150$ GHz (Draine, 2004; Draine & Fraisse, 2009). Second, polarised synchrotron emission arises from electrons spiralling in the galactic magnetic field and is dominant at lower frequencies (Kogut et al., 2007). At first order the Galactic synchrotron radiation’s spectral behaviour can be described by a power law

[TABLE]

where $\beta_{\mathrm{sync}}$ is the sky mean spectral index and the current best fit values are $-3.09\pm 0.05$ at 23 GHz (Kogut, 2012) for temperature and $\beta_{s}=-3.13\pm 0.13$ for polarisation (Gold et al., 2011; Planck Collaboration et al., 2018). However, there is strong evidence that the synchrotron spectral index, for both polarised and temperature emission, varies across the sky (Miville-Deschênes et al., 2008; Fuskeland et al., 2014; Jew, 2017; Jew & Grumitt, 2019; Krachmalnicoff et al., 2018) and shows spectral steepening (Kogut, 2012; Dickinson et al., 2018). Third, gravitational lensing B modes are generated as light propagates from the surface of last scattering to the observer. The light is lensed by intervening matter and the lensing shears the E mode signal to B mode signal (Seljak, 1996; Blanchard & Schneider, 1987). The final contribution, patchy tau, arises due to inhomegenous reionization, where spatial variation in the optical depth generates T, E and B mode anisotropies (Weller, 1999; Hu, 2000). However as this effect is small compared to current and upcoming B mode bounds (Roy et al., 2018; Namikawa, 2018), we ignore this effect in this work.

The primordial B mode signal distinguishes itself from these other sources as it has a blackbody spectrum with temperature $2.726\pm 0.010$ K (Mather et al., 1994) and is thought to have highly Gaussian fluctuations. Dust and synchrotron emission have distinctly different spectra and whilst gravitational lensing signals have an identical spectra, they are highly non-Gaussian. Many different approaches to foreground cleaning have been proposed, see Ichiki (2014) for a review of some of these approaches; several of which have successfully been applied to Planck data (Planck Collaboration IX, 2016; Planck Collaboration XII, 2014). Methods to remove lensing, known as delensing, are described in Seljak & Hirata (2004); Sehgal et al. (2017). The properties of galactic foregrounds are less well constrained than the lensing foregrounds and are the focus of this work.

To remove foregrounds with high precision the foregrounds need to be accurately characterised. Recent measurements have started to constrain the spatial and spectral properties of the dust and synchrotron. Planck measurements have characterised the distribution of polarisation fractions, their correlation with intensity measurements and studied how these properties related to the galactic magnetic field (Planck Collaboration Int. XIX, 2015). Further they found that the galactic dust E mode power has roughly a factor of two more power than the B mode and that the dust T-B power-spectrum is non-zero (Planck Collaboration et al., 2018). Work by Planck Collaboration Int. XIX (2015); Rotti & Huffenberger (2016) have shown that the dust polarised emission varies significantly across the sky. Outside the galactic plane polarised synchrotron emission is mainly from large filaments with high polarisation fractions and these have been extensively studied (e.g. Kogut et al., 2007; Planck Collaboration XXV, 2016) including in the context of CMB B mode foregrounds (Vidal et al., 2015). Building on the work of Kogut et al. (2007); Page et al. (2007), Choi & Page (2015) explored the correlation between the synchrotron emission and the dust emission. Most recently Rana et al. (2018) explored the properties of the temperature bispectrum of synchrotron emission at 408 MHz.

Complimenting these measurements there has been extensive theoretical work on the properties of the foregrounds. Due to shocks and anisotropic motions in the interstellar medium, it is expected that the dust properties will be highly non-Gaussian and anisotropic (see e.g. Elmegreen & Scalo, 2004, for a review). Theoretical models (Caldwell et al., 2017) and numerical magnetohydrodynamics (Kritsuk et al., 2017; Vansyngel et al., 2017; Kandel et al., 2018) have explored the physics behind the observed ratio of dust EE to BB power. Work by Kamionkowski & Kovetz (2014) suggested that the anisotropy of the galactic dust could be in the form of local hexadecapolar type deviations. This was further explored in Philcox et al. (2018), where this anisotropy was measured in simulations and the detectability of the local hexadecapole with future surveys was discussed.

Once foreground cleaned CMB maps are obtained tools are need to validate that any remaining signal is not residual foregrounds. The importance of validating cleaned maps was highlighted in recent work (e.g. Madhavacheril & Hill, 2018) that found the Planck foreground cleaned temperature maps contain significant residual thermal Sunyaev Zel’dovich signal. Traditional approaches rely on cross correlations with foreground dominated maps to constrain residual foregrounds. Whilst cross correlations are a very useful tool, alone they may not be sufficient to validate a potential primordial signal, particularly if the power spectra have been used during the foreground removal process or if foregrounds de-correlate with frequency, though there is currently no evidence to show that they do (Sheehy & Slosar, 2018). Kamionkowski & Kovetz (2014) and Rotti & Huffenberger (2016) suggested using measurements of the anisotropy to characterise remaining foregrounds, as the primordial signal is expected to be isotropic but residual foregrounds are not. Recently von Hausegger et al. (2018) explored how skewness and kurtosis measurements can be used to constrain the residual foreground contamination. In Planck Collaboration XVI (2016) various additional methods (skewness, Minkowski functionals and N-point functions) were used to search for deviations from statistical isotropy and homogeneity, which could be indicative of residual foregrounds.

Our work explores how the bispectrum can be used to characterise the foregrounds and test for residual signals, exploiting the fact that the primordial signal is expected to be highly Gaussian. The bispectrum is the harmonic equivalent of the three point function and vanishes for Gaussian signals. The bispectrum is the lowest order non-Gaussian statistic. We first seek to characterise the bispectrum of the galactic foregrounds using data from the Planck satellite. Then we demonstrate how measurements of the bispectrum can be used to test for residual foregrounds in foreground-cleaned maps. This work builds on previous bispectrum measurements of galactic foregrounds (Komatsu et al., 2002; Renzi, 2012), which used bispectrum measurements of foregrounds to avoid biases in primordial non-Gaussianity estimators. Our examination of the bispectrum extends the work in Planck Collaboration XVI (2016) by looking at the full bispectrum (as opposed to two subsets of the three point function), by including polarized maps and by correlating cleaned maps with foreground dominated maps. Understanding the bispectrum from galactic foregrounds, and especially any residual foregrounds in cleaned maps, is very important for constraints on primordial non-Gaussianity. As this work was being prepared for publication, Jung et al. (2018) published their work, in which they examine the temperature bispectrum of galactic foregrounds and show that these foregrounds can significantly bias primordial non-Gaussianity measurements. Our work is highly complementary to theirs, with this work including polarized bispectra and discussing their use in verifying cleaned maps.

We use a binned bispectrum approach (Bucher et al., 2010, 2016). The binned bispectrum is blind approach that requires no theoretical model of the signal and can constrain non-smooth signals. The cost of this broad approach is that it can be less optimal constraints that more targeted approaches. The details of our implementation of the binned bispectrum are described in Section 2. In Section 3 we describe the data sets used in this work and briefly outline our analysis pipeline. Our results are presented in Section 4 and then are discussed, along with our conclusions, in Section 5. In Appendix A we discuss the details of the estimator variance and elaborate on some of the analysis choices used in this work.

2 Bispectrum Estimator

In this section we will briefly overview the binned bispectrum estimator for the parity even and odd cases. For both estimators we must decompose the maps into spherical harmonic components. The measurements of the stokes I component, $\Delta T(\mathbf{n})$, are decomposed using spin-zero harmonics as

[TABLE]

When no mask is applied, the stokes Q and U polarisation components can be decomposed into spin two spherical harmonics

[TABLE]

which can then be further decomposed into two scalar fields

[TABLE]

When a sky mask is applied the decomposition mixes E and B modes together. In this work we use pure E and B estimators, as described in Smith et al. (2007) and Grain et al. (2012), to obtain maps which are free of this leakage.

2.1 Parity Even Estimator

The bispectrum is given by the ensemble average of three spherical harmonic coefficients

[TABLE]

where $B^{X,Y,Z}(\ell_{1},m_{1},\ell_{2},m_{2},\ell_{3},m_{3})$ is the bispectrum between maps $X$,$Y$ and $Z$. Under the assumptions that this signal is statistically homogeneous and isotropic, the bispectrum can be decomposed as (Komatsu & Spergel, 2001)

[TABLE]

where the terms in brackets are Wigner 3j symbols and $b_{\ell_{1},\ell_{2},\ell_{3}}$ is the reduced bispectrum. The Wigner 3j symbol is related to the Clebcsh-Gordan coefficients that describe the coupling of two angular momenta. The Wigner 3j symbol vanishes unless the triangle conditions are satisfied. They require that

[TABLE]

This estimator is zero unless $\ell_{1}+\ell_{2}+\ell_{3}=\mathrm{even}$ and so is only sensitive to even parity configurations. For this reason we will refer to this estimator as the parity-even estimator. The foregrounds are neither isotropic nor homogeneous and by focusing on the reduced bispectrum we lose any anisotropic information (this is also an issue for most power spectra approaches). It is left to future work to explore the anisotropic contributions to the bispectra.

We use a binned bispectrum method similar to that described in Bucher et al. (2016). Using the fact that the Wigner 3j symbols can be evaluated using the Gaunt integral

[TABLE]

we can estimate the reduced bispectrum as

[TABLE]

where $N_{\ell_{1},\ell_{2},\ell_{3}}$ is a normalisation constant. The binned estimator is a simple modification of the above formula. The maps are filtered in harmonic space to contain only modes with $\ell$ within the bin. The filtered map with $\ell$ satisfying $\ell_{i}<\ell\leq\ell_{i+1}$ is denoted as $W^{X}_{i}(\mathbf{n})$ and is given explicitly by

[TABLE]

The binned bispectrum estimator is then given by

[TABLE]

The estimator normalisation is given by

[TABLE]

The Gaussian part of the estimator’s variance is given by

[TABLE]

where $C^{X,Y}$ is the power spectrum between map $X$ and map $Y$ and $g_{\ell_{1},\ell_{2},\ell_{3}},$ is 6 if all its arguments are equal, 2 if any two of the arguments are equal and 1 if none of them are equal. When masks are applied to the data the variance of the estimator is altered by a factor of $f_{\rm{sky}}$.

2.2 Parity Odd Estimator

Planck Collaboration et al. (2018) recently measured a non-zero BT power spectrum for the dust. This is an odd-parity signal and motivated us to consider odd parity bispectra, those which have $\ell_{1}+\ell_{2}+\ell_{3}=\mathrm{odd}$. Odd parity bispectra were originally investigated in the context of parity violating inflationary models and searches for primordial magnetic fields (Shiraishi, 2012; Shiraishi et al., 2013). If the foreground signals are isotropic, homogeneous and parity invariant signals, then we would expect to find only parity even T and E bispectra (Shiraishi et al., 2011; Kamionkowski & Souradeep, 2011). However, even in the absence of parity violating signals, bispectra involving odd numbers of B modes naturally have odd parity (Meerburg et al., 2016).

The reduced bispectrum estimator described in Section 2.1 is only sensitive to bispectra with even parity, those that satisfy $\ell_{1}+\ell_{2}+\ell_{3}=\mathrm{even}$. To measure parity odd bispectra we introduce a second estimator, hereafter called the parity-odd estimator, based on the work of (Shiraishi et al., 2014, 2015). For notational convenience we define the following function

[TABLE]

We define the odd parity reduced bispectrum as

[TABLE]

We note that this definition is not unique (as was discussed in Shiraishi et al., 2014) and different definitions would result in scalings of the reduced bispectrum by factors of $\sim\sqrt{\ell}$. Now we can write down an estimator for the full odd parity bispectrum as

[TABLE]

We symmetrise with respect to $h^{m_{1},m_{2},m_{3}}_{\ell_{1},\ell_{2},\ell_{3}}$ so that our definition of $\hat{b}^{odd,X,X,X}_{\ell_{1},\ell_{2},\ell_{3}}$ is symmetric under the permutation of $\ell_{1}$, $\ell_{2}$ and $\ell_{3}$. Our odd parity estimator vanishes for bispectra that do not satisfy the triangle conditions or have $\ell_{1}+\ell_{2}+\ell_{3}=\mathrm{even}$. In a similar manner to the parity even bispectrum we define a set of filtered maps. The filtered maps are obtained via a spin weighted spherical harmonic transform and are defined by

[TABLE]

Then the binned odd parity bispectrum estimator is given by

[TABLE]

In analogy to the parity even case, the normalization is

[TABLE]

and the variance is

[TABLE]

2.3 Linear Term

The variance of the estimators defined above are only correct in the case that the data is homogenous and isotropic. In the case of a real analysis this condition is broken by the detector noise, masking and the galactic foreground signals. To account for these effects the estimator must be altered to include a linear correction term. As this term has been discussed thoroughly in the literature Creminelli et al. (2006); Yadav et al. (2008); Bucher et al. (2016) here we just summarise the results. The parity even estimator given in Eq. 13 is altered to

[TABLE]

and the parity odd estimator given in Eq. 2.2 is altered in an identical manner.

Computation of the linear term requires a large number of simulations. In this work we explore many different configurations of maps and this would require a very large number of simulations. However, the cuts we use, which are described in Section 3.2, mean that the linear term from the noise and mask cuts has only a small effect on our estimator, decreasing the variance by $\sim 10\%$. This is because the linear term from the mask an anisotropic noise most significantly affect squeezed configurations (Creminelli et al., 2006; Bucher et al., 2016) and these configurations are masked with our low $\ell$ cut. For this reason we do not calculate the linear term for the results discussed below. There are two an important caveats here. Firstly the foregrounds are anisotropic and their potential impact on the linear term is discussed in Section 3.3. Secondly, we tested the importance of the linear term using Planck single frequency sky simulations (Planck Collaboration XII, 2016) and using component separate simulations. The component separation method used was the Spectral Matching Independent Component Analysis SMICA method. The SMICA technique is a component separation method that is capable of minimising the variance of a linear combination of frequency channels in the spatial, harmonic, or hybrid domains. In extends beyond a standard ILC by fitting a model for the frequency map weights (as opposed to using the observed covariance). However in our analysis we also use other component separated maps. The noise in these maps is likely very different from the single frequency and SMICA maps and without noise simulations from the component separation pipelines we are unable to accurately asses the importance of this noise. With that said, all the maps used here passed our null test and we are satisfied that the linear term for our configurations can be ignored even for these caveated cases, though in future work will further investigate this.

3 Data Sets and Pipeline

3.1 Data sets

For this work we used data from the Planck 2015 data release (Planck Collaboration X, 2016; Planck Collaboration IX, 2016). For our fiducial studies of the foregrounds we use the Commander synchrotron map and the 353 GHz map as a tracer of the dust. For a detailed description of the Commander component separation method we refer the reader to Eriksen et al. (2006, 2008) and here we very briefly overview the method. Commander is a Bayesian, parametric map-based method which separates a set of different frequency maps into their components. The map components are modelled by an amplitude for each pixel and a spectral parameterisation. The amplitudes and spectral indexes in each pixel are then obtained by an MCMC Gibbs sampling algorithm. The Commander method uses temperature data from the Planck Satellite, WMAP and a reprocessed version of the 408 MHz radio continuum all-sky map HASLAM map and only Planck data for the polarization maps (Haslam et al., 1982; Bennett et al., 2013; Remazeilles et al., 2015). The effective beam beam full width half-maximum (FWHM )of the 353 GHz maps is 4.86 arcmin.

We also investigate the frequency dependence of the dust by investigating correlations between the Planck 353 GHz and the Planck Collaboration’s reprocessing of the IRAS and IRIS maps (Planck Collaboration XI, 2014), hereafter we refer to this as the IRIS map. This temperature map is obtain from combining the processing of the Infrared Astronomical Satellite (IRAS) and COBE/DIRBE from Schlegel et al. (1998) with a reprocessing of this data-set by Miville-Deschênes & Lagache (2005). This combined map provides the best measure of the galactic dust at $100\mu$m by utilizing more accurate large scale features from the Schlegel et al. (1998) with the high resolution small scale modes from Miville-Deschênes & Lagache (2005). The effective beam FWHM of the IRIS map is 4.3 arcmin.

For our synchrotron analysis we used the data from the S-band Polarization All Sky Survey (S-PASS) (Carretti et al., 2019). The SPASS observed the southern sky at Dec$<-1^{\circ}$ at 2.3 GHz in temperature and polarization. This survey has low noise (with a signal-to-noise ratio larger than 3 on around 94$\%$ of the Q and U pixels) and high resolution, with a 8.9 arcmin FWHM, and so provides a matching data-set to the dust study. After masking the galaxy we can only probe a relatively small fraction of the sky and so we complement this analysis with the higher noise and lower resolution Commander maximum likelihood synchrotron maps. The Commander temperature map has a reference frequency of 408 MHz and is smoothed such that the effective beam (FWHM) is 60 arcmin, and the Commander synchrotron polarization maps have reference frequencies of 30GHz and are smoothed to have an effective FWHM of 40 arcmin.

In Section 4.4, where we investigate cleaned CMB maps, we used maps produced by Planck’s four different cleaning methods: Commander, SMICA, Spectral Estimation Via Expectation Minimization (SEVEM) and Needlet Internal Linear Combination (NILC). The Needlet method works in the needlet domain to reconstruct the CMB anisotropies, whilst minimising the variance, and the SEVEM method is a template fitting method that removes foreground templates to construct a CMB map, with the templates constructed from combinations of the individual frequency channels. These four methods approach the task of foreground cleaning from four different directions that range from data orientated, minimal assumption methods, such as NILC and SEVEM, to parametric methods, such as Commander. For the purpose of this work it is sufficient to note that the different assumptions, approaches and spaces (i.e. wavelet, pixel and harmonic space) in which these methods work mean that these methods are a good representation of the different cleaning methods available in the literature. For more detailed information on the foreground separation methods we refer the reader to Planck Collaboration X (2016) for a general overview, Delabrouille et al. (2009) for NILC, Fernández-Cobos et al. (2012) for SEVEM and Cardoso et al. (2008) for SMICA.

3.2 Pipeline

Before estimating the bispectrum from the maps we preform several pre-processing steps. Firstly we apply a galactic sky mask. Our fiducial mask is the $60\%$ sky mask provided by the Planck experiment but in Section 4 we investigate the effect of the sky cut by using the $20\%$, $40\%$ and $60\%$ Planck masks (Planck Collaboration I, 2016). The mask is apodized with a $3$ degree Gaussian to reduce the effect of mode coupling. We apply the mask as without it our results would be entirely dominated by the bright signal from the galactic plane and because we are interested in characterizing the foregrounds in sky regions of interest to cosmology. Next we apply the Planck point source mask and a second galactic mask that masks the brightest $2\%$ of the sky within in the 60% sky mask. We discuss the motivation behind this mask in Section 3.3. We then iteratively infill these masks to apodize them. Gruetjen et al. (2015) showed that infilling point sources allowed nearly optimal (i.e. minimum variance) efficient bispectrum estimation, see Bucher et al. (2016) for a more detailed discussion of the effect of infilling on the binned bispectrum. We then use the libsharp (Reinecke & Seljebotn, 2013) and HEALPIX (Górski et al., 2005) packages combined with the pure E/B decomposition method from Grain et al. (2012) to obtain the spherical harmonic coefficients. Finally as the maps released as part of the Planck 2015 data release have had the modes with $\ell<30$ in polarization masked due to residual systematics (Planck Collaboration VII, 2016), we filter our maps for both polarization and temperature to remove all modes with $\ell<40$.

The error bars on our measurements are computed using Gaussian simulations that have processed with the same pipeline, including same masks and cuts. These error-bars show good agreement with error-bars computed using Eqs. 2.1 and 2.2. We validated this procedure with two tests; firstly we applied our method to a set of end-to-end noise simulations from the Planck satellite, secondly we examined bispectra constructed from the difference of the half ring maps; these maps should have the same variance properties as the data without any bispectrum signal.

3.3 Signal Inhomogeneity and non-Gaussianity

We have made three suboptimal analysis choices: our low $\ell$ cut in temperature, the second dust mask described in Section 3.2 and choosing to use the 353 GHz maps rather than the Commander dust map, which should be a cleaner map of the dust. These choices were necessary as without them we found we were unable to accurately model the estimator’s variance. The dust is highly anisotropic and highly non-Gaussian and we found without these cuts the estimator variance was larger than the Gaussian prediction. There are two possible causes, firstly as the dust is anisotropic it could contribute to the linear term, see Section 2.3. Currently we are not able to calculate the dust’s contributions to the linear term as this requires simulating the anisotropic distribution of the dust. The second possibility is that, as the dust is highly non-Gaussian, there could be significant contributions to the variance from non-Gaussian contributions to the six point function, which we are currently unable to model. Accurately calculating the linear term and non-Gaussian variance is the subject of on-going work and is necessary to obtain optimal results. The 2% dust mask was chosen as it was found that with this mask, and the low $\ell$ cut, we can model the variance of our estimator, using the empirically measured power spectra, to within $\sim 10\%$. We examine this in more detail, and present a potential alternative to these cuts, in Appendix A. Similarly we found that the high level of non-Gaussianity in the synchrotron temperature map results in a significantly larger estimator variance than expected. We found this was not easily suppressed via masking or $\ell$ cuts and so we opt to only use configurations involving one temperature leg (with the exception of the TTT bispectrum) as these configurations can be well modelled.

4 Foreground Bispectra

The bispectrum is a three dimensional function and for visualisation we use two projections of the data for ease of visualisation. The first is the two dimensional bispectrum defined as

[TABLE]

and its associated two dimensional chi-squared

[TABLE]

where $\ell_{i}$ is the $\ell$ of the centre of the ith bin and N is the number of modes in each bin. In our notation the capital subscript indices indicate the scale of the averaged leg and the superscript indices in brackets $(Y,Z)$ indicate which fields have been average over. The two dimensional bispectrum enables us to explore some of the shape dependence of the bispectra. In this work we will consider many different bispectra and to compress the data further we use a second projection, to a one dimensional bispectrum, to allow an overview of the many bispectra. The one dimensional bispectrum defined as

[TABLE]

and the one dimensional chi-squared

[TABLE]

While the 1-D bispectrum is useful for visualisation purposes, it can obscure the signal as the variance for each bin is set by the noisiest bin included in the sum. The one dimensional $\chi^{2}$ is based on the Gaussian covariance and assumes that the bispectra bins are independent. The application of a mask introduces mode coupling between the bispectra bins violating the independence assumptions, however we choose the width of our bins such that the bin to bin coupling constrained to be $<5\%$ such that this assumption is approximately true. For the largest mask, $f_{\mathrm{sky}}=60\%$, we use bin widths starting from $\Delta\ell=35$ and for the smallest mask, $f_{\mathrm{sky}}=20\%$, we use bin widths starting from $\Delta\ell=45$. Note that we deconvolve the beams for all the bispectra presented in this work.

Before presenting our results, we first present a simple application of the 1-D statistics presented above in order to aid the interpretation of our results and as a simple verification of our pipieline. We apply these 1-D statistics to a set of 100 Gaussian realizations of the CMB and 100 maps with Gaussian CMB and with Poisson point sources. In Figure 2 we plot the average 1D bispectrum from the 100 Poisson maps and we plot the theoretical expectation (which is just a constant). In Figure 2 we plot the average 1D chi-squared for the 100 Poisson maps and for the 100 Gaussian maps. As expected the Gaussian maps $\chi^{2}$ is consistent with 1 across all the bins, which says there is no evidence of non-Gaussianity. This is not the case for the Poisson maps which deviate significantly from 1.

4.1 Dust Bispectra

4.1.1 353 GHz Bispectra

It is useful to first review the 353 GHz power spectrum. In Figure 3 we plot TT, EE and BB pseudo-Cl power-spectra for the 353 GHz maps. We see that the BB mode power is a factor of $\sim 2$ less than the EE power and than the power spectra are well described by power laws $C_{\ell}\propto\ell^{\alpha}$ with $\alpha\sim-2.5$ for all three spectra. Now we move to the bispectrum results. In Figure 4a we present the 1-D parity even binned bispectrum of the *Planck *353 GHz maps. In this projection, we see a strong signal for several bispectra combinations, notably TTT, ETT, BBT, BBE, EET and EEE. The remaining combinations seem to show no clear signal. By examining the 1-D parity even chi-squared, shown in Figure 4b, we see that the other configurations are consistent with zero non-Gaussianity.

Having observed non-zero bispectrum it is necessary to verify that is a sky signal. To do this we use the difference of the two-half ring data splits Planck Collaboration VIII (2016). This combination should have no signal and only noise. In Figure 6 we plot the 1-D chi-squared for the noise splits. We see no strong evidence of non-Gaussianity in any of the splits. We also have analysed the Planck 353 GHz noise simulations and found that the levels of non-Gaussianity in these maps are consistent with zero. Secondly we wish to verify this is galactic dust and not another source of non-Gaussianity, such as extra-galactic point sources (Argüeso et al., 2003; Lacasa et al., 2014; Crawford et al., 2014; Coulton et al., 2017). To do this we plot the signal as a function of sky fraction. The masks are constructed from the 353 GHz maps such that the cleanest parts of the sky are selected, thus as we reduce the used fraction of sky we should see a reduction in the size of the bispectrum signal. In Figure 6 we plot the BBT and EET bispectrum signals as a function of the sky cut and find that as the fraction of the sky is reduced the signals are reduced in amplitude.

Having verified that this the signal from galactic dust we now explore the shape dependence for TTT, TBB and EEE. In Figure 7 we present the 2D bispectra for these configurations as well as the corresponding 2D chi-squared. Generally, we find that the dust bispectrum is strongest in the squeezed and folded configurations. It should be noted that the noise in the polarized maps limits our ability to probe these bispectra to significantly higher $\ell$ as well as reduces the configurations that we can probe (equilateral configurations at similar scales to squeezed configurations will be significantly noisier).

In Planck Collaboration Int. XXX (2016) it was found that the dust EE power was twice that of the BB power. In Figure 6 we find that the amplitude of the EET bispectrum is greater than the BBT bispectrum for all sky cuts by a factor of $\sim 2$. Planck found evidence for non zero BT power-spectrum signal, a parity invariance violating signal, thus it is of interest to apply our parity even bispectra estimator to the configurations BTT, BTE, BEE and BBB, which should be zero in the case of no parity violation. We find no evidence for non-zero BTT, BTE,BEE or BBB bispectra with our parity even estimator, thus, for these sky cuts, there is no evidence of parity violation at the bispectrum level. We then consider the parity odd 353 GHz bispectrum and the resulting 1-D chi-squared are shown in Figure 9. We find measurable levels of parity odd non-Gaussianity for many of the bispectrum configurations, most significantly the BET, BEE and BTT. We find no evidence for parity violating non-Gaussianity in our parity-odd estimator measurements, which would appear as non-zero TTT, TTE etc bispectra. In Figure 10 we explore the shape dependence of these parity odd bispectra. These bispectra are also strongly peaked in the squeezed limit and show a strongly negative signal.

The dust power spectra are well characterised by a power laws ($C_{\ell}\propto\ell^{\alpha}$ with $\alpha\sim-2.5$ Planck Collaboration et al., 2018) and we explore the scale dependence of the bispectrum in Figure 9 where we plot the squeezed slice of three parity odd dust bispectra. The 1-D bispectrum plotted elsewhere merges many different configurations and can hide the scale dependence of the bispectrum, and thus for a cleaner probe of the scale dependence we plot the squeezed slice. We see that the BTT, BET and BEE bispectra are roughly described by $b(\ell,\ell,\ell)\propto\ell^{\alpha}$ where $\alpha\sim-2$, higher signal to noise measurements are necessary to obtain a more precise relation.

Beyond characterizing these bispectra, we would like to understand the features that cause them. Fully understanding their origin is a complex task that should be explored with simulations and so here we only briefly explore some of the possible implications of these observed bispectra. The filamentary structure seen in the temperature and polarization maps likely leads to the folded and equilateral bispectrum shapes seen in TTT and EEE as these bispectrum configurations corresponds to elongated over-dense regions. As was discussed in Planck Collaboration Int. XXX (2016), the amplitude of the dust polarisation power spectrum at $\ell=80$ is correlated with the dust intensity. This should result in squeezed TEE and TBB bispectra and we observe a strong squeezed signal for these configurations. Physically this arises as lines of sight with high levels of dust intensity also have high levels of polarization. We observe a strong positive squeezed ETT bispectrum, which implies that the small scale temperature power is correlated with the large scale E mode power. Finally we also observe a strong negative parity odd BTT (and a weaker but still negative BTE) bispectrum. To understand these parity odd bispectra we consider these bispectra in the flat-sky limit, where building an intuition for these bispectra is simpler. In the flat-sky limit, assuming isotropy, homogeneity and that the underlying physics is parity conserving, bispectra between temperature and E mode fields can be written as

[TABLE]

where $X,Y,Z\in\{T,E\}$ and the reduced bispectrum $b(\ell_{1},\ell_{2},\ell_{3})$ which contains the physical information, depends only on the magnitude of the Fourier modes. The parity odd bispectrum between one B mode and two T/E fields has the following form Meerburg et al. (2016)

[TABLE]

Unlike the T and E mode only bispectra, bispectra involving B fields depend on the cross product between the B field and the T / E mode field. This means these bispectra explicitly probe correlations between B modes and perpendicular T/E modes. For the BTT (as well as BTE and BEE bispectra), this implies that there is a significant correlation between the variation of the small scale T/E power with the gradient of the perpendicular B modes.

4.1.2 Frequency Evolution

We then examine the frequency evolution of the dust bispectrum by investigating the IRIS map, described in Section 3, along with correlations with the 353 GHz map. We restrict this analysis to temperature as the IRIS map is a temperature only map. We find that the IRIS map has a strong bispectrum signal and in Figures 11a and 11b we examine the bispectrum. We find a qualitatively similar features in the IRIS bispectrum and 353 GHz bispectrum; a signal with a strong squeezed limit and similar scale dependence in the equilateral configuration.

To quantitatively examine the degree of correlation we introduce a bispectrum cross-correlation coefficient as

[TABLE]

This can be used on the full bispectrum or on the 1D or 2D quantities. In order to evaluate the errors on this quantity we propagate the errors using the method described in Appendix B of Villaescusa-Navarro et al. (2018). In Figure 12 we plot the bispectrum correlation coefficient for three different sky masks. We see a generally strong correlation between the two maps. This correlation is not perfect as there is spatial variation in the dust temperature spectral index (see e.g. Planck Collaboration XI, 2014) This results in a de-correlation of the bispectra, as the relative brightness of regions changes as a function of frequency. In Figure 12 we see there is a weak evidence for variation in the correlation coefficient as a function of sky, with the smallest sky region seeming to show more de-correlation than the largest regions. This increased de-correlation is likely due to the increasing importance of anisotropies in the cosmic infrared background (CIB) in low dust regions. The CIB de-correlates with frequency as this is composed of the emission from galaxies with similar rest frame spectra, but located at different redshifts and thus having different observed frequency dependencies. Planck Collaboration Int. XVII (2014) showed that the CIB was only $29\%$ correlated between 353 GHz and 100 $\mu$m. The impact of the CIB becomes more important at smaller scales and perhaps explains the slight decrease in the correlation coefficient at smaller scales.

4.2 Synchrotron Bispectra

4.2.1 SPASS Synchrotron Bispectra

To analyse the SPASS data we make four modifications to the pipeline described in Section 3.2. Firstly we alter the Planck sky mask to be restricted to the region observed by SPASS. Secondly we introduce a mask to mask Centaurus A, the * Large Magellanic Cloud* and *Fornax A * as described in Krachmalnicoff et al. (2018). Next we restrict our multipole range to $40\leq\ell\leq 200$. This is motivated by the results from Krachmalnicoff et al. (2018) that the maps are dominated by point sources above $\ell>200$; this is also seen in our analyses. The high significance of the measurements of these point sources means they dominate the bispectrum and have significant higher order correlations, that contribute to the bispectrum errors via the connected four and six point functions.

The final change is to mask the point sources detected in the SPASS polarization maps. Without masking we see strong evidence of all the bispectra configurations, including parity violating bispectra such as BBB. These bispectra are all dominated by point sources. Uncorrelated polarised point sources should have no bispectra in configurations with an odd number of polarization fields, such as TTE or EEE, as the point sources should cancel each other. However, that is only true for averaged quantities, any single realisation can have a non-zero bispectrum, particularly when the source number counts has an extended high flux tail. For example, this arises if only a small number of bright polarized sources are dominating the bispectrum. In order to study the bispectrum of the diffuse galactic synchrotron emission we need to mask the point sources. For the intensity map this is done using an iterative in-painting as discussed in Section 3.2. For polarization this in-painting procedure leads to mixing of E and B power and has the potential to mix E and B bispectra. To avoid this we just mask these point sources and then apodize them to prevent leakage of power. This apodization is performed with 3 degree Gaussian (as is used for the Galactic mask) and if many point sources were masked this would result a significant fraction of the map being masked. Thus we can only mask the brightest sources found in the Lamee et al. (2016) polarization catalog and map. Even so with the apodized masking of the three extended sources and the bright sources we have a sizeable mask and potentially suppress non-Gaussianity. In future work we will explore more efficient methods to mask polarized point sources without apodization. In temperature we mask all sources with flux greater than 200mJy and which are detected in the maps or in the Meyers et al. (2017) catalog at more than $4\sigma$.

In Figures 13a and 13b we plot the parity even 1D bispectrum and chi-sqared for the SPASS maps. We find strong evidence for a temperature bispecturm which exhibits a weak scale dependence. This likely arises from a large scale bispectrum from the diffuse emission combined with emission from point sources on smaller scales, as is seen in the power-spectrum (Krachmalnicoff et al., 2018). We see no evidence, after masking the brightest sources, for polarized bispectra or cross correlations.

4.2.2 Planck synchrotron maps

The Planck synchrotron maps are provided at lower resolution than the dust maps and so we consider a reduced $\ell$ range ($\ell\leq 400$) for this analysis. In Figure 14a and 14b we present the 2-D parity even binned bispectrum and 2-D chi-squared of the Planck Commander temperature synchrotron map. We find that there is strong evidence for the TTT bispectrum. As was discussed in Section 3.3, we find that the estimator variance for configurations involving multiple synchrotron temperature maps is larger than expected. This is attributed to the large non-Gaussianity and inhomogeneity of the synchrotron emission and the origins of the extra variance are discuss in Appendix A. In light of this, the chi-squared shown in Figure 14b can be taken as measure of how non-Gaussian the synchrotron is, but not a measure of the significance of the detection of the bispectrum. To avoid this increased variance we only considered polarization configurations involving one temperature leg and find no evidence for bispectrum involving the synchrotron E and B maps. We also search for, and find no evidence for any odd-parity non-Gaussianity. Examining the shape of the synchrotron temperature bispectrum we see that it has most significance in the equilateral limit. As was noted in Jung et al. (2018) and Planck Collaboration X (2016), this is caused by the contribution of unresolved point sources and is evidenced by our examination of the higher resolution SPASS maps.

4.2.3 SPASS - Commander cross correlations

Given the two measures of the synchrotron emission we can investigate the evolution of the bispectrum as a function of frequency. Using the bispectrum correlation coefficient defined in Eq. 30 we investigate the degree of correlation between the SPASS and Commander synchrotron maps. To do this we apply the SPASS mask to the Planck synchrotron map and use the same point source masking level of 1Jy - dictated by the Planck resolution. The resulting correlation coefficient is plotted in Figure 16. We see that these two maps are relatively strongly correlated at the bispectrum level. The synchrotron bispectrum has contributions from the diffuse emission, which has a variable spectral index (see e.g. Guzmán, A. E. et al., 2011), and from radio point sources, which can be broadly divided into two populations -flat spectrum and steep spectrum radio sources (see e.g. Massardi et al., 2010). The different spectral indices of these components, combined with the spatial variations, explains the imperfect correlation between these bispectra and is a possible explanation for the difference seen in the two bispectra configurations.

4.3 Dust - Synchrotron Cross Bispectra

In Choi & Page (2015) they studied the correlation between the dust and synchrotron polarisation signals. Physically this correlation arises as both effects are influenced by the galactic magnetic field. Analogously to that we investigated the dust-synchrotron bispectrum using the 353 GHz map and either the Commander synchrotron map or the SPASS maps. We find that the dust-synchrotron bispectrum is consistent with zero for all of the masks. The 1-D chi-squared for the parity even E mode bispectra are shown in Figure 16. We find no evidence for a 353-synchrotron bispectrum. We note that revisiting this correlation with full-sky, SPASS quality-maps at frequencies used in CMB analyses would be interesting. Though, as was described above, a new method for treating polarized point sources will be required.

4.4 Cleaned Map Bispectra

Motivated by the signals seen above, we investigated whether bispectra can be used to search for residual foregrounds. Before considering foreground cleaned maps it is necessary to see if we see a bispectrum in the Planck single frequency maps. The bispectrum in these maps will be significantly reduced when compared to the results in Section 4.1. In Figures 17a and 17b we examined the 1-D chi-squared for cross bispectra between the 353 GHz T, E and B and the B mode maps from the 100 GHz, 143 GHz and 217 GHz channels. Using the strong foregrounds in the 353 GHz maps we look for correlations in the B map. We find strong evidence for the parity even BTB bispectrum and odd parity EEB bispectrum at 143 and 217 GHz, with hints of a signal in the 100 GHz map.

The signal seen in Figures 17a and 17b means that the bispectrum could be used to test for residual foregrounds. In Figures 18a and 18b we plot the cross bispectra between the 353 GHz T,E and B maps with the SMICA B mode map. As can be seen there is no strong evidence of a residual bispectrum. As discussed in Appendix A we have neglected the linear term in this analysis and this results in an underestimation of the variance of $~{}10\%$ and so our error bars are slightly underestimated. In Figures 19a and 19b we explore the residual correlations for the four different component separation methods. We find high levels of consistency between the foreground cleaning methods, with none of the methods showing evidence of residual foregrounds. Similarly we can test for residuals foregrounds in the temperature and E mode polarisation maps and the results are shown in Figures 20a and 20b. We see no evidence for residual foregrounds.

5 Discussion and Conclusions

The results presented in Section 4 show that there are large levels of non-Gaussianity in the galactic foregrounds, particularly the galactic dust. It also suggests that these bispectra could be useful for cross checking residual foregrounds in cleaned temperature maps. Whilst we found no strong evidence for residual foreground non-Gaussianity in the component separated maps, the work presented here is the first step. To further investigate the consistency with the Gaussian expectation we need to more accurately model the signal variance. Ideally we would also like to push to larger sky areas and to lower $\ell$, however further work is required to understand the excess variance described in detail in Appendix A.

The strong dust and synchrotron bispectra observed here need to be accounted for when searching for primordial non-Gaussianity. This is particularly important as these bispectra peak in the squeezed configurations and so could bias measurements of local type non-Gaussianity. The large polarization signals means that, unlike extragalactic terms discuss in Hill (2018), these biases cannot be avoided by solely using polarization data. The lack of evidence for bispectra in the cleaned maps means that the Planck non-Gaussianity results (Planck Collaboration XVII, 2016) should be unaffected by residual Galactic foregrounds (which, for the temperature foregrounds, has been thoroughly examined in Jung et al. (2018)). In future searches for scalar-tensor non-Gaussianity (Meerburg et al., 2016) these contaminants will need to be tightly controlled and the tools discussed here will be very useful for validating those results.

In this paper we have focused on characterising the bispectra signals and discussing the utility of the bispectrum to constrain residual foregrounds, however there is far more information available. Planck Collaboration Int. XIX (2015) examined the physical origins of the polarized dust signal. Bispectrum measurements can also be used to study the physical properties of the dust and complement the results of Planck Collaboration Int. XIX (2015). In particular, in Burkhart et al. (2009) they demonstrate first using simulations, and then in Burkhart et al. (2010) using measurements of the Small Magellanic Cloud, that column density bispectra can be used to constrain Alfvén and sonic Mach number. Thus in the future bispectrum measurements could potentially be used to constrain the magnetic and kinetic properties of the interstellar medium.

6 Acknowledgements

The authors would like to thank Anthony Challinor, Adri Duivenvoorden, Joanna Dunkley, Daan Meerburg and Blake Sherwin for useful discussions. This work has made use of S-band Polarisation All Sky Survey (S-PASS) data. The authors would also like to thank the referee for their useful comments and great suggestions.

Appendix A An examination of the Bispectrum Estimator Variance

In this Appendix we explore more the comments made in Section 3.3. To begin with we review the contributions to the estimator variances before examining the larger than naively expected observed variance. We then discuss our approach to mitigate this problem, before finally presenting an alternative method to alleviate the problem.

A.1 Review of Estimator Variance

In this section we limit the discussion to the parity even estimator, for conciseness, however this discussion generalises trivially to the parity odd estimator. The formulae for the estimator variance given in Section 2 are valid either for a full sky measurement of a homogenous field or for a measurement of an inhomogeneous field with the linear used in the estimator (Babich, 2005; Komatsu et al., 2002). In the case when the linear term is neglected, or cannot be calculated, there are the additional terms, which have the schematic form

[TABLE]

For the homogeneous and isotropic case $\langle a_{\ell_{1},m_{1}}a_{\ell_{2},-m_{2}}\rangle=(-1)^{m_{2}}C_{\ell_{1}}\delta_{\ell_{1},\ell_{2}}\delta_{m_{1},m_{2}}$, the following Wigner 3j property

[TABLE]

means that these terms are proportional to the monopole, which is set to zero. For anisotropic fields it is no longer the case that $\langle a_{\ell_{1},m_{1}}a_{\ell_{2},-m_{2}}\rangle=(-1)^{m_{2}}C_{\ell_{1}}\delta_{\ell_{1},\ell_{2}}\delta_{m_{1},m_{2}}$, there will be off diagonal terms. These new contributions to the variance can be significant (Bucher et al., 2016).

Secondly, when the non-Gaussianty is strong there can be non-Gaussian contributions to the covariance matrix. So far the contributions we have considered are the Gaussian (or disconnected) components of the variance, which can be calculated via Wick’s Theorem. However for non-Gaussian fields there can also be other components (the connected or non-Gaussian contributions), see for example Kayo et al. (2013) for an overview of the non-Gaussian contributions to the bispectrum covariance matrix. These have the schematic form of

[TABLE]

where $T_{\ell_{a},\ell_{b},\ell_{c},\ell_{d}}$ is the connected trispectrum, $B_{\ell_{a},\ell_{b},\ell_{c}}$ is the bispectrum and $S_{\ell_{a},\ell_{b},\ell_{c},\ell_{d},\ell_{e},\ell_{f}}$ is the connected six point function. For a weakly non-Gaussian field these will be subdominant to the connected terms, but that is not necessarily the case for a strongly non-Gaussian field.

A.2 Application to dust bispectra

In Figure 22 we plot the results from a simple null test. We compute the bispectrum, and the 1D chi-squared, between two 545 GHz temperature data maps and a simulated Planck 353 GHz noise map. As there is no correlation between the simulated map and the data we expect a vanishing bispectrum and a 1D chi-squared consistent with 1. For the 545 GHz map we do not see this, instead we see a non zero signal. The covariance used to calculate the 1D chi-squared was obtained from Gaussian simulations and matches the analytical value calculated from Eq. 2.1. As there is no bispectrum signal, the failure of this test implies that the variance is larger than is given by Eq. 2.1. In Figure 22 we also plot the result from using the Commander dust temperature map and find that it also fails this test. There are two possible causes of this, either the linear term is insufficient or there are non-Gaussian contributions to the covariance matrix.

To calculate the linear term we need to be able to evaluate the off-diagonal power spectrum term, $\langle a_{\ell_{1},m_{1}}a_{\ell_{2},m_{2}}\rangle.$ For even modest $\ell$ this matrix is very large and impractical to fully calculate. Instead the linear term is typically calculated through simulations as described in Smith & Zaldarriaga (2011) (or Bucher et al. (2016) for the binned equivalent). In analyses of extragalactic sources the main inhomogeneous terms arise from instrument noise and galactic sky cuts, the extragalactic signals tend to be homogeneous. As these sources of inhomogeneity can be easily simulated the linear term can be accurately calculated. In our case we have noise simulations of the 353 GHz maps and simulations of the 545 GHz noise. However the galactic signal itself is very anisotropic and the anisotropy of this signal is not included in our calculations as it is very difficult to simulate. This means the residual signal seen in Figure 22 could arise as we are missing the terms given in Eq. A.1 for the dust signal itself.

There is a second possible explanation: as can be seen in Figure 4b the dust is very non-Gaussian. It would not be surprising then for there to be significant non zero higher point correlation functions. These would result in extra contributions to the variance of the form given in Eq. A.1. It is very challenging to calculate these contributions to the covariance matrix.

A.3 Mitigation strategies

The approach we used to overcome this was to use the Planck 353 GHz maps and to mask the data. The Planck 353 GHz map has a slightly lower signal to noise than the 545 GHz or Commander dust map and so the contribution of the non-Gaussian variance terms should be suppressed. By applying a real space mask (as described in Section 3.2) and an harmonic space cut of low $\ell$ modes we reduce the inhomogeneity and non-Gaussianity. In Figure 22 we plot the same null test as described above for this masked 353 GHz map, and we see that it passes this (and our other) null tests. In fact with these cuts (particularly the low $\ell$ cut) we find that the linear term has only a very small effect on our estimator and so we neglect it in this analysis. This can be seen in Figure 22 where we plot the average 2D chi-squared for 30 simulations of Planck SMICA noise without including the linear term. We find that mean is $\sim 10\%$ higher than the expectation (of 1), which corresponds to a $\sim 10\%$ underestimation of the variance and is sufficiently accurate for this work. We also see that the variance of all the configurations seems to be equally increased and thus we do not expect that neglecting the linear term will bias any of our results.

As masking the signal is counter productive (as we wish to study the signal), we briefly note that an alternative method that will be explored more in future work. Instead of masking the signal we take an approach to reduce the inhomogeneity and higher point functions. As the Commander dust map is positive definite this can be achieved by taking the log of this map. In Figure 22 we plot the same null test as described above for this transformed map, we find that it also passes this null test. In Figure 23 we plot the 1D chi-squared from correlations between the log Commander dust map with the Commander dust E and B maps. We see strong evidence of non-Gaussianity, at a similar level to that seen in using the masked 353 GHz map, suggesting that the log operation has not destroyed all the correlations and retains a high level of information. In future work we will further examine these correlations.