© ESO 2020

1. Introduction

Gravitational lensing distorts our view of the last-scattering surface, generating new non-Gaussian signals and B-mode polarization, as well as smoothing the shape of the observed power spectra. The large distance to recombination means that each photon is effectively independently lensed many times, boosting the signal compared to other second- and higher-order effects. The sharply-defined acoustic scale in the unlensed cosmic microwave background (CMB) perturbation power also makes small magnification and shear distortions easily detectable, allowing us to use observations of the lensed sky to reconstruct the lensing deflections and hence learn about the large-scale structure and geometry of the Universe between recombination and today (Blanchard & Schneider 1987; Hu & Okamoto 2002; Lewis & Challinor 2006). In this paper we present the final Planck¹ lensing reconstruction analysis, giving the most significant detection of lensing to date over 70% of the sky. We also give new results for polarization-based reconstructions, a combination of Planck’s lensing reconstruction and measurements of the cosmic infrared background (CIB), and delensing of the CMB temperature and polarization fields.

In the Planck 2013 analysis (Planck Collaboration XVII 2014; hereafter PL2013) we produced the first nearly full-sky lensing reconstruction based on the nominal-mission temperature data. In the 2015 analysis (Planck Collaboration XV 2016; hereafter PL2015) this was updated to include the full-mission temperature data, as well as polarization, along with a variety of analysis improvements. The final full-mission analysis presented here uses essentially the same data as PL2015: most of the map-level improvements discussed by Planck Collaboration Int. XLVI (2016) and Planck Collaboration III (2020) are focussed on large scales (especially the low-ℓ polarization), which have almost no impact on lensing reconstruction (the lensing analysis does not include multipoles ℓ <  100). Instead, we focus on improvements in the simulations, optimality of the lensing reconstruction, foreground masking, and new results such as the polarization-only reconstruction, joint analysis with the CIB, and delensing. We highlight the following main results.

– The most significant measurement of the CMB lensing power spectrum to date is performed, 9σ from polarization alone, and 40σ using the minimum-variance estimate combining temperature and polarization data on 67% of the sky, over the (conservative) multipole range 8 ≤ L ≤ 400.

– A new best estimate is made of the lensing potential over 58% of the sky by combining information from the Planck CMB lensing reconstruction and high-frequency maps as a probe of the CIB. The CIB is expected to be highly correlated with the CMB lensing potential, and although the CIB does not provide robust independent information on the lensing power spectrum, the map can provide an improved estimate of the actual realization of lensing modes down to small scales. The joint estimate gives the best picture we currently have of the lensing potential.

– Using the lensing-reconstruction maps, we demonstrate that the CMB acoustic peaks can be delensed, detecting peak sharpening at 11σ from the Planck reconstruction alone and 15σ upon further combination with the CIB (corresponding to removal of about 40% of the lensing effect). We also detect at 9σ a decrease in power of the B-mode polarization after delensing.

– Using the Planck lensing likelihood alone we place a 3.5% constraint on the parameter combination ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}$. This has comparable statistical power to current constraints from galaxy lensing, but the high-redshift CMB source plane gives a different degeneracy direction compared to the ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.5}$ combination from galaxy lensing at lower redshift. Combining a baryon density prior with measurements of baryon acoustic oscillations (BAOs) in the galaxy distribution gives a competitive measurement of σ₈, Ω_m, and H₀. We can also break the degeneracy by combining our likelihood with the first-year lensing results from the Dark Energy Survey (DES; Troxel et al. 2018), giving the tightest lensing-only constraints on these parameters.

Our baseline lensing reconstruction map is shown in Fig. 1. In Sect. 2 we explain how this was obtained, and the changes compared to our analysis in PL2015. We also describe the new optimal filtering approach used in our best polarization analysis. In Sect. 3 we present our main results, including power-spectrum estimates, cosmological parameter constraints, and a joint estimation of the lensing potential using the CIB. We end the section by using the estimates of the lensing map to delens the CMB, reducing the B-mode polarization power and sharpening the acoustic peaks. In Sect. 4 we describe in detail a number of null and consistency tests, explaining the motivation for our data cuts and the limits of our understanding of the data. We also discuss possible contaminating signals, and assess whether they are potentially important for our results. In Sect. 5 we briefly describe the various data products that are made available to the community, and we end with conclusions in Sect. 6. A series of appendices describe some technical details of the calculation of various biases that are subtracted, and derive the error model for the Monte Carlo estimates.

2. Data and methodology

This final Planck lensing analysis is based on the 2018 Planck HFI maps as described in detail in Planck Collaboration III (2020). Our baseline analysis uses the SMICA foreground-cleaned CMB map described in Planck Collaboration IV (2020), and includes both temperature and polarization information. We use the Planck Full Focal Plane (FFP10) simulations, described in detail in Planck Collaboration III (2020), to remove a number of bias terms and correctly normalize the lensing power-spectrum estimates. Our analysis methodology is based on the previous Planck analyses, as described in PL2013 and PL2015. After a summary of the methodology, Sect. 2.1 also lists the changes and improvements with respect to PL2015. Some details of the covariance matrix are discussed in Sect. 2.2, and details of the filtering in Sect. 2.3. The main set of codes applying the quadratic estimators is made public as the python package PLANCKLENS².

2.1. Lensing reconstruction

The five main steps of the lensing reconstruction are as follows.

Filtering of the CMB maps. The observed sky maps are cut by a Galactic mask and have noise, so filtering is applied to remove the mask and approximately optimally weight for the noise. The lensing quadratic estimators use as input optimal Wiener-filtered X = T, E, and B CMB multipoles, as well as inverse-variance-weighted CMB maps. The latter maps can be obtained easily from the Wiener-filtered multipoles by dividing by the fiducial CMB power spectra $C_{\mathcal{l}}^{\text{fid}}$ before projecting onto maps. We write the observed temperature T and polarization (written as the spin ±2 combination of Stokes parameters _±2P ≡ Q ± iU) pixelized data as

$$\begin{array}{c}\hfill \left(\begin{array}{c}{T}^{\mathrm{dat}}\\ {}_2{P}^{\mathrm{dat}}\\ {}_{-2}{P}^{\mathrm{dat}}\end{array}\right)=\mathcal{B}\mathcal{Y}\left(\begin{array}{c}T\\ E\\ B\end{array}\right)+\mathrm{noise},\end{array}$$(1)

where T, E, and B on the right-hand side are the multipole coefficients of the true temperature and E- and B-mode polarization. The matrix 𝒴 contains the appropriate (spin-weighted) spherical harmonic functions to map from multipoles to the sky, and the matrix ℬ accounts for the real-space operations of beam and pixel convolution. We further use the notation 𝒯 ≡ ℬ𝒴 for the complete transfer function from multipoles to the pixelized sky. The Wiener-filtered multipoles are obtained from the pixelized data as

$$\begin{array}{c}\hfill \left(\begin{array}{c}{T}^{\mathrm{WF}}\\ E^{\mathrm{WF}}\\ B^{\mathrm{WF}}\end{array}\right)\equiv C^{\mathrm{fid}}{\mathcal{T}}^{†}{\mathrm{Cov}}^{-1}\left(\begin{array}{c}{T}^{\mathrm{dat}}\\ {}_2{P}^{\mathrm{dat}}\\ {}_{-2}{P}^{\mathrm{dat}}\end{array}\right),\end{array}$$(2)

where the pixel-space covariance is Cov = 𝒯C^fid𝒯^† + N. Here, C^fid is a fiducial set of CMB spectra and N is the pixel-space noise covariance matrix, which we approximate as diagonal. As in previous releases, our baseline results use independently-filtered temperature and polarization maps (i.e., we always neglect $C_{\mathcal{l}}^{TE}$ in Cov⁻¹ in Eq. (2)) at the cost of a 3% increase in reconstruction noise on our conservative multipole range (L ≤ 400). The large matrix inversion is performed with a multigrid-preconditioned conjugate-gradient search (Smith et al. 2007). The temperature monopole and dipole are projected out, being assigned formally infinite noise. As in PL2015, we use only CMB multipoles 100 ≤ ℓ ≤ 2048 from these filtered maps. Our baseline analysis approximates the noise as isotropic in the filtering, which has the advantage of making the lensing estimator normalization roughly isotropic across the sky at the expense of some loss of optimality. In this case we also slightly rescale the filtered multipoles so that the effective full-sky transfer function matches the one seen empirically on the filtered simulations, with a minimal impact on the band powers. We also present new more optimally-filtered results, as discussed in Sect. 2.3.

Construction of the quadratic lensing estimator. We determine $\widehat{\varphi }$ from pairs of filtered maps, and our implementation now follows Carron & Lewis (2017). This differs slightly from PL2015, allowing us to produce minimum-variance (MV) estimators from filtered maps much faster, which is useful given the variety of tests performed for this release. We calculate a spin-1 real-space (unnormalized) lensing displacement estimate

$$\begin{array}{c}\hfill {}_1\widehat{d}(\widehat{\mathit{n}})=-{\displaystyle \sum _{s=0,\pm 2}{}_{-s}\overline{X}(\widehat{\mathit{n}})\left[{ð}_s{X}^{\mathrm{WF}}\right](\widehat{\mathit{n}}),}\end{array}$$(3)

where ð is the spin-raising operator, and the pre-subscript s on a field denotes the spin. The quadratic estimator involves products of the real-space inverse-variance filtered maps

$$\begin{array}{c}\hfill \overline{X}(\widehat{\mathit{n}})\equiv \left[{\mathcal{B}}^{†}{\mathrm{Cov}}^{-1}{X}^{\mathrm{dat}}\right](\widehat{\mathit{n}}),\end{array}$$(4)

and the gradients of the Wiener-filtered maps

$$\begin{array}{cc}& \left[{ð}_0{X}^{\mathrm{WF}}\right](\widehat{\mathit{n}})\equiv {\displaystyle \sum _{\ell m}}\sqrt{\ell (\ell +1)}{T}_{\ell m}^{\mathrm{WF}}{}_1{Y}_{\ell m}(\widehat{\mathit{n}}),\hfill \\ \hfill & \left[{ð}_{-2}{X}^{\mathrm{WF}}\right](\widehat{\mathit{n}})\equiv -{\displaystyle \sum _{\ell m}}\sqrt{(\ell +2)(\ell -1)}[E_{\ell m}^{\mathrm{WF}}-i{B}_{\ell m}^{\mathrm{WF}}]{}_{-1}{Y}_{\ell m}(\widehat{\mathit{n}}),\hfill \\ \hfill & \left[{ð}_2{X}^{\mathrm{WF}}\right](\widehat{\mathit{n}})\equiv -{\displaystyle \sum _{\ell m}}\sqrt{(\ell -2)(\ell +3)}[E_{\ell m}^{\mathrm{WF}}+i{B}_{\ell m}^{\mathrm{WF}}]{}_3{Y}_{\ell m}(\widehat{\mathit{n}}).\hfill \end{array}$$(5)

The deflection estimate ${}_1\widehat{d}(\widehat{\mathit{n}})$ is decomposed directly into gradient (g) and curl (c) components by using a spin-1 harmonic transform, where the gradient piece contains the information on the lensing potential and the curl component is expected to be zero:³

$$\begin{array}{c}\hfill {}_{\pm 1}\widehat{d}(\widehat{\mathit{n}})\equiv \mp {\displaystyle \sum _{\mathit{LM}}\left(\frac{{\widehat{g}}_{\mathit{LM}}\pm i{\widehat{c}}_{\mathit{LM}}}{\sqrt{L(L+1)}}\right){}_{\pm 1}{Y}_{\mathit{LM}}(\widehat{\mathit{n}}).}\end{array}$$(6)

By default we produce three estimators, namely temperature-only (s = 0), polarization-only (s = ±2), and MV (s = 0, ±2), rather than the traditional full set TT, TE, TB, EE, and EB estimators of Okamoto & Hu (2003). The temperature-polarization coupling $C_{\mathcal{l}}^{TE}$ is neglected in the C^fid factor that appears in the Wiener-filter of Eq. (2) for temperature- and polarization-only estimators, but is included in the MV reconstruction. We use lensed CMB spectra in Eq. (2) to make the estimator nearly unbiased to non-perturbative order (Hanson et al. 2011; Lewis et al. 2011). When producing the full set of individual quadratic estimators, we simply use the same equations after setting to zero the appropriate set of filtered maps entering Eq. (3).

The estimators described above only differ from the implementation described in PL2015 by the presence of the filtered B modes, B^WF, in Eq. (5). This affects the TB and EB estimators, and introduces a BB component in the polarization and MV estimators, which yields no lensing information to leading order in a cosmology with only lensing B modes. However, these modifications have a negligible impact on the reconstruction band powers and their covariance (a maximal fractional change of 0.4% for polarization only, 0.05% for MV) and we make no attempt at further optimization here.

Mean-field subtraction and normalization. This involves modification of the lensing deflection estimators in Eq. (3). Masking and other anisotropies bias the reconstruction and complicate the estimator’s response to the lensing potential, which is diagonal in the harmonic domain only under idealized conditions. The mean field is the map-level signal expected from mask, noise, and other anisotropic features of the map in the absence of lensing; we subtract this mean-field bias after estimating it using the quadratic estimator mean over our most faithful set of simulations. As in PL2015, we apply an approximate isotropic normalization at the map level, calculated analytically for the full sky following Okamoto & Hu (2003), using isotropic effective beams and noise levels in the filters. Our lensing map estimate becomes

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{\mathit{LM}}\equiv \frac{1}{{\mathcal{R}}_L^{\varphi }}({\widehat{g}}_{\mathit{LM}}-{\langle {\widehat{g}}_{\mathit{LM}}\rangle }_{\mathrm{MC}}),\end{array}$$(7)

and similarly for the lensing curl. With the notational conventions adopted above, the responses are identical to those defined in PL2015. The isotropic normalization is fairly accurate on average: cross-spectra between reconstructions from masked simulations and the true input lensing realizations match expectations to sub-percent levels on all but the largest scales. For the released lensing maps, the subtracted mean field is calculated across the entire available set of 300 simulations (see below), and is also provided. For power-spectrum estimation, we use cross-spectrum estimators of maps with independent Monte Carlo noise on the mean-field subtraction, obtained using 30 independent simulations for the mean field subtracted from each map (60 simulations in total). This number is motivated by a nearly optimal trade-off between uncertainties in the Monte Carlo estimates of the mean field and biases that affect the reconstruction band-power covariance matrix (see Sect. 2.2).

Calculation of the power spectrum of the lensing map and subtraction of additional biases. More specifically we need to perform subtraction of the so-called N⁽⁰⁾ and N⁽¹⁾ lensing biases, as well as point-source contamination. For a pair of lensing map estimates ${\widehat{\varphi }}_1$ and ${\widehat{\varphi }}_2$, we use the same simple cross-spectrum estimator as in PL2015,

$$\begin{array}{c}\hfill {\widehat{C}}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}\equiv \frac{1}{(2L+1)f_{\mathrm{sky}}}{\displaystyle \sum _{M=-L}^L{\widehat{\varphi }}_{1,LM}{\widehat{\varphi }}_{2,LM}^{\ast },}\end{array}$$(8)

from which biases are subtracted:

$$\begin{array}{c}\hfill {\widehat{C}}_L^{\varphi \varphi }\equiv {\widehat{C}}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}-{\mathrm{\Delta }{C}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}|}_{\mathrm{RDN}0}-{\mathrm{\Delta }{C}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}|}_{\mathrm{N}1}-{\mathrm{\Delta }{C}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}|}_{\mathrm{PS}}.\end{array}$$(9)

Lensing power-spectrum estimation is designed to probe the connected 4-point function of the data that is induced by lensing. The combination of the mean-field subtraction and the first bias term $\mathrm{\Delta }{C}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}{|}_{\mathrm{RDN}0}$ (N⁽⁰⁾) in Eq. (9) subtracts the disconnected signal expected from Gaussian fluctuations even in the absence of lensing, and is calculated using the same realization-dependent N⁽⁰⁾ (RD-N⁽⁰⁾) estimator described in PL2015 (and summarized in Appendix A for the baseline cases). The $\mathrm{\Delta }{C}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}{|}_{\mathrm{N}1}$ term subtracts an $\mathcal{O}(C_L^{\varphi \varphi })$ signal term (N⁽¹⁾) coming from non-primary couplings of the connected 4-point function (Kesden et al. 2003), and in our baseline analysis this is calculated using a full-sky analytic approximation in the fiducial model as described in Appendix A. The signal-dependence of this term is handled consistently in the likelihood as described in Sect. 3.2. We tested an alternative simulation-based N⁽¹⁾ calculation and a deconvolution technique as described in Sect. 4.8. The point-source (PS) bias term $\mathrm{\Delta }{C}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}{|}_{\mathrm{PS}}$ subtracts the contribution from the connected 4-point function of unclustered point sources. The amplitude of this correction is estimated from the data, as described in PL2015.

The MV estimator empirically has very slightly larger reconstruction noise N⁽⁰⁾ than the TT estimator for lensing multipoles L  ≃  1000 and beyond, as was also the case in PL2015: the simple combination of the various quadratic estimators in Eq. (3) is slightly suboptimal if the fiducial covariance matrix does not exactly match that of the data. This is at most a 2% effect at the highest multipoles of the reconstruction, and is sourced by our choice of independently filtering the temperature and polarization data (i.e., the neglect of $C_{\mathcal{l}}^{TE}$ in Cov in Eq. (2)). Therefore, we have not attempted further optimization in Eq. (3).

Binning, and application of a multiplicative correction. This final correction is obtained through Monte Carlo simulations to account for various approximations made in the previous steps, including correcting the approximate isotropic normalization assumed. After converting the lensing potential spectra to convergence (κ) spectra ($C_L^{\kappa \kappa }=L^2{(L+1)}^2{C}_L^{\varphi \varphi }/4$), we define our band-power estimates

$$\begin{array}{c}\hfill {\widehat{C}}_{L_{\mathrm{b}}}^{\kappa \kappa }\equiv {\displaystyle \left(\sum _L{\mathcal{B}}_{\mathrm{b}}^L{\widehat{C}}_L^{\kappa \kappa }\right)}\left(\frac{{\sum }_L{\mathcal{B}}_{\mathrm{b}}^L{C}_L^{\kappa \kappa ,\mathrm{fid}}}{{\sum }_L{\mathcal{B}}_{\mathrm{b}}^L{\langle {\widehat{C}}_L^{\kappa \kappa }\rangle }_{\mathrm{MC}}}\right).\end{array}$$(10)

The binning functions ${\mathcal{B}}_{\mathrm{b}}^L$ use an approximate inverse-variance weighting $V_L^{-1}\propto (2L+1)f_{\mathrm{sky}}{\mathcal{R}}_{\mathcal{L}}^{\mathcal{2}}/[2L^4{(L+1)}^4]$ to produce roughly optimal signal amplitudes:

$$\begin{array}{c}\hfill {\mathcal{B}}_{\mathrm{b}}^L=C_{L_{\mathrm{b}}}^{\kappa \kappa ,\mathrm{f}id}\frac{C_L^{\kappa \kappa ,\mathrm{fid}}{V}_L^{-1}}{{\sum }_{L^{\prime }}{\left(C_{L^{\prime }}^{\kappa \kappa ,\mathrm{fid}}\right)}^2{V}_{L^{\prime }}^{-1}},L_{\min }^{\mathrm{b}}\le L\le L_{\max }^{\mathrm{b}}.\end{array}$$(11)

Equation (11) rescales the amplitude measurements by the fiducial convergence spectrum interpolated to the bin multipole L = L_b, where the bin multipoles are the weighted means,

$$\begin{array}{c}\hfill L_{\mathrm{b}}\equiv \frac{{\sum }_{L}L{\mathcal{B}}_{\mathrm{b}}^L}{{\sum }_{L^{\prime }}{\mathcal{B}}_{\mathrm{b}}^{L^{\prime }}}·\end{array}$$(12)

This choice ensures that the binned fiducial spectrum goes exactly through the fiducial model at L = L_b, so plotting band-power bin values at L = L_b against the unbinned fiducial model gives a fair visual comparison of whether the observed band power is higher or lower than the fiducial one (assuming that the true spectrum shape is close to the fiducial shape). For a flat convergence spectrum, Eq. (12) gives the centre of mass L of the bin.

In a change to the earlier analyses, we no longer subtract a Monte Carlo correction from the estimated lensing power spectrum, instead making a multiplicative correction. The ratio on the right-hand side of Eq. (10) is our multiplicative Monte Carlo correction, which corrects for the various isotropic and simplifying approximations we make in constructing the unbinned power-spectrum estimator. The simulation-averaged band powers ${\langle {\widehat{C}}_L^{\kappa \kappa }\rangle }_{\mathrm{MC}}$ in Eq. (10) are built from simulations as from the data according to Eq. (9), but with a cheaper Monte Carlo N⁽⁰⁾ (MC-N⁽⁰⁾) estimation described in Appendix A, and no point-source correction (since the simulations are free of point sources). The (reciprocal of the) Monte Carlo correction for our baseline MV band powers is illustrated later in Sect. 2.3 (see Fig. 3 there).

Other differences to the analysis of PL2015 include the following points.

– An improved mask, with reduced point-source contamination for the same sky fraction. The amplitude of the point-source correction decreased by a factor of 1.9, and the detection of this point-source contamination is now marginal at 1.7σ. The 2013 and 2015 lensing analyses used essentially the same mask, constructed as described in PL2013. This is now updated using a combination of unapodized masks: a SMICA-based confidence mask;⁴ the 2015 70% Galactic mask; and the point-source masks at 143 GHz and 217 GHz. We also consider a mask targeted at the resolved Sunyaev–Zeldovich (SZ) clusters with S/N >  5 listed in the 2015 SZ catalogue⁵. This has little impact on the results, but is included in the baseline analysis, leaving a total unmasked sky fraction f_sky = 0.671. A reconstruction map without the SZ mask is also made available for use in SZ studies.

– We continue to use the foreground-cleaned SMICA maps for our baseline analysis; however, the details of the SMICA processing have changed, as described in Planck Collaboration IV (2020). Specifically, the SMICA weights at high ℓ relevant for lensing are now optimized over a region of the sky away from the Galaxy (but larger than the area included in the lensing mask), significantly changing the relative weighting of the frequency channels on small scales. This changes the noise and residual foreground realization in the SMICA maps compared to the 2015 analysis, and hence the lensing reconstruction data points scatter with respect to 2015 by more than would be expected from individual frequencies. The 2018 SMICA maps also correct an LFI map calibration issue in the 2015 maps that affected the amplitude around the first peak; however, this had little impact on the 2015 lensing analysis, since the great majority of the signal comes from smaller scales.

– Monte Carlo evaluation of the mean field and bias terms are now based on Planck FFP10 simulations, described in detail in Planck Collaboration III (2020). In addition to many processing changes, the simulations fix an error in the FFP8 simulation pipeline used for PL2015, which led to aberration (due to the motion of the Solar System relative to the CMB rest frame) not being simulated; the new simulations include the expected level of aberration and the associated modulation (Planck Collaboration XXVII 2014), although this has minimal impact on the lensing analysis. There are only 300 noise FFP10 simulations, so we now include various sources of Monte Carlo error from the finite number of simulations as additional contributions to the covariance matrix. The noise simulations were generated using a single fiducial foreground and CMB realization, which is subtracted before adding to the signal simulations. Small nonlinearities in the processing cause a weakly correlated residual between simulations. This residual can be detected both in the temperature and polarization simulation mean fields at very high lensing multipoles (see Appendix B), where it can be seen that the impact on our band powers is completely negligible compared to the error bars.

There is a roughly 3% mismatch in power at multipoles ℓ ≃ 2000 between the data and the FFP10 temperature simulations (which have no variance from residual foregrounds). We account for this by adding isotropic Gaussian noise to the simulations, with a spectrum given by the power difference. The size of this component is roughly 5 μK-arcmin with a weak scale dependence. In polarization, as discussed in Planck Collaboration IV (2020), the simulation power can be slightly larger than the data power. In this case, for consistency, we add a small additional noise component to the data maps.

– The lensing maps that we release are provided to higher L_max  =  4096 than in 2015. Multipoles at L ≫ 60 become increasingly noise dominated, but some residual signal is present at L >  2048, so we increase the range, following requests related to cross-correlation and cluster analyses (Geach & Peacock 2017; Singh et al. 2017). The reconstruction with the full multipole range is made publicly available, but in this paper we only show results for the power spectrum at multipoles up to L_max = 2048; we have not studied the reliability of reconstructions at higher multipoles, so we recommend they be used with caution.

– The treatment of the Monte Carlo (MC) correction differs, being now multiplicative instead of additive. After subtraction of the lensing biases and formation of the band powers, we calculate the MC correction by taking the ratio to the appropriately binned fiducial $C_L^{\text{fid},\varphi \varphi }$. The choice of a multiplicative correction is more appropriate for mode-mixing effects, where corrections are expected to scale with the signal, and for calibration of the quadratic estimator responses when using inhomogeneous filtering. Our baseline reconstruction MC correction is most important on large scales (where it is around 10%), but only has a small impact on the band-power errors.

– The lensing likelihood is constructed as before, following Appendix C of PL2015. We now include L ≤ 4096 in the calculation of the fiducial N⁽¹⁾ bias that we subtract, and include L ≤ 2500 in the linear correction to account for the model-dependence of N⁽¹⁾ relative to the fiducial lensing power. The 2015 MV likelihood contained an almost inconsequential error in the calculation of the response of N⁽¹⁾ to the polarization power that has now been corrected.

The construction of the covariance matrix also differs slightly, with additional small terms to take into account uncertainties in several factors that are calibrated in simulations (see Sect. 2.2). For “lensing-only” parameter results we now adopt slightly tighter priors, and marginalize out the dependence on the CMB spectra given the observed Planck data, as described in Sect. 3.2.1.

– Our fiducial model, the same as used to generate the FFP10 simulations, is now a spatially-flat ΛCDM cosmology with: baryon density ω_b ≡ Ω_bh² = 0.02216; cold dark matter density ω_c ≡ Ω_ch² = 0.1203; two massless neutrinos and one massive with mass 0.06 eV; Hubble constant H₀ = 100h km s⁻¹ Mpc⁻¹ with h = 0.670; spectral index of the power spectrum of the primordial curvature perturbation n_s = 0.964; amplitude of the primordial power spectrum (at k = 0.05 Mpc⁻¹) A_s = 2.119 × 10⁻⁹; and Thomson optical depth through reionization τ = 0.060.

2.2. Covariance matrix

Our band-power covariance matrix is obtained from the FFP10 simulation suite. Out of the 300 simulations, 60 are used for the mean-field subtraction (30 for each of the quadratic reconstructions that are correlated to form the power spectrum), and 240 for estimation of the lensing biases, Monte Carlo correction, and band-power covariance matrix. It is impractical to perform the same, expensive, realization-dependent N⁽⁰⁾ (RD-N⁽⁰⁾) subtraction on all these simulations for evaluation of the covariance matrix, and therefore, as in previous releases, we use a cheaper semi-analytic calculation, as detailed in PL2015, which only requires empirical spectra of the CMB data. This semi-analytic calculation is only accurate to 1–2%, which is not enough for debiasing where sub-percent accuracy is required to recover the lensing signal at high lensing multipoles; however, it is sufficient for the covariance matrix calculation.

New to this release are two corrections to the covariance matrix that slightly increase the error bars. First, we take into account Monte Carlo uncertainties in the mean-field, RDN0, and MC corrections. As detailed in Appendix C, for Planck noise levels the additional variance ${\sigma }_{\text{MC}}^2$ caused by the finite number of simulations can be written to a good approximation in terms of the band-power statistical⁶ errors ${\sigma }_{\text{BP}}^2$ as

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{MC}}^2\simeq (\frac{2}{N_{\mathrm{MF}}}+\frac{9}{N_{\mathrm{Bias}}}){\sigma }_{\text{BP}}^2.\end{array}$$(13)

Here, N_MF is the number of simulations entering the mean-field subtraction, and N_Bias the number used for the noise biases and MC correction. Our choice (N_MF = 60 and N_Bias = 240) is close to optimal, given the 300 simulations at our disposal and our choice of N⁽⁰⁾ estimator. To account for the finite number of simulations, we have simply rescaled the entire covariance matrix by this factor, a 7% increase in covariance, irrespective of binning. Second, we also rescale our inverse covariance matrix by the factor (Hartlap et al. 2007)

$$\begin{array}{c}\hfill {\alpha }_{\mathrm{cov}}=\frac{N_{\mathrm{var}}-N_{\mathrm{bins}}-2}{N_{\mathrm{var}}-1},\end{array}$$(14)

where N_var (which equals N_bias in our analysis) is the number of simulations used to estimate the covariance matrix, to correct for the bias that would otherwise be present in the inverse covariance matrix that is used in the likelihood. We construct two likelihoods: one based on the conservative multipole range 8 ≤ L ≤ 400, for which the number of band-power bins N_bins = 9 and 1/α_cov = 1.035 (i.e., effectively a 3.5% increase in band-power covariance); and one on the aggressive multipole range 8 ≤ L ≤ 2048, for which N_bins = 16 and 1/α_cov = 1.06. After our semi-analytical realization-dependent debiasing, the covariance matrix shows no obvious structure on either multipole range. We find all individual cross-correlation coefficients to be smaller than 10% and consistent with zero, a constraint limited by the number of simulations available. We choose to include the off-diagonal elements in the likelihood.

Following PL2015 we subtract a point-source template correction from our band powers, with an internally measured amplitude (the point-source shot-noise trispectrum ${\widehat{S}}_4$). We neglect the contribution to the error from point-source subtraction uncertainty, since for this release the estimated error on ${\widehat{S}}_4$ would formally inflate band-power errors by at most 0.4% at L ≃ 300, where the correction is strongest, and much less elsewhere.

2.3. Inhomogeneous filtering

Approximating the noise as isotropic for filtering is suboptimal because the Planck scanning results in significant noise anisotropy with a dynamic range of 10 for polarization and 5 for temperature, after allowing for residual foregrounds (see Fig. 2). In this section we describe a new polarization-only reconstruction using inhomogeneous filtering, demonstrating a large improvement over polarization results that use homogeneous filtering. However, inhomogeneous filtering is not used for our main cosmology results including temperature, where it makes little difference but would complicate the interpretation.

The SMICA CMB map is constructed from Planck frequency maps using isotropic weights $w_{\mathcal{l}}^{X,f}$ per frequency channel f. To construct the noise variance map used in the inhomogeneous filtering, we first combine the variance maps from the individual frequency maps with the SMICA weights to obtain the total noise variance in each pixel of the SMICA map. More specifically, in polarization, neglecting Q, U noise correlations, defining the pixel noise variance ${\sigma }_P^2\equiv \left({\sigma }_Q^2+{\sigma }_U^2\right)/2$, and neglecting differences between ${\sigma }_Q^2$ and ${\sigma }_U^2$, we have

$$\begin{array}{c}\hfill {\sigma }_P^2(\widehat{\mathit{n}})={\displaystyle \sum _{\mathrm{freq}.f}{\int }_{S_2}\mathrm{d}\widehat{\mathit{n}}\prime {\sigma }_{P,f}^2(\widehat{\mathit{n}}\prime )\sum _{s=\pm 2}\frac{1}{4}{[{\xi }_{2,s}^{E,f}(\widehat{\mathit{n}}·\widehat{\mathit{n}}\prime )\pm {\xi }_{2,s}^{B,f}(\widehat{\mathit{n}}·\widehat{\mathit{n}}\prime )]}^2,}\end{array}$$(15)

where

$$\begin{array}{c}\hfill {\xi }_{s,s\prime }^{X,f}(\mu )\equiv {\displaystyle \sum _{\ell }\left(\frac{2\ell +1}{4\pi }\right)w_{\ell }^{X,f}{d}_{s,s\prime }^{\ell }(\mu ),}\end{array}$$(16)

and $d_{s,s^{\prime }}^{\mathcal{l}}$ are reduced Wigner d-matrices. This equation follows from transforming the noise maps at each frequency f, which have pixel variance ${\sigma }_{P,f}^2(\widehat{\mathit{n}})$ and are further assumed uncorrelated across frequencies, into their E- and B-modes, applying the SMICA weights $w_{\mathcal{l}}^{X,E}$ and $w_{\mathcal{l}}^{B,f}$, and transforming back to Q and U in pixel space. Averaging the pixel variances of the resulting Q and U noise maps yields ${\sigma }_P^2$. Physically, the variances of the frequency maps are combined non-locally with kernels that derive from the convolutions implied by the SMICA weights. In temperature,

$$\begin{array}{c}\hfill {\sigma }_T^2(\widehat{\mathit{n}})={\displaystyle \sum _{\mathrm{freq}.f}{\int }_{S_2}\mathrm{d}\widehat{\mathit{n}}\prime {\sigma }_{T,f}^2(\widehat{\mathit{n}}\prime ){\left[{\xi }_{0,0}^{T,f}(\widehat{\mathit{n}}·\widehat{\mathit{n}}\prime )\right]}^2·}\end{array}$$(17)

The pixel noise is expected to be correlated to some degree, both due to the mapmaking process and the SMICA weighting. We have seen no evidence that this is relevant to the reconstruction and for simplicity neglect it for the filter. SMICA uses a hybrid method with two different set of weights in temperature. Only the second set, used for high multipoles and well away from the Galactic plane, is relevant for the sky area and multipoles used by the lensing analysis. The small mismatch in the power in the FFP10 simulations and the data is corrected as described in Sect. 2.1. The noise variance maps are shown in Fig. 2.

Using the variance maps in our filtering slightly slows down convergence to the solution, especially in temperature. Empirically, we found that using smoothed variance maps produces reconstructions with the same signal-to-noise ratio (S/N), as long as the mostly quadrupolar structure of the map is resolved. We use the variance maps of Fig. 2 smoothed with a Gaussian width of 10° for our quoted results; the impact on the execution time compared to homogeneous filtering is negligible, with the outputs virtually identical to filtering with the high-resolution variance maps. Using the noise anisotropy in the filter downweights modes in more noisy regions of the map, and hence improves the optimality of the estimator, especially in polarization, where the Planck data are noise dominated on most scales. The expected cross-correlation coefficient of the lensing potential estimate to the true signal improves by 20% over the conservative multipole range 8 ≤ L ≤ 400, and the S/N of the band powers increases by 30%. The filter has less impact on the temperature, which is signal dominated, and the reconstruction is of essentially the same quality with or without inhomogeneous filtering.

One disadvantage of the anisotropic filter is that it complicates the estimator’s response: the correct normalization of the estimator becomes position dependent. We have not attempted to perform a full map-level normalization, instead simply making the additional correction as part of the Monte Carlo correction we apply to our band powers. At Planck’s noise levels we expect this procedure to be very close to optimal for polarization, and for most (but not all) scales for temperature. Approximating the sky as a collection of independent patches with roughly constant noise within a patch, a full-sky optimal lensing spectrum estimation is obtained by inverse-variance weighting correctly-normalized spectra in each patch. Since the estimators are optimal in each patch, the estimator normalization in each patch is identical to the reconstruction noise level N⁽⁰⁾. Therefore, whenever the reconstruction noise dominates the lens cosmic variance in the band-power errors, inverse-variance weighting of the patch spectra is equivalent to uniform weighting of the unnormalized estimators’ spectra.

We can predict the lensing spectrum’s Monte Carlo correction fairly accurately using the simple independent-patch approximation. Let ℛ_L denote the response of the quadratic estimator to the lensing signal, such that under idealized conditions the properly normalized lensing map estimate is (as in Eq. (7))

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{\mathit{LM}}\equiv \frac{{\widehat{g}}_{\mathit{LM}}}{{\mathcal{R}}_L}·\end{array}$$(18)

Applying a single fiducial response ${\mathcal{R}}_L^{\mathrm{fid}}$ on the full-sky estimate, the local estimate in a patch centred on $\widehat{\mathit{n}}$ is biased by a factor ${\mathcal{R}}_L(\widehat{\mathit{n}})/{\mathcal{R}}_L^{\mathrm{fid}}$, where ${\mathcal{R}}_L(\widehat{\mathit{n}})$ is the true response according to the local temperature and polarization filtering noise levels. We may then write the multipoles of the full-sky lensing map as a sum of multipoles extracted over the patches,

$$\begin{array}{c}\hfill {\widehat{\varphi }}_{\mathit{LM}}\simeq {\displaystyle \sum _{\text{patches p}}\frac{{\mathcal{R}}_L({\widehat{\mathit{n}}}^{\mathrm{p}})}{{\mathcal{R}}_L^{\mathrm{fid}}}{\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{p}},}\end{array}$$(19)

where each unbiased component ${\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{p}}$ is obtained from the patch p. Using a large number of patches, neglecting correlations between patches, and turning the sum into an integral gives the following useful approximate result for the correlation of the estimator with the input

$$\begin{array}{c}\hfill \frac{\langle {\widehat{C}}_L^{\widehat{\varphi }{\varphi }_{\mathrm{in}}}\rangle }{C_L^{\varphi \varphi ,\mathrm{f}id}}\simeq {\displaystyle \int \frac{\mathrm{d}\widehat{\mathit{n}}}{4\pi }\left(\frac{{\mathcal{R}}_L(\widehat{\mathit{n}})}{{\mathcal{R}}_L^{\mathrm{fid}}}\right),}\end{array}$$(20)

and equivalently for the estimator power spectrum,

$$\begin{array}{c}\hfill \frac{\langle {\widehat{C}}_L^{\varphi \varphi }\rangle }{C_L^{\varphi \varphi ,\mathrm{f}id}}\simeq {\displaystyle \int \frac{\mathrm{d}\widehat{\mathit{n}}}{4\pi }{\left(\frac{{\mathcal{R}}_L(\widehat{\mathit{n}})}{{\mathcal{R}}_L^{\mathrm{fid}}}\right)}^2.}\end{array}$$(21)

The spectrum-level correction of Eq. (21) is only close to the squared map-level correction of Eq. (20) (which can be made close to unity using a refined choice of fiducial response) if the true responses do not vary strongly across the sky. This is the case for the signal-dominated temperature map; however, the responses vary by almost an order of magnitude in the polarization map. The blue points in Fig. 3 show the empirical Monte Carlo correction we apply to our inhomogeneously-filtered polarization band powers, together with the prediction from Eq. (21). The agreement is visually very good, with a residual at low-L that originates from masking, also found on our baseline, homogeneously-filtered MV band powers (orange points). This large-scale MC correction has a significant dependence on the sky cut, but little dependence on other analysis choices; the lensing reconstruction is close to local in real space, but this breaks down near the mask boundaries. Finally, while our power-spectrum estimator in Eq. (8) does not attempt to remove any mode-mixing effect of masking on the lensing estimate, we note that the large-scale MC correction is not simply just a ϕ mode-mixing effect: using a pseudo-C_ℓ inversion to construct the band powers from the masked ϕ map makes almost no difference to the MC correction.

As part of the Planck 2018 release, we make available products containing both versions of lensing maps. Optimally-weighted maps can be used if S/N is critical, while the isotropically-weighted maps can be used to simplify cross-correlation analyses if required. Our baseline results and likelihoods are still derived from the simpler isotropic filtering, but there is a substantial improvement in the polarization-only reconstruction when using the more optimal weighting. In all cases the reconstruction noise properties of the maps are best assessed using the corresponding set of released simulations.

3. Results

3.1. Lensing-reconstruction map and power spectrum

In Fig. 1 we show our baseline Wiener-filtered minimum-variance lensing deflection estimate from the Planck temperature and polarization SMICA CMB maps. This is shown as a map of

$$\begin{array}{c}\hfill {\widehat{\alpha }}_{\mathit{LM}}^{\mathrm{WF}}=\sqrt{L(L+1)}\frac{C_L^{\varphi \varphi ,\mathrm{fid}}}{C_L^{\varphi \varphi ,\mathrm{fid}}+N_L^{\varphi \varphi }}{\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{MV}},\end{array}$$(22)

where $C_L^{\varphi \varphi ,\mathrm{fid}}$ is the lensing potential power spectrum in our fiducial model and $N_L^{\varphi \varphi }$ is the noise power spectrum of the reconstruction. The quantity ${\alpha }_{\mathit{LM}}^{\mathrm{WF}}=\sqrt{L(L+1)}{\varphi }_{\mathit{LM}}^{\mathrm{WF}}$ is equivalent to the Wiener-filtered gradient mode (or E mode) of the lensing deflection angle. For power-spectrum estimates we plot ${[L(L+1)]}^2{C}_L^{\varphi \varphi }/2\pi =L(L+1)C_L^{\alpha \alpha }/2\pi$, so that a map of α has the same relation to the plotted power spectrum as the CMB temperature map does to $\mathcal{l}(\mathcal{l}+1)C_{\mathcal{l}}^{TT}/2\pi$. As in 2015 we exclude L <  8 due to the high sensitivity to the mean-field subtraction there. The characteristic scale of the lensing modes visible in the reconstruction is L ≃ 60, corresponding to the peak of the deflection power spectrum, where the S/N is of order 1. The left panel in Fig. 4 shows our corresponding baseline MV reconstruction power spectra over the conservative and aggressive multipole ranges.

In addition to the MV reconstruction, we also provide temperature-only results, as well as two variants of the polarization-only reconstruction: the first polarization reconstruction uses the same homogeneous filtering as the temperature and MV results; the second uses the more optimal filter, based on the variance maps shown in Fig. 2. The inhomogeneous filtering gives a large improvement in the precision of the polarization-only reconstruction, as shown in the right panel in Fig. 4. No significant improvement with inhomogeneous filtering is expected (or found) for the temperature and MV reconstructions, so we do not give results for them.

Table 1 lists our band-power measurements. For each spectrum, we provide the amplitude relative to the fiducial band powers in the first column, and the fiducial band powers in the second. The binning function was given in Eq. (11) and uses an approximate analytic inverse-variance weighting of the unbinned spectra. The fiducial band powers can therefore show slight variations because the noise varies between the different reconstructions.

Section 3.2.1 introduces our new lensing-only likelihood, marginalizing over the CMB spectra. Using this likelihood to obtain lensing amplitude summary statistics $\widehat{A}$, (with $\widehat{A}=1$ for ${\widehat{C}}_L^{\varphi \varphi }$ equal to the best-fit ΛCDM model to the Planck temperature and polarization power spectra and the reconstructed lensing power⁷), we obtain

$$\begin{array}{cc}\hfill {}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{8\to 400}^{\varphi ,\mathrm{MV}}& =1.011\pm 0.028\text{(CMB marginalized)}\hfill \end{array}$$(23)

over the conservative multipole range L = 8–400. This corresponds to a slightly higher value of the lensing spectrum than PL2015 (for which ${}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{40\to 400}^{\varphi ,\mathrm{MV}}=0.995\pm 0.026$) with a similar significance. The shift is mostly driven by the temperature reconstruction, whose amplitude is higher than 2015 by 0.8σ. This shift is consistent with that expected from the change in methodology and data: the bin 8 ≤ L ≤ 40 was not included in the 2015 lensing likelihood, and is around 1σ high in temperature (causing a 0.3σ amplitude shift), and the mask and SMICA weights have changed. We have evaluated the expected deviation for these two changes by comparison to reconstructions on the new SMICA maps with the 2015 mask, and on the new mask, using the 2015 SMICA weights. Using the observed shifts in the simulated reconstructions after the indicated changes, we find expected amplitude differences of 0.18σ and 0.33σ, respectively. Discarding any additional changes in the data processing, adding these in quadrature results in the observed total amplitude shift of 1.3σ. Over the aggressive multipole range, the measured amplitude is

$$\begin{array}{c}\hfill {}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{8\to 2048}^{\varphi ,\mathrm{MV}}=0.995\pm 0.026\text{(CMB marginalized)}.\end{array}$$(24)

As discussed in detail in Sect. 4, the high-L range fails a pair of consistency tests and we advise against using the full range for parameter constraints.

The temperature reconstruction still largely dominates our MV estimate, with amplitudes

$$\begin{array}{cc}& {}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{8\to 400}^{\varphi ,TT}=1.026\pm 0.035\text{(CMB marginalized)},\hfill \end{array}$$(25)

$$\begin{array}{cc}& {}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{8\to 2048}^{\varphi ,TT}=1.004\pm 0.033\text{(CMB marginalized)}.\hfill \end{array}$$(26)

Amplitude statistics for the polarization-only reconstructions are as follows (neglecting the CMB marginalization and other very small likelihood linear corrections):

$$\begin{array}{cc}& {}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{8\to 400}^{\varphi ,PP}=0.85\pm 0.16\text{(homogeneous filtering)};\hfill \end{array}$$(27)

$$\begin{array}{cc}& {}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{8\to 400}^{\varphi ,PP}=0.95\pm 0.11\text{(inhomogeneous filtering)}.\hfill \end{array}$$(28)

This is formally a 5σ measurement for our baseline filtering, and roughly 9σ with the optimized filtering. As can be seen in Fig. 4, the improvement of the polarization reconstruction is scale-dependent, with most gain achieved on small scales. This behaviour is consistent with analytic expectations, calculated using the independent-patch approximation introduced in Sect. 2.3.

All reconstruction band powers are consistent with a ΛCDM cosmology fit to the Planck CMB power spectra. Figure 5 presents a summary plot of our new MV band powers together with a compilation of other recent measurements, and the previous results from PL2015.

3.2. Likelihood and parameter constraints

We construct a lensing likelihood from the power-spectrum reconstruction following the same method as PL2015. We now expand the default conservative multipole range to 8 ≤ L ≤ 400, but exclude the higher multipoles to be conservative, given marginal evidence for null-test failures in the lensing curl and frequency consistency at L >  400. Multipoles L <  8 are very sensitive to the fidelity of the simulations due to the large mean field there, and we continue to exclude them for robustness (though only L = 2 looks clearly anomalous). The likelihood for the band powers over 8 ≤ L ≤ 400 is approximated as Gaussian, with a fixed covariance estimated from simulations, but the power-spectrum band powers are corrected perturbatively for changes in normalization and N⁽¹⁾ due to parameter-dependent deviations from the fiducial model. We neglect a possible dependence on cosmology of the small Monte Carlo correction.

The final likelihood is of the form⁸

$$\begin{array}{c}\hfill -2log{\mathcal{L}}_{\varphi }={\mathcal{B}}_i^L({\widehat{C}}_L^{\varphi \varphi }-C_L^{\varphi \varphi ,\mathrm{th}}){\left[{\mathrm{\Sigma }}^{-1}\right]}^{\mathit{ij}}{\mathcal{B}}_j^{L\prime }({\widehat{C}}_{L\prime }^{\varphi \varphi }-C_{L\prime }^{\varphi \varphi ,\mathrm{th}}),\end{array}$$(29)

where Σ is the covariance matrix and the binning functions ${\mathcal{B}}_i^L$ are defined in Eq. (11). The binned “theory” power spectrum for cosmological parameters θ is given in the linear approximation by

$$\begin{array}{c}\hfill {\mathcal{B}}_i^L{C}_L^{\varphi \varphi ,\mathrm{th}}\simeq {\mathcal{B}}_i^L{C_L^{\varphi \varphi }|}_{\mathit{\theta }}+M_i^{a,\ell \prime }({C_{\ell \prime }^a|}_{\mathit{\theta }}-{C_{\ell \prime }^a|}_{\mathrm{fid}}),\end{array}$$(30)

where a sums over both the $C_L^{\varphi \varphi }$ and CMB power-spectra terms, and the linear correction matrix $M_i^{a,{\mathcal{l}}^{\prime }}$ can be pre-computed in the fiducial model. The linear correction accounts for the N⁽¹⁾ dependence on $C_L^{\varphi \varphi }$, and the dependence of the lensing response and N⁽¹⁾ on the CMB power spectra; explicitly,

$$\begin{array}{cc}& M_i^{\varphi ,L\prime }={\mathcal{B}}_i^L\frac{\partial }{\partial C_{L\prime }^{\varphi }}{\mathrm{\Delta }{C}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}|}_{\mathrm{N}1}\hfill \end{array}$$(31)

$$\begin{array}{cc}& M_i^{X,\ell \prime }={\mathcal{B}}_i^L\frac{\partial }{\partial C_{\ell \prime }^X}({\mathrm{\Delta }{C}_L^{{\widehat{\varphi }}_1{\widehat{\varphi }}_2}|}_{\mathrm{N}1}+ln\left({[{\mathcal{R}}_L^{\varphi }]}^2\right){C_L^{\varphi }|}_{\mathrm{fid}}),\hfill \end{array}$$(32)

where $C_L^{\varphi }$ derivatives are understood not to act on the lensing contribution to lensed power spectra, X is one of the CMB power spectra, and $C_{\mathcal{l}}^X$ derivatives do not act on the fiducial power spectra in the estimator weights.

The 2015 CMB likelihoods were based on an LFI polarization likelihood at low multipoles, but the new 2018 low-ℓ likelihood now uses the HFI data for the low-ℓ polarization and gives a constraint on the optical depth with a considerably smaller uncertainty (Planck Collaboration V 2020). Since uncertainty in the optical depth is the main limitation to inferring the perturbation power-spectrum amplitude from CMB power-spectrum measurements, this means that the 2018 CMB likelihoods can constrain amplitude-related parameters significantly better than in 2015, reducing the relative impact of the information coming from the lensing likelihood. However, it is still important to check consistency using the lensing power spectrum. The lensing spectrum also contains some shape information and can probe extensions to ΛCDM in some directions of parameter space that cannot be constrained with primary CMB power-spectrum measurements alone (e.g., due to the geometric degeneracy). The new 8 ≤ L <  40 bin only has a small impact on ΛCDM parameter constraints, but is valuable for some extended models. For example, some modified gravity models, or the presence of compensated isocurvature modes, can give substantial changes to the lensing spectrum at low multipoles as discussed in PCP18 and Planck Collaboration X (2020).

3.2.1. Constraints from lensing alone and comparison with CMB

To do a “lensing-only” analysis, in PL2015 we fixed the theoretical CMB power spectra, which are required for the linear correction in Eq. (30), to a ΛCDM fit to the CMB power-spectrum data. We now remove any dependence on the theoretical model of the CMB power spectra by marginalizing out the theoretical $C_{\mathcal{l}}^{\text{CMB}}$ by approximating their distribution as Gaussian, where up to a constant

$$\begin{array}{c}\hfill -2lnP(C^{\mathrm{CMB}}|{\widehat{C}}^{\mathrm{CMB}})\simeq (C^{\mathrm{CMB}}-{\widehat{C}}^{\mathrm{CMB}}){\mathrm{cov}}_{\mathrm{CMB}}^{-1}(C^{\mathrm{CMB}}-{\widehat{C}}^{\mathrm{CMB}}).\end{array}$$(33)

Here, C^CMB and ${\widehat{C}}^{\mathrm{CMB}}$ are vectors of CMB TT, TE, and EE power-spectrum values at each multipole, with ${\widehat{C}}^{\mathrm{CMB}}$ a data estimate of the CMB power spectra without foregrounds (or noise), which could be measured in various ways. The covariance matrix of the CMB power spectra is cov_CMB. Integrating out C^CMB, the likelihood then takes the form of Eq. (29) with covariance increased to account for the uncertainty in the CMB power spectra,

$$\begin{array}{c}\hfill {\overline{\mathrm{\Sigma }}}_{\mathit{ij}}={\mathrm{\Sigma }}_{\mathit{ij}}+M_i^{X,\ell }{\mathrm{cov}}_{\mathrm{CMB}}^{X\ell ;Y\ell \prime }{M}_j^{Y,\ell \prime },\end{array}$$(34)

and the theory spectrum shifted by the linear correction for the observed CMB power,

$$\begin{array}{cc}\hfill {\mathcal{B}}_i^L{\overline{C}}_L^{\varphi \varphi ,\mathrm{th}}& \simeq ({\mathcal{B}}_i^L+M_i^{\varphi ,L}){C_L^{\varphi \varphi }|}_{\mathit{\theta }}\hfill \\ \hfill & -M_i^{\varphi ,L}{C_L^{\varphi \varphi }|}_{\mathrm{fid}}+M_i^{X,\ell }({\widehat{C}}_{\ell }^X-{C_{\ell }^X|}_{\mathrm{fid}}),\hfill \end{array}$$(35)

where X, Y are summed only over CMB spectra. The combined term on the second line is now a constant, so the likelihood only depends on cosmological parameters via $C_L^{\varphi \varphi }{|}_{\mathit{\theta }}$. We evaluate the CMB power correction using the plik_lite band powers, which are calculated from the full plik high-ℓ likelihood by marginalizing over the foreground model without any further assumptions about cosmology (Planck Collaboration XI 2016; Planck Collaboration V 2020). To relate plik_lite bins to the $M_i^{X,{\mathcal{l}}^{\prime }}$ bins, we assume that the underlying CMB power spectra are represented only by modes that are smooth over Δℓ = 50. The plik_lite bandpower covariance cov_CMB is similarly used to calculate Eq. (34). The increase in the diagonal of the covariance is about 6% at its largest, and the linear correction shifts lensing amplitude estimates slightly compared to using a ΛCDM best fit. The shift is largely explained because, over the ℓ range that the lensing reconstruction is sensitive to, the CMB TT data are somewhat less sharply peaked than the ΛCDM model (which also shows up in a preference for the phenomenological lensing amplitude parameter A_L >  1 when fitting just CMB power-spectrum data, as discussed in Planck Collaboration VI 2020); smaller dC_ℓ/dℓ between the acoustic peaks leads to a smaller lensing signal response, so the theory model value ${\overline{C}}_L^{\varphi \varphi ,\mathrm{th}}$ is decreased (by approximately 1.5% compared to the ΛCDM best fit).

We follow PL2015 in adopting some weak priors for constraining parameters from the lensing likelihood without using the Planck CMB power-spectrum data. Specifically, we fix the optical depth to reionization to be τ = 0.055, put a prior on the spectral index of n_s = 0.96 ± 0.02, and limit the range of the reduced Hubble constant to 0.4 <  h <  1. We also place a prior on the baryon density of Ω_bh² = 0.0222 ± 0.0005, motivated by D/H measurements in quasar absorption-line systems combined with the predictions of big-bang nucleosynthesis (BBN)⁹. The exact choice of Ω_bh² prior has very little effect on lensing-only constraints, but the prior is useful to constrain the sound horizon (since this has a weak but important dependence on Ω_bh²) for joint combination with baryon oscillation (BAO) data. We adopt the same methodology and other priors as Planck Collaboration VI (2020; hereafter PCP18), using camb (Lewis et al. 2000) to calculate theoretical predictions with HMcode to correct for nonlinear growth (Mead et al. 2016). Our CosmoMC (Lewis 2013) parameter chains are available on the Planck Legacy Archive¹⁰, where for comparison we also provide alternative results with a different set of cosmological priors consistent with those used by the DES collaboration (DES Collaboration 2018b). Parameter limits, confidence contours and marginalized constraints are calculated from the chains using the GetDist package (Lewis 2019), following the same conventions as Planck Collaboration XIII (2016).

Figure 6 shows the lensing-only ΛCDM constraint on σ₈ and Ω_m. As discussed in detail in PL2015, the lensing data constrain a narrow band in the 3D σ₈ − Ω_m − H₀ parameter space, corresponding to

$$\begin{array}{c}\hfill \frac{{\sigma }_8}{0.8}{\left(\frac{h}{0.67}\right)}^{-1}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{0.3}\right)}^{-0.27}=0.999\pm 0.026(68\%,\mathrm{lensing}\mathrm{only}),\end{array}$$(36)

or the tighter and slightly less prior-dependent 2% constraint

$$\begin{array}{c}\hfill \frac{{\sigma }_8}{0.8}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}}}{0.3}\right)}^{0.23}{\left(\frac{{\mathrm{\Omega }}_{\mathrm{m}}{h}^2}{0.13}\right)}^{-0.32}=0.986\pm 0.020(68\%,\mathrm{lensing}\mathrm{only}).\end{array}$$(37)

The allowed region projects into a band in the Ω_m–σ₈ plane with

$$\begin{array}{c}\hfill {\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}=0.589\pm 0.020(68\%,\mathrm{lensing}\mathrm{only}).\end{array}$$(38)

The corresponding result using a fixed CMB fit for the CMB power spectrum is ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}=0.586\pm 0.020$, which is consistent with the similar constraint, ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}=0.591\pm 0.021$, found in PL2015. The roughly 0.25σ shift down in this parameter is consistent with the slight increase in ${\widehat{A}}^{\varphi }$ because of the anti-correlation of ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}$ with the lensing deflection power, as discussed in PL2015. While the tight three-parameter constraints of Eqs. (36) and (37) depend on our priors (for example weakening by a factor of 2–4 if the baryon density prior and other priors are weakened substantially) the 2D projection of Eq. (38) is much more stable (see Table 2 for examples of prior sensitivity). The Planck 2018 power-spectrum constraints give slightly lower values of σ₈ compared to the 2015 analysis due to the lower optical depth, which increases the overlap between the lensing-only and CMB power-spectrum contours, making them very consistent within the ΛCDM model.

Combining CMB lensing with BAO data (Beutler et al. 2011; Ross et al. 2015; Alam et al. 2017), and recalling that we are placing a prior on Ω_bh² so that the sound horizon is fairly well constrained, we can break the main degeneracy and constrain individual parameters, giving the ΛCDM constraints

$$\begin{array}{c}\hfill \begin{array}{c}{H}_0=67.9_{-1.3}^{+1.2}\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}\hfill \\ {\sigma }_8=0.811\pm 0.019\hfill \\ {\mathrm{\Omega }}_{\mathrm{m}}=0.{303}_{-0.018}^{+0.016}\hfill \end{array}\}68\%,\text{lensing+BAO}.\end{array}$$(39)

The value of the Hubble constant inferred here assuming ΛCDM is in good agreement with other inverse distance-ladder measurements (Aubourg et al. 2015; DES Collaboration 2018a), and the ΛCDM result from Planck power spectra in PCP18, but is somewhat in tension with (i.e., lower than) more model-independent values obtained using distance-ladder measurements (Riess et al. 2018).

Massive neutrinos suppress the growth of structure on scales smaller than the neutrino free-streaming scale. The combination of CMB lensing and BAO data is expected to be a particularly clean way to measure the absolute neutrino mass scale via this effect. Allowing for a varying neutrino mass, the constraints from lensing with BAO are very broad and peak away from the base ∑m_ν = 0.06 eV we assumed for ΛCDM, though not at a significant level (see Table 2). Remaining degeneracies can be broken by using the acoustic-scale measurement from the Planck CMB power spectra. The acoustic scale parameter θ_*, the ratio of the sound horizon at recombination to the angular diameter distance, is very robustly measured almost independently of the cosmological model (since many acoustic peaks are measured by Planck at high precision). Using θ_* is equivalent to using an additional high-precision BAO measurement at the recombination redshift. For convenience we use the θ_MC parameter, which is an accurate approximation to θ_*, and conservatively take 100θ_MC = 1.0409 ± 0.0006 (consistent with the Planck data in a wide range of non-ΛCDM models). Using this, we have a neutrino mass constraint based only on lensing and geometric measurements combined with our priors:

$$\begin{array}{c}\hfill {\displaystyle \sum m_{\nu }<0.60\mathrm{eV}(95\%,\text{lensing+BAO+}{\theta }_{\mathrm{MC}}).}\end{array}$$(40)

The other parameters determining the background geometry are very tightly constrained by the inverse distance ladder, and the amplitude parameter is still well measured, though with lower mean value, due to the effect of neutrinos suppressing structure growth:

$$\begin{array}{c}\hfill \begin{array}{c}{H}_0=(67.4\pm 0.8)\mathrm{km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}\hfill \\ {\sigma }_8=0.{786}_{-0.023}^{+0.028}\hfill \\ {\mathrm{\Omega }}_{\mathrm{m}}=0.306\pm 0.009\hfill \end{array}\}68\%,\text{lensing+BAO+}{\theta }_{\mathrm{MC}}.\end{array}$$(41)

3.2.2. Joint Planck parameter constraints

The CMB lensing power spectrum is consistent with expectations from the Planck CMB power-spectrum results, and combining the two can tighten constraints on the amplitude parameters and geometrical parameters that are limited by geometrical degeneracies when measured from the anisotropy power spectrum alone. Figure 7 shows the combined constraint in the σ₈–Ω_m plane, giving the tight (sub-percent) amplitude ΛCDM result

$$\begin{array}{c}\hfill {\sigma }_8=0.811\pm 0.006(68\%,\mathit{Planck}\mathrm{TT},\mathrm{TE},\mathrm{EE}+\text{low}\mathrm{E}+\text{lensing}).\end{array}$$(42)

This result uses the temperature and E-mode polarization CMB power spectrum likelihoods, including both at low-ℓ, which we denote by Planck TT,TE,EE+lowE following PCP18. Constraints on some other parameters, also in combination with BAO, are shown in Table 2; for many more results and further discussion see PCP18. As in the previous analyses, the Planck high-ℓ CMB power spectra continue to prefer larger fluctuation amplitudes (related to the continuing preference for high A_L at a level nearing 3σ; see PCP18), with the low-ℓ optical depth constraint from E-mode polarization pulling the amplitude back to values that are more consistent with the lensing analyses. Even with the low-ℓ polarization, the power-spectrum data prefer around 1σ higher σ₈ than the lensing data alone, with the joint constraint lying in between. Values of σ₈ inferred from the CMB power-spectrum data in a given theoretical model cannot be arbitrarily low because the reionization optical depth cannot be arbitrarily small; observations of high-redshift quasars (Fan et al. 2006) indicate that reionization was largely complete by redshift z ≃ 6.5, and this additional constraint is shown in the dotted line in Fig. 7.

3.2.3. Joint CMB-lensing and galaxy-lensing constraints

Cosmic shear of galaxies can be used to measure the lensing potential with lower-redshift sources than the CMB. Since the source galaxies and lines of sight to the CMB partly overlap, in general the signals are correlated due to both correlated lensing shear and the intrinsic alignment of the source galaxies in the tidal shear field probed by CMB lensing. Since our CMB lensing reconstruction covers approximately 70% of the sky, the area will overlap with most surveys. The cross-correlation signal has been detected with a variety of lensing data (Hand et al. 2015; Liu & Hill 2015; Kirk et al. 2016; Harnois-Déraps et al. 2017; Singh et al. 2017) and may ultimately be a useful way to improve parameter constraints and constrain galaxy-lensing systematic effects (Vallinotto 2012; Das et al. 2013; Larsen & Challinor 2016; Schaan et al. 2017).

Here we do not study the correlation directly, but simply consider constraints from combining the CMB lensing likelihood with the cosmic shear likelihood from the Dark Energy Survey (DES; Troxel et al. 2018). In principle the likelihoods are not independent because of the cross-correlation; however, since the fractional overlap of the full CMB lensing map with the DES survey area is relatively small, and since the Planck lensing reconstruction is noise dominated on most scales, the correlation should be a small correction in practice and we neglect it here. We use the DES lensing (cosmic shear) likelihood, data cuts, nuisance parameters, and nuisance parameter priors as described by Troxel et al. (2018), DES Collaboration (2018b), and Krause et al. (2017). However we use cosmological parameter priors consistent with our own CMB lensing-only analysis described in Sect. 3.2.1¹¹.

Figure 8 shows the Planck and DES lensing-only ΛCDM constraints in the Ω_m–σ₈ plane, together with the joint constraint, compared to the result from the Planck CMB power spectra. The DES lensing constraint is of comparable statistical power to CMB lensing, but due to the significantly lower mean source redshift the degeneracy directions are different (with DES cosmic shear approximately constraining ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.5}$ and CMB lensing constraining ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}$). The combination of the two lensing results therefore breaks a large part of the degeneracy, giving a substantially tighter constraint than either alone. The lensing results separately, and jointly, are both consistent with the main Planck power-spectrum results. The joint result in the Ω_m–σ₈ plane constrains the combined direction

$$\begin{array}{c}\hfill {\sigma }_8{({\mathrm{\Omega }}_{\mathrm{m}}/0.3)}^{0.35}=0.{798}_{-0.019}^{+0.024}[68\%,(\mathrm{DES}+Planck)\mathrm{lensing}],\end{array}$$(43)

although the posterior is not very Gaussian due to a non-trivial DES lensing-intrinsic-alignment parameter degeneracy. If instead we adopt the cosmological parameter priors of Troxel et al. (2018), but fixing the neutrino mass, then the lensing-only joint result is more Gaussian, with ${\sigma }_8{({\mathrm{\Omega }}_{\mathrm{m}}/0.3)}^{0.4}=0.{797}_{-0.018}^{+0.022}$; this is tighter than the constraint obtained by Troxel et al. (2018) in combination with galaxy clustering data¹².

Cosmological parameter degeneracies (and degeneracies with intrinsic-alignment and other nuisance parameters) limit the precision of the DES lensing-only results. Results can be tightened by using different priors (as in the DES analysis DES Collaboration 2018b, see also the analysis in PCP18 using the DES priors). Adding additional data also substantially reduces the degeneracy. For example, adding BAO data gives the tighter contours shown in Fig. 9. DES lensing+ BAO gives the marginalized constraints on individual parameters

$$\begin{array}{c}\hfill \begin{array}{c}{H}_0=70.8_{-2.8}^{+2.1}{\mathrm{km}\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}\hfill \\ {\sigma }_8=0.728\pm 0.052\hfill \\ {\mathrm{\Omega }}_{\mathrm{m}}=0.{348}_{-0.041}^{+0.034}\hfill \end{array}\}68\%,\text{DES lensing+BAO},\end{array}$$(44)

also consistent with the main Planck power-spectrum parameter analysis. However, as shown in Fig. 9, even after combining with BAO, the DES lensing results have a substantial remaining σ₈–Ω_m–H₀ degeneracy that limits the individual constraints. For CMB lensing the situation is much better, since the CMB lensing comes from a single well-defined source redshift plane that leads to a tight σ₈–Ω_m–H₀ degeneracy that intersects with the BAO Ω_m–H₀ degeneracy in a much narrower region of parameter space (giving the parameter constraints of Eq. (39)). Adding CMB lensing data to the DES result therefore gives much tighter constraints on individual parameters:

$$\begin{array}{c}\hfill \begin{array}{c}{H}_0=(67.6\pm 1.0){\mathrm{km}\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}\hfill \\ {\sigma }_8=0.805\pm 0.014\hfill \\ {\mathrm{\Omega }}_{\mathrm{m}}=0.295\pm 0.011\hfill \end{array}\}\begin{array}{c}68\%,\text{DES lensing}\hfill \\ +\mathit{Planck}\text{lensing+BAO}.\hfill \end{array}\end{array}$$(45)

The constraining power on individual parameters is dominated by CMB lensing+BAO, but the DES lensing data do help partly to break the remaining CMB lensing degeneracy, giving tighter constraints than from CMB lensing+BAO alone. The tight joint constraint on the Hubble constant here is consistent with the result of DES Collaboration (2018a), but without the use of galaxy clustering data (which are sensitive to details of bias modeling). Our combined result with the full DES joint likelihood, including galaxy clustering and galaxy–galaxy lensing (DES Collaboration 2018b), is given in Table 2.

3.2.4. Parameters from likelihood variations

The baseline conservative likelihood, restricted to 8 ≤ L ≤ 400, was chosen to be robust; however, we cannot rule out the possibility that the moderate null-test failures at higher L that motivate our choice L ≤ 400 are simply statistical fluctuations, in which case it would be perfectly legitimate to use the full multipole range for parameter constraints. Comparing different likelihood variations also allows us to assess how changes in the spectrum propagate into changes in parameters, and hence robustness of the parameter constraints.

Table 2 compares parameter constraints from the baseline likelihood with a number of variations with different binning, multipole ranges, and with and without using polarization in the reconstruction. As a general trend the higher multipoles are below the best-fit to the conservative range, so including lensing multipoles up L ≤ 2048 tends to pull amplitude parameters to lower values, although for the MV reconstruction only to an extent that is compatible with expected shifts when including more data. The shift between TT and MV results is, however, more anomalous: over the full multipole range the mean value of ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}$ shifts by more than 1σ despite the error bar only shrinking by a bit more than 10% (assuming Gaussian statistics, the shift is greater than 2σ unusual). Over the conservative multipole range the shift is more acceptable, although still somewhat conspicuous. Even over the conservative range the band powers are slightly tilted with respect to a ΛCDM best fit, which gives a mild (1–2)σ preference for non-zero neutrino mass when combined with BAO. We discuss the consistency of the different data ranges in more detail at the data level in Sect. 4 below.

3.3. Joint CIB–CMB lensing potential reconstruction

The Planck lensing reconstruction has S/N peaking at unity around the peak of the deflection-angle power spectrum at L ≃ 60, but is noise dominated on smaller scales. However, the cosmic infrared background (CIB) is a high-redshift tracer of the matter distribution and known to be correlated at the roughly 80% level to the CMB lensing potential and hence, potentially, is a good proxy for it (Song et al. 2003; Planck Collaboration XVIII 2014; Sherwin & Schmittfull 2015). The CIB is well measured by Planck’s higher frequency channels (Planck Collaboration XVIII 2011; Planck Collaboration XXX 2014; Planck Collaboration Int. XLVIII 2016; Mak et al. 2017), with high S/N to much smaller scales than probed directly by the lensing reconstruction. It has already been demonstrated that Planck’s CIB measurement can be used to delens the CMB acoustic peaks with about the same efficiency as Planck’s internal measurement (Larsen et al. 2016; Carron et al. 2017), with the CIB carrying substantially more information on small scales. For the purpose of delensing degree-scale B modes, most of the lensing signal required is from lensing multipoles L ≃ 500, where the Planck MV lensing reconstruction map is fully noise dominated, making the CIB especially valuable until higher-sensitivity internal CMB-polarization-based reconstructions are available. The main difficulty in using the CIB as a tracer is contamination by Galactic dust and modeling of the cross-correlation coefficient. The Galactic dust is more of a problem on large scales, just where the Planck lensing reconstruction S/N peaks. This suggests that a joint analysis could potentially give a substantial improvement to the lensing potential determination, and hence also improve the efficiency of delensing based on Planck data (Sherwin & Schmittfull 2015). A combined tracer map, including galaxy number counts from WISE in addition to the CIB and the 2015 Planck lensing map, has recently been presented over 43% of the sky by Yu et al. (2017). Here we focus on the Planck-based combination of the CIB and lensing reconstructions on 60% of the sky, with similar maximum delensing efficiency, and use them to delens all four Planck CMB spectra.

We combine the (internal) MV quadratic estimator reconstruction with the CIB map provided by the GNILC (Generalized Needlet Internal Linear Combination) component-separation method at 353, 545, and 857 GHz (Planck Collaboration Int. XLVIII 2016, based on the 2015 Planck data release and available at the Planck Legacy Archive). We use a mask given by the union of the lensing mask and the GNILC mask, leaving 60% of the sky unmasked (58% after apodization). Given the strong coherence of the CIB across this range of frequencies, we see very little gain in using combinations of the three frequency maps, and our baseline results use the 545-GHz CIB map only. Figure 10 shows the cross-spectra of the lensing temperature-only, polarization-only and MV reconstructions with the GNILC 545-GHz map. The bispectrum between residual foregrounds in the CMB temperature maps that enter the quadratic estimator and the CIB can potentially bias the CIB-lensing cross-spectra. Figure 10 shows two tests of such bias, neither of which shows any evidence for systematic contamination. The first performs lens reconstruction with a SMICA CMB map constructed to deproject the thermal-SZ effect (i.e., the weighting across frequencies nulls any component with the frequency spectrum of the thermal SZ effect). The second uses only polarization in the lensing reconstruction, which is expected to be essentially free of extragalactic foregrounds, at the cost of significantly larger reconstruction noise for Planck. Contamination from residual CIB in the SMICA temperature map is assessed below.

We build lensing estimates as follows. Consider a set of tracers Iⁱ of the lensing potential, with auto-power spectra $C_L^{I_i{I}_i}$, cross-correlation coefficient matrix ρ_L, and cross-correlation coefficients to the true lensing potential ${\rho }_L^{i\varphi }$. Assuming that we can treat the lensing and tracer fields as being approximately jointly Gaussian, it is straightforward to obtain a maximum a posteriori (MAP) estimate for the true lensing given the tracers, yielding an optimally-filtered potential map ${\widehat{\varphi }}^{\mathrm{MAP}}$. We can write this in terms of (isotropic) weights w_L acting in harmonic space on renormalized maps with unit spectra, through

$$\begin{array}{c}\hfill \frac{{\widehat{\varphi }}_{\mathit{LM}}^{\mathrm{MAP}}}{\sqrt{C_L^{\varphi \varphi ,\mathrm{fid}}}}\equiv {\displaystyle \sum _i{w}_L^i\frac{I_{\mathit{LM}}^i}{\sqrt{C_L^{I_i{I}_i}}},}\end{array}$$(46)

where the optimal weights are given at each lensing multipole by

$$\begin{array}{c}\hfill w_L^j={\displaystyle \sum _i{\rho }_L^{\varphi i}{({\mathit{\rho }}_L^{-1})}_{\mathit{ij}}.}\end{array}$$(47)

The expected resulting squared cross-correlation coefficient ${\rho }_L^2$ is the weighted sum ${\sum }_i{w}_L^i{\rho }_L^{\varphi i}$. To compare the noise in the joint reconstruction to the quadratic-estimator reconstruction noise $N_L^{(0)}$, we can also define an effective noise level N_L so that

$$\begin{array}{c}\hfill \frac{C_L^{\varphi \varphi ,\mathrm{fid}}}{C_L^{\varphi \varphi ,\mathrm{fid}}+N_L}={\rho }_L^2={\displaystyle \sum _i{w}_L^i{\rho }_L^{\varphi i}.}\end{array}$$(48)

We adopt a purely empirical approach to obtain the weights, together with the auto-spectra on the right-hand side of Eq. (46), and take $C_L^{\varphi \varphi ,\mathrm{fid}}$ as our FFP10 fiducial lensing spectrum. For each pair of tracers, we perform a bicubic spline fit across 12 bins of the observed cross-correlation coefficient and auto-spectra, with errors for each bin determined from the observed scatter. The spectra and cross-spectra entering the cross-correlation coefficient are calculated by deconvolving the effects of the mask from pseudo-C_ℓ spectra (Wandelt et al. 2001), using a 12′ apodization window. To estimate the cross-correlation of the tracers to the true lensing potential, we assume that the quadratic estimator is unbiased (which, according to the FFP10 simulations, is a very good approximation except at the very lowest multipoles); the cross-spectra of the CIB maps to the quadratic estimator are then used as proxies to the true C^CIB ϕ at the corresponding frequency.

These CIB-lensing cross-spectra get a contribution from the CIB bispectrum from residual CIB present in the CMB maps. However, this is very small: Sect. 4.5 discusses foreground contamination to our lensing estimates more generally using a dedicated simulation set. Using these simulations, we can cross-correlate the quadratic estimator applied to the expected residual CIB in the SMICA maps to the CIB component at 545 GHz. We find a bias of at most 1% at ℓ ≃ 2000, which we can safely neglect when building the tracers.

The cross-correlation coefficients ${\rho }_L^{\varphi \mathrm{GNILC}}$ built in this way have significant noise, but can barely be distinguished from each other across frequencies. We choose to use the same fiducial ρ^ϕCIB at all frequencies, which we obtain by averaging the estimates at 353 and 545 GHz. Figure 11 shows (smoothed) cross-correlation coefficients of the true lensing map to the MV estimator (blue), together with that of the 545-GHz GNILC map (orange), and their combination (green). Delensing efficiencies, ${ϵ}_L\equiv {\rho }_L^2$, are a good measure of how much lensing can be removed at each lensing multipole if the maps are used for a delensing analysis (see Sect. 3.4). Effective noise levels from Eq. (48) are shown in the lower plot for the MV reconstruction, the 545-GHz GNILC lensing estimate, and their combination, showing how the CIB particularly improves the S/N at smaller scales.

We apply the weights given by Eq. (47) to the pseudo-a_LMs of the tracers, obtained using the apodized mask, keeping afterwards the multipole range 10 ≤ L ≤ 2000. We do not perform any further low-L cuts to the GNILC maps, since, for our purposes at least, they do not display obvious signs of dust contamination, and they make only a small contribution to the combination with the internal lensing reconstruction on large scales. We note that our delensing analysis below is robust to the low-L quadratic estimator mean field, as demonstrated by Carron et al. (2017).

Our joint CIB and internal lensing reconstruction maps are shown in Fig. 12. As expected, adding the CIB tracer dramatically improves the resolution of the lensing reconstruction, effectively filling in the small-scale modes by using the phase information about the sky realization from the CIB, with amplitude determined from cross-correlation of the CIB with the lensing reconstruction. The joint map is available in the 2018 Planck public release. It could be used for delensing future low-ℓ polarization data, and is independent of the CMB at ℓ <  100 (where we cut the input CMB maps), except for possibly a small contamination in the CIB map. We have also constructed a version of the map that uses no input B modes, to avoid delensing biases (see Sect. 3.4.1 below for further discussion).

The joint CIB-lensing map is an excellent tracer of the realization of the lensing field; however, the lensing amplitude information that it contains all comes from the lensing reconstruction, since we do not have an accurate predictive model for the CIB signal. It cannot therefore be used directly to improve parameter constraints from its power spectrum.

3.4. Delensing Planck power spectra

CMB B-mode polarization from lensing is present on all scales, and for small tensor-to-scalar ratios could become an important source of confusion for the signal from primordial gravitational waves. However, using a lensing reconstruction it is possible to estimate and subtract most of the B-mode signal, a process known as delensing (Knox & Song 2002; Kesden et al. 2002; Hirata & Seljak 2003). This may ultimately be crucial for future observations to detect low levels of primordial gravitational waves. Delensing can also be applied to the temperature and E-mode polarization maps, so that the corresponding peak-smoothing effect in the power spectra can be partly removed, sharpening the acoustic peaks and restoring more of the information to the power spectrum (Green et al. 2017). Peak sharpening by delensing has been convincingly demonstrated by Larsen et al. (2016, using CIB maps as a lensing tracer) and Carron et al. (2017, using the internal Planck 2015 lensing reconstruction). Carron et al. (2017) and Manzotti et al. (2017) (the latter using the CIB as measured by Herschel at 500 μm to delens SPT data) have also successfully demonstrated that delensing can reduce the power in the B-mode polarization, although currently instrumental noise rather than lensing is still limiting the detectable level of primordial B modes.

The three black lines in the upper panel of Fig. 11 show the leading eigenvector of the matrix $\mathcal{l}\partial C_{\mathcal{l}}^{TT,EE,BB}/\partial \ln C_L^{\varphi \varphi }$, for 100 ≤ ℓ ≤ 2048, showing the scales over which lensing modes are relevant for lensing of the power spectra. Peak sharpening in T, E requires good delensing efficiency at lensing multipoles L ≲ 250, while removal of BB lensing power requires smaller-scale lensing reconstruction. As discussed in Sect. 3.3, a joint CIB/internal delensing analysis is expected to be significantly better for delensing all signals because of the complementarity of scales.

For characterization of our delensing analysis, we need to extend the FFP10 simulation suite to include CIB components. We perform this in the simplest manner, simulating the CIB as an isotropic Gaussian field, obeying our estimates of the various cross-correlations. To produce the CIB tracers I_i at each frequency ν_i, we first rescale the FFP10 input lensing potential harmonic coefficients according to ρ^ϕi, then add independent Gaussian noise ϵ_i:

$$\begin{array}{c}\hfill I_{\mathit{LM}}^i=\sqrt{\frac{{\widehat{C}}_L^{I^i{I}^i}}{C_L^{\varphi \varphi ,\mathrm{f}id}}}{\rho }_L^{\varphi i}{\varphi }_{\mathit{LM}}^{\mathrm{input}}+\sqrt{{\widehat{C}}_L^{I^i{I}^i}}{ϵ}_{\mathit{LM}}^i,\end{array}$$(49)

where the noise is generated according to the reduced covariance

$$\begin{array}{c}\hfill \langle {ϵ}_{\mathit{LM}}^i{ϵ}_{\mathit{LM}}^{j\ast }\rangle ={\rho }_L^{\mathit{ij}}-{\rho }_L^{\varphi i}{\rho }_L^{\varphi j}.\end{array}$$(50)

After convolving with the GNILC 5′ FWHM Gaussian beam, the maps are masked and analysed in the same way as the data maps. We separately use a simulation suite that includes no lensing effects, but are generated with lensed CMB spectra. In this case the simulated CIB maps are just the noise components described by Eq. (49), with ${\rho }_L^{\varphi i}=0$. The simulated harmonic coefficients are conservatively computed over a range of scales ℓ_max ≤ 2500, slightly larger than the one used for the data combination to avoid edge effects due to masking.

We first delens the BB power spectrum using a template subtraction method in Sect. 3.4.1. We then use a debiasing technique to isolate the delensing signature on the CMB signal for the TT, TE, EE and BB spectra, with the help of a direct remapping method in Sect. 3.4.2.

3.4.1. C${}_{\ell }^{\mathit{BB}}$-delensing

Using a filtered version of the E-mode map, and a tracer of the lensing ϕ field, we can build a template for the lensed B modes that we can subtract from the data to reduce the B-mode power. To reduce the B-mode power as much as possible, the fields must be Wiener filtered (Smith et al. 2009, 2012; Sherwin & Schmittfull 2015). The extent by which we can then reduce the B-mode power depends on the quality of the tracer and the noise in the E-mode map. The correlations of the lensing tracer maps to the true lensing field are shown in Fig. 11. For the relevant lensing multipoles (L ≃ 500) about 0.6², that is, 35% of the lensing can be removed using GNILC, with a moderate improvement in combination with the Planck internal reconstruction. The Planck small-scale polarization data are noisy, and at CMB E-mode scale ℓ ≃ 500, from where the large-scale B-power takes significant contributions, the E-mode filter quality $C_{\mathcal{l}}^{EE}/(C_{\mathcal{l}}^{EE}+N_{\mathcal{l}}^{EE}/b_{\mathcal{l}}^2)$ (where b_ℓ is the combined beam and pixel window function) is no greater than 0.6, so we cannot expect to remove more than about 20% of the lensing B-mode power. This is in contrast to current ground-based experiments, for which the much lower polarization noise level makes the lensing map fidelity the main limiting factor, as demonstrated by Manzotti et al. (2017).

We now describe our template in more detail. We start with the Wiener-filtered E-mode map, the same as produced by our pipeline for input to the quadratic estimators (E^WF in Eq. (2)). We use the inhomogeneously-filtered version, since this visibly increases the local fidelity of the template and our delensing efficiencies. We build the polarization (Stokes parameters) P^(E) ≡ Q^(E) + iU^(E) from this E-mode map, and simply remap the polarization according to the deflection angle determined by the lensing tracer, that is,

$$\begin{array}{c}\hfill P^{\mathrm{template}}(\widehat{\mathit{n}})=P^{(E)}(\widehat{\mathit{n}}+\mathbf{\nabla }{\widehat{\varphi }}_{\mathrm{WF}}).\end{array}$$(51)

After projecting into the B-mode component, this forms our template for the lensed B modes. The simplified notation $\widehat{\mathit{n}}+\mathbf{\nabla }{\widehat{\varphi }}_{\mathrm{WF}}$ denotes the displacement of length $|\mathbf{\nabla }{\widehat{\varphi }}_{\mathrm{WF}}|$ on the sphere along the geodesic defined by $\mathbf{\nabla }{\widehat{\varphi }}_{\mathrm{WF}}$, including the small rotation of the spherical-polar polarization components induced by the parallel transport of the polarization tensor (Challinor & Chon 2002; Lewis 2005). We perform the interpolation using the Python curved-sky lensing tools from LENSIT¹³, with a bicubic spline interpolation algorithm. The observed E-mode map we use to build the template is itself lensed, while ideally one would use an unlensed E-mode map. In principle, we could try and improve on that by obtaining an optimally-filtered E map including our deflection tracer in the Wiener filtering, using the algorithm developed by Carron & Lewis (2017); however, the difference is a second-order effect in the deflection, and we neglect it here.

The cross-correlations of these templates to the B-mode CMB map provide strong indirect detections of the lensing B-mode power in the data. The cross-spectra are displayed in Fig. 13, for templates constructed from the lensing quadratic estimator (blue), the CIB at 545 GHz (orange), and in combination (green). These templates are cross-correlated to our B-mode map, obtained from the Wiener-filtered $B_{\ell m}^{\mathrm{WF}}$ after dividing by $C_{\mathcal{l}}^{BB,\text{fid}}/\left(C_{\mathcal{l}}^{BB,\text{fid}}+N_{\mathcal{l}}^{BB}/b_{\mathcal{l}}^2\right)$, where b_ℓ is the combined beam and pixel window function. This is an unbiased estimate of the B-mode map away from the mask boundaries. The black line shows the fiducial lensing B-power of roughly (5 μK-arcmin)² at low multipoles. We use the quadratic TT+TE+EE estimator, discarding the EB and TB parts, so that all tracers are statistically independent from the B-mode map. The raw cross-spectra actually provide only a filtered version of the B-mode power, since the Wiener-filtering of E and ϕ maps reduces the power in the template. We follow Hanson et al. (2013), assuming isotropy of the filtering, and simply rescale the cross-spectra by the ratio of the result calculated on the FFP10 simulation suite to the simulated B-mode power. In all cases, the data cross-correlation amplitudes are consistent with the fiducial model within 1σ, with the lensing B-mode power detected at 13, 18, and 20σ for templates constructed from the lensing quadratic estimator, the CIB, and in combination, respectively.

To delens, the lensed B-mode template is simply subtracted from our B-mode map estimate, obtained as described above. Figure 14 shows the result of this procedure for different choices of tracer, after building the delensed difference spectrum that partly cancels cosmic variance and noise:

$$\begin{array}{c}\hfill \mathrm{\Delta }{\widehat{C}}_{\ell }\equiv {\widehat{C}}_{\ell }^{\mathrm{del}}-{\widehat{C}}_{\ell }^{\mathrm{dat}}.\end{array}$$(52)

The spectra are obtained from the position-space B-mode maps before and after the B-mode template subtraction, using (scalar) deconvolution of the pseudo-C_ℓ spectra to account for the effect of the mask. Table 3 lists the results of fitting a single delensing efficiency parameter $\mathrm{\Delta }{C}_{\ell }^{\mathit{BB}}=-ϵC_{\ell }^{BB,\mathrm{fid}}$ across all the bins shown in Fig. 14. The efficiency defined in this way reaches 12% for the combined tracers. Large-scale B modes are better delensed, but with larger statistical errors. For instance, for 100 ≤ ℓ ≤ 300 we find a delensing efficiency of 0.217 ± 0.047 (${\widehat{\varphi }}^{\mathrm{MV}}$+CIB at 545 GHz), with a predicted value of 0.214, in accordance with the naive expectations laid out at the beginning of this section based on the E and ϕ map quality.

The presence of common B modes in the quadratic estimator reconstruction noise and the field being delensed would create a spurious delensing signature coming from leading disconnected correlators (Teng et al. 2011; Carron et al. 2017). As we did when estimating the lensing B-mode power spectrum above, we avoid this by using the quadratic estimator involving only TT+TE+EE, discarding the EB and TB parts. The lensing tracers used in Fig. 14 are then statistically independent from the B-mode map, at negligible cost in tracer fidelity.

3.4.2. Acoustic peak sharpening (de-smoothing)

We now test the expected acoustic peak re-sharpening of the T and E power spectra after delensing. For the case of the quadratic estimator, internal delensing biases originate from the statistical dependence between the lensing reconstruction noise and the CMB maps. Here, the biases are more difficult to avoid than for BB-delensing in the previous section, and we refer the reader to Carron et al. (2017) for a detailed discussion of these biases. Proposed methods to correct for such biases include using a non-overlapping set of modes (Sehgal et al. 2017), at the cost of some loss of S/N, or a much more elaborate realization-dependent higher-point function estimation (Namikawa 2017). We follow the procedure of Carron et al. (2017), which we summarize below, to correct for these biases and evaluate our delensing efficiencies. This procedure has no impact on the covariance of the biased delensed spectra. This covariance has not been well studied, and this procedure will probably not be optimal for ambitious future experiments aiming to transfer information from the lensing higher-point statistics back to the CMB power spectrum, as advocated by Green et al. (2017).

The procedure is as follows: using our lensing tracers, we remap an estimate of the CMB temperature and polarization fields, filtered in the multipole range 100 ≤ ℓ ≤ 2048, according to the estimate of the inverse displacement $-\mathbf{\nabla }{\widehat{\varphi }}_{\mathrm{WF}}$. The CMB maps are built from the inverse-variance, homogeneously-filtered T, E, and B CMB harmonic coefficients, rescaled by the isotropic limit of the filter (since we neglect the TE correlation in the filtering, these limits are simply given by $C_{\ell }^{TT,EE,BB,\mathrm{fid}}+N_{\ell }^{TT,EE,BB}/b_{\ell }^2$); the resulting filtered maps are expected to be unbiased away from the mask boundaries. As before, we calculate differences of delensed and raw T, E, and B spectra after deconvolution of their pseudo-C_ℓ power spectra for the effects of the mask. To assess the delensing effects and subtract the internal delensing biases, we subtract the same spectral differences obtained using the same delensing procedure, but on Gaussian simulations (generated using lensed CMB spectra) and uncorrelated CIB:

$$\begin{array}{c}\hfill B_{\ell }^{\mathrm{Gauss}}\equiv {\langle {\widehat{C}}_{\ell }^{\mathrm{del}}-{\widehat{C}}_{\ell }^{\mathrm{dat}}\rangle }_{\mathrm{Gauss}}.\end{array}$$(53)

The combination ${\widehat{C}}_{\ell }^{\mathrm{del}}-{\widehat{C}}_{\ell }^{\mathrm{dat}}-B_{\ell }^{\mathrm{Gauss}}$ should vanish in the mean by construction in the absence of lensing, and we use these spectra to quantify the significance of the delensing that we detect. However, to assess how much delensing was really performed on the CMB, further corrections are required. First, $B_{\mathcal{l}}^{\text{Gauss}}$ subtracts not only the spurious delensing bias, but also the CMB peak smoothing effect due to independent reconstruction noise of the tracer. Second, the spectral difference ${\widehat{C}}_{\ell }^{\mathrm{del}}-{\widehat{C}}_{\ell }^{\mathrm{dat}}-B_{\ell }^{\mathrm{Gauss}}$ still contains the difference of the noise spectra caused by the signal part of the tracer. To isolate the CMB peak-sharpening, we therefore build the combination

$$\begin{array}{c}\hfill \mathrm{\Delta }{\widehat{C}}_{\ell ,\mathrm{debias}}\equiv {\widehat{C}}_{\ell }^{\mathrm{del}}-{\widehat{C}}_{\ell }^{\mathrm{dat}}-B_{\ell }^{\mathrm{Gauss}}-B_{\ell }^{\mathrm{Noise}}+B_{\ell }^{\mathrm{CMB}},\end{array}$$(54)

with

$$\begin{array}{c}\hfill B_{\ell }^{\mathrm{Noise}}\equiv {\langle {\widehat{C}}_{\ell }^{\mathrm{del}}-{\widehat{C}}_{\ell }^{\mathrm{dat}}\rangle }_{\mathrm{signal}{\widehat{\varphi }}_{\mathrm{WF}}},\end{array}$$(55)

and

$$\begin{array}{c}\hfill B_{\ell }^{\mathrm{CMB}}\equiv {\langle {\widehat{C}}_{\ell }^{\mathrm{del}}-{\widehat{C}}_{\ell }^{\mathrm{dat}}\rangle }_{\mathrm{independent}{\widehat{\varphi }}_{\mathrm{WF}}}.\end{array}$$(56)

In Eq. (55), the lensing tracer only contains its signal part, and in Eq. (56) the lensing tracer is constructed using a Gaussian simulation independent of the one being delensed. Only the noise part of the simulation is considered in Eq. (55), and only the CMB part in Eq. (56).

Figure 15 shows the resulting spectral differences and Table 4 collects summary statistics. We list delensing efficiencies ϵ, obtained as a straightforward χ²-minimization of the observed binned spectral differences

$$\begin{array}{c}\hfill {\chi }^2(ϵ)={\displaystyle \sum _{\mathrm{bins}b}\frac{1}{{\sigma }_{{\ell }_{\mathrm{b}}}^2}{[\mathrm{\Delta }{\widehat{C}}_{{\ell }_{\mathrm{b}},\mathrm{debias}}-ϵ(C_{{\ell }_{\mathrm{b}}}^{\mathrm{unl},\mathrm{fid}}-C_{{\ell }_{\mathrm{b}}}^{\mathrm{len},\mathrm{fid}})]}^2,}\end{array}$$(57)

together with predictions obtained from our simulations and the best-fit reduced χ². Here, $C_{{\ell }_{\mathrm{b}}}^{\mathrm{unl},\mathrm{fid}}-C_{{\ell }_{\mathrm{b}}}^{\mathrm{len},\mathrm{fid}}$ is the binned difference of the unlensed and lensed power spectra in the fiducial model, and ${\sigma }_{{\mathcal{l}}_{\text{b}}}^2$ is the variance of the binned $\mathrm{\Delta }{\widehat{C}}_{{\ell }_{\mathrm{b}},\mathrm{debias}}$ estimated from simulations. We detect delensing at high significance in the difference between the delensed and raw spectra, which removes most of the cosmic variance and noise common to both. When combining the tracers we also list the improvement achieved (defined with respect to the most powerful single tracer). The improvement is always significant except for the B-mode power, where the quadratic estimator provides only a little additional information on the required scales.

The last column of Table 4 and the lower panel of Fig. 14 show the delensing of the B-mode power with the remapping method. The procedure we follow in this section, using Eq. (54), is designed isolate the change in the CMB signal power spectrum due to the delensing remapping. This is not, by construction, the same as the actual change in the power spectrum of the maps: for example, for B-mode delensing, remapping by the reconstruction noise producesB modes, and remapping the instrument noise also changes the B-mode power. The delensing efficiency figures for BB that are quoted in Table 4 are therefore not directly comparable to those reported in Table 3 for the template-delensing method, which are based on the actual change in B-mode power after delensing. The template method is also optimized to account for instrument noise (through the Wiener-filtering of the E-mode map in the construction of the template). The template-based efficiencies are, therefore, the relevant ones for assessing improvements in primordial B-mode limits from delensing. However, since the PlanckB-mode measurement is noise dominated, there would be a negligible improvement in primordial B-mode limits due to the small decrease in the lensing component.

The peak sharpening in the other CMB spectra could in principle slightly improve parameter constraints (Green et al. 2017); however, the likelihood model is complicated and the forecast improvement for Planck is very small, so we do not use the delensed spectra for cosmological parameter analyses. Our demonstration of peak sharpening is an important proof of principle, and also a useful consistency check on the lensing analysis.

4. Null and consistency tests

We now turn to various consistency tests of the lensing and lensing curl reconstruction band powers. In Sect. 4.1 we first discuss the empirical distribution of the band powers in simulations and the presence of features. In Sect. 4.2 we focus on the temperature reconstruction, and present detailed consistency tests based on comparisons of different reconstructions; Sect. 4.3 then compares results from various polarization estimators with the temperature reconstruction. In Sect. 4.4 we present cross-half-mission data splits and noise null tests. In Sect. 4.5 we test for the impact of correlations between the foregrounds and lensing signal. Our entire analysis uses the FFP10 cosmology as the fiducial model; we test in Sect. 4.6 that this does not affect our results. Section 4.7 discusses the impact of the non-Gaussianity of the lensing field, and Sect. 4.8 presents alternative calculations of the N⁽¹⁾ bias.

4.1. Band-power distribution and features

We first discuss the empirical statistical behaviour of our baseline data band powers compared to the result expected from FFP10 fiducial simulations. PL2015 already hinted at two possible features in the temperature-only lensing reconstruction (a sharp dip in the lensing spectrum at L ≃ 700, and a low curl spectrum over a broad range of scales). These two features are still present, as displayed in Fig. 16. The upper panel shows the relative deviation of ${\widehat{C}}_L^{\varphi \varphi }$ to the Planck TT,TE,EE+lowE+lensing best-fit cosmology and the lower panel the curl spectrum, both on the aggressive multipole range.

The dip in the measured ${\widehat{C}}_L^{\varphi \varphi }$ band powers in the upper panel of Fig. 16 is slightly less pronounced than in the previous release; the bin covering 638 ≤ L ≤ 762 shows now a 2.8σ deviation from the best-fit spectrum, compared to 3.8σ in 2015 for the same best-fit spectrum. This significance is unchanged after marginalization of the CMB spectra uncertainties. Despite the large errors at high-L, the measurement can affect results if these multipoles are used for parameter inference. Can the dip be a statistical fluke? Empirically, we find that 5% of our simulated temperature reconstructions show equal or stronger outliers (the most extreme TT outlier across all aggressive bins and simulations is −4σ at L ≃ 50). The feature in the data is therefore consistent with a fluctuation that is only slightly unusual. A statistical fluctuation is also consistent with Fig. 10, which shows that the cross-spectrum with the GNILC CIB tracer has no significant dip at these multipoles.

The lower panel of Fig. 16 shows our lensing curl reconstruction, which is expected to be zero to within tiny post-Born corrections that can be safely neglected at Planck noise levels (Pratten & Lewis 2016; Fabbian et al. 2018). Our curl amplitudes ${}_{\mathrm{fid}}{\widehat{A}}^{\mathrm{\Omega },\mathit{TT}}$ are built using the same binning procedure as for the lensing gradient in Eq. (10), using a reference flat spectrum

$$\begin{array}{c}\hfill \frac{L^2{(L+1)}^2}{2\pi }{C}_L^{\mathrm{\Omega }\mathrm{\Omega },\mathrm{f}id}\equiv {10}^{-7}.\end{array}$$(58)

Orange points use a coarser binning and emphasise the low curl found over a fairly wide range of scales at 300 ≤ L ≤ 900, very similar to the result found in PL2015. None of our original aggressive band powers (blue) are anomalous individually, but the resulting coarse central bin lies formally 4.2σ away from zero. This binning was custom-made in order to emphasise the feature, and any assessment of its true significance must take this into account. We have done so as follows: for each of our FFP10 simulations, we identified within the aggressive binning the longest sequence of adjacent positive or negative curl amplitudes, and inverse-variance weighted these points to produce the largest outlier possible¹⁴. While there are plenty of longer sequences of consistently low or high amplitudes, only one of the 240 simulations shows a larger combined deviation. Hence we may assign an approximate revised PTE of 1/240 ≃ 0.4% to the feature. This PTE roughly matches the analytic prediction from Gaussian statistics for our band powers (0.6%). This is still very suggestive of the presence of an uncontrolled component affecting our curl band powers; the following section presents additional curl consistency checks in more detail.

4.2. Temperature lensing and lensing curl consistency and stability tests

We now turn to how the band powers differ between various reconstructions. We focus this subsection on temperature-only reconstructions, which dominate the signal and show stronger features than the MV estimator. For each one of these reconstructions, the methodology follows that described in Sect. 2, and uses the same set of 300 FFP10 simulations. All lensing biases, the point-source and Monte-Carlo corrections are recalculated consistently.

Figure 17 presents a visual summary of our consistency tests based on gradient reconstructions, and Fig. 18 for tests based on the curl. We now describe these in detail. For each pair of reconstructions, as listed below, we use the band powers to fit an amplitude relative to the fiducial model spectrum (or the flat spectrum of Eq. (58) for the curl), and build the difference of the amplitudes between the two reconstruction, both on data and FFP10 simulations. In the simulations, this difference is very well described by a Gaussian centred on zero. The matrices of differences shown in Figs. 17 and 18 are colour-coded by the value of the difference in the data in units of the standard deviation of the difference from simulations. The lower triangle shows the results over our conservative range 8 ≤ L ≤ 400, and the upper triangle the high-L range 401 ≤ L ≤ 2048. Reading the lower triangle from the left, the reconstruction labelled on the left-hand side has an inferred lensing amplitude larger than the one labelled on the top if it has a red colour, and a smaller amplitude if it has a blue colour. Likewise, reading the upper triangle from the top, the reconstruction labelled on the top is larger for a red colour. If the amplitude shift is larger than twice the expected standard deviation, the precise value is indicated, normally in black. A small number of cells show the deviation written in white. This indicates that even though the deviation is anomalously large, it is still small in absolute terms (less than one fifth of the standard deviation of the left or upper reconstruction amplitude, whichever is larger), hence not relevant for practical purposes. The black solid lines separate for clarity the following classes of test (boldface refers to the labels used to identify the various reconstructions in Fig. 17).

Component-separation methods. In addition to using our baseline SMICA map, we perform reconstructions on the CMB maps from the three other component-separation methods, NILC, SEVEM, and COMMANDER, as described in Planck Collaboration IV (2020), and the version of the SMICA CMB map that deprojects the thermal SZ signal (SMICA noSZ). In each case, the lensing reconstruction is performed on the union of the temperature and polarization confidence masks recommended by Planck Collaboration IV (2020), together with the same f_sky = 70% Galactic mask, CO mask, point-source mask, and SZ-cluster-targeted mask we use for our baseline reconstruction.

Sky cuts. We test a series of sky cuts, with corresponding masks shown in Fig. 19. Conservative Galactic dust masks with f_sky = 57% (du57) and f_sky = 44% (du44) test for Galactic foreground contamination. A reconstruction using a more aggressive galactic mask with f_sky = 80% is also shown (f80). Masking to keep only the ecliptic poles (ecl.pol.) or the ecliptic equator (ecl.eq.) tests effects related to the scanning. Masking, in Galactic coordinates, to keep only the south-east and northwest quadrants (SENW), or the southwest and northeast quadrants (SWNE), is relevant for effects related to HEALPix (Gorski et al. 2005) pixelization. We have used further two additional masks (not shown in Fig. 19), testing for contamination from SZ clusters. The SZ-unm. reconstruction uses our baseline mask but without masking SZ clusters, and SZ-cons. includes a more conservative SZ-cluster mask. This conservative cluster mask was constructed by extending our SZ-cluster mask to at least three times the virial radius (θ₅₀₀) of each cluster, and enlarging the cluster catalogue: instead of using only the internally-detected Planck SZ catalogue, we added X-ray clusters listed in the Meta-Catalog of X-ray detected Clusters of galaxies (MCXC; Piffaretti et al. 2011). This contains additional objects that are not individually detected within Planck data but whose SZ decrement is clearly visible after stacking.

CMB multipole cuts and other analysis choices. Most of the signal in the lensing reconstruction comes from CMB multipoles centred around ℓ ≃ 1450, with the Gaussian reconstruction noise N⁽⁰⁾ over the range 8 ≤ L ≤ 400 taking roughly equal contributions from either side of ℓ = 1450. The ℓ_max = 1450 reconstruction only includes CMB multipoles 100 ≤ ℓ ≤ 1450, and ℓ_min = 1451 only uses 1451 ≤ ℓ ≤ 2048. Further, ℓ_max = 1900 excludes the contribution from the smallest scales by considering 100 ≤ ℓ ≤ 1900. We also test consistency between two variations of our pipeline. In our baseline analysis, power spectra are evaluated from harmonic coefficients (either of the Wiener-filtered CMB maps, or the lensing potentials) using the simple estimator in Eq. (8). With the PCL variation, we first remap the harmonic coefficients of the convergence field estimate $\widehat{\kappa }$ to real space, compute pseudo-C_ℓ power spectra after further masking with a 20′ apodized mask, and deconvolve the effects of the apodized mask to obtain our final power spectrum estimates. Finally, with the W2500 reconstruction we test details of our filtering code, by reconstructing the Wiener-filtered map over a larger range of scales (2 ≤ ℓ ≤ 2500), and cutting the map to our baseline range (100 ≤ ℓ ≤ 2048) afterwards. This tests for leakage very close to the boundary (ℓ = 2048 in our baseline analysis).

Bias-hardened estimators. This block tests a series of temperature-only lensing estimators orthogonalized to potential contaminants, called “bias-hardened” estimators, following Namikawa et al. (2013). Sources of statistical anisotropy in the data that are not properly accounted for in the simulations can contaminate the lensing signal. Any such source, s, with a simple parametrization of its impact on the CMB covariance can be assigned an optimal quadratic estimator ${\widehat{g}}^s$. We can then define a new lensing estimate (before mean-field subtraction and normalization)

$$\begin{array}{c}\hfill {\widehat{g}}_{\mathit{LM}}^{\varphi -s}\equiv {\widehat{g}}_{\mathit{LM}}^{\varphi }-\frac{{\mathcal{R}}_L^{\varphi s}}{{\mathcal{R}}_L^{\mathit{ss}}}{\widehat{g}}_{\mathit{LM}}^s,\end{array}$$(59)

which has vanishing response to s, at the cost of an increase in variance. Here, the response functions ${\mathcal{R}}_L^{\varphi s}$ and ${\mathcal{R}}_L^{ss}$ describe the responses to s of the lensing estimator ${\widehat{g}}^{\varphi }$ and source estimator ${\widehat{g}}^s$, respectively. We give results bias hardened against point-source (S²) contamination (Osborne et al. 2014; S-hard.), where a spatially-varying point source signal $S^2(\widehat{\mathit{n}})$ is sought in the CMB data covariance, $\langle T^{\mathrm{dat}}({\widehat{\mathit{n}}}_i)T^{\mathrm{dat}}({\widehat{\mathit{n}}}_j)\rangle \ni \int {\mathrm{d}}^2\widehat{\mathit{n}}{S}^2(\widehat{\mathit{n}})\mathcal{B}({\widehat{\mathit{n}}}_i,\widehat{\mathit{n}}){\mathcal{B}}^{\ast }({\widehat{\mathit{n}}}_j,\widehat{\mathit{n}})$. The corresponding unnormalized optimal estimator is

$$\begin{array}{c}\hfill {\widehat{g}}^{S^2}(\widehat{\mathit{n}})\equiv {\overline{T}}^2(\widehat{\mathit{n}}).\end{array}$$(60)

We also harden against anisotropy in the variance of the instrumental noise (N-hard.). In this case the anisotropy is parametrized as $\langle T^{\mathrm{dat}}({\widehat{\mathit{n}}}_i)T^{\mathrm{dat}}({\widehat{\mathit{n}}}_j)\rangle \ni {\sigma }^2({\widehat{\mathit{n}}}_i){\delta }_{{\widehat{\mathit{n}}}_i{\widehat{\mathit{n}}}_j}$, with quadratic esimator

$$\begin{array}{c}\hfill {\widehat{g}}^{{\sigma }^2}(\widehat{\mathit{n}})\equiv {{\displaystyle \left(\sum _{\ell m}\frac{{\overline{T}}_{\ell m}}{b_{\ell }}{}_0{Y}_{\ell m}(\widehat{\mathit{n}})\right)}}^2,\end{array}$$(61)

where b_ℓ is an effective, isotropic approximation to the beam and pixel transfer function. We further tested hardening against spatial modulation of the CMB signal by a field $m(\widehat{\mathit{n}})$ (for example as a check against an effective spatial variation in calibration), which we denote by M-hard. The anisotropy is of the form $T^{\mathrm{dat}}({\widehat{\mathit{n}}}_i)\ni \int {\mathrm{d}}^2\widehat{\mathit{n}}\mathcal{B}({\widehat{\mathit{n}}}_i,\widehat{\mathit{n}})T(\widehat{\mathit{n}})[1+m(\widehat{\mathit{n}})]$, with resulting estimator

$$\begin{array}{c}\hfill {\widehat{g}}^m(\widehat{\mathit{n}})\equiv \overline{T}(\widehat{\mathit{n}})T^{\mathrm{WF}}(\widehat{\mathit{n}}).\end{array}$$(62)

We have additionally attempted hardening against a residual dust amplitude (D-hard.). The motivation in this case is to look for a dust component consisting of a statistically-isotropic field $d(\widehat{\mathit{n}})$ with spectral shape ${\mathcal{D}}_{\ell }^{\mathit{dd}}\equiv \ell (\ell +1)C_{\ell }^{\mathit{dd}}/(2\pi )$ modulated by a varying amplitude $A_d(\widehat{\mathit{n}})$:

$$\begin{array}{c}\hfill T^{\mathrm{dust}}(\widehat{\mathit{n}})\equiv [1+A_d(\widehat{\mathit{n}})]d(\widehat{\mathit{n}});{\mathcal{D}}_{\ell }^d\equiv \frac{{(100/\ell )}^{0.4}}{{[1+{(\ell /160)}^2]}^{0.085}}·\end{array}$$(63)

The power spectrum shape ${\mathcal{D}}_{\mathcal{l}}^{dd}$ is similar to that used to model unmasked dust power in the Planck CMB likelihood (Planck Collaboration V 2020) or in Mak et al. (2017), based on measured spectral differences between various Galactic masks. The quadratic estimator for the squared dust amplitude has the form

$$\begin{array}{c}\hfill {\widehat{g}}^{A_d}(\widehat{\mathit{n}})\equiv \overline{T}(\widehat{\mathit{n}}){\displaystyle \left(\sum _{\ell m}{C}_{\ell }^{\mathit{dd}}{\overline{T}}_{\ell m}{}_0{Y}_{\ell m}(\widehat{\mathit{n}})\right).}\end{array}$$(64)

We also applied this hardened estimator to the 217-GHz channel (217 GHz D-hard.). All these contaminant estimators are orthogonal to the lensing curl mode and do not appear in Fig. 18.

Frequency maps. We test temperature reconstructions using the HFI 143- and 217-GHz frequency channels. When using frequency maps, we perform a simple foreground cleaning of the maps by projecting out dust and CIB templates in the filtering step. This ensures that these two modes of the data are not used in the reconstruction, but, of course, the extent to which they remove dust and CIB depends on the quality of the templates. We take the 857-GHz channel as a dust template, and the GNILC CIB map at 353 GHz for the CIB. Figure 17 lists several reconstruction variations, one using our baseline f_sky = 67% mask (143 GHz, 217 GHz), and another using an f_sky = 44% mask specifically built to minimize dust contamination (143 GHz du44, 217 GHz du44). We also list a reconstruction using the 143-GHz channel with more aggressive SZ-cluster masking (143 GHz SZ-cons.) or no SZ-cluster masking at all (143 GHz SZ-unm.).

The most striking inconsistencies appearing in Fig. 17 are from the frequency reconstructions. These reconstructions use a simple template-cleaning method that is cruder than that used for the main component-separated products. On the conservative multipole range only the 217 GHz reconstruction has possible anomalies, but on the high-L range, both 143- and 217-GHz reconstructions display a striking series of inconsistencies. This demonstrates the importance of precise temperature cleaning for lensing reconstruction. Specifically, this figure strongly suggests that the 217 GHz reconstruction contains substantial Galactic foreground contamination: the dust-masked 217 GHz du44 shows very large, anomalous shifts compared to 217 GHz, on both multipole ranges displayed. While the 217 GHz row and column are predominantly red, so this reconstruction has larger lensing amplitude than most of the others, the situation is reversed for 217 GHz du44, which has slightly lower amplitudes. This is also visible on the lower panel, where 217 GHz shows a much larger lensing amplitude, which is substantially reduced on the f_sky = 44% dust mask.

All five component-separation methods show good consistency over the conservative range. However, we do see some anomalous behaviour at high-L. Notably, we find

$$\begin{array}{c}\hfill {}_{\mathrm{fid}}{\widehat{A}}_{401\to 2048}^{\varphi ,\mathit{TT},\mathrm{NILC}}-{}_{\mathrm{fid}}{\widehat{A}}_{401\to 2048}^{\varphi ,\mathit{TT},\mathrm{SEVEM}}=0.088\pm 0.029,\end{array}$$(65)

formally a 3σ deviation from zero. The actual shift in the high-L amplitude between these two reconstructions is not that large however: about 0.5σ of our high-L baseline measurement. None of the component-separation methods show large shifts in an absolute sense. Our baseline band powers using SMICA appear to be stable with respect to choice of mask, giving shifts consistent with expectations for the change in sky area. The reconstruction using the aggressive galactic mask f80 does not display suspicious deviations on either lensing multipole range, nor does it exacerbate the band-power features even slightly. SMICA noSZ relies more on the 217-GHz channel, and gives a slightly larger amplitude. While in itself not suspicious, du57 shows a significant decrease in amplitude from SMICA noSZ, and this could indicate a higher level of Galactic dust contamination in the SMICA noSZ map. We find generally no evidence for thermal SZ (tSZ) contamination: we obtain almost identical band powers using our baseline tSZ-cluster mask, the conservative tSZ-cluster mask, or no cluster mask at all, both from the 143-GHz channel or on the SMICA maps. At 143-GHz, we do however detect and subtract a much larger point-source correction without the SZ-cluster mask. This is consistent with the predictions of further tSZ-dedicated tests provided in Sect. 4.5.

The 143-GHz channel shows some consistency issues on the high-L range significant at around the 2.5σ level, where it displays consistently lower amplitudes than all other reconstructions. This is true independently of sky-cuts and bias-hardening tests. The cause of these issues is currently not understood and this further motivates us to use the conservative multipole range for our baseline cosmology results.

Figure 18 shows the same tests for the lensing curl spectrum estimated from temperature-only reconstructions. The curl estimator has less sensitivity to many sources of systematic effects; apart from the lensing gradient mode, it is also orthogonal (i.e., has zero linear response) to residual point-source contamination, CMB signal modulation, or to the dust amplitude estimator discussed above. Hence, the curl reconstruction is a stringent test of the noise properties of the simulations versus those of the data maps. Figure 18 shows that our lensing curl power is very stable across reconstructions, and that we find the same extended region of negative curl in all tests. The curl spectrum on the conservative multipole range shows sensitivity to small-scale modes of the maps used for the reconstruction (around ℓ ≃ 2000). The shifts between W2500 and ℓ_max = 1900 compared to our baseline are formally significant. As shown in the lower panel, the absolute value of the W2500 shift is very small, so this probably indicates a slight mischaracterization of the covariance rather than a systematic effect that is important for the cosmological interpretation of the baseline analysis. ℓ_max = 1900 shifts by a standard deviation on the conservative range, giving mild evidence of a tension. However, a similar problem is not seen in our lensing reconstruction. On the high-L range, the only potentially problematic test is ℓ_min = 1451, which gives a lower curl. This reconstruction is noisy and only weakly constraining on the high-L range, and shows no suspicious offset on the conservative L-range. The largest shifts overall are observed on the ecliptic poles compared to the ecliptic equators, and to a similar extent on SWNE versus SENW, for both ranges. These two cuts define regions with substantial overlap, with the curl power being higher (i.e., more consistent with zero) on the ecliptic poles, where the noise is considerably lower. However, the two reconstructions being compared share no common sky area, so differences are expected, and these shifts are not anomalous compared to simulations. The curl feature is negative, which makes it hard to explain by small inaccuracies in the simulations: after RD-N⁽⁰⁾ subtraction, which already corrects for errors in the power in the simulations to linear order, any error in the reconstructed lensing power from inaccuracies in the simulations should be quadratic in the difference and positive, so a negative curl power cannot easily be explained by a mismatch between the data and the simulations alone.

There is also a shift between SWNE and SENW (and ecliptic poles and equator) in the ϕ reconstruction, with the regions around the ecliptic poles (SENW) having higher lensing amplitude. However the variance from splitting the sky in two is large, and as before these shifts are not clearly significant compared to simulations. It is interesting to note that the SWNE region gives a larger A_L than SENW when fitting for the amplitude of the lensing smoothing effect in the temperature power spectrum (see PCP18 for discussion of A_L more generally). This is the opposite to what one would expect if the difference were due to lensing, but could be consistent with a statistical fluctuations.

4.3. Individual estimator crosses and temperature-polarization consistency

The baseline MV reconstruction optimally combines all the information in the temperature and polarization maps (in the approximation of neglecting T − E correlations in the filtering). We can also construct estimators separately for each pair of T, E, and B maps, giving five¹⁵ distinct lensing estimators. These can easily be constructed from the MV pipeline simply by zeroing the input maps that are not required for each estimator. Figure 20 shows a comparison of the MV power spectrum and all 15 distinct auto- and cross-spectra of the separate reconstructions. The lensing spectrum is measured with high S/N by many of the estimators. The polarization EB estimator is nearly independent of the TT estimator, but their cross-spectrum shows very good agreement with the baseline reconstruction.

Figure 21 shows the reconstructed lensing curl spectra. Curl reconstruction summary amplitude statistics are very consistent on the conservative range between all pairs of estimators; however, this is less true at higher lensing multipoles. Of particular interest is the range 264 ≤ L ≤ 901, built to emphasise the temperature curl feature seen in Fig. 16. Using the fiducial lensing response functions and N⁽¹⁾ subtraction, we find

$$\begin{array}{c}\hfill {}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },\mathit{TT}}=-0.072\pm 0.017.\end{array}$$(66)

The feature is specific to the temperature-only reconstruction; the next most constraining estimator over this multipole range is TT × EB, with

$$\begin{array}{c}\hfill {}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },TT\times EB}=0.016\pm 0.039.\end{array}$$(67)

The polarization-only estimator is much noisier but does not show anomalies on this range either:

$$\begin{array}{c}\hfill {}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },\mathrm{P}}=0.07\pm 0.17.\end{array}$$(68)

The amplitude of the curl feature is reduced by almost a factor of 2 in the minimum-variance reconstruction compared to the temperature-only reconstruction. This reduction is larger than expected from the FFP10 simulation behaviour. We find

$$\begin{array}{c}\hfill {}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },\mathrm{MV}}-{}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },\mathit{TT}}=0.0327\pm 0.0094,\end{array}$$(69)

a shift that is unexpected at the 3.5σ level (neglecting look-elsewhere effects). There is no similarly surprising shift in the lensing spectrum, however.

4.4. Noise tests

We now present tests of the noise by considering reconstructions using the two SMICA half-mission maps (HM1 and HM2) and their difference. First, we perform gradient and curl lensing reconstructions on their half-difference, nulling the CMB signal. We use only the noise part of the simulation suite for this analysis. The result is shown in Fig. 22 as the blue points labelled (HM₁−HM₂)⁴ (rescaled by a factor of 20 for visibility). The band powers are consistent with zero to within 2σ, as expected, with overall amplitude

$$\begin{array}{c}\hfill \begin{array}{c}{}_{\mathrm{fid}}{\widehat{A}}_{8\to 2048}^{\varphi ,\mathit{TT}}=0.002\pm 0.004,\hfill \\ {}_{\mathrm{fid}}{\widehat{A}}_{8\to 2048}^{\mathrm{\Omega },\mathit{TT}}=-0.00046\pm 0.00024.\hfill \end{array}\end{array}$$(70)

Second, we compare our baseline reconstruction to that obtained from half-mission maps. We consider two types of reconstruction, labelled (HM₁×HM₂)² and ${\text{HM}}_1^2\times {\text{HM}}_2^2$. In the first reconstruction, we estimate the lensing spectrum using a single lensing map, built using the first half-mission and second half-mission maps. In this way, there is no noise mean field in the estimate, and the power spectrum reconstruction noise only has one noise contraction. In the second reconstruction, we cross-correlate the lensing map built using the first half-mission data to another lensing map built using the second half-mission data. Both lensing maps now have a noise mean field, but the lensing reconstruction noise comes purely from CMB fluctuations. One caveat to these consistency tests is that there are no half-mission FFP10 CMB simulations accompanying the 2018 Planck release (only noise simulations). Instead, we use FFP9 half-mission CMB simulations built for the 2015 release, co-adding the frequency CMB simulations with the SMICA weighting of the 2018 release. The FFP9 and FFP10 CMB simulations share the same random seed, and apart from the weights the main differences lie in the details of the effective beam across the sky, which is important at low lensing multipoles for the lensing mean field. Using beam-processed CMB simulations also has a slight impact on the RD-N⁽⁰⁾ bias. Hence, these half-mission reconstructions should not be considered as reliable as our baseline ones. The gradient and curl half-mission spectra are displayed in the upper and lower plots, respectively, in Fig. 22, and show overall consistency with our baseline. The curl feature is less pronounced in (HM₁×HM₂)², but only slightly so. Using the FFP10 fiducial model to build summary amplitudes over the multipole range L = 264–901 covering the curl feature, we find

$$\begin{array}{c}\hfill \begin{array}{c}{}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },\mathit{TT}}=-0.058\pm 0.018\left[{({\mathrm{HM}}_1\times {\mathrm{HM}}_2)}^2\right],\hfill \\ {}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },\mathit{TT}}=-0.071\pm 0.018[{\mathrm{HM}}_1^2\times {\mathrm{HM}}_2^2],\hfill \end{array}\end{array}$$(71)

compared to −0.072 ± 0.017 for our baseline. The shifts are consistent within 1σ with those obtained from the simulation differences. On the other hand, the “dip” in the gradient spectrum around L ≃ 700 is more pronounced. We have

$$\begin{array}{c}\hfill \begin{array}{c}{}_{\mathrm{fid}}{\widehat{A}}_{638\to 762}^{\varphi ,\mathit{TT}}=0.11\pm 0.25\left[{({\mathrm{HM}}_1\times {\mathrm{HM}}_2)}^2\right],\hfill \\ {}_{\mathrm{fid}}{\widehat{A}}_{638\to 762}^{\varphi ,\mathit{TT}}=0.20\pm 0.25[{\mathrm{HM}}_1^2\times {\mathrm{HM}}_2^2],\hfill \end{array}\end{array}$$(72)

compared to 0.29 ± 0.24 for our baseline. These shifts are again within 1σ of the expected differences, but lower the spectrum further with respect to the ΛCDM prediction. The low value is driven by the first half-mission map as we now discuss. Taking amplitudes with respect to the Planck best-fit Planck TT,TE,EE+lowE+lensing, and marginalizing over CMB uncertainties in order to quantify better the discrepancy, we have

$$\begin{array}{c}\hfill \begin{array}{c}{}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{638\to 762}^{\varphi ,TT}=-0.51\pm 0.35[{({\mathrm{HM}}_1)}^4;\text{CMB marginalized}]\hfill \\ {}_{\mathrm{b}.\mathrm{f}.}{\widehat{A}}_{638\to 762}^{\varphi ,TT}=0.43\pm 0.32[{({\mathrm{HM}}_2)}^4;\text{CMB marginalized}].\hfill \end{array}\end{array}$$(73)

The (HM₁)⁴ result is formally a 4.5σ deviation from the best-fit lensing spectrum value. The remaining band powers of the HM₁ reconstruction show no inconsistencies with the best-fit, however, and the observed shift between the two amplitudes in Eq. (73) is within 2σ of the difference expected from simulations. The origin of these curious dip differences is not currently understood, but the dip is not included in our conservative multipole range, so it is not directly relevant for our baseline results.

The MV baseline lensing spectrum and the half-mission spectra are still low on these multipole ranges for both curl and gradient, but less than the temperature-only spectra. We find

$$\begin{array}{c}\hfill \begin{array}{c}{}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },\mathrm{MV}}=-0.029\pm 0.017\left[{({\mathrm{HM}}_1\times {\mathrm{HM}}_2)}^2\right],\hfill \\ {}_{\mathrm{fid}}{\widehat{A}}_{264\to 901}^{\mathrm{\Omega },\mathrm{MV}}=-0.038\pm 0.018[{\mathrm{HM}}_1^2\times {\mathrm{HM}}_2^2],\hfill \end{array}\end{array}$$(74)

compared to −0.039 ± 0.015 for our baseline, and

$$\begin{array}{c}\hfill \begin{array}{c}{}_{\mathrm{fid}}{\widehat{A}}_{638\to 762}^{\varphi ,\mathrm{MV}}=0.43\pm 0.24\left[{({\mathrm{HM}}_1\times {\mathrm{HM}}_2)}^2\right],\hfill \\ {}_{\mathrm{fid}}{\widehat{A}}_{638\to 762}^{\varphi ,\mathrm{MV}}=0.30\pm 0.27[{\mathrm{HM}}_1^2\times {\mathrm{HM}}_2^2],\hfill \end{array}\end{array}$$(75)

compared to 0.45 ± 0.23 for our baseline.

4.5. Tests of CMB lensing foreground correlations

We add residual extragalactic foreground power to our simulations as an independent Gaussian component. In reality, extragalactic foregrounds are non-Gaussian and are correlated with the lensing signal since both are affected by the same matter perturbations. There is therefore a concern that residual foregrounds could lead to additional contributions to the quadratic estimator power that are not corrected for in our pipeline. In particular SZ, CIB, and clustered point-source contributions (and/or their local power) could be directly correlated with the lensing potential. The polarization reconstruction is expected to be essentially free from these contaminants due to the much lower level of foreground power in the small-scale polarization maps. Our use of foreground-cleaned SMICA maps is also expected to remove the bulk of the CIB foreground signal in temperature; however, SZ and unresolved extragalactic radio sources are largely unaffected, and foreground residuals remain due to instrumental noise. Semi-analytic estimates suggest that the possible bias on the temperature reconstruction power should not be large at the sensitivity level of Planck (Osborne et al. 2014; van Engelen et al. 2014; Ferraro & Hill 2018), but it is important to test this more directly.

To assess the possible impact of correlated foregrounds we have used a small number of non-Gaussian foreground simulations that are constructed to have approximately the correct correlation structure. The Planck Sky Model (PSM; Delabrouille et al. 2013; Planck Collaboration XII 2016) software is extended to consistently generate maps of the CMB lensing potential, CIB, and SZ components, including their correlations. To do this, the Boltzmann code class (Blas et al. 2011) is used to calculate the correlated angular power spectra of the unlensed CMB, lensing potential, and matter distribution (neglecting matter-CMB correlations) over 64 concentric shells between redshifts z = 0.01 and z = 6 up to a maximum multipole l_max = 4096. The power spectra for the FFP10 fiducial model are calculated using the nonlinear halofit model from Takahashi et al. (2012) using modes up to k_max = 1 h Mpc⁻¹. The lensing potential is assumed to be Gaussian, but the matter density fields ${\rho }_i(\widehat{\mathit{n}})$ in each shell are simulated as log-normal fields; the covariance from class is therefore mapped into the covariance of the log fields, so that a Gaussian realization of $s_i(\widehat{\mathit{n}})$ can be exponentiated to obtain log-normal matter density fields with the correct covariance structure (Greiner & Enßlin 2015). The CMB component is lensed with the lensing potential using LensPix (Lewis 2005). We have checked that the impact of lensing of the CIB is negligible, as expected (Schaan et al. 2018).

The CIB galaxies are grouped into three different populations according to their spectral energy distributions (SEDs): proto-spheroid, spiral, and starburst. The flux density in each population is randomly distributed according to the redshift-dependent number counts from Planck ERCSC (Negrello et al. 2013), JCMT/SCUBA-2 (Chen et al. 2013), AzTEC/ASTE (Scott et al. 2012), and Herschel/SPIRE (Béthermin et al. 2013) observations. Each redshift shell is populated with CIB galaxies and galaxy clusters with probabilities proportional to the density contrast distribution, with a population-dependent bias. The CIB maps constructed in this way have power spectra that agree with measurements from Planck data at high frequencies using several different methods for dealing with Galactic dust and CMB contamination (Planck Collaboration XXX 2014; Mak et al. 2017; Lenz et al. 2019).

Galaxy clusters are simulated by drawing a catalogue of halos from a Poisson distribution of the Tinker mass function (Tinker et al. 2008). The thermal SZ emission of each halo at a given redshift and mass is modeled according to the temperature (Arnaud et al. 2005) and the pressure profiles given by Arnaud et al. (2011). The simulated SZ clusters are then distributed over the redshift shells with probabilities proportional to matter density fields to have consistent correlations between SZ, CMB, and CIB. The resulting cluster counts are consistent with those observed by Planck. The kinetic SZ component is simulated by using the electron density of the clusters and a Gaussian realization of the 3D cluster velocity field derived from the density fluctuations (following Peebles 1993).

We first use the simulations to assess the expected impact of the CIB component on our lensing reconstruction. We test the 217-GHz channel component (without any cleaning), and a cleaned map constructed by combining the different frequency simulations using the SMICA frequency weights. We use the same multipole cuts as our baseline analysis, but include the full sky for simplicity. We see a CIB signature in the lensing spectrum consistent with previous work (van Engelen et al. 2014), with two principal biases: a positive bias from the CIB trispectrum (labelled CIB⁴ in the following); and a larger contribution from the bispectrum (CIB² × κ). We isolate the first contribution by performing direct lensing reconstruction on the simulated CIB map, then estimating the power spectrum subtracting (disconnected) biases with Gaussian isotropic realizations from the CIB power spectrum. To estimate the bispectrum contribution, we cross-correlate the lensing estimator $\widehat{\varphi }$ as applied on the CIB map to the input lensing potential. The resulting bias is twice this cross-spectrum. Figure 23 shows the biases that are expected. For the uncleaned 217-GHz channel (blue) there is at most about a 2% bias over the conservative multipole range. This dominant (bispectrum) contribution is effectively reduced by one order of magnitude after SMICA cleaning (orange), making it negligible. It is larger, but still sub-percent level, in the case of the tSZ-deprojected SMICA weighting (green). The trispectrum biases are even smaller (red and purple). A similar analysis for the lensing curl shows that all these terms are expected to be completely negligible for the curl power.

We now turn to the thermal SZ component, focussing on the 143-GHz channel since it is small in the 217-GHz channel where the detector bandpasses are centred on the null of the tSZ spectrum. We perform several reconstructions on the same CMB simulation, with and without thermal SZ, and with and without masking to emulate the effect of the cluster mask that we apply in our baseline analysis. These tests allow us both to demonstrate the good performance of our point-source correction procedure, even in the absence of masking, as well as the expected robustness of our band-powers. To build the simulation cluster mask, we mask the same number of objects as our baseline mask (rescaled by the respective observed area), starting from the objects with the strongest integrated Compton Y₅₀₀ emission, and using the same criteria to define the masking radii for each object. This procedure results in a total of 1148 objects masked (out of 377 563) for a resulting masked sky fraction of roughly 0.3%. Some effectively pure point-source signal remains in the map, which would be detected by the Planck source-detection methodology and also masked in our actual data lensing mask. However, these residual signals have no impact on these tests and we perform no further cleaning.

Figure 24 collects the different reconstructions, showing for clarity the conservative multipole range on the left panel, and the remaining high-L multipoles on the right panel. We show deviations of the reconstruction from the fiducial input lensing power spectrum, with the blue points showing the reconstruction for the reference simulation without the SZ mask. This full-sky reference simulation includes instrumental noise in the form of Gaussian noise, independent between pixels but with amplitude in each pixel following the 143-GHz channel variance map. Owing to the larger sky area, the blue error bars are slightly tighter than those on our 143-GHz channel data reconstruction. The tSZ signal appears largely as a point-source-like contamination, and we distinguish results before (orange) and after (green) point-source correction, calculated and subtracted from the band powers in the same way as in the main data analysis. The correction is very effective at correcting the power over the conservative range, but residual additional power is clearly visible at high L. Finally, the red points show the result after application of the cluster mask. At this point, the contamination (i.e., the difference from the reconstructed power for the reference simulation) becomes a negligible fraction of a standard deviation. We also perform the same tests for the curl reconstruction band powers; the curl estimator does not respond to point-source signals, and we find a negligible impact, even without any masking.

Applying a cluster mask systematically removes small regions of the sky with strong lensing convergence power, and could bias the lensing spectrum estimate. We assess the size of this effect on simulations by comparing lensing spectra before and after masking. The cluster masks are produced for each simulation in the way described above, and we obtain spectra after deconvolving from pseudo-$C_L^{\kappa \kappa }$ power spectra the effects of the (slightly enlarged and then apodized) mask. We find this bias to be at most 0.25% in amplitude in all bins (limited by sampling noise), and we do not consider it further.

Finally, we perform test reconstructions using the entire set of foregrounds available (CIB, thermal and kinetic SZ, and point sources) using the SMICA frequency weighting, and incorporating the cluster mask described above. The foregrounds contribute an additional 5% power around ℓ = 2000 (from the CIB and point sources), in good agreement with the mismatch observed in the data compared to the (clean) FFP10 simulation suite. We perform a full lensing reconstruction on this map, adding Gaussian isotropic power to the Monte Carlo simulations to account for the difference in power, in the same way as for our analysis of the actual data, and compare to a reference reconstruction that includes the CMB component only. We find shifts in all bins that are typically 0.1(0.13)σ of our baseline reconstructions in the conservative (aggressive) multipole range. All deviations seem consistent with the scatter expected from the addition of the extra power at high CMB multipoles.

4.6. Test of dependence on fiducial model

The lensing pipeline assumes a fiducial model for the true CMB and lensing power spectra, with the likelihood perturbatively correcting for differences between the fiducial model and the actual model spectra at each point in parameter space. For our results to be robust, they should be almost independent of which fiducial model was assumed. To test this, we carry out a simplified analysis using a different fiducial model: a non-flat ΛCDM model best-fit to the CMB power spectra with Ω_K ≃ −0.04. The test model has 15–20% more lensing power than the FFP10 fiducial model, and also has CMB power spectra that differ by approximately 2% at high ℓ. We note that we do not need to test large deviations of the CMB spectra because they are now empirically measured to good accuracy: the main variations still allowed in the spectral shape are due to residual cosmic variance and foreground uncertainties, both of which are a small fraction of the total CMB signal on scales relevant for lensing reconstruction.

Running a full set of FFP10 simulations with a different fiducial model would be numerically very expensive, so we instead generate a set of simpler idealized isotropic-beam simulations with the FFP10 fiducial model and the test Ω_K ≠ 0 model, and check for differences in the lensing reconstruction based on these two sets of simulations and theoretical spectra. The spectra enter into the filtering of the data, the analytic estimate of the estimator response, and via the simulations used for evaluation of the mean field, RD-N⁽⁰⁾ and MC corrections, and via the likelihood in the perturbative correction functions and fiducial value of N⁽¹⁾. The simpler set of simulations we use here lacks the detailed beam shape and scanning model that affects the signal at low-L. For this reason we compare the reconstructions not on the actual Planck data, but on one of the simulations, generated with the FFP10 CMB spectra as input. We show in Sect. 4.8 below that the N⁽¹⁾ contribution on the cut sky is accurately modeled by full-sky analytic results, so we use the full sky for convenience in this test.

We test robustness to the choice of fiducial model by analysing a single FFP10 (Ω_K = 0) simulation with the consistent fiducial simulation set, and then by analysing it using the inconsistent (Ω_K ≠ 0) fiducial simulation set. With the consistent fiducial model, using the full lensing multipole range to fit for a lensing amplitude we find

$$\begin{array}{c}\hfill {}_{\mathrm{fid}}{\widehat{A}}_{8\to 2048}^{\widehat{\varphi },\mathit{TT}}=1.032\pm 0.025\text{(FFP10 fiducial sim.).}\end{array}$$(76)

Over this same multipole range, fitting the same reconstructed power with a lensing amplitude relative to the Ω_K fiducial model results in 0.88 ± 0.021, consistent with the increased lensing power in this model and excluding the model at more than 5σ. If we now perform lensing reconstruction on the same data, but with the Ω_K model replacing the fiducial FFP10 model throughout the analysis, and then fit a lensing amplitude relative to the FFP10 model in the likelihood, we get

$$\begin{array}{c}\hfill \begin{array}{c}{}_{\mathrm{fid}}{A}_{8\to 2048}^{\widehat{\varphi },\mathit{TT}}=1.009\pm 0.026\text{(no lin. corrections)}\hfill \\ {}_{\mathrm{fid}}{A}_{8\to 2048}^{\widehat{\varphi },\mathit{TT}}=1.026\pm 0.026\text{(with lin. corrections)}\hfill \\ ({\mathrm{\Omega }}_K\mathrm{fiducial}\mathrm{simulations}).\hfill \end{array}\end{array}$$(77)

While the estimate without linear corrections differs by one standard deviation from that in Eq. (76), the linear corrections are effective at reducing the discrepancy to a small fraction of a standard devitation, demonstrating consistency of our estimator and likelihood methodology. Note that we do not expect perfect agreement, since the estimator’s filtering and quadratic estimator weights do differ slightly. For the most part, the linear correction in the likelihood occurs at high lensing multipoles for the N⁽¹⁾ bias subtraction, where the lensing power spectrum is up to 20% higher in the Ω_K ≠ 0 model than the flat model (see Fig. 25).

4.7. Lensing Gaussianity assumption

The quadratic estimator formalism and FFP10 simulations assume that the lensing potential is Gaussian, but it is expected to be non-Gaussian at some level due to nonlinear growth of large-scale structure (LSS) and post-Born lensing (Pratten & Lewis 2016). A non-vanishing LSS bispectrum will source contributions to the lensing reconstruction spectrum involving three powers of ϕ, potentially giving rise to an additional $N_L^{(3/2)}$ lensing bias as well as the N⁽⁰⁾ and N⁽¹⁾ biases that we already model (Böhm et al. 2016). A subset of the contractions leading to this bias were first studied by Böhm et al. (2016) assuming tree-level perturbation theory for the LSS bispectrum and neglecting post-Born contributions, finding that the terms are negligible at Planck noise levels. Beck et al. (2018) and Böhm et al. (2018) have shown on simulations that including post-Born effects reduces the signal even further, as expected from the opposite sign of post-Born and LSS contributions to most configurations of the lensing bispectrum (Pratten & Lewis 2016). We therefore consider the bias on the lensing reconstruction power spectrum to be negligible and neglect it. The bias on cross-correlations between the lensing spectrum and low-redshift large-scale structure tracers could, however, be larger due to larger low-redshift LSS non-Gaussianity, and smaller opposite-sign post-Born contributions at low redshift. However, for Planck noise levels, following Fabbian et al. (2019) we estimate that the cross-correlation bias remains below 1% for tracers at z ≳ 0.2 where there is a useful cross-correlation signal, and hence should also be negligible compared to errors.

Note that although the lensing signal is expected to be close to Gaussian, the lensing reconstruction noise is expected to be non-Gaussian, since it is a nonlinear function of the maps, for example the reconstruction 1-point function is skewed (Liu et al. 2016). However, this non-Gaussianity does not affect our analysis, and the power spectrum band powers are well approximated as Gaussian to the required level of accuracy.

4.8. Tests of the N⁽¹⁾ lensing bias

The N⁽¹⁾ bias originates from secondary contractions of the lensed CMB trispectrum (Kesden et al. 2003; Hanson et al. 2011), and affects both gradient and curl deflection reconstructions. We follow PL2015 in that our baseline band powers simply subtract the N⁽¹⁾ bias evaluated in a fiducial model, with the likelihood correcting for the model dependence perturbatively. This section outlines some tests to show that this is not significantly suboptimal, nor a source of bias. We first discuss two ways to make N⁽¹⁾-subtracted band powers without assuming a fiducial lensing spectrum. We then test a simulation-based N⁽¹⁾ estimator, assessing the impact of sky cuts and sky curvature; this also forms our baseline N⁽¹⁾ calculation in the case of inhomogeneous filtering.

The N⁽¹⁾ correction captures the estimator’s linear response to C^ϕϕ from non-primary contractions. As such, we can write (neglecting point-source and MC corrections, and complications due to masking in the following discussion)

$$\begin{array}{c}\hfill \langle |{\widehat{\varphi }}_{\mathit{LM}}{|}^2\rangle =N_L^{(0)}+{\displaystyle \sum _{L\prime }({\delta }_{LL\prime }+N_{LL\prime }^{(1)})C_{L\prime }^{\varphi \varphi }.}\end{array}$$(78)

Instead of subtracting the right-most term using a fiducial $C_L^{\varphi \varphi }$, we can invert this matrix relation to obtain the lensing spectrum. This can also be generalized to calculate the curl N⁽¹⁾, which is used for the curl null test. With X, Y each standing for the gradient or curl deflection, let N^(1),XY be the Y-induced N⁽¹⁾ contribution to the X-spectrum estimate, we can write

$$\begin{array}{cc}\hfill \left(\begin{array}{c}\langle |{\widehat{\varphi }}_{\mathit{LM}}{|}^2\rangle \\ \langle |{\widehat{\mathrm{\Omega }}}_{\mathit{LM}}{|}^2\rangle \end{array}\right)& =\left(\begin{array}{c}{N}_L^{(0),\varphi \varphi }\\ N_L^{(0),\mathrm{\Omega }\mathrm{\Omega }}\end{array}\right)\hfill \\ \hfill & +{\displaystyle \sum _{L\prime }\left(\begin{array}{cc}{\delta }_{LL\prime }+N_{LL\prime }^{(1),\varphi \varphi }& N_{LL\prime }^{(1),\varphi \mathrm{\Omega }}\\ N_{LL\prime }^{(1),\mathrm{\Omega }\varphi }& {\delta }_{LL\prime }+N_{LL\prime }^{(1),\mathrm{\Omega }\mathrm{\Omega }}\end{array}\right)\left(\begin{array}{c}{C}_{L\prime }^{\varphi \varphi }\\ C_{L\prime }^{\mathrm{\Omega }\mathrm{\Omega }}\end{array}\right).}\hfill \end{array}$$(79)

In this equation $C_L^{\mathrm{\Omega }\mathrm{\Omega }}$ represents a hypothetical curl-like deflection caused, for example, by an instrumental systematic effect that we assume is uncorrelated with the lensing signal. Sub-pixel effects or pointing errors could create such a signal (see Appendix B).

Inverting Eq. (79) provides slightly different lensing gradient and curl band powers, labelled “xpRD-N⁽¹⁾” in Fig. 26, while inverting Eq. (78) provides yet another set of lensing band powers, labelled “RD-N⁽¹⁾.” The results for the temperature-only reconstructions are shown in Fig. 26. Lensing band powers are shown in the top panel, after taking the ratio to the FFP10 fiducial spectrum; the lower panel shows the curl band powers. In both cases, the points with error bars show our baseline reconstructions, using the fiducial N⁽¹⁾ subtraction. For comparison, the red curves display the case of no N⁽¹⁾ subtraction, and the black dashed curves the N⁽¹⁾ bias itself. In implementing the direct inversion of the linear relations above, we have truncated L′ at 3000, and the matrix elements are evaluated numerically according to the flat-sky isotropic analytic form given in Appendix A. These different methods for dealing with the N⁽¹⁾ bias have very little impact on the band powers, with a change of at most a small fraction of 1σ on the smallest scales, and are consistent with results from our fiducial model. The N⁽¹⁾ subtraction method does not substantially affect our curl null-test results either. A deconvolution of the N⁽¹⁾ bias might also be desirable in order to reduce off diagonal contributions to the band powers’ covariance (Peloton et al. 2017), but we found the gain to be negligible.

We also test a simulation-based N⁽¹⁾ estimation, labelled “MC-N⁽¹⁾” in Fig. 26. We use pairs of simulations that share the same lensing deflections but have different CMB realizations, to build an estimate of the bias. The precise way to do this follows Story et al. (2015), and is also reproduced for completeness in Appendix A. In contrast to the other methods described above, this way of estimating the bias takes into account our exact sky cuts, and also includes sky curvature (which is neglected in our analytic N⁽¹⁾ calculations). We use simplified, effectively isotropic CMB simulations that do not contain the full anisotropic Planck beam model or SMICA processing, since a full re-simulation would be too numerically expensive. The agreement between the simulation-based N⁽¹⁾ estimator and the other methods is good, to the point that it does not visibly affect the band powers.

We use the simulation-based MC-N⁽¹⁾ as the baseline for our results using inhomogeneous filtering. Figure 27 shows MC-N⁽¹⁾ (blue points) for the polarization reconstruction using the inhomogeneously-filtered maps. The large dynamic range of the noise-variance map used for filtering makes the naive analytic prediction (the orange curve) fail to reproduce the MC-N⁽¹⁾ shape. However, it is straightforward to understand the simulation result analytically by using the same toy model discussed in Sect. 2.3. Let $N_L^{(1)}(\widehat{\mathit{n}})$ be the analytic (full-sky) N⁽¹⁾ bias calculated using the noise levels at angular position $\widehat{\mathit{n}}$, as specified by the variance map used in the filtering. The green curve in Fig. 27 shows the weighted average

$$\begin{array}{c}\hfill N_L^{(1)}\simeq {\displaystyle \int \frac{\mathrm{d}\widehat{\mathit{n}}}{4\pi }{\left(\frac{{\mathcal{R}}_L(\widehat{\mathit{n}})}{{\mathcal{R}}_L^{\mathrm{fid}}}\right)}^2{N}_L^{(1)}(\widehat{\mathit{n}}).}\end{array}$$(80)

This is evaluated by splitting the sky into 15 patches with roughly constant noise levels, and successfully reproduces the simulation-based N⁽¹⁾.

5. Data products

The final Planck lensing data products are available on the Planck Legacy Archive¹⁶, and described in more detail in the Explanatory Supplement (Planck Collaboration 2018). The lensing analysis data release consists of the following products:

– the baseline MV convergence (κ) reconstruction map up to L_max = 4096 based on the SMICA CMB map, along with the corresponding simulations, mask, mean field, lensing biases and response functions;

– variations using only temperature information, only polarization information, and with no SZ-cluster mask;

– temperature, polarization and MV lensing maps using optimal filtering for the noise inhomogeneity, together with simulations.

– temperature lensing maps and simulations built from the SMICA tSZ-deprojected maps (up to L_max = 2048);

– joint MV+CIB and TT+TE+EE+CIB reconstruction κ maps, primarily for use with delensing and for plotting the current best-estimate lensing potential, together with the corresponding simulation suite.

– lensing B-mode templates and simulations built from our best E-mode estimate and the TT+TE+EE lensing reconstruction, and in combination with the CIB;

– power spectrum band powers, covariance, and linear correction matrices for the various analyses;

– likelihood code using the band powers for the conservative and aggressive multipole ranges from the baseline analysis; and

– cosmological parameter tables and MCMC chains.

6. Conclusions

We have presented the final official Planck lensing analysis, and described in detail the limits of our understanding of the data. The baseline lensing reconstruction, over nearly 70% of the sky and using lensing multipoles 8 ≤ L ≤ 400, is robust to a wide variety of tests. The reconstruction S/N is of order 1 at L ≲ 100, but retains significant statistical power to smaller scales. It gives interesting constraints on cosmological parameters on its own, yielding percent-level estimates of ${\sigma }_8{\mathrm{\Omega }}_{\mathrm{m}}^{0.25}$, and tight constraints on individual parameters when combined with BAO and a baryon density prior. The Planck lensing results are currently competitive with galaxy survey constraints, and in the case of σ₈ substantially more powerful due to the weaker degeneracy with Ω_m. In combination with lower-redshift tracers, degeneracies can be further broken, and CMB lensing provides a powerful high-redshift baseline for joint constraints.

We showed that a joint analysis with the Planck CIB map can be used to provide lensing estimates out to much smaller scales than with lensing reconstruction alone. The combined map provides our current best estimate of the integrated mass in the Universe between today and recombination, and can be used for delensing analyses. We demonstrated that delensing can already be achieved by Planck, giving peak-sharpening and a reduction in B-mode power in line with expectations. For future B-mode polarization observations, delensing will probably be essential.

We studied in detail the limits of our understanding of the data. In particular, the deficit of curl power on small scales persists and appears to be very robust to analysis choices. The roughly 3σ significance of this signal is at a level where it starts to be concerning, and may be correlated with sky direction in a way that suggests an unknown systematic issue. However, the changes could also be purely statistical, and we cannot at this point give any likely origin of a systematic signal that could explain it. The reconstructed power on smaller scales, L >  400, also shows less stability to changes in foreground modeling and map choices, so we restrict attention to the conservative multipole range 8 ≤ L ≤ 400 for our main cosmology results. Our likelihoods are made publicly available for both the conservative and aggressive (8 ≤ L ≤ 2048) ranges, and users may opt to use the full multipole range at their discretion. For cross-correlation studies, where the detailed modeling of the auto-spectrum bias is not required, the full multipole range may be more reliable, but we have not performed detailed consistency checks for that application.

In addition to the baseline analysis, we have also provided a number of substantial analysis improvements. In particular, we gave the first demonstration that anisotropic filtering of the polarization can significantly improve the S/N in the reconstruction, giving our best polarization-only band power lensing estimates.

Our baseline lensing reconstruction maps and simulations are made publicly available, along with the joint CIB map and variations with and without SZ-cluster masking. No planned experiments will be able to provide comparable quality maps over the full sky for many years. Ongoing and forthcoming ground-based observations are expected to improve greatly the reconstruction over patches of the sky, but the Planck reconstructions are likely to remain the only tracers of the very largest-scale lensing modes for the immediate future.

---

Acknowledgments

We thank Duncan Hanson for all his work on the previous Planck analyses without which the current work would not have been possible. We thank the DES collaboration for sharing their data. Some of the results in this paper have been derived using the HEALPix package. Support is acknowledged from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement No. [616170], and from the Science and Technology Facilities Council [grant numbers ST/L000652/1 and ST/N000927/1, respectively]. The Planck Collaboration acknowledges the support of: ESA; CNES, and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at https://www.cosmos.esa.int/web/planck/planck-collaboration.

References

Appendix A: Power spectrum biases

This appendix describes in more detail the calculation of the lensing power spectrum biases. For brevity, here we only present results for our baseline analysis, where all four input maps entering the power spectrum estimator come from the same data and simulation sets. PL2015 contains a thorough description of the most general case with mixtures of input maps. Figure A.1 shows the size of the various biases for our baseline analysis.

We introduce three types of lensing power spectrum estimators built using the simulations: ${\widehat{C}}_L^{\mathit{ii}}$ is the power spectrum of the quadratic estimator with both legs¹⁷ being simulation i; ${\widehat{C}}_L^{\mathit{di}}$ is the analogous spectrum where one leg of the quadratic estimator is always the data, and the second leg is simulation i; and ${\widehat{C}}_L^{\mathit{ij}}$ is where simulation i is the first leg and simulation j the second leg and i ≠ j. We note that ${\widehat{C}}_L^{\mathit{ii}}$ contains a mean-field subtraction, but the other two spectra do not (the mean field would vanish, since the different maps are independent to a very good approximation). The Gaussian lensing reconstruction noise biases MC-N⁽⁰⁾ and RD-N⁽⁰⁾ are then defined as

$$\begin{array}{c}\hfill \begin{array}{c}\text{MC-}{N}_L^{(0)}\equiv {\langle 2{\widehat{C}}_L^{\mathit{ij}}\rangle }_{N_{\mathrm{bias}}},\hfill \\ \text{RD-}{N}_L^{(0)}\equiv {\langle 4{\widehat{C}}_L^{\mathit{di}}-2{\widehat{C}}_L^{\mathit{ij}}\rangle }_{N_{\mathrm{bias}}},\hfill \end{array}\end{array}$$(A.1)

where angle brackets denote an average over simulations, with index i, and in our implementation j is always equal to i + 1 (cyclically). The numerical factors in front of ${\widehat{C}}_L^{\mathit{ij}}$ (or ${\widehat{C}}_L^{\mathit{di}}$) in Eq. (A.1) account for the fact that the quadratic estimators entering these spectra are built with two independent CMB maps, hence only half of the Gaussian contractions are captured. In principle, many more pairs could be used to estimate $\langle 2{\widehat{C}}_L^{\mathit{ij}}\rangle$; this would require considerably more resources, and the Monte Carlo error on this term is already a small correction, as demonstrated in Appendix C. The Monte Carlo reconstruction ${\langle {\widehat{C}}_L^{\varphi \varphi }\rangle }_{\mathrm{MC}}$ (used to build the Monte Carlo correction, Eq. (10), after binning) is

$$\begin{array}{c}\hfill {\langle {\widehat{C}}_L^{\varphi \varphi }\rangle }_{\mathrm{MC}}\equiv {\langle {\widehat{C}}_L^{\mathit{ii}}-2{\widehat{C}}_L^{\mathit{ij}}-N_L^{(1)}\rangle }_{N_{\mathrm{bias}}}.\end{array}$$(A.2)

In the last definition, the analytic expression for N⁽¹⁾ is used, with the exception of our inhomogeneously-filtered reconstructions, for which a Monte Carlo estimate MC-N⁽¹⁾ is used instead. In general, for a source of anisotropy s, the analytic s-induced N⁽¹⁾ is calculated according to the flat-sky approximation using (Kesden et al. 2003)

$$\begin{array}{cc}\hfill N_L^{(1)XYIJ,s}& =\frac{1}{{\mathcal{R}}_L^{\mathit{XY}}}\frac{1}{{\mathcal{R}}_L^{\mathit{IJ}}}\int \frac{{\mathrm{d}}^2{\mathit{\ell }}_1{d}^2{\mathit{\ell }}_{1^{\prime }}}{{(2\pi )}^4}{F}_{l_1}^X{F}_{l_2}^Y{F}_{l_{1^{\prime }}}^I{F}_{l_{2^{\prime }}}^J\hfill \\ \hfill & \times W^{\mathit{XY}}({\mathit{\ell }}_1,{\mathit{\ell }}_2)W^{\mathit{IJ}}({\mathit{\ell }}_{1^{\prime }},{\mathit{\ell }}_{2^{\prime }})\hfill \\ \hfill & \times [C_{|{\mathit{\ell }}_1+{\mathit{\ell }}_1^{\prime }|}^{\mathit{ss}}{f}^{XI,s}({\mathit{\ell }}_1,{\mathit{\ell }}_1^{\prime })f^{YJ,s}({\mathit{\ell }}_2,{\mathit{\ell }}_2^{\prime })\hfill \\ \hfill & +C_{|{\mathit{\ell }}_1+{\mathit{\ell }}_2^{\prime }|}^{\mathit{ss}}{f}^{XJ,s}({\mathit{\ell }}_1,{\mathit{\ell }}_2^{\prime })f^{YI,s}({\mathit{\ell }}_2,{\mathit{\ell }}_1^{\prime })],\hfill \end{array}$$(A.3)

with ${\mathit{\ell }}_1+{\mathit{\ell }}_2=\mathbf{L}=-({\mathit{\ell }}_{\mathbf{1}\prime }+{\mathit{\ell }}_{\mathbf{2}\prime })$, and the labels X, Y, I, J stand for the T, E, or B map. Further, $F_l^X$ is the analytic inverse-variance filter 𝒯^†Cov⁻¹𝒯 in the isotropic limit for field X (diagonal in T, E, B space since we perform temperature and polarization filtering separately), W^XY are the flat-sky estimator weights as applied to the inverse-variance filtered X and Y fields, f^XY, s is the CMB XY flat-sky covariance response to anisotropy source s, and $C_L^{ss}$ is the spectrum of that source of anisotropy. Results in this paper use lensing as source of anisotropy s almost exclusively, but we have also built curl-induced and point-source induced N⁽¹⁾ biases for robustness tests.

The Monte Carlo MC-N⁽¹⁾ (Story et al. 2015) is built as follows. We generate additional noiseless simulations in pairs, where each pair shares the same input deflection field, but have independent unlensed CMB. These simulations are not propagated through the Planck beam processing, since this would be computationally too expensive, but are generated with the effective, isotropic transfer function instead (the accuracy requirements on N⁽¹⁾ are much lower than on N⁽⁰⁾, and Sect. 4.8 demonstrates that the impact of the CMB non-idealities on N⁽¹⁾ are negligible). These simulations also neglect $C_L^{T\varphi }$ and $C_L^{E\varphi }$, which are negligible on the relevant scales. We again form spectra mixing quadratic estimates ${\widehat{C}}_L^{\mathit{ij}}$ and ${\widehat{C}}_L^{ii\prime }$ from these simulations, where ${\widehat{C}}_L^{\mathit{ij}}$ is defined as above on the first pair member of these new simulations, and ${\widehat{C}}_L^{ii\prime }$ takes the quadratic estimate built from the first pair member on one leg and the second pair member on the second leg. Then

$$\begin{array}{c}\hfill \text{MC-}{N}_L^{(1)}\equiv \langle 2{\widehat{C}}_L^{i{i}^{\prime }}-2{\widehat{C}}_L^{\mathit{ij}}\rangle .\end{array}$$(A.4)

Appendix B: Mean fields

In this appendix we characterize the main components of the mean field, including very high lensing multipoles L >  2048 that are not used in this paper, but are included in the lensing maps that we release. The mean field is obtained by averaging lensing reconstruction maps constructed from the FFP10 simulations, and for many effects we rely on the fidelity of these simulations for accurate subtraction of the mean field from the data.

Figure B.1 shows the different contributions to the mean field of the temperature-only reconstruction. The left-hand panel focuses on the low-L range (L ≤ 400) and the right-hand panel on higher multipoles. The total mean field is shown as the blue curve. The orange curve shows the contribution to the mean field from the noise simulations. This captures various effects, including the large inhomogeneity of the pixel noise variance (see below and Fig. 2), and residual noise correlations from the mapmaking procedure. The green curve shows the contribution from statistical anisotropy of the observed CMB, through the effects of beam anisotropies and pixelization. The orange and green curves were obtained from full-sky reconstructions to remove the effect of masking, and used idealized statistically-isotropic CMB and noise components, respectively (which do not contribute to the mean field). The mask contribution itself is shown as the red curve, which we obtain by differencing the mean-field spectra obtained from the full, anisotropic, and masked simulations to the mean field from the same simulations without any masking.

As in previous releases, the mean field is very large at low multipoles, where it is mainly due to masking, with the noise and beam-anisotropy contribution being the next most important. The pink curve in Fig. B.1 shows the expected contribution from the inhomogeneities in the pixel noise variance only (using the variance map in Fig. 2 and assuming there are no pixel-pixel correlations). The expected noise mean-field contribution is, from Eq. (3),

$$\begin{array}{c}\hfill \langle {}_1\widehat{d}(\widehat{\mathit{n}})\rangle \text{(noise inhom. only)}={\displaystyle \sum _{\mathit{LM}}{h}_L{\sigma }_{T,LM}^2{}_1{Y}_{\mathit{LM}}(\widehat{\mathit{n}}),}\end{array}$$(B.1)

and hence, according to Eq. (6), is a pure gradient field with ${\mathit{g}}_{\mathit{LM}}=-\sqrt{L(L+1)}{h}_L{\sigma }_{T,LM}^2$. Here, ${\sigma }_{T,LM}^2$ is the spin-0 harmonic transform of the temperature noise variance map, and

$$\begin{array}{c}\hfill \begin{array}{c}{h}_L\equiv 2\pi {\int }_{-1}^{1}d\mu {\xi }_{0,0}^{\overline{T}}(\mu ){\xi }_{0,1}^{T^{\mathrm{WF}}}(\mu )d_{-1,0}^L(\mu ),\hfill \\ {\xi }_{0,0}^{\overline{T}}(\mu )\equiv {\sum }_{\ell =100}^{2048}\left(\frac{2\ell +1}{4\pi }\right){\left(\frac{b_{\ell }}{b_{\ell }^2{C}_{\ell }^{\mathit{TT}}+N_{\ell }}\right)}_{\mathrm{fid}}{d}_{0,0}^{\ell }(\mu ),\hfill \\ {\xi }_{0,1}^{T^{\mathrm{WF}}}(\mu )\equiv {\sum }_{\ell =100}^{2048}\left(\frac{2\ell +1}{4\pi }\right){\left(\frac{b_{\ell }{C}_{\ell }^{\mathit{TT}}}{b_{\ell }^2{C}_{\ell }^{\mathit{TT}}+N_{\ell }}\right)}_{\mathrm{fid}}\sqrt{\ell (\ell +1)}{d}_{1,0}^{\ell }(\mu ).\hfill \end{array}\end{array}$$(B.2)

Here, b_ℓ, $C_{\mathcal{l}}^{TT}$, and N_ℓ are the fiducial combined beam and pixel window function, fiducial CMB spectrum, and noise spectrum used in the filtering (we use a flat noise level of 35 μK-arcmin). The contribution to the mean field from inhomogeneous pixel noise is only visible at low lensing multipoles.

At higher multipoles, 500 ≤ L ≤ 2500, the noise mean field drops sharply, but the CMB mean field remains clearly visible, at least partly because of pixelization effects. To some approximation, the mapmaking procedure simply bins individual time-ordered data into the assigned pixels. Approximating the temperature field locally as a gradient, this gives rise to an effective deflection field, which can be predicted given the full pointing information for each frequency channel. Since this sub-pixel mis-centring is the same in all simulations, it appears rather directly in the mean-field estimate (with a possible bias, since, unlike lensing, the sub-pixel deflections act on the sky signal after beam convolution). The prediction for the 217-GHz channel sub-pixel deflection spectrum is shown as the purple curve in Fig. B.1, and has the approximate power spectrum of a white-noise displacement component $C_L^{\kappa \kappa ,\text{sub-pixel}}\propto L^2$ (as expected for a large number of independent hits) with an amplitude of around 0.​​′05. The cross-spectrum of the mean-field estimate for the SMICA map with the predicted sub-pixel displacements for the 217 GHz channel recovers the expected amplitude in the range 500 ≤ L ≤ 2500 fairly well.

At yet higher multipoles, L ≥ 2500, while the CMB mean field can still be seen in cross-correlation, the most striking feature is a very sharp rise in the mean field sourced by the noise simulations. We caution that this rise is likely to be largely an artefact of the way that the simulations were constructed: as discussed in Sect. 2, the noise simulations are not exactly independent, due to small nonlinearities in the data processing causing correlations after subtraction of the common fiducial CMB and foreground realizations used to make the simulations. This introduces a contribution to the mean field given by the result of applying the lensing quadratic estimator to the common component of the noise simulations. We estimate this common component from empirical averages of the first 150 noise simulations and the last 150 noise simulations. By applying the lensing quadratic estimator to each of these simulation averages in turn, and cross-correlating the results, we attempt to isolate the mean-field power sourced by the common signal. The result is shown as the red curve in Fig. B.1; this agrees well with the sharp increase in the mean field (blue curve) seen at L ≥ 2500. The red curve contains a residual Monte Carlo noise level of 150⁻² in the mean field, which is visible at low L. On these scales the mean field of the common component is too small to be separated from the noise. In polarization, the mean field of the common component is slightly larger than in temperature; however, its contribution to the spectrum is a negligibly small fraction of the standard deviation of the band powers.

Appendix C: Covariance matrix corrections

As described in Sect. 2.2, our covariance matrix takes into account Monte Carlo errors in the various quantities obtained from simulations, increasing the covariance by roughly 10% over our conservative multipole range. These errors come mainly from the mean-field estimate and the lensing N⁽⁰⁾ bias (see Appendix A). We use N_bias = 60 simulations for the mean field and N_bias = 240 independent simulations for N⁽⁰⁾. This section describes how the impact of the Monte Carlo error is calculated. For simplicity we use the following form of the Monte Carlo-corrected final power spectrum estimator, where the Monte Carlo correction is applied additively instead of multiplicatively:

$$\begin{array}{c}\hfill \begin{array}{cc}{\widehat{C}}_L^{\varphi \varphi }\hfill & ={\widehat{C}}_L^{\mathit{dd}}-\text{RD-}{N}_L^{(0)}-N_L^{(1)}-{\langle {\widehat{C}}_L^{\varphi \varphi }\rangle }_{\mathrm{MC}}+C_L^{\varphi \varphi ,\mathrm{fid}}\hfill \\ \hfill & ={\widehat{C}}_L^{\mathit{dd}}-{\langle 4{\widehat{C}}_L^{\mathit{di}}-4{\widehat{C}}_L^{\mathit{ij}}+{\widehat{C}}_L^{\mathit{ii}}\rangle }_{N_{\mathrm{bias}}}+C_L^{\varphi \varphi ,\mathrm{fid}}\hfill \end{array}\end{array}$$(C.1)

(and neglecting the point-source correction), obtained combining Eqs. (A.1) and (A.2). Here, ${\widehat{C}}_L^{\mathit{dd}}$ is the lensing reconstruction power spectrum estimate from the data, including mean-field subtraction (i.e., it is equivalent to ${\widehat{C}}_L^{\widehat{\varphi }\widehat{\varphi }}$ defined in the main text in Eq. (8), but here we wish to emphasise that it is constructed from the data).

C.1. Monte Carlo errors on the mean field

The mean field is defined as an ensemble average ${\varphi }_{\mathit{LM}}^{\mathrm{MF}}\equiv \langle {\widehat{\mathit{g}}}_{\mathit{LM}}\rangle /{\mathcal{R}}_L^{\varphi }$, but evaluated with a finite number of simulations N_bias. The mean field only enters two terms in the reconstructed spectrum of Eq. (C.1): the data lensing reconstruction spectrum ${\widehat{C}}_L^{\mathit{dd}}$; and similarly ${\widehat{C}}_L^{\mathit{ii}}$ on the simulations used for the Monte Carlo correction calculation. To avoid Monte Carlo noise bias in these spectrum estimates, we use two independent sets of N_bias/2 simulations to subtract the mean field from each quadratic estimator before forming their cross-spectrum. Denoting the (small and independent) error on the estimated mean field for each quadratic estimate by $\delta {\varphi }_{1,2}^{\mathrm{MF}}$, all relevant terms are of the form

$$\begin{array}{c}\hfill {\widehat{C}}_L^{\varphi \varphi }\supset (\frac{{\widehat{g}}_{\mathit{LM}}}{{\mathcal{R}}_L^{\varphi }}-{\varphi }_{\mathit{LM}}^{\mathrm{MF}}-\delta {\varphi }_{1,LM}^{\mathrm{MF}}){(\frac{{\widehat{g}}_{\mathit{LM}}}{{\mathcal{R}}_L^{\varphi }}-{\varphi }_{\mathit{LM}}^{\mathrm{MF}}-\delta {\varphi }_{2,LM}^{\mathrm{MF}})}^{\ast },\end{array}$$(C.2)

where ${\widehat{g}}_{\mathit{LM}}$ is either the data lensing map ${\widehat{g}}_{\mathit{LM}}^d$ (for ${\widehat{C}}_L^{\mathit{dd}}$) or that of simulation i (${\widehat{g}}_{\mathit{LM}}^i$ for ${\widehat{C}}_L^{\mathit{ii}}$). Summing the terms from ${\widehat{C}}_L^{\mathit{dd}}$ and ${\widehat{C}}_L^{\mathit{ii}}$ at fixed L and M, the contribution to the error on the spectrum is

$$\begin{array}{c}\hfill \delta {\widehat{C}}_L^{\varphi \varphi }\supset -\delta {\varphi }_{1,LM}^{\mathrm{MF}}{[(\frac{{\widehat{g}}_{\mathit{LM}}^d}{{\mathcal{R}}_L^{\varphi }}-{\varphi }_{\mathit{LM}}^{\mathrm{MF}})-{\langle \frac{{\widehat{g}}_{\mathit{LM}}^i}{{\mathcal{R}}_L^{\varphi }}-{\varphi }_{\mathit{LM}}^{\mathrm{MF}}\rangle }_{N_{\mathrm{bias}}}]}^{\ast }+{(1\leftrightarrow 2)}^{\ast }.\end{array}$$(C.3)

We note that terms of the form $\delta {\varphi }_{1,LM}^{\mathrm{MF}}{(\delta {\varphi }_{2,LM}^{\mathrm{MF}})}^{*}$ cancel between $\delta {\widehat{C}}_L^{\mathit{dd}}$ and $\delta {\widehat{C}}_L^{\mathit{ii}}$. We proceed to calculate the averaged squared error on the spectrum. The spectrum of the mean-field error is

$$\begin{array}{c}\hfill \langle |\delta {\varphi }_{\mathit{LM}}^{\mathrm{MF}}{|}^2\rangle \simeq \frac{2}{N_{\mathrm{MF}}}(C_L^{\varphi \varphi }+N_L^{(0)}),\end{array}$$(C.4)

neglecting the N⁽¹⁾ contribution. We note that the factor of two is present since we use only half of the N_MF simulations on each leg. We also have

$$\begin{array}{cc}\hfill \langle {|(\frac{{\widehat{g}}_{\mathit{LM}}^d}{{\mathcal{R}}_L^{\varphi }}-{\varphi }_{\mathit{LM}}^{\mathrm{MF}})-{\langle \frac{{\widehat{g}}_{\mathit{LM}}^i}{{\mathcal{R}}_L^{\varphi }}-{\varphi }_{\mathit{LM}}^{\mathrm{MF}}\rangle }_{N_{\mathrm{bias}}}|}^2\rangle & \simeq (C_L^{\varphi \varphi }+N_L^{(0)})\hfill \\ \hfill & \times (1+\frac{1}{N_{\mathrm{bias}}}).\hfill \end{array}$$(C.5)

Finally, averaging over M, allowing crudely for an effective number of modes (2L + 1)f_sky, and neglecting the correction scaling with 1/N_bias, the total error induced in the reconstruction spectrum due to Monte Carlo errors in the mean-field evaluation has variance:

$$\begin{array}{c}\hfill \mathrm{var}\left({\delta {\widehat{C}}_L^{\varphi \varphi }|}_{\mathrm{MF}}\right)\simeq \frac{4}{N_{\mathrm{MF}}}\frac{{(C_L^{\varphi \varphi }+N_L^{(0)})}^2}{(2L+1)f_{\mathrm{sky}}}·\end{array}$$(C.6)

C.2. N⁽⁰⁾ Monte Carlo errors

We now consider the error on the spectrum due to the finite number N_bias of simulations in the average in Eq. (C.1),

$$\begin{array}{c}\hfill {\delta {\widehat{C}}_L^{\varphi \varphi }|}_{\mathrm{Biases}}\equiv -\delta {\langle 4{\widehat{C}}_L^{\mathit{di}}-4{\widehat{C}}_L^{\mathit{ij}}+{\widehat{C}}_L^{\mathit{ii}}\rangle }_{N_{\mathrm{bias}}}.\end{array}$$(C.7)

We proceed by approximating the reconstructions ${\widehat{g}}_{\mathit{LM}}$ as Gaussian, in which case correlations between the errors in ${\widehat{C}}_L^{\mathit{ii}}$, ${\widehat{C}}_L^{\mathit{di}}$, and ${\widehat{C}}_L^{\mathit{ij}}$ vanish. Empirical verification confirms that to a good approximation, errors in these spectra are independent, with approximate Gaussian variance

$$\begin{array}{cc}& \mathrm{var}\left(\delta {\langle {\widehat{C}}_L^{\mathit{ii}}\rangle }_{N_{\mathrm{bias}}}\right)\simeq \frac{1}{N_{\mathrm{bias}}}\frac{2{(C_L^{\varphi \varphi }+N_L^{(0)})}^2}{(2L+1)f_{\mathrm{sky}}},\hfill \end{array}$$(C.8)

$$\begin{array}{cc}& \mathrm{var}\left(\delta {\langle 2{\widehat{C}}_L^{\mathit{di}}\rangle }_{N_{\mathrm{bias}}}\right)=\mathrm{var}\left(\delta {\langle 2{\widehat{C}}_L^{\mathit{ij}}\rangle }_{N_{\mathrm{bias}}}\right)\simeq \frac{1}{N_{\mathrm{bias}}}\frac{2{\left(N_L^{(0)}\right)}^2}{(2L+1)f_{\mathrm{sky}}}.\hfill \end{array}$$(C.9)

C.3. Total Monte Carlo errors

Combining Eqs. ((C.6)–(C.9)) gives the total variance of the Monte Carlo error on the power spectrum of the reconstruction. Neglecting the absence of $C_L^{\varphi \varphi }$ in Eq. (C.9), which is accurate on almost all scales, we obtain

$$\begin{array}{c}\hfill {\sigma }_{\mathrm{MC},L}^2\simeq (\frac{2}{N_{\mathrm{MF}}}+\frac{9}{N_{\mathrm{Bias}}})\frac{2{(C_L^{\varphi \varphi }+N_L^{(0)})}^2}{(2L+1)f_{\mathrm{sky}}}·\end{array}$$(C.10)

This gives Eq. (13) in the main text, after identifying the diagonal band-power Gaussian variance with ${\sigma }_{BP}^2$. Figure C.1 shows a rough empirical estimate of the Monte Carlo errors on our MV band powers compared with the analytic estimate of Eq. (C.10). The points were estimated from the scatter of five data band powers, each constructed using 60 independent simulations (with N_bias = 12 and N_bias = 48) out of the full set of 300 simulations. The figure shows the empirical variance (divided by 5) normalized to the statistical variance of the band powers, ${\sigma }_{BP}^2$; this agrees reasonably well with the analytic estimate (the correction is small, so it does not need to be calculated accurately).

All Tables

All Figures