that of the likelihood model).

Figure 14:  Lensing potential estimate pseudo-power spectra for QE (left
panel) and MAP (right panel) as obtained from the CMB-S4 LAT simulations. The
solid dark line is the fiducial lensing potential. We show the result for all
three simulation sets: $\rm{M}_{00}$ in dot dashed, $\rm{M}_{07}$ solid, and
$\rm{M}_{09}$ dashed. However, their spectra often overlap and are barely
distinguishable, particularly so in the QE case. Pink (blue) shows the QE or
MAP lensing potential before (after) mean-field subtraction, gray (green)
shows the QE mean-field or MAP residual mean-field with (without) Monte-Carlo
noise. The mean-field drops off quickly for multipoles larger than $L\simeq
10$. The QE mean-field (green and gray) contains a large signal at large
angular scales, and is the lensing potential’s (pink) main component on that
scale. On large scales, the residual MAP mean-field signal is of the order of
the lensing potential, many orders of magnitude smaller than its QE counter
part (green, left panel). This means that at each iteration, the QE mean-field
is a good proxy for the true mean-field term.

Figure 15:  Absolute difference maps for the QE (left column) and MAP (right
column) deflection-mean-fields, ${\rm MF}_{07}-{\rm MF}_{00}$ (top panels),
and ${\rm MF}_{09}-{\rm MF}_{00}$ (bottom panels). There are clear features
seen in every panel, which originate at least in part from the foreground
residuals.

We now turn to three quantitative tests on the impact of the mean-field on the
delensing efficiency. The changes applied are the following.

  * (I)

A posteriori mean-field subtraction: we obtain the average of the MAP
reconstructions obtained from Eq.(58) and subtract it from each lensing
potential estimate (a procedure similar to the QE mean-field subtraction).

  * (II)

Swapped mean-field: we use the simplest simulation set available (the Gaussian
foreground model $\rm{M}_{00}{}$) to build the mean-field template for the
reconstructions on all simulation sets. Hence, the mean-field subtraction will
not remove the non-Gaussianity of the foregrounds (and also be wrong in a
broader sense, since the different models have different power).

  * (III)

Mean-$B$-lensing template subtraction (similar to point (I), but at the level
of the $B$-lensing templates): we subtract from each $B$-template the average
of templates obtained on a set of reconstructions.

In all of these tests, we calculate the residual lensing amplitude on the
large, green contoured sky patch $p_{50}$ (see Fig. 9). This includes some
mean-field signals in the outskirt of the SPDP, yielding conservative results
due to the small weight of the SPDP edge for the $r$-inference. For
convenience, we only use 100 simulations, with almost exactly same
configuration as described in Sec. IV.2 (here, we boosted the numerical
accuracy parameters to obtain more precise large scale lenses). We obtain
$A_{\rm lens}$ from the ratio of power, following Eq. (56), along with the
ensemble standard deviation. We summarize our findings in terms of a single
value $\frac{\Delta A_{\rm lens}}{A_{\rm lens}}$ in Table 4 that compares the
tests to the baseline configuration. Here, $\Delta A_{\rm lens}$ is the
difference in residual lensing amplitude between the baseline and modified
analyses.

(I): Mean-field subtraction:
For test (I), we subtract the MAP mean-field from the lensing-potential
estimates prior to the lensing operation, and calculate the $B$-lensing
templates, denoted $B^{\rm LT}_{\rm MF^{-}}$. Explicitly,

$\displaystyle B^{\rm LT,MAP}_{\rm MF^{-}}$
$\displaystyle=\left(\hat{}\bm{\alpha}^{\rm MAP}-\hat{}\bm{\alpha}^{\rm
MF}\right)\circ\hat{E}^{\rm unl}.$ (61)

The result is shown in the first row of Table 4. The impact is always small,
and mostly consistent with zero (at $2\sigma$ for simulation set
$\rm{M}_{00}$). The given error bars are the one-simulation-equivalents, and
therefore, averaging across all simulations, we would detect a tiny
improvement in all cases. Nevertheless, this confirms that simply subtracting
the QE mean-field as calculated for the starting point is sufficient.

(II): Swapped mean-field:
We now perform the lensing reconstructions with a mean-field different to the
baseline. For simulation sets $\rm{M}_{07}$ and $\rm{M}_{09}$, we perform
reconstructions where the mean-field gradient term is set to that of model
$\rm{M}_{00}$ (doing so also changes the starting point for the quadratic
estimator but this alone does not affect the converged solution).

Thus, the lensing reconstruction is performed on an inhomogeneous and non-
Gaussian foreground simulation set, but its mean-field is calculated from a
simpler simulation set (with differences shown in Fig. 15), with homogeneous
and Gaussian foregrounds. The result is shown in Table 4 in the second row,
and is consistent with zero. This demonstrates that the choice of mean-field
in the lensing reconstruction does only have a small impact. The mean-field
from $\rm{M}_{00}$ only traces the patch mask and noise variance map, yet it
nevertheless yields similar results for all simulation sets.

(III): Mean-$B$-lensing Template Subtraction:
Figure 16 shows an example $B$-lensing template for $\rm{M}_{07}$ and
$\rm{M}_{09}$. The MAP mean-$B$-lensing templates for $\rm{M}_{07}$ and
$\rm{M}_{09}$, together with the MAP difference maps,
$\rm{M}_{07}{}-\rm{M}_{00}{}$, and $\rm{M}_{09}{}-\rm{M}_{00}{}$ are shown in
Fig. 17.

Figure 16:  $B$-lensing templates for $\rm{M}_{07}$ (left panel) and
$\rm{M}_{09}$ (right panel). All plots are band-passed to $30<\ell<200$.

Figure 17:  _Left column_ : MAP mean $B$-lensing templates for $\rm{M}_{07}$
(top) and $\rm{M}_{09}$ (bottom). _Right column_ : Difference maps between the
mean $B$-lensing templates $\rm{M}_{07}$ and $\rm{M}_{09}$ with $\rm{M}_{00}$.
All plots are band-passed to $30<\ell<200$.

For $\rm{M}_{09}$ (bottom left panel) we see in Fig. 17 clear features along
the edge of the patch, hinting at contamination from foreground residuals.
$\rm{M}_{07}$ does not show any dominant features. The difference plots (right
panels) show that for $\rm{M}_{07}$ minus $\rm{M}_{00}$, the mean-$B$-lensing
template follows the foreground amplitude modulation of $\rm{M}_{07}$. For
$\rm{M}_{09}$ minus $\rm{M}_{00}$, we see bright spots at about the same areas
where foreground contamination appears in the mean-field of $\rm{M}_{09}$,
also clearly visible in the mean-$B$-lensing template in Fig. 16. We also see
a small amplitude ringing effect along the latitude of the SPDP, which we
believe is an artefact of the map-building procedure, introduced when modes
below $30$ were excised from the maps.

We test for the impact of these features by subtracting the mean-$B$-lensing
template from each $B$-lensing template,

$B_{\ell m}^{\rm LT,MAP}\rightarrow B_{\ell m}^{\rm LT,MAP}-\left\langle
B_{\ell m}^{\rm LT,MAP}\right\rangle_{\rm MC},$ (62)

where in the average on the right-hand side we always exclude the simulation
index for which $B^{\rm LT,MAP}$ is calculated.

However, similar to the QE and MAP mean-field, subtracting does not only
remove contaminants, but also adds some level of Monte-Carlo noise. For this
reason, we calculate power spectra by dividing the simulation sets into two
splits and taking cross-spectra. The result is shown in the third column of
Table 4. We do not find any improvement, and the error bars are rather large.

Table 4:  Mean-field test results. We give the relative change in terms of residual lensing amplitude, compared to the ‘baseline’ configuration, $\frac{\Delta A_{\rm lens}}{A_{\rm lens}}$, a negative number indicating an improvement. See text for the definition of the individual tests. Test |  $\rm{M}_{00}$ |  $\rm{M}_{07}$ |  $\rm{M}_{09}$
---|---|---|---
I: | $-0.0028\pm 0.0014$ | $-0.0005\pm 0.0014$ | $-0.0009\pm 0.0014$
II: |  | $\ \ \ 0.0001\pm 0.004$ | $\ \ \ 0.0039\pm 0.01$
III: | $-0.0005\pm 0.02$ | $-0.0002\pm 0.015$ | $-0.0028\pm 0.03$

Summarizing these tests, they confirm that our results are robust against
changes to the mean-field treatment, and more generally that the impact of the
mean-field is small more generally.

###  IV.3 Chile Deep Patch Delensing

As a last test, we repeat the calculation of the $B$-lensing template for a
different CMB-S4 configuration. Specifically, we assume that the Deep Field
would be observed from Chile, referring to this as the Chile Deep Patch (CDP).
Observing the sky from Chile would, among other things, change the scan
strategy of the telescopes, therefore impacting the noise levels and observed
sky fractions. Figure 18 shows the derived noise-level map after component-
separation, in equatorial coordinates. The central noise level is
$\hat{n}_{\rm lev}=0.6\ \mu$K-arcmin and therefore about $50\%$ higher than
for SPDP, and the scan area is much broader.

We only use the foreground models $\rm{M}_{00}$ and $\rm{M}_{07}$, and repeat
lensing reconstruction for 100 realizations of the CDP component-separated
maps. We update our noise model and overlapping $B$-mode deprojection matrix
to this case, and leave other ingredients of our algorithm essentially the
same.

Figure 18:  This is a full sky Noise-level map for the Chile Configuration and
in equatorial coordinates. The contours show edges that contain the masks A
(blue), B (red), C (green), and D (yellow), and the numbers in the color bar
refer to the contour noise level. The lowest (highest) noise of about 0.65
($13.19$)$\ \mu$K-arcmin is indicated by the blue (dark gray) contour and all
are derived from the Chile hits-count map.

To improve numerical stability for the iterative lensing reconstruction step,
we have increased the mask by excluding areas for which the noise is more than
$100$ times larger than the central noise value. These are mostly pixels near
the galactic plane with high foreground residuals that can safely be
discarded. Iterative lensing reconstruction succeeds on the entire resulting
area.

The bulk of the constraining power comes from the red region in Fig. 18. The
larger sky patches indicated by the yellow and green contours partly cover the
Galactic plane. Although lensing reconstruction does converge there as well,
the foreground residuals are too large in these areas, leading to
uninteresting delensing results.

In the red region and for $\rm{M}_{07}$, we find303030The results for QE are
calculated using non-perturbative remapping of the $E$-mode, improving
performance slightly, as before. $A_{\rm lens}^{\rm QE}=0.30\pm 0.02$, $A_{\rm
lens}^{\rm MAP}=0.121\pm 0.007$. For $\rm{M}_{00}$, we find $A_{\rm lens}^{\rm
QE}=0.29\pm 0.02$, $A_{\rm lens}^{\rm MAP}=0.115\pm 0.006$.

These empirical results for iterative delensing are in reasonable agreement
with our theoretical predictions, ${A}^{\rm MAP}_{\rm lens}=0.11$ for
$\rm{M}_{07}$, and ${A}^{\rm MAP}_{\rm lens}=0.10$ for $\rm{M}_{00}$, matching
as well those expected from the CMB-S4 $r$-forecast paper Abazajian  _et al._
(2022).

###  IV.4 Constraints on $r$

In our previous CMB-S4 $r$ constraint paper Abazajian  _et al._ (2022) one of
the two re-analysis methods used was the parametric multi-component cross-
spectral likelihood that was first introduced in Ref. (BICEP2/Keck
Collaboration  _et al._ , 2015), and which has been used in all BICEP/Keck
analyses up to the most recent (Ade  _et al._ , 2021). Here we take this
likelihood and add the lensing templates derived in the previous section as an
additional pseudo-frequency band, in a similar manner to the joint analysis of
BICEP/Keck and SPT data presented in Ref. BICEP/Keck Collaboration  _et al._
(2021).

To obtain sufficient $E/B$ separation purity we calculate the auto- and cross-
spectra between the maps using the $\mathtt{S^{2}HAT}$ library. The resulting
bandpowers are then compared against a likelihood model built using the
Hamimeche-Lewis Hamimeche and Lewis (2008) approximation with the following
parameters: the tensor-to-scalar ratio $r$; the amplitudes of the dust and
synchrotron $B$-mode power spectra, their frequency spectral indices, and
their spatial power law indices; the frequency-independent spatial correlation
between dust and synchrotron; and the dust correlation between 217 and 353
GHz. See appendix G of Ref. Ade  _et al._ (2018) for a full description of the
parameterization.

The likelihood requires statistical characterization of the signal and noise
parts of the lensing template auto-spectrum and of their cross-spectra to the
SAT frequency maps. In this first analysis, we do this in the following crude
(but effective) manner. Making the assumption that the template may be written
as a filtered version of the true CMB lensing $B$-mode plus an additive noise
term, we compute an effective, isotropic filter function from 100 cross-
spectra to their known simulation CMB inputs (neglecting the precise position
dependence of the filter). This allows us to build predictions for the cross-
spectra to the SAT maps. Further, we obtain the lensing template auto-spectrum
noise power by subtracting from the templates spectra the signal term as
predicted from the isotropic filter model.

Figure 19:  Results from reanalysis of map-based simulations of CMB-S4 observing the South Pole Deep Patch. The mean and standard deviation of the maximum likelihood $r$ values evaluated over ensembles of simulation realizations are shown as the central point and horizontal error bar, respectively. Blue and orange points show the results for $r=0$ and $r=0.003$ realizations, with the vertical lines indicating these values. The darker points are results without the lensing template, while the lighter points show the substantial benefit of including the iteratively-reconstructed templates described in this paper.  Table 5:  Results from reanalysis of map-based simulations of CMB-S4 observing the South Pole Deep Patch, including the reconstructed lensing templates described in this paper.  $r$ value | Sky model | $\sigma(r)\times 10^{4}$ |  $r$ bias $\times 10^{4}$
---|---|---|---
0 | 0 | 4.3 | -1.0
| 7 | 4.7 | -1.3
| 9 | 5.8 | 1.1
0.003 | 0 | 6.4 | 1.3
| 7 | 6.9 | 1.4
| 9 | 7.8 | 4.7

For each of the 500 simulation realizations, and each of the three foreground
models, we perform a maximum-likelihood (ML) search, both with and without the
lensing template. We use the mean and the standard deviation of a given
parameter value to assess the bias and experimental uncertainty of that
parameter—where the parameter $r$ is of course by far of the greatest
interest. The SPDP results are shown in Fig. 19 for the realizations with
$r=0$ and those with $r=0.003$. As expected, the addition of the lensing
template (effective delensing) dramatically reduces the uncertainty $\sigma$
on $r$, here denoted $\sigma(r)$. In the case of the non-Gaussian model
$\rm{M}_{09}$, a significant bias is also removed. Table 5 gives the with-
template results in numerical form. These can be compared to the results in
the rightmost column of Table 2 of our previous paper Abazajian  _et al._
(2022). The results for $\rm{M}_{00}$ are most directly comparable. Previously
for $r=0$ we had $\sigma(r)=5.7\times 10^{-4}$, to be compared to the present
$4.3\times 10^{-4}$, the latter being below the $\sigma(r)=5\times 10^{-4}$
science goal of CMB-S4 Abazajian  _et al._ (2019).

We have preliminary ML search runs for the CDP simulations as well. The
uniform amplitude $\rm{M}_{00}$ model is a gross oversimplification for this
larger sky patch. Nevertheless, for input $r=0$ this model gives
$\sigma(r)=5.2\times 10^{-4}$, which is consistent with the expectation that
the best sensitivity will be obtained with the most concentrated sky coverage.
For $\rm{M}_{07}$ the effective foreground amplitude averaged over the patch
is extremely high and the likelihoods result is not robust. We find that
$\sigma(r)$ is only increased to $7.0\times 10^{-4}$, but there is significant
bias—further investigations are ongoing. No results are available for
$\rm{M}_{09}$ at the time of writing.

##  V Conclusion

In this paper, we have presented an optimal lensing reconstruction software on
the curved sky, usable under realistic conditions, inclusive of inhomogeneous
noise, and on sky patches of arbitrary size. With it, we achieve the CMB-S4
science goal on $r$ by removing $92\%$ to $93\%$ of the $B$-lensing power,
depending on the simulation set.

The method mostly generalizes that of Ref. Carron and Lewis (2017) to curved-
sky geometry. We discussed carefully all ingredients, and showed that the few
approximations that are made are perfectly under control given the moderate
size of $\Lambda$CDM deflection angles (as is often the case, the sky
curvature can be neglected on the scale of the deflection), and that the
reconstruction should be optimal provided that the CMB likelihood model is a
fair representation of the data.

We found most difficulties of practical nature in the reconstruction process
at low lensing multipoles. There, the mean-field is prominent, and more
generally, the behavior of the reconstruction of the largest scale lenses is
less stable, and also often requires higher accuracy in the most important
steps (lensing remappings, recovery of the delensed $E$-mode) than for
smaller-scale lenses. We have used here as the baseline a compromise, where
the mean-field and largest lenses are not precisely reconstructed, which for
most applications is perfectly satisfactory: these large lenses contribute
little to $B$-mode delensing on degree scales, and, for the purpose of
lensing-spectrum reconstruction, cosmic variance dominates there anyway, so
that the quadratic estimator solution can be satisfactory. Recent work
Reinecke  _et al._ (2023) allowed us to solve most of these issues, making
upcoming analyses easier as well as slightly improving the lensing
reconstruction further.

New to this paper is the characterization of the noise anisotropies induced by
delensing as the only source of the mean-field, which has no direct analog in
quadratic-estimator theory. We showed that it is a negligible source of
complication, in particular for situations where uncertainties in the
statistical properties of the data overshadow its contribution. One is
naturally lead to drastic simplifications in the treatment of the mean-field,
which becomes as simple as for the standard quadratic estimator.

The iterative reconstruction usually converges in a reasonable number of steps
(rarely did we use more than 15 iterations), but is by no means a smooth
process. In fact, the non-Gaussian nature of the likelihood is such that for
low noise levels full Newton steps are much too strong, resulting generically
in unrealistic, non-invertible deflection-field estimates. We found that
reducing the step in a scale-dependent manner can correct that problem, but
this requires undesirable ad-hoc tuning, with little improvements each time,
except for the first few iterations; calculation of the gradients is costly
and generally improvement of the convergence after a few steps would be highly
desirable. This is left for future work. For the purpose of delensing, it
helps that recovery of the lensing $B$-mode appears easier than that of the
lensing deflection field and $E$-mode map individually: errors in both tracers
display statistical dependence in order to still produce the right level of
$B$-modes seen in the data.

Also new to this paper is the joint reconstruction of the lensing post-Born
curl together with the standard lensing-gradient mode. The joint
reconstruction brought no complications with it, and the output is correctly
described by our predictions of the iterative $N^{(0)}$ and $N^{(1)}$ biases.
At the current time the lensing curl is really only useful as a sanity check
of the lensing reconstruction, but at the CMB-S4 depth probed here, it will be
also be a signal, whose detectability is boosted by large factors with our
optimal method.

We focused mostly on simulations of a very deep patch, using semi-realistic
inhomogeneous CMB-S4 noise maps, together with two non-Gaussian foreground
models. Our baseline reconstructions on ILC-cleaned maps are not perfect near
the edges and on large scales. Nevertheless this has little impact overall,
and the tracer is able to remove close to 93% of the lensing over the relevant
area, as expected. Of course, available Galactic foreground models are still
crude on small scales, but we find encouragingly that the main effect of the
presence of the foregrounds on the lensing templates appears to be the
increased effective noise of the maps. This extends cautious optimism of Ref.
Beck  _et al._ (2020) to the case of iterative delensing. For the most complex
model $\rm{M}_{09}{}$ a non-Gaussian feature is present close to the edges,
which is visible when averaging the templates.

We built predictions for the delensing capabilities and lensing bias
calculations, using a slightly modified version of the well-known algorithm
introduced in Ref. Smith  _et al._ (2012). It is quite noteworthy how well
these predictions (which do not exactly follow from first principles), seem to
work even under non-idealized conditions, and it is reassuring to see that
forecasts with these tools can probably be trusted.

To summarize, Ref. Kamionkowski and Kovetz (2016) writes that $B$-mode
delensing is an ‘ambitious, sophisticated, and challenging endeavor’: we
believe that this paper makes the prospect of an eventually successful
delensing much more realistic. Of course, much work remains to be done. Our
noise simulations were inhomogeneous but otherwise did not have the full
complexity expected from a ground-based experiment like CMB-S4. While at least
some of such non-idealities are expected to impact predominantly scales larger
than those most relevant for lensing reconstruction, it is difficult to assess
the importance of these effects for delensing without more realistic
simulations. This is left for future work. Other areas that need development
include the study of systematic effects, better foreground models and cleaning
methods, non-linear large-scale structure simulation inputs, and further code
improvements.

###### Acknowledgements.

The authors thank Victor Buza, Benjamin Racine for previous work relevant to
this paper. CMB-S4 is supported by the Director, Office of Science, Office of
High Energy Physics of the U.S. Department of Energy under Contract No.DE-
AC02-05CH11231; by the National Energy Research Scientific Computing Center, a
DOE Office of Science User Facility under the same contract; and by the
Divisions of Physics and Astronomical Sciences and the Office of Polar
Programs of the U.S. National Science Foundation under Mid-Scale Research
Infrastructure award OPP-1935892. Considerable additional support is provided
by the many CMB-S4 team members and their institutions. This research used
resources of the National Energy Research Scientific Computing Center, a DOE
Office of Science User Facility supported by the Office of Science of the U.S.
Department of Energy under Contract No. DE-AC02-05CH11231. JC, SB, and LL
acknowledge support from a SNSF Eccellenza Professorial Fellowship (No.
186879). CB acknowledges support from the U.S. Department of Energy Office of
Science (DE-SC0011784). This work was supported by a grant from the Swiss
National Supercomputing Centre (CSCS) under project ID s1203 and sm80. Some
computations in this paper were run on the Odyssey/Cannon cluster supported by
the FAS Science Division Research Computing Group at Harvard University.

## References

  * Seljak and Zaldarriaga (1997) Uros Seljak and Matias Zaldarriaga, “Signature of gravity waves in polarization of the microwave background,” Phys. Rev. Lett. 78, 2054–2057 (1997), arXiv:astro-ph/9609169 .
  * Seljak (1997) Uros Seljak, “Measuring polarization in cosmic microwave background,” Astrophys. J. 482, 6 (1997), arXiv:astro-ph/9608131 .
  * Kamionkowski  _et al._ (1997) Marc Kamionkowski, Arthur Kosowsky, and Albert Stebbins, “Statistics of cosmic microwave background polarization,” Phys. Rev. D 55, 7368–7388 (1997), arXiv:astro-ph/9611125 .
  * Akrami  _et al._ (2020) Y. Akrami _et al._ (Planck), “Planck 2018 results. X. Constraints on inflation,” Astron. Astrophys. 641, A10 (2020), arXiv:1807.06211 [astro-ph.CO] .
  * Kamionkowski and Kovetz (2016) Marc Kamionkowski and Ely D. Kovetz, “The Quest for B Modes from Inflationary Gravitational Waves,” Ann. Rev. Astron. Astrophys. 54, 227–269 (2016), arXiv:1510.06042 [astro-ph.CO] .
  * Achúcarro  _et al._ (2022) Ana Achúcarro _et al._ , “Inflation: Theory and Observations,” (2022), arXiv:2203.08128 [astro-ph.CO] .
  * Lyth and Riotto (1999) David H. Lyth and Antonio Riotto, “Particle physics models of inflation and the cosmological density perturbation,” Phys. Rept. 314, 1–146 (1999), arXiv:hep-ph/9807278 .
  * Baumann and McAllister (2015) Daniel Baumann and Liam McAllister,  _Inflation and String Theory_, Cambridge Monographs on Mathematical Physics (Cambridge University Press, 2015) arXiv:1404.2601 [hep-th] .
  * Hu and White (1997) Wayne Hu and Martin J. White, “A CMB polarization primer,” New Astron. 2, 323 (1997), arXiv:astro-ph/9706147 .
  * Ade  _et al._ (2021) P. A. R. Ade _et al._ (BICEP, Keck), “Improved Constraints on Primordial Gravitational Waves using Planck, WMAP, and BICEP/Keck Observations through the 2018 Observing Season,” Phys. Rev. Lett. 127, 151301 (2021), arXiv:2110.00483 [astro-ph.CO] .
  * Tristram  _et al._ (2022) M. Tristram _et al._ , “Improved limits on the tensor-to-scalar ratio using BICEP and Planck data,” Phys. Rev. D 105, 083524 (2022), arXiv:2112.07961 [astro-ph.CO] .
  * Lewis and Challinor (2006) Antony Lewis and Anthony Challinor, “Weak gravitational lensing of the cmb,” Phys. Rept. 429, 1–65 (2006), arXiv:astro-ph/0601594 [astro-ph] .
  * Zaldarriaga and Seljak (1998) Matias Zaldarriaga and Uros Seljak, “Gravitational lensing effect on cosmic microwave background polarization,” Phys. Rev. D 58, 023003 (1998), arXiv:astro-ph/9803150 .
  * Knox and Song (2002) Lloyd Knox and Yong-Seon Song, “A Limit on the detectability of the energy scale of inflation,” Phys. Rev. Lett. 89, 011303 (2002), arXiv:astro-ph/0202286 .
  * Kesden  _et al._ (2002) Michael Kesden, Asantha Cooray, and Marc Kamionkowski, “Separation of gravitational wave and cosmic shear contributions to cosmic microwave background polarization,” Phys. Rev. Lett. 89, 011304 (2002), arXiv:astro-ph/0202434 .
  * Abazajian  _et al._ (2016) Kevork N. Abazajian _et al._ (CMB-S4), “CMB-S4 Science Book, First Edition,” (2016), arXiv:1610.02743 [astro-ph.CO] .
  * Abazajian  _et al._ (2022) Kevork Abazajian _et al._ (CMB-S4), “CMB-S4: Forecasting Constraints on Primordial Gravitational Waves,” Astrophys. J. 926, 54 (2022), arXiv:2008.12619 [astro-ph.CO] .
  * Challinor  _et al._ (2018) Anthony Challinor _et al._ (CORE), “Exploring cosmic origins with CORE: Gravitational lensing of the CMB,” JCAP 04, 018 (2018), arXiv:1707.02259 [astro-ph.CO] .
  * Baleato Lizancos  _et al._ (2021a) Antón Baleato Lizancos, Anthony Challinor, and Julien Carron, “Limitations of CMB $B$-mode template delensing,” Phys. Rev. D 103, 023518 (2021a), arXiv:2010.14286 [astro-ph.CO] .
  * Carron (2019) Julien Carron, “Optimal constraints on primordial gravitational waves from the lensed CMB,” Phys. Rev. D 99, 043518 (2019), arXiv:1808.10349 [astro-ph.CO] .
  * Manzotti (2018) A. Manzotti, “Future cosmic microwave background delensing with galaxy surveys,” Phys. Rev. D 97, 043527 (2018), arXiv:1710.11038 [astro-ph.CO] .
  * Ade  _et al._ (2019) Peter Ade _et al._ (Simons Observatory), “The Simons Observatory: Science goals and forecasts,” JCAP 02, 056 (2019), arXiv:1808.07445 [astro-ph.CO] .
  * Namikawa  _et al._ (2022) Toshiya Namikawa _et al._ , “Simons Observatory: Constraining inflationary gravitational waves with multitracer B-mode delensing,” Phys. Rev. D 105, 023511 (2022), arXiv:2110.09730 [astro-ph.CO] .
  * Hu and Okamoto (2002) Wayne Hu and Takemi Okamoto, “Mass reconstruction with cmb polarization,” Astrophys. J. 574, 566–574 (2002), arXiv:astro-ph/0111606 .
  * Okamoto and Hu (2003) Takemi Okamoto and Wayne Hu, “CMB lensing reconstruction on the full sky,” Phys. Rev. D67, 083002 (2003), arXiv:astro-ph/0301031 [astro-ph] .
  * Maniyar  _et al._ (2021) Abhishek S. Maniyar, Yacine Ali-Haïmoud, Julien Carron, Antony Lewis, and Mathew S. Madhavacheril, “Quadratic estimators for CMB weak lensing,” Phys. Rev. D 103, 083524 (2021), arXiv:2101.12193 [astro-ph.CO] .
  * Hirata and Seljak (2003a) Christopher M. Hirata and Uros Seljak, “Reconstruction of lensing from the cosmic microwave background polarization,” Phys. Rev. D68, 083002 (2003a), arXiv:astro-ph/0306354 [astro-ph] .
  * Hirata and Seljak (2003b) Christopher M. Hirata and Uros Seljak, “Analyzing weak lensing of the cosmic microwave background using the likelihood function,” Phys. Rev. D67, 043001 (2003b), arXiv:astro-ph/0209489 [astro-ph] .
  * Carron and Lewis (2017) Julien Carron and Antony Lewis, “Maximum a posteriori CMB lensing reconstruction,” Phys. Rev. D 96, 063510 (2017), arXiv:1704.08230 [astro-ph.CO] .
  * Adachi  _et al._ (2020) S. Adachi _et al._ (POLARBEAR), “Internal delensing of Cosmic Microwave Background polarization $B$-modes with the POLARBEAR experiment,” Phys. Rev. Lett. 124, 131301 (2020), arXiv:1909.13832 [astro-ph.CO] .
  * Legrand and Carron (2022) Louis Legrand and Julien Carron, “Lensing power spectrum of the cosmic microwave background with deep polarization experiments,” Phys. Rev. D 105, 123519 (2022), arXiv:2112.05764 [astro-ph.CO] .
  * Legrand and Carron (2023) Louis Legrand and Julien Carron, “Robust and efficient CMB lensing power spectrum from polarization surveys,” (2023), arXiv:2304.02584 [astro-ph.CO] .
  * Aurlien  _et al._ (2022) Ragnhild Aurlien _et al._ , “Foreground Separation and Constraints on Primordial Gravitational Waves with the PICO Space Mission,” (2022), arXiv:2211.14342 [astro-ph.CO] .
  * Hanany  _et al._ (2019) Shaul Hanany _et al._ (NASA PICO), “PICO: Probe of Inflation and Cosmic Origins,” (2019), arXiv:1902.10541 [astro-ph.IM] .
  * Osborne  _et al._ (2014) Stephen J. Osborne, Duncan Hanson, and Olivier Doré, “Extragalactic Foreground Contamination in Temperature-based CMB Lens Reconstruction,” JCAP 03, 024 (2014), arXiv:1310.7547 [astro-ph.CO] .
  * van Engelen  _et al._ (2014) A. van Engelen, S. Bhattacharya, N. Sehgal, G. P. Holder, O. Zahn, and D. Nagai, “CMB Lensing Power Spectrum Biases from Galaxies and Clusters using High-angular Resolution Temperature Maps,” Astrophys. J. 786, 13 (2014), arXiv:1310.7023 [astro-ph.CO] .
  * Ferraro and Hill (2018) Simone Ferraro and J. Colin Hill, “Bias to CMB Lensing Reconstruction from Temperature Anisotropies due to Large-Scale Galaxy Motions,” Phys. Rev. D97, 023512 (2018), arXiv:1705.06751 [astro-ph.CO] .
  * Darwish  _et al._ (2023) Omar Darwish, Blake D. Sherwin, Noah Sailer, Emmanuel Schaan, and Simone Ferraro, “Optimizing foreground mitigation for CMB lensing with combined multifrequency and geometric methods,” Phys. Rev. D 107, 043519 (2023), arXiv:2111.00462 [astro-ph.CO] .
  * Baleato Lizancos and Ferraro (2022) Antón Baleato Lizancos and Simone Ferraro, “Impact of extragalactic foregrounds on internal delensing of the CMB B-mode polarization,” Phys. Rev. D 106, 063534 (2022), arXiv:2205.09000 [astro-ph.CO] .
  * Omori (2022) Yuuki Omori, “Agora: Multi-Component Simulation for Cross-Survey Science,” (2022), arXiv:2212.07420 [astro-ph.CO] .
  * MacCrann  _et al._ (2023) Niall MacCrann _et al._ (ACT), “The Atacama Cosmology Telescope: Mitigating the impact of extragalactic foregrounds for the DR6 CMB lensing analysis,” (2023), arXiv:2304.05196 [astro-ph.CO] .
  * Anderes  _et al._ (2015) Ethan Anderes, Benjamin Wandelt, and Guilhem Lavaux, “Bayesian inference of CMB gravitational lensing,” Astrophys. J. 808, 152 (2015), arXiv:1412.4079 [astro-ph.CO] .
  * Millea  _et al._ (2020) Marius Millea, Ethan Anderes, and Benjamin D. Wandelt, “Sampling-based inference of the primordial CMB and gravitational lensing,” Phys. Rev. D 102, 123542 (2020), arXiv:2002.00965 [astro-ph.CO] .
  * Millea  _et al._ (2021) M. Millea _et al._ , “Optimal Cosmic Microwave Background Lensing Reconstruction and Parameter Estimation with SPTpol Data,” Astrophys. J. 922, 259 (2021), arXiv:2012.01709 [astro-ph.CO] .
  * Millea  _et al._ (2019) Marius Millea, Ethan Anderes, and Benjamin D. Wandelt, “Bayesian delensing of CMB temperature and polarization,” Phys. Rev. D 100, 023509 (2019), arXiv:1708.06753 [astro-ph.CO] .
  * Millea and Seljak (2022) Marius Millea and Uros Seljak, “Marginal unbiased score expansion and application to CMB lensing,” Phys. Rev. D 105, 103531 (2022), arXiv:2112.09354 [astro-ph.CO] .
  * Bianchini and Millea (2023) Federico Bianchini and Marius Millea, “Inference of gravitational lensing and patchy reionization with future CMB data,” Phys. Rev. D 107, 043521 (2023), arXiv:2210.10893 [astro-ph.CO] .
  * Hanson  _et al._ (2010) Duncan Hanson, Anthony Challinor, and Antony Lewis, “Weak lensing of the CMB,” Gen. Rel. Grav. 42, 2197–2218 (2010), arXiv:0911.0612 [astro-ph.CO] .
  * Hanson  _et al._ (2011) Duncan Hanson, Anthony Challinor, George Efstathiou, and Pawel Bielewicz, “CMB temperature lensing power reconstruction,” Phys. Rev. D 83, 043005 (2011), arXiv:1008.4403 [astro-ph.CO] .
  * Fabbian  _et al._ (2019) Giulio Fabbian, Antony Lewis, and Dominic Beck, “CMB lensing reconstruction biases in cross-correlation with large-scale structure probes,” JCAP 10, 057 (2019), arXiv:1906.08760 [astro-ph.CO] .
  * Hanson  _et al._ (2009) Duncan Hanson, Graca Rocha, and Krzysztof Gorski, “Lensing reconstruction from PLANCK sky maps: inhomogeneous noise,” Mon. Not. Roy. Astron. Soc. 400, 2169–2173 (2009), arXiv:0907.1927 [astro-ph.CO] .
  * Benoit-Levy  _et al._ (2013) Aurelien Benoit-Levy, Typhaine Dechelette, Karim Benabed, Jean-Francois Cardoso, Duncan Hanson, and Simon Prunet, “Full-sky CMB lensing reconstruction in presence of sky-cuts,” Astron. Astrophys. 555, A37 (2013), arXiv:1301.4145 [astro-ph.CO] .
  * Aghanim  _et al._ (2020a) N. Aghanim _et al._ (Planck), “Planck 2018 results. VIII. Gravitational lensing,” Astron. Astrophys. 641, A8 (2020a), arXiv:1807.06210 [astro-ph.CO] .
  * Nocedal (1980) Jorge Nocedal, “Updating quasi-newton matrices with limited storage,” Mathematics of Computation 35, 773–782 (1980).
  * Lewis  _et al._ (2017) Antony Lewis, Alex Hall, and Anthony Challinor, “Emission-angle and polarization-rotation effects in the lensed CMB,” JCAP 08, 023 (2017), arXiv:1706.02673 [astro-ph.CO] .
  * Lewis (2005) Antony Lewis, “Lensed CMB simulation and parameter estimation,” Phys. Rev. D71, 083008 (2005), arXiv:astro-ph/0502469 [astro-ph] .
  * Lavaux and Wandelt (2010) Guilhem Lavaux and Benjamin D. Wandelt, “Fast CMB lensing using statistical interpolation on the sphere,” Astrophys. J. Suppl. 191, 32–42 (2010), arXiv:1003.4984 [astro-ph.CO] .
  * Challinor and Chon (2002) Anthony Challinor and Gayoung Chon, “Geometry of weak lensing of CMB polarization,” Phys. Rev. D 66, 127301 (2002), arXiv:astro-ph/0301064 .
  * Aghanim  _et al._ (2020b) N. Aghanim _et al._ (Planck), “Planck 2018 results. VIII. Gravitational lensing,” Astron. Astrophys. 641, A8 (2020b), arXiv:1807.06210 [astro-ph.CO] .
  * Carron  _et al._ (2022) Julien Carron, Mark Mirmelstein, and Antony Lewis, “CMB lensing from Planck PR4 maps,” JCAP 09, 039 (2022), arXiv:2206.07773 [astro-ph.CO] .
  * Qu  _et al._ (2023) Frank J. Qu _et al._ (ACT), “The Atacama Cosmology Telescope: A Measurement of the DR6 CMB Lensing Power Spectrum and its Implications for Structure Growth,” (2023), arXiv:2304.05202 [astro-ph.CO] .
  * Teng  _et al._ (2011) Wei-Hsiang Teng, Chao-Lin Kuo, and Jiun-Huei Proty Wu, “Cosmic Microwave Background Delensing Revisited: Residual Biases and a Simple Fix,” (2011), arXiv:1102.5729 [astro-ph.CO] .
  * Carron  _et al._ (2017) Julien Carron, Antony Lewis, and Anthony Challinor, “Internal delensing of Planck CMB temperature and polarization,” JCAP 1705, 035 (2017), arXiv:1701.01712 [astro-ph.CO] .
  * Namikawa (2017) Toshiya Namikawa, “CMB internal delensing with general optimal estimator for higher-order correlations,” Phys. Rev. D 95, 103514 (2017), arXiv:1703.00169 [astro-ph.CO] .
  * Baleato Lizancos  _et al._ (2021b) Antón Baleato Lizancos, Anthony Challinor, and Julien Carron, “Impact of internal-delensing biases on searches for primordial $\mathsf{B}$-modes of CMB polarisation,” JCAP 03, 016 (2021b), arXiv:2007.01622 [astro-ph.CO] .
  * Hager (1989) William W. Hager, “Updating the inverse of a matrix,” SIAM Review 31, 221–239 (1989), https://doi.org/10.1137/1031049 .
  * Smith  _et al._ (2007) Kendrick M. Smith, Oliver Zahn, and Olivier Dore, “Detection of Gravitational Lensing in the Cosmic Microwave Background,” Phys. Rev. D76, 043510 (2007), arXiv:0705.3980 [astro-ph] .
  * Smith  _et al._ (2012) Kendrick M. Smith, Duncan Hanson, Marilena LoVerde, Christopher M. Hirata, and Oliver Zahn, “Delensing CMB Polarization with External Datasets,” JCAP 06, 014 (2012), arXiv:1010.0048 [astro-ph.CO] .
  * Aghanim  _et al._ (2014) N. Aghanim _et al._ (Planck), “Planck 2013 results. XXVII. Doppler boosting of the CMB: Eppur si muove,” Astron. Astrophys. 571, A27 (2014), arXiv:1303.5087 [astro-ph.CO] .
  * Cooray  _et al._ (2005) Asantha Cooray, Marc Kamionkowski, and Robert R. Caldwell, “Cosmic shear of the microwave background: The Curl diagnostic,” Phys. Rev. D 71, 123527 (2005), arXiv:astro-ph/0503002 .
  * Pratten and Lewis (2016) Geraint Pratten and Antony Lewis, “Impact of post-Born lensing on the CMB,” JCAP 1608, 047 (2016), arXiv:1605.05662 [astro-ph.CO] .
  * Robertson and Lewis (2023) Mathew Robertson and Antony Lewis, “How to detect lensing rotation,” JCAP 08, 048 (2023), arXiv:2303.13313 [astro-ph.CO] .
  * Namikawa  _et al._ (2012) Toshiya Namikawa, Daisuke Yamauchi, and Atsushi Taruya, “Full-sky lensing reconstruction of gradient and curl modes from CMB maps,” JCAP 01, 007 (2012), arXiv:1110.1718 [astro-ph.CO] .
  * Adam  _et al._ (2016) R. Adam _et al._ (Planck), “Planck intermediate results. XXX. The angular power spectrum of polarized dust emission at intermediate and high Galactic latitudes,” Astron. Astrophys. 586, A133 (2016), arXiv:1409.5738 [astro-ph.CO] .
  * Vansyngel  _et al._ (2017) Flavien Vansyngel, François Boulanger, Tuhin Ghosh, Benjamin D. Wandelt, Jonathan Aumont, Andrea Bracco, François Levrier, Peter G. Martin, and Ludovic Montier, “Statistical simulations of the dust foreground to cosmic microwave background polarization,” Astron. Astrophys. 603, A62 (2017), [Erratum: Astron. Astrophys.618,C4(2018)], arXiv:1611.02577 [astro-ph.CO] .
  * Tegmark  _et al._ (2003) Max Tegmark, Angelica de Oliveira-Costa, and Andrew Hamilton, “A high resolution foreground cleaned CMB map from WMAP,” Phys. Rev. D 68, 123523 (2003), arXiv:astro-ph/0302496 .
  * Górski  _et al._ (2005) K. M. Górski, E. Hivon, A. J. Banday, B. D. Wand elt, F. K. Hansen, M. Reinecke, and M. Bartelmann, “HEALPix: A Framework for High-Resolution Discretization and Fast Analysis of Data Distributed on the Sphere,” Astrophys. J. 622, 759–771 (2005), arXiv:astro-ph/0409513 [astro-ph] .
  * Chon  _et al._ (2004) Gayoung Chon, Anthony Challinor, Simon Prunet, Eric Hivon, and Istvan Szapudi, “Fast estimation of polarization power spectra using correlation functions,” Mon. Not. Roy. Astron. Soc. 350, 914 (2004), arXiv:astro-ph/0303414 .
  * Mirmelstein  _et al._ (2019) Mark Mirmelstein, Julien Carron, and Antony Lewis, “Optimal filtering for CMB lensing reconstruction,” Phys. Rev. D 100, 123509 (2019), arXiv:1909.02653 [astro-ph.CO] .
  * Reinecke and Seljebotn (2013) M. Reinecke and D. S. Seljebotn, “Libsharp - spherical harmonic transforms revisited,” Astronomy & Astrophysics 554, A112 (2013).
  * Reinecke  _et al._ (2023) Martin Reinecke, Sebastian Belkner, and Julien Carron, “Improved CMB (de-)lensing using general spherical harmonic transforms,” (2023), arXiv:2304.10431 [astro-ph.CO] .
  * BICEP2/Keck Collaboration  _et al._ (2015) BICEP2/Keck Collaboration, Planck Collaboration, _et al._ , “Joint Analysis of BICEP2/Keck Array and Planck Data,” Phys. Rev. Lett. 114, 101301 (2015), arXiv:1502.00612 [astro-ph.CO] .
  * BICEP/Keck Collaboration  _et al._ (2021) BICEP/Keck Collaboration, SPTpol Collaboration, _et al._ , “A demonstration of improved constraints on primordial gravitational waves with delensing,” Phys. Rev. D 103, 022004 (2021), arXiv:2011.08163 [astro-ph.CO] .
  * Hamimeche and Lewis (2008) Samira Hamimeche and Antony Lewis, “Likelihood Analysis of CMB Temperature and Polarization Power Spectra,” Phys. Rev. D 77, 103013 (2008), arXiv:0801.0554 [astro-ph] .
  * Ade  _et al._ (2018) P. A. R. Ade _et al._ (BICEP2, Keck Array), “BICEP2 / Keck Array X: Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season,” Phys. Rev. Lett. 121, 221301 (2018), arXiv:1810.05216 [astro-ph.CO] .
  * Abazajian  _et al._ (2019) Kevork Abazajian _et al._ , “CMB-S4 Science Case, Reference Design, and Project Plan,” (2019), arXiv:1907.04473 [astro-ph.IM] .
  * Beck  _et al._ (2020) Dominic Beck, Josquin Errard, and Radek Stompor, “Impact of Polarized Galactic Foreground Emission on CMB Lensing Reconstruction and Delensing of B-Modes,” JCAP 06, 030 (2020), arXiv:2001.02641 [astro-ph.CO] .
  * Newman and Penrose (1966) E. T. Newman and R. Penrose, “Note on the Bondi-Metzner-Sachs group,” J. Math. Phys. 7, 863–870 (1966).
  * Goldberg  _et al._ (1967) J. N. Goldberg, A. J. MacFarlane, E. T. Newman, F. Rohrlich, and E. C. G. Sudarshan, “Spin s spherical harmonics and edth,” J. Math. Phys. 8, 2155 (1967).
  * Lewis  _et al._ (2002) Antony Lewis, Anthony Challinor, and Neil Turok, “Analysis of CMB polarization on an incomplete sky,” Phys. Rev. D 65, 023505 (2002), arXiv:astro-ph/0106536 .

##  Appendix A Spin-weighted fields and harmonics

Throughout this paper, we make heavy use of spin-weighted spherical harmonic
Newman and Penrose (1966); Goldberg  _et al._ (1967); Lewis  _et al._ (2002)
transforms with spins between $0$ and $3$. A spin-weighted field
$\prescript{}{s}{f}(\hat{n})$ ($s\geq 0)$ on the sphere and parametrized with
colatitude and longitude $(\theta,\varphi)$ is defined with reference to the
local axes, $\bm{e}_{\theta}$ and $\bm{e}_{\varphi}$. Our conventions are that
$\bm{e}_{\theta}$ and $\bm{e}_{\varphi}$ point southwards and eastwards,
respectively.

Under a clockwise rotation at $\hat{n}$ by some angle $\psi$, the field
transforms, by definition, as

$\prescript{}{s}{f}(\hat{n})\rightarrow\>e^{is\psi(\hat{n})}\>\prescript{}{s}{f}(\hat{n}).$
(63)

A spin-$s$ field can conveniently be decomposed into its gradient $(G)$ and
curl $(C)$ harmonics by using the spin-weighted spherical harmonics,
${}_{s}Y_{\ell m}(\hat{n})$. The expansion has the form

$\prescript{}{\pm s}{f}(\hat{n})=-(\pm 1)^{s}\sum_{\ell m}\left(G_{\ell m}\pm
iC_{\ell m}\right)\>\prescript{}{\pm s}{Y_{\ell m}(\hat{n})}.$ (64)

Further, let $\partial^{+}$ and $\partial^{-}$ be the spin-raising and spin-
lowering operators Newman and Penrose (1966); Goldberg  _et al._ (1967),

$\partial^{\pm}_{s}f=-\left(\sin\theta\right)^{\pm
s}\left(\partial_{\theta}\pm\frac{i}{\sin\theta}\partial_{\varphi}\right)\left(\sin\theta\right)^{\mp
s}\>_{s}f\,.$ (65)

We follow generally the conventions of Ref. Lewis  _et al._ (2002), and refer
to the appendices of that paper for a discussion in the context of the CMB.
With the help of these operators, one can introduce orthonormal basis
functions (the spin-weighted spherical harmonics ${}_{s}Y_{\ell m}$) for spin-
weighted fields. Explicitly,

$\begin{split}\partial^{\pm}_{s}{Y_{\ell m}}&=\pm\sqrt{(\ell\mp s)(\ell\pm
s+1)}\>\prescript{}{s\pm 1}{Y_{\ell m}}.\end{split}$ (66)

A lensing remapping on the sphere can be described by a vector field
(concretely, a field of spin-$1$), and this vector field can, in turn, be
described by two scalar (spin-$0$) lensing potentials, $\phi(\hat{n})$ and
$\Omega(\hat{n})$, with the spin-$1$ deflection vector defined as follows:

$\begin{split}{}_{1}\alpha(\hat{n})&=-\eth^{+}\left(\phi(\hat{n})+i\>\Omega(\hat{n})\right)\\\
&=-\sum_{\ell m}\sqrt{\ell(\ell+1)}\left(\phi_{\ell m}+i\Omega_{\ell
m}\right)\>\prescript{}{1}{Y_{\ell m}}(\hat{n})\,.\end{split}$ (67)

The second line follows from Eq. (66). As can be seen by comparing this to Eq.
(64), $\phi$ defines the gradient and $\Omega$ the curl component of the
deflection vector. For a pure gradient deflection field, its real and
imaginary components at $\hat{n}$ are simply $\partial_{\theta}\phi(\hat{n})$
and $\csc(\theta)\partial_{\varphi}\phi(\hat{n})$, and are generally denoted
$\alpha_{\theta}$ and $\alpha_{\varphi}$. If convenient, we also make use of
the polar representation of the vector, with $\alpha$ the magnitude of the
field, and $\beta$ the angle to the $\theta$ coordinate direction, such that,
at each point on the sphere, the following relation holds:

$\begin{split}{}_{\pm 1}\alpha(\hat{n})\equiv\alpha_{\theta}(\hat{n})\pm
i\alpha_{\varphi}(\hat{n})\equiv\alpha(\hat{n})e^{\pm
i\beta(\hat{n})}.\end{split}$ (68)

##  Appendix B Delensed-noise mean-field

This section gives details on the delensed-noise mean-field (the variation of
the covariance matrix log-determinant with $\bm{\alpha}$) in the case of an
otherwise isotropic configuration. Using the Woodbury determinant relation,
the relevant log-determinant is

$\frac{1}{2}\ln\det\left[C^{EE,\rm
unl,-1}+N_{\bm{\alpha}}^{-1}\right]\equiv\frac{1}{2}\ln D\,,$ (69)

with the inverse delensed $E$-mode noise matrix

$N_{\bm{\alpha}}^{-1}=\,_{2}\mathcal{Y}^{\dagger}\mathcal{D}_{\bm{\alpha}}^{\dagger}\mathcal{B}^{\dagger}N^{-1}\mathcal{B}\mathcal{D}_{\bm{\alpha}}\,{}_{2}\mathcal{Y}.$
(70)

There are two variations to be considered. The first variation of the log-
determinant can be written in the following way, using simple matrix algebra,

$\begin{split}\delta\ln D&=\textrm{Tr}\>C^{EE,\rm unl}(C^{EE,\rm
unl}+N_{\bm{\alpha}})^{-1}\delta\ln N_{\bm{\alpha}}^{-1}\\\
&\equiv\textrm{Tr}\>\mathcal{W}_{\bm{\alpha}}\delta\ln
N_{\bm{\alpha}}^{-1}\,.\end{split}$ (71)

For convenience, we have used the notation $\delta\ln N_{\bm{\alpha}}^{-1}$
for the matrix $N_{\bm{\alpha}}\delta N_{\bm{\alpha}}^{-1}$. This matrix has
roughly constant contributions on all scales. It is contracted with the
Wiener-filtering matrix (defined as $\mathcal{W}_{\bm{\alpha}}$ in the second
line), thereby suppressing contributions from noisy $E$-modes. The impact of
$\bm{\alpha}$ in the Wiener-filter matrix is, however, only significant when
the noise is comparable to $C^{EE,\rm unl}$. Hence, for low noise experiments,
the isotropic Wiener filter $\mathcal{W}_{\ell}=C_{\ell}/(C_{\ell}+N_{\ell})$,
which is the leading-order contribution, already gives an almost exact result.
We now discuss $\delta\ln N_{\bm{\alpha}}^{-1}$ in more detail, and then look
at the second variation as well.

The action of lenses much larger than the beam commutes with the beam to good
accuracy. It is useful to consider this approximation first, which results in
the correct low-$L$ behavior, and which also renders the results more
transparent. This assumption also corresponds to the exact result in the case
of zero beamwidth (but arbitrary large white noise level, so that this still
corresponds to a somewhat physically relevant situation). Under this
assumption, we can swap the order of appearance of the beam and deflection
operations in Eq. (70). Assuming homogeneous noise with full sky coverage, and
using the exact relation
$[\mathcal{D}_{\bm{\alpha}}^{\dagger}\mathcal{D}_{\bm{\alpha}}][\hat{n},\hat{n}^{\prime}]=|A_{\bm{\alpha}^{-1}}|(\hat{n})\delta^{D}(\hat{n},\hat{n}^{\prime})$,
the noise matrix greatly simplifies.

For convenience, we now adopt notation convenient for intensity measurements
by using spin-$0$ spherical harmonics. However, the simulations discussed in
this paper include polarized data and have no primordial $B$-modes, and we
generally must use spin-$2$ harmonics to correctly describe it. These two
differences lead to slight complications in the argument below. They are,
nevertheless, only formal, and lead to identical leading-order and similar
second-order results, and will be discussed at the end.

The inverse noise matrix, in the approximation discussed above, becomes

$\left[N_{\bm{\alpha}}^{-1}\right]_{\ell m,\ell^{\prime}m^{\prime}}=\int
d\hat{n}\frac{|A_{\bm{\alpha}^{-1}}|(\hat{n})}{N_{\ell}^{1/2}N_{\ell^{\prime}}^{1/2}}Y^{\dagger}_{\ell
m}(\hat{n})Y_{\ell^{\prime}m^{\prime}}(\hat{n})\,.$ (72)

Assuming that the bandlimits give a sufficiently large range, the spherical
harmonics project onto a complete set, and the matrix has the explicit inverse

$\left[N_{\bm{\alpha}}\right]_{\ell m,\ell^{\prime}m^{\prime}}=\int
d^{2}\hat{n}\frac{N_{\ell}^{1/2}N_{\ell^{\prime}}^{1/2}}{|A_{\bm{\alpha}^{-1}}|(\hat{n})}Y^{\dagger}_{\ell
m}(\hat{n})Y_{\ell^{\prime}m^{\prime}}(\hat{n})\,.$ (73)

With this, we find that the variation $N_{\bm{\alpha}}\delta
N_{\bm{\alpha}}^{-1}$ in Eq. (71) is

$\begin{split}&\left[\delta\ln N_{\bm{\alpha}}^{-1}\right]_{\ell
m,\ell^{\prime}m^{\prime}}\\\
&=\frac{N^{1/2}_{\ell}}{N^{1/2}_{\ell^{\prime}}}\int
d^{2}\hat{n}\>\delta\ln|A_{\bm{\alpha}^{-1}}|(\hat{n})Y^{\dagger}_{\ell
m}(\hat{n})Y_{\ell^{\prime}m^{\prime}}(\hat{n}).\end{split}$ (74)

Further, we can then write the leading term of Eq. (71) as

$\delta\ln D=\left(\delta\int\frac{d^{2}\hat{n}^{\prime}}{4\pi}\>\ln
A_{\bm{\alpha}^{-1}}(\hat{n}^{\prime})\right)\left(\sum_{\ell}(2\ell+1)\mathcal{W}_{\ell}\right)\,.$
(75)

By using $\ln A_{\bm{\alpha}^{-1}}(\hat{n}^{\prime})=-\ln
A_{\bm{\alpha}}(\hat{n})$, where $\hat{n}^{\prime}$ and $\hat{n}$ are the
deflected and undeflected positions, respectively, and performing the change
of parameterization in the integral, the prefactor is

$-\delta\int\frac{d^{2}\hat{n}}{4\pi}|A_{\bm{\alpha}}|(\hat{n})\ln|A_{\bm{\alpha}}|(\hat{n})\sim-\delta\int\frac{d^{2}\hat{n}}{4\pi}2\kappa^{2}(\hat{n})\,.$
(76)

By definition, $g^{\rm MF,\kappa}_{LM}=\frac{1}{2}\delta\ln
D/\delta\kappa_{LM}$, and we recover

$g^{\rm
MF,\kappa}_{LM}=-2\kappa_{LM}\left(\sum_{\ell}\frac{2\ell+1}{4\pi}\mathcal{W}_{\ell}\right)\,,$
(77)

as given in the text in Eq. (24).

Let us now discuss the second variation. It is natural to continue working
with $\ln A_{\bm{\alpha}^{-1}}$, since the second variation $\delta^{2}\ln
N_{\bm{\alpha}}^{-1}$ vanishes (see Eq. LABEL:eq:lnMvar). For this reason, the
only contribution is from the change of the Wiener filter. This may be written
as

$\delta^{2}\ln
D=\left[\textrm{Tr}\>\mathcal{\mathcal{W}_{\bm{\alpha}}}\delta\ln
N_{\bm{\alpha}}^{-1}(1-\mathcal{\mathcal{W}_{\bm{\alpha}}})\delta\ln
N_{\bm{\alpha}}^{-1}\right]\,.$ (78)

The full result following from this is given in Eq. (80). The squeezed limit,
however, (low $L$, high $\ell$, and $\ell\sim\ell^{\prime}$) is easily
extracted:

$\left(\int\frac{d^{2}\hat{n}}{4\pi}|A_{\bm{\alpha}}|(\ln|A_{\bm{\alpha}}|)^{2}\right)\left(\sum_{\ell}(2\ell+1)\mathcal{W}_{\ell}(1-\mathcal{W}_{\ell})\right)\,.$
(79)

The prefactor is $\int(2\kappa)^{2}+$ higher orders. This is of the same order
in $\kappa$ as the first variation, but with much lower relevance due to
$\mathcal{W}_{\ell}(1-\mathcal{W}_{\ell})$ and taking only contributions when
$C_{\ell}\sim N_{\ell}$. The exact second variation result for spin-$0$
intensity is

$\begin{split}\frac{1}{2}\ln D\ni&\frac{1}{2}\sum_{LM}|2\kappa_{LM}|^{2}\\\
&\cdot
2\pi\int_{-1}^{1}d\mu\>d^{L}_{00}(\mu)\xi_{00}^{\mathcal{W}}(\mu)\xi_{00}^{\mathcal{(}1-\mathcal{W})}(\mu)\,,\end{split}$
(80)

with, for any spins $a$, $b$, and filter $F$,

$\xi_{ab}^{F}(\mu)=\sum_{\ell}\frac{2\ell+1}{4\pi}F_{\ell}\>d^{\ell}_{ab}(\mu).$
(81)

For polarization, and in the limit of vanishing $C_{\ell}^{BB,\rm unl}$, the
differences are the following. First, the Wiener-filter
$\mathcal{W}_{\bm{\alpha}}$ becomes pure $EE$ to leading order, but
$(1-\mathcal{W}_{\bm{\alpha}})$ has both an $EE$ and a $BB$ component. This
last term originates from the contribution of the $EB$ part of the quadratic
estimator. Second, we must use the spin-2 harmonics. The first variation is
unchanged, but the second becomes

$\begin{split}&\frac{1}{2}\ln
D\ni\frac{1}{2}\sum_{LM}|2\kappa_{LM}|^{2}\times\\\
&\left\\{2\pi\int_{-1}^{1}\>d^{L}_{00}(\mu)\frac{1}{2}\left[\xi_{2,2}^{\mathcal{W}}(\mu)\xi_{-2,-2}^{1-\mathcal{W}}(\mu)+\xi^{\mathcal{W}}_{-2,2}(\mu)\xi^{1-\mathcal{W}}_{2,-2}(\mu)\right]\right.\\\
&\left.+2\pi\int_{-1}^{1}\>d^{L}_{00}(\mu)\frac{1}{2}\left[\xi_{2,2}^{\mathcal{W}}(\mu)\xi_{-2,-2}^{1}(\mu)-\xi^{\mathcal{W}}_{-2,2}(\mu)\xi^{1}_{2,-2}(\mu)\right]\right\\},\\\
\end{split}$

where the second line is the $EB$ contribution. Combining this with Eq. (77)
gives the final result for the mean-field for tiny beams,

$\begin{split}&g^{\rm
MF,\kappa}_{LM}=-2\kappa_{LM}\left\\{\left[\sum_{\ell}\frac{2\ell+1}{4\pi}\mathcal{W}_{\ell}\right]\right.\\\
&-2\pi\int_{-1}^{1}\>d^{L}_{00}(\mu)\frac{1}{2}\left[\xi_{2,2}^{\mathcal{W}}(\mu)\xi_{-2,-2}^{1-\mathcal{W}}(\mu)+\xi^{\mathcal{W}}_{-2,2}(\mu)\xi^{1-\mathcal{W}}_{2,-2}(\mu)\right]\\\
&\left.-2\pi\int_{-1}^{1}\>d^{L}_{00}(\mu)\frac{1}{2}\left[\xi_{2,2}^{\mathcal{W}}(\mu)\xi_{-2,-2}^{1}(\mu)-\xi^{\mathcal{W}}_{-2,2}(\mu)\xi^{1}_{2,-2}(\mu)\right]\right\\}\,.\\\
\end{split}$ (82)

The first two lines are the leading and sub-leading $EE$ contributions, while
the third line is the one from the $EB$-part of the quadratic estimator, which
is basically zero at low $L$, as expected.

##  Appendix C Curved-sky likelihood gradients

In this section, we give details on the derivation of the curved-sky
likelihood gradients, Eq. (18). For this, we consider a small variation
$\bm{\epsilon}$ to the deflection vector $\bm{\alpha}(\hat{n})$ in the plane
tangent to $\hat{n}$ in Fig. 2. The deflection vector $\bm{\alpha}$ gives the
coordinates of the deflected position $\hat{n}^{\prime}$. If we were to
slightly change the deflection vector, for example to
$\bm{\alpha}+\bm{\epsilon}$, this would lead us to a position
$\hat{n}^{\prime}_{\epsilon}$. We can picture this as a squeezed spherical
triangle, with two large sides, one of length $|\bm{\alpha}|$ joining
$\hat{n}$ to $\hat{n}^{\prime}$, the other of length
$|\bm{\alpha}(\hat{n})+\bm{\epsilon}(\hat{n})|$ joining $\hat{n}$ to
$\hat{n}^{\prime}_{\epsilon}$. The remaining, squeezed side of the triangle
joining $\hat{n}^{\prime}$ to $\hat{n}^{\prime}_{\epsilon}$ lies, for
practical purposes, in the plane tangent to $\hat{n}^{\prime}$ and can be
treated to first order in $\bm{\epsilon}$.

Our strategy is then as follows. Working in spin-weight components, we can
associate to the change in deflection $\>{}_{1}\epsilon(\hat{n})$ at $\hat{n}$
a deflection field $\>{}_{1}\epsilon(\hat{n}^{\prime})$ at $\hat{n}^{\prime}$,
and the remapping from $\hat{n}^{\prime}$ to $\hat{n}^{\prime}_{\epsilon}$
that it induces can be treated perturbatively. Doing so, and using Eq. (2) for
the ‘large’ deflection and Eq. (6) for the small deflection, we can write

$\begin{split}&2\frac{\delta}{\delta_{\pm
1}\alpha(\hat{n})}\left[\mathcal{D}_{\bm{\alpha}}\>\,{}_{s}\mathcal{T}\right](\hat{n})\\\
&=-e^{is(\beta-\beta^{\prime})}\left[\frac{\partial\>_{-1}\epsilon(\hat{n}^{\prime})}{\partial\>_{\pm
1}\epsilon(\hat{n})}\partial^{+}_{s}\mathcal{T}(\hat{n}^{\prime})+\frac{\partial\>_{1}\epsilon(\hat{n}^{\prime})}{\partial\>_{\pm
1}\epsilon(\hat{n})}\partial^{-}_{s}\mathcal{T}(\hat{n}^{\prime})\right]\,.\end{split}$
(83)

In this way, we only need to obtain the linear, explicit dependency of the
components $\>{}_{1}\epsilon(\hat{n}^{\prime})$ to
$\>{}_{1}\epsilon(\hat{n})$.

It is convenient to split
$\bm{\epsilon}(\hat{n})=\bm{\epsilon}_{\parallel}(\hat{n})+\bm{\epsilon}_{\perp}(\hat{n})$
in the components parallel and perpendicular to $\bm{\alpha}(\hat{n})$. Using,
for example, simple spherical trigonometry, it is not too difficult to see
that the parallel component gives a contribution
$\bm{\epsilon}_{\parallel}(\hat{n})$ to the squeezed side of the spherical
triangle, but the perpendicular contribution is reduced by the sky curvature
to $\bm{\epsilon}_{\perp}(\hat{n})\left(\sin\alpha/\alpha\right)$.

The spin-weight components of $\bm{\epsilon}(\hat{n})$ at $\hat{n}$ are given
by

$\>{}_{1}\epsilon(\hat{n})=e^{i\beta}(\epsilon_{\parallel}(\hat{n})+i\epsilon_{\perp}(\hat{n}))\,.$
(84)

At $\hat{n}^{\prime}$, for the reason just discussed above, the components of
the deflection vector pointing to $\hat{n}^{\prime}_{\epsilon}$ become instead

$\>{}_{1}\epsilon(\hat{n}^{\prime})=e^{i\beta^{\prime}}\left(\epsilon_{\parallel}(\hat{n})+i\left(\frac{\sin\alpha}{\alpha}\right)\epsilon_{\perp}(\hat{n})\right)\,.$
(85)

We can rearrange these two relations, and find one piece that is linear in the
deflection angle $\alpha$, and a second piece that is quadratic:

$\begin{split}&{}_{1}\epsilon(\hat{n}^{\prime})=\>e^{-i(\beta-\beta^{\prime})}\>_{1}\epsilon(\hat{n})\\\
&-e^{-i(\beta-\beta^{\prime})}\frac{1}{2}\left(1-\frac{\sin\alpha}{\alpha}\right)\left[\>\>{}_{1}\epsilon(\hat{n})-e^{2i\beta}\>_{-1}\epsilon(\hat{n})\right].\end{split}$
(86)

Plugging this into Eq. (83), and using the definition of the remapping by
$\bm{\alpha}$, Eq. (2), recovers Eq. (18), as given in the text, after
identification of $e^{\pm 2i\beta}$ to $\left({}_{\pm
1}\alpha(\hat{n})/\alpha(\hat{n})\right)^{2}$.

For $\alpha=\pi$, the perpendicular component is reduced to zero, since all
great circles starting at $\hat{n}$ meet at the antipodal point. This provides
a simple sanity check of these expressions: consider varying the deflection
phase $\beta$, then the gradient of the deflected fields must exactly vanish.
Since, by definition, ${}_{\pm 1}\alpha(\hat{n})=\alpha(\hat{n})e^{\pm
i\beta(\hat{n})}$, we may write in general

$\begin{split}&\frac{1}{i}\frac{\delta\left[\mathcal{D}_{\bm{\alpha}}\>\,{}_{s}\mathcal{T}\right](\hat{n})}{\delta\beta(\hat{n})}=\\\
&{}_{1}\alpha(\hat{n})\frac{\delta\left[\mathcal{D}_{\bm{\alpha}}\>\,{}_{s}\mathcal{T}\right](\hat{n})}{\delta_{1}\alpha(\hat{n})}-\,_{-1}\alpha(\hat{n})\frac{\delta\left[\mathcal{D}_{\bm{\alpha}}\>\,{}_{s}\mathcal{T}\right](\hat{n})}{\delta_{-1}\alpha(\hat{n})}\,.\end{split}$
(87)

The spin-weighted gradients (Eq. 18) themselves become

$\begin{split}&\frac{2\delta}{\delta_{\pm
1}\alpha(\hat{n})}\left[\mathcal{D}_{\bm{\alpha}}\>\,{}_{s}\mathcal{T}\right](\hat{n})\\\
&\overbrace{=}^{\alpha(\hat{n})=\pi}-\frac{1}{2}[\mathcal{D}_{\bm{\alpha}}\eth^{\mp}_{s}\mathcal{T}](\hat{n})-\frac{1}{2}\left(\frac{{}_{\mp
1}\alpha(\hat{n})}{\alpha(\hat{n})}\right)^{2}[\mathcal{D}_{\bm{\alpha}}\eth^{\pm}_{s}\mathcal{T}](\hat{n})\,.\end{split}$
(88)

It is not difficult to see that the $\beta$ gradient Eq. (87) always vanishes,
as it should.

##  Appendix D Magnification determinant

In this section, we describe a series of steps to derive the exact form of the
remapping Jacobian determinant in Eq. (15).
In standard polar coordinates, we use the local orthonormal basis at
$\hat{n}$, yielding $\bm{e}_{r}=\hat{n},\bm{e}_{\theta}$, and
$\bm{e}_{\varphi}$. Variations along $\bm{e}_{\theta}$ and $\bm{e}_{\varphi}$
induce at the deflected position $\hat{n}^{\prime}$ shifts
$\partial_{\theta}\hat{n}^{\prime}$ and
$\frac{1}{\sin\theta}\partial_{\varphi}\hat{n}^{\prime}$ in the plane tangent
to the sphere at $\hat{n}^{\prime}$. We can write these shifts using
$\bm{e}_{r},\bm{e}_{\theta}$, and $\bm{e}_{\varphi}$.

With this, the Jacobian determinant can be written as

$|A|=\left|\partial_{\theta}\>\hat{n}^{\prime}\times\partial_{\tilde{}\varphi}\>\hat{n}^{\prime}\right|,$
(89)

where we have used the notation $\partial_{\tilde{}\varphi}$ for
$\frac{1}{\sin\theta}\partial_{\varphi}$. The third deflection vector
$\bm{d}=\hat{n}^{\prime}-\hat{n}$ can be written in this basis as

$\hat{n}^{\prime}=\hat{n}+\bm{d}\equiv(1+\tilde{}\alpha_{r})\bm{e}_{r}+\tilde{}\alpha_{\theta}\bm{e}_{\theta}+\tilde{}\alpha_{\varphi}\bm{e}_{\varphi}\,,$
(90)

with

$\begin{split}1+\tilde{}\alpha_{r}=\cos\alpha,\quad\tilde{}\alpha_{\theta}=\frac{\sin\alpha}{\alpha}\alpha_{\theta},\quad\tilde{}\alpha_{\varphi}=\frac{\sin\alpha}{\alpha}\alpha_{\varphi}.\end{split}$
(91)

The variation of $\hat{n}^{\prime}$ can be calculated from the variations of
${}_{1}\tilde{\alpha}$, and those of $\bm{e}_{r},\bm{e}_{\theta}$, and
$\bm{e}_{\varphi}$. A straightforward, albeit somewhat lengthy calculation,
gives

$\begin{split}&\partial_{\theta}\hat{n}^{\prime}\times\partial_{\tilde{}\varphi}\hat{n}^{\prime}\\\
&=\begin{pmatrix}\partial_{\theta}\tilde{}\alpha_{r}-\tilde{}\alpha_{\theta}\\\
1+\tilde{}\alpha_{r}+\partial_{\theta}\tilde{}\alpha_{\theta}\\\
\partial_{\theta}\tilde{}\alpha_{\varphi}\end{pmatrix}\times\begin{pmatrix}\partial_{\tilde{}\varphi}\tilde{}\alpha_{r}-\tilde{}\alpha_{\varphi}\\\
\partial_{\tilde{}\varphi}\tilde{}\alpha_{\theta}-\tilde{}\alpha_{\varphi}\cot\theta\\\
1+\tilde{}\alpha_{r}+\tilde{}\alpha_{\theta}\cot\theta+\partial_{\tilde{}\varphi}\tilde{}\alpha_{\varphi}\end{pmatrix}\,.\end{split}$
(92)

With the definitions of the convergence and field rotation,
$\tilde{\kappa}+i\tilde{\omega}=\frac{1}{2}\tilde{\eth}\>_{1}\tilde{\alpha}$,
shear $\tilde{\gamma}=\frac{1}{2}\eth\>_{1}\tilde{\alpha}$, as well as a new
spin-$1$ field
${}_{1}f=f_{1}+if_{2}\equiv\>_{1}\tilde{\alpha}+\eth\tilde{}\alpha_{r}$, we
may write

$\begin{split}\partial_{\theta}\hat{n}^{\prime}\times\partial_{\tilde{}\varphi}\hat{n}^{\prime}=\begin{pmatrix}(1+\tilde{}\alpha_{r}-\tilde{\kappa})^{2}+\tilde{\omega}^{2}-\tilde{\gamma}^{2}\\\
-f_{1}(1+\tilde{}\alpha_{r}-\tilde{\kappa}+\tilde{\gamma}_{1})+f_{2}(-\tilde{\gamma}_{2}-\tilde{\omega})\\\
-f_{1}(-\tilde{\gamma}_{2}+\tilde{\omega})+f_{2}(1+\tilde{}\alpha_{r}-\tilde{\kappa}-\tilde{\gamma}_{1})\end{pmatrix}\,.\end{split}$
(93)

The first ($\bm{e}_{r}$) component reduces to the usual magnification matrix
determinant $1-2\kappa$ for infinitesimal deflections. The $\bm{e}_{\theta}$
and $\bm{e}_{\varphi}$ components form the real and imaginary parts of a
spin-$1$ field, giving a correction quadratic in the deflection angle.
To simplify all the calculations, we introduce and calculate with spin-$0$
quantities $\tilde{\eta},\eta$ and $\tilde{\xi},\xi$, built from the
3-dimensional and tangential deflection , respectively,

$\begin{split}\tilde{\eta}&\equiv
1+\tilde{}\alpha_{r}-\tilde{\kappa}+i\tilde{\omega},\quad\eta\equiv
1-\kappa+i\omega\,,\\\ \tilde{\xi}&\equiv
e^{-2i\beta}\tilde{\gamma},\quad\xi\equiv e^{-2i\beta}\gamma\,.\end{split}$
(94)

A short calculation, applying $\eth$ on $\tilde{\eta}$ and $\tilde{\xi}$ shows
that they are quite simply related to their tangential counterparts:

$\begin{split}\begin{pmatrix}\tilde{\eta}\\\
\tilde{\xi}\end{pmatrix}=&\begin{pmatrix}j_{0}(\alpha)-\frac{\alpha}{2}j_{1}(\alpha)&\frac{\alpha}{2}j_{1}(\alpha)\\\
\frac{\alpha}{2}j_{1}(\alpha)&j_{0}-\frac{\alpha}{2}j_{1}(\alpha)\end{pmatrix}\begin{pmatrix}\eta\\\
\xi\end{pmatrix}\\\ &-\frac{\alpha}{2}j_{1}(\alpha)\begin{pmatrix}1\\\
1\end{pmatrix}\,,\end{split}$ (95)

with $j_{0}(\alpha)=\sin(\alpha)/\alpha$ and $\alpha
j_{1}(\alpha)=-\cos(\alpha)+\sin(\alpha)/\alpha$ the first two spherical
Bessel functions. With this notation in place, the $\bm{e}_{r}$ component of
the vector product is

$\begin{split}&(1+\tilde{}\alpha_{r}-\tilde{\kappa})^{2}+\tilde{\omega}^{2}-\tilde{\gamma}^{2}=|\tilde{\eta}|^{2}-|\tilde{\xi}|^{2}\\\
&=\cos\alpha\left(j_{0}(\alpha)\left(|\eta|^{2}-|\xi|^{2}\right)-\frac{\alpha}{2}j_{1}(\alpha)(\eta+\eta^{*}-\xi-\xi^{*})\right)\,.\end{split}$
(96)

The contributions of the $\bm{e}_{\theta}$ and $\bm{e}_{\varphi}$ can also be
expressed easily in terms of $\eta$ and $\xi$, with the exact same result as
Eq. (96), but with a $\sin\alpha$ prefactor replaced by $\cos\alpha$. The full
Jacobian determinant then collects a factor $\cos^{2}\alpha$ and
$\sin^{2}\alpha$ from the component parallel and perpendicular to $\hat{n}$,
respectively, with the result

$|A|=\left|j_{0}(\alpha)\left(|\eta|^{2}-|\xi|^{2}\right)-\frac{\alpha}{2}j_{1}(\alpha)(\eta+\eta^{*}-\xi-\xi^{*})\right|,$
(97)

equivalent to the expression given in Eq. (15).

The first line of Eq. (15) looks familiar in the context of lensing, arguably
less so the second. We can illustrate its role with the case of a locally
constant deflection field, for which $\kappa=\omega=\gamma=0$ locally, with
the result $A(\hat{n})=\cos\alpha$. Consider a tiny disk centred at the North
Pole, with a locally constant deflection field there. If we increase the
magnitude of the deflection field there gradually from zero, the circle starts
to move towards the equator. The diameter of the disk parallel to the
deflection vector stays unchanged, since the points are lying on the same
great circle. The diameter transverse to the deflection gets smaller owing to
focussing. Upon reaching the equator ($\alpha=\pi/2$), the disk has become a
line, with the component parallel to the equator being squeezed to zero,
consistent with $A=0$. Increasing $\alpha$ further, at the South Pole
($\alpha=\pi$) the line has becomes a disk again, but with $A=-1$, since the
transverse geodesics did cross each other on the equator.

##  Appendix E Other miscellaneous formulae

Generally, the undeflected angles can easily be obtained from the relation

$\hat{n}^{\prime}=\cos\alpha\>\hat{n}+\frac{\sin\alpha}{\alpha}\left(\alpha_{\theta}\bm{e}_{\theta}+\alpha_{\varphi}\bm{e}_{\varphi}\right)\,.$
(98)

The explicit relations are (Lewis, 2005, for example)

$\begin{split}\cos\theta^{\prime}=-&\alpha_{\theta}\frac{\sin\alpha}{\alpha}\sin\theta+\cos\alpha\cos\theta\,,\\\
\sin\theta^{\prime}\sin\left(\varphi^{\prime}-\varphi\right)=&\alpha_{\varphi}\frac{\sin\alpha}{\alpha}\,,\\\
\sin\theta^{\prime}\cos\left(\varphi^{\prime}-\varphi\right)=&\alpha_{\theta}\frac{\sin\alpha}{\alpha}\cos\theta+\cos\alpha\sin\theta\,.\\\
\end{split}$ (99)

The explicit dependence on the deflection amplitude is always second order and
weak, since $\alpha\simeq 2\textrm{ arcmin}=5.8\times 10^{-4}$. Conversely,
the components of a deflection field can be obtained from the angles using

$\begin{split}\alpha_{\varphi}\left(\frac{\sin\alpha}{\alpha}\right)&=\sin(\varphi^{\prime}-\varphi)\sin\theta^{\prime}\,,\\\
\alpha_{\theta}\left(\frac{\sin\alpha}{\alpha}\right)&=\sin(\theta^{\prime}-\theta)-2\sin^{2}\left(\frac{1}{2}(\varphi^{\prime}-\varphi)\right)\cos\theta\sin\theta^{\prime}\,.\end{split}\
$ (100)

To leading order the phase of the lensing remapping is given by313131Ref.
Challinor and Chon (2002) lacks a factor $\textrm{csc}(\theta)$ in their Eq.
8.

$\beta-\beta^{\prime}=\alpha_{\varphi}\left[\textrm{cot}~{}\theta-\frac{1}{2}\textrm{csc}^{2}\theta(1+\cos^{2}\theta)\alpha_{\theta}+\cdots\right]\,,$
(101)

and can also easily be calculated in general from

$\beta^{\prime}=\arctan_{2}(\alpha_{\varphi},\alpha\sin\alpha\cot\theta+\alpha_{\theta}\cos\alpha)\,.$
(102)