© E. Camphuis et al. 2022

1 Introduction

One of the most powerful probes of cosmology is the observation of the cosmic microwave background (CMB) anisotropies. The NASA WMAP and the ESA Planck satellite CMB measurements marked the entry into the era of precision cosmology, in which many Λ cold dark matter (ΛCDM) matter cosmological parameters are measured with uncertainties smaller than 1% (Hinshaw et al. 2013; Planck Collaboration V 2020). Ongoing and upcoming ground-based and satellite experiments such as the Atacama Cosmology Telescope (ACT; Aiola et al. 2020), the South Pole Telescope (SPT; Dutcher et al. 2021), the Simons Observatory (SO; Ade et al. 2019), CMB-Stage 4 (CMB-S4; Gallardo et al. 2022), and Litebird (Hazumi et al. 2012) will provide yet more information about the nature of our Universe.

Because primary CMB anisotropies in intensity and polarization are distributed as a Gaussian random field, most of the cosmological information is contained in the angular power spectrum of the CMB anisotropies. As the evolution of the primary anisotropies is linear, the multipoles of the angular power spectrum are uncorrelated when the full sky is observed. However, any realistic experiment requires masking parts of the sky, either to avoid regions that are highly contaminated by foregrounds (e.g., galactic emission or point sources) or because the scanning strategy is designed to observe specific regions of the sky. The estimation of the power spectrum on the masked sky, the so-called pseudo-power spectrum, is biased, and different multipoles become correlated (Hivon et al. 2002). An unbiased estimator of the spectra can then be obtained through the MASTER approach (Hivon et al. 2002), as implemented in the PolSpice¹ software, for instance (Szapudi et al. 2001; Chon et al. 2004). A robust inference of cosmological parameters requires accurate covariance matrices that describe the variance of the spectra along their diagonal, as well as the correlations between multipoles in the off-diagonal terms. Pseudo-C_ℓ covariance matrices are corrected for the effect of the mask using MASTER to obtain the covariance matrices for the unbiased C_ℓ estimator. Inaccuracies in the covariance matrix estimation can lead to the misestimation of cosmological parameters and of their uncertainties (Dodelson & Schneider 2013; Sellentin & Starck 2019).

Covariance matrices can be calculated through the use of simulations. The number of simulations determines the accuracy of the estimator. As the simulations are expensive to produce, the obtained noisy realization of the covariance has to be regularized (Balkenhol & Reichardt 2022)². Alternatively, it is possible to calculate pseudo-C_ℓ covariance matrices analytically. However, these depend on integrals whose exact numerical implementation is computationally expensive. Thus, approximations have been proposed in previous works to make these calculations efficient; see, for example, Efstathiou (2004), Nicola et al. (2021) and Friedrich et al. (2021).

We analyze the problem of computing accurate analytical covariance matrices. We take the specific case of the South Pole Telescope third-generation (SPT-3G) experiment, which observes the CMB anisotropies in temperature and polarization on a small sky patch that corresponds to about 4% of the sky. On such a small sky region, the calculated power spectra has strong correlations between multipoles. The existing approximations of the covariance matrix can be less accurate in these conditions. Considering this particular case is thus a particularly stringent test of the validity of analytical algorithms.

We implement for the first time the exact computationally expensive calculation of the covariance matrices, which we find to be numerically feasible at multipoles smaller than ℓ ≲ ℓ_max,ex ≡ 1000 through a new algorithm that gains one order of numerical complexity over the brute-force approach, resulting in a thousandfold speed improvement. Then, we test the existing approximations, and find that they are accurate at the 5% level in the case of the SPT-3G footprint when the power spectrum is averaged in wide bandpowers. We then propose a new approximation that improves the existing algorithms to attain an accuracy of 1% in the same case. Finally, we describe how the covariance matrix of the PolSpice C_ℓ estimator can be calculated from the pseudo-C_ℓ covariance matrix.

While in this work, we focused on the specific case of the SPT-3G CMB experiment, many considerations can be applied more broadly to any probe relying on the calculation of power spectra and covariance matrices. We also highlight that in this work we only consider the signal-signal part of the covariance matrix. Analytical approaches to modeling the noise contribution have already been developed and used in CMB experiments, such as Planck; see Appendix C.1 of Planck Collaboration XI (2016). They rely on assumptions on the noise properties, such as isotropy or whiteness. When these assumptions do not hold, the noise-noise and signal-noise components can be obtained directly from the data, as was done in the SPT-3G analysis (Lueker et al. 2010; Dutcher et al. 2021). The integration of the noise part in our exact formalism and new approximation is beyond the scope of this paper.

The paper is organized as follows. In Sect. 2 we introduce the pseudo-power spectrum estimator and its covariance. In Sect. 3 we perform the exact calculation of the covariance matrix. In Sect. 4 we describe the existing approximations for the calculation of the covariance matrix, and we test their accuracy against the exact computation. Section 5 presents our new approximation, which is more accurate. Section 6 describes how the covariance matrix of the PolSpice estimator can be calculated. We conclude in Sect. 7, and some detailed calculations are given in the appendices.

2 Covariance of the pseudo-power spectrum

2.1 Pseudo-power spectrum

Cosmic microwave background anisotropies in intensity and polarization can be described as maps of Stokes parameters $T\left(\widehat{n}\right)$, $Q\left(\widehat{n}\right)$, $U\left(\widehat{n}\right)$ for a direction $\widehat{n}$ of the sky. They are Gaussian random fields, fully characterized by their angular power spectra $\left(C_{\ell }^{\text{TT}},C_{\ell }^{\text{EE}},C_{\ell }^{\text{BB}},C_{\ell }^{\text{TE}}\right)$, which are the variances of the harmonic coefficients $a_{\ell m}^{\text{T}}$, $a_{\ell m}^{\text{E}}$, $a_{\ell m}^{\text{B}}$ obtained by spherical harmonic decomposition of the maps. Cosmological models allow the computation of the expectation of the different power spectra in an ideal full-sky case. However, data only ever cover a part of the sky. We describe the partial coverage with the weight map $W\left(\widehat{n}\right)$. The power spectrum of masked maps, labeled ${\tilde{C}}_{\ell }^{\text{XY}}$, is usually called throughout pseudo-power spectrum. Its expression for temperature is given in Eq. (A.6). It can be computed from the masked harmonic coefficients ${\tilde{a}}_{\ell m}^{\text{X}}$, which are directly related to the unmasked ones, $a_{{\ell }^{\prime }{m}^{\prime }}^{\text{X}}$, by the application of the mode-coupling kernels _sI_{ℓmℓ′m′} [W]. In the case of temperature, we write $${\tilde{a}}_{\ell m}^{\text{T}}={\displaystyle \sum _{\ell \prime m\prime }{a}_{\ell \prime m\prime }^{\text{T}}{}_0{I}_{\ell m{\ell }^{\prime }{m}^{\prime }}\left[W\right],}$$(1)

where we have defined the mode coupling kernels³ $${}_s{I}_{\ell m{\ell }^{\prime }{m}^{\prime }}\left[W\right]\equiv {\displaystyle \int \text{d}{\widehat{u}}_s}{Y}_{\ell m}\left(\widehat{u}\right)W{\left(\widehat{u}\right)}_s{Y}_{{\ell }^{\prime }{m}^{\prime }}^{*}\left(\widehat{u}\right).$$(2)

These coupling kernels are an important component in the following discussions⁴. In the full-sky case, the orthonormality properties of spin-weighted spherical harmonics ensure that _sI_{ℓmℓ′m′} [1] = δ_ℓℓ′δ_mm′. We recall in Appendix A some summation properties of products of coupling matrices that appear in the computation of pseudo-power spectra and their covariance. In particular, they are related to the symmetric coupling kernel acting on a power spectrum A, labeled Ξ[A], with $${\text{Ξ}}_{\ell {\ell }^{\prime }}^{s{s}^{\prime }}\left[A\right]\equiv {\displaystyle \sum _L\frac{2L+1}{4\pi }{A}_L}\left(\begin{array}{lll}\ell \hfill & {\ell }^{\prime }\hfill & L\hfill \\ s\hfill & -s\hfill & 0\hfill \end{array}\right)\left(\begin{array}{lll}\ell \hfill & {\ell }^{\prime }\hfill & L\hfill \\ s^{\prime }\hfill & -s^{\prime }\hfill & 0\hfill \end{array}\right).$$(3)

This operator, introduced in Efstathiou (2004), can also be seen as acting on a map A with power spectrum A_ℓ. In the following, we use the notation Ξ[A] ≡ Ξ[A]. We recall that the average of the pseudo-spectrum is related to the underlying power spectrum by the application of the asymmetric coupling kernel computed for the mask W, also known as the MASTER mode-coupling matrix M. In the case of temperature, we have $$\langle {\tilde{C}}_{\ell }^{TT}\rangle ={\displaystyle \sum _{{\ell }^{\prime }}{}_0{M}_{\ell {\ell }^{\prime }}\left[W\right]C_{{\ell }^{\prime }}^{\text{TT}}}$$(4) $${}_0{M}_{\ell {\ell }^{\prime }}\left[W\right]\equiv \left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[W\right].$$(5)

In this work, without loss of generality, we develop the computations for the intensity case (i.e., s = s′ = 0), the polarization case being similar. We also assume that a single mask is used for both temperature and polarization. When required, we highlight the differences between the temperature and polarization cases and give insight into the importance of the single-mask assumption.

2.2 Covariance

Estimating the covariance of the measured power-spectrum is crucial to assess the agreement between data and model predictions and to constrain cosmological parameters from CMB maps. As discussed in Hivon et al. (2002) and demonstrated in Appendix A, masking breaks the statistical isotropy and induces correlations between the modes of the pseudo-spectrum. The details of the derivation of the analytical expression of the pseudo-spectrum covariance can be found in Appendix B. We give here the expression in terms of the coupling matrices ₀I and the true underlying intensity power spectrum C_ℓ, for the temperature case, $$\begin{array}{lll}{\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}\hfill & \equiv \hfill & \mathrm{cov}\left({\tilde{C}}_{\ell },{\tilde{C}}_{{\ell }^{\prime }}\right),\hfill \\ \hfill & =\hfill & \frac{2}{\left(2\ell +1\right)\left(2{\ell }^{\prime }+1\right)}{\displaystyle \sum _{m{m}^{\prime }}{\displaystyle \sum _{{\ell }_1{m}_1}{\displaystyle \sum _{{\ell }_2{m}_2}{C}_{{\ell }_1}{C}_{{\ell }_2}}}}\hfill \\ \hfill & \hfill & {\times }_0{I}_{\ell m{\ell }_{1}m{}_1}\left[W\right]{}_0{I}_{{\ell }^{\prime }{m}^{\prime }{\ell }_1{m}_1}^{*}\left[W\right]{}_0{I}_{\ell m{\ell }_2{m}_2}^{*}\left[W\right]{}_0{I}_{{\ell }^{\prime }{m}^{\prime }{\ell }_{2}m{}_2}\left[W\right].\hfill \end{array}$$(6)

As shown in Eq. (1), the mode-coupling coefficient ₀I kernels relate the underlying harmonic coefficients to the harmonic coefficients measured on the sky through the mask. In the analytic expression of the covariance, they represent the coupling between modes due to partial sky coverage. An expression similar to Eq. (6) can be written for polarization, using spin-2 spherical harmonics, that is, s, s′ = ±2. These expressions mix the EE and BB power spectra.

The expression in Eq. (6) involves several convolutions, and its evaluation is computationally expensive. The full computation scales as $O\left({\ell }_{\max }^6\right)$, ℓ_maxbeing the largest multipole, making the exact computation of this covariance a daunting task given the currently available computation power. We have developed an algorithm that allows the computation of the covariance matrix at low multipoles with a gain of an order of magnitude in computational time. We discuss this result in Sect. 3.

This approach was previously unavailable, therefore existing work has relied on approximations of Eq. (6). In Sect. 4 we present different approximations that have been proposed in previous works, and we then validate them against our full computation. This validation was performed for a small survey footprint, where spectral modes are highly correlated. These correlations can challenge the assumptions made in the different approximations. Throughout this work, we use a test-case inspired by SPT-3G. The footprint of the first year of the survey presented in Dutcher et al. (2021) covers roughly 4% of the sky and is displayed in Fig. 1 along with the mask power spectrum W_ℓ. We apodized the mask with a Gaussian window function of 30 arcmin full width at half maximum, using an algorithm similar to the one used in Planck (Planck Collaboration VI 2020). We also show in Fig. 1 the power spectrum of one of the masks used in the Planck cosmological analysis, which covers a much larger patch of the sky, around 70% before apodization. The precision of the standard approximation of the covariance was validated in the latter case, but it needs to be assessed for a smaller survey area.

3 Exact computation

An exact calculation of the pseudo-power spectrum covariance matrix can be obtained by integrating Eq. (6). We propose an algorithm that performs the computation in $O\left({\ell }_{\max }^5\right)$, typically gaining a thousandfold speed-up compared to the direct implementation in $O\left({\ell }_{\max }^6\right)$. This is achieved with the fast harmonic transform tools implemented in the HEALPix library⁵. It enables the exact computation of the covariance matrix, albeit on a limited range of multipoles. In this work, we have computed the full covariance up to ℓ_max,ex ≡ 1000, and calculated a few rows of the matrix at ℓ > ℓ_max,ex. This allows the direct comparison of the various analytic covariance approximation formulae with the exact calculation.

In the following, we describe the algorithm we developed to perform this computation. We validate it with Monte-Carlo estimates of the covariance for the reference SPT-3G survey.

Fig. 1 Survey area and the mask power spectrum. Top panel: CMB temperature anisotropies on the SPT-3G patch in galactic coordinates. The dark green line delimits the survey footprint. The vertical and horizontal bold black lines are the zero-longitude and zero-latitude coordinates, respectively. The SPT-3G patch covers roughly 4% of the sky. Bottom panel: mask power spectra as defined in Eq. (A.4) for SPT-3G and for the 143 GHz map used in the Planck cosmological analysis, which covers around 70% of the sky. The spectra have been renormalized by their first value for comparison purposes. Masks corresponding to small sky fractions, such as that of SPT-3G, have a shallower power spectrum than large ones. We emphasize that the mask used here does not include point-source masking.

3.1 Algorithm description

We focus on the computation of a given row of the covariance matrix $\widetilde{\text{Σ}}$. This allows us either to compute a full covariance matrix at low multipoles or to test our approximations on a selection of rows. We first start by describing the computation of the covariance of the intensity spectrum. A diagrammatic implementation of our calculation is presented in Fig. 2. For the polarization spectra, the calculation follows a similar pattern, and the difference between the two cases are discussed in the next section.

We derive in Appendix B the expression of the covariance of the pseudo-spectrum that leads to Eq. (6). In particular, we express the covariance as a sum over m and m’ of the square of the correlation $\langle {\tilde{a}}_{\ell m}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{*}\rangle$, $${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}=\frac{2}{\left(2\ell +1\right)\left(2{\ell }^{\prime }+1\right)}{\displaystyle \sum _{m{m}^{\prime }}{\left|\langle {\tilde{a}}_{\ell m}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{*}\rangle \right|}^2.}$$(7)

The harmonic coefficients correlation can be written as $$\begin{array}{lll}\langle {\tilde{a}}_{\ell m}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{*}\rangle \hfill & =\hfill & {\displaystyle \sum _{LM}{C}_L{}_0{I}_{\ell mLM}{}_0{I}_{{\ell }^{\prime }{m}^{\prime }LM}^{*},}\hfill \\ \hfill & =\hfill & {\displaystyle \int \text{d}\widehat{u}{}_0{Y}_{\ell m}\left(\widehat{u}\right)W\left(\widehat{u}\right){\displaystyle \sum _{LM}{\left\{C_L{}_0{I}_{{\ell }^{\prime }{m}^{\prime }LM}^{*}\right\}}_0{Y}_{LM}^{*}\left(\widehat{u}\right),}}\hfill \end{array}$$(8)

where we have used Eq. (2) to expand one of the ₀I kernels, reorganized the equation, and used the fact that the power spectrum C_L and the mask W are real quantities (we also dropped the explicit W dependences of the kernel to simplify notations). For fixed ℓ′ and m′, the rightmost part of the equation can be seen as the complex conjugate of the backward spherical harmonic transform of a set of spherical harmonic coefficients into a map X^ℓ′m′, defined as $$X^{{\ell }^{\prime }{m}^{\prime }}\left(\widehat{u}\right)\equiv {\displaystyle \sum _{LM}{x}_{LM}^{{\ell }^{\prime }{m}^{\prime }}{Y}_{LM}\left(\widehat{u}\right),}$$(9)

where we defined the spherical harmonic coefficients with $$x_{LM}^{{\ell }^{\prime }{m}^{\prime }}\equiv C_L{}_0{I}_{{\ell }^{\prime }{m}^{\prime }LM}.$$(10)

Here, we emphasize that the map X^ℓ′m′ is a complex map, thus it needs special care when it is decomposed into harmonic coefficients. With these definitions, Eq. (8) reduces to $$\langle {\tilde{a}}_{\ell m}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{*}\rangle ={\displaystyle \int \text{d}\widehat{u}{}_0{Y}_{\ell m}\left(\widehat{u}\right)W\left(\widehat{u}\right)X^{{\ell }^{\prime }{m}^{\prime }*}\left(\widehat{u}\right),}$$(11)

where we recognize the forward harmonic transform of the map X^ℓ′m′, masked by W. Thus, a spherical harmonic transform of a masked map can produce the correlation $\langle {\tilde{a}}_{\ell m}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{*}\rangle$ for all ℓ, m and a fixed pair of ℓ′, m′. As we discussed, this X^ℓ′m′ map is defined by a set of spherical harmonic coefficients whose expression is given in Eq. (10).

The computation of the $x_{LM}^{{\ell }^{\prime }{m}^{\prime }}$ coefficients requires the evaluation of the ₀I kernel. When Eq. (2) is used for a fixed ℓ′, m′, the ₀I_LMℓ′m′ kernel can be computed as the spherical harmonic transform of a masked ₀Y_ℓ′m′ map for all the L, M indices. When everything is added, we see that for a choice of ℓ′, m′, the computation of $\langle {\tilde{a}}_{\ell m}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{*}\rangle$ for all ℓ, m reduces to two forward and one backward spherical harmonic transforms, as summarized in Fig. 2.

In practice, these decompositions can be performed with HEALPix, which takes advantage of a specific pixelation scheme to make the computation more efficient. This is where the gain announced at the beginning of this section comes from, and it allows us to implement the exact computation. HEALPix decompositions typically scale as O(ℓ′^’3)⁶. We repeated the decompositions resulting in $\langle {\tilde{a}}_{\ell m}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{*}\rangle$ for all m′ ∊[−ℓ′,ℓ′] to perform the summation in Eq. (7). Thus, the computation of a single row ${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}$, for all ℓ and fixed ℓ′ scales as O(ℓ′⁴). Finally, the computation of a full covariance matrix for all ℓ′ scales as $O\left({\ell }_{\max ,\text{ex}}^5\right)$.

Additional optimizations can be implemented in the algorithm by degrading maps and running the algorithm at a lower HEALPix resolution, n_side, for small multipoles. HEALPix computations are precise up to ℓ ~ 2n_side, hence choosing a map resolution on the order of the multipole is sufficient to precisely compute the close-to-diagonal elements of the covariance. This operation requires a degraded version of the mask, which must be computed while avoiding aliasing from small-scale features. This can be done by implementing a hard cutoff of the mask harmonic coefficients before degrading its resolution. This allowed us to compute the exact covariance up to multipole ℓ_max,ex = 1000. The algorithm requires 300 h of CPU-time to compute a row of the intensity (TTTT) and polarization (EEEE) matrices at multipole ℓ = 950 with map resolution n_side = 1024. It is also well suited to a potential GPU implementation, which could lead to more speed-ups.

Fig. 2 Algorithm used to compute one row of the covariance ${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}$, using the HEALPix tools for a fixed ℓ′ and varying ℓ. Square, diamond, or circle boxes are arrays representing maps, power spectra, or spherical harmonic coefficients, respectively. Operations are symbolized by arrows and described alongside. The indices of the arrays are indicated on the side of the corresponding boxes. For example, near the top of the diagram, the HEALPIX.map2alm operation applied on the product Yl’m′ W produces the array 0ILMl′m′, with indices L, M. At the bottom of the diagram, the operations before the final summation produce the array $\langle {\tilde{a}}_{\ell m}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{*}\rangle$ with indices ℓ,m for a fixed ℓ’,m’ pair. This part of the algorithm scales as O(ℓ’3) because this is the scaling of HEALPix operations (which we applied three times) when the resolution of the map is comparable to the maximum multipole index considered. The last operation, summing over the indices m, m′, requires repeating the precedent steps for all m′ ∊ [−ℓ′, ℓ′], thus repeating them 2ℓ′ + 1 times. Therefore, the final complexity for producing a single row with fixed multipole index ℓ′ is O(ℓ′4). Computing the covariance matrix for all ℓ′ then increases the computational time to O(ℓ′5).

3.2 Polarization

The polarized case is very similar to the intensity case detailed in the previous subsection. We only describe the EEEE case, which gives a general template to the other polarization and temperaturexpolarization cases.

The polarized version of Eq. (8) is given by Eq. (6) of Challinor & Chon (2005), which reads $$\begin{array}{lll}\langle {\tilde{a}}_{\ell m}^{\text{E}}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{\text{E}*}\rangle \hfill & =\hfill & {\displaystyle \sum _{LM}[C_L^{\text{EE}}+I_{\ell mLM+}{I}_{{\ell }^{\prime }{m}^{\prime }LM}^{*}}\hfill \\ \hfill & \hfill & C_L^{\text{BB}}-I_{\ell mLM-}{I}_{{\ell }^{\prime }{m}^{\prime }LM}^{*}],\hfill \end{array}$$(12)

where we defined the Hermitian coupling coefficients $${}_{\pm }{I}_{\ell mLM}=\frac{1}{2}\left({}_{+2}{I}_{\ell mLM}\pm {}_{-2}{I}_{\ell mLM}\right),$$(13)

with the spin-weighted coupling coefficients ±₂I defined as for intensity, see Eq. (2). Reordering the terms and repeating the same operations as in the previous section, the final harmonic coefficient correlation can be seen as the masked forward spherical harmonic decomposition of two maps $Z_1^{{\ell }^{\prime }{m}^{\prime }},Z_2^{{\ell }^{\prime }{m}^{\prime }}$, $$\begin{array}{lll}\langle {\tilde{a}}_{\ell m}^{\text{E}}{\tilde{a}}_{{\ell }^{\prime }{m}^{\prime }}^{\text{E*}}\rangle \hfill & =\hfill & \frac{1}{2}[{\displaystyle \int \text{d}\widehat{u}W\left(\widehat{u}\right)Z_1^{{\ell }^{\prime }{m}^{\prime }}\left(\widehat{u}\right)\left({}_{+2}{Y}_{\ell m}^{*}{+}_{-2}{Y}_{\ell m}^{*}\right)\left(\widehat{u}\right)}\hfill \\ \hfill & \hfill & -i{\displaystyle \int \text{d}\widehat{u}W\left(\widehat{u}\right)Z_1^{{\ell }^{\prime }{m}^{\prime }}\left(\widehat{u}\right)\left({}_{+2}{Y}_{\ell m}^{*}-{}_{-2}{Y}_{\ell m}^{*}\right)\left(\widehat{u}\right)}].\hfill \end{array}$$(14)

The maps $Z_1^{{\ell }^{\prime }{m}^{\prime }},Z_2^{{\ell }^{\prime }{m}^{\prime }}$ are obtained using a backward spherical harmonic decomposition of the coefficients $x_{LM}^{{\ell }^{\prime }{m}^{\prime };\text{E,B}}$, $$\left(Z_1^{{\ell }^{\prime }{m}^{\prime }}-i{Z}_2^{{\ell }^{\prime }{m}^{\prime }}\right)\left(\widehat{u}\right)\equiv {\displaystyle \sum _{LM}\frac{x_{LM}^{{\ell }^{\prime }{m}^{\prime };\text{E}}+x_{LM}^{{\ell }^{\prime }{m}^{\prime };\text{B}}}{2}{}_{+2}{Y}_{LM}}\left(\widehat{u}\right),$$(15) $$\left(Z_1^{{\ell }^{\prime }{m}^{\prime }}+i{Z}_2^{{\ell }^{\prime }{m}^{\prime }}\right)\left(\widehat{u}\right)\equiv {\displaystyle \sum _{LM}\frac{x_{LM}^{{\ell }^{\prime }{m}^{\prime };\text{E}}-x_{LM}^{{\ell }^{\prime }{m}^{\prime };\text{B}}}{2}{}_{-2}{Y}_{LM}}\left(\widehat{u}\right).$$(16)

The set of harmonic coefficients $x_{LM}^{\text{E,B}}$ is defined similarly to the temperature case (Eq. (10)), and it was obtained by filtering the coefficients computed with a masked forward harmonic decomposition of the spin-2 spherical harmonics ±₂Y_{ℓ′, m′}, $$\{\begin{array}{l}{x}_{LM}^{{\ell }^{\prime }{m}^{\prime };\text{E}}=C_L^{\text{EE}}{}_{+}{I}_{{\ell }^{\prime }{m}^{\prime }LM,}\hfill \\ x_{LM}^{{\ell }^{\prime }{m}^{\prime };\text{B}}=C_L^{\text{BB}}{}_{-}{I}_{{\ell }^{\prime }{m}^{\prime }LM.}\hfill \end{array}$$(17)

This algorithm can be extended for any combination of spectra for the other polarization cases, including the cross-correlation between temperature and polarization that we do not treat here.

3.3 Validation on simulations

We compared the results of our implementation of the exact computation with a MC estimate of the covariance, obtained with N_sim simulations. The MC covariance terms are expected to be Wishart distributed with N_sim degrees of freedom, as explained in Lueker et al. (2010). We can estimate their variance to be $$\langle {\left({\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}^{\text{sim}}-\langle {\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}^{\text{sim}}\rangle \right)}^2\rangle =\frac{{\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}^2+{\widetilde{\text{Σ}}}_{\ell \ell }{\widetilde{\text{Σ}}}_{{\ell }^{\prime }{\ell }^{\prime }}}{N_{\text{sim}}}.$$(18)

N_sim = 10 000 allows us to reach a percent-level accuracy on the diagonal. This is the number of realizations that we use for the MC covariance. For this validation, we used the mask shown in Fig. 1. In this idealized setting, we did not include a point source mask.

We performed an exact computation of the TTTT and EEEE covariance up to ℓ_max,ex = 1000, using our algorithm and degrading the mask to smaller resolutions. Furthermore, we computed the rows every 25 multipoles of the matrix up to ℓ_max = 1500, as well as at a few well-chosen multipoles that correspond to peaks and troughs of the spectra up to ℓ_max = 2000.

We first focus on the diagonal of the covariance. Figure 3 presents the comparison between the exact computation of the diagonal, obtained by selecting the corresponding value in the rows we computed, and the MC evaluation. The two agree within the MC noise expected for N_sim.

For the off-diagonal terms, we show in Fig. 4 a few rows of the exact and MC covariance. The rows agree within the MC noise. The correlation between the modes falls to the percent level within a distance |ℓ – ℓ′| ~ 25 bands around the diagonal. The correlation matrix is defined as the covariance renormalized by its diagonal, $${\sigma }_{\ell {\ell }^{\prime }}\equiv \frac{{\sum }_{\ell {\ell }^{\prime }}}{\sqrt{{\sum }_{\ell \ell }{\sum }_{{\ell }^{\prime }{\ell }^{\prime }}}}.$$(19)

We display the exact and MC correlation matrix in the same multipole range in Fig. 5. Our tests demonstrate that our implementation of the exact computation of the pseudo-spectrum covariance is correct, at least at the level of accuracy that can be reached by MC estimation.

While massive MC estimates as we performed here in this idealized case can produce accurate of-diagonal term estimates, performing numerous MC estimates can be more challenging in the case of a realistic experiment. It requires regularization approaches, where the number of possible simulations is limited by the computational cost of mock-observations.

We also stress that our covariance matrix cannot be directly compared to the one used in Dutcher et al. (2021) because the presence of a point source mask (which we did not include in our simple example), the complications brought about by introducing realistic noise and scanning strategy, and projection effects (the analysis in Dutcher et al. (2021) is performed using a flat-sky approach) can all yield different levels of correlations between the modes.

Fig. 3 Relative difference of diagonals ${\widetilde{\text{Σ}}}_{\ell \ell }/{\widetilde{\text{Σ}}}_{\ell \ell }^{\text{sim}}-1$ for temperature TTTT (top) and polarization EEEE (bottom). In red we plot the relative differences every 25 multipoles until ℓ = 1500 and a few well-chosen ones (at the locations of peaks and troughs of the spectra) up to ℓ = 2000. In gray, the same quantity is plotted for all multipoles for ℓ ∊ [ℓcut = 200, ℓmax,ex = 1000]. We are able to compute the covariance exactly only for a limited number of rows, and it is computationally cheaper for lower multipoles, justifying our choice of full calculation at ℓ < ℓmax,ex and partial calculation for larger multipoles. This plot shows the agreement between the two approaches and validates our exact calculation.

Fig. 4 Rows of the exact covariance (solid red) and of the simulated one (dashed blue) for multipole ℓ′ = 1000 (left-hand side) and ℓ′ = 1800 (right-hand side). The dotted black line shows the expected MC uncertainty on the simulated covariance from Eq. (18), using Nsim = 10 000 realizations. The two approaches agree very well. The MC standard deviation is as large as the covariance off-diagonal terms at ∣ℓ − ℓ′∣ ~ 25.

Fig. 5 EEEE correlation matrix (see Eq. (19)) obtained from the exact (left) and MC (right) calculations of the covariance matrix. The exact computation allows us to have the full correlation matrix, but the MC approach is limited by numerical noise. All terms below 10−4 are plotted in dark blue. Results for TTTT are comparable.

4 Existing approximations and their accuracy on a small patch of the sky

Our algorithm allows us to obtain the exact covariance only for ℓ < 1000 or for a few rows at higher ℓ′s due to the expensive computing resources required. The usual analytical approach consists of using approximations of Eq. (6) to decrease the computational cost. In this section, we introduce a new framework to express the approximations of the covariance matrix and use it to list the different methods proposed in the literature. Then, we test and discuss their accuracy against our exact computation.

4.1 General framework

Before discussing the approximations of the covariance, we define a few quantities that help relate the various approximations to each other. We rewrite Eq. (6) as $${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}=\frac{2}{\left(2\ell +1\right)\left(2{\ell }^{\prime }+1\right)}{\displaystyle \sum _{{\ell }_1{\ell }_2}{C}_{{\ell }_1}{\text{Θ}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}}\left[W\right]C_{{\ell }_2},$$(20)

introducing ${\text{Θ}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\left[W\right]$. The covariance coupling kernel, defined as the sum over the multipole orders (m, m′,…) of the coupling coefficients, reads $${\text{Θ}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\left[W\right]\equiv {\displaystyle \sum _{m{m}_1{m}^{\prime }{m}_2}\left({}_0{I}_{\ell m{\ell }_1{m}_1}{}_0{I}_{{\ell }^{\prime }{m}^{\prime }{\ell }_1{m}_1}^{*}{}_0{I}_{{\ell }^{\prime }{m}^{\prime }{\ell }_2{m}_2}{}_0{I}_{\ell m{\ell }_2{m}_2}^{*}\right)\left[W\right].}$$(21)

The covariance coupling kernel Θ represents the coupling between the modes of the theoretical underlying power spectrum C_ℓ′ for i = 1, 2 depending on the index of the pseudo-covariance (ℓ, ℓ′ <). We chose to show the two indices related to the covariance (ℓ, ℓ′) as subscripts and the two summing indices (ℓ₁, ℓ₂) as superscripts. We considered a single-mask temperature case, for which the coupling kernel is symmetric with respect to the exchange of multipole indices ℓ ↔ ℓ′ or ℓ₁ ↔ ℓ₂. In the following, we write our results in this case for the sake of simplicity, but they are valid regardless of the choice of single or multiple masks. In the case of spectra obtained from maps with different masks, or in the case of cross-spectra, the kernel is not symmetric. While the results of this work apply to both cases, considering multiple masks increases computing cost.

Using the completeness relations of spherical harmonics, given in Eqs. (A.10) and (A.11), we can write $${\displaystyle \sum _{{\ell }_1{\ell }_2}{\text{Θ}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}}\left[W\right]=\left(2\ell +1\right)\left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[W^2\right].$$(22)

We can now define the reduced covariance coupling kernel as $${\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\left[W\right]\equiv \frac{{\text{Θ}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\left[W\right]}{\left(2\ell +1\right)\left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[W\right]},$$(23)

for which $${\displaystyle \sum _{{\ell }_1{\ell }_2}{\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}}\left[W\right]=1.$$(24)

With these notations, we can rewrite the covariance as $${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}=2{\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[W^2\right]{\displaystyle \sum _{{\ell }_1{\ell }_2}{C}_{{\ell }_1}{\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}}\left[W\right]C_{{\ell }_2}.$$(25)

The symmetric mode-coupling kernel Ξ[W²] provides the purely geometric coupling due to sky masking and is common to all approximations of the covariance. It only depends on the power spectrum of the squared mask, and it is easy to be convinced that up to a normalization, it corresponds to the exact covariance in the case of a constant power spectrum. Its computation scales as $O\left({\ell }_{\max }^3\right)$. This could be improved by noting, as was done by Louis et al. (2020), that at small enough scales, Ξ[W²] is close to a Toeplitz matrix, allowing us to further reduce the scaling to $O\left({\ell }_{\max }^2\right)$ for a wide range of modes. The sum on the right-hand side of Eq. (25) describes the contribution of the signal power spectrum modulated by the kernel $\overline{\text{Θ}}$, which depends on the mask. This is the sum that all approximations try to simplify, replacing the kernel $\overline{\text{Θ}}$ with a simpler ansatz. In the following, we describe all approximations in terms of this redefinition of the covariance matrix. Since the reduced covariance coupling kernel is normalized, every approximation formulated in this formalism yields the exact covariance for a constant underlying power spectrum C_ℓ = N.

4.2 Approximations

4.2.1 Narrow-kernel approximation

Based on the observation that the coupling coefficients ₀I in Eq. (1) are narrow and peak at their first multipole indices ℓ or ℓ′, Efstathiou (2004) introduced the following approximation of Eq. (25), taking the convolving spectra $C_{{\ell }_i,}i=1,2$ out of the sum, and replacing them by the power spectrum evaluated at the first multipole index of the coupling coefficients (i.e., the covariance indices of $\overline{\text{Θ}}$), C_ℓ or C_ℓ′. Following the notation introduced in García-García et al. (2019), we refer to this approximation of the covariance as the narrow-kernel approximation (NKA), $$\begin{array}{lll}{\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}\hfill & \approx \hfill & 2C_{\ell }{C}_{{\ell }^{\prime }}{\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[W^2\right]{\displaystyle \sum _{{\ell }_1{\ell }_2}{\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\left[W\right]}\hfill \\ \hfill & =\hfill & 2C_{\ell }{C}_{{\ell }^{\prime }}{\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[W^2\right]\equiv {\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}^{\text{NKA}}.\hfill \end{array}$$(26)

In terms of the reduced covariance coupling kernel, the NKA uses $${\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\left[W\right]\approx {\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2;\text{NKA}}\left[W\right]\equiv \frac{{\delta }_{\ell {\ell }_1}{\delta }_{{\ell }^{\prime }{\ell }_2}+{\delta }_{{\ell }^{\prime }{\ell }_1}{\delta }_{\ell {\ell }_2}}{2}.$$(27)

The approximation is exact for the full sky. It provides an accurate estimator whenever the underlying power spectrum C_ℓ varies slowly as a function of ℓ compared to the typical size of the operators ₀I. This condition is fulfilled when the amplitude of the mask power spectrum drops quickly with multipole ℓ, which is the case for large sky fractions observed with a mask that contains no small-scale features. This is shown in Fig. 1, for example, where we plot the power spectrum of one of the masks used in the Planck analysis. In this case, the above approximation holds for multipoles much larger than those for which the mask spectrum contains power.

The NKA was first introduced in intensity by Efstathiou (2004) and extended to polarization in Challinor & Chon (2005). As in the temperature case, the approximated covariances in polarization are expressed as a function of the polarization spectra EE and BB and the symmetric coupling kernels ${\text{Ξ}}_{\ell {\ell }^{\prime }}^{{}_{\pm }2{,}_{\pm }2}$. The expressions of the approximated polarization covariances mix EE and BB due to leakage that appears because the sky is masked; see Eqs. (25)–(27) of Challinor & Chon (2005).

The NKA has been widely used, for instance, in the Planck cosmological analysis, which masked only small portions of the full sky; see Planck Collaboration XI (2016). However, it has never been thoroughly tested on small sky fractions. As shown in Fig. 1, the mask power spectrum in the case of the small survey footprint of SPT-3G drops much more slowly than the large Planck one. From this observation, we expect the mode-coupling kernels _sI to be wider, as can be deduced from Eq. (2). As a result, the theoretical underlying spectrum C_ℓ might not be treated as constant compared to the covariance coupling kernels in the sums of Eq. (25), and SPT-3G may be outside the regime of validity of the NKA assumption. This is tested at the end of this section. We now list some proposed improvements to the NKA.

4.2.2 Friedrich approximation

A straightforward extension of the NKA has been proposed in Friedrich et al. (2021). It is based on the observation that the reduced covariance-coupling kernel $\overline{\text{Θ}}$ has four maxima at [ℓ = ℓ₁,ℓ′ = ℓ₁], [ℓ = ℓ₁, ℓ′ = ℓ₂], [ℓ = ℓ₂, ℓ′ = ℓ₁], and [ℓ = ℓ₂, ℓ′ = ℓ₂]. This suggests the following form of the reduced covariance coupling matrix: $${\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\left[W\right]\approx {\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2;\text{FRI}}\left[W\right]\equiv \frac{{\delta }_{\ell {\ell }_1}+{\delta }_{{\ell }^{\prime }{\ell }_1}}{2}\frac{{\delta }_{\ell {\ell }_2}+{\delta }_{{\ell }^{\prime }{\ell }_2}}{2}.$$(28)

Thus, the approximated covariance is $${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}^{\text{FRI}}\equiv 2{\Xi }_{\ell {\ell }^{\prime }}^{00}\left[W^2\right]{\left(\frac{C_{\ell }+C_{{\ell }^{\prime }}}{2}\right)}^2.$$(29)

We refer to this approximation as FRI in the rest of the article.

Fig. 6 Relative differences of binned approximations with respect to the exact binned covariance: ${\widetilde{\text{Σ}}}_{b{b}^{\prime }}^{\text{APP}}/{\widetilde{\text{Σ}}}_{b{b}^{\prime }}-1$ for TTTT (left) and EEEE (right), with binning Δℓ = 50. In the first row, we plot the relative differences for the diagonal, i.e., b = b′, while in the second row, we plot those of the first off-diagonal, i.e., b′ = b + 1. The NKA (light blue dashed), FRI (purple dashed-double-dot) and INKA (dark blue dashed-dot) approximations are accurate at the 5% level, whereas the ACC approximation (solid red) reaches the 1% level, both in intensity and polarization for multipoles larger than ℓcut = 200. The relative differences are plotted for bins that include multipoles up to ℓmax,ex = 1000 because it is the maximum multipole for which we computed all the rows of the exact covariance. The third row displays the corresponding binned underlying renormalized spectrum TT or EE to show that the difference of the covariances lies in the peaks and troughs of the spectra.

4.2.3 Improved narrow-kernel approximation

Nicola et al. (2021) proposed an improved version of the NKA, the improved narrow-kernel approximation (INKA). In this approximation, the Dirac functions in Eq. (27) are replaced by ₀M, the renormalized MASTER mode-coupling kernel, as defined in Appendix A.3. It reads $${\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\left[W\right]\approx {\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2;\text{FRI}}\left[W\right]\equiv \frac{{}_0{\overline{M}}_{\ell {\ell }_1}{}_0{\overline{M}}_{{\ell }^{\prime }{\ell }_2}+{}_0{\overline{M}}_{{\ell }^{\prime }{\ell }_1}{}_0{\overline{M}}_{\ell {\ell }_2}}{2}.$$(30)

The convolution in Eq. (6) indeed averages the power spectra $C_{{\ell }_i,}i=1,2$ over multipoles close to ℓ and ℓ′. We can take advantage of this by replacing the convolution by a multiplication with a smoothed power spectrum. When $\overline{C}\equiv \overline{M}C$, the resulting covariance can be written as $${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}^{\text{INKA}}\equiv 2{\overline{C}}_{\ell }{\overline{C}}_{{\ell }^{\prime }}{\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[W^2\right].$$(31)

All the NKA, FRI, and INKA scale as $O\left({\ell }_{\max }^3\right)$, which are the computing resources needed to obtain the coupling kernels Ξ and $\overline{M}$. This is a significant improvement over the $O\left({\ell }_{\max }^5\right)$ scaling of our full computation. We now validate the approximations in the case of small surveys using the expensive exact computation of the covariance matrix.

4.3 Accuracy

We tested the accuracy of the NKA, FRI and INKA using the exact computation in the case of the SPT-3G small survey footprint shown in Fig. 1. For this mask, the correlations between modes are significant, as we showed in Fig. 4. In this case, it is customary to bin the individual multipoles into wider bandpowers. For this reason, we performed all of our tests on a binned version of the covariance. Given the shape of the power spectrum of the mask and the correlations that we expect from it, we adopted a Δℓ = 50 binning with ℓ(ℓ + ℓ)/(2π) weights to flatten the dynamics of the spectra in each bin. With this bin size, we expect that most of the correlations between bandpowers are concentrated in the first off-diagonal bin. We also conservatively excluded the first ℓ_cut = 200 multipoles from our analysis. They are more challenging to measure on a small survey footprint as they can suffer from leakage from the super-survey scales. We restrict our comparison to the multipoles between ℓ_cut and ℓ_max,ex = 1000, where we have carried out the exact computation of all the matrix rows.

We present in Fig. 6 a comparison between the exact computation and the NKA, FRI, and INKA for the diagonal and first off-diagonal of the TTTT and EEEE binned covariances. We discuss the performance of our new accurate covariance coupling (ACC) approximation, also shown in the figure, in the next section. The existing approximations provide good estimates of these elements of the covariance as they fall within 5% of the accuracy. The amplitude of the errors varies at different multipoles. Even though the FRI and INKA schemes were implemented to improve upon the simple NKA, their errors are of similar amplitude for this choice of binning. However, all approximations fail to recover the binned covariance at the percent level. On the third row, we show that the difference of the covariances lies in the peaks and troughs of the spectra. As the covariance coupling kernel acts as a symmetric convolution on the spectrum (see Eq. (25)), it is sensitive to the second-order derivative of the spectrum rather than the first. An error on the covariance coupling kernels leads to a larger relative difference in the covariance at multipoles where the curvature of the spectrum is maximum. In Sect. 5 we use the knowledge gained from the exact computation to propose an improved approximation scheme.

Because we cannot easily compute the full matrix to present binned results, we only compare some unbinned rows in Fig. 7 at higher multipoles. The shaded regions in this figure give the lowest values of the relative difference for the approximation within multipoles ℓ ∊ [ℓ_cut = 200, ℓ_max,ex = 1000] and ℓ′ ∊ [ℓ − 2Δℓ, ℓ + 2Δℓ]. These are the covariance terms for which we can calculate the full binned covariance; their accuracy is shown in Fig. 6.

Furthermore, the lines in Fig. 7 show the same quantity as the shaded regions, that is, the maximum relative difference, but for multipoles ℓ ∊ [ℓ_max,ex, 2000], estimated over the sparse number of rows for which we computed the matrix exactly. The difference with the exact covariance for all approximations at large multipoles is always within the same error range as for lower multipoles. This shows that the approximations still work with the same precision at higher multipoles, both for temperature and polarization, and that the accuracy of the approximations in the unbinned case quickly falls below 20% when Δ = 50.

Fig. 7 Relative difference between the unbinned approximated covariance matrices compared to the exact calculation, ${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}^{\text{APP}}/{\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}-1$, as a function of Δ = ℓ − ℓ′. We show the NKA (light blue), INKA (dark blue), and FRI (purple) approximations for TTTT (top) and EEEE (bottom). We only show the relative differences for the ℓ at which the deviation is largest, with ℓ ∊ [ℓcut = 200,ℓmax,ex = 1000] (shaded regions) or ℓ ∊ [ℓmax,ex, 2000] (lines).

Fig. 8 Algorithm for computing the reduced covariance -coupling kernels using HEALPix tools. We use the same notation as in the diagram of Fig. 2. The HEALPix functions require $O\left(n_{\text{side}}^3\right)$, with nside the chosen resolution of maps W and Yℓm. As we choose the resolution nside to be on the order of the multipoles indices ℓ,ℓ′, it is equivalent to say that they require O((ℓ + ℓ′)3). As operations are done O(ℓ + ℓ′) times, the whole operation of computing ${\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}$, is O((ℓ + ℓ′)4). Finally, it is clear in this diagram that the computing time of this kernel is at least doubled when multiple masks are used. Different masks would be used as inputs in the first line. As a result, the coupling coefficients would need to be computed for each of the masks, as shown in Eq. (D.3).

4.4 Structure of the reduced covariance-coupling kernel

Our expression of the covariance matrix approximations in terms of the normalized coupling kernel $\overline{\text{Θ}}$ in Eq. (27) gives us a very efficient tool for examining the validity of each approximation and for better understanding their differences. We designed an algorithm to calculate this kernel exactly, similar to the one described in Sect. 3 for the exact calculation of the matrix. We show a diagram of the algorithm in Fig. 8.

The reduced covariance-coupling kernel is then displayed in Fig. 9 for the INKA, NKA, and FRI approximations compared to the exact computation for a fiducial multipole ℓ = 200. The kernels are represented as matrices as a function of ℓ₁, ℓ₂ for different fixed choices of the indices ℓ, ℓ′. Columns show the results for different choices of (ℓ′ = ℓ − Δ, with Δ = 0 (i.e., the kernels for the diagonal terms of the covariance matrix, e.g., ${\widetilde{\text{Σ}}}_{200,200}$) or Δ = 10, 50 (i.e., the kernels for the off-diagonal terms separated by 10 or 50 multipoles, e.g., ${\widetilde{\text{Σ}}}_{200,190}$). We recall that the reduced kernel is multiplied to $C_{{\ell }_1}$, $C_{{\ell }_2}$ and summed over the indices ℓ,ℓ₂ in Eq. (25). Hence, Fig. 9 directly shows the weight of the ℓ₁, ℓ₂ power spectra that contribute to the ${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}$ element of the covariance matrix.

We first focus on the kernels for the diagonal of the covariance matrix, Δ = 0, shown in the first column of Fig. 9. All kernels peak at ℓ₁ = ℓ₂ = ℓ, as expected. However, it is clear from the exact calculation that the kernel has a significant width compared to the CMB power spectrum. This is more clearly shown in Fig. 10, where we plot a slice of the coupling kernels for ℓ₁ = 200. The width of the kernel cannot be neglected compared to the slope of the CMB power spectrum. This justifies the INKA, which replaces the Dirac δ functions of NKA and FRI by renormalized mode coupling kernels, see Eq. (31). However, as shown in this figure, the INKA kernel is slightly broader than the exact calculation and its amplitude is smaller. This explains why INKA underestimates the covariance diagonal in the peaks of the power spectrum and overestimates it on the troughs, as shown in Fig. 6: it averages the underlying power spectrum in a wider range of multipoles. Conversely, the NKA/FRI kernels are much thinner than in the exact computation, and so they overestimate the diagonal in the peaks and underestimate it in the troughs of the power spectrum.

Second, we focus on the off-diagonal terms, Δ = 10, 50, shown in the second and third column of Fig. 9. The difference between the exact computation and all the existing approximations is striking, and it is clear that the kernel has more structure than the simple approximated forms. For close off-diagonal terms such as Δ = 10, the true kernel peaks at its central index ${\ell }_1={\ell }_2=\overline{\ell }\equiv \left(\ell +{\ell }^{\prime }\right)/2$. For far off-diagonal terms such as Δ = 50, there are four maxima, as predicted by the FRI approximation, which are partially missed by the INKA. Moreover, the true coupling has more dynamics and also covers negative values. Therefore the different approximations, even the INKA approximation, fail to correctly represent the off-diagonal terms of the covariance, as observed in Fig. 7.

Fig. 9 Reduced covariance coupling kernels ${\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}$ for ℓ = 200 and ℓ’′= ℓ + Δ, with Δ = [0, 10, 50] shown in the three different columns. All kernels in one column share the same color map. All of the matrices are shown as a function of ℓ1 and ℓ2 and are centered on = (ℓ + ℓ′)/2. Each of the displayed kernels is properly normalized according to Eq. (22). The plots are restricted to the elements that have a significant value. Top row: approximated INKA kernels and the positions at which the delta distributions of the NKA (circles) and FRI (crosses) approximation peak. Bottom row: exact kernels. The comparison between the two highlights how much of the structure of the exact kernel is missed by the different approximations.

Fig. 10 Slices of the exact covariance coupling kernel (solid red) vs. the approximated NKA and FRI (dashed light blue) and INKA (dot-dashed dark blue) ones (right scale). We show for comparison the CMB intensity power spectrum (left scale, solid gray line).

Fig. 11 Diagonal (top) and row (bottom) of the reduced covariance coupling kernels ${\overline{\text{Θ}}}_{\ell ,{\ell }^{\prime }}^{{\ell }_1,{\ell }_2}$ for ℓ ∈ [150, 206, 300, 650, 750] and ∆ = ℓ′ − ℓ = 50 as a function of ∆2 = ℓ2 − ℓ. The row (bottom) is plotted for ℓ1 = ℓ, i.e. ∆1 = ℓ1 − ℓ = 0. The plots show that for different I but the same ∆, the kernels are very similar, differing only at the 5% level. A similar result can be shown for other values of ∆. This leads us to formulate our new approximation, where we assume $\overline{\text{Θ}}$ to depend only on the multipole separations ∆, ∆1, ∆2.

5 New approximation for the covariance

5.1 Improved approach: Approximated covariance coupling

Our ability to calculate the exact reduced covariance coupling matrix $\overline{\text{Θ}}$, described in Sect. 4, allows us to introduce a new approximation for the computation of the pseudo-power spectrum covariance matrix. We note that for a fixed ∆ = ℓ′ − ℓ, the structure of ${\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}$ appears to be invariant. In other words, the coupling matrices contributing to the ${\widetilde{\text{Σ}}}_{\ell {\ell }^{\prime }}$, term of the covariance matrix only depend on the distance ∆ from the diagonal. This is demonstrated in Fig. 11, where we plot diagonal and rows of the exact calculation of ${\overline{\text{Θ}}}_{\ell ,\ell +\text{Δ}}\text{for\hspace{0.17em}Δ}=50$ and different ℓ, for ℓ₁ = ℓ. The plot reveals that the kernels are nearly identical when they are plotted as a function of ∆₂ = ℓ₂ − ℓ. We thus infer that in general, when the 0 matrices are written as a function of ∆₁ = ℓ₁ − ℓ and ∆₂ = ℓ₂ − ℓ, they only depend on ∆ = ℓ′ − ℓ for any ℓ and ℓ′. The difference between kernels computed at different ℓ_S for same ∆ is small, at the 5% percent level. We can thus assume that for any choice of multipole ℓ, λ, $$\forall \left(\ell ,\lambda \right),{\overline{\text{Θ}}}_{\ell \left(\ell +\text{Δ}\right)}^{\left(\ell +{\text{Δ}}_1\right)\left(\ell +{\text{Δ}}_2\right)}\approx {\overline{\text{Θ}}}_{\lambda \left(\lambda +\text{Δ}\right)}^{\left(\lambda +{\text{Δ}}_1\right)\left(\lambda +{\text{Δ}}_2\right)}.$$(32)

An analytical justification of this approximation is provided in Appendix E using the asymptotic expansion of the Wigner-3j symbols when ℓ is large.

This suggests a new approximation where the coupling kernel just has to be computed at a given fiducial ℓ for all relevant values of ∆, ∆₁, and ∆₂. We call this the new approximated covariance coupling method, ACC. More precisely, the ACC kernel is given by $${\overline{\text{Θ}}}_{\ell {\ell }^{\prime }}^{{\ell }_1{\ell }_2}\approx {\overline{\text{Θ}}}_{\ell \left(\ell +\text{Δ}\right);\text{ACC}}^{\left(\ell +{\text{Δ}}_1\right)\left(\ell +{\text{Δ}}_2\right)}\equiv {\overline{\text{Θ}}}_{\ell *\left(\ell *+\text{Δ}\right)}^{\left(\ell *+{\text{Δ}}_1\right)\left(\ell *+{\text{Δ}}_2\right)},$$(33)

where we perform the exact and costly computation only for the ${\overline{\text{Θ}}}_{\ell *\left(\ell *+\text{Δ}\right)}^{\left(\ell *+{\text{Δ}}_1\right)\left(\ell *+{\text{Δ}}_2\right)}$. We are free to choose the reference ℓ* multipole. Because there are no significant long-range correlations in our case, we can pick a low⁷ ℓ* and use a low n_side map resolution. We have to ensure, however, that ℓ* is larger than ℓ_cut, the low-ℓ cutoff that was introduced to avoid issues with large-scale leakages. Close to ℓ_cut, the exact computation can be used. With a small mask, large-scale modes are difficult to measure and are usually excluded from the cosmological analysis. We can also restrict the range of ∆ to the number of off-diagonal terms of interest in the covariance matrix. The correlation falls quickly (see Figs. 4 and 5), and in practice, we can restrict it to |∆| < d_max, with d_max being on the order of a few times the correlation length. Similarly, the kernels fall quickly in ∆₁, ∆₂, so that we can also restrict ourselves to a small region of a similar order, and in the case of the single-mask analysis, use the symmetry around ∆ ↔ −∆ to reduce the computational cost.

While we only presented temperature coupling kernels in Fig. 9, the situation is identical in polarization, and a similar approximation can be built; see Appendix E. We used this approximation with ℓ* = 300, n_side = 512, and d_max = 100 to compute the ACC results in Fig. 6.

5.2 Accuracy and scaling

We validate the accuracy of our ACC approximation and compare it to the other approximations in Fig. 6. Our new approximation succeeds at estimating the covariance within an error of 1% for all multipoles larger than ℓ_cut in intensity and polarization, which is an improvement of a factor of ~4 over previous approximations. This is also shown in Figs. 7 and 12, which is just the same as Fig. 7, but focused on the EEEE ACC residuals with respect to INKA. This figure shows that the ACC approximation estimates both the diagonal and off-diagonal terms of the covariance matrix far better.

Figure 8 shows the computations needed to obtain the covariance-coupling kernels. Following the same argumentation as for the exact computation of the Sect. 3, we can show that the computation of a single kernel 0 scales as O((ℓ + ℓ′)⁴). As a result, because we need to compute one for each diagonal index ∆ ∈ [0, d_max], the final ACC approximation scales as $O\left(n_{\text{side}}^4{d}_{\max }\right)$, where n_side is the map resolution chosen to compute the kernels. The computing resources needed to obtain $\overline{\text{Θ}}$ for all approximations are summarized in Table 1. They add up to the resources needed to compute the symmetric coupling kernel Ξ in Eq. (25). In practice, the kernel $\overline{M}$ needed to build INKA is often already known for the sky analysis because it has the same structure as Ξ, so that the effective complexity for this approximation is O(1).

5.3 Point-source mask

We did not include a point-source mask in the survey footprint. Point-source masks significantly complicate the problem as the power spectrum of the mask will have power at large multipoles, hence it will extend the correlation length. This has been an issue for all analysis thus far. Apodizing the point-source masks helps to alleviate the problem, but at the price of discarding a significant area of the usable sky. Even in the case of a large survey footprint, such as Planck, the point-source masks have been shown to break the NKA. In this case, the issue was mitigated using a simulation-based correction. For FRI, INKA, and ACC, preliminary work that we have performed also suggests that they fail when sources are included. We expect these approximations to perform poorly because the improvements over NKA are focused on the central shape of the reduced covariance coupling matrix $\overline{\text{Θ}}$, while the sources tend to affect the far off-diagonal terms, which increases the correlation between distant multipoles. Particularly for the ACC approximation, the stronger correlations at distant multipoles will break the asymptotic behavior of the Wigner-3j symbols shown in Appendix E. We can expect that this reduces the validity of the approximated invariance by translation along the covariance diagonal, which is at the core of the ACC approximation. More work is required to assess the accuracy of ACC and other approximations in this case. However, different approaches can be adopted to mitigate the effect of a point-source mask. For example, analytical solutions might be found (Gratton et al., in prep.), the maps might be inpainted (Benoit-Lévy et al. 2013), or a MC correction such as the one used in Planck might be employed (Planck Collaboration XI 2016).

Fig. 12 Zoom of Fig. 7. We focus on the relative differences of the ACC and INKA covariance matrices with respect to the exact covariance for EEEE. The deviations found at ℓ < 1000 (shaded regions) are similar to those found at the higher multipoles (lines), showing that the approximations work at the same level of accuracy in the two cases.

6 Covariance of the PolSpice estimator

The pseudo-power spectrum is a biased estimator of the true underlying spectrum of the masked CMB maps. To recover an unbiased estimator, we can apply the MASTER (Hivon et al. 2002) formalism, which inverts the mode-coupling matrix and applies it to the biased estimator. Similarly, we can use the PolSpice (Szapudi et al. 2001; Chon et al. 2004) algorithm, which corrects the two-point correlation function for the effect of the mask in real space and then converts the result back into harmonic space. However, when the sky footprint is small and no large angular scales are observed, the unbinned mode-coupling matrix becomes singular. Analogously, the PolSpice conversion of the two-point correlation function into a power spectrum cannot be performed. In the first case, the mode-coupling matrix must be binned to allow the inversion. In the second case, we must apodize (i.e., cut gradually) the large angular scales of the two-point correlation function before calculating the corresponding power spectrum. This introduces a small bias in the final estimator that cannot be corrected for.

In this section, we explain in detail how the covariance matrix is calculated for the PolSpice estimator starting from a pseudo-power spectrum covariance matrix, which we produce through our ACC approximation. We show how the effect of the correction of the mask is included, as well as the small bias introduced by the PolSpice apodization of the two-point correlation function. In particular, we show that this apodization can be expressed in harmonic space, allowing us to relate the PolSpice spectrum covariance matrix to that of the pseudo-spectrum with a convolution.

6.1 MASTER equation

The pseudo-power spectrum is related to the true spectrum through the well-known MASTER equation introduced in Hivon et al. (2002), $$\langle {\tilde{C}}_{\ell }^{\text{TT}}\rangle ={\displaystyle \sum _{{\ell }^{\prime }}{}_0{M}_{\ell {\ell }^{\prime }}{C}_{{\ell }^{\prime }}^{\text{TT}},}$$(34)

with similar equations for polarization; see Appendix A.2. This bias comes from the information that is missing due to the masked sky. Given the weighted mask $W\left(\widehat{n}\right)$, we can compute ₀M using Eq. (A.15). Provided that ₀M is invertible, an unbiased estimator can be constructed. These relations can be expressed in real space using the two-point correlation functions ξ, which for a statistically isotropic sky depend only on the relative angle between two directions, $$\langle T\left({\widehat{n}}_1\right)T\left({\widehat{n}}_2\right)\rangle =\xi \left(\arccos \left({\widehat{n}}_1·{\widehat{n}}_2\right)\right).$$(35)

They can be related to the power spectrum C_ℓ using a Legendre series, with $$\xi \left(\theta \right)={\displaystyle \sum _{\ell }\frac{2\ell +1}{4\pi }{C}_{\ell }{P}_{\ell }\left(\cos \theta \right),}$$(36) $$C_{\ell }=2\pi {\displaystyle {\int }_{-1}^1\text{d}\cos \theta \xi \left(\theta \right)P_{\ell }\left(\cos \theta \right).}$$(37)

When we define the correlation function $\tilde{\xi }$ of the masked sky in the same manner, associated with the pseudo-spectrum ${\tilde{C}}_{\ell }$, we obtain from Eq. (34) by applying the decomposition in a Legendre series the following relation: $$\begin{array}{cc}\langle \tilde{\xi }\left(\theta \right)\rangle =w\left(\theta \right)\xi \left(\theta \right),& \forall \theta \in \left[0,\pi \right],\end{array}$$(38)

where w(θ) is the mask angular correlation function (more details can be found in Appendix C). From this relation, we can establish that the MASTER mode-coupling matrix ₀M is invertible only when the correlation function of the mask w(θ) is nonzero for all θ ∈ [0, π], which implies that the survey area explores all angular separations on the sky. While this is valid for an almost full-sky analysis such as Planck, it does not hold for experiments observing small patches, such as SPT-3G, where angular scales larger than θ ~ 30 deg are unexplored. As a result, ₀M is not invertible. Binning allows the regularization of the MASTER matrix and thus to build a nearly unbiased estimator of the bandpowers. This approach is described in Hivon et al. (2002) and is adopted in NaMaster⁸ (Alonso et al. 2019). Similarly, we show in the next section how the unobserved large angular scales are handled in the PolSpice estimator.

6.2 Regularizing with PolSpice

6.2.1 Temperature

The pseudo-power spectrum estimator can be regularized in real space following the PolSpice approach in Szapudi et al. (2001). The pseudo-correlation function ξ is smoothed with a scalar apodizing function f^apo(θ), which cuts out large θ and then corrects for the bias coming from the weighted mask described in Eq. (38). The scalar apodizing function goes smoothly from f^apo(0) = 1 to f^apo(θ_max) = 0 to avoid Fourier ringing. θ_max should be chosen as the maximum angular size of the weighted mask. A new correlation function estimator ξ(θ) is defined as $$\widehat{\xi }\left(\theta \right)\equiv g\left(\theta \right)\tilde{\xi }\left(\theta \right),$$(39)

with $$g\left(\theta \right)=\{\begin{array}{ll}{f}^{\text{apo}}\left(\theta \right)/w\left(\theta \right)\hfill & \forall \theta \in \left[0,{\theta }_{,\max }\right),\hfill \\ 0\hfill & \forall \theta \in \left[{\theta }_{\max },\pi \right].\hfill \end{array}$$(40)

The function g is well defined and smooth for all angles through the apodization f^apo. As a consequence of Eqs. (38) and (39), the PolSpice estimator of the correlation function can be related on average to the true underlying correlation function with $$\begin{array}{cc}\langle \widehat{\xi }\left(\theta \right)\rangle =f^{\text{apo}}\left(\theta \right)\xi \left(\theta \right)& \forall \theta \in \left[0,\pi \right].\end{array}$$(41)

Returning to harmonic space using a Legendre transform, this operation can be expressed as $$\langle {\widehat{C}}_{\ell }\rangle ={\displaystyle \sum _{{\ell }^{\prime }}{}_0{K}_{\ell {\ell }^{\prime }}{C}_{{\ell }^{\prime }}.}$$(42) $${}_0{K}_{\ell {\ell }^{\prime }}=\left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[f^{\text{apo}}\right].$$(43)

The PolSpice kernel ₀K is obtained from the scalar apodizing function f^apo with an extended definition of the operators Ξ (see Appendix C and Eq. (C.7) for more details). The operator acts on the Legendre transform of f^apo.

The advantage of PolSpice, which performs the regularization in real space rather than in harmonic space, is that it replaces an t-space convolution by an integration and a multiplication, which are faster and numerically more stable. This produces an estimator for all multipoles ℓ. We denote this estimator with a hat, for instance, ${\widehat{C}}_{\ell }^{\text{XY}}$. This regularization (which is only required for small sky patches) introduces a bias in the PolSpice estimator that cannot be corrected for because the coupling is not invertible generally as the apodizing function f^apo reaches zero. The bias is small because ₀K is properly normalized, that is, ∑_ℓ0K_ℓℓ′ = 1. Furthermore, the regularization increases the correlations between unbinned modes. The PolSpice kernels behave as window functions, mixing multipoles of the pseudo-power spectrum. The lack of information at large scales induces the inability to distinguish multipoles that are close to each other. For this reason, the spectrum estimator is binned in ranges larger than the typical correlation between multipoles for a cosmological analysis.

6.2.2 Polarization

PolSpice allows correcting for the bias introduced by the cut sky in the same manner as for the polarized spectra. It also allows decoupling the EE and BB estimator; see Challinor & Chon (2005) or Appendix C. Similarly to the intensity case, we can express the effect of the PolSpice real-space regularization in spherical harmonics by defining the polarized PolSpice kernels, ${}_{\pm 2}{K}_{\ell {\ell }^{\prime }}=\left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{2\pm 2}\left[f^{\text{apo}}\right]$. The PolSpice estimator follows for X ∈ [E, B] $$\langle {\widehat{C}}_{\ell }^{\text{XX}}\rangle ={\displaystyle \sum _{{\ell }^{\prime }}{}_{-2}{K}_{\ell {\ell }^{\prime }}{C}_{{\ell }^{\prime }}^{\text{XX}}}.$$(44)

For the temperature×polarization case, we can show that $$\langle {\widehat{C}}_{\ell }^{\text{TE}}\rangle ={\displaystyle \sum _{{\ell }^{\prime }}{}_{\times }{K}_{\ell {\ell }^{\prime }}{C}_{{\ell }^{\prime }}^{\text{TE}}},$$(45)

with $${}_{\times }{K}_{\ell {\ell }^{\prime }}=\left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{20}\left[f^{\text{apo}}\right].$$(46)

6.3 Relating PolSpice and MASTER in harmonic space

We can translate the relations Eq. (39) into harmonic space in temperature and polarization to obtain the PolSpice estimator as a harmonic convolution of the pseudo-power spectrum estimator, $${\widehat{C}}_{\ell }^{\text{TT}}={\displaystyle \sum _{{\ell }^{\prime }}{}_0{G}_{\ell {\ell }^{\prime }}{\tilde{C}}_{{\ell }^{\prime }}^{\text{TT}}},$$(47) $${\widehat{C}}_{\ell }^{\text{EE}}-{\widehat{C}}_{\ell }^{\text{BB}}={\displaystyle \sum _{{\ell }^{\prime }}{}_{-2}{G}_{\ell {\ell }^{\prime }}\left({\tilde{C}}_{{\ell }^{\prime }}^{\text{EE}}-{\tilde{C}}_{{\ell }^{\prime }}^{\text{BB}}\right).}$$(48) $${\widehat{C}}_{\ell }^{\text{EE}}+{\widehat{C}}_{\ell }^{\text{BB}}={\displaystyle \sum _{{\ell }^{\prime }}{}_{\text{dec}}{G}_{\ell {\ell }^{\prime }}\left({\tilde{C}}_{{\ell }^{\prime }}^{\text{EE}}+{\tilde{C}}_{{\ell }^{\prime }}^{\text{BB}}\right),}$$(49) $${\widehat{C}}_{\ell }^{\text{TE}}={\displaystyle \sum _{{\ell }^{\prime }}{}_{\times }{G}_{\ell {\ell }^{\prime }}{\tilde{C}}_{{\ell }^{\prime }}^{\text{TE}}.}$$(50)

The G kernels are constructed in the same manner as the PolSpice kernels, with the operator Ξ acting on the function g = f^apo/w according to Eq. (C.7) (or equivalently, on the associated power spectrum of g via Legendre transform). They are given by $${}_0{G}_{\ell {\ell }^{\prime }}=\left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{00}\left[g\right],$$(51) $${}_{-2}{G}_{\ell {\ell }^{\prime }}=\left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{2-2}\left[g\right],$$(52) $${}_{\times }{G}_{\ell {\ell }^{\prime }}=\left(2{\ell }^{\prime }+1\right){\text{Ξ}}_{\ell {\ell }^{\prime }}^{20}\left[g\right],$$(53) $${}_{\text{dec}}{G}_{\ell {\ell }^{\prime }}=\frac{2{\ell }^{\prime }+1}{2}{\displaystyle {\int }_{-1}^{1}g\left(\theta \right)d_{22}^{\ell }\left(\theta \right)d_{2-2}^{{\ell }^{\prime }}\left(\theta \right)\text{d\hspace{0.17em}}\cos \left(\theta \right)}.$$(54)

The first three equations above reduce to the inverse of the master kernels when the PolSpice apodization function is set to 1, that is, no apodization. The last kernel, referred to as dec, is the kernel that allows the decoupling of the EE and BB spectra. Appendix C gives more details on this point. _decG is associated with integral relations in real space, thus its harmonic expression is not straightforward. This expression requires more numerical resources to be computed because no closed relations exist for the Wigner d-matrix with different multipole indices. It can be obtained with PolSpice for all ℓ by setting ${\tilde{C}}_{{\ell }^{\prime }}^{\text{EE}}={\tilde{C}}_{{\ell }^{\prime }}^{\text{BB}}={\delta }_{\ell {\ell }^{\prime }}$ as input.

6.4 Covariance of the PolSpice estimator

Given the previous relations in Eqs. (47)–(50), the covariance of the PolSpice estimator can be written as a convolution of the covariance of the pseudo-power spectrum, with $$\begin{array}{lll}{\widehat{\text{Σ}}}_{L{L}^{\prime }}^{\text{TTTT}}{}_{}\hfill & \equiv \hfill & \mathrm{cov}\left({\widehat{C}}_{\ell }^{\text{TT}},{\widehat{C}}_{{\ell }^{\prime }}^{\text{TT}}\right)\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{L{L}^{\prime }}{}_0{G}_{lL}{\widetilde{\text{Σ}}}_{L{L}^{\prime }}^{\text{TTTT}}{}_0{G}_{l^{\prime }{L}^{\prime }}.}\hfill \end{array}$$(55)

For polarization, there is mixing between the EE and BB components in the covariance. We write the polarized EEEE PolSpice covariance after defining ${}_{\pm }G\equiv \frac{1}{2}\left({}_{\text{dec}}G{\pm }_{-2}G\right)$ as $$\begin{array}{lll}{\widehat{\text{Σ}}}^{\text{EEEE}}\hfill & =\hfill & {}_{+}G{\widetilde{\text{Σ}}}^{\text{EEEE}}{}_{+}{G}^{\top }{+}_{-}G{\widetilde{\text{Σ}}}^{\text{BBEE}}{}_{+}{G}^{\top }\hfill \\ \hfill & \hfill & {+}_{+}G{\widetilde{\text{Σ}}}^{\text{EEBB}}{}_{-}{G}^{\top }{+}_{-}G{\widetilde{\text{Σ}}}^{\text{BBBB}}{}_{-}{G}^{\top }.\hfill \end{array}$$(56)

The polarized PolSpice covariance is built on the polarized pseudo-covariance, mixing the components EE and BB, thanks to the kernel _±G; see Fig. 13. This figure displays a row of the kernels that were computed on the mask SPT-3G used in our analysis. It shows the window functions that are applied to the pseudo-power spectra to produce the PolSpice spectra. They correct for the bias due to the mask, but introduce a small bias due to the lack of information at large scales. The temperature kernel ₀G and the polarization ₊G kernel are almost identical. The leakage kernel ₋G (all negative), which accounts for the mixing of the E and B polarization pseudo-spectra in the PolSpice spectrum, is orders of magnitudes smaller than the other two. Hence, the BB covariance terms do not affect the EEEE covariance. On the other hand, the EE terms affect the BBBB covariance because the EE spectrum is a few orders of magnitude larger than BB. The PolSpice apodizing function of the correlation function we used is $$f^{\text{apo}}\left(\theta \right)=\{\begin{array}{ll}\frac{1}{2}\left(1+\cos \frac{\pi \theta }{{\theta }_{\max }}\right)\hfill & \forall \theta <{\theta }_{\max }\hfill \\ 0\hfill & \text{otherwise.}\hfill \end{array}$$(57)

Here we have, θ_max = π/6. Without apodization, but with partial sky, as for Planck, the decoupling kernel _decG is not null, still resulting in a nonzero ₋G kernel. However, it is orders of magnitude smaller than ₊G, hence we can compute EEEE ignoring the leakage from covariance terms that include BB.

Fig. 13 Amplitude of the PolSpice convolution kernels G for the SPT-3G footprint. The negative terms are plotted with thinner lines.

6.5 Accuracy of the covariance approximations for the PolSpice estimator

We can build estimates of the PolSpice spectrum covariance by convolving the pseudo-spectrum covariance with the appropriate kernels following Eqs. (55) and (56). To calculate the pseudo-spectrum covariance, we can use the NKA, INKA, FRI, and ACC approximations or the exact computation. Figure 14 shows the accuracy of the binned PolSpice covariance calculated with the approximations compared to the exact calculation. The results are similar to those we found for the accuracy of the pseudo-spectrum covariances shown in Fig. 6. The NKA, INKA, and FRI approaches provide a good estimate of the PolSpice covariance. However, the ACC approach improves the existing approximations dramatically. This shows that our results for the accuracy of the pseudo-covariance also hold for the PolSpice.

7 Summary and conclusions

One of the key ingredients of cosmological analysis based on power spectra are covariance matrices. Accurate covariance matrices ensure precise error bars and an unbiased estimation of cosmological parameters. The analytical estimation of these matrices can be difficult when small sky fractions are observed because existing approximations might fail. We have considered the specific example of estimating accurate analytical signal-signal covariance matrices for the SPT-3G CMB experiment, whose survey covers about 4%, without masking the contribution of point sources. We considered the cases of estimating the matrix for pseudo-power spectrum and for the PolSpice power spectrum estimator.

First, we implemented in Sect. 3 an expensive exact calculation of the covariance of the pseudo-power spectrum in intensity and polarization for the first time. We used a map-based algorithm that is accelerated thanks to the HEALPix pixelation tools. We were thus able to compute the entire covariance matrix up to ℓ_max,ex = 1000 exactly. We also obtained a selection of rows of the covariance of particular interest up to ℓ = 2000.

Based on this result, we were able to estimate the accuracy of the existing approximations in Sect. 4 precisely by comparing them to the binned exact covariances of the pseudo-power spectra measured on the SPT-3G patch. The approximations were found to be precise to the 5% level.

Then, using the code we developed for an exact computation of the covariance matrix, we estimated the covariance-coupling kernel θ, which determines how the CMB power spectrum couples into the covariance matrix. We were able to understand why the existing approximations in the literature fail to achieve a precision better than 5%. We then proposed a new approximation in Sect. 5, the ACC, which is more computationally expensive than the existing approximations, but allows a more precise estimation of the covariance matrix at the 1% level.

Finally, we showed in Sect. 6 that we are able to build the covariance of the PolSpice power spectrum in both temperature and polarization using a harmonic correction. This computation is exact and based on the PolSpice algorithm real-space corrections that we translated into harmonic space. Through this correction, we produced estimates of the PolSpice covariance matrix based on the previous approximations of the pseudo-power spectrum covariance. The accuracy of the resulting PolSpice covariance approximations is the same as for the pseudo-power spectrum.

While this paper considered the particular example of the SPT-3G experiment, the results can be extended to a non-CMB power spectrum analysis such as that of weak-lensing shear or photometric catalogs. We would also like to stress that the accuracy of any of the approximations presented in this paper (existing or new) is reduced when a point-source mask is included in the sky footprint. Nevertheless, the exact computation of the covariance matrices still holds in this particular case. While previous experiments have included the effect of point-source masks through the use of simulations (see, e.g., Planck Collaboration XI 2016) or by inpainting the holes with constrained realizations (see, e.g., Benoit-Lévy et al. 2013), additional work is required to find an analytical calculation of this contribution.

Fig. 14 Relative differences of binned PolSpice covariance matrices calculated using approximations of the pseudo-spectrum covariance with respect to the exact computation: ${\widehat{\text{Σ}}}_{b{b}^{\prime }}^{\text{APP}}/{\widehat{\text{Σ}}}_{b{b}^{\prime }}-1$, for TTTT (left-hand side) and EEEE (right-hand side), with binning Δℓ = 50. On the first row we plot the relative differences for the diagonal, i.e. b = b′, while on the second row we plot the differences for the first off-diagonal, i.e. b′ = b + 1. Similarly to the case of pseudo-covariances in Fig. 6, we find acceptable accuracy for the NKA (dash light blue), INKA (dash-dot dark blue) and FRI (dash-double-dot purple) approximations, while our ACC approximation (solid red) improves overall the others. The PolSpice covariance matrices have been calculated using Eqs. (55) and (56). The last row displays the corresponding binned underlying renormalized spectrum TT or EE, to highlight the fact that the differences in the covariances are on the peaks and in the troughs of the spectra, i.e. where the curvature of the spectrum is maximal.

Acknowledgements

We are grateful to the SPT collaboration for discussions and suggestions. We thank Lennart Balkenhol for discussions and insightful comments on the manuscript. This work has received funding from the French Centre National d’Etudes Spatiales (CNES). This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 101001897).

Appendix A Analysis of the curved sky

This appendix describes the mathematical tools that are used on curved sky for a CMB analysis. We make use of the spherical harmonic decomposition of Gaussian fields. We introduce various operators that allow us to express the couplings and the covariance of the power spectra. We make use of some geometrical relations of spherical harmonics to obtain our results. In this appendix, we introduce formulae that can either be used in the temperature or in the polarization case.

Appendix A.1 Temperature

We first consider the case of a map of the CMB intensity anisotropies $T\left(\widehat{n}\right)$ observed in direction $\widehat{n}$. The anisotropies are distributed as a Gaussian random field with a corresponding power spectrum $C_{\ell }^{\text{TT}}$, observed through a mask $W\left(\widehat{n}\right)$.

Harmonic coefficients and underlying power spectrum

The intensity map can be decomposed with spin-0 spherical harmonics to obtain the harmonic coefficients and their variance, the intensity power spectrum, which fully characterizes the physical properties of the field, $$a_{\ell m}^{\text{T}}={\displaystyle \int \text{d}\widehat{u}}T{\left(\widehat{u}\right)}_0{Y}_{\ell m}^{*}\left(\widehat{u}\right),$$(A.1) $$\langle a_{\ell m}^{\text{T}}{a}_{{\ell }^{\prime }{m}^{\prime }}^{\text{T*}}\rangle =C_{\ell }^{\text{TT}}{\delta }_{\ell {\ell }^{\prime }}{\delta }_{m{m}^{\prime }}.$$(A.2)

Here the brackets 〈〉 indicate an average over many realizations of the maps.

Weighted mask $W\left(\widehat{n}\right)$

The weighted mask is a real map with weights from 0 to 1 that is used to taper the data on the border of the survey area in order to reduce Fourier ringing when harmonic decomposition is used. We define the mask harmonic coefficients and its power spectrum as $$w_{\ell m}={\displaystyle \int \text{d}\widehat{u}}W{\left(\widehat{u}\right)}_0{Y}_{\ell m}^{*}\left(\widehat{u}\right),$$(A.3) $$W_{\ell }\equiv \frac{1}{2\ell +1}{\displaystyle \sum _m{w}_{\ell m}{w}_{\ell m}^{*}.}$$(A.4)

Pseudo-power spectrum estimator

One way to obtain a biased estimator of the power spectrum in CMB experiments is to define the pseudo-harmonic coefficients and the pseudo-power spectrum, $${\tilde{a}}_{\ell m}^{\text{T}}\equiv {\displaystyle \int \text{d}\widehat{u}W\left(\widehat{u}\right)T{\left(\widehat{u}\right)}_0{Y}_{\ell m}^{*}\left(\widehat{u}\right)},$$(A.5) $${\tilde{C}}_{\ell }^{\text{TT}}\equiv \frac{1}{2\ell +1}{\displaystyle \sum _{m=-\ell }^{\ell }{\left|{\tilde{a}}_{\ell m}^{\text{T}}\right|}^2.}$$(A.6)

Relation between harmonic coefficients

We relate the masked pseudo-harmonic coefficients to the unmasked coefficient with $${\tilde{a}}_{\ell m}^{\text{T}}={\displaystyle \sum _{{\ell }^{\prime }{m}^{\prime }}{a}_{{\ell }^{\prime }{m}^{\prime }}^{\text{T}}{}_0{I}_{\ell m{\ell }^{\prime }{m}^{\prime }}\left[W\right].}$$(A.7)

The I[W] couplings are defined below and can be expressed in terms of sums over Wigner-3j symbols and the w_ℓm, with $$\begin{array}{lll}{}_s{I}_{\ell m{\ell }^{\prime }{m}^{\prime }}\left[W\right]\hfill & =\hfill & {\displaystyle \int \text{d}{\widehat{u}}_s{Y}_{\ell m}\left(\widehat{u}\right)W{\left(\widehat{u}\right)}_s{Y}_{{\ell }^{\prime }{m}^{\prime }}^{*}\left(\widehat{u}\right),}\hfill \\ \hfill & =\hfill & {\displaystyle \sum _{LM}{w}_{LM}{\left(-1\right)}^{m^{\prime }}{\left[\frac{\left(2\ell +1\right)\left(2{\ell }^{\prime }+1\right)\left(2L+1\right)}{4\pi }\right]}^{1/2}}\hfill \end{array}$$(A.8) $$\left(\begin{array}{lll}\ell \hfill & {\ell }^{\prime }\hfill & L\hfill \\ -s\hfill & s\hfill & 0\hfill \end{array}\right)\times \left(\begin{array}{lll}\ell \hfill & {\ell }^{\prime }\hfill & L\hfill \\ m\hfill & -m\hfill & M\hfill \end{array}\right).$$(A.9)

Here, we anticipated the extension of this notation to the polarized case, which deals with spin-2 fields. The mask-dependent _sI_ℓmℓ′mv [W] coupling coefficients relate the underlying harmonic coefficients to the measured pseudo-harmonic coefficients. In the full-sky case, we obtain _sI_{ℓmℓ′m′} [1] = δ_ℓℓ′δ_mm′ through the orthonormality properties of the spin-weighted spherical harmonics.

We introduce some useful relations that are demonstrated in Hivon et al. (2002), $${\displaystyle \sum _{\ell m}{}_s{I}_{{\ell }_1{m}_1{\ell }_m}{\left[W\right]}_s{I}_{{\ell }_2{m}_2\ell m}^{*}\left[W\right]{=}_s{I}_{{\ell }_1{m}_1{\ell }_2{m}_2}\left[W^2\right]},$$(A.10) $${\displaystyle \sum _{m_1{m}_2}\frac{{}_s{I}_{{\ell }_1{m}_1{\ell }_2{m}_2}{\left[W\right]}_{s^{\prime }}{I}_{{\ell }_1{m}_1{\ell }_2{m}_2}^{*}\left[W\right]}{\left(2{\ell }_1+1\right)\left(2{\ell }_2+1\right)}={\text{Ξ}}_{{\ell }_1{\ell }_2}^{s{s}^{\prime }}\left[W\right].}$$(A.11)

Here we introduced the symmetric operator ${\text{Ξ}}^{s{s}^{\prime }}$ acting on a power spectrum A_ℓ, $${\text{Ξ}}_{\ell {\ell }^{\prime }}^{s{s}^{\prime }}\left|A\right|\equiv {\displaystyle \sum _L\frac{2L+1}{4\pi }{A}_L\left(\begin{array}{lll}\ell \hfill & {\ell }^{\prime }\hfill & L\hfill \\ s\hfill & -s\hfill & 0\hfill \end{array}\right)\left(\begin{array}{lll}\ell \hfill & {\ell }^{\prime }\hfill & L\hfill \\ s^{\prime }\hfill & -s^{\prime }\hfill & 0\hfill \end{array}\right).}$$(A.12)

We extend this definition to an operator acting on a map $A\left(\widehat{n}\right)$, with $${\text{Ξ}}_{\ell {\ell }^{\prime }}^{s{s}^{\prime }}\left|A\right|\equiv {\text{Ξ}}_{\ell {\ell }^{\prime }}^{s{s}^{\prime }}\left|A\right|,$$(A.13)