Maximum a posteriori CMB lensing reconstruction

Gravitational lensing of the cosmic microwave background (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, 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.

Figure 1

A demonstration of how our Wiener filter produces optimal (maximum a posteriori) estimates on the unlensed CMB maps from masked data. The simulated temperature map is comparable to a Planck configuration, and we use the exact input simulated deflection and exact input unlensed $C_{\ell }^{TT}$ spectrum in the filter. The upper-left panel shows the simulated masked temperature data map, with a homogeneous noise level of $35\mu \mathrm{K}$-arcmin and a beam FWHM of 7-arcmin. The (unapodized) mask is built out of a portion of the public Planck lensing mask, to which we have added a band surrounding the patch on all sides. The upper right panel shows the reconstructed unlensed map $T^{\mathrm{WF}}$. The residual to the true input CMB map ($T^{\mathrm{WF}}-T^{\text{input}}$) is shown on the lower-right panel. The lower-left panel shows the residual (on the unmasked pixels) of the result obtained when the Wiener filter instead uses no deflection but the lensed CMB spectrum in place of the unlensed spectrum (as in the standard quadratic estimator). These residuals are several times larger in magnitude (the same color scale is sometimes saturated), and display the anisotropic swirly patterns generated by the pure gradient deflection field.

Figure 2

Power spectra of the three gradients ${\mathit{g}}^{\mathrm{QD}}$ (blue), ${\mathit{g}}^{\mathrm{MF}}$ (orange) and ${\mathit{g}}^{\mathrm{PR}}$ (green) for lensing reconstruction from polarization on a masked patch, together with the total gradient spectrum (red). The algorithm works by reducing the red curve as much possible to find the most probable lensing map. The mask is shown in Fig. 1, and causes the large contribution from the mean field at low multipoles. All curves are normalized to the quadratic estimator normalization $N_L^0$. The upper panel shows the gradient spectra at the first iteration step, where the deflection is the Wiener-filtered quadratic estimator, and the lower panel shows the result after 20 iterations. At this point, the gradient has hit the mean-field MC noise floor (purple on the lower panel) on all scales and the solution cannot be improved by more iterations. The MC noise floor is mean-field estimator dependent and inversely proportional to the number of simulations used (here 511 per iteration). The dot-dashed line shows predictions for the MC noise floor built from an isotropic likelihood, which are inaccurate at low multipoles because of the sky cuts. At low multipoles, the improved reconstruction relies on accurately cancelling the mean field contribution, but on intermediate scales the decrease in the quadratic estimate is immediately visible.

Figure 3

The top-left panel shows a simulated lensing field used as input to a lensing reconstruction analysis. The input Stokes $Q$ and $U$ maps are masked with the same mask shown on Fig. 1, and we assume a polarization-sensitive experiment having polarization noise level of $1.5·\sqrt{2}\mu \mathrm{K}\text{−}\mathrm{arcmin}$ and 3 arcmin FWHM beam. The middle panel shows on the same color scale the Wiener-filtered quadratic estimate, which is the starting point of the iterative solution for the maximum a posteriori solution, Eq. (3.10). The top-right panel shows the converged solution. The top three panels show the displacement-like scalar field with transform $|\mathbf{\ell }|{\widehat{\varphi }}_{\mathbf{\ell }}$, Eq. (4.5). The lower two panels show the convergence maps $\kappa (\mathbf{\ell })=-\frac{1}{2}{\ell }^2\widehat{\varphi }(\mathbf{\ell })$ (only the central regions covering one fourth of the map) for the quadratic estimator (left) and iterative solution (right). The iterative solution can resolve structure down to smaller scales, and improvement can also be seen in the masked regions.

Figure 4

The delensing efficiencies as function of lensing $L$ multipole with either the quadratic estimator (dashed colored lines) or the iterative solution (solid colored lines), for current and futuristic noise levels. The curves were obtained from 128 simulated spectra and cross-spectra of the input lensing with idealized quadratic and iterative reconstructions. The green points show for comparison an estimate of the efficiency from the reconstruction on the masked sky described in Sec. 4b, which has identical noise level to the corresponding green curve. Only data far away from the mask edges were used to produce these points. Also shown are estimates of the efficiency using CIB maps, obtained as discussed in the main text from the public Planck GNILC maps at 545 GHz. The brown points show the efficiency obtained from 60% of the sky, while the purple points were obtained on the cleanest (according to the GNILC dust map) 4% of the sky. The black line shows the contribution per log-multipole bin of $C_L^{\varphi \varphi }$ to the total $B$-mode power on the scales relevant for a primordial $B$-mode measurement, see Eq. (4.8). The average of the colored curves weighted by the black line gives the approximate delensing efficiency relevant to each observation. The residual lensing $B$ power is listed in Table 1 together with the noise levels and expected improvement on tensor-to-scalar ratio constraints.

Figure 5

Schematic aid to the curved-sky lensed CMB likelihood gradient calculation. The variation of a lensed tensor $\mathcal{T}(\mathit{n}+\mathit{\alpha }(\mathit{n}))$ with respect to the deflection must be calculated at the deflected position $\mathit{n}+\mathit{\alpha }(\mathit{n})$. This point is defined by following the geodesic from $\mathit{n}$ in the direction $\mathit{\alpha }$ for length $|\mathit{\alpha }|$. A variation $\delta \mathit{\alpha }(\mathit{n})$ in the deflection vector at $\mathit{n}$ shifts the geodesic slightly. The exact first order change in position is given by the vector field $\mathbf{J}$ (solid arrows), proportional to $\delta \mathit{\alpha }$ infinitesimally close to $\mathit{n}$ and evolving along the original geodesic according to the Jacobi equations set by the curvature tensor. For simplicity, we instead evaluate the variation using the parallel-transported $\delta \mathit{\alpha }$ (dashed arrows) instead of $\mathbf{J}$. This neglects the focusing effect of the sphere curvature, slightly overestimating the geodesic deviation. The relative error in the gradient normal component is quadratic in the deflection angle and equal to $(1-\sin \alpha /\alpha )\sim 6\times {10}^{-8}$ for 2 arcmin deflections. This is completely negligible and of similar order as the small-angle approximation and neglected physical effects such as polarization rotation [15]: it is safe to neglect sky curvature on the scale of the deflection angles.

Figure 6

The expected contribution of the deflection-induced mean field for the first step of iterative reconstruction from polarization with noise level of $1.5·\sqrt{2}\mu \mathrm{K}\text{−}\mathrm{arcmin}$ and no sky cut. The prediction (orange line) is calculated using an input deflection with spectrum equal to that expected for the first iteration estimate of the maximum a posteriori solution [i.e., the spectrum of the Wiener-filtered deflection, Eq. (b8)]. Also shown as the (barely visible) blue curve is the measured mean field, with MC noise visible on small scales. The accuracy of the prediction is at least percent-level. Both curves were normalized by $N^0$, as the gradients on Fig 2, and can be directly compared. This shows that the $\varphi$-induced mean field plays very little role in this configuration on the first iteration, where the bulk of the reconstruction improvement is performed. The dotted-dashed black line shows (with the same normalization) the single-simulation Monte-Carlo noise of the mean-field estimator used in this work, Eq. (b12). For this configuration, it improves upon the naive $N_0$ noise (solid black) by more than an order of magnitude.