CMB internal delensing with general optimal estimator for higher-order correlations

TL;DR

This paper introduces a new realization-dependent estimator for internal delensing of CMB B modes, effectively reducing biases from correlations and improving sensitivity to primordial signals.

Contribution

It develops a general optimal estimator for higher-order correlations that corrects biases in CMB internal delensing without significantly reducing delensing efficiency.

Findings

The RD estimator effectively corrects biases in CMB internal delensing.

The method reduces sensitivity to simulation uncertainties.

Delensing efficiency is maintained with bias correction.

Abstract

We present a new method for delensing $B$ modes of the cosmic microwave background (CMB) using a lensing potential reconstructed from the same realization of the CMB polarization (CMB internal delensing). The $B$-mode delensing is required to improve sensitivity to primary $B$ modes generated by, e.g., the inflationary gravitational waves, axion-like particles, modified gravity, primordial magnetic fields, and topological defects such as cosmic strings. However, the CMB internal delensing suffers from substantial biases due to correlations between observed CMB maps to be delensed and that used for reconstructing a lensing potential. We construct a realization-dependent (RD) estimator for correcting these biases by deriving a general optimal estimator for higher-order correlations. The RD estimator is less sensitive to simulation uncertainties. Compared to the previous methods, we find that the RD estimator corrects the biases without substantial degradation of the delensing efficiency.

Equations

$$\tilde{X}(\hat{\mathbf{n}})=X(\hat{\mathbf{n}}+\nabla \varphi (\hat{\mathbf{n}})).$$

$${\Theta }_{\mathbf{\ell }}$$

$$E_{\mathbf{\ell }}\pm \mathrm{i}{B}_{\mathbf{\ell }}$$

$${\tilde{\Theta }}_{\mathbf{\ell }}$$

$${\tilde{E}}_{\mathbf{\ell }}$$

$${\tilde{B}}_{\mathbf{\ell }}$$

$$z_{\mathbf{\ell }\mathbf{L}}$$

$$w_{\mathbf{\ell }\mathbf{L}}$$

$$⟨{\tilde{E}}_{\mathbf{\ell }}{\tilde{B}}_{\mathbf{L}-\mathbf{\ell }}{⟩}_{\mathrm{CMB}}=-{\tilde{C}}_{\ell }^{\mathrm{EE}}{w}_{\mathbf{L}-\mathbf{\ell },\mathbf{L}}{\varphi }_{\mathbf{L}},$$

$${\hat{\varphi }}_{\mathbf{L}}^{EB}=A_L\int \frac{{\mathrm{d}}^2\mathbf{\ell }}{(2\pi {)}^2}\frac{-{\tilde{C}}_{\ell }^{\mathrm{EE}}{w}_{\mathbf{L}-\mathbf{\ell },\mathbf{L}}}{{\hat{C}}_{\ell }^{\mathrm{EE}}{\hat{C}}_{|\mathbf{L}-\mathbf{\ell }|}^{\mathrm{BB}}}{\hat{E}}_{\mathbf{\ell }}{\hat{B}}_{\mathbf{L}-\mathbf{\ell }}.$$

$$(\hat{E}\star \hat{\varphi }{)}_{\mathbf{\ell }}=-\int \frac{{\mathrm{d}}^2{\mathbf{\ell }}^{\prime }}{(2\pi {)}^2}\frac{C_{{\ell }^{\prime }}^{\mathrm{EE}}}{{\hat{C}}_{{\ell }^{\prime }}^{\mathrm{EE}}}\frac{C_{|\mathbf{\ell }-{\mathbf{\ell }}^{\prime }|}^{\varphi \varphi }}{{\hat{C}}_{|\mathbf{\ell }-{\mathbf{\ell }}^{\prime }|}^{\varphi \varphi }}{w}_{\mathbf{\ell },\mathbf{\ell }-{\mathbf{\ell }}^{\prime }}{\hat{E}}_{{\mathbf{\ell }}^{\prime }}{\hat{\varphi }}_{\mathbf{\ell }-{\mathbf{\ell }}^{\prime }}.$$

$$C_{\ell }^{{\hat{B}}^{\mathrm{d}}{\hat{B}}^{\mathrm{d}}}=\int \frac{{\mathrm{d}}^2{\mathbf{\ell }}^{\prime }}{(2\pi {)}^2}{w}_{\mathbf{\ell },\mathbf{\ell }-{\mathbf{\ell }}^{\prime }}^2\left(1-\frac{(C_{{\ell }^{\prime }}^{\mathrm{EE}}{)}^2}{{\hat{C}}_{{\ell }^{\prime }}^{\mathrm{EE}}}\frac{(C_{|\mathbf{\ell }-{\mathbf{\ell }}^{\prime }|}^{\varphi \varphi }{)}^2}{{\hat{C}}_{|\mathbf{\ell }-{\mathbf{\ell }}^{\prime }|}^{\varphi \varphi }}\right).$$

$$⟨|{\hat{B}}^{\mathrm{d}}{|}^2⟩=⟨|\hat{B}-(\hat{E}\star {\hat{\varphi }}^{\prime }){|}^2⟩=(2\pi {)}^2{\delta }_0^{\mathrm{D}}{C}^{{\hat{B}}^{\mathrm{d}}{\hat{B}}^{\mathrm{d}}}.$$

$$N^{\mathrm{bias}}$$

$$=-2⟨\hat{B}(\hat{E}\star (\hat{\varphi }-{\hat{\varphi }}^{\prime }){)}^{\ast }⟩$$

$$+⟨|(\hat{E}\star \hat{\varphi }){|}^2-|(\hat{E}\star {\hat{\varphi }}^{\prime }){|}^2⟩$$

$$=N^{(4)}+N^{(6)},$$

$$N^{(4)}$$

$$N^{(6)}$$

$$N^{(4)}=-2⟨{\hat{B}}^{\ast }(\hat{E}\star \hat{\varphi }){⟩}_{\mathrm{d}}-2⟨{\hat{B}}^{\ast }(\hat{E}\star (\hat{\varphi }-\varphi )){⟩}_{\mathrm{c}},$$

$$(\hat{E}\star \hat{\varphi }{)}_{\mathbf{\ell }}$$

$$W_{\mathbf{\ell }{\mathbf{\ell }}^{\prime }{\mathbf{\ell }}^{\prime \prime }}\equiv \frac{C_{{\ell }^{\prime }}^{\mathrm{EE}}{C}_{|\mathbf{\ell }-{\mathbf{\ell }}^{\prime }|}^{\varphi \varphi }{w}_{\mathbf{\ell },\mathbf{\ell }-{\mathbf{\ell }}^{\prime }}}{{\hat{C}}_{{\ell }^{\prime }}^{\mathrm{EE}}{\hat{C}}_{|\mathbf{\ell }-{\mathbf{\ell }}^{\prime }|}^{\varphi \varphi }}\frac{A_{|\mathbf{\ell }-{\mathbf{\ell }}^{\prime }|}{\tilde{C}}_{{\ell }^{\prime \prime }}^{\mathrm{EE}}{w}_{\mathbf{\ell }-{\mathbf{\ell }}^{\prime }-{\mathbf{\ell }}^{\prime \prime },\mathbf{\ell }-{\mathbf{\ell }}^{\prime }}}{{\hat{C}}_{{\ell }^{\prime \prime }}^{\mathrm{EE}}{\hat{C}}_{|\mathbf{\ell }-{\mathbf{\ell }}^{\prime }-{\mathbf{\ell }}^{\prime \prime }|}^{\mathrm{BB}}}.$$

$$⟨|(\hat{E}\star \hat{\varphi }){|}_{\mathbf{\ell }}^2⟩=\int \frac{{\mathrm{d}}^2{\mathbf{\ell }}_1^{\prime }}{(2\pi {)}^2}\int \frac{{\mathrm{d}}^2{\mathbf{\ell }}_1^{\prime \prime }}{(2\pi {)}^2}\int \frac{{\mathrm{d}}^2{\mathbf{\ell }}_2^{\prime }}{(2\pi {)}^2}\int \frac{{\mathrm{d}}^2{\mathbf{\ell }}_2^{\prime \prime }}{(2\pi {)}^2}$$

$$\times W_{\mathbf{\ell }{\mathbf{\ell }}_1^{\prime }{\mathbf{\ell }}_1^{\prime \prime }}{W}_{\mathbf{\ell }{\mathbf{\ell }}_2^{\prime }{\mathbf{\ell }}_2^{\prime \prime }}⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_1^{\prime \prime }}{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩.$$

$$⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_1^{\prime \prime }}{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{⟩}_{\mathrm{d}}$$

$$=2⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_1^{\prime \prime }}⟩⟨{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{⟩}_{\mathrm{c}}$$

$$+2⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩⟨{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{⟩}_{\mathrm{c}}$$

$$+⟨{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}^{\ast }{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{⟩}_{\mathrm{c}}$$

$$+⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩⟨{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{⟩}_{\mathrm{c}}$$

$$+⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_1^{\prime \prime }}⟩⟨{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩⟨{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩$$

$$+⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }⟩⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩⟨{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩,$$

$$⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime }}^{\ast }⟩⟨{\hat{E}}_{{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩⟨{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{B}}_{\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }⟩$$

$$=⟨{\hat{E}}_{1,{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{1,{\mathbf{\ell }}_2^{\prime }}^{\ast }{⟩}_1⟨{\hat{E}}_{2,{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{2,{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{⟩}_2⟨{\hat{B}}_{3,\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{B}}_{3,\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{⟩}_3$$

$$=⟨{\hat{E}}_{1,{\mathbf{\ell }}_1^{\prime }}{\hat{E}}_{1,{\mathbf{\ell }}_2^{\prime }}^{\ast }{\hat{E}}_{2,{\mathbf{\ell }}_1^{\prime \prime }}{\hat{E}}_{2,{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{\hat{B}}_{3,\mathbf{\ell }-{\mathbf{\ell }}_1^{\prime }-{\mathbf{\ell }}_1^{\prime \prime }}{\hat{B}}_{3,\mathbf{\ell }-{\mathbf{\ell }}_2^{\prime }-{\mathbf{\ell }}_2^{\prime \prime }}^{\ast }{⟩}_{1,2,3},$$

$$⟨|(\hat{E}\star \hat{\varphi }){|}^2⟩\supset ⟨(E_1\star {\hat{\varphi }}^{E_2{B}_3})(E_2\star {\hat{\varphi }}^{E_1{B}_3}{)}^{\ast }{⟩}_{1,2,3}.$$

$$N^{(6)}$$

$$+2⟨⟨E_1\star {\hat{\varphi }}^{E_2{B}_2}(E_2\star {\hat{\varphi }}^{E_1{B}_2}{)}^{\ast }{⟩}_1{⟩}_{2,\mathrm{c}}$$

$$+⟨⟨E_2\star {\hat{\varphi }}^{E_2{B}_1}(E_2\star {\hat{\varphi }}^{E_2{B}_1}{)}^{\ast }{⟩}_1{⟩}_{2,\mathrm{c}}$$

$$+⟨⟨E_2\star {\hat{\varphi }}^{E_1{B}_2}(E_2\star {\hat{\varphi }}^{E_1{B}_2}{)}^{\ast }{⟩}_1{⟩}_{2,\mathrm{c}}$$

$$+⟨(E_1\star {\hat{\varphi }}^{E_1{B}_3})(E_2\star {\hat{\varphi }}^{E_2{B}_3}{)}^{\ast }{⟩}_{1,2,3}$$

$$+⟨(E_1\star {\hat{\varphi }}^{E_2{B}_3})(E_2\star {\hat{\varphi }}^{E_1{B}_3}{)}^{\ast }{⟩}_{1,2,3}$$

Full text

I Introduction

The cosmic microwave background (CMB) anisotropies are one of the best probe in cosmology. $B$ modes of the CMB polarization could be generated by the inflationary gravitational waves, axion-like particles, primordial magnetic fields, modified gravity, and topological defects such as cosmic strings (e.g., Refs. Polnarev (1985); Lue et al. (1999); Kamionkowski et al. (1997); Seljak and Zaldarriaga (1997); Seljak and Slosar (2006); Shaw and Lewis (2010); Pogosian et al. (2011); Teng et al. (2011). In addition, the gravitational lensing converts some of the $E$ modes into $B$ modes Zaldarriaga and Seljak (1998), and this effect has been measured by multiple CMB experiments (e.g., Refs. SPTpol Collaboration (2013) (Hanson, D. et al.); POLARBEAR Collaboration (2014); ACT Collaboration (2016) (Sherwin, Blake D. et. al.); Planck Collaboration (2016a); Bicep2 / Keck Array Collaboration (2016)). The recent precise measurement of the CMB polarization by the BICEP2/Keck Array experiment indicates that the observed $B$ modes are dominated by the gravitational lensing $B$ modes Bicep2 / Keck Array Collaboration (2016). Thus, removal of the lensing-induced $B$ modes, usually referred to as delensing, will be soon required to improve sensitivity to the primary $B$ modes Kesden et al. (2002); Knox and Song (2002). The temperature and $E$-mode delensing are also important to constrain, e.g., the effective number of the relativistic species Green et al. (2016).

Several works have discussed methods to remove the lensing contributions in the CMB anisotropies. Refs. Seljak and Hirata (2004); Smith et al. (2012) provide a method to construct a template of the lensed CMB fluctuations and to remove it from the observed CMB fluctuations. LABEL:Green:2016 proposes a delensing technique to remove the lensing-induced CMB anisotropies including the higher-order terms of the lensing potential. Refs. Simard et al. (2015); Sherwin and Schmittfull (2015) show the delensing method using the cosmic infrared background (CIB) (and this technique is demonstrated by Refs. Larsen et al. (2016); SPTpol Collaboration (2017) (Manzotti, A. et al.)).

The delensing methods, however, suffer from biases if the CMB maps are delensed with a lensing potential derived from the CMB maps (CMB internal delensing). In the CMB internal delensing, we use the same realization of the observed CMB anisotropies for the delensing and reconstruction of the lensing potential. Therefore, the CMB anisotropies to be delensed correlates with that used in the lensing potential reconstruction. As we discuss latter in Sec. II, this correlation substantially biases the power and scatter of the delensed CMB spectrum (delensing bias).

Methods to correct the delensing bias have been explored in several works. The delensing bias is pointed out by LABEL:Teng:2011xc and is characterized by LABEL:Namikawa:2014a (see also LABEL:Sehgal:2016 for the CMB temperature and LABEL:Carron:2017 for the Planck cases). Refs. Teng et al. (2011); Sehgal et al. (2016) proposed a method to mitigate the delensing bias by splitting CMB anisotropies into annuli in multipole space so that CMB multipoles to be delensed and that used for the lensing reconstruction are uncorrelated. In this method, however, the efficiency of the delensing is degraded because we lose the CMB multipoles in the lensing reconstruction and precision of the reconstructed lensing potential decreases.

In this paper, we study an alternative way of correcting the delensing bias to avoid degradation of the delensing efficiency and to estimate the delensed spectrum more accurately than using the simulation alone. To construct an estimator for the delensing bias, we use the similar analogy of the “realization-dependent” (RD) method in the lensing reconstruction Namikawa et al. (2013); Planck Collaboration (2014, 2015). In the reconstruction of the lensing potential, the leading-order bias in estimating the lensing-potential power spectrum comes from the disconnected part of the observed CMB trispectrum. In the RD method, the bias is estimated by combining both simulated and observed CMB fluctuations so that the estimator is optimal and less sensitive to the error in the simulated CMB covariance. The RD method in the lensing reconstruction is motivated by the optimal trispectrum estimator Regan et al. (2010).

Compared to the lensing reconstruction, the delensing bias contains higher order correlations such as six point correlations. Thus, we first generalize the optimal trispectrum estimator to the case with high-order correlations (polyspectra). Then we employ the optimal polyspectra estimator to construct RD estimators for the delensing bias. We demonstrate that the RD method removes the delensing bias without degradation of the delensing efficiency. In this paper, for simplicity, we derive the RD method in the flat-sky approximation but the extension to the full sky approach is straightforward.

This paper is organized as follows. Sec. II reviews the $B$-mode delensing and the delensing bias. Sec. IV derives an optimal estimator for polyspectra based on the Edgeworth expansion of the CMB likelihood. Sec. V provides the RD estimators for the delensing bias and shows results of numerical simulation. Sec. VI is devoted to summary and discussion.

II CMB internal delensing

Here we briefly summarize notations in this paper and methods for the lensing reconstruction and delensing.

II.1 CMB polarization

The lensing effect on the CMB anisotropies is expressed by a remapping of the primary CMB anisotropies. Denoting the primary CMB fluctuations at position $\hat{\bm{n}}$ on the last scattering surface as $X(\hat{\bm{n}})$, where $X$ is either the CMB temperature $\Theta$, or polarization anisotropies $Q\pm{\rm i}\hskip 0.50003ptU$, the lensed CMB anisotropies in a direction $\hat{\bm{n}}$, are given by (e.g., LABEL:Lewis:2006fu)

[TABLE]

Here the two-dimensional vector, $\bm{\nabla}\phi(\hat{\bm{n}})$, is the deflection angle. The Fourier modes of the temperature anisotropies, and the $E$ and $B$ modes are defined as

[TABLE]

where $\varphi_{\bm{\ell}}$ is the angle of $\bm{\ell}$ measured from the $x$-axis. Substituting Eq. (1) into the above equation gives the following expression of the lensed CMB anisotropies up to first order of the lensing potential Hu and Okamoto (2002):

[TABLE]

where the weight functions are given by

[TABLE]

If the primordial $B$ modes are ignored, the $B$ modes induced by gravitational lensing are expressed as a convolution of the lensing potential and primary $E$ modes.

II.2 Lensing reconstruction

Precision of the lensing reconstruction in the future will be almost determined by the $EB$ quadratic estimator, and we focus on the case using the $EB$ estimator. For a fixed realization of the lensing potential in the universe, the lensing effect on CMB maps produces off-diagonal elements of the CMB covariance. From Eq. (6), the off-diagonal elements become Hu and Okamoto (2002)

[TABLE]

where $\langle\cdots\rangle_{\rm CMB}$ denotes the ensemble average over unlensed CMB anisotropies with a fixed realization of the lensing potential $\phi$. We ignore the higher-order terms of the lensing potential and contributions from primary $B$ modes 111 Note that the presence of the primordial gravitational-wave (GW) B modes does not affect the precision of the lensing reconstruction and delensing efficiency. This is because the GW B modes are significant only at the large scale ($\ell<100$) and the lensing potential is reconstructed from small scale B modes. . $\widetilde{C}^{\rm EE}_{\ell}$ is the lensed $E$-mode spectrum to mitigate the higher-order terms of the lensing-potential power spectrum (Hanson et al., 2011; Lewis et al., 2011). The lensing estimators are then described as Hu and Okamoto (2002)

[TABLE]

Here $\widehat{E}$ ($\widehat{B}$) is the observed $E$ ($B$) modes, $\widehat{C}^{\rm EE}$ ($\widehat{C}^{\rm BB}$) is the observed $E$ ($B$) mode spectrum, and $A_{L}$ is the quadratic estimator normalization given by LABEL:Hu:2001kj. In practical analysis, the filtering functions, $1/\widehat{C}^{\rm EE}_{\ell}$ and $1/\widehat{C}^{\rm BB}_{\ell}$, are generalized to the inverse of the covariance matrix of $E$ and $B$ modes, and are determined by simulation (e.g., see LABEL:Story:2014hni). Unless otherwise stated, $\widehat{\phi}$ is equal to $\widehat{\phi}^{EB}$ with the correction of the mean-field bias, $\langle\widehat{\phi}\rangle$.

II.3 Delensing

A method to remove lensing $B$ modes from observed $B$ modes is to make a template of the lensing $B$ modes with a measured lensing potential and $E$ modes. Hereafter, we refer to this technique as the template method. The template of the lensing $B$ modes up to first order of the lensing potential is defined as Smith et al. (2012)

[TABLE]

Here the optimal weight is given as a product of the Wiener filters of the $E$ modes and lensing potential. In practical analysis, the denominator of the Wiener filter is determined by simulations. This expression is derived by assuming that the reconstructed lensing potential is uncorrelated with the CMB anisotropies Smith et al. (2012). While it would be possible to make a template including this correlation, Eq. (11) is numerically tractable and useful for delensing analysis in practical cases. A lensing $B$-mode template to remove the higher-order terms of $\phi$ is also possible to construct but Eq. (11) is enough to reproduce most of the lensing $B$ modes Smith et al. (2012). The delensed $B$-mode power spectrum is then given by Smith et al. (2012)

[TABLE]

Alternatively, in high sensitivity polarization experiments, the delensed $B$ modes are obtained by the inverse remapping of the observed CMB map using the Wiener-filtered lensing potential Green et al. (2016). Hereafter, we call this approach the remapping method. While the remapping method is a simple way to remove higher order terms of the lensing potential, the inverse-remapping procedure simultaneously remaps noise fluctuations which should be taken into account for a noisy experiment Carron et al. (2017).

III Delensing bias

To delens $B$ modes, the observed $B$ modes to be delensed and that used in the lensing reconstruction are correlated and the delensed $B$ modes have significant biases. Here we discuss the bias in the delensed $B$-mode power spectrum.

III.1 CMB internal delensing with simulation

To study the CMB internal delensing, we simulate the lensed CMB, noise and reconstructed lensing potential maps as follows. We first generate $1200$ realizations of unlensed CMB and lensing potential maps in $10$deg$\times 10$deg region as random Gaussian fields using CAMB output spectra Lewis et al. (2000) 222 The first $300$ realizations are used as “observed” CMB fluctuations while the other $3\times 300$ realizations are used for estimating the (RD) delensing biases. . We assume the cosmological parameters consistent with LABEL:P15:main without the inflationary B modes. The lensed CMB maps are obtained by remapping the unlensed CMB maps with the lensing potential Louis et al. (2013). The instrumental noise spectrum is computed with the formula of LABEL:Knox:1999 assuming $\sqrt{2}\times 3\,\mu\text{K-arcmin}$ white noise in polarization map and $1$ arcmin Gaussian beam, and the noise maps are generated with this noise spectrum. The lensing potential is reconstructed from the lensed+noise (observed) maps with the $EB$ quadratic estimator using CMB multipoles in the range between $200$ and $3000$. We do not use the multipoles at $\ell<200$ since the large-scale B modes are dominated by the Galactic foregrounds and are also not important in the lensing potential reconstruction. To characterize the delensing bias, we also use the input lensing potential with a random Gaussian noise whose power corresponds to the reconstruction noise 333 The leading order reconstruction noise is expressed as the disconnected bias, but the higher-order biases are not negligible in the following discussion. Thus we denote the reconstruction noise spectrum as $C_{L}^{\widehat{\phi}\widehat{\phi}}-C_{L}^{\phi\phi}$. .

We then perform delensing with the above lensing potentials and observed CMB maps. In the template method, the lensing $B$ modes are estimated from Eq. (11) and are subtracted from the observed $B$ modes. In the remapping method, the observed $Q$/$U$ maps are remapped by the Wiener-filtered lensing potential and we obtain the delensed $B$ modes from the delensed $Q$/$U$ maps. We do not apply Wiener filter for the CMB polarizations since its impact is negligible in our case.

III.2 Comparison between template and remapping methods

Fig. 1 shows delensed $B$ modes with the template and remapping methods. For simplicity, we use the input lensing potential with the random Gaussian reconstruction noise since it does not create the delensing bias. The delensed $B$-mode power spectra using the template and remapping methods are similar since the lensing $B$ modes are approximately expressed as Eq. (6) Challinor and Lewis (2005). In the following sections, we focus on the method of the (linear order) template method.

III.3 Delensing bias in $B$ mode power spectrum

Next we use the reconstructed lensing potential for delensing. Fig. 2 shows the delensed $B$-mode spectrum using the input lensing potential with the random Gaussian reconstruction noise or the reconstructed lensing potential. 444 The instrumental noise power spectrum is subtracted from the delensed $B$-mode spectrum and thus the debiased spectra can be negative.

Since the delensed $B$-mode spectrum with the random Gaussian reconstruction noise does not contain the delensing bias, the difference of the $B$-mode spectra between the two cases describes the delensing bias. The discrepancy appears at $L\gtrsim 200$ and the scatter of the power spectrum is also reduced.

To understand the delensing bias, we expand the delensed $B$-mode spectrum into the four and six point correlations. If the lensing potential is estimated from an external data or the reconstruction noise is uncorrelated with the CMB anisotropies, the delensed $B$-mode power spectrum is simply equivalent to the power spectrum given in Eq. (12):

[TABLE]

Here $\widehat{\phi}^{\prime}$ is given as a sum of the true lensing potential and the random Gaussian reconstruction noise, $\delta^{\rm D}_{\bm{0}}$ is the Delta function, and we omit the multipole $\bm{\ell}$. In the CMB internal delensing, the delensed spectrum has additional contributions. The additional terms are given as

[TABLE]

where we define

[TABLE]

The first term, $N^{(4)}$, is a four point correlation and leads to the following bias terms,

[TABLE]

where $\langle\cdots\rangle_{\rm d}$ and $\langle\cdots\rangle_{\rm c}$ are the disconnected and connected parts, respectively.

The second term, $N^{(6)}$, is a six point correlation and is decomposed into a connected term, products of the power spectrum and trispectrum, and products of three power spectra. To see this, from Eqs. 11 and 10, we explicitly rewrite the lensing $B$-mode template as

[TABLE]

where we define

[TABLE]

The power spectrum of the above quantity is the six point correlation, and we obtain

[TABLE]

The disconnected part of the six point correlation is expanded as

[TABLE]

where we omit the terms irrelevant to the delensing bias (i.e., the terms contained in $|(\widehat{E}\star\widehat{\phi}^{\prime})|^{2}$) and use the fact that the correlation between $E$ and $B$ modes vanishes.

The evaluation of the above equation is simplified by introducing an operator, $\langle\cdots\rangle_{\bm{i}}$, which is the ensemble average over quantities in $i$th set of realizations. For example, denoting $E_{i}$ and $B_{i}$ as the $E$ and $B$ modes in $i$th set of simulation, we can rewrite the last term of Eq. (21) as

[TABLE]

where the index ($i=1,2,3$) denotes independent set of realizations. Substituting the above term into Eq. (20), we obtain

[TABLE]

Here we again omit the multipole. Substituting Eq. (21) into Eq. (20), and including the connected part of the six point correlation, we obtain

[TABLE]

Here $\langle\cdots\rangle_{i,\rm c}$ extracts the connected part of the ensemble average. This expression is convenient to evaluate the bias using simulation. In the above equation, the first four terms are the product of the trispectrum and power spectrum, and the next two terms are the product of the three power spectra. The last term contains the pentaspectrum.

Fig. 3 shows the power spectra of each contribution in the four point and six point delensing bias. The most significant contribution comes from the disconnected part of the four point correlation $N^{(4),\rm d}$. Some terms involved in $N^{(6)}$ are also significant. The last term in Eq. (24) is negligible.

Fig. 4 plots the dominant contribution of the delensing bias which is given by

[TABLE]

$N^{(4),\rm c}$ cancels with a term involved in $2\langle\langle E_{1}\star\widehat{\phi}^{E_{2}B_{2}}(E_{2}\star\widehat{\phi}^{E_{1}B_{2}})^{*}\rangle_{\bm{1}}\rangle_{2,\rm c}$. The above expression can be derived by rewriting the delensed $B$ modes as Namikawa and Nagata (2014)

[TABLE]

Because the third term is proportional to $\widehat{B}$ Namikawa and Nagata (2014), the bias depends on the realization of the observed $B$ mode. The primary $B$ modes which we want to detect are also simultaneously subtracted.

III.4 Mitigating delensing bias

Let us discuss how to mitigate the delensing bias.

Several works have discussed a method to avoid the delensing bias by splitting multipole range of CMB into annuli in $\ell$ space so that CMB maps to be delensed and that used for the lensing reconstruction are uncorrelated Teng et al. (2011); Sehgal et al. (2016). Fig. 5 shows the delensed $B$-mode spectra using the $\ell$-splitting method of LABEL:Teng:2011xc where CMB multipole is split into four annuli with equal spacing between $200$ and $3000$. We also show the delensing efficiency using the standard and the $\ell$-splitting methods where we use the random Gaussian reconstruction noise to avoid the delensing bias. The efficiency of delensing is defined as the ratio of the residual lensing $B$-mode spectrum to the original lensing-$B$ mode spectrum, indicating that how significantly the lensing $B$ mode is removed. In the $\ell$-splitting method, the delensing bias (difference between cyan and magenta points) is not significant. However, comparing the standard method (black) with the $\ell$-splitting technique (blue) in the random $\phi$ noise cases, the use of the $\ell$-splitting method decreases the delensing efficiency. This is because the CMB multipoles used for the lensing reconstruction decrease and the precision of the lensing potential reconstruction is degraded.

As an alternative method, in the following sections, we discuss the RD method which utilizes observed CMB maps to estimate the delensed $B$-mode power spectrum without degradation of the delensing efficiency and is more reliable than using simulation alone. This method emerges naturally as the delensing bias depends on observed $B$ mode. To construct an estimator, we employ the similar analogy of the RD method for reconstructing the lensing potential Namikawa et al. (2013) which is motivated by the optimal trispectrum estimator. We extend this technique to the delensing bias by deriving the optimal estimator for the polyspectra in Sec. IV, and applying the optimal polyspectra estimator to the disconnected part of the delensing bias which has most significant contribution in the bias in Sec. V.

IV Optimal estimator for polyspectra

IV.1 Derivation

We consider a set of variables $\bm{x}=(x_{1},x_{2},\dots,x_{N})$. For example, the index of $x_{i}$ denotes the multipoles ($\bm{\ell}$) in the flat sky analysis and the harmonic coefficients ($\ell m$) in the full sky analysis, and specifies the temperature, $E$ or $B$ modes. Expanding the likelihood of these variables around the Gaussian likelihood, a term which contains $n$th order of the correlation is given by (e.g., Refs. Amendola (1996); Regan et al. (2010))

[TABLE]

where $k_{i_{1}\dots i_{n}}$ is the coefficients of the $n$th order derivative, $i_{1},\cdots,i_{n}\in[1,N]$, and the operator $\sum_{i_{j}}$ sums over all possible combinations of $i_{1},\cdots,i_{n}$. For example, if $n=6$, $k$ is the sum of the six point correlation and products of the three point correlations. The Gaussian likelihood is given as

[TABLE]

where $c_{ij}=(C^{-1})^{x_{i}x_{j}}$ with $C$ denoting the covariance between $\bm{x}$.

To derive an explicit expression of the derivative of the Gaussian likelihood with respect to $x_{i_{1}},x_{i_{2}},\dots,x_{i_{n}}$, we first define the following quantities:

[TABLE]

Here $m$ is a non-negative integer and $p^{n}_{m}=0$ if $m>[n/2]$. Again, to simplify the numerical calculation, we introduce the operator, $\langle\cdots\rangle_{\bm{i}}$, which takes the ensemble average over $i$th set of realizations, $\overline{x}^{\bm{i}}$, and $\overline{x}$ is an inverse-variance filtered variable fields as

[TABLE]

For example, if $n=6$, we obtain

[TABLE]

and $p^{6}_{m}=0$ for $m>3$. Using

[TABLE]

the quantities $p^{n}_{m}$ satisfy

[TABLE]

With $p^{n}_{m}$, we find that the $n$th order derivative of the Gaussian likelihood is simply expressed as

[TABLE]

To show this, we first consider the case with $n=1$. The left-hand side of the above equation becomes

[TABLE]

On the other hand, the quantity defined in Eq. (29) is given by $p^{1}_{0}=\overline{x}_{i_{1}}$ and the right-hand side of Eq. (37) coincides with Eq. (38). To show the case with $n$ greater than $1$, we differentiate Eq. (37) with respect to $x_{i_{n+1}}$ and obtain

[TABLE]

Replacing $m$ with $m+1$ in the second term, we obtain

[TABLE]

Using Eq. (36), the first term becomes

[TABLE]

Substituting Eq. (41) into the first term of Eq. (40) leads to Eq. (37) and Eq. (37) is satisfied for any $n$.

The optimal estimator of the $n$th order polyspectra, $k$, is proportional to the derivative of the likelihood with respect to $k$ as similar to the trispectrum Namikawa et al. (2013). From Eq. (37), the estimator of $k$ is given as

[TABLE]

where $A_{i_{1}\cdots i_{n}}$ is a normalization of the estimator so that the estimator is unbiased and $\overline{k}_{i_{1}\cdots i_{n}}$ is the unnormalized estimator. The second term in the parenthesis is the estimator of the disconnected part of $k$:

[TABLE]

If the products of the polyspectra in $k$ are negligible (e.g., the product of the three-point correlations in $k^{i_{1}\cdots i_{6}}$), Eq. (42) becomes the estimator for the connected part of the $n$-point correlation.

Note that an optimal estimator for the amplitude of the $n$-point correlation is defined as

[TABLE]

where the variance of the $n$-point correlation, $\sigma^{2}_{i_{1}\cdots i_{n}}$, is usually assumed to be the product of the two point-correlations.

IV.2 Suboptimal estimator

Here we discuss the accuracy of the estimator given in Eq. (42). The estimator requires a set of simulated data, $\overline{x}^{\bm{i}}$. In the trispectrum case ($n=4$), when we use an incorrect covariance $\langle\overline{x}_{i}\overline{x}_{j}\rangle=\langle\overline{x}_{i}\overline{x}_{j}\rangle^{\rm true}+\delta\langle\overline{x}_{i}\overline{x}_{j}\rangle$ in the simulation, the estimator of Eq. (42) has no contribution at linear order of $\delta\langle\overline{x}_{i}\overline{x}_{j}\rangle$ Namikawa et al. (2013). Thus, the estimator is less sensitive to a mismatch between the simulated and observed covariance and more accurate than correcting the disconnected bias with the simulation alone.

In general cases, however, the estimator has contributions from linear order of $\delta\langle\overline{x}_{i}\overline{x}_{j}\rangle$. For example, in $n=6$, the following terms involved in Eq. (32) produce biases:

[TABLE]

A possible way to mitigate these terms is to redefine $p^{6}_{1}$ as

[TABLE]

The second term is added so that its ensemble average cancels with Eq. (46), while the third term is introduced such that the additional terms have zero mean. Combining $p^{6}_{m}$ in Eqs. 32 and 33 and (34), the estimator for the disconnected part of the polyspectra is given by

[TABLE]

Note that we discuss an alternative way to derive the above estimator in Appendix B by generalizing the discussion in LABEL:Carron:2017 to the polyspectra. The above estimator has no contributions from linear order of $\delta\langle\overline{x}_{i}\overline{x}_{j}\rangle$ and also from error of the trispectrum at linear order, $\delta\langle\overline{x}_{i}\overline{x}_{j}\overline{x}_{k}\overline{x}_{l}\rangle_{\rm c}$. Thus, the suboptimal estimator of Eq. (48) is more accurate than the estimator defined in Eq. (44).

V Realization-dependent delensing

We here discuss the RD estimators for the removal of the delensing bias. Since the dominant contribution to the delensing bias comes from the disconnected part, we apply the above optimal polyspectra estimators given in Eqs. 44 and 48 to the disconnected part of the delensing bias. Since the derivation requires a lengthy calculation, here we only show the results and the details of the calculation are given in Appendix C.

From Eqs. 14 and 16, the bias is described as

[TABLE]

$N^{(4)}$ is a four point correlation and its disconnected part is the most dominant contribution in the delensing bias. The disconnected part of $N^{(4)}$ can be replaced with the following RD estimator (see Appendix C.1);

[TABLE]

Here we define $X_{i,\pm j\dots}=X_{i}\pm X_{j}+\dots$ and $i=0$ means the real data. $N^{(6)}$ is a six point correlation and its disconnected part is the dominant source of the bias. According to the optimal polyspectra estimator, the RD estimator for the disconnected part of the six point correlation is given by $\widehat{N}^{(6)}\propto p^{6}_{1}-p^{6}_{2}+p^{6}_{3}$ where $p^{6}_{i}$ is defined in Eq. (29) with four $E$ and two $B$ modes. We find that the RD estimator for the disconnected part of the six point correlation becomes (see Appendix C.2)

[TABLE]

where we define $2\widehat{\phi}_{ij}=\widehat{\phi}^{E_{i}B_{j}}+\widehat{\phi}^{E_{j}B_{i}}$ and

[TABLE]

As discussed in the previous section, the suboptimal estimator of Eq. (48) is more accurate than the estimator of Eq. (51). The suboptimal RD estimator for the six point correlation is given by (see Appendix C.3)

[TABLE]

where we define

[TABLE]

As shown in Fig. 4, the dominant terms of the delensing bias are given by Eq. (25). The RD estimator for Eq. (25) can be constructed by extracting the corresponding terms in the suboptimal estimator. The details are given in Appendix C and the result is

[TABLE]

The above RD estimators are more reliable than the use of the delensing bias estimated from simulation alone. In the followings, we check how the use of $\widehat{N}^{(6),\rm d,sub}$ or $\widehat{N}^{(6),\rm d,dom}$ degrades the performance of the delensing.

Fig. 6 shows the difference between the $B$-mode power spectra using the RD estimators and with the random Gaussian reconstruction noise, divided by the $1\sigma$ simulation error. After the correction of the delensing bias with the RD estimators, the delensed spectrum becomes consistent with that in the case with the random Gaussian reconstruction noise. Fig. 7 shows the $1\sigma$ statistical error of the delensed $B$-mode spectrum after the bias correction. The statistical errors are almost unchanged between the three RD estimators. We conclude that the use of the suboptimal RD estimators ($\widehat{N}^{\rm sub}$ or $\widehat{N}^{\rm dom}$) does not degrade the performance of the delensing.

VI Summary and discussion

We have derived the RD estimators to correct the delensing bias. We first showed that the delensing bias is significant in the delensed $B$-mode power spectrum. To correct the bias, we derived the optimal polyspectra estimator in general cases. We then formulated the RD estimators for the delensing bias with the polyspectra estimator. Unlike the $\ell$ splitting method proposed in previous works, the RD method corrects the delensing bias without substantial loss of the delensing efficiency.

A similar bias on $B$ modes discussed in this paper would be appeared when we remove, for example, the spatially varying components of foregrounds estimated from the same realization of CMB maps with a quadratic estimator. The RD method presented in this paper would be also useful to mitigate such biases.

We have focused on the $B$-mode delensing but the delensing bias is also significant in the CMB temperature and $E$-mode delensing Sehgal et al. (2016); Carron et al. (2017) (see also Appendix A). Unlike the $B$-mode delensing focused in this paper, the higher-order terms of the lensing potential is not negligible in the temperature an $E$-mode delensing, and the remapping method or the template method including the higher-order terms of the lensing potential is more efficient than the linear-order template method. The delensing bias contains non-negligible contributions from the higher-order correlations (e.g. eight-point correlation of the observed CMB maps). To apply the RD method described in this paper to the temperature and $E$-mode delensing, we need to include the optimal polyspectra up to more than six point correlations.

In this paper, we have focused on the application of the optimal polyspectra estimator to the delensing. The estimator is also useful to estimate, e.g., the bispectrum of the reconstructed lensing potential which is a six point correlation of the observed CMB anisotropies, and higher-order statistics of the kinetic SZ effects. The lensing bispectrum is generated by the nonlinear growth of the large-scale structure and the post-Born corrections. Since measurements of the lensing bispectrum will be useful to extract additional cosmological information in future CMB experiments Namikawa (2016); Boehm et al. (2016); Pratten and Lewis (2016); Liu et al. (2016), we will investigate the feasibility of the RD method for a lensing bispectrum measurement in our future work.

Acknowledgements.

We would like to thank Chao-Lin Kuo for support of this work and Antony Lewis for helpful comments. This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231.

Appendix A Delensing bias in other cases

In this section, we show the delensing biases in the delensed temperature, $E$, and $B$ mode auto power spectra and temperature-$E$-mode cross power spectrum. We reconstruct lensing potentials from $\Theta\Theta$, $\Theta E$, $\Theta B$, $EE$, and $EB$ quadratic estimators of LABEL:Hu:2001kj. The temperature, $E$, and $B$ modes are delensed by the lensing potentials with the remapping method.

Fig. 8 shows the delensing biases using a lensing potential reconstructed from these five estimators. The delensing biases are defined as $\Delta C_{\ell}=C_{\ell}^{\rm res,\phi^{\rm inp}}-C_{\ell}^{\rm res}$ where $C_{\ell}^{\rm res}$ and $C_{\ell}^{\rm res,\phi^{\rm inp}}$ are the delensed spectra using a reconstructed lensing potential and the input lensing potential with the random Gaussian reconstruction noise, respectively. We also show the difference between the lensed and unlensed power spectra. In the delensed $C^{\Theta\Theta}_{\ell}$, $C^{\rm\Theta E}_{\ell}$, $C^{\rm EE}_{\ell}$, and $C^{\rm BB}_{\ell}$, the most significant delensing biases come from $\Theta\Theta$, $\Theta E$, $EE$, and $EB$ quadratic estimators, respectively.

Appendix B Derivation of the suboptimal polyspectra estimator

Here we derive Eq. (48) by generalizing the discussion of LABEL:Carron:2017 to the case with the six point correlation. In the six point correlation, errors in the covariance and trispectrum could lead to errors in the estimate of the delensing bias. To remove the errors, we consider the following estimator:

[TABLE]

where we define $c_{ij}=\langle\overline{x}_{i}\overline{x}_{j}\rangle$, $t_{jklm}=\langle\overline{x}_{j}\overline{x}_{k}\overline{x}_{l}\overline{x}_{m}\rangle_{\rm c}$, and

[TABLE]

The sum of $K^{(2,4)}$ and $K^{(2,2,2)}$ is the disconnected part of the six point correlation. The expectation of the above estimator is equal to the first term, $K^{(2,4)}+K^{(2,2,2)}$. If there are errors in the covariance and trispectrum, $\delta c$ and $\delta t$, the above estimator does not contain biases from the linear order of $\delta c$ and $\delta t$.

To simplify the above estimator, we use

[TABLE]

Substituting the above equations into Eq. (56), we find

[TABLE]

This estimator coincides with Eq. (48).

Appendix C Derivation of the RD estimator for delensing bias

In this section, we derive the RD estimator in general cases of the CMB internal delensing. We consider delensed CMB anisotropies $\widehat{X}$ using a reconstructed lensing potential $\widehat{\phi}^{YZ}$ obtained from observed CMB anisotropies, $\widehat{Y}$ and $\widehat{Z}$.

C.1 Four point correlation

We first discuss the RD estimator for the disconnected part of the four point correlation $N^{(4),d}$. In the CMB internal delensing, the delensing bias from the disconnected part of the four point correlation is given by

[TABLE]

From Eq. (44), the optimal trispectrum estimator for the disconnected part of $-2(\widehat{X}\star\widehat{\phi}^{YZ})\widehat{W}^{*}$ is given as $p^{4}_{1}-p^{4}_{2}$ where $p^{4}_{m}$ is defined in Eq. (29). The quantity, $p^{4}_{1}$, becomes

[TABLE]

Here we define a symmetric estimator,

[TABLE]

with $i=0$ meaning the real data, and ignore the mean-field bias because its contribution is absorbed into the mean-field corrected estimator $\widehat{\phi}=\widehat{\phi}^{YZ}-\langle\widehat{\phi}\rangle$ after combining $p_{m}$. We then replace $\widehat{\phi}^{YZ}$ with $\widehat{\phi}$. The last term is the correction for the mean-field bias in $(\widehat{X}\star\widehat{\phi})\widehat{W}^{*}$. The term, $p^{4}_{2}$, becomes

[TABLE]

Using Eqs. (66) and (68), the optimal estimator for $N^{(4),d}$ is written as

[TABLE]

where we define $X_{i,\pm j\dots}=X_{i}\pm X_{j}+\dots$.

The RD estimator of Eq. (50) is derived by substituting $X=Y=E$, $Z=W=B$ into Eq. (69). For the temperature delensing with the temperature quadratic estimator, the RD estimator is obtained by substituting $X=Y=Z=W=\Theta$ into Eq. (69) and the result is

[TABLE]

where $\widehat{\phi}^{\Theta\Theta}$ is the temperature quadratic estimator Hu (2001).

C.2 Six point correlation

Next we consider the estimator for the six point correlation. For clarity, we omit the operation, $\langle\cdots\rangle_{\bm{i}}$. In the CMB internal delensing using the template method, the six point correlation is written as $|\widehat{X}\star\widehat{\phi}^{YZ}|^{2}$. We first derive the estimator based on Eq. (42) and the suboptimal estimator of Eq. (48) is derived in the next section.

According to the optimal estimator of Eq. (42), the RD estimator for the disconnected part of the six point correlation is given as $p^{6}_{1}-p^{6}_{2}+p^{6}_{3}$. The quantity, $p^{6}_{1}$, has $15$ terms and becomes

[TABLE]

Using the symmetric estimator of Eq. (67), the above quantity is reduced to

[TABLE]

The quantity, $p^{6}_{2}$, consists of $45$ terms. The terms which do not have $\langle\widehat{\phi}\rangle$ are given by

[TABLE]

To reduce the above terms, we first rewrite the above equation as

[TABLE]

Using the symmetric estimator of Eq. (67), we obtain

[TABLE]

where we use $\widehat{X}\star\widehat{\phi}^{Y_{1}Z_{2}}(\widehat{X}\star(\widehat{\phi}^{Y_{2}Z_{1}}+\widehat{\phi}^{Y_{1}Z_{2}}))^{*}=2|\widehat{X}\star\widehat{\phi}_{12}|^{2}$. The quantity, $p^{6}_{3}$, is given by

[TABLE]

Combining the above quantities, we obtain

[TABLE]

Using $\widehat{X}_{i,\pm j,\dots}=\widehat{X}_{i}\pm\widehat{X}_{j}+\dots$ with $i=0$ denoting the real data, we find

[TABLE]

Here we use the fact that the left-hand side of the above fourth equation is real, and we also exchange the index ($1\leftrightarrow 3$) in the right-hand side of the last equation. The sum of the above terms becomes

[TABLE]

The other terms in Eq. (77) except $|\widehat{X}_{1}\star\widehat{\phi}|^{2}$ are summarized as

[TABLE]

The sum of Eqs. 79 and 80 becomes

[TABLE]

where we define a realization dependent $B$ mode as

[TABLE]

The estimator for the disconnected part of the six point correlation is recast as

[TABLE]

To obtain the delensed power spectrum of Eq. (12), we need to add the term $|\widehat{X}_{2}\star\widehat{\phi}_{11}|^{2}$ to the above estimator because $\widehat{N}$ is the estimator for the all disconnected parts of the six point correlation and $|\widehat{X}_{2}\star\widehat{\phi}_{11}|^{2}$ is involved in $\widehat{N}$. To recover this term, we add

[TABLE]

This correction is motivated by the following estimator Carron et al. (2017)

[TABLE]

where $x$ is any quantity and $f$ is a function of $x$. Subtracting Eq. (84) from $\widehat{N}$, we obtain

[TABLE]

The RD estimator of Eq. (51) is derived by substituting $X=Y=E$, $Z=B$ into Eq. (86).

C.3 Six point correlation with suboptimal estimator

Next we derive the suboptimal estimator for the six point correlation. From Eq. (48), the suboptimal estimator is given by

[TABLE]

Using Eqs. 72, 75 and 76, we obtain

[TABLE]

The last term of Eq. (87) is given by

[TABLE]

Combining Eqs. 88 and 89, we obtain

[TABLE]

where $|\widehat{X}_{1}\star\widehat{\phi}|^{2}+|\widehat{X}\star\widehat{\phi}_{11}|^{2}-|\widehat{X}_{1}\star\widehat{\phi}_{22}|^{2}$ is subtracted from the above equation.

To simplify the above equation, we use

[TABLE]

Using Eq. (78) and the above equation (with exchanging the index, $1\leftrightarrow 2$), we obtain

[TABLE]

The sum of the second and third lines of the above equation is given by

[TABLE]

Substituting the above equation into Eq. (92), we find

[TABLE]

where we define a realization-dependent $B$ mode as

[TABLE]

Eq. (53) is derived with $X=Y=E$ and $Z=B$ in Eq. (94).

C.4 Six point correlation relevant to the dominant term of the delensing bias

Here we consider the RD estimator for Eq. (25). To extract the corresponding terms of Eq. (25) in Eq. (48), we first rewrite Eq. (48) as

[TABLE]

Here $\overline{k}_{\i_{3}i_{4}i_{5}i_{6}}$ is the unnormalized estimator of the trispectrum defined in Eq. (43).

We first consider the term containing the connected four point correlation in Eq. (25). The corresponding terms in Eq. (96) are given by

[TABLE]

The permutations of the above quantity are not relevant to Eq. (25). The above equation is recast as

[TABLE]

From the above equation, the RD estimator for the second term of Eq. (25) is obtained by replacing $\overline{x}_{i_{1}}\to\widehat{E}$, $\overline{x}_{i_{2}}\overline{x}_{i_{3}}\to\widehat{\phi}^{EB}$, $\overline{x}_{i_{4}}\to\widehat{E}$ and $\overline{x}_{i_{5}}\overline{x}_{i_{6}}\to\widehat{\phi}^{EB}$. The result is

[TABLE]

To simplify the above equation, we use

[TABLE]

Note that the ensemble average of the above first and third equations is zero while that of the second equation gives the connected part in Eq. (25). With the above equations, Eq. (100) becomes

[TABLE]

Next we consider the full disconnected part of the six point correlation (i.e., the third term) involved in Eq. (25). The corresponding terms in Eq. (96) are described by

[TABLE]

Replacing $\overline{x}_{i_{1}}\to\widehat{E}$, $\overline{x}_{i_{2}}\overline{x}_{i_{3}}\to\widehat{\phi}^{EB}$, $\overline{x}_{i_{4}}\to\widehat{E}$ and $\overline{x}_{i_{5}}\overline{x}_{i_{6}}\to\widehat{\phi}^{EB}$, we obtain the RD estimator for the third term of Eq. (25) as

[TABLE]

Combining $\widehat{N}^{\rm dom,1}$ and $\widehat{N}^{\rm dom,2}$, the RD estimator of Eq. (25) is given as

[TABLE]