General modal estimation for cross-bispectra

TL;DR

This paper introduces a fast, optimal estimator for measuring angular bispectra between two weakly non-Gaussian fields, applicable to various shapes and domains, demonstrated on Planck CMB data.

Contribution

It presents a versatile, separable modal bispectrum estimator capable of analyzing multiple bispectrum shapes and domains, with implementation for Planck CMB data analysis.

Findings

Effective measurement of angular bispectra from observational data.

Applicable to both factorizable and non-factorizable bispectrum templates.

Used in the official Planck satellite analysis of CMB maps.

Abstract

We describe a fast, optimal estimator for measuring angular bispectra between two correlated weakly non-Gaussian fields ($Y$ and $Z$) from observational datasets, based on a separable modal bispectrum expansion. Our methodology is applicable to (1) any shape of the input theoretical bispectrum templates (factorizable or not), (2) both even and odd $\ell_1 + \ell_2 + \ell_3$ multipole domains and (3) both amplitude ($f_{\rm NL})$ bispectrum estimation and full bispectrum reconstruction, considering either joint estimation of ($YYY$, $ZZZ$, $YYZ$ and $ZZY$) shapes, or independent estimation of auto-bispectra ($YYY$ or $ZZZ$) and cross-bispectra ($YYZ$ or $ZZY$); hence, it has quite high versatility. The methodology described here was implemented and used for the official analysis of temperature and polarization cosmic microwave background maps from the $Planck$ satellite.

Tables

Table 1. Table 1 : Averages of three types of estimators ( 3.19 ), ( 5.10 ) and ( 5.14 ) computed from 50 CMB maps sourced by the parity-even scalar equilateral NG (left two columns) and the parity-odd tensor equilateral one (right two columns). The maps used in the columns “ f NL = A subscript 𝑓 NL 𝐴 f_{\rm NL}=A ” and “ f NL T ​ T ​ E = A superscript subscript 𝑓 NL 𝑇 𝑇 𝐸 𝐴 f_{\rm NL}^{TTE}=A ” are created assuming ( f NL T ​ T ​ T , f NL E ​ E ​ E , f NL T ​ T ​ E , f NL E ​ E ​ T ) = ( A , A , A , A ) superscript subscript 𝑓 NL 𝑇 𝑇 𝑇 superscript subscript 𝑓 NL 𝐸 𝐸 𝐸 superscript subscript 𝑓 NL 𝑇 𝑇 𝐸 superscript subscript 𝑓 NL 𝐸 𝐸 𝑇 𝐴 𝐴 𝐴 𝐴 (f_{\rm NL}^{TTT},f_{\rm NL}^{EEE},f_{\rm NL}^{TTE},f_{\rm NL}^{EET})=(A,A,A,A) and ( 0 , 0 , A , 0 ) 0 0 𝐴 0 (0,0,A,0) , respectively. If the estimators are optimal, the bold numbers are expected to reach A 𝐴 A .

	equilateral NG		tensor NG
	(${\ell }_1+{\ell }_2+{\ell }_3=\mathrm{even}$)		(${\ell }_1+{\ell }_2+{\ell }_3=\mathrm{odd}$)
	$f_{\mathrm{NL}}=300$	$f_{\mathrm{NL}}^{TTE}=300$	$f_{\mathrm{NL}}=2000$	$f_{\mathrm{NL}}^{TTE}=2000$
$⟨ℰ⟩$ (3.19)	292	123	2030	1450
$⟨{ℰ}_{TTT}⟩$ (5.10)	301	1.66	3430	1510
$⟨{ℰ}_{EEE}⟩$ (5.10)	290	$-1.67$	1670	$-234$
$⟨{ℰ}_{TTE}⟩$ (5.10)	290	298	2030	2020
$⟨{ℰ}_{EET}⟩$ (5.10)	297	3.75	2000	$-21.7$
$⟨{ℰ}_{TTE}⟩$ (5.14)	64.3	295	2090	2030

Equations

$$\displaystyle\begin{split}\Braket{a_{\ell_{1}m_{1}}^{Y}a_{\ell_{2}m_{2}}^{Y}a_{\ell_{3}m_{3}}^{Y}}&\equiv f_{\rm NL}\,B_{\ell_{1}\ell_{2}\ell_{3}}^{YYY}\left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right),\\ \Braket{a_{\ell_{1}m_{1}}^{Z}a_{\ell_{2}m_{2}}^{Z}a_{\ell_{3}m_{3}}^{Z}}&\equiv f_{\rm NL}\,B_{\ell_{1}\ell_{2}\ell_{3}}^{ZZZ}\left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right),\\ \Braket{a_{\ell_{1}m_{1}}^{Y}a_{\ell_{2}m_{2}}^{Y}a_{\ell_{3}m_{3}}^{Z}}&\equiv f_{\rm NL}\,B_{\ell_{1}\ell_{2}\ell_{3}}^{YYZ}\left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right),\\ \Braket{a_{\ell_{1}m_{1}}^{Z}a_{\ell_{2}m_{2}}^{Z}a_{\ell_{3}m_{3}}^{Y}}&\equiv f_{\rm NL}\,B_{\ell_{1}\ell_{2}\ell_{3}}^{ZZY}\left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right),\end{split}$$

$$\mathcal{E}=\frac{1}{6F}\sum _{X_i{X}_i^{\prime }}\sum _{{\ell }_i}(-1{)}^{{\ell }_1+{\ell }_2+{\ell }_3}{B}_{{\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3}(C_{{\ell }_1}^{-1}{)}^{X_1{X}_1^{\prime }}(C_{{\ell }_2}^{-1}{)}^{X_2{X}_2^{\prime }}(C_{{\ell }_3}^{-1}{)}^{X_3{X}_3^{\prime }}{\mathcal{B}}_{{\ell }_1{\ell }_2{\ell }_3}^{X_1^{\prime }{X}_2^{\prime }{X}_3^{\prime }},$$

$${\cal B}_{\ell_{1}\ell_{2}\ell_{3}}^{X_{1}^{\prime}X_{2}^{\prime}X_{3}^{\prime}}\equiv\sum_{m_{i}}\left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right)\left[a_{\ell_{1}m_{1}}^{X_{1}^{\prime}}a_{\ell_{2}m_{2}}^{X_{2}^{\prime}}a_{\ell_{3}m_{3}}^{X_{3}^{\prime}}-\left(a_{\ell_{1}m_{1}}^{X_{1}^{\prime}}\Braket{a_{\ell_{2}m_{2}}^{X_{2}^{\prime}}a_{\ell_{3}m_{3}}^{X_{3}^{\prime}}}+2\ {\rm perms}\right)\right],$$

$$F=\frac{1}{6}\sum _{X_i{X}_i^{\prime }}\sum _{{\ell }_i}(-1{)}^{{\ell }_1+{\ell }_2+{\ell }_3}{B}_{{\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3}(C_{{\ell }_1}^{-1}{)}^{X_1{X}_1^{\prime }}(C_{{\ell }_2}^{-1}{)}^{X_2{X}_2^{\prime }}(C_{{\ell }_3}^{-1}{)}^{X_3{X}_3^{\prime }}{B}_{{\ell }_1{\ell }_2{\ell }_3}^{X_1^{\prime }{X}_2^{\prime }{X}_3^{\prime }}.$$

$$\left(\begin{array}[]{c}a_{\ell m}^{Y}\\ a_{\ell m}^{Z}\end{array}\right)=L_{\ell}^{-1}\left(\begin{array}[]{c}g_{\ell m}^{Y}+A_{\ell m}^{YYY}\\ g_{\ell m}^{Z}+3A_{\ell m}^{ZYY}+3A_{\ell m}^{ZZY}+A_{\ell m}^{ZZZ}\end{array}\right),$$

$$L_{\ell }$$

$$A_{{\ell }_1{m}_1}^{X_1{X}_2{X}_3}$$

$${\ell }_1+{\ell }_2+{\ell }_3=\mathrm{even},|{\ell }_1-{\ell }_2|\le {\ell }_3\le |{\ell }_1+{\ell }_2|.$$

$$B_{(e){\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3}\equiv h_{{\ell }_1{\ell }_2{\ell }_3}{b}_{(e){\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3},$$

$$h_{\ell_{1}\ell_{2}\ell_{3}}\equiv\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}\\ 0&0&0\end{array}\right),$$

$$Q_{ijk}^{X_1{X}_2{X}_3}({\ell }_1,{\ell }_2,{\ell }_3)\equiv q_{[i}^{X_1}({\ell }_1)q_j^{X_2}({\ell }_2)q_{k]}^{X_3}({\ell }_3),$$

$$\begin{array}{rl}{q}_{[i}^Y({\ell }_1)q_j^Y({\ell }_2)q_{k]}^Y({\ell }_3)& =\frac{1}{6}\left[q_i^Y({\ell }_1)q_j^Y({\ell }_2)q_k^Y({\ell }_3)+5\mathrm{perms}\mathrm{in}i,j,k\right],\\ q_{[i}^Y({\ell }_1)q_j^Y({\ell }_2)q_{k]}^Z({\ell }_3)& =\frac{1}{2}\left[q_i^Y({\ell }_1)q_j^Y({\ell }_2)+q_j^Y({\ell }_1)q_i^Y({\ell }_2)\right]q_k^Z({\ell }_3).\end{array}$$

$$\frac{h_{{\ell }_1{\ell }_2{\ell }_3}{b}_{(e){\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3}}{\sqrt{C_{{\ell }_1}^{X_1{X}_1}{C}_{{\ell }_2}^{X_2{X}_2}{C}_{{\ell }_3}^{X_3{X}_3}}}.$$

$$\frac{v_{{\ell }_1}{v}_{{\ell }_2}{v}_{{\ell }_3}{b}_{(e){\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3}}{\sqrt{C_{{\ell }_1}^{X_1{X}_1}{C}_{{\ell }_2}^{X_2{X}_2}{C}_{{\ell }_3}^{X_3{X}_3}}}=\sum _{ijk}{\alpha }_{ijk}^{X_1{X}_2{X}_3}{Q}_{ijk}^{X_1{X}_2{X}_3}({\ell }_1,{\ell }_2,{\ell }_3).$$

$${\gamma }_{np}^{X_1{X}_2{X}_3}\equiv {⟨Q_n^{X_1{X}_2{X}_3}({\ell }_1,{\ell }_2,{\ell }_3)Q_p^{X_1{X}_2{X}_3}({\ell }_1,{\ell }_2,{\ell }_3)⟩}_e,$$

$${⟨f_{{\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3{X}_1^{\prime }{X}_2^{\prime }{X}_3^{\prime }}⟩}_e\equiv \sum _{{\ell }_i}{\left(\frac{h_{{\ell }_1{\ell }_2{\ell }_3}}{v_{{\ell }_1}{v}_{{\ell }_2}{v}_{{\ell }_3}}\right)}^2{f}_{{\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3{X}_1^{\prime }{X}_2^{\prime }{X}_3^{\prime }},$$

$${\alpha }_n^{X_1{X}_2{X}_3}=\sum _p{\left({\gamma }^{X_1{X}_2{X}_2}\right)}_{np}^{-1}{⟨\frac{v_{{\ell }_1}{v}_{{\ell }_2}{v}_{{\ell }_3}{b}_{(e){\ell }_1{\ell }_2{\ell }_3}^{X_1{X}_2{X}_3}}{\sqrt{C_{{\ell }_1}^{X_1{X}_1}{C}_{{\ell }_2}^{X_2{X}_2}{C}_{{\ell }_3}^{X_3{X}_3}}}{Q}_p^{X_1{X}_2{X}_3}({\ell }_1,{\ell }_2,{\ell }_3)⟩}_e.$$

$$R_n^{X_1{X}_2{X}_3}({\ell }_1,{\ell }_2,{\ell }_3)=\sum _p{\lambda }_{np}^{X_1{X}_2{X}_3}{Q}_p^{X_1{X}_2{X}_3}({\ell }_1,{\ell }_2,{\ell }_3).$$

$${\left({\gamma }^{X_1{X}_2{X}_3}\right)}^{-1}={\left({\lambda }^{X_1{X}_2{X}_3}\right)}^{\top }{\lambda }^{X_1{X}_2{X}_3}.$$

$${\alpha }_n^{X_1{X}_2{X}_3}=\sum _p{\lambda }_{pn}^{X_1{X}_2{X}_3}{\alpha }_p^{R{X}_1{X}_2{X}_3}.$$

$$\begin{array}{rl}{\gamma }_{np}^{X_1{X}_2{X}_3}& =8{\pi }^2{\int }_{-1}^{1}d\mu {\zeta }_{[i[i^{\prime }}^{X_1}(\mu ){\zeta }_{j{j}^{\prime }}^{X_2}(\mu ){\zeta }_{k]k^{\prime }]}^{X_3}(\mu ),\\ {\zeta }_{i{i}^{\prime }}^X(\mu )& \equiv \sum _{\ell }\frac{2\ell +1}{4\pi }\frac{q_i^X(\ell )q_{i^{\prime }}^X(\ell )}{v_{\ell }^2}{P}_{\ell }(\mu ).\end{array}$$

$$h_{{\ell }_1{\ell }_2{\ell }_3}^2=\frac{(2{\ell }_1+1)(2{\ell }_2+1)(2{\ell }_3+1)}{8\pi }{\int }_{-1}^{1}d\mu P_{{\ell }_1}(\mu )P_{{\ell }_2}(\mu )P_{{\ell }_3}(\mu ),$$

$$h_{\ell_{1}\ell_{2}\ell_{3}}\left(\begin{array}[]{ccc}\ell_{1}&\ell_{2}&\ell_{3}\\ m_{1}&m_{2}&m_{3}\end{array}\right)=\int d^{2}\hat{\bf n}\,Y_{\ell_{1}m_{1}}(\hat{\bf n})Y_{\ell_{2}m_{2}}(\hat{\bf n})Y_{\ell_{3}m_{3}}(\hat{\bf n}).$$

$$M_i^X(\hat{\mathbf{n}})$$

$${\overline{a}}_{\ell m}^X$$

$$\mathcal{E}$$

$$\begin{array}{rl}{\beta }_{\mathrm{cub}ijk}^{X_1{X}_2{X}_3}& =\int d^2\hat{\mathbf{n}}{M}_{[i}^{X_1}(\hat{\mathbf{n}})M_j^{X_2}(\hat{\mathbf{n}})M_{k]}^{X_3}(\hat{\mathbf{n}}),\\ {\beta }_{\mathrm{lin}ijk}^{X_1{X}_2{X}_3}& =\int d^2\hat{\mathbf{n}}[⟨M_{[i}^{X_1}(\hat{\mathbf{n}})M_j^{X_2}(\hat{\mathbf{n}})⟩M_{k]}^{X_3}(\hat{\mathbf{n}})+⟨M_{[j}^{X_2}(\hat{\mathbf{n}})M_k^{X_3}(\hat{\mathbf{n}})⟩M_{i]}^{X_1}(\hat{\mathbf{n}})\\ & +⟨M_{[k}^{X_3}(\hat{\mathbf{n}})M_i^{X_1}(\hat{\mathbf{n}})⟩M_{j]}^{X_2}(\hat{\mathbf{n}})].\end{array}$$

$$\left(\begin{array}[]{c}\overline{a}_{\ell m}^{Y}\\ \overline{a}_{\ell m}^{Z}\end{array}\right)\equiv\left(\begin{array}[]{cc}1&\frac{-C_{\ell}^{YZ}}{\sqrt{C_{\ell}^{YY}C_{\ell}^{ZZ}-(C_{\ell}^{YZ})^{2}}}\\ 0&\frac{\sqrt{C_{\ell}^{YY}C_{\ell}^{ZZ}}}{\sqrt{C_{\ell}^{YY}C_{\ell}^{ZZ}-(C_{\ell}^{YZ})^{2}}}\end{array}\right)\left(\begin{array}[]{c}\hat{a}_{\ell m}^{Y}\\ \hat{a}_{\ell m}^{Z}\end{array}\right)$$

$$F_{YYY}=\frac{1}{6}\sum _{n=0}^{n_{\max }}{\left({\alpha }_n^{RYYY}\right)}^2,$$

$$F=\frac{1}{6}\sum _{X_i{X}_i^{\prime }}\sum _{np}{\alpha }_n^{X_1{X}_2{X}_3}{\tau }_{np}^{X_1{X}_2{X}_3{X}_1^{\prime }{X}_2^{\prime }{X}_3^{\prime }}{\alpha }_p^{X_1^{\prime }{X}_2^{\prime }{X}_3^{\prime }},$$

$${\tau }_{np}^{X_1{X}_2{X}_3{X}_1^{\prime }{X}_2^{\prime }{X}_3^{\prime }}\equiv {⟨Q_n^{X_1{X}_2{X}_3}{Q}_p^{X_1^{\prime }{X}_2^{\prime }{X}_3^{\prime }}\left[\prod _{i=1}^3\sqrt{C_{{\ell }_i}^{X_i{X}_i}}(C_{{\ell }_i}^{-1}{)}^{X_i{X}_i^{\prime }}\sqrt{C_{{\ell }_i}^{X_i^{\prime }{X}_i^{\prime }}}\right]⟩}_e.$$

Full text

1 Introduction

Non-Gaussianity (NG) in cosmic observables is an important indicator of the statistical and physical properties of the cosmological fluctuation field, carrying information both about primordial perturbations and later processes, such as reionization and structure formation.

A typical and useful approach to test cosmological NG is the study of beyond two-point correlation functions. Among them, the bispectrum, i.e., the Fourier transform of the three-point correlation, often plays a special role. For example, if we focus on the primordial fluctuation field generated by inflation, the bispectrum of Cosmic Microwave Background (CMB) anisotropies is generally the most powerful test of models beyond standard single-field slow roll. Stringent constraints on primordial NG have been initially obtained using the auto-bispectrum of the temperature field ($T$) in WMAP and Planck data [1, 2]. These were further refined by including the CMB polarization field in the latest Planck data analysis [3]. Optimal estimation of primordial NG from CMB temperature and polarization data requires constraining jointly two auto-bispectra ($TTT$ and $EEE$) and two mixed ones ($TTE$ and $EET$). There are many other interesting examples, in the literature, of cross-bispectra containing relevant cosmological information, see e.g. refs. [4, 5, 6, 7, 8, 9, 10] for the correlators containing CMB B-mode polarization and refs. [11, 12] for those with CMB spectral distortion anisotropies.

A main purpose of this paper is therefore to establish a general methodology for optimal estimation of angular bispectra composed of two correlated weakly NG fields. In order for it to be as general as possible, as well as to encompass all the examples mentioned in the list above, we want our approach to be able to deal with both parity-odd and parity-even bispectra, eventually written in non-separable form (i.e. angular bispectra that cannot be explicitly written as factorized products of the three multipole moments $\ell_{1}$, $\ell_{2}$ and $\ell_{3}$; these models pose significant computational challenges, as explained in detail in the main text below). We tackle these problems by specifically resorting to a modal estimation approach, in which a general bispectrum is written as a linear combination of factorizable basis templates, and the latter are fit to the data. The work developed here will therefore provide a generalization of the modal bispectrum estimation pipeline developed for auto-bispectra in refs. [13, 14, 15, 16, 17]. Another extension of the modal estimation approach, aimed at joint estimation of temperature and polarization CMB bispectra, is developed and discussed in ref. [18]. There, the author orthogonalizes the input bispectrum vector $(B_{\ell_{1}\ell_{2}\ell_{3}}^{YYY},B_{\ell_{1}\ell_{2}\ell_{3}}^{ZZZ},B_{\ell_{1}\ell_{2}\ell_{3}}^{YYZ},B_{\ell_{1}\ell_{2}\ell_{3}}^{ZZY})$ (where $Y$ and $Z$ are two general weakly NG random fields; in the case of CMB Planck analysis, $Y$ and $Z$ represent the temperature and E-mode polarization fields). After rotation, the problem is reduced to that of estimating four independent auto-bispectra in the new variables.

In this paper we take a different approach, namely we directly decompose the original bispectra, without rotation. As we show in the main text, this approach can also yield a factorized form of the estimator (2.2), allowing for fast bispectrum estimation, although with a slightly worse performance than the orthogonalized method. On the other hand, the approach discussed here allows us to very easily obtain a fast, optimal modal estimator for the (physical, un-rotated) mixed bispectra ($YYZ$ or $ZZY$) alone. This can actually be applied to validation and consistency tests in a joint bispectrum analysis, as well as to the estimation of specific cross-bispectra, such as, for example, CMB $TTB$ bispectra mentioned above, which have recently received significant attention in the literature. The methodology we develop in this paper is therefore different and complementary to the orthogonalized modal estimator of ref. [18]. We note that both methods have already been implemented and used for a large range of applications in the Planck data analysis. They are dubbed, respectively, Modal1 (estimator discussed in this paper) and Modal2 (orthogonalized estimator) methods in the Planck 2015 analysis of primordial NG [3].

This paper is organized as follows. In section 2, we briefly recall some general results in the theory of bispectrum estimation, and report the factorized forms of the optimal angular bispectrum estimator (2.2). In section 3 we describe our modal decomposition methodology for even-parity mixed bispectra [13, 14, 15]. In section 4, we generalize the method to odd-parity mixed bispectra. In section 5, we discuss how our modal approach can be directly used not only to build fast estimators of bispectrum amplitudes ($f_{\rm NL}$) from a joint analysis of the four auto- and cross-bispectra, but also to independently reconstruct their shapes. In section 6 we discuss the implementation and validation of our method. Our conclusions are then summarized in the final section.

2 Optimal $f_{\rm NL}$ estimation and non-Gaussian map generation

We consider a pair of weakly NG fields, $Y$ and $Z$, and define their angular auto- and cross-bispectra as follows:

[TABLE]

where $f_{\rm NL}$ is the usual NG-amplitude parameter, defining the strength of the NG signal. For example, in CMB studies of primordial NG, $Y$ and $Z$ can be any two between the temperature, E-mode polarization and B-mode polarization fields. Here we note that while in this paper we demonstrate the method for fields on a 2 sphere, with the specific example of the CMB, the method can be trivially extended to any general set of $M$ correlated fields in $N$-dimensions.111For example, a complete CMB bispectrum estimation requires the method for $M=3$ (in order to deal with $TEB$ correlator). It can also be achieved by a simple extension of our $M=2$ method.

We start here by focusing on the problem of optimal $f_{\rm NL}$ estimation, using all the available auto- and cross-bispectra between $Y$ and $Z$. In other words, we consider the problem of estimating the bispectrum amplitude by selecting a specific NG model and fitting the theoretically predicted shape to the data. The optimal joint-estimator for $f_{\rm NL}$ was originally described in refs. [19, 20] and it can be written as:222 This formula is perturbatively derived assuming the weakness of NG and hence is not suitable for analyzing strong NG templates such as some secondary bispectra and bispectra estimated using fields in the nonlinear regime. Moreover, this is justified only when the covariance matrix is approximately diagonalized. Such a situation has actually been realized in the Planck 2015 data analysis owing to a recursive inpainting process [3]. We explore this case in the following and derive the fast modal estimator (3.19). If discarding the diagonal covariance matrix approximation, our derivation will be slightly complicated by replacing $\sum_{X^{\prime}}(C_{\ell}^{-1})^{XX^{\prime}}a_{\ell m}^{X^{\prime}}$ with $\sum_{X^{\prime}}\sum_{\ell^{\prime}m^{\prime}}({\cal C}^{-1})_{\ell m\ell^{\prime}m^{\prime}}^{XX^{\prime}}a_{\ell^{\prime}m^{\prime}}^{X^{\prime}}$, where ${\cal C}_{\ell m\ell^{\prime}m^{\prime}}^{XX^{\prime}}\equiv\Braket{a_{\ell m}^{X*}a_{\ell^{\prime}m^{\prime}}^{X^{\prime}}}$. It is then expected that the inverse of the covariance matrix of the $\beta$ coefficients are inserted into the estimator (3.19) as Ref. [18] reported.

[TABLE]

where $X_{i}$ and $X_{i}^{\prime}$ span over $Y$ and $Z$, and $(C_{\ell}^{-1})^{XX^{\prime}}$ is the $(X,X^{\prime})$ element of the inverse 2D angular power spectrum matrix $\left(\begin{smallmatrix}C_{\ell}^{YY}&C_{\ell}^{YZ}\\ C_{\ell}^{YZ}&C_{\ell}^{ZZ}\end{smallmatrix}\right)^{-1}$, ${\cal B}_{\ell_{1}\ell_{2}\ell_{3}}^{X_{1}^{\prime}X_{2}^{\prime}X_{3}^{\prime}}$ are the observed angle-averaged bispectra given by

[TABLE]

and $F$ is the Fisher matrix:

[TABLE]

Note that we reported here explicitly the $(-1)^{\ell_{1}+\ell_{2}+\ell_{3}}$ term, which is usually dropped in the literature, since we are interesting here in a general estimator dealing with both even and odd $\ell_{1}+\ell_{2}+\ell_{3}$ configurations.

Our goal in the next few sections will be that of defining a general modal approach to implement this optimal joint-estimator of $f_{\rm NL}$. We will however point out in section 5 that, with some simple, direct rearrangement of the output and no further computational cost, the same modal estimation pipeline also produces model-independent reconstruction of the full auto- and cross-bispectra from the data, as well as fast and optimal $f_{\rm NL}$ measurements, obtained by fitting $YYZ$, $ZZY$, $YYY$ and $ZZZ$ independently, and not jointly with all other combinations.

For estimator testing purposes, it is useful to have a procedure to generate NG maps with given power spectrum and bispectrum. A fast algorithm achieving this goal was originally proposed in ref. [21] for auto-spectra, and generalized to multiple fields with non-vanishing cross-spectra in ref. [18].

A set of multipoles $(a_{\ell m}^{Y},a_{\ell m}^{Z})$, with specified two-point correlation $\braket{a_{\ell_{1}m_{1}}^{X_{1}}a_{\ell_{2}m_{2}}^{X_{2}}}=C_{\ell_{1}}^{X_{1}X_{2}}(-1)^{m_{1}}\delta_{\ell_{1},\ell_{2}}\delta_{m_{1},-m_{2}}$ and the 3-point one (2.1) can be generated according to the following prescription:

[TABLE]

where $g_{\ell m}^{X}$ are orthonormal Gaussian realizations obeying $\braket{g_{\ell_{1}m_{1}}^{X_{1}}g_{\ell_{2}m_{2}}^{X_{2}}}=(-1)^{m_{1}}\delta_{\ell_{1},\ell_{2}}\delta_{m_{1},-m_{2}}\delta_{X_{1},X_{2}}$ and

[TABLE]

We will show later how the generation of NG maps can be made computationally efficient using the modal approach.

3 Modal estimator for parity-even bispectra

The well-known, general problem with formula (2.2) is that the brute-force computation of this estimator demands ${\cal O}(\ell_{\rm max}^{5})$ operations, where $\ell_{\rm max}$ is the maximum number of available multipoles. When $\ell_{\rm max}\gtrsim 10^{3}$ (comparable to Planck resolution, in the CMB case [2, 3]), this becomes a completely unfeasible task. In order to reduce the computational cost in bispectrum estimation, a variety of computational methodologies have been developed. We here work with one of them, the so-called Modal methodology, originally developed in refs. [14, 15]. For other approaches, the reader can refer to, e.g., ref. [22] (KSW estimator), ref. [23] (Binned estimator), and ref. [24] (Wavelet method).333 In fact both the binned and wavelet approaches, and indeed all optimal methods for bispectrum estimation, can be encompassed by the modal approach simply by changing the set of basis functions to the appropriate form, e.g. top hat functions in $\ell$-space for bins or the harmonic transform of the wavelet for wavelet approaches. This is a reflection of the fact that the only solution to evaluating the optimal bispectrum estimator is to represent the bispectrum in a separable form and the modal method is simply a generic procedure for doing this.

It was first realized by the authors of ref. [22] that the number of operations required for NG estimation can be reduced to ${\cal O}(\ell_{\rm max}^{3})$, with a massive ${\cal O}(\ell_{\rm max}^{2})$ gain, if the theoretical bispectrum template in input is written in separable form. This is however often not the case. The point of the modal methodology is that of tightly approximating the input bispectrum, by expanding it into finite sets of separable eigenfunctions in bispectrum space. This approach was initially developed for parity-even, temperature-only, CMB bispectrum measurements [14, 15] and was later extended to the parity-odd temperature-only case [16, 17]. In the following, we consider the direct extension of the modal methodology for auto-bispectra, with the goal to perform joint estimation of CMB temperature and polarization bispectra, or any other general combination of auto- and cross-bispectra.

In this section, we consider parity-even angle-averaged bispectra, the nonvanishing signal of which is confined to

[TABLE]

Such bispectra can be expressed as

[TABLE]

where a geometrical factor:

[TABLE]

filters the multipoles satisfying eq. (3.1). The reduced bispectra $b_{(e)\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$ contain all physical information.

3.1 Modal decomposition

Our initial goal is that of finding a factorizable representation for a general mixed reduced bispectrum, $b_{(e)\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$. This is done here via a straightforward generalization of the modal approach for auto-bispectra, described in ref. [14]. To start with, we consider a set of factorizable basis templates, $Q_{ijk}^{X_{1}X_{2}X_{3}}$, defined as:

[TABLE]

where $q_{i}^{X}(\ell)$ is a general set of eigenfunctions, which can be in principle arbitrarily chosen. Different choices can however have a strong impact on the rate of convergence for specific input theoretical shapes. Common choices are a polynomial set or a set of oscillatory trigonometric functions; see refs. [14, 18] for more details on this issue. The notation $[i,j,k]$ represents permutations over the indices $i$, $j$ and $k$, in order to obtain suitably symmetrized representations, for either auto- or cross-bispectra, i.e.

[TABLE]

Now we want to approximate the reduced bispectrum in a way that ensures the minimum error is introduced to the estimator (2.2). To do this we must decompose the signal-to-noise weighted reduced bispectrum which we can easily determine from the Fisher matrix (2.4), to be:

[TABLE]

Unfortunately $h_{\ell_{1}\ell_{2}\ell_{3}}$ is not easily separable and so we instead mimic it with the good separable approximation $v_{\ell_{1}}v_{\ell_{2}}v_{\ell_{3}}$ where $v_{\ell}=(2\ell+1)^{1/6}$; see again [14] for a more detailed discussion on this specific point. We can now decompose the approximately signal-to-noise weighted bispectrum into our separable modal basis $Q_{ijk}^{X_{1}X_{2}X_{3}}(\ell_{1},\ell_{2},\ell_{3})$, as

[TABLE]

The level of convergence of the expansion to the original shape is measured in terms of the following shape correlator, which naturally defines a (parity-even) inner product in bispectrum space:

[TABLE]

where

[TABLE]

and the squared term in brackets corrects for using our approximation to $h_{\ell_{1}\ell_{2}\ell_{3}}$ in the decomposition. Note also that in this last expression we have labeled the triples $ijk$ in $Q_{ijk}^{X_{1}X_{2}X_{3}}$ by means of a single index $n$. Using this, the modal coefficients are expressed as

[TABLE]

Note that, in general, $\gamma_{np}^{X_{1}X_{2}X_{3}}\neq\delta_{np}$, and the $Q_{n}^{X_{1}X_{2}X_{3}}$ templates then form a non-orthonormal basis. We can however always produce an orthonormal basis $R_{n}$, starting from $Q_{n}$, via a suitable rotation:

[TABLE]

In the formula above, the lower triangular, rotation matrix $\lambda_{np}^{X_{1}X_{2}X_{3}}$ is related to $\gamma_{np}^{X_{1}X_{2}X_{3}}$ via a Cholesky decomposition:

[TABLE]

In total analogy with what done above, one can expand the input bispectrum replacing the $Q_{n}^{X_{1}X_{2}X_{3}}$ basis templates with their orthonormal counterparts $R_{n}^{X_{1}X_{2}X_{3}}$, extracting the new expansion coefficients $\alpha_{n}^{RX_{1}X_{2}X_{3}}$ via the inner product $\braket{b,R}_{e}$. The relation between $\alpha_{n}^{X_{1}X_{2}X_{3}}$ and $\alpha_{n}^{RX_{1}X_{2}X_{3}}$ is then obtained as

[TABLE]

It must however be noted that the orthonormal $R$ templates are generally not separable. Since separability is crucial for fast estimation, the actual estimator employs the non-orthonormal $Q$-basis, and the rotation is then often applied to the measured template amplitudes, in order to simplify and make the interpretation of the results more transparent.

In the modal decomposition process, supposing that the input bispectrum, $b_{(e)\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$, is precomputed, the most time-consuming task comes from the computation of the inner products $\gamma\sim\Braket{Q,Q}_{e}$ and $\alpha\sim\Braket{b,Q}_{e}$. This computation requires ${\cal O}(\ell_{\rm max}^{3})$ operations in a brute-force approach and it is inevitable when obtaining $\alpha_{n}^{X_{1}X_{2}X_{3}}$, unless the starting bispectrum is already written in separable form; note, however, that this calculation is required only once, for a given shape. The expansion coefficients $\alpha$ can be computed and stored away for each shape under study, without constituting a bottleneck for the data estimation pipeline. Since the $Q$-basis is separable, we can always speed up the calculation of $\gamma_{np}^{X_{1}X_{2}X_{3}}$, by using the separable formula [21]:

[TABLE]

The derivation of this formula is based on the following integral representation of the Wigner $3j$ symbols:

[TABLE]

where $P_{\ell}(\mu)$ is the Legendre polynomial. Employing eq. (3.14) allows for a massive speed up of the algorithm, since it reduces the total number of operations by a factor $\ell_{\rm max}$.

3.2 Separable estimator

As firstly pointed out in ref. [22], if we can factorize the bispectrum, it is always possible to write the estimator (2.2) in separable form, by resorting on the identity (Gaunt integral):

[TABLE]

If we consider our factorized decomposition for auto- and cross-bispectra, we can generalize the modal estimator derivation of ref. [14], obtaining that the joint auto- and cross-bispectrum estimator can be rewritten into an integral of the products of the following filtered maps [note how separability is made possible also by the choice of weights in eq. (3.5)]:

[TABLE]

Using this expression, eq. (2.2) reduces to

[TABLE]

where the new coefficients $\beta_{n}^{X_{1}X_{2}X_{3}}\equiv\beta_{{\rm cub}\,n}^{X_{1}X_{2}X_{3}}-\beta_{{\rm lin}\,n}^{X_{1}X_{2}X_{3}}$ contain all the information about the observed bispectra ${\cal B}_{\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$. Their cubic and linear parts are represented, respectively, as

[TABLE]

To summarize, the simplicity of the estimator in eq. (3.19) has been created by the linear transformation of the observed $a^{X}_{\ell m}$ described by eq. (3.18). It is this transformation which removed all the cross terms from the estimator. This is what will allow us to simply estimate bispectra both jointly and individually using this method. This is in contrast to other methods like that in ref. [18] that instead removes the cross terms in the estimator by orthonormalizing the $a^{X}_{\ell m}$ which entangles the individual bispectra so they can only be jointly constrained. The two approaches can be simply related by:

[TABLE]

where $\hat{a}_{\ell m}^{X}\equiv\sum_{X^{\prime}}L_{\ell}^{XX^{\prime}}a_{\ell m}^{X^{\prime}}$, with $L_{\ell}^{XX^{\prime}}$ the rotational matrix (2.10), are the orthonormalized multipoles of ref. [18].

In order to get a complete expression for a fast estimator, we need to ensure that also the normalization can be computed fast. This requires fast numerical evaluation of the Fisher matrix for our expanded templates. If we deal only with auto-bispectra ($YYY$), then this is straightforward, since it is easy to verify, in such case, that:

[TABLE]

where $n_{\rm max}$ is the total number of basis templates used in the expansion. The situation becomes however more complex in the joint auto- and cross-bispectrum estimation case which we are considering here, since now the different bispectra are correlated, and we need to evaluate their full covariance. The straightforward, but slow approach would be to explicitly compute our approximating bispectrum template by summing over the modes ($b\sim\sum_{n}\alpha_{n}Q_{n}$). We can then insert the re-summed $b_{(e)\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$ into eq. (2.4) and numerically evaluate this expression. This would scale like ${\cal O}(\ell_{\rm max}^{3})$, even when the input bispectrum shape is explicitly separable. Furthermore, the whole process needs to be repeated not just for any new shape, but also any time we change the noise properties in the input data, since the noise power spectrum appears in the weights. This is very inconvenient, especially when a large number of data validation and consistency checks are performed (as it was the case, for example, in the Planck primordial NG data validation campaign, for which the currently described pipeline was intensively used).

Fortunately, the separability of the modal basis allows also in this case for a large speed up of the algorithm. We start by expanding the Fisher matrix (2.4) in the $Q$-basis:

[TABLE]

where $\tau^{X_{1}X_{2}X_{3}X_{1}^{\prime}X_{2}^{\prime}X_{3}^{\prime}}_{np}$ is a $8(n_{\rm max}+1)\times 8(n_{\rm max}+1)$ real symmetric matrix, reading

[TABLE]

Exploiting eq. (3.15), we can then apply the same factorization trick as done in eq. (3.14)

[TABLE]

The number of required numerical operations falls again by ${\cal O}(\ell_{\rm max})$ times, in comparison with the brute force computation. Note how we are computing here the covariance of the basis modes, which we indicated as $\tau_{np}$. This is completely model independent. The computation of $\tau_{np}$ can therefore always be carried efficiently, even for non-separable input bispectra. Finally, since $\tau_{np}$ represents essentially the mode covariance matrix between different bispectra, it also plays a key role when the estimator is used for full, model-independent bispectrum reconstruction, rather than for measuring $f_{\rm NL}$; see section 5.3 for further explanations on this issue.

3.3 Non-Gaussian simulations

The procedure described in the previous section can be similarly applied to produce a factorized expression of the NG simulation formula, eq. (2.13), leading to:

[TABLE]

where

[TABLE]

This allows creating a NG map through ${\cal O}(\ell_{\rm max}^{4})$ operations.

4 Spin-weighted modal estimator for parity-odd bispectra

In several important and interesting cases, the angular bispectrum takes values on a parity-odd domain:

[TABLE]

This happens not only in presence of parity-breaking primordial models, but also e.g. when evaluating cross-bispectra between CMB temperature and polarization B-modes (TTB), in a parity-conserving Universe. It is therefore useful to extend the previously described algorithm for fast estimation, to include the parity-odd case.

The derivation procedure is essentially the same as in the parity-even case, aside from the use of a spin-weighted geometrical factor:

[TABLE]

Here, the bispectrum templates obeying the parity-odd condition (4.1), $B_{(o)\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$, are decomposed into

[TABLE]

where $b_{(o)\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$ are the parity-odd reduced bispectra and

[TABLE]

Hence $b_{(o)\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}=b_{(o)\ell_{2}\ell_{1}\ell_{3}}^{X_{2}X_{1}X_{3}}=b_{(o)\ell_{1}\ell_{3}\ell_{2}}^{X_{1}X_{3}X_{2}}=b_{(o)\ell_{3}\ell_{1}\ell_{2}}^{X_{3}X_{1}X_{2}}=b_{(o)\ell_{2}\ell_{3}\ell_{1}}^{X_{2}X_{3}X_{1}}=b_{(o)\ell_{3}\ell_{2}\ell_{1}}^{X_{3}X_{2}X_{1}}$ holds. In the zero-spin limit, $x=y=z=0$, ${}_{[xyz]}H_{\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}=h_{\ell_{1}\ell_{2}\ell_{3}}^{0~{}0~{}0}=h_{\ell_{1}\ell_{2}\ell_{3}}$ just selects the even $\ell_{1}+\ell_{2}+\ell_{3}$ signal. In other words, in order to pick up all the nonvanishing signal satisfying both eq. (4.1) and $\ell_{1},\ell_{2},\ell_{3}\geq 2$, we here have to choose nonzero $x$, $y$ and $z$. Then, the selection rules of $h_{\ell_{1}\ell_{2}\ell_{3}}^{x~{}y~{}z}$ force us to obey $x+y+z=0$ and $-2\leq x,y,z\leq 2$. There are some possible combinations, and one of them is $(x,y,z)=(1,1,-2)$ [16].

In the next section, we will examine this case, and extend the parity-odd auto-bispectrum method developed in refs. [16, 17].

4.1 Modal decomposition

In analogy to parity-even bispectra, parity-odd reduced bispectra can be decomposed as

[TABLE]

where we include the imaginary unit in the weighting, and $\alpha_{n}^{X_{1}X_{2}X_{3}}\in\mathbb{R}$, since $b_{(o)\ell_{1}\ell_{2}\ell_{3}}^{X_{1}X_{2}X_{3}}$ take pure imaginary values.

The inner product in odd $\ell_{1}+\ell_{2}+\ell_{3}$ space is defined with a spin-dependent weighting; namely,

[TABLE]

Unlike the parity-even case (3.9), we here have to limit the summation range to $\ell_{1}+\ell_{2}+\ell_{3}={\rm odd}$ by hand. By employing it, the modal coefficients are expressed as

[TABLE]

The $\gamma_{np}^{X_{1}X_{2}X_{3}}$ matrix is defined and factorized as:

[TABLE]

where the maps depending on spin and parity are given by

[TABLE]

and $\sum_{a+b+c=o}$ denotes a summation over parity-odd pairs, i.e., $\sum_{a+b+c=o}{}^{(a)}\zeta{}^{(b)}\zeta{}^{(c)}\zeta={}^{(o)}\zeta{}^{(e)}\zeta{}^{(e)}\zeta+{}^{(e)}\zeta{}^{(o)}\zeta{}^{(e)}\zeta+{}^{(e)}\zeta{}^{(e)}\zeta{}^{(o)}\zeta+{}^{(o)}\zeta{}^{(o)}\zeta{}^{(o)}\zeta$. To derive the above, we have used the following formula:

[TABLE]

where ${}_{s}\lambda_{\ell m}(\mu)\equiv{}_{s}Y_{\ell m}(\mu,\phi)e^{-im\phi}$.

Similarly to the parity-even case, we can move to the orthonormal basis via the rotation (3.11). The decomposition into the orthonormal basis is written as

[TABLE]

where the coefficients, given by

[TABLE]

are related to $\alpha_{n}^{X_{1}X_{2}X_{3}}$ through eq. (3.13).

4.2 Separable estimator

The parity-odd estimator can be factorized in the same way as the parity-even case, only by replacing the Gaunt integral (3.16) with its spin-weighted version, reading

[TABLE]

We can now generalize eqs. (3.17) and (3.26) to include information on spin and parity, obtaining new filtered maps:

[TABLE]

where $\overline{a}_{\ell m}^{X}$ is given by eq. (3.18). The parity-odd estimator can now be written exactly in the same form as eq. (3.19); namely

[TABLE]

however, the coefficients $\beta_{ijk}^{X_{1}X_{2}X_{3}}=\beta_{{\rm cub}\,ijk}^{X_{1}X_{2}X_{3}}-\beta_{{\rm lin}\,ijk}^{X_{1}X_{2}X_{3}}$ and the matrix $\tau_{np}^{X_{1}X_{2}X_{3}X_{1}^{\prime}X_{2}^{\prime}X_{3}^{\prime}}$ are now composed of the new filtered maps:

[TABLE]

and the mode expansion of the bispectrum covariance now reads:

[TABLE]

This concludes the construction of our parity-odd fast estimator.

4.3 Non-Gaussian simulations

The algorithm to produce NG simulations can also be generalized to the parity-odd case. The NG part of the simulated $a_{\ell m}^{X}$, given by eq. (2.13), is now rewritten into the separable form:

[TABLE]

where

[TABLE]

This enables fast numerical computation, as in the parity even case.

5 General bispectrum estimation and bispectrum reconstruction

So far we have described the implementation of a general modal pipeline for joint estimation of the $f_{\rm NL}$ bispectrum amplitude parameter, using all possible bispectra built from the combination of two weakly NG field, $Y$ and $Z$. The starting fast estimator expression, eq. (2.2), was firstly obtained by the authors of ref. [20], who started from a $\chi^{2}$ statistic for the bispectrum, and managed to express the full covariance between different auto- and mixed-bispectra in factorizable form. This approach yields the cubic part of the estimator. The linear (“mean field”) term in $a_{\ell m}$ can then be computed by minimizing the variance of the cubic part, in presence of realistic isotropy-breaking effects, such as sky masking or an inhomogenous spatial noise distribution.

Since the derivation of the separable estimator formula, eq. (2.2), was only briefly sketched in ref. [20], we provide here a detailed derivation, starting from a different, yet customary, approach, based on an Edgeworth expansion of the $a_{\ell m}$ likelihood in presence of mild NG. We show in particular that starting from the likelihood allows us to obtain directly the factorizable estimator expression, without further manipulations of the bispectrum covariance. We can then straightforwardly derive expression for the fast, factorized estimator of $f_{\rm NL}$ from a given cross-bispectrum – say $B_{\ell_{1}\ell_{2}\ell_{3}}^{YYZ}$ – alone (not jointly with other cross-bispectra and auto-bispectra). Finally, and importantly, this approach also allows deriving a formula for full reconstruction of auto- and cross-bispectra, in a model-independent way, and at no computational extra-cost; we will discuss this further in sections 5.2 and 5.3.

The Edgeworth expansion of the $f_{\rm NL}$ likelihood, generalized from the auto-bispectrum case [25], is:

[TABLE]

where ${\cal C}_{\ell m\ell^{\prime}m^{\prime}}^{XX^{\prime}}\equiv\Braket{a_{\ell m}^{X*}a_{\ell^{\prime}m^{\prime}}^{X^{\prime}}}$. Then, because of the weakness of NG,

[TABLE]

holds, where, one can find, by explicitly taking the derivatives:

[TABLE]

with $C_{\ell m\ell^{\prime}m^{\prime}}^{XX^{\prime}}\equiv\Braket{a_{\ell m}^{X}a_{\ell^{\prime}m^{\prime}}^{X^{\prime}}}=(-1)^{m}{\cal C}_{\ell-m\ell^{\prime}m^{\prime}}^{XX^{\prime}}$.

This expression is the starting point for the derivation of the estimators, described in the following.

5.1 Amplitude estimator

We start by considering a general case in which each of the possible bispectra out of the combination of $Y$ and $Z$ has an independent amplitude, i.e., we now deal with the set of parameters $\{f_{\rm NL}^{YYY},f_{\rm NL}^{ZZZ},f_{\rm NL}^{YYZ},f_{\rm NL}^{ZZY}\}$. The usual expression (2.2) for the $f_{\rm NL}$ estimator is of course immediately derived from the results in this section, imposing that the four amplitudes are the same. The slightly more general case discussed here can be directly useful, for example, to test the presence of foreground residuals in primordial NG analysis of CMB temperature and polarization data. In principle, these residuals can in fact generate spurious bispectra with different amplitudes [26]. We note that, due to correlations between temperature and polarization data, a test of this kind would probably be robust only after preliminary projecting temperature out of the polarization component, in order to make the measured amplitudes independent. Moreover, the algebraic steps below, using independent amplitudes, can be straightforwardly generalized to obtain a (model-independent) mode- and bispectrum-reconstruction estimator, see section 5.3.

An optimal, unbiased estimator ${\cal E}_{I}$ must satisfy the following conditions (the former defining unbiasedness, and the latter the saturation of the Cramér-Rao bound):

[TABLE]

where $I,J\in\{YYY,ZZZ,YYZ,ZZY\}$ and the $4\times 4$ Fisher matrix is defined as

[TABLE]

With eqs. (5.2) and (5.3), the score function is computed as

[TABLE]

where $\displaystyle\sum_{X_{1}X_{2}X_{3}}^{\{YYY\}}\equiv\sum_{X_{1}X_{2}X_{3}}^{YYY}$ and $\displaystyle\sum_{X_{1}X_{2}X_{3}}^{\{YYZ\}}\equiv\sum_{X_{1}X_{2}X_{3}}^{\begin{subarray}{c}YYZ\\ YZY\\ ZYY\end{subarray}}$. For derivation, we have used the diagonal covariance matrix approximation:

[TABLE]

Each element of the Fisher matrix is thus given by

[TABLE]

Imposing the conditions (5.4) then leads to the following expression for the general estimator:

[TABLE]

We can now apply the modal decomposition to the expression above, leading to the final form for the modal estimator

[TABLE]

where $\alpha_{n}^{X_{1}X_{2}X_{3}}$, $\beta_{n}^{X_{1}X_{2}X_{3}}$ and $\tau_{np}^{X_{1}X_{2}X_{3}X_{1}^{\prime}X_{2}^{\prime}X_{3}^{\prime}}$ for even and odd $\ell_{1}+\ell_{2}+\ell_{3}$ modes are computed as described in sections 3 and 4, respectively. Notice that this expression is not the same as eq. (3.19), due to different normalization factors. Of course it reduces to eq. (3.19) when $f_{\rm NL}^{YYY}=f_{\rm NL}^{ZZZ}=f_{\rm NL}^{YYZ}=f_{\rm NL}^{ZZY}$.

5.2 Separable estimator for $f_{\rm NL}^{YYZ}$

Here we consider the estimation of $f_{\rm NL}^{YYZ}$ alone. Obviously, eqs. (5.9) and (5.10) directly yield the $f_{\rm NL}^{YYZ}$ estimator for models in which the other three bispectra vanish, i.e. $f_{\rm NL}^{YYY}=f_{\rm NL}^{ZZZ}=f_{\rm NL}^{ZZY}=0$ (for such situations see e.g., refs. [11, 12, 8, 9, 10]). The optimal estimator, obeying

[TABLE]

takes then the expected form:

[TABLE]

Via the modal decomposition, the estimator reduces to:

[TABLE]

5.3 Bispectrum reconstruction

In a previous work [14] it was shown that the modal estimator does not just directly allows us to estimate the amplitude $f_{\rm NL}$, but also to reconstruct the full auto-bispectrum of the map. To this purpose, it suffices to consider e.g. the estimated $\beta^{YYY}$ coefficients from eqs. (3.20) and (4.18), and build the linear combination:

[TABLE]

This can be generalized to the mixed bispectrum case. Therefore, as it is always the case with the modal approach, even if we started with an amplitude $f_{\rm NL}$ estimator, what we obtain is actually also a general method for full bispectrum reconstruction (with a finite resolution, imposed by the truncation of the modal expansion series, as it would be for any binning scheme). The bispectrum reconstruction formula for the full auto- and cross-bispectrum case can be obtained by simply taking the expected value of the $\beta$ coefficients, defined in eqs. (3.20) and (4.18). This yields:

[TABLE]

where the $\tau$ matrix was given by eqs. (3.25) and (4.19). Note that our previous derivation of the $f_{\rm NL}$ estimator, in the case with four independent amplitudes, can also be used to obtain this reconstruction expression, via estimation of the expansion coefficients $\alpha$. To get the optimal $\hat{\alpha}$ estimator, simply consider eq. (5.3) and replace $\braket{a_{\ell_{1}m_{1}}^{X_{1}}a_{\ell_{2}m_{2}}^{X_{2}}a_{\ell_{3}m_{3}}^{X_{3}}}$ with the modal expansion:

[TABLE]

Then, maximize with respect to the mode amplitudes $\alpha^{X_{1}X_{2}X_{3}}_{n}$, following completely analogous algebraic steps to what done for $f_{\rm NL}$ in the previous sections. The minimum variance estimate $\hat{\alpha}_{n}^{X_{1}X_{2}X_{3}}$ can now be re-inserted in eq. (5.19) to obtain eq. (5.16), which is then the optimal modal bispectrum reconstruction formula. The approach is of course valid both in the parity-even and parity-odd case.

6 Pipeline implementation and validation

In this section, we test the modal estimators developed above via the measurements of CMB temperature and E-mode polarization bispectra, i.e., $B_{\ell_{1}\ell_{2}\ell_{3}}^{TTT}$, $B_{\ell_{1}\ell_{2}\ell_{3}}^{EEE}$, $B_{\ell_{1}\ell_{2}\ell_{3}}^{TTE}$ and $B_{\ell_{1}\ell_{2}\ell_{3}}^{EET}$, for both parity-even and parity-odd cases.

For theoretical templates of the bispectra, we adopt the parity-even one from the usual scalar-mode equilateral NG [27] and the parity-odd one from the tensor-mode equilateral NG [7, 17, 3]. For these two cases, we produce NG simulations, using eqs. (3.27) and (4.20) respectively.

We assume the amplitude of the bispectra to be $f_{\rm NL}^{\rm equil}=300$ and $f_{\rm NL}^{\rm tens}=2000$. We start by computing the averages of our estimators from $50$ NG simulations, considering both joint and independent estimation. These are described respectively by eq. (2.2) – with $\beta$ adapted to either the parity-even or parity-odd case – and by eqs. (5.10) and (5.14).

We then extract the standard deviations, $\sqrt{\Braket{{\cal E}^{2}}}$ and $\sqrt{\Braket{{\cal E}_{I}^{2}}}$, from a separate set of 500 Gaussian realizations. For simplicity, we assume an ideal full-sky noiseless survey and hence all types of secondary uncertainties due to, e.g., finite beam widths, sky-cuts by masking, instrumental noises and foreground contaminations need not be taken into account (we note however again that this pipeline was thoroughly tested and applied to the primordial NG analysis of Planck data, where all these complications are obviously present).

In both cases, the modal decomposition is truncated when we reach a $95\%$ level of correlation between the original and the expanded shape. To speed up the computation, the multipole ranges are restricted to $2\leq\ell_{T}\leq 500$ and $35\leq\ell_{E}\leq 500$ in the parity-even analysis, and $2\leq\ell_{T},\ell_{E}\leq 100$ in the parity-odd one.

The averages and the standard deviations of each estimator are summarized in tables 2 and 2, respectively. We can see how the input $f_{\rm NL}$ values and the Cramér-Rao bounds are correctly recovered, both in the joint and independent estimation case; thus, we conclude that our estimators work properly.444One may feel at a glance that, in the parity-odd tensor NG case, $\braket{{\cal E}_{TTT}}$ and $\braket{{\cal E}_{EEE}}$ are a bit far from the input $f_{\rm NL}$ values, however, note that the error bars are also relatively large as shown in table 2, therefore the results are completely consistent.

7 Conclusions

In this paper, we have proposed and described in detail a general modal methodology for reconstructing and estimating amplitudes of auto- and cross-angular bispectra, considering two correlated weakly-NG fields. As usual in the modal bispectrum approach, any theoretically-predicted shape is decomposed using suitable, separable modal basis templates, which are then fit to the data. The freedom in the choice of the modal basis, and the possibility to use it to get factorizable expansions of any input shape, make this approach very versatile.

In a previous study [18], the modal decomposition was performed after the four input bispectrum templates ($YYY$, $ZZZ$, $YYZ$ and $ZZY$) were orthogonalized via a proper rotational transformation. In contrast, in the present approach, original bispectrum templates are decomposed as they are without rotation, and their covariance is dealt with later, again via a separable modal expansion. Both methodologies, of course, yield the same results when jointly constraining the amplitude of the four bispectra, $f_{\rm NL}$, given in eq. (2.1). The orthogonalized method is slightly faster for joint estimation, whereas the method developed here directly enables us to independently reconstruct and estimate amplitudes of single cross-bispectra, $f_{\rm NL}^{YYZ}$ or $f_{\rm NL}^{ZZY}$ (5.14), since those are decomposed from the start without rotation. The methodology presented here generalizes our previous work also by including the odd $\ell_{1}+\ell_{2}+\ell_{3}$ domain. This significantly extends the range of shapes which can be investigated. The pipeline presented here is an ideal starting point to include B-mode information in future modal CMB bispectrum analysis, since some bispectra have in this case odd-parity and we are explicitly interested in independent estimation of the TTB cross-bispectrum. [4, 5, 6, 7, 8, 9, 10]. However, while $TTB$ estimation will be one of the main future applications for our pipeline, it is by no means the only one. Other interesting cross-bispectra arise for example in NG studies of spectral distortion anisotropies [11, 12]. More in general, this methodology can be a flexible starting point for estimation of angular auto- and cross-bispectra for general random fields, including e.g. galaxy clustering and gravitational lensing studies. We showed in section 6 some applications and tests of our pipeline, including NG maps generation and bispectrum estimation, with both parity-even and parity-odd examples. We also extensively validated and applied this same pipeline, as one of the main methods for primordial NG analysis of Planck temperature and polarization data [3]. Further applications will be discussed in future studies.

Acknowledgments

M. S. was supported by JSPS Grant-in-Aid for Research Activity Start-up Grant Number 17H07319 and JSPS Grant-in-Aid for Early-Career Scientists Grant Number 19K14718. Numerical computations by M. S. were in part carried out on Cray XC50 at Center for Computational Astrophysics, National Astronomical Observatory of Japan. Part of this work was undertaken on the STFC COSMOS@DiRAC HPC Facility at the University of Cambridge, funded by UK BIS NEI grants