# Quadratic estimators for CMB weak lensing

Abhishek S. Maniyar<EMAIL_ADDRESS>Center for Cosmology and
Particle Physics, Department of Physics, New York University, New York, NY
10003, USA  Yacine Ali-Haïmoud  Center for Cosmology and Particle Physics,
Department of Physics, New York University, New York, NY 10003, USA  Julien
Carron  Université de Genève, Département de Physique Théorique et CAP, 24
Quai Ansermet, CH-1211 Genève 4, Switzerland  Antony Lewis  Department of
Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK  Mathew S.
Madhavacheril  Perimeter Institute for Theoretical Physics, Waterloo, ON,
Canada N2L 2Y5

###### Abstract

In recent years, weak lensing of the cosmic microwave background (CMB) has
emerged as a powerful tool to probe fundamental physics, such as neutrino
masses, primordial non-Gaussianity, dark energy, and modified gravity. The
prime target of CMB lensing surveys is the lensing potential, which is
reconstructed from the observed CMB temperature $T$ and polarization $E$ and
$B$ fields. Until very recently, this reconstruction has been performed with
quadratic estimators (QEs), which, although known to be suboptimal for high-
sensitivity experiments, are numerically efficient, and useful to make
forecasts and cross-check the results of more sophisticated likelihood-based
methods. It is expected that ongoing and near-future CMB experiments such as
AdvACT, SPT-3G and the Simons Observatory (SO), will also rely on QEs. In this
work, we review different QEs, and clarify and quantify their differences. In
particular, we show that the Hu-Okamoto (HO02) estimator is not the absolute
optimal lensing estimator that can be constructed out of quadratic
combinations of $T,E$ and $B$ fields. Instead, we derive the global-minimum-
variance (GMV) lensing quadratic estimator. Although this estimator can be
found elsewhere in the literature, it was erroneously described as equivalent
to the HO02 estimator, and has never been used in real data analyses. Here, we
show explicitly that the HO02 estimator is suboptimal to the GMV estimator,
with a reconstruction noise larger by up to $\sim 9\%$ for a SO-like
experiment. We further show that the QE used in the Planck, and recent SPT
lensing analysis is suboptimal to both the HO02 and GMV estimator, and would
have a reconstruction noise up to $\sim 11\%$ larger than that of the GMV
estimator for a SO-like experiment. In addition to clarifying differences
between different QEs, this work should thus provide motivation to implement
the GMV estimator in future lensing analyses relying on QEs.

##  I Introduction

Weak gravitational lensing of the cosmic microwave background (CMB) arises
from the deflection of CMB photons as they travel to us from the last
scattering surface, through the inhomogeneous Universe Blanchard and Schneider
(1987); see e.g. Ref. Lewis and Challinor (2006) for a review. The deflection
angle is proportional to the gradient of the lensing potential $\phi$, which
is determined by the projected mass distribution along the line of sight.
Reconstructing $\phi$ is therefore a powerful cosmological tool, as it gives
direct access to the projected distribution of the _total_ matter – baryonic
and dark – without relying on biased tracers Seljak and Zaldarriaga (1999).
Among other applications, the power spectrum of the lensing potential and its
cross-correlation with other tracers of large-scale structure are a sensitive
probe of the growth of matter fluctuations, primordial non-Gaussianity,
neutrino masses, dark energy, and modified gravity Lewis and Challinor (2006);
Allison  _et al._ (2015); Schmittfull and Seljak (2018). CMB lensing has been
successfully measured by ACT, SPT, Planck, BICEP and POLARBEAR Das  _et al._
(2011); Sherwin  _et al._ (2017); van Engelen  _et al._ (2012); Ade  _et al._
(2016a); Aghanim  _et al._ (2020); Omori  _et al._ (2017); Story  _et al._
(2015); Wu  _et al._ (2019); Millea  _et al._ (2020); Ade  _et al._ (2016b,
2014); Adachi  _et al._ (2020). Current and upcoming wide-field CMB
experiments such as AdvACT Henderson  _et al._ (2016), SPT-3G Benson  _et al._
(2014) and the Simons Observatory (SO) Ade  _et al._ (2019) will measure the
lensing potential with even higher signal-to-noise ratio. Looking ahead, next-
generation “Stage-4” instrumental concepts with unprecedented depth and
angular resolution are currently under development Abazajian  _et al._ (2016),
with CMB lensing as one of their main science goals Abazajian  _et al._
(2019).

One of the main signatures of weak lensing is the induced correlations between
unequal Fourier modes of the CMB temperature and polarization fields. It is
therefore natural to seek to estimate $\phi$ out of linear combinations of
terms quadratic in different modes of the observed fields Zaldarriaga and
Seljak (1999); Hu (2001); and indeed, almost all CMB lensing analyses thus far
have relied on such quadratic estimators. For the next-generation Stage-4-like
CMB experiments (CMBS4), quadratic estimators are known to be suboptimal,
especially for polarization Hirata and Seljak (2003a). More elaborate
algorithms are being developed, such as the gradient-inversion method
Hadzhiyska  _et al._ (2019), or likelihood-based methods Hirata and Seljak
(2003b, a); Carron and Lewis (2017); Millea  _et al._ (2020). Meanwhile,
quadratic estimators remain the workhorse tool for current and near-future CMB
experiments like AdvACT, SPT-3G, and SO. They have the advantages of being
very simple to implement and computationally efficient, and will serve as
useful cross-checks even when more accurate and computationally demanding
methods are employed with future data.

The main goal of this paper is to clarify and quantify the differences between
several quadratic estimators commonly used for CMB lensing reconstruction. Our
most important point is that the well-known Hu and Okamoto Hu and Okamoto
(2002) (hereafter, HO02) estimator is _not_ the optimal quadratic estimator
that can be constructed from temperature and polarization maps, even if
generalized to the full sky, and even when using non-perturbative response
functions Lewis  _et al._ (2011); Fabbian  _et al._ (2019). Instead, we derive
the global-minimum-variance (hereafter GMV) quadratic estimator built out of
all possible quadratic combination of $T,E$ and $B$. The GMV estimator was in
fact first derived in Hirata and Seljak Hirata and Seljak (2003a), as the
weak-signal limit of their likelihood-based method. Nevertheless, it was
stated there and in subsequent works that this estimator is equivalent to that
of HO02. We explicitly show that this is not the case, and that the
reconstruction noise of the GMV estimator can be up to $\sim 9\%$ lower than
that of the HO02 estimator on large angular scales. We also generalize it to
be accurately unbiased accounting for higher-order lensing effects.
Furthermore, we show that the quadratic estimator used in the Planck
collaboration Ade  _et al._ (2016a); Aghanim  _et al._ (2020) and SPT
collaboration Wu  _et al._ (2019) lensing analyses, obtained by neglecting
$C^{TE}_{\ell}$ in the inverse filter matrix, is suboptimal to _both_ the GMV
and HO02 estimators. For a SO-like experiment, this suboptimal estimator is up
to $\sim 11\%$ noisier that the GMV estimator. This may motivate implementing
the GMV estimator in future analyses, despite the possible added complexity of
jointly filtering temperature and polarization maps.

The remainder of this paper is organized as follows. After introducing our
notation and convention in Sec. II, we review the HO02 estimator and its close
cousin, the Okamoto-Hu Okamoto and Hu (2003) (hereafter OH03) estimator in
Sec. III. We then derive the GMV estimator in Sec. IV and explicitly show how
it differs from the HO02 estimator. We describe the suboptimal lensing
estimator of Ref. Aghanim  _et al._ (2020) in Sec. V. Finally, we compare
estimators in Sec. VI for different instrumental setups, and conclude in Sec.
VII.

##  II Notation and conventions

We denote by capital letters $X,Y=T,E,B$ the observed (lensed and noisy) CMB
temperature and polarization fields, and by $\phi$ the projected lensing
potential. Throughout we work in the flat-sky approximation; we denote two-
dimensional Fourier wavenumbers by $\boldsymbol{l}$ for CMB fields and
$\boldsymbol{L}$ for the lensing potential.

The power spectra of the observed temperature and polarization fields are
defined as

$\displaystyle\langle
X(\boldsymbol{l})Y(\boldsymbol{l}^{\prime})\rangle=(2\pi)^{2}\delta(\boldsymbol{l}+\boldsymbol{l}^{\prime})C_{l}^{XY}\,,$
(1)

where $C_{l}^{XY}$ is the total cross-power spectrum of the lensed fields,
including detector noise added in quadrature (for $X=Y$). It can also include
contributions from other sources of variance such as residual foreground
contamination. In this expression the angular brackets denote taking ensemble
averages over the primordial CMB, detector noise, along with the underlying
large scale structure.

Gravitational lensing affects the auto- and cross-power spectra of CMB fields,
and moreover produces correlations between non-opposite $\boldsymbol{l}$
modes, proportional to the projected lensing potential. The response of the
non-opposite correlations to lensing can be quantified by non-perturbative
response functions $f_{XY}$ defined by

$\left\langle\frac{\delta}{\delta\phi(\boldsymbol{L})}\left(X(\boldsymbol{l})Y(\boldsymbol{l}^{\prime})\right)\right\rangle=\delta(\boldsymbol{l}+\boldsymbol{l}^{\prime}-\boldsymbol{L})f_{XY}(\boldsymbol{l},\boldsymbol{l}^{\prime}).$
(2)

The coupling coefficients $f_{XY}$ appearing in Eq. (2) are given explicitly
in Table 1. They depend on the lensed gradient spectra
$\widetilde{C}_{l}^{X\nabla Y}$ defined in Refs. Lewis  _et al._ (2011);
Fabbian  _et al._ (2019), which generalize the unlensed spectra used in the
original work of HO02 so that the response function for each lensing mode
includes the important higher-order effect of other lensing modes. The BB term
has negligible contribution to the signal-to-noise ratio of the reconstructed
$\phi$ field and thus we neglect it in our analysis. Note that different
foregrounds can also contribute to off-diagonal correlations (e.g. Ferraro and
Hill, 2018; Schaan and Ferraro, 2019), but we do not include them in this
work.

$\alpha$ |  $f_{\alpha}({\boldsymbol{l}_{1},\boldsymbol{l}_{2}})$ [
---|---
$TT$ |  $\widetilde{C}_{l_{1}}^{T\nabla T}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{1})+\widetilde{C}_{l_{2}}^{T\nabla T}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{2})$[
$TE$ |  $\widetilde{C}_{l_{1}}^{T\nabla E}\cos 2\varphi_{\boldsymbol{l}_{1}\boldsymbol{l}_{2}}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{1})+\widetilde{C}_{l_{2}}^{T\nabla E}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{2})$[
$EE$ |  $[\widetilde{C}_{l_{1}}^{E\nabla E}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{1})+\widetilde{C}_{l_{2}}^{E\nabla E}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{2})]\cos 2\varphi_{\boldsymbol{l}_{1}\boldsymbol{l}_{2}}$ [
$TB$ |  $\widetilde{C}_{l_{1}}^{T\nabla E}\sin 2\varphi_{\boldsymbol{l}_{1}\boldsymbol{l}_{2}}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{1})$[
$EB$ |  $[\widetilde{C}_{l_{1}}^{E\nabla E}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{1})+\widetilde{C}_{l_{2}}^{B\nabla B}({\boldsymbol{L}}\cdot{\boldsymbol{l}}_{2})]\sin 2\varphi_{\boldsymbol{l}_{1}\boldsymbol{l}_{2}}$ [
Table 1: CMB lensing correlation coefficients.
$\varphi_{\boldsymbol{l}_{1}\boldsymbol{l}_{2}}$ is the angle between
$\boldsymbol{l}_{1}$ and $\boldsymbol{l}_{2}$. The quantity
$\widetilde{C}^{X\nabla Y}$ is the lensed gradient spectrum, defined in Refs.
Lewis  _et al._ (2011); Fabbian  _et al._ (2019). Note that we do not include
curl-like terms $\widetilde{C}_{l}^{TP_{\bot}},\widetilde{C}_{l}^{PP_{\bot}}$,
which are always subdominant Fabbian  _et al._ (2019).

In the remainder of this work we will describe different estimators
$\hat{\phi}_{\alpha}$ for the lensing potential. All these estimators are
required to be unbiased, i.e. such that
$\langle\hat{\phi}_{\alpha}\rangle=\phi$. They are however noisy, and we
define their variance (or reconstruction noise) $N_{\alpha}(L)$ through

$\langle(\hat{\phi}_{\alpha}-\phi)(\boldsymbol{L})(\hat{\phi}_{\alpha}-\phi)(\boldsymbol{L}^{\prime})\rangle=(2\pi)^{2}\delta_{\rm
D}(\boldsymbol{L}+\boldsymbol{L}^{\prime})N_{\alpha}(L).$ (3)

Here, for optimizing the signal to noise, we only consider the primary
Gaussian disconnected contractions of the lensed fields,
$N_{\alpha}^{(0)}(L)$; in Appendix. A we give an explicit form for the
$N_{\alpha}^{(1)}(L)$ contractions Kesden  _et al._ (2003) that should also be
included in any full data likelihood analysis. The superscript values 0 and 1
represent the order to which the variance $N_{\alpha}(L)$ explicitly depends
on $C^{\phi\phi}_{L}$.

We will often deal with convolutions in Fourier space, and for brevity,
introduce the compact notation

$\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}...\equiv\iint\frac{d^{2}l_{1}d^{2}l_{2}}{(2\pi)^{2}}\delta_{\rm
D}(\boldsymbol{l}_{1}+\boldsymbol{l}_{2}-\boldsymbol{L})...$ (4)

We define the Fourier transform of a configuration-space function
$A(\boldsymbol{\hat{n}})$ as

$\mathcal{F}[A(\boldsymbol{\hat{n}})](\boldsymbol{l})\equiv\int
d^{2}\boldsymbol{\hat{n}}~{}A(\boldsymbol{\hat{n}})\textrm{e}^{-i\boldsymbol{l}\cdot\boldsymbol{\hat{n}}},$
(5)

and the inverse-Fourier transform of a harmonic-space function
$B(\boldsymbol{l})$ as

$\mathcal{F}^{-1}[B(\boldsymbol{l})](\boldsymbol{\hat{n}})\equiv\int\frac{d^{2}\boldsymbol{l}}{(2\pi)^{2}}B(\boldsymbol{l})\textrm{e}^{i\boldsymbol{l}\cdot\boldsymbol{\hat{n}}}.$
(6)

##  III Hu and Okamoto estimators

We now briefly rederive the HO02 and OH03 quadratic estimators for the lensing
potential, setting the stage for our subsequent derivation of the global-
minimum-variance estimator.

The approach of HO02 consists in constructing the single-pair estimators
$\hat{\phi}_{\alpha}({\boldsymbol{L}})$ separately for each pair
$\alpha=TT,TE,EE,TB,EB$, and then combining them together to form their
minimum variance estimator.

OH03 moreover derive efficient full-sky single-pair estimators in
configuration space. They are identical to the HO02 estimators, except for the
$TE$ estimator, which is slightly sub-optimal. We give explicit expressions
for these estimators in the flat-sky limit in Section III.2. Here again, the
final OH03 estimator is obtained by combining these single-pair estimators.
Note that HO02 and OH03 used unlensed spectra in the response functions rather
than the lensed gradient spectra. We still refer to the estimators constructed
with the lensed gradient spectra as the HO02 and OH03 estimators, given that
the procedure is identical.

###  III.1 Single-pair minimum-variance quadratic estimators in harmonic space

We start by constructing quadratic estimators out of a single pair $XY$, of
the form

$\hat{\phi}_{XY}(\boldsymbol{L})\equiv\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}X(\boldsymbol{l}_{1})Y(\boldsymbol{l}_{2})F_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}).$
(7)

For the estimator to be unbiased, the weights $F_{XY}$ must satisfy the
constraint

$\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})F_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=1.$
(8)

The noise of this estimator (defined as in Eq. (3)) is

$\displaystyle
N_{XY}(L)=\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}F_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\Big{(}F_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})C_{l_{1}}^{XX}C_{l_{2}}^{YY}$
$\displaystyle+F_{XY}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})C_{l_{1}}^{XY}C_{l_{2}}^{XY}\Big{)}.~{}~{}~{}$
(9)

####  III.1.1 All pairs except $TE$

For all pairs except $TE$, either $X=Y$ or $C_{l}^{XY}$ = 0. As a consequence,
the variance of the estimator takes the form

$\displaystyle
N_{XY}(L)=(1+\delta_{XY})\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}C_{l_{1}}^{XX}C_{l_{2}}^{YY}[F_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})]^{2}.~{}~{}~{}~{}$
(10)

Minimizing this variance under the constraint (8) results in the following
coefficients

$\displaystyle F_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ $\displaystyle=$
$\displaystyle\lambda_{XY}(L)~{}\frac{f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})}{(1+\delta_{XY})C_{l_{1}}^{XX}C_{l_{2}}^{YY}},$
(11) $\displaystyle\lambda_{XY}(L)$ $\displaystyle\equiv$
$\displaystyle\left[\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\frac{[f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})]^{2}}{(1+\delta_{XY})C_{l_{1}}^{XX}C_{l_{2}}^{YY}}\right]^{-1}.$
(12)

Inserting back into Eq. (10), we find the corresponding minimum variance
$N_{XY}(L)=\lambda_{XY}(L)$.

####  III.1.2 Special case of $XY=TE$

We may decompose $F_{TE}$ into a symmetric and antisymmetric piece:

$\displaystyle F_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ $\displaystyle=$
$\displaystyle
F_{TE}^{+}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})+F_{TE}^{-}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}),$
(13) $\displaystyle F_{TE}^{\pm}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$
$\displaystyle\equiv$
$\displaystyle\frac{1}{2}\left(F_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\pm
F_{TE}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})\right).$ (14)

For each pair $(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$, we further define the
2-dimensional vector

$\boldsymbol{F}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\equiv\left(F_{TE}^{+}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}),F_{TE}^{-}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\right).$
(15)

After some algebra, and only keeping even functions of
$(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ in the integral, Eq. (9) can be
rewritten as

$N_{TE}(L)=\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\boldsymbol{F}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\cdot\boldsymbol{M}(l_{1},l_{2})\cdot\boldsymbol{F}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}),$
(16)

where for each pair $(l_{1},l_{2})$, the 2 by 2 matrix
$\boldsymbol{M}(l_{1},l_{2})$ is given by

$\displaystyle\boldsymbol{M}(l_{1},l_{2})=$
$\displaystyle\begin{pmatrix}C_{(l_{1}}^{TT}C_{l_{2})}^{EE}+C_{l_{1}}^{TE}C_{l_{2}}^{TE}&C_{[l_{1}}^{TT}C_{l_{2}]}^{EE}\\\
C_{[l_{1}}^{TT}C_{l_{2}]}^{EE}&C_{(l_{1}}^{TT}C_{l_{2})}^{EE}-C_{l_{1}}^{TE}C_{l_{2}}^{TE}\end{pmatrix},~{}~{}~{}$
(17)

where $A_{(l_{1}l_{2})}\equiv(A_{l_{1}l_{2}}+A_{l_{2}l_{1}})/2$ and
$A_{[l_{1}l_{2}]}\equiv(A_{l_{1}l_{2}}-A_{l_{2}l_{1}})/2$ are the symmetric
and antisymmetric parts of $A_{l_{1}l_{2}}$.

Similarly, we may define the symmetric and antisymmetric parts of the
correlation coefficients $f_{TE}^{\pm}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$
and the two-dimensional vector $\boldsymbol{f}=(f_{TE}^{+},f_{TE}^{-})$, for
each pair $(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$, and rewrite the
constraint (8) as

$\displaystyle\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\boldsymbol{F}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\cdot\boldsymbol{f}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=1.$
(18)

Minimizing the variance (16) under this constraint leads to the solution

$\boldsymbol{F}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=\lambda(L)~{}\boldsymbol{M}^{-1}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\cdot\boldsymbol{f}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}),$
(19)

where the Lagrange multiplier $\lambda$ is obtained from the constraint (18):

$\lambda(L)=\left(\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\boldsymbol{f}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\cdot\boldsymbol{M}^{-1}(l_{1},l_{2})\cdot\boldsymbol{f}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\right)^{-1}.$
(20)

The 2 by 2 matrix $\boldsymbol{M}(l_{1},l_{2})$ is easily invertible, and
after re-expressing Eq. (19) in terms of the original
$F_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ and
$f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$, one recovers the HO02 optimal
weights for $TE$, namely, with our notation,

$\displaystyle F_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ $\displaystyle=$
$\displaystyle\lambda_{TE}(L)\frac{C_{l_{1}}^{EE}C_{l_{2}}^{TT}f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})-C_{l_{1}}^{TE}C_{l_{2}}^{TE}f_{TE}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})}{C_{l_{1}}^{TT}C_{l_{2}}^{EE}C_{l_{1}}^{EE}C_{l_{2}}^{TT}-\left(C_{l_{1}}^{TE}C_{l_{2}}^{TE}\right)^{2}},$
(21) $\displaystyle\lambda_{TE}(L)$ $\displaystyle\equiv$
$\displaystyle\Bigg{[}\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\frac{C_{l_{1}}^{EE}C_{l_{2}}^{TT}f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})-C_{l_{1}}^{TE}C_{l_{2}}^{TE}f_{TE}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})}{C_{l_{1}}^{TT}C_{l_{2}}^{EE}C_{l_{1}}^{EE}C_{l_{2}}^{TT}-\left(C_{l_{1}}^{TE}C_{l_{2}}^{TE}\right)^{2}}\Bigg{]}^{-1}\,.$
(22)

Inserting Eq. (19) into Eq. (16), we see that the noise of the minimum-
variance estimator is just $N_{TE}(L)=\lambda_{TE}(L)$.

###  III.2 Single-pair efficient configuration-space estimators

####  III.2.1 All pairs except $TE$

The response coefficients $f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ can
all be written as linear combinations of products of functions of
$\boldsymbol{l}_{1}$ with functions of $\boldsymbol{l}_{2}$, with coefficients
depending on $L$. From Eq. (11), we see that this property transfers to the
optimal weights $F_{XY}$ for all pairs except $TE$. As a consequence, all
single-pair estimators except $TE$ can be written as sums of convolutions of
functions of $\boldsymbol{l}_{1}$ with functions of $\boldsymbol{l}_{2}$. This
implies that they can be written as a sum of products of functions of
configuration space – they are “separable” in configuration space. This allows
to use Fast Fourier Transforms (FFTs) (or fast harmonic transforms for full-
sky expressions Okamoto and Hu (2003)) to compute them efficiently.

Similar to OH03, we define the following bilinear operator of harmonic-space
functions:

$\boldsymbol{\mathcal{P}}[A(\boldsymbol{l}_{1}),B(\boldsymbol{l}_{2})](\boldsymbol{\hat{n}})\equiv\boldsymbol{\nabla}\mathcal{F}^{-1}[A(\boldsymbol{l}_{1})](\boldsymbol{\hat{n}})\times\mathcal{F}^{-1}[B(\boldsymbol{l}_{2})](\boldsymbol{\hat{n}}).$
(23)

The single-pair estimators can all be written in the form

$\hat{\phi}_{XY}(\boldsymbol{\hat{n}})=-\boldsymbol{\nabla}\cdot\mathcal{F}^{-1}\left[\lambda_{XY}(L)\mathcal{F}\left[\boldsymbol{\psi}_{XY}(\boldsymbol{\hat{n}})\right]\right],$
(24)

with

$\displaystyle\boldsymbol{\psi}_{TT}$ $\displaystyle=$
$\displaystyle\boldsymbol{\mathcal{P}}\left[\frac{\tilde{C}_{l_{1}}^{T\nabla
T}}{C_{l_{1}}^{TT}}T(\boldsymbol{l}_{1}),\frac{T(\boldsymbol{l}_{2})}{C_{l_{2}}^{TT}}\right],$
(25) $\displaystyle\boldsymbol{\psi}_{EE}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\epsilon=\pm
1}\boldsymbol{\mathcal{P}}\left[\textrm{e}^{2i\epsilon\varphi_{\boldsymbol{l}_{1}}}\frac{\tilde{C}_{l_{1}}^{E\nabla
E}}{C_{l_{1}}^{EE}}E(\boldsymbol{l}_{1}),\textrm{e}^{-2i\epsilon\varphi_{\boldsymbol{l}_{2}}}\frac{E(\boldsymbol{l}_{2})}{C_{l_{2}}^{EE}}\right],~{}~{}$
(26) $\displaystyle\boldsymbol{\psi}_{TB}$ $\displaystyle=$
$\displaystyle\frac{1}{2i}\sum_{\epsilon=\pm
1}\epsilon\boldsymbol{\mathcal{P}}\left[\textrm{e}^{2i\epsilon\varphi_{\boldsymbol{l}_{1}}}\frac{\tilde{C}_{l_{1}}^{T\nabla
E}}{C_{l_{1}}^{TT}}T(\boldsymbol{l}_{1}),\textrm{e}^{-2i\epsilon\varphi_{\boldsymbol{l}_{2}}}\frac{B(\boldsymbol{l}_{2})}{C_{l_{2}}^{BB}}\right],~{}~{}~{}$
(27) $\displaystyle\boldsymbol{\psi}_{EB}$ $\displaystyle=$
$\displaystyle\frac{1}{2i}\sum_{\epsilon=\pm
1}\epsilon\boldsymbol{\mathcal{P}}\left[\textrm{e}^{2i\epsilon\varphi_{\boldsymbol{l}_{1}}}\frac{\tilde{C}_{l_{1}}^{E\nabla
E}}{C_{l_{1}}^{EE}}E(\boldsymbol{l}_{1}),\textrm{e}^{-2i\epsilon\varphi_{\boldsymbol{l}_{2}}}\frac{B(\boldsymbol{l}_{2})}{C_{l_{2}}^{BB}}\right].~{}~{}~{}$
(28)

These expressions are the flat-sky limit of the OH03 full-sky expressions.

####  III.2.2 Case of $XY=TE$

The separability property is not satisfied by the $TE$ estimator, due to the
non-factorizable term in the denominator of $F_{TE}$ in Eq. (21). Instead of
the optimal $F_{TE}$, one can use a slightly suboptimal coefficient, obtained
by setting $C_{l}^{TE}=0$ in Eq. (21), namely

$\displaystyle F^{\rm eff}_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$
$\displaystyle=$ $\displaystyle\lambda^{\rm
eff}_{TE}(L)\frac{f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})}{C_{l_{1}}^{TT}C_{l_{2}}^{EE}},$
(29) $\displaystyle\lambda^{\rm eff}_{TE}(L)$ $\displaystyle\equiv$
$\displaystyle\Bigg{[}\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\frac{[f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})]^{2}}{C_{l_{1}}^{TT}C_{l_{2}}^{EE}}\Bigg{]}^{-1}\,.$
(30)

The resulting suboptimal estimator $\hat{\phi}_{TE}^{\rm eff}$ also takes the
form Eq. (24), with

$\displaystyle\boldsymbol{\psi}_{TE}^{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{1}{2}\sum_{\epsilon=\pm
1}\boldsymbol{\mathcal{P}}\left[\textrm{e}^{2i\epsilon\varphi_{\boldsymbol{l}_{1}}}\frac{\tilde{C}_{l_{1}}^{T\nabla
E}}{C_{l_{1}}^{TT}}T(\boldsymbol{l}_{1}),\textrm{e}^{-2i\epsilon\varphi_{\boldsymbol{l}_{2}}}\frac{E(\boldsymbol{l}_{2})}{C_{l_{2}}^{EE}}\right]$
(31)
$\displaystyle+\boldsymbol{\mathcal{P}}\left[\frac{\tilde{C}_{l_{2}}^{T\nabla
E}}{C_{l_{2}}^{EE}}E(\boldsymbol{l}_{2}),\frac{T(\boldsymbol{l}_{1})}{C_{l_{1}}^{TT}}\right].$

###  III.3 Optimal combination of single-pair estimators

Given the five single-pair estimators $\hat{\phi}_{\alpha}({\boldsymbol{L}})$
constructed for each $\alpha\in\\{TT,TE,EE,TB,EB\\}$, HO02 combine them to
form the estimator

$\hat{\phi}_{\rm
HO02}({\boldsymbol{L}})=\sum_{\alpha}w_{\alpha}(L)\hat{\phi}_{\alpha}({\boldsymbol{L}})\,,$
(32)

where the optimal weights $w_{\alpha}(L)$ are obtained by minimizing the
variance of the linear combination with the constraint that they sum up to
unity i.e. $\sum_{\alpha}w_{\alpha}=1$. Subject to this constraint, one gets

$w_{\alpha}=\frac{\sum_{\beta}({\boldsymbol{N}}^{-1})_{\alpha\beta}}{\sum_{\beta\gamma}({\boldsymbol{N}}^{-1})_{\beta\gamma}},$
(33)

where for each $L$, $\boldsymbol{N}_{\alpha\beta}(L)$ is the covariance matrix
of the separate estimators $\hat{\phi}_{\alpha}$, whose elements are obtained
by the generalization of Eq. (3) to the cross-correlation of two estimators.
The overall noise of this estimator is then $N_{\rm
HO}\equiv\left(\sum_{\alpha\beta}({\boldsymbol{N}}^{-1})_{\alpha\beta}\right)^{-1}$.

The final HO02 estimator thus takes the form

$\hat{\phi}_{\rm
HO02}(\boldsymbol{L})=\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\sum_{XY}F_{XY}^{\rm
HO02}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})X(\boldsymbol{l}_{1})Y(\boldsymbol{l}_{2}),$
(34)

where the sum runs over the five unique pairs $XY=TT,TE,EE,TB,EB$, and the
weights are proportional to the single-pair optimal weights, each with a
different proportionality coefficient:

$F_{XY}^{\rm
HO02}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=w_{XY}(L)F_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}).$
(35)

The same procedure can be carried with the single-pair separable estimators of
OH03. These estimators are all identical to the minimum-variance estimators of
HO02, except for $\hat{\phi}_{TE}^{\rm eff}$, which is slightly suboptimal
relative to $\hat{\phi}_{TE}$. Upon combining all five estimators, the OH03
estimator also takes the form of Eq. (34), with weights

$F_{XY}^{\rm OH03}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=w_{XY}^{\rm
eff}(L)\lambda_{XY}^{\rm
eff}(L)\frac{f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})}{(1+\delta_{XY})C_{l_{1}}^{XX}C_{l_{2}}^{YY}}.$
(36)

##  IV Global minimum-variance quadratic estimator

It is easy to see that the final HO02 estimator is a linear combination of
terms quadratic in $T,E,B$. Rather than splitting the optimization process in
two steps, we instead directly seek the global minimum variance quadratic
estimator, in one single step. By doing so, we can account for the
correlations between different $XY$ pairs _for each
$(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$_, rather than only after integrating
over $(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$, as done in the HO02 estimator.
The estimator built this way is therefore necessarily less noisy than the HO02
estimator, as we will show explicitly.

###  IV.1 Harmonic-space expression

We start by deriving the global-minimum-variance (hereafter GMV) estimator in
harmonic space, following the steps of Appendix A of Hirata & Seljak Hirata
and Seljak (2003a).

For each Fourier mode $\boldsymbol{l}$, we define the three-dimensional vector
$\boldsymbol{X}(\boldsymbol{l})=[T(\boldsymbol{l}),E(\boldsymbol{l}),B(\boldsymbol{l})]$.
We seek an estimator of the form

$\hat{\phi}(\boldsymbol{L})=\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}X^{i}(\boldsymbol{l}_{1})\Xi_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})X^{j}(\boldsymbol{l}_{2}),$
(37)

where we use the Einstein summation convention. Without loss of generality, we
may assume
$\Xi_{ji}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})=\Xi_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$,
as only the part of the integrand symmetric under exchange of
$(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ contributes to the integral.

For $\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\neq\boldsymbol{0}$, we define
$f_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ through

$\left\langle\frac{\delta}{\delta\phi(\boldsymbol{L})}\left(X^{i}(\boldsymbol{l}_{1})X^{j}(\boldsymbol{l}_{2})\right)\right\rangle=\delta(\boldsymbol{l}_{1}+\boldsymbol{l}_{2}-\boldsymbol{L})f_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}).$
(38)

In other words, if $i=1,2,3$ correspond to $X^{i}=T,E,B$, we have
$f_{11}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=f_{TT}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}),f_{12}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=f_{21}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})$,
etc… Here again, we have
$f_{ji}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})=f_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$.

Requiring the estimator to be unbiased thus leads the constraint equation

$\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\Xi_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})f_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=1.$
(39)

Using the symmetry properties of $\Xi$, the variance of this estimator is then

$\displaystyle N(\boldsymbol{L})$ $\displaystyle=$ $\displaystyle
2\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\Xi_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\Xi_{pq}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})C_{l_{1}}^{ip}C_{l_{2}}^{jq}.$
(40)

Minimizing this variance under the constraint (39) leads to

$C_{l_{1}}^{ip}\Xi_{pq}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})C_{l_{2}}^{jq}=\frac{\lambda(L)}{2}f_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2}),$
(41)

where $\lambda(L)$ is a Lagrange multiplier. This equation is more easily
solved in matrix form. For each $l$, we define the three by three symmetric
matrix $[\boldsymbol{C}_{l}]$ with elements $C_{l}^{ij}$; similarly, for each
pair $(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$, we define the three by three
matrices $[\boldsymbol{\Xi}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})]$ and
$[\boldsymbol{f}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})]$. Equation (41) then
has the solution

$[\boldsymbol{\Xi}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})]=\frac{\lambda(L)}{2}[\boldsymbol{C}_{l_{1}}]^{-1}[\boldsymbol{f}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})][\boldsymbol{C}_{l_{2}}]^{-1}.$
(42)

Inserting back into the constraint equation, we obtain

$\lambda(L)^{-1}=\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\frac{1}{2}\textrm{Tr}\left([\boldsymbol{C}_{l_{1}}]^{-1}[\boldsymbol{f}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})][\boldsymbol{C}_{l_{2}}]^{-1}[\boldsymbol{f}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})]\right).$
(43)

The noise of the minimum-variance estimator is then simply $N(L)=N_{\rm
GMV}(L)=\lambda(L)$.

Putting everything together, the GMV estimator takes the final form

$\displaystyle\hat{\phi}_{\rm GMV}(\boldsymbol{L})$ $\displaystyle=$
$\displaystyle\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\sum_{XY}F_{XY}^{\rm
GMV}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})X(\boldsymbol{l}_{1})Y(\boldsymbol{l}_{2}),$
(44)

where the sum runs over the five unique pairs $XY=TT,TE,EE,TB,EB$. Explicitly,
the weights are

$\displaystyle F_{TT}^{\rm
GMV}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=\frac{N_{\rm
GMV}(L)}{2D_{l_{1}}D_{l_{2}}}\times$
$\displaystyle~{}~{}~{}\Big{[}C_{l_{1}}^{EE}C_{l_{2}}^{EE}f_{TT}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})+C_{l_{1}}^{TE}C_{l_{2}}^{TE}f_{EE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$
$\displaystyle~{}~{}~{}~{}~{}-C_{l_{1}}^{EE}C_{l_{2}}^{TE}f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})-C_{l_{2}}^{EE}C_{l_{1}}^{TE}f_{TE}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})\Big{]},~{}~{}~{}~{}~{}$
(45) $\displaystyle F_{EE}^{\rm
GMV}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=\frac{N_{\rm
GMV}(L)}{2D_{l_{1}}D_{l_{2}}}\times$
$\displaystyle~{}~{}~{}\Big{[}C_{l_{1}}^{TE}C_{l_{2}}^{TE}f_{TT}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})+C_{l_{1}}^{TT}C_{l_{2}}^{TT}f_{EE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$
$\displaystyle~{}~{}~{}~{}~{}-C_{l_{1}}^{TE}C_{l_{2}}^{TT}f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})-C_{l_{2}}^{TE}C_{l_{1}}^{TT}f_{TE}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})\Big{]},$
(46) $\displaystyle F_{TE}^{\rm
GMV}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=\frac{N_{\rm
GMV}(L)}{D_{l_{1}}D_{l_{2}}}\times$
$\displaystyle~{}~{}~{}\Big{[}-C_{l_{1}}^{TE}C_{l_{2}}^{EE}f_{TT}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})-C_{l_{1}}^{TE}C_{l_{2}}^{TT}f_{EE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$
$\displaystyle~{}~{}~{}~{}~{}+C_{l_{1}}^{EE}C_{l_{2}}^{TT}f_{TE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})+C_{l_{1}}^{TE}C_{l_{2}}^{TE}f_{TE}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})\Big{]},$
(47) $\displaystyle F_{TB}^{\rm
GMV}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=\frac{N_{\rm
GMV}(L)}{D_{l_{1}}C_{l_{2}}^{BB}}\times$
$\displaystyle~{}~{}~{}\Big{[}C_{l_{1}}^{EE}f_{TB}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})-C_{l_{1}}^{TE}f_{EB}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\Big{]},$
(48) $\displaystyle F_{EB}^{\rm
GMV}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})=\frac{N_{\rm
GMV}(L)}{D_{l_{1}}C_{l_{2}}^{BB}}\times$
$\displaystyle~{}~{}~{}\Big{[}-C_{l_{1}}^{TE}f_{TB}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})+C_{l_{1}}^{TT}f_{EB}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\Big{]}\,,$
(49)

where $D_{l}\equiv C^{TT}_{l}C^{EE}_{l}-[C^{TE}_{l}]^{2}$.

These explicit expressions should make it very clear that the GMV estimator is
different from the HO02 estimator. Put differently, the first step of the
likelihood-based iterative technique of Ref. Hirata and Seljak (2003a) is
_not_ equivalent to the HO02 estimator. Indeed, in the HO02 estimator, the
weight $F_{XY}$ of each pair $XY$ is proportional to $f_{XY}$ only (times a
function of $L$), even after combining all single-pair estimators; in
contrast, for the GMV estimator, the weight of each pair is a linear
combination of the response coefficients from _all_ pairs. The weights in the
GMV estimator would not separately minimize the variance of an individual $XY$
pair, but they provide a global optimum when combining all the pairs together.

These expressions also show that the weights are all sums of products of
function of $\boldsymbol{l}_{1}$ with functions of $\boldsymbol{l}_{2}$,
including $F_{TE}^{\rm GMV}$. In other words, the GMV estimator is separable
without requiring any additional approximation, making it well adapted for
efficient computations, as we now discuss.

As a side note, let us point out that the GMV estimator (just like the HO02 an
OH03 estimators) can be split into two pieces, built from $\\{TT,TE,EE\\}$ and
$\\{TB,EB\\}$, respectively, which are uncorrelated as
$C^{TB}_{\ell}=C^{EB}_{\ell}=0$. Our publicly available Python code
GlobalLensQuest first computes these two separate uncorrelated estimators
$\hat{\phi}_{a}\equiv\hat{\phi}_{\\{TT,TE,EE\\}}$ and
$\hat{\phi}_{b}\equiv\hat{\phi}_{\\{TB,EB\\}}$, and then combines them with
inverse variance weighting to obtain the GMV estimator.

###  IV.2 Compact configuration-space expression

We may write the GMV estimator in even more compact form by defining the
inverse-covariance-weighted fields

$\overline{\boldsymbol{X}}(\boldsymbol{l})\equiv[\boldsymbol{C}_{l}]^{-1}\boldsymbol{X}(\boldsymbol{l}),$
(50)

and write

$\overline{T}(\boldsymbol{l}_{1})\equiv\overline{X}^{1}(\boldsymbol{l}_{1})=[C_{l_{1}}^{EE}T(\boldsymbol{l}_{1})-C_{l_{1}}^{TE}E(\boldsymbol{l}_{1})]/D_{l_{1}},$
(51)

and similarly
$\overline{E}(\boldsymbol{l}_{1})=\overline{X}^{2},\overline{B}(\boldsymbol{l}_{1})=\overline{X}^{3}$.
We then have

$\hat{\phi}_{\rm
GMV}(L)=\frac{\lambda(L)}{2}\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\overline{X}^{i}(\boldsymbol{l})f_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\overline{X}^{j}(\boldsymbol{l}_{2}).$
(52)

We moreover define the Wiener-filtered fields

$X^{i}_{\rm
WF}(\boldsymbol{l})\equiv\widetilde{C}_{l}^{ij}~{}\overline{X}^{j},$ (53)

where $\widetilde{C}_{l}^{11}\equiv\widetilde{C}_{l}^{T\nabla
T},\widetilde{C}_{l}^{12}=\widetilde{C}_{l}^{21}\equiv\widetilde{C}_{l}^{T\nabla
E}$, etc., and write

$T_{\rm WF}(\boldsymbol{l})\equiv X^{1}_{\rm
WF}(\boldsymbol{l})\equiv\widetilde{C}_{l}^{T\nabla
T~{}}\overline{T}(\boldsymbol{l})+\widetilde{C}_{l}^{T\nabla
E}~{}\overline{E}(\boldsymbol{l}),$ (54)

and similarly for $E_{\rm WF}(\boldsymbol{l})\equiv X^{2}_{\rm
WF}(\boldsymbol{l})$.

In terms of these fields, the GMV estimator takes the particularly simple form

$\displaystyle\hat{\phi}_{\rm
GMV}(\boldsymbol{L})=\frac{\lambda(L)}{2}\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}(\boldsymbol{L}\cdot\boldsymbol{l}_{1})\Big{[}2T_{\rm
WF}(\boldsymbol{l}_{1})\overline{T}(\boldsymbol{l}_{2})$
$\displaystyle+\sum_{\epsilon=\pm 1}E_{\rm
WF}(\boldsymbol{l}_{1})\textrm{e}^{-2i\epsilon\varphi_{\boldsymbol{l}_{1}}}\left(\overline{E}(\boldsymbol{l}_{2})+i\epsilon\overline{B}(\boldsymbol{l}_{2})\right)\textrm{e}^{2i\epsilon\varphi_{\boldsymbol{l}_{2}}}\Big{]}.~{}$
(55)

This form is well adapted for efficient evaluation, as it is the sum of
convolutions of functions of $\boldsymbol{l}_{1}$ with functions of
$\boldsymbol{l}_{2}$. To see this, let us define

$\displaystyle~{}_{\pm 2}E_{\rm WF}(\boldsymbol{l})$ $\displaystyle\equiv$
$\displaystyle E_{\rm WF}(\boldsymbol{l})\textrm{e}^{\pm
2i\varphi_{\boldsymbol{l}}},$ (56) $\displaystyle~{}_{\pm
2}\overline{P}(\boldsymbol{l})$ $\displaystyle\equiv$
$\displaystyle{\frac{1}{2}}\left(\overline{E}(\boldsymbol{l})\pm
i\overline{B}(\boldsymbol{l})\right)\textrm{e}^{\pm
2i\varphi_{\boldsymbol{l}}}.$ (57)

We may then express the GMV estimator in terms of the configuration-space
versions of these fields (i.e. their inverse-Fourier transforms):

$\hat{\phi}_{\rm
GMV}(\boldsymbol{\hat{n}})=-\boldsymbol{\nabla}\cdot\mathcal{F}^{-1}\left[\lambda(L)\mathcal{F}[\boldsymbol{\psi}_{\rm
GMV}(\boldsymbol{\hat{n}})]\right],\\\ $ (58)

where

$\displaystyle\boldsymbol{\psi}_{\rm GMV}(\boldsymbol{\hat{n}})$
$\displaystyle=$ $\displaystyle\boldsymbol{\nabla}T_{\rm
WF}(\boldsymbol{\hat{n}})~{}\overline{T}(\boldsymbol{\hat{n}})$ (59)
$\displaystyle+\sum_{s=\pm 2}\boldsymbol{\nabla}[_{-s}E_{\rm
WF}(\boldsymbol{\hat{n}})]~{}_{s}\overline{P}(\boldsymbol{\hat{n}}).$

This expression is the flat-sky equivalent of Eq. (3) in Ref. Aghanim  _et
al._ (2020), derived in Ref. Carron (2019). A similar expression is derived in
Ref. Peloton  _et al._ (2017), in terms of $(T,Q,U)$ rather than $(T,E,B)$;
nevertheless, it is also incorrectly stated in that paper that this estimator
is identical to an estimator built out of single-pair estimators, i.e. the
HO02 estimator.

##  V Suboptimal quadratic estimator (SQE) used in recent data analyses

While the full expression for the configuration-space GMV estimator was
already known (although it was not known that it differs from the HO02
estimator) Aghanim  _et al._ (2020); Carron (2019), in practice only an
approximate version was used for data analyses thus far. Instead of using the
full covariance matrix $[\boldsymbol{C}_{l}]$ in Eq. (50), the Planck
collaboration Ade  _et al._ (2016a); Aghanim  _et al._ (2020) and SPT
collaboration Wu  _et al._ (2019) approximate it as diagonal by setting
$C_{l}^{TE}=0$ – note that Aghanim  _et al._ (2020) still use the exact
response coefficients $f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$. This
simplification allows to deal with a cut-sky setup with a lower computational
cost; it moreover preserves the configuration space separability. We denote
this suboptimal quadratic estimator SQE. Explicitly, the weights of this
estimator are

$\displaystyle F_{XY}^{\rm SQE}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$
$\displaystyle=$ $\displaystyle\lambda_{\rm
SQE}(L)\frac{f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})}{(1+\delta_{XY})C_{l_{1}}^{XX}C_{l_{2}}^{YY}},$
(60) $\displaystyle\lambda_{\rm SQE}(L)$ $\displaystyle\equiv$
$\displaystyle\left(\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\sum_{XY}\frac{f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})^{2}}{(1+\delta_{XY})C_{l_{1}}^{XX}C_{l_{2}}^{YY}}\right)^{-1},~{}~{}$
(61)

where again the sum runs of the five distinct pairs $XY$.

By definition, this estimator is suboptimal relative to the GMV estimator.
Furthermore, it should be clear from Eq. (36) that it is also noisier than the
OH03 estimator (and as we will see, also noisier than the HO02 estimator).
Indeed, the OH03 estimator accounts for the covariances between different
single-pair estimators, which depend on $C_{l}^{TE}$; as a consequence, the
$L$-dependent proportionality constant in Eq. (36) is different for each pair.
In contrast, the SQE amounts to neglecting correlations between single-pair
estimators, and simply using their inverse-variance combination, leading to
the same coefficient $\lambda_{\rm SQE}(L)$ for all weights in Eq. (60). Given
that the SQE estimator is effectively a linear combination of single-pair
estimator, and that the OH03 weights represent the optimal linear combination
of single-pair estimators, we conclude that the SQE estimator must be
suboptimal to the OH03 estimator.

To compute the noise of this estimator, we must account for $C_{l}^{TE}\neq
0$. We find

$\displaystyle N_{\rm SQE}(L)=\lambda_{\rm
SQE}^{2}(L)\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\sum_{XY}\frac{f_{XY}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})}{(1+\delta_{XY})C_{l_{1}}^{XX}C_{l_{2}}^{YY}}$
$\displaystyle\sum_{UV}\left[\frac{f_{UV}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})C_{l_{1}}^{XU}C_{l_{2}}^{YV}}{(1+\delta_{UV})C_{l_{1}}^{UU}C_{l_{2}}^{VV}}+\frac{f_{UV}(\boldsymbol{l}_{2},\boldsymbol{l}_{1})C_{l_{1}}^{XV}C_{l_{2}}^{YU}}{(1+\delta_{UV})C_{l_{2}}^{UU}C_{l_{1}}^{VV}}\right].~{}~{}$
(62)

The SQE estimator enables faster evaluation with a cut-sky, at the cost of
only $\sim 3\%$ increase in the reconstruction noise for Planck (Aghanim  _et
al._ , 2020), as we confirm in Fig. 2. However, we will see that for more
sensitive experimental setups, the reconstruction noise penalty can be more
than $10\%$ and thus a full joint filtering analysis of the temperature and
polarization maps would be beneficial in future.

##  VI Quantitative comparison of different quadratic estimators

###  VI.1 Experimental setups

In order to evaluate the variance of different quadratic estimators, we use
three different setups which correspond to the Planck, SO-like, and CMBS4-like
experiments. In Table 2, we provide the adopted specifications for these
setups. The Gaussian random noise of the detector is calculated as

$C^{T,E/B}_{\ell}|_{\mathrm{noise}}=(\Delta_{T,P})^{2}e^{\ell(\ell+1)\sigma^{2}/8\ln{2}}$
(63)

where $\Delta_{T,P}$ denote the white-noise of the detector in $\mu$K-radian,
and $\sigma$ is the full width at half maximum (FWHM) of the beam in arcmin.

Experiment | $\ell_{\mathrm{max}}$ | $\Delta_{T}$ | $\Delta_{P}$ |  $\sigma$ [
---|---|---|---|---
|  |  $\mu\mathrm{K}$-arcmin |  $\mu\mathrm{K}$-arcmin | arcmin [
Planck | 3000 | 35.0 | 60.0 | 5.0 [
SO | 3000 | 8.0 | 8.0$\sqrt{2}$ | 1.4 [
CMBS4 | 3000 | 1.0 | 1.0$\sqrt{2}$ | 1.0 [
Table 2: Experimental specifications used in this work.

It has been shown that extra-galactic foregrounds can bias the CMB lensing
reconstruction from temperature maps and different strategies have been
proposed to mitigate this issue (e.g. Ferraro and Hill, 2018; Schaan and
Ferraro, 2019; Madhavacheril and Hill, 2018; Darwish  _et al._ , 2021; van
Engelen  _et al._ , 2014). We neglect the foregrounds in this study, and
simply choose $\ell_{\mathrm{max}}=3000$ in both temperature and polarization.
Although it is possible to go for a much higher $\ell_{\mathrm{max}}$ in
polarization than in temperature due to lack of strongly polarized
foregrounds, for simplicity we take
$\ell^{T}_{\mathrm{max}}=\ell^{P}_{\mathrm{max}}$. We have checked our results
with different $\ell_{\mathrm{max}}$ ranges and found no drastic difference in
our results.

###  VI.2 Comparison of reconstruction noises

We start by comparing the OH03 and HO02 quadratic estimators in Figure 1. We
find that the HO02 estimator is systematically less noisy than the OH03
estimator for large angular scales (by less than 0.5%). Interestingly, the
OH03 estimator becomes slightly less noisy than the HO02 estimator for $L$
larger than several hundreds. We have checked that our numerical integrals are
converged to better than 0.01% relative accuracy up to $L\approx 2000$, and to
better than 0.03% for $L\lesssim 3000$ which we also show in Fig. 4. This
gives us confidence that the $\sim 0.08\%$ improvement of OH03 estimator over
HO02 seen for a CMBS4-like experiment is real and not a numerical artifact.

The lower noise of OH03 at small angular scales may appear surprising at
first, given that this estimator uses a TE estimator suboptimal to that of
HO02. However, the suboptimal TE estimator does not guarantee that the overall
OH03 estimator (obtained by optimally combining the five single-pair
estimators) is noisier than the overall HO02 estimator: indeed, if the OH03 TE
estimator happens to be more correlated with the TT and EE estimators than the
HO02 TE estimator which is the case here, the overall combination of OH03
estimators can be less noisy than that of HO02.

Figure 1: Fractional difference between the variance of the OH03 and HO02
estimators. Different colors correspond to the different experimental setups
described in Sec. VI.1.

In Fig. 2, we show the ratios of the reconstruction noise of the GMV and HO02
estimators to that of the SQE estimator for different experimental setups. As
expected, we find that the variance of the GMV quadratic estimator is lower
than that of the HO02 and SQE estimators. Also, the variance of HO02 estimator
is smaller than the SQE estimator. For a Planck-like experimental setup, on
large angular scales ($L\lesssim 500$), the difference between the SQE and GMV
estimators is of the order of 3% and can reach $\sim 6\%$ around $L\sim 2000$.
However, for more sensitive experiments, this difference reaches $\sim 11\%$
and $\sim 12\%$ at $L\lesssim 100$ and $L\sim 2000$ for SO- or CMBS4-like
experiments, respectively. This result may motivate using the full covariance
matrix $[\boldsymbol{C}_{l}]$ in Eq. (50) instead of assuming $C_{l}^{TE}=0$,
in order to obtain more precise results in future data analyses.

Figure 2: Ratio of the minimum variance reconstruction noise of the GMV and
SQE estimators, and HO02 and SQE estimators for different experimental setups.

From Fig. 2, we can also see that for a Planck-like experimental setup, at
small angular scales ($L\gtrsim 1000$) the GMV and HO02 estimators perform
almost equally well, while for SO- and CMBS4-like experiment GMV outperforms
HO02 everywhere. On large angular scales, the difference between the GMV and
HO02 is much more significant for all the experiments considered here. For
Planck, polarization noise is significant, so the difference in lensing
estimators is quite small on all scales. However, although CMBS4-like
experiments are $EB$-dominated for the purpose of lensing reconstruction, the
improvement of the GMV over HO02 is driven by a significant difference in the
TT-TE-EE part of the minimum-variance estimator (rather than improved
filtering of $E$ improving $EB$; see Appendix C). The effect on the combined
MV estimator is largest at $L\lesssim 200$ where the MV estimator has the
largest contribution from TT-TE-EE.

##  VII Conclusions

Quadratic estimators (QEs) are widely used to reconstruct the CMB lensing
potential from CMB temperature and polarization maps. In fact, up until the
very recent POLARBEAR Adachi  _et al._ (2020) and SPTPol Millea  _et al._
(2020) results, all maps of the CMB lensing potential have been constructed
using QEs. In this work, we present a clear comparison between different QEs,
both in terms of explicit equations, and quantitatively, by comparing their
reconstruction noise. Importantly, we show that the Hu-Okamoto Hu and Okamoto
(2002) (HO02) optimization method, consisting in first constructing optimal
single-pair quadratic estimators, and then their optimal linear combination,
does not lead to the absolute minimum-variance QE. Instead, we derive the
global-minimum-variance (GMV) QE, which minimizes the variance of quadratic
temperature and polarization combinations in one single step.

Interestingly, the GMV estimator derived here had been hiding in plain sight
in previous works. It is equivalent to the first step (or weak-signal limit)
of likelihood-based methods Hirata and Seljak (2003a); Peloton  _et al._
(2017); Carron and Lewis (2017), which is therefore _not_ equivalent to the
HO02 estimator, contrary to what was previously thought (although technically
the GMV estimator with non-perturbative lensed gradient weights presented here
is a modification to the first-step of likelihood based estimators so that the
result is non-perturbatively unbiased). Our work is the first to note that the
HO02 estimator is not the global-optimum quadratic estimator, and make this
point sharply clear through explicit expressions, as well as numerical
comparisons. Indeed, we show that the reconstruction noise of the GMV
estimator is lower than that of the HO02 estimator by up to $\sim 9\%$ for a
SO-like experiment.

We also study the suboptimal QE used in the 2018 Planck Ade  _et al._ (2016a);
Carron and Lewis (2017); Aghanim  _et al._ (2020) and recent SPT Wu  _et al._
(2019) lensing analyses (SQE), which is obtained from the GMV quadratic
estimator (appropriately generalized to account for beam and pixel
convolution), with the additional approximation of neglecting $C_{l}^{TE}=0$
in the inverse filter matrix. We show that this approximation makes the SQE
suboptimal not only relative to the GMV estimator, but also relative to the
HO02 and OH03 estimators. We evaluate the reconstruction noise of the
different estimators for ongoing and planned CMB experimental setups and find
that while the improvement in the reconstruction noise between the SQE and
HO02 estimator is of order $\sim 1-8\%$, the difference between the SQE and
GMV estimators is $\sim 9-12\%$ for more sensitive experiments, especially on
large angular scales $L\lesssim 10^{2}$ and scales around $L\sim 2000$. This
improvement amounts to achieving a better sensitivity for the same experiment
at no additional cost. This should motivate overcoming the added complexity
associated with joint filtering of cut-sky temperature and polarization maps,
in order to be able to use the GMV estimator in future lensing data analyses.
Our Python code
GlobalLensQuest111https://github.com/abhimaniyar/GlobalLensQuest to compare
and compute the noise variances of the HO02, OH03, GMV, and SQE estimators and
Julien Carron’s codes LensIt222https://github.com/carronj/LensIt and
plancklens333https://github.com/carronj/plancklens which can perform the
optimal GMV operation with anisotropic noise and cut-sky are publicly
available.

While in this work we have chosen to present relevant equations in the flat-
sky limit for conciseness, it is straightforward to generalize our results to
the full-sky case. The flat-sky fields $X(\boldsymbol{l})$ are to be replaced
by the full-sky harmonic coefficients $X_{\ell m}$; the generalization of Eq.
(2) and of the coupling coefficients
$f_{XY}(\boldsymbol{l},\boldsymbol{l}^{\prime})$ are then provided in Ref.
Okamoto and Hu (2003). While Okamoto and Hu (2003) provide the full-sky
expressions for the HO02 estimator, Carron (2019) provide the full-sky
expressions for the GMV estimator incorporating the instrumental beam response
and anisotropic noise. We expect a comparable improvement over the full-sky
version of the HO02 estimator Okamoto and Hu (2003) when using the GMV
estimator Carron (2019).

The approach presented here for the GMV estimator would also apply to any
other joint estimator constructed from similar linear combinations of other
estimators. For example, the foreground-immune hybrid QE of Ref. Schaan and
Ferraro (2019) splits the $TT$ lensing estimator into magnification-only and
shear-only estimators, and then form a hybrid estimator through a minimum-
variance linear combination of these two estimators. Their hybrid estimator
can be further optimized by following the logic presented here, i.e. searching
for the global-minimum-variance shear and magnification estimator, accounting
for correlations for each $(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$, rather
than after integration over $\boldsymbol{l}_{1},\boldsymbol{l}_{2}$.

For SO- and CMBS4-like experiments, on large angular scales ($L<100$) the
reconstruction is expected to be signal dominated Ade  _et al._ (2019);
Abazajian  _et al._ (2019). The power spectrum $C^{\phi\phi}_{L}$ uncertainty
is therefore dominated by cosmic variance and using the GMV estimator rather
than the HO02 or SQE estimators would not drastically affect the measurement
of $C^{\phi\phi}_{L}$ on these large angular scales. However, the reduction in
the noise of the reconstructed $\phi$ field on the signal-dominated large
angular scales will be beneficial for science goals which involve cross-
correlation of the $\phi$ field with other tracers of large-scale structure
Ade  _et al._ (2019); Schmittfull and Seljak (2018), utilizing the sample
variance cancellation through cross-correlations. Also, lensing induced
B-modes act as a source of noise and limit the measurement of the primordial
B-modes Zaldarriaga and Seljak (1998) which is a major scientific goal for CMB
experiments. These modes can be removed using map-level estimates of both the
primordial E-modes and lensing potential $\phi$ with a technique called
delensing and depend on the estimate of the particular realization of $\phi$
in the given patch of the sky Seljak and Hirata (2004); Ade  _et al._ (2019);
Abazajian  _et al._ (2019). Lower noise estimates of the $\phi$ field will
therefore be crucial for such operations and motivate using the GMV estimator
instead of other QEs.

Even if future CMBS4 experiments will likely use likelihood-based iterative
methods to reconstruct the lensing potential, QEs will remain very useful as a
forecasting and cross-checking tool. More immediately, QEs are still the
primary lensing reconstruction tool for current and near-future CMB
experiments. The GMV estimator will therefore be a useful tool to harvest even
more information out of the CMB data.

## Acknowledgements

We acknowledge useful discussions with Colin Hill, Emmanuel Schaan, and Simone
Ferraro. YAH is supported by NSF award No 1820861 and NASA grant No
80NSSC20K0532. JC acknowledges support from a SNSF Eccellenza Professorial
Fellowship (No. 186879). AL has support from the UK STFC grant ST/T000473/1.

##  Appendix A Explicit form for the GMV $N^{(1)}(L)$

As pointed out in Sec. II, in the main text we only optimize relative to the
reconstruction noise $N^{(0)}(L)$. Here we give an explicit expression for the
additional $N^{(1)}(L)$ bias Kesden  _et al._ (2003) for the GMV estimator,
for which an explicit expression has not been provided in the literature.

The covariance of the GMV estimator given by Eq. (37) becomes

$\displaystyle\langle\hat{\phi}(\boldsymbol{L})\hat{\phi}(\boldsymbol{L}^{\prime})\rangle=\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\int_{\begin{subarray}{c}\boldsymbol{l}^{\prime}_{1}+\boldsymbol{l}^{\prime}_{2}\\\
=\boldsymbol{L}^{\prime}\end{subarray}}$ $\displaystyle\langle
X^{i}(\boldsymbol{l}_{1})X^{j}(\boldsymbol{l}_{2})X^{p}(\boldsymbol{l}^{\prime}_{1})X^{q}(\boldsymbol{l}^{\prime}_{2})\rangle\Xi_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\Xi_{pq}(\boldsymbol{l}^{\prime}_{1},\boldsymbol{l}^{\prime}_{2})\,.$
(64)

Evaluating the expectation value in the brackets

$\displaystyle\langle
X^{i}(\boldsymbol{l}_{1})X^{j}(\boldsymbol{l}_{2})X^{p}(\boldsymbol{l}^{\prime}_{1})X^{q}(\boldsymbol{l}^{\prime}_{2})\rangle$
$\displaystyle=$
$\displaystyle(2\pi)^{4}\big{[}C^{ij}_{l_{1}}C^{pq}_{l^{\prime}_{1}}\delta_{\rm
D}(\boldsymbol{L})\delta_{\rm
D}(\boldsymbol{L}^{\prime})+C^{ip}_{l_{1}}C^{jq}_{l^{\prime}_{1}}\delta_{\rm
D}(\boldsymbol{l}_{1}+\boldsymbol{l}^{\prime}_{1})\delta_{\rm
D}(\boldsymbol{l}_{2}+\boldsymbol{l}^{\prime}_{2})$ (65)
$\displaystyle+C^{iq}_{l_{1}}C^{jp}_{l^{\prime}_{1}}\delta_{\rm
D}(\boldsymbol{l}_{1}+\boldsymbol{l}^{\prime}_{2})\delta_{\rm
D}(\boldsymbol{l}_{2}+\boldsymbol{l}^{\prime}_{1})\big{]}+(2\pi)^{2}T^{ijpq}(\boldsymbol{l}_{1},\boldsymbol{l}_{2},\boldsymbol{l}^{\prime}_{1},\boldsymbol{l}^{\prime}_{2})\delta_{\rm
D}(\boldsymbol{L}+\boldsymbol{L}^{\prime})\,,$

where the first term in the square brackets disappears because
$\boldsymbol{L}\neq 0$ and rest of the terms in the square bracket represent
$N^{(0)}(L)$ and
$T^{ijpq}(\boldsymbol{l}_{1},\boldsymbol{l}_{2},\boldsymbol{l}^{\prime}_{1},\boldsymbol{l}^{\prime}_{2})$,
the trispectrum containing terms that contribute to the lensing power spectrum
signal and the signal-dependent $N^{(1)}(L)$ bias. Following Kesden  _et al._
(2003), the trispectrum term can be written in terms of
$f_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$ to first order in explicit
$C^{\phi\phi}_{L}$ as

$\displaystyle
T^{ijpq}(\boldsymbol{l}_{1},\boldsymbol{l}_{2},\boldsymbol{l}^{\prime}_{1},\boldsymbol{l}^{\prime}_{2})=C^{\phi\phi}_{|\boldsymbol{l}_{1}+\boldsymbol{l}_{2}|}f_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})f_{pq}(\boldsymbol{l}^{\prime}_{1},\boldsymbol{l}^{\prime}_{2})+C^{\phi\phi}_{|\boldsymbol{l}_{1}+\boldsymbol{l}^{\prime}_{1}|}f_{ip}(\boldsymbol{l}_{1},\boldsymbol{l}^{\prime}_{1})f_{jq}(\boldsymbol{l}_{2},\boldsymbol{l}^{\prime}_{2})+C^{\phi\phi}_{|\boldsymbol{l}_{1}+\boldsymbol{l}^{\prime}_{2}|}f_{iq}(\boldsymbol{l}_{1},\boldsymbol{l}^{\prime}_{2})f_{jp}(\boldsymbol{l}_{2},\boldsymbol{l}^{\prime}_{1}).~{}~{}~{}~{}$
(66)

Substituting Eqs. (65) and (66) in Eq. (64) we have

$\langle\hat{\phi}(\boldsymbol{L})\hat{\phi}(\boldsymbol{L}^{\prime})\rangle=(2\pi)^{2}\delta_{\rm
D}(\boldsymbol{L}+\boldsymbol{L}^{\prime})[C^{\phi\phi}_{L}+N^{(0)}(L)+N^{(1)}(L)]\,,$
(67)

where

$\displaystyle N^{(1)}(L)$ $\displaystyle\equiv$
$\displaystyle\int_{\begin{subarray}{c}\boldsymbol{l}_{1}+\boldsymbol{l}_{2}\\\
=\boldsymbol{L}\end{subarray}}\int_{\begin{subarray}{c}\boldsymbol{l}^{\prime}_{1}+\boldsymbol{l}^{\prime}_{2}\\\
=\boldsymbol{L}^{\prime}\end{subarray}}\Xi_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})\Xi_{pq}(\boldsymbol{l}^{\prime}_{1},\boldsymbol{l}^{\prime}_{2})\times\big{[}C^{\phi\phi}_{|\boldsymbol{l}_{1}+\boldsymbol{l}^{\prime}_{1}|}f_{ip}(\boldsymbol{l}_{1},\boldsymbol{l}^{\prime}_{1})f_{jq}(\boldsymbol{l}_{2},\boldsymbol{l}^{\prime}_{2})+C^{\phi\phi}_{|\boldsymbol{l}_{1}+\boldsymbol{l}^{\prime}_{2}|}f_{iq}(\boldsymbol{l}_{1},\boldsymbol{l}^{\prime}_{2})f_{jp}(\boldsymbol{l}_{2},\boldsymbol{l}^{\prime}_{1})\big{]}\,,$
(68)

The optimal weight matrix $\Xi_{ij}(\boldsymbol{l}_{1},\boldsymbol{l}_{2})$
was determined in order to minimize the variance while only considering
$N^{(0)}(L)$.

Figure 3: Comparison of the $N^{(0)}(L)$ and $N^{(1)}(L)$ curves for the GMV
estimator (left) and ratio of the $N^{(1)}(L)$ for SQE and GMV estimators
(right) for given experimental configurations.

In the left panel of Fig. 3, we show the comparison between the $N^{(0)}(L)$
and $N^{(1)}(L)$ curves for the GMV estimator for our given experimental
configurations. We find that for Planck-like experiments $N^{(1)}(L)$ is a
couple of orders of magnitude smaller than $N^{(0)}(L)$, while for more
sensitive (less noisy) SO- and CMBS4-like experiments, $N^{(1)}(L)$ is a
factor of few to an order of magnitude smaller than $N^{(0)}(L)$. On small
scales it can however become comparable to the signal spectrum and is
important to model in any likelihood analysis. It would be straightforward to
apply the perturbative likelihood approximation of Ade  _et al._ (2016a);
Aghanim  _et al._ (2020) (which accounts consistently for the signal
dependence of $N^{(1)}(L)$) to the GMV estimator. In the right panel of Fig.
3, we show the ratio of the $N^{(1)}(L)$ for SQE and GMV estimators for
different experiments considered here. For SO- and CMBS4-like experiments,
$N^{(1)}(L)$ bias for the GMV estimator is smaller than for the SQE estimator
for $L\lesssim 1800$. For Planck-like experiments, apart from the large
angular scales $L\lesssim 200$ where GMV estimator gives smaller $N^{(1)}(L)$
bias than the SQE estimator, both the estimators have almost the same
$N^{(1)}(L)$ bias.

##  Appendix B Numerical convergence

In Fig. 4, we show the result of the convergence test we perform for our
Python code. It shows the % change in the noise calculation of HO02 (dashed
curves) and OH03 (solid curves) estimators when we double the number of steps
in the angular part of the integration for the given $\ell_{\rm max}$ and
$L\lesssim 3000$. The dashed HO02 curves mostly overlap with the solid OH03
curves and thus are not distinctly visible. This shows that our numerical
integrals are converged to better than 0.01% relative accuracy up to $L\approx
2000$, and to better than 0.03% for $L\lesssim 3000$, and thus makes us
confident that the improvement observed for OH03 over HO02 for small angular
scales is not a numerical artifact.

Figure 4: % change in the noise curves for HO02 and OH03 estimators when we
double the number of integration steps in our Python code.

##  Appendix C GMV to HO02 comparison

We perform the following exercise to compare the GMV and HO02 estimators. As
mentioned at the end of Sec. IV.1, the GMV estimator can be split into two
independent estimators $\hat{\phi}^{\rm GMV}_{\\{TT,TE,EE\\}}$ and
$\hat{\phi}^{\rm GMV}_{\\{TB,EB\\}}$. We compare the minimum variance
reconstruction noise of these individual estimators with their HO02
counterparts i.e. $\hat{\phi}^{\rm HO02}_{\\{TT,TE,EE\\}}$ and
$\hat{\phi}^{\rm HO02}_{\\{TB,EB\\}}$. This is shown in Fig. 5, where we plot
the ratio of the reconstruction noises for the GMV and HO02 versions of these
two estimators. As we can see, HO02 performs almost equally well as the GMV
estimator for $\\{TB,EB\\}$ pair. The overall improvement of the GMV estimator
over HO02 estimator is thus mainly driven by $\\{TT,TE,EE\\}$, especially for
more sensitive SO- and CMBS4-like experiments. We also show an ideal case
setup in Fig. 5 which corresponds to a noise-less experiment with the same
multipole ranges as other experiments considered. HO02 estimator for a
CMBS4-like setup considered here performs almost as well as the GMV estimator
for $\\{TB,EB\\}$ pair and very slightly under-performs for the
$\\{TT,TE,EE\\}$ set when compared to the ideal case setup.

Figure 5: Ratio of the minimum variance reconstruction noise of the GMV and
HO02 estimators for $\hat{\phi}_{\\{TT,TE,EE\\}}$ and
$\hat{\phi}_{\\{TB,EB\\}}$ estimators for different experimental setups. The
ideal setup corresponds to a noise-less case.

## References

  * Blanchard and Schneider (1987) A. Blanchard and J. Schneider, Astron. Astrophys. 184, 1 (1987).
  * Lewis and Challinor (2006) A. Lewis and A. Challinor, Phys. Rep. 429, 1 (2006), arXiv:astro-ph/0601594 [astro-ph.CO] .
  * Seljak and Zaldarriaga (1999) U. Seljak and M. Zaldarriaga, Phys. Rev. Lett. 82, 2636 (1999), arXiv:astro-ph/9810092 [astro-ph] .
  * Allison  _et al._ (2015) R. Allison, P. Caucal, E. Calabrese, J. Dunkley, and T. Louis, Phys. Rev. D 92, 123535 (2015), arXiv:1509.07471 [astro-ph.CO] .
  * Schmittfull and Seljak (2018) M. Schmittfull and U. Seljak, Phys. Rev. D 97, 123540 (2018), arXiv:1710.09465 [astro-ph.CO] .
  * Das  _et al._ (2011) S. Das _et al._ , Phys. Rev. Lett. 107, 021301 (2011), arXiv:1103.2124 [astro-ph.CO] .
  * Sherwin  _et al._ (2017) B. D. Sherwin _et al._ , Phys. Rev. D 95, 123529 (2017), arXiv:1611.09753 [astro-ph.CO] .
  * van Engelen  _et al._ (2012) A. van Engelen _et al._ , Astrophys. J. 756, 142 (2012), arXiv:1202.0546 [astro-ph.CO] .
  * Ade  _et al._ (2016a) P. A. R. Ade _et al._ (Planck Collaboration), Astron. Astrophys. 594, A15 (2016a), arXiv:1502.01591 [astro-ph.CO] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck Collaboration), Astron. Astrophys. 641, A8 (2020), arXiv:1807.06210 [astro-ph.CO] .
  * Omori  _et al._ (2017) Y. Omori _et al._ , Astrophys. J. 849, 124 (2017), arXiv:1705.00743 [astro-ph.CO] .
  * Story  _et al._ (2015) K. T. Story _et al._ , Astrophys. J. 810, 50 (2015), arXiv:1412.4760 [astro-ph.CO] .
  * Wu  _et al._ (2019) W. L. K. Wu _et al._ , Astrophys. J. 884, 70 (2019), arXiv:1905.05777 [astro-ph.CO] .
  * Millea  _et al._ (2020) M. Millea _et al._ , (2020), arXiv:2012.01709 [astro-ph.CO] .
  * Ade  _et al._ (2016b) P. A. R. Ade _et al._ (BICEP2 Collaboration and Keck Array Collaboration), Astrophys. J. 833, 228 (2016b), arXiv:1606.01968 [astro-ph.CO] .
  * Ade  _et al._ (2014) P. A. R. Ade _et al._ (Polarbear Collaboration), Phys. Rev. Lett. 113, 021301 (2014), arXiv:1312.6646 [astro-ph.CO] .
  * Adachi  _et al._ (2020) S. Adachi _et al._ (Polarbear Collaboration), Phys. Rev. Lett. 124, 131301 (2020), arXiv:1909.13832 [astro-ph.CO] .
  * Henderson  _et al._ (2016) S. W. Henderson _et al._ , Journal of Low Temperature Physics 184, 772 (2016), arXiv:1510.02809 [astro-ph.IM] .
  * Benson  _et al._ (2014) B. A. Benson _et al._ , in  _Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VII_, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, Vol. 9153, edited by W. S. Holland and J. Zmuidzinas (2014) p. 91531P, arXiv:1407.2973 [astro-ph.IM] .
  * Ade  _et al._ (2019) P. Ade _et al._ (Simons Observatory Collaboration), JCAP 2019, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * Abazajian  _et al._ (2016) K. N. Abazajian _et al._ , arXiv e-prints (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Abazajian  _et al._ (2019) K. Abazajian _et al._ (CMB-S4 Collaboration), arXiv e-prints (2019), arXiv:1907.04473 [astro-ph.IM] .
  * Zaldarriaga and Seljak (1999) M. Zaldarriaga and U. Seljak, Phys. Rev. D 59, 123507 (1999), arXiv:astro-ph/9810257 [astro-ph] .
  * Hu (2001) W. Hu, Astrophys. J. Lett. 557, L79 (2001), arXiv:astro-ph/0105424 [astro-ph] .
  * Hirata and Seljak (2003a) C. M. Hirata and U. Seljak, Phys. Rev. D 68, 083002 (2003a), arXiv:astro-ph/0306354 [astro-ph] .
  * Hadzhiyska  _et al._ (2019) B. Hadzhiyska, B. D. Sherwin, M. Madhavacheril, and S. Ferraro, Phys. Rev. D 100, 023547 (2019), arXiv:1905.04217 [astro-ph.CO] .
  * Hirata and Seljak (2003b) C. M. Hirata and U. Seljak, Phys. Rev. D 67, 043001 (2003b), arXiv:astro-ph/0209489 [astro-ph] .
  * Carron and Lewis (2017) J. Carron and A. Lewis, Phys. Rev. D 96, 063510 (2017), arXiv:1704.08230 [astro-ph.CO] .
  * Millea  _et al._ (2020) M. Millea, E. Anderes, and B. D. Wandelt, Phys. Rev. D 102, 123542 (2020), arXiv:2002.00965 [astro-ph.CO] .
  * Hu and Okamoto (2002) W. Hu and T. Okamoto, Astrophys. J. 574, 566 (2002), arXiv:astro-ph/0111606 [astro-ph] .
  * Lewis  _et al._ (2011) A. Lewis, A. Challinor, and D. Hanson, JCAP 2011, 018 (2011), arXiv:1101.2234 [astro-ph.CO] .
  * Fabbian  _et al._ (2019) G. Fabbian, A. Lewis, and D. Beck, JCAP 2019, 057 (2019), arXiv:1906.08760 [astro-ph.CO] .
  * Okamoto and Hu (2003) T. Okamoto and W. Hu, Phys. Rev. D 67, 083002 (2003), arXiv:astro-ph/0301031 [astro-ph] .
  * Ferraro and Hill (2018) S. Ferraro and J. C. Hill, Phys. Rev. D 97, 023512 (2018), arXiv:1705.06751 [astro-ph.CO] .
  * Schaan and Ferraro (2019) E. Schaan and S. Ferraro, Phys. Rev. Lett. 122, 181301 (2019), arXiv:1804.06403 [astro-ph.CO] .
  * Kesden  _et al._ (2003) M. Kesden, A. Cooray, and M. Kamionkowski, Phys. Rev. D 67, 123507 (2003), arXiv:astro-ph/0302536 [astro-ph] .
  * Carron (2019) J. Carron, (2019), arXiv:1908.02016 [astro-ph.CO] .
  * Peloton  _et al._ (2017) J. Peloton, M. Schmittfull, A. Lewis, J. Carron, and O. Zahn, Phys. Rev. D 95, 043508 (2017), arXiv:1611.01446 [astro-ph.CO] .
  * Madhavacheril and Hill (2018) M. S. Madhavacheril and J. C. Hill, Phys. Rev. D 98, 023534 (2018), arXiv:1802.08230 [astro-ph.CO] .
  * Darwish  _et al._ (2021) O. Darwish _et al._ , Mon. Not. Roy. Ast. Soc. 500, 2250 (2021), arXiv:2004.01139 [astro-ph.CO] .
  * van Engelen  _et al._ (2014) A. van Engelen, S. Bhattacharya, N. Sehgal, G. P. Holder, O. Zahn, and D. Nagai, Astrophys. J. 786, 13 (2014), arXiv:1310.7023 [astro-ph.CO] .
  * Zaldarriaga and Seljak (1998) M. Zaldarriaga and U. Seljak, Phys. Rev. D 58, 023003 (1998), astro-ph/9803150 .
  * Seljak and Hirata (2004) U. Seljak and C. M. Hirata, Phys. Rev. D 69, 043005 (2004), arXiv:astro-ph/0310163 [astro-ph] .