Exact joint likelihood of pseudo-$C_\ell$ estimates from correlated Gaussian cosmological fields

TL;DR

This paper derives an exact joint likelihood for pseudo-$C_5$ power spectrum estimates from Gaussian cosmological fields, applicable to various fields and masks, enhancing the robustness of cosmological parameter inference.

Contribution

It introduces a method to compute the exact joint likelihood of pseudo-$C_5$ estimates for correlated Gaussian fields with arbitrary masks, without assumptions on mask geometry.

Findings

Likelihood accurately recovers the multivariate distribution of $EE$, $BB$, and $EB$ spectra.

Method is applicable to CMB, weak lensing, and galaxy clustering analyses.

Simulations confirm the method's effectiveness in realistic scenarios.

Abstract

We present the exact joint likelihood of pseudo-$C_\ell$ power spectrum estimates measured from an arbitrary number of Gaussian cosmological fields. Our method is applicable to both spin-0 fields and spin-2 fields, including a mixture of the two, and is relevant to Cosmic Microwave Background, weak lensing and galaxy clustering analyses. We show that Gaussian cosmological fields are mixed by a mask in such a way that retains their Gaussianity, without making any assumptions about the mask geometry. We then show that each auto- or cross-pseudo-$C_\ell$ estimator can be written as a quadratic form, and apply the known joint distribution of quadratic forms to obtain the exact joint likelihood of a set of pseudo-$C_\ell$ estimates in the presence of an arbitrary mask. Considering the polarisation of the Cosmic Microwave Background as an example, we show using simulations that our likelihood recovers the full, exact multivariate distribution of $EE$, $BB$ and $EB$ pseudo-$C_\ell$ power spectra. Our method provides a route to robust cosmological constraints from future Cosmic Microwave Background and large-scale structure surveys in an era of ever-increasing statistical precision.

Equations

$$\alpha \left(\Omega \right)=\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{a}_{\ell m}^{\left(\alpha \right)}{Y}_{\ell m}\left(\Omega \right),$$

$$\begin{array}{rl}\alpha \left(\Omega \right)& =\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }\left(a_{\ell m}^{E\left(\alpha \right)}+i{a}_{\ell m}^{B\left(\alpha \right)}\right){}_2{Y}_{\ell m}\left(\Omega \right),\\ {\alpha }^{\ast }\left(\Omega \right)& =\sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }\left(a_{\ell m}^{E\left(\alpha \right)}-i{a}_{\ell m}^{B\left(\alpha \right)}\right){}_{-2}{Y}_{\ell m}\left(\Omega \right),\end{array}$$

$$a_{\ell -m}^{\left(\alpha \right)}={\left(-1\right)}^m{\left(a_{\ell m}^{\left(\alpha \right)}\right)}^{\ast },$$

$${\text{C}}_{\ell }^{\alpha \beta }=C_{\ell }^{\alpha \beta }.$$

$$\tilde{\alpha }\left(\Omega \right)=W_{\alpha }\left(\Omega \right)\alpha \left(\Omega \right).$$

$${\tilde{a}}_{\ell m}^{\left(\alpha \right)}=\sum _{{\ell }^{\prime }=0}^{\infty }\sum _{m^{\prime }=-{\ell }^{\prime }}^{{\ell }^{\prime }}{}_0{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}{a}_{{\ell }^{\prime }{m}^{\prime }}^{\left(\alpha \right)};$$

$${\tilde{a}}_{\ell m}^{E\left(\alpha \right)}=\sum _{{\ell }^{\prime }=0}^{\infty }\sum _{m^{\prime }=-{\ell }^{\prime }}^{{\ell }^{\prime }}\left(W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{+}{a}_{{\ell }^{\prime }{m}^{\prime }}^{E\left(\alpha \right)}+W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{-}{a}_{{\ell }^{\prime }{m}^{\prime }}^{B\left(\alpha \right)}\right);$$

$${\tilde{a}}_{\ell m}^{B\left(\alpha \right)}=\sum _{{\ell }^{\prime }=0}^{\infty }\sum _{m^{\prime }=-{\ell }^{\prime }}^{{\ell }^{\prime }}\left(W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{+}{a}_{{\ell }^{\prime }{m}^{\prime }}^{B\left(\alpha \right)}-W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{-}{a}_{{\ell }^{\prime }{m}^{\prime }}^{E\left(\alpha \right)}\right).$$

$${}_s{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}=\int d{\Omega }_s{Y}_{{\ell }^{\prime }{m}^{\prime }}\left(\Omega \right)W_{\alpha }{\left(\Omega \right)}_s{Y}_{\ell m}^{\ast }\left(\Omega \right),$$

$$W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{+}=\frac{1}{2}\left({}_2{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}{+}_{-2}{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}\right);W_{\ell {\ell }^{\prime }m{m}^{\prime }}^{-}=\frac{i}{2}\left({}_2{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}{-}_{-2}{W}_{\ell {\ell }^{\prime }}^{m{m}^{\prime }}\right),$$

$${\tilde{a}}_{\ell m}^{\left(\alpha \right)}=\sum _{{\ell }^{\prime }{m}^{\prime }}\sum _{\beta }\frac{\partial {\tilde{a}}_{\ell m}^{\left(\alpha \right)}}{\partial a_{{\ell }^{\prime }{m}^{\prime }}^{\left(\beta \right)}}{a}_{{\ell }^{\prime }{m}^{\prime }}^{\left(\beta \right)},$$

$$\begin{split}\operatorname{Re}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right)=\sum_{\beta,\ell^{\prime}}\Bigg{[}&\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}0}^{\left(\beta\right)}}\right)\operatorname{Re}\left(a_{\ell^{\prime}0}^{\left(\beta\right)}\right)\\ &+\sum_{m^{\prime}>0}\Bigg{(}\left[\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}\right)+\left(-1\right)^{m^{\prime}}\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\beta\right)}}\right)\right]\operatorname{Re}\left(a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}\right)-\left[\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}\right)-\left(-1\right)^{m^{\prime}}\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\beta\right)}}\right)\right]\operatorname{Im}\left(a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}\right)\Bigg{)}\Bigg{]};\end{split}$$

$$\begin{split}\operatorname{Im}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right)=\sum_{\beta,\ell^{\prime}}\Bigg{[}&\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}0}^{\left(\beta\right)}}\right)\operatorname{Re}\left(a_{\ell^{\prime}0}^{\left(\beta\right)}\right)\\ &+\sum_{m^{\prime}>0}\Bigg{(}\left[\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}\right)+\left(-1\right)^{m^{\prime}}\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\beta\right)}}\right)\right]\operatorname{Re}\left(a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}\right)+\left[\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}\right)-\left(-1\right)^{m^{\prime}}\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\beta\right)}}\right)\right]\operatorname{Im}\left(a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}\right)\Bigg{)}\Bigg{]}.\end{split}$$

$$Y=\sum _i{a}_i{X}_i,X_i\sim \mathcal{N}\left({\mu }_i,{\sigma }_i^2\right)⟹Y\sim \mathcal{N}\left(\sum _i{a}_i{\mu }_i,\sum _{i{i}^{\prime }}{a}_i{a}_{i^{\prime }}\mathrm{Cov}\left(X_i,X_{i^{\prime }}\right)\right),$$

$$\begin{split}\mathrm{Var}\left(\operatorname{Re}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right)\right)=\sum_{\beta,\gamma}\sum_{\ell^{\prime}}C_{\ell^{\prime}}^{\beta\gamma}\Bigg{[}\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}0}^{\left(\beta\right)}}\right)\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}0}^{\left(\gamma\right)}}\right)+\frac{1}{2}\sum_{m^{\prime}>0}\Bigg{(}&\operatorname{Re}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\gamma\right)}}\right)+\operatorname{Re}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\beta\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\gamma\right)}}\right)\\ &+\left(-1\right)^{m^{\prime}}\left[\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\gamma\right)}}\right)+\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\beta\right)}}\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\gamma\right)}}\right)\right]\Bigg{)}\Bigg{]};\end{split}$$

$$\begin{split}\mathrm{Var}\left(\operatorname{Im}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right)\right)=\sum_{\beta,\gamma}\sum_{\ell^{\prime}}C_{\ell^{\prime}}^{\beta\gamma}\Bigg{[}\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}0}^{\left(\beta\right)}}\right)\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}0}^{\left(\gamma\right)}}\right)+\frac{1}{2}\sum_{m^{\prime}>0}\Bigg{(}&\operatorname{Re}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\gamma\right)}}\right)+\operatorname{Re}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\beta\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\gamma\right)}}\right)\\ &-\left(-1\right)^{m^{\prime}}\left[\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\gamma\right)}}\right)+\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}-m^{\prime}}^{\left(\beta\right)}}\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime}m^{\prime}}^{\left(\gamma\right)}}\right)\right]\Bigg{)}\Bigg{]},\end{split}$$

$$\tilde{\text{a}}\sim \mathcal{N}\left(\text{0},\mathbf{\Sigma }\right),$$

$$\mathrm{Cov}\left(\sum _i{\alpha }_i{X}_i,\sum _j{\beta }_j{Y}_j\right)=\sum _{ij}{\alpha }_i{\beta }_j\mathrm{Cov}\left(X_i,Y_j\right),$$

$$\begin{split}\mathrm{Cov}\left(\operatorname{Re}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right),\operatorname{Re}\left(\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}\right)\right)=\sum_{\gamma,\varepsilon}\sum_{\ell^{\prime\prime}}C_{\ell^{\prime\prime}}^{\gamma\varepsilon}\Bigg{[}&\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}0}^{\left(\gamma\right)}}\right)\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}0}^{\left(\varepsilon\right)}}\right)\\ &+\frac{1}{2}\sum_{m^{\prime\prime}>0}\Bigg{(}\operatorname{Re}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\gamma\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)+\operatorname{Re}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\gamma\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(-1\right)^{m^{\prime\prime}}\left[\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\gamma\right)}}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)+\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\gamma\right)}}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)\right]\Bigg{)}\Bigg{]};\end{split}$$

$$\begin{split}\mathrm{Cov}\left(\operatorname{Im}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right),\operatorname{Im}\left(\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}\right)\right)=\sum_{\gamma,\varepsilon}\sum_{\ell^{\prime\prime}}C_{\ell^{\prime\prime}}^{\gamma\varepsilon}\Bigg{[}&\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}0}^{\left(\gamma\right)}}\right)\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}0}^{\left(\varepsilon\right)}}\right)\\ &+\frac{1}{2}\sum_{m^{\prime\prime}>0}\Bigg{(}\operatorname{Re}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\gamma\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)+\operatorname{Re}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\gamma\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}-\left(-1\right)^{m^{\prime\prime}}\left[\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\gamma\right)}}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)+\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\gamma\right)}}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)\right]\Bigg{)}\Bigg{]};\end{split}$$

$$\begin{split}\mathrm{Cov}\left(\operatorname{Re}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right),\operatorname{Im}\left(\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}\right)\right)=\sum_{\gamma,\varepsilon}\sum_{\ell^{\prime\prime}}C_{\ell^{\prime\prime}}^{\gamma\varepsilon}\Bigg{[}&\operatorname{Re}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}0}^{\left(\gamma\right)}}\right)\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}0}^{\left(\varepsilon\right)}}\right)\\ &+\frac{1}{2}\sum_{m^{\prime\prime}>0}\Bigg{(}\operatorname{Im}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\gamma\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)+\operatorname{Im}\left(\left[\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\gamma\right)}}\right]^{*}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)\\ &~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}+\left(-1\right)^{m^{\prime\prime}}\left[\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\gamma\right)}}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)+\operatorname{Im}\left(\frac{\partial\widetilde{a}_{\ell m}^{\left(\alpha\right)}}{\partial a_{\ell^{\prime\prime}-m^{\prime\prime}}^{\left(\gamma\right)}}\frac{\partial\widetilde{a}_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}}{\partial a_{\ell^{\prime\prime}m^{\prime\prime}}^{\left(\varepsilon\right)}}\right)\right]\Bigg{)}\Bigg{]}.\end{split}$$

$${\hat{C}}_{\ell }^{\alpha \beta }=\frac{1}{2\ell +1}\left[\mathrm{Re}\left(a_{\ell 0}^{\left(\alpha \right)}\right)\mathrm{Re}\left(a_{\ell 0}^{\left(\beta \right)}\right)+2\sum _{m=1}^{\ell }\left[\mathrm{Re}\left(a_{\ell m}^{\left(\alpha \right)}\right)\mathrm{Re}\left(a_{\ell m}^{\left(\beta \right)}\right)+\mathrm{Im}\left(a_{\ell m}^{\left(\alpha \right)}\right)\mathrm{Im}\left(a_{\ell m}^{\left(\beta \right)}\right)\right]\right]=\frac{1}{2\ell +1}\sum _{m=-\ell }^{\ell }{a}_{\ell m}^{\left(\alpha \right)}{\left(a_{\ell m}^{\left(\beta \right)}\right)}^{\ast },$$

$${\tilde{C}}_{\ell }^{\alpha \beta }=\frac{1}{2\ell +1}\left[\mathrm{Re}\left({\tilde{a}}_{\ell 0}^{\left(\alpha \right)}\right)\mathrm{Re}\left({\tilde{a}}_{\ell 0}^{\left(\beta \right)}\right)+2\sum _{m=1}^{\ell }\left[\mathrm{Re}\left({\tilde{a}}_{\ell m}^{\left(\alpha \right)}\right)\mathrm{Re}\left({\tilde{a}}_{\ell m}^{\left(\beta \right)}\right)+\mathrm{Im}\left({\tilde{a}}_{\ell m}^{\left(\alpha \right)}\right)\mathrm{Im}\left({\tilde{a}}_{\ell m}^{\left(\beta \right)}\right)\right]\right].$$

$${\tilde{C}}_{\ell }^{\alpha \beta }={\tilde{\text{a}}}^T{\text{M}}_{\ell }^{\alpha \beta }\tilde{\text{a}},$$

$$\tilde{\text{a}}=\left(\begin{array}{c}{\tilde{\text{a}}}^{\left(\mathrm{field}1\right)}\\ {\tilde{\text{a}}}^{\left(\mathrm{field}2\right)}\\ ⋮\\ {\tilde{\text{a}}}^{\left(\mathrm{field}N\right)}\end{array}\right);{\tilde{\text{a}}}^{\left(\mathrm{field}\alpha \right)}=\left(\begin{array}{c}{\tilde{\text{a}}}_{\ell =0}^{\left(\alpha \right)}\\ {\tilde{\text{a}}}_{\ell =1}^{\left(\alpha \right)}\\ ⋮\\ {\tilde{\text{a}}}_{{\ell }_{\max }}^{\left(\alpha \right)}\end{array}\right);{\tilde{\text{a}}}_{\ell }^{\left(\alpha \right)}=\left(\begin{array}{c}\mathrm{Re}\left({\tilde{a}}_{\ell 0}^{\left(\alpha \right)}\right)\\ \mathrm{Re}\left({\tilde{a}}_{\ell 1}^{\left(\alpha \right)}\right)\\ \mathrm{Im}\left({\tilde{a}}_{\ell 1}^{\left(\alpha \right)}\right)\\ \mathrm{Re}\left({\tilde{a}}_{\ell 2}^{\left(\alpha \right)}\right)\\ \mathrm{Im}\left({\tilde{a}}_{\ell 2}^{\left(\alpha \right)}\right)\\ ⋮\\ \mathrm{Re}\left({\tilde{a}}_{\ell \ell }^{\left(\alpha \right)}\right)\\ \mathrm{Im}\left({\tilde{a}}_{\ell \ell }^{\left(\alpha \right)}\right)\end{array}\right).$$

$$\displaystyle\textbf{{M}}_{\ell}^{\alpha\beta}=\begin{array}[]{r l}\mathit{Field:}&\begin{array}[]{c c c c c c c c}\hphantom{~{}\,}1&2&\ldots&\hphantom{~{}\,}\alpha&\hphantom{\,}\ldots&\hphantom{~{}\,}\beta&\hphantom{\,}\ldots&N\end{array}\\ \begin{array}[]{r}1\\ 2\\ \vdots\\ \alpha\\ \vdots\\ \beta\\ \vdots\\ N\end{array}&\left(\begin{array}[]{c c c c c c c c}\textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}\\ \textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}\\ \vdots&\vdots&\ddots&\vdots&\ddots&\vdots&\ddots&\vdots\\ \textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{M}}_{\ell}&\ldots&\textbf{{0}}\\ \vdots&\vdots&\ddots&\vdots&\ddots&\vdots&\ddots&\vdots\\ \textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{M}}_{\ell}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}\\ \vdots&\vdots&\ddots&\vdots&\ddots&\vdots&\ddots&\vdots\\ \textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}\end{array}\right)\end{array};\quad\textbf{{M}}_{\ell}=\begin{array}[]{c l}\ell^{\prime}:&\begin{array}[]{c c c c c c c}\hphantom{~{}\,}0&1&\ldots&\hphantom{~{}}\ell&\hphantom{\,}\ldots&\ell_{\mathrm{max}}\end{array}\\ \begin{array}[]{c}0\\ 1\\ \vdots\\ \ell\\ \vdots\\ \ell_{\mathrm{max}}\end{array}&\left(\begin{array}[]{c c c c c c c}\textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}\\ \textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}\\ \vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\ \textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{M}}&\ldots&\textbf{{0}}\\ \vdots&\vdots&\ddots&\vdots&\ddots&\vdots\\ \textbf{{0}}&\textbf{{0}}&\ldots&\textbf{{0}}&\ldots&\textbf{{0}}\\ \end{array}\right)\end{array};$$

$$\text{M}=\left(\frac{1+{\delta }_{\alpha \beta }}{2\ell +1}\right)\mathrm{diag}\left(\frac{1}{2},1,1,1,\dots \right).$$

$$\varphi\left(\textbf{{t}}\right)=\Bigg{\lvert}\textbf{{I}}-2i\sum_{\ell}\sum_{\alpha\beta}t_{\ell}^{\alpha\beta}\textbf{{M}}_{\ell}^{\alpha\beta}\bm{\Sigma}\Bigg{\rvert}^{-1/2},$$

$$\phi \left(\text{t}\right)=\int d\tilde{\text{C}}\exp \left(i\text{t}\cdot \tilde{\text{C}}\right)f\left(\tilde{\text{C}}\right),$$

$$f\left(\widetilde{\textbf{{C}}}\right)=\frac{1}{\left(2\pi\right)^{\nu}}\int d\textbf{{t}}\exp\left(-i\widetilde{\textbf{{C}}}\cdot\textbf{{t}}\right)\Bigg{\lvert}\textbf{{I}}-2i\sum_{\ell}\sum_{\alpha\beta}t_{\ell}^{\alpha\beta}\textbf{{M}}_{\ell}^{\alpha\beta}\bm{\Sigma}\Bigg{\rvert}^{-1/2}.$$

$$\tilde{\text{a}}={\left(\mathrm{Re}{\tilde{a}}_{20},\mathrm{Re}{\tilde{a}}_{21},\mathrm{Im}{\tilde{a}}_{21},\mathrm{Re}{\tilde{a}}_{22},\mathrm{Im}{\tilde{a}}_{22},\mathrm{Re}{\tilde{a}}_{40},\mathrm{Re}{\tilde{a}}_{41},\mathrm{Im}{\tilde{a}}_{41},\mathrm{Re}{\tilde{a}}_{42},\mathrm{Im}{\tilde{a}}_{42},\mathrm{Re}{\tilde{a}}_{43},\mathrm{Im}{\tilde{a}}_{43},\mathrm{Re}{\tilde{a}}_{44},\mathrm{Im}{\tilde{a}}_{44}\right)}^T,$$

$$y_{\ell }^{\alpha \beta }={\text{x}}^T{\text{E}}_{\ell }^{\alpha \beta }\text{x}\text{with }{\text{E}}_{\ell }^{\alpha \beta }=\frac{1}{2}{\text{C}}^{-1}\frac{\partial \text{C}}{\partial C_{\ell }^{\alpha \beta }}{\text{C}}^{-1},$$

Full text

††pubyear: 2019††pagerange: Exact joint likelihood of pseudo-$C_{\ell}$ estimates from correlated Gaussian cosmological fields–B

1 Introduction

In the age of precision cosmology, all possible sources of bias in the process of cosmological inference must be examined and controlled. One such source of bias that has recently been receiving an increasing amount of attention is the choice of likelihood function, which is an essential ingredient of Bayesian inference. The traditional, convenient choice is a multivariate Gaussian likelihood. Depending on the cosmological observable in question, this assumption may be anything from (near-)exact to severely wrong. In the latter case, this error can potentially propagate through to biased cosmological constraints.

Here we consider the likelihood of observed cosmological power spectra. Even in contemporary analyses these are routinely assumed to be Gaussian distributed (e.g. Hikage et al., 2019; Liu & Madhavacheril, 2019; Planck Collaboration et al., 2019), and the choice of covariance matrix in a Gaussian likelihood is an extremely active area of research (e.g. Kodwani et al., 2019; Harnois-Déraps et al., 2019). However, it is well known that the true power spectrum likelihood is strongly non-Gaussian, especially on large scales, both for Cosmic Microwave Background (Sun et al., 2013) and weak lensing observations (Hartlap et al., 2009; Sellentin & Heavens, 2018; Sellentin et al., 2018).

For full sky observations with isotropic noise and in the absence of systematic effects, the exact likelihood of observed power spectra is known (Percival & Brown, 2006). However, the situation is more complex for realistic observations, and specifically in the presence of a mask. Wandelt et al. (2001) presented the exact distribution of power in a single multipole of a single spin-[math] field in the presence of an azimuthally symmetric mask. Here we extend this to obtain the exact likelihood of an arbitrary number of multipoles of an arbitrary number of correlated spin-0 and spin-2 fields, each observed with a mask of any geometry. This is the first mathematically exact approach to such a generalisation. Previously, many approximate forms have been suggested, some of which attempt to model other realistic effects beyond a cut sky (Percival & Brown, 2006; Hamimeche & Lewis, 2008; Mangilli et al., 2015; Kalus et al., 2016). Alternative approaches include Gaussianisation of the data vector (e.g. Wang et al., 2019) or likelihood-free inference (Alsing & Wandelt, 2019; Taylor et al., 2019).

In this paper we will consider the exact distribution of power spectrum estimates measured on a cut sky using the pseudo-$C_{\ell}$ estimator. This estimator has the advantage of speed as it relies only on spherical harmonic transforms, which have fast implementations in software packages such as HEALPix (Górski et al., 2005). It has been used in the Hyper Suprime-Cam first-year weak lensing analysis in Hikage et al. (2019) (under the flat-sky approximation using Fourier transforms in place of spherical harmonic transforms) and is also planned to be used for the weak lensing analysis of future Euclid data. With the statistical power offered by next-generation experiments such as Euclid, it is essential to understand the likelihood of all estimators in order to ensure reliable cosmological constraints are obtained. Alternatives to the pseudo-$C_{\ell}$ estimator include the Quadratic Maximum Likelihood (QML) estimator (Tegmark, 1997), which is designed to return minimum-variance estimates provided a suitable choice of fiducial model is made. In Section 2 we derive the distribution of observed spherical harmonic coefficients, the pseudo-$a_{\ell m}$s, show that their Gaussianity is preserved in the presence of a mask—regardless of its geometry— and present exact expressions for the covariance of any pair of pseudo-$a_{\ell m}$s. In Section 3 we show that each pseudo-$C_{\ell}$ estimator may be written as a quadratic form in the pseudo-$a_{\ell m}$s, and use the known joint distribution of such quadratic forms to derive the full, exact joint likelihood of pseudo-$C_{\ell}$ estimates. We also show that the same formalism can be applied to obtain the exact joint likelihood of QML estimates. We consider the polarisation of the Cosmic Microwave Background (CMB) as an example in Sections 4 and 5, where we demonstrate that our pseudo-$C_{\ell}$ likelihood reproduces the exact distribution of $EE$, $BB$ and $EB$ pseudo-$C_{\ell}$ power spectra in the presence of a general mask. We discuss the implications of our work and present our conclusions in Section 6.

2 Pseudo-$a_{\ell m}$ distribution

In this section we derive the distribution of spherical harmonic coefficients on a general cut sky, the pseudo-$a_{\ell m}$s.

2.1 Full-sky $a_{\ell m}$ distribution

We begin by considering the distribution of $a_{\ell m}$s on the full sky. We assume correlated, isotropic Gaussian spin-0 or spin-2 cosmological fields $\alpha\left(\Omega\right)$, where $\Omega$ represents sky coordinates. We will use lower case Greek characters to represent cosmological fields throughout. Each spin-0 field can be decomposed in terms of spherical harmonics as

[TABLE]

where $a_{\ell m}^{\left(\alpha\right)}$ are the spherical harmonic coefficients for field $\alpha$. Each (complex) spin-2 field can be decomposed in terms of spin-weighted spherical harmonics as

[TABLE]

where the superscript $*$ denotes complex conjugation, and $a_{\ell m}^{E\left(\alpha\right)}$ and $a_{\ell m}^{B\left(\alpha\right)}$ are the spherical harmonic coefficients of the $E$- and $B$-mode components of field $\alpha$, respectively.

For spin-0 and spin-2 spherical harmonics, the $a_{\ell m}$s have the property that

[TABLE]

which further implies that the $m=0$ $a_{\ell m}$s must be real. For a given multipole $\ell$, the real and imaginary parts of $a_{\ell m}^{\left(\alpha\right)}$ for all $m>0$ are independently and identically distributed as multivariate Gaussian with mean 0 and covariance matrix $\frac{1}{2}\textbf{{C}}_{\ell}$, where the elements of $\textbf{{C}}_{\ell}$ are given by

[TABLE]

Here, each of $\alpha$ and $\beta$ can represent either a spin-0 field or the $E$- or $B$-mode component of a spin-2 field, and $C_{\ell}^{\alpha\beta}$ is the underlying angular cross-power spectrum between fields $\alpha$ and $\beta$ (which may be 0). $\textbf{{C}}_{\ell}$ is symmetric such that $C_{\ell}^{\alpha\beta}=C_{\ell}^{\beta\alpha}$, and $\alpha=\beta$ gives the auto-power spectrum $C_{\ell}^{\alpha\alpha}$. In the case of $m=0$, the lack of an imaginary part means that all of the underlying power is in the real part, and hence the $a_{\ell 0}^{\left(\alpha\right)}$ are multivariate Gaussian distributed with mean 0 and covariance matrix $\textbf{{C}}_{\ell}$. The real part of each $a_{\ell m}^{\left(\alpha\right)}$ is independent from all imaginary parts, and vice versa. All $a_{\ell m}^{\left(\alpha\right)}$ are independent between different $\ell$ and $m$; the only cross-correlation is between $a_{\ell m}^{\left(\alpha\right)}$s having the same $\ell$ and the same $m$, but different fields $\alpha$. The $m<0$ $a_{\ell m}$s can be regarded as deterministic functions of their positive-$m$ counterparts following Equation (3), rather than separate random variables.

2.2 Effect of a cut sky

In real space, the effect of a mask is to multiply each field by a window function $W_{\alpha}\left(\Omega\right)$, which can in general be unique to each field:

[TABLE]

The multiplication in real space is equivalent to a convolution in spherical harmonic space, which has the effect of mixing the $a_{\ell m}$s between different $\ell$, $m$, and between $E$- and $B$-modes in the case of spin-2 fields, to give the pseudo-$a_{\ell m}$s (Lewis et al., 2001; Brown et al., 2005):

[TABLE]

The spin-weighted spherical harmonic space window functions are given by

[TABLE]

where we note that the optional field-dependence of the mask reflected in $W_{\alpha}\left(\Omega\right)$ also means that each ${}_{s}W_{\ell\ell^{\prime}}^{mm^{\prime}}$ is implicitly field-specific, but we drop the $\alpha$ on the left-hand side to limit the number of indices. Following Brown et al. (2005) we have defined

[TABLE]

which may also be specific to each field.

2.3 Pseudo-$a_{\ell m}$ distribution

Equations (6)–(8) describe the effects of a mask on spin-0 and spin-2 fields. More generally, any observed pseudo-$a_{\ell m}$ can be written as a linear combination of full-sky $a_{\ell m}$s as

[TABLE]

where $\alpha$ and $\beta$ may each be either a single spin-0 field or the $E$- or $B$-mode component of a spin-2 field. We can expand this into real and imaginary parts as $\widetilde{a}_{\ell m}^{\left(\alpha\right)}=\operatorname{Re}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right)+i\operatorname{Im}\left(\widetilde{a}_{\ell m}^{\left(\alpha\right)}\right)$, and we show in Appendix A that the respective parts are given by

[TABLE]

[TABLE]

Each derivative in these equations is a constant weight that is uniquely specified by the geometry of the mask, and is independent of each $a_{\ell^{\prime}m^{\prime}}^{\left(\beta\right)}$. Therefore, Equations (12) and (13) are simply linear combinations of Gaussian random variables. Any linear combination of Gaussian variables is itself Gaussian distributed, with mean and variance adding linearly such that

[TABLE]

where we have followed the convention that $\sim$ means “distributed as” and $\mathcal{N}\left(\mu,\sigma^{2}\right)$ represents a univariate Gaussian distribution with mean $\mu$ and variance $\sigma^{2}$. $\mathrm{Cov}\left(\cdot,\cdot\right)$ represents the covariance of two random variables, and the covariance of any variable with itself is its variance. From Equation (14) we find that the real and imaginary parts of $\widetilde{a}_{\ell m}^{\left(\alpha\right)}$ each follow a Gaussian distribution with zero mean, since the full-sky $a_{\ell m}$s each have zero mean. We show in Appendix A that the variance of the real and imaginary part of each pseudo-$a_{\ell m}$ is

[TABLE]

[TABLE]

where $\beta$ and $\gamma$ are all fields correlated with field $\alpha$, each of which may be either spin-0 or the $E$- or $B$-mode component of a spin-2 field.

Since both the real and imaginary parts of each pseudo-$a_{\ell m}$ are Gaussian distributed, the joint distribution for both the real and imaginary parts of all pseudo-$a_{\ell m}$s for all fields, contained in the vector $\widetilde{\textbf{{a}}}$, can be described by a multivariate Gaussian distribution,

[TABLE]

with covariance matrix $\bm{\Sigma}$ whose diagonal elements are given by Equations (15) and (16). The off-diagonal elements can be calculated using the rule—which is not specific to the Gaussian distribution—that the covariance of two linear combinations of random variables is given by

[TABLE]

where $\alpha_{i}$ and $\beta_{j}$ are constant weights. We show in Appendix A that this leads to the following expressions for the elements of $\bm{\Sigma}$:

[TABLE]

[TABLE]

[TABLE]

This is our first key result: that the spherical harmonic coefficients of correlated Gaussian fields measured on a cut sky, which follow a multivariate Gaussian distribution, have mean 0 and covariance matrix $\bm{\Sigma}$ whose elements are given by Equations (19)–(21). For a given cosmological model and survey mask, these expressions can then be applied to a particular set of observational probes by specifying the fields in question ($\alpha,\beta$), the cosmological signals ($C_{\ell}^{\gamma\epsilon}$), and all non-zero derivative terms, which once again are completely specified by the geometry of the survey mask. (We demonstrate this explicitly for the case of the CMB later in Section 4.1.)

3 Pseudo-$C_{\ell}$ distribution

In this section we derive the distribution of pseudo-$C_{\ell}$ estimates for auto- and cross-power spectra of an arbitrary number of correlated spin-0 and spin-2 fields, for a general cut sky.

3.1 The pseudo-$C_{\ell}$ estimator

On the full sky, the covariance between either real parts or imaginary parts of two $a_{\ell m}$s having the same $(\ell,m)$ is equal to the corresponding underlying power spectrum, modulo a factor of $\frac{1}{2}$ for $m>0$ (see Section 2.1). Therefore, an unbiased estimator of the power spectrum is given by an appropriately weighted sample covariance of observed $a_{\ell m}$s:

[TABLE]

where we implicitly take the real part of the result. The auto-spectrum estimator is given by the special case where $\alpha=\beta$. Under a cut sky, this becomes the pseudo-$C_{\ell}$ estimator (Wandelt et al., 2001; Hivon et al., 2002; Brown et al., 2005), written in terms of the pseudo-$a_{\ell m}$s:

[TABLE]

We can write this as a matrix equation involving the vector of all pseudo-$a_{\ell m}$s, $\widetilde{\textbf{{a}}}$:

[TABLE]

where $\textbf{{M}}_{\ell}^{\alpha\beta}$ is a real symmetric matrix chosen to pick out the correct elements of $\widetilde{\textbf{{a}}}$ to match the expression in Equation (23). This may not be obvious, so to demonstrate we will choose an order for $\widetilde{\textbf{{a}}}$ and present the corresponding form for $\textbf{{M}}_{\ell}^{\alpha\beta}$, but we point out that this order is arbitrary and could be chosen, for example, to optimise computation of the likelihood. For $N$ correlated fields, each of which is either a spin-0 field or the $E$- or $B$-mode component of a spin-2 field, and a maximum measured multipole of $\ell_{\mathrm{max}}$, we can decompose $\widetilde{\textbf{{a}}}$ into three hierarchical levels as

[TABLE]

We can then similarly decompose $\textbf{{M}}_{\ell}^{\alpha\beta}$ into blocks to pick out the correct elements:

[TABLE]

The Kronecker delta $\delta_{\alpha\beta}$ accounts for the auto-spectrum case where $\alpha=\beta$. M has $2\ell+1$ diagonal elements, so the size of M is dependent on $\ell$ but its form is not. Each block in $\textbf{{M}}_{\ell}$ has $2\ell^{\prime}+1$ elements in each dimension, so $\textbf{{M}}_{\ell}$ has a total of $\sum_{\ell^{\prime}=0}^{\ell_{\mathrm{max}}}\left(2\ell^{\prime}+1\right)=\left(\ell_{\mathrm{max}}+1\right)^{2}$ elements along each side. This is also true of every zero block in $\textbf{{M}}_{\ell}^{\alpha\beta}$, so $\textbf{{M}}_{\ell}^{\alpha\beta}$ has a total of $N\times\left(\ell_{\mathrm{max}}+1\right)^{2}$ elements along each side. We note that the choice of $\ell_{\mathrm{\min}}=0$ is arbitrary, and any other $\ell_{\mathrm{\min}}$ could be chosen as long as this choice is consistent between $\textbf{{M}}_{\ell}^{\alpha\beta}$ and $\widetilde{\textbf{{a}}}$.

3.2 The joint distribution of quadratic forms

Each cut-sky cross- or auto-spectrum estimator can be written in the form of Equation (24), each with the same multivariate Gaussian vector $\widetilde{\textbf{{a}}}$ and different symmetric matrix $\textbf{{M}}_{\ell}^{\alpha\beta}$ depending on the field(s) and multipole in question. The joint characteristic function (CF) of an arbitrary number of such quadratic forms, $\varphi\left(\textbf{{t}}\right)$, has a known analytic form (Good, 1963):

[TABLE]

where $\lvert\cdot\rvert$ denotes the matrix determinant. The joint CF is defined as the expectation value $\langle e^{i\textbf{{t}}\cdot\widetilde{\textbf{{C}}}}\rangle$, which can also be written as a Fourier integral, i.e.

[TABLE]

where each $\left\{\widetilde{C}_{\ell}^{\alpha\beta},t_{\ell}^{\alpha\beta}\right\}$ are a Fourier pair111This should not be confused with the Fourier pair of $\left\{\ell,\theta\right\}$ which relate the power spectrum $C_{\ell}$ to the correlation function $\xi\left(\theta\right)$. and $\widetilde{\textbf{{C}}}$ and t are the respective vectors of all $\widetilde{C}_{\ell}^{\alpha\beta}$ and $t_{\ell}^{\alpha\beta}$ for all fields and multipoles in question. $f\left(\widetilde{\textbf{{C}}}\right)$ is the joint probability density function of all observed pseudo-$C_{\ell}$s. Equation (63) may then be inverted to yield the exact joint distribution of the pseudo-$C_{\ell}$s:

[TABLE]

Here $\nu$ is the length of the data vector, which is equal to $N\times\left(\ell_{\mathrm{max}}+1\right)^{2}$ when considering $N$ fields and all multipoles from [math] to $\ell_{\mathrm{max}}$. Equation (64) is normalised by construction, such that it integrates to 1 over the range of possible $\widetilde{\textbf{{C}}}$.

Equation (64) implicitly depends on the cosmological model via the covariance matrix of pseudo-$a_{\ell m}$s, $\bm{\Sigma}$, whose elements are given in Equations (19)–(21). These elements depend on both the underlying power spectra and the mask for each field. When regarded as a function of the model parameters $\bm{\theta}$ for fixed observed pseudo-$C_{\ell}$ estimates, Equation (64) is hence equivalent to the likelihood of the parameters, $\mathcal{L}\left(\bm{\theta}\right)$. Therefore, Equation (64) is our second key result: the exact joint likelihood of an arbitrary number of auto- and cross-pseudo-$C_{\ell}$ estimates from an arbitrary number of correlated spin-0 and spin-2 fields. The posterior distribution of the parameters $p\left(\bm{\theta}\right)$ may then be obtained via Bayes’ theorem, $p\left(\bm{\theta}\right)\propto\pi\left(\bm{\theta}\right)\mathcal{L}\left(\bm{\theta}\right)$, where $\pi\left(\bm{\theta}\right)$ is the prior distribution of the model parameters $\bm{\theta}$, and the proportionality constant is equal to the Bayesian evidence such that the posterior distribution is normalised.

We note the resemblance of Equation (64) to Equation (16) of Wandelt et al. (2001) and to Equation (C6) of Hamimeche & Lewis (2008), as all three are instances of the general distribution of quadratic forms. However, our Equation (64) is a much more generally applicable result than Equation (16) of Wandelt et al. (2001), which is restricted to a single multipole of a single spin-0 power spectrum observed with an azimuthally symmetric mask. Equation (C6) of Hamimeche & Lewis (2008), meanwhile, describes the distribution of a single multipole of a full-sky weighted cross-spectrum estimator.

3.2.1 The likelihood of a subset of pseudo-$C_{\ell}$ estimates

In Equations (25)–(61) we have presented a formalism for obtaining the full multivariate likelihood for all correlated fields and multipoles from $\ell_{\mathrm{min}}$ to $\ell_{\mathrm{max}}$. However, this likelihood may be adapted depending on exactly which pseudo-$C_{\ell}$s are required, and only the elements of $\bm{\Sigma}$ that contribute to the relevant estimators must be calculated. For example, if we wanted to calculate the joint likelihood of $\widetilde{C}_{\ell=2}^{\alpha\alpha}$ and $\widetilde{C}_{\ell=4}^{\alpha\alpha}$, we would have a pseudo-$a_{\ell m}$ vector of

[TABLE]

where we have omitted the $\left(\alpha\right)$ subscripts. We would then choose $\textbf{{M}}_{2}^{\alpha\alpha}$ and $\textbf{{M}}_{4}^{\alpha\alpha}$ to pick out the appropriate elements of $\widetilde{\textbf{{a}}}$.

It is mathematically straightforward to extend this formalism to obtain the exact distribution of deconvolved power estimates, as this is a linear operation in the pseudo-$C_{\ell}$s and hence a quadratic form in the pseudo-$a_{\ell m}$s. However, this approach offers no additional constraining power and introduces additional possible sources of error, so we choose not to discuss this further in this paper. Similarly, the extension to obtain the distribution of bandpowers through linear binning of multipoles—which is often necessary if one wishes to obtain deconvolved power estimates—is straightforward, requiring only appropriate changes to the selection matrices $\textbf{{M}}_{\ell}^{\alpha\beta}$.

3.2.2 The likelihood of Quadratic Maximum Likelihood estimates

At low $\ell$ one may prefer to use a Quadratic Maximum Likelihood (QML) power spectrum estimator (Tegmark, 1997), due to its optimality. The price is an increased computational cost, and diminishing returns at higher multipoles, compared to the pseudo-$C_{\ell}$ estimator (Efstathiou, 2004). By design, the QML estimator is also a quadratic form:

[TABLE]

where x is the vector of all pixel values, which is multivariate Gaussian distributed with covariance C. The $y_{\ell}$ may be scaled and linearly transformed to provide unbiased estimates of the underlying power spectrum, but—as with pseudo-$C_{\ell}$ estimates—this step is unnecessary for cosmological inference and here we will only discuss the distribution of $y_{\ell}$ values.

From Equation (66) and the general distribution of quadratic forms, it is straightforward to write down the joint distribution of a set of QML estimates y,

[TABLE]

where $\nu$ remains the length of the data vector, which is now equal to the total number of pixels across all maps. As with Equation (64) for pseudo-$C_{\ell}$ estimates, Equation (67) gives the likelihood of model parameters $\mathcal{L}\left(\bm{\theta}\right)$ when considered as a function of the model for fixed y. A potentially useful property of Equation (67) is that the matrices needed to evaluate it—C and the set of $\textbf{{E}}_{\ell}^{\alpha\beta}$—are the same matrices needed to evaluate the estimator itself, reducing the amount of additional work needed to evaluate the likelihood.

Theoretically, one could also obtain the joint distribution of both pseudo-$C_{\ell}$ and QML estimates. This would require writing both as quadratic forms in the same underlying quantity, such as the full-sky $a_{\ell m}$s. This is mathematically possible because the transforms from full-sky $a_{\ell m}$s to both pixels and cut-sky pseudo-$a_{\ell m}$s are linear, meaning that a quadratic form in one can also be written as a quadratic form in the other. These transforms could then be encoded in the selection matrix in place of $\textbf{{M}}_{\ell}^{\alpha\beta}$ in Equation (64) or $\textbf{{E}}_{\ell}^{\alpha\beta}$ in Equation (67), with the comparatively simple full-sky $a_{\ell m}$ covariance matrix in place of $\bm{\Sigma}$ or C respectively. However, in practice this is unlikely to be an attractive or computationally tractable option. An alternative would be to form the approximate joint distribution using a Gaussian copula, which would avoid the need to write both sets of estimates in terms of a common basis. Instead it would require the correlation coefficient or covariance between QML and pseudo-$C_{\ell}$ estimates, expressions for which are given in Efstathiou (2004).

4 Application to Cosmic Microwave Background polarisation

In this section we will demonstrate that the pseudo-$C_{\ell}$ likelihood presented in Section 3 is correct by calculating the full joint distribution of three multipoles of the CMB polarisation power spectra, and comparing to the observed distribution obtained from simulations. We will show that our likelihood exactly reproduces correlations between spectra. It also naturally models correlations between any number of multipoles, but here we consider three-dimensional distributions only to limit the computational resources needed to calculate each full distribution.

4.1 CMB pseudo-$\widetilde{a}_{\ell m}$ covariance

In the CMB case, the relevant fields are the temperature field, the $E$-mode polarisation and the $B$-mode polarisation. We therefore identify

[TABLE]

The fields are mixed as (Lewis et al., 2001; Brown et al., 2005)

[TABLE]

giving the relevant derivatives

[TABLE]

Here we consider only the polarisation fields, as CMB polarisation is the focus of many current and future experiments searching for evidence of inflation (e.g. Hui et al., 2018; Ade et al., 2019; Abazajian et al., 2016; Hazumi et al., 2019). However, our formalism naturally extends to include the temperature field, as well as including cross-correlation between fields observed by different detectors. Inserting the derivatives in Equation (70) into the general pseudo-$a_{\ell m}$ covariance matrix elements given in Equations (19)–(21) gives the elements of the CMB polarisation pseudo-$a_{\ell m}$ covariance matrix. Here we assume that all underlying $C_{\ell}^{EB}$ vanish, as is the case in any model in which parity is conserved, but Equations (19)–(21) naturally allow for non-zero underlying $C_{\ell}^{EB}$.

4.2 Implementation

The calculation of the pseudo-$a_{\ell m}$ covariance matrix requires theory power spectra and harmonic space window functions. We use camb (Lewis et al., 2000; Howlett et al., 2012) through the CosmoSIS222https://bitbucket.org/joezuntz/cosmosis interface (Zuntz et al., 2015) to generate CMB power spectra using the $\Lambda$CDM$+r$ model with Planck 2018 best fit parameters from Planck Collaboration et al. (2018) and tensor-to-scalar ratio $r=0.01$. We manually implement the harmonic space window function transform in Equation (9) using spin-weighted spherical harmonics provided by the spherical_functions Python package.333https://github.com/moble/spherical_functions However, we found that these only work reliably up to $\ell\approx 35$, so for our tests here we impose $\ell_{\mathrm{max}}=30$. In future work, we intend to calculate the harmonic space window functions by expressing them in terms of Wigner $3j$ symbols (Hivon et al., 2002) and making use of their recursion relations (Lewis et al. 2001; see also the Appendices of Hamimeche & Lewis 2008, 2009).

When calculating the joint characteristic function in Equation (62), we find it is more numerically stable to use an alternative form that avoids the need to evaluate the determinant of a complex matrix:

[TABLE]

i.e., the product is over all eigenvalues of the real matrix $\sum_{\ell}\sum_{\alpha\beta}t_{\ell}^{\alpha\beta}\textbf{{M}}_{\ell}^{\alpha\beta}\bm{\Sigma}$. We show in Appendix B that this form is mathematically equivalent to the form presented in Equation (62). To calculate the full joint likelihood distribution from the joint CF we implement Equation (64) by writing it in terms of a Fast Fourier Transform (FFT):

[TABLE]

where we have used the index $k$ as a proxy to encapsulate all summation variables ($\ell$, $\alpha$, $\beta$) in Equation (71). We use the numpy444https://numpy.org $n$-dimensional FFT (Oliphant, 2007):

[TABLE]

where each $t_{k}$ is discretised as $t_{k}=t_{k0}+m_{k}\Delta t_{k}$.

4.3 Simulations

We generated 276 million realisations of the CMB polarisation field from the fixed theory power spectra using the synfast routine in HEALPix555https://healpix.sourceforge.io (Górski et al., 2005) with a resolution of nside $=128$. To match the assumptions made in the theoretical distribution, we also impose $\ell_{\text{max}}=30$ in the input power. We note that the accuracy of a finite $\ell_{\text{max}}$ requires either a band-limited signal (such as is the case for the CMB), a well-behaved mask or apodisation to strongly suppress long-range mode mixing. For each realisation, we applied the polarisation field mask used in the WMAP 9-year analysis, corresponding to a sky fraction of $73.2\,\%$. The mask is described further in Bennett et al. (2013). We then measured the $EE$, $BB$ and $EB$ pseudo-$C_{\ell}$ power spectra from the masked maps using anafast. We formed marginal histograms of the power at $\ell=2$, 5 and 10 from each spectrum and a three-dimensional histogram of the measured $\ell=2$ power from all three spectra, which we selected to demonstrate the ability of our likelihood to exactly describe correlations between spectra. It also naturally (and exactly) describes correlations between multipoles, though we do not show this here to limit the computational complexity and number of required simulations. We have instead chosen to focus on the fact that our likelihood exactly models the mixing of power between $E$- and $B$-modes, as this is a particular difficulty for approximate methods due to the large discrepancy in the underlying $E$- and $B$-mode power. We have separately tested and confirmed the ability of our likelihood to exactly describe correlations between multipoles.

4.4 Comparison to approximation

We also consider an approximation to our exact likelihood, to demonstrate the importance of accurately modelling the full multidimensional distribution. On the full sky, the joint distribution of $\widehat{C}_{\ell}^{EE}$, $\widehat{C}_{\ell}^{BB}$ and $\widehat{C}_{\ell}^{EB}$ for a single fixed multipole $\ell$ follows a Wishart distribution (e.g. Percival & Brown, 2006):

[TABLE]

where $W_{p}\left(d,\textbf{{S}}\right)$ denotes the rank-$p$ Wishart distribution with $d$ degrees of freedom and scale matrix S. We adapt this for the cut sky using “effective" parameters of $\nu_{\mathrm{eff}}$ and $\textbf{{W}}_{\ell}^{\mathrm{eff}}$. We have chosen these effective parameters to be those that best fit the marginal distributions of $\widetilde{C}_{\ell}^{EE}$, $\widetilde{C}_{\ell}^{BB}$ and $\widetilde{C}_{\ell}^{EB}$, as calculated with our likelihood. We do not consider this to be a viable or attractive alternative to the exact likelihood; however, we include it to demonstrate that even an approximation that fits the marginals well fails to capture the full multidimensional likelihood. In addition to this, the Wishart approximation can only model correlations between different spectra for the same $\ell$, and not correlations between multipoles. Our likelihood, in contrast, exactly models correlations not only between spectra but also between multipoles, a behaviour which we have verified separately to the tests presented in this paper. However, as discussed above, we have chosen to focus on the former behaviour in this example.

A common alternative is to use a multivariate Gaussian approximation to the likelihood. With a suitable choice of mean, the problem then reduces to finding the most appropriate choice of covariance matrix, which can be calibrated with simulations. However, at low multipoles the true likelihood is so skewed that a Gaussian is an extremely poor approximation, which can lead to biased results and, if not modified, can assign non-zero probability to unphysical results such as negative auto-power spectra. For this reason, the 2018 Planck analysis uses a pixel-based likelihood below $\ell=30$ as described in Planck Collaboration et al. (2019).

5 Results

5.1 Marginal distributions

In Fig. 1 we show the marginal distributions of $\widetilde{C}_{\ell}^{EE}$, $\widetilde{C}_{\ell}^{BB}$ and $\widetilde{C}_{\ell}^{EB}$ for each of $\ell=2$, 5 and 10, to compare the prediction of our exact likelihood to the distributions observed in our simulations. Each histogram uses 300 bins. There is an excellent fit between the predicted and observed distribution, with no visible noise in the histograms due to the large number of events in each marginal distribution. The predicted likelihood exactly reproduces both the shape and amplitude of the observed distributions, including the considerable skewness in the auto-spectra. This skewness is reduced for higher multipoles, which is consistent with the full-sky behaviour of the likelihood. We have scaled each auto-$\widetilde{C}_{\ell}$ in Fig. 1 by the relevant theory $C_{\ell}$ used to generate both the theoretical likelihood and the simulations. In the case of the cross-spectrum $\widetilde{C}_{\ell}^{EB}$, there is no input $C_{\ell}^{EB}$ so we instead normalise by $\sqrt{C_{\ell}^{EE}C_{\ell}^{BB}}$. This scaling allows us to observe that the $E$-mode power is reduced by the sky cut while the $B$-mode power is increased. This is a result of $E$–$B$ mixing: the $EE$ power spectrum is much larger in magnitude than the $BB$ power spectrum (by a factor $\sim 200$ at $\ell=2$), meaning that the $E$–$B$ mixing induced by the mask leads to a relative increase in $B$-mode power at the expense of $E$-mode power.

5.2 Correlation between spectra

As well as exactly reproducing marginal distributions, our likelihood naturally describes correlations both between multipoles of the same spectrum and between spectra. As described in Section 4, we have formed the three-dimensional joint likelihood of $\widetilde{C}_{\ell}^{EE}$, $\widetilde{C}_{\ell}^{BB}$ and $\widetilde{C}_{\ell}^{EB}$ for $\ell=2$. We formed the corresponding simulated distribution by binning events in three dimensions, using 100 bins in each dimension. We then integrated our exact likelihood over the volume of each histogram bin to allow for comparison between theory and simulations. In Figs. 2–4 we show two-dimensional slices of this three-dimensional likelihood. Each slice corresponds to fixing a single histogram bin in one dimension and shows all bins in the other two dimensions. Our likelihood appears to accurately match the observed distributions in all six slices to within pixel noise that arises from the finite number of realisations in the simulations. The right-hand panel for each slice shows the logarithmic fractional residual, defined as

[TABLE]

Bins with no sampled events have $r=0$ and appear as white in Figs. 2–4. These areas were not explicitly excluded, but their probability is very low (albeit non-zero). We do not otherwise see any clear evidence of structure in these residuals, indicating that these bins contain only noise. This is mostly at the level of $r\approx-4$ to $-2$, except for a small number of outlying bins whose probability density is so low that the expected number of events in each bin from our 276 million simulations is significantly less than 1, leading to fractional residual values up to $r\approx 2$ in those bins in which an event was observed.

5.2.1 Comparison to approximation

In this section we compare to a Wishart distribution with fitted parameters $\nu_{\mathrm{eff}}$ and $\textbf{{W}}_{\ell}^{\mathrm{eff}}$, as described in Section 4.4. We are not advocating the use of this approximation; on the contrary, we aim to demonstrate the merits of using our exact likelihood. For this reason, we are not concerned with whether we can obtain appropriate values of $\nu_{\mathrm{eff}}$ and $\textbf{{W}}_{\ell}^{\mathrm{eff}}$ in practice. We simultaneously fitted the three marginals of a $p=2$ Wishart distribution to the exact marginals, and obtained the following best-fitting values:

[TABLE]

We have no reason to expect that these values would also be the best-fitting values for a different $\ell$ or for a different input cosmology, and they would certainly be different for another mask. We show the resulting marginal distributions in Fig. 5, where we have omitted the simulated histogram for clarity. The fit is almost perfect for $\widetilde{C}_{\ell}^{EE}$, indicating that this marginal distribution closely follows a gamma distribution as in the full-sky case. There is slightly more deviation for $\widetilde{C}_{\ell}^{BB}$ and $\widetilde{C}_{\ell}^{EB}$.

We integrated the Wishart probability density over each histogram bin in three-dimensional space, as we did with our exact likelihood. In Fig. 6 we show one-dimensional slices through this three-dimensional likelihood. For each slice, we have fixed a single bin for two dimensions and show the distribution across the third dimension. We see that while the marginal distributions suggest a near-exact fit for the Wishart approximation, the one-dimensional slices reveal that the approximation incorrectly distributes probability in some parts of the three-dimensional space relative to other parts. Our exact likelihood, on the other hand, faithfully reproduces the observed distribution throughout. This is also seen in Fig. 7, which shows three one-dimensional slices in different diagonal directions across the three-dimensional space. In each of these slices, the Wishart approximation overestimates the probability density relative to the true distribution, implying that it must underestimate it in other parts of the distribution, given that the whole distribution is normalised to integrate to 1.

As well as failing to accurately reproduce the full distribution for a fixed $\ell$, the Wishart distribution cannot naturally be extended to include correlations between multipoles, whereas our exact likelihood automatically produces the full joint distribution between all multipoles of all power spectra. In some cases an approximation may work better over the joint distribution of many multipoles than for a single $\ell$ (Hamimeche & Lewis, 2008), but the inverse can also be true (Elsner & Wandelt, 2012). In any case, the choice of approximation for the purposes of this comparison is unimportant: no approximation can completely match the exact distribution.

6 Conclusions

We have presented the exact joint likelihood of an arbitrary number of pseudo-$C_{\ell}$ estimates from correlated Gaussian fields, valid for both auto- and cross-power spectra and for any mask geometry. The likelihood—given in Equation (64)—naturally models both intrinsic correlations between spin-0 and spin-2 fields, and correlations induced by a cut sky which result in the mixing between spherical harmonic coefficients. The pseudo-$a_{\ell m}$s follow a multivariate Gaussian distribution with covariance matrix elements given by Equations (19)–(21). In Equation (67) we have also presented the exact joint likelihood for QML power spectrum estimates. An accurate likelihood function is an essential companion to any estimator for unbiased cosmological inference, but until now a complete likelihood for either the pseudo-$C_{\ell}$ or QML estimator on an arbitrary sky has not been known.

We have shown how our likelihood can be applied to observations of the polarisation of the Cosmic Microwave Background. This is especially relevant given current and future experiments aimed at detecting primordial $B$-modes, which require exquisite control of all possible sources of systematic bias. One such source of bias is an inexact likelihood function, so knowledge of the exact likelihood could play an important role in extracting cosmological information from polarisation measurements in an unbiased manner. In particular, it exactly models the leakage of $E$-mode power into the much smaller $B$-mode signal. Our likelihood also extends naturally to include correlations between temperature anisotropy and polarisation, including cross-correlations between any number of detectors, where each observed temperature or polarisation field may have its own mask. It does not account for weak gravitational lensing of the CMB, which breaks the assumption of Gaussianity at higher multipoles.

The exact likelihood will also be extremely useful for weak lensing observations. It will perhaps be most valuable at relatively low multipoles, as this is the regime where the common assumption of a Gaussian likelihood for power spectrum estimates is least applicable due to the considerable skewness in the true likelihood (e.g. Sellentin et al., 2018). These low multipoles correspond to large physical scales, for which it is an excellent approximation to describe the spin-2 cosmic shear field as a Gaussian field. At higher multipoles there will be significant deviations from Gaussianity, so our likelihood cannot be considered exact on small scales. The likelihood naturally extends to describe the full distribution of auto- and cross-power spectra between an arbitrary number of redshift bins, each with its own mask. It is thus well suited for extracting robust cosmological constraints from tomographic galaxy clustering and weak lensing shear power spectrum measurements in multiple redshift bins, such as those that will result from future analyses of the Euclid (Laureijs et al., 2011), LSST (LSST Science Collaboration et al., 2009) and SKA (Square Kilometre Array Cosmology Science Working Group et al., 2018) surveys.

Use of the exact pseudo-$C_{\ell}$ likelihood is likely to be competitive in terms of speed when considering a small number of bandpower estimates. Its main strength compared to an exact pixel-based likelihood is that for a single power estimate the pseudo-$C_{\ell}$ likelihood requires the determinant evaluation of a matrix of size $\sim\ell$ compared to $\sim\ell^{2}$ for a pixel-based method. This means that it may be evaluated at much higher $\ell$ than is possible for a pixel-based method. However, once many power estimates are considered the scaling is less competitive. The pixel-based method would offer all additional multipoles $\ell^{\prime}<\ell$ without significant additional computational cost, while the exact pseudo-$C_{\ell}$ likelihood scales as $\sim k^{N}$, where $N$ is the number of (band)power estimates and $k\sim 200$. The cost is driven by the need to evaluate the characteristic function in Equation (62) for every value of the vector t. The range of t must cover a wide enough space for the integral in the likelihood expression to converge, while at a sufficiently high resolution so that its curvature is accurately represented. Each point in t-space carries its own determinant or eigenvalue calculation (which is not the case for a single power estimate, due to the simple scaling of the determinant of a single matrix). We therefore recommend use of the exact likelihood only in the case of a small number of bandpower estimates, and are exploring alternative (necessarily approximate) approaches for the joint distribution of many power estimates, including the use of a copula with the exact marginal distributions. Copula methods have previously been described in a cosmological context, but only with approximate marginal distributions (Benabed et al., 2009; Sato et al., 2010, 2011). Alternative approaches previously explored in the literature include approximate extensions to the Wishart distribution to model correlations between multipoles (Hamimeche & Lewis, 2008; Mangilli et al., 2015). Computational limitations in the likelihood calculation may also be mitigated to some extent by potential speed increases at other levels in the inference process such as neural net-assisted sampling (Manrique-Yus & Sellentin, 2019).

Despite the limitations of its direct use, knowledge of the exact pseudo-$C_{\ell}$ and QML likelihood is extremely useful as a starting point and testing benchmark for developing fast, accurate approximations. The current best practice for cosmological inference (e.g. Planck Collaboration et al., 2019) is to use an exact pixel-based likelihood at low multipoles and to switch to an approximate Gaussian power spectrum likelihood for higher multipoles, at the point at which a pixel-based likelihood becomes computationally unfeasible (at $\ell=29$ in the case of Planck, while for weak lensing analyses with many redshift bins an exact pixel-based method may not be feasible at all). Methods derived from this exact likelihood may fill an important niche between these two regimes, allowing the use of an exact or near-exact likelihood up to higher multipoles than is currently possible. This may be a powerful tool for interpreting future observations, given the increased statistical precision that they will offer. Our likelihood also has the advantage that it can naturally describe the cross-correlation power spectrum measured between two different maps, in contrast to exact pixel-based methods which are not readily adapted to extracting just the cross-correlation information. Considering only cross-spectra in this way makes the cosmological analysis insensitive to the details of the noise bias, which will be especially relevant for cosmic shear observations, for which shape noise is an important and uncertain factor.

More work is required to adapt our techniques so that they are useful in practice. In future work, we plan to investigate the feasibility of implementing approximations to the exact likelihood as part of a realistic cosmological inference pipeline, and to benchmark their performance against alternatives such as an exact pixel-based likelihood. It also remains to be seen up to what maximum multipole these methods can be applied—both computationally, and also regarding the effects of non-Gaussianity at higher multipoles. A cut sky mixes power between multipoles, so even relatively low $\ell$ modes may be contaminated by non-Gaussian power in the presence of a severe mask. We plan to assess the impact of this, and possibly even incorporate it into a more advanced weak lensing likelihood. We have also not considered noise or systematic errors at this stage.

Acknowledgements

We thank Andrew Jaffe, Benjamin Joachimi and Andy Taylor for useful discussions, and the referee for helpful comments that improved the manuscript. REU is supported by a studentship from the UK Science and Technology Facilities Council. LW is supported by a UK Space Agency grant. Some of the results in this paper have been derived using the HEALPix (Górski et al., 2005) package.

Appendix A General pseudo-$a_{\ell m}$ covariance derivation

In this appendix we derive the general covariance matrix elements of the pseudo-$a_{\ell m}$s. First we present the derivation of the real and imaginary parts of the general pseudo-$a_{\ell m}$s given in Equations (12) and (13). We begin with the general form of a pseudo-$a_{\ell m}$ as a weighted sum of full-sky $a_{\ell m}s$ as given in Equation (11). We expand the sum over $m^{\prime}$ into separate sums for $m^{\prime}<0$ and $m^{\prime}>0$, and a term for $m^{\prime}=0$:

[TABLE]

Writing the $m^{\prime}<0$ contributions in terms of their $m^{\prime}>0$ counterparts using Equation (3) gives

[TABLE]

We now expand this into real and imaginary parts as

[TABLE]

which, after some algebra, gives the real and imaginary parts in Equations (12) and (13).

We now present the derivation of the covariance between two of the real parts, as given in Equation (19). The derivations of the imaginary–imaginary and real–imaginary covariance given in Equations (20) and (21) follow analogously. We begin by inserting the expression for the real part of a pseudo-$a_{\ell m}$ from Equation (12):

[TABLE]

We expand this into a linear sum of full-sky $a_{\ell m}$ covariances, noting that the terms involving covariance between real and imaginary full-sky $a_{\ell m}$s immediately vanish, as do terms involving covariance between $m=0$ and $m>0$ $a_{\ell m}$s:

[TABLE]

The full-sky covariances are given by

[TABLE]

and for $m>0$,

[TABLE]

Inserting these into Equation (81), evaluating the delta functions and rearranging, we obtain

[TABLE]

Finally, we use the identities that for complex $A$ and $B$,

[TABLE]

to obtain the final covariance given in Equation (19). The variance of a single real part is a special case of Equation (19) having $\alpha=\beta$, $\ell=\ell^{\prime}$ and $m=m^{\prime}$.

Appendix B Equivalence of the two forms of the joint characteristic function

Here we demonstrate that the expression for the joint characteristic function used in our implementation of the likelihood, given in Equation (71) and denoted here by $\varphi_{\mathrm{alt}}$, is mathematically equivalent to the known analytic form from Good (1963), given in Equation (62) and denoted here by $\varphi_{\mathrm{Good}}$.

The form given in Good (1963) is written in terms of the determinant of the complex square matrix which we will denote by X:

[TABLE]

where I is the identity matrix. The determinant of any matrix is equal to the product of its eigenvalues, so we may write $\varphi_{\mathrm{Good}}$ as

[TABLE]

The eigenvalues of X, $\lambda\left(\textbf{{X}}\right)$, are defined by the eigenvalue equation,

[TABLE]

Now let us define another matrix Y such that $\textbf{{X}}=\textbf{{I}}-\textbf{{Y}}$. Inserting this into Equation (88), we obtain

[TABLE]

For any matrix A, the determinant of $-\textbf{{A}}$ is given by

[TABLE]

depending on whether the rank of A is even or odd. Using this fact in Equation (89) implies that

[TABLE]

Therefore, the eigenvalues of X and Y are related as

[TABLE]

This allows us to write $\varphi_{\mathrm{Good}}$ in terms of the eigenvalues of Y,

[TABLE]

Finally, we use the fact that the eigenvalues of a scalar multiple of a matrix are equal to the scalar multiplied by the original matrix,

[TABLE]

to extract the multiple of $2i$ from Y:

[TABLE]

which we insert into Equation (93) to obtain our alternative form,

[TABLE]