© ESO, 2014

1. Introduction

This paper, one of a set of papers associated with the 2013 release of data from the Planck¹ mission (Planck Collaboration I 2014), describes our reconstruction of the cosmic microwave background (CMB) lensing potential based on 15 months of data, estimation of the lensing potential power spectrum, and a first set of associated science results.

When Blanchard and Schneider first considered the effect of gravitational lensing on the CMB anisotropies in 1987, they wrote with guarded optimism that although “such an observation is far from present possibilities [...] such an effect will not be impossible to find and to identify in the future” (Blanchard & Schneider 1987). In the proceeding years, and with the emergence of the concordance ΛCDM cosmology, a standard theoretical picture has emerged, in which the large-scale, linear structures of the Universe which intercede between ourselves and the CMB last-scattering surface induce small but coherent (Cole & Efstathiou 1989) deflections of the observed CMB temperature and polarisation anisotropies, with a typical magnitude of 2′. These deflections blur the acoustic peaks (Seljak 1996), generate small-scale power (Linder 1990; Metcalf & Silk 1997), non-Gaussianity (Bernardeau 1997), and convert a portion of the dominant E-mode polarisation to B-mode (Zaldarriaga & Seljak 1998). Gravitational lensing of the CMB is both a nuisance, in that it obscures the primordial fluctuations (Knox & Song 2002), as well as a potentially useful source of information; the characteristic signatures of lensing provide a measure of the distribution of mass in the Universe at intermediate redshifts (typically 0.1 <z< 5). In the ΛCDM framework, there exist accurate methods to calculate the effects of lensing on the CMB power spectra (Challinor & Lewis 2005), as well as optimal estimators for the distinct statistical signatures of lensing (Hu & Okamoto 2002; Hirata & Seljak 2003a).

In recent years there have been a number of increasingly sensitive experimental measurements of CMB lensing. Lensing has been measured in the data of the WMAP satellite both in cross-correlation with large-scale-structure probed by galaxy surveys (Hirata et al. 2004, 2008; Smith et al. 2007; Feng et al. 2012a), as well as internally at lower signal-to-noise (Smidt et al. 2011; Feng et al. 2012b). The current generation of low-noise, high-resolution ground-based experiments has done even better; the Atacama Cosmology Telescope (ACT) has provided an internal detection of lensing at 4.6σ (Das et al. 2011, 2013), and the South Pole Telescope detects lensing at 8.1σ in the temperature power spectrum, and 6.3σ from a direct reconstruction of the lensing potential (Keisler et al. 2011; van Engelen et al. 2012; Story et al. 2013). Significant measurements of the correlation between the reconstructed lensing potential and other tracers of large-scale structure have also been observed (Bleem et al. 2012; Sherwin et al. 2012).

Planck enters this field with unique full-sky, multi-frequency coverage. Nominal map noise levels for the first data release (approximately 105, 45, and 60 μK arcmin for the three CMB channels at 100, 143, and 217 GHz respectively) are approximately five times lower than those of WMAP (or twenty five times lower in power), and the Planck beams (approximately 10′, 7′ and 5′ at 100, 143, and 217 GHz), are small enough to probe the deflections typical of lensing. Full sky coverage is particularly beneficial for the statistical analysis of lensing effects, as much of the “noise” in temperature lens reconstruction comes from CMB fluctuations themselves, which can only be beaten down by averaging over many modes.

Lensing performs a remapping of the CMB fluctuations, such that the observed temperature anisotropy in direction $\hat{{n}}$ is given in terms of the unlensed, “primordial” temperature anisotropy as (e.g. Lewis & Challinor 2006) $\begin{array}{ccc}\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{\right)}& \mathrm{=}& {\mathit{T}}^{\mathrm{unl}}\mathrm{(}\hat{{n}}\mathrm{+}\mathrm{\nabla }\mathit{\phi }\mathrm{(}\hat{{n}}\mathrm{\left)}\mathrm{\right)}\mathit{,}\\ & \mathrm{=}& \end{array}$(1)where $\mathit{\phi }\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ is the CMB lensing potential, defined by $\mathit{\phi }\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}\mathrm{-}\mathrm{2}{\mathrm{\int }}_{\mathrm{0}}^{{\mathit{\chi }}_{\mathrm{\ast }}}\mathrm{d}\mathit{\chi }\frac{{\mathit{f}}_{\mathit{K}}\mathrm{(}{\mathit{\chi }}_{\mathrm{\ast }}\mathrm{-}\mathit{\chi }\mathrm{)}}{{\mathit{f}}_{\mathit{K}}\mathrm{\left(}{\mathit{\chi }}_{\mathrm{\ast }}\mathrm{\right)}{\mathit{f}}_{\mathit{K}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}}\mathrm{\Psi }\mathrm{(}\mathit{\chi }\hat{{n}}\mathrm{;}{\mathit{\eta }}_{\mathrm{0}}\mathrm{-}\mathit{\chi }\mathrm{)}\mathit{.}$(2)Here χ is conformal distance (with χ_∗ ≈ 14 000 Mpc denoting the distance to the CMB last-scattering surface) and $\mathrm{\Psi }\mathrm{\left(}\mathit{\chi }\hat{{n}}\mathit{,\eta }\mathrm{\right)}$ is the (Weyl) gravitational potential at conformal distance χ along the direction $\hat{{n}}$ at conformal time η (the conformal time today is denoted as η₀). The angular-diameter distance f_K(χ) depends on the curvature of the Universe, and is given by ${\mathit{f}}_{\mathit{K}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\mathrm{=}\{\begin{array}{c}\\ & \mathrm{forK}\mathrm{>}\mathrm{0}\mathrm{\left(}\mathrm{closed}\mathrm{\right)}\mathit{,}\\ & \mathrm{forK}\mathrm{=}\mathrm{0}\mathrm{\left(}\mathrm{flat}\mathrm{\right)}\mathit{,}\\ & \mathrm{forK}\mathrm{<}\mathrm{0}\mathrm{\left(}\mathrm{open}\mathrm{\right)}\mathit{.}\end{array}$(3)The lensing potential is a measure of the integrated mass distribution back to the last-scattering surface. To first order, its effect on the CMB is to introduce a correlation between the lensed temperature and the gradient of the unlensed temperature, a property which can be exploited to make a (noisy) reconstruction of the lensing potential itself.

In Fig. 1 we plot the noise power spectrum ${\mathit{N}}_{\mathit{L}}^{\mathit{\phi \phi }}$ for reconstruction of the lensing potential using the three Planck frequencies which are most sensitive to the CMB anisotropies on the arcminute angular scales at which lensing effects become apparent.

The angular size of the Planck beams (5′ FWHM and greater) does not allow a high signal-to-noise ratio (S/N) reconstruction of the lensing potential for any individual mode (our highest S/N on an individual mode is approximately 2/3 for the 143 and 217 GHz channels, or 3/4 for a minimum-variance combination of both channels), however with full-sky coverage the large number of modes that are probed provides considerable statistical power. To provide a feeling for the statistical weight of different regions of the lensing measurement, in Fig. 2 we plot (forecasted) contributions to the total detection significance for the potential power spectrum ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ as a function of lensing multipole L. In addition to the power spectrum of the lensing potential, there is tremendous statistical power in cross-correlation of the Planck lensing potential with other tracers of the matter distribution. In Fig. 2 we also plot forecasted S/N contributions for several representative tracers.

This paper describes the production, characterization, and first science results for two Planck-derived lensing products:

• 1.

  A map of the CMB lensing potential $\mathit{\phi }\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ over a large fraction of the sky (approximately 70%). This represents an integrated measure of mass in the entire visible Universe, with a peak sensitivity to redshifts of z ~ 2. At the resolution of Planck, this map provides an estimate of the lensing potential down to angular scales of 5′ at L = 2048, corresponding to structures on the order of 3 Mpc in size at z = 2.

• 2.

  An estimate of the lensing potential power spectrum ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ and an associated likelihood, which is used in the cosmological parameter analysis of Planck Collaboration XVI (2014). Our likelihood is based on the lensing multipole range 40 ≤ L ≤ 400. This multipole range (highlighted as a dark grey band in Fig. 2), was chosen as the range in which Planck has the greatest sensitivity to lensing power, encapsulating over 90% of the anticipated signal-to-noise, while conservatively avoiding the low-L multipoles where mean-field corrections due to survey anisotropy (discussed in Appendix C) are large, and the high-L multipoles where there are large corrections to the power spectra from Gaussian (disconnected) noise bias. Distilled to a single amplitude, our likelihood corresponds to a 4% measurement of the amplitude of the fiducial ΛCDM lensing power spectrum, or a 2% measurement of the amplitude of the matter fluctuations (neglecting parameter degeneracies).

Our efforts to validate these products are aided by the frequency coverage of the three Planck channels that we employ, which span a wide range of foreground, beam, and noise properties. For the mask levels that we use, the root-mean-squared (RMS) foreground contamination predicted by the Planck sky model (Delabrouille et al. 2013) has an amplitude of 14, 22, and 70 μK at 100, 143, and 217 GHz, which can be compared to a CMB RMS for the Planck best-fitting ΛCDM power spectrum of approximately 110 μK. The dominant foreground component at all three CMB frequencies is dust emission, both from our Galaxy as well as the cosmic infrared background (CIB), although at 100 GHz free-free emission is thought to constitute approximately 15% of the foreground RMS. Contamination from the thermal Sunyaev-Zeldovich (tSZ) effect is a potential worry at 100 and 143 GHz, but negligible at 217 GHz (Sunyaev & Zeldovich 1980). On the instrumental side, these frequency channels also span a wide range of beam asymmetry, with typical ellipticities of 19%, 4%, and 18% at 100, 143, and 217 GHz. The magnitude of correlated noise on small scales (due to deconvolution of the bolometer time response) also varies significantly. The ratio of the noise power (before beam deconvolution) at ℓ = 1500 to that at ℓ = 500 is a factor of 1.5, 1.1, and 1.0 at 100, 143, and 217 GHz. The agreement of lens reconstructions based on combinations of these three channels allows a powerful suite of consistency tests for both foreground and instrumental biases. We will further validate the robustness of our result to foreground contamination using the component-separated maps from the Planck consortium (Planck Collaboration XII 2014).

At face value, the 4% measurement of ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ in our fiducial likelihood corresponds to a 25σ detection of gravitational lensing effects. In fact, a significant fraction (approximately 25% of our error bar) is due to sample variance of the lenses themselves, and so the actual “detection” of lensing effects (under the null hypothesis of no lensing) is significantly higher. We have also been conservative in terms of mask and multipole range in the construction of our fiducial lensing likelihood. As we will show in Sect. 7.1, we obtain consistent results on sky fractions larger than our fiducial 70% sky mask.

The Planck lensing potential is part of a significant shift for CMB lensing science from the detection regime to that of precision cosmological probe. The NVSS quasar catalogue, for example, has been a focus of previous lensing cross-correlation studies with WMAP (Hirata et al. 2004; Smith et al. 2007; Hirata et al. 2008), where evidence for cross-correlation was found at approximately 3.5σ. As we will see in Sect. 6.3, the significance for this correlation with Planck is now 20σ. Notably, this is less than the significance with which lensing may be detected internally with Planck. The lensing potential measured by Planck now has sufficient signal-to-noise that shot noise of the NVSS quasar catalogue is the limiting source of noise in the cross-correlation.

The majority of this paper is dedicated to the production and testing of the Planck lensing map and power spectrum estimate. Our focus here is on extracting the non-Gaussian signatures of lensing, although we note that lensing effects are also apparent at high significance (10σ) as a smoothing effect in the Planck temperature power spectra (Planck Collaboration XV 2014). We begin in Sect. 2, where we describe and motivate our methodology for producing unbiased estimates of the lensing potential and its power spectrum. The Planck maps and data cuts that are used for this purpose are described in Sect. 3, and the simulations that we use to characterize our reconstruction and its uncertainties are described in Sect. 4. In Sect. 5 we give an overview of our error budget, and discuss the various sources of systematic and statistical uncertainty for our lensing estimates. In Sect. 6 we present our main results: the first Planck lensing map and a corresponding estimate of the lensing potential power spectrum. The likelihood based on this power spectrum is combined with the likelihood for the temperature anisotropy power spectrum (Planck Collaboration XV 2014) to derive parameter constraints in Planck Collaboration XVI (2014). In Sect. 6.1 we highlight a subset of parameter results where the information provided by the lensing likelihood has proven particularly useful. In the concordance ΛCDM cosmology, there is believed to be a correlation between the CMB lensing potential and the low-ℓ temperature anisotropies, driven by the effects of dark energy. We also present a measurement of this correlation in Sect. 6.2. Finally, we connect our lensing potential map to other tracers of large-scale structure with several illustrative cross-correlations using galaxy, quasar, cluster and infrared source catalogues in Sect. 6.3. These main results are followed in Sect. 7 by a large suite of systematic and consistency tests, where we perform null tests against a variety of different data cuts and processing. We conclude in Sect. 8. A series of appendices provide further details on some technical aspects of our methodology and lensing potential estimates.

Throughout this paper, when we refer to the concordance or fiducial ΛCDM cosmology we are referring to a model with baryon density ω_b = Ω_bh² = 0.0221, cold dark matter density ω_c = Ω_ch² = 0.1199, density parameter for the cosmological constant Ω_Λ = 0.6910, Hubble parameter H₀ = 100 h km s^-1 Mpc^-1 with h = 0.6778, spectral index of the power spectrum of the primordial curvature perturbation n_s = 0.96, amplitude of the primordial power spectrum (at k = 0.05 Mpc^-1) A_s = 2.21 × 10^-9, and Thomson optical depth through reionization τ = 0.093. These values were determined from a pre-publication analysis of the Planck temperature power spectrum, but are consistent with the best-fit values quoted in Planck Collaboration XVI (2014). Throughout this work we shall frequently quote lensing bandpower amplitudes relative to this fiducial model for ease of comparison, although as discussed in Sect. 2 our lensing likelihood itself is designed to be insensitive to the choice of fiducial model.

2. Methodology

In this section, we detail our methodology for reconstructing the lensing potential and estimating its angular power spectrum. These are both obtained by exploiting the distinctive statistical properties of the lensed CMB.

• 1.

  If we consider a fixed lensing potential applied to multiple realizations of the CMB temperature anisotropies, then lensing introduces statistical anisotropy into the observed CMB; the fluctuations are still Gaussian, however the covariance varies as a function of position and orientation on the sky. We use this idea to obtain a (noisy) estimate of $\mathit{\phi }\mathrm{\left(}\hat{{n}}\mathrm{\right)}$. The noise of this map is a combination of instrumental noise and statistical noise due to the fact that we only have a single realization of the CMB to observe, analogous to shape noise in galaxy lensing.

• 2.

  If we consider averaging over realizations of both the lensing potential and the CMB fluctuations, then lensing introduces non-Gaussianity into the observed CMB. This appears at lowest order in the connected part of the CMB 4-point function, or trispectrum². We use this to measure the lensing power spectrum ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$.

The estimators that we use are derived from maximizing the likelihood function of the lensed CMB under the assumption that the instrumental noise is Gaussian and the lensed CMB is perturbatively Gaussian and statistically isotropic. These estimators are optimal (in the minimum-variance sense). In cases where we have made suboptimal analysis choices, we provide estimates of the loss of signal-to-noise.

2.1. Lens reconstruction

To gain intuition for the process of lens reconstruction, it is useful to consider the effect of lensing on a small patch of the sky. Lensing remaps the temperature fluctuations by a deflection field $\mathrm{\nabla }\mathit{\phi }\mathrm{\left(}\hat{{n}}\mathrm{\right)}$. The part of $\mathrm{\nabla }\mathit{\phi }\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ that is constant over our patch is not an observable effect; it describes only a re-centering of the map. The variation of the deflection field across the patch is observable, however. This can be usefully decomposed into convergence (κ) and shear modes (γ₊, γ₋) as $\mathrm{-}{\mathrm{\nabla }}_{\mathit{i}}{\mathrm{\nabla }}_{\mathit{j}}\mathit{\phi }\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}\left[\begin{array}{c}\\ \mathit{\kappa }\mathrm{+}{\mathit{\gamma }}_{\mathrm{+}}& {\mathit{\gamma }}_{\mathrm{-}}\\ {\mathit{\gamma }}_{\mathrm{-}}& \mathit{\kappa }\mathrm{-}{\mathit{\gamma }}_{\mathrm{+}}\end{array}\right]\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{.}$(4)If we observe a patch that is small enough that these quantities can be taken as constant, then the observational consequences are simple. The convergence mode causes a local change of scale, either magnifying or demagnifying the fluctuations. Taking the local power spectrum of our small patch, we would find that the CMB peaks would shift to larger or smaller scales, relative to the full-sky average. The shear modes also describe changes of scale, however they are now orientation dependent. On a small patch, convergence and shear estimators can be constructed from local estimates of the (orientation-dependent) power spectrum and then stitched together to recover the lensing potential φ (Zaldarriaga & Seljak 1998; Bucher et al. 2012). This procedure describes a quadratic estimator for the local convergence and shear.

From the description above, it is not immediately clear how to go about stitching together estimates of convergence and shear in different regions of the sky, or what weight to give the local power spectrum estimates as a function of scale. These questions can be resolved by considering a generic form for the quadratic estimator, and optimizing its weight function for sensitivity to lensing (Okamoto & Hu 2003). To first order in the lensing potential, the statistical anisotropy introduced by lensing appears as an off-diagonal contribution to the covariance matrix of the CMB: $\mathrm{\Delta }\mathrm{⟨}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{⟩}\mathrm{=}\sum _{\mathit{LM}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& \mathit{L}\\ {\mathit{m}}_{\mathrm{1}}& {\mathit{m}}_{\mathrm{2}}& \mathrm{-}\mathit{M}\end{array}\right){\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}{\mathit{\phi }}_{\mathit{LM}}\mathit{,}$(5)where the average ⟨⟩ is taken over CMB realizations with a fixed lensing potential. Here the bracketed term is a Wigner 3j symbol, ${\mathit{\phi }}_{\mathit{LM}}\mathrm{=}{}^{\mathrm{\int }}\mathrm{d}^{\mathrm{2}}\hat{{n}}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{\phi }\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ is the harmonic transform of the lensing potential, and the weight function ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}$ is given by $\begin{array}{ccc}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}\mathrm{=}\mathrm{-}\sqrt{\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}}\sqrt{\mathit{L}\mathrm{(}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\ell }}_{\mathrm{1}}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}}& & \\ \mathrm{\times }{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{TT}}\left(\frac{\mathrm{1}\mathrm{+}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathit{L}}}{\mathrm{2}}\right)\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\mathrm{+}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{\leftrightarrow }{\mathit{\ell }}_{\mathrm{2}}\mathrm{)}\mathit{.}& & \end{array}$(6)Here ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$ is the ensemble-average power spectrum of the lensed CMB. In our analysis, we will use the fiducial model described at the end of Sect. 1 to determine ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$, however our lensing likelihood can be renormalized to account for uncertainties in this model. In this approach, errors in the fiducial model do not bias our lensing bandpower estimates, they only result in slight sub-optimality. Note that we use the lensed power spectrum here, rather than the unlensed spectrum that is sometimes used in the literature, as this is accurate to higher order in φ (Lewis et al. 2011), an improvement which is necessary at Planck sensitivity (Hanson et al. 2011). Use of the unlensed spectrum would lead to biases on the order of 15% at L< 200.

Now we construct a quadratic estimator to search for the covariance which is introduced by lensing. We will use several different estimators for the lensing potential, as well as to probe possible point-source contamination, and so it will be useful to keep this discussion as general as possible. A completely generic quadratic estimator for the lensing potential can be written as $\mathit{\phi ̂}\begin{array}{c}\mathit{x}\\ \mathit{LM}\end{array}\mathrm{=}\sum _{{\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}}{\mathrm{\left[}{\mathrm{ℛ}}^{\mathit{x\phi }}\mathrm{\right]}}_{\mathit{LM,}{\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}}^{-1}\mathrm{[}\mathit{x̅}{\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}\mathrm{-}\mathit{x̅}{\begin{array}{c}\mathrm{MF}\\ {\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}\end{array}}^{\mathrm{]}}\mathit{,}$(7)where ℛ^xφ is a normalization matrix, and is a quadratic “building block” which takes in a pair of filtered sky maps $\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\\ \mathit{\ell m}\end{array}$ and $\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ \mathit{\ell m}\end{array}$, and sums over their empirical covariance matrix with a weight function ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{x}}$: $\mathit{x̅}\mathit{LM}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& \mathit{L}\\ {\mathit{m}}_{\mathrm{1}}& {\mathit{m}}_{\mathrm{2}}& \mathrm{-}\mathit{M}\end{array}\right){\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{x}}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\\ {\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\end{array}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ {\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}\end{array}\mathit{.}$(8)The “mean-field” term $\mathit{x̅}\begin{array}{c}\mathrm{MF}\\ \mathit{LM}\end{array}$ accounts for all known sources of statistical anisotropy in the map, which could otherwise bias the lensing estimate. It is given by $\mathit{x̅}\begin{array}{c}\mathrm{MF}\\ \mathit{LM}\end{array}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\mathit{,}{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& \mathit{L}\\ {\mathit{m}}_{\mathrm{1}}& {\mathit{m}}_{\mathrm{2}}& \mathrm{-}\mathit{M}\end{array}\right){\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{x}}\mathrm{⟨}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\\ {\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}\end{array}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ {\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}\end{array}\mathrm{⟩}\mathit{,}$(9)where the ensemble average here is taken over realizations of the CMB and noise.

We may now optimize the generic quadratic estimator above. If the primordial CMB fluctuations and instrumental noise are Gaussian and the lensing potential is fixed, then the likelihood for the observed CMB fluctuations is still a Gaussian, which may be maximized with respect to the lensing potential modes φ_LM (Hirata & Seljak 2003a). The optimal quadratic estimator is the first step of an iterative maximization of this likelihood, and it has been shown that additional iterations of the estimator are not necessary for temperature lens reconstruction (Hirata & Seljak 2003a; Okamoto & Hu 2003). The optimal quadratic estimator has the following choices for the weight function and filtering.

• 1.

  The weight function W^x should be a matched filter for the covariance induced by lensing (i.e., one should use , with weight function given by Eq. (6)). We shall use this weight function for all of our fiducial results, although for consistency tests we will also use “bias-hardened” estimators, which have weight functions constructed to be orthogonal to certain systematic effects (Namikawa et al. 2013). This is discussed further in Sect. 7.4.

• 2.

  The filtered temperature multipoles should be given by , where T is a beam-deconvolved sky map and C is its total signal+noise covariance matrix. We describe our approximate implementation of this filtering in Appendix B. When combining multiple frequencies for our minimum-variance estimator, all of the available data are combined into a single map which is then filtered and used for both input multipoles of the quadratic estimator. It can be desirable to use different pairs of maps however, and we use this for several consistency tests. For example, we feed maps with independent noise realizations into the quadratic estimator to avoid possible noise biases in Sect. 7.3.

In the quadratic maximum-likelihood estimator, the mean-field correction emerges from the determinant term in the likelihood function, and it can be seen that the normalization matrix ℛ is the Fisher matrix for the φ_LM; this means that the normalization is the same as the covariance matrix of the lens reconstruction, and so the unnormalized lensing estimate is equivalent to an inverse-variance-weighted lens reconstruction, which is precisely the quantity needed for most statistical analysis. This is why we have denoted it with an overbar, in analogy to .

We choose to treat the map noise as if it were homogeneous when constructing the filtered , and do not account for variation with hit count across the sky. This is a slightly suboptimal filtering choice; in Appendix. B we estimate that it leads to a 5% loss of total signal-to-noise when constraining the power spectrum of the lensing potential. The advantage of this approach, however, is that far from the mask boundaries our filtering asymptotes to a simple form, given by $\mathit{T̅}\mathit{\ell m}\mathrm{\approx }{\mathrm{[}{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\mathrm{+}{{\mathit{C}}_{\mathit{\ell }}^{\mathit{NN}}}^{\mathrm{]}}}^{-1}{\mathit{T}}_{\mathit{\ell m}}\mathrm{\equiv }{\mathit{F}}_{\mathit{\ell }}{\mathit{T}}_{\mathit{\ell m}}\mathit{,}$(10)where ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$ is the temperature power spectrum and ${\mathit{C}}_{\mathit{\ell }}^{\mathit{NN}}$ is the power spectrum of the homogeneous noise level that we use in our filtering. For the purposes of compact notation, in the following equations we combine both of these elements in the “filter function” F_ℓ. The asymptotic form of our filtering, Eq. (10), will prove useful, as it means that the normalization of our estimator, as well as its variance and response to various systematic effects, may be accurately modelled analytically. It allows us to propagate uncertainties in the beam transfer function and CMB power spectrum, for example, directly to our lens reconstruction. This filtering choice also means that the normalization does not vary as a function of position on the sky, which simplifies the analysis of cross-correlations between the lensing potential map and external tracers. Under the approximation of Eq. (10), the estimator normalization is given by ${\mathrm{ℛ}}_{\mathit{LM,}{\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}}^{\mathit{x\phi }}\mathrm{=}{\mathit{\delta }}_{\mathit{L}{\mathit{L}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{M}{\mathit{M}}^{\mathrm{\prime }}}{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{x\phi }}\mathit{,}$(11)where the response function ${\mathrm{ℛ}}_{\mathit{L}}^{\mathit{x\phi }}$ for filtered maps and is ${\mathrm{ℛ}}_{\mathit{L}}^{\mathit{x\phi ,}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{x}}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}\mathit{.}$(12)This can be read as “the response of estimator x to lensing on scale L”. The filter functions F_ℓ are those used for and respectively. In cases where the filter functions are obvious, we will drop the indices above.

Putting all of the above together, for a chosen quadratic estimator we obtain normalized, mean-field-debiased estimates of the lensing potential φ as $\mathit{\phi ̂}\begin{array}{c}\mathit{x}\\ \mathit{LM}\end{array}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{x\phi }}}\mathrm{(}\mathit{x̅}\mathit{LM}\mathrm{-}\mathit{x̅}{\begin{array}{c}\mathrm{MF}\\ \mathit{LM}\end{array}}^{\mathrm{)}}\mathit{.}$(13)Note that our normalization function ${\mathrm{ℛ}}_{\mathit{L}}^{\mathit{x\phi }}$ is only approximate, but we will verify its accuracy in Sect. 4. For the standard lensing estimator of Okamoto & Hu (2003; which uses the weight function of Eq. (6)), we use x = φ. This estimator is denoted simply as .

2.2. Lensing power spectrum estimation

We form estimates for the power spectrum of the lensing potential by taking spectra of the lensing estimates from Sect. 2.1, using a simple pseudo-C_ℓ estimator. In order to reduce mode coupling, as well as to downweight regions near the analysis boundary where the mean-field due to masking can be large, we take the power spectrum from an apodized version of our lensing estimate, given by $\begin{array}{c}\mathrm{􏽥}\\ {\mathit{\phi }}_{\mathit{LM}}^{\mathit{x}}\end{array}\mathrm{=}{\mathit{P}}_{\mathit{L}}^{-1}\mathrm{\int }\mathrm{d}\hat{{n}}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\begin{array}{c}\mathrm{􏽥}\\ \mathit{M}\end{array}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\left[\sum _{{\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}}{\mathit{Y}}_{{\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\mathit{P}}_{{\mathit{L}}^{\mathrm{\prime }}}\mathit{\phi ̂}\begin{array}{c}\mathit{x}\\ {\mathit{L}}^{\mathrm{\prime }}{\mathit{M}}^{\mathrm{\prime }}\end{array}\right]\mathit{,}$(14)where $\begin{array}{c}{}_{\mathrm{􏽥}}\\ \mathit{M}\end{array}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ is an apodized version of the analysis mask $\mathit{M}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ used in our filtering and P_L ≡ L(L + 1) is an approximate pre-whitening operation. The construction of $\begin{array}{c}{}_{\mathrm{􏽥}}\\ \mathit{M}\end{array}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ is described in Sect. 3. Our fiducial apodization occurs over a band of approximately 5deg, and effectively reduces the sky fraction by 9%.

The power spectrum of $\begin{array}{c}{\mathrm{􏽥}}^{}\\ \mathit{\phi }\end{array}$ probes the 4-point function of the observed CMB, which contains both disconnected and connected parts. We model it as being due to a combination of Gaussian CMB fluctuations, lensing effects and unresolved point-source shot noise, and estimate the power spectrum of the lensing potential with $\begin{array}{ccc}\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L,x}\end{array}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{sky}\mathit{,}\mathrm{2}}^{-1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{\mathit{M}}\mathrm{|}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{\phi }}_{\mathit{LM}}^{\mathit{x}}\end{array}{\mathrm{|}}^{\mathrm{2}}\mathrm{-}{\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{0}}& & \\ \mathrm{-}{\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{1}}\mathrm{-}{\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{PS}}\mathrm{-}{\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{MC}}\mathit{,}& & \end{array}$(15)where ${\mathit{f}}_{\mathrm{sky}\mathit{,}\mathrm{2}}\mathrm{=}{}^{\mathrm{\int }}\mathrm{d}\hat{{n}}\begin{array}{c}{}_{\mathrm{􏽥}}\\ {\mathit{M}}^{\mathrm{2}}\end{array}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{/}\mathrm{4}\mathit{\pi }$ is the average value of the square of the apodizing mask. The first line of Eq. (15)isolates the connected part of the CMB 4-point function, or trispectrum, which would be zero for Gaussian fluctuations. The second line contains corrections which isolate the part of the trispectrum which is directly proportional to the non-Gaussianity induced by ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$. In the following paragraphs, we explain these terms in more detail.

The first correction term ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{0}}$ subtracts the (large) disconnected contribution to the power spectrum of $\begin{array}{c}{\mathrm{􏽥}}^{}\\ \mathit{\phi }\end{array}$. To determine this term, we use the data-dependent subtraction which emerges for maximum-likelihood estimators of the CMB trispectrum (Regan et al. 2010; see also Appendix D). For lensing, this procedure has the additional advantage of reducing the correlation between different multipoles L ≠ L′ of the lens reconstruction to a level which is negligible at Planck resolution and noise levels (Schmittfull et al. 2013), as well as reducing sensitivity to uncertainties in our model of the CMB and noise covariance matrices (Namikawa et al. 2013). Writing the power spectrum of $\begin{array}{c}{\mathrm{􏽥}}^{}\\ {\mathit{\phi }}_{\mathit{LM}}\end{array}$ explicitly as a function of the four inverse-variance filtered temperature maps ${\mathit{C}}_{\mathit{L,x}}^{\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}}\mathrm{\left[}\mathit{T̅}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathit{,}\mathit{T̅}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\mathit{,}\mathit{T̅}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathit{,}\mathit{T̅}\mathrm{\left(}\mathrm{4}\mathrm{\right)}\mathrm{\right]}\mathrm{\equiv }\frac{{\mathit{f}}_{\mathrm{sky}\mathit{,}\mathrm{2}}^{-1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{\mathit{M}}\mathrm{|}\begin{array}{c}\mathrm{􏽥}\\ {\mathit{\phi }}_{\mathit{LM}}^{\mathit{x}}\end{array}{\mathrm{|}}^{\mathrm{2}}\mathit{,}$(16)the disconnected contribution reads $\begin{array}{ccc}{\mathrm{\Delta }{\mathit{C}}_{\mathit{L,x}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{0}}\mathrm{=}⟨\mathrm{-}{\mathit{C}}_{\mathit{L,x}}^{\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}}\mathrm{[}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ {\mathit{mc}}^{\mathrm{\prime }}\end{array}\mathit{,}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\\ {\mathit{mc}}^{\mathrm{\prime }}\end{array}\mathit{,}\mathit{T̅}{\begin{array}{c}\mathrm{\left(}\mathrm{4}\mathrm{\right)}\\ \mathit{mc}\end{array}}^{\mathrm{]}}& & \\ \mathrm{+}{\mathit{C}}_{\mathit{L,x}}^{\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}}\mathrm{[}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\mathit{,}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}^{\mathrm{]}}\mathrm{+}{\mathit{C}}_{\mathit{L,x}}^{\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}}\mathrm{[}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\mathit{,}\mathit{T̅}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathit{,}\mathit{T̅}{\begin{array}{c}\mathrm{\left(}\mathrm{4}\mathrm{\right)}\\ \mathit{mc}\end{array}}^{\mathrm{]}}& & \\ \mathrm{+}{\mathit{C}}_{\mathit{L,x}}^{\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}}\mathrm{[}\mathit{T̅}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathit{,}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathit{,}\mathit{T̅}{\begin{array}{c}\mathrm{\left(}\mathrm{4}\mathrm{\right)}\\ \mathit{mc}\end{array}}^{\mathrm{]}}\mathrm{+}{\mathit{C}}_{\mathit{L,x}}^{\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}}\mathrm{[}\mathit{T̅}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathit{,}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}^{\mathrm{]}}& & \\ \mathrm{-}{\mathit{C}}_{\mathit{L,x}}^{\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}\begin{array}{c}\mathrm{􏽥}\\ \mathit{\phi }\end{array}}\mathrm{[}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\\ {\mathit{mc}}^{\mathrm{\prime }}\end{array}\mathit{,}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\\ \mathit{mc}\end{array}\mathit{,}\mathit{T̅}{\begin{array}{c}\mathrm{\left(}\mathrm{4}\mathrm{\right)}\\ {\mathit{mc}}^{\mathrm{\prime }}\end{array}}^{\mathrm{]}}{⟩}_{\mathit{mc}\mathit{,}{\mathit{mc}}^{\mathrm{\prime }}}\mathit{,}& & \end{array}$(17)where indicates a Monte-Carlo simulation of the corresponding map. The ensemble average is taken over two sets of independent realizations mc and mc′. Note that because of the way we have used pairs of Monte-Carlo simulations and data with independent CMB and noise realizations, the mean-field correction is zero for all of the terms above.

The term $\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}{\mathrm{|}}_{\mathit{N}\mathrm{1}}$ corrects for the “N⁽¹⁾” bias due to secondary contractions of the lensing trispectrum (Hu 2001; Kesden et al. 2003). It is only a large effect at L> 100, and so we calculate it using a flat-sky expression as $\begin{array}{ccc}{\mathrm{\Delta }{\mathit{C}}_{\mathrm{\left|}{L}\mathrm{\right|}\mathit{,x}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{1}}& \mathrm{=}& \frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathrm{\left|}{L}\mathrm{\right|}}^{\mathit{x\phi ,}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{\left(}\mathrm{2}\mathrm{\right)}}{\mathrm{ℛ}}_{\mathrm{\left|}{L}\mathrm{\right|}}^{\mathit{x\phi ,}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathrm{\left(}\mathrm{4}\mathrm{\right)}}}\mathrm{\int }\frac{{\mathrm{d}}^{\mathrm{2}}{{l}}_{\mathrm{1}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}\mathrm{\int }\frac{{\mathrm{d}}^{\mathrm{2}}{{l}}_{\mathrm{3}}}{\mathrm{(}\mathrm{2}\mathit{\pi }{\mathrm{)}}^{\mathrm{2}}}\\ & & \mathrm{\times }\hspace{0.17em}{\mathit{F}}_{\mathrm{\left|}{{l}}_{\mathrm{1}}\mathrm{\right|}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathit{F}}_{\mathrm{\left|}{{l}}_{\mathrm{2}}\mathrm{\right|}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}{\mathit{F}}_{\mathrm{\left|}{{l}}_{\mathrm{3}}\mathrm{\right|}}^{\mathrm{\left(}\mathrm{3}\mathrm{\right)}}{\mathit{F}}_{\mathrm{\left|}{{l}}_{\mathrm{4}}\mathrm{\right|}}^{\mathrm{\left(}\mathrm{4}\mathrm{\right)}}{\mathit{W}}^{\mathit{x}}\mathrm{\left(}{{l}}_{\mathrm{1}}\mathit{,}{{l}}_{\mathrm{2}}\mathrm{\right)}{\mathit{W}}^{\mathit{x}}\mathrm{\left(}{{l}}_{\mathrm{3}}\mathit{,}{{l}}_{\mathrm{4}}\mathrm{\right)}\\ & & \mathrm{\times }\hspace{0.17em}[{\mathit{C}}_{\mathrm{|}{{l}}_{\mathrm{1}}\mathrm{-}{{l}}_{\mathrm{3}}\mathrm{|}}^{\mathit{\phi \phi ,}\mathrm{fid}\mathit{.}}{\mathit{W}}^{\mathit{\phi }}\mathrm{(}\mathrm{-}{{l}}_{\mathrm{1}}\mathit{,}{{l}}_{\mathrm{3}}\mathrm{)}{\mathit{W}}^{\mathit{\phi }}\mathrm{(}\mathrm{-}{{l}}_{\mathrm{2}}\mathit{,}{{l}}_{\mathrm{4}}\mathrm{)}\\ & & \end{array}$(18)where l₁ + l₂ = l₃ + l₄ = L and ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi ,}\hspace{0.17em}\mathrm{fid}\mathit{.}}$ is a fiducial model for the lensing potential power spectrum. The W(l,l′) are flat-sky analogues of the full-sky weight functions. The flat-sky lensing weight function, for example, is ${\mathit{W}}^{\mathit{\phi }}\mathrm{\left(}{{l}}_{\mathrm{1}}\mathit{,}{{l}}_{\mathrm{2}}\mathrm{\right)}\mathrm{=}{\mathit{C}}_{\mathrm{\left|}{{l}}_{\mathrm{1}}\mathrm{\right|}}^{\mathit{TT}}{{l}}_{\mathrm{1}}\mathrm{·}{L}\mathrm{+}{\mathit{C}}_{\mathrm{\left|}{{l}}_{\mathrm{2}}\mathrm{\right|}}^{\mathit{TT}}{{l}}_{\mathrm{2}}\mathrm{·}{L}\mathit{.}$(19)The N⁽¹⁾ term is proportional to the lensing potential power spectrum, and so in principle it should be used to improve our constraints on ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ rather than subtracted as an additive bias. However the statistical power of this term is relatively small at Planck noise levels. From a Fisher matrix calculation, the trispectrum contractions which source the N⁽¹⁾ term are only detectable in the Planck data at 4σ significance, compared to approximately 25σ for the primary contractions. We choose simply to subtract the N⁽¹⁾ term from our power spectrum estimates. There is a small cosmological uncertainty in the N⁽¹⁾ correction due to uncertainty in the ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ power spectrum, which we discuss in Sect. 5.3.

The ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{PS}}$ term is a correction for the bias induced by the shot-noise of unresolved point sources (including SZ clusters), and will be discussed in more detail in Sect. 2.4.

Finally, the ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{MC}}$ term is a small correction that we obtain by estimating $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$ following the procedure above on a number of lensed CMB realizations, and then subtracting the input power spectrum. This term can be non-zero due to pseudo-C_ℓ leakage effects from masking, which we have not accounted for other than apodization, errors in our calculation of the N⁽¹⁾ term, or errors in the normalization at the power spectrum level. We will find that ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{MC}}$ is sufficiently small that in practice it does not matter whether we account for it as a renormalization or an additive offset, and we choose to treat it as an offset for simplicity.

We will ultimately characterize the uncertainty of $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$ by Monte-Carlo, however the following analytical expression is a useful approximation $\mathrm{Var}\mathrm{\left(}\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L,x}\end{array}\mathrm{\right)}\mathrm{\approx }{\mathit{V}}_{\mathit{L}}\mathrm{\equiv }\frac{\mathrm{1}}{{\mathit{f}}_{\mathrm{sky}\mathit{,}\mathrm{2}}}\frac{\mathrm{2}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}{\mathrm{[}{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{+}{{\mathit{N}}_{\mathit{L,x}}^{\mathit{\phi \phi }}}^{\mathrm{]}}}^{\mathrm{2}}$(20)where ${\mathit{N}}_{\mathit{L,x}}^{\mathit{\phi \phi }}$ is the reconstruction noise level. We take it to be $\begin{array}{ccc}{\mathit{N}}_{\mathit{L,x}}^{\mathit{\phi \phi }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}\frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{x\phi ,}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{\left(}\mathrm{2}\mathrm{\right)}}{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{x\phi ,}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathrm{\left(}\mathrm{4}\mathrm{\right)}}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{4}}{\left|{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{x}}\right|}^{\mathrm{2}}& & \\ \mathrm{\times }\mathrm{(}\mathit{Ĉ}\begin{array}{c}\mathit{TT,}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\\ {\mathit{\ell }}_{\mathrm{1}}\end{array}\mathit{Ĉ}\begin{array}{c}\mathit{TT,}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\mathrm{\left(}\mathrm{4}\mathrm{\right)}\\ {\mathit{\ell }}_{\mathrm{2}}\end{array}\mathrm{+}\mathit{Ĉ}\begin{array}{c}\mathit{TT,}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{\left(}\mathrm{4}\mathrm{\right)}\\ {\mathit{\ell }}_{\mathrm{1}}\end{array}\mathit{Ĉ}{\begin{array}{c}\mathit{TT,}\mathrm{\left(}\mathrm{2}\mathrm{\right)}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\\ {\mathit{\ell }}_{\mathrm{2}}\end{array}}^{\mathrm{)}}\mathit{,}& & \end{array}$(21)where $\mathit{Ĉ}\begin{array}{c}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{\left(}\mathit{j}\mathrm{\right)}\\ \mathit{L}\end{array}$ is the ensemble-average cross-spectrum between and , given by $\mathit{Ĉ}\begin{array}{c}\mathit{TT,}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{\left(}\mathit{j}\mathrm{\right)}\\ \mathit{L}\end{array}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{sky}\mathit{,}\mathrm{\left(}\mathit{i}\mathrm{\right)}\mathrm{\left(}\mathit{j}\mathrm{\right)}}^{-1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{\mathit{M}}\mathrm{⟨}\mathit{T̅}\begin{array}{c}\mathrm{\left(}\mathit{i}\mathrm{\right)}\\ \mathit{LM}\end{array}\mathit{T̅}{\begin{array}{c}\mathrm{\left(}\mathit{j}\mathrm{\right)}\mathrm{\ast }\\ \mathit{LM}\end{array}}^{\mathrm{⟩}}\mathit{,}$(22)where f_sky,(i)(j) is the fraction of sky common to both and . We will use this analytical estimate of the variance for weighting our ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ estimates in statistical analyses.

2.3. Lensing power spectrum likelihood

Based on our measurements of the lensing potential power spectrum, we construct a Gaussian likelihood based on bins in ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$. Our likelihood has the form $\mathrm{-}\mathrm{2}\ln {\mathrm{ℒ}}_{\mathit{\phi }}\mathrm{\left(}{\mathit{C}}^{\mathit{\phi \phi }}\mathrm{\right)}\mathrm{=}{\mathrm{ℬ}}_{\mathit{i}}^{\mathit{L}}\mathrm{(}\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}\mathrm{-}{{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}}^{\mathrm{)}}{\mathrm{\left[}{\mathrm{\Sigma }}^{-1}\mathrm{\right]}}^{\mathit{ij}}{\mathrm{ℬ}}_{\mathit{j}}^{{\mathit{L}}^{\mathrm{\prime }}}\mathrm{(}\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ {\mathit{L}}^{\mathrm{\prime }}\end{array}\mathrm{-}{{\mathit{C}}_{{\mathit{L}}^{\mathrm{\prime }}}^{\mathit{\phi \phi }}}^{\mathrm{)}}\mathit{,}$(23)where ℬ represents a binning function, Σ is the covariance matrix between bins, and sums are performed over paired upper/lower indices. As the shape and amplitude of ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ are constrained strongly in the concordance ΛCDM model, with which our results are broadly consistent, we choose a binning function designed to maximize our sensitivity to small departures from the fiducial ΛCDM expectation. This is given by ${\mathrm{ℬ}}_{\mathit{i}}^{\mathit{L}}\mathrm{=}\frac{{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi ,}\hspace{0.17em}\mathrm{fid}\mathit{.}}{\mathit{V}}_{\mathit{L}}^{-1}}{{\sum }_{{\mathit{L}}^{\mathrm{\prime }}\mathrm{=}{\mathit{L}}_{\min }^{\mathit{i}}}^{{\mathit{L}}_{\max }^{\mathit{i}}}\mathrm{(}{\mathit{C}}_{{\mathit{L}}^{\mathrm{\prime }}}^{\mathit{\phi \phi ,}\hspace{0.17em}\mathrm{fid}\mathit{.}}{\mathrm{)}}^{\mathrm{2}}{\mathit{V}}_{{\mathit{L}}^{\mathrm{\prime }}}^{-1}} \mathrm{if}{\mathit{L}}_{\min }^{\mathit{i}}\mathrm{\le }\mathit{L}\mathrm{\le }{\mathit{L}}_{\max }^{\mathit{i}}\mathit{,}$(24)where V_L is defined in Eq. (20)and ${\mathit{L}}_{\min }^{\mathit{i}}$ and ${\mathit{L}}_{\max }^{\mathit{i}}$ define the multipole range of the bin. Our binned results correspond to an estimate of the amplitude of a fiducial ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ power spectrum in a given multipole range, normalized to unity for the fiducial model. We denote these amplitude estimates explicitly for a given lens reconstruction $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$ as $\mathit{Â}\mathit{i}\mathrm{=}{\mathrm{ℬ}}_{\mathit{i}}^{\mathit{L}}\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}\mathit{.}$(25)In principle there are several reasons why a likelihood approach such as the one we have described above could fail; the usual issues with a pseudo-C_ℓ likelihood, such as the non-Gaussianity of the ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ are compounded by the fact that the φ estimates themselves are non-Gaussian, and derived from the temperature data itself, which means that the measurement uncertainties on the lens reconstruction are potentially correlated with those of the temperature power spectrum. In Appendix D, we validate the above approach to the likelihood by considering these issues in more detail.

2.4. Unresolved point-source correction

At the high multipoles (ℓ > 1000) that provide most of the modes which are useful for lens reconstruction, the contribution from unresolved (and therefore undetected and unmasked) extragalactic foregrounds becomes apparent in the Planck power spectra. For measurements of the lensing potential power spectrum, we should therefore be concerned with possible biases from the trispectrum non-Gaussianity of these unresolved sources. The work of Osborne et al. (2014) indicates that at Planck resolution and sensitivity, the contribution of unresolved foregrounds to lens reconstruction is small but potentially non-negligible, with the largest contribution coming from the shot noise trispectrum of unresolved radio and tSZ sources. Given uncertainties in modelling the non-Gaussianity of unresolved point sources, our approach is heavily data dependent. In addition to tests for the consistency of lens reconstructions at 143 and 217 GHz, we use two additional methods to measure and correct for point-source contamination in our analysis.

In our fiducial lensing power spectrum analysis, we measure the amplitude of the shot-noise contribution to the data trispectrum and correct the measured lensing spectrum accordingly. We measure the shot-noise amplitude using the power spectrum of a quadratic estimator designed to detect the “noise” due to sources, similar to the approach advocated in Munshi et al. (2011a). The trispectrum (or connected four-point function) for point-source shot noise is defined in position space as $\begin{array}{ccc}\mathrm{⟨}\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{1}\mathrm{\right)}\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{2}\mathrm{\right)}\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{3}\mathrm{\right)}\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{4}\mathrm{\right)}\mathrm{⟩}\mathrm{=}& & \\ {\mathit{S}}^{\mathrm{4}}\mathit{\delta }\mathrm{(}\hat{{n}}\mathrm{1}\mathrm{-}\hat{{n}}\mathrm{2}\mathrm{)}\mathit{\delta }\mathrm{(}\hat{{n}}\mathrm{2}\mathrm{-}\hat{{n}}\mathrm{3}\mathrm{)}\mathit{\delta }\mathrm{(}\hat{{n}}\mathrm{3}\mathrm{-}\hat{{n}}\mathrm{4}\mathrm{)}\mathit{,}& & \end{array}$(26)or in harmonic space: $\mathrm{⟨}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{4}}{\mathit{m}}_{\mathrm{4}}}\mathrm{⟩}\mathrm{=}{\mathit{S}}^{\mathrm{4}}\mathrm{\int }\mathrm{d}\hat{{n}}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{1}}}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{3}}{\mathit{m}}_{\mathrm{3}}}{\mathit{Y}}_{{\mathit{\ell }}_{\mathrm{4}}{\mathit{m}}_{\mathrm{4}}}\mathit{.}$(27)The quadratic estimator that we use to measure the amplitude S⁴ is denoted as , and is defined by the weight function ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{s}}\mathrm{=}\sqrt{\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& \mathit{L}\\ \mathrm{0}& \mathrm{0}& \mathrm{0}\end{array}\right)\mathit{.}$(28)We measure the amplitude of the shot-noise trispectrum S⁴ as $\begin{array}{c}\mathrm{􏽣}\\ {\mathit{S}}^{\mathrm{4}}\end{array}\mathrm{=}\left(\sum _{\mathit{L}}\mathit{Ĉ}\begin{array}{c}\mathit{ss}\\ \mathit{L}\end{array}{\mathit{S}}_{\mathit{L}}^{-1}\right)/\left(\sum _{{\mathit{L}}^{\mathrm{\prime }}}{\mathit{S}}_{{\mathit{L}}^{\mathrm{\prime }}}^{-1}\right)\mathit{,}$(29)where $\mathit{Ĉ}\begin{array}{c}\mathit{ss}\\ \mathit{L}\end{array}$ is defined analogously to Eq. (15)and the sums are taken over $\mathit{L,}{\mathit{L}}^{\mathrm{\prime }}\mathrm{\in }\mathrm{\left[}{\mathit{L}}_{\min }^{\mathit{s}}\mathit{,}{\mathit{L}}_{\max }^{\mathit{s}}\mathrm{\right]}$. The power spectrum weighting S_L is given by ${\mathit{S}}_{\mathit{L}}\mathrm{\equiv }\frac{\mathrm{1}}{{\mathit{f}}_{\mathrm{sky}\mathit{,}\mathrm{2}}}\frac{\mathrm{2}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{ss,}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{\left(}\mathrm{2}\mathrm{\right)}}{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{ss,}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathrm{\left(}\mathrm{4}\mathrm{\right)}}}\mathrm{·}$(30)In the limit that $\mathrm{\left[}{\mathit{L}}_{\min }^{\mathit{s}}\mathit{,}{\mathit{L}}_{\max }^{\mathit{s}}\mathrm{\right]}\mathrm{=}\mathrm{\left[}\mathrm{0}\mathit{,}\mathrm{\infty }\mathrm{\right]}$ this estimator is equivalent to measuring the shot-noise trispectrum using the 1-point kurtosis of the map. The advantage of the trispectrum-related power spectrum approach is that it allows us to separate out regions of the trispectrum measurement that are contaminated by lensing, as well as to look for evidence of clustering, which would appear as a deviation of the measured spectrum $\mathit{Ĉ}\begin{array}{c}\mathit{ss}\\ \mathit{L}\end{array}$ from the shape expected for unclustered sources. We then calculate and remove an estimated bias to the measured lensing potential given by $\mathrm{\Delta }\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L,x}\end{array}{\mathrm{|}}_{\mathrm{PS}}\mathrm{=}\begin{array}{c}\mathrm{􏽣}\\ {\mathit{S}}^{\mathrm{4}}\end{array}{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{xs,}\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{\left(}\mathrm{2}\mathrm{\right)}}{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{xs,}\mathrm{\left(}\mathrm{3}\mathrm{\right)}\mathrm{\left(}\mathrm{4}\mathrm{\right)}}\mathit{.}$(31)The shot-noise correction described above does not take into account the correlation of sources with the dark matter distribution. And hence the lensing potential. In Sect. 7.4 we therefore additionally perform tests using a “point source bias-hardened estimator”, constructed using the weight function of Eq. (28). This bias-hardened estimator has zero response to both the point-source shot-noise trispectrum (S⁴), as well as to the primary trispectrum contractions involving the correlation between point sources and the lensing potential (S²φ).

3. Data and cuts

Planck sky maps:

The majority of the results in this paper are based on the Planck nominal-mission frequency maps at 100, 143, and 217 GHz, built from the first 15.5 months of data. These are in HEALPix³ format, with resolution parameter N_side = 2048, corresponding to pixels with a typical width of . These have effective noise levels of approximately 105 μK arcmin at 100 GHz, 45 μK arcmin at 143 GHz, and 60 μK arcmin at 216 GHz. The beam widths shrink with frequency, and are 10′ at 100 GHz, 7′ at 143 GHz, and 5′ at 217 GHz. Our primary products – a map of the lensing potential and an associated power spectrum likelihood – are based on a minimum-variance (hereafter MV) combination of the 143 and 217 GHz maps. Although lensing may be detected at a significance of approximately 10σ at 100 GHz, the CMB modes used in this reconstruction are already sample-variance dominated at 143 and 217 GHz, and so adding it to the MV combination would lead to negligible improvement. In addition to these maps, we also use the 857 GHz Planck map as a dust template, which is projected out of both maps independently in our filtering procedure (described further in Appendix B). This projection is primarily intended to remove diffuse Galactic dust contamination, although it also removes a portion of the CIB fluctuations, which have a similar spectral index to that of Galactic dust over these frequency bands (Planck Collaboration XXIV 2011; Planck Collaboration XVIII 2011). As a simple approach, this template projection overlooks several potential difficulties, including variation of Galactic dust spectral indices across the sky as well as mismatch between the beams at 100/143/217 and 857 GHz, although we find it is adequate for our purposes. As lens reconstruction is most sensitive to small-scale modes, the coupling to large-scale Galactic foregrounds is relatively weak. In Sect. 7 we will also perform lens reconstruction using the more rigorously component-separated maps of Planck Collaboration XII (2014), finding good agreement with our baseline results.

In analyzing these maps, we use the fiducial beam transfer functions described in Planck Collaboration VII (2014) and Planck Collaboration (2013). There are uncertainties associated with these transfer functions, which we propagate to an uncertainty in the lensing estimator normalization in Sect. 5.2.

Galaxy mask:

We avoid the majority of Galactic foreground power using the temperature analysis masks described in Planck Collaboration (2013). These are constructed using a combination of the 30 GHz and 353 GHz maps, corrected for an estimate of the CMB contribution, smoothed to 5′ and thresholded until a desired sky fraction is obtained (Planck Collaboration XII 2014). For our baseline results, we use the 70% masks (which remove 30% of the sky), although in Sect. 7 we will show that we obtain consistent results with both larger and smaller masks. When computing the power spectrum, we additionally multiply our φ estimates by an apodized version of the Galaxy mask. Each pixel outside of the masked region is multiplied by an apodized weight varying between zero at the mask boundary and unity at a distance greater than 5deg from the closest masked pixel. We use a sinusoid weight function, similar to the one used in Namikawa et al. (2013) and Benoit-Levy et al. (2013). Note that due to this apodization, the effective sky fraction used in our power spectrum analysis is approximately 9% lower than the sky fraction of our reconstruction.

CO and extended-object masks:

We mask regions believed to be contaminated by carbon-monoxide (CO) emission lines at 100 and 217 GHz using the type 2 CO map described in Planck Collaboration XIII (2014). We mask all pixels above 60 μK in the map after degrading to N_side = 256, and then restore by hand all isolated pixels that are removed by this cut but do not appear to be extended CO regions. This low-resolution mask is prograded back to N_side = 2048 and smoothed at 20′, slightly larger than the map resolution of 15′. We additionally exclude extended nearby objects (the two Magellanic Clouds, M31, M33, and M81) by cutting out disks centred on their locations. The radii of these disks range from 250′ for the Large Magellanic Cloud to 30′ for M81.

Point-source masks:

We remove detected point sources (or otherwise compact objects) using a mask constructed from a combination of sources identified in the Planck Early Release Compact Source Catalogue (ERCSC; Planck Collaboration Int. VII 2013), the SZ clusters from the Planck Cluster Catalogue (PCC; Planck Collaboration XXIX 2014), and the Planck Catalogue of Compact Sources (PCCS; Planck Collaboration XXVIII 2014). We produce individual masks for 100, 143, and 217 GHz, cutting out detected sources with disks having radii of either 3σ or 5σ depending on their flux level, where σ describes the Gaussian beam-width of the given channel and is taken to be 4.1′ at 100 GHz, 3.1′ at 143 GHz, and 2.1′ at 217 GHz. SZ clusters are not masked at 217 GHz. For the PCCS we mask all sources with S/N ≥ 5 in a given band, or S/N ≥ 10 in either adjacent frequency band, using a 3σ disk, extending to a 5σ disk for any sources with S/N ≥ 10 in the target band. For the ERCSC, we make a 3σ cut for every source detected at the target frequency, as well as for every source that is detected with a flux greater than 1 Jy at 70 GHz, 0.7 Jy at 100 GHz, 0.4 Jy at 143 GHz, 0.5 Jy at 217 GHz, or 0.9 Jy at 353 GHz. For sources that satisfy the brightness criterion at the target frequency, we also enlarge the cut to 5σ. For SZ clusters we make a cut between 3σ and 5σ depending on the cluster angular size θ₅₀₀ for all sources in the MMF1 catalogue with S/N ≥ 5. We further mask cool cores (CC) using the Early Cold Cores catalogue (Planck Collaboration VII 2011) in a similar fashion, using a maximum 15′ excluding radius.

In Sect. 7.2 we will test our sensitivity to this SNR cut, varying it from S/N ≥ 5 to S/N ≥ 4.2 and S/N ≥ 10.

4. Simulations

We require simulations of the data described in Sect. 3 for several aspects of our analysis: to determine the lensing-mean field φ^MF, to determine the ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{0}}$ and ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{MC}}$ correction terms of Eq. (15), and to determine our measurement error bars, as well as to validate our reconstruction methodology. Here we outline the basic algorithm we use to generate these simulations.

• 1.

  Simulate the unlensed CMB and lensing potential. Using the fiducial ΛCDM theoretical unlensed power spectra for the cosmology given at end of Sect. 1, we generate realizations ${\mathit{T}}_{\mathit{\ell m}}^{\mathrm{unl}}$ and φ_LM of the unlensed temperature and lensing potential with appropriate correlations. By default we simulate these fields up to ℓ_max = 2560. In addition to these stochastic components, we also include a dipolar lensing mode in our simulations, to account for the Doppler aberration due to our motion with respect to the CMB (Kamionkowski & Knox 2003; Challinor & van Leeuwen 2002): ${\mathit{\phi }}_{\mathit{\beta }}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{\approx }{\beta }\mathrm{·}\hat{{n}}\mathit{,}$(32)where β is in the direction of Galactic coordinates (l,b) = (264.4,48.4), with amplitude | β | = 0.00123, corresponding to a velocity of 369 km s^-1 (Hinshaw et al. 2009).

• 2.

  Lens the temperature field. We use a fast spherical harmonic transform to compute the temperature and deflection fields $\begin{array}{ccc}{\mathit{T}}^{\mathrm{unl}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{\ell m}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\mathit{T}}_{\mathit{\ell m}}^{\mathrm{unl}}\mathit{,}\\ {d}\mathrm{\left(}\hat{{n}}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{M}}\mathrm{\nabla }{\mathit{Y}}_{\mathit{LM}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\mathit{\phi }}_{\mathit{LM}}\mathit{.}\end{array}$(33)This computation is performed on an equicylindrical pixelization (ECP) grid with (N_θ,N_φ) = (16 384,32 768) equally-spaced points in (θ,φ). We determine the lensed temperature on this ECP grid as $\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}{\mathit{T}}^{\mathrm{unl}}\mathrm{(}\hat{{n}}\mathrm{+}{d}\mathrm{(}\hat{{n}}\mathrm{\left)}\mathrm{\right)}\mathit{.}$(34)This equation is evaluated following the same algorithm as the LENSPix code (Lewis 2005), with 2D cubic Lagrange interpolation on the ECP grid, which has conveniently uniform pixel spacing. We then obtain the lensed temperature multipoles ${\mathit{T}}_{\mathit{\ell m}}^{\mathit{\phi }}$ by a reverse harmonic transform. The average power spectrum of these lensed simulations agrees with the implementation in CAMB (Challinor & Lewis 2005) to better than 0.1% at l< 2048. We have also used the algorithm of Basak et al. (2009), which gives consistent results. In addition to the lensing dipole, our motion with respect to the CMB also causes a small modulation of the observed fluctuations, which we determine from the lensed temperature field as ${\mathit{T}}^{\mathrm{mod}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}{\mathit{f}}_{\mathit{\nu }}{\beta }\mathrm{·}\hat{{n}}\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{,}$(35)where f_ν is a frequency dependent boost factor that we approximate as 1.5, 2, and 3 at 100, 143, and 217 GHz, respectively. There is an extended discussion of this factor in Planck Collaboration XXVII (2014).

• 3.

  Generate a Gaussian realization of extragalactic foreground power ${\mathit{T}}_{\mathit{\ell m}}^{\mathrm{fg}}$. We simulate Gaussian foreground power at 100, 143, and 217 GHz such that the auto- and cross-spectra of the foreground power at these frequencies has the form ${\mathit{C}}_{\mathit{\ell }}^{\mathit{A}\mathrm{\times }\mathit{B}}\mathrm{=}{\mathit{A}}_{\mathit{A}\mathrm{\times }\mathit{B}}^{\mathrm{PS}}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{PS}}\mathrm{+}{\mathit{A}}_{\mathit{A}\mathrm{\times }\mathit{B}}^{\mathrm{CIB}}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}}\mathit{,}$(36)where A and B label frequency bands. The template power spectra are based on the modelling of Planck Collaboration XV (2014); Planck Collaboration XVI (2014), and are given by $\begin{array}{ccc}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{PS}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{{\mathrm{3000}}^{\mathrm{2}}}\mathit{,}{\mathit{C}}_{\mathit{\ell }}^{\mathrm{CIB}}\mathrm{=}\frac{\mathrm{2}\mathit{\pi }}{\mathit{\ell }\mathrm{(}\mathit{\ell }\mathrm{+}\mathrm{1}\mathrm{)}}{\left(\frac{\mathit{\ell }}{\mathrm{3000}}\right)}^{\mathrm{0.8}}\mathrm{·}& & \end{array}$(37)We use the parameters $\begin{array}{cc}{\mathit{A}}_{\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{100}}^{\mathrm{PS}}& \mathrm{=}\mathrm{208}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}\\ \mathrm{=}\mathrm{72}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}& {\mathit{A}}_{\mathrm{100}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}}^{\mathrm{PS}}\\ {\mathit{A}}_{\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{143}}^{\mathrm{PS}}& \mathrm{=}\mathrm{64}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}\\ \mathrm{=}\mathrm{43}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}& {\mathit{A}}_{\mathrm{217}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}}^{\mathrm{PS}}\\ {\mathit{A}}_{\mathrm{143}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{143}}^{\mathrm{CIB}}& \mathrm{=}\mathrm{4}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}\\ \mathrm{=}\mathrm{14}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{2}}& {\mathit{A}}_{\mathrm{217}\hspace{0.17em}\mathrm{\times }\hspace{0.17em}\mathrm{217}}^{\mathrm{CIB}}\end{array}$The parameters above were obtained by fitting the auto- and cross-spectra of our inverse-variance-filtered temperature maps to the templates above. Our point source masks differ from those used in Planck Collaboration XVI (2014), and so these numbers cannot be directly compared to the values there. The purpose of including a Gaussian foreground component in our simulations is primarily so that our measurement error bars will include the additional scatter due to the foreground power (we will see in Sect. 5 that foreground power constitutes a few percent of our total error budget). Our lensing power spectrum results are otherwise unaffected by this component of our simulations. Our lensing estimates do include an overall correction for point-source shot-noise, although we do not include the non-Gaussian covariance of the shot-noise in our error budget, given that point sources are a small fraction of the total power at the multipoles probed by Planck. We test this point-source treatment further in Sect. 7.2.

• 4.

  Convolve T^φ + T^fg + T^mod with the (asymmetric) instrumental beam using Eq. (C.7). This results in a HEALPix map of the beam-convolved CMB ${\mathit{T}}^{\mathit{b}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$, from which we extract beam-convolved multipoles ${\mathit{T}}_{\mathit{\ell m}}^{\mathit{b}}$.

• 5.

  Make a final map of the observed sky $\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ with the HEALPix N_side = 2048 pixelization, with $\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\ell m}}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}{\mathit{B}}_{\mathit{\ell }}^{\mathrm{pix}}{\mathit{T}}_{\mathit{\ell m}}\mathrm{+}\mathit{n}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{,}$(38)where ${\mathit{B}}_{\mathit{\ell }}^{\mathrm{pix}}$ is the HEALPix N_side = 2048 pixel window function, which accounts for the effective smearing of the beam due to the binning of hits into pixels. Although this smearing is an inadequate description of the pixelization-smearing effect for any individual pixel, it is a good approximation when averaged over all pixels on the sky. We will discuss pixelization effects further in Sect. C.4. The quantity $\mathit{n}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ is a simulated realization of the instrumental noise. We use noise realizations from the FFP6 simulation set (Planck Collaboration I 2014). In order to make the power spectrum of our simulations more closely match the power spectrum of the data, we find that it is necessary to include an additional white noise component with an amplitude of 20, 10, and 10 μK arcmin at 100, 143, and 217 GHz respectively. Indeed, the FFP6 simulation are produced on the basis of the half-ring estimates of the noise which are known to underestimate the map noise (Planck Collaboration VI 2014).

The simulations described above are used throughout this work both to debias our lensing estimates, as well as to characterize their uncertainty. Here we briefly use them to validate the normalization of our lensing estimates. There are two relevant normalizations: the map-level normalization which describes how our reconstruction traces the input φ realization, as well as the analogous normalization for the lensing potential power spectrum. Throughout this work, we use the analytical normalization estimates of Eq. (12)to normalize our map estimates. In Fig. 3 we plot the average cross-correlation between the input and reconstructed φ_LM for our minimum-variance reconstruction, divided by a simple f_sky factor to account for missing power in the mask. We can see that our analytical map-level normalization is accurate at the 1% level. This is small enough that we do not further renormalize our φ map. At the power spectrum level, we find similar accuracy of our normalization in the 40 <L< 400 region which forms the basis of our fiducial lensing likelihood. The power spectrum estimates in Fig. 3 are obtained by averaging the first line of Eq. (15)over a set of lensed but otherwise Gaussian simulations. This 1% error in the power spectrum normalization is accounted for in our likelihood by the ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{MC}}$ term. Given the good agreement of the ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ power spectrum used in our simulations with that measured from data, it does not matter whether we incorporate the Monte-Carlo correction as a multiplicative or additive term. We have chosen to incorporate it as additive for simplicity. At higher multipoles, there is additional scatter in our power spectrum normalization test (up to 10% for the plotted bins). It is difficult to know from this test whether it is due to actual scatter in our power spectrum normalization (which seems unlikely, given that the map-level normalization is quite accurate), a small mis-estimation of the N⁽¹⁾ term (which begins to dominate over ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ above L = 500), or simply pseudo-C_ℓ leakage issues. In any case, it is still small, and is absorbed into our ΔC_L|_MC term when estimating power spectra at these multipoles. Although we have only presented results for our minimum-variance reconstruction here, those for individual frequency results are similar.

As a visual illustration, and a preview of our data results in Sect. 6, in Fig. 4, we show a simulated lens reconstruction as well as the input φ map. This gives a visual impression of the signal-to-noise in our lens reconstruction.

5. Error budget

In this section, we describe the measurement and systematic error budget for our estimation of the lensing potential power spectrum. This is broken down into three sections; in Sect. 5.1 we describe our measurement (or “statistical”) error bars, which are due to the fact that we have only a single noisy sky with a finite number of modes to observe. In Sect. 5.2 we consider uncertainty in the instrumental beam transfer function, which we will see propagates to a normalization uncertainty for our lensing estimates. Finally, in Sect. 5.3 we discuss the effect of cosmological uncertainty; possible errors in the fiducial model for ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$ result in a normalization uncertainty for our lensing estimates, and uncertainties in the fiducial ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ power spectrum lead to uncertainties in the ${\mathit{N}}_{\mathit{L}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ correction. As a guide to the relative size and scale dependence of these terms, in Fig. 5 we summarize the error budget for our fiducial minimum-variance lens reconstruction, based on 143 and 217 GHz. Individual frequency bands, as well as 100 GHz are qualitatively similar.

5.1. Measurement

Although our measurement uncertainties are ultimately assigned by Monte Carlo, we can use the analytical expression of Eq. (21)to gain intuition for how they are sourced by various components. Our simple model of the sky after masking and dust cleaning is that it consists of three uncorrelated signals: CMB, instrumental noise, and unresolved isotropic foreground power. The noise variance of the lens reconstruction in Eq. (21)involves two power spectra, and so we can think of the noise contribution as the sum of six possible terms involving pairs of the CMB, noise, and foreground power spectra. In Fig. 6 we combine these contributions into three representative contributions to the reconstruction noise: “pure CMB” in which both spectra are due to CMB fluctuations; the “noise” contribution in which either both spectra are those for noise power, or one is noise and one is CMB; and, finally, the “foreground” contribution in which either one or both of the spectra are due to unresolved foreground power. We can see that for most reconstruction multipoles, the pure CMB contribution constitutes the largest part of the reconstruction noise, followed by noise. The unresolved foreground power is a fairly small contribution to our measurement error. Note that the dominant terms for both the “noise” and “foreground” contributions are the ones in which one of the spectra is a CMB fluctuation. For this reason, we will focus less on the use of cross-spectra to avoid noise biases than is done for the usual CMB power spectra (Planck Collaboration XV 2014), although we will perform consistency tests using cross-spectra of data to avoid noise biases. Note that our realization-dependent method for removing the disconnected noise bias (Eq. (17)) means that the majority of this contribution is estimated directly from the data itself, reducing our sensitivity to uncertainty in the noise and foreground power.

5.2. Beam transfer function

Errors in the effective beam transfer function appear as an error in the normalization of our lensing estimates. For simplicity here we will describe the case for a single standard quadratic lensing estimator that uses the same map for both of its inputs, although when dealing with combinations of channels for our actual results we account for differences in the beam transfer function and errors between “legs” of the estimator (or pair of estimators, in the case of the lensing potential power spectrum).

We model the effect on from a discrepancy between the “measured” and “true” effective beam transfer function as $\mathrm{⟨}\mathit{\phi ̂}\mathit{LM}\mathrm{⟩}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi \phi }}}\left[\frac{\mathrm{1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{2}}{\left|{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}\right|}^{\mathrm{2}}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{2}}}\frac{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{true}}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{meas}}}\frac{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{true}}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{meas}}}\right]{\mathit{\phi }}_{\mathit{LM}}\mathit{,}$(39)where the average is over CMB fluctuations for fixed lenses. The statistics of B^true/B^meas are discussed in Planck Collaboration VII (2014), where it is found that the beam uncertainty may be described by a small number N_eig of eigenmodes (typically five), $\ln \left(\frac{{\mathit{B}}_{\mathit{\ell }}^{\mathrm{meas}}}{{\mathit{B}}_{\mathit{\ell }}^{\mathrm{true}}}\right)\mathrm{=}\sum _{\mathit{i}\mathrm{=}\mathrm{1}}^{{\mathit{N}}_{\mathrm{eig}}}{\mathit{g}}_{\mathit{i}}{\mathit{E}}_{\mathit{\ell }}^{{\mathit{b}}_{\mathit{i}}}\mathit{.}$(40)Here, ${\mathit{E}}_{\mathit{\ell }}^{{\mathit{b}}_{\mathit{i}}}$ are the uncertainty eigenmodes (indexed by i) for the beam determination, and g_i are independent Gaussian random variables with unit variance. The frequency-band eigenmodes are small for 100−217 GHz (no more than 1% at l< 2048), and so we can equate $\frac{{\mathit{B}}_{\mathit{\ell }}^{\mathrm{true}}}{{\mathit{B}}_{\mathit{\ell }}^{\mathrm{meas}}}\mathrm{\approx }\mathrm{1}\mathrm{-}\ln \left(\frac{{\mathit{B}}_{\mathit{\ell }}^{\mathrm{meas}}}{{\mathit{B}}_{\mathit{\ell }}^{\mathrm{true}}}\right)\mathrm{·}$(41)We can also ignore terms of order E² and propagate the independent beam uncertainties directly to a corresponding set of independent eigenmodes for φ. These eigenmodes are given as ${\mathit{E}}_{\mathit{L}}^{{\mathit{\phi }}_{\mathit{i}}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi \phi }}}\left[\frac{\mathrm{1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{2}}{\left|{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}\right|}^{\mathrm{2}}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{2}}}\sum _{\mathit{i}}\mathrm{(}{\mathit{E}}_{{\mathit{\ell }}_{\mathrm{1}}}^{{\mathit{b}}_{\mathit{i}}}\mathrm{+}{{\mathit{E}}_{{\mathit{\ell }}_{\mathrm{2}}}^{{\mathit{b}}_{\mathit{i}}}}^{\mathrm{)}}\right]\mathit{,}$(42)with the model that $\frac{\mathrm{⟨}\mathit{\phi ̂}\mathit{LM}\mathrm{⟩}}{{\mathit{\phi }}_{\mathit{LM}}}\mathrm{=}\mathrm{1}\mathrm{+}\sum _{\mathit{i}\mathrm{=}\mathrm{0}}^{{\mathit{N}}_{\mathrm{eig}}}{\mathit{g}}_{\mathit{i}}{\mathit{E}}_{\mathit{L}}^{{\mathit{\phi }}_{\mathit{i}}}\mathit{.}$(43)

In Fig. 7 we plot the beam uncertainty eigenmodes as well as the corresponding normalization eigenmodes. The normalization eigenmodes are flatter than the beam uncertainty eigenmodes, and do not converge to zero at low-L, reflecting the fact that at all scales the lensing estimate takes most of its weight from modes of ℓ₁ and ℓ₂ close to the noise/beam cutoff (at ℓ ≈ 1000 for 100 GHz and ℓ ≈ 1500 for 143 and 217 GHz). As can be seen in Fig. 7, the beam normalization eigenmode is essentially twice the value of the beam uncertainty eigenmode on these scales.

Owing to the small size of the eigenmodes, the propagation to a fractional uncertainty on ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ is accomplished essentially by applying a factor of two to the eigenmodes of Eq. (42). The normalization uncertainty due to beams for the lensing potential power spectrum is therefore approximately 1% at 100 GHz and 0.5% at 143 and 217 GHz. When computing the beam uncertainty for our minimum-variance result, we take the 143 and 217 GHz beam eigenmodes to be uncorrelated.

5.3. Cosmological

Similarly to the case with the beam transfer function of the previous subsection, uncertainty in the cosmological model may also be propagated to an uncertainty in the estimator normalization. There are two aspects in which cosmological uncertainty can enter our analysis. First, when estimating the ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{1}}$ contribution of Eq. (15)we assume a fiducial ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$, obtained from the parameters of the fiducial ΛCDM cosmology described in Sect. 1. We take a random subset of samples from Planck MCMC chains for the ΛCDM cosmology and calculate the corresponding ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{1}}$ bias. Compared to our fiducial correction, we see small-scale fluctuations with an amplitude of 0.2% and a typical scale of ΔL = 150, superimposed on a larger overall amplitude scatter with a standard deviation of approximately 6%. This result is reasonable: in the ΛCDM model the PlanckTT measurements have the power to detect the smoothing effects of gravitational lensing on the acoustic peaks at a significance level of 10σ, corresponding to a 10% constraint on the overall amplitude of ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$, consistent with the result from analysis of individual samples above. To account for cosmological uncertainty in the ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{1}}$ correction, we therefore assign it a 10% overall amplitude uncertainty. This correction is small for the L range considered in our fiducial likelihood, and so this uncertainty is essentially negligible compared to our statistical error bars (although we include it in our error budget nevertheless).

The second point at which cosmological uncertainty enters our analysis is in the normalization of our estimates. In Eq. (12), as well as in our tests of the normalization in Sect. 4, we have assumed that the statistically-anisotropic covariance induced by lensing is equal to that assumed when weighting the estimator, given by Eq. (6). If the true cosmological power spectrum differs from the fiducial one used in Eq. (6)then we will misestimate φ as $\frac{\mathrm{⟨}\mathit{\phi ̂}\mathit{LM}\mathrm{⟩}}{{\mathit{\phi }}_{\mathit{LM}}}\mathrm{=}\mathrm{1}\mathrm{+}{\mathrm{\Delta }}_{\mathit{L}}^{\mathit{TT}}\mathrm{=}\frac{\mathrm{1}}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi \phi }}}\left[\frac{\mathrm{1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{2}}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi ,}\hspace{0.17em}\mathrm{true}}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{1}}}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{2}}}\right]\mathit{,}$(44)where ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi ,}\mathrm{true}}$ is given by Eq. (6), with the modification that it is evaluated with the “true” ensemble-average lensed CMB power spectrum rather than the fiducial one ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$. We have chosen to account for this normalization uncertainty coherently with the temperature likelihood. Given a proposed ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$, ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ pair, we renormalize ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{\to }\mathrm{(}\mathrm{1}\mathrm{+}{\mathrm{\Delta }}_{\mathit{L}}^{\mathit{TT}}{\mathrm{)}}^{\mathrm{2}}{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ in our likelihood (Eq. (23)) so that it is directly comparable to our measurement.

6. Results

In this section, we present the primary results of this paper; a measurement of the CMB lensing potential over approximately 70% of the sky and its associated power spectrum estimate $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$. Our lensing potential map and power spectrum have several science implications, which we will discuss in proceeding subsections.

• 1.

  The lensing power spectrum probes the matter power spectrum integrated back to the last-scattering surface, and provides additional low-redshift leverage for CMB-alone parameter constraints with a broad maximum at z ~ 2. Highlights of the additional information that the Planck lensing likelihood contributes to the parameter constraints of Planck Collaboration XVI (2014) are presented in Sect. 6.1.

• 2.

  Our lensing potential map correlates with the CMB temperature via the late-time integrated Sachs-Wolfe (ISW) effect sourced by dark energy. We measure the cross-spectrum between our lens reconstruction and the low-ℓ temperature multipoles in Sect. 6.2, which we find is present at the level predicted in ΛCDM, and discrepant with the null hypothesis of no correlation at approximately 2.5σ.

• 3.

  Our lensing potential map correlates with the fluctuations of the cosmic infrared background, which also trace large-scale structure and have good redshift overlap with the CMB lensing potential. We do not discuss this correlation here, but instead refer to the detailed analysis of Planck Collaboration XVIII (2014).

• 4.

  Finally, our lensing potential map correlates with external tracers of large-scale structure. In Sect. 6.3 we present cross-correlations with several representative galaxy catalogues.

Our primary results are based on a minimum-variance (MV) combination of the 143 and 217 GHz Planck maps, although as an initial consistency check we will also present most of our results in this section for the individual frequency maps as well. In Sect. 7 we will present a more extensive suite of consistency tests on our power spectrum estimates and lensing map. In particular in Sect. 7.1, we will compare with results obtained on different masks and on the foreground-cleaned Planck maps described in Planck Collaboration XII (2014).

Our MV lensing map is shown in Fig. 8, where we plot the Wiener-filtered reconstruction of the lensing potential. At the noise levels of our reconstruction, approximately half of the modes in this map are noise. At a visual level, the lensing map of Fig. 8 has comparable signal-to-noise to the COBE-DMR maps of the CMB temperature anisotropy (Bennett et al. 1996). Future measurements, including the Planck full mission release which contains approximately twice the amount of data used here, will improve on this first full-sky map of the CMB lensing potential. As is illustrated in the simulated reconstruction of Fig. 4, there will be clear visual correlations between this map and future measurements. For comparison with the full MV result, in Fig. 9 we also plot the individual 143  and 217 GHz reconstructions for the southern hemisphere.

In Fig. 10 we plot the power spectra of our individual 100, 143, and 217 GHz reconstructions as well as the minimum-variance reconstruction. The agreement of all four spectra is striking.

Overall, our power spectrum measurement is reasonably consistent with the ΛCDM prediction, given our measurement error bars. Dividing the L ∈ [ 1,2048 ] multipole range into bins of ΔL = 64 and binning uniformly in $\mathrm{\left[}\mathit{L}\mathrm{\right(}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$, we obtain a reduced χ² for the difference between our power spectrum estimate and the model of 40.7 with 32 degrees of freedom. The associated probability to exceed is 14%. On a detailed level, there are some discrepancies between the shape and amplitude of our power spectrum and the fiducial model however. Our likelihood is based on the multipole range 40 ≤ L ≤ 400, which captures 90% of the forecasted ensemble-average signal-to-noise for an amplitude constraint on ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ in our fiducial cosmological model (see Fig. 2). This range was chosen as the region of our spectrum least likely to be contaminated by systematic effects (primarily uncertainties in the mean-field corrections at low-L, and uncertainties in the Gaussian and point-source bias corrections at high-L). Estimating an average amplitude for the fiducial lensing power spectrum for a single bin over this multipole range using Eq. (25)we find an amplitude of Â_{40 → 400} = 0.94 ± 0.04 relative to the fiducial model (which has A = 1). The power in this region is consistent with the fiducial model, although 1.5σ low (the corresponding probability-to-exceed for the χ² of this difference is 15%). The low- and high-L extent of our likelihood were deliberately chosen to have enough expected lensing signal to enable a 10σ detection of lensing on either side, bookending our likelihood with two additional consistency tests. On the low-L side, we have a good agreement with the expected power. As will be discussed in Sect. 7.4, our measurement at L< 10 fails some consistency tests at a level comparable to the expected signal. The L< 10 modes, which we suspect are somewhat contaminated by errors in the mean-field subtraction, are nevertheless consistent with the fiducial expectation, as can be seen in Fig. 10; we measure Â_{1 → 10} = 0.44 ± 0.54. Extending to the lower limit of our likelihood, with a single bin from 10 ≤ L ≤ 40 we measure Â_{10 → 40} = 1.02 ± 0.12. On the high-L side of our fiducial likelihood, there is tension however. Extending from the final likelihood multipole at L = 400 to the maximum multipole of our reconstruction, we find Â_{400 → 2048} = 0.68 ± 0.13, which is in tension with A = 1 at a level of just over 2.4σ. The relatively low power in our reconstruction is driven by a dip relative to the ΛCDM model spectrum between 500 <L< 750, as can be seen in Fig. 10. We show this feature more clearly in the residual plot of Fig. 11. This deficit of power is in turn driven by the 143 GHz data. For an estimate of the power spectrum using only 143 GHz, we measure $\mathit{Â}\begin{array}{c}\mathrm{143}\\ \mathrm{400}\mathrm{\to }\mathrm{2048}\end{array}\mathrm{=}\mathrm{0.37}\mathrm{\pm }\mathrm{0.18}$. The 217 GHz reconstruction is more consistent with the model, having $\mathit{Â}\begin{array}{c}\mathrm{217}\\ \mathrm{400}\mathrm{\to }\mathrm{2048}\end{array}\mathrm{=}\mathrm{0.82}\mathrm{\pm }\mathrm{0.17}$. These two measurements are in tension; we have $\mathit{Â}\begin{array}{c}\mathrm{217}\mathrm{-}\mathrm{143}\\ \mathrm{400}\mathrm{\to }\mathrm{2048}\end{array}\mathrm{=}\mathrm{0.45}\mathrm{\pm }\mathrm{0.18}$, which is a 2.5σ discrepancy. The error bar on this difference accounts for the expected correlation between the two channels due to the fact that they see the same CMB sky. A larger set of consistency tests will be presented in Sect. 7. We note for now that the bins from 40 <L< 400 used in our likelihood pass all consistency tests, and show better agreement between 143 and 217 GHz. Although L< 40 and L> 400 are not included in our nominal likelihood, when discussing the use of the lensing likelihood for cosmological parameter constraints in the following section we will perform additional cross-checks using these bins to ascertain whether they would have any significant implications for cosmology.

In addition to the Planck power spectrum measurements, in Fig. 11 we have overplotted the ACT and SPT measurements of the lensing potential power spectrum (Das et al. 2013; van Engelen et al. 2012). It is clear that all are very consistent. The Planck data provides the largest signal-to-noise of these measurements; as we have already discussed the 40 <L< 400 lensing likelihood provides a 4% constraint on the amplitude of the lensing potential power spectrum, while the constraint from current ACT and SPT measurements are 32% and 16% respectively. These measurements are nevertheless quite complementary. As a function of angular scale, the full-sky Planck power spectrum estimate has the smallest uncertainty per multipole of all three experiments at L< 500, at which point the additional small-scale modes up to ℓ_max = 3000 used in the SPT lensing analysis lead to smaller error bars. The good agreement in these estimates of ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ is reassuring; in addition to the fact that the experiments and analyses are completely independent, these measurements are sourced from fairly independent angular scales in the temperature map, with ℓ ≲ 1600 in the case of Planck, ℓ < 2300 in the case of ACT, and ℓ < 3000 in the case of SPT. Cross-correlation of the Planck lensing map with these independent measures of the lensing potential will provide an additional cross-check on their consistency, however at the power spectrum level they are already in good agreement.

6.1. Parameters

Weak gravitational lensing of the CMB provides sensitivity to cosmological parameters affecting the late-time growth of structure that are otherwise degenerate in the primary CMB anisotropies imprinted around recombination. Examples include the dark energy density in models with spatial curvature and the mass of neutrinos that are light enough (m_ν< 0.5 eV) still to have been relativistic at recombination.

To connect our measurement of the lensing power spectrum to parameters, we construct a lensing likelihood nominally based on the multipole range 40 ≤ L ≤ 400, cut into eight equal-width bins with ΔL = 45 to maintain parameter leverage from shape information in addition to our overall amplitude constraint. In Table 1 we present bandpowers for these eight bins using the individual 100, 143, and 217 GHz reconstructions as well as the MV reconstruction that is the basis for our nominal likelihood. The bandpower estimates and their uncertainties are broken down into constituent parts as discussed in Sect. 2. Based on these bandpowers, we form a likelihood following Eq. (23). The measurement errors on each bin are measured by Monte-Carlo using 1000 simulations, and the bins are sufficiently wide that we can neglect any small covariance between them (this is discussed further in Appendix D). We analytically marginalize over uncertainties that are correlated between bins, including them in the measurement covariance matrix. This includes beam transfer function uncertainties (as described in Sect. 5.2), uncertainties in the point source correction (Sect. 7.2) and uncertainty in the N⁽¹⁾ correction.

As the lensing likelihood is always used in conjunction with the PlanckTT power spectrum likelihood, we coherently account for uncertainty in ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$ by renormalizing our lensing potential measurement for each sample, as described in Sect. 5.3.

The lensing likelihood is combined with the main PlanckTT likelihood (Planck Collaboration XV 2014) – constructed from the temperature (pseudo) cross-spectra between detector sets at intermediate and high multipoles, and an exact approach for Gaussian temperature anisotropies at low multipoles – in Planck Collaboration XVI (2014) to derive parameter constraints for the six-parameter ΛCDM model and well-motivated extensions. Lensing also affects the power spectrum, or 2-point function, of the CMB anisotropies, and this effect is accounted for routinely in all Planck results. On the angular scales relevant for Planck, the main effect is a smoothing of the acoustic peaks and this is detected at around 10σ in the Planck temperature power spectrum (Planck Collaboration XVI 2014). The information about ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ that is contained in the lensed temperature power spectrum for multipoles ℓ ≲ 3000 is limited to the amplitude of a single eigenmode (Smith et al. 2006). In extensions of ΛCDM with a single additional late-time parameter, lensing of the power spectrum itself can therefore break the geometric degeneracy (Stompor & Efstathiou 1999; Sherwin et al. 2011; van Engelen et al. 2012; Planck Collaboration XVI 2014). As discussed in Appendix D and Schmittfull et al. (2013), cosmic variance of the lenses produces weak correlations between the CMB 2-point function and our estimates of ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$, but they are small enough that ignoring the correlations in combining the two likelihoods should produce only sub-percent underestimates of the errors in physical cosmological parameters.

In the following, we illustrate the additional constraining power of our ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ measurements in ΛCDM models and one-parameter extensions, highlighting those results from Planck Collaboration XVI (2014) where the lensing likelihood is influential.

6.1.1. Six-parameter ΛCDM model

In the six-parameter ΛCDM model, the matter densities, Hubble constant and spectral index of the primordial curvature perturbations are tightly constrained by the Planck temperature power spectrum alone. However, in the absence of lensing the amplitude A_s of the primordial power spectrum and the reionization optical depth τ are degenerate, with only the combination A_se^{− 2τ}, which directly controls the amplitude of the anisotropy power spectrum on intermediate and small scales, being well determined. This degeneracy is broken by large-angle polarization since the power from scattering at reionization depends on the combination A_sτ². In this first release of Planck data, we use the WMAP nine-year polarization maps (Bennett et al. 2013) in combination with Planck temperature data. With this data combination, ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ is rather tightly constrained in the ΛCDM model (see Fig. 12) and the direct measurements reported here provide a non-trivial consistency test of the model.

The eight ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ bandpowers used in the lensing likelihood are compared to the expected spectrum in Fig. 12 (upper-left panel). For the latter, we have used parameter values determined from the main Planck likelihood in combination with WMAP polarization (hereafter denoted WP) and small-scale power spectrum measurements (hereafter highL) from ACT (Das et al. 2013) and SPT (Reichardt et al. 2012)⁴. In this plot, we have renormalized the measurements and their error bars (rather than the theory) using the best-fit model with a variant of the procedure described in Sect. 5.3. Since the lensed temperature power spectrum in the best-fit model is very close to that in the fiducial model used to normalise the power spectrum estimates throughout this paper, the power spectrum renormalisation factor $\mathrm{(}\mathrm{1}\mathrm{+}{\mathrm{\Delta }}_{\mathit{L}}^{\mathit{TT}}{\mathrm{)}}^{\mathrm{2}}$ of Eq. (44) differs from unity by less than 0.5%. The predicted ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ in the best-fit model differs from the fiducial model by less than 2.5% for L< 1000. The best-fit model is a good fit to the measurements, with χ² = 10.9 and the corresponding probability to exceed equal to 21%⁵. Significantly, we see that the ΛCDM model, calibrated with the CMB fluctuations imprinted around z = 1100, correctly predicts the evolution of structure and geometry at much lower redshifts. The 68% uncertainty in the ΛCDM prediction of ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ is shown by the dashed lines in the upper-left panel of Fig. 12. We can assess consistency with the direct measurements, properly accounting for this uncertainty, by introducing an additional parameter ${\mathit{A}}_{\mathrm{L}}^{\mathit{\phi \phi }}$ that scales the theory ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ in the lensing likelihood. (Note that we choose not to alter the lensing effect in ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$.) As reported in Planck Collaboration XVI (2014) we find

${\mathit{A}}_{\mathrm{L}}^{\mathit{\phi \phi }}\mathrm{=}\mathrm{0.99}\mathrm{\pm }\mathrm{0.05} \mathrm{\left(}\mathrm{68}\mathrm{\%;}\mathit{Planck}\mathrm{+}\mathrm{lensing}\mathrm{+}\mathrm{WP}\mathrm{+}\mathrm{highL}\mathrm{\right)}\mathit{,}$in excellent agreement with ${\mathit{A}}_{\mathrm{L}}^{\mathit{\phi \phi }}\mathrm{=}\mathrm{1}$. This analysis differs in detail from the simple one-parameter analysis with the fiducial template, reported in Table 1 as A = 0.943 ± 0.04 for the MV reconstruction. The latter does not account for uncertainty in the amplitude or shape of the template between bandpowers due to variations in the six ΛCDM parameters, nor small variations in the lensed temperature power spectrum (which affect the normalisation of the reconstruction power spectrum).

An alternative route to breaking the A_s-τ degeneracy is possible for the first time with Planck. Since ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ is directly proportional to A_s, the lensing power spectrum measurements and the smoothing effect of lensing in ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$ (which at leading order varies as ${\mathit{A}}_{\mathrm{s}}^{\mathrm{2}}{\mathrm{e}}^{\mathrm{-}\mathrm{2}\mathit{\tau }}$) can separately constrain A_s and τwithout large-angle polarization data. The variation of ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ with τ in ΛCDM models constrained only by the Planck temperature power spectrum is illustrated in the upper-right panel of Fig. 12, and suggests that the direct ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ measurements may be able to improve constraints on τ further. This is indeed the case, as shown in Fig. 13 where we compare the posterior distribution of τ for the Planck temperature likelihood alone with that including the lensing likelihood. We find $\begin{array}{ccc}\mathit{\tau }& \mathrm{=}& \mathrm{0.097}\mathrm{\pm }\mathrm{0.038} \mathrm{\left(}\mathrm{68}\mathrm{\%;}\mathit{Planck}\mathrm{\right)}\\ \mathit{\tau }& \mathrm{=}& \mathrm{0.089}\mathrm{\pm }\mathrm{0.032} \mathrm{\left(}\mathrm{68}\mathrm{\%;}\mathit{Planck}\mathrm{+}\mathrm{lensing}\mathrm{\right)}\mathit{.}\end{array}$At 95% confidence, we can place a lower limit on the optical depth of 0.04 (Planck+lensing). This is very close to the optical depth for instantaneous reionization at z = 6, providing further support for reionization being an extended process.

The τ constraints via the lensing route are consistent with, though weaker, than those from WMAP polarization. However, since the latter measurement requires very aggressive cleaning of Galactic emission (see e.g., Fig. 17 of Page et al. 2007), the lensing constraints are an important cross-check.

6.1.2. Effect of the large and small scales on the six-parameter ΛCDM model

Before exploring the further parameters that can be constrained with the lensing likelihood, we test the effect on the ΛCDM model of adding the large-scale (10 ≤ L ≤ 40) and small-scale (400 ≤ L ≤ 2048) lensing data to our likelihood. Adding additional data will produce random shifts in the posterior distributions of parameters, but these should be small here since the multipole range 40 ≤ L ≤ 400 is designed to capture over 90% of the signal-to-noise (on an amplitude measurement). If the additional data are expected to have little statistical power, i.e., the error bars on parameters do not change greatly, but their addition produces large shifts in the posteriors, this would be symptomatic either of internal tensions between the data or an incorrect model.

In Fig. 14, we compare the posterior distributions of the ΛCDM parameters for Planck+WP+highL alone with those after combining with various lensing likelihoods. Adding our fiducial lensing likelihood (second column) reduces the errors on parameters by a small amount and the median values shift by rather less than 1σ for all parameters. The largest gain is for Ω_ch² (and H₀) where the errors improve by 20%. Adding further large- and small-scale data produces no significant reduction in error bars, as expected. Parameter medians also change very little except for Ω_ch², which is dragged low by a further 0.3σ on adding the small-scale lensing information. (The shift in H₀ is due to the anti-correlation between H₀ and Ω_ch² caused by the acoustic-scale degeneracy in the temperature power spectrum; see Planck Collaboration XVI 2014.) These findings are consistent with the power spectrum amplitude measurements discussed in Sect. 6: we can lower the lensing power by reducing the matter density, and this is favoured by the lower amplitudes measured from the small-scale lensing power spectrum.

The tension between the small-scale power and the power over the L = 40–400 range included in our fiducial likelihood, coupled with our lower confidence in the accuracy of the bias removal on small scales, is the reason that we do not include these smaller scales at this stage in the Planck lensing likelihood.

6.1.3. Spatial curvature and dark energy

Inflation models with sufficient number of e-folds of expansion naturally predict that the Universe should be very close to spatially flat. Constraining any departures from flatness is therefore a critical test of inflationary cosmology. However, the primary CMB anisotropies alone suffer from a geometric degeneracy, whereby models with identical primordial power spectra, physical matter densities, and angular-diameter distance to last-scattering have almost identical power spectra (Efstathiou & Bond 1999). The degeneracy is partly broken by lensing (Stompor & Efstathiou 1999), with small additional contributions from the late-ISW effect (on large scales) and by projection effects in curved models (Howlett et al. 2012). In ΛCDM models with curvature, the geometric degeneracy is two-dimensional, involving the curvature and dark energy density, and this limits the precision with which either can be determined from the CMB alone.

With the high-significance detection of lensing by Planck in the temperature power spectrum (Planck Collaboration XVI 2014), and via the lens reconstruction reported here, the geometric degeneracy is partially broken, as shown in Fig. 15. The long tail of closed models with low dark energy density (and expansion rate at low redshift) allowed by the geometric degeneracy have too much lensing power to be consistent with Planck’s measured temperature and lensing power spectra (see also Fig. 12). We find marginalised constraints on the curvature parameter of $\begin{array}{ccc}{\mathrm{\Omega }}_{\mathit{K}}& \mathrm{=}& \mathrm{-}\mathrm{0.04}{\mathrm{2}}_{-0.018}^{\mathrm{+}\mathrm{0.027}} \mathrm{(}\mathrm{68}\mathrm{\%;}\mathit{Planck}\mathrm{+}\mathrm{WP}\mathrm{+}\mathrm{highL}\mathrm{)}\\ {\mathrm{\Omega }}_{\mathit{K}}& \mathrm{=}& \mathrm{-}\mathrm{0.009}{\mathrm{6}}_{-0.0082}^{\mathrm{+}\mathrm{0.010}} \mathrm{(}\mathrm{68}\mathrm{\%;}\mathit{Planck}\mathrm{+}\mathrm{lensing}\mathrm{+}\mathrm{WP}\mathrm{+}\mathrm{highL}\mathrm{)}\mathit{,}\end{array}$so that lensing reconstruction reduces the uncertainty on Ω_K by more than a factor of two over limits driven by the smoothing effect on the acoustic peaks of ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$. This improvement is consistent with the spread in ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ in curved models constrained by the temperature power spectrum, relative to the errors on the reconstruction power spectrum; see Fig. 12. Note that the mean value of Ω_K also moves towards zero with the inclusion of the ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ measurements. Adding the high-L and low-L data to the likelihood brings no more than a percent-level improvement on the constraint. We see that the CMB alone now constrains the geometry to be flat at the percent level. Previous constraints on curvature via CMB lensing have been reported by SPT in combination with the WMAP-7 data: ${\mathrm{\Omega }}_{\mathit{K}}\mathrm{=}\mathrm{-}\mathrm{0.00}{\mathrm{3}}_{-0.018}^{\mathrm{+}\mathrm{0.014}}$ (68%; Story et al. 2013). This constraint is consistent, though almost a factor of two weaker, than that from Planck. Tighter constraints on curvature result from combining the Planck data with other astrophysical data, such as baryon acoustic oscillations, as discussed in Planck Collaboration XVI (2014).

Lensing effects provide evidence for dark energy from the CMB alone, independent of other astrophysical data (Sherwin et al. 2011; van Engelen et al. 2012). In curved ΛCDM models, we find marginalised constraints on Ω_Λ of $\begin{array}{ccc}{\mathrm{\Omega }}_{\mathrm{\Lambda }}& \mathrm{=}& \mathrm{0.5}{\mathrm{7}}_{-0.055}^{\mathrm{+}\mathrm{0.073}} \mathrm{(}\mathrm{68}\mathrm{\%;}\mathit{Planck}\mathrm{+}\mathrm{WP}\mathrm{+}\mathrm{highL}\mathrm{)}\\ {\mathrm{\Omega }}_{\mathrm{\Lambda }}& \mathrm{=}& \mathrm{0.6}{\mathrm{7}}_{-0.023}^{\mathrm{+}\mathrm{0.027}} \mathrm{(}\mathrm{68}\mathrm{\%;}\mathit{Planck}\mathrm{+}\mathrm{lensing}\mathrm{+}\mathrm{WP}\mathrm{+}\mathrm{highL}\mathrm{)}\mathit{.}\end{array}$Again, lensing reconstruction improves the errors by more than a factor of two over those from the temperature power spectrum alone. Note that part of the preference for closed models from the Planck + WP + highL data combination (i.e., without the 4-point function) is related to the tendency of the temperature power spectrum, considered over the full multipole range, to favour models with more lensing power than in the best-fit ΛCDM model. This effect is discussed in some detail in Planck Collaboration XVI (2014).

6.1.4. Neutrino masses

The unique effect in the unlensed temperature power spectrum of massive neutrinos that are still relativistic at recombination is small. With the angular scale of the acoustic peaks fixed from measurements of the temperature power spectrum, neutrino masses increase the expansion rate at z> 1 where dark energy is negligible⁶ and so suppress clustering on scales smaller than the horizon size at the non-relativistic transition (Kaplinghat et al. 2003). This effect reduces ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ for L> 10 (see Fig. 12) and gives less smoothing of the acoustic peaks in ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$. As discussed in Planck Collaboration XVI (2014), the constraint on ∑ m_ν from the Planck temperature power spectrum (and WMAP low-ℓ polarization) is driven by the smoothing effect of lensing: ∑ m_ν< 0.66 eV (95%; Planck+WP+highL). Curiously, this constraint is weakened by additionally including the lensing likelihood to $\begin{array}{c}\sum \end{array}{\mathit{m}}_{\mathit{\nu }}\mathit{<}\mathrm{0.85}\hspace{0.17em}\mathrm{eV}\mathit{,} \mathrm{(}\mathrm{95}\mathrm{\%;}\mathit{Planck}\mathrm{+}\mathrm{lensing}\mathrm{+}\mathrm{WP}\mathrm{+}\mathrm{highL}\mathrm{)}\mathit{,}$(45)reflecting mild tensions between the measured lensing and temperature power spectra, with the former preferring larger neutrino masses than the latter. Possible origins of this tension are explored further in Planck Collaboration XVI (2014) and are thought to involve both the ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ measurements and features in the measured ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$ on large scales (ℓ < 40) and small scales ℓ > 2000 that are not fit well by the ΛCDM+foreground model. In regard to ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$, Fisher estimates show that the bandpowers in the range 130 <L< 309 carry most of the statistical weight in determining the marginal error on ∑ m_ν, and Fig. 12 reveals a preference for high ∑ m_ν from this part of the spectrum. We have checked that removing the first bandpower from the lensing likelihood, which is the least stable to data cuts and the details of foreground cleaning as discussed in Sect. 7, has little impact on our neutrino mass constraints. We also note that a similar trend for lower lensing power than the ΛCDM expectation on intermediate scales is seen in the ACT and SPT measurements (Fig. 11). Adding the high-L information to the likelihood weakens the constraint further, pushing the 95% limit to 1.07 eV. This is consistent with our small-scale measurement having a significantly lower amplitude. At this stage it is unclear what to make of this mild tension between neutrino mass constraints from the 4-point function and those from the 2-point function, and we caution over-interpreting the results. We expect to be able to say more on this issue with the further data that will be made available in future Planck data releases.

6.2. Correlation with the ISW Effect

As CMB photons travel to us from the last-scattering surface, the gravitational potentials they traverse may undergo a non-negligible amount of evolution. This produces a net redshift or blueshift of the photons concerned, as they fall into and then escape from the evolving potentials. The overall result is a contribution to the CMB temperature anisotropy known as the late-time ISW effect, or the Rees-Sciama (R-S) effect depending on whether the evolution of the potentials concerned is in the linear (ISW) or non-linear (R-S) regime of structure formation (Sachs & Wolfe 1967; Rees & Sciama 1968). In the epoch of dark energy domination, which occurs after z ~ 0.5 for the concordance ΛCDM cosmology, large-scale potentials tend to decay over time as space expands, resulting in a net blueshifting of the CMB photons which traverse these potentials.

In the concordance ΛCDM model, there is significant overlap between the large-scale structure which sources the CMB lensing potential φ and the ISW effect (greater than 90% at L< 100), although it should be kept in mind that we cannot observe the ISW component by itself, and so the effective correlation with the total CMB temperature is much smaller, on the order of 20%.

The correlation between the lensing potential and the ISW effect results in a non-zero bispectrum or three-point function for the observed CMB fluctuations. This bispectrum is peaked for “squeezed” configurations, in which one short leg at low-ℓ supported by the ISW contribution is matched to the lensing-induced correlation between two small-scale modes at high-ℓ. Constraints on the amplitude of the lensing-ISW bispectrum using several different estimators are presented in Planck Collaboration XXIV (2014). Here we will present an additional constraint, in which the bispectrum measurement is recast as an estimate for the amplitude of the cross-spectrum ${\mathit{C}}_{\mathit{L}}^{\mathit{T\phi }}$, using the filtering and frequency map combinations of our baseline lensing reconstruction. Our measurements are in good agreement with those made in Planck Collaboration XXIV (2014); a detailed comparison of several lensing-ISW bispectrum estimators, including the one used here, is presented in Planck Collaboration XIX (2014).

Following Lewis et al. (2011), we begin with an estimator for the cross-spectrum of the lensing potential and the ISW effect as $\begin{array}{ccc}\mathit{Ĉ}\begin{array}{c}\mathit{T\phi }\\ \mathit{L}\end{array}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{sky}}^{-1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{\mathit{M}}\mathit{T̂}\mathit{LM}\mathit{\phi ̂}\begin{array}{c}\mathrm{\ast }\\ \mathit{LM}\end{array}\mathit{,}& & \end{array}$(46)where $\mathit{T̂}{}_{\mathit{\ell m}}\mathrm{=}{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\mathit{T̅}{}_{\mathit{\ell m}}$ is the Wiener-filtered temperature map and is given in Eq. (13). In Fig. 16 we plot the measured cross-spectra for our individual frequency reconstructions at 100, 143, and 217 GHz as well as the MV reconstruction. We also plot the mean and scatter expected in the fiducial ΛCDM model.

To compare quantitatively the overall level of the measured ${\mathit{C}}_{\mathit{L}}^{\mathit{T\phi }}$ correlation to the value in ΛCDM, we estimate an overall amplitude for the cross-spectrum as $\mathit{Â}\mathit{T\phi }\mathrm{=}{\mathrm{𝒩}}^{\mathit{T\phi }}\sum _{\mathit{L}\mathrm{=}{\mathit{L}}_{\min }}^{{\mathit{L}}_{\max }}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{C}}_{\mathit{L}}^{\mathit{T\phi ,}\mathrm{fid}\mathit{.}}\mathit{Ĉ}\begin{array}{c}\mathit{T\phi }\\ \mathit{L}\end{array}\mathit{/}\mathrm{\left(}{\mathit{C}}_{\mathit{L}}^{\mathit{TT}}{\mathit{N}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{\right)}\mathit{.}$(47)The overall normalization is determined from Monte-Carlo simulations. For our processing of the data, we find that it is well approximated (at the 5% level) by the analytical approximation ${\mathrm{𝒩}}^{\mathit{T\phi }}\mathrm{\approx }{\left[\sum _{\mathit{L}\mathrm{=}{\mathit{L}}_{\min }}^{{\mathit{L}}_{\max }}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathrm{(}{{\mathit{C}}_{\mathit{L}}^{\mathit{T\phi ,}\mathrm{fid}\mathit{.}}}^{\mathrm{)}}}^{\mathrm{2}}\mathit{/}\mathrm{\left(}{\mathit{C}}_{\mathit{L}}^{\mathit{TT}}{\mathit{N}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{\right)}\right]}^{-1}\mathit{.}$(48)The estimator above is equivalent to the KSW and skew-C_ℓ estimators of Komatsu et al. (2005); Munshi et al. (2011b) for the lensing-ISW bispectrum that are used in Planck Collaboration XXIV (2014; up to implementation details such as filtering). The mean-field subtraction performed when computing can be identified with the linear term of Creminelli et al. (2006), which is necessary to minimize the estimator variance. The contribution to the total S/N of this estimator as a function of the short leg L is plotted in Fig. 2, where it can be seen that the constraining power for the fiducial correlation is almost entirely at L< 100.

In Table 2 we present measured values for the amplitude of the lensing-ISW bispectrum using Eq. (46). The uncertainties on Â^Tφ are determined by Monte-Carlo. We use the multipole range 10 <L< 100, given some of the potential systematic issues with L< 10 identified in Sect. 7.4, although as can be seen from Fig. 16, the inclusion of lower multipoles does not significantly affect our results. Note that for the ISW-lensing measurements, inaccuracies in the mean-field subtraction do not bias the estimator although they may degrade the statistical errors on large scales. The differences between the different amplitude fits are well within the expected scatter, as we show in Table 3.

As a point of interest, we have also split our amplitude constraint into the contribution from even and odd multipoles. There are well known odd/even-multipole power asymmetries in the temperature anisotropies on large angular scales, the study of which is somewhat limited by the small number of available modes (Land & Magueijo 2005; Kim & Naselsky 2010; Gruppuso et al. 2011; Bennett et al. 2011). The lensing potential gives a potentially new window on these power asymmetries, as a third somewhat independent measurement of power on large angular scales. As we can see in Table 2, there is a large difference between the odd/even-L contributions in our lensing-ISW bispectrum estimate, related to the odd/even temperature asymmetry.

In Table 3 we present differences between the amplitude measurements of Table 2, as well as the expected scatter, accounting for correlations between the estimates due to common CMB and (in the case of MV vs. 143 or 217 GHz) noise fluctuations. We see that our estimates are all very consistent. The difference in measured odd/even multipole power, while striking, is not statistically significant at greater than 2σ.

To conclude, we see a correlation between our lensing potential measurement and the large-scale temperature anisotropies that is consistent with the level expected due to the present-day dark energy domination in our fiducial ΛCDM cosmology. Our amplitude measurements are in agreement with those presented in Planck Collaboration XXIV (2014); these independent measurements are compared in detail in Planck Collaboration XIX (2014).

6.3. Correlation with galaxy catalogues

The CMB lensing potential is an integrated measure of all matter in the Universe back to the last-scattering surface, with a broad kernel which peaks at z ~ 2, but has significant contributions from both lower and higher redshifts. There are therefore significant correlations between the CMB lensing potential and other tracers of large-scale structure. Such correlations have already been observed between the lensing potential reconstructed from WMAP data and the NVSS quasar catalogue (Smith et al. 2007; Hirata et al. 2008, on the order of 3σ in significance), as well as between the lensing potential from SPT with optical and infrared catalogues from BCS, WISE, and Spitzer/IRAC (Bleem et al. 2012, between 4 and 5σ), and the ACT lensing potential with quasars from SDSS (Sherwin et al. 2012, at 3.8σ). Here we will show several representative examples of correlations with the Planck lensing potential, between 7σ and 20σ in significance. Our goal is not to perform an in depth study, but rather to demonstrate the power of our public lensing map. In addition to the correlations presented here, there is a powerful correlation with, the cosmic infrared background (CIB). Correlation of the Planck lensing potential with the CIB fluctuations as probed by the highest frequency Planck channels is observed at greater than 40σ in significance, and has been subjected to a more detailed analysis and modelling, which is presented in an accompanying paper, Planck Collaboration XVIII (2014).

To predict expected levels of correlation for a given galaxy catalogue, we use the Limber approximation (Limber 1954) with a simple linear bias model, in which the cross-spectrum between two mass tracers is given by ${\mathit{C}}_{\mathit{L}}^{\mathit{g\phi }}\mathrm{=}\mathrm{\int }\mathrm{d}\mathit{\chi }{\mathit{K}}_{\mathit{L}}^{\mathit{g}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}{\mathit{K}}_{\mathit{L}}^{\mathit{\phi }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}\mathit{P}\mathrm{(}\mathit{k}\mathrm{=}\mathit{L}\mathit{/}\mathit{\chi ,\chi }\mathrm{)}\mathit{,}$(49)where χ is conformal distance and P(k,χ) is the 3D matter power spectrum for wavenumber k at conformal lookback time χ. K^φ and K^g are kernels associated with lensing and with the galaxy catalogue of interest respectively, and are given by $\begin{array}{ccc}{\mathit{K}}_{\mathit{L}}^{\mathit{\phi }}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}& \mathrm{=}& \mathrm{-}\frac{\mathrm{3}{\mathrm{\Omega }}_{\mathrm{m}}{\mathit{H}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{L}}^{\mathrm{2}}}\frac{\mathit{\chi }}{\mathit{a}}\left(\frac{{\mathit{\chi }}_{\mathrm{\ast }}\mathrm{-}\mathit{\chi }}{{\mathit{\chi }}_{\mathrm{\ast }}\mathit{\chi }}\right)\mathit{,}\\ {\mathit{K}}_{\mathit{L}}^{\mathit{g}}\mathrm{\left(}\mathit{\chi }\mathrm{\right)}& \mathrm{=}& \frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{z}}\frac{\mathrm{d}\mathit{z}}{\mathrm{d}\mathit{\chi }}\frac{\mathit{b}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathit{\chi }}\mathrm{·}\end{array}$(50)Here b(z) is a redshift-dependent linear bias parameter, and dN/dz describes the redshift distribution of the galaxy population. As with the lensing-ISW bispectrum, we can construct a simple pseudo-C_L estimator for the cross correlation using $\mathit{Ĉ}\begin{array}{c}\mathit{g\phi }\\ \mathit{L}\end{array}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{sky}\mathit{,}\mathrm{j}}^{-1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{\mathit{M}}{\mathit{g}}_{\mathit{LM}}\mathit{\phi ̂}\begin{array}{c}\mathrm{\ast }\\ \mathit{LM}\end{array}\mathit{,}$(51)where f_sky,j is the sky fraction common to both the lens reconstruction and the catalogue, and g_LM is the harmonic transform of the galaxy fractional overdensity. Denoting the positions of the objects p in the catalogue as , the transform is given by ${\mathit{g}}_{\mathit{LM}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{N}}^{\mathit{gg}}}\sum _{\mathit{p}}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{n̂}\mathit{p}\mathrm{\right)}\mathit{,}$(52)where N^gg is the surface density of objects in gal/steradian. With a fiducial model ${\mathit{C}}_{\mathit{L}}^{\mathit{g\phi ,}\mathrm{fid}\mathit{.}}$ for the cross-spectrum obtained using Eq. (49), the minimum-variance estimator for its overall amplitude is $\mathit{Â}\mathit{g\phi }\mathrm{=}{\mathrm{𝒩}}^{\mathit{g\phi }}\sum _{\mathit{L}}\frac{\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{C}}_{\mathit{L}}^{\mathit{g\phi ,}\mathrm{fid}\mathit{.}}\mathit{Ĉ}\begin{array}{c}\mathit{g\phi }\\ \mathit{L}\end{array}}{\mathrm{(}{\mathit{C}}_{\mathit{L}}^{\mathit{gg}}\mathrm{+}{\mathit{N}}_{\mathit{gg}}^{}\mathrm{)}\mathrm{(}{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{+}{\mathit{N}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{)}}\mathrm{\equiv }\sum _{\mathit{L}}{\mathit{A}}_{\mathit{L}}^{\mathit{g\phi }}\mathit{Ĉ}\begin{array}{c}\mathit{g\phi }\\ \mathit{L}\end{array}\mathit{,}$(53)where ${\mathit{C}}_{\mathit{L}}^{\mathit{gg}}$ is the signal power spectrum of the catalogue, which can be estimated using Eq. (49)with K^g for both weight functions. Here we have defined the spectrum ${\mathit{A}}_{\mathit{L}}^{\mathit{g\phi }}$ as a scaling at each multipole L of the cross-correlation power spectrum. The normalization is given by ${\mathrm{𝒩}}^{\mathit{g\phi }}\mathrm{=}{\left[\sum _{\mathit{L}}\frac{\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathrm{(}{{\mathit{C}}_{\mathit{L}}^{\mathit{g\phi ,}\mathrm{fid}\mathit{.}}}^{\mathrm{)}}}^{\mathrm{2}}}{\mathrm{(}{\mathit{C}}_{\mathit{L}}^{\mathit{gg}}\mathrm{+}{\mathit{N}}_{\mathit{gg}}^{}\mathrm{)}\mathrm{(}{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{+}{\mathit{N}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{)}}\right]}^{-1}\mathit{.}$(54)In Fig. 17 we plot the contributions to Â^gφ as a function of L for several surveys that have significant correlations with the Planck MV lensing potential: the NVSS quasar catalogue, the MaxBCG cluster catalogue, an SDSS LRG catalogue, and an infrared catalogue from the WISE satellite. The error bars for each correlation are measured from the scatter of simulated lens reconstructions correlated with each catalogue map, and are in generally good agreement (at the 20% level) with analytical expectations. These catalogues are discussed in more detail below.

• 1.

  NVSS Quasars: the NRAO VLA Sky Survey (NVSS; Condonet al. 1998) is a catalogue ofapproximately two million sources north ofδ = −40deg which is 50% complete at 2.5 mJy. Most of the bright sources are AGN-powered radio galaxies and quasars. We process this catalogue following Smith et al. (2007), pixelizing the catalogue at HEALPix N_side = 256 and projecting out the azimuthally-symmetric modes of the galaxy distribution in ecliptic coordinates to avoid systematic striping effects in the NVSS dataset. We model the expected cross-correlation for this catalogue using a constant b(z) = 1.7 and a redshift distribution centered at z₀ = 1.1 given by $\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{z}}\mathrm{\propto }\{\begin{array}{c}\\ \exp \left(\mathrm{-}\frac{\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}\mathrm{(}\mathrm{0.8}{\mathrm{)}}^{\mathrm{2}}}\right)& \mathrm{(}\mathit{z}\mathrm{\le }{\mathit{z}}_{\mathrm{0}}\mathrm{)}\\ \exp \left(\mathrm{-}\frac{\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}\mathrm{(}\mathrm{0.3}{\mathrm{)}}^{\mathrm{2}}}\right)& \mathrm{(}\mathit{z}\mathrm{\ge }{\mathit{z}}_{\mathrm{0}}\mathrm{)}\mathit{.}\end{array}$(55)For this model, in the correlation with the MV lens reconstruction we measure an amplitude of $\mathit{Â}\begin{array}{c}\mathit{g\phi }\\ \mathrm{NVSS}\end{array}\mathrm{=}\mathrm{1.03}\mathrm{\pm }\mathrm{0.05}$.

• 2.

  SDSS LRGs: We use the LRG catalogue of Ross et al. (2011); Ho et al. (2012) based on Sloan Digital Sky Survey Data Release 8 (SDSS DR8), which covers 25% of the sky. After cutting to select all sources with photometric redshift 0.4 ≤ z ≤ 0.8, and p_gal> 0.2, we are left with approximately 1.4 × 10⁶ objects with a mean redshift of z = 0.55 and a scatter of ± 0.07. Apart from the cut above, we do not perform any additional weighting on p_gal. We model this catalogue using dN/dz taken from the histogram of photometric redshifts, and take b(z) = 2. We measure $\mathit{Â}\begin{array}{c}\mathit{g\phi }\\ \mathrm{LRGs}\end{array}\mathrm{=}\mathrm{0.96}\mathrm{\pm }\mathrm{0.10}$, very consistent with expectation.

• 3.

  MaxBCG Clusters: The MaxBCG cluster catalogue (Koester et al. 2007) is a collection of 13,823 clusters over approximately 20% of the sky selected from the SDSS photometric data, covering a redshift range 0.1 ≤ z ≤ 0.3. It is believed to be 90% pure and more than 85% complete for clusters with M ≥ 1 × 10¹⁴M_⊙. To simplify the sky coverage, we have discarded the three southern SDSS stripes in the catalogue, which reduces the overall sky coverage to approximately 17%. There are accurate photometric redshifts (Δ_z ~ 0.01) for all objects in the catalogue, and so we can construct dN/dz directly from the histogram of the redshift distribution. Although these clusters are at very low redshift compared to the typical structures which source the CMB lensing potential, they are strong tracers of dark matter, with an effective bias parameter of b(z) = 3 (Hütsi 2010). We obtain a similar average bias parameter ⟨ b(M,z) ⟩ for the MaxBCG clusters if we combine the mass-richness relation of Bauer et al. (2012) and the halo bias prescription of Tinker et al. (2010). Here we measure a correlation with the Planck lensing potential of $\mathit{Â}\begin{array}{c}\mathit{g\phi }\\ \mathrm{MaxBCG}\end{array}\mathrm{=}\mathrm{1.54}\mathrm{\pm }\mathrm{0.21}$. This is significantly larger than expected given the simple model above, although as can be seen in Fig. 17 the shape of the correlation is in reasonable agreement.

• 4.

  WISE Catalogue: The Wide Field Survey Infrared Explorer (WISE) satellite (Wright et al. 2010) has mapped the full sky in four frequency bands W1 − W4 at 3.4, 4.6, 12, and 22 μm respectively. We start from the full mission catalogue, which contains over 560 × 10⁶ objects. To obtain a catalogue with roughly uniform sensitivity over the full sky and to eliminate stellar contamination we follow Kovacs et al. (2013), selecting all sources with W1 magnitudes less than 15.2 at galactic latitudes greater than 10°, and require W1 − W2 > 0.2 and W2 − W3 > 2.9. We cut all sources that are flagged as potentially spurious, or for which more than 30% of frames observing a given source are marked as contaminated by moon-glow. These cuts leave us with 2 308 751 sources. We additionally remove all sources that lie outside an f_sky = 0.6 Galaxy mask. Following Goto et al. (2012), we model the redshift distribution of these sources using $\frac{\mathrm{d}\mathit{N}}{\mathrm{d}\mathit{z}}\mathrm{\propto }\mathit{z}\exp \left(\mathrm{-}\frac{\mathrm{(}\mathit{z}\mathrm{-}{\mathit{z}}_{\mathrm{0}}{\mathrm{)}}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_{\mathit{z}}^{\mathrm{2}}}\right)\mathit{,}$(56)with z₀ = 0.1, σ_z = 0.1, and where the proportionality constant is chosen such that ^∫dzdN/ dz = 1. The bias is taken to be b = 1. For this model we measure an amplitude of $\mathit{Â}\begin{array}{c}\mathit{g\phi }\\ \mathrm{WISE}\end{array}\mathrm{=}\mathrm{0.97}\mathrm{\pm }\mathrm{0.13}$. Note that the cuts we have used above are fairly conservative – in particular, relaxing the magnitude cut above we can obtain much larger cross-correlation signals (on the order of 40σ in significance).

It can be seen from Fig. 17 that we have significant correlations with the objects in all four of these catalogues, with amplitudes and shapes roughly at the expected level. The correlations here are only a demonstrative sample of what is possible with the Planck lensing potential map. Note that although the lensing potential has a peak sensitivity to high-redshift tracers (such as NVSS), significant correlations may also be observed with relatively low redshift catalogues. From the perspective of cross-correlation, the Planck lensing potential represents a mass survey with a well understood redshift kernel over the full-sky, and we anticipate that correlations with it will prove useful as a direct probe of the correlation between luminous and dark matter.

7. Consistency and systematic tests

In this section we perform a number of consistency and systematic tests. Most of our tests will focus on the robustness of the MV reconstruction that forms the basis for our fiducial lensing likelihood. In Fig. 18 we summarize several of our main consistency checks: foreground tests obtained by comparing individual frequencies, adjusting the Galactic mask level, or comparing to results obtained analyzing component-separated maps (discussed in Sect. 7.1), tests of our correction at the power spectrum level for contamination by the shot-noise from unresolved point-source (in Sect. 7.2), tests for possible noise bias obtained by taking appropriate correlations of maps with independent noise realizations (discussed further in Sect. 7.3), and tests using alternative estimators that are bias hardened against various systematics (Sect. 7.4). In all cases, we find that our MV power spectrum is generally consistent with these alternative reconstructions.

Our baseline reconstruction method uses the methodology outlined in Sect. 2, as well as a specific choice of filtering to produce the that are the inputs to our lensing estimators. We test our robustness to these choices in Sect. 7.5, as well as the implementation of our main reconstruction pipeline, by comparing to three alternative and largely independent pipelines.

In addition to these targeted tests, we have also performed several generic null tests for unexpected signals in our data. These are deferred to Appendix A, where we study survey-to-survey discrepancies, lensing curl modes, and the map-level statistics of our lens reconstruction.

7.1. Foreground contamination

In this section, we test the consistency of our results with respect to the treatment of Galactic foreground contamination. We use several different tests to gain perspective on possible issues:

• 1.

  We compare individual frequency reconstructions. The agreement of the lensing power spectra estimated at 100, 143, and 217 GHz is apparent within the error bars plotted in Fig. 10. However, because a large fraction of the uncertainty in these estimates comes from CMB fluctuations themselves they are significantly correlated. In Fig. 18 we plot the difference between the 143 and 217 GHz reconstructions (the scatter for 100 GHz is too large to provide a useful test of the MV reconstruction so we have not plotted it, although it is also consistent). Both are in reasonable agreement, although we can see that the somewhat large amplitude for our first likelihood bin, just above L = 40 is being driven by 217 GHz.

• 2.

  We perform lens reconstruction using both more aggressive and more conservative Galaxy masks than our fiducial f_sky = 0.7 analysis mask, constructed following the description in Sec. 3. The results of this test are also shown in Fig. 18, for masks with f_sky = 0.6 and f_sky = 0.8. We can see a highly significant discrepancy at very low-L between our fiducial f_sky = 0.7 mask and the more aggressive f_sky = 0.8 mask, however other bins (in particular in the 40 <L< 400 range) are untouched. Our fiducial results are completely consistent with the more conservative f_sky = 0.6 results.

• 3.

  We perform lens reconstruction on the component-separated CMB maps from Planck Collaboration XII (2014). These are the product of four different foreground-removal algorithms: Ruler, Nilc, Sevem, and Smica. Each of these methods combines the full set of nine Planck frequency maps from 30 → 857 GHz to obtain a best-estimate CMB map. This makes a detailed characterization of these maps more difficult. However analysis of these maps does provide a valuable consistency test. In Fig. 18 we compare the reconstructions obtained on these component-separated maps to our MV reconstruction, analyzed using the same Galaxy mask. The foreground-cleaned maps are generally in good agreement with our fiducial reconstruction, although the first bin above L = 40 is consistently about one standard-deviation lower than our MV result, which as we have already pointed out is driven high by 217 GHz. The confidence masks produced in Planck Collaboration XII (2014) cover significantly more sky than the f_sky = 0.7 mask used in our baseline analysis.

The remainder of this subsection describes our analysis of the component-separated maps in more detail. For the power spectrum tests of Fig. 18, we analyze these maps following the same methodology as for our fiducial results. As a mask, we use the union of our fiducial f_sky = 0.7 Galaxy mask, point source masks at 100, 143, and 217 GHz, as well as the “method” masks provided for each component separation algorithm. We characterize the bias terms of Eq. (15)as well as the statistical uncertainties of the lens reconstruction on each map using the same lensed CMB and FFP6 noise realizations that are discussed in Sect. 4. These have been run through each of the component-separation algorithms above, resulting in a set of simulated foreground-cleaned maps (note that no foreground signal was included in the simulations which were run through the separation algorithms, and so these simulations only reproduce the mixing employed by each algorithm of the individual frequency maps). To match the power spectrum of these simulations to the power spectrum of the data maps, we find it is necessary to add extragalactic foreground power following the model in Sect. 4, with A_cib = 18 μK² and A_src = 28 μK². The resulting simulations have a power spectrum that agrees with that of the CMB map estimate based on the data to better than 2% at l< 2048. This could be improved slightly by tailoring a specific correction for each map. We also add homogeneous pixel noise with a level of 12 μK arcmin. If we neglected this power, the agreement would be only at the 8% level, primarily due to the noise term (the A_cib and A_src contributions are each at the level of 1−2%). Due to the procedure that we use to subtract the disconnected noise bias (Eq. (17)) from our lensing power spectrum estimates, the inclusion of these components does not significantly affect our results, but comparison with the values used for our single-frequency simulations in Sect. 4 are a useful indicator of the extent to which the component-separation algorithms are able to remove extragalactic foreground power in the high-ℓ regime.

As already discussed, our results on the component-separated CMB maps are presented in Fig. 18. Because the CMB and FFP6 noise components of the foreground-cleaned map simulations are the same as those used to characterize our fiducial lens reconstruction, we can measure the expected scatter between the component-separated maps and our fiducial reconstruction. This scatter will be slightly overestimated because we have not attempted to model coherently the contribution to the reconstruction noise from residual diffuse extragalactic foreground power. For the eight bins in 40 ≤ L ≤ 400 on which our fiducial likelihood is based, we measure a χ² for the difference between our fiducial reconstruction and the corresponding foreground-cleaned reconstruction of χ² = (3.14,4.3,2.5,14.7) for Nilc, Smica, Sevem, and Ruler respectively. These χ² values associated have probability-to-exceed (PTE) values of (79%,64%,86%,2%). At the level that we are able to test, the Nilc, Smica, and Sevem foreground-cleaned maps give results that are quantitatively consistent with our fiducial reconstruction. There is more scatter between our fiducial reconstruction and the Ruler map than expected from simulations, as evidenced by a very high χ² for the difference, however as can be seen in Fig. 18, there are not any clear systematic differences. Indeed, the discrepancy for the bins plotted in Fig. 18 (which differ somewhat from the linear bins used in our likelihood) is much less significant than for the bins of our fiducial likelihood.

When using the component-separated maps above, we have used the same f_sky = 0.7 Galactic mask as for our MV result, although the confidence regions associated with each foreground-cleaned map allow more sky, ranging up to f_sky = 0.94 for the Nilc method. We have used the metis pipeline (described later in Sect. 7.5) to test whether this improved sky coverage could benefit our lens reconstruction. The same method has been used in Planck Collaboration XII (2014) to evaluate possible biases to lens reconstruction induced by these methods using the FFP6 simulated CMB realization, described in Planck Collaboration I (2014), indicating that the different component-separation algorithms do not alter significantly the lensing signal (at the level which can be tested on a single simulation). Analyzing the Nilc map, which has the largest confidence region, we find that we can increase the usable sky area up to f_sky = 0.87 without encountering significant Galactic contamination. In Fig. 19 we show the striking improvement in sky coverage on the Nilc map. Smica and Sevem are very similar; we have not considered Ruler because of its larger noise level.

Power spectrum estimates at this mask level show consistency with the MV reconstruction within two standard deviations of the measurement uncertainty. The increased sky coverage does not bring significant improvements in the error-bars of the power spectrum, however. Using Eq. (20) as an estimate of the power spectrum variance, the larger sky coverage yields only a 3.5% improvement at L< 40 over the MV result, decreasing down to 0 at L = 400. This could be due to the different weighting used in the component separation compared to the one of the MV map, which results in slightly noisier maps for our purpose. While the component-separated maps allow for a reduced mask, they lead to a marginal improvement of the power spectrum uncertainties. Nevertheless, their agreement with the MV result is reassuring.

7.2. Point source correction

As can be seen in Table 1, the unresolved point-source shot-noise correction in any individual band for our MV likelihood is on the order of a few percent, reaching up to 6% for the highest multipole bands. Averaged over the 40 ≤ L ≤ 400 band, the shot noise correction amounts to a 2% shift in the amplitude of $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$, which is small but non-negligible compared to our statistical uncertainty of 4%. Physically, the amplitudes of our source corrections are reasonable; at 143 GHz we measure $\mathit{Ŝ}\begin{array}{c}\mathrm{4}\\ \mathrm{143}\end{array}\mathrm{=}\mathrm{(}\mathrm{1.3}\mathrm{\pm }\mathrm{0.6}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-12}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{4}}$. From the radio point-source model of De Zotti et al. (2010), this corresponds to an effective flux cut of approximately 150 mJy at this frequency, roughly comparable to that expected for the S/N> 5 cut we make when masking sources in our fiducial analysis (Planck Collaboration XXVIII 2014). The shot noise measured at 217 GHz is lower, as expected given the smaller contribution from radio sources, with $\mathit{Ŝ}\begin{array}{c}\mathrm{4}\\ \mathrm{217}\end{array}\mathrm{=}\mathrm{(}\mathrm{0.4}\mathrm{\pm }\mathrm{0.4}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-12}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{4}}$. The shot-noise level that we measure for the MV reconstruction lies between these two, with $\mathit{Ŝ}\begin{array}{c}\mathrm{4}\\ \mathrm{MV}\end{array}\mathrm{=}\mathrm{(}\mathrm{0.6}\mathrm{\pm }\mathrm{0.3}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-12}\hspace{0.17em}\mathit{\mu }{\mathrm{K}}^{\mathrm{4}}$.

As a consistency test of our S/N> 5 source masks, we construct two new sets of source masks that only exclude S/N> 10 sources (to check the stability of our reconstruction under significantly higher source contamination), or that exclude all detected sources with S/N> 4.2 (to test the conservatism of our fiducial flux cut). The results of this test are shown in the masks panel of Fig. 18. Although not shown explicitly in the figure, the source correction term with the S/N> 10 mask is considerably larger than our fiducial correction. For this mask we measure $\mathit{Ŝ}\begin{array}{c}\mathrm{4}\\ \mathrm{MV}\end{array}\mathrm{=}\mathrm{(}\mathrm{2.8}\mathrm{\pm }\mathrm{0.3}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-12}\hspace{0.17em}\mathit{\mu }$K⁴, now a 9σ detection of shot noise in the trispectrum. The importance of the correction grows correspondingly, becoming a 9% overall correction to the ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ amplitude in 40 ≤ L ≤ 400. With the source correction, we measure $\mathit{Â}\begin{array}{c}\mathrm{MV}\\ \mathrm{40}\mathrm{\to }\mathrm{400}\end{array}\mathrm{=}\mathrm{0.95}\mathrm{\pm }\mathrm{0.04}$, consistent with the $\mathit{Â}\begin{array}{c}\mathrm{MV}\\ \mathrm{40}\mathrm{\to }\mathrm{400}\end{array}\mathrm{=}\mathrm{0.94}\mathrm{\pm }\mathrm{0.04}$ measurement for our fiducial source mask. Without the source correction we would have measured $\mathit{Â}\begin{array}{c}\mathrm{MV}\\ \mathrm{40}\mathrm{\to }\mathrm{400}\end{array}\mathrm{=}\mathrm{1.04}\mathrm{\pm }\mathrm{0.04}$, significantly discrepant. We can see that even for this extremely permissive source mask, our correction for the point-source shot-noise trispectrum does a good job of rendering $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$ insensitive to source contamination. Turning to the S/N> 4.2 mask, we measure $\mathit{Ŝ}\begin{array}{c}\mathrm{4}\\ \mathrm{MV}\mathit{,}\mathrm{SNR}\mathrm{4}\end{array}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{0.1}\mathrm{\pm }\mathrm{0.3}\mathrm{)}\mathrm{\times }{\mathrm{10}}^{-12}\hspace{0.17em}\mathit{\mu }$K⁴, reduced from a 2σ measurement of shot noise to one very consistent with zero. For this reconstruction we measure Â_{40 → 400} = 0.93 ± 0.04, in good agreement with the MV result, although there does appear to be an overall trend for lower power above L = 400 in the S/N> 4.2 mask. At a detailed level, both of these results do fail χ² tests on the difference with the MV reconstruction. As the change in sky fraction between these three reconstructions is small (the difference between the fiducial and S/N> 4.2 mask is 1%, and even smaller for the S/N> 10 mask), even the small differences quoted above are significantly larger than they should be if the source correction were perfect. This can be seen in Fig. 18, where the error bars on the difference between each of these reconstructions and the MV result is too small to be visible for the plotted points. For the 40 ≤ L ≤ 400 range on which our likelihood is based, we believe that the discrepancy between the power spectra estimated with these source masks is acceptable. Particularly given the results of the S/N> 10 cut, the source correction that we perform appears to behave reasonably, and reduces the sensitivity of our estimates to the shot-noise trispectrum of unmasked point sources. In Sect. 7.4, we will also use a modified lensing estimator that has zero response to the shot-noise trispectrum and therefore does not require a source correction. As with the tests above, we will find good overall agreement with the standard estimator used for the MV reconstruction.

7.3. Instrumental noise bias

Errors in the instrumental noise model are a source of concern for our lensing power spectrum estimates, in which all four legs of the trispectrum estimate come from a single map. In this situation there are three Wick contractions of the noise which contribute to the power spectrum estimate in two distinct ways, depending on whether they couple between the two quadratic estimators in the power spectrum estimate or not:

• 1.

  The Wick contraction that couples noise contributions inside the same estimator leads to a mean-field term as in Eq. (9). This term is sourced by the statistical anisotropy of the noise (due to the uneven hit distribution across the sky, or the correlations along the scan direction induced by deconvolution of the bolometer time constant). The shape of this mean-field is dictated by the Planck scan strategy, which is coherent over large scales, and so this mean-field is primarily at low-L, below the L_min = 40 cutoff of our likelihood. It can be avoided by forming a quadratic estimator in which the two input maps have different noise realizations. We perform a null test against this contribution in Fig. 18, where the (143 × 217) × (143 × 217) panel shows a comparison of our MV power spectrum estimate to that obtained using a quadratic estimator in which one leg comes from 143 GHz and one from 217 GHz. We see generally good agreement with the MV result. For the eight bins in 40 ≤ L ≤ 400 of our fiducial likelihood, we find a χ² for the difference (accounting for correlations between the two estimates) of 4.2, which has a PTE of 83%. The noise mean-field contribution can also be reduced by using an appropriate bias-hardened estimator. We will perform additional tests with such estimators in Sect. 7.4.

• 2.

  The two Wick contractions that couple noise terms across-estimators contribute to the disconnected bias ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{0}}$. This term describes the noise contribution to the reconstruction variance. Our use of the data covariance matrix in Eq. (17)makes our estimate of the disconnected noise bias insensitive to small errors in the noise model. The noise contribution can also be avoided by taking the cross-spectrum of two quadratic estimators with independent noise. We show the result of such a correlation compared to our MV reconstruction in Fig. 18, where the (143 × 143) × (217 × 217) panel shows the cross-spectrum between the individual 143 and 217 GHz reconstructions. Again, we see generally good agreement with the MV result. For the eight bins in 40 ≤ L ≤ 400 of our fiducial likelihood, we find a χ² for the difference (accounting for correlations between the two estimates) of 6.5, which has a PTE of 59%.

Based on these tests, we believe that noise bias cannot be a significant contaminant for our power spectrum estimates.

7.4. Bias-hardened estimators

Most of the results in this paper use the standard weight function ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}$ to form quadratic estimates for the lensing potential. This weight function is a matched filter for lensing, and results in estimators with minimal variance. As we discuss in Appendix C, however, the standard estimator has large mean-field contributions at low-L, primary due to masking and noise inhomogeneity. As pointed out in Namikawa et al. (2013), these mean-fields can be mitigated with the use of “bias-hardened” lensing estimators, which use weight functions ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{x}}$ specially constructed to project out these effects. Following the notation of Appendix C, we consider an uncertain bias field z_LM with associated weight function ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{z}}$ that gives a contribution to the off-diagonal elements of the CMB covariance matrix (and therefore the estimator mean-field) of $\mathrm{\Delta }\mathrm{⟨}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}\mathrm{⟩}\mathrm{=}\sum _{\mathit{LM}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& \mathit{L}\\ {\mathit{m}}_{\mathrm{1}}& {\mathit{m}}_{\mathrm{2}}& \mathit{M}\end{array}\right){\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{z}}{\mathit{z}}_{\mathit{LM}}\mathit{.}$(57)We can then form an estimator which is insensitive to z_LM simply as ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathrm{(}\mathit{\phi }\mathrm{-}\mathit{z}\mathrm{)}}\mathrm{=}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}\mathrm{-}\frac{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi z}}}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{zz}}}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{z}}\mathit{.}$(58)The quadratic estimator (φ − z) that uses the weight function above has a significantly reduced mean-field contribution from z_LM. This bias-hardening procedure may be repeated iteratively to produce weight functions that are insensitive to several sources of mean-field simultaneously. In Fig. 20 we plot a comparison at low-L between estimates with the standard estimator and estimators that are bias hardened against noise and mask mean-fields. We see generally good agreement between these sets of estimates for L ≥ 10. Below L = 10, the pseudo-spectrum of the difference between the standard and bias-hardened estimators is unexpectedly large. For this reason, we band-limit our lensing potential map to L ≥ 10, and do not consider L< 10 for any of the quantitative results in this paper.

In addition to low-L consistency tests, we may also use these bias-hardened estimators to estimate ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$. The resulting power spectra are plotted in Fig. 18. Taking the difference with the eight bins in 40 ≤ L ≤ 400 of our fiducial MV reconstruction, with scatter estimated from simulations, we find χ² values of 6.6 and 4.3 for the mask/noise-hardened estimators respectively. These have corresponding PTE values of 58% and 83%. The error bars on ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ obtained with the bias hardened estimators used here are generally between 10 and 20% larger than the error bars of the standard estimator. We also construct an estimator that is bias hardened against the point source weight function ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{s}}$ of Eq. (28). This estimator has the distinction of having zero response to point-source shot-noise, as well as any correlation between the point-source shot noise power and the lensing potential. Again, we find consistent results with the standard estimator, with χ² = 5 and a PTE of 76%.

7.5. Alternative methods

All of the primary results in this paper use a lens reconstruction pipeline based on the methodology outlined in Sect. 2. As a robustness test, both of this methodology and of its implementation, we have implemented three independent reconstruction pipelines. These independent pipelines make significantly different choices than our baseline approach, primarily in the choice of the inverse-variance filter function and the calculation of the correction terms in Eq. (15). Our main result of this section is Fig. 21, where we compare the lens reconstruction power spectra for these alternative pipelines to our baseline results at 143 and 217 GHz, as well as the MV reconstruction. The agreement is excellent. The only aspect of implementation common to these four results (apart from the data that was analyzed) is the set of simulations, described in Sect. 4, which are used to estimate the various bias terms and characterize the estimator scatter.

We now proceed to describe these alternative pipelines in somewhat more detail. We note first two common aspects of all three alternative pipelines, which differ from our baseline results. First, rather than accounting for the Galaxy and point-source mask in the filter function as is done for our baseline method, these methods use a completely diagonal filter, which is the same as the form (given in Eq. (10)) that our baseline filtering asymptotes to far from the mask boundaries. Instead, as we will discuss, the alternative pipelines deal with the mask using different choices of data preparation and selection. A second difference common to all three pipelines with our baseline approach is that when computing the disconnected noise bias ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{0}}$, these methods do not use the two-point expression of Eq. (17), but rather the approximation to it based on Eq. (21), evaluated using an estimate of the data power spectrum. Because this correction is quite large, the agreement of results that use this alternative method to calculate it is reassuring. We now proceed to describe the individual methods in more detail; the common motivation in the development of each has been the reduction of the sharp gradients induced when masking, which can induce a mean-field several orders of magnitude larger than the lensing signal at low multipoles (as discussed in Appendix C). Each method takes a different approach to mitigating this mask effect, as discussed below:

• 1.

  The method iso consists of applying the standard quadratic lensing estimator to the sky map after multiplying by an apodized Galactic mask, and filling point source holes using local constrained Gaussian realizations of the CMB signal+noise. The mask mean-field is proportional to the power spectrum of the mask, and so as apodization smooths the mask boundary (suppressing its power spectrum on small scales), it correspondingly reduces the mean-field significantly. The combination of apodization and source filling makes this estimator very fast to apply to simulations but does require an involved set of correction terms and f_sky factors. Our implementation and calculation of this method is described in detail in Benoit-Levy et al. (2013).

• 2.

  The metis method consists of inpainting the Galactic mask as well as the point-source holes, using the sparse-inpainting algorithm described in Abrial et al. (2007, 2008). The resulting map resembles a full-sky CMB map, and therefore has no mask mean-field contribution. The inhomogeneous noise and beam-induced mean-fields do still have to be corrected however. Our implementation is based on that described in Perotto et al. (2010), with several improvements. In Perotto et al. (2010) lens reconstruction was performed on the inpainted map and then analyzed on the full sky. However, further inspection has revealed that there are some spurious features in the lens reconstruction, localized to the inpainted region inside the Galactic mask. This is likely due to the inhomogeneous noise in Planck that was ignored in previous work and cannot be reproduced by the inpainter. We therefore remask the full-sky lens reconstruction with an apodized Galactic mask (as in Eq. (14)) to remove these regions from our analysis. We follow the same procedure when evaluating the analytical expression for the ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{N}\mathrm{0}}$ bias, prewhitening and then applying an apodized Galactic mask to the inpainted temperature multipoles to estimate their power spectrum. Small residual biases are corrected using the same ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{MC}}$ procedure used in the main method.

• 3.

  The patches method avoids the Galactic mask completely by cutting the sky into a collection of 410 small overlapping 10deg × 10deg patches centered on the locations of N_side = 8 HEALPix pixels, which are then analyzed under the flat-sky approximation. Our implementation of this method is described in Plaszczynski et al. (2012). As with the iso method, point-source holes are filled using constrained Gaussian realizations. The patches are extracted from a pre-whitened CMB map, and apodized with a Kaiser-Bessel window function. The Fourier modes in each patch are fitted in real space using a fast Fourier-Toeplitz algorithm. No mean-field correction is applied. Residual biases due to noise inhomogeneity are removed using a ${\mathrm{\Delta }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}|}_{\mathit{MC}}$ correction that is found to be small. The patches method has been particularly useful in the early stages of our analysis, to identify outliers caused by unmasked point sources.

As can be seen in Fig. 21, all three of these methods are in good agreement with the results of our baseline method, providing reassurance that our results are insensitive to the precise details of our data filtering and reconstruction methodology.

8. Conclusions

The Planck maps have unprecedented sensitivity to gravitational lensing effects. We see significant and consistent measurements of lensing for each of the high-resolution CMB channels at 100, 143 and 217 GHz. Even the noisiest channel that we have considered, 100 GHz, provides a 10σ detection of lensing, which is greater than all previous detections. Our fiducial lens reconstruction, based on a minimum-variance combination of the 143 and 217 GHz channels does even better, with a detection of lensing (relative to the null hypothesis of no lensing) at a significance of greater than 25σ. Notably, the noise on our reconstruction is low enough that it is no longer the limiting source of noise for many correlations with large-scale structure catalogues (several examples of which we have given in Sect. 6.3). This marks a shift for CMB lensing, from the detection regime into that of standard cosmological probe. Our lensing potential map is publicly available, and we look forward to the uses which may be found fo it.

The percent-level Planck lensing potential measurement pushes into the realm of precision cosmology, and requires careful validation tests that we have performed in Sect. 7. Our fiducial likelihood, based on the 40 ≤ L ≤ 400 range that is most sensitive to lensing, passes all of the tests that we have performed at an acceptable level. Most importantly, it is consistent with individual frequency reconstructions, with more aggressive or more conservative masking, and with results obtained on more rigorously component-separated maps. Our measurement of the lensing potential power spectrum ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ is broadly consistent with the ΛCDM expectations of Planck Collaboration XVI (2014), although there are some shape and amplitude tensions that lead to smaller-than-forecasted improvements on parameters such as the sum of neutrino masses ∑ m_ν. As we saw in Sect. 6.1, our estimate of the potential power spectrum has direct cosmological implications, where it significantly improves CMB power-spectrum-only constraints on curvature (or alternatively, dark energy). Our lensing likelihood also breaks the A_s − τ degeneracy, allowing a measurement of the optical depth to reionization which is a useful independent cross-check on the values determined with large angular-scale polarization measurements.

CMB lensing science has undergone a remarkable development in the previous two years (Das et al. 2011; Sherwin et al. 2011; van Engelen et al. 2012). The next year is expected to bring even greater improvements. The 2500 deg² lensing survey of the South Pole Telescope is expected have comparable statistical power to the one reported here, with a highly complementary weighting in L-space that is more sensitive to smaller angular scale lenses. There is also additional Planck data, which will improve on the measurement reported here considerably. The results of this paper are based on the nominal-mission Planck maps, collected during the first fifteen months of science operations. The cryogens used to maintain the 100 mK cooling stage for the Planck HFI instrument lasted considerably longer than the nominal-mission lifetime, to a full-mission duration of approximately thirty months. This additional data, set to be released in 2014, results in map noise levels that are roughly a factor two lower in power than those used here. As instrumental noise constitutes approximately half of our error budget, we expect a corresponding decrease of approximately 25% in the variance of our lensing potential map as well as the uncertainties on ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$. The Planck polarization data, also set for release in 2014, will provide an additional powerful probe of lensing.

---

Acknowledgments

The development of Planck has been supported by: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MICINN, 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); and PRACE (EU). A description of the Planck Collaboration and a list of its members, including the technical or scientific activities in which they have been involved, can be found at http://www.sciops.esa.int/index.php?project=planck&page=Planck_Collaboration. Some of the results in this paper have been derived using the HEALPix package (Górski et al. 2005). This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the US Department of Energy under Contract No. DE-AC02-05CH11231, as well as of the IN2P3 Computer Center (http://cc.in2p3.fr) and the Planck-HFI data processing center infrastructures hosted at the Institut d’Astrophysique de Paris (France) and financially supported by CNES. We acknowledge support from the Science and Technology Facilities Council [grant number ST/I000976/1].

References

Appendix A: Null tests

In this appendix we present several generic null tests of our lens reconstruction. In Appendix A.1 we perform tests using lens reconstructions on the two individual six-month full-sky surveys contained in the first Planck data release. In Appendix. A.2 we perform the curl-mode null test which is standard practice in Galaxy lensing. Finally, in Appendix A.3 we will perform several tests on the isotropy and Gaussianity of our MV lensing potential map.

Appendix A.1: Survey consistency

The Planck nominal mission consists of approximately two full-sky six-month surveys. Our fiducial analysis is based on maps made from the full data range, however we may also analyze the first two surveys separately to increase our sensitivity to effects such as beam asymmetry, correlated noise along the scan direction, and zodiacal light. In Fig. A.1 we plot several consistency tests exploiting the breakdown of the data into separate surveys, for both 143 GHz and 217 GHz. We test for gross contamination of either survey by estimating the lensing power spectrum on each individually and then difference the resulting $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$, finding results consistent with the scatter expected from simple signal+noise simulations. To test the noise model more thoroughly, we also run lens reconstruction on a map constructed by taking the half-difference of the two surveys, filtering it in the same way as the full nominal-mission map. As can be seen in Fig. A.1, the power spectrum of this reconstruction null test is small compared to the measurement errors on $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$. Finally, to isolate the effect of systematics which couple to the CMB and which flip sign between surveys (such as odd moments of the beam asymmetry), we take the cross-spectrum of lens reconstruction estimates obtained from the half-difference of the two survey maps with those obtained from the half-sum. Again, this null test shows residuals that are negligibly small compared to our measurement errors.

Appendix A.2: Curl-mode reconstruction

In addition to the standard gradient lensing potential φ sourced by large-scale structure, we may also use the Planck maps to reconstruct a curl-mode lensing field ψ. The weight function associated with the curl field is given by (Namikawa et al. 2012) $\begin{array}{ccc}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\psi }}\mathrm{=}\mathrm{-}\sqrt{\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}}\sqrt{\mathit{L}\mathrm{(}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathit{\ell }}_{\mathrm{1}}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}}& & \\ \mathrm{\times }{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{TT}}\left(\frac{\mathrm{1}\mathrm{-}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathit{L}}}{\mathrm{2}}\right)\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\mathrm{+}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{\leftrightarrow }{\mathit{\ell }}_{\mathrm{2}}\mathrm{)}\mathit{.}& & \end{array}$(A.1)Using this weight function, we may construct an estimator for the curl-mode power spectrum, completely analogous to that for φ. The response function of the curl estimator to point source shot noise is zero, and so we neglect this correction. Consistency of the reconstructed curl mode with zero is a standard systematic test in weak lensing, as the curl field due to gravitational lensing should be negligible (Hirata & Seljak 2003b). One caveat for this test in the CMB lensing context is that the variance of the curl mode reconstruction is still affected by the presence of gravitational lenses, by means of an “N⁽¹⁾”-type (Kesden et al. 2003) contraction in the trispectrum (van Engelen et al. 2012). Our curl-mode power spectrum estimates are plotted in Fig. A.2, after correction for the “N⁽¹⁾”-type bias using Eq. (18). This bias is also given explicitly for the curl mode in Benoit-Levy et al. (2013). We do not see evidence for anomalous power in the curl reconstruction, although we do find in Sect. C.4 that the curl estimate is contaminated by pixelization effects at a level too small to be seen in this plot.

Appendix A.3: Statistical tests at the map level

In this section we present several sanity tests of our lensing estimates at the map level.

1-point PDF of :

in the left panel of Fig. A.3 we plot the histogram, or 1-point function, of our unnormalized potential estimate for the MV reconstruction, bandlimited to 10 ≤ L ≤ 2048. As can be seen from the histogram, the map is slightly non-Gaussian, as might be expected for an estimator that is based on the product of two Gaussian random fields. In the right panel of Fig. A.3 we plot the same histogram for a simulated MV reconstruction, which shows a comparable level of non-Gaussian structure. In Appendix D we show that this small non-Gaussianity is negligible for our lensing likelihood.

Power spectra on patches:

to test for gross contamination of our reconstruction from any localized region of the sky, we use the patches method (described in Sect 7.5) to build an estimate of the lensing power spectra on a collection of approximately 400 10deg × 10deg tiles. We fill point-source holes with local constrained Gaussian realizations of the CMB signal. The patches were chosen to have no overlap with the Galaxy mask.

The distribution of the power spectra estimated for each patch can be used as a diagnostic for deficiency in the mask and dust correction. Significantly contaminated but localized regions would appear as outliers. This procedure successfully detected several point sources that were overlooked in early versions of our source masks. In Fig. A.4 we present results of this analysis using our final Galaxy and point source masks for 143 GHz and 217 GHz.

Isotropy on patches:

we can also use the same data to build a coarse test of the isotropy of the reconstructed lensing map. We perform lens reconstruction on patches which are oriented along the local spherical basis (e_θ,e_φ), and then stack the corresponding 2D power spectrum estimates to build a bi-dimensional power spectrum for the 143 GHz and 217 GHz maps. This is shown in Fig. A.5, where we have additionally applied smoothing based on an undecimated wavelet transform to reduce the visual noise. The overplotted contour lines are placed at levels given by the azimuthally averaged power spectrum for L = 100, 200, 300 and 400.

Appendix B: Filtering

In this section, we discuss our filtering methodology, that takes us from a set of maps ${\mathit{T}}_{\mathit{p}}^{\mathrm{obs}\mathit{,\nu }}$ indexed by frequency ν and pixel p to the multipoles that are fed into the quadratic estimators used to measure lensing and possible contaminants. For an observed CMB that consists of Gaussian signal and noise, the optimal filtering (in the minimum-variance sense) for the construction of any quadratic anisotropy estimator is “C^-1” filtering, in which the observed map is beam-deconvolved and then multiplied by the inverse of a signal+noise covariance matrix. For the purposes of lens reconstruction, the Gaussian assumption above is quite accurate for Planck. The Planck map noise is well approximated as Gaussian (Planck Collaboration VI 2014). The observed CMB signal is less-obviously well approximated as Gaussian, with the non-Gaussianity due to lensing detected at high significance, however for lens reconstruction at Planck noise levels the difference in signal-to-noise using a filter constructed in the Gaussian limit versus a more complicated “delensing” filter that properly accounts for lensing non-Gaussianity has been shown to be negligible (Hirata & Seljak 2003a). Other sources of non-Gaussianity, such as unresolved point sources, are completely subdominant to either primary CMB or noise at all of the multipoles probed by our lens reconstruction, and so we neglect them in the construction of our filters.

For a single map, with a symmetric beam transfer function, full-sky coverage and homogeneous map noise levels, the C^-1 filter is given by $\begin{array}{ccc}\mathit{T̅}\mathit{\ell m}\mathrm{=}\frac{{\mathit{F}}_{\mathit{l}}}{{\mathit{B}}_{\mathit{l}}}\sum _{\mathit{p}\mathrm{=}\mathrm{1}}^{{\mathit{n}}_{\mathrm{pix}}}\frac{\mathrm{4}\mathit{\pi }}{{\mathit{n}}_{\mathrm{pix}}}{\mathit{Y}}_{\mathit{\ell m}}^{\mathrm{\ast }}\mathrm{\left(}{{n̂}}_{\mathit{p}}\mathrm{\right)}{\mathit{T}}_{\mathit{p}}^{\mathrm{obs}}\mathit{,}& & \end{array}$(B.1)where gives the value of the spherical harmonic at the center of pixel p, B_l gives the beam transfer function, and F_l is given by ${\mathit{F}}_{\mathit{l}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\mathrm{+}{\mathrm{(}\frac{\mathit{\pi }}{\mathrm{180}\mathrm{\times }\mathrm{60}}\mathit{\sigma }{{\mathit{B}}_{\mathit{\ell }}^{-1}}^{\mathrm{)}}}^{\mathrm{2}}}\mathit{,}$(B.2)where σ is the map noise level in units of μK arcmin. The maps we use must always be masked to remove bright Galactic foregrounds and point sources, and so we do not ever use the full-sky filter above, although as we will discuss the filter that we do use asymptotes to this form in regions far from the mask, and the simple form F_l of the filter function given in Eq. (B.2)will prove useful in many analytical calculations.

For the more realistic case of inhomogeneous map noise levels, with pixel-space noise correlations and a sky cut, the construction of is more involved. To account for the sky cut, we can take the noise covariance to be infinite for masked pixels. Reformulating the C^-1 filter using the inverse of the noise covariance matrix ${\mathit{N}}_{\mathit{p}{\mathit{p}}^{\mathrm{\prime }}}^{-1}$ in pixel space (defined to be zero if either p or p′ is a masked pixel), and allowing for the joint-filtering of multiple input maps (with no noise correlations between them) we have $\mathit{T̅}\mathit{\ell m}\mathrm{=}{\mathrm{(}{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}}^{\mathrm{)}}}^{-1}\sum _{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\sum _{\mathit{p}{\mathit{p}}^{\mathrm{\prime }}}\sum _{\mathit{\nu }}{\mathrm{𝒞}}_{\mathit{\ell m,}{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{-1}{\mathrm{𝒴}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{\mathrm{\ast }\mathit{p,\nu }}{\mathit{N}}_{\mathit{p}{\mathit{p}}^{\mathrm{\prime }}}^{\mathrm{-}\mathrm{1}\mathit{,\nu }}{\mathit{T}}_{{\mathit{p}}^{\mathrm{\prime }}}^{\mathrm{obs}\mathit{,\nu }}\mathit{.}$(B.3)where ν indexes the maps which are being combined and the pointing matrix ${\mathrm{𝒴}}_{\mathit{\ell m}}^{\mathit{p,\nu }}\mathrm{\equiv }{\mathit{B}}_{\mathit{\ell }}^{\mathit{\nu }}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\hat{{n}}{}_{\mathit{p}}\mathrm{\right)}\mathrm{4}\mathit{\pi }\mathit{/}{\mathit{n}}_{\mathrm{pix}}$ gives the value of the spherical harmonic at the center of pixel p, convolved with the appropriate beam+pixel transfer function ${\mathit{B}}_{\mathit{\ell }}^{\mathit{\nu }}$. The matrix is given by ${\mathrm{𝒞}}_{\mathit{\ell m,}{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{\equiv }{\mathrm{(}{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}}^{\mathrm{)}}}^{-1}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{+}\sum _{\mathit{p,}{\mathit{p}}^{\mathrm{\prime }}}\sum _{\mathit{\nu }}{\mathrm{𝒴}}_{\mathit{\ell m}}^{\mathrm{\ast }\mathit{p,\nu }}{\mathit{N}}_{\mathit{p}{\mathit{p}}^{\mathrm{\prime }}}^{\mathrm{-}\mathrm{1}\mathit{,\nu }}{\mathrm{𝒴}}_{{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{{\mathit{p}}^{\mathrm{\prime }}\mathit{,\nu }}\mathit{.}$(B.4)The filter of Eq. (B.3)is the one that we will use, with a slightly suboptimal choice for the noise covariance matrix N_pp′ that is described below.

The Planck map noise has several sources of pixel-pixel correlations for p ≠ p′. On small scales, there are correlations along the scan direction induced by the deconvolution of the bolometer time response (Planck Collaboration VII 2014). On large scales, there are correlations between scan rings, due to residual 1 /f noise fluctuations not cancelled by the baseline subtraction procedure (Planck Collaboration VI 2014). Accounting for these pixel correlations would be computationally intractable, as the full covariance matrix N_pp′ (an ${\mathit{n}}_{\mathrm{pix}}^{\mathrm{4}}$ sized object) is too large to work with on a modern computer, and so we neglect them in our noise filter. The remaining structure in the noise matrix is due to the variable coverage in different regions of the sky. Noise levels are lowest near the Ecliptic poles, which are visited more frequently than the region around the Ecliptic plane due to the Planck scan strategy. Although the variable noise level is tractable to include in our filtering procedure (and is indeed included in the C^-1 filtering of e.g., Planck Collaboration XXIV 2014), we have chosen to neglect this aspect of the noise, using instead a fixed effective noise level. We use an inverse-noise map in pixel space given by ${\mathit{N}}_{\mathit{p}}^{\mathrm{-}\mathrm{1}\mathit{,\nu }}\mathrm{=}{\left(\frac{\mathrm{180}\mathrm{\times }\mathrm{60}\mathrm{\times }{\mathit{M}}_{\mathit{p}}}{\mathit{\pi }\mathrm{\times }{\mathit{\sigma }}_{\mathit{\nu }}}\right)}^{\mathrm{2}}\left(\frac{\mathrm{4}\mathit{\pi }}{{\mathit{n}}_{\mathrm{pix}}}\right)\mathit{,}$(B.5)where M_p is a sky mask and σ_ν is an average map noise level in units of μK arcmin. In addition to this noise weighting, we also project out the five modes corresponding to the monopole, dipole, and 857 GHz Planck map (as a simple dust template). The noise level for each of these modes is taken to infinity using the Woodbury formula; putting each of these modes into a map ${\mathit{T}}_{\mathit{p}}^{\mathit{i}}$ indexed by i, our total noise covariance matrix is given in terms of the inverse noise map above by ${\mathit{N}}_{\mathit{p}{\mathit{p}}^{\mathrm{\prime }}}^{\mathrm{-}\mathrm{1}\mathit{,\nu }}\mathrm{=}{\mathit{\delta }}_{\mathit{p}{\mathit{p}}^{\mathrm{\prime }}}{\mathit{N}}_{\mathit{p}}^{\mathrm{-}\mathrm{1}\mathit{,\nu }}\mathrm{-}\sum _{\mathit{ij}}{\mathrm{𝒯}}_{\mathit{ij}}^{-1}{\mathit{T}}_{\mathit{p}}^{\mathit{i}}{\mathit{T}}_{{\mathit{p}}^{\mathrm{\prime }}}^{\mathit{j}}{\mathit{N}}_{\mathit{p}}^{\mathrm{-}\mathrm{1}\mathit{,\nu }}{\mathit{N}}_{{\mathit{p}}^{\mathrm{\prime }}}^{\mathrm{-}\mathrm{1}\mathit{,\nu }}\mathit{,}$(B.6)where is the overlap matrix between templates given by ${\mathrm{𝒯}}_{\mathit{ij}}\mathrm{=}\sum _{\mathit{p}}{\mathit{T}}_{\mathit{p}}^{\mathit{i}}{\mathit{T}}_{\mathit{p}}^{\mathit{j}}{\mathit{N}}_{\mathit{p}}^{\mathrm{-}\mathrm{1}\mathit{,\nu }}\mathit{.}$(B.7)Our motivation for taking a fixed noise level is that with this approach, in regions sufficiently far from the mask boundary, our filter asymptotes to the diagonal form of Eq. (B.1). This means that the normalization of our lensing estimates can be well-approximated analytically, which is very useful for the propagation of systematic effects, and also that the normalization of our lensing estimates does not vary across the sky with noise level, which simplifies cross-correlation analysis. Our C^-1 filter is therefore optimally accounting for masking effects, but not for noise correlations and inhomogeneity. We estimate the suboptimality of neglecting these noise properties by calculating the quantity $\frac{\mathrm{(}\mathit{S}\mathit{/}\mathit{N}{\mathrm{)}}^{\mathit{use}}}{\mathrm{(}\mathit{S}\mathit{/}\mathit{N}{\mathrm{)}}^{\mathit{opt}}}\mathrm{=}\frac{{\mathrm{(}{{\mathit{R}}_{\mathit{L}}^{\mathit{\phi \phi ,}\mathit{use}}}^{\mathrm{)}}}^{\mathrm{2}}}{{\mathrm{(}{{\mathit{R}}_{\mathit{L}}^{\mathit{\phi \phi ,}\mathit{opt}}}^{\mathrm{)}}}^{}}{\left(\sum _{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}}\frac{\mathrm{1}}{\mathrm{2}}{\left|{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{\phi }}\right|}^{\mathrm{2}}\frac{{\mathrm{(}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{use}}{{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{use}}}^{\mathrm{)}}}^{\mathrm{2}}}{{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{opt}}{\mathit{F}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{opt}}}\right)}^{-1}\mathit{,}$(B.8)where ${\mathit{F}}_{\mathit{\ell }}^{\mathit{opt}}$ is the optimal filter and ${\mathit{F}}_{\mathit{\ell }}^{\mathit{use}}$ is the suboptimal filter that we have actually used. This equation gives the S/N loss as a function of lens multipole L, however in practice we find that the L dependence is small enough that it suffices to quote a single average loss. To estimate the degradation due to ignoring noise correlations we set ${\mathit{F}}_{\mathit{\ell }}^{\mathit{opt}}\mathrm{=}\frac{\mathrm{1}}{{\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}\mathrm{+}{\mathit{B}}_{\mathit{\ell }}^{\mathrm{-}\mathrm{2}\mathit{,\nu }}{\mathit{N}}_{\mathit{L}}^{\mathit{TT}}}\mathit{,}$(B.9)where ${\mathit{N}}_{\mathit{L}}^{\mathit{TT}}$ is the power spectrum of the map noise. We find that the degradation due to neglect of noise correlations is small; less than 2% for all L ≤ 2048 at 100 GHz, and less than 0.1% at 143 and 217 GHz. To calculate the degradation due to ignoring noise inhomogeneity, we determine the map noise level in the 3072 regions corresponding to N_side = 16 HEALPix pixels, take ${\mathit{F}}_{\mathit{\ell }}^{\mathit{opt}}$ using Eq. (B.2)with the local noise level, and estimate a resulting S/N degradation using Eq. (B.8). The neglect of noise inhomogeneity is the dominant suboptimality of our filtering, although it is still small. We find an average S/N loss (averaged over the entire sky) of approximately 4% at 100, 143, and 217 GHz, consistent with the simulation-based results of Hanson et al. (2009). We take this loss as justified, given the simpler normalization properties of our lensing estimates when neglecting variations in the map noise level.

Appendix C: Mean-fields

As discussed in Sect. 2, the quadratic lensing estimators that we use are designed to detect statistical anisotropy induced by lensing. There are a number of non-lensing sources of statistical anisotropy that can mimic the lensing signal to some extent. In our analysis, the effects we consider are

• 1.

  The application of a sky mask, which introduces sharp gradients that may be misinterpreted as lensing.

• 2.

  Noise inhomogeneity, which causes the overall power to fluctuate across the sky and can resemble the convergence component of lensing.

• 3.

  Beam asymmetry, which smears the fluctuations more along one direction than another and can mimic the shear component of lensing.

• 4.

  Pixelization, in which detector samples are accumulated into pixels, introduces a spurious deflection field on the pixel scale because the centroid of the hit distribution in each pixel does not necessarily lie at the pixel center.

In our analysis, we account for most of these effects with a corrective mean-field term, given by Eq. (9), which is determined using Monte-Carlo simulations. In this appendix, we will break this mean-field down into its constituent parts and discuss each in more detail. As an overview of the results in this section, in Fig. C.1 we plot estimates for the three largest mean-fields, due to masking, noise inhomogeneity, and beam asymmetry at 143 GHz (100 and 217 GHz are qualitatively similar). These mean-fields all have most of their contributions on very large scales, dictated by the coherency of the scan strategy in the case of beam asymmetry and noise inhomogeneity, and of the large-scale nature of the Galactic foregrounds in the case of the sky mask. We note also that the mask mean-field is concentrated near the edges of the Galaxy cut, and so for our power spectrum analysis (where we apodize the estimated φ map, removing the contributions near the edges of the mask), the effective mean-field is significantly reduced.

Our discussion will focus on constructing simple models for each source of mean-field. Following Hanson et al. (2010), we will identify each of the individual contributions to the mean-field with a tracer z_LM that sources a contribution to the CMB covariance matrix given by $\mathrm{\Delta }\mathrm{⟨}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{m}}_{\mathrm{2}}}{\mathit{T}}_{{\mathit{\ell }}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}}^{\mathrm{\ast }}\mathrm{⟩}\mathrm{=}\sum _{\mathit{LM}}{\mathit{z}}_{\mathit{LM}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}\left(\begin{array}{c}\\ {\mathit{\ell }}_{\mathrm{1}}& {\mathit{\ell }}_{\mathrm{2}}& \mathit{L}\\ {\mathit{m}}_{\mathrm{1}}& {\mathit{m}}_{\mathrm{2}}& \mathrm{-}\mathit{M}\end{array}\right){\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{z}}\mathit{,}$(C.1)where ${\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{z}}$ is a weight function describing how z_LM couples multipoles. Such a contaminant leads to a bias for the standard lensing estimator given by $\mathit{\phi ̂}\begin{array}{c}\mathrm{MF}\\ \mathit{LM}\end{array}\mathrm{=}\frac{{\mathrm{ℛ}}_{\mathit{LM}}^{\mathit{\phi z}}}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi \phi }}}{\mathit{z}}_{\mathit{LM}}\mathit{,}$(C.2)where the response function ${\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi z}}$ is defined in Eq. (12). The analytical forms for the mean-fields that we present here are used in Sect. 7.4 to construct “bias hardened” estimators that have smaller mean-field corrections, as proposed by Namikawa et al. (2013).

Appendix C.1: Mask

The application of a mask $\mathit{M}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ to the beam-convolved sky, followed by subsequent beam deconvolution, can be modelled as a statistical anisotropy “M”, with ${\mathit{M}}_{\mathit{LM}}\mathrm{=}{}^{\mathrm{\int }}\mathrm{d}\hat{{n}}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{M}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ and weight function (Namikawa et al. 2013) $\begin{array}{ccc}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{M}}\mathrm{=}\sqrt{\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}}\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)& & \\ \mathrm{\times }\left({\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{TT}}\frac{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}\mathrm{+}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{TT}}\frac{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}}\right)\mathrm{·}& & \end{array}$(C.3)The mask mean-field using our fiducial Galaxy and point-source analysis mask (described in Sect. 3) is plotted in Fig. C.1. At very low-L, the mean field induced by this mask can be much larger than the reconstruction noise level. At L = 2 for example, it can be up to 10⁴ times greater than the reconstruction noise level, although the amplitude of the mean-field is a very steep function of L.

The mask mean-field is proportional to ${\mathit{C}}_{\mathit{\ell }}^{\mathit{TT}}$, and so with percent-level agreement between the power spectrum of our simulations we use to determine the mean-field and the data power spectrum (discussed in Sect. 4), we can expect errors on the power spectra of the mask mean-field at the level of one part in 10⁴. At L = 2 for 143 and 217 GHz this corresponds to a level of $\mathrm{\left[}\mathit{L}\mathrm{\right(}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{\approx }\mathrm{4}\mathrm{\times }{\mathrm{10}}^{-7}$, comparable to the lensing potential power spectrum. There are several good ways to mitigate the effect of the mask mean-field. By apodizing the mask before lens reconstruction, the gradients at the mask boundary can be greatly reduced, which significantly suppresses the mean-field (van Engelen et al. 2012; Benoit-Levy et al. 2013). At the map level, the only issue with such an apodization for our purposes is that it results in a position-dependent normalization which can complicate cross-correlation analyses. A related approach, which we have taken in our power spectrum analysis, is to apodize the lens reconstruction after performing the reconstruction. This amounts to removing the regions of near the mask boundary where the mean-field has all of its features. This “avoidance” approach can also be used easily in cross-correlation analyses, as in Hirata et al. (2004). A third approach is that of Namikawa et al. (2013), who show that it is possible to construct “bias-hardened” lensing estimators which have smaller mean-field contributions. We have already shown consistency tests based on this approach in Sect. 7.4. In the low-L bias-hardened comparisons of Fig. 20, we saw discrepancies with the standard estimator of $\mathrm{\left[}\mathit{L}\mathrm{\right(}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}{\mathrm{]}}^{\mathrm{2}}{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}\mathrm{\approx }\mathrm{4}\mathrm{\times }{\mathrm{10}}^{-7}$ and 7 × 10^-7, precisely the level expected above.

For our publicly released lensing map, we have chosen simply to high-pass filter our map to L> 10, above which we believe that mean-field errors are subdominant to the reconstruction noise.

Appendix C.2: Instrumental noise

The instrumental noise of Planck is rich in statistically-anisotropic structure, as discussed in Planck Collaboration VI (2014), which is sourced by the following

• 1.

  Striping effects, due to residual noise correlations on long time scales after destriping.

• 2.

  Correlations along the scan direction induced by the deconvolution of the detector time response and the low-pass rolloff filter.

• 3.

  Large variations in the pixel noise level due to the larger density of pixel hits near the Ecliptic poles.

All of these effects can be avoided by forming estimators using pairs of data with independent noise realizations, as we have done in Sect. 7.3, although our baseline results do not use such cross-correlations. Here we discuss the size of the noise mean-field. Items (1) and (2) are difficult to study analytically. For white noise with a variance in each pixel given by $\mathit{N}\mathrm{\left(}\hat{{n}}{}_{\mathit{p}}\mathrm{\right)}$, item (3) can be identified with a statistical anisotropy “N” defined by Hanson et al. (2009)$\begin{array}{ccc}{\mathit{N}}_{\mathit{LM}}& \mathrm{=}& \int \mathrm{d}\hat{{n}}\mathit{p}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\hat{{n}}\mathit{p}\mathrm{\right)}\mathit{N}\mathrm{\left(}\hat{{n}}\mathit{p}\mathrm{\right)}\\ {\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{N}}& \mathrm{=}& \sqrt{\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}}\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\frac{\mathrm{1}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}}\frac{\mathrm{1}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}}\mathit{,}\end{array}$(C.4)Note that the shape of this bias is dictated by the hit count distribution, which is known quite well, rather than the noise level, which is more uncertain. The shape of the noise mean-field is plotted in Fig. C.1. As with the mask, it is large compared to our signal and noise but falls off quickly with multipole. For L< 10, given that there are known issues with our noise modelling at the 5% level in power (as discussed in Sect. 4), it can be important to construct quadratic estimators from correlations of maps with independent noise realizations (apart from our bias-hardened estimator tests in Sect. 7.4, this is also done for L = 1 in Planck Collaboration XXVII 2014), however at the L> 10 multipoles of our lensing map and power spectrum, uncertainties in the noise mean-field should be negligible.

Appendix C.3: Beam asymmetry

Planck observes the sky after convolution with a “scanning beam”, which captures its effective response to the sky as a function of displacement from the nominal pointing direction. Asymmetry in the shape of this beam (for example ellipticity) means that the fluctuations in the map are smeared depending on the beam orientation, which can mimic the shearing effects of gravitational lensing. Here we present semi-analytical estimates of this beam asymmetry mean-field, and discuss how simulations of beam asymmetry are included in our simulations to correct for it. As with the noise inhomogeneity, we will see that the large-scale coherency of the Planck scan strategy means that beam asymmetries only mimic large-scale modes of φ_LM.

Decomposing the scanning beam into harmonic coefficients B_ℓm, each time-ordered data (TOD) sample can be modelled as (neglecting the contribution from instrumental noise, which is independent of beam asymmetry) ${\mathit{T}}_{\mathit{i}}\mathrm{=}\sum _{\mathit{\ell ms}}{\mathrm{e}}^{\mathrm{-}\mathit{is}{\mathit{\alpha }}_{\mathit{i}}}{\mathit{B}}_{\mathit{\ell s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}{\mathit{s}}_{}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathit{,}{\mathit{\phi }}_{\mathit{i}}\mathrm{\right)}\mathit{,}$(C.5)where the TOD samples are indexed by i, and $\begin{array}{c}{\mathrm{˜}}_{}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}$ is the underlying sky signal. The spin spherical harmonic ${\mathit{s}}_{}{\mathit{Y}}_{\mathit{\ell m}}$ rotates the scanning beam to the pointing location (θ,φ), while the e^{− isα_i} factor gives it the correct orientation. On the small scales that lensing is sensitive to, it is a good approximation to assume that the procedure of mapmaking from TOD samples is just a process of binning: $\mathit{T}\mathrm{\left(}\hat{{n}}\mathit{p}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{\in }\mathit{p}}{\mathit{T}}_{\mathit{i}}\mathit{/}{\mathit{H}}_{\mathit{p}}\mathit{,}$(C.6)where $\hat{{n}}{}_{\mathit{p}}$ is the centroid of pixel p, and the sum is taken over all hits assigned to pixel p. H_p is the number of hits in pixel p. Here we will assume that the pointing for each observation of pixel p is identical to the pixel center $\hat{{n}}{}_{\mathit{p}}$. The purpose of this assumption is simply to isolate the effect of beam asymmetry from the (small) effect of pixelization. We will discuss the pixelization effect in Sect. C.4. Placing all hits at pixel centers, combining Eq. (C.5)and Eq. (C.6)we have $\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{s}}\mathit{w}\mathrm{(}\hat{{n}}\mathit{,}\mathrm{-}\mathit{s}\mathrm{)}\left[\sum _{\mathit{\ell m}}{\mathit{B}}_{\mathit{\ell s}}\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}{\mathit{s}}_{}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\right]\mathit{,}$(C.7)where the “scan strategy maps” $\mathit{w}\mathrm{\left(}\hat{{n}}\mathit{,s}\mathrm{\right)}$ are given by $\mathit{w}\mathrm{(}\hat{{n}}\mathit{p}\mathit{,}\mathrm{-}\mathit{s}\mathrm{)}\mathrm{=}\sum _{\mathit{i}\mathrm{\in }\mathit{p}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{s}{\mathit{\alpha }}_{\mathit{i}}}\mathit{/}{\mathit{H}}_{\mathit{p}}\mathit{.}$(C.8)The scan strategy maps are spin-s objects, with corresponding harmonic decomposition ${\mathit{s}}_{}{\mathit{w}}_{\mathit{LM}}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\hat{{n}}{\mathit{s}}_{}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{w}\mathrm{\left(}\hat{{n}}\mathit{,s}\mathrm{\right)}\mathit{.}$(C.9)We use Eq. (C.7)to include the effect of asymmetric beams in the simulations that we use to debias our lensing estimates. We can also propagate the effect of beam asymmetry to the lensing mean-field directly. Each mode s of the scan strategy can be identified with a statistical anisotropy “b^s” of the type defined in Eq. (C.1), with ${\mathit{b}}_{\mathit{LM}}^{\mathit{s}}\mathrm{=}{\mathit{s}}_{}{\mathit{w}}_{\mathit{LM}}$ and $\begin{array}{ccc}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{{\mathit{b}}^{\mathit{s}}}\mathrm{=}\sqrt{\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}}& & \\ \mathrm{\times }\left[\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\frac{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}\mathit{s}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}\mathrm{\ast }}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{TT}}\mathrm{+}\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\frac{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}\mathit{s}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}\mathrm{\ast }}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{2}}}^{\mathit{TT}}\right]\mathit{.}& & \end{array}$(C.10)In Fig. C.1 we have plotted the expected full-sky bias for the best-fit Planck beams, which are measured from planet observations in Planck Collaboration VII (2014). There are considerable differences in the mean-field between frequencies, as well as the steep scale dependence, which is dictated by the large-scale coherency of the scan strategy.

Appendix C.4: Pixelization

In the process of mapmaking, TOD samples are accumulated into bins defined by HEALPix pixels. The Planck HFI maps are generated at HEALPix resolution 11, with a typical size of , which is sufficiently close to the scale of deflections expected from gravitational lensing effects that this remapping of the observed CMB is a source of concern. To model this effect, we begin from Eq. (C.5), but treat the beam as symmetric to isolate the effect of pixelization from that of beam asymmetry. Binning hits into pixels, we then have $\begin{array}{ccc}\mathit{T}\mathrm{\left(}\hat{{n}}\mathrm{\right)}& \mathrm{=}& \sum _{\mathit{i}\mathrm{\in }\mathit{p}}\sum _{\mathit{\ell m}}{\mathit{B}}_{\mathit{\ell }\mathrm{0}}\begin{array}{c}\mathrm{˜}\\ {\mathit{T}}_{\mathit{\ell m}}\end{array}{\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}{\mathit{\theta }}_{\mathit{i}}\mathit{,}{\mathit{\phi }}_{\mathit{i}}\mathrm{\right)}\mathit{/}\mathit{H}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\\ & \mathrm{\approx }& \end{array}$(C.11)where $\hat{{n}}$ is the center of pixel p and in the second line we have made a gradient approximation to ${\mathit{Y}}_{\mathit{\ell m}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}$ about the pixel center, and used the effective deflection vector ${{d}}^{\mathrm{pix}}\mathrm{\left(}\hat{{n}}\mathit{p}\mathrm{\right)}\mathrm{=}\sum _{\mathit{i}\mathrm{\in }\mathit{p}}\mathrm{(}\hat{{n}}\mathit{i}\mathrm{-}\hat{{n}}{\mathit{p}}^{\mathrm{)}}\mathit{/}{\mathit{H}}_{\mathit{p}}\mathit{,}$(C.12)where $\hat{{n}}{}_{\mathit{i}}$ is the pointing for observation i, which is assigned to pixel p (with center at $\hat{{n}}{}_{\mathit{p}}$ and number of hits H_p). Implicit in this model is an assumption that on the arcminute pixelization scale the beam-convolved sky fluctuations are quite smooth, and can be modelled as a gradient. Pixelization therefore introduces an effective deflection field, given by d^pix, which biases our reconstruction of the lensing potential. We could in principle correct for the spurious lensing field which this introduces by including the remapping of Eq. (C.11)in the simulations used to estimate the mean-field correction, however as we will see it is quite small. We can decompose the spin-1 deflection field as ${\mathit{d}}_{\mathit{LM}}^{\mathrm{pix}}\mathrm{=}\mathrm{\int }{\mathrm{d}}^{\mathrm{2}}\hat{{n}}\hspace{0.17em}{\mathrm{1}}_{}^{}{\mathit{Y}}_{\mathit{LM}}^{\mathrm{\ast }}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathrm{[}{\mathit{d}}_{\mathit{\theta }}^{\mathrm{pix}}\mathrm{+}\mathrm{i}{{\mathit{d}}_{\mathit{\phi }}^{\mathrm{pix}}}^{\mathrm{]}}\mathrm{\left(}\hat{{n}}\mathrm{\right)}\mathit{,}$(C.13)where d_θ and d_φ are the components of d^pix along the θ and φ directions respectively. The associated weight function, following Eq. (C.1), is given by

$\begin{array}{ccc}{\mathit{W}}_{{\mathit{\ell }}_{\mathrm{1}}{\mathit{\ell }}_{\mathrm{2}}\mathit{L}}^{\mathit{d}}\mathrm{=}\mathrm{-}\sqrt{\frac{\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}{\mathit{\ell }}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{)}\mathrm{(}\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathrm{4}\mathit{\pi }}}{\mathit{C}}_{{\mathit{\ell }}_{\mathrm{1}}}^{\mathit{TT}}& & \\ \mathrm{\times }\sqrt{\frac{{\mathit{\ell }}_{\mathrm{1}}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{+}\mathrm{1}\mathrm{)}}{\mathit{L}\mathrm{(}\mathit{L}\mathrm{+}\mathrm{1}\mathrm{)}}}\frac{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{1}}}}{{\mathit{B}}_{{\mathit{\ell }}_{\mathrm{2}}}}\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\mathrm{+}\mathrm{(}{\mathit{\ell }}_{\mathrm{1}}\mathrm{\leftrightarrow }{\mathit{\ell }}_{\mathrm{2}}\mathrm{)}\mathit{.}& & \end{array}$(C.14)The pixelization deflection field can also be usefully split into gradient- and curl-type parts, as ${\mathit{d}}_{\mathit{LM}}^{\mathrm{pix}}\mathrm{=}{\mathit{d}}_{\mathit{LM}}^{\mathit{G,}\hspace{0.17em}\mathrm{pix}}\mathrm{+}\mathrm{i}{\mathit{d}}_{\mathit{LM}}^{\mathit{C,}\hspace{0.17em}\mathrm{pix}}\mathit{,}$(C.15)where the gradient and curl parts both satisfy the reality condition ${\mathit{d}}_{\mathit{LM}}^{\mathit{G}\mathit{/}\mathit{C,}\hspace{0.17em}\mathrm{pix}}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}{\mathit{d}}_{\mathit{L}\mathrm{-}\mathit{M}}^{\mathrm{\ast }\mathit{G}\mathit{/}\mathit{C,}\hspace{0.17em}\mathrm{pix}}\mathit{.}$(C.16)Only the gradient mode couples to the lensing estimator, while the curl part only couples to the curl-mode lensing estimator (discussed in Sect. A.2). The statistical anisotropy introduced by ${\mathit{d}}_{\mathit{G,}\hspace{0.17em}\mathrm{pix}}^{}$ is almost completely degenerate with that of CMB lensing, differing only by the fact that it operates on the beam-convolved CMB sky rather than the primordial one. In Fig. C.2 we plot the expected bias to the lensing power spectrum estimator, which is given by $\begin{array}{ccc}{\mathrm{\Delta }\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}|}_{\mathrm{pix}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}{\mathrm{(}{{\mathrm{ℛ}}_{\mathit{L}}^{\mathit{\phi d}}}^{\mathrm{)}}}^{\mathrm{2}}\sum _{\mathit{M}}{\left|{\mathit{d}}_{\mathit{LM}}^{\mathit{G}}\right|}^{\mathrm{2}}\mathit{,}& & \end{array}$(C.17)as well as the analogous bias ${\mathrm{\Delta }\mathit{Ĉ}\begin{array}{c}\mathit{\psi \psi }\\ \mathit{L}\end{array}|}_{\mathrm{pix}}$ for the curl null-test estimator. To test for the presence of the pixelization lensing field in our reconstructions, we cross-correlate the deflection fields ${\mathit{d}}_{\mathit{LM}}^{\mathit{G}\mathit{/}\mathit{C,}\hspace{0.17em}\mathrm{pix}}$ with our lensing estimates , . In Fig. C.2 we plot $\mathit{Ĉ}\begin{array}{c}\mathit{\phi d}\\ \mathit{L}\end{array}\mathrm{=}\frac{{\mathit{f}}_{\mathrm{sky}}^{-1}}{\mathrm{2}\mathit{L}\mathrm{+}\mathrm{1}}\sum _{\mathit{M}}\mathit{\phi ̂}\mathit{LM}{\mathit{d}}_{\mathit{LM}}^{\mathit{G}\mathrm{\ast }}\mathit{,}$(C.18)and similarly for ψ. We detect both ${\mathit{d}}_{\mathit{LM}}^{\mathit{G}}$ and ${\mathit{d}}_{\mathit{LM}}^{\mathit{C}}$ at greater than 15σ significance for all of our three main frequency bands, at the expected level. This detection is only possible because we know d_LM precisely from the satellite pointing. Except at the very highest multipoles, the gradient mode of the deflection lensing field is always significantly smaller than the CMB lensing potential. At L< 250, which provides most of the weight for the Planck lensing likelihood, the bias to the lensing potential is significantly less than 0.1%, while our statistical measurement uncertainty is at the 4% level.

Appendix D: Likelihood validation

In the construction of our lensing likelihood in Sect. 2.3 we have used several simplifying approximations to the statistics of the reconstructed lensing potential, which we discuss briefly here.

Uncorrelated binned power spectrum estimates: in Fig. D.1 we plot Monte-Carlo estimates of the bandpower correlation. These plots do not include the effect of beam uncertainties, errors in the point-source correction or cosmological uncertainty in the subtraction of the N⁽¹⁾ bias, which we include in our likelihood covariance matrix analytically. Even for adjacent bandpowers, the correlation coefficient is constrained to be less than 10% and we ignore such correlations in constructing the likelihood. As discussed in Sect. 5, however, we do include bandpower correlations from beam uncertainties and the point-source and N⁽¹⁾ correction in the likelihood.

Gaussianity of the binned power spectrum estimates: our estimator for the lensing potential power spectrum is a sum over the trispectrum of the observed CMB. There is no simple, rigorous argument to expect our bandpower estimates to be Gaussian distributed, however keeping in mind that we are collapsing the information in a multi-million pixel map down to eight eigenmodes, it is quite reasonable to expect that the central limit theorem will drive their distribution in this direction. With a set of simulated lens reconstructions, it is straightforward to assess the Gaussianity of our amplitude estimates. In Fig. D.2 we plot histograms of the bandpower coefficients from our set of fiducial simulations showing that the Gaussian approximation is good for the marginal distributions. Further tests of the assumption of a Gaussian likelihood are made in Schmittfull et al. (2013) for isotropic surveys. There it is shown that parameter constraints for a simple two-parameter model of ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ scatter across simulations in a manner consistent with the likelihood width in any realization.

No correlations between the lensing potential and temperature power spectra: lensing may be detected at a significance of 10σ as a smearing of the acoustic peaks in the temperature power spectrum; see Planck Collaboration XVI (2014). A potential worry when simply combining the temperature likelihood with the lensing power spectrum likelihood is that we are double counting the same information (Lewis 2005; Perotto et al. 2006; Smith et al. 2006) since we do not account for the correlated effect of sample variance of the lenses in the estimated $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$ and $\mathit{Ĉ}\begin{array}{c}\mathit{TT}\\ \mathit{\ell }\end{array}$ spectra. This effect is explored in detail in Schmittfull et al. (2013) where it is shown to derive from the connected part of the 6-point function for the lensed temperature anisotropies. They estimate for Planck that the maximum correlation between $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$ and $\mathit{Ĉ}\begin{array}{c}\mathit{TT}\\ \mathit{\ell }\end{array}$ for unbinned spectra is around 0.05%, and that ignoring this leads to a misestimation of the errors on the lensing amplitude A_L (see Sect. 6.1) in a joint analysis of only 2%. A more pressing concern is that the statistical errors on the reconstruction power are correlated with those for the temperature power spectrum due to sample variance of the CMB anisotropies. This effect is also discussed in detail in Schmittfull et al. (2013). The correlation derives from the fully disconnected part of the CMB 6-point function and leads to correlations of at most 0.2% between the unbinned spectra. However, these correlations are removed exactly by the realization-dependent aspect of the N⁽⁰⁾ debiasing procedure discussed in Sect. 2.2.

To see that the disconnected contribution to $\mathrm{cov}\mathrm{\left(}\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}\mathit{,}\mathit{Ĉ}\begin{array}{c}\mathit{TT}\\ \mathit{\ell }\end{array}\mathrm{\right)}$ vanishes, we note that with our data-dependent N⁽⁰⁾ subtraction the data enters the $\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}$ in the form⁷$\begin{array}{ccc}\mathit{Ĉ}\begin{array}{c}\mathit{\phi \phi }\\ \mathit{L}\end{array}& \mathrm{~}& \frac{\mathit{\partial }}{\mathit{\partial }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}}\mathrm{⟨}{\mathit{T}}_{\mathit{i}}{\mathit{T}}_{\mathit{j}}{\mathit{T}}_{\mathit{k}}{\mathit{T}}_{\mathit{l}}{\mathrm{⟩}}_{\mathrm{c}}(\mathit{T̅}\mathit{i}\mathit{T̅}\mathit{j}\mathit{T̅}\mathit{k}\mathit{T̅}\mathit{l}\\ & & \mathrm{-}\mathrm{[}{\mathit{C}}_{\mathit{ij}}^{-1}\mathit{T̅}\mathit{k}\mathit{T̅}\mathit{l}\mathrm{+}{\mathit{C}}_{\mathit{kl}}^{-1}\mathit{T̅}\mathit{i}\mathit{T̅}\mathit{j}\mathrm{-}\mathit{C̅}\begin{array}{c}-1\\ \mathit{ij}\end{array}\mathit{C̅}{\begin{array}{c}-1\\ \mathit{kl}\end{array}}^{\mathrm{]}}\\ & & \mathrm{-}\mathrm{[}{\mathit{C}}_{\mathit{ik}}^{-1}\mathit{T̅}\mathit{j}\mathit{T̅}\mathit{l}\mathrm{+}{\mathit{C}}_{\mathit{jl}}^{-1}\mathit{T̅}\mathit{i}\mathit{T̅}\mathit{k}\mathrm{-}\mathit{C̅}\begin{array}{c}-1\\ \mathit{ik}\end{array}\mathit{C̅}{\begin{array}{c}-1\\ \mathit{jl}\end{array}}^{\mathrm{]}}\\ & & \end{array}$(D.1)Here, is the inverse-variance filtered temperature, with covariance ${\mathit{C}}_{\mathit{ij}}^{-1}$, and summation over repeated indices is assumed. Equation (D.1) arises naturally for optimal trispectrum estimation of weakly non-Gaussian fields (e.g., Regan et al. 2010). The expectation value of the data combination on the right-hand side is simply the connected 4-point function, ⟨ T_iT_jT_kT_l ⟩ _c, but the data-dependent quadratic terms improve the variance of the estimator. The right-hand side of Eq. (D.1) is explicitly symmetric under regroupings of pairs of indices so that only the primary coupling of the trispectrum need be included. For lensing, the primary coupling factors such that $\frac{\mathit{\partial }}{\mathit{\partial }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}}\mathrm{⟨}{\mathit{T}}_{\mathit{i}}{\mathit{T}}_{\mathit{j}}{\mathit{T}}_{\mathit{k}}{\mathit{T}}_{\mathit{l}}{\mathrm{⟩}}_{\mathrm{c}}^{\mathrm{primary}}\mathrm{=}\sum _{\mathit{M}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}{\mathit{X}}_{\mathit{LM}}^{\mathit{ij}}{\mathit{X}}_{\mathit{L}\mathrm{-}\mathit{M}}^{\mathit{kl}}\mathrm{;}$(D.2)see Hu (2001); Hanson et al. (2011) for explicit expressions. The first set of terms in Eq. (D.1) combine to give $\begin{array}{ccc}& & \frac{\mathit{\partial }}{\mathit{\partial }{\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}}\mathrm{⟨}{\mathit{T}}_{\mathit{i}}{\mathit{T}}_{\mathit{j}}{\mathit{T}}_{\mathit{k}}{\mathit{T}}_{\mathit{l}}{\mathrm{⟩}}_{\mathrm{c}}^{\mathrm{primary}}\\ & & \mathrm{\times }\mathrm{(}\mathit{T̅}\mathit{i}\mathit{T̅}\mathit{j}\mathit{T̅}\mathit{k}\mathit{T̅}\mathit{l}\mathrm{-}{\mathrm{[}{\mathit{C}}_{\mathit{ij}}^{-1}\mathit{T̅}\mathit{k}\mathit{T̅}\mathit{l}\mathrm{+}{\mathit{C}}_{\mathit{kl}}^{-1}\mathit{T̅}\mathit{i}\mathit{T̅}\mathit{j}\mathrm{-}\mathit{C̅}\begin{array}{c}-1\\ \mathit{ij}\end{array}\mathit{C̅}{\begin{array}{c}-1\\ \mathit{kl}\end{array}}^{\mathrm{]}}}^{\mathrm{)}}\\ & & \mathrm{=}\sum _{\mathit{M}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{M}}\mathrm{[}{\mathit{X}}_{\mathit{LM}}^{\mathit{ij}}{\mathrm{(}\mathit{T̅}\mathit{i}\mathit{T̅}\mathit{j}\mathrm{-}{{\mathit{C}}_{\mathit{ij}}^{-1}}^{\mathrm{)}}}^{\mathrm{]}}\mathrm{[}{\mathit{X}}_{\mathit{L}\mathrm{-}\mathit{M}}^{\mathit{kl}}{\mathrm{(}\mathit{T̅}\mathit{k}\mathit{T̅}\mathit{k}\mathrm{-}{{\mathit{C}}_{\mathit{kl}}^{-1}}^{\mathrm{)}}}^{\mathrm{]}}\mathit{,}\end{array}$(D.3)where the terms in square brackets are exactly of the form of mean-field-subtracted, unnormalized estimates of φ_LM and φ_{L − M}. The remaining terms are equivalent to the N⁽⁰⁾ subtraction in Eq. (17). Any quadratic estimate of the temperature power spectrum will involve a linear combination of terms of the form ${\mathit{T}}_{\mathit{p}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathit{T}}_{\mathit{q}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$, where ${\mathit{T}}_{\mathit{p}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}$ and ${\mathit{T}}_{\mathit{q}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ may be different noisy estimates of the temperature anisotropies but are correlated with themselves and T_i through (at least) the common temperature anisotropies. If we form the covariance between ${\mathit{T}}_{\mathit{p}}^{\mathrm{\left(}\mathrm{1}\mathrm{\right)}}{\mathit{T}}_{\mathit{q}}^{\mathrm{\left(}\mathrm{2}\mathrm{\right)}}$ and the data combination in Eq. (D.1) that enters the ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ estimate, the disconnected terms give the contribution to the covariance between the estimated lensing and temperature power spectra due to sample variance of the temperature anisotropies. However, it is straightforward to show that by including the additional terms in Eq. (D.1) that are quadratic in the disconnected contribution vanishes. We see that, as well as mitigating against errors in the CMB and noise covariance matrices (see Sect. 2), and reducing the covariance between ${\mathit{C}}_{\mathit{L}}^{\mathit{\phi \phi }}$ estimates (Hanson et al. 2011), the data-dependent N⁽⁰⁾ subtraction that we adopt removes much of the covariance with the temperature power spectrum.

The final effect that correlates the lensing and temperature power spectrum estimates is the correlation between T and φ due to the late-time ISW effect in the former (Sect. 6.2). This produces only local correlations between the spectra and they fall rapidly with L being less than 0.5% at L = 40. Their impact on the errors on the lensing amplitude A_L in a joint estimate from the temperature and lensing likelihoods is negligible.

All Tables

All Figures