<figure><img src="image_1.png" /></figure>

# A Measurement of the Cosmic Microwave Background Lensing Potential and Power Spectrum from 500 deg² of SPTpol Temperature and Polarization Data

W. L. K. Wu¹, L. M. Mocanu¹,²,³, P. A. R. Ade⁴, A. J. Anderson⁵, J. E. Austermann⁶, J. S. Avva⁷, J. A. Beall⁶, A. N. Bender¹,⁸, B. A. Benson¹,²,⁵, F. Bianchini⁹, L. E. Bleem¹,⁸, J. E. Carlstrom¹,²,⁸,¹⁰,¹¹, C. L. Chang¹,²,⁸, H. C. Chiang¹²,¹³, R. Citron¹, C. Corbett Moran¹⁴, T. M. Crawford¹,², A. T. Crites¹,²,¹⁵, T. de Haan⁷,¹⁶, M. A. Dobbs¹²,¹⁷, W. Everett¹⁸, J. Gallicchio¹,¹⁹, E. M. George⁷,²⁰, A. Gilbert¹², N. Gupta⁹, N. W. Halverson¹⁸,²¹, N. Harrington⁷, J. W. Henning¹,⁸, G. C. Hilton⁶, G. P. Holder¹⁷,²²,²³, W. L. Holzapfel⁷, Z. Hou¹, J. D. Hrubes²⁴, N. Huang⁷, J. Hubmayr⁶, K. D. Irwin²⁵,²⁶, L. Knox²⁷, A. T. Lee⁷,¹⁶, D. Li⁶,²⁵, A. Lowitz¹, A. Manzotti¹,²⁸, J. J. McMahon²⁹, S. S. Meyer¹,²,¹⁰,¹¹, M. Millea⁷, J. Montgomery¹²,¹⁷, A. Nadolski²²,²³, T. Natoli¹,², J. P. Nibarger⁶, G. I. Noble¹², V. Novosad³⁰, Y. Omori²⁶,³¹, S. Padin¹,²,¹⁵, S. Patil⁹, C. Pryke³², C. L. Reichardt⁹, J. E. Ruhl³³, B. R. Saliwanchik¹³,³⁴, J. T. Sayre¹⁸,²¹, K. K. Schaffer¹,¹¹,³⁵, C. Sievers¹, G. Simard¹², G. Smecher¹²,³⁶, A. A. Stark³⁷, K. T. Story²⁶,³¹, C. Tucker⁴, K. Vanderlinde³⁸,³⁹, T. Veach⁴⁰, J. D. Vieira²²,²³, G. Wang⁸, N. Whitehorn⁴¹, and V. Yefremenko⁸

¹ Kavli Institute for Cosmological Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA; wlwu@kicp.uchicago.edu

² Department of Astronomy and Astrophysics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA

³ Institute of Theoretical Astrophysics, University of Oslo, P.O.Box 1029 Blindern, NO-0315 Oslo, Norway

⁴ Cardiff University, Cardiff CF10 3XQ, UK

⁵ Fermi National Accelerator Laboratory, MS209, P.O. Box 500, Batavia, IL 60510, USA

⁶ NIST Quantum Devices Group, 325 Broadway Mailcode 817.03, Boulder, CO 80305, USA

⁷ Department of Physics, University of California, Berkeley, CA 94720, USA

⁸ High Energy Physics Division, Argonne National Laboratory, 9700 S. Cass Avenue, Argonne, IL 60439, USA

⁹ School of Physics, University of Melbourne, Parkville, VIC 3010, Australia

¹⁰ Department of Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA

¹¹ Enrico Fermi Institute, University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA

¹² Department of Physics, McGill University, 3600 Rue University, Montreal, Quebec H3A 2T8, Canada

¹³ School of Mathematics, Statistics & Computer Science, University of KwaZulu-Natal, Durban, South Africa

¹⁴ TAPIR, Walter Burke Institute for Theoretical Physics, California Institute of Technology, 1200 E California Blvd, Pasadena, CA 91125, USA

¹⁵ California Institute of Technology, MS 249-17, 1216 E. California Blvd., Pasadena, CA 91125, USA

¹⁶ Physics Division, Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA

¹⁷ Canadian Institute for Advanced Research, CIFAR Program in Cosmology and Gravity, Toronto, ON, M5G 1Z8, Canada

¹⁸ Department of Astrophysical and Planetary Sciences, University of Colorado, Boulder, CO 80309, USA

¹⁹ Harvey Mudd College, 301 Platt Blvd., Claremont, CA 91711, USA

²⁰ European Southern Observatory, Karl-Schwarzschild-Str. 2, D-85748 Garching bei München, Germany

²¹ Department of Physics, University of Colorado, Boulder, CO 80309, USA

²² Astronomy Department, University of Illinois at Urbana-Champaign, 1002 W. Green Street, Urbana, IL 61801, USA

²³ Department of Physics, University of Illinois at Urbana-Champaign, 1110 W. Green Street, Urbana, IL 61801, USA

²⁴ University of Chicago, 5640 South Ellis Avenue, Chicago, IL 60637, USA

²⁵ SLAC National Accelerator Laboratory, 2575 Sand Hill Road, Menlo Park, CA 94025, USA

²⁶ Dept. of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, CA 94305, USA

²⁷ Department of Physics, University of California, One Shields Avenue, Davis, CA 95616, USA

²⁸ Institut d'Astrophysique de Paris, 98 bis boulevard Arago, F-75014 Paris, France

²⁹ Department of Physics, University of Michigan, 450 Church Street, Ann Arbor, MI 48109, USA

³⁰ Materials Sciences Division, Argonne National Laboratory, 9700 S. Cass Avenue, Argonne, IL 60439, USA

³¹ Kavli Institute for Particle Astrophysics and Cosmology, Stanford University, 452 Lomita Mall, Stanford, CA 94305, USA

³² School of Physics and Astronomy, University of Minnesota, 116 Church Street S.E. Minneapolis, MN 55455, USA

³³ Physics Department, Center for Education and Research in Cosmology and Astrophysics, Case Western Reserve University, Cleveland, OH 44106, USA

³⁴ Department of Physics, Yale University, 217 Prospect Street, New Haven, CT 06511, USA

³⁵ Liberal Arts Department, School of the Art Institute of Chicago, 112 S Michigan Ave, Chicago, IL 60603, USA

³⁶ Three-Speed Logic, Inc., Vancouver, B.C., V6A 2J8, Canada

³⁷ Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA

³⁸ Dunlap Institute for Astronomy & Astrophysics, University of Toronto, 50 St George St, Toronto, ON, M5S 3H4, Canada

³⁹ Department of Astronomy & Astrophysics, University of Toronto, 50 St George St, Toronto, ON, M5S 3H4, Canada

⁴⁰ Department of Astronomy, University of Maryland College Park, MD 20742, USA

⁴¹ Department of Physics and Astronomy, University of California, Los Angeles, CA 90095, USA

Received 2019 May 14; revised 2019 July 31; accepted 2019 September 3; published 2019 October 14

## Abstract

We present a measurement of the cosmic microwave background lensing potential using 500 deg² of 150 GHz data from the SPTpol receiver on the South Pole Telescope. The lensing potential is reconstructed with signal-to-noise per mode greater than unity at lensing multipoles $L \lesssim 250$, using a quadratic estimator on a combination of cosmic microwave background temperature and polarization maps. We report measurements of the lensing potential power spectrum in the multipole range of $100 < L < 2000$ from sets of temperature-only ($T$), polarization-only (POL), and minimum-variance (MV) estimators. We measure the lensing amplitude by taking the ratio of the measured spectrum to the expected spectrum from the best-fit $\Lambda$ cold dark matter model to the Planck 2015 TT + low $P$ + lensing data set. For the minimum-variance estimator, we find $A_{MV} = 0.944 \pm 0.058(\text{Stat.}) \pm 0.025(\text{Sys.})$; restricting to only polarization data, we find $A_{POL} = 0.906 \pm 0.090(\text{Stat.}) \pm 0.040(\text{Sys.})$. Considering statistical
uncertainties alone, this is the most precise polarization-only lensing amplitude constraint to date (10.1σ) and is
more precise than our temperature-only constraint. We perform null tests and consistency checks and find no
evidence for significant contamination.

***Key words:*** cosmology: cosmic background radiation – large-scale structure of the universe

**1. Introduction**

Gravitational potentials associated with large-scale structure deflect the paths of cosmic microwave background (CMB) photons as they propagate from the surface of last scattering—a process called gravitational lensing (Blanchard & Schneider 1987). Gravitational lensing breaks the statistical isotropy of the CMB and introduces correlations across CMB temperature and polarization fluctuations on different angular scales. These correlations are proportional to the projected gravitational potentials integrated along the line of sight and therefore can be used to reconstruct the lensing potential (Lewis & Challinor 2006). The lensing potential is a probe of the growth of large-scale structure and the geometry of the universe between the epoch of recombination and today. Thus, from CMB observations alone, we can extract information about the universe at both the redshift of last scattering ($z \simeq 1100$) and the redshifts of structure formation ($z \lesssim 3$) and dark energy domination. This makes CMB lensing a powerful tool for pursuing some of the most ambitious goals in cosmology and particle physics today (e.g., CMB-S4 Collaboration et al. 2016), including constraining the sum of neutrino masses and the amplitude of primordial gravitational waves (Lesgourgues & Pastor 2006; Kamionkowski & Kovetz 2016).

While using the CMB lensing measurements from this work, we will not detect the sum of neutrino masses or significantly constrain primordial gravitational waves through delensing (Manzotti et al. 2017), these measurements will nevertheless provide relevant constraints to parameters of the standard cosmological model lambda cold dark matter (ΛCDM). This is of particular interest currently because some optical probes of gravitational lensing are in mild tension with Planck's CMB + lensing constraints on matter density and fluctuations (Abbott et al. 2018; Hildebrandt et al. 2018; Hikage et al. 2019). The optical lensing measurements use the effect of gravitational lensing on the apparent shapes of background galaxies to measure the intervening gravitational potentials. Compared with using galaxies, an attractive characteristic of CMB lensing is that the source plane has nearly Gaussian statistics with a well-characterized angular power spectrum and is at a high and well-known redshift $z = 1089.8 \pm 0.2$ (Planck Collaboration et al. 2016b, 2018a). Owing to the high-redshift source plane, CMB lensing probes the integrated matter fluctuations to redshifts beyond optical surveys. Furthermore, CMB lensing measurements have different instrumental and astrophysical systematics compared with those from optical surveys. The use of independent probes will therefore help us investigate the source of the tension. It is thus an important goal for us to identify potential sources of instrumental and/or astrophysical systematics as the field advances the precision of CMB lensing measurements.

In the last few years, CMB lensing has entered the era of precision measurements. This lensing effect was first detected through cross correlations with radio sources (Smith et al. 2007). Since then, the CMB lensing potential power spectrum has been measured by multiple groups using temperature data only (*T*, Das et al. 2011; van Engelen et al. 2012; Planck

Collaboration et al. 2014; Omori et al. 2017), polarization data only (POL, POLARBEAR Collaboration et al. 2014; BICEP2 Collaboration et al. 2016), and combinations of temperature and polarization data (Story et al. 2015; Sherwin et al. 2017; Planck Collaboration et al. 2018b). The most precise lensing amplitude measurement, at $40\sigma$, comes from Planck's minimum-variance (MV) estimator that combines both temperature and polarization estimators; in that measurement, the temperature reconstruction contributes most of the signal-to-noise ratio (S/N). More generally, in prior lensing measurements that used both temperature and polarization maps, the *T* estimator has always dominated the overall measurement precision. In this work, for the first time, the polarization-only (POL) measurement is more constraining than the *T* measurement. Furthermore, if we consider only the statistical uncertainty, we have a 10% constraint (10.1σ) on the lensing amplitude using polarization data alone—the tightest constraint of its kind to date.

In this work, we extend the lensing measurement to a 500 deg² field from the 100 deg² field of Story et al. (2015, hereafter S15). We build on the lensing pipeline presented in S15 with two main modifications. First, instead of including the Monte Carlo (MC) bias—the difference between the recovered lensing spectrum from simulation and the input spectrum—as a systematic uncertainty, we identify its main contributor and correct for the bias using a multiplicative correction factor. Second, instead of treating extragalactic foreground biases as negligible, we subtract an expected foreground bias from the *T* and MV lensing spectra.

Because the input CMB maps are of similar depths as those used in S15, we also measure lensing modes with S/N better than unity for $L \lesssim 250$. However, in this analysis, we have ~5 times more sky area and therefore are able to make a more precise measurement of the lensing spectrum. With both the *T* and POL lensing amplitudes well constrained, we improve the precision of the MV lensing amplitude measurement from S15's 14% to 6%. This is approaching the ~3% precision of the lensing amplitude measurement from Planck (Planck Collaboration et al. 2018b). These two measurements arrive at their respective precisions from very different regimes: while the Planck measurement covers 67% of the sky, each lensing mode is measured with low S/N; our measurement covers only 1% of the sky, but many lensing modes are measured with high S/N. When our measurements are projected to cosmological parameter space, the constraint in the $\sigma_8-\Omega_M$ plane is only slightly weaker than those from Planck (Planck Collaboration et al. 2018b). This will be useful in illuminating the aforementioned mild tension between optical probes and Planck. In a future paper, we will present cosmological parameter constraints and comparisons with other lensing probes.

This paper is organized as follows: in Section 2, we describe the data set used and the simulated skies generated for this analysis. In Section 3, we summarize the lensing analysis pipeline and describe new aspects. We present the lensing measurements in Section 4 and show that our measurements are robust against systematics in Section 5.1. In Section 5.2, we
account for systematic uncertainties from sources that can bias
the lensing measurement. We conclude in Section 6.

## 2. Data and Simulations

In this section, we describe the SPTpol survey, the data
processing steps taken to generate the data maps, and the
simulated skies generated for this analysis.

### 2.1. SPTpol 500 deg² Survey

The South Pole Telescope (SPT, Padin et al. 2008; Carlstrom et al. 2011) is a 10 m diameter off-axis Gregorian telescope located at the Amundsen-Scott South Pole Station in Antarctica. In this work, we use data from the 150 GHz detectors from the first polarization-sensitive receiver on SPT, SPTpol. The SPTpol 500 deg² survey spans 4 hours in right ascension (R.A.), from 22h to 2h, and 15° in declination, from -65° to -50°. We use observations conducted between 2013 April 30 and 2015 October 27 (after the 100 deg² survey for S15 was finished), which include 3491 independent maps of the 500 deg² survey footprint. The field was observed using two strategies. Initially, we used a “lead-trail” scanning strategy, where the field was divided into two halves in R.A. The telescope first scanned the lead half and then switched to the trail half such that they were both observed over the same range in azimuth, and therefore the same patch of ground. In 2014 May, the scanning strategy was switched to a “full-field” strategy, in which the full R.A. range of the field was covered in a single observation.

### 2.2. Data Processing

The data reduction for this set of maps follows that applied
to the TE/EE power spectrum analysis of the same
field (Henning et al. 2018, hereafter H18). We therefore
highlight here only aspects that are different or are particularly
relevant for this analysis.

An observation is built from a collection of constant-elevation
scans, where the telescope moves from one end of the R.A.
range to the other. For every scan, we filter the time streams,
which corresponds in map space to mode removal along the scan
direction. We subtract from each detector's time stream a
Legendre polynomial as an effective high-pass filter. We choose
a third-order polynomial for the lead-trail scans and a fifth-order
polynomial for the full-field scans. We combine this with a high-
pass filter with a cutoff frequency corresponding to angular
multipole $l = 100$ in the scan direction. During the polynomial
and high-pass filtering step, regions within 5' of point sources
brighter than 50 mJy at 150 GHz are masked in the time streams
to avoid ringing artifacts. Finally, we low-pass the time streams
at an effective multipole of $l = 7500$ in the scan direction to
avoid high-frequency noise aliasing into the signal due to the
adopted pixel resolution of 1' in the maps.

Electrical cross talk between detectors can bias our measure-
ment. In S15, we accounted for this bias as a systematic
uncertainty of 5% on the MV lensing amplitude. In this analysis,
we correct the cross talk between detectors at the time stream
level as described in H18. With this correction, cross talk is
suppressed by more than an order of magnitude. With the 5%
uncertainty from S15 being an upper limit before this
suppression, we conclude that the residual cross talk introduces
negligible bias to our lensing amplitude measurement.

Before binning into maps, we calibrate the time streams
relative to each other using observations of the HII region

RCW38 and an internal chopped thermal source. The per-
detector polarization angles are calibrated based on measure-
ments from a polarized thermal source (Crites et al. 2015).
After that, the time stream data are binned into T/Q/U maps
with 1' × 1' pixels using the oblique Lambert azimuthal equal-
area projection.

We apply the following corrections to the coadded map of the
individual observations to obtain the final map: $T \rightarrow P$
monopole subtraction, global polarization rotation, absolute
calibration, and source and boundary masking. We measure
$T \rightarrow P$ leakage by taking a weighted average over multipole
space of the cross correlation of the temperature map with either
the Q or the U map. We find leakage factors of $\epsilon^Q = 0.018$
and $\epsilon^U = 0.008$. We obtain their uncertainties from the spread of
leakage factors derived from 100 different half-splits of the data
and find the fractional uncertainties to be 0.1% and 0.3%,
respectively. We deproject a monopole leakage term from the
maps by subtracting a copy of the temperature map scaled by
these factors from the Q and U maps.

Assuming ΛCDM cosmology, we expect the cross spectra
TB and EB to be consistent with zero. Therefore, we can
estimate and apply the global polarization angle rotation
needed to minimize the TB and EB correlations (Keating et al.
2013). We find the rotation angle to be 0.63 ± 0.04 and rotate
the Q and U maps accordingly.

The final absolute calibration is obtained by comparing the
final coadded SPTpol *T* map with *Planck* over the angular
multipole range 600 < *l* < 1000. Specifically, we take the
ratio of the cross spectrum of two half-depth SPTpol maps to
the cross spectrum of full-depth SPTpol and *Planck* maps and
require the ratio to be consistent with 1 to determine the
calibration factor. We estimate the polarization efficiency (or
polarization calibration factor) similarly, by comparing a full-
depth SPTpolE-mode map with the *Planck* E-mode map. The
temperature and polarization calibration factors ($T_{cal}$ and $P_{cal}$) as
derived are 0.9088 and 1.06, respectively (see H18 for details).
We apply $T_{cal}$ to the temperature map and $T_{cal} \times P_{cal}$ to the
polarization maps to obtain calibrated temperature and
polarization maps. In the legacy *Planck* release, their polariza-
tion efficiency estimates are found to be potentially biased at
the 1%–2% level (see Table 9 of Planck Collaboration et al.
2018c). To circumvent potential biases when we tie our
polarization map calibration to *Planck's* E-mode map, we
adjust $P_{cal}$ by a multiplicative factor obtained from H18 without
using *Planck's* polarization maps. In H18, $T_{cal}$ and $P_{cal}$ were free
parameters in the fits of the PlanckTT+SPTpolEETE data
set to the ΛCDM + foregrounds + nuisance parameters model,
with priors set by the $T_{cal}$ and $P_{cal}$ values and uncertainties
derived from the comparison against *Planck* described above.
In this work, we multiply $P_{cal}$ by the best-fit $P_{cal}$ parameter
(1.01) from H18 (see Table 5 of H18). In Section 5.2, we use
the posterior uncertainties of the $T_{cal}$ and $P_{cal}$ parameters to
quantify their contributions to the systematic uncertainties of
the lensing amplitude measurements.

We apply a mask that defines the boundary of the T/Q/U
maps to downweight the noisy field edges. In addition, we
mask bright point sources in the map. We use a 5' radius to
mask point sources from the SPT-SZ catalog by W. Everett
et al. (in preparation) with flux density above 6 mJy at either 95
or 150 GHz that are in the 500 deg² footprint. We use a 10'
radius for point sources with flux density greater than 90 mJy at
either 95 or 150 GHz. Clusters detected in Bleem et al. (2015)
within the 500 deg² survey footprint are masked using a 10' radius. The mask has a top hat profile for both the boundary and the sources, and mode mixing due to the mask is suppressed by the inverse-variance map-filtering step (Section 3.2). We compare this lensing measurement with one using a cosine-tapered apodization on the mask edges in Section 5.2 and find them to be consistent.

The noise levels in the coadded maps are 11.8 µK–arcmin in temperature and 8.3 µK–arcmin in polarization over the multipole range $1000 < l < 3000$.<fnref n="42" /> The map depth is similar to that of the 100 deg² field used in S15.<fnref n="43" /> Because this map covers five times the sky area, the sample variance of the lensing power spectrum is reduced and hence the lensing spectrum measurement precision is improved.

## 2.3. Simulations

We use simulated skies to estimate the mean-field (MF) bias $\bar{\phi}^{\text{MF}}$ to the lensing potential map (Section 3.1), to correct the analytical response of the lensing estimator (Section 3.1), to correct the expected biases ($N_L^{(0),\text{RD}}$ and $N_L^{(1)}$, see Section 3.3) in the raw lensing power spectrum, and to estimate the uncertainties of the lensing measurements (Section 4).

The simulated skies contain the CMB, foregrounds, and instrumental noise. The input cosmology for generating the CMB is the best-fit ΛCDM model to the Planck plikHM_TT_lowTEB_lensing data set (second column in Table 4 of Planck Collaboration et al. 2016a). Using the best-fit cosmology, we use CAMB<fnref n="44" /> to generate theory spectra, which are then fed into HEALPIX (Górski et al. 2005) to generate the spherical harmonic coefficients $a_{\ell m}$ for T, E, B, and the lensing potential $\phi$. The CMB $a_{\ell m}$ are projected to maps and lensed by $\phi$ using LENSPIX (Lewis 2005). The lensed maps are then converted back to $a_{\ell m}$, at which point foreground fluctuations are added, and the $a_{\ell m}$ are multiplied by the $l$-space instrument beam window function. The $a_{\ell m}$ are subsequently projected on an equidistant cylindrical projection grid for mock observation, which produces mock skies that are processed identically as the real sky.

The foreground components are generated as Gaussian realizations of model angular power spectra. We include power from thermal and kinematic Sunyaev–Zel’dovich effects (tSZ and kSZ), the cosmic infrared background (CIB), radio sources, and galactic dust. The amplitudes for tSZ, kSZ, and CIB components are taken from George et al. (2015), which has the same source masking threshold as this work. We use an amplitude $D_{\ell=3000}^{\text{tSZ}+\text{kSZ}} = 5.66 \mu\text{K}^2$ with a model shape from Shaw et al. (2010) for the sum of tSZ and kSZ components.<fnref n="45" /> We use $D_{\ell=3000}^{\text{radio}} = 1.06 \mu\text{K}^2$ with $D_{\ell} \propto l^2$ for the radio source component. We use $D_{\ell=3000}^{\text{CIB},P} = 9.16 \mu\text{K}^2$ with $D_{\ell} \propto l^2$ for the

unclustered CIB component and $D_{\ell=3000}^{\text{CIB,cl}} = 3.46 \mu\text{K}^2$ with $D_{\ell} \propto l^{0.8}$ for the clustered CIB. By modeling these terms as Gaussian, we have neglected the potential bias their non-Gaussianities can introduce to the lensing spectrum. We treat this potential bias explicitly later in the analysis (see Section 3.3). We assume a 2% polarization fraction for all the Poisson-distributed (unclustered) components to model extragalactic polarized emission (Seiffert et al. 2007). We model galactic dust power in temperature and polarization as a power law with $D_{\ell} \propto l^{-0.42}$, with $D_{\ell=80}^{\text{TT,dust}} = 1.15 \mu\text{K}^2$, $D_{\ell=80}^{\text{EE,dust}} = 0.0236 \mu\text{K}^2$, and $D_{\ell=80}^{\text{BB,dust}} = 0.0118 \mu\text{K}^2$ (Keisler et al. 2015).

The instrumental noise realizations are generated by subtracting half of all the observations from the other half. By randomly grouping the observations into one of the halves, we create 500 different realizations of noise from the data themselves. The noise realizations are added to the simulated skies after mock observation.

We generate 500 lensed skies including CMB, foregrounds, and noise. All 500 skies are used for estimating the lensing spectrum bias $N_L^{(0),\text{RD}}$. We use 100 skies to estimate the MF ($\bar{\phi}^{\text{MF}}$, see Section 3.1): 50 for each of the two lensing potential $\phi$ estimates that enter the lensing spectrum calculation ($\hat{\phi}_L^{\text{UV}}$ and $\hat{\phi}_L^{\text{XY}}$ in Equation (6)). We use the other 400 to correct the analytical response of the lensing estimator and to obtain the statistical uncertainty of the lensing spectrum. In addition, we generate 500 unlensed skies with all other inputs being identical. We use this set of simulations for testing the lensing pipeline and assessing the probability of detecting lensing from unlensed skies.

The $N_L^{(1)}$ bias is estimated using a different set of 200 lensed skies. We start with 200 realizations of unlensed skies and divide them into two groups. We then lens one sky from each group by the same $\phi$ (Section 3.3). We do not add instrumental noise or foregrounds to this set of simulations, because $N_L^{(1)}$ comes from correlations between the CMB and the lensing modes. At the map-filtering step, we filter this set of simulations assuming the same level of foregrounds and noise as the other set of lensed skies.

## 3. Lensing Analysis

In this section, we describe the lensing pipeline that produces an unbiased estimate of the lensing potential and of the lensing power spectrum. This analysis follows that of S15 except for the treatment of the MC bias and foreground subtraction in the lensing spectrum. We summarize the steps here and refer the reader to S15 for details.

### 3.1. Estimating the Lensing Potential

The unlensed CMB sky is well approximated by a statistically isotropic, Gaussian random field. Gravitational lensing breaks the statistical isotropy of the fluctuations and introduces correlations between the otherwise uncorrelated Fourier modes of the CMB temperature and polarization maps. We use these correlations to estimate the underlying lensing potential with pairs of filtered maps using the quadratic estimator of Hu & Okamoto (2002):

$$ \bar{\phi}_L^{XY} = \int d^2\ell W_{\ell,\ell-L}^{XY} \bar{X}_\ell \bar{Y}_*_{\ell-L}, \quad (1) $$

<footnote n="42">⁴² The temperature map noise is higher than the polarization map noise in this multipole range because of atmospheric noise contributions.</footnote>

<footnote n="43">⁴³ The noise levels provided in S15 were estimated without the $T_{\text{cal}}$ and $P_{\text{cal}}$ corrections and are thus higher.</footnote>

<footnote n="44">⁴⁴ http://camb.info</footnote>

<footnote n="45">⁴⁵ Given that we mask all clusters in the Bleem et al. (2015) catalog, this tSZ level is high by ~2 $\mu\text{K}^2$. As a result, we nonoptimally downweight high-$l$ modes in the temperature data map and slightly overestimate the $N_L^{(0)}$ noise (Section 3.3). We compare the analytic $N_L^{(0)}$ between $D_{\ell=3000}^{\text{tSZ}+\text{kSZ}} = 5.66 \mu\text{K}^2$ and $D_{\ell=3000}^{\text{tSZ}+\text{kSZ}} = 3.61 \mu\text{K}^2$ (keeping all other inputs equal) and find that the latter is ~2% lower. Because there is both signal and noise variance in the uncertainty of the lensing amplitude measurement, and the T lensing amplitude uncertainty is 12% (see Section 4), the T lensing amplitude would at most be reduced by ~2% of 12%, which is negligible.</footnote>
where $\bar{X}$ and $\bar{Y}$ are filtered $T$, $E$, or $B$ fields as outlined in Section 3.2, $W$ is a weight function unique to each input pair of $XY$ maps, $\ell$ are modes of the CMB, and $L$ are modes of the lensing potential. We form lensing potentials from $XY \in [TT, TE, EE, EB, TB]$ and modify $W$ from Hu & Okamoto (2002) by replacing the un lensed CMB spectra with lensed spectra to reduce a higher-order bias ($N_L^{(2)}$, Hanson et al. 2011). As written, $\bar{\phi}$ is a biased estimate of $\phi$. To arrive at an unbiased estimate of $\phi$, we remove an additive bias (the MF) and normalize the estimator by the response (defined below).

The MF can arise from masking and inhomogeneous noise— sources that introduce mode coupling across angular scales. These mode couplings persist even when the CMB and $\phi$ realizations are different. We can therefore estimate the MF by averaging $\bar{\phi}$ from many realizations of input lensed CMB maps. For the $L$ range considered in this work, the MF is subdominant compared with the lensing signal spectrum. Specifically, it is $\lesssim 30\%$ of the lensing power spectrum in the first $L$ bin. It grows larger toward larger angular scales, dominated by the effect of the mask.

We construct the response using an analytical calculation corrected by simulations. The analytic response of the estimator $R_{|L|}^{XY,Analytic}$ is

$$R_{|L|}^{XY,Analytic} = \int d^2\ell W_{\ell,\ell-L}^{XY} \times W_{\ell,\ell-L}^{XY} \mathcal{F}_{\ell}^{X} \mathcal{F}_{\ell-L}^{Y}, \quad (2)$$

where $\mathcal{F}_{\ell}^{X} = (C_{\ell}^{XX} + N_{\ell}^{XX})^{-1}$ describes the diagonal approximation to the filter applied to the input CMB maps (Section 3.2). Because this filter assumes spatial stationarity of the statistics of the signal and the noise, it does not account for nonstationary effects from, for example, the boundary and source mask and causes the response to be slightly misestimated. To account for this, we use simulations to estimate the correction to the analytic response. The total response is thus $R_{|L|}^{XY} = R_{|L|}^{XY,Analytic} \times R_{L}^{XY,MC}$, with $R_{L}^{XY,MC}$ denoting the correction estimated from simulations. We extract $R_{L}^{XY,MC}$ by first taking the cross spectrum between the input $\phi^{in}$ and the intermediate $\phi$ estimate $\hat{\phi}^{\prime XY}$, which has been MF subtracted and normalized by the analytic response: $\hat{\phi}_{L}^{\prime XY} = (\bar{\phi}_{L}^{XY} - \bar{\phi}_{L}^{XY,MF})/R_{|L|}^{XY,Analytic}$. We then take the ratio of the average of the cross spectra $\langle \phi_{L}^{in*} \hat{\phi}_{L}^{\prime XY} \rangle$ averaged within each $L$ annulus with the input spectrum over 400 simulation realizations. Similar to S15, we find the $R_{L}^{XY,MC}$ to be a $\le 10\%$ correction.

The normalized and MF-corrected $\hat{\phi}$ is

$$\hat{\phi}_{L}^{XY} = \frac{1}{R_{|L|}^{XY}} (\bar{\phi}_{L}^{XY} - \bar{\phi}_{L}^{XY,MF}) \quad (3)$$

for the individual estimators where $XY \in [TT, TE, EE, EB, TB]$. In this work, we use $T$ to denote the lensing potential and spectrum constructed using only temperature data, that is, $\hat{\phi}_{L}^{TT}$ and its spectrum. We use POL to denote the potential and spectrum constructed using only polarization data, and MV to denote the potential and spectrum constructed using both temperature and polarization data.

To construct the combined MV and POL estimators, we first form a weighted average of the input estimators to form the intermediate $\hat{\phi}_{L}^{\prime}$. We use the inverse-noise variance of the input

estimators as weights, which are approximated by their analytical responses. We then extract $R_{L}^{MC}$ by forming cross spectra between $\hat{\phi}_{L}^{\prime MV}$ or $\hat{\phi}_{L}^{\prime POL}$ with the input $\phi_{L}^{in}$. Putting it all together, the unbiased MV and POL lensing potentials are constructed as

$$\hat{\phi}_{L} = \frac{1}{R_{L}^{MC}} \frac{\sum_{XY} \bar{\phi}_{L}^{XY} - \bar{\phi}_{L}^{XY,MF}}{\sum_{XY} R_{|L|}^{XY,Analytic}}, \quad (4)$$

where $XY \in [TT, TE, EE, EB, TB]$ for the MV estimator and $XY \in [EE, EB]$ for the POL estimator.

## 3.2. Input CMB Map Filtering

We filter our input maps with an inverse-variance ($C^{-1}$) filter, which is derived such that the variance of a lensing field reconstructed by a quadratic estimator is minimized. The filter is constructed identically as in S15: we assume the data maps to be composed of three components: sky signal, “sky noise,” and pixel domain noise. The sky signal and noise are modeled in the Fourier domain and they include the CMB, astrophysical foregrounds, and atmospheric noise. The pixel domain noise is modeled as white, uncorrelated, and spatially nonvarying inside the mask. Concretely, we solve for the inverse-variance filtered Fourier modes $\bar{X}$ in the following expression using conjugate-gradient-descent:

$$[S^{-1} + P^{\dagger}n^{-1}P]S\bar{X} = P^{\dagger}n^{-1}d. \quad (5)$$

Here $S \equiv C_{\ell}^{XX} + N_{\ell}^{XX}$ describes the sky signal and noise components, $n^{-1}$ is the inverse of the map noise variance and is zero for masked pixels, $P$ applies the filter transfer function, and inverse-Fourier transforms the map from pixel space to Fourier space. $P$ additionally transforms $Q/U$ to $E/B$ for polarization maps, and $d$ denotes the input pixel-space $T/Q/U$ maps.

This filter approaches the simple form $\bar{X} = S^{-1}P^{\dagger}d = (C_{\ell}^{XX} + N_{\ell}^{XX})^{-1}\bar{X}$ in regions far away from the mask boundary. To see that, one can absorb the pixel-domain noise into the $N_{\ell}^{XX}$ term (which is valid under our model of the data for pixels far from the mask boundary); in the limit where $n \to 0$, the diagonal form is exact.

The inputs to this filter are as follows. For the sky signal component of $S$, we use the same lensed CMB and astrophysical foreground spectra used to generate simulations (Section 2.3). For the sky noise component, we take the averaged power spectra from the temperature and polarization noise realizations and from them subtract noise floors of 7 µK–arcmin to form $N_{\ell}^{XX}$. We set 7 µK–arcmin as the white noise level of the pixel-domain noise of the $T, Q$, and $U$ maps. The filter transfer function $P$ includes time stream filtering, beam, and pixel-window function. We approximate the time stream filtering transfer function using simplified simulations that capture the lost modes along the $\ell_x = 0$ axis due to our scan strategy. We measure the SPT beam using the response of the detectors to Venus, as was done in H18. In addition, we fit the measured beam with a function of form $B(\ell) = \sum_{i=1,2} A_i \exp(-0.5(\ell/\ell_i)^{p_i})$ to obtain a smooth profile at $\ell < 400$. We incorporate the pixel-window function for the map pixel size of $1' \times 1'$. Our input CMB map $\ell$ range for this analysis is $|\ell_x| > 100$ set by the time stream $\ell$-space high-pass filter and $|\ell| < 3000$ set by concerns of foreground contamination (Section 5.2).
### 3.3. Estimating the Lensing Potential Power Spectrum

Upon obtaining the unbiased estimate of the lensing potential $\hat{\phi}$, we calculate the raw power spectrum of the lensing potential $C_L^{\hat{\phi}^{UV}\hat{\phi}^{XY}}$ by forming cross spectra of $\hat{\phi}_L^{XY}$ and $\hat{\phi}_L^{UV}$ where $UV, XY \in [TT, TE, EE, EB, TB]$:

$$C_L^{\hat{\phi}^{UV}\hat{\phi}^{XY}} \equiv f_{\text{mask}}^{-1} \sum_{|L|=L} \langle \hat{\phi}_L^{UV} \hat{\phi}^{*XY} \rangle, \quad (6)$$

where $f_{\text{mask}}$ is the average value of the fourth power of the mask.

This power spectrum is biased, and we correct four sources of bias in this analysis. Two of them, the $N_L^{(0)}$ and $N_L^{(1)}$ biases, arise from spurious correlations of the input fields. The third source of bias comes from foregrounds, $\Delta C_L^{\phi\phi,FG}$. Finally, we correct for a multiplicative bias due to higher-order coupling of the source mask $f_{\text{PS}}$. Therefore, the unbiased lensing spectrum is estimated as

$$\hat{C}_L^{\hat{\phi}\hat{\phi}} = f_{\text{PS}} [C_L^{\hat{\phi}\hat{\phi}} - N_L^{(0),\text{RD}} - N_L^{(1)}] - \Delta C_L^{\phi\phi,\text{FG}}. \quad (7)$$

Figure 1 shows $N_L^{(0),\text{RD}}$, $N_L^{(1)}$, and $\Delta C_L^{\phi\phi,\text{FG}}$ for the MV reconstruction. $N_L^{(0),\text{RD}}$ and $N_L^{(1)}$ are similar to those in S15 because the CMB map noise levels for these two works are about the same. We describe how each bias term is estimated in the following paragraphs.

The lensing spectrum is a 4-point function, or trispectrum, of the observed fields and thus contains, in addition to the connected term caused by lensing, a disconnected term from correlations of Gaussian fields ($N_L^{(0)}$). Secondary contractions of the trispectrum give rise to connected terms that also bias the lensing spectrum ($N_L^{(1)}$). It is called $N_L^{(1)}$ because it is first order in $C_L^{\phi\phi}$ (Kesden et al. 2003; Hanson et al. 2011). We estimate the power of the $N_L^{(0)}$ and $N_L^{(1)}$ terms using simulations and subtract them from the raw lensing spectrum.

We estimate the disconnected term using a realization-dependent method ($N_L^{(0),\text{RD}}$, Namikawa et al. 2013). The standard $N_L^{(0)}$ is estimated by computing the lensing spectra with two input maps that form $\hat{\phi}$ from different simulation realizations. Because input maps from different realizations do not share the lensing potential, the resultant power in the lensing spectra comes from spurious correlations of the maps. With the realization-dependent method, we use a combination of simulation and data maps. Specifically, on top of the mix of simulation realizations, we also construct lensing potentials using data in one of the input maps. Accounting for all the different combinations, we have

$$\begin{aligned} N_L^{(0),\text{RD}} = & \langle + C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_d, \bar{V}_{\text{MC}}, \bar{X}_d, \bar{Y}_{\text{MC}}] \\ & + C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_{\text{MC}}, \bar{V}_d, \bar{X}_d, \bar{Y}_{\text{MC}}] \\ & + C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_d, \bar{V}_{\text{MC}}, \bar{X}_{\text{MC}}, \bar{Y}_d] \\ & + C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_{\text{MC}}, \bar{V}_d, \bar{X}_{\text{MC}}, \bar{Y}_d] \\ & - C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_{\text{MC}}, \bar{V}_{\text{MC}'}, \bar{X}_{\text{MC}'}, \bar{Y}_{\text{MC}'}] \\ & - C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_{\text{MC}}, \bar{V}_{\text{MC}'}, \bar{X}_{\text{MC}'}, \bar{Y}_{\text{MC}}] \rangle_{\text{MC},\text{MC}'}, \quad (8) \end{aligned}$$

where $C_L^{\hat{\phi}^{UV}\hat{\phi}^{XY}}$ is expressed as $C_L^{\hat{\phi}\hat{\phi}} [\bar{U}, \bar{V}, \bar{X}, \bar{Y}]$ to more clearly indicate the sources of the input maps. Here MC and MC' denote simulation skies from different realizations, and d

<figure><img src="image_2.png" /><figcaption>Figure 1. Bias terms subtracted from the raw MV lensing spectrum: $N_L^{(0),\text{RD}}$, $N_L^{(1)}$ (converted to $C_L^{\kappa\kappa}$ units, Equation (12)), and foreground bias $\Delta C_L^{\kappa\kappa,\text{FG}}$. The theoretical lensing convergence spectrum for the fiducial cosmology is shown as a black line. The disconnected term in the CMB 4-point function $N_L^{(0),\text{RD}}$ is shown as a dashed yellow line (Equation (8)). The spurious correlated power between the CMB and the lensing potential $N_L^{(1)}$ is shown as a dashed pink line (Equation (9)). The sum of these two terms is labeled “total noise bias” in the figure as a solid green line. The foreground bias from the tSZ trispectrum, CIB trispectrum, and the tSZ and CIB correlation with $\kappa$ is shown as a dashed dark blue line (Section 3.3). The total noise bias is also an estimate of the noise in the reconstructed $\phi$ map. Lensing modes with values of L for which the total bias is less than the signal spectrum are measured at S/N greater than unity. For this measurement, that includes all modes with $L \lesssim 250$, or angular scales of roughly a degree and larger.</figcaption></figure>

denotes the data map. The realization-dependent method produces a better estimate of the disconnected term because it reduces the bias from the mismatch between the fiducial cosmology used for generating the simulations and that in the data. In addition, this method reduces the covariance between the lensing potential bandpowers.

The second bias term $N_L^{(1)}$ comes from spurious correlations between the CMB and the lensing potential and is proportional to $C_L^{\phi\phi}$. We estimate it using pairs of simulated skies that have different realizations of unlensed CMB lensed with the same $\phi$:

$$\begin{aligned} N_L^{(1)} = & \langle + C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_{\phi^1,\text{MC}}, \bar{V}_{\phi^1,\text{MC}'}, \bar{X}_{\phi^1,\text{MC}}, \bar{Y}_{\phi^1,\text{MC}'}] \\ & + C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_{\phi^1,\text{MC}}, \bar{V}_{\phi^1,\text{MC}'}, \bar{X}_{\phi^1,\text{MC}'}, \bar{Y}_{\phi^1,\text{MC}'}] \\ & - C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_{\text{MC}}, \bar{V}_{\text{MC}'}, \bar{X}_{\text{MC}'}, \bar{Y}_{\text{MC}'}] \\ & - C_L^{\hat{\phi}\hat{\phi}} [\bar{U}_{\text{MC}}, \bar{V}_{\text{MC}'}, \bar{X}_{\text{MC}'}, \bar{Y}_{\text{MC}}] \rangle_{\text{MC},\text{MC}'}, \quad (9) \end{aligned}$$

where the subscript $\phi^1$ denotes CMB maps lensed by the same $\phi$ realization.

We account for biases to the lensing spectrum due to foregrounds in the temperature map. As studied in van Engelen et al. (2014), both tSZ and CIB have a trispectrum, which leads to a response in the lensing power spectrum. This bias enters the lensing spectrum through the 4-point function of the temperature map, thus modifying the T spectrum and the MV
spectrum. Additionally, Because both tSZ and CIB trace the same large-scale structure as the lensing field, the non-Gaussianities of both fields can mimic lensing and couple through the $\phi$ estimator in a coherent way that correlates with $\phi$, forming a nonzero bispectrum of the matter density field (denoted by $tSZ^{2-\kappa}$ and CIB$^{2-\kappa}$). This effect biases the MV lensing spectrum and spectra from pairs of estimators of the form $\langle\phi^{TT}\phi^{UV}\rangle$, where $UV \in [TT, TE, EE, EB, TB]$. The level of foreground bias is scale dependent. It is negative and close to flat at $\sim 2\%$ for $L < 1300$ for the total bias contributed from tSZ and CIB through their trispectra and correlations with $\phi$. As $L$ increases, the magnitude of this bias decreases and reaches a null at about $L$ of 2000. We subtract the relevant terms from the $T$ and MV spectra by using the bias estimates from van Engelen et al. (2014). Specifically, we subtract the total foreground bias coming from the tSZ trispectra, CIB trispectra, $tSZ^{2-\kappa}$, and CIB$^{2-\kappa}$ from the $T$ lensing spectrum. For the MV spectrum, we subtract the total foreground bias that enters through $\langle\phi^{TT}\phi^{TT}\rangle$ and $tSZ^{2-\kappa}$ and CIB$^{2-\kappa}$ biases that enter through $\langle\phi^{TT}\phi^{UV}\rangle$ with $UV \in [TE, EE, EB, TB]$. We compute the bias fraction by forming a weighted average of the foreground biases to $\langle\phi^{TT}\phi^{TT}\rangle$ and $\langle\phi^{TT}\phi^{UV}\rangle$. We use as weights the fractional contribution of the estimators to the MV estimator. The size of this bias in the MV spectrum is shown in Figure 1.

MC bias describes the difference between the recovered amplitude in simulations and the input spectrum. It is typically found to be small (e.g., Sherwin et al. 2017; Planck Collaboration et al. 2018b). In S15, we found that the mean MV amplitude $A_{MV}$ from simulations was 3% below unity and treated this discrepancy as a source of systematic uncertainty. In this work, we find that the main source of our MC bias comes from the inclusion of the point source and cluster mask. When analyzing the set of simulations with the source mask removed (while keeping the boundary mask), we are able to recover $A_{MV} = 1$ to within $1\sigma$ of the standard error of 400 sky realizations (0.058/√400 = 0.0029) for the multipole range we report in this work. Furthermore, comparing the mean recovered spectrum of the MV, POL, and $T$ estimators with and without the source mask applied to the input maps, we observe a relatively constant multiplicative offset in the range of $100 < L \lesssim 800$. We therefore conclude that the main source of our MC bias is due to higher-order coupling generated from the presence of source masks in the map that is not accounted for by $f_{\text{mask}}$. We construct a multiplicative correction for this MC bias $f_{\text{PS}}$ as an inverse of the mean of the simulation lensing amplitudes $A_{\text{lens}}$ estimated from $100 < L < 602$:

$$f_{\text{PS}}^{XY,UV} = \langle A_{\text{lens}}^{XY,UV} \rangle^{-1}. \qquad (10)$$

To check the stability of this estimate, we vary the maximum $L$ range used between 446 and 813 and find the resultant simulation $A_{\text{lens}}$ from all pairs of estimators to be consistent with unity to within $2\sigma$ of their standard errors. The correction from this step is about 5%. Specifically, $f_{\text{PS}}$ for the MV, POL, and $T$ are $f_{\text{PS}}^{\text{MV}} = 1.05$, $f_{\text{PS}}^{\text{POL}} = 1.05$, and $f_{\text{PS}}^{\text{T}} = 1.07$. For $L < 100$, the mean of the simulation spectra with the source mask removed is 5%–10% below the input spectrum. Therefore, this MC bias correction is not applicable for multipoles below 100, and we do not report results below $L$ of 100.

In constructing our simulations and applying the quadratic estimator to recover the lensing potential, we have assumed $\phi$

to be Gaussian. However, nonlinear structure growth and post-Born lensing would introduce non-Gaussianities to $\phi$ (Böhm et al. 2016; Pratten & Lewis 2016). These non-Gaussianities produce the so-called $N_L^{(3/2)}$ bias, as studied in Böhm et al. (2016, 2018) and Beck et al. (2018). For the $L$ range considered, the size of this bias is $\sim 0.5\%$ for the temperature reconstruction and is negligible for the EB reconstruction for input CMB maps similar in noise levels and multipole range to those in this work (Beck et al. 2018; Böhm et al. 2018). These are small compared with our other sources of uncertainties, and we neglect them in our results.

We report our results in binned bandpowers. First, we derive the per-bin amplitude as the ratio of the unbiased lensing spectrum to the input theory spectrum:

$$A_b^{UVXY} = \frac{C_b^{\phi^{UV}\phi^{XY}}}{C_b^{\phi\phi,\text{theory}}}, \qquad (11)$$

with the subscript $b$ denoting a binned quantity; $C_b$ is a weighted average of the $C_L$ inputs for $L$ inside the boundaries of the bin, with the weights $w$ designed to maximize S/N: $w_L^{UV,XY} = C_L^{\phi\phi,\text{theory}}/\text{Var}(C_L^{\hat{\phi}^{UV}\hat{\phi}^{XY}})$. We obtain the variance from the corresponding set of simulation cross spectra.

We report the bandpowers in lensing convergence ($\kappa$) instead of the lensing potential $\phi$. The convergence field is $-1/2$ of the divergence of the deflection field, which is the gradient of the lensing potential $\phi$ (Lewis & Challinor 2006):

$$\kappa = -\frac{1}{2} \nabla^2 \phi. \qquad (12)$$

In Fourier space, they are related by $\kappa_L = (L(L+1))\phi_L/2$. The reported bandpowers are derived as the product of the data amplitude $A_b$ and the input theory spectrum at bin center $L_b$,

$$\hat{C}_{L_b}^{\kappa\kappa} \equiv \frac{(L_b(L_b+1))^2}{4} A_b C_{L_b}^{\phi\phi,\text{theory}}. \qquad (13)$$

The overall lensing amplitude for each estimator is calculated identically as the per-bin amplitude in Equation (11) using the whole reported $L$ range.

Differences between the fiducial cosmology and the cosmology of the SPTpol patch would produce different measured lensing amplitudes. The cosmology dependence enters through the $N_L^{(1)}$ bias. In this work, because we choose a fiducial cosmology that is consistent with the data, we expect the difference in the $N_L^{(1)}$ bias to be small. To test this, we sample the lensing amplitude given the fiducial cosmology with and without corrections to $N_L^{(1)}$. We find the difference in the lensing amplitude for the MV estimator to be 0.007 (0.1$\sigma$). We therefore neglect this correction in this work.

## 4. Results

In the following, we discuss our main results: the lensing convergence map, the lensing spectrum, and the lensing amplitude measurement. We compare our MV map to S15's. We discuss the relative weights of the $T$ and the POL results, compare the sizes of the systematic uncertainties to the statistical uncertainties, and put our measurement in the context of other lensing measurements.

In Figure 2, we show the MV, $T$, and POL lensing convergence. The maps are smoothed by a 1° FWHM Gaussian to highlight the higher S/N modes at the larger scales. At the
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Lensing $\kappa$ map reconstructed from the SPTpol 500 deg$^2$ field data, smoothed by a 1° FWHM Gaussian to highlight the signal-dominated modes. Top: the $\kappa$ map from the MV lensing estimator, which combines all temperature and polarization information. Bottom: $\kappa$ maps from the temperature estimator ($T$, left) and from combining the polarization estimators EE and EB (POL, right). At the angular scales shown, the POL estimator recovers the lensing potential with slightly higher S/N than the $T$ estimator. Therefore, it has higher weight in the MV combination and traces the fluctuations of the MV $\kappa$ map with higher fidelity.</figcaption></figure>

angular scales shown, the POL reconstruction has higher S/N (lower $N_L^{(0)}$) than the $T$ reconstruction. Therefore, while the $T$ and the POL $\kappa$ map fluctuations both trace those in the MV map in broad strokes, we see that the POL modes trace the MV modes more faithfully. The lensing modes in the MV map are reconstructed with S/N better than unity for $L \lesssim 250$. It is the largest lensing map reconstructed with this S/N level from the CMB to date.

In Figure 3, we compare the MV $\kappa$ map from S15 and the one from this work over the same region of the sky. We observe that the two convergence maps show nearly identical degree-scale structure. The S15 data were taken before the SPTpol 500 deg$^2$ survey data (see Section 2.1) and thus constitute an independent data set from that used in this work. In addition, because the two analyses have similar per-mode reconstruction noise, the visual agreement of the modes between the two is a useful consistency check.

We present the lensing power spectrum measurement in logarithmically spaced bins in the range $100 < L < 2000$. We

list the MV lensing bandpowers and their uncertainties in Table 1. The lower bound of the $L$ range is chosen by the region of validity of the MC bias correction. The upper bound is set by computing the uncertainties on the lensing amplitudes as we include higher multipoles and seeing no gain in S/N of the amplitude going beyond $L$ of 2000. In Figure 4, we show the bandpowers from the MV spectrum, the POL spectrum, and the $T$ spectrum. We see that the error bars of the POL spectrum are smaller than those of the $T$ spectrum at $L \lesssim 600$ and vice versa on smaller angular scales. This is consistent with the $N_L^{(0)}$ of the $T$ reconstruction at small angular scales being lower than that of the POL reconstruction—the small angular scale modes are better reconstructed by the $T$ estimator.

We measure the overall lensing amplitude for each estimator and find

$$A_{\text{MV}} = 0.944 \pm 0.058 \text{ (Stat.)} \pm 0.025 \text{ (Sys.),}$$

$$A_{\text{POL}} = 0.906 \pm 0.090 \text{ (Stat.)} \pm 0.040 \text{ (Sys.),}$$

$$A_{\text{T}} = 0.835 \pm 0.103 \text{ (Stat.)} \pm 0.012 \text{ (Sys.).}$$
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison with the S15 $\kappa$ map: Left: the $\kappa$ map reconstructed by S15 using the MV estimator with maps from SPTpol 100 deg$^2$ field observations (the sign is opposite from that shown in S15 because of a sign error in plotting in S15). Right: cutout of the 100 deg$^2$ field from the MV $\kappa$ map in this analysis, overlaid with contours from the S15 MV $\kappa$ at [-0.048, -0.032, -0.016, 0, 0.016, 0.032, 0.048]. Both maps are smoothed by a 1° FWHM Gaussian. Data that enter the S15 reconstruction were taken before the 500 deg$^2$ survey and thus are independent from the data used in this analysis. Because the modes shown are measured to S/N greater than unity, the fluctuations of the two maps visually tracing each other serve as a consistency check.</figcaption></figure>

<table><caption>**Table 1**
MV Lensing Bandpowers</caption><thead><tr><th>[L<sub>min</sub></th><th>L<sub>max</sub></th><th>L<sub>b</sub></th><th>10<sup>7</sup> C&#770;<sup>&kappa;&kappa;</sup><sub>L<sub>b</sub></sub></th></tr></thead><tbody><tr><td>[100]</td><td>133]</td><td>117</td><td>1.144 &plusmn; 0.230</td></tr><tr><td>[134]</td><td>181]</td><td>158</td><td>1.041 &plusmn; 0.161</td></tr><tr><td>[182]</td><td>244]</td><td>213</td><td>0.876 &plusmn; 0.115</td></tr><tr><td>[245]</td><td>330]</td><td>288</td><td>0.509 &plusmn; 0.080</td></tr><tr><td>[331]</td><td>446]</td><td>389</td><td>0.423 &plusmn; 0.058</td></tr><tr><td>[447]</td><td>602]</td><td>525</td><td>0.260 &plusmn; 0.042</td></tr><tr><td>[603]</td><td>813]</td><td>708</td><td>0.198 &plusmn; 0.032</td></tr><tr><td>[814]</td><td>1097]</td><td>956</td><td>0.086 &plusmn; 0.025</td></tr><tr><td>[1098]</td><td>1481]</td><td>1290</td><td>0.081 &plusmn; 0.023</td></tr><tr><td>[1482]</td><td>1998]</td><td>1741</td><td>0.053 &plusmn; 0.025</td></tr></tbody></table>

**Note.** The bandpowers for the MV spectrum as defined in Equation (13) and shown in Figure 4. Bins are evenly spaced in log(L), and bandpowers are reported at the center of each bin.

We derive the statistical uncertainties from the standard deviations of the lensing amplitudes from simulations. We detail the sources of systematic uncertainties in Section 5.2 and break them down in Table 3. Considering statistical uncertainties alone, we measure the lensing amplitude with 6% uncertainty using the MV estimator and with 10% uncertainty with the POL estimator. For the *T* estimator, we measure the lensing amplitude with 12% uncertainty. Having chosen the same cuts in multipole space for both the input temperature and polarization maps, this shows that the S/N per mode in the input polarization maps are now high enough that the POL estimators give more stringent measurements of the lensing amplitude than the *T* estimator. In future analyses, the *T* lensing spectrum is sample variance limited and cannot be improved by lowering the temperature map noise levels. Instead, it can be improved by including information from higher multipoles and/or more sky area. However, lowering the noise levels of the polarization maps can still improve the lensing measurement from polarization estimators. Specifically, unlike the temperature estimator, the $N_L^{(0)}$ of the EB estimator is not limited by unlensed power in the map, because there is little unlensed *B* mode power to contribute to $N_L^{(0)}$ in the multipole

range important for lensing reconstruction. In addition to surpassing the measurement uncertainty of the *T* lensing amplitude, considering statistical uncertainties alone, our POL lensing amplitude is the most precise amplitude measurement (10.1σ) using polarization data alone to date.

The systematic uncertainties for the MV and the POL estimators are ~40% of their respective statistical uncertainties, whereas the systematic uncertainty is subdominant for the *T* estimator compared with its statistical uncertainty. For both the MV and the POL estimators, the systematic uncertainty budget is dominated by the $P_{cal}$ uncertainty (Section 5.2). Including the systematic uncertainties in the MV amplitude measurement, we measure $A_{MV}$ with 7% uncertainty.

We detect lensing at very high significance. From reconstructing $\phi$ using 400 unlensed simulations, the standard deviation of $A_{MV}^{\text{unl}}$ is 0.024. The observed amplitude of $A_{MV} = 0.944$ would thus correspond to a 39σ fluctuation.

Compared with other ground-based measurements, our result has the tightest constraint on the lensing amplitude. In Figure 5, we show our lensing power spectrum measurement against previous measurements. Our measurement is consistent with the measurement by Omori et al. (2017). In that work, they reconstruct lensing using a combined temperature map from SPT-SZ and Planck over the common 2500 deg$^2$ of sky. They measure the lensing amplitude to be 0.95 ± 0.06 relative to the best-fit ΛCDM model to the Planck 2015 plikHM_TT_low-TEB_lensing data set (same as the fiducial cosmology used in this work). The most recent lensing analysis of all-sky Planck data found the best-fit lensing amplitude to be 1.011 ± 0.028 against the Planck 2018 TTTEE lowE lensing cosmology (Planck Collaboration et al. 2018b). To compare our measurement with this model, we refit our MV bandpowers and get $A_{MV} = 0.946 \pm 0.058$ (Stat.) ± 0.025 (Sys.), consistent with Planck’s lensing measurement.

## 5. Null Tests, Consistency Checks, and Systematic Uncertainties

In this section, we summarize the null and consistency tests we have performed on the data and account for the systematic uncertainties in our lensing amplitude measurements. We report test results from the MV, POL, and *T* estimators.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Lensing convergence bandpowers estimated from SPTpol 500 deg² field data. We show bandpowers from the MV, POL, and *T* estimators. The POL and *T* bandpowers are shifted in *L* for clarity. The *T* and POL bandpowers are consistent with each other given the error bars of the bandpowers. The reconstruction noise of the POL estimator is lower than that of *T* for *L* ≲ 600 and vice versa on smaller angular scales. This provides a sense of the angular scales at which each estimator gives better S/N. The black solid line shows the lensing convergence spectrum from the fiducial cosmology.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Lensing power spectrum measurements from this work and S15 (MV) (Story et al. 2015), POLARBEAR (POL) (POLARBEAR Collaboration et al. 2014), BICEP2/Keck (POL) (BICEP2 Collaboration et al. 2016), ACTPol (MV) (Sherwin et al. 2017), SPT-SZ + *Planck* (*T*) (Omori et al. 2017), and *Planck* (MV) (Planck Collaboration et al. 2018b). The gray solid line is the lensing spectrum from the best-fit ΛCDM model to the *Planck* plikHM_TT_lowTEB_lensing data set. This work contains data from about five times more area than S15, and the sizes of the error bars of the signal-dominated multipoles reflect the decrease in sample variance. In addition, it is the tightest lensing amplitude measurement of all the lensing measurements made using only ground-based telescope data to date.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. Results of the power spectrum consistency tests and curl null test. The figure shows $C_{L}^{\kappa\kappa}$ bandpowers from different analysis choices and the curl test bandpowers. The MV bandpowers and error bars from the baseline analysis are plotted as boxes. The bandpowers calculated with different analysis choices are consistent with the baseline MV bandpowers. The bandpowers from the curl null test are consistent with zero.</figcaption></figure>

## 5.1. Null Tests and Consistency Checks

We quantify the results of our tests using summary statistics comparing data bandpower differences with simulation bandpower-difference distributions, where the differences are taken between bandpowers obtained from the baseline analysis pipeline and from a pipeline with one analysis change. We plot these bandpowers in Figure 6 to provide an absolute sense of how much the bandpowers and their error bars change given the various analysis variations.

Quantitatively, we calculate the $\chi^2$ of the data-difference spectrum ($\Delta C_{b,\text{data}}^{\phi\phi}$) against the mean of the simulation-difference spectrum ($\langle\Delta C_{b,\text{sim}}^{\phi\phi}\rangle$) using the variance of the simulation-difference spectra ($\sigma_{b,\text{sys}}^2$), as the difference bandpowers are largely uncorrelated:

$$ \chi_{\text{sys}}^2 = \sum_b \frac{(\Delta C_{b,\text{data}}^{\phi\phi} - \langle\Delta C_{b,\text{sim}}^{\phi\phi}\rangle)^2}{\sigma_{b,\text{sys}}^2}. \quad (14) $$

The probability-to-exceed (PTE) is calculated from a $\chi^2$ distribution with 10 degrees of freedom, as we have 10 bandpowers. For the curl test, $\sigma_{b,\text{sys}}$ comes from the distribution of the simulation curl spectra. Figure 7 provides a visual summary of these tests. It shows the data-difference bandpowers for each test for the three estimators. The error bars are generated from the distributions of simulation-difference bandpowers. The data points and the error bars from each estimator in each L bin are scaled by the $1\sigma$ lensing spectrum uncertainties from that estimator in that bin. We list the $\chi^2$ and PTE values for each test and each estimator in Table 2.

We compare the differences in overall lensing amplitude similarly. We calculate the difference in the lensing amplitudes

between the baseline and the alternate analysis choice setup: $\Delta A_{\text{data}}$. We form the $\chi^2$ by comparing the data-difference amplitude against the distribution of the difference amplitudes in simulations.

**Varying $\ell_{\text{xmin}}$ and $\ell_{\text{max}}$:** We vary the multipole range of the input CMB maps used for lensing reconstruction to check the following: (1) consistency of bandpowers as we include more or fewer CMB modes and (2) impact of foregrounds as we increase the maximum multipole. As we increase the maximum multipole used in the input temperature map, the contamination of the CMB by foregrounds like tSZ and CIB increases. As discussed in Section 3.3, both of these inputs can bias the lensing spectrum. On the low $\ell$ side, we remove $|\ell_x| < \ell_{\text{xmin}}$ modes because the combination of our observing strategy and time stream filtering removes modes approximately along the $\ell_x = 0$ axis. We run two cases for the $\ell_{\text{xmin}}$ cut: $\ell_{\text{xmin}} = 50$ and $\ell_{\text{xmin}} = 200$. From Table 2, we see that the differences of the $\ell_{\text{xmin}}$ cases are consistent with the expected simulation distribution. We can in principle set the $\ell_{\text{xmin}}$ in the baseline analysis to be 50. However, because the applied time stream filtering removes $\ell \lesssim 100$ modes, and the number of modes available between $\ell_{\text{xmin}}$ of 50 and 100 is small compared with the number of modes at the high-$\ell$ end, we do not expect there to be significant improvement in S/N and therefore choose to set $\ell_{\text{xmin}} = 100$ for the baseline analysis. For the $\ell_{\text{max}}$ tests, we vary the $\ell_{\text{max}}$ cut between $\ell = 2500$ and $\ell = 3500$ in steps of 500. We observe from Table 2 that moving between $\ell_{\text{max}}$ of 3000 and 2500 and between $\ell_{\text{max}}$ of 3000 and 3500, the data bandpower differences are consistent with bandpower differences from simulations for all three sets of estimators. We set the baseline $\ell_{\text{max}}$ to be 3000 to keep the systematic uncertainty due to the subtraction of foreground biases subdominant to the statistical uncertainty of the measurements. For future analyses, estimators designed to reduce foreground biases (Osborne et al. 2014; Madhavacheril & Hill 2018; Schaan & Ferraro 2019) can be employed, potentially sacrificing some statistical power for reduced systematic uncertainty.

**Source masking threshold:** We check the impact of extragalactic foregrounds on the lensing measurements by varying the masking thresholds of point sources and galaxy clusters. We raise the thresholds of the point source and the cluster masks so that about 20% more sources are kept in the input maps. The modified cuts correspond to a flux threshold of 7.5 mJy and a cluster detection significance threshold of 4.7, whereas the baseline thresholds (see Section 2.2) are 6 mJy and 4.5 (i.e., the threshold for inclusion in the Bleem et al. 2015 catalog). From the bandpower PTEs and amplitude PTEs listed in Table 2, we see that the differences in bandpowers and amplitudes are consistent with expected variations from simulations. Therefore, we conclude that the masking thresholds we apply in the baseline analysis reduce foreground biases to the lensing spectrum sufficiently.

**Apodization:** We test for the effects of using a top hat function for the boundary and source mask by redoing the entire analysis using a cosine taper instead. The radius of the cosine taper is 10' for the boundary mask and 5' for the sources. The data difference is typical given the simulation-difference distribution for all three estimators.

**Curl test:** The deflection field $d(\hat{n})$ that remaps the primordial CMB anisotropies can be generically decomposed into a gradient and a curl component:

$$ d(\hat{n}) = \nabla\phi(\hat{n}) + \star\nabla\Omega(\hat{n}), \quad (15) $$
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The difference bandpowers ($\Delta C_L^{KK}$) and their uncertainties between the baseline analysis and analyses with one change as indicated in the labels, scaled by the uncertainties of the respective MV, *T*, and POL bandpowers. The $1\sigma$ bands of the MV, *T*, and POL estimators are shaded in gray. This is a visual summary of the $\chi^2$ values listed in Table 2. We find the analysis-variation bandpowers and best-fit lensing amplitudes to be consistent with the baseline setup for all three estimators.</figcaption></figure>

<table><caption>**Table 2**
Consistency and Null Tests $\chi^2$ and PTE</caption><thead><tr><th>Test Name</th><th>χ<sup>2</sup><sub>MV</sub></th><th>(PTE)</th><th>ΔA<sub>MV</sub><br>±var(ΔA<sub>MV</sub>)</th><th>(PTE)</th><th>χ<sup>2</sup><sub>T</sub></th><th>(PTE)</th><th>ΔA<sub>T</sub><br>±var(ΔA<sub>T</sub>)</th><th>(PTE)</th><th>χ<sup>2</sup><sub>POL</sub></th><th>(PTE)</th><th>ΔA<sub>POL</sub><br>±var(ΔA<sub>POL</sub>)</th><th>(PTE)</th></tr></thead><tbody><tr><td>ℓ<sub>xmin</sub> = 50</td><td>7.2</td><td>(0.71)</td><td>-0.002 ± 0.013</td><td>(0.76)</td><td>5.8</td><td>(0.83)</td><td>0.005 ± 0.019</td><td>(0.86)</td><td>9.7</td><td>(0.47)</td><td>-0.015 ± 0.025</td><td>(0.55)</td></tr><tr><td>ℓ<sub>xmin</sub> = 200</td><td>4.3</td><td>(0.93)</td><td>0.017 ± 0.021</td><td>(0.44)</td><td>10.4</td><td>(0.40)</td><td>0.035 ± 0.026</td><td>(0.17)</td><td>3.3</td><td>(0.97)</td><td>-0.005 ± 0.037</td><td>(0.90)</td></tr><tr><td>ℓ<sub>max</sub> = 2500</td><td>14.2</td><td>(0.17)</td><td>-0.000 ± 0.035</td><td>(0.98)</td><td>12.0</td><td>(0.28)</td><td>0.040 ± 0.093</td><td>(0.65)</td><td>5.9</td><td>(0.83)</td><td>-0.021 ± 0.018</td><td>(0.25)</td></tr><tr><td>ℓ<sub>max</sub> = 3500</td><td>17.2</td><td>(0.07)</td><td>0.014 ± 0.023</td><td>(0.66)</td><td>8.2</td><td>(0.61)</td><td>0.085 ± 0.055</td><td>(0.17)</td><td>11.0</td><td>(0.36)</td><td>-0.006 ± 0.008</td><td>(0.48)</td></tr><tr><td>Mask thres.</td><td>13.2</td><td>(0.21)</td><td>0.006 ± 0.008</td><td>(0.43)</td><td>15.6</td><td>(0.11)</td><td>0.021 ± 0.014</td><td>(0.10)</td><td>14.6</td><td>(0.15)</td><td>-0.002 ± 0.013</td><td>(0.88)</td></tr><tr><td>Apod. mask</td><td>4.8</td><td>(0.91)</td><td>0.004 ± 0.012</td><td>(0.90)</td><td>13.3</td><td>(0.20)</td><td>0.035 ± 0.022</td><td>(0.13)</td><td>3.7</td><td>(0.96)</td><td>-0.002 ± 0.022</td><td>(0.85)</td></tr><tr><td>Curl</td><td>9.3</td><td>(0.50)</td><td>-0.007 ± 0.005</td><td>(0.20)</td><td>7.5</td><td>(0.67)</td><td>-0.013 ± 0.013</td><td>(0.29)</td><td>9.4</td><td>(0.49)</td><td>-0.010 ± 0.009</td><td>(0.28)</td></tr></tbody></table>

**Note.** Results of systematics tests. For each test, we show the $\chi^2$ and PTE of the difference bandpowers and the difference amplitudes and its associated PTEs for each of the MV, *T*, and POL estimators.

where ⋆ is a 90° rotation operator and $\Omega(\hat{n})$ is a pseudo scalar field that sources the curl component (see e.g., Hirata & Seljak 2003; Namikawa et al. 2012). In the Born approximation, the lensing potential only sources the gradient component of the deflection field. At current noise levels and barring unknown physics, the curl component is expected to be consistent with zero. Therefore, estimating the curl of the deflection field serves as a test of the existence of any components in the data generated from non-Gaussian secondary effects or foregrounds (Cooray et al. 2005). Analogous to estimating the lensing $\phi$ field, for the curl component, we estimate the pseudo scalar field $\Omega$. The curl estimator is orthogonal to the lensing (gradient) estimator and has the same form as the lensing estimator Equation (1). The weights of the curl estimator are designed to have response to the curl instead of the gradient. We implement the weights as presented in Namikawa et al. (2012). The curl spectrum is derived

analogous to the lensing spectrum with two differences: (1) the theory input is set to a flat spectrum $C_L^{\Omega\Omega} = 10^{-7}$, where it is used for uniformly weighting the modes when binning and as a reference spectrum for the amplitude calculation; (2) no response correction from simulations is applied to the pseudo scalar field as the expected signal is zero. The bottom panel of Figure 7 shows the power spectra for the MV, POL, and *T* curl estimators and their uncertainties scaled by the respective uncertainties of their lensing spectra. The bandpower PTEs and the amplitude PTEs for the three estimators are listed in Table 2. They are all consistent with zero. The bandpower PTE and the amplitude PTE are 0.50 and 0.20 for the MV curl estimator. We thus see no evidence of contamination to our lensing estimate from non-Gaussian secondaries.

**Consistency across estimators:** In Figure 8, we show the reconstructed lensing spectrum from all 15 pairs of the five estimators *TT*, *TE*, *EE*, *EB*, and *TB* and the MV spectrum. We
<figure><img src="image_9.png" /><figcaption>Figure 8. Comparison of cross spectra of individual estimators and the MV combination MV × MV. The lensing spectrum from the fiducial cosmology, the Planck plikHM_TT_lowTEB_lensing ΛCDM best fit, is plotted in solid black lines. The best fit from each estimator is plotted in dashed gray lines. The amplitudes are computed relative to the fiducial lensing spectrum. The data points and error bars are scaled by 0.5 for the TB × TB spectrum to allow all the points to be shown. We find a PTE = 0.44 when comparing the 150 bandpowers against the best fit of the MV spectrum using a covariance matrix constructed from simulations.</figcaption></figure>

test for consistency of the lensing spectra from the 15 pairs of estimators with the best fit from the MV estimator by calculating χ² and PTE. We calculate the χ² of the 150 bandpowers against the binned theory spectrum scaled by the best-fit MV amplitude:

$$ \chi^2 = (d - m) \text{Cov}^{-1}(d - m)^{\dagger}, \qquad (16) $$

where d is the 150 data bandpowers, m is 15 copies of the scaled binned theory spectrum, and Cov denotes the covariance matrix. We construct the covariance matrix of the 150 bandpowers of 15 estimators from 400 simulations. We set the off-diagonal terms of each subblock off the main diagonal to zero as we expect there to be and have verified that there is little to no correlation across different L bins between estimators. For the main diagonal, the first and second bins are correlated at the ~10% level for 12 of the estimators, and we keep them while setting the rest of the covariance elements to zero. The PTE compared with the set of 400 sets of

simulated bandpowers is 0.44. Comparing instead the best-fit amplitudes of the 15 pairs of estimators to the MV best-fit amplitude, we note that the EB × TB and EB × EB pairs are low by 2.4σ and 2.1σ, respectively, while the rest of the pairs are within 2σ of the MV's amplitude. With the number of tests we have performed, it is not unusual to see 2σ outliers.

## 5.2. Systematic Uncertainty

In this section, we summarize the sources of systematic uncertainty and our accounting of them in the lensing amplitude measurement. We quantify the uncertainties in the lensing amplitude measurement due to uncertainties in the beam measurement, temperature and polarization calibrations, $T \to P$ leakage correction, global polarization angle rotation, and the applied foreground template. We address the potential impact of non-Gaussian polarized foregrounds. The sources of systematic uncertainty and their respective impact on the lensing amplitude measurements are listed in Table 3.
<table><caption>**Table 3**
Systematic Uncertainties</caption><thead><tr><th>Type</th><th>ΔA<sub>MV</sub></th><th>ΔA<sub>POL</sub></th><th>ΔA<sub>T</sub></th></tr></thead><tbody><tr><td>ΔA<sub>beam</sub></td><td>0.008</td><td>0.010</td><td>0.005</td></tr><tr><td>ΔA<sub>cal</sub></td><td>0.023</td><td>0.039</td><td>0.008</td></tr><tr><td>ΔA<sub>T→P</sub></td><td>≤0.001</td><td>≤0.001</td><td>N/A</td></tr><tr><td>ΔA<sub>pol.rot.</sub></td><td>&lt;0.001</td><td>&lt;0.001</td><td>N/A</td></tr><tr><td>ΔA<sub>fg</sub></td><td>0.004</td><td>N/A</td><td>0.008</td></tr><tr><td>ΔA<sub>tot</sub></td><td>0.025</td><td>0.040</td><td>0.012</td></tr></tbody></table>

**Note.** The contributions to the systematic uncertainty budget. The total systematic uncertainties are evaluated by taking the quadrature sum of the individual contributions.

**Beam uncertainty:** We take the beam measurement uncertainty $\Delta B_{\ell}$ derived from the beam covariance matrix in H18 and convolve $(1+\Delta B_{\ell})$ with the input $\bar{X}$ maps while keeping all the simulation maps the same. We analyze the data maps as if their beam were the mean measured beam $B_{\ell}$ as opposed to the modified $B_{\ell}$ $(1+\Delta B_{\ell})$. The difference in the lensing amplitude between this analysis and the baseline analysis quantifies the effect of an underestimation of the beam profile by $1\sigma$ of the beam measurement uncertainty across the entire angular multipole range. We find $\Delta A_{MV} = 0.008$, $\Delta A_{POL} = 0.010$, and $\Delta A_{T} = 0.005$. The systematic uncertainties induced by the beam uncertainty are small $(\sim 0.1\sigma)$ compared with the respective statistical uncertainties of the measured lensing amplitudes in T, POL, and MV.

**Temperature and polarization calibrations:** We apply the temperature and polarization calibration factors derived in H18 to our data $T/Q/U$ maps as described in Section 2.2. From the posterior distributions of H18, we obtain the uncertainties of the $T_{cal}$ factor $\delta T_{cal}$ and the $P_{cal}$ factor $\delta P_{cal}$, which are 0.3% and 0.6%, respectively. The $T_{cal}$ and $P_{cal}$ are applied to the raw temperature and polarization maps through $T = T^{\text{raw}} \times T_{cal}$ and $(Q/U) = (Q/U)^{\text{raw}} \times T_{cal} \times P_{cal}$. Keeping the simulated maps fixed to the baseline analysis, we scale the data maps by $(1+\delta T_{cal})$ for the temperature map and $(1+\delta T_{cal})(1+\delta P_{cal})$ for the polarization maps and we recalculate the data lensing amplitudes.46 The two pieces in the amplitude calculation that change because of the modified data maps are $C_{L}^{\hat{\phi}\hat{\phi}}$ and $N_{L}^{(0),RD}$. To illustrate how they contribute to the overall change in the lensing amplitude, we consider here the temperature-only estimator. The term $C_{L}^{\hat{\phi}\hat{\phi},T}$, which contains four copies of the data map, would shift by $(1+\delta T_{cal})^4$. For $N_{L}^{(0),RD}$, the four terms where two of the four input maps are replaced by data maps get a multiplicative correction of $(1+\delta T_{cal})^2$. The overall shift in $\hat{C}_{L}^{\hat{\phi}\hat{\phi},T}$ is therefore

$$ \Delta \hat{C}_{L}^{\hat{\phi}\hat{\phi}} = \Delta C_{L}^{\hat{\phi}\hat{\phi}} - \Delta N_{L}^{(0),RD}, \quad (17) $$

where

$$ \Delta C_{L}^{\hat{\phi}\hat{\phi}} = [(1 + \delta T_{cal})^4 - 1] C_{L}^{\hat{\phi}\hat{\phi}} \approx 4\delta T_{cal} C_{L}^{\hat{\phi}\hat{\phi}} \quad (18) $$

and

$$ \begin{aligned} \Delta N_{L}^{(0),RD} \approx & 2\delta T_{cal} \langle C_{L}^{\hat{\phi}\hat{\phi}}[\bar{T}_{d}, \bar{T}_{MC}, \bar{T}_{d}, \bar{T}_{MC}] \\ & + C_{L}^{\hat{\phi}\hat{\phi}}[\bar{T}_{MC}, \bar{T}_{d}, \bar{T}_{d}, \bar{T}_{MC}] \\ & + C_{L}^{\hat{\phi}\hat{\phi}}[\bar{T}_{d}, \bar{T}_{MC}, \bar{T}_{MC}, \bar{T}_{d}] \\ & + C_{L}^{\hat{\phi}\hat{\phi}}[\bar{T}_{MC}, \bar{T}_{d}, \bar{T}_{MC}, \bar{T}_{d}] \rangle_{\text{MC}}, \end{aligned} \quad (19) $$

with MC denoting simulation realizations and d denoting data. The POL case has a similar form as the temperature-only case with $(1+\delta T_{cal}) \rightarrow (1+\delta T_{cal})(1+\delta P_{cal})$. The difference in the resultant measured amplitudes due to an offset in the calibration factors depends on the relative amplitudes of $C_{L}^{\hat{\phi}\hat{\phi}}$ and $N_{L}^{(0),RD}$ at each multipole L. Therefore, we quantify the difference by running the baseline analysis with the temperature and polarization calibration of the data maps shifted by $1\sigma$. We find that $\Delta A_{MV} = 0.023$, $\Delta A_{POL} = 0.039$, and $\Delta A_{T} = 0.008$. The shifts in the lensing amplitudes are dominated by the polarization calibration uncertainty. Furthermore, the calibration systematic uncertainty is almost half of the statistical uncertainty of the POL amplitude. Because the S/N of future lensing measurements will be driven by the polarization estimators, we will need more precise polarization calibration for the measurements to remain dominated by statistical uncertainties. This can be achieved by cross-calibrating with deeper CMB maps or over larger areas of sky, assuming an external CMB map exists with more accurate polarization calibration than Planck. Relatedly, lensing measures mode coupling, which in principle can be extracted irrespective of the input maps calibration. An example of circumventing the systematic uncertainty contribution caused by calibration uncertainties in the input maps is to use the measured power spectra in one of the weight functions in the analytic response in Equation (2). In this approach, the response moves together with the calibration of the input maps and therefore eliminates this systematic uncertainty.

**Temperature-to-polarization leakage:** We estimate the bias to the lensing amplitude measurement caused by misestimating the $T \rightarrow P$ leakage factors. Similar to quantifying the systematic offsets by the beam uncertainty and the $T_{cal}/P_{cal}$ uncertainties, here we modify the data polarization maps by over-subtracting a $\epsilon^{(Q/U)}$-scaled copy of the temperature map by $1\sigma$ in $\epsilon^Q$ and $\epsilon^U$, while keeping the rest of the analysis the same. We find the change in $A_{MV}$ and $A_{POL}$ to be negligible, at $0.001\sigma$ and $0.002\sigma$ respectively.

**Global polarization angle rotation:** There is a 6% uncertainty in the global rotation angle measured through the minimization of the TB and EB spectra discussed in Section 2.2. How much would an offset in the polarization angle rotation bias the lensing amplitude measurement? We run the baseline analysis with the data polarization maps rotated an extra 6% on top of the measured angle from the minimization procedure. We find that $A_{MV}$ and $A_{POL}$ change by less than $0.01\sigma$ of their respective statistical uncertainties.

**Extragalactic foregrounds:** As discussed in Section 3.3, we subtract templates of expected foreground biases from the T and MV lensing power spectra $\hat{C}_{L}^{\hat{\phi}\hat{\phi}}$ given models of CIB and tSZ from van Engelen et al. (2014). The templates are taken from the mean of the range of models. We estimate the

<footnote n="46">46 In S15, we accounted for the systematic uncertainties from calibration uncertainty as $4 \times \delta T_{cal}$ and $4 \times \delta P_{cal}$. This underestimates the systematic uncertainty because the object that gets scaled by $\delta T_{cal}$ and $\delta P_{cal}$ is $C_{L}^{\hat{\phi}\hat{\phi}}$ instead of $\hat{C}_{L}^{\hat{\phi}\hat{\phi}}$. This means that the calibration offset applies to the noise bias terms as well as the signal term and can be a factor of a couple larger than $\hat{C}_{L}^{\hat{\phi}\hat{\phi}}$ itself.</footnote>
systematic uncertainty from the template subtraction step by
measuring the lensing amplitude difference using the maximum
and minimum of template values allowed by the model space.
We find the lensing amplitudes $A_{\text{MV}}$ to shift by $\pm 0.004$ and $A_{T}$
to shift by $\pm 0.008$, both less than $0.1\sigma$ of their statistical
uncertainties. For the temperature estimator, this source of bias
is $\sim$1% of the lensing amplitude and is of the same magnitude
as the bias from the temperature calibration. For polarization,
we expect the only significant source of extragalactic fore-
grounds to come from the unclustered point-source component,
because SZ and the clustered CIB are negligibly polarized. The
SZ effects are expected to be polarized at the less than 1% level
for the cluster masking threshold of this work (Birkinshaw
1999; Carlstrom et al. 2002). The clustered CIB component is a
modulation of the mean power from all sources (Scott &
White 1999) and is thus effectively unpolarized. For the
unclustered point sources, we measure the mean squared
polarization fraction to be $\langle p^{2} \rangle = (0.95 \pm 0.11) \times 10^{-3}$ ($\sim$3%
polarized) at 150 GHz for sources above 6 mJy in the SPTpol
field (Gupta et al. 2019). Assuming the polarization fraction is
constant with flux density, the extragalactic foregrounds are
thus on average significantly less polarized than the CMB.
Because the level of bias introduced from temperature
foregrounds is at the percent level, we conclude that the bias
from polarized sources is negligible for this work.

Galactic foregrounds: Galactic foregrounds are non-Gaussian and can contribute to both the lensing estimator, imparting a bias to the lensing spectrum measurement, and the curl estimator, causing the curl null test to fail. We mitigate these potential effects by (1) observing a low-foreground patch and (2) suppressing power from angular scales that are foreground dominated. We use the curl null test to check for potential contamination from galactic foregrounds. The SPTpol 500 deg² field is one of the lowest galactic foreground regions (e.g., Planck Collaboration et al. 2018d). Galactic foreground power is fainter than the CMB in both temperature and in *E*-mode polarization for the angular scales considered. However, polarized thermal dust dominates the *B*-mode power on angular scales below multipole *l* of 200 (BICEP2 Collaboration et al. 2018). Challinor et al. (2018) showed that by including the foreground power spectrum as a noise term in the $C^{-1}$ filter of the input CMB fields, the polarization estimator can recover unbiased lensing amplitude to subpercent accuracy. In this work, galactic dust power in *EE* and *BB* at the same level as the simulation inputs are included as noise terms in the map filters. Finally, if galactic foregrounds had been a significant source of bias, we would fail the curl null tests. The data pass the curl null test for the *T*, POL and MV estimators, suggesting that biases from galactic foregrounds are subdominant at present.

We summarize the systematic uncertainties contributed from
each source addressed in this section in Table 3. We add the
differences in lensing amplitude from each source in quadrature
and include the result as a systematic uncertainty in our lensing
amplitude measurements. The total systematic uncertainties are
$A_{\text{MV}}^{\text{sys}} = 0.025$, $A_{\text{POL}}^{\text{sys}} = 0.040$, and $A_{T}^{\text{sys}} = 0.012$, respectively.
The *T* systematic uncertainty is subdominant given the
statistical uncertainties of the *T* lensing amplitude, while those
for the MV and POL reach ~0.4σ of their respective statistical
uncertainties.

6. Conclusions

In this work, we present a new measurement of the CMB lensing potential from the 500 deg² SPTpol survey. We use a quadratic estimator to extract the lensing potential $\phi$ from combinations of temperature and polarization maps. Our per-mode measurement has S/N similar to that in S15 because the noise levels of the input maps are similar—we measure $\phi$ with S/N better than unity in the multipole range $100 < L \lesssim 250$.

We measure the lensing amplitude with our MV combina-
tion of estimators to be $A_{\text{MV}} = 0.944 \pm 0.058$ (Stat.) $\pm$
0.025 (Sys.). Of the published lensing measurements using
data only from ground-based telescopes, this is the tightest
measurement of the lensing amplitude to date. Using unlensed
simulations to estimate the probability of the measurement, we
find that in the absence of true CMB lensing our measurement
would be a 39$\sigma$ outlier—that is, we detect lensing at very high
significance. We make individual combinations of temperature
and polarization estimators. Considering statistical uncertain-
ties alone, our POL lensing constraint is the most precise
among measurements of its kind to date. We perform null tests
and consistency checks on our results and find no evidence for
significant contamination. We estimate the size of the
systematic uncertainty in our lensing measurement from
uncertainties in calibration, beam measurement, $T \rightarrow P$ leakage
correction, global polarization rotation, and subtraction of
extragalactic foreground bias. We find that our systematic
uncertainty is nearly half as large as our statistical uncertainty
and dominated by the uncertainty in the polarization calibra-
tion. Looking forward, we will have to either improve the
precision of the polarization calibration or find solutions to self-
calibrate the lensing estimator in order to reduce this systematic
uncertainty.

The lensing spectrum measurement in this work is at
sufficiently high precision to provide relevant independent
constraints on cosmological parameters like $\Omega_M$, $\sigma_8$, and the
sum of neutrino masses. We will report cosmological parameter
constraints in a future paper. While this work represents the
first analysis in which the lensing amplitude is better
constrained by polarization, rather than temperature data, this
will become typical for future experiments as map noise levels
continue to decrease. The *EB* estimator will eventually become
the most constraining CMB lensing estimator, at least for
lensing multipoles below several hundred. For the expected
survey map depths of the currently operating SPT-3G
experiment (3.0, 2.2, and 8.8 $\mu$K-arcmin at 95, 150, and
220 GHz, respectively, Bender et al. 2018), the *EB* estimator is
projected to provide the highest S/N lensing amplitude
measurement of all pairs of input CMB maps.

The advantages of using polarization maps for CMB lensing
measurements are as follows. First, foregrounds in CMB
temperature maps, for example, tSZ, are the main sources of
bias in cross correlations with other dark matter tracers (Abbott
et al. 2019; Namikawa et al. 2019). Polarization maps are less
contaminated by extragalactic foregrounds, and thus more
modes on small angular scales (high *l*) can be used for lensing
reconstruction and for cross correlation, improving the S/N of
the measurement in addition to reducing foreground biases.
Second, lensing estimators are typically limited by the $N_L^{(0)}$
noise contributed by the unlensed power in the input maps. The
*EB* estimator, however, is not as limited because unlensed *B*
power is subdominant to lensing *B* power on subdegree scales.
Furthermore, to first order in $\phi$, lensing *B* modes are sourced
completely by *E* modes. These factors make methods beyond
the quadratic estimator that approach the optimal solution
(Millea et al. 2019) particularly efficient in improving the S/N
of the *EB* estimator.

High-S/N lensing measurements are essential for constrain-
ing the sum of neutrino masses and the amplitude of primordial
gravitational waves, parameterized by the tensor-to-scalar ratio
*r*. In the coming few years, the uncertainty on *r* from upcoming
CMB experiments will be limited by the variance of the lensing
*B* modes. To make progress in constraining *r*, we need lensing
measurements with high S/N per mode for delensing (Man-
zotti et al. 2017). For example, delensing with the lensing map
from SPT-3G will be crucial for the BICEP Array experi-
ment (Hui et al. 2018) to reach its projected uncertainty on *r* of
*σ*(*r*) ~ 0.003, in which delensing improves *σ*(*r*) by a factor of
about 2.5. For the next-generation ground-based CMB
experiment CMB-S4, delensing is even more crucial for
achieving the projected *r* constraint of *σ*(*r*) = 5 × 10⁻⁴
(CMB-S4 Collaboration et al. 2016).

The authors would like to acknowledge helpful comments
from Chang Feng and Srinivasan Raghunathan on the manu-
script. The South Pole Telescope program is supported by the
National Science Foundation through grant PLR-1248097.
Partial support is also provided by the NSF Physics Frontier
Center grant PHY-0114422 to the Kavli Institute of Cosmo-
logical Physics at the University of Chicago, the Kavli
Foundation, and the Gordon and Betty Moore Foundation
through grant GBMF No. 947 to the University of Chicago.
This work is also supported by the U.S. Department of Energy.
W.L.K.W. is supported in part by the Kavli Institute for
Cosmological Physics at the University of Chicago through
grant NSF PHY-1125897 and an endowment from the Kavli
Foundation and its founder Fred Kavli. J.W.H. is supported by
the National Science Foundation under award No. AST-
1402161. C.R. acknowledges support from an Australian
Research Council Future Fellowship (FT150100074). B.B. is
supported by the Fermi Research Alliance LLC under contract
No. De-AC02-07CH11359 with the U.S. Department of
Energy. The Cardiff authors acknowledge support from the
UK Science and Technologies Facilities Council (STFC). The
CU Boulder group acknowledges support from NSF AST-
0956135. The McGill authors acknowledge funding from the
Natural Sciences and Engineering Research Council of Canada,
Canadian Institute for Advanced Research, and the Fonds de
Recherche du Québec—Nature et technologies. The UCLA
authors acknowledge support from NSF AST-1716965 and
CSSI-1835865. Work at Argonne National Lab is supported by
UChicago Argonne LLC, Operator of Argonne National
Laboratory (Argonne). Argonne, a U.S. Department of Energy
Office of Science Laboratory, is operated under contract No.
DE-AC02-06CH11357. We also acknowledge support from
the Argonne Center for Nanoscale Materials. This research
used resources of the National Energy Research Scientific
Computing Center (NERSC), a U.S. Department of Energy
Office of Science User Facility operated under contract No.
DE-AC02-05CH11231. The data analysis pipeline also uses the
scientific python stack (Jones et al. 2001; Hunter 2007; van der
Walt et al. 2011) and the HDF5 file format (The HDF Group
1997). G. Longhi’s color scheme suggestions are gratefully
acknowledged.

ORCID iDs

W. L. K. Wu https://orcid.org/0000-0001-5411-6920

P. A. R. Ade https://orcid.org/0000-0002-5127-0401

A. N. Bender https://orcid.org/0000-0001-5868-0748

B. A. Benson https://orcid.org/0000-0002-5108-6823

F. Bianchini https://orcid.org/0000-0003-4847-3483

R. Citron https://orcid.org/0000-0001-8920-0356

T. M. Crawford https://orcid.org/0000-0001-9000-5013

T. de Haan https://orcid.org/0000-0001-5105-9473

G. P. Holder https://orcid.org/0000-0002-0463-6394

A. Manzotti https://orcid.org/0000-0003-0057-877X

Y. Omori https://orcid.org/0000-0002-0963-7310

C. L. Reichardt https://orcid.org/0000-0003-2226-9169

B. R. Saliwanchik https://orcid.org/0000-0002-5089-7472

G. Simard https://orcid.org/0000-0001-5800-3990

J. D. Vieira https://orcid.org/0000-0001-7192-3871

N. Whitehorn https://orcid.org/0000-0002-3157-0407

References

Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2018, PhRvD, 98, 043526

Abbott, T. M. C., Abdalla, F. B., Alarcon, A., et al. 2019, PhRvD, 100, 023541

Beck, D., Fabbian, G., & Errard, J. 2018, PhRvD, 98, 043512

Bender, A. N., Ade, P. A. R., Ahmed, Z., et al. 2018, Proc. SPIE, 10708, 1070803

BICEP2 Collaboration, Keck Array Collaboration, Ade, P. A. R., et al. 2016, ApJ, 833, 228

BICEP2 Collaboration, Keck Array Collaboration, Ade, P. A. R., et al. 2018, PhRvL, 121, 221301

Birkinshaw, M. 1999, PhR, 310, 97

Blanchard, A., & Schneider, J. 1987, A&A, 184, 1

Bleem, L. E., Stalder, B., de Haan, T., et al. 2015, ApJS, 216, 27

Böhm, V., Schmittfull, M., & Sherwin, B. D. 2016, PhRvD, 94, 043519

Böhm, V., Sherwin, B. D., Liu, J., et al. 2018, PhRvD, 98, 123510

Carlstrom, J. E., Ade, P. A. R., Aird, K. A., et al. 2011, PASP, 123, 568

Carlstrom, J. E., Holder, G. P., & Reese, E. D. 2002, ARA&A, 40, 643

Challinor, A., Allison, R., Carron, J., et al. 2018, JCAP, 2018, 018

CMB-S4 Collaboration, Abazajian, K. N., Adshead, P., et al. 2016, arXiv:1610.02743

Cooray, A., Kamionkowski, M., & Caldwell, R. R. 2005, PhRvD, 71, 123527

Crites, A. T., Henning, J. W., Ade, P. A. R., et al. 2015, ApJ, 805, 36

Das, S., Sherwin, B. D., Aguirre, P., et al. 2011, PhRvL, 107, 021301

George, E. M., Reichardt, C. L., Aird, K. A., et al. 2015, ApJ, 799, 177

Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759

Gupta, N., Reichardt, C. L., Ade, P. A. R., et al. 2019, arXiv:1907.02156

Hanson, D., Challinor, A., Efstathiou, G., & Bielewicz, P. 2011, PhRvD, 83, 043005

Henning, J. W., Sayre, J. T., Reichardt, C. L., et al. 2018, ApJ, 852, 97

Hikage, C., Oguri, M., Hamana, T., et al. 2019, PASJ, 71, 43

Hildebrandt, H., Köhlinger, F., van den Busch, J. L., et al. 2018, A&A, submitted (arXiv:1812.06076)

Hirata, C. M., & Seljak, U. 2003, PhRvD, 68, 083002

Hu, W., & Okamoto, T. 2002, ApJ, 574, 566

Hui, H., Ade, P. A. R., Ahmed, Z., et al. 2018, Proc. SPIE, 10708, 1070807

Hunter, J. D. 2007, CSE, 9, 90

Jones, E., Oliphant, T., Peterson, P., et al. 2001, SciPy: Open Source Scientific Tools for Python, http://www.scipy.org/

Kamionkowski, M., & Kovetz, E. D. 2016, ARA&A, 54, 227

Keating, B. G., Shimon, M., & Yadav, A. P. S. 2013, ApJL, 762, L23

Keisler, R., Hoover, S., Harrington, N., et al. 2015, ApJ, 807, 151

Kesden, M., Cooray, A., & Kamionkowski, M. 2003, PhRvD, 67, 123507

Lesgourgues, J., & Pastor, S. 2006, PhR, 429, 307

Lewis, A. 2005, PhRvD, 71, 083008

Lewis, A., & Challinor, A. 2006, PhR, 429, 1

Madhavacheril, M. S., & Hill, J. C. 2018, PhRvD, 98, 023534

Manzotti, A., Story, K. T., Wu, W. L. K., et al. 2017, ApJ, 846, 45

Millea, M., Anderes, E., & Wandelt, B. D. 2019, PhRvD, 100, 023509

Namikawa, T., Chinone, Y., Miyatake, H., et al. 2019, ApJ, 882, 62

Namikawa, T., Hanson, D., & Takahashi, R. 2013, MNRAS, 431, 609

Namikawa, T., Yamauchi, D., & Taruya, A. 2012, JCAP, 1, 7

Omori, Y., Chown, R., Simard, G., et al. 2017, ApJ, 849, 124

Osborne, S. J., Hanson, D., & Doré, O. 2014, JCAP, 3, 24
Padin, S., Staniszewski, Z., Keisler, R., et al. 2008, ApOpt, 47, 4418

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2014, A&A, 571, 17

Planck Collaboration, Ade, P. A. R., Aghanim, N., et al. 2016a, A&A, 594, A13

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018a, arXiv:1807.06209

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018b, arXiv:1807.06210

Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018c, arXiv:1807.06207

Planck Collaboration, Aghanim, N., Arnaud, M., et al. 2016b, A&A, 594, A11

Planck Collaboration, Akrami, Y., Ashdown, M., et al. 2018d, arXiv:1801.04945

POLARBEAR Collaboration, Ade, P. A. R., Akiba, Y., et al. 2014, PhRvL, 113, 021301

Pratten, G., & Lewis, A. 2016, JCAP, 2016, 047

Schaan, E., & Ferraro, S. 2019, PhRvL, 122, 181301

Scott, D., & White, M. 1999, A&A, 346, 1

Seiffert, M., Borys, C., Scott, D., & Halpern, M. 2007, MNRAS, 374, 409

Shaw, L. D., Nagai, D., Bhattacharya, S., & Lau, E. T. 2010, ApJ, 725, 1452

Sherwin, B. D., van Engelen, A., Sehgal, N., et al. 2017, PhRvD, 95, 123529

Smith, K. M., Zahn, O., & Doré, O. 2007, PhRvD, 76, 043510

Story, K. T., Hanson, D., Ade, P. A. R., et al. 2015, ApJ, 810, 50

The HDF Group 1997, Hierarchical Data Format v5 (Champaign, IL: HDF Group)

van der Walt, S., Colbert, S., & Varoquaux, G. 2011, CSE, 13, 22

van Engelen, A., Bhattacharya, S., Sehgal, N., et al. 2014, ApJ, 786, 13

van Engelen, A., Keisler, R., Zahn, O., et al. 2012, ApJ, 756, 142