Issue 
A&A
Volume 641, September 2020
Planck 2018 results



Article Number  A8  
Number of page(s)  42  
Section  Cosmology (including clusters of galaxies)  
DOI  https://doi.org/10.1051/00046361/201833886  
Published online  11 September 2020 
Planck 2018 results
VIII. Gravitational lensing
^{1}
AIM, CEA, CNRS, Université ParisSaclay, Université ParisDiderot, Sorbonne Paris Cité, 91191 GifsurYvette, France
^{2}
APC, AstroParticule et Cosmologie, Université Paris Diderot, CNRS/IN2P3, CEA/lrfu, Observatoire de Paris, Sorbonne Paris Cité 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13, France
^{3}
African Institute for Mathematical Sciences, 68 Melrose Road, Muizenberg, Cape Town, South Africa
^{4}
Aix Marseille Univ, CNRS, CNES, LAM, Marseille, France
^{5}
Astrophysics Group, Cavendish Laboratory, University of Cambridge, J J Thomson Avenue, Cambridge CB3 0HE, UK
^{6}
Astrophysics & Cosmology Research Unit, School of Mathematics, Statistics & Computer Science, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban 4000, South Africa
^{7}
CITA, University of Toronto, 60 St. George St., Toronto, ON M5S3H8, Canada
^{8}
CNRS, IRAP, 9 Av. Colonel Roche, BP 44346, 31028 Toulouse Cedex 4, France
^{9}
Cahill Center for Astronomy and Astrophysics, California Institute of Technology, Pasadena, CA 91125, USA
^{10}
California Institute of Technology, Pasadena, California, USA
^{11}
Centre for Theoretical Cosmology, DAMTP, University of Cambridge, Wilberforce Road, Cambridge CB3 0WA, UK
^{12}
Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley, CA, USA
^{13}
Département de Physique Théorique, Université de Genève, 24, Quai E. Ansermet, 1211 Genève 4, Switzerland
^{14}
Département de Physique, École normale supérieure, PSL Research University, CNRS, 24 rue Lhomond, 75005 Paris, France
^{15}
Departamento de Astrofísica, Universidad de La Laguna (ULL), 38206 La Laguna, Tenerife, Spain
^{16}
Departamento de Física, Universidad de Oviedo, C/ Federico García Lorca, 18, Oviedo, Spain
^{17}
Department of Astrophysics/IMAPP, Radboud University, PO Box 9010, 6500 Nijmegen, The Netherlands
^{18}
Department of Mathematics, University of Stellenbosch, Stellenbosch 7602, South Africa
^{19}
Department of Physics & Astronomy, University of British Columbia, 6224 Agricultural Road, Vancouver, BC, Canada
^{20}
Department of Physics & Astronomy, University of the Western Cape, Cape Town 7535, South Africa
^{21}
Department of Physics and Astronomy, University of Sussex, Brighton BN1 9QH, UK
^{22}
Department of Physics, University of Helsinki, Gustaf Hällströmin katu 2a, Helsinki, Finland
^{23}
Department of Physics, Princeton University, Princeton, NJ, USA
^{24}
Department of Physics, University of California, Berkeley, CA, USA
^{25}
Department of Physics, University of California, One Shields Avenue, Davis, CA, USA
^{26}
Department of Physics, University of California, Santa Barbara, CA, USA
^{27}
Dipartimento di Fisica e Astronomia G. Galilei, Università degli Studi di Padova, Via Marzolo 8, 35131 Padova, Italy
^{28}
Dipartimento di Fisica e Scienze della Terra, Università di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{29}
Dipartimento di Fisica, Università La Sapienza, P.le A. Moro 2, Roma, Italy
^{30}
Dipartimento di Fisica, Università degli Studi di Milano, Via Celoria, 16, Milano, Italy
^{31}
Dipartimento di Fisica, Università degli Studi di Trieste, Via A. Valerio 2, Trieste, Italy
^{32}
Dipartimento di Fisica, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy
^{33}
European Space Agency, ESAC, Planck Science Office, Camino bajo del Castillo, s/n, Urbanización Villafranca del Castillo, Villanueva de la Cañada, Madrid, Spain
^{34}
European Space Agency, ESTEC, Keplerlaan 1, 2201 Noordwijk, The Netherlands
^{35}
Gran Sasso Science Institute, INFN, Viale F. Crispi 7, 67100 L’Aquila, Italy
^{36}
HEP Division, Argonne National Laboratory, Lemont, IL 60439, USA
^{37}
Haverford College Astronomy Department, 370 Lancaster Avenue, Haverford, PA, USA
^{38}
Helsinki Institute of Physics, University of Helsinki, Gustaf Hällströmin katu 2, Helsinki, Finland
^{39}
INAF – OAS Bologna, Istituto Nazionale di Astrofisica – Osservatorio di Astrofisica e Scienza dello Spazio di Bologna, Area della Ricerca del CNR, Via Gobetti 101, 40129 Bologna, Italy
^{40}
INAF – Osservatorio Astronomico di Padova, Vicolo dell’Osservatorio 5, Padova, Italy
^{41}
INAF – Osservatorio Astronomico di Trieste, Via G.B. Tiepolo 11, Trieste, Italy
^{42}
INAF, Istituto di Radioastronomia, Via Piero Gobetti 101, 40129 Bologna, Italy
^{43}
INAF/IASF Milano, Via E. Bassini 15, Milano, Italy
^{44}
INFN – CNAF, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{45}
INFN, Sezione di Bologna, Viale Berti Pichat 6/2, 40127 Bologna, Italy
^{46}
INFN, Sezione di Ferrara, Via Saragat 1, 44122 Ferrara, Italy
^{47}
INFN, Sezione di Milano, Via Celoria 16, Milano, Italy
^{48}
INFN, Sezione di Roma 1, Università di Roma Sapienza, Piazzale Aldo Moro 2, 00185 Roma, Italy
^{49}
INFN, Sezione di Roma 2, Università di Roma Tor Vergata, Via della Ricerca Scientifica, 1, Roma, Italy
^{50}
Imperial College London, Astrophysics group, Blackett Laboratory, Prince Consort Road, London SW7 2AZ, UK
^{51}
Institut d’Astrophysique Spatiale, CNRS, Univ. ParisSud, Université ParisSaclay, Bât. 121, 91405 Orsay Cedex, France
^{52}
Institut d’Astrophysique de Paris, CNRS (UMR7095), 98bis boulevard Arago, 75014 Paris, France
^{53}
Institute Lorentz, Leiden University, PO Box 9506, Leiden 2300, The Netherlands
^{54}
Institute of Astronomy, University of Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{55}
Institute of Theoretical Astrophysics, University of Oslo, Blindern, Oslo, Norway
^{56}
Instituto de Astrofísica de Canarias, C/Vía Láctea s/n, La Laguna, Tenerife, Spain
^{57}
Instituto de Física de Cantabria (CSICUniversidad de Cantabria), Avda. de los Castros s/n, Santander, Spain
^{58}
Istituto Nazionale di Fisica Nucleare, Sezione di Padova, Via Marzolo 8, 35131 Padova, Italy
^{59}
Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena, CA, USA
^{60}
Jodrell Bank Centre for Astrophysics, Alan Turing Building, School of Physics and Astronomy, The University of Manchester, Oxford Road, Manchester M13 9PL, UK
^{61}
Kavli Institute for Cosmology Cambridge, Madingley Road, Cambridge CB3 0HA, UK
^{62}
Kavli Institute for the Physics and Mathematics of the Universe (Kavli IPMU, WPI), UTIAS, The University of Tokyo, Chiba 2778583, Japan
^{63}
Laboratoire de Physique Subatomique et Cosmologie, Université GrenobleAlpes, CNRS/IN2P3, 53 rue des Martyrs, 38026 Grenoble Cedex, France
^{64}
Laboratoire de Physique Théorique, Université ParisSud 11 & CNRS, Bâtiment 210, 91405 Orsay, France
^{65}
Low Temperature Laboratory, Department of Applied Physics, Aalto University, Espoo 00076, Aalto, Finland
^{66}
MaxPlanckInstitut für Astrophysik, KarlSchwarzschildStr. 1, 85741 Garching, Germany
^{67}
Mullard Space Science Laboratory, University College London, Surrey RH5 6NT, UK
^{68}
NAOCUKZN Computational Astrophysics Centre (NUCAC), University of KwaZuluNatal, Durban 4000, South Africa
^{69}
National Centre for Nuclear Research, ul. L. Pasteura 7, 02093 Warsaw, Poland
^{70}
Purple Mountain Observatory, No. 8 Yuan Hua Road, 210034 Nanjing, PR China
^{71}
SISSA, Astrophysics Sector, via Bonomea 265, 34136, Trieste, Italy
^{72}
San Diego Supercomputer Center, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093, USA
^{73}
School of Chemistry and Physics, University of KwaZuluNatal, Westville Campus, Private Bag X54001, Durban, 4000, South Africa
^{74}
School of Physical Sciences, National Institute of Science Education and Research, HBNI, Jatni, 752050 Odissa, India
^{75}
School of Physics and Astronomy, Cardiff University, Queens Buildings, The Parade, Cardiff CF24 3AA, UK
^{76}
School of Physics and Astronomy, Sun Yatsen University, 2 Daxue Rd, Tangjia, Zhuhai, PR China
^{77}
School of Physics and Astronomy, University of Nottingham, Nottingham NG7 2RD, UK
^{78}
School of Physics, Indian Institute of Science Education and Research Thiruvananthapuram, Maruthamala PO, Vithura, Thiruvananthapuram, 695551 Kerala, India
^{79}
School of Physics, The University of New South Wales, Sydney, NSW 2052, Australia
^{80}
Simon Fraser University, Department of Physics, 8888 University Drive, Burnaby, BC, Canada
^{81}
Sorbonne Université, Observatoire de Paris, Université PSL, École normale supérieure, CNRS, LERMA, 75005 Paris, France
^{82}
Sorbonne Université, UMR7095, Institut d’Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
^{83}
Space Research Institute (IKI), Russian Academy of Sciences, Profsoyuznaya Str, 84/32, Moscow 117997, Russia
^{84}
Space Science Data Center – Agenzia Spaziale Italiana, Via del Politecnico snc, 00133 Roma, Italy
^{85}
Space Sciences Laboratory, University of California, Berkeley, California, USA
^{86}
The Oskar Klein Centre for Cosmoparticle Physics, Department of Physics, Stockholm University, AlbaNova, 106 91 Stockholm, Sweden
^{87}
Université de Toulouse, UPSOMP, IRAP, 31028 Toulouse Cedex4, France
^{88}
Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
Received:
17
July
2018
Accepted:
18
July
2019
We present measurements of the cosmic microwave background (CMB) lensing potential using the final Planck 2018 temperature and polarization data. Using polarization maps filtered to account for the noise anisotropy, we increase the significance of the detection of lensing in the polarization maps from 5σ to 9σ. Combined with temperature, lensing is detected at 40σ. We present an extensive set of tests of the robustness of the lensingpotential power spectrum, and construct a minimumvariance estimator likelihood over lensing multipoles 8 ≤ L ≤ 400 (extending the range to lower L compared to 2015), which we use to constrain cosmological parameters. We find good consistency between lensing constraints and the results from the Planck CMB power spectra within the ΛCDM model. Combined with baryon density and other weak priors, the lensing analysis alone constrains σ_{8}Ω_{m}^{0.25} = 0.589 ± 0.020 (1σ errors). Also combining with baryon acoustic oscillation data, we find tight individual parameter constraints, σ_{8} = 0.811 ± 0.019, H_{0} = 67.9_{−1.3}^{+1.2} km s^{−1} Mpc^{−1}, and Ω_{m} = 0.303_{−0.018}^{+0.016}. Combining with Planck CMB power spectrum data, we measure σ_{8} to better than 1% precision, finding σ_{8} = 0.811 ± 0.006. CMB lensing reconstruction data are complementary to galaxy lensing data at lower redshift, having a different degeneracy direction in σ_{8} − Ω_{m} space; we find consistency with the lensing results from the Dark Energy Survey, and give combined lensingonly parameter constraints that are tighter than joint results using galaxy clustering. Using the Planck cosmic infrared background (CIB) maps as an additional tracer of highredshift matter, we make a combined Planckonly estimate of the lensing potential over 60% of the sky with considerably more smallscale signal. We additionally demonstrate delensing of the Planck power spectra using the joint and individual lensing potential estimates, detecting a maximum removal of 40% of the lensinginduced power in all spectra. The improvement in the sharpening of the acoustic peaks by including both CIB and the quadratic lensing reconstruction is detected at high significance.
Key words: gravitational lensing: weak / cosmological parameters / cosmic background radiation / largescale structure of Universe / cosmology: observations
© ESO 2020
1. Introduction
Gravitational lensing distorts our view of the lastscattering surface, generating new nonGaussian signals and Bmode 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 higherorder effects. The sharplydefined 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 largescale 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^{1} lensing reconstruction analysis, giving the most significant detection of lensing to date over 70% of the sky. We also give new results for polarizationbased 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 fullsky lensing reconstruction based on the nominalmission temperature data. In the 2015 analysis (Planck Collaboration XV 2016; hereafter PL2015) this was updated to include the fullmission temperature data, as well as polarization, along with a variety of analysis improvements. The final fullmission analysis presented here uses essentially the same data as PL2015: most of the maplevel 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 polarizationonly 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 minimumvariance 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 highfrequency 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 lensingreconstruction 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 Bmode polarization after delensing.
– Using the Planck lensing likelihood alone we place a 3.5% constraint on the parameter combination . This has comparable statistical power to current constraints from galaxy lensing, but the highredshift CMB source plane gives a different degeneracy direction compared to the 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 σ_{8}, Ω_{m}, and H_{0}. We can also break the degeneracy by combining our likelihood with the firstyear lensing results from the Dark Energy Survey (DES; Troxel et al. 2018), giving the tightest lensingonly 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 powerspectrum 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 Bmode 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.
Fig. 1. Mollweide projection in Galactic coordinates of the lensingdeflection reconstruction map from our baseline minimumvariance (MV) analysis. We show the Wienerfiltered displacementlike scalar field with multipoles , corresponding to the gradient mode (or E mode) of the lensing deflection angle. Modes with L < 8 have been filtered out. 
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 foregroundcleaned 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 powerspectrum 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}.
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 Wienerfiltered X = T, E, and B CMB multipoles, as well as inversevarianceweighted CMB maps. The latter maps can be obtained easily from the Wienerfiltered multipoles by dividing by the fiducial CMB power spectra before projecting onto maps. We write the observed temperature T and polarization (written as the spin ±2 combination of Stokes parameters _{±2}P ≡ Q ± iU) pixelized data as
where T, E, and B on the righthand side are the multipole coefficients of the true temperature and E and Bmode polarization. The matrix 𝒴 contains the appropriate (spinweighted) spherical harmonic functions to map from multipoles to the sky, and the matrix ℬ accounts for the realspace operations of beam and pixel convolution. We further use the notation 𝒯 ≡ ℬ𝒴 for the complete transfer function from multipoles to the pixelized sky. The Wienerfiltered multipoles are obtained from the pixelized data as
where the pixelspace covariance is Cov = 𝒯C^{fid}𝒯^{†} + N. Here, C^{fid} is a fiducial set of CMB spectra and N is the pixelspace noise covariance matrix, which we approximate as diagonal. As in previous releases, our baseline results use independentlyfiltered temperature and polarization maps (i.e., we always neglect in Cov^{−1} 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 multigridpreconditioned conjugategradient 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 fullsky transfer function matches the one seen empirically on the filtered simulations, with a minimal impact on the band powers. We also present new more optimallyfiltered results, as discussed in Sect. 2.3.
Construction of the quadratic lensing estimator. We determine from pairs of filtered maps, and our implementation now follows Carron & Lewis (2017). This differs slightly from PL2015, allowing us to produce minimumvariance (MV) estimators from filtered maps much faster, which is useful given the variety of tests performed for this release. We calculate a spin1 realspace (unnormalized) lensing displacement estimate
where ð is the spinraising operator, and the presubscript s on a field denotes the spin. The quadratic estimator involves products of the realspace inversevariance filtered maps
and the gradients of the Wienerfiltered maps
The deflection estimate is decomposed directly into gradient (g) and curl (c) components by using a spin1 harmonic transform, where the gradient piece contains the information on the lensing potential and the curl component is expected to be zero:^{3}
By default we produce three estimators, namely temperatureonly (s = 0), polarizationonly (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 temperaturepolarization coupling is neglected in the C^{fid} factor that appears in the Wienerfilter of Eq. (2) for temperature and polarizationonly estimators, but is included in the MV reconstruction. We use lensed CMB spectra in Eq. (2) to make the estimator nearly unbiased to nonperturbative 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.
Meanfield 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 maplevel signal expected from mask, noise, and other anisotropic features of the map in the absence of lensing; we subtract this meanfield 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
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: crossspectra between reconstructions from masked simulations and the true input lensing realizations match expectations to subpercent 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 powerspectrum estimation, we use crossspectrum estimators of maps with independent Monte Carlo noise on the meanfield 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 tradeoff between uncertainties in the Monte Carlo estimates of the mean field and biases that affect the reconstruction bandpower 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 socalled N^{(0)} and N^{(1)} lensing biases, as well as pointsource contamination. For a pair of lensing map estimates and , we use the same simple crossspectrum estimator as in PL2015,
from which biases are subtracted:
Lensing powerspectrum estimation is designed to probe the connected 4point function of the data that is induced by lensing. The combination of the meanfield subtraction and the first bias term (N^{(0)}) in Eq. (9) subtracts the disconnected signal expected from Gaussian fluctuations even in the absence of lensing, and is calculated using the same realizationdependent N^{(0)} (RDN^{(0)}) estimator described in PL2015 (and summarized in Appendix A for the baseline cases). The term subtracts an signal term (N^{(1)}) coming from nonprimary couplings of the connected 4point function (Kesden et al. 2003), and in our baseline analysis this is calculated using a fullsky analytic approximation in the fiducial model as described in Appendix A. The signaldependence of this term is handled consistently in the likelihood as described in Sect. 3.2. We tested an alternative simulationbased N^{(1)} calculation and a deconvolution technique as described in Sect. 4.8. The pointsource (PS) bias term subtracts the contribution from the connected 4point 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^{(0)} 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 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 (), we define our bandpower estimates
The binning functions use an approximate inversevariance weighting to produce roughly optimal signal amplitudes:
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,
This choice ensures that the binned fiducial spectrum goes exactly through the fiducial model at L = L_{b}, so plotting bandpower 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 righthand 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 powerspectrum estimator. The simulationaveraged band powers in Eq. (10) are built from simulations as from the data according to Eq. (9), but with a cheaper Monte Carlo N^{(0)} (MCN^{(0)}) estimation described in Appendix A, and no pointsource 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 pointsource contamination for the same sky fraction. The amplitude of the pointsource correction decreased by a factor of 1.9, and the detection of this pointsource 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 SMICAbased confidence mask;^{4} the 2015 70% Galactic mask; and the pointsource 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^{5}. 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 foregroundcleaned 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 μKarcmin 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 crosscorrelation 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 . The choice of a multiplicative correction is more appropriate for modemixing 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 bandpower 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^{(1)} bias that we subtract, and include L ≤ 2500 in the linear correction to account for the modeldependence of N^{(1)} relative to the fiducial lensing power. The 2015 MV likelihood contained an almost inconsequential error in the calculation of the response of N^{(1)} 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 “lensingonly” 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 spatiallyflat ΛCDM cosmology with: baryon density ω_{b} ≡ Ω_{b}h^{2} = 0.02216; cold dark matter density ω_{c} ≡ Ω_{c}h^{2} = 0.1203; two massless neutrinos and one massive with mass 0.06 eV; Hubble constant H_{0} = 100h km s^{−1} Mpc^{−1} 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^{−1}) A_{s} = 2.119 × 10^{−9}; and Thomson optical depth through reionization τ = 0.060.
2.2. Covariance matrix
Our bandpower covariance matrix is obtained from the FFP10 simulation suite. Out of the 300 simulations, 60 are used for the meanfield 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 bandpower covariance matrix. It is impractical to perform the same, expensive, realizationdependent N^{(0)} (RDN^{(0)}) subtraction on all these simulations for evaluation of the covariance matrix, and therefore, as in previous releases, we use a cheaper semianalytic calculation, as detailed in PL2015, which only requires empirical spectra of the CMB data. This semianalytic calculation is only accurate to 1–2%, which is not enough for debiasing where subpercent 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 meanfield, RDN0, and MC corrections. As detailed in Appendix C, for Planck noise levels the additional variance caused by the finite number of simulations can be written to a good approximation in terms of the bandpower statistical^{6} errors as
Here, N_{MF} is the number of simulations entering the meanfield 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^{(0)} 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)
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 bandpower bins N_{bins} = 9 and 1/α_{cov} = 1.035 (i.e., effectively a 3.5% increase in bandpower covariance); and one on the aggressive multipole range 8 ≤ L ≤ 2048, for which N_{bins} = 16 and 1/α_{cov} = 1.06. After our semianalytical realizationdependent debiasing, the covariance matrix shows no obvious structure on either multipole range. We find all individual crosscorrelation coefficients to be smaller than 10% and consistent with zero, a constraint limited by the number of simulations available. We choose to include the offdiagonal elements in the likelihood.
Following PL2015 we subtract a pointsource template correction from our band powers, with an internally measured amplitude (the pointsource shotnoise trispectrum ). We neglect the contribution to the error from pointsource subtraction uncertainty, since for this release the estimated error on would formally inflate bandpower 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 polarizationonly 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.
Fig. 2. Noisevariance maps (shown as noise rms in μKarcmin) that we use to filter the SMICA CMB maps that are fed into the quadratic estimators when performing inhomogeneous filtering. The upper panel shows the temperature noise map, with median over our unmasked sky area of 27 μKarcmin. We use a common noise map for Q and U polarization, also neglecting Q and U noise correlations, shown in the lower panel, which spans an entire order of magnitude, with median 52 μKarcmin (larger than times the temperature noise because not all the Planck detectors are polarized). In temperature, the variance map has a homogeneous (approximately 5 μKarcmin) contribution from the isotropic additional Gaussian power that we add to the simulations to account for residual foreground contamination. 
The SMICA CMB map is constructed from Planck frequency maps using isotropic weights 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 , and neglecting differences between and , we have
where
and are reduced Wigner dmatrices. This equation follows from transforming the noise maps at each frequency f, which have pixel variance and are further assumed uncorrelated across frequencies, into their E and Bmodes, applying the SMICA weights and , and transforming back to Q and U in pixel space. Averaging the pixel variances of the resulting Q and U noise maps yields . Physically, the variances of the frequency maps are combined nonlocally with kernels that derive from the convolutions implied by the SMICA weights. In temperature,
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 signaltonoise 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 highresolution 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 crosscorrelation 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 maplevel 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 fullsky optimal lensing spectrum estimation is obtained by inversevariance weighting correctlynormalized 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^{(0)}. Therefore, whenever the reconstruction noise dominates the lens cosmic variance in the bandpower errors, inversevariance 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 independentpatch 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))
Applying a single fiducial response on the fullsky estimate, the local estimate in a patch centred on is biased by a factor , where is the true response according to the local temperature and polarization filtering noise levels. We may then write the multipoles of the fullsky lensing map as a sum of multipoles extracted over the patches,
where each unbiased component 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
and equivalently for the estimator power spectrum,
The spectrumlevel correction of Eq. (21) is only close to the squared maplevel 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 signaldominated 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 inhomogeneouslyfiltered polarization band powers, together with the prediction from Eq. (21). The agreement is visually very good, with a residual at lowL that originates from masking, also found on our baseline, homogeneouslyfiltered MV band powers (orange points). This largescale 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 powerspectrum estimator in Eq. (8) does not attempt to remove any modemixing effect of masking on the lensing estimate, we note that the largescale MC correction is not simply just a ϕ modemixing effect: using a pseudoC_{ℓ} inversion to construct the band powers from the masked ϕ map makes almost no difference to the MC correction.
Fig. 3. Monte Carloderived multiplicative normalization corrections for the polarization reconstruction, using inhomogeneous filtering (blue points). Band powers are divided by these numbers to provide our final estimates. This correction is not small, and is sourced by the large spatial variation of the estimator response; however, it is very well reproduced by the approximate analytic model of Eq. (21), shown as the dashed blue line. The correction for our baseline MV band powers using homogeneous filtering is shown as the orange points. 
As part of the Planck 2018 release, we make available products containing both versions of lensing maps. Optimallyweighted maps can be used if S/N is critical, while the isotropicallyweighted maps can be used to simplify crosscorrelation analyses if required. Our baseline results and likelihoods are still derived from the simpler isotropic filtering, but there is a substantial improvement in the polarizationonly 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. Lensingreconstruction map and power spectrum
In Fig. 1 we show our baseline Wienerfiltered minimumvariance lensing deflection estimate from the Planck temperature and polarization SMICA CMB maps. This is shown as a map of
where is the lensing potential power spectrum in our fiducial model and is the noise power spectrum of the reconstruction. The quantity is equivalent to the Wienerfiltered gradient mode (or E mode) of the lensing deflection angle. For powerspectrum estimates we plot , so that a map of α has the same relation to the plotted power spectrum as the CMB temperature map does to . As in 2015 we exclude L < 8 due to the high sensitivity to the meanfield 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.
Fig. 4. Planck 2018 lensing reconstruction band powers (values and multipole ranges are listed in Table 1). Left: minimumvariance (MV) lensing band powers, shown here using the aggressive (blue, 8 ≤ L ≤ 2048) and conservative (orange, 8 ≤ L ≤ 400) multipole ranges. The dots show the weighted bin centres and the fiducial lensing power spectrum is shown as the black line. Right: comparison of polarizationonly band powers using homogeneous map filtering (blue boxes, with dots showing the weighted bin centres) and the more optimal inhomogeneous filtering (orange error bars). The inhomogeneous filtering gives a scaledependent increase in S/N, amounting to a reduction of 30% in the error on the amplitude of the power spectrum over the conservative multipole range shown. The black line is the fiducial lensing power spectrum. 
In addition to the MV reconstruction, we also provide temperatureonly results, as well as two variants of the polarizationonly 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 polarizationonly 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 bandpower 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 inversevariance weighting of the unbinned spectra. The fiducial band powers can therefore show slight variations because the noise varies between the different reconstructions.
Lensingreconstruction powerspectrum bandpower amplitudes and errors for the temperatureonly (TT), combined temperatureandpolarization minimumvariance (MV), and polarization estimators (PP).
Section 3.2.1 introduces our new lensingonly likelihood, marginalizing over the CMB spectra. Using this likelihood to obtain lensing amplitude summary statistics , (with for equal to the bestfit ΛCDM model to the Planck temperature and polarization power spectra and the reconstructed lensing power^{7}), we obtain
over the conservative multipole range L = 8–400. This corresponds to a slightly higher value of the lensing spectrum than PL2015 (for which ) 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
As discussed in detail in Sect. 4, the highL 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
Amplitude statistics for the polarizationonly reconstructions are as follows (neglecting the CMB marginalization and other very small likelihood linear corrections):
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 scaledependent, with most gain achieved on small scales. This behaviour is consistent with analytic expectations, calculated using the independentpatch 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.
Fig. 5. Planck 2018 lensing powerspectrum band powers (pink boxes) over the aggressive multipole range. The 2015 analysis band powers (green) were calculated assuming a slightly different fiducial model and have not been (linearly) corrected to the 2018 model. Also shown are recent measurements by the ACTPol (Sherwin et al. 2017), SPTpol (Story et al. 2015), and SPTSZ (Simard et al. 2017) collaborations. The SPTSZ measurement is not completely independent, since the SPTSZ reconstruction also uses temperature data from Planck, but with subdominant weight over the smaller sky area used. The black line shows the lensing potential power spectrum for the ΛCDM bestfit parameters to the Planck 2018 likelihoods (Planck TT,TE,EE+lowE, which excludes the lensing reconstruction). 
3.2. Likelihood and parameter constraints
We construct a lensing likelihood from the powerspectrum 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 nulltest 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 powerspectrum band powers are corrected perturbatively for changes in normalization and N^{(1)} due to parameterdependent 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^{8}
where Σ is the covariance matrix and the binning functions are defined in Eq. (11). The binned “theory” power spectrum for cosmological parameters θ is given in the linear approximation by
where a sums over both the and CMB powerspectra terms, and the linear correction matrix can be precomputed in the fiducial model. The linear correction accounts for the N^{(1)} dependence on , and the dependence of the lensing response and N^{(1)} on the CMB power spectra; explicitly,
where derivatives are understood not to act on the lensing contribution to lensed power spectra, X is one of the CMB power spectra, and 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 powerspectrum amplitude from CMB powerspectrum measurements, this means that the 2018 CMB likelihoods can constrain amplituderelated 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 powerspectrum 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 “lensingonly” 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 powerspectrum data. We now remove any dependence on the theoretical model of the CMB power spectra by marginalizing out the theoretical by approximating their distribution as Gaussian, where up to a constant
Here, C^{CMB} and are vectors of CMB TT, TE, and EE powerspectrum values at each multipole, with 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,
and the theory spectrum shifted by the linear correction for the observed CMB power,
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 . 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 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 powerspectrum 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 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 powerspectrum 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 Ω_{b}h^{2} = 0.0222 ± 0.0005, motivated by D/H measurements in quasar absorptionline systems combined with the predictions of bigbang nucleosynthesis (BBN)^{9}. The exact choice of Ω_{b}h^{2} prior has very little effect on lensingonly constraints, but the prior is useful to constrain the sound horizon (since this has a weak but important dependence on Ω_{b}h^{2}) 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^{10}, 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 lensingonly ΛCDM constraint on σ_{8} and Ω_{m}. As discussed in detail in PL2015, the lensing data constrain a narrow band in the 3D σ_{8} − Ω_{m} − H_{0} parameter space, corresponding to
Fig. 6. Constraints on Ω_{m} and σ_{8} in the baseΛCDM model from CMB lensing alone (posterior sample points coloured by the value of the Hubble constant in units of km s^{−1} Mpc^{−1}) using the priors described in the text. Grey bands give corresponding 1σ and 2σ lensingonly constraints using the approximate fit . The joint 68% and 95% constraints from CMB lensing with the addition of BAO data (Beutler et al. 2011; Ross et al. 2015; Alam et al. 2017) are shown as the dashed contours, and the constraint from the Planck CMB powerspectrum data is shown for comparison as the solid contours. 
or the tighter and slightly less priordependent 2% constraint
The allowed region projects into a band in the Ω_{m}–σ_{8} plane with
The corresponding result using a fixed CMB fit for the CMB power spectrum is , which is consistent with the similar constraint, , found in PL2015. The roughly 0.25σ shift down in this parameter is consistent with the slight increase in because of the anticorrelation of with the lensing deflection power, as discussed in PL2015. While the tight threeparameter 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 powerspectrum constraints give slightly lower values of σ_{8} compared to the 2015 analysis due to the lower optical depth, which increases the overlap between the lensingonly and CMB powerspectrum contours, making them very consistent within the ΛCDM model.
Parameter constraints for different lensing datasets and priors, with and without galaxy BAO data.
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 Ω_{b}h^{2} so that the sound horizon is fairly well constrained, we can break the main degeneracy and constrain individual parameters, giving the ΛCDM constraints
The value of the Hubble constant inferred here assuming ΛCDM is in good agreement with other inverse distanceladder 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 modelindependent values obtained using distanceladder measurements (Riess et al. 2018).
Massive neutrinos suppress the growth of structure on scales smaller than the neutrino freestreaming 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 acousticscale 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 highprecision 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:
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:
3.2.2. Joint Planck parameter constraints
The CMB lensing power spectrum is consistent with expectations from the Planck CMB powerspectrum 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 σ_{8}–Ω_{m} plane, giving the tight (subpercent) amplitude ΛCDM result
Fig. 7. Constraints on Ω_{m} and σ_{8} in the baseΛCDM model from Planck temperature and polarization power spectra (red), and the tighter combined constraint with CMB lensing (blue). The dashed line shows the joint result when the reionization redshift is restricted to z_{re} > 6.5 to be consistent with observations of highredshift quasars (Fan et al. 2006). Contours contain 68% and 95% of the probability. 
This result uses the temperature and Emode 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 Emode polarization pulling the amplitude back to values that are more consistent with the lensing analyses. Even with the lowℓ polarization, the powerspectrum data prefer around 1σ higher σ_{8} than the lensing data alone, with the joint constraint lying in between. Values of σ_{8} inferred from the CMB powerspectrum data in a given theoretical model cannot be arbitrarily low because the reionization optical depth cannot be arbitrarily small; observations of highredshift 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 CMBlensing and galaxylensing constraints
Cosmic shear of galaxies can be used to measure the lensing potential with lowerredshift 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 crosscorrelation signal has been detected with a variety of lensing data (Hand et al. 2015; Liu & Hill 2015; Kirk et al. 2016; HarnoisDéraps et al. 2017; Singh et al. 2017) and may ultimately be a useful way to improve parameter constraints and constrain galaxylensing 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 crosscorrelation; 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 lensingonly analysis described in Sect. 3.2.1^{11}.
Figure 8 shows the Planck and DES lensingonly ΛCDM constraints in the Ω_{m}–σ_{8} 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 and CMB lensing constraining ). 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 powerspectrum results. The joint result in the Ω_{m}–σ_{8} plane constrains the combined direction
Fig. 8. Constraints on Ω_{m} and σ_{8} in the base ΛCDM model from DES galaxy lensing (green), Planck CMB lensing (grey), and the joint constraint (red). The Planck powerspectrum constraint is shown in blue. Here, we adopt cosmological parameter priors consistent with the CMB lensingonly analysis, which differ from the priors assumed by the DES collaboration (Troxel et al. 2018). The odd shape of the DES lensing and joint contours is due to a nontrivial degeneracy with the intrinsic alignment parameters, giving a region of parameter space with large negative intrinsic alignment amplitudes that cannot be excluded by current lensing data alone (but is reduced by different choices of cosmological parameter priors; see PCP18). Contours contain 68% and 95% of the probability. 
although the posterior is not very Gaussian due to a nontrivial DES lensingintrinsicalignment parameter degeneracy. If instead we adopt the cosmological parameter priors of Troxel et al. (2018), but fixing the neutrino mass, then the lensingonly joint result is more Gaussian, with ; this is tighter than the constraint obtained by Troxel et al. (2018) in combination with galaxy clustering data^{12}.
Cosmological parameter degeneracies (and degeneracies with intrinsicalignment and other nuisance parameters) limit the precision of the DES lensingonly 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
Fig. 9. ΛCDM model constraints on Ω_{m}, σ_{8}, and H_{0} (in units of km s^{−1} Mpc^{−1}) from DES galaxy lensing+BAO (green), Planck CMB lensing+BAO (grey), and the joint constraint of DES lensing, CMB lensing, and galaxy BAO (red). All these results use the standard lensing priors described in the text, including the baryon density prior Ω_{b}h^{2} = 0.0222 ± 0.0005. The Planck powerspectrum constraints are shown in blue. Contours contain 68% and 95% of the probability. 
also consistent with the main Planck powerspectrum parameter analysis. However, as shown in Fig. 9, even after combining with BAO, the DES lensing results have a substantial remaining σ_{8}–Ω_{m}–H_{0} degeneracy that limits the individual constraints. For CMB lensing the situation is much better, since the CMB lensing comes from a single welldefined source redshift plane that leads to a tight σ_{8}–Ω_{m}–H_{0} degeneracy that intersects with the BAO Ω_{m}–H_{0} 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:
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 nulltest 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 bestfit 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 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 nonzero 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 deflectionangle power spectrum at L ≃ 60, but is noise dominated on smaller scales. However, the cosmic infrared background (CIB) is a highredshift 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 degreescale 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 highersensitivity internal CMBpolarizationbased reconstructions are available. The main difficulty in using the CIB as a tracer is contamination by Galactic dust and modeling of the crosscorrelation 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 Planckbased 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) componentseparation 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 545GHz CIB map only. Figure 10 shows the crossspectra of the lensing temperatureonly, polarizationonly and MV reconstructions with the GNILC 545GHz map. The bispectrum between residual foregrounds in the CMB temperature maps that enter the quadratic estimator and the CIB can potentially bias the CIBlensing crossspectra. 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 thermalSZ 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.
Fig. 10. Empirical crossspectra between our lensing reconstructions (temperatureonly in blue, polarizationonly in orange, and MV in green) and the GNILC 545GHz CIB map. The red points use the temperature lensing reconstruction from the thermalSZdeprojected SMICA CMB map. The dashed line shows the smooth spline fit to the green points, which we use to weight the lensing tracers when forming our combined lensing potential estimate. 
We build lensing estimates as follows. Consider a set of tracers I^{i} of the lensing potential, with autopower spectra , crosscorrelation coefficient matrix ρ_{L}, and crosscorrelation coefficients to the true lensing potential . 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 optimallyfiltered potential map . We can write this in terms of (isotropic) weights w_{L} acting in harmonic space on renormalized maps with unit spectra, through
where the optimal weights are given at each lensing multipole by
The expected resulting squared crosscorrelation coefficient is the weighted sum . To compare the noise in the joint reconstruction to the quadraticestimator reconstruction noise , we can also define an effective noise level N_{L} so that
We adopt a purely empirical approach to obtain the weights, together with the autospectra on the righthand side of Eq. (46), and take as our FFP10 fiducial lensing spectrum. For each pair of tracers, we perform a bicubic spline fit across 12 bins of the observed crosscorrelation coefficient and autospectra, with errors for each bin determined from the observed scatter. The spectra and crossspectra entering the crosscorrelation coefficient are calculated by deconvolving the effects of the mask from pseudoC_{ℓ} spectra (Wandelt et al. 2001), using a 12′ apodization window. To estimate the crosscorrelation 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 crossspectra of the CIB maps to the quadratic estimator are then used as proxies to the true C^{CIB ϕ} at the corresponding frequency.
These CIBlensing crossspectra 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 crosscorrelate 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 crosscorrelation coefficients 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) crosscorrelation coefficients of the true lensing map to the MV estimator (blue), together with that of the 545GHz GNILC map (orange), and their combination (green). Delensing efficiencies, , 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 545GHz GNILC lensing estimate, and their combination, showing how the CIB particularly improves the S/N at smaller scales.
Fig. 11. Top: expected crosscorrelation coefficient of the true lensing potential to the minimumvariance (MV) estimator, the 545GHz GNILC CIB map (orange), and their combination (green). To build these curves, the crossspectrum of the CIB map to the true lensing potential is approximated by its crossspectrum to the MV (quadraticestimator) lensing reconstruction, displayed in Fig. 10, and GNILC CIB curves are smooth spline fits. The red curve shows for comparison the crosscorrelation coefficient of the CIB to the MV quadratic estimator. This crosscorrelation is suppressed by instrument noise, foreground residuals, and shot noise in the CIB map and reconstruction noise in the lensing quadratic estimator. The black curves show the lensing kernels that contribute to the lensing of the temperature (solid), Emode polarization (dashed), and Bmode polarization (dashdotted) power spectra, as described in the text. Bottom: effective reconstruction noise levels N_{L} for each of these tracers, as defined by Eq. (48), and, for comparison, the theoretical lensing spectrum . The quadratic estimator reconstruction noise is slightly underestimated at low L in this figure, since we have neglected Monte Carlo corrections when combining the tracers. 
We apply the weights given by Eq. (47) to the pseudoa_{LM}s of the tracers, obtained using the apodized mask, keeping afterwards the multipole range 10 ≤ L ≤ 2000. We do not perform any further lowL 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 lowL 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 smallscale modes by using the phase information about the sky realization from the CIB, with amplitude determined from crosscorrelation 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).
Fig. 12. Comparison of lensing maps constructed from the minimumvariance quadratic estimator alone (upper panel) and in combination with the CIB, as traced by the 545GHz GNILC frequency map (lower panel). The combination is performed on 60% of the sky, defined as the union of the lensing mask and the GNILC mask. Maps show the orthographic projection of the Wienerfiltered displacement E mode, the scalar field with multipoles , with 10 ≤ L ≤ 2000. The left and right panels are centred on the north and south Galactic poles, respectively. While the two reconstruction maps are clearly strongly correlated, the combined map has substantially more smallscale power, due to the higher S/N of the CIB on small scales. 
The joint CIBlensing 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 Bmode polarization from lensing is present on all scales, and for small tensortoscalar 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 Bmode 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 Emode polarization maps, so that the corresponding peaksmoothing 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 Bmode 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 , 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 smallerscale 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 crosscorrelations. 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}:
where the noise is generated according to the reduced covariance
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 . 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. Cdelensing
Using a filtered version of the Emode 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 Bmode power. To reduce the Bmode 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 Bmode power depends on the quality of the tracer and the noise in the Emode 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^{2}, that is, 35% of the lensing can be removed using GNILC, with a moderate improvement in combination with the Planck internal reconstruction. The Planck smallscale polarization data are noisy, and at CMB Emode scale ℓ ≃ 500, from where the largescale Bpower takes significant contributions, the Emode filter quality (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 Bmode power. This is in contrast to current groundbased 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 Wienerfiltered Emode map, the same as produced by our pipeline for input to the quadratic estimators (E^{WF} in Eq. (2)). We use the inhomogeneouslyfiltered 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 Emode map, and simply remap the polarization according to the deflection angle determined by the lensing tracer, that is,
After projecting into the Bmode component, this forms our template for the lensed B modes. The simplified notation denotes the displacement of length on the sphere along the geodesic defined by , including the small rotation of the sphericalpolar polarization components induced by the parallel transport of the polarization tensor (Challinor & Chon 2002; Lewis 2005). We perform the interpolation using the Python curvedsky lensing tools from LENSIT^{13}, with a bicubic spline interpolation algorithm. The observed Emode map we use to build the template is itself lensed, while ideally one would use an unlensed Emode map. In principle, we could try and improve on that by obtaining an optimallyfiltered E map including our deflection tracer in the Wiener filtering, using the algorithm developed by Carron & Lewis (2017); however, the difference is a secondorder effect in the deflection, and we neglect it here.
The crosscorrelations of these templates to the Bmode CMB map provide strong indirect detections of the lensing Bmode power in the data. The crossspectra 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 crosscorrelated to our Bmode map, obtained from the Wienerfiltered after dividing by , where b_{ℓ} is the combined beam and pixel window function. This is an unbiased estimate of the Bmode map away from the mask boundaries. The black line shows the fiducial lensing Bpower of roughly (5 μKarcmin)^{2} 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 Bmode map. The raw crossspectra actually provide only a filtered version of the Bmode power, since the Wienerfiltering 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 crossspectra by the ratio of the result calculated on the FFP10 simulation suite to the simulated Bmode power. In all cases, the data crosscorrelation amplitudes are consistent with the fiducial model within 1σ, with the lensing Bmode power detected at 13, 18, and 20σ for templates constructed from the lensing quadratic estimator, the CIB, and in combination, respectively.
Fig. 13. Estimates of the lensing Bmode power spectrum from crossspectra between the Bmode polarization map and various Bmode templates constructed using different lensing tracers, as measured on 58% of the sky: the TT+TE+EE quadratic estimator (blue); the CIB map at 545 GHz (orange); and their combination (green). The raw crossspectra are a filtered version of the Bmode spectrum and have been rescaled using the expected scaling computed from simulations assuming an isotropic filter. The dashed black line is the FFP10 fiducial lensing Bmode spectrum. 
To delens, the lensed Bmode template is simply subtracted from our Bmode 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:
Fig. 14. Difference between the delensed and original Bmode power spectrum for the SMICA CMB polarization maps. In the upper panel, the delensed Bmode map is obtained by template subtraction. The templates are constructed using the labelled lensing tracers, and the Wienerfiltered Emode map, as discussed in Sect. 3.4.1. The dashed lines show predictions obtained by repeating these operations on the FFP10 set of Planck simulations. The black curve, for our fiducial model, shows the difference expected for perfect delensing. Summary statistics are presented in Table 3. In the lower panel, a remapping method is used, and the delensing signature is obtained after subtraction of biases, as described in Sect. 3.4.2, with summary statistics in Table 4. In this case, the actual total Bmode power is not decreased after delensing. 
The spectra are obtained from the positionspace Bmode maps before and after the Bmode template subtraction, using (scalar) deconvolution of the pseudoC_{ℓ} spectra to account for the effect of the mask. Table 3 lists the results of fitting a single delensing efficiency parameter across all the bins shown in Fig. 14. The efficiency defined in this way reaches 12% for the combined tracers. Largescale 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 (+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.
Reduction in the BB power spectrum after lensing Bmode template subtraction, for combinations of different lensing tracers.
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 Bmode 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 Bmode map, at negligible cost in tracer fidelity.
3.4.2. Acoustic peak sharpening (desmoothing)
We now test the expected acoustic peak resharpening 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 BBdelensing 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 nonoverlapping set of modes (Sehgal et al. 2017), at the cost of some loss of S/N, or a much more elaborate realizationdependent higherpoint 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 higherpoint 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 . The CMB maps are built from the inversevariance, homogeneouslyfiltered 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 ); 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 pseudoC_{ℓ} 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:
The combination 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, 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 still contains the difference of the noise spectra caused by the signal part of the tracer. To isolate the CMB peaksharpening, we therefore build the combination
with
and
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 χ^{2}minimization of the observed binned spectral differences
Fig. 15. Acoustic peak resharpening after delensing of the SMICA CMB maps. From top to bottom: difference between the delensed and original TT, TE, and EE spectra. Data points show delensing with our internal, minimumvariance, quadratic lensing estimate (blue), with the GNILC CIB map (orange), and with their optimal combination (green), as described in the text. For reference, the black line shows half the difference of the unlensed to lensed CMB spectra in the fiducial cosmology. In all cases, the solid coloured lines show the FFP10 simulation predictions, obtained after performing the same operations on each simulation as for the data. The dashed lines show the noise delensing bias defined by Eq. (55) (only relevant here for EE at high multipoles), while the dotdashed lines show the CMB peak smoothing caused by the tracer reconstruction noise, Eq. (56). 
together with predictions obtained from our simulations and the bestfit reduced χ^{2}. Here, is the binned difference of the unlensed and lensed power spectra in the fiducial model, and is the variance of the binned 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 Bmode 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 Bmode 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 Bmode delensing, remapping by the reconstruction noise producesB modes, and remapping the instrument noise also changes the Bmode 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 templatedelensing method, which are based on the actual change in Bmode power after delensing. The template method is also optimized to account for instrument noise (through the Wienerfiltering of the Emode map in the construction of the template). The templatebased efficiencies are, therefore, the relevant ones for assessing improvements in primordial Bmode limits from delensing. However, since the PlanckBmode measurement is noise dominated, there would be a negligible improvement in primordial Bmode 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 crosshalfmission 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 nonGaussianity of the lensing field, and Sect. 4.8 presents alternative calculations of the N^{(1)} bias.
4.1. Bandpower 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 temperatureonly 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 to the Planck TT,TE,EE+lowE+lensing bestfit cosmology and the lower panel the curl spectrum, both on the aggressive multipole range.
Fig. 16. Top: fractional differences between the lensing power spectrum band powers from the ϕ^{TT} reconstruction over the aggressive multipole range and the Planck 2018 bestfit cosmology using Planck TT,TE,EE+lowE+lensing (blue points). Bottom: temperatureonly curl reconstruction bandpower amplitudes. In both cases, the orange points show the same data with a coarser binning, where the construction of the binning function and midpoints follows the same procedure as the blue points (see Sect. 2.1). The curl bin at 264 ≤ L ≤ 901 deviates from zero with a formal significance of 4.3σ. This bin was chosen to maximize the deviation; after simple consideration of “lookelsewhere” effects as described in the text, we evaluate the PTE of that deviation to be 0.41% (or roughly 2.9 Gaussian σ), still suggestive of a problem in the curl reconstruction. 
The dip in the measured 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 bestfit spectrum, compared to 3.8σ in 2015 for the same bestfit spectrum. This significance is unchanged after marginalization of the CMB spectra uncertainties. Despite the large errors at highL, 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 crossspectrum 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 postBorn corrections that can be safely neglected at Planck noise levels (Pratten & Lewis 2016; Fabbian et al. 2018). Our curl amplitudes are built using the same binning procedure as for the lensing gradient in Eq. (10), using a reference flat spectrum
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 custommade 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 inversevariance weighted these points to produce the largest outlier possible^{14}. 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 temperatureonly 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 pointsource and MonteCarlo 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 colourcoded 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 highL range 401 ≤ L ≤ 2048. Reading the lower triangle from the left, the reconstruction labelled on the lefthand 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).
Fig. 17. Summary of lensingamplitudemeasurement consistency tests using temperatureonly reconstructions. Upper panel: difference between the measured lensing amplitude from pairs of reconstructions in units of the standard deviation of the difference from simulations. This tests whether the observed shifts are compatible with the behaviour expected from the FFP10 simulations. The lower triangle shows the amplitude measured from the conservative range 8 ≤ L ≤ 400, with a red colour indicating a larger amplitude of the reconstruction labelled on the left; the upper triangle uses the highL range (401 ≤ L ≤ 2048), with red indicating a larger amplitude of the reconstruction labelled at the top. Numerical significances are quoted for those tests that show differences greater than twice the expected standard deviation. The first entry is our baseline temperatureonly reconstruction; see the text for a description of all the other reconstructions. Lower panel: actual amplitudes (relative to our FFP10 fiducial model), with the result using the conservative range in blue and the highL range in orange. The dashed lines show for reference the values in our baseline reconstruction. 
Fig. 18. Same as Fig. 17, but for the curl reconstruction. Here the lower panel shows the curl amplitude over the conservative multipole range (8 ≤ L ≤ 400; blue) and over the range covering the negative curl feature shown in Fig. 16 (264 ≤ L ≤ 901; orange). We use the same binning scheme as for the lensing gradient (Eq. (10)), with an amplitude of unity corresponding to a flat curl spectrum. 
Componentseparation methods. In addition to using our baseline SMICA map, we perform reconstructions on the CMB maps from the three other componentseparation 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, pointsource mask, and SZclustertargeted 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 southeast 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 SZunm. reconstruction uses our baseline mask but without masking SZ clusters, and SZcons. includes a more conservative SZcluster mask. This conservative cluster mask was constructed by extending our SZcluster mask to at least three times the virial radius (θ_{500}) of each cluster, and enlarging the cluster catalogue: instead of using only the internallydetected Planck SZ catalogue, we added Xray clusters listed in the MetaCatalog of Xray 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.
Fig. 19. Various masks used for the lensing analysis and consistency tests. The topleft mask is the default lensing mask, including cutting out point sources, resolved SZ clusters, and the Galaxy; the topmiddle mask shows the difference between the new mask and the one used for the 2015 analysis. The remaining masks are used for consistency tests, as described in Sect. 4.2. The dust masks are constructed by thresholding the smoothed pixel variance obtained from the product of the two halfmission 545GHz Planck maps after filtering to include only multipoles 1000 ≤ ℓ ≤ 2000. 
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^{(0)} 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 Wienerfiltered 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 to real space, compute pseudoC_{ℓ} 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 Wienerfiltered 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).
Biashardened estimators. This block tests a series of temperatureonly lensing estimators orthogonalized to potential contaminants, called “biashardened” 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 . We can then define a new lensing estimate (before meanfield subtraction and normalization)
which has vanishing response to s, at the cost of an increase in variance. Here, the response functions and describe the responses to s of the lensing estimator and source estimator , respectively. We give results bias hardened against pointsource (S^{2}) contamination (Osborne et al. 2014; Shard.), where a spatiallyvarying point source signal is sought in the CMB data covariance, . The corresponding unnormalized optimal estimator is
We also harden against anisotropy in the variance of the instrumental noise (Nhard.). In this case the anisotropy is parametrized as , with quadratic esimator
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 (for example as a check against an effective spatial variation in calibration), which we denote by Mhard. The anisotropy is of the form , with resulting estimator
We have additionally attempted hardening against a residual dust amplitude (Dhard.). The motivation in this case is to look for a dust component consisting of a statisticallyisotropic field with spectral shape modulated by a varying amplitude :
The power spectrum shape 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
We also applied this hardened estimator to the 217GHz channel (217 GHz Dhard.). 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 217GHz 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 857GHz 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 143GHz channel with more aggressive SZcluster masking (143 GHz SZcons.) or no SZcluster masking at all (143 GHz SZunm.).
The most striking inconsistencies appearing in Fig. 17 are from the frequency reconstructions. These reconstructions use a simple templatecleaning method that is cruder than that used for the main componentseparated products. On the conservative multipole range only the 217 GHz reconstruction has possible anomalies, but on the highL range, both 143 and 217GHz 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 dustmasked 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 componentseparation methods show good consistency over the conservative range. However, we do see some anomalous behaviour at highL. Notably, we find
formally a 3σ deviation from zero. The actual shift in the highL amplitude between these two reconstructions is not that large however: about 0.5σ of our highL baseline measurement. None of the componentseparation 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 bandpower features even slightly. SMICA noSZ relies more on the 217GHz 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 tSZcluster mask, the conservative tSZcluster mask, or no cluster mask at all, both from the 143GHz channel or on the SMICA maps. At 143GHz, we do however detect and subtract a much larger pointsource correction without the SZcluster mask. This is consistent with the predictions of further tSZdedicated tests provided in Sect. 4.5.
The 143GHz channel shows some consistency issues on the highL range significant at around the 2.5σ level, where it displays consistently lower amplitudes than all other reconstructions. This is true independently of skycuts and biashardening 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 temperatureonly 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 pointsource 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 smallscale 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 highL range, the only potentially problematic test is ℓ_{min} = 1451, which gives a lower curl. This reconstruction is noisy and only weakly constraining on the highL range, and shows no suspicious offset on the conservative Lrange. 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 RDN^{(0)} 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 temperaturepolarization 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^{15} 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 crossspectra 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 crossspectrum shows very good agreement with the baseline reconstruction.
Fig. 20. Lensing reconstruction band powers over the full aggressive multipole range for the minimumvariance estimator (MV), and separately for all auto and crossspectra of the five quadratic lensing estimators that can be formed from the temperature (T) and polarization (E and B) SMICA CMB maps. The black line shows the Planck TT,TE,EE+lowE+lensing bestfit spectrum. In each lower left corner the χ^{2} probability to exceed (PTE) is given for the reconstruction (evaluated at the fiducial lensing spectrum), calculated across our conservative multipole range. 
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^{(1)} subtraction, we find
The feature is specific to the temperatureonly reconstruction; the next most constraining estimator over this multipole range is TT × EB, with
The polarizationonly estimator is much noisier but does not show anomalies on this range either:
The amplitude of the curl feature is reduced by almost a factor of 2 in the minimumvariance reconstruction compared to the temperatureonly reconstruction. This reduction is larger than expected from the FFP10 simulation behaviour. We find
a shift that is unexpected at the 3.5σ level (neglecting lookelsewhere 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 halfmission maps (HM1 and HM2) and their difference. First, we perform gradient and curl lensing reconstructions on their halfdifference, 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_{1}−HM_{2})^{4} (rescaled by a factor of 20 for visibility). The band powers are consistent with zero to within 2σ, as expected, with overall amplitude
Fig. 22. Halfmission lensing reconstructions tests for the gradient (ϕ; upper) and curl (Ω; lower). Blue points show reconstruction power spectra from SMICA halfmission difference maps, after multiplication by a factor of 20. The orange points show the autospectrum of the reconstruction built from HM_{1} and HM_{2} (with no noise mean field), and the green points the crossspectrum of the reconstruction built from HM_{1} with that built from HM_{2} (with no noise contribution to the reconstruction noise N^{(0)}). The black line in the upper plot shows the Planck TT,TE,EE+lowE+lensing bestfit spectrum. 
Second, we compare our baseline reconstruction to that obtained from halfmission maps. We consider two types of reconstruction, labelled (HM_{1}×HM_{2})^{2} and . In the first reconstruction, we estimate the lensing spectrum using a single lensing map, built using the first halfmission and second halfmission 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 crosscorrelate the lensing map built using the first halfmission data to another lensing map built using the second halfmission 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 halfmission FFP10 CMB simulations accompanying the 2018 Planck release (only noise simulations). Instead, we use FFP9 halfmission CMB simulations built for the 2015 release, coadding 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 beamprocessed CMB simulations also has a slight impact on the RDN^{(0)} bias. Hence, these halfmission reconstructions should not be considered as reliable as our baseline ones. The gradient and curl halfmission 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_{1}×HM_{2})^{2}, 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
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
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 halfmission map as we now discuss. Taking amplitudes with respect to the Planck bestfit Planck TT,TE,EE+lowE+lensing, and marginalizing over CMB uncertainties in order to quantify better the discrepancy, we have
The (HM_{1})^{4} result is formally a 4.5σ deviation from the bestfit lensing spectrum value. The remaining band powers of the HM_{1} reconstruction show no inconsistencies with the bestfit, 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 halfmission spectra are still low on these multipole ranges for both curl and gradient, but less than the temperatureonly spectra. We find
compared to −0.039 ± 0.015 for our baseline, and
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 nonGaussian 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 pointsource 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 smallscale polarization maps. Our use of foregroundcleaned 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. Semianalytic 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 nonGaussian 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 matterCMB 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^{−1}. The lensing potential is assumed to be Gaussian, but the matter density fields in each shell are simulated as lognormal fields; the covariance from class is therefore mapped into the covariance of the log fields, so that a Gaussian realization of can be exponentiated to obtain lognormal 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): protospheroid, spiral, and starburst. The flux density in each population is randomly distributed according to the redshiftdependent number counts from Planck ERCSC (Negrello et al. 2013), JCMT/SCUBA2 (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 populationdependent 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 217GHz 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^{4} in the following); and a larger contribution from the bispectrum (CIB^{2} × κ). 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 crosscorrelate the lensing estimator as applied on the CIB map to the input lensing potential. The resulting bias is twice this crossspectrum. Figure 23 shows the biases that are expected. For the uncleaned 217GHz 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 subpercent level, in the case of the tSZdeprojected 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.
Fig. 23. Estimates of the expected CIBinduced biases (as a fraction of the lensing potential power spectrum) to the temperatureonly lensing spectrum reconstruction at 217 GHz (without any cleaning), and for the SMICA frequency weighting. The largescalestructure bispectrum causes the quadratic estimator applied to the CIB map to correlate with the lensing signal, sourcing a 2% negative bias at 217 GHz (blue) on our conservative multipole range. This is reduced by an order of magnitude after SMICA cleaning (orange). The bias remains subpercent in the tSZdeprojected SMICA weighting (green). The CIB trispectrum sources a smaller contribution (red and purple at 217 GHz and for SMICA frequency weighting, respectively). Biases to the curl reconstruction are negligible. Error bars on this figure are estimated from the scatter of the unbinned power. 
We now turn to the thermal SZ component, focussing on the 143GHz channel since it is small in the 217GHz 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 pointsource correction procedure, even in the absence of masking, as well as the expected robustness of our bandpowers. 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_{500} 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 pointsource signal remains in the map, which would be detected by the Planck sourcedetection 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 highL 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 fullsky reference simulation includes instrumental noise in the form of Gaussian noise, independent between pixels but with amplitude in each pixel following the 143GHz channel variance map. Owing to the larger sky area, the blue error bars are slightly tighter than those on our 143GHz channel data reconstruction. The tSZ signal appears largely as a pointsourcelike contamination, and we distinguish results before (orange) and after (green) pointsource 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 pointsource signals, and we find a negligible impact, even without any masking.
Fig. 24. Impact of the thermal SZ (tSZ) on the temperatureonly lensing reconstruction spectrum at 143 GHz. We show the fractional difference between several reconstruction power spectra based on a 143 GHz fullsky simulation and the input lensing potential power spectrum (where the simulation includes tSZ and CMB, and Gaussian noise generated according to the 143GHz channel pixel variance map). Without any masking, the contamination appears predominantly as a large pointsource trispectrum contribution (orange). The pointsource correction is quite effective at reducing this term below L ≃ 400, with visible residual contamination at higher L (green). After applying a f_{sky} = 0.997 cluster mask (red), the differences from the power reconstructed in the reference CMBonly simulation with no masking (blue) are further reduced to a small fraction of a standard deviation in all bin. The left panel shows the conservative multipole range and the right panel the higher lensing multipoles. 
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 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 nonflat ΛCDM model bestfit 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 isotropicbeam 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, RDN^{(0)} and MC corrections, and via the likelihood in the perturbative correction functions and fiducial value of N^{(1)}. The simpler set of simulations we use here lacks the detailed beam shape and scanning model that affects the signal at lowL. 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^{(1)} contribution on the cut sky is accurately modeled by fullsky 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
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
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^{(1)} bias subtraction, where the lensing power spectrum is up to 20% higher in the Ω_{K} ≠ 0 model than the flat model (see Fig. 25).
Fig. 25. Comparison of reconstructions over the same highL multipole range on the same simulated map, using different fiducial model assumptions. The blue points use a fiducial model, set of simulations, and estimator ingredients based on the FFP10 input cosmology. The orange points use a cosmology with Ω_{K} ≃ −0.04, where the lensing power spectrum is higher by 10–20% (black line), and the CMB spectra differ by roughly 2% at high ℓ. The linear corrections to the likelihood account for the differences in fiducial model, and are mostly important at highL in this case; they are effective (green) at bringing the reconstructions with the different fiducial models into good agreement. 
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 nonGaussian at some level due to nonlinear growth of largescale structure (LSS) and postBorn lensing (Pratten & Lewis 2016). A nonvanishing LSS bispectrum will source contributions to the lensing reconstruction spectrum involving three powers of ϕ, potentially giving rise to an additional lensing bias as well as the N^{(0)} and N^{(1)} 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 treelevel perturbation theory for the LSS bispectrum and neglecting postBorn 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 postBorn effects reduces the signal even further, as expected from the opposite sign of postBorn 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 crosscorrelations between the lensing spectrum and lowredshift largescale structure tracers could, however, be larger due to larger lowredshift LSS nonGaussianity, and smaller oppositesign postBorn contributions at low redshift. However, for Planck noise levels, following Fabbian et al. (2019) we estimate that the crosscorrelation bias remains below 1% for tracers at z ≳ 0.2 where there is a useful crosscorrelation 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 nonGaussian, since it is a nonlinear function of the maps, for example the reconstruction 1point function is skewed (Liu et al. 2016). However, this nonGaussianity 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^{(1)} lensing bias
The N^{(1)} 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^{(1)} 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^{(1)}subtracted band powers without assuming a fiducial lensing spectrum. We then test a simulationbased N^{(1)} estimator, assessing the impact of sky cuts and sky curvature; this also forms our baseline N^{(1)} calculation in the case of inhomogeneous filtering.
The N^{(1)} correction captures the estimator’s linear response to C^{ϕϕ} from nonprimary contractions. As such, we can write (neglecting pointsource and MC corrections, and complications due to masking in the following discussion)
Instead of subtracting the rightmost term using a fiducial , we can invert this matrix relation to obtain the lensing spectrum. This can also be generalized to calculate the curl N^{(1)}, 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 Yinduced N^{(1)} contribution to the Xspectrum estimate, we can write
In this equation represents a hypothetical curllike deflection caused, for example, by an instrumental systematic effect that we assume is uncorrelated with the lensing signal. Subpixel 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 “xpRDN^{(1)}” in Fig. 26, while inverting Eq. (78) provides yet another set of lensing band powers, labelled “RDN^{(1)}.” The results for the temperatureonly 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^{(1)} subtraction. For comparison, the red curves display the case of no N^{(1)} subtraction, and the black dashed curves the N^{(1)} 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 flatsky isotropic analytic form given in Appendix A. These different methods for dealing with the N^{(1)} 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^{(1)} subtraction method does not substantially affect our curl nulltest results either. A deconvolution of the N^{(1)} 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.
Fig. 26. Comparison of reconstruction band powers for various methods to remove the N^{(1)} lensing bias. Our baseline band powers are built by subtracting the bias calculated analytically with the FFP10 fiducial model, and are shown as the blue points with error bars. The curves labelled “RDN^{(1)}” and “xpRDN^{(1)}” invert Eqs. (78) and (79), respectively, and make no assumption about the theory lensing spectrum. The curves labelled “MCN^{(1)}” subtract a simulationbased estimate of N^{(1)}. For comparison the red lines show the band powers without any N^{(1)} subtraction, and the dashed black lines the N^{(1)} bias itself. The top plot shows the lensing reconstruction from temperature only, and is normalized to the FFP10 fiducial . The lower plot shows the temperatureonly curl reconstruction (which on small scales appears to be coincidentally close to the result one would expect with no N^{(1)} at all). 
We also test a simulationbased N^{(1)} estimation, labelled “MCN^{(1)}” 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^{(1)} calculations). We use simplified, effectively isotropic CMB simulations that do not contain the full anisotropic Planck beam model or SMICA processing, since a full resimulation would be too numerically expensive. The agreement between the simulationbased N^{(1)} estimator and the other methods is good, to the point that it does not visibly affect the band powers.
We use the simulationbased MCN^{(1)} as the baseline for our results using inhomogeneous filtering. Figure 27 shows MCN^{(1)} (blue points) for the polarization reconstruction using the inhomogeneouslyfiltered maps. The large dynamic range of the noisevariance map used for filtering makes the naive analytic prediction (the orange curve) fail to reproduce the MCN^{(1)} shape. However, it is straightforward to understand the simulation result analytically by using the same toy model discussed in Sect. 2.3. Let be the analytic (fullsky) N^{(1)} bias calculated using the noise levels at angular position , as specified by the variance map used in the filtering. The green curve in Fig. 27 shows the weighted average
Fig. 27. Comparison of the simulationbased N^{(1)} lensing bias estimates (blue points) to analytic estimates in the case of our inhomogeneouslyfiltered polarization reconstruction. The orange curve shows the naive prediction using a single mapaveraged noise level, which does not capture the variation of the N^{(1)} bias across the sky. The green curve shows the weighted average Eq. (80). The black line shows the fiducial lensing power spectrum. 
This is evaluated by splitting the sky into 15 patches with roughly constant noise levels, and successfully reproduces the simulationbased N^{(1)}.
5. Data products
The final Planck lensing data products are available on the Planck Legacy Archive^{16}, 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 SZcluster 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 tSZdeprojected 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 bestestimate lensing potential, together with the corresponding simulation suite.
– lensing Bmode templates and simulations built from our best Emode 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 percentlevel estimates of , 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 σ_{8} substantially more powerful due to the weaker degeneracy with Ω_{m}. In combination with lowerredshift tracers, degeneracies can be further broken, and CMB lensing provides a powerful highredshift 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 peaksharpening and a reduction in Bmode power in line with expectations. For future Bmode 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 crosscorrelation studies, where the detailed modeling of the autospectrum 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 polarizationonly 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 SZcluster masking. No planned experiments will be able to provide comparable quality maps over the full sky for many years. Ongoing and forthcoming groundbased 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 largestscale lensing modes for the immediate future.
Planck (https://www.esa.int/Planck) is a project of the European Space Agency (ESA) with instruments provided by two scientific consortia funded by ESA member states and led by Principal Investigators from France and Italy, telescope reflectors provided through a collaboration between ESA and a scientific consortium led and funded by Denmark, and additional contributions from NASA (USA).
The SMICA mask was a preliminary mask constructed for the SMICA 2018 analysis; it differs from the final 2018 componentseparation mask described in Planck Collaboration IV (2020), since this was finalized later. We make the mask used for the lensing analysis available with the other lensing products, and show results using the final componentseparation mask in Table 2 for comparison.
This likelihood combination corresponds to that denoted Planck TT,TE,EE+lowE+lensing in Planck Collaboration VI (2020); this is the baseline combination advocated there for parameter constraints.
Planck Collaboration XVII (2014) lists in Appendix C several arguments and tests (performed with more simulations than we are using in this paper) that justifies our use of a Gaussian likelihood. These tests performed on the updated simulations do not show any qualitative difference.
From a set of seven quasar absorptionline observations, Cooke et al. (2018) estimate a primordial deuterium ratio 10^{5}D/H = 2.527 ± 0.030. Assuming that standard BBN can be solved exactly, the D/H measurement can be converted into an Ω_{b}h^{2} measurement with notional 1σ statistical error of 1.6 × 10^{−4}. However, as discussed in PCP18, the central value depends on various nuclear rate parameters that are uncertain at this level of accuracy. For example, adopting the theoretical rate of Marcucci et al. (2016), rather than the defaults in the PArthENoPE code (Pisanti et al. 2008), results in a central value shifted to Ω_{b}h^{2} = 0.02198 compared to Ω_{b}h^{2} = 0.02270, while Cooke et al. (2018) quote a central value of Ω_{b}h^{2} = 0.02166 using Marcucci et al. (2016) but a different BBN code. Our conservative BBN prior is centred at the midpoint of these two differences, with error bar increased so that the different results (and other rate uncertainties) lie within approximately 1σ of each other.
Chains at https://pla.esac.esa.int, description and parameter tables in Planck Collaboration (2018).
In particular, our flat parameter prior 0.4 < h < 1 on the Hubble constant 100h km s^{−1} Mpc^{−1} is wider than the 0.55 < h < 0.91 range assumed by Troxel et al. (2018), and our flat prior on Ω_{c}h^{2} shifts some of the probability weight with respect to the flat prior on Ω_{m} used by DES. We also approximate the contribution of massive neutrinos as a single mass eigenstate with m_{ν} = 0.06 eV in our baseΛCDM model, rather than marginalizing over the sum of the neutrino masses with a flat prior, and use HMcode (Mead et al. 2016) for nonlinear corrections.
For fixed neutrino mass, the DES joint constraint using the DES priors described in the caption of Table 2 is σ_{8}(Ω_{m}/0.3)^{0.5} = 0.793 ± 0.024. Marginalizing over neutrino mass, the joint DESPlanck lensingonly result using DES priors is .
Due to the lack of a sufficiently precise realizationdependent debiaser for all simulations, for this analysis we have used the covariance matrix built with a realizationindependent Monte Carlo N^{(0)} (defined similarly to the gradientmode MCN^{(0)} in Appendix A) consistently on simulations and data.
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/20072013)/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/INSUIN2P3INP (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/planckcollaboration.
References
 Alam, S., Ata, M., Bailey, S., et al. 2017, MNRAS, 470, 2617 [NASA ADS] [CrossRef] [Google Scholar]
 Arnaud, M., Pointecouteau, E., & Pratt, G. W. 2005, A&A, 441, 893 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Arnaud, M., Pratt, G. W., Piffaretti, R., et al. 2011, A&A, 517, A92 [Google Scholar]
 Aubourg, É., Bailey, S., Bautista, J. E., et al. 2015, Phys. Rev. D, 92, 123516 [NASA ADS] [CrossRef] [Google Scholar]
 Beck, D., Fabbian, G., & Errard, J. 2018, Phys. Rev. D, 98, 043512 [CrossRef] [Google Scholar]
 Béthermin, M., Wang, L., Doré, O., et al. 2013, A&A, 557, A66 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Beutler, F., Blake, C., Colless, M., et al. 2011, MNRAS, 416, 3017 [NASA ADS] [CrossRef] [Google Scholar]
 Blanchard, A., & Schneider, J. 1987, A&A, 184, 1 [NASA ADS] [Google Scholar]
 Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 7, 034 [NASA ADS] [CrossRef] [Google Scholar]
 Böhm, V., Schmittfull, M., & Sherwin, B. D. 2016, Phys. Rev. D, 94, 043519 [CrossRef] [Google Scholar]
 Böhm, V., Sherwin, B. D., Liu, J., et al. 2018, Phys. Rev. D, 98, 123510 [CrossRef] [Google Scholar]
 Carron, J., & Lewis, A. 2017, Phys. Rev. D, 96, 063510 [CrossRef] [Google Scholar]
 Carron, J., Lewis, A., & Challinor, A. 2017, JCAP, 1705, 035 [Google Scholar]
 Challinor, A., & Chon, G. 2002, Phys. Rev. D, 66, 127301 [CrossRef] [Google Scholar]
 Chen, C.C., Cowie, L. L., Barger, A. J., et al. 2013, ApJ, 762, 81 [NASA ADS] [CrossRef] [Google Scholar]
 Cooke, R. J., Pettini, M., & Steidel, C. C. 2018, ApJ, 855, 102 [NASA ADS] [CrossRef] [Google Scholar]
 Das, S., Errard, J., & Spergel, D. 2013, ArXiv eprints [arXiv:1311.2338] [Google Scholar]
 Delabrouille, J., Betoule, M., Melin, J.B., et al. 2013, A&A, 553, A96 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 DES Collaboration 2018a, MNRAS, 480, 3879 [NASA ADS] [CrossRef] [Google Scholar]
 DES Collaboration 2018b, Phys. Rev. D, 98, 043526 [Google Scholar]
 Fabbian, G., Calabrese, M., & Carbone, C. 2018, JCAP, 1802, 050 [CrossRef] [Google Scholar]
 Fabbian, G., Lewis, A., & Beck, D. 2019, JCAP, 10, 057 [CrossRef] [Google Scholar]
 Fan, X.H., Carilli, C. L., & Keating, B. 2006, ARA&A, 44, 415 [NASA ADS] [CrossRef] [Google Scholar]
 Ferraro, S., & Hill, J. C. 2018, Phys. Rev. D, 97, 023512 [CrossRef] [Google Scholar]
 Geach, J. E., & Peacock, J. A. 2017, Nat. Astron., 1, 795 [NASA ADS] [CrossRef] [Google Scholar]
 Gorski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
 Green, D., Meyers, J., & van Engelen, A. 2017, JCAP, 1712, 005 [CrossRef] [Google Scholar]
 Greiner, M., & Enßlin, T. A. 2015, A&A, 574, A86 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hand, N., Leauthaud, A., Das, S., et al. 2015, Phys. Rev. D, 91, 062001 [CrossRef] [Google Scholar]
 Hanson, D., Challinor, A., Efstathiou, G., & Bielewicz, P. 2011, Phys. Rev. D, 83, 043005 [NASA ADS] [CrossRef] [Google Scholar]
 Hanson, D., Hoover, S., Crites, A., et al. 2013, Phys. Rev. Lett., 111, 141301 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 HarnoisDéraps, J., Tröster, T., Chisari, N. E., et al. 2017, MNRAS, 471, 1619 [NASA ADS] [CrossRef] [Google Scholar]
 Hartlap, J., Simon, P., & Schneider, P. 2007, A&A, 464, 399 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hirata, C. M., & Seljak, U. 2003, Phys. Rev. D, 68, 083002 [NASA ADS] [CrossRef] [Google Scholar]
 Hu, W., & Okamoto, T. 2002, ApJ, 574, 566 [Google Scholar]
 Kesden, M., Cooray, A., & Kamionkowski, M. 2002, Phys. Rev. Lett., 89, 011304 [NASA ADS] [CrossRef] [Google Scholar]
 Kesden, M., Cooray, A., & Kamionkowski, M. 2003, Phys. Rev. D, 67, 123507 [NASA ADS] [CrossRef] [Google Scholar]
 Kirk, D., Omori, Y., BenoitLévy, A., et al. 2016, MNRAS, 459, 21 [CrossRef] [Google Scholar]
 Knox, L., & Song, Y.S. 2002, Phys. Rev. Lett., 89, 011303 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 Krause, E., Eifler, T. F., Zuntz, J., et al. 2017, ArXiv eprints [arXiv:1706.09359] [Google Scholar]
 Larsen, P., & Challinor, A. 2016, MNRAS, 461, 4343 [CrossRef] [Google Scholar]
 Larsen, P., Challinor, A., Sherwin, B. D., & Mak, D. 2016, Phys. Rev. Lett., 117, 151102 [CrossRef] [Google Scholar]
 Lenz, D., Doré, O., & Lagache, G. 2019, ApJ, 883, 75 [CrossRef] [Google Scholar]
 Lewis, A. 2005, Phys. Rev. D, 71, 083008 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A. 2013, Phys. Rev. D, 87, 103529 [NASA ADS] [CrossRef] [Google Scholar]
 Lewis, A. 2019, GetDist: a Python Package for Analysing Monte Carlo Samples, https://getdist.readthedocs.io [Google Scholar]
 Lewis, A., & Challinor, A. 2006, Phys. Rep., 429, 1 [Google Scholar]
 Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473 [Google Scholar]
 Lewis, A., Challinor, A., & Hanson, D. 2011, JCAP, 1103, 018 [NASA ADS] [CrossRef] [Google Scholar]
 Liu, J., & Hill, J. C. 2015, Phys. Rev. D, 92, 063517 [CrossRef] [Google Scholar]
 Liu, J., Hill, J. C., Sherwin, B. D., et al. 2016, Phys. Rev. D, 94, 103501 [CrossRef] [Google Scholar]
 Mak, D. S. Y., Challinor, A., Efstathiou, G., Lagache, G., & Lagache, G. 2017, MNRAS, 466, 286 [NASA ADS] [CrossRef] [Google Scholar]
 Manzotti, A., Story, K. T., Wu, W. L. K., et al. 2017, ApJ, 846, 45 [CrossRef] [Google Scholar]
 Marcucci, L. E., Mangano, G., Kievsky, A., & Viviani, M. 2016, Phys. Rev. Lett., 116, 102501; erratum: Phys. Rev. Lett., 117, 049901 (2016) [CrossRef] [PubMed] [Google Scholar]
 Mead, A., Heymans, C., Lombriser, L., et al. 2016, MNRAS, 459, 1468 [Google Scholar]
 Namikawa, T. 2017, Phys. Rev. D, 95, 103514 [CrossRef] [Google Scholar]
 Namikawa, T., Hanson, D., & Takahashi, R. 2013, MNRAS, 431, 609 [NASA ADS] [CrossRef] [Google Scholar]
 Negrello, M., Clemens, M., GonzalezNuevo, J., et al. 2013, MNRAS, 429, 1309 [NASA ADS] [CrossRef] [Google Scholar]
 Okamoto, T., & Hu, W. 2003, Phys. Rev. D, 67, 083002 [NASA ADS] [CrossRef] [Google Scholar]
 Osborne, S. J., Hanson, D., & Doré, O. 2014, JCAP, 1403, 024 [NASA ADS] [CrossRef] [Google Scholar]
 Peebles, P. J. E. 1993, Principles of Physical Cosmology (Princeton: Princeton University Press) [Google Scholar]
 Peloton, J., Schmittfull, M., Lewis, A., Carron, J., & Zahn, O. 2017, Phys. Rev. D, 95, 043508 [CrossRef] [Google Scholar]
 Piffaretti, R., Arnaud, M., Pratt, G. W., Pointecouteau, E., & Melin, J. B. 2011, A&A, 534, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pisanti, O., Cirillo, A., Esposito, S., et al. 2008, Comput. Phys. Commun., 178, 956 [NASA ADS] [CrossRef] [Google Scholar]
 Planck Collaboration. 2018, The Legacy Explanatory Supplement, http://wiki.cosmos.esa.int/plancklegacyarchive (ESA) [Google Scholar]
 Planck Collaboration XVIII. 2011, A&A, 536, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVII. 2014, A&A, 571, A17 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XVIII. 2014, A&A, 571, A18 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXVII. 2014, A&A, 571, A27 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XXX. 2014, A&A, 571, A30 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2016, A&A, 594, A11 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2016, A&A, 594, A12 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XIII. 2016, A&A, 594, A13 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XV. 2016, A&A, 594, A15 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration I. 2020, A&A, 641, A1 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration II. 2020, A&A, 641, A2 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration III. 2020, A&A, 641, A3 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IV. 2020, A&A, 641, A4 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration V. 2020, A&A, 641, A5 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VI. 2020, A&A, 641, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VII. 2020, A&A, 641, A7 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration VIII. 2020, A&A, 641, A8 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration IX. 2020, A&A, 641, A9 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration X. 2020, A&A, 641, A10 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XI. 2020, A&A, 641, A11 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration XII. 2020, A&A, 641, A12 [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVI. 2016, A&A, 596, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Planck Collaboration Int. XLVIII. 2016, A&A, 596, A109 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Pratten, G., & Lewis, A. 2016, JCAP, 1608, 047 [CrossRef] [Google Scholar]
 Riess, A. G., Casertano, S., Yuan, W., et al. 2018, ApJ, 855, 136 [NASA ADS] [CrossRef] [Google Scholar]
 Ross, A. J., Samushia, L., Howlett, C., et al. 2015, MNRAS, 449, 835 [NASA ADS] [CrossRef] [Google Scholar]
 Schaan, E., Krause, E., Eifler, T., et al. 2017, Phys. Rev. D, 95, 123512 [NASA ADS] [CrossRef] [Google Scholar]
 Schaan, E., Ferraro, S., & Spergel, D. N. 2018, Phys. Rev. D, 97, 123539 [CrossRef] [Google Scholar]
 Scott, K. S., Wilson, G. W., Aretxaga, I., et al. 2012, MNRAS, 423, 575 [NASA ADS] [CrossRef] [Google Scholar]
 Sehgal, N., Madhavacheril, M. S., Sherwin, B., & van Engelen, A. 2017, Phys. Rev. D, 95, 103512 [CrossRef] [Google Scholar]
 Sherwin, B. D., & Schmittfull, M. 2015, Phys. Rev. D, 92, 043005 [NASA ADS] [CrossRef] [Google Scholar]
 Sherwin, B. D., van Engelen, A., Sehgal, N., et al. 2017, Phys. Rev. D, 95, 123529 [NASA ADS] [CrossRef] [Google Scholar]
 Simard, G., Omori, Y., Aylor, K., et al. 2017, ApJ, 860, 137 [CrossRef] [Google Scholar]
 Singh, S., Mandelbaum, R., & Brownstein, J. R. 2017, MNRAS, 464, 2120 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, R. E., Peacock, J. A., Jenkins, A., et al. 2003, MNRAS, 341, 1311 [Google Scholar]
 Smith, K. M., Hanson, D., LoVerde, M., Hirata, C. M., & Zahn, O. 2012, JCAP, 1206, 014 [Google Scholar]
 Song, Y.S., Cooray, A., Knox, L., & Zaldarriaga, M. 2003, ApJ, 590, 664 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. M., Zahn, O., & Dore, O. 2007, Phys. Rev. D, 76, 043510 [NASA ADS] [CrossRef] [Google Scholar]
 Smith, K. M., Cooray, A., Das, S., et al. 2009, AIP Conf. Proc., 1141, 121 [NASA ADS] [CrossRef] [Google Scholar]
 Story, K. T., Hanson, D., Ade, P. A. R., et al. 2015, ApJ, 810, 50 [NASA ADS] [CrossRef] [Google Scholar]
 Takahashi, R., Sato, M., Nishimichi, T., Taruya, A., & Oguri, M. 2012, ApJ, 761, 152 [NASA ADS] [CrossRef] [Google Scholar]
 Teng, W.H., Kuo, C.L., & Wu, J.H. P. 2011, ArXiv eprints [arXiv:1102.5729] [Google Scholar]
 Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709 [NASA ADS] [CrossRef] [Google Scholar]
 Troxel, M. A., MacCrann, N., Zuntz, J., et al. 2018, Phys. Rev. D, 98, 043528 [NASA ADS] [CrossRef] [Google Scholar]
 Vallinotto, A. 2012, ApJ, 759, 32 [CrossRef] [Google Scholar]
 van Engelen, A., Bhattacharya, S., Sehgal, N., et al. 2014, ApJ, 786, 13 [NASA ADS] [CrossRef] [Google Scholar]
 Wandelt, B. D., Hivon, E., & Gorski, K. M. 2001, Phys. Rev. D, 64, 083003 [NASA ADS] [CrossRef] [Google Scholar]
 Yu, B., Hill, J. C., & Sherwin, B. D. 2017, Phys. Rev. D, 96, 123511 [CrossRef] [Google Scholar]
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.
Fig. A.1. Compilation of the power spectrum reconstruction biases for our baseline MV analysis. The blue and orange lines show RDN^{(0)} and the analytic N^{(1)} estimate. The green curve shows the pointsource correction, and the red points the additive Monte Carlo correction (shown here after a coarse binning was applied) defined by the difference of Eq. (A.2) to the simulation input fiducial lensing spectrum. This correction is most prominent at low lensing multipoles, where it is negative (dashed) and due to masking. The error bars are obtained from the 240 simulations used to calculate the correction. The black line shows the fiducial lensing power spectrum. The black dot at L = 1 shows the almost pure dipole effective deflection caused by our motion with respect to the CMB frame (aberration). This dipole is included in our simulations, hence subtracted from our lensing reconstructions together with other sources of anisotropy through the mean field. 
We introduce three types of lensing power spectrum estimators built using the simulations: is the power spectrum of the quadratic estimator with both legs^{17} being simulation i; is the analogous spectrum where one leg of the quadratic estimator is always the data, and the second leg is simulation i; and is where simulation i is the first leg and simulation j the second leg and i ≠ j. We note that contains a meanfield 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 MCN^{(0)} and RDN^{(0)} are then defined as
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 (or ) 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 ; 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 (used to build the Monte Carlo correction, Eq. (10), after binning) is
In the last definition, the analytic expression for N^{(1)} is used, with the exception of our inhomogeneouslyfiltered reconstructions, for which a Monte Carlo estimate MCN^{(1)} is used instead. In general, for a source of anisotropy s, the analytic sinduced N^{(1)} is calculated according to the flatsky approximation using (Kesden et al. 2003)
with , and the labels X, Y, I, J stand for the T, E, or B map. Further, is the analytic inversevariance filter 𝒯^{†}Cov^{−1}𝒯 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 flatsky estimator weights as applied to the inversevariance filtered X and Y fields, f^{XY, s} is the CMB XY flatsky covariance response to anisotropy source s, and 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 curlinduced and pointsource induced N^{(1)} biases for robustness tests.
The Monte Carlo MCN^{(1)} (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^{(1)} are much lower than on N^{(0)}, and Sect. 4.8 demonstrates that the impact of the CMB nonidealities on N^{(1)} are negligible). These simulations also neglect and , which are negligible on the relevant scales. We again form spectra mixing quadratic estimates and from these simulations, where is defined as above on the first pair member of these new simulations, and takes the quadratic estimate built from the first pair member on one leg and the second pair member on the second leg. Then
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 temperatureonly reconstruction. The lefthand panel focuses on the lowL range (L ≤ 400) and the righthand 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 fullsky reconstructions to remove the effect of masking, and used idealized statisticallyisotropic 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 meanfield spectra obtained from the full, anisotropic, and masked simulations to the mean field from the same simulations without any masking.
Fig. B.1. Characterization of the meanfield power spectrum and its components in our temperatureonly reconstruction. The lefthand part shows the convergencelike power spectrum on large scales, while the righthand part shows the deflection power spectrum on smaller scales. The blue line shows the total mean field as obtained from the SMICA FFP10 simulation set. Lines showing the contributions from statistical anisotropy of the mask (red), noise (orange), and the beamconvolved and pixelized CMB anisotropies (green) are plotted with rough 68% confidence regions shown shaded on the right panel (the uncertainty arises from the finite number of simulations used to estimate the mean fields). The mask mean field dominates at low multipoles and is barely distinguishable from the total on the lefthand plot. The orange and green curves were obtained from fullsky lensing reconstructions using idealized, statisticallyisotropic CMB or noise components, respectively. The red curve was obtained by differencing the meanfield spectra on the full anisotropic simulations, with and without masking. Pixelization effects due to subpixel pointing offsets are expected to appear as an approximately whitenoise lensing deflection component of amplitude 0.′05 (the brown line shows the prediction for the 217GHz channel), and is clearly detected by crosscorrelating the mean field to the subpixel deflection prediction (pink curve). The rise of the mean field at L ≃ 3000 originating in the noise maps is mostly a simulation artefact, sourced by nonlinearities in the data processing that cause slight correlations between the simulations (due to the fixed fiducial CMB and foregrounds used to make the noise simulations); the purple curve shows the lensing quadratic estimator applied to the empirical average of the noise simulations. The predicted mean field due to inhomogeneity of the variance of idealized uncorrelated pixel noise (grey curve, as derived from Eq. (B.1)) is much weaker on smaller scales and is only visible at low lensing multipoles in this figure. 
As in previous releases, the mean field is very large at low multipoles, where it is mainly due to masking, with the noise and beamanisotropy 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 pixelpixel correlations). The expected noise meanfield contribution is, from Eq. (3),
and hence, according to Eq. (6), is a pure gradient field with . Here, is the spin0 harmonic transform of the temperature noise variance map, and
Here, b_{ℓ}, , 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 μKarcmin). 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 timeordered 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 subpixel miscentring is the same in all simulations, it appears rather directly in the meanfield estimate (with a possible bias, since, unlike lensing, the subpixel deflections act on the sky signal after beam convolution). The prediction for the 217GHz channel subpixel deflection spectrum is shown as the purple curve in Fig. B.1, and has the approximate power spectrum of a whitenoise displacement component (as expected for a large number of independent hits) with an amplitude of around 0.′05. The crossspectrum of the meanfield estimate for the SMICA map with the predicted subpixel 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 crosscorrelation, 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 crosscorrelating the results, we attempt to isolate the meanfield 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^{−2} 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 meanfield estimate and the lensing N^{(0)} bias (see Appendix A). We use N_{bias} = 60 simulations for the mean field and N_{bias} = 240 independent simulations for N^{(0)}. This section describes how the impact of the Monte Carlo error is calculated. For simplicity we use the following form of the Monte Carlocorrected final power spectrum estimator, where the Monte Carlo correction is applied additively instead of multiplicatively:
(and neglecting the pointsource correction), obtained combining Eqs. (A.1) and (A.2). Here, is the lensing reconstruction power spectrum estimate from the data, including meanfield subtraction (i.e., it is equivalent to 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 , 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 ; and similarly 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 crossspectrum. Denoting the (small and independent) error on the estimated mean field for each quadratic estimate by , all relevant terms are of the form
where is either the data lensing map (for ) or that of simulation i ( for ). Summing the terms from and at fixed L and M, the contribution to the error on the spectrum is
We note that terms of the form cancel between and . We proceed to calculate the averaged squared error on the spectrum. The spectrum of the meanfield error is
neglecting the N^{(1)} 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
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 meanfield evaluation has variance:
C.2. N^{(0)} 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),
We proceed by approximating the reconstructions as Gaussian, in which case correlations between the errors in , , and vanish. Empirical verification confirms that to a good approximation, errors in these spectra are independent, with approximate Gaussian variance
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 in Eq. (C.9), which is accurate on almost all scales, we obtain
This gives Eq. (13) in the main text, after identifying the diagonal bandpower Gaussian variance with . 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, ; this agrees reasonably well with the analytic estimate (the correction is small, so it does not need to be calculated accurately).
Fig. C.1. Crude empirical estimates of the additional variance on the MV reconstruction band powers (on the conservative multipole range) due to the finite number of simulations used for meanfield and bias correction. The additional variance is normalized by the baseline statistical variance of the band powers. These empirical estimates were obtained by splitting the simulations into five independent subsets, building reconstruction band powers using each subset of the simulations, and measuring the scatter between the five sets of band powers. The results are scaled by 1/5 to approximate the full simulation set used in the main analysis. The horizontal line shows the simple constant relative correction applied to the covariance matrix, Eq. (13), as obtained in Appendix C. The error bars on the empirical variance (identical for each bin, after rescaling by ) are those expected in the fiducial model of Eq. (13), assuming Gaussian statistics of the lensing maps. 
All Tables
Lensingreconstruction powerspectrum bandpower amplitudes and errors for the temperatureonly (TT), combined temperatureandpolarization minimumvariance (MV), and polarization estimators (PP).
Parameter constraints for different lensing datasets and priors, with and without galaxy BAO data.
Reduction in the BB power spectrum after lensing Bmode template subtraction, for combinations of different lensing tracers.
All Figures
Fig. 1. Mollweide projection in Galactic coordinates of the lensingdeflection reconstruction map from our baseline minimumvariance (MV) analysis. We show the Wienerfiltered displacementlike scalar field with multipoles , corresponding to the gradient mode (or E mode) of the lensing deflection angle. Modes with L < 8 have been filtered out. 

In the text 
Fig. 2. Noisevariance maps (shown as noise rms in μKarcmin) that we use to filter the SMICA CMB maps that are fed into the quadratic estimators when performing inhomogeneous filtering. The upper panel shows the temperature noise map, with median over our unmasked sky area of 27 μKarcmin. We use a common noise map for Q and U polarization, also neglecting Q and U noise correlations, shown in the lower panel, which spans an entire order of magnitude, with median 52 μKarcmin (larger than times the temperature noise because not all the Planck detectors are polarized). In temperature, the variance map has a homogeneous (approximately 5 μKarcmin) contribution from the isotropic additional Gaussian power that we add to the simulations to account for residual foreground contamination. 

In the text 
Fig. 3. Monte Carloderived multiplicative normalization corrections for the polarization reconstruction, using inhomogeneous filtering (blue points). Band powers are divided by these numbers to provide our final estimates. This correction is not small, and is sourced by the large spatial variation of the estimator response; however, it is very well reproduced by the approximate analytic model of Eq. (21), shown as the dashed blue line. The correction for our baseline MV band powers using homogeneous filtering is shown as the orange points. 

In the text 
Fig. 4. Planck 2018 lensing reconstruction band powers (values and multipole ranges are listed in Table 1). Left: minimumvariance (MV) lensing band powers, shown here using the aggressive (blue, 8 ≤ L ≤ 2048) and conservative (orange, 8 ≤ L ≤ 400) multipole ranges. The dots show the weighted bin centres and the fiducial lensing power spectrum is shown as the black line. Right: comparison of polarizationonly band powers using homogeneous map filtering (blue boxes, with dots showing the weighted bin centres) and the more optimal inhomogeneous filtering (orange error bars). The inhomogeneous filtering gives a scaledependent increase in S/N, amounting to a reduction of 30% in the error on the amplitude of the power spectrum over the conservative multipole range shown. The black line is the fiducial lensing power spectrum. 

In the text 
Fig. 5. Planck 2018 lensing powerspectrum band powers (pink boxes) over the aggressive multipole range. The 2015 analysis band powers (green) were calculated assuming a slightly different fiducial model and have not been (linearly) corrected to the 2018 model. Also shown are recent measurements by the ACTPol (Sherwin et al. 2017), SPTpol (Story et al. 2015), and SPTSZ (Simard et al. 2017) collaborations. The SPTSZ measurement is not completely independent, since the SPTSZ reconstruction also uses temperature data from Planck, but with subdominant weight over the smaller sky area used. The black line shows the lensing potential power spectrum for the ΛCDM bestfit parameters to the Planck 2018 likelihoods (Planck TT,TE,EE+lowE, which excludes the lensing reconstruction). 

In the text 
Fig. 6. Constraints on Ω_{m} and σ_{8} in the baseΛCDM model from CMB lensing alone (posterior sample points coloured by the value of the Hubble constant in units of km s^{−1} Mpc^{−1}) using the priors described in the text. Grey bands give corresponding 1σ and 2σ lensingonly constraints using the approximate fit . The joint 68% and 95% constraints from CMB lensing with the addition of BAO data (Beutler et al. 2011; Ross et al. 2015; Alam et al. 2017) are shown as the dashed contours, and the constraint from the Planck CMB powerspectrum data is shown for comparison as the solid contours. 

In the text 
Fig. 7. Constraints on Ω_{m} and σ_{8} in the baseΛCDM model from Planck temperature and polarization power spectra (red), and the tighter combined constraint with CMB lensing (blue). The dashed line shows the joint result when the reionization redshift is restricted to z_{re} > 6.5 to be consistent with observations of highredshift quasars (Fan et al. 2006). Contours contain 68% and 95% of the probability. 

In the text 
Fig. 8. Constraints on Ω_{m} and σ_{8} in the base ΛCDM model from DES galaxy lensing (green), Planck CMB lensing (grey), and the joint constraint (red). The Planck powerspectrum constraint is shown in blue. Here, we adopt cosmological parameter priors consistent with the CMB lensingonly analysis, which differ from the priors assumed by the DES collaboration (Troxel et al. 2018). The odd shape of the DES lensing and joint contours is due to a nontrivial degeneracy with the intrinsic alignment parameters, giving a region of parameter space with large negative intrinsic alignment amplitudes that cannot be excluded by current lensing data alone (but is reduced by different choices of cosmological parameter priors; see PCP18). Contours contain 68% and 95% of the probability. 

In the text 
Fig. 9. ΛCDM model constraints on Ω_{m}, σ_{8}, and H_{0} (in units of km s^{−1} Mpc^{−1}) from DES galaxy lensing+BAO (green), Planck CMB lensing+BAO (grey), and the joint constraint of DES lensing, CMB lensing, and galaxy BAO (red). All these results use the standard lensing priors described in the text, including the baryon density prior Ω_{b}h^{2} = 0.0222 ± 0.0005. The Planck powerspectrum constraints are shown in blue. Contours contain 68% and 95% of the probability. 

In the text 
Fig. 10. Empirical crossspectra between our lensing reconstructions (temperatureonly in blue, polarizationonly in orange, and MV in green) and the GNILC 545GHz CIB map. The red points use the temperature lensing reconstruction from the thermalSZdeprojected SMICA CMB map. The dashed line shows the smooth spline fit to the green points, which we use to weight the lensing tracers when forming our combined lensing potential estimate. 

In the text 
Fig. 11. Top: expected crosscorrelation coefficient of the true lensing potential to the minimumvariance (MV) estimator, the 545GHz GNILC CIB map (orange), and their combination (green). To build these curves, the crossspectrum of the CIB map to the true lensing potential is approximated by its crossspectrum to the MV (quadraticestimator) lensing reconstruction, displayed in Fig. 10, and GNILC CIB curves are smooth spline fits. The red curve shows for comparison the crosscorrelation coefficient of the CIB to the MV quadratic estimator. This crosscorrelation is suppressed by instrument noise, foreground residuals, and shot noise in the CIB map and reconstruction noise in the lensing quadratic estimator. The black curves show the lensing kernels that contribute to the lensing of the temperature (solid), Emode polarization (dashed), and Bmode polarization (dashdotted) power spectra, as described in the text. Bottom: effective reconstruction noise levels N_{L} for each of these tracers, as defined by Eq. (48), and, for comparison, the theoretical lensing spectrum . The quadratic estimator reconstruction noise is slightly underestimated at low L in this figure, since we have neglected Monte Carlo corrections when combining the tracers. 

In the text 
Fig. 12. Comparison of lensing maps constructed from the minimumvariance quadratic estimator alone (upper panel) and in combination with the CIB, as traced by the 545GHz GNILC frequency map (lower panel). The combination is performed on 60% of the sky, defined as the union of the lensing mask and the GNILC mask. Maps show the orthographic projection of the Wienerfiltered displacement E mode, the scalar field with multipoles , with 10 ≤ L ≤ 2000. The left and right panels are centred on the north and south Galactic poles, respectively. While the two reconstruction maps are clearly strongly correlated, the combined map has substantially more smallscale power, due to the higher S/N of the CIB on small scales. 

In the text 
Fig. 13. Estimates of the lensing Bmode power spectrum from crossspectra between the Bmode polarization map and various Bmode templates constructed using different lensing tracers, as measured on 58% of the sky: the TT+TE+EE quadratic estimator (blue); the CIB map at 545 GHz (orange); and their combination (green). The raw crossspectra are a filtered version of the Bmode spectrum and have been rescaled using the expected scaling computed from simulations assuming an isotropic filter. The dashed black line is the FFP10 fiducial lensing Bmode spectrum. 

In the text 
Fig. 14. Difference between the delensed and original Bmode power spectrum for the SMICA CMB polarization maps. In the upper panel, the delensed Bmode map is obtained by template subtraction. The templates are constructed using the labelled lensing tracers, and the Wienerfiltered Emode map, as discussed in Sect. 3.4.1. The dashed lines show predictions obtained by repeating these operations on the FFP10 set of Planck simulations. The black curve, for our fiducial model, shows the difference expected for perfect delensing. Summary statistics are presented in Table 3. In the lower panel, a remapping method is used, and the delensing signature is obtained after subtraction of biases, as described in Sect. 3.4.2, with summary statistics in Table 4. In this case, the actual total Bmode power is not decreased after delensing. 

In the text 
Fig. 15. Acoustic peak resharpening after delensing of the SMICA CMB maps. From top to bottom: difference between the delensed and original TT, TE, and EE spectra. Data points show delensing with our internal, minimumvariance, quadratic lensing estimate (blue), with the GNILC CIB map (orange), and with their optimal combination (green), as described in the text. For reference, the black line shows half the difference of the unlensed to lensed CMB spectra in the fiducial cosmology. In all cases, the solid coloured lines show the FFP10 simulation predictions, obtained after performing the same operations on each simulation as for the data. The dashed lines show the noise delensing bias defined by Eq. (55) (only relevant here for EE at high multipoles), while the dotdashed lines show the CMB peak smoothing caused by the tracer reconstruction noise, Eq. (56). 

In the text 
Fig. 16. Top: fractional differences between the lensing power spectrum band powers from the ϕ^{TT} reconstruction over the aggressive multipole range and the Planck 2018 bestfit cosmology using Planck TT,TE,EE+lowE+lensing (blue points). Bottom: temperatureonly curl reconstruction bandpower amplitudes. In both cases, the orange points show the same data with a coarser binning, where the construction of the binning function and midpoints follows the same procedure as the blue points (see Sect. 2.1). The curl bin at 264 ≤ L ≤ 901 deviates from zero with a formal significance of 4.3σ. This bin was chosen to maximize the deviation; after simple consideration of “lookelsewhere” effects as described in the text, we evaluate the PTE of that deviation to be 0.41% (or roughly 2.9 Gaussian σ), still suggestive of a problem in the curl reconstruction. 

In the text 
Fig. 17. Summary of lensingamplitudemeasurement consistency tests using temperatureonly reconstructions. Upper panel: difference between the measured lensing amplitude from pairs of reconstructions in units of the standard deviation of the difference from simulations. This tests whether the observed shifts are compatible with the behaviour expected from the FFP10 simulations. The lower triangle shows the amplitude measured from the conservative range 8 ≤ L ≤ 400, with a red colour indicating a larger amplitude of the reconstruction labelled on the left; the upper triangle uses the highL range (401 ≤ L ≤ 2048), with red indicating a larger amplitude of the reconstruction labelled at the top. Numerical significances are quoted for those tests that show differences greater than twice the expected standard deviation. The first entry is our baseline temperatureonly reconstruction; see the text for a description of all the other reconstructions. Lower panel: actual amplitudes (relative to our FFP10 fiducial model), with the result using the conservative range in blue and the highL range in orange. The dashed lines show for reference the values in our baseline reconstruction. 

In the text 
Fig. 18. Same as Fig. 17, but for the curl reconstruction. Here the lower panel shows the curl amplitude over the conservative multipole range (8 ≤ L ≤ 400; blue) and over the range covering the negative curl feature shown in Fig. 16 (264 ≤ L ≤ 901; orange). We use the same binning scheme as for the lensing gradient (Eq. (10)), with an amplitude of unity corresponding to a flat curl spectrum. 

In the text 
Fig. 19. Various masks used for the lensing analysis and consistency tests. The topleft mask is the default lensing mask, including cutting out point sources, resolved SZ clusters, and the Galaxy; the topmiddle mask shows the difference between the new mask and the one used for the 2015 analysis. The remaining masks are used for consistency tests, as described in Sect. 4.2. The dust masks are constructed by thresholding the smoothed pixel variance obtained from the product of the two halfmission 545GHz Planck maps after filtering to include only multipoles 1000 ≤ ℓ ≤ 2000. 

In the text 
Fig. 20. Lensing reconstruction band powers over the full aggressive multipole range for the minimumvariance estimator (MV), and separately for all auto and crossspectra of the five quadratic lensing estimators that can be formed from the temperature (T) and polarization (E and B) SMICA CMB maps. The black line shows the Planck TT,TE,EE+lowE+lensing bestfit spectrum. In each lower left corner the χ^{2} probability to exceed (PTE) is given for the reconstruction (evaluated at the fiducial lensing spectrum), calculated across our conservative multipole range. 

In the text 
Fig. 21. Same as Fig. 20, but for the lensing curl reconstruction. 

In the text 
Fig. 22. Halfmission lensing reconstructions tests for the gradient (ϕ; upper) and curl (Ω; lower). Blue points show reconstruction power spectra from SMICA halfmission difference maps, after multiplication by a factor of 20. The orange points show the autospectrum of the reconstruction built from HM_{1} and HM_{2} (with no noise mean field), and the green points the crossspectrum of the reconstruction built from HM_{1} with that built from HM_{2} (with no noise contribution to the reconstruction noise N^{(0)}). The black line in the upper plot shows the Planck TT,TE,EE+lowE+lensing bestfit spectrum. 

In the text 
Fig. 23. Estimates of the expected CIBinduced biases (as a fraction of the lensing potential power spectrum) to the temperatureonly lensing spectrum reconstruction at 217 GHz (without any cleaning), and for the SMICA frequency weighting. The largescalestructure bispectrum causes the quadratic estimator applied to the CIB map to correlate with the lensing signal, sourcing a 2% negative bias at 217 GHz (blue) on our conservative multipole range. This is reduced by an order of magnitude after SMICA cleaning (orange). The bias remains subpercent in the tSZdeprojected SMICA weighting (green). The CIB trispectrum sources a smaller contribution (red and purple at 217 GHz and for SMICA frequency weighting, respectively). Biases to the curl reconstruction are negligible. Error bars on this figure are estimated from the scatter of the unbinned power. 

In the text 
Fig. 24. Impact of the thermal SZ (tSZ) on the temperatureonly lensing reconstruction spectrum at 143 GHz. We show the fractional difference between several reconstruction power spectra based on a 143 GHz fullsky simulation and the input lensing potential power spectrum (where the simulation includes tSZ and CMB, and Gaussian noise generated according to the 143GHz channel pixel variance map). Without any masking, the contamination appears predominantly as a large pointsource trispectrum contribution (orange). The pointsource correction is quite effective at reducing this term below L ≃ 400, with visible residual contamination at higher L (green). After applying a f_{sky} = 0.997 cluster mask (red), the differences from the power reconstructed in the reference CMBonly simulation with no masking (blue) are further reduced to a small fraction of a standard deviation in all bin. The left panel shows the conservative multipole range and the right panel the higher lensing multipoles. 

In the text 
Fig. 25. Comparison of reconstructions over the same highL multipole range on the same simulated map, using different fiducial model assumptions. The blue points use a fiducial model, set of simulations, and estimator ingredients based on the FFP10 input cosmology. The orange points use a cosmology with Ω_{K} ≃ −0.04, where the lensing power spectrum is higher by 10–20% (black line), and the CMB spectra differ by roughly 2% at high ℓ. The linear corrections to the likelihood account for the differences in fiducial model, and are mostly important at highL in this case; they are effective (green) at bringing the reconstructions with the different fiducial models into good agreement. 

In the text 
Fig. 26. Comparison of reconstruction band powers for various methods to remove the N^{(1)} lensing bias. Our baseline band powers are built by subtracting the bias calculated analytically with the FFP10 fiducial model, and are shown as the blue points with error bars. The curves labelled “RDN^{(1)}” and “xpRDN^{(1)}” invert Eqs. (78) and (79), respectively, and make no assumption about the theory lensing spectrum. The curves labelled “MCN^{(1)}” subtract a simulationbased estimate of N^{(1)}. For comparison the red lines show the band powers without any N^{(1)} subtraction, and the dashed black lines the N^{(1)} bias itself. The top plot shows the lensing reconstruction from temperature only, and is normalized to the FFP10 fiducial . The lower plot shows the temperatureonly curl reconstruction (which on small scales appears to be coincidentally close to the result one would expect with no N^{(1)} at all). 

In the text 
Fig. 27. Comparison of the simulationbased N^{(1)} lensing bias estimates (blue points) to analytic estimates in the case of our inhomogeneouslyfiltered polarization reconstruction. The orange curve shows the naive prediction using a single mapaveraged noise level, which does not capture the variation of the N^{(1)} bias across the sky. The green curve shows the weighted average Eq. (80). The black line shows the fiducial lensing power spectrum. 

In the text 
Fig. A.1. Compilation of the power spectrum reconstruction biases for our baseline MV analysis. The blue and orange lines show RDN^{(0)} and the analytic N^{(1)} estimate. The green curve shows the pointsource correction, and the red points the additive Monte Carlo correction (shown here after a coarse binning was applied) defined by the difference of Eq. (A.2) to the simulation input fiducial lensing spectrum. This correction is most prominent at low lensing multipoles, where it is negative (dashed) and due to masking. The error bars are obtained from the 240 simulations used to calculate the correction. The black line shows the fiducial lensing power spectrum. The black dot at L = 1 shows the almost pure dipole effective deflection caused by our motion with respect to the CMB frame (aberration). This dipole is included in our simulations, hence subtracted from our lensing reconstructions together with other sources of anisotropy through the mean field. 

In the text 
Fig. B.1. Characterization of the meanfield power spectrum and its components in our temperatureonly reconstruction. The lefthand part shows the convergencelike power spectrum on large scales, while the righthand part shows the deflection power spectrum on smaller scales. The blue line shows the total mean field as obtained from the SMICA FFP10 simulation set. Lines showing the contributions from statistical anisotropy of the mask (red), noise (orange), and the beamconvolved and pixelized CMB anisotropies (green) are plotted with rough 68% confidence regions shown shaded on the right panel (the uncertainty arises from the finite number of simulations used to estimate the mean fields). The mask mean field dominates at low multipoles and is barely distinguishable from the total on the lefthand plot. The orange and green curves were obtained from fullsky lensing reconstructions using idealized, statisticallyisotropic CMB or noise components, respectively. The red curve was obtained by differencing the meanfield spectra on the full anisotropic simulations, with and without masking. Pixelization effects due to subpixel pointing offsets are expected to appear as an approximately whitenoise lensing deflection component of amplitude 0.′05 (the brown line shows the prediction for the 217GHz channel), and is clearly detected by crosscorrelating the mean field to the subpixel deflection prediction (pink curve). The rise of the mean field at L ≃ 3000 originating in the noise maps is mostly a simulation artefact, sourced by nonlinearities in the data processing that cause slight correlations between the simulations (due to the fixed fiducial CMB and foregrounds used to make the noise simulations); the purple curve shows the lensing quadratic estimator applied to the empirical average of the noise simulations. The predicted mean field due to inhomogeneity of the variance of idealized uncorrelated pixel noise (grey curve, as derived from Eq. (B.1)) is much weaker on smaller scales and is only visible at low lensing multipoles in this figure. 

In the text 
Fig. C.1. Crude empirical estimates of the additional variance on the MV reconstruction band powers (on the conservative multipole range) due to the finite number of simulations used for meanfield and bias correction. The additional variance is normalized by the baseline statistical variance of the band powers. These empirical estimates were obtained by splitting the simulations into five independent subsets, building reconstruction band powers using each subset of the simulations, and measuring the scatter between the five sets of band powers. The results are scaled by 1/5 to approximate the full simulation set used in the main analysis. The horizontal line shows the simple constant relative correction applied to the covariance matrix, Eq. (13), as obtained in Appendix C. The error bars on the empirical variance (identical for each bin, after rescaling by ) are those expected in the fiducial model of Eq. (13), assuming Gaussian statistics of the lensing maps. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.