© ESO, 2014

1. Introduction

The measurement and study of high-order moments of the cosmic microwave background (CMB) and the large scale structure (LSS) has emerged as an important field in observational cosmology because it probes deviations from Gaussianity. For example, observing non-Gaussianity (NG) from the CMB or the large scale structure would probe primordial NG dating back to the generation of cosmological perturbations, and this in turn would constrain severely the primordial models generating these fluctuations (e.g. inflation). For example, a NG detection with the sensitivities of today’s experiments would rule out a whole class of inflation models, among which the standard single-field slow roll inflation (Acquaviva et al. 2003; Maldacena 2003; Creminelli & Zaldarriaga 2004). However, the latest measurement by the Planck Collaboration has put stringent constraints on this possibility (Planck Collaboration XXIV 2014).

The most popular form for this primordial non-Gaussianity is the so-called “local” type, parametrised by a factor f_NL, for which the Bardeen potential takes the form $\mathrm{\Phi }\mathrm{\left(}{x}\mathrm{\right)}\mathrm{=}{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\left(}{x}\mathrm{\right)}\mathrm{+}{\mathit{f}}_{\mathrm{NL}}\mathrm{\left(}{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{\right(}{x}{\mathrm{)}}^{\mathrm{2}}\mathrm{-}\mathrm{⟨}{\mathrm{\Phi }}_{\mathrm{G}}\mathrm{(}{x}{\mathrm{)}}^{\mathrm{2}}{\mathrm{⟩}}^{\mathrm{)}}\mathit{,}$(1)where Φ_G is the Gaussian part of the potential. It is on this form of NG that most efforts have focused and that the tightest constraints were obtained. Indeed a fast estimators for f_NL was developed and applied to the CMB data (Komatsu et al. 2005, 2009, 2011; Planck Collaboration XXIV 2014). However, at CMB frequencies, many contaminating signals are present, which typically have a non-Gaussian distribution. While the galactic emission is mainly confined to a defined area, the galactic plane, and may thus be masked, extragalactic point sources are present all over the sky. The brightest sources can be masked but the faintest, below detection limits and confusion levels, can neither be detected nor masked. These unresolved sources are therefore one potential source of secondary non-Gaussianity that needs to be estimated.

Essentially, at CMB frequencies the extragalactic point sources split into two main populations defined with respect to the characteristic frequency domain they emit at, namely radio and infrared (IR). We thus consider radio loud sources with strong magnetic fields that yield synchrotron or free-free emission in the radio domain and IR galaxies with thermal emission from dust heated by young stars in the IR domain.

At the CMB frequencies, radio galaxies can be considered as randomly distributed on the sky (Toffolatti et al. 1998; González-Nuevo et al. 2005). They thus have constant polyspectra, in particular their power spectrum and bispectrum. They are hence characterised entirely by their number counts, with, say, the latest model by Tucci et al. (2011). In this context, their bispectrum can be written as ${\mathit{b}}_{\mathrm{PS}}\mathrm{=}\mathrm{\int }{\mathit{S}}^{\mathrm{3}}\hspace{0.17em}\frac{\mathrm{d}\mathit{n}}{\mathrm{d}\mathit{S}}\hspace{0.17em}\mathrm{d}\mathit{S,}$(2)where S is the source flux density in Jy. (b_PS may then be converted to CMB temperature elevation through Planck’s law, but is usually quoted in dimensionless units ΔT/T in the literature).

Unlike radio galaxies, the cosmic infrared background (CIB), produced by IR galaxies, is strongly clustered. Not only are galaxy number counts important for modelling this case, but also the distribution of the CIB anisotropies. Models for the IR galaxies number counts (Lapi et al. 2011; Somerville et al. 2012; Béthermin et al. 2012) and their power spectrum (Pénin et al. 2012; Xia et al. 2012; Thacker et al. 2013) have recently been developed. Furthermore, the most recent CMB experiments, such as Planck (Planck Collaboration I 2014; Planck Collaboration I 2011), the Atacama Cosmology Telescope (Das et al. 2011) and the South Pole Telescope (Keisler et al. 2011), that are probing high frequencies have permitted the detection of CIB on several scales (e.g. Hall et al. 2010; Planck Collaboration XVIII 2011; Viero et al. 2013; Planck Collaboration XXX 2014). At present, extragalactic point sources are thus well characterised at the 1- and 2-point levels, but the study of their higher orders is still in its early stages. Theoretical predictions were performed by Argüeso et al. (2003) and Lacasa et al. (2012), among others. The estimate of contamination of the primordial NG parameter f_NL by the point sources was performed for example by Babich & Pierpaoli (2008), Lacasa et al. (2012) in temperature, and by Curto et al. (2013) accounting for polarisation. They show that radio and IR populations yield different biases since their bispectra have different shapes. Statistical methods of estimating the non-Gaussianity from radio point sources, based on wavelets (Argüeso et al. 2006) or on the bispectrum (Komatsu et al. 2005), have been developed, and first measurements were made on WMAP data (Komatsu et al. 2009, 2011). Recently, detection of the bispectrum of the CIB anisotropies was performed by Planck Collaboration XXX (2014) and Crawford et al. (2014).

The purpose of this article is to propose a fast estimator, taking advantage of a prescription developed in Lacasa et al. (2012) for measuring the amplitude of the non-Gaussianity (bispectrum) of clustered point sources. This estimator can be used in non-Gaussianity studies to disentangle the primordial and foreground contributions, and among foregrounds it can be used to disentangle the contributions from clustered and unclustered point sources, which may help for putting constraints on their respective models.

The article is organised as follows. In Sect. 2 we briefly recall the prescription for the bispectrum of IR sources and introduce the estimator for its amplitude. In Sect. 3 we tackle the problems of partial sky coverage and anisotropic noise. In Sect. 4 we present the expected signal-to-noise ratio (S/N) of detection at several frequencies. In Sect. 5 we show how to combine our new estimator with the ones used for the CMB and radio galaxies to produce a joint constraint of non-Gaussianity, and we conclude in Sect. 6.

2. Bispectrum of clustered sources

2.1. Bispectrum

Given a full-sky map of the temperature fluctuations ΔT(n) of some signal in direction n, it can be decomposed in the spherical harmonic basis ${\mathit{a}}_{\mathit{\ell m}}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}{n}{\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{\left(}{n}\mathrm{\right)}\mathrm{\Delta }\mathit{T}\mathrm{\left(}{n}\mathrm{\right)}$(3)with the usual orthonormal spherical harmonics Y_ℓm$\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}{n}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}{n}\mathrm{\right)}{\mathit{Y}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{\left(}{n}\mathrm{\right)}\hspace{0.17em}\mathrm{=}\hspace{0.17em}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\mathit{.}$(4)The power spectrum C_ℓ of the signal is the 2-point correlation function in harmonic space given as $\mathrm{⟨}{\mathit{a}}_{\mathit{\ell m}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{⟩}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}\hspace{0.17em}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}\hspace{0.17em}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\mathit{.}$(5)The bispectrum b_ℓ1ℓ2ℓ3 is the three-point correlation function in harmonic space: $\mathrm{⟨}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathrm{⟩}\mathrm{=}{\mathit{G}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{{\mathit{m}}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{3}}}\mathrm{\times }{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}$(6)with the Gaunt coefficient $\begin{array}{ccc}{\mathit{G}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}& \mathrm{=}& \int {\mathrm{d}}^{\mathrm{2}}{n}\hspace{0.17em}{\mathit{Y}}_{\mathrm{123}}\mathrm{\left(}{n}\mathrm{\right)}\\ & \mathrm{=}& \sqrt{\frac{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}{\mathrm{)}}_{\mathrm{123}}}{\mathrm{4}\mathit{\pi }}}\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\end{array}$where Y_i = Y_ℓimi. In the following the subscript 123 denotes the product of the corresponding variables e.g. X₁₂₃ ≡ X₁X₂X₃. G_1,2,3 is zero unless the triplet (ℓ₁,ℓ₂,ℓ₃) follows the triangle inequalities and m₁ + m₂ + m₃ = 0.

The bispectrum estimator is $\mathit{b̂}\mathrm{123}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}_{\mathrm{123}}}\sum _{{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{m}}_{\mathrm{2}}\mathit{,}{\mathit{m}}_{\mathrm{3}}}{\mathit{G}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}{\mathit{a}}_{\mathrm{1}}\hspace{0.17em}{\mathit{a}}_{\mathrm{2}}\hspace{0.17em}\mathit{a}3$(9)with ${\mathit{N}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{=}\frac{\mathrm{(}\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}{\mathrm{)}}_{\mathrm{123}}}{\mathrm{4}\mathit{\pi }}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& {\mathit{\ell }}_{\mathrm{3}}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)2$(10)being the “number of modes” for the (ℓ₁,ℓ₂,ℓ₃) triangle.

In the following, we note in bold the bispectrum that accounts for the beam effect B_ℓ$\begin{array}{ccc}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}\mathrm{=}{\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{3}}}\mathit{.}& & \end{array}$We also note that $\begin{array}{ccc}{\mathrm{𝒞}}_{\mathit{\ell }}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}{\mathit{B}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{noise}}& & \end{array}$is the power spectrum accounting for both the noise and the beam effects.

2.2. Optimal estimator of the amplitude of the bispectrum

The non-Gaussianity from unclustered sources at CMB frequencies is entirely characterised by a single parameter b_PS that can be estimated on a map (Komatsu et al. 2005). As a result their contamination to f_NL can be marginalised over. By contrast, the non-Gaussianity from clustered sources is more complicated and is characterised by the full bispectrum, whose estimation is computationally intensive (with $\mathrm{𝒪}\mathrm{\left(}{\mathit{N}}_{\mathrm{pix}}^{\mathrm{5}\mathit{/}\mathrm{2}}\mathrm{\right)}$ operations).

In this context and following work by Argüeso et al. (2003), Lacasa et al. (2012) proposed a full-sky analytical prescription to compute the bispectrum of extragalactic sources from their number counts and power spectra. In this prescription, the bispectrum arising from a general population of point sources writes as ${\mathit{b}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{PS}}\mathrm{=}\mathit{\alpha }\sqrt{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{3}}}}$(11)with C_ℓ being the power spectrum of the considered point sources and $\mathit{\alpha }\mathrm{=}\frac{\mathrm{\int }{\mathit{S}}^{\mathrm{3}}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\mathrm{d}\mathit{S}}{{{\mathrm{(}}^{\mathrm{\int }}{\mathit{S}}^{\mathrm{2}}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\mathrm{d}{\mathit{S}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}\mathrm{=}\frac{{\mathit{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{white}}}{{\mathrm{(}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{white}}}^{\mathrm{)}}}^{\mathrm{3}\mathit{/}\mathrm{2}}}$(12)where S is the source flux density. The integral runs up to the detection limit of the resolved sources indicated by a flux cut S_cut, $\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}$ are the number counts per steradian, and the superscript “white” denotes the contribution from a white-noise with the same number counts. Equation (11) reduces to the usual formula ${\mathit{b}}_{\mathrm{PS}}\mathrm{=}{}^{\mathrm{\int }}\mathit{S}^{\mathrm{3}}\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{S}}\mathrm{d}\mathit{S}$ in the unclustered case (${\mathit{C}}_{\mathit{\ell }}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{white}}$). In that sense Eq. (11) is a generalisation that includes the clustering effect.

Lacasa et al. (2012) have tested the validity of the prescription against numerical simulations developed by Sehgal et al. (2010). They have shown that the prescription agrees with the full bispectrum measured on the simulated maps.

The level of non-Gaussianity induced by point sources can be estimated by measuring the bispectrum directly, or it can be quantified by the amplitude of their bispectra compared to a template bispectrum. Such an estimator has already been proposed for the unclustered sources (Komatsu et al. 2005). In the following, we construct a similar estimator of the bispectrum amplitude for the clustered sources. To do so, we define the template bispectrum as the prescription given in Eq. (11). In practice, the amplitude of the bispectrum is measured by minimising the χ² of the map bispectrum to the template: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}\mathit{A}\mathrm{\right)}\mathrm{=}\sum _{{\mathit{\ell }}_{\mathrm{1}}\hspace{0.17em}\mathrm{\le }\hspace{0.17em}{\mathit{\ell }}_{\mathrm{2}}\hspace{0.17em}\mathrm{\le }\hspace{0.17em}{\mathit{\ell }}_{\mathrm{3}}}\frac{{\mathrm{(}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}\mathrm{-}\mathit{A}\hspace{0.17em}{{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{clust}}}^{\mathrm{)}}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}\mathrm{\right)}}$(13)with ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}\mathrm{\right)}\mathrm{=}\frac{{\mathrm{𝒞}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{tot}}{\mathrm{𝒞}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{tot}}{\mathrm{𝒞}}_{{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{tot}}}{{\mathit{N}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}}\mathrm{\times }{\mathrm{\Delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}$(14)being the bispectrum variance in the weak non-Gaussianity limit¹. In this expression, ${\mathrm{𝒞}}_{\mathit{\ell }}^{\mathrm{tot}}$ is the map power spectrum, and we have $\begin{array}{ccc}{\mathrm{\Delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}& \mathrm{=}& \mathrm{1}\mathrm{+}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\mathrm{+}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{+}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathrm{+}\mathrm{2}\hspace{0.17em}{\mathit{\delta }}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\\ & \mathrm{=}& \{\begin{array}{c}\\ & \\ & \\ & \end{array}\end{array}$(15)Then the maximum likelihood estimator is $\mathit{Â}\mathrm{=}\sum _{{\mathit{\ell }}_{\mathrm{1}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{2}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{3}}}\frac{{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{clust}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}\mathrm{\right)}}\mathrm{\times }{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{Â}\mathrm{\right)}$(16)with ${\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}\mathit{Â}\mathrm{\right)}\mathrm{=}{\left(\sum _{{\mathit{\ell }}_{\mathrm{1}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{2}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{3}}}\frac{{\mathrm{(}{{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{clust}}}^{\mathrm{)}}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}\mathrm{\right)}}\right)}^{-1}$(17)being a normalisation factor giving the error bar in the assumption of a Gaussian likelihood, which is valid in the weak NG limit when σ²(ℓ₁,ℓ₂,ℓ₃) can be considered independent of A.

The theoretical bispectrum, Eq. (11), is separable²; hence a simpler computation of the numerator of Eq. (16) can be devised, inspired by Komatsu et al. (2005). Indeed, defining the filtered map, $\mathit{F}\mathrm{\left(}{n}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell m}}\frac{{\mathit{\alpha }}^{\mathrm{1}\mathit{/}\mathrm{3}}\hspace{0.17em}\sqrt{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{clust}}}\hspace{0.17em}{\mathit{b}}_{\mathit{\ell }}}{{\mathrm{𝒞}}_{\mathit{\ell }}^{\mathrm{tot}}}\hspace{0.17em}{\mathit{a}}_{\mathit{\ell m}}\hspace{0.17em}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}{n}\mathrm{\right)}$(18)and ${\mathrm{𝒮}}_{\mathrm{clust}}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}{n}\mathit{F}\mathrm{(}{n}{\mathrm{)}}^{\mathrm{3}}\mathit{,}$(19)we have ${\mathrm{𝒮}}_{\mathrm{clust}}\mathrm{=}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\frac{{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{clust}}\hat{{b}}\mathrm{123}}{{\mathrm{𝒞}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{tot}}{\mathrm{𝒞}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{tot}}{\mathrm{𝒞}}_{{\mathit{\ell }}_{\mathrm{3}}}^{\mathrm{tot}}}{\mathit{N}}_{\mathrm{123}}\mathrm{=}\mathrm{6}\mathrm{\times }\sum _{{\mathit{\ell }}_{\mathrm{1}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{2}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{3}}}\frac{{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{clust}}\hat{{b}}\mathrm{123}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}\mathrm{\right)}}\mathrm{·}$(20)Furthermore, it is worth noting that this approach is faster than a full bispectrum analysis, with $\mathrm{𝒪}\mathrm{\left(}{\mathit{N}}_{\mathrm{pix}}^{\mathrm{3}\mathit{/}\mathrm{2}}\mathrm{\right)}$ operations.

3. Masked sky and inhomogeneous noise

In realistic cases of CMB analysis, the statistical isotropy of the signal is broken by inhomogeneous noise, e.g. due to scanning strategy, or by masking large areas of the sky, e.g. those contaminated by galactic emissions. In these cases, the estimator of the bispectrum amplitude for clustered sources is no longer optimal and is biased in a non-trivial way. Nevertheless and similar to the case of f_NL, the bias and the lack of optimality can both be tackled through adapted modifications to the estimator.

3.1. Debiasing

When a covariance matrix of the map(s) can be estimated, Creminelli et al. (2006) have shown that Wiener filtering of the map(s) will debias the non-Gaussian estimator from anisotropic contaminants or noise. Specifically, if C is the estimated covariance matrix in harmonic space, the estimator Eq. (16) can be debiased by applying the change $\frac{{\mathit{a}}_{\mathit{\ell m}}}{{\mathrm{𝒞}}_{\mathit{\ell }}^{\mathrm{tot}}}\mathrm{\to }{\mathit{a}}_{\mathit{\ell m}}^{\mathrm{\prime }}\mathrm{=}{\mathrm{(}{{C}}^{-1}\mathrm{·}{{a}}^{\mathrm{)}}}_{\mathit{\ell m}}\mathrm{=}\sum _{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}{\mathit{C}}_{\mathit{\ell m,}{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{-1}{\mathit{a}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}$(21)in the filtered map Eq. (18), where ${a}\mathrm{=}{{}^{\mathrm{(}}\mathit{a}_{\mathit{\ell m}}{}^{\mathrm{)}}}_{\mathit{\ell }\mathrm{=}\mathrm{2..}{\mathit{\ell }}_{\max }\hspace{0.17em}\mathit{,}\hspace{0.17em}\mathit{m}\mathrm{=}\mathrm{-}\mathit{\ell }\mathit{..}\mathit{\ell }}$ is the vector of harmonic coefficients.

This step is sometimes called C^-1 prefiltering, and an efficient algorithm for it has been designed by Elsner & Wandelt (2013).

3.2. Reducing variance and reaching optimality

When isotropy is broken, the three-point correlation function used to define the bispectrum in Eq. (6) no longer has minimal variance. It must be replaced by the Wick product of the three harmonic coefficients given by $\begin{array}{ccc}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}& \mathrm{\to }& {\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathrm{-}\mathrm{⟨}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{⟩}\hspace{0.17em}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\\ & & \mathrm{-}\mathrm{⟨}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathrm{⟩}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{-}\mathrm{⟨}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}\mathrm{⟩}{\mathit{a}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}\mathit{.}\end{array}$(22)This expression exhibits the same mean as in the isotropic case but it has a lower variance. Furthermore, when applied to the isotropic case Eq. (22) gives the same bispectrum estimator as in Eq. (6) (e.g. for ℓ₁ ≠ ℓ₂ ≠ ℓ₃ all the expectation values vanish).

Only the two-point correlation function is to be considered for the expectation values. Therefore in practice, the latter are obtained from a sufficiently large number of Gaussian realisations with same power spectrum as that of the signal considered.

Inputting the linear corrections given by Eq. (22) into the expression of S_IR in Eq. (20) yields ${\mathit{S}}_{\mathrm{clust}}\mathrm{\to }\begin{array}{c}\mathrm{˜}\\ {\mathit{S}}_{\mathrm{clust}}\end{array}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}{n}\mathit{F}\mathrm{(}{n}{\mathrm{)}}^{\mathrm{3}}\mathrm{-}\mathrm{3}\mathrm{\times }\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}{n}\mathit{F}\mathrm{(}{n}\mathrm{)}\hspace{0.17em}\mathit{G}\mathrm{(}{n}\mathrm{)}$(23)with $\begin{array}{ccc}\mathit{G}\mathrm{\left(}{n}\mathrm{\right)}& \mathrm{=}& \sum _{\mathrm{12}}{\mathit{\alpha }}^{\mathrm{2}\mathit{/}\mathrm{3}}\sqrt{{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{IR}}\hspace{0.17em}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{IR}}}\mathrm{⟨}{\mathit{a}}_{\mathrm{1}}^{\mathrm{\prime }}\hspace{0.17em}{\mathit{a}}_{\mathrm{2}}^{\mathrm{\prime }}{\mathrm{⟩}}_{\mathrm{MC}}{\mathit{Y}}_{\mathrm{1}}\mathrm{\left(}{n}\mathrm{\right)}\hspace{0.17em}{\mathit{Y}}_{\mathrm{2}}\mathrm{\left(}{n}\mathrm{\right)}\\ & \mathrm{=}& \mathrm{⟨}\mathit{F}\mathrm{(}{n}{\mathrm{)}}^{\mathrm{2}}{\mathrm{⟩}}_{\mathrm{MC}}\end{array}$where the brackets ⟨ ⟩ _MC stand for the average of Gaussian simulations.

4. Expected detection significance for the CIB

We now apply our formalism for clustered sources to the case of the CIB. We consider here both lensed and unlensed IR galaxies as a single population of dusty star-forming galaxies.

We present the expected detection significance in terms of the S/N for the CIB bispectrum amplitude estimator. To this end, we make use of the values of α and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}}$ that we measure on the simulations by Sehgal et al. (2010). These publicly available simulations developed for the Atacama Cosmology Telescope experiment cover, in particular, the frequencies 150, 220, 280 and 350 GHz. Furthermore, these values of α and ${\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}}$ were shown to adequately reproduce the CIB bispectrum from 150 to 350 GHz frequencies (Lacasa et al. 2012). The instrumental characteristics (beams and noise levels) used are specified in Table 2.

4.1. Ideal cosmic variance-limited case

We first focus on an ideal case where the signal is made solely of CIB galaxies, i.e. noise and other astrophysical contributions are neglected. We will compute the S/N of the detection of the non-Gaussian signal as defined by the amplitude of the bispectrum. Each triplet ℓ₁ ≤ ℓ₂ ≤ ℓ₃ contains information proportional to its number of configurations: ${\mathit{\sigma }}^{-2}\mathrm{\left(}\mathit{Â}\mathrm{\right)}\mathrm{=}{\mathit{\alpha }}^{\mathrm{2}}\sum _{{\mathit{\ell }}_{\mathrm{1}}\hspace{0.17em}\mathrm{\le }\hspace{0.17em}{\mathit{\ell }}_{\mathrm{2}}\hspace{0.17em}\mathrm{\le }\hspace{0.17em}{\mathit{\ell }}_{\mathrm{3}}}{\mathit{N}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathit{.}$(26)In other terms, the S/N of the detection increases with the total number of bispectrum configurations: $\mathit{S}\mathit{/}\mathit{N}\mathrm{=}\frac{\mathrm{1}}{\mathit{\sigma }\mathrm{\left(}\mathit{Â}\mathrm{\right)}}\mathrm{\propto }\sqrt{{\mathit{N}}_{\mathrm{tot}}\mathrm{\left(}{\mathit{\ell }}_{\max }\mathrm{\right)}}$(27)with ${\mathit{N}}_{\mathrm{tot}}\mathrm{\left(}{\mathit{\ell }}_{\max }\mathrm{\right)}\mathrm{=}\sum _{{\mathit{\ell }}_{\min }\hspace{0.17em}\mathrm{\le }\hspace{0.17em}{\mathit{\ell }}_{\mathrm{1}}\hspace{0.17em}\mathrm{\le }\hspace{0.17em}{\mathit{\ell }}_{\mathrm{2}}\hspace{0.17em}\mathrm{\le }\hspace{0.17em}{\mathit{\ell }}_{\mathrm{3}}\hspace{0.17em}\mathrm{\le }\hspace{0.17em}{\mathit{\ell }}_{\max }}{\mathit{N}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}{\mathit{\ell }}_{\mathrm{3}}}\mathit{.}$(28)For example in the case of full-sky maps of the CIB at frequencies 150, 220, 280, and 350 GHz with angular scales ranging from ℓ_min = 2 to ℓ_max = 2048, the obtained S/N are ~1000, and they do not vary much from frequency to frequency (see first row of Table 1).

4.2. Realistic case (CMB+CIB)

In a real experiment, the CIB map is convolved by the instrumental beam, and the signal is contaminated by noise and other astrophysical/cosmological signals. The most prominent contaminants are the CMB and the galaxy, where the latter can be masked. The brightest point sources are also masked. Hence the S/N is significantly decreased compared to the ideal case. Table 1 summarises the expected S/N with the instrumental characteristics specified in Table 2. We show how the detection significance decreases from the ideal case to the case of a convolved but "perfect" full-sky CIB map, then including CMB contamination but still full-sky, and finally with a 50% sky-fraction mimicking the mask that may be applied to avoid galactic dust emission. We also include the case where 90% of the CMB can be removed by component separation methods. As a matter of fact, component separation methods can efficiently estimate the CMB signal from multifrequency observation with an error in amplitude below 10% (see Delabrouille & Cardoso 2007; Remazeilles et al. 2011; Bobin et al. 2013 and references therein). The estimated CMB map can then be substracted from frequency maps to reveal the weaker CIB signal from unresolved sources as performed for example by Planck Collaboration XVIII (2011) and Planck Collaboration XXX (2014). For the partial-sky coverage cases, we assume that the optimisation and the debiasing described in Sects. 3.1 and 3.2 have been applied with the S/N scaling as ${\mathit{f}}_{\mathrm{SKY}}^{\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}}$.

The S/N of the detection is severely reduced from ideal to more realistic cases especially at lower frequencies. The signal is still detectable at all frequencies when including beam and noise effects. However, when the CMB is included the signal becomes undetectable at 150 GHz. Above this frequency the detection significance is S/N ~ 5σ at 220 GHz and reaches 40 to 300σ at 280 and 350 GHz where the CIB emission is dominant even in presence of CMB contribution.

The removal of the CMB component from the frequency maps allows us to further improve the detection at 220 GHz by increasing the S/N to ~46σ. In this case, the CIB non-Gaussianity at 150 GHz seems marginally detectable. However at this frequency the non-Gaussianity from radio sources becomes a problem (while at higher frequencies IR sources are forecasted to dominate the unresolved source population).

5. Joint NG estimation

When several non-Gaussian signals are present, a joint estimation of their amplitudes that takes their covariances into account is necessary. In the present study we focus on the main extragalactic non-Gaussian signals: unclustered sources (in particular radio sources, but also IR late-type galaxies and, more generally, any population exhibiting constant polyspectra), the CIB, and the CMB, with potential local-type non-Gaussianity. A joint estimation requires minimisation of the chi-square: ${\mathit{\chi }}^{\mathrm{2}}\mathrm{\left(}{\mathit{f}}_{\mathrm{NL}}\mathit{,}{\mathit{b}}_{\mathrm{PS}}\mathit{,}{\mathit{A}}_{\mathrm{CIB}}\mathrm{\right)}\mathrm{=}\sum _{{\mathit{\ell }}_{\mathrm{1}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{2}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{3}}}\frac{{\mathrm{(}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{obs}}\mathrm{-}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{model}}\mathrm{(}{\mathit{f}}_{\mathrm{NL}}\mathit{,}{\mathit{b}}_{\mathrm{PS}}\mathit{,}{\mathit{A}}_{\mathrm{CIB}}{\mathrm{\left)}}^{\mathrm{\right)}}}^{\mathrm{2}}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}\mathrm{\right)}}$(29)with ${{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{model}}\mathrm{\left(}{\mathit{f}}_{\mathrm{NL}}\mathit{,}{\mathit{b}}_{\mathrm{PS}}\mathit{,}{\mathit{A}}_{\mathrm{CIB}}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathrm{NL}}\hspace{0.17em}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{CMB}}\mathrm{+}{\mathit{b}}_{\mathrm{PS}}\hspace{0.17em}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{unclust}}\mathrm{+}{\mathit{A}}_{\mathrm{CIB}}\hspace{0.17em}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathrm{CIB}}$(30)and (involved in the definition of σ²(ℓ₁,ℓ₂,ℓ₃), see Eq. (14)) ${\mathrm{𝒞}}_{\mathit{\ell }}^{\mathrm{tot}}\mathrm{=}\mathrm{(}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CMB}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{unclust}}\mathrm{+}{{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}}}^{\mathrm{)}}{\mathit{B}}_{\mathit{\ell }}^{\mathrm{2}}\mathrm{+}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{noise}}\mathit{.}$(31)We use the standard parameter b_PS for compatibility with the literature (instead of noting A_unclust).

Then, if we define the scalar product between two bispectra b^α and b^β: $\mathrm{⟨}{{b}}^{\mathit{\alpha }}\mathit{,}{{b}}^{\mathit{\beta }}\mathrm{⟩}\mathrm{=}\sum _{{\mathit{\ell }}_{\mathrm{1}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{2}}\mathrm{\le }{\mathit{\ell }}_{\mathrm{3}}}\frac{{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathit{\alpha }}{{b}}_{\mathrm{1}\mathit{,}\mathrm{2}\mathit{,}\mathrm{3}}^{\mathit{\beta }}}{{\mathit{\sigma }}^{\mathrm{2}}\mathrm{\left(}{\mathit{\ell }}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}\mathit{,}{\mathit{\ell }}_{\mathrm{3}}\mathrm{\right)}}\mathrm{·}$(32)Minimising Eq. (29) corresponds to solving the linear system: $\begin{array}{ccc}& \left(\begin{array}{c}\\ & & \\ & & \\ & & \end{array}\right)\mathrm{·}\left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathrm{=}& \\ & \left(\begin{array}{c}\\ \\ \\ \end{array}\right)\mathit{.}& \end{array}$(33)If we introduce-noted with upper tilde-the estimators which consider only one source of non-Gaussianity, e.g. $\begin{array}{c}\mathrm{˜}\\ {\mathit{A}}_{\mathrm{CIB}}\end{array}\mathrm{=}\frac{\mathrm{⟨}{{b}}_{\mathrm{obs}}\mathit{,}{{b}}_{\mathrm{CIB}}\mathrm{⟩}}{\mathrm{⟨}{{b}}_{\mathrm{CIB}}\mathit{,}{{b}}_{\mathrm{CIB}}\mathrm{⟩}}\mathit{,}$(34)then Eq. (33) can be rewritten to define the joint NG estimators – noted with upper hat – as $\left(\begin{array}{c}\\ \mathit{f̂}\mathrm{NL}\\ \mathit{b̂}\mathrm{PS}\\ \mathit{Â}\mathrm{CIB}\end{array}\right)\mathrm{=}{\mathrm{ℳ}}^{-1}\mathrm{·}\left(\begin{array}{c}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{f}}_{\mathrm{NL}}\end{array}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{b}}_{\mathrm{PS}}\end{array}\\ \begin{array}{c}\mathrm{˜}\\ {\mathit{A}}_{\mathrm{CIB}}\end{array}\end{array}\right)$(35)with $\mathrm{ℳ}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{1}& \mathrm{\Delta }{\mathit{f}}_{\mathrm{NL}}^{\mathrm{unclust}}& \mathrm{\Delta }{\mathit{f}}_{\mathrm{NL}}^{\mathrm{CIB}}\\ \mathrm{\Delta }{\mathit{b}}_{\mathrm{PS}}^{\mathrm{CMB}}& \mathrm{1}& \mathrm{\Delta }{\mathit{b}}_{\mathrm{PS}}^{\mathrm{CIB}}\\ \mathrm{\Delta }{\mathit{A}}_{\mathrm{CIB}}^{\mathrm{CMB}}& \mathrm{\Delta }{\mathit{A}}_{\mathrm{CIB}}^{\mathrm{unclust}}& \mathrm{1}\end{array}\right)\mathit{,}$(36)and $\begin{array}{ccc}\mathrm{\Delta }{\mathit{f}}_{\mathrm{NL}}^{\mathrm{unclust}}\mathrm{=}\frac{\mathrm{⟨}{{b}}_{\mathrm{unclust}}\mathit{,}{{b}}_{\mathrm{CMB}}\mathrm{⟩}}{\mathrm{⟨}{{b}}_{\mathrm{CMB}}\mathit{,}{{b}}_{\mathrm{CMB}}\mathrm{⟩}} \mathrm{\Delta }{\mathit{f}}_{\mathrm{NL}}^{\mathrm{CIB}}\mathrm{=}\frac{\mathrm{⟨}{{b}}_{\mathrm{CIB}}\mathit{,}{{b}}_{\mathrm{CMB}}\mathrm{⟩}}{\mathrm{⟨}{{b}}_{\mathrm{CMB}}\mathit{,}{{b}}_{\mathrm{CMB}}\mathrm{⟩}}& & \\ \mathrm{\Delta }{\mathit{b}}_{\mathrm{PS}}^{\mathrm{CMB}}\mathrm{=}\frac{\mathrm{⟨}{{b}}_{\mathrm{CMB}}\mathit{,}{{b}}_{\mathrm{unclust}}\mathrm{⟩}}{\mathrm{⟨}{{b}}_{\mathrm{unclust}}\mathit{,}{{b}}_{\mathrm{unclust}}\mathrm{⟩}} \mathrm{\Delta }{\mathit{b}}_{\mathrm{PS}}^{\mathrm{IR}}\mathrm{=}\frac{\mathrm{⟨}{{b}}_{\mathrm{CIB}}\mathit{,}{{b}}_{\mathrm{unclust}}\mathrm{⟩}}{\mathrm{⟨}{{b}}_{\mathrm{unclust}}\mathit{,}{{b}}_{\mathrm{unclust}}\mathrm{⟩}}& & \\ \mathrm{\Delta }{\mathit{A}}_{\mathrm{CIB}}^{\mathrm{CMB}}\mathrm{=}\frac{\mathrm{⟨}{{b}}_{\mathrm{CMB}}\mathit{,}{{b}}_{\mathrm{CIB}}\mathrm{⟩}}{\mathrm{⟨}{{b}}_{\mathrm{CIB}}\mathit{,}{{b}}_{\mathrm{CIB}}\mathrm{⟩}} \mathrm{\Delta }{\mathit{A}}_{\mathrm{CIB}}^{\mathrm{unclust}}\mathrm{=}\frac{\mathrm{⟨}{{b}}_{\mathrm{unclust}}\mathit{,}{{b}}_{\mathrm{CIB}}\mathrm{⟩}}{\mathrm{⟨}{{b}}_{\mathrm{CIB}}\mathit{,}{{b}}_{\mathrm{CIB}}\mathrm{⟩}}\mathrm{·}& & \end{array}$We note in particular that the estimators $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{f}}_{\mathrm{NL}}\end{array}$ and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{b}}_{\mathrm{PS}}\end{array}$ are the ones defined in Komatsu et al. (2005), while $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{A}}_{\mathrm{CIB}}\end{array}$ is the estimator proposed in the present work.

We compute, for illustration, the mixing matrix at 220 GHz. Neglecting the noise and beam effect and assuming the ERCSC flux cuts (Planck Collaboration VII 2011), we find for ℓ_max = 2048, $\mathrm{ℳ}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{1}& \mathrm{0.124}& \mathrm{0.631}\\ \mathrm{2.08}\mathrm{\times }{\mathrm{10}}^{-5}& \mathrm{1}& \mathrm{1.79}\\ \mathrm{3.16}\mathrm{\times }{\mathrm{10}}^{-5}& \mathrm{0.535}& \mathrm{1}\end{array}\right)$(37)where we renormalised the unclustered sources amplitude to get numbers of order 1; i.e., the mixing matrix must be understood to be applied to ${\mathit{b}}_{\mathrm{PS}}\mathit{/}{\mathit{b}}_{\mathrm{PS}}^{\mathrm{expected}}$, where ${\mathit{b}}_{\mathrm{PS}}^{\mathrm{expected}}\mathrm{=}\mathrm{1.8}\mathrm{\times }{\mathrm{10}}^{-28}$ is the expected unclustered sources bispectrum that we derived from the simulations of Sehgal et al. (2010) masking all sources above the ERCSC flux cut (S_cut = 0.28 Jy).

The matrix is mostly bloc diagonal, with the f_NL part decoupling from the b_PS and A_CIB parts. We have checked that unless b_PS ≫ 10^-27, f_NL estimation can be considered unaffected by extragalactic foregrounds at this frequency. Conversely, estimation of NG from unclustered and clustered sources is unaffected by f_NL. However, we see that the two contributions affect each other heavily, with the b_PS-A_CIB submatrix being poorly conditioned (condition number ~180).

When we include instrumental effects, Table 2, (same ℓ_max and flux cut), we find the mixing matrix: $\mathrm{ℳ}\mathrm{=}\left(\begin{array}{c}\\ \mathrm{1}& \mathrm{0.0461}& \mathrm{0.245}\\ \mathrm{1.21}\mathrm{\times }{\mathrm{10}}^{-4}& \mathrm{1}& \mathrm{2.12}\\ \mathrm{1.31}\mathrm{\times }{\mathrm{10}}^{-4}& \mathrm{0.435}& \mathrm{1}\end{array}\right)\mathit{.}$(38)The relative importance of CMB NG compared to point source NG has increased. Indeed, beam and noise effects mostly decrease the S/N in configurations with high multipoles. This down-weights NG from point sources which is dominant at high ℓ, while CMB NG dominates at lower multipoles. However CMB and point source NG estimation can still be considered decoupled unless b_PS ≫ 10^-27. In the b_PS-A_CIB submatrix, we see that the relative importance of b_PS has decreased compared to A_CIB. Indeed the CIB signal is most important on large angular scales, while the unclustered sources S/N comes mostly from small angular scales that are affected by the beam smoothing and the noise.

The b_PS-A_CIB submatrix is still ill-conditioned when we account for beam and noise effect, and this is due to a high correlation between the CIB and unclustered sources’ bispectra templates³. This correlation can be quantified by the quantity r = cosθ ∈ [ − 1,1 ]: ${\mathit{r}}_{\mathit{\alpha ,\beta }}\mathrm{=}\frac{\mathrm{⟨}{{b}}_{\mathit{\alpha }}\mathit{,}{{b}}_{\mathit{\beta }}\mathrm{⟩}}{\sqrt{\mathrm{⟨}{{b}}_{\mathit{\alpha }}\mathit{,}{{b}}_{\mathit{\alpha }}\mathrm{⟩}\hspace{0.17em}\mathrm{⟨}{{b}}_{\mathit{\beta }}\mathit{,}{{b}}_{\mathit{\beta }}\mathrm{⟩}}}$(39)with α,β = unclust, CIB, or CMB. In other terms θ_α,β is the angle between the vectors b_α,b_β in the vector space of bispectra with the scalar product defined in Eq. (32). At ℓ_max = 2048, with the instrumental configuration of Table 2, the correlation coefficients of the different bispectra are listed in Table 3.

In Fig. 1, we plot these coefficients as a function of the maximum multipole used in the analysis, with the same instrumental specifications.

Fig. 1 Correlation coefficients between the radio, CIB, and CMB bispectra at 220 GHz with 5 arcmin Gaussian beam and 10-8 ΔT/T noise. Upper panel: correlation between IR and CMB in red, between radio and CMB in blue, as a function of ℓmax. Lower panel: correlation between IR and radio.

When analysing only large angular scales (small ℓ_max), the unclustered and CIB bispectrum templates are anti-correlated with the CMB bispectrum. This is because both unclustered and CIB bispectra estimators are positive, while the CMB is negative due to the Sachs-Wolfe effect term (Bucher et al. 2010). However, when going to higher resolution, the first acoustic peak makes the CMB bispectrum positive with high values in some configurations (e.g. equilateral). This makes the correlation between the unclustered and CMB bispectrum positive for ℓ ≥ 168. The CIB/CMB correlation in turn becomes positive at slightly higher multipoles, ℓ = 190, because the CIB bispectrum peaks in the squeezed configurations where the CMB bispectrum is still negative. At larger multipoles, because of the complex pattern of the CMB bispectrum and the changes of sign with acoustic peaks, the correlation coefficients tend asymptotically to zero.

The CIB and unclustered sources bispectrum templates are positively correlated since they are both positive, and this correlation is weak when analysing only the lowest multipoles, but it increases rapidly over 80% with a small “dip” at ℓ_max = 200−300 at the location of the first CMB acoustic peak. Indeed the CMB spectrum is involved in the denominator of Eq. (32), so at ℓ_max = 200−300 it typically down-weights non-squeezed configurations where unclustered sources and CIB bispectra are the most correlated. Nevertheless, the correlation between unclustered sources and CIB bispectra asymptotically tends to one at high angular resolutions, at the same time as the CIB bispectrum flattens at high multipoles.

Fig. 2 1σ and 2σ confidence levels for (ACIB, bPS) in a joint estimation at 220 GHz ℓmax = 2048 with a 5 arcmin Gaussian beam and 10-8 ΔT/T noise.

Finally, in Fig. 2 we forecast the likelihood contours in the (A_CIB, b_PS) plane, which we expect with the aforementioned instrumental specifications. The A_CIB and b_PS estimations are quite degenerate, as expected from the high correlations of their templates, which significantly degrades the constraint that can be put on each of them independently. However, we can still have an order-of-magnitude estimation for both of them. Also, including the prior b_PS ≥ 0 (and also to some extent A_CIB ≥ 0) would improve the constraints. Finally, CMB removal through component separation may also help by decreasing the variance of both estimators (although not their correlation).

6. Conclusions

Statistical methods of estimating the non-Gaussianity from unclustered extragalactic point sources have been developed in the past (Komatsu et al. 2005; Argüeso et al. 2006); however, the subject has not been tackled previously for the case of clustered sources. The prescription proposed in Lacasa et al. (2012) yields a separable form for the bispectrum of clustered sources. This allowed us to build a fast estimator for the clustered sources NG amplitude in this article, à-la-KSW (Komatsu et al. 2005). The estimator can account for a lensed population of sources if the number counts and power spectrum entering as ingredients account for it. We showed how this estimator can be optimised to account for partial sky coverage and inhomogeneous noise. We computed the detection significance that can be expected at several high frequencies dominated by the CIB. Finally we showed how this new estimator can be used with estimators of the unclustered sources and CMB non-Gaussianity to build a joint estimation of sources of non-Gaussianity. We find that the estimation of CMB NG is barely correlated with the estimation of point source NG, so that unless the latter is too important –e.g. if few sources are masked– both estimations can be considered decoupled. In particular, this is the case for the cosmological channels (100–220 GHz) with flux cuts lower or equal to that of the Planck Early Release Compact Source Catalog. On the contrary, for point sources, clustered and unclustered NG have similar shapes especially at high multipoles that make their amplitude estimation degenerate. We show how this hampers the detectability of the CIB NG at 220 GHz. Increasing this detectability can be achieved for example by subtracting the CMB or by including priors on b_PS in the analysis (e.g. positivity and/or upper limits given by the flux cut and radio-source models).

This study opens up the possibility of robust measurements of non-Gaussianity where contaminations are treated in a systematic way. It also opens up the possibility of detecting point source NG, which may allow their respective models to be constrained.

Acknowledgments

The authors thank an anonymous referee for the comments and suggestions. They thank B. Wandelt for interesting discussions. They acknowledge the use of the HEALpix package (Gorski et al. 2005).

References

All Tables

All Figures

	Fig. 1 Correlation coefficients between the radio, CIB, and CMB bispectra at 220 GHz with 5 arcmin Gaussian beam and 10-8 ΔT/T noise. Upper panel: correlation between IR and CMB in red, between radio and CMB in blue, as a function of ℓmax. Lower panel: correlation between IR and radio.
In the text

	Fig. 2 1σ and 2σ confidence levels for (ACIB, bPS) in a joint estimation at 220 GHz ℓmax = 2048 with a 5 arcmin Gaussian beam and 10-8 ΔT/T noise.
In the text