Maximum a posteriori CMB lensing reconstruction

TL;DR

This paper introduces LensIt, an efficient iterative method for maximum a posteriori CMB lensing reconstruction that significantly improves signal-to-noise ratio over quadratic estimators, especially for polarization data.

Contribution

The authors develop and implement LensIt, a novel algorithm that finds the most probable lensing map without significant approximations, enhancing the accuracy of CMB lensing reconstruction.

Findings

LensIt achieves maximal theoretical improvement in signal-to-noise ratio.

Delensing with LensIt reduces tensor-to-scalar ratio errors by about a factor of 2.

The method is efficient, converging quickly with linear computational cost.

Abstract

Gravitational lensing of the CMB is a valuable cosmological signal that correlates to tracers of large-scale structure and acts as a important source of confusion for primordial $B$-mode polarization. State-of-the-art lensing reconstruction analyses use quadratic estimators, which are easily applicable to data. However, these estimators are known to be suboptimal, in particular for polarization, and large improvements are expected to be possible for high signal-to-noise polarization experiments. We develop a method and numerical code, $\rm{LensIt}$, that is able to find efficiently the most probable lensing map, introducing no significant approximations to the lensed CMB likelihood, and applicable to beamed and masked data with inhomogeneous noise. It works by iteratively reconstructing the primordial unlensed CMB using a deflection estimate and its inverse, and removing residual lensing from these maps with quadratic estimator techniques. Roughly linear computational cost is maintained due to fast convergence of iterative searches, combined with the local nature of lensing. The method achieves the maximal improvement in signal to noise expected from analytical considerations on the unmasked parts of the sky. Delensing with this optimal map leads to forecast tensor-to-scalar ratio parameter errors improved by a factor $\simeq 2 $ compared to the quadratic estimator in a CMB stage IV configuration.

Tables

Table 1. Table 1: The first three rows give the temperature and polarization noise levels and beam width input to the simulations used in Sec.  IV.3 to obtain the delensing efficiencies shown on Fig.  4 . The next two rows show the effective B 𝐵 B -mode lensing power achievable on the scales relevant for primordial B 𝐵 B -mode measurement, using the quadratic (MV) estimator and the iterated solution respectively. The latter results accurately match predictions using an iterated Gaussian noise level as described in the main text, shown on the sixth row. The last two rows show the fractional improvement on the error bar of the tensor to scalar ratio r 𝑟 r (compared to the case of no delensing, assuming r = 0 𝑟 0 r=0 ). These results are for the idealized case of no foreground or mean-field contamination, and show some sensitivity to the largest multipole L 𝐿 L that can be delensed. Ratios calculated using a lens multipole cut at L min = 100 subscript 𝐿 min 100 L_{\rm{min}}=100 instead of L min = 40 subscript 𝐿 min 40 L_{\rm{min}}=40 are shown in parentheses.

	Planck	S.O.	S4-wide	S4-deep
$N_{\mathrm{lev}}^T$ / ($\mu$K arcmin)	35	3.0	1.5	0.38
$N_{\mathrm{lev}}^P$ / ($\mu$K arcmin)	55	4.2	2.1	0.53
Beam FWHM $/$ (arcmin)	6.5	3.0	3.0	3.0
$B^{\mathrm{len}}$ / ($\mu$K arcmin) (quadratic)	4.5	3.2	2.8	2.3
$B^{\mathrm{len}}$ / ($\mu$K arcmin) (iterative)	4.5	3.0	2.4	1.4
pred. (iter. $N_L^0$)	4.5	3.0	2.4	1.4
$\sigma \left(r\right)$ impr. (quadr.)	1.0	1.7(1.5)	2.5(2.2)	4.1(3.3)
$\sigma \left(r\right)$ impr. (iter.)	1.0	1.8(1.6)	3.1(2.6)	10.6(6.3)

Equations

$$⟨X_{\mathbf{\ell }}{X}_{{\mathbf{\ell }}^{\prime }}^{†}⟩={\delta }_{\mathbf{\ell }{\mathbf{\ell }}^{\prime }}{C}_{\ell }^{\mathrm{unl}}.$$

$$D_{\mathbf{\ell }}^{TT}(\mathbf{x})\equiv \frac{1}{\sqrt{V}}{e}^{i\mathbf{\ell }\cdot \left(\mathbf{x}+\mathbf{\alpha }\left(\mathbf{x}\right)\right)},D_{\mathbf{\ell }}^{TE}(\mathbf{x})=D_{\mathbf{\ell }}^{TB}(\mathbf{x})=0.$$

$$X^{\mathrm{dat}}=BDX+n.$$

$${\text{Cov}}_{\mathbf{\alpha }}\equiv ⟨X^{\mathrm{dat}}{X}^{\mathrm{dat},†}⟩=BD{C}^{\mathrm{unl}}{D}^{†}{B}^{†}+N,$$

$$\begin{array}{rl}\ln p(X^{\mathrm{dat}}|\mathbf{\alpha })& =-\frac{1}{2}{X}^{\mathrm{dat}}\cdot {\text{Cov}}_{\mathbf{\alpha }}^{-1}{X}^{\mathrm{dat}}-\frac{1}{2}\det {\text{Cov}}_{\mathbf{\alpha }}.\end{array}$$

$$\ln p(\varphi |X^{\mathrm{dat}})=\ln p(X^{\mathrm{dat}}|\varphi )-\frac{1}{2}\sum _{\mathbf{L}}\frac{{\varphi }_{\mathbf{L}}^2}{C_{\mathbf{L}}^{\varphi \varphi }},$$

$$g_a^{\mathrm{tot}}\equiv \frac{\delta \ln p(\mathbf{\alpha }|X^{\mathrm{dat}})}{\delta {\alpha }_a(\mathbf{n})}=g_a^{\mathrm{QD}}-g_a^{\mathrm{MF}}+g_a^{\mathrm{PR}},$$

$$\frac{\delta \ln p(X^{\mathrm{dat}}|\mathbf{\alpha })}{\delta {\alpha }_a(\mathbf{n})}=g_a^{\text{QD}}(\mathbf{n})-g_a^{\mathrm{MF}}(\mathbf{n}).$$

$$g_a^{\text{QD}}(\mathbf{n})={\left[V_{\mathbf{\alpha }}{X}^{\mathrm{dat}}\right]}^i(\mathbf{n}){\left[W_{\mathbf{\alpha }}^a{X}^{\mathrm{dat}}\right]}_i(\mathbf{n}),$$

$$V_{\mathbf{\alpha }}=B^{†}{\text{Cov}}_{\mathbf{\alpha }}^{-1},W_{\mathbf{\alpha }}^a=D{\nabla }_a{C}^{\mathrm{unl}}{D}^{†}{B}^{†}{\text{Cov}}_{\mathbf{\alpha }}^{-1}.$$

$$X_{\mathbf{\alpha }}^{\mathrm{WF}}\equiv C^{\mathrm{unl}}{D}^{†}{B}^{†}{\text{Cov}}_{\mathbf{\alpha }}^{-1}{X}^{\mathrm{dat}}.$$

$$W_{\mathbf{\alpha }}^a{X}^{\mathrm{dat}}(\mathbf{x})=D{\nabla }_a{X}_{\mathbf{\alpha }}^{\mathrm{WF}}(\mathbf{x}).$$

$$V_{\mathbf{\alpha }}{X}^{\mathrm{dat}}$$

$$\begin{array}{rl}{g}_a^{\mathrm{MF}}(\mathbf{n})& =\frac{1}{2}\frac{\delta \ln \det {\text{Cov}}_{\mathbf{\alpha }}}{\delta {\alpha }_a(\mathbf{x})}=⟨g_a^{\text{QD}}(\mathbf{n})⟩.\end{array}$$

$$X_{\mathbf{\alpha }}^{\mathrm{WF}}={\left[{\left(C^{\mathrm{unl}}\right)}^{-1}+D^{†}{B}^{†}{N}^{-1}BD\right]}^{-1}{D}^{†}{B}^{†}{N}^{-1}{X}^{\mathrm{dat}},$$

$$\begin{array}{rl}{\left[D^{†}T\right]}_{\mathbf{\ell }}& =\frac{1}{\sqrt{V}}\int d^{2}x{e}^{-i\mathbf{\ell }\cdot (\mathbf{x}+\mathbf{\alpha }(\mathbf{x}))}T(\mathbf{x})\\ & =\frac{1}{\sqrt{V}}\int d^{2}x{e}^{-i\mathbf{\ell }\cdot \mathbf{x}}|M_{{\mathbf{\alpha }}^{-1}}|(\mathbf{x})T(\mathbf{x}+{\mathbf{\alpha }}^{-1}(\mathbf{x})).\end{array}$$

$${\left[M_{\mathbf{\alpha }}\right]}_{ab}(\mathbf{x})={\delta }_{ab}+\frac{\partial {\alpha }_a}{\partial x_b}(\mathbf{x}).$$

$$\mathbf{x}+\mathbf{\alpha }(\mathbf{x})+{\mathbf{\alpha }}^{-1}(\mathbf{x}+\mathbf{\alpha }(\mathbf{x}))\equiv \mathbf{x}.$$

$$\begin{array}{rl}& {{\mathbf{\alpha }}^{-1}}_{N+1}(\mathbf{n})={{\mathbf{\alpha }}^{-1}}_N(\mathbf{n})\\ & -M_{\mathbf{\alpha }}^{-1}(\mathbf{n}+{{\mathbf{\alpha }}^{-1}}_N(\mathbf{n}))\cdot ({{\mathbf{\alpha }}^{-1}}_N(\mathbf{n})+\mathbf{\alpha }(\mathbf{n}+{{\mathbf{\alpha }}^{-1}}_N(\mathbf{n}))).\end{array}$$

$${\mathbf{g}}^{\mathrm{MF}}(\mathbf{n})=⟨{\mathbf{g}}^{\mathrm{QD}}(\mathbf{n})⟩,$$

$${\left[H^{-1}\right]}_{\mathbf{L}{\mathbf{L}}^{\prime }}^{ab}\equiv -\frac{{\delta }^2\ln p(X^{\mathrm{dat}}|\mathbf{\alpha })}{\delta {\alpha }_{\mathbf{L}}^a\delta {\alpha }_{{\mathbf{L}}^{\prime }}^{b,\ast }}$$

$$H_{N+1}=(1+\rho s{y}^t)H_N(1+\rho y{s}^t)-\rho s{s}^t,$$

$$s(\mathbf{n})={\varphi }_{N+1}(\mathbf{n})-{\varphi }_N(\mathbf{n}),y(\mathbf{n})=g_{N+1}^{\mathrm{tot}}(\mathbf{n})-g_N^{\mathrm{tot}}(\mathbf{n}),$$

$${\mathbf{\alpha }}_0(\mathbf{L})=\frac{C_L^{\varphi \varphi }}{C_L^{\varphi \varphi }+N_L^{0,\mathrm{len}}}i\mathbf{L}{\hat{\varphi }}^{\mathrm{qest}}(\mathbf{L}),$$

$${\mathbf{\alpha }}_{N+1}={\mathbf{\alpha }}_N+\lambda H_N{\mathbf{g}}_N.$$

$${\hat{C}}_L^{\mathbf{g}}=\frac{1}{f_{\mathrm{patch}}}\frac{1}{n_L}\sum _{\mathbf{\ell }\text{ in }L\text{ bin }}|{\mathbf{g}}_{\mathbf{\ell }}{|}^2.$$

$${\left(N_L^0\right)}^2{\hat{C}}_L^{\mathbf{g}},\text{for }\mathbf{g}={\mathbf{g}}^{\mathrm{QD}},{\mathbf{g}}^{\mathrm{MF}},{\mathbf{g}}^{\mathrm{PR}}\text{ and }{\mathbf{g}}^{\mathrm{tot}}.$$

$$H_L^0={\left(\frac{1}{N_L^{0,\mathrm{unl}}}+\frac{1}{C_L^{\varphi \varphi }}\right)}^{-1}.$$

$${\hat{C}}_L^{\mathrm{MF}}=\frac{N_1{N}_2}{{\left(N_1+N_2\right)}^2}{\hat{C}}_L^{{\mathrm{MF}}_1-{\mathrm{MF}}_2}.$$

$$d(\mathbf{x})\equiv \frac{1}{\sqrt{V}}\sum _{\mathbf{L}}L{\hat{\varphi }}_{\mathbf{L}}{e}^{i\mathbf{L}\cdot \mathbf{x}}.$$

Full text

I Introduction

The large-scale structure of the Universe deflects Cosmic Microwave Background (CMB) photons by a few arcminutes, introducing a characteristic signature in the fluctuations in the CMB temperature and polarization Lewis and Challinor (2006). The statistical homogeneity and isotropy of the CMB gets distorted locally, and sizeable higher-order statistics are produced. Lensing estimators use these higher-order statistics to construct an integrated measure of the linear mass fluctuations in the Universe that cross-correlates to all traditional large-scale structure tracers. After the first direct detection of lensing in the CMB by the ACT team Das et al. (2011), the SPT  van Engelen et al. (2012), POLARBEAR Ade et al. (2014a), SPTpol Story et al. (2015), BICEP2-KECK Ade et al. (2016a) and Planck Ade et al. (2014b) collaborations have also reported the detection of the lensing signal and published band-powers. The most decisive detection yet was by the Planck satellite Planck Collaboration et al. (2016): its full-sky coverage comes with a statistical power that simply cannot be matched at the present time.

Current measurements all use quadratic estimator techniques, first devised in optimized form by Refs. Hu and Okamoto (2002); Okamoto and Hu (2003). The quadratic estimator uses optimally-weighted two-point statistics of the data maps to reconstruct the deflection field. At current noise levels, this estimator is nearly optimal. The science returns from the use of more sophisticated techniques are expected to be small, and no other type of estimator has been applied to data so far. However, the situation will have changed by the time of CMB stage IV (CMB-S4) Abazajian et al. (2016), if not before. At this point the polarization instrumental noise is expected to become smaller than the $\sim 5\mu K$ arcmin lensing $B$ mode. Barring welcome detections, lensing will become the most relevant cosmological source of confusion in the search for primordial $B$ modes Seljak and Hirata (2004), and more optimal delensing methods will become critical.

If the noise and primordial polarization $B$ mode is negligible, a well-known variable-counting argument Hirata and Seljak (2003a) suggests that as long as the lensing is fully described by a gradient deflection remapping of the unlensed fields, the observed lensed $E$ and $B$ fields should contain enough information to reconstruct essentially perfectly both the lensing potential and the unlensed $E$ field. Fundamental limits are well below near-future sensitivities, including corrections to the remapping approximation from emission angle, time delay and polarization rotation Hu and Cooray (2001); Lewis et al. (2017), lensing curl modes from second order post-Born lens-lens couplings Hirata and Seljak (2003a); Pratten and Lewis (2016), intrinsic non-linearities of the CMB at recombination Mollerach et al. (2004); Fidler et al. (2014), and second-order sourced vector and tensor modes Mollerach et al. (2004); Baumann et al. (2007); Sarkar et al. (2008); Adamek et al. (2016). With the last science release from the Planck team in sight, it therefore seems timely to revisit alternative, more optimal CMB lensing estimation. This paper presents and discusses a new implementation of a maximum a posteriori estimate of the lensing potential from CMB data.

Motivation for this work is not limited to primordial $B$ modes. The CMB lensing kernel peaks at $z\sim 2$ and overlaps the galaxy and weak lensing surveys targeting the dark sector of the Universe. The correlated information is expected to contribute to breaking important degeneracies and to help with systematics, so optimal CMB lensing mass maps will also be useful for use with large-scale structure observations. Iterative estimates may also be useful even at higher noise levels, in particular in the presence of sky cuts or wildly inhomogeneous noise maps where the analytic response of the quadratic estimator is inaccurate. An optimal estimate of the potential map might prove better than the current simpler practice of sweeping these deviations into Monte-Carlo (MC) corrections to the spectrum estimate. Finally, also looking a bit ahead, successful exploration of the lensed CMB likelihood may prove useful more generally, opening a path towards optimal joint estimation of the primary CMB and lensing potential.

It is clear how iterative estimates should work in practice, at least intuitively: delens the data using the quadratic deflection estimate, then again apply a quadratic estimator on the resulting maps, with possibly modified weights, and iterate until convergence Smith et al. (2012). Of course, a world of potential complications lurks in the details, and no canonically best approach is known at present. Formally, the code we present finds a maximum of the posterior probability density function (PDF) of the lensing potential. As such, it is similar in spirit to the first iterative estimator proposed for temperature reconstruction Hirata and Seljak (2003b) and polarization Hirata and Seljak (2003a). A maximum likelihood approach to lensing reconstruction is also discussed, without implementation, in the review article Ref. Hanson et al. (2010). In contrast to Refs. Hirata and Seljak (2003b, a), our implementation can be considered exact, in the sense that it maximizes the relevant functions without introducing approximations, under the assumption of Gaussian unlensed CMB, noise and deflection fields. It can also account for beams, sky cuts and other non-ideal effects.

The quadratic estimator has the convenient property of being relatively straightforward. It can be implemented using a small number of harmonic transforms Dvorkin and Smith (2009); Planck Collaboration et al. (2016), keeping the overall numerical cost under control (dominated in the Planck analysis by the cost of the inverse-variance filtering step). It seems unavoidable that alternative more-optimal approaches must be substantially more costly, and our implementation is no exception. At each iteration step, maximum a posteriori unlensed CMB maps are produced, under the assumption that the current deflection estimate is the correct one. This operation in effect solves a large $N_{\rm{pix}}\times N_{\rm{pix}}$ set of linear equations, and must be itself performed via an iterative method, each step involving a fair number of lensing operations, even in the absence of sky cuts or other non-ideal effects.

Nevertheless, lensing and lensing reconstruction have the advantage of being very local in position space. All operations scale linearly with the number of resolution elements, or follow the cost of an harmonic transform, and the good convergence properties of the iterative searches proposed here keep the total computational burden under control. We also provide GPU implementations of the most expensive steps.

We use the flat sky approximation throughout the main text. Appendix A describes the implementation on the curved sky, using the machinery of spin-weight spherical harmonics. The implementation is otherwise identical in all respects, though we have so far only thoroughly tested everything on the flat sky where the numerical implementation is faster. We expect the same convergence properties of iterative estimator on the curved sky: empirically, the only effect we observe increasing the area is to rescale the total execution time, which is reasonable given that lensing distortions are very much localized. Furthermore, iterative delensing will probably initially be most useful on deep observations of a small patch of sky where the flat sky approximation is accurate (Abazajian et al., 2016).

Sections II and III describe the algorithm and details of its implementation respectively. Section IV provides tests and applications. We summarize and conclude in Section V.

II Description

Let us first establish some notation. Let $\mathbf{x},\mathbf{y}$ be points on a patch of the flat sky of area $V$, and $\mathbf{r}=\mathbf{x}-\mathbf{y}$ be the separation vector. The primary, unlensed Stokes CMB fields $T,E$ or $B$ are written as $X(\mathbf{x})$, and $X^{\rm dat}$ denotes the observed Stokes data $T,Q$ and $U$ on the data pixels, inclusive of noise and transfer function. We use $a,b$ in $(0,1)$ to denote the two cartesian axes of the flat sky, and use the symmetric Fourier convention, which is closest to the traditional curved sky normalization. We use the notation ${\bm{\ell}}$ for multipoles of the CMB maps and $\mathbf{L}$ for the lensing maps.

We denote the primordial, unlensed CMB modes $\{T,E,B\}$ as a column matrix $X$, with primordial spectral matrix $C^{\rm unl}_{\ell}$

[TABLE]

This matrix is diagonal with respect to multipole index, but not necessarily across $T,E,B$ indices. Also let $D$ be the deflection operation that maps these unlensed CMB modes to the real space, lensed, Stokes parameters. For instance, in temperature we may write explicitly on the flat sky

[TABLE]

The polarization components are similar but involves the spin-2 flat-sky harmonics. Here, and throughout, we the use approximation that the lensed fields are entirely defined by a remapping of the unlensed fields, where ${\bm{\alpha}}(\mathbf{x})$ is the lensing deflection angle that relates the observed lensed field at $\mathbf{x}$ to the unlensed fields at $\mathbf{x}+{\bm{\alpha}}(\mathbf{x})$.

II.1 Model

The model for the CMB data $X^{\rm dat}$ on the observed pixels is given by a linear response matrix $B$ operating on the lensed sky (which includes, for example, the effect of the instrumental beam and pixel window function), plus independent noise $n$, so that

[TABLE]

The pixel-pixel covariance can be written in compact notation using a series of linear operators as follows:

[TABLE]

where $N$ is the noise covariance matrix, which we assume is diagonal in pixel space. Unlensed CMB fields and the noise on each pixel are assumed to obey Gaussian statistics, so the likelihood is also Gaussian. The log-likelihood is then

[TABLE]

We need to use a prior on the statistics of the deflection field to regularize the large number of poorly-constrained small-scale modes. The $\Lambda$CDM CMB lensing potential $\phi$ is expected to be nearly linear, so choosing Gaussian field statistics for $\phi$ is a natural choice, and will likely remain accurate in the foreseeable future on the scales where the lensing potential can accurately be reconstructed. Using a Gaussian prior on the signal does not prevent reconstruction of any non-Gaussian signal that may actually be present (as expected from non-linear structure growth and post-Born lensing Namikawa (2016); Pratten and Lewis (2016)).

We assume pure gradient lensing deflections, in which case the log-posterior becomes, up to irrelevant constants,

[TABLE]

where the likelihood is given by Eq. 5 with ${\bm{\alpha}}=\nabla\phi$. Curl and joint curl / gradient reconstruction is very analogous, but should be of limited physical relevance in the foreseeable future, mainly serving as a consistency check on the gradient reconstruction analysis. An interesting first prospect would be the detection of the post-Born curl signal, forecast to be marginally detectable in the bispectrum with CMB-S4 (Pratten and Lewis, 2016), to which this methodology could also be applied.

II.2 Gradients

To maximize the log-posterior we consider the derivative of the log-posterior with respect to the deflection. The total gradient, $g$, splits naturally into three pieces:

[TABLE]

one from the quadratic part of the likelihood ($\bm{g}^{\rm QD}$), one from the likelihood covariance determinant ($\bm{g}^{\rm MF}$, the mean-field), and one ($\bm{g}^{\rm PR}$) from the prior. The choice of the odd sign of $\bm{g}^{\rm MF}$ is more natural and becomes clear later on. The prior gradient is straightforward to evaluate assuming Gaussian statistics.

The gradients of the likelihood are first calculated in real space, with the gradients with respect to the two cartesian components of the deflection giving

[TABLE]

These are then rotated to harmonic space to give the gradient and curl components. The piece quadratic in the data

[TABLE]

is made up of two legs with data weights

[TABLE]

The gradient matrix $\nabla_{a}C^{\rm unl}$ is block diagonal in harmonic space with blocks $i{\bm{\ell}}_{a}C^{\rm unl}_{\bm{\ell}}$. These weights are identical, in the absence of deflection, to the (unnormalized) traditional Minimum Variance (MV) lensing quadratic estimators evaluated with unlensed spectra.

Both legs of the quadratic estimator (the two terms in Eq. (9)) can be written in terms of reconstructed unlensed CMB modes, as follows. Consider the most probable primordial CMB modes $X^{\rm WF}_{{\bm{\alpha}}}$ given the data, under the assumption that ${\bm{\alpha}}$ is the true deflection field, and that they are Gaussian fields with power $C_{\ell}^{\rm unl}$. The maximum a posteriori (MAP) unlensed CMB maps are formally given by the Wiener-filtered data111For a signal seen under a linear response $s^{\rm dat}=Rs^{\rm true}+n$ the maximum a posteriori reconstructed signal $s$ is given assuming Gaussian statistics by $s^{\rm WF}=(S^{-1}+R^{t}N^{-1}R)^{-1}R^{t}N^{-1}s^{\rm dat}=SR^{t}(RSR^{t}+N)^{-1}s^{\rm dat}$

[TABLE]

The leg $W^{a}_{\bm{\alpha}}X^{\rm dat}$ of the quadratic piece is then simply the deflected gradient of these maps

[TABLE]

The other leg can be written as the inverse-noise-weighted residual between the data and how the inferred primordial modes are predicted to appear:

[TABLE]

The calculation of these two terms is simple once the maps $X^{\rm WF}$ are reconstructed. Our implementation is discussed in Sec. III.1.

The second part of the likelihood gradient is the contribution from the mean field

[TABLE]

The average here is over realizations of the data, with displacement ${\bm{\alpha}}$ held fixed. The second equality follows from observing that the first variation of a log-likelihood always vanish in the mean, and that $g^{\rm MF}_{a}(\bm{n})$ itself is independent of the data (and hence is equal to its expectation). The mean field serves the same purpose here as for the traditional quadratic estimator: to subtract the known sources of anisotropy from the quadratic estimate. It depends on the current estimate of the deflection, because ${\bm{\alpha}}$ at this iteration acts as a known source of anisotropy when measuring residual lensing at the next iteration. Implementations are discussed in Sec. III.3 and Appendix B.

III Implementation

This section presents some details of our implementation. The main numerical difficulty lies in the calculation of the Wiener-filtered modes $X^{\rm WF}_{\bm{\alpha}}$ in the presence of the deflection ${\bm{\alpha}}$ (and anisotropic noise, masks, etc.). This is discussed in Sec. III.1. The Wiener-filtering operation requires the inversion of the deflection field, described in Sec. III.2, and the mean-field evaluation is discussed in Sec. III.3. We also make use of curvature information to improve convergence, as discussed in Sec. III.4. Finally, we describe our choice of starting point in Sec. III.5, and summarize the workflow of the method in Sec. III.6.

III.1 Reconstruction of the unlensed CMB

The MAP estimate of the unlensed CMB given a deflection field ${\bm{\alpha}}$ is formally given by Eq. 11, and typically requires solving a large system of linear equations. This is not straightforward even in the absence of sky cuts or other non-ideal effects, since the deflection field breaks isotropy so that the harmonic transforms do not diagonalize the system. However, for all realistic situations we have investigated, we found that the additional complication of the deflection field was minor in comparison to (typically highly anisotropic) sky cuts.

Our implementation is as follows. We first transform Eq. 11 to the following form

[TABLE]

and solve for the large inverse in brackets with conjugate gradient descent. Here one should understand the bracketed matrix to act on the space of non-zero unlensed CMB modes. There is no ambiguity regarding the unlensed spectra $(C^{\rm unl})^{-1}$, since all modes in $X^{\rm WF}$ are exactly zero when they correspond to fiducial $C_{\ell}^{\rm unl}$ that are zero: the Wiener filter builds the maximum a posteriori $X$ maps, hence these vanish whenever the prior variance $C^{\rm unl}$ does. We found Eq. (15) to be more efficient than other possible ways to perform the mask deconvolution, and it also is more efficient when the lensing operations are the only source of complications.

For the noise matrix $N$ we use an input variance map that is diagonal in pixel space. The noise can be inhomogeneous, and we can also add to the noise matrix a set of templates which are projected out. This is useful for instance to account for poorly understood low-$\ell$ noise, or to project out any templates for galactic dust. Using conjugate gradient descent requires a reasonably fast way to apply $\left(C^{\rm unl}\right)^{-1}+D^{\dagger}B^{\dagger}N^{-1}BD$ to vectors $X$. The first term is diagonal in harmonic space and poses no problem. The second term requires application of the lensing operators $D$ and $D^{\dagger}$, and of the inverse noise matrix. From Eq. 16, applying the lensing operator to a map is simply achieved by harmonic transforms followed by lensing of the resulting map. As discussed in more detail in Sec. III.2, $D^{\dagger}$ also involves the inverse harmonic transform and a delensing operation (lensing with the inverse deflection field). It therefore has the same complexity as forward lensing, provided the inverse deflection field has been precomputed. The inverse noise matrix is simple under the assumption that it is diagonal in pixel space. The inclusion of templates is only a minor complication as long as there are only a reasonable number of them. Finally, assuming isotropic beams, the beam operations are fast in harmonic space.

All in all, application of the bracketed matrix in Eq. 15 requires 4 harmonic transform and 2 lensing operations, multiplied by 1, 2 or 3 for temperature only, polarization only or joint reconstruction respectively. Lensing of maps (or the displacement inversion) is not a cheap operation, but the cost scales linearly with the number of pixels and can easily be parallelized. Our implementations, including a GPU implementation, are discussed in Sec III.2.

The use of a good preconditioner is mandatory for convergence in acceptable time, especially when dealing with masked maps. We use a multigrid preconditioner following Ref. Smith et al. (2007), where a set of working resolutions is set up so that lower-resolution inverses are used to precondition those at higher resolution. Specifically, we extend the qcinv package222https://github.com/dhanson/qcinv by Duncan Hanson to include the lensing operations. At the lowest resolution stage, we use a dense preconditioner. We offer no unique recipe of a good multigrid chain as performance appears to depend substantially on the specific configuration. The solution $X^{\rm WF}_{{\bm{\alpha}}_{N}}$ obtained at iteration ${\bm{\alpha}}_{N}$ can however be used as starting point for iteration $N+1$, which significantly speeds whole process as ${\bm{\alpha}}_{N}$ settles down to the converged estimate.

Finally, we note that in ideal situations where ${\bm{\alpha}}$ is the only source of anisotropy, the use of a simple diagonal preconditioner is much faster, with no need to resort to a multigrid solution, at least up to noise levels of a CMB-S4 configuration we have been testing.

III.2 Lensing and delensing operations

Lensing of maps is done at a resolution of 0.7 arcminutes, using a standard bicubic spline interpolation. Lensing is an expensive operation, even if easily computed in parallel, and there are a large number of maps to process until convergence is reached, and this can dominate the overall computational cost in typical runs. We found that porting the lensing on GPU, using a GPU-optimized implementation (Ruijters et al., 2008; Ruijters and Thévenaz, 2012) can provide substantial speed-up. This is one of the implementations that we provide.

In addition to the forward lensing operation, the filtering step also requires applying $D^{\dagger}$. This is equivalent to applying the inverse deflection together with multiplication by the magnification. To see this, consider the temperature part only. From the explicit form of the operator $D$ in Eq. 16 we have

[TABLE]

The second line follows from the first after the obvious change of variable $\mathbf{x}\rightarrow\mathbf{x}+{\bm{\alpha}}(\mathbf{x})$, where $|M|$ is the magnification matrix determinant that accounts for the change of volume element in these new coordinates:

[TABLE]

The inverse deflection ${\bm{\alpha}^{-1}}$ is defined by the condition that points deflected by ${\bm{\alpha}}$ are remapped to themselves

[TABLE]

Eq. 16 has the simple form of the harmonic transform of the delensed temperature map, multiplied by the magnification of the inverse deflection. The generalization to polarization is immediate.

We therefore need to obtain the inverse deflection field. While some approximation to the inverse deflection is possible given the noise levels of current data (Anderes et al., 2015; Carron et al., 2017), we found the exact inversion is always well-behaved for a $\Lambda$CDM displacement, is always in the weak-lensing regime, and is not a bottleneck for our reconstruction. The inversion is a very localized operation which is easily parallelized, for which we use a simple real space Newton-Raphson scheme on a high-resolution grid. Specifically, following Ref. Carron et al. (2017), we solve iteratively for ${\bm{\alpha}^{-1}}(\bm{n})$ using

[TABLE]

In practice, we use the same $0.7$ arcmin grid spacing that we use for the lensing operations, in which case 3 iterations starting from ${\bm{\alpha}^{-1}}=0$ are enough for essentially exact inversion of a typical $\Lambda$CDM deflection field. Typical resulting r.m.s. fractional residuals on the deflection amplitude are as low as $2\cdot 10^{-5}$. For lensing reconstruction in a realistic situation, the forward deflection is much smoother in comparison to a typical $\Lambda$CDM deflection owing to the prior effectively filtering out many small-scale modes, and coarser resolutions may also be used.

III.3 Mean field evaluation

Provided with a large number of data simulations, the mean field may be evaluated using

[TABLE]

i.e., by repeating the quadratic estimate on a number $N_{\rm MC}$ of independent simulations of the data maps and averaging to get the mean field. In practice, each of these naive estimates of $\bm{g}^{\rm MF}$ has spectrum $|\bm{g}^{\rm QD}|^{2}$ and a large resulting Monte-Carlo (MC) noise, containing the signal and noise parts of $\bm{g}^{\rm QD}$. On small scales this noise is typically much larger than $\bm{g}^{\rm MF}$, so for a reasonable number of MC simulations the mean field subtraction would effectively be adding noise with power $|\bm{g}^{\rm QD}|^{2}/N_{\rm MC}$ to the estimation of the gradient on these scales. It is therefore desirable to obtain better ways to estimate the mean field. We suggest two types of trick to accelerate convergence of the mean-field estimation.

The first simply subtracts some of the Monte Carlo noise by subtracting a mean field calculated using an isotropic approximation to the data likelihood, using the same random phases for the simulations. This introduces no bias, since by isotropy the correction vanishes in the mean, but has the virtue of cancelling part of the MC noise where the anisotropy is mild. For example, the isotropic approximation could consist of recalculating the same quadratic estimate but setting ${\bm{\alpha}}$ to zero in the weights (in the absence of other non-ideal effects), or using a simulation extended to full sky in the presence of sky cuts.

The second trick is to modify the weights of the quadratic estimator, in a way that keeps its expectation value (i.e., the mean field) constant. This is discussed in more detail in Appendix B. Combined, these tricks can lead to orders of magnitude decrease of the MC noise on the mean-field estimator, drastically reducing the number of simulations required for the same target accuracy.

In principle, different random phases must be used for the simulations at each iteration step: usage of the same phases at each step causes artificial convergence of the iteration towards what is an approximation to the true posterior. This approximation might still be fairly good, however, if enough simulations are used. At any given scale it is the mean field MC noise that sets the accuracy at which the true maximum a posteriori deflection solution can be determined. We refer to this later on as the MC noise floor.

III.4 Curvature

Finally, to perform an efficient search for the optimal point we need the curvature of the likelihood as well as the gradient. Specifically, to perform efficient Newton-type iteration across parameter space, the inverse curvature is needed. Curvature matrices such as

[TABLE]

can never be evaluated exactly in reasonable time. We proceed as follows: starting with an initial isotropic guess, $H_{0}$, we perform a rank two update to $H$ every time we move across parameter space. At each iteration, two maps are saved to disk and can be used to apply recursively the inverse curvature matrix to any vector. We use the limited memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) update (Nocedal, 1980), for which

[TABLE]

with

[TABLE]

and $\rho=1/y^{t}s$. Built in this way, the inverse curvature takes into account non-Gaussian and realization-dependent aspects of the likelihood, and we found can dramatically improve the convergence properties of the iterative search. Since the inverse curvature approximates the covariance, it can also be used to assess the width of the posterior density function, giving us approximate confidence regions for free at the end of the iterative process.

III.5 Starting point

If the posterior density were exactly Gaussian, a single Newton step starting from ${\bm{\alpha}}\equiv 0$ would bring us directly to the optimal solution. This solution matches (neglecting the difference between the realization-dependent curvature and its average) the Wiener-filtered quadratic estimator calculated with unlensed weights Okamoto and Hu (2003). However, we use lensed weights since they provide a better quadratic reconstruction Hanson et al. (2011), especially on large scales. Explicitly, we use

[TABLE]

where $N^{0,\rm len}_{L}$ is the Gaussian reconstruction noise of the quadratic estimator, calculated with the lensed CMB spectra. $N^{0,\rm len}_{L}$ is calculated using its real space flat-sky representation (see Appendix B). The quadratic estimator is the minimum variance (MV) estimator built with the set of maps considered: $T$ alone, $Q$ and $U$ polarization alone, or the three Stokes maps in combination. The filtering step is described in Sec.III.1, with the $({\bm{\alpha}}=0,C_{\ell}^{\rm len})$ Wiener filter using an input noise variance map and fiducial beam transfer function. In temperature, there are known ways to optimize the weights further Peloton et al. (2017), but Eq. (24) works well for our purposes.

III.6 Summary of the workflow

We are now in position to summarize and describe the workflow of the iterative search for the maximum a posteriori point. Initially, the displacement ${\bm{\alpha}}_{0}$ is set at the Wiener-filtered quadratic estimator as described in Sec. III.5. To get the optimal reconstruction we apply the following steps recursively until satisfactory convergence is reached:

The displacement ${\bm{\alpha}}_{N}$ is inverted to give ${\bm{\alpha}}^{-1}_{N}$ and cached. 2. 2.

With the deflection and its inverse, the delensed maps $X^{\rm WF}_{\bm{\alpha}}$ are obtained with the $({\bm{\alpha}}_{N},C_{\ell}^{\rm unl})$ Wiener filter using multigrid-preconditioned conjugate gradient inversion. This is the most expensive step of the whole process by some margin. 3. 3.

With the delensed CMB at hand, the quadratic part of gradients are calculated from Eqs. (12) and (13). 4. 4.

In parallel (unless neglected, or some other approximation scheme is used), the mean field contributions to the gradient are calculated by repeating steps 2-3 on a number of simulated maps, using the tricks discussed in Sec. B. 5. 5.

The total gradient $\bm{g}_{N}$ is then obtained, the inverse curvature updated according to the BFGS scheme of (22), and the displacement ${\bm{\alpha}}_{N+1}$ is found along the Newton descent direction:

[TABLE]

The parameter $\lambda$ helps improve convergence. For CMB-S4-like configurations, we picked $\lambda=1/2$ at all steps, whereas the full Newton step $\lambda=1$ can safely be used at higher noise levels.

IV Results

For this preliminary investigation we report three tests of our reconstruction method. First, using a simulated lensed map similar current Planck public data, we demonstrate use of the Wiener-filtering procedure to extract the maximum a posteriori estimate of the unlensed CMB. Second, we simulate a lensing reconstruction from polarization on a masked field, with noise levels of next-generation CMB experiments. In its simplest incarnation, the algorithm does not directly provide a delensed primordial $B$-mode map, but the iteratively reconstructed potential map can still be used to delens the observed $B$ modes. Third and finally, while explicit delensing of $B$-modes with these improved lensing maps will be demonstrated in upcoming work, we discuss the increase in correlation coefficient (and hence delensing efficiency) that the method can achieve with upcoming CMB data.

IV.1 Wiener filtering

Before turning to lensing reconstruction, we first demonstrate our Wiener filtering technique to accurately estimate the unlensed CMB from a simulated masked temperature CMB map. We convolved an input simulated lensed CMB sky with a beam of $7$-arcmin FWHM, and added homogeneous $35\mu$K-arcmin noise. Using the exact input lensing potential and input unlensed spectrum, we obtained the reconstructed unlensed CMB modes $T^{\rm WF}$, up to $\ell\leq 3500$. The mask was chosen arbitrarily as a piece of the public Planck 2015 analysis lensing mask, and we have further excised a band along all sides of the $600\>\rm{deg}^{2}$ patch so that the boundary of the unmasked patch is non-periodic. This leaves $f_{\rm patch}=60\%$ of the patch unmasked.

The upper panels of Fig. 1 show a comparison of simulated lensed data and the estimated unlensed CMB ($T^{\rm WF}$) using the known input deflection field. The lower panels compare the residuals $T^{\rm WF}-T^{\rm input}$, using either the deflection-dependent Wiener filter (${\bm{\alpha}}=\nabla\phi^{\rm input},C^{\rm unl}_{\ell}$, right), or the usual quadratic estimate filter that we use for the initial estimate (${\bm{\alpha}}=0,C^{\rm len}_{\ell}$, left). The former has small nearly-isotropic, near-uniform residuals set by the noise map, the latter shows the swirly patterns characteristic of lensing, with residuals of much larger amplitude.

It takes a couple of minutes on a modern laptop to reconstruct these modes, up to the point where the residual norm of the solution to the linear system of equations has decreased across the full patch by 5 orders of magnitude. At this point the values inside the masked regions have converged to sub-percent level.

IV.2 Lensing reconstruction on a masked patch

Next, we demonstrate lensing potential reconstruction from a simulated futuristic polarization-based experiment, including masking. We consider, on the same mask, the iterative reconstruction of the lensing potential from Stokes $Q,U$ simulated maps with noise level of $1.5\cdot\sqrt{2}\>\mu$K-arcmin, and a beam of $3$ FWHM. This corresponds roughly to noise levels expected for the baseline, widefield CMB-S4 configuration. We use a vanishing fiducial gravitational wave amplitude, $C^{BB,\rm unl}_{\ell}=0$, in which case the Wiener filter reconstructs the optimal unlensed $E$ mode map from $Q$ and $U$, assuming the $B$ mode is due to lensing and noise exclusively. We reconstruct multipoles of the $E$-mode map up to $\ell\leq 3500$, and the lensing potential over the same multipole range. The search converges just as well to higher multipoles, but there is little information there since the reconstruction becomes completely noise dominated. At no point do we apply low multipole cuts to the lensing potential, demonstrating that the iterative search can handle masking and the resulting large mean field adequately.

Convergence of the iterative search towards the optimal solution can be explicitly checked on all scales by monitoring the change in the posterior gradient and its components. The upper panel of Fig. 2 show the power spectra of the different pieces of the posterior gradient $g$ at the starting point of the iteration, for which the lensing reconstruction is given by the MV estimate. These gradient maps are smooth across the mask, so to estimate the spectra we simply rescale naive spectra estimates by $1/{f_{\rm patch}}$ for the purpose of this figure. Assigning Fourier modes to $L$-bin according to $L=|{\bm{\ell}}|-1/2$, resulting in $n_{L}\sim(2L+1)V/4\pi$ modes per bin, we build

[TABLE]

Shown are the quadratic piece (blue), the mean field (orange) and the prior (green) spectrum. The total gradient spectrum is shown as the red line. For easier comparison, all gradients have been normalized with the quadratic estimator normalization $N^{0}_{L}$, so the quantity plotted is actually

[TABLE]

With this normalization, the reconstruction noise in the quadratic estimate before any iteration is $N^{0}_{L}$ itself, shown as the solid black line. On scales with large signal to noise ($N^{0}\ll C_{L}^{\phi\phi}$) the inverse curvature is (to a crude approximation) $\sim N^{0}_{L}$. Thus, on these scales, $\left(N_{L}^{0}\right)^{2}C_{L}^{\bm{g}\rm tot}$ is also roughly the spectrum of the Newton increment in Eq. 25 added to the potential estimate at the corresponding iteration step.

The quadratic part of the gradient is pure reconstruction noise at small scales and is cancelled by the prior gradient. On large scales, the quadratic piece is dominated by the mean-field contamination from the mask. Since the mask mean field is basically independent from $\hat{\phi}$, it will vary little from iteration to iteration, and improving the potential estimate (by reducing the total gradient) at low multipoles demands precise evaluation of this term. On the other hand, on intermediate scales we can see that $\bm{g}^{\rm QD}$ is the dominant contribution to the total gradient at the start of the iterations. Hence, on these scales the lensing map can be improved without relying on cancellation of the mean field.

We used 511 simulations at each step to estimate the mean field. The upturn of the orange curve at $L\simeq 1000$ shows the onset of the MC noise dominated regime, where the mean-field estimate becomes pure MC noise. The dot-dashed black line shows an analytic prediction for the expected MC noise neglecting sky cuts, calculated with the tools from Appendix B. The MC noise spectrum is smaller than $N^{0}_{L}$ by four orders of magnitude, so iterations should be able to reduce the gradient amplitude by a similar amount. On large scales, the MC noise stays smaller than the prior gradient, and thus the iterative procedure will have exhausted information from the data before it hits the MC noise floor. However, we will see below that the isotropic prediction for the MC noise is inaccurate on large scales, where there is a substantial contribution from sky cuts. A small contribution also comes from the deflection field. It is possible to predict this contribution perturbatively, since the mean field just follows the spatial distribution of the deflection field at each step: see Appendix B. This contribution from the deflection is comparatively larger for temperature reconstruction, and also for temperature in combination with polarization.

We start the posterior inverse curvature $H^{0}$ with the isotropic estimate

[TABLE]

The second term is the prior curvature, and for the first term (the likelihood curvature) we used the unlensed weights. The choice of initial curvature is not critical as long as the BFGS scheme is used to update it. However, using the lensed weights can lead to the algorithm taking steps that are too large and hence give poor convergence, so we use unlensed weights instead (which gives a curvature that is slightly too large, but works well): steps that are too large should be avoided, as the search relies on the displacement being invertible at each step. More optimal curvature estimates might be built using partially lensed weights.

Convergence is acceptably quick, and after $\sim 9$ iterations the bulk of the improvement has been gained, with the gradient spectrum reduced by $2$–$3$ orders of magnitude. Only small variations on the large-scale modes are visible in the deflection maps after this point. The lower panel of Fig 2 shows the spectra after 20 iterations. The purple curve show an empirical estimate of the mean-field MC noise. This is calculated by splitting our set of simulations into two independent sets of size $N_{1}$ and $N_{2}$ with corresponding mean-field prediction $\bm{g}^{\rm{MF}_{1}}$ and $\bm{g}^{\rm{MF}_{2}}$, and building

[TABLE]

This is much larger than the analytic isotropic prediction on large scales because of the mask contribution to the MC noise. The total gradient closely follows the MC noise curve, and no further improvement can be achieved after the 20 iterations. At intermediate scales, the quadratic gradient is visibly much reduced and is now in equilibrium with the prior.

Finally, we show the reconstructed lensing map in Fig. 3. From left to right in the top row we show the input lensing map, the quadratic estimate, and the converged iterative solution. Here we plot the displacement-like but isotropic spin-0 transforms

[TABLE]

The reconstruction is visibly improved, both by large-scale modes filling in the masked regions and by the presence of finer-grained structure well inside the patch. The bottom row of Fig. 3 shows a zoom in of the central area, showing instead the lensing convergence where the improvement on small scales is more clearly visible.

IV.3 Delensing efficiency

How can our lensing reconstruction method help with measurement of primordial tensor modes? In the absence of a fiducial non-vanishing $C_{\ell}^{BB}$, for which there is at present no preferred choice, no delensed $B$-mode map is directly produced by the algorithm as the prior sets it to zero. However, it is well known that the lensing map can be used to remove some of the lensing signal in the observed $B$-mode map. Reduction of lensing signal in the $B$-mode map will result in some degree of improvement on tensor constraints, since the lensing $B$ modes act as a source of noise for any primordial signal.

Delensing of $B$-mode polarization has recently been demonstrated on Planck data by remapping the Stokes maps Carron et al. (2017), and by the SPT team Manzotti et al. (2017) using a template subtraction method. In both cases, the expected reduction of lensing-like power is approximately set by the squared cross-correlation coefficient of the measured lensing map to the true lensing map, which we call the delensing efficiency:

[TABLE]

Fig. 4 shows this cross-correlation coefficient for simulated reconstructions. We built these curves using 128 idealized simulations, with homogeneous input noise maps and no sky cuts. In this case, the deflection field is the only source of mean field at each step, which is sufficiently well described by the perturbative predictions derived and discussed in the Appendix B. In all cases, we have considered joint temperature and polarization (MV) reconstruction, with sharp multipole cuts $10<\ell\leq 3000$. For the Planck curve, we cut at $2048$, following the public analysis Ade et al. (2016b). Shown are the delensing efficiencies expected for the MV quadratic estimator (dashed colour) and the iterated, converged solution (solid). Besides Planck, we also show curves for the Simons Observatory333www.simonsobservatory.org, and two distinct CMB-S4-like configurations: one for a wide but shallow coverage, and one for a deep survey with sensitivity increased by a factor of about four. The assumed beam and noise levels are shown in Table LABEL:table:eps, ignoring all experimental complications, and are not meant to be necessarily very accurate representations of the experiment label. We have used Gaussian beams of 3 arcmin FWHM in all cases, again with the exception of Planck  where we used 6.5 arcmin. The S4-wide configuration is identical to those of the optimal reconstruction performed on the masked sky in Sec. IV.2. We show an estimate of the efficiency of this reconstruction as the green data points. The points are obtained from Eq. 31, using pseudo-$C_{\mathbf{L}}$ estimates after enlarging the mask conservatively near the mask boundaries, leaving $f_{\rm patch}\sim 35\%$, in order to avoid any edge effects. The points stand in very good agreement to expectations, demonstrating that masking does not substantially affect the reconstruction quality away from the mask boundaries.

Fig. 4 also shows the delensing efficiency reachable using the publicly available444http://pla.esac.esa.int GNILC Cosmic Infrared Background (CIB) map Aghanim et al. (2016) as an external tracer of the lensing map (purple points). We used the GNILC reconstruction at $545$ GHz, covering $60\%$ of the sky, and Planck 2015 lensing potential map to build these points. We estimated the CIB auto spectrum and the lensing-CIB cross-spectrum on the union of their released masks after apodization on a scale of $12\,{\rm arcmin}$, deconvolving the pseudo-$C_{\ell}$ estimates from the mask coupling matrix. We show the CIB efficiency estimate

[TABLE]

as the brown data points. The fiducial lensing spectrum is based on the Planck 2015 cosmology. This estimate is justified in so far as the lensing map is an unbiased tracer of the true lensing, and in the absence of spurious cross-correlation between the two maps. Comparison to previous works on CIB delensing (Fig. 2 of Ref. Larsen et al. (2016), using Planck cleaned $545$ GHz map, and Fig. 1 of Ref. Manzotti et al. (2017) from the SPT team, using Herschel 500 $\mu m$ map) shows good consistency555Note that both references show the cross-correlation coefficient $\rho_{L}$, while we show the efficiency $\rho_{L}^{2}$ over the relevant scales.

On this large sky fraction ($60\%$), contamination by galactic dust may reduce somewhat the cross-correlation to the lensing in the GNILC CIB map, and the brown points are slightly lower than expected for clean maps Ade et al. (2014c); Sherwin and Schmittfull (2015). Larger CIB efficiencies might be possible on cleaner regions of the sky, or with improved dust cleaning from future observations. For comparison we also show as the purple data points the efficiency on a smaller but cleaner area, using a mask built by thresholding the GNILC dust map at 545 GHz, keeping only $4\%$ of the sky unmasked.

Not all multipoles are equally important for the purpose of $B$-mode delensing. To a good approximation, the $B$ power $C_{\ell}^{B}$ depends linearly on the lensing deflection spectrum, hence we may write the delensed $B$-mode power as

[TABLE]

If the very-low $\ell$ reionization peak cannot be probed or is discarded, the tensor-mode recombination peak (at roughly $40\leq\ell\leq 100$) determines the scale where delensing is most important. The black line on Fig. 4 shows $LdC^{BB}_{\ell}/d\ln C_{L}^{\phi\phi}$, after averaging over this $\ell$-multipole range, and normalized such that it $L$-integrates to unity. By construction, weighting the efficiency curves on this figure against this line gives the delensing efficiency relevant for primordial $B$ modes around the recombination peak.

Effective residual delensed $B$-mode noise amplitudes are listed on the second set of rows of Table LABEL:table:eps. The delensed $B$-mode lensing power is calculated from the unlensed $E$ spectrum and a reduced lensing spectrum given by

[TABLE]

where the efficiencies are those shown in Fig. 4. The numbers in the table are the mean power over $40\leq\ell\leq 100$. The partially lensed spectra as well as the coupling matrix $dC^{BB}_{\ell}/d\ln C_{L}^{\phi\phi}$ are obtained with the Python camb package666camb.readthedocs.io. The Planck number matches well the result of the $B$-mode delensing analysis performed by Ref. Carron et al. (2017) on data. We also give predictions for the iterated solution. These predictions are obtained following Ref. Smith et al. (2012) by iteratively producing delensed power spectra and MV reconstruction noises $N_{L}^{0}$, using at each step the reduced lensing power in Eq. 34 with efficiencies

[TABLE]

to calculate the partially delensed $B$-mode power used when calculating $N_{L}^{0}$ for the next step. The predictions stand in excellent agreement with our simulated reconstructions.

Table LABEL:table:eps also shows the expected delensing improvement of constraints on the tensor-to-scalar ratio $r$, comparing results using the iterative lensing estimator to those using the quadratic estimator. These are calculated from a toy $r$ estimator variance estimate, assuming $r=0$,

[TABLE]

The delensed $B$ power is the one calculated according to Eq. 34. More realistic forecasts, for example including foreground cleaning, are well beyond the scope of this paper, but we note that these ratios stand in agreement with expectations Abazajian et al. (2016).

V Summary

We presented an iterative method for CMB lensing reconstruction, and showed that for future high-sensitivity observations it can produce substantially better results than the quadratic estimator. Even with non-trivial masking, the method remains numerically tractable and produces results in agreement with naive expectations. For low noise levels the large-scale lensing modes are all reconstructed with high signal to noise, even by the quadratic estimator, so the cosmological information is limited by cosmic variance. The main information gain from the iterative estimator comes on smaller scales where the quadratic estimator reconstruction noise starts to be substantial. However, for $B$-mode delensing, even small errors on the reconstruction of the large-scale lensing realization can lead to residual lensing $B$-mode power, so the improvement in signal to noise is important on all scales.

The algorithm works by extracting residual lensing from optimally reconstructed unlensed CMB maps. As such, it produces both estimates of the lensing potential and the delensed CMB maps. Note, however, that the method does not directly produce a delensed $B$-mode map, unless a prior spectrum is adopted for the unlensed $B$-mode spectrum. Nevertheless, the resulting deflection estimate, alone or in combination with the delensed $E$ map, may be used to delens the observed polarization map, giving improved delensing efficiency compared to using a quadratic estimator reconstruction.

The algorithm maximizes the posterior probability for the lensing potential, assuming Gaussianity of the unlensed maps and noise. A solution to this same problem was first attempted in Refs. Hirata and Seljak (2003b, a). These references, working in the absence of non-ideal effects, introduced several approximations to reduce the computational burden, avoiding in particular the anisotropic inverse variance filtering step. For similar reasons involving the difficulty of a global analysis, Ref. Anderes et al. (2011) introduced a local likelihood reconstruction method, where the lensing map is approximated as quadratic in small neighborhoods, and large wavelengths are ignored. We have demonstrated how a conjugate gradient inversion can handle the global inversion very efficiently, and, crucially, can also successfully be applied in the presence of sky cuts and other realistic non-idealities. Our solution is the first that does not rely on approximations once the fiducial ingredients of the likelihood and prior have been chosen. This means that, given enough computational resources, the resulting lensing potential map is optimal and cannot be improved upon.

In practice, one limiting factor is the mean-field calculation. Unless some approximation is used, the mean field is calculated with a finite number of simulations, and this sets a Monte-Carlo noise floor that cannot be improved upon by further iterations. However, we demonstrated that for realistic situations reconstructions can be successfully performed on masked data for current and next-generation CMB experiments. We also showed how a perturbative approximation to the mean field is adequate in the absence of non-ideal effects, allowing very fast iterative reconstructions in this case. We expect the methods and codes described and tested here to be useful for the planning and execution of future CMB lensing analyses.

The modular, fully parallelized pipeline (using MPI) is written in Python, internally calling parts written in C, and/or sending these to a GPU device using the pyCUDA interface (Klöckner et al., 2012). The flat-sky code is publicly available777https://github.com/carronj/LensIt. We also described the curved-sky algorithm; this will be tested on data and reported elsewhere.

Our successful exploration of the lensed CMB likelihood suggests several interesting possibilities for future investigation and improvement. The iterative estimate takes as an input the fiducial unlensed CMB spectra, which we have taken to include no primordial $B$ modes so that our posterior (MAP) estimate of the unlensed $B$ modes is exactly zero. This prevents us directly obtaining an optimal measurement of a delensed gravitational wave signal, which must be obtained afterwards using a more standard template subtraction or point remapping method. By allowing for non-zero unlensed $B$ modes in the prior, we could also allow direct joint estimation of the lensing together with the primordial signal. Exactly how best to do this, given the unknown amplitude of the primordial signal and complications with delensing biases, is worth careful future consideration. Another important future direction is to go beyond estimation of the lensing map to also provide optimal lensing power spectrum estimates (and estimates of the delensed CMB power spectra). Within the maximum posterior density framework, building a posterior density for the lensing power spectrum formally requires an intractable marginalization over the deflection field, though approximations can certainly be built Hirata and Seljak (2003a) and are worth further study. Finally, this also opens exciting prospects for cluster CMB lensing Baxter et al. (2015); Madhavacheril et al. (2015), by allowing non-parametric cluster mass profile measurements from the full likelihood.

Acknowledgements.

We thank Duncan Hanson for making his qcinv code publicly available, on which our filtering code is partly based. The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement No. [616170]. 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. Part of this paper is based observations obtained with Planck, an ESA science mission with instruments and contributions directly funded by ESA Member States, NASA, and Canada.

Appendix A Curved sky gradients

We give here the curved sky likelihood gradients, analogous to the flat sky version given in the main text. We first state the results and describe the implementation, and then provide a derivation. This requires only repeated use of the gradient and curl decomposition of a complex spin $s$ field, which are readily available in widespread packages. With real part $\mathcal{R}$ and imaginary part $\mathcal{I}$, we may write ${}_{\pm|s|}f$ as

[TABLE]

Then the gradient curl component are defined as ($|s|>0$)

[TABLE]

We follow here for convenience the sign conventions adopted e.g. by the relevant spin harmonic transform routines of the widespread HEALpix package Gorski et al. (2005). The definitions of the spin harmonics follow e.g. Ref. Lewis et al. (2002). The $E,B$ decomposition of the spin-$2$ polarization field is

[TABLE]

We give the results for the joint $T,Q,U$ analysis, the restriction to temperature only or polarization only is straightforward.

We aim to obtain the likelihood gradients with respect to the displacement modes, in analogy to the flat sky derivation in the main text. We only need to derive the quadratic part of the gradient: the mean field is, as before, its average, and the Gaussian prior is straightforward. We define it by projecting the gradient vector onto the spin basis {\mbox{\boldmath{e}}}_{\pm}. This basis is associated to the cartesian orthonormal frame {\mbox{\boldmath{e}}}_{1},{\mbox{\boldmath{e}}}_{2} orthogonal to $n$, given by {\mbox{\boldmath{e}}}_{\pm}\equiv{\mbox{\boldmath{e}}}_{1}\pm i{\mbox{\boldmath{e}}}_{2}. The spin $\pm 1$ quadratic gradient is then defined as

[TABLE]

The unnormalized potential ($\phi$) and curl potential ($\Omega$) quadratic estimators (as could have been obtained directly by taking gradients with respect to $\phi$ and $\Omega$) are then simply given by the harmonic expansion of ${}_{\pm 1}g$:

[TABLE]

Postponing the derivation, the end result is as follows: with $\bm{n}$ and $\bm{n}^{\prime}$ the undeflected and deflected points,

[TABLE]

In this equation, the left leg of the quadratic product is the inverse noise weighted residual

[TABLE]

and the right leg is given by deflected gradients of the Wiener-filtered maps. Explicitly,

[TABLE]

The only difference between Eq. 42 and traditional position-space curved-sky implementation of the quadratic estimator (as stated above, without the $N_{0}$ normalization) are the use of the unlensed spectra instead of the lensed spectra in producing $X^{\rm WF}$, together with the presence of the deflection operations, both in the filter and explicitly in Eq. 42.

We now justify Eq. 42. The model for the observed signal with noise $n$ is

[TABLE]

where on the curved sky the operator $D$ sends the unlensed $T,E$ and $B$ CMB modes to the deflected Stokes map with definite spin $0,\pm 2$ (i.e. ${}_{\pm 2}P$ and not $Q,U$):

[TABLE]

Similarly, $B$ projects the spin maps $T,_{\pm 2}P$ to the observed $X^{\rm dat}=T^{\rm dat},Q^{\rm dat}$ and $U^{\rm dat}$. From the definitions given in Eqs. 40 and 43, we have

[TABLE]

On the second line we used the spin $s$ as the index for the different components of the residual and gradient maps. How to make sense of and evaluate the variations of $D$? From Eq. 46, all elements are spin-weighted harmonics at the deflected position, hence we need to understand how this position changes under a variation of the deflection. The geometry is sketched on Fig. 5. On the curved sky, the notation $\bm{n}\rightarrow\bm{n}+{\bm{\alpha}}(\bm{n})$ indicates displacement of length $|{\bm{\alpha}}(\bm{n})|$ from $\bm{n}$ along the geodesic in direction ${\bm{\alpha}}$. The polarization axes are parallel transported along the geodesic, leading to some small change in $Q,U$ from the resulting misalignment with the coordinate vectors at the new point (Challinor and Chon, 2002; Lewis, 2005). Varying ${\bm{\alpha}}(\bm{n})$ by a small amount $\delta{\bm{\alpha}}(\bm{n})$ give rises to a slightly different geodesic. On the flat sky, the end separation vector between the points will be $\delta{\bm{\alpha}}(\bm{n})$, but this is not so on the curved sphere. Since $\alpha$ is typically a few arcminutes, the difference is very small, and we will neglect it. We now justify this, by deriving the exact but less practical result.

On any Riemannian manifold, the separation between the geodesics is described by the Jacobi vector $\bf J$. $\bf J$ is initially zero, has initial velocity $\delta{\bm{\alpha}}(\bm{n})$, and its acceleration is set by the Riemann curvature tensor through the Jacobi equations. The positive curvature of the sphere will reduce the separation vector compared to the flat sky. The covariant, first order change in a tensor $\mathcal{T}$ on the manifold is

[TABLE]

On the sphere, using a parallel orthonormal frame, with one vector $\bm{e_{\parallel}}$ initially aligned with ${\bm{\alpha}}(\bm{n})$, $\bf J$ is given by

[TABLE]

Hence, in this frame,

[TABLE]

This differs by $\sin\alpha/\alpha$ in the perpendicular component from the approximation that we use, where instead we use the parallel-transported $\delta{\bm{\alpha}}(\bm{n})$ (and not $\delta{\bm{\alpha}}(\bm{n}^{\prime})$),

[TABLE]

This is extremely accurate, since $1-\sin\alpha/\alpha\sim 10^{-7}$ for $\sim 2$ arcminutes deflections. With this approximation, we can make use of the spin lowering and raising form of covariant derivatives for spin weighted functions Lewis et al. (2002); Challinor and Chon (2002). The equivalent of Eq. (50) for spin-weight quantities ${}_{s}\mathcal{T}$ is

[TABLE]

Using repeatedly

[TABLE]

on all the $D$ matrix entries of Eq. (46) gives the result in Eqs. 42 and 44.

Appendix B Mean field

In this appendix we first discuss a perturbative analytic expression for the deflection-induced contribution to the mean field and then the tricks we used to reduce the number of simulations needed to calculate the mean field. We use the $\star$ operator for multiplication of infinite dimensional matrices across the survey area (continuous sky indices). In the isotropic limit this multiplication reduces to a standard convolution. Sums over discrete indices (CMB pixels, or Stoke fields) are indicated by juxtaposition. In particular, the weight matrices $W$ of the quadratic estimators act with field and pixel indices on the right, and with field and sky indices on the left. The beam operation $B$ maps the Stokes sky onto the data Stokes pixelization. Hence it acts on field and sky indices on the right, and field and pixel indices on the left.

B.1 Deflection-induced mean-field contribution

A good handle on the deflection-induced contribution to the mean field can be obtained perturbatively. We start from

[TABLE]

where unlike the previous appendix, we do not need to distinguish between upper and lower indices on the flat sky.

Recall that $\mathbf{x}$ refers to an arbitrary point on the sky, unrelated to the pixelization: we view each element of the covariance matrix as a functional of the deflection field. The linear response of the mean field is the second variation of the log-determinant functional. Hence,

[TABLE]

The first term (identical to the likelihood Fisher matrix), when evaluated at zero displacement, is simply minus the inverse $N^{0}_{L}$ lensing quadratic estimator response (evaluated with unlensed weights). It is convenient to introduce $\xi$, the real-space two-point function of the unlensed CMB fields, and $\xi_{,a}$, its derivative with respect to coordinate axis $a$. For each element of the covariance matrix, we may then write

[TABLE]

In this equation a sum over Stokes-field indices is implicit, and the pixel and further field indices are omitted on both sides in order to prevent visual cluttering. The second variation becomes

[TABLE]

where $\delta^{D}$ is the Dirac delta function. Introducing the isotropic operator

[TABLE]

all explicit dependence on the pixelization has disappeared, and a short calculation gives

[TABLE]

In harmonic space, the last term is a constant, and ensures the response to the unobservable deflection monopole $R_{ab}({\bm{\ell}}=0)$ vanishes as it should. The inverse of minus the first two terms is the usual $N^{0}$ lensing bias (here displayed with unlensed weights, and before projection onto gradient and curl components). For low noise experiments, $\xi K\xi\sim\xi$ on most scales, causing large cancellations between the second and third terms.

By design, we have thus

[TABLE]

where all terms can easily be calculated with a series of Fourier transforms.

Fig. 6 shows the expected contribution of the deflection mean field for the polarization reconstruction considered in Sec. IV.2, but with no sky cuts so that the deflection is the only mean-field source. The blue line shows the mean-field spectrum estimate obtained by averaging 500 simulations. As a test case we used an input deflection field $\nabla\phi_{0}$, with spectrum

[TABLE]

equivalent, from Eq. 24, to the spectrum of the reconstruction expected at the first iteration step. The orange curve shows the predicted spectrum from Eq. 59, which is in very good agreement given the MC noise of the mean-field estimation, visible at high multipoles.

B.2 Mean field tricks

To improve the calculation of the mean field, we introduce two tricks that together can reduce the MC noise power in the mean-field estimate by orders of magnitude. The first, and often the most powerful, is to modify the weights used in the quadratic estimator that is averaged to calculate the mean field from the simulations. The usefulness of this trick, however, depends on the adequacy of $\textrm{Cov}_{\bm{\alpha}}$ to represent the true covariance of the data. Specifically, it makes use of

[TABLE]

In practice, $\textrm{Cov}_{\bm{\alpha}}$ input to the likelihood is always going to be only an approximation to the true unknown data covariance. Accurate timestream simulations $X^{\rm dat}$ can sometimes be used to quantify the mean field (and other biases) for the standard quadratic estimator. But moderate numbers of simulations do not directly provide the means to accurately calculate $\textrm{Cov}_{\bm{\alpha}}$ or apply its inverse. The simulations themselves may of course also only capture only parts of the complexity entering the relevant systematics and the data processing. These are difficulties affecting the quadratic estimator as well, not specifically the iterative scheme proposed in this paper, that requires some amount of testing in a realistic situation.

From Eq. (14), the mean field satisfies

[TABLE]

Performing the average using the normal weights of Eq. (10) and assuming we can use Eq. 62 results in

[TABLE]

where the trace is over field indices. Only the diagonal of $H^{a}_{\bm{\alpha}}$ is relevant for the mean field, both in respect to sky and field indices, and since $H^{a}$ is a vector it vanishes for isotropic fields where there is no mean field.

It is possible to construct mean-field estimators with modified weights, or acting on different maps: as long as the expectation of the quadratic estimator remains the same they will produce unbiased estimates of the mean field. We choose weights to apply to independent unit variance Gaussian variables $s$, with $\left\langle ss^{\dagger}\right\rangle={\rm diag}(1)$ on the unmasked pixels, and use the following pair of weights

[TABLE]

The matrix inverse can be performed in the same way as the usual filtering. We now proceed to justify this choice and explain why it has lower variance (though is not exactly minimum variance).

Ideally, we would like to choose a pair of weights $W_{1},W^{a}_{2}$ to minimize the Gaussian MC noise on the mean field estimate while keeping the constraint $[W_{1}\left\langle ss^{\dagger}\right\rangle W^{a}_{2}{}^{\dagger}](\mathbf{x},\mathbf{x})=\mathrm{Tr}H_{\bm{\alpha}}^{a}(\mathbf{x},\mathbf{x})$. The real-space MC noise covariance is

[TABLE]

For simplicity, we use the constraint $W_{1}\left\langle ss^{\dagger}\right\rangle W^{a}_{2}{}^{\dagger}=W_{1}W^{a}_{2}{}^{\dagger}=H^{a}$, which is much more stringent than only matching the diagonal but makes things more tractable. We then consider minimizing the variance in the isotropic limit, in which case $W_{1}W^{a}_{2}{}^{\dagger}$ becomes a real-space convolution and hence the constraint means that in harmonic space $W^{a}_{2}({\bm{\ell}})=H^{a}_{\bm{\alpha}}({\bm{\ell}})/W_{1}({\bm{\ell}})$, where $W_{1}({\bm{\ell}})$ is a free function of multipole. A natural measure to minimize is the integrated variance from the above equation

[TABLE]

Using the constraint equation, minimizing $N_{0}(\mathbf{r}=0)$ gives

[TABLE]

so the two weight functions are proportional. From the constraint this implies that $W_{1}({\bm{\ell}})\propto\sqrt{H^{a}_{\bm{\alpha}}({\bm{\ell}})H^{a}_{\bm{\alpha}}(\mathbf{r}=0)}$.

The key message of this calculation is that to get low MC noise the scale dependence of the weight functions should be similar, though once isotropy is broken by the deflection in practice it is never possible to obtain exact square roots even in the absence of non-ideal effects. Other aspects matter as well, such as how well the unavoidable matrix inversion behaves with the chosen weights, and the contribution of the non-ideal effects to the final MC-noise floor. The form given in Eq. 65 is an empirical compromise between these considerations, that we found works well in practice. It equilibrates the weights in a very crude way, simply by having the same powers of the signal on each leg (though not the exact scale dependence). Despite being crude, it reduces the MC noise floor by more than one order of magnitude for the polarization reconstruction performed in the main text. This is shown on Fig. 6. The dash-dotted black line the MC noise of the estimator given by Eq. 65, and should be compared to the naive estimate MC noise given by the $N_{0}$ curve (black, solid). Both curves were calculated with the lensed spectra weights, which is a slightly conservative estimate choice as the iterative search converges towards the optimal solution. We found empirically that further modifying the weights to make their scale-dependence closer (for example multiplying $W_{1}$ by $\sqrt{\ell}$ and $W_{2}$ by $1/\sqrt{\ell}$ to equilibrate power in the no-beam, no-lensing, no-noise limit) did not give substantial further improvements.

The second trick is to subtract from each MC estimate the same estimate but based on an isotropic approximation to the posterior that is a close as possible to the true one, as described in Sec. III.3.