CMB lensing with shear-only reconstruction on the full sky

Reconstruction of gravitational lensing effects in the CMB from current and upcoming surveys is still dominated by temperature anisotropies. Extragalactic foregrounds in temperature maps can induce significant biases in the lensing power spectrum obtained with the standard quadratic estimators. Techniques such as masking cannot remove these foregrounds fully, and the residuals can still lead to large biases if unaccounted for. In this paper, we study the"shear-only"estimator, an example of a class of geometric methods that suppress extragalactic foreground contamination while making only minimal assumptions about foreground properties. The shear-only estimator has only been formulated in the flat-sky limit and so is not easily applied to wide surveys. Here, we derive the full-sky version of the shear-only estimator and its generalisation to an $m=2$ multipole estimator that has improved performance for lensing reconstruction on smaller scales. The multipole estimator is generally not separable, and so is expensive to compute. We explore separable approximations based on a singular-value decomposition, which allow efficient evaluation of the estimator with real-space methods. Finally, we apply these estimators to simulations that include extragalactic foregrounds and verify their efficacy in suppressing foreground biases.


I. INTRODUCTION
Gravitational lensing of the CMB encodes a wealth of information about our Universe.Observing the deflections produced by the intervening large-scale structure on the paths of CMB photons allows us to make integrated measurements of the projected matter distribution to high redshifts [1].CMB lensing provides us with a powerful probe to constrain parameters such as neutrino masses [2,3] and dark energy [4].Analyses with Planck and AdvACT data, building on earlier work by the Atacama Cosmology Telescope (ACT) and South Pole Telescope [5,6], have demonstrated the great potential of this approach; see [7][8][9] for the most recent results.Current and future surveys, such as AdvACT [10], SPT-3G [11], Simons Observatory [12] and CMB-S4 [13], will improve the precision of CMB lensing measurements significantly, making the identification and reduction of systematic biases increasingly important.While one expects polarisation information to dominate in the reconstruction of lensing from future surveys, many current and upcoming CMB surveys will still rely heavily on temperature.In this regime, extragalactic foreground contamination from the cosmic infrared background (CIB), the thermal Sunyaev-Zel'dovich effect (tSZ), the kinematic Sunyaev-Zel'dovich effect (kSZ) and radio point sources (PS) can leak into the lensing estimator producing significant biases if unaccounted for [14].Several mitigation methods for these biases have been proposed.For example, masking out sources from a known catalogue can decrease this bias, and techniques such as bias hardening [15,16], which involves reconstructing and projecting out foregrounds, are useful for cases when the statistical properties of the foregrounds are known.Another method is multi-frequency component separation [17], which can reduce or null specific foregrounds, but it was found that simultaneously reducing the CIB and tSZ increases the noise by a large factor.An improved technique, building upon [17], was introduced in [18] to eliminate foregrounds from the tSZ while preserving most of the signal-to-noise.Finally, [19,20] explore which combinations of multi-frequency cleaning and geometric methods (bias hardening and shear) are most effective in controlling lensing biases with only modest reduction in signal-to-noise.
In this paper, we focus on the shear estimator introduced in [21], which built on earlier work exploring the role of magnification and shear in CMB lensing reconstruction [22,23].The idea behind the shear estimator is to exploit the different geometric effects on the local CMB two-point function of lensing and extragalactic foregrounds to separate them.In the limit where large-scale lenses are reconstructed from small-scale temperature anisotropies, weak lensing produces local distortions in the 2D CMB power spectrum with an isotropic part (i.e., monopole) due to lensing convergence and a quadrupolar part due to lensing shear.
The quadratic estimators usually employed in lensing reconstruction [24,25] optimally combine convergence and shear in the large-scale-lens limit.However, extragalactic foreground power predominantly biases the local monopole power spectrum, leaving the shear-only estimator much less affected by foregrounds than the standard quadratic estimator.Moreover, with the shear-only estimator, one can include smaller-scale temperature modes in the reconstruction without introducing unacceptable levels of bias, thus mitigating the loss of signal-to-noise from discarding convergence information in the case of high-resolution, low-noise observations [21].
For the reconstruction of smaller-scale lenses, it is no longer true that the lensing convergence and shear can be considered constant over the coherence scale of the CMB.In this limit, lensing not only introduces monopole (m = 0) and quadrupole (m = 2) couplings in the local two-point function of the lensed CMB but also higher-order couplings.Furthermore, the dependencies of the m = 0 and m = 2 couplings on the angular scale of the CMB fluctuations and lenses deviate significantly from their large-lens limits.For reconstruction on smaller scales, one can formulate a set of multipole estimators, each extracting information from a specific m [21].Most of the reconstruction signal-to-noise is still contained in the m = 0 and m = 2 estimators, and extragalactic foregrounds are expected still to bias mainly m = 0.However, it is necessary to use the correct scale dependence of the m = 2 estimator to avoid the poor performance of the shear estimator when it is extended directly to reconstruction of smaller-scale lenses.This makes efficient evaluation with real-space methods difficult since the m = 2 estimator is not generally separable.
The shear reconstruction discussed in [21] is based on the flat-sky approximation.Since current and future highresolution CMB experiments will cover a significant fraction of the sky, a full-sky formulation of the shear estimator is required.In this paper, we derive the full-sky version of the shear estimator and show how it can be evaluated efficiently in real space with spin-weighted spherical harmonic transforms.We generalise further to an m = 2 multipole estimator to avoid the sub-optimal performance of the shear estimator on smaller scales.We suggest a simple separable approximation for this estimator using singular-value decomposition, which allows efficient evaluation with real-space methods.[Although extending the shear formalism to polarization-based estimators is straightforward, as shown very recently in the flat-sky limit by [26], our focus in this paper is solely on the shear formalism for the temperature part of the lensing estimator.This is motivated by the fact that the dominant sources of extragalactic foregrounds from the CIB or tSZ are only weakly polarized compared to the CMB fluctuations on the scales that dominate lensing reconstruction, meaning that the biases they induce are relatively smaller in polarization-based reconstructions than for temperature.Moreover, temperature-based estimators still make a significant contribution to lensing reconstruction at the sensitivity of current surveys.While lensing with future surveys, such as CMB-S4 [13], will instead be dominated by polarization and will achieve much higher signal-to-noise, the main contributor -the EB quadratic estimator -is effectively a shear estimator and therefore naturally suppresses foreground biases.] This paper is organised as follows.In Sec.II we review CMB lensing reconstruction and multipole estimators in the flat-sky approximation and show how to generalise these to the curved sky.For the full-sky shear estimator, we show that it provides an unbiased recovery of the lensing power spectrum using lensed CMB maps.Section III then explores ways of improving the shear estimator by constructing a separable approximation to the m = 2 estimator using singular-value decomposition.Finally, in Sec.IV we test the efficacy of the full-sky shear and m = 2 estimators in suppressing extragalactic foreground contamination by measuring the lensing power spectrum using CMB simulations injected with foregrounds from the Websky simulation [27]. in Appendix A we discuss the estimator normalisation and reconstruction noise power.In Appendix B we derive the form of the full-sky shear estimator in spherical-harmonic space.

II. THEORY: SHEAR AND MULTIPOLE ESTIMATORS
In this section, we introduce the standard lensing quadratic estimator and decompose it into a series of multipole estimators in the flat-sky limit.Starting from this, we then derive the full-sky form of the multipole (m = 2) and shear estimators.
A. Multipole and shear estimators in the flat-sky limit Lensing produces local distortions in the CMB two-point function, breaking the statistical isotropy of the unlensed CMB field and hence introducing new correlations between different Fourier modes over a range of wavenumbers defined by the lensing potential.Averaging over an ensemble of temperature fields for a fixed lensing potential ϕ results in off-diagonal correlations in the observed temperature field T : with and C T T ℓ the temperature power spectrum. 1The standard quadratic estimator [24] exploits the coupling of otherwise-independent temperature modes to reconstruct the lensing field φ by combining pairs of appropriately filtered temperature fields: where A QE L is a multipole-dependent normalisation to make the reconstructed field unbiased and C total ℓ denotes the total temperature power spectrum, including residual foregrounds and instrumental noise.We follow the standard practice of using upper-case L to denote a lensing multipole and lower-case ℓ to refer to CMB multipoles.
The angular dependence of the lensing response function can be expanded in a Fourier series in θ L,ℓ , the angle between L and ℓ: (3) The expansion only involves even multipoles m ∈ 2N because f ϕ (ℓ + L/2, L/2 − ℓ) is invariant under L → −L, i.e., θ L,ℓ → θ L,ℓ + π.In the limit that the lenses are on much larger scales than the CMB fluctuations they are lensing, L ≪ ℓ, the expansion (3) is dominated by the m = 0 (monopole) and m = 2 (quadrupole) moments.The former corresponds to isotropic magnification or demagnification and the latter to shear.Expanding in x ≡ L/ℓ, we have which involves only m = 0 and m = 2 terms at leading order.The multipole expansion of the response function gives rise to a family of lensing estimators characterised by the multipole m, generally of the form where the normalisation A m L is chosen so that the estimator is unbiased: The multipole weight functions g m L,ℓ may be chosen in various ways.In Ref. [21], the minimum-variance estimator at each multipole is constructed, in which case the ℓ-dependence of g m L,ℓ follows where the integral here depends only on the magnitudes ℓ and L. We use a simpler form in this paper (with a relatively minor impact on optimality) whereby we replace C total |ℓ±L/2| appearing explicitly in the weight function in Eq. ( 7) with C total ℓ , which is correct for their product to O(x2 ).Reference [21] show that most of the information in the lensing reconstruction is captured by the m = 0 and m = 2 multipole estimators, even for smaller-scale lenses where the squeezed limit L ≪ ℓ does not apply.
It is convenient to split the QE estimator into this family of multipole estimators because some multipoles are more affected by foregrounds than others.The m = 2 estimator, for instance, is expected to be more robust to extragalactic foregrounds since they primarily bias the m = 0 estimator [21].We discuss this further in Sec.II B below.
The above multipole estimators are generally not easy to implement efficiently because they are non-separable expressions of L and ℓ.To allow for fast evaluation with real-space methods, we first consider the squeezed-limit of the m = 2 estimator.In this case, we use the approximate form of f 2 L,ℓ given by the leading-order quadrupole part of Eq. ( 4): which is clearly separable in L and ℓ.We make a further simplication [21], replacing T (ℓ + L/2)T (L/2 − ℓ) with T (ℓ)T (L − ℓ) to allow fast evaluation of the estimator.Note that this is not simply a variable transformation as we do not change the arguments of the weight function or the angle θ L,ℓ .The foreground deprojection argument still holds at leading order in this case (see Sec. II B).With these modifications, we obtain the shear estimator where the shear weight function is The shear normalisation is obtained from Eq. ( 6).Equation ( 9) can be evaluated efficiently by first noting that the angular term cos(2θ L,ℓ ) can be expressed in terms of the contraction of the symmetric, trace-free tensors L⟨i Lj⟩ and l⟨i lj⟩ (where overhats denote unit vectors in this context and the angular brackets denote the symmetric, trace-free part of the tensor): We then simplify the shear estimator as follows: where the filtered temperature field is The last line of Eq. (12) shows that the shear estimator is equivalent to extracting the E-mode part of the product of the temperature field and the symmetric, trace-free derivative of the filtered temperature field.Expressing the estimator in this form makes its translation to the curved sky rather straightforward, as we discuss in Sec.II C. The shear estimator can be evaluated very efficiently, and has excellent immunity to extragalactic foregrounds.However, approximating f 2 L,ℓ with its squeezed-limit f shear L,ℓ in the weight function results in poor performance at high L, where L ≪ ℓ is not a good approximation [21].In Sec.III, we suggest an alternative, separable approximation to f 2 L,ℓ that performs better than the shear estimator at high L.The approximation we develop there takes as its starting point the asymmetric form of the standard quadratic estimator: We then replace C total |L−ℓ| → C total ℓ and expand the asymmetric lensing response function in Fourier series which now involves both even and odd multipoles.In the squeezed limit, the expansion is still dominated by m = 0 and m = 2 and the expressions for f 0 L,ℓ and f 2 L,ℓ reduce to their symmetric counterparts f 0 L,ℓ and f 2 L,ℓ , respectively.We expect the majority of the signal-to-noise to remain in these multipoles even for smaller-scale lenses (as was the case for the symmetric estimator above).The argument for foreground immunity of the m = 2 estimator still holds in this asymmetric case, also (see Sec. II B).
We end this section by noting that an alternative to considering the m = 2 estimator (plus higher-multipole estimators) is to remove the m = 0 contribution from the standard QE, as was done in [29].However, it is difficult to find an efficient implementation of this scheme since it requires a very accurate separable approximation to the monopole estimator as any error in the monopole removal will cause leakage of foreground contamination.

B. Foreground immunity
To see why extragalactic foregrounds predominantly affect the m = 0 multipole estimator, we consider a simple toy model of extragalactic sources all with the same circularly symmetric angular profile F (x) in the temperature map.Assume that these are Poisson sampled from a fluctuating-mean source density n(x) = n[1 + bδ(x)] where n is the global mean source density, δ(x) is a projected density field correlated with the CMB lensing potential, and b is a linear bias.If we average over the Poisson fluctuations at fixed δ(x), the two-point function of the source contribution to the temperature map, f (x), satisfies to first order in δ and for ℓ 1 ̸ = 0 and ℓ 2 ̸ = 0.If we consider applying a multipole estimator φm (L) to a map composed of CMB, instrument noise and foregrounds f , the foreground bias in the correlation of the reconstructed ϕ with the true ϕ is ⟨ φm L (f, f )ϕ(L ′ )⟩, where φm L (f, f ) is the multipole estimator applied to the field f .At leading order, and for L ̸ = 0, this bias is so the bias predominantly appears in the m = 0 estimator.Note that this holds true even at relatively large L provided that the source profile F (x) is very compact.It also remains true if we replace to allow fast evaluation of the estimator, as discussed above.In this case, Eq. ( 17) becomes Although the expansion of F (ℓ)F (L − ℓ) now introduces terms suppressed by only one power of x = L/ℓ, these have m = 1 angular dependence and so do not bias the m = 2 estimator.The m = 2 foreground terms are still suppressed in the integrand by (L/ℓ) 2 .
If instead, we consider the power spectrum of the reconstructed lensing potential, foreground biases similar to those above arise from terms like ⟨ φm L (f, f ) φm L ′ (T, T )⟩ since averaging over the unlensed CMB reduces the second quadratic estimator to ϕ(L ′ ).These terms are similarly suppressed for m ̸ = 0.However, additional "secondary bispectrum" terms, of the form ⟨ φm L (T, f ) φm L ′ (T, f )⟩, where the CMB and Poisson fluctuations are averaged across quadratic estimators, are not suppressed.However, these biases are generally small; see [21] and Sec.IV.Foreground trispectrum biases also arise that involve the connected four-point function of the foregrounds ⟨ φm L (f, f ) φm L ′ (f, f )⟩ c .These are also suppressed in the limit where shot noise dominates.In this case, our toy model gives and the trispectrum bias separates to involve the same integral as in Eq. ( 17): C. Full-sky formalism In this section, we construct the full-sky versions of the shear and m = 2 estimators.We start with the real-space form of the shear estimator in the flat-sky limit, Eq. ( 12).It is easy to extend this to the curved sky, using sphericalharmonic functions as a basis, converting the partial derivatives into covariant derivatives on the sphere ∂ → ∇, and the integration measure from d 2 x → d 2 n.These changes give This expression can be evaluated efficiently by writing the symmetric, trace-free derivatives of the spherical-harmonic functions in terms of spin-±2 spherical harmonics and then using (fast) spherical-harmonic transforms; see Appendix B. Evaluating the resulting Gaunt integral in terms of Wigner-3j symbols allows us to derive the following harmonic-space form of the full-sky shear estimator: (−1) M g shear ℓ1,ℓ2 (L) Here, the full-sky shear weight function g shear ℓ1,ℓ2 (L) has the same structure as its flat-sky counterpart in Eq. ( 10): where ω 2 ℓ ≡ ℓ(ℓ + 1).The 3j symbol here enforces mode coupling, the spherical analogue of ℓ 1 + ℓ 2 = L, and the expression in square brackets accounts for cos(2θ L,ℓ1 ) noting that, by the cosine rule, With the usual curvature correction, ℓ → ω ℓ , and assuming that all multipoles are much larger than √ 2, we recover the correspondence to Eq. ( 9).The explicit form of the spherical normalization, A shear L , is given in Appendix A. One can similarly derive the full-sky version of the m = 2 estimator, obtained by replacing f shear L,ℓ by f 2 L,ℓ (or the asymmetric form f 2 L,ℓ ) in the reconstruction weight function.We now have where the non-separable filtered temperature field T F,m=2 L (n) is given by The harmonic-space form is where the weight function g m=2 ℓ1,ℓ2 (L) is no longer a separable function of L, ℓ 1 and ℓ 2 , and is given by This m = 2 estimator performs better than the shear estimator on small scales but is computationally expensive to evaluate due to the non-separability of the weight function.Reconstruction at each L requires evaluation of a separate filtered temperature field, T F,m=2 L (n), making the entire reconstruction a factor of L max slower than for the separable shear estimator.The reconstruction would therefore scale as O(L 4 max ), which is generally infeasible, particularly as the reconstruction typically has to be run many times for simulation-based removal of biases in the reconstructed lensing power spectrum.In Sec.III, we suggest a simple way of writing the non-separable f 2 L,ℓ (or, actually, f 2 L,ℓ ) as a sum of separable terms, allowing reconstruction to be performed with far fewer filtered fields.
Finally, for completeness, we point out that the same prescription can be extended to calculate higher-multipole estimators even though, as noted above, most of the signal-to-noise is contained in the monopole and quadrupole.As an example, consider the m = 4 estimator.On the flat sky, the real-space form of the estimator follows from noting that cos(4θ L,ℓ ) = 8 L⟨i Lj Lk Ll⟩ l⟨i lj lk ll⟩ .
This allows us to express the m = 4 response function directly in terms of ℓ and L resulting in the following real-space estimator on the full sky (for L ≥ 4): where the filtered field is now Here, we have simply replaced the flat-sky ℓ 4 , which comes from converting multiplication by the unit vector l to a derivative, with ω 4 ℓ and similarly for L 4 .The harmonic-space form is with weight function Note that in going from the flat-sky to full-sky, we could have chosen to use (ℓ + 4)!/(ℓ − 4)! instead of ω 4 ℓ for the terms that arise with the 4th derivatives.The fractional differences are O(1/ℓ 2 ) and would produce only small changes in optimality (at intermediate and small scales) while simplifying the above weight functions.
We test the full-sky shear estimator using full-sky lensed CMB simulations with the specifications of an upcoming Stage-3 experiment, with 1.4 arcmin beam width (full-width at half maximum) and 7 µK-arcmin instrument white noise at a frequency of 148 GHz.We first test the estimator using full-sky temperature maps without foreground contamination.The cross-power spectrum between the reconstructed convergence field κ (related to the lensing field via κ = −∇ 2 ϕ/2) and the input agrees with the theory lensing spectrum at the percent level.This shows that our estimator is correctly normalised (we use the full-sky normalisation given in Appendix A).We further verify that we can recover an unbiased estimate of the convergence power spectrum, C κκ L = L 2 (L + 1) 2 C ϕϕ L /4, by subtracting several biases from the empirical auto-spectrum of the reconstruction, C κκ L : L .
Here, N L is the Gaussian bias produced from the disconnected part of the CMB four-point function that enters ⟨C κκ L ⟩.This can be thought of as the power spectrum of the statistical reconstruction noise sourced by chance, Gaussian fluctuations in the CMB and instrument noise that mimic the effects of lensing.We estimate this bias by forming different pairings of the simulation that is being treated as data and independent simulations following the method described in [15].We use 100 different realisations of simulation pairs to obtain an average over simulations.The N (1) bias, which arises from the connected part of the CMB four-point function and at leading order is linear in the lensing power spectrum [30], is estimated using 200 pairs of simulations, with the same lensing realisation for each member of the pair but different unlensed CMB realisations, based on [31].The de-biased bandpowers of the shear reconstruction are shown in Fig. 1.
FIG. 1. Lensing convergence reconstructed on the full sky with the shear-only estimator (Eq.21) applied to a simulated, foreground-free, noisy CMB temperature map.The cross-power spectrum of the reconstructed and input convergence is shown in green, and bandpowers with ∆L = 60 of the auto-power spectrum of the reconstruction are in red.The latter is corrected for N (0) and N (1) bias terms (shown in blue and orange, dashed, respectively).Both the auto-and cross-spectra agree well with the expected power (black).

III. BEYOND THE LARGE-SCALE-LENS REGIME: SVD EXPANSION OF THE MULTIPOLE KERNELS
The shear estimator in Eq. ( 21) is separable and so efficient to evaluate, but this comes at the cost of increased noise in the reconstruction on small scales.The sub-optimality of the shear estimator is apparent from Fig. 2, which shows that its disconnected noise bias N (0) L has a spike at small scales (see also Fig. 2 in [21]).This spike arises because the shear estimator has zero response to lenses at this particular scale.The noise biases in the figure are computed on the full sky using Eq.(A8).
The full m = 2 estimator, which in this paper we approximate by Eq. ( 25) with non-separable weights (26), has better noise performance than the shear estimator for L > 100 (for the survey specifications adopted here) as shown in Fig. 2. In particular, if the weights are constructed from the m = 2 component of the asymmetric lensing response function, f 2 L,ℓ , the noise spike is eliminated.However, such m = 2 estimators are inefficient to evaluate since the weights are not separable.
A simple work-around, which we have found to perform quite well, is to retain the squeezed-limit, separable approximation (i.e., the shear estimator) on large scales where its performance is similar to the full m = 2 estimator, but to approximate the f 2 L,ℓ as a sum of separable terms on smaller scales.We find these separable terms by singularvalue decomposition (SVD) [32].
In detail, we construct a hybrid approximation to f 2 L,ℓ as follows.For L < L * , we use the shear approximation; for L > L * , we perform a singular-value decomposition of the block of the asymmetric, convergence response function, 2 f 2 L,ℓ /ω 2 L , with L * ≤ L ≤ L max and 2 ≤ ℓ ≤ ℓ max , and approximate this with the first n largest singular values.We found that keeping n = 20 SVD terms gives a reasonable balance between computational efficiency and optimality.We chose the lensing multipole L * at which to switch such that the reconstruction noise on the SVD-based estimator is lower than that of the shear, which for our survey parameters is L * = 1000.(Note that there is a range of L in which the above condition is true, and L * was chosen empirically based on good SVD convergence.Other metrics could certainly be used to determine n and L * .)The complete form of our hybrid-SVD response function is Here Λ i corresponds to the ith singular value, U L,i are the components of the ith left singular vector and V ℓ,i are the components of the ith right singular vector.We see that the SVD naturally decomposes the response into a sum of separable terms and hence the reconstruction can be performed efficiently for each component.The separable SVD estimator is given explicitly by (for where the filtered field for the ith singular value is The normalisation A SVD L is chosen, as usual, to ensure the estimator has the correct response to lensing.We show in Sec.IV that the separable approximation is still very effective at suppressing foregrounds, as expected since we have not altered the m = 2 geometric structure of the estimator.As can be seen from Fig. 2, the noise performance is very close to that of the full (asymmetric) m = 2 estimator for L ≥ L * , and around a factor of two better than the shear estimator.

IV. TESTING THE SENSITIVITY TO FOREGROUNDS USING SIMULATIONS
To test the sensitivity of the estimators to extragalactic foregrounds, we use the component maps of the Websky extragalactic foreground simulations [27].These include CIB, tSZ and kSZ at 143 GHz.The power spectra of these foregrounds are shown in Fig. 3.In a real analysis, bright galaxy clusters and sources would be dealt with either by masking (i.e., excising regions around them) or in-painting (masking, but with the resulting holes filled with constrained realizations).We mimic this for point sources in our analysis, without introducing the complications of having to deal with masked or in-painted maps, as follows.We apply a matched-filter, with the profile corresponding to the instrumental beam and noise power given by the sum of instrumental noise and foreground power, to maps including the full lensed CMB plus foregrounds plus white noise.Sources with recovered flux density greater than 5 mJy are catalogued and regions around them are removed from the foreground maps only.These masked foreground maps are then combined with lensed CMB and 7 µK-arcmin white noise to form the final temperature map given by T total = T CMB + T f + T noise , where we have written explicitly the contributions to the observed temperature map FIG. 2. Lensing (convergence) reconstruction noise power spectra for the standard quadratic estimator (QE; black), the shear estimator (green), the m = 2 symmetric estimator (symmetric f 2 ; orange), the m = 2 asymmetric estimator (asymmetric f 2 ; blue), and the SVD approximation to the latter (asymmetric f 2 -SVD; red).The lensing convergence power spectrum is also plotted (dashed black).Note that the convergence spectra are related to the spectra of the lensing potential by, e.g., The survey specifications are the same as in Sec.II C.
from the lensed CMB T CMB , the extragalactic foregrounds T f and the detector noise T noise .We do not mimic the masking of bright galaxy clusters.As noted below, this means that our results for the bias should be considered rather extreme, particularly for the trispectrum bias.
To assess the bias induced by foregrounds in the auto-power spectrum of the reconstruction, Ĉκκ L , we evaluate the primary foreground bispectrum 2⟨Q[T f , T f ]κ⟩ and the foreground trispectrum term from which the disconnected (Gaussian) contribution is subtracted using simulations.Here Q[T A , T B ] represents a quadratic estimator (we consider the standard quadratic, shear and hybrid-SVD estimators) applied to maps T A and T B .We do not consider the secondary bispectrum bias discussed in [21], as it was found to be subdominant to the primary bispectrum and the trispectrum biases (and we expect the same to hold for our estimator variants).
The primary bispectrum bias on the lensing power spectrum can be seen in Fig. 4 for three choices of the maximum CMB multipole used in the reconstruction: ℓ max = 3000, 3500 and 5000.Power spectra are binned with ∆L = 60.For all the ℓ max choices, significant biases are observed in the standard quadratic estimator, while the shear estimator can remove the bias very effectively, in agreement with the flat-sky results of [21].Furthermore, one can see that the bias induced in the hybrid-SVD estimator is smaller than that of the shear on small scales and comparable to that of the shear on large scales.The improvement in noise performance of the hybrid-SVD estimator compared to the shear estimator can also be appreciated, where the shaded 1 σ bandpower errors (which include reconstruction noise and lensing sample variance) for the hybrid-SVD estimator (red) lie between the shear (green) and the standard quadratic estimator (blue), the latter having the lowest variance but a large foreground bias.Similar improvements can be seen in Fig. 5, where the trispectrum bias reduces significantly when switching from the standard quadratic estimator to the shear or hybrid-SVD estimator.Although for ℓ max = 3500 and ℓ max = 5000, the bias is no longer significantly smaller than the statistical error, the improvement compared with the standard quadratic estimator is still large.Furthermore, it should be noted that the trispectrum bias is particularly sensitive to the most massive galaxy clusters, which can be straightforwardly detected and removed by masking or inpainting in a real analysis.As we have not carried this out here, we expect a significant reduction in trispectrum bias for all estimators in practice.
FIG. 3. Power spectra of the tSZ (orange), kSZ (green) and CIB (red) extragalactic foregrounds at 148 GHz from the Websky simulation [27].The lensed CMB power spectrum is also shown (in blue).The combined total map is matched-filtered and masked for point sources above a flux density of 5 mJy in our tests.FIG. 4. Relative primary bispectrum bias on the CMB convergence power spectrum due to the combined effect of foregrounds for the standard quadratic estimator (QE; blue), the shear estimator (green) and the hybrid-SVD estimator based on SVD of the asymmetric f 2 L,ℓ response (SVD asymmetric; red).From top to bottom, we vary the maximum CMB temperature multipole used in the reconstruction, ℓmax = {3000, 3500, 5000}.The shaded bands indicate the 1 σ statistical reconstruction error for the different estimators for bandpowers of width ∆L = 60.The bias is well contained within the statistical errors for the shear and hybrid-SVD estimators but there is significant bias in the standard quadratic estimator.FIG. 5.As Fig. 4 but for the relative trispectrum bias.Similar to the bispectrum biases, the shear and hybrid-SVD estimators suppress the trispectrum bias significantly across most of the lensing multipole range.Note that the bias for the standard quadratic estimator with ℓmax = 5000 is off the scale of the plot.It is worth noting that the trispectrum bias is particularly sensitive to bright tSZ clusters and further improvement can be obtained via cluster masking.

V. DISCUSSION AND CONCLUSIONS
We showed how to formulate foreground-immune multipole estimators for CMB lensing reconstruction, particularly the m = 2 estimator that contains most of the signal-to-noise, on the spherical sky.This allows the straightforward application of the estimators proposed in [21] to large-area surveys such as Planck, AdvACT and the forthcoming Simons Observatory.Generally, these estimators are not separable and so cannot easily be evaluated efficiently.Previous separable approximations -the shear estimator introduced in [21] -have sub-optimal reconstruction noise when reconstructing small-scale lenses.We presented a simple, first attempt at producing a separable approximation to the full m = 2 estimator based on singular-value decomposition of the part of its response function at intermediate and large lensing multipoles.We tested the performance of this hybrid-SVD estimator, along with the shear approximation and the standard quadratic estimator, on the Websky [27] foreground simulation.As in the flat-sky tests considered in [21], we found the shear estimator to be very effective in suppressing foreground biases even on single-frequency maps.The same is true of the hybrid-SVD estimator, but it has the advantage of higher signal-to-noise on small scales.
The field of CMB lensing has experienced a fast transition from first detection [5,6] to precision measurements in the last 15 years.With current and upcoming surveys of the CMB, such as AdvACT, SPT-3G and Simons Observatory, probing the millimetre sky with increasing resolution and sky coverage, we can expect further rapid improvements in the quality of lensing products reconstructed from the CMB.However, improvements in statistical noise must be met with more stringent control of systematic effects, such as those from extragalactic foregrounds in temperature maps.The methods explored in this paper provide a robust way to measure CMB lensing, which is largely immune to the effect of these foregrounds.They can be added to the existing repertoire of methods to mitigate foregrounds, such as multi-frequency cleaning and bias hardening, and can be used in combination with these to improve optimality further.For example, it was found in [19] that a robust estimator to reduce foreground biases while having a low impact on signal-to-noise tends to consist of a combination of bias hardening (for point sources and tSZ cluster profiles), explicit tSZ deprojection in multi-frequency foreground cleaning and the shear estimator.[Furthermore, Ref. [33] pointed out that extragalactic foregrounds in temperature can also bias B-mode delensing efforts for upcoming experiments like Simons Observatory, leading to bias in the inferred amplitude of the power spectrum of primordial gravitational waves if the non-Gaussianity of extragalactic foregrounds is not accounted for.These biases can also be effectively mitigated without significantly compromising the delensing efficiency if the lensing map is obtained using the shear-only estimator described here.] [In this work, we did not explore the impact of polarized extragalactic foregrounds.This is because for current and next-generation surveys, like ACT and Simons Observatory, the temperature-based (T T ) estimator carries significant weight (about 2/3 of the statistical weight for an ACT-like survey) and the foreground bias in temperature-based reconstruction is larger than for polarization since the degree of polarization of extragalactic foregounds is small compared to the CMB fluctuations.However, for future deep surveys, such as CMB-S4, the polarization channel will play a dominant role in lensing reconstruction when noise and foreground levels become below the lensing B-mode levels of around 5 µK-arcmin [34], giving significant improvements in the signal-to-noise of the reconstruction.Although extragalactic-foreground biases to polarization-based lensing reconstructions are small compared to temperature, particularly for the EB estimator, which is naturally a shear estimator and will carry most of the statistical weight for future surveys, the lower statistical reconstruction noise means the absolute requirements on foreground bias will be much more stringent.Extending the shear estimator for polarization-based estimators would enable a tighter control of potential polarized extragalactic-foreground contamination (coming mainly from bright radio and infrared point sources) as shown in [26] for the flat-sky case.There it was found there that using the minimum-variance combination of shear estimators incurs only a modest 20% noise penalty compared to the standard minimum-variance estimator for a CMB-S4-like survey.An extension of such polarization-based shear estimators to the full-sky case would be interesting but is deferred to future work.]polar coordinates (θ, ϕ).Expanding the filtered field in (B1) in terms of spherical harmonics, where, as discussed in the main text, the term in square brackets on the right accounts for the cos(2θ L,ℓ1 ) weighting in the flat-sky limit.With this simplification, the shear weight function reduces to g shear ℓ1,ℓ2 (L) = 1 2 (2L + 1)(2ℓ 1 + 1)(2ℓ 2 + 1) 16π given as Eq. ( 23) in the main text.