# A demonstration of improved constraints on primordial gravitational waves with delensing

P.A.R. Ade, Z. Ahmed, M. Amiri, A.J. Anderson, J.E. Austermann, J.S.

Avva, D. Barkats, R. Basu Thakur, J.A. Beall, A.N. Bender, et al.

► To cite this version:

P.A.R. Ade, Z. Ahmed, M. Amiri, A.J. Anderson, J.E. Austermann, et al. A demonstration of improved constraints on primordial gravitational waves with delensing. Phys.Rev.D, 2021, 103 (2), pp.022004. 10.1103/PhysRevD.103.022004. hal-03047658

HAL Id: hal-03047658

https://hal.science/hal-03047658

Submitted on 3 Jul 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
A demonstration of improved constraints
on primordial gravitational waves with delensing

P. A. R. Ade,¹ Z. Ahmed,²,³ M. Amiri,⁴ A. J. Anderson,⁵ J. E. Austermann,⁶,⁷ J. S. Avva,⁸ D. Barkats,⁹ R. Basu Thakur,¹⁰

J. A. Beall,⁶ A. N. Bender,¹¹,¹² B. A. Benson,⁵,¹²,¹³ F. Bianchini,¹⁴ C. A. Bischoff,¹⁵ L. E. Bleem,¹¹,¹² J. J. Bock,¹⁰,¹⁶

H. Boenish,¹⁷ E. Bullock,¹⁸ V. Buza,⁹,¹⁷ J. E. Carlstrom,¹²,¹⁹,¹¹,¹³,²⁰ C. L. Chang,¹²,¹¹,¹³ J. R. Cheshire IV,¹⁸

H. C. Chiang,²¹,²² T.-L. Chou,¹²,¹⁹ R. Citron,²³ J. Connors,⁹,⁶ C. Corbett Moran,¹⁶ J. Cornelison,⁹ T. M. Crawford,¹²,¹³

A. T. Crites,¹²,¹³,²⁴,²⁵ M. Crumrine,²⁶ A. Cukierman,³,²,²⁷ T. de Haan,²⁸ M. Dierickx,⁹ M. A. Dobbs,²¹,²⁹ L. Duband,³⁰

W. Everett,³¹ S. Fatigoni,⁴ J. P. Filippini,³²,³³ S. Fliescher,²⁶ J. Gallicchio,¹²,³⁴ E. M. George,³⁵,⁸ T. St. Germaine,⁹

N. Goeckner-Wald,²⁷ D. C. Goldfinger,⁹,²⁷ J. Grayson,²⁷ N. Gupta,¹⁴ G. Hall,²⁶ M. Halpern,⁴ N. W. Halverson,³¹,⁷

S. Harrison,⁹ S. Henderson,³,² J. W. Henning,¹¹,¹² S. R. Hildebrandt,¹⁰,¹⁶ G. C. Hilton,⁶ G. P. Holder,³³,³²,²⁹

W. L. Holzapfel,⁸ J. D. Hrubes,²³ N. Huang,⁸ J. Hubmayr,⁶ H. Hui,¹⁰ K. D. Irwin,²,²⁷ J. Kang,²⁷ K. S. Karkare,⁹,¹²

E. Karpel,²⁷ S. Kefeli,¹⁰ S. A. Kernasovskiy,²⁷ L. Knox,³⁶ J. M. Kovac,⁹,¹⁷ C. L. Kuo,²⁷,³,² K. Lau,²⁶ A. T. Lee,⁸,³⁷

E. M. Leitch,¹²,¹³ D. Li,⁶,² A. Lowitz,¹³ A. Manzotti,¹²,³⁸ J. J. McMahon,¹²,¹⁹,¹³ K. G. Megerian,¹⁶ S. S. Meyer,¹²,¹⁹,¹³,²⁰

M. Millea,⁸ L. M. Mocanu,¹²,¹³,³⁹ L. Moncelsi,¹⁰ J. Montgomery,²¹ A. Nadolski,³³,³² T. Namikawa,⁴⁰ T. Natoli,¹³,¹²,²⁴

C. B. Netterfield,⁴¹,²⁹ H. T. Nguyen,¹⁶ J. P. Nibarger,⁶ G. Noble,²¹ V. Novosad,⁴² R. O’Brien,¹⁰,¹⁶ R. W. Ogburn IV,²⁷,³,²

Y. Omori,³,²⁷ S. Padin,¹²,¹³,¹⁰ S. Palladino,¹⁵ S. Patil,¹⁴ T. Prouve,³⁰ C. Pryke,²⁶,¹⁸ B. Racine,⁹ C. L. Reichardt,¹⁴

C. D. Reintsema,⁶ S. Richter,⁹ J. E. Ruhl,⁴³ B. R. Saliwanchik,⁴³,⁴⁴ K. K. Schaffer,¹²,²⁰,⁴⁵ A. Schillaci,¹⁰ B. L. Schmitt,⁹

R. Schwarz,²⁶ C. D. Sheehy,⁴⁶ C. Sievers,²³ G. Smecher,²¹,⁴⁷ A. Soliman,¹⁰ A. A. Stark,⁹ B. Steinbach,¹⁰

R. V. Sudiwala,¹ G. P. Teply,¹⁰,⁴⁸ K. L. Thompson,²⁷,³ J. E. Tolan,²⁷ C. Tucker,¹ A. D. Turner,¹⁶ C. Umiltà,³² T. Veach,⁴⁹

J. D. Vieira,³³,³² A. G. Vieregg,²⁰,¹² A. Wandui,¹⁰ G. Wang,¹¹ A. C. Weber,¹⁶ N. Whitehorn,⁵⁰ D. V. Wiebe,⁴

J. Willmert,²⁶ C. L. Wong,⁹,¹⁷ W. L. K. Wu,¹²,²,³,* H. Yang,²⁷ V. Yefremenko,¹¹ K. W. Yoon,²⁷,³,² E. Young,³,²,²⁷

C. Yu,²⁷ L. Zeng,⁹,¹⁷ and C. Zhang¹⁰

(The BICEP/Keck and SPTpol Collaborations)

¹Cardiff University, Cardiff CF10 3XQ, United Kingdom

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

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

⁴Department of Physics and Astronomy, University of British Columbia,
Vancouver, British Columbia, V6T 1Z1, Canada

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

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

⁷Department of Physics, University of Colorado, Boulder, Colorado 80309, USA

⁸Department of Physics, University of California, Berkeley, California 94720, USA

⁹Harvard-Smithsonian Center for Astrophysics, 60 Garden Street, Cambridge, Massachusetts 02138, USA

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

¹¹High Energy Physics Division, Argonne National Laboratory,
9700 S. Cass Avenue, Argonne, Illinois 60439, USA

¹²Kavli Institute for Cosmological Physics, University of Chicago,
5640 South Ellis Avenue, Chicago, Illinois 60637, USA

¹³Department of Astronomy and Astrophysics, University of Chicago,
5640 South Ellis Avenue, Chicago, Illinois 60637, USA

¹⁴School of Physics, University of Melbourne, Parkville, Victoria 3010, Australia

¹⁵Department of Physics, University of Cincinnati, Cincinnati, Ohio 45221, USA

¹⁶Jet Propulsion Laboratory, Pasadena, California 91109, USA

¹⁷Department of Physics, Harvard University, Cambridge, Massachusetts 02138, USA

¹⁸Minnesota Institute for Astrophysics, University of Minnesota,
Minneapolis, Minnesota 55455, USA

¹⁹Department of Physics, University of Chicago, 5640 South Ellis Avenue, Chicago, Illinois 60637, USA

²⁰Enrico Fermi Institute, University of Chicago, 5640 South Ellis Avenue, Chicago, Illinois 60637, USA

²¹Department of Physics, McGill University, 3600 Rue University, Montreal, Quebec H3A 2T8, Canada
<sup>22</sup>School of Mathematics, Statistics and Computer Science, University of KwaZulu-Natal,
Durban, South Africa

<sup>23</sup>University of Chicago, 5640 South Ellis Avenue, Chicago, Illinois 60637, USA

<sup>24</sup>Dunlap Institute for Astronomy and Astrophysics, University of Toronto,
50 St George St, Toronto, Ontario, M5S 3H4, Canada

<sup>25</sup>Department of Astronomy and Astrophysics, University of Toronto,
50 St George St, Toronto, Ontario, M5S 3H4, Canada

<sup>26</sup>School of Physics and Astronomy, University of Minnesota,
116 Church Street S.E. Minneapolis, Minnesota 55455, USA

<sup>27</sup>Department of Physics, Stanford University, 382 Via Pueblo Mall, Stanford, California 94305, USA

<sup>28</sup>High Energy Accelerator Research Organization (KEK), Tsukuba, Ibaraki 305-0801, Japan

<sup>29</sup>Canadian Institute for Advanced Research, CIFAR Program in Gravity and the Extreme Universe,
Toronto, Ontario, M5G 1Z8, Canada

<sup>30</sup>Service des Basses Températures, Commissariat à l'Energie Atomique, 38054 Grenoble, France

<sup>31</sup>Department of Astrophysical and Planetary Sciences, University of Colorado,
Boulder, Colorado 80309, USA

<sup>32</sup>Department of Physics, University of Illinois Urbana-Champaign,
1110 W. Green Street, Urbana, Illinois 61801, USA

<sup>33</sup>Astronomy Department, University of Illinois at Urbana-Champaign,
1002 W. Green Street, Urbana, Illinois 61801, USA

<sup>34</sup>Harvey Mudd College, 301 Platt Blvd., Claremont, California 91711, USA

<sup>35</sup>European Southern Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei München, Germany

<sup>36</sup>Department of Physics, University of California, One Shields Avenue, Davis, California 95616, USA

<sup>37</sup>Physics Division, Lawrence Berkeley National Laboratory, Berkeley, California 94720, USA

<sup>38</sup>Institut d'Astrophysique de Paris, 98 bis boulevard Arago, 75014 Paris, France

<sup>39</sup>Institute of Theoretical Astrophysics, University of Oslo, P.O.Box 1029 Blindern, N-0315 Oslo, Norway

<sup>40</sup>Department of Applied Mathematics and Theoretical Physics, University of Cambridge,
Cambridge, CB3 0WA, United Kingdom

<sup>41</sup>Department of Physics, University of Toronto, Toronto, Ontario, M5S 1A7, Canada

<sup>42</sup>Materials Sciences Division, Argonne National Laboratory,
9700S. Cass Avenue, Argonne, Illinois 60439, USA

<sup>43</sup>Physics Department, Center for Education and Research in Cosmology and Astrophysics,
Case Western Reserve University, Cleveland, Ohio 44106, USA

<sup>44</sup>Department of Physics, Yale University, P.O. Box 208120, New Haven, Connecticut 06520-8120, USA

<sup>45</sup>Liberal Arts Department, School of the Art Institute of Chicago,
112S Michigan Ave, Chicago, Illinois 60603, USA

<sup>46</sup>Physics Department, Brookhaven National Laboratory, Upton, New York 11973, USA

<sup>47</sup>Three-Speed Logic, Inc., Victoria, British Columbia, V8S 3Z5, Canada

<sup>48</sup>Department of Physics, University of California at San Diego, La Jolla, California 92093, USA

<sup>49</sup>Space Science and Engineering Division, Southwest Research Institute, San Antonio, Texas 78238, USA

<sup>50</sup>Department of Physics and Astronomy, Michigan State University,
567 Wilson Road, East Lansing, Michigan 48824, USA

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

(Received 18 November 2020; accepted 23 December 2020; published 26 January 2021)

We present a constraint on the tensor-to-scalar ratio, *r*, derived from measurements of cosmic microwave background (CMB) polarization *B*-modes with “delensing,” whereby the uncertainty on *r* contributed by the sample variance of the gravitational lensing *B*-modes is reduced by cross-correlating against a lensing *B*-mode template. This template is constructed by combining an estimate of the polarized CMB with a tracer of the projected large-scale structure. The large-scale-structure tracer used is a map of the cosmic infrared background derived from *Planck* satellite data, while the polarized CMB map comes from a combination of South Pole Telescope, BICEP/Keck, and *Planck* data. We expand the BICEP/Keck likelihood analysis framework to accept a lensing template and apply it to the BICEP/Keck dataset collected through 2014 using the same parametric foreground modeling as in the previous analysis. From simulations, we find that the uncertainty on *r* is reduced by ~10%, from $\sigma(r) = 0.024$ to 0.022, which can be compared with a ~26% reduction obtained when using a perfect lensing template or if there were zero lensing *B*-modes. Applying the

<footnote>* Corresponding author: W. L. K. Wu.
wlwu@slac.stanford.edu</footnote>
technique to the real data, the constraint on $r$ is improved from $r_{0.05} < 0.090$ to $r_{0.05} < 0.082$ (95% C.L.).
This is the first demonstration of improvement in an $r$ constraint through delensing.

DOI: 10.1103/PhysRevD.103.022004

# I. INTRODUCTION

Inflation describes a period of near-exponential expansion during the earliest moments of the Universe. The inflationary paradigm provides conceptual solutions to problems arising from the big bang description of the early Universe including the horizon problem and the flatness problem. Furthermore, inflationary models make testable predictions about perturbations away from perfect homogeneity and isotropy [1]. These predictions have been confirmed in observations of the cosmic microwave background (CMB) temperature and polarization anisotropies. They include the Gaussianity, phase-synchronicity, and near-scale-invariance of the scalar density fluctuations, and superhorizon correlation of the CMB anisotropies [2]. However, one prediction from inflation that has yet to be confirmed is the existence of a stochastic primordial gravitational wave (PGW) background.

PGWs are generically predicted in many inflationary models. Their amplitude is parametrized by $r$, the ratio of the amplitudes of the tensor and scalar perturbation spectra at a pivot scale ($k_* = 0.05$ Mpc$^{-1}$ in this work). If PGWs exist, they would imprint a specific divergence-free ($B$-mode) signature in the polarization of the CMB [3,4]. This makes CMB polarization a promising avenue in the search for PGWs.

However, PGWs are not the only source of $B$-modes. Thermal dust and synchrotron emission within our Galaxy produce polarized foreground patterns which contain $B$-modes [5,6]. Additionally, there is a source of $B$-modes, called the “lensing $B$-mode,” produced by gravitational lensing of the CMB [7]. If there were no inhomogeneities in the matter between us and the last scattering surface, then scalar perturbations from inflation would produce a purely curl-free ($E$-mode) CMB polarization pattern. However, during their propagation to us, the polarized CMB photons undergo small gravitational deflections by the forming large-scale structure along the line of sight. This produces a $B$-mode component which is small compared to the source $E$-modes, and which has already been detected by a number of experiments [8–14].

The BICEP/Keck experiments have deployed CMB polarization telescopes optimized for measurements at the “recombination bump” in the predicted PGW-generated $B$-mode spectrum (harmonic multipoles $\ell \sim 80$, or angular scales of $\sim 2$ deg). To separate out the Galactic dust and synchrotron components, which have different frequency spectral shapes than the blackbody emission of the CMB, BICEP/Keck observes in several frequency bands, and the $r$ analyses also incorporate maps at additional frequencies

from the WMAP and Planck satellites. The existing analysis pipeline takes all possible auto- and cross-spectra of the maps at different frequencies and compares these against a parametric model of CMB and foregrounds [12,14] to set constraints on $r$ which are close to optimal given the available data. Alternative approaches involving “cleaning coefficient” subtraction of a dust template map (as measured at a higher frequency) would in general be less powerful [e.g., [15]].

In contrast to the foregrounds, the lensing component has the same frequency spectral shape as the PGW component and thus cannot be constrained using multifrequency observations. Given an estimate of the projected gravitational potential responsible for CMB lensing and the observed CMB $E$-mode pattern, one can estimate the $B$-modes which have been produced by the lensing effect. Subtracting these from the observed $B$-modes has been demonstrated to reduce $B$-mode power in several recent works [16–20]. However, none of these works have demonstrated a reduction in the $B$-mode measurement uncertainties at large angular scales—a necessary step to achieve improved constraints on PGWs. This subtraction process is usually referred to as “delensing.” But in this work, we take a different approach and therefore broaden the meaning of delensing to include any process which reduces the effective lensing sample variance in the $B$-mode measurements. Specifically, we extend the BICEP/Keck analysis pipeline to accept an estimate of the lensing $B$-modes as a “lensing template”—an additional pseudo-frequency band against which cross-spectra are taken. This (optimally) reduces the effective sample variance of the lensing $B$-mode component and hence reduces the uncertainty of the PGW contribution.

The lensing potential $\phi$ can be computed using higher-order statistics of the CMB pattern itself [21]. However, since the lensing potential is a weighted integral of the mass distribution along the line of sight between us and the last scattering surface, we may also approximate it by other tracers of this mass distribution. At the noise levels of current CMB observations, it turns out to be better to use a cosmic infrared background (CIB) [22,23] map rather than one of the available CMB lensing reconstructions [18,24] directly. To use an alternate tracer of $\phi$, we need to know the degree of correlation between it and the true CMB lensing potential—if this were misestimated it could potentially lead to a false detection of PGW. This correlation may be found empirically from the cross-correlation of the tracer with a reconstruction of the CMB lensing potential. In this paper, we use a CIB map from Planck
generated using the generalized needlet internal linear combination (GNILC) component separation algorithm [25] as the $\phi$ tracer and estimate its correlation with the lensing potential using a Planck minimum-variance lensing map [26].

To estimate the lensing template, in addition to the tracer of the lensing potential, one also needs the best available estimate of the observed CMB polarization pattern. Since the lensing operation mixes modes over a wide range of angular scales, the inclusion of small-scale E-modes is important for precise estimation of the lensing B-modes at the angular scales of interest ($\ell \sim 80$). Therefore, we use arcminute-resolution maps from the South Pole Telescope (SPT) second-generation camera SPTpol, augmenting these with polarization measurements from BICEP/Keck and Planck.

In this paper, we add the CIB-derived lensing template to the previous “BK14” analysis [12] which utilizes data from BICEP/Keck through the 2014 observation season. With the addition of the lensing template, we demonstrate a ~10% reduction in the uncertainty on $r$ for the BK14 dataset, to be compared with a ~26% reduction in the uncertainty on $r$ when using a perfect lensing template or if there were zero lensing B-modes. This shows that the lensing sample variance is a subdominant fraction of the uncertainty on $r$ for BK14. However, it will be an increasingly limiting factor going forward. Therefore, this analysis serves as a proof of principle, and a first step towards future analyses where delensing will more significantly improve $\sigma(r)$.

This paper is organized as follows: In Sec. II, we describe the construction of the lensing template and the extension to the BICEP/Keck pipeline to include the lensing template. In Sec. III, we describe the data and simulation sets of the CMB maps, how we combine the $Q/U$ maps

from SPTpol, BICEP/Keck, and Planck, and the data and simulations of the $\phi$ tracer. We validate our simulations and pipeline in Sec. IV and test for systematics in Sec. V. We present our results in Sec. VI and conclude in Sec. VII.

## II. METHOD

In this section, we describe new elements added to the BICEP/Keck analysis framework to incorporate information on the lensing B-modes in the BICEP/Keck patch, with the aim of reducing the effective uncertainty of the observed B-modes, and thereby reducing the uncertainty on $r$. We illustrate the incorporation of the lensing template into the BICEP/Keck likelihood analysis framework schematically in Fig. 1. There are two main areas of new development: (1) constructing a lensing template and (2) extending the BICEP/Keck pipeline to include the lensing template. We will describe each aspect in the following subsections.

### A. Constructing the lensing template

The key element to constraining the lensing B-modes in the BICEP/Keck patch is making an estimate of these modes. To do this, we use two inputs: (1) a tracer of the CMB lensing potential $\phi$ from large-scale structure observations and (2) observed $Q/U$ polarization maps. We construct the lensing template using an “undeflect-and-difference” method in which we undeflect the observed $Q/U$ maps using the $\phi$ tracer and subtract the undeflected maps from the input.

Formally, we take the lensed, polarized CMB fields $\tilde{X}$, which are related to the un lensed CMB fields $X$ by

$$ \tilde{X}_{\pm}(\hat{n}) = X_{\pm}(\hat{n} + \nabla\phi(\hat{n})), \quad (1) $$

<figure><img src="image_3.png" /><figcaption>FIG. 1. Schematic of the analysis flow in this work. The rectangular blocks denote input maps; the blocks with rounded corners denote operations on maps. The teal-colored region highlights the inputs to, and processes involved in generating, the lensing template. The input maps include the SPTpol, BICEP/Keck, and Planck $Q/U$ maps, the Planck GNILC CIB map, and the Planck minimum-variance (MV) reconstruction of $\phi$. The Planck MV $\phi$ is in parentheses because, instead of using it as a $\phi$ tracer, we use it to filter and normalize the CIB map and for generating simulations. The unshaded region denotes the standard BICEP/Keck $r$ analysis, where auto- and cross-spectra of multifrequency maps from BICEP/Keck, Planck, and WMAP form the input data for computing likelihoods to extract parameter constraints. The lensing template is injected into the standard analysis as an additional pseudofrequency band.</figcaption></figure>
where $X_{\pm} = (Q \pm iU)$ and $\nabla\phi$ denotes the deflection field [27]. We undo the deflection by remapping the $\tilde{Q}$ and $\tilde{U}$ polarization fields by $-\nabla\phi$, evaluated at the lensed positions $\hat{n}'$ (not delensed positions $\hat{n}$),

$$X_{\pm}^{d}(\hat{n}) = \tilde{X}_{\pm}(\hat{n}'), \qquad (2)$$

where $\hat{n}' = \hat{n} - \nabla\phi(\hat{n}')$ and $X^d$ denotes the undeflected field. Therefore,

$$X_{\pm}(\hat{n}) = X_{\pm}^{d}(\hat{n}) = \tilde{X}_{\pm}(\hat{n} - \nabla\phi(\hat{n} - \nabla\phi(\hat{n} - \dots))), \quad (3)$$

where the last expression is used for the practical implementation, and ... denotes the recursion that locates the position to which the value at $\hat{n}$ was deflected from the un lensed plane.

Specifically, the undeflection is implemented by first computing the amount of deflection at the lensed positions, denoted by $(dx, dy)$, on the delensed map pixel grid $(x, y)$. To do that we first evaluate $\nabla\phi$ at $(x, y)$ to get $(dx', dy')$, and then we evaluate $\nabla\phi$ at $(x - dx', y - dy')$, and so on. We find that the solution converges after 1 recursion, which means that with the notation given, $(dx, dy)$ is $\nabla\phi$ at $(x - dx', y - dy')$. The evaluation of $\nabla\phi$ at any grid point is done by interpolating $\nabla\phi$ values in HEALPix format using first-order Taylor expansion. We then remap the $\tilde{Q}/\tilde{U}$ map pixels at $(x - dx, y - dy)$ to $(x, y)$ via cubic interpolation. We note that by evaluating the deflection field at the lensed positions, we do not incur the small $O(\nabla\phi \cdot \nabla)\nabla\phi$ error found in similar algorithms that evaluate $\nabla\phi$ at the delensed positions [17,28,29].

The lensing templates $Q^i/U^i$ are then derived by subtracting the obtained undeflected map from the observed (lensed) one,

$$Q^i(\hat{n}) = \tilde{Q}(\hat{n}) - Q^d(\hat{n}) \qquad (4)$$

$$U^i(\hat{n}) = \tilde{U}(\hat{n}) - U^d(\hat{n}). \qquad (5)$$

We test the algorithm on noiseless lensed simulations using the $\phi$ maps which were used to lens them. The correlation of the resulting lensing B-mode template with the difference of the lensed and un lensed input skies is $\geq 95\%$ for the angular scales used in this analysis. This is sufficiently accurate at the current noise levels. In other work, lensing templates have also been constructed after transforming to harmonic space, converting to $E/B$, and lensing by a $\phi$ tracer using expressions derived from the first-order Taylor expansion of Eq. (1) [16,18,19]. At the noise levels of the current analysis, the two approaches perform similarly in constraining the lensing B-mode contribution to the observed B-modes. We discuss in more detail the differences of the two approaches in the Appendix.

Since the undeflect-and-difference operation corresponds to an all-with-all mixing in Fourier space, to obtain

the lowest possible lensing template noise in the $\ell$ range of interest we first Wiener filter [e.g., [29]] the $Q/U$ and $\phi$-tracer maps. We filter the $Q/U$ maps by a 2D Wiener filter in Fourier space,

$$\tilde{Q}(\ell) \rightarrow \frac{C_{\ell}^{EE}}{C_{\ell}^{EE} + N_{\ell}^{EE}} \tilde{Q}(\ell), \qquad (6)$$

$$\tilde{U}(\ell) \rightarrow \frac{C_{\ell}^{EE}}{C_{\ell}^{EE} + N_{\ell}^{EE}} \tilde{U}(\ell), \qquad (7)$$

to account for anisotropic noise and mode-loss due to filtering. $C_{\ell}^{EE}$ and $N_{\ell}^{EE}$ are 2D power spectra of the E-mode signal and noise components, constructed from a weighted combination of $Q/U$ maps from the three experiments SPTpol, BICEP/Keck, and Planck. We describe the procedure to combine the $Q/U$ maps and the details of the Wiener filter in Sec. III A 3.

We filter and normalize the $\phi$ tracer in spherical harmonic space according to

$$\phi_{LM}^{T} = \left( \frac{C_{L}^{T\phi'}}{C_{L}^{TT}} \right) T_{LM}, \qquad (8)$$

where $T$ denotes the tracer and $\phi'$ is an unbiased, but noisy map of the true CMB lensing potential [22,23]. To see that this weighting is a joint normalization and Wiener filter, we write $T_{LM} = g_L \phi_{LM} + n_{LM}$, where $g_L$ is the relative normalization factor (and unit conversion), $\phi$ is the true (noiseless) lensing potential, and $n$ is the effective noise in the tracer pattern (the part which does not correlate with $\phi$) with power spectrum $N_L^{TT}$. Expanding and taking the expectation value, we get

$$\frac{C_L^{T\phi'}}{C_L^{TT}} = \frac{g_L C_L^{\phi\phi}}{g_L^2 C_L^{\phi\phi} + N_L^{TT}} = g_L^{-1} \frac{C_L^{\phi\phi}}{C_L^{\phi\phi} + N_L^{TT}/g_L^2}, \qquad (9)$$

fulfilling its role of normalization and filtering. In this paper the tracer $T$ is a CIB map from Planck and $\phi'$ is a lensing reconstruction map also from Planck—see Sec. III B below for further details.

With the $Q/U$ lensing templates constructed, we then take them as an additional pseudofrequency band for input into the existing BICEP/Keck analysis.

## B. Adding the lensing template to the existing analysis framework

The development of the existing BICEP/Keck r analysis framework has been described in a series of papers [12,14,15,30]. Briefly, we take all possible auto- and cross power spectra between the available frequency bands, and then compare the resulting set of bandpowers to their expectation values under a parametric model of lensed- $\Lambda$CDM + dust + synchrotron + r using an expansion of
the Hamimeche-Lewis likelihood approximation [31]. It is a straightforward extension to this framework to include the lensing template as an additional pseudofrequency band. To do this we require reliable simulations of the signal and noise content of the lensing template so that we can (1) debias its autospectrum, (2) determine the expectation values of the auto- and cross-spectra involving the lensing template, and (3) determine the variance of these band-powers, and their covariance with other bandpowers. These simulations are described in Sec. III below. Here we describe a few complications with respect to the normal procedure which arise in the steps above.

The lensing template is formed from two kinds of input maps (the $Q/U$ maps and the $\phi$ tracer) which both contain relevant amounts of noise. The Planck CIB map has very high signal-to-instrumental-noise. However, the integrated dust emission from star-forming galaxies back to the last scattering surface weights differently over redshift than the deflection of CMB photons, and these galaxies do not perfectly trace the underlying mass density field. This means that the CIB only partially correlates with the true lensing potential. For the purposes of this paper, the $\phi$ tracer signal is the portion of the CIB that is correlated with the true lensing potential $\phi$; the $\phi$ tracer noise corresponds to the uncorrelated portion. We detail our $\phi$ tracer simulations in Sec. III B 2.

We remove the noise bias of the lensing template autospectrum by subtracting the noise autospectrum estimated from simulations. Schematically, the lensing template B-mode autospectrum is

$$ \begin{aligned} \langle L_B^2 \rangle &= \langle \langle (s_{QU} + n_{QU}) * (s_\phi + n_\phi) \rangle^2 \rangle \\ &= \langle (s_{QU} * s_\phi)^2 \rangle + \langle (s_{QU} * n_\phi)^2 \rangle \\ &\quad + \langle (n_{QU} * s_\phi)^2 \rangle + \langle (n_{QU} * n_\phi)^2 \rangle, \end{aligned} \quad (10) $$

where $s_X$ and $n_X$ denote the signal and noise components of field $X \in [QU, \phi]$ and $*$ denotes the following steps: undeflect-and-difference, Fourier transform, and convert from $Q/U$ to B-modes. In writing the second line, we have assumed all the cross terms have zero expectation value. We estimate the noise autospectrum from simulations as

$$ \langle (n_{QU} * (s_{\phi} + n_{\phi}))^2 \rangle + \langle (s_{QU} * n_{\phi})^2 \rangle, \quad (11) $$

averaged over all simulation realizations and subtract it from Eq. (10). The $Q/U$ and $\phi$ input signal and noise maps are Wiener filtered in the same way as the data maps (and the simulation signal + noise maps). Empirically, when adding this inferred noise bias to the mean of the signal-only simulation spectra ($Q/U$ signal undeflected with $\phi$ tracer signal) one obtains the mean of the signal + noise simulation spectra ($Q/U$ signal + noise undeflected with $\phi$ tracer signal + noise) to high fractional precision.

In the BICEP/Keck standard procedure, the filter/beam suppression of the bandpower values is computed using sets of maps which each contains power at only a single multipole $l$ passed through the “observing matrix” as described in Sec. VI. C of [30]. However, since the lensing template is derived in a very different manner to the standard BICEP/Keck maps, the usual observing matrix is not applicable, and we fall back to a simulation-based approach. We rescale both the data and simulation lensing template auto- and cross-spectra by the ratio of the input lensing spectrum $C_{\ell}^{BB}$ to the average of the signal-only simulation bandpowers. This step overrides the normalization part of Eq. (9) applied to the $\phi$ tracer. However, accurate knowledge of the degree of correlation between the lensing tracer and the true lensing potential is still required to avoid bias on $r$ (see Sec. IV).

In the standard BICEP/Keck procedure, the bandpower covariance matrix is constructed by taking the auto- and cross-spectra of the signal and noise components of the simulations as described in Appendix H of [14]. Since the lensing template is formed from two maps which both have signal and noise components we expand the usual procedure to form additional cross-spectra and combine the results appropriately.

With this extended analysis framework, we can now incorporate lensing templates constructed using simulations and data to the BICEP/Keck likelihood and constrain the model parameters.

### III. DATA AND SIMULATIONS

The BICEP/Keck analysis pipeline relies on signal-only, noise-only, and signal + noise simulations, the construction of which is described in Sec. V of [30]. We reuse the data maps and simulations including Gaussian realizations of Galactic dust from the BK14 analysis unchanged. The data maps include the WMAP and Planck bands with BICEP/Keck filtering applied (as described in Sec. II. A of [15]). To add the lensing template as an additional pseudofrequency band, we need data maps and corresponding simulations of it. Since the lensing template is constructed from $Q/U$ CMB maps and a CIB $\phi$ tracer, we in turn need data maps and simulations of both of these. As a prestep, we combine the SPTpol, BICEP/Keck, and Planck $Q/U$ maps to generate a synthetic map which has the best possible signal-to-noise at all points in the 2D Fourier plane. Figure 1 gives a schematic view of the steps involved in generating the various maps.

#### A. $Q/U$ CMB maps

Below we describe the data processing of the SPTpol, BICEP/Keck, and Planck $Q/U$ maps that are relevant in the construction of the combined $Q/U$ maps and their Wiener filter. The combined, Wiener filtered $Q/U$ maps are the
inputs to the undeflect-and-difference step which is used to
construct the lensing template.

1. Data CMB maps

SPTpol maps: We use SPTpol maps made specifically for this analysis using 150 GHz observations taken between 2013 and 2015 by the SPTpol camera [32] on the South Pole Telescope [33]. The SPTpol 500 deg² survey field is centered at RA 0h and Dec. -57.5°, matching the BICEP/Keck field. The polarization map depth is ~10 µK arcmin in the multipole range of 300 $\lesssim \ell \lesssim 2000$. The time stream processing is identical to that in [34], except for the polynomial-filter order and the low-pass filter. We fit and subtract a third-order/sixth-order polynomial from the time stream of each detector over the RA extent of the lead-trail/full-field observations. We choose the low-pass filter based on the pixel size. This set of SPTpol maps is binned into 5 arcminute-sized pixels, which are a ×3 resolution superset of the BICEP/Keck map pixels. To reduce aliasing given the pixel size, we apply a low-pass filter to the time stream that corresponds to $\ell \sim 1900$. The polarization maps, in addition to the calibration factors included through calibrating the temperature map against Planck, have an extra polarization calibration factor ($P_{cal}$) applied. The polarization calibration factor is taken from [34] and is obtained by forming a cross-spectrum between the SPTpol E-mode map and an E-mode map from Planck.<fnref n="1" /> We discuss impacts on $r$ from biases in $P_{cal}$ in Sec. V.

BICEP/Keck maps: We use the BICEP2/Keck 150 GHz band $Q/U$ maps from BK14. These have noise of ~3 µK arcmin over an effective area of 395 deg² centered at RA 0h, Dec. -57.5°. The BICEP2 and Keck Array telescopes have ~30 arcminute resolution at 150 GHz. This limits the highest angular multipole to which they are sensitive to $\ell$ of hundreds. As described in Secs. III and IV of [30] the construction of the maps involves time-stream filtering. Specifically, a third-order polynomial was subtracted from the time streams of each detector over each scan. Across the ~30° scan throw on the sky, this approximately corresponds to removing $\ell_x < 20$ modes. These maps are binned in 0.25° rectangular pixels in RA and Dec, and calibrated by forming cross-spectra with the Planck temperature map.

Planck maps: We use the 143 GHz $Q/U$ “full mission” maps from Planck public release 2 as the input to the combined three-experiment $Q/U$ maps. We convert the map to $a_{\ell m}$, apply an antialiasing filter by low-pass filtering at $\ell = 2100$, render a $N_{\text{side}} = 2048$ HEALPix map, and interpolate to the same 5-arcminute pixel grid as used for the SPTpol maps.

2. Simulated CMB maps

We reuse the BK14 simulated maps unchanged. We thus need to make corresponding simulations for the SPTpol and Planck maps. The BICEP/Keck CMB sky realizations have remained the same since originally described in Sec. V of [30]. These are the unlensed $a_{\ell m}$, Gaussian realizations of $\phi$ given the input cosmology, and lensed $a_{\ell m}$ generated using LensPix [35]. The BICEP/Keck simulations were originally generated with a maximum $\ell$ of 1536 which is adequate given the beam sizes of the telescopes. To match more closely the pixel scale of the SPTpol and Planck maps in this analysis, we generate additional higher-$\ell a_{\ell m}$, graft these onto the existing unlensed values, pass through LensPix, and graft the output onto the existing lensed values. Since lensing to some degree mixes angular scales this is clearly only approximately correct, but we note that the amount of lensing B-modes below $\ell$ of 350 sourced from E-modes between $\ell$ of 1536 and 2100 (the pixel-scale) is negligible [see e.g., Fig. 2 of [22]]. We refer to this set of input lensed and unlensed $a_{\ell m}$ as the extended set and the original set as the standard set.

We generate SPTpol simulations for this analysis using the extended set of $a_{\ell m}$. In a procedure similar to that used to generate the existing BICEP/Keck simulations, we multiply the input $a_{\ell m}$ by the instrument beam, “mock-observe” the skies by creating time-stream samples given the pointing information of each detector, apply the same time-stream level filters as applied to data, and bin to maps in the pixelization used for the real data. Corresponding noise realizations are generated by the standard method used in both BICEP/Keck and SPT analyses—differencing combinations of halves of data maps, where the halves are defined so that the weights of each half are close to equal.

We generate simulated maps for Planck 143 GHz by first taking the $a_{\ell m}$ from the extended set and multiplying them by the Planck 143 GHz beam. We then low-pass filter and process as for the real Planck map. Corresponding noise realizations are taken from the Planck FFP8 simulations and processed identically to the real Planck map. We generate 499 realizations of signal and noise skies for each experiment as is the BICEP/Keck standard.

3. Combining and filtering the Q/U maps
from SPTpol, BICEP/Keck, and Planck

A factor that impacts the delensing efficiency (the recovery of the lensing B-modes) is the per-mode noise of the input $Q/U$ maps. The lower the noise per mode, the better the lensing templates trace the true lensing B-modes. The lensing B-modes at multipole $\ell$ are mostly sourced by E-modes from a range of multipoles slightly higher in $\ell$ (smaller in angular scale) [see e.g., Fig. 2 of [22]]. BICEP/Keck does not image these smaller-scale E-modes very well because of its large beam size. Therefore, it is advantageous to combine with polarization measurements from other,

<footnote n="1">¹Planck COMMANDER maps: COM CMB IQU-commander 1024 R2.02 full.</footnote>
<figure><img src="image_4.png" /><figcaption>FIG. 2. Simulated 2D E-mode signal power spectra of BICEP/Keck, SPTpol, and Planck. The axis scales for the BICEP/Keck Fourier plane are zoomed in compared with the rest of the panels to focus on the modes accessible by BICEP/Keck's small apertures. The color stretch in all four panels is identical. For BICEP/Keck and SPTpol because the observations are being made at the South Pole with scans along the azimuth direction, scanwise filtering leads to modes along the $\ell_y$ axis being suppressed. These filtered modes along the $\ell_x$ axis can be partially filled in using measurements from Planck. To generate the combined, Wiener filtered 2D E-mode signal power spectra on the rightmost panel, the three sets of modes to the left are corrected for beam and filtering, combined using inverse-noise weighting, and Wiener filtered to suppress modes which remain noisy in the combined set (as described in Sec. III A 3). We see that some modes remain unavailable for lensing template construction at $|\ell_x| < 100$ and $|\ell_y| > 500$.</figcaption></figure>

higher-resolution experiments such as SPTpol and Planck to increase the signal-to-noise ratio of the input $Q/U$ maps, and thus the E-modes.

We combine the three maps in Fourier space. We divide from the 2D mode sets of the three experiments their respective 2D transfer functions, taken as the square-root of the mean of the 2D E-mode power spectra of the signal-only simulations divided by the mean of the corresponding spectra of the (unfiltered) input maps. We also divide the 2D noise power spectra by the same ratio (without the square root). We then combine the SPTpol, BICEP/Keck, and Planck $Q/U$ modes using an inverse-variance weighting taken from the mean of the 2D noise power spectra. Specifically, the combined $Q/U$ mode sets are

$$ \mathcal{X}(\boldsymbol{\ell}) = \sum_i w_i(\boldsymbol{\ell}) \mathcal{X}_i(\boldsymbol{\ell}), \quad (12) $$

where $\mathcal{X} \in [Q, U]$, $i \in [\text{SPTpol}, \text{BICEP/Keck}, \text{Planck}]$, and $w_i$ denotes the weight

$$ w_i(\boldsymbol{\ell}) = \frac{N_i^{-1}(\boldsymbol{\ell})}{\sum_i N_i^{-1}(\boldsymbol{\ell})}. \quad (13) $$

Here, $N_i(\boldsymbol{\ell})$ denotes the mean of the transfer-function-divided 2D angular power spectra of the E-mode noise realizations from experiment i. We additionally impose $\ell_x$ cuts by artificially increasing the noise below some $\ell_x$ to remove modes that are empirically found to be unrecoverable due to the scan-wise time-stream filtering. We set $\ell_x$ to 25 for BICEP/Keck and $\ell_x$ to 50 for SPTpol.

Before passing the combined $Q/U$ map to the lensing template construction step, we apply a Wiener filter as described in Sec. II A above. The $C_{\ell}^{EE}$ in Eq. (6) is the 2D

input E-mode power spectrum and the $N_{\ell}^{EE}$ of the combined modes is

$$ N_{\ell}^{EE} = \sum_{i} w_{i}^{2}(\boldsymbol{\ell}) N_{i}(\boldsymbol{\ell}), \quad (14) $$

with $w_i(\boldsymbol{\ell})$ given by Eq. (13).

In the above, we transform to Fourier space and back again, and hence need to choose an apodization mask. Due to the small instantaneous field of view of the SPTpol camera as compared to the size of the observation region, the integration time map (inverse noise variance map) is a near uniform rectangular box tapering to zero over a few degrees at the edges. In contrast, the BICEP/Keck integration time map has no uniform central region and tapers smoothly and continuously from a peak in the middle (see for instance Fig. 1 of [15]), with nonzero coverage extending well outside the SPTpol region. (Planck observes the full sky and has close to uniform coverage across the sky region in question.) To perform the map combination we need to pick a single apodization function for all three input maps. We choose to use the one built from the SPTpol integration time map, with a cosine taper with a radius of 1 deg. This is because SPTpol is the experiment with the most restrictive sky coverage but the best mode coverage. This means that the resulting lensing template does not cover the full BICEP/Keck sky region. In addition, because of the chosen spatial weighting of pixels, we introduce suboptimality in the combination.

Figure 2 illustrates the process. The left three panels show the 2D E-mode signal power spectra for the three experiments. We see the $\Lambda$CDM E-mode spectrum rolled-off by the beam window function of each telescope. Because of their scan strategies and the applied scanwise
<figure><img src="image_5.png" /><figcaption>FIG. 3. The top two panels show the experiment-combined and Wiener filtered $Q/U$ maps. The middle two panels show the $x$ and $y$ derivatives of the normalized and Wiener filtered Planck CIB map. Signal and noise are approximately equal in these maps. Due to the foreshortening effect the RA deflections are larger and increase towards more negative Dec. The $Q/U$ maps in the top panels are undeflected by the angles shown in the middle panels and differenced with the initial maps to form the lensing template $Q/U$ maps shown in the bottom panels.</figcaption></figure>

filtering, BICEP/Keck and SPTpol have filtered out the modes along the $\ell_y$ axis; while Planck has isotropic mode coverage. The right panel shows the combined mode set after the final Wiener filter step, so only modes measured with good signal-to-noise are retained. At $\ell > 500$ Planck does not have good per-mode signal-to-noise so the modes along the $\ell_y$ axis beyond this multipole cannot be filled in.

We next proceed to inverse-Fourier transform the combined and Wiener filtered $Q/U$ modes back to image space where they are ready to be undeflected by the gradient of the $\phi$ tracer. The $Q/U$ maps at this stage are shown in the top panels of Fig. 3.

## B. CIB map

With the combined and filtered $Q/U$ maps in hand, we next need a $\phi$ tracer map. In the following, we describe the characteristics of the CIB map used in this analysis, and how we generate simulations of it in the BICEP/Keck patch.

## 1. CIB data

It is possible to reconstruct the lensing potential field $\phi$ from the CMB temperature and polarization patterns [18], and in the future this will become the best $\phi$ estimate for delensing [36]. However, at the currently available noise levels the most effective available $\phi$ tracer is the CIB, even though it is only partially correlated with $\phi$ [22]. Specifically, we use the 545 GHz CIB map from Planck generated using the GNILC algorithm [25].<fnref n="2" /> We also considered using the CIB maps generated by [37] and will discuss that later in this section. To determine the degree to which the GNILC CIB map is correlated with $\phi$, we use the Planck 2015 minimum-variance lensing reconstruction map<fnref n="3" /> [26] and make the assumption that this is an unbiased (although noisy) representation of the true $\phi$ pattern.

<footnote n="2">²CIB map: COM_CompMap_CIB-GNILC-F545_2048_R2.00.</footnote>

<footnote n="3">³Planck lensing map: COM_CompMap_Lensing_2048_R2.00.</footnote>
<figure><img src="image_6.png" /></figure>

FIG. 4. The light red regions denote the *Planck* lensing mask, used for computing the “full sky” average of GNILC CIB and $\phi'$ map cross-correlation. The yellow regions are the eight patches with similar size and unpolarized dust amplitudes as the BICEP/Keck patch. These patches used for measuring the mean and scatter of the CIB autospectra and CIB $\times \phi'$, which are used as inputs to simulating CIB and filtering the CIB map. The overlaps between the yellow patches are small and apodization is applied when calculating the auto- and cross-spectra. The BICEP/Keck patch is shaded in blue. The background is the *Planck* dust intensity map.

We refer the reader to [38] for a detailed discussion of the *Planck* CIB map. Briefly, the GNILC component-separation technique [39] disentangles different components of emission using both frequency and spatial (angular-scale dependence) information. In this case, the GNILC algorithm was applied to *Planck* data to disentangle Galactic dust emission and CIB anisotropies. Even though both components share similar frequency spectral signatures, they have distinct angular power spectra. Thus by using priors on the angular power spectra of the CIB, Galactic dust, the CMB, and the instrumental noise, these components can be (partially) separated. We note that the algorithm was developed mainly for extracting Galactic dust, and regions with different levels of Galactic dust can be expected to have different efficiencies of CIB recovery [e.g., [37,40]]. Therefore, in the following, we quantify the GNILC CIB map correlation with the *Planck* estimate of $\phi$ empirically in selected parts of the sky.

To select patches for estimating the CIB-$\phi$ correlations, we measure the mean amplitude in a *Planck* dust temperature map<fnref n="4" /> of ~500 deg²-sized circles throughout the sky. Amongst the eight selected patches, as shown in Fig. 4, the ratios of the mean amplitudes in the patches vs that in the BICEP/Keck patch range from 0.6 to 1.7. These are thus similar to the BICEP/Keck region in terms of their unpolarized dust intensities.

In these patches, we compute the autospectra and cross-spectra using PolSpice [41] [42]. We show the correlations

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

FIG. 5. The binned correlation factor $\rho_L = C_L^{\phi'} / \sqrt{C_L^{\text{II}} C_L^{\phi\phi}}$ for the eight patches, the full sky, and the BICEP/Keck patch. “Full sky” corresponds to the overlap area between the *Planck* lensing mask and the GNILC CIB map. The $\rho_L$ in the BICEP/Keck patch is consistent with those measured across the eight patches. The yellow band denoted by “eight patches” is the mean and standard deviation of $\rho_L$ across the eight patches. The error bars for the red and blue points are computed by taking the standard error of $\rho_L$ within each $\Delta L = 100$ bin. The red and blue points are shifted for clarity.

$\rho_L$, defined as $C_L^{\phi'} / \sqrt{C_L^{\text{II}} C_L^{\phi\phi}}$, for a few different regions of the sky in Fig. 5. Here *I* denotes the CIB map, $\phi'$ denotes the *Planck* lensing estimate, and $C_L^{\phi\phi}$ is the theory spectrum from the fiducial cosmology used in [43]. Comparing the correlations of the CIB map and the lensing map in the selected patches with that from the full overlap between the two maps (labeled “full sky”), we observe that the correlations within the patches are higher than the correlation in the larger region that includes lower Galactic latitudes and hence higher dust levels. The full-overlap region correlation is ~62% for *L* between 150 and 550, whereas the mean correlation in the patches is ~69% over the same *L* range. Figure 5 also shows the cross-correlation in the BICEP/Keck patch, which appears to be consistent with the eight circular patches.

As a cross-check, we compare within the BICEP/Keck patch the cross-spectrum of the GNILC CIB map and the *Planck* lensing map against the cross-spectrum of a CIB map produced by [37] and the *Planck* lensing map. This CIB map has been cleaned using neutral hydrogen (HI) as a Galactic foreground tracer, with an HI column density threshold of $2.5 \times 10^{20} \text{ cm}^{-2}$. We find the lensing correlation in the two CIB maps to be consistent with each other, thus providing additional evidence that in the BICEP/Keck map region, the GNILC CIB map does not show the reduced correlation which is expected, and seen, in regions closer to the Galactic plane.

The filter and normalization of the $\phi$ tracer is given in general form in Eq. (8). In this case, we take it as the average over the eight patches of the cross-spectra of the CIB and the lensing map divided by the average of the CIB autospectra,

<footnote n="4">⁴Thermal dust emission map: ThermalDust-commander_2048_R2.00/index.html.</footnote>
$$ \phi_{LM}^{\text{CIB}} = \left( \frac{\langle C_L^{\text{I}\phi'} \rangle_{\text{patches}}}{\langle C_L^{\text{II}} \rangle_{\text{patches}}} \right) I_{LM}. \qquad (15) $$

To further prevent Galactic dust contamination, we additionally impose a $L_{\min} = 100$ cut. This filter and normalization is applied to the real data as well as the simulated CIB realizations which are described in the next section. We render the normalized and Wiener filtered CIB $a_{\ell m}$ and its associated gradients to HEALPix maps of $N_{\text{side}} = 512$, and then interpolate and convert the gradient maps to derivatives with respect to our pixel grid. The derivatives are shown in the middle panels of Fig. 3.

## 2. CIB simulations

We use CIB simulations to estimate the expected level of lensing B-modes in the lensing template, to form the bandpower covariance in the likelihood analysis, and as inputs to null tests.

We generate CIB simulations based on the input Gaussian $\phi$ fields of the BICEP/Keck simulation set described in Sec. III A 2. To convert the $\phi$ fields to CIB fields, we use the autospectrum of the CIB, $C_L^{\text{II}}$, and the cross-spectrum of the CIB and the Planck lensing estimate $\phi'$, $C_L^{\text{I}\phi'}$. We construct each CIB field by rescaling each input $\phi$ field so that the cross-spectrum of the rescaled field with the input $\phi$ is $C_L^{\text{I}\phi'}$. We then add to the rescaled $\phi$ field Gaussian noise so that its autospectrum is $C_L^{\text{II}}$. Formally, we construct the signal part of the CIB simulations, $I_{LM}^{\text{S}}$, as

$$ I_{LM}^{\text{S}} = \frac{C_L^{\text{I}\phi'}}{C_L^{\phi\phi}} \phi_{LM}, \qquad (16) $$

where $\phi_{LM}$ are the spherical harmonic coefficients of the input $\phi$ fields. We construct the noise part of the CIB simulations, $I_{LM}^{\text{N}}$, by generating Gaussian random fields with power spectrum described by $C_L^{\text{II}} - (C_L^{\text{I}\phi'})^2/C_L^{\phi\phi}$. The total CIB field is the sum of the two terms $I_{LM} = I_{LM}^{\text{S}} + I_{LM}^{\text{N}}$.

We have 499 realizations of $\phi_{LM}$. For each $\phi_{LM}$, we form $I_{LM}$ as described in the previous paragraph with $C_L^{\phi\phi}$ from the input theory<fnref n="5" /> and $C_L^{\text{II}}$ and $C_L^{\text{I}\phi'}$ sampled from the measured mean and covariance of $C_L^{\text{II}}$ and $C_L^{\text{I}\phi'}$ from the eight patches selected in Sec. III B 1. In the limit of many realizations, the simulated $I_{LM}$ will have the same covariance structure in $C_L^{\text{II}}$ and $C_L^{\text{I}\phi'}$ as measured from the eight patches. The advantage of sampling $C_L^{\text{II}}$ and $C_L^{\text{I}\phi'}$ as opposed to using the measured mean from the eight patches is that

the potential patch-to-patch variation of the CIB autospectrum, and the cross-spectrum between CIB and $\phi'$, is built into the simulations. Therefore, the uncertainties in the CIB measurements are propagated to the uncertainty in the r measurement.

At this point we use the method described in Sec. II A to undeflect the combined $Q/U$ data and simulation maps with the data and simulation CIB maps to form the real and simulated lensing templates. The lensing templates from the real data are shown in the bottom panels of Fig. 3.

We have now laid out the lensing template construction, the extension of the BICEP/Keck analysis framework, and the input simulations and data used in this paper. The next steps include demonstrating the robustness of these extensions to potential biases and misestimations of inputs.

## IV. PIPELINE AND SIMULATION VALIDATION

In this section, we demonstrate the robustness of the pipeline in the limit of perfect delensing, quantify the level of bias to our inference of r given potential misestimations in the inputs to our simulations, and estimate the impact on $\sigma(r)$ given variations in the simulation setup. To do that we use the set of lensed-$\Lambda$CDM + dust + noise simulations ($r=0$) from the BK14 paper. We run maximum-likelihood searches of the baseline lensed-$\Lambda$CDM + dust + synchrotron + r model as described in Appendix E. 3 of the BK14 paper, in this case adding a lensing template.

### A. r recovery with perfect delensing

To validate the addition of the lensing template as a pseudoband in the BICEP/Keck analysis framework, we run maximum-likelihood searches in two configurations—unlensed input CMB skies without lensing templates and lensed input CMB skies with perfect lensing templates. The perfect lensing templates are constructed by differencing the filtered, noiseless, lensed and unlensed $Q/U$ skies. If the likelihood works as intended, we expect the recovered r values from the two sets of simulations to be extremely close to each other on a realization-by-realization basis. We find that the differences between the recovered r values $|\Delta r| \lesssim 0.002$. We also find that at our current noise level, even if we have perfect knowledge of the lensing B-modes in our patch, the uncertainty on r is reduced only by 26% from $\sigma(r) = 0.024$ to $\sigma(r) = 0.018$. This means that lensing uncertainty is subdominant compared to uncertainties from foregrounds and instrument noise in the BK14 dataset.

### B. Biases to r from misestimations of inputs

We investigate the bias to r from the following: (1) misestimation of the correlation between the CIB map and $\phi$, and (2) biases in polarization efficiency in the input $Q/U$ maps.

<footnote n="5">⁵Here we have taken $C_L^{\phi\phi}$ as the Planck 2013 cosmology used to generate the BICEP/Keck simulations introduced in Sec. III A 2. This is slightly different than the latest Planck cosmology which is implicit in the $C_L^{\text{I}\phi'}$ of Eq. (16). Arguably it would be more self-consistent to use the latest $C_L^{\phi\phi}$ here. However, we have checked that this makes no practical difference at the current sensitivity level.</footnote>
**Misestimation of $C_L^{1\phi}$:** As discussed in Sec. III B 1, we compute $C_L^{\text{II}}$ and $C_L^{1\phi'}$ used in the normalization and Wiener filter of the CIB map as the mean of the $C_L^{\text{II}}$ and $C_L^{1\phi'}$ from eight patches. In addition, we generate simulations of CIB based on the mean and scatter of the $C_L^{\text{II}}$ and $C_L^{1\phi'}$ spectra measured from the patches. Here, we consider the case in which the actual CIB cross-spectrum with $\phi$ is offset from the measured mean $C_L^{1\phi'}$. A plausible way in which the measured $C_L^{1\phi'}$ might be biased from the true $C_L^{1\phi}$ is through a bias in the $\phi'$ reconstruction due to CIB in the input temperature maps. In that case, the measured $C_L^{1\phi'}$ would contain a term that comes from CIB $\times \phi(\text{CIB, CIB})$, where $\phi(\text{CIB, CIB})$ denotes the CIB power that is leaked through the $\phi$ estimator applied to the CMB maps.

We construct a test for this bias, which proceeds as follows: using the measured mean $C_L^{\text{II}}$ and $C_L^{1\phi'}$, we generate simulated CIB skies as described in Sec. III B 2. This set of simulations is the assumed truth. We then generate two sets of CIB skies whose $C_L^{1\phi'}$ is either half a $\sigma_{\text{sp}}$ above or below the measured mean $C_L^{1\phi'}$, where $\sigma_{\text{sp}}$ is measured from the spread across the 8 patches. We process these CIB skies as if they had the mean $C_L^{1\phi'}$; i.e., we normalize and Wiener filter these maps using the mean $C_L^{1\phi'}$ and $C_L^{\text{II}}$. We then proceed to construct lensing templates and calculate auto- and cross-spectra with the rest of the BK14 maps, exactly as in the baseline analysis. The bandpower covariance matrix is derived using lensing templates constructed from the nominal, unbiased set of CIB skies. We then run maximum-likelihood searches on these two sets of simulations for the model parameters $r, A_d, A_{\text{sync}}, \beta_d$, and $\beta_s$. We determine the bias on $r$ by comparing the means of the maximum likelihood $r$ values from the nominal set and the half-$\sigma_{\text{sp}}$ offset sets. We observe that the mean $r$ is biased by $0.2\sigma$, where $\sigma$ denotes the uncertainty of the $r$ measurement (i.e., the width of the $r$ distribution of the nominal set).

To get a sense of how relevant this bias is, we compare the half-$\sigma_{\text{sp}}$ offset we introduce into the simulations to a worst-case scenario of CIB leakage in the reconstructed $\phi'$ map. Reference [44] estimated the term CIB $\times \phi(\text{CIB, CIB})$ using $\phi'$ reconstructed from the Planck 545 GHz maps without foreground cleaning and found the bias to be below ~5% for $L < 1024$, the $L$ range used in this work. A 5% bias is smaller than the half $\sigma_{\text{sp}}$ shift considered. Furthermore, the Planck lensing map used to calculate $C_L^{1\phi'}$ was constructed using the SMICA input maps that are foreground-suppressed. Therefore, we expect the $0.2\sigma$ bias to be an overestimate of potential biases from misestimating $C_L^{1\phi'}$.

**Misestimation of polarization efficiency:** The Planck Collaboration has found that their polarization efficiency calibration could potentially be biased at the 1%–2% level [see e.g., Table 9 of [45]]. The SPTpol Q/U maps are

calibrated using a Planck E-mode map [34]. Therefore, it is reasonable to ask by how much $r$ would be biased if the calibration of the input Q/U maps is biased.

We construct the test by artificially scaling the SPTpol Q/U simulated maps low by 1.7% and analyzing the maps as if they had the original amplitudes. In other words, similar to the half-$\sigma_{\text{sp}}$ $C_L^{1\phi'}$ shift test above, the rest of the pipeline is held identical and the only change is the input SPTpol Q/U maps. For simplicity, instead of using the combined Q/U map, we use only SPTpol simulated maps for this test. Comparing the mean of the recovered maximum likelihood $r$ values for the nominal set with that for the biased set, we find negligible differences. We conclude that biases at this level in the polarization efficiency are not an issue for this analysis.

### C. Impact on $\sigma(r)$ from variations of inputs

We investigate the impact on $\sigma(r)$ from two effects: (1) non-Gaussianities in the input CIB map, and (2) inclusion of patch-to-patch variation in $C_L^{\text{II}}$ and $C_L^{1\phi'}$ in the generation of the CIB realizations.

**Non-Gaussianity of the CIB:** As discussed in Sec. III B 2, we generate our CIB simulations based on the $\phi$ realizations used to lens the simulated CMB input skies. While the $\phi$ realizations are Gaussian, the true $\phi$ has some non-Gaussianities due to nonlinear growth of structure [46]. However, the contribution to lensing B-modes from non-Gaussian $\phi$ is subdominant over the angular scales considered [47]. It is thus sufficient to model $\phi$ and the portion of CIB that correlates with $\phi$, the signal term (Eq. (16), as Gaussian. In addition to the signal term, we simulate the noise term of the CIB $I_{\text{LM}}^N$—the portion of the CIB that does not correlate with $\phi$—as Gaussian realizations given the measured $C_L^{\text{II}}$, $C_L^{1\phi'}$, and the input $C_L^{\phi\phi}$. However, the CIB is known to be quite non-Gaussian; its bispectra at the angular scales relevant to this work have been measured by [48] with high signal-to-noise. Therefore, one could imagine that simulating the CIB $I_{\text{LM}}^N$ as Gaussian fluctuations would cause the lensing template fluctuation to be underestimated. With the underestimation of the lensing template fluctuation, $\sigma(r)$ would be underestimated. Here we estimate the impact on $\sigma(r)$ when we increase the lensing template fluctuation.

To get a handle on how much to increase the lensing template fluctuation, we build lensing templates using a simulated CIB sky from Websky mocks, [49] which are built based on an approximation to full N-body halo catalogs [50,51]. From the full-sky CIB realization, we make 80 cutouts of size similar to the BICEP/Keck patch, undeflect-and-difference the Q/U maps, and compute the lensing template bandpower variances. We then generate matching Gaussian realizations of CIB using the $C_L^{\text{II}}$ and $C_L^{1\phi}$ measured between the simulated CIB map and the corresponding $\phi$ map (provided as a $\kappa$ map, where
$\kappa = -\nabla^2\phi/2$). Using these Gaussian CIB realizations, we generate lensing templates and calculate their bandpower variances. For the $L$ range considered in this analysis, the ratio of the lensing template $1\sigma$ uncertainties between templates generated from Gaussian CIB and those from N-body based CIB is $0.97 \pm 0.07$. This suggests that the lensing template bandpower variance is sufficiently modeled using Gaussian simulations. Furthermore, [52] performed a similar test using galaxy densities as $\phi$ tracer and found that the difference in the lensing template covariance between the Gaussian and their simulations is within the Monte Carlo uncertainty of the number of simulations considered.

Since the above tests could still be limited by the number of non-Gaussian simulated CIB skies, we ask how much $\sigma(r)$ could be impacted because of some low level of unmodeled non-Gaussianity in the $\phi$ tracer. To do that, we increase the values in the lensing template autospectrum subblock of the bandpower covariance matrix by 10% and perform maximum-likelihood searches on the baseline set of simulations. The resultant $\sigma(r)$ estimated from the width of the $r$ value distribution is negligibly different to the baseline case. Therefore, we conclude that at the current level of noise, unmodeled non-Gaussianities of the CIB have negligible impact on the uncertainty of the $r$ measurement.

Patch-to-patch variation in $C_L^{II}$ and $C_L^{I\phi'}$: We construct the CIB realizations using samples of $C_L^{II}$ and $C_L^{I\phi'}$ drawn from the measured covariance of $C_L^{II}$ and $C_L^{I\phi'}$ across eight patches. By doing this we incorporate the patch-to-patch variation in the CIB auto- and cross-spectrum with $\phi'$ into the uncertainty on $r$. Here we check how large this effect is by comparing the $\sigma(r)$ estimated from a set of CIB simulations generated with fixed $C_L^{II}$ and $C_L^{I\phi'}$ with that estimated from a set of CIB simulations generated from a distribution of $C_L^{II}$ and $C_L^{I\phi'}$. We find that the $\sigma(r)$ from these two sets of simulations are compatible to within MC uncertainty. This means that the uncertainty on $r$ introduced by the uncertainties in $C_L^{II}$ and $C_L^{I\phi'}$ is subdominant compared with the noise and sample variance of the lensing templates.

Having estimated the biases to $r$ caused by possible biases in the CIB and $Q/U$ maps and found them to be small, and having shown the impact on $\sigma(r)$ due to unmodeled non-Gaussianity of the CIB to be minimal, we now turn to testing the robustness of the simulations against unmodeled Galactic foregrounds using the data themselves.

## V. SYSTEMATICS CHECKS

Previous BICEP/Keck papers include “jackknife” internal consistency tests on the 95 and 150 GHz maps used here [12,30,53]. In this section, we provide similar tests of the

auto- and cross-spectra of the newly introduced lensing template. We consider the following ways in which the simulations can fail to sufficiently describe the statistics of the lensing template:

(1) Galactic dust in the input 150 GHz $Q/U$ maps leaks into the lensing template,

(2) low-$\ell$ systematic residuals in the Planck polarization maps leak into the lensing template,

(3) non-Gaussian Galactic dust residuals in the CIB map introduce extra power in the lensing template beyond that described by Gaussian modeling of uncorrelated power.

All of the above would (i) increase the power of the lensing template autospectrum, and (ii) introduce potential chance coupling with the observed B-modes.

Galactic dust power is subdominant to E-mode power over the angular scales relevant to producing the lensing B-mode template, and the simulated $Q/U$ maps used in Sec. III A 2 do not include a dust component. However, we would still like to check that the Galactic dust component in the $Q/U$ data maps does not significantly contribute to the lensing template auto-spectrum. For the CIB map, any components that contribute to the CIB autospectrum but are uncorrelated with $\phi'$ are modeled as Gaussian fluctuations. Therefore, the unmodeled non-Gaussian Galactic foregrounds could contribute extra fluctuation in the lensing templates when used to deflect the CMB maps. In addition, they could contribute extra template power when deflecting the unmodeled Galactic foregrounds in the $Q/U$ maps.

To address the question of whether the simulations are a sufficient description of the data given these unmodeled effects, we test the consistency of the lensing template auto- and cross-spectra against simulations. Specifically, we perform spectrum-difference tests where we compare the difference spectrum of data between the baseline $\ell$ and $L$ ranges and variant $\ell$ and $L$ ranges against the corresponding differences in simulation. We calculate two quantities, $\chi^2$ and $\chi$, as follows. Firstly

$$ \chi_{\text{sys}}^2 = \Delta C_{\ell}^{\dagger} \text{Cov}^{-1} \Delta C_{\ell}, \quad (17) $$

where $\Delta C_{\ell}$ denotes the binned data difference spectrum and Cov is the bandpower covariance matrix formed from the difference spectra of the corresponding simulations. And secondly,

$$ \chi_{\text{sys}} = \sum_{\ell} \Delta C_{\ell} / \sigma_{\ell, \text{diff}}, \quad (18) $$

where $\sigma_{\ell, \text{diff}}$ denotes the standard deviation from the simulation difference spectra. Figure 6 shows the difference spectra for the lensing template autospectrum (LT $\times$ LT), lensing template cross-spectrum with the BK14 95 GHz map (LT $\times$ BK14$_{95}$), and lensing template cross-spectrum
<figure><img src="image_8.png" /><figcaption>FIG. 6. Difference bandpowers ($\Delta C_{\ell}$, see definition in text) between the baseline analysis and analyses with one parameter changed, and the uncertainties on those difference bandpowers, both scaled by the statistical uncertainties on the baseline analysis bandpowers. The label at the top left-hand corner of each row indicates which parameter has been modified and how it is modified. The left to right columns show the difference bandpowers from the lensing template autospectrum, lensing template cross-spectrum with the BK14 95 GHz map, and lensing template cross-spectrum with the BK14 150 GHz map. The gray bands indicate the $0.5\sigma$ statistical uncertainty of the baseline spectra. The $\chi^2_{sys}$ and PTE of the difference bandpowers are listed in Table I. We find the data difference bandpowers to be consistent with the spread in the simulation difference bandpowers.</figcaption></figure>

with the BK14 150 GHz map (LT × BK14$_{150}$). The PTE values from $\chi^2_{sys}$ and $\chi_{sys}$ are listed in Table I.

### A. L-cuts on CIB map

At large angular scales the CIB map could be contaminated by Galactic dust and thus a test in which the $L_{min}$ for the CIB map is varied could be sensitive to its impact. The unmodeled non-Gaussianity of residual Galactic foregrounds in the CIB map would cause the lensing template to have larger variance than it would otherwise. We test the hypothesis that the simulations are sufficient descriptions of the real data by differencing the lensing template auto- and cross-spectra generated using the baseline $L_{min} = 100$ for the CIB map and those generated with $L_{min} = 50$

<table><caption>TABLE I. The PTE values from $\chi^2_{sys}$ and $\chi_{sys}$ (separated by a comma) with different CIB input $L_{min}$ and $Q/U$ map input $\ell_{min}$, compared with the baseline setup. LT × LT, LT × 95, and LT × 150 denote the lensing template autospectrum, lensing template cross-spectrum with the BK14 95 GHz map, and with the BK14 150 GHz map, respectively.</caption><thead><tr><th>Variation /spectrum</th><th>LT × LT</th><th>LT × 95</th><th>LT × 150</th></tr></thead><tbody><tr><td>CIB L<sub>min</sub> = 50</td><td>0.36, 0.12</td><td>0.80, 0.23</td><td>0.66, 0.09</td></tr><tr><td>CIB L<sub>min</sub> = 150</td><td>0.91, 0.67</td><td>0.68, 0.63</td><td>0.25, 0.88</td></tr><tr><td>Q/U l<sub>min</sub> = 100</td><td>0.76, 0.70</td><td>0.09, 0.84</td><td>0.34, 0.52</td></tr><tr><td>Q/U l<sub>min</sub> = 200</td><td>0.76, 0.75</td><td>0.36, 0.57</td><td>0.28, 0.62</td></tr></tbody></table>

and $L_{min} = 150$. The PTEs from the $L_{min}$ difference spectra show that the data differences are sufficiently described by the simulation-difference distributions.

### B. $\ell$-cuts on $Q/U$ maps

Galactic dust contributes a fraction of the total power in the $Q/U$ maps on the largest scales. Additionally, there could be low levels of unmodeled systematic residuals [54] in the Planck $Q/U$ maps that could leak power to the lensing templates. Similar to the test done with the CIB map, we set the $\ell_{min}$ of the input $Q/U$ map to two different levels compared to the baseline (no explicit $\ell_{min}$ set) and compute difference spectra between the variant $\ell_{min}$ and the baseline. For $\ell_{min} = 100$ and $\ell_{min} = 200$, we find the PTEs from the difference spectra to be consistent with the simulation-difference distributions.

We thus conclude that at the current level of noise, the lensing template auto-spectrum and the cross-spectra with the 95 GHz and 150 GHz maps do not contain unmodeled systematics from large angular scales of the input $Q/U$ and CIB maps large enough to be incompatible with the simulation distributions.

## VI. RESULTS

We now proceed to repeat the parameter constraint analysis from the BK14 paper [12] including the lensing template extension described and validated above. We
<figure><img src="image_9.png" /><figcaption>FIG. 7. BB auto- and cross-spectra calculated using BICEP2/Keck 95 and 150 GHz maps, the Planck 353 GHz map, and the lensing template developed in this paper. The black lines show the model expectation values for lensed-ΛCDM, while the red lines show the expectation values of the baseline lensed-ΛCDM + dust model from the BK14 analysis ($r = 0, A_d = 4.3 \ \mu\text{K}^2, \beta_d = 1.6, \alpha_d = -0.4$), and the error bars are scaled to that model. Compared to the BK14 BB spectra, which contain both foregrounds and lensing components, the lensing template represents an alternate way to estimate the lensing B-modes which is largely foreground-immune, and, as we see here, provides good signal-to-noise in the resulting auto- and cross-spectra.</figcaption></figure>

present two main results in this work. First, we estimate $\sigma(r)$ with delensing by running maximum-likelihood searches on the set of lensed-ΛCDM + dust + noise simulations from BK14. Second, we explore the likelihood space of the real data and provide constraints on $r$ and the foreground model parameters.

In Sec. III, we described the construction of a lensing template using the Planck GNILC CIB map and the combined Q/U maps from SPTpol, BICEP/Keck, and Planck. Figure 7 shows the auto- and cross-spectra of this lensing template with the maps that most significantly constrain the model parameters—the BICEP/Keck 95 and 150 GHz maps, and the Planck 353 GHz map. The lensing template auto- and cross-spectra shown in Fig. 7, plus the

<table><caption>TABLE II. Priors imposed on each parameter for both maximum-likelihood search and posterior sampling for the baseline analysis. $U(a, b)$ denotes uniform distribution between $[a, b]$. $N(\mu, \sigma^2)$ denotes normal distribution with mean $\mu$ and variance $\sigma^2$.</caption><thead><tr><th>Parameter</th><th>ML search</th><th>Sampling</th></tr></thead><tbody><tr><td>r</td><td>U(-0.5, 0.5)</td><td>U(0, 0.5)</td></tr><tr><td>A<sub>d</sub></td><td>U(-2, 15)</td><td>U(0, 15)</td></tr><tr><td>A<sub>sync</sub></td><td>U(-2, 15)</td><td>U(0, 50)</td></tr><tr><td>β<sub>d</sub></td><td>N(1.6, 0.11<sup>2</sup>)</td><td>N(1.59, 0.11<sup>2</sup>)</td></tr><tr><td>β<sub>s</sub></td><td>N(-3.1, 0.3<sup>2</sup>)</td><td>N(-3.1, 0.3<sup>2</sup>)</td></tr><tr><td>α<sub>d</sub></td><td>U(-1, 0)</td><td>U(-1, 0)</td></tr><tr><td>α<sub>s</sub></td><td>Fixed</td><td>U(-1, 0)</td></tr><tr><td>ε</td><td>Fixed</td><td>U(0, 1)</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>FIG. 8. Histograms of maximum-likelihood values of $r$, $A_d$, and $A_{sync}$ from 499 realizations of BK14 + LT (blue) and BK14 (gray) lensed-ΛCDM + dust + noise simulations in the baseline model with six free parameters: $r$, $A_d$, $A_{sync}$, $\beta_d$, $\beta_s$ and $\alpha_d$. The red lines mark the means of the distributions for the BK14 + LT simulation set, and the gray dashed lines mark the input values. $\sigma(r)$ from the BK14 + LT (BK14) simulation set is 0.022 (0.024) from the leftmost panel.</figcaption></figure>

additional cross-spectra with the other bands of WMAP and Planck, are the new additions to the bandpower data vector input to the likelihood analysis. It is interesting to note that the error bars are much smaller at low $\ell$ for $LT \times LT$ than for $LT \times 150$. This is because, although the 150 GHz map noise is very small, the dust sample variance is large.

## A. Reduction in $\sigma(r)$

The inclusion of the lensing template cross-spectra reduces the effective sample variance of the lensing component of the observed B-modes. This is the reason that the uncertainty of the *r* component can be reduced when we add a lensing template to the likelihood.

In BK14, we introduced $\sigma(r)$ as a measure of the intrinsic constraining power of a given set of experimental data. In contrast to the width of the 68% highest posterior density interval as derived from the real data this measure is not subject to noise fluctuation within that single realization. To compare the $\sigma(r)$ from the BK14 dataset and the BK14 dataset with lensing template included, we repeat the analysis of Appendix E.3 of the BK14 paper. We run maximum-likelihood searches with the baseline lensed-ΛCDM + dust + synchrotron + *r* model on the lensed-ΛCDM + dust + noise simulations for the two cases. The parameters and priors are the same as in BK14 and are summarized in Table II. The amplitudes at $\ell = 80$ of the dust and synchrotron BB spectra defined at 353 GHz and 23 GHz are denoted by $A_d$ and $A_{sync}$, respectively; $\beta$ and $\alpha$ denote the frequency and spatial spectral indices, with subscripts *d* and *s* referring to dust and synchrotron respectively; $\epsilon$ denotes the dust-synchrotron correlation. Flat priors are applied to *r*, $A_d$, $A_{sync}$ & $\alpha_d$, and Gaussian priors are applied to $\beta_d$ and $\beta_s$. Figure 8 shows the distributions of maximum likelihood *r*, $A_d$, and $A_{sync}$ values. With the inclusion of the lensing

template, we reduce $\sigma(r)$ from 0.024 to 0.022, a ~10% reduction.<fnref n="6" />

We also generate simulated lensing templates using only one of SPTpol, BICEP/Keck, and Planck for the input $Q/U$ maps. We add the single-experiment lensing template to the BK14 simulation set and perform maximum-likelihood searches. We find that the $\sigma(r)$ from LT$_{SPTpol}$, LT$_{BICEP/Keck}$, and LT$_{Planck}$ to be 0.0223, 0.0230, and 0.0236 respectively.<fnref n="7" /> This shows that the SPTpol $Q/U$ maps contribute most to recovering the lensing B-modes. The fact that LT$_{BICEP/Keck}$ contributes more than LT$_{Planck}$ shows that the signal-to-noise per mode at low $\ell$ is more important than having a wider range in $\ell$ for the particular combination of the $\ell$ range and noise levels between BICEP/Keck and Planck.

## B. Parameter posteriors of BK14 with delensing

We now repeat the eight-parameter likelihood evaluation of the real data as in the BK14 paper. We again use COSMOMC [55] and the lensed-ΛCDM + dust + synchrotron + *r* model with parameters and priors summarized in Table II. Figure 9 shows the posterior distributions of the baseline analysis compared with the BK14 result. The peak and 68% credible regions of the marginalized *r* distribution are shifted down from the BK14 values of $0.028_{-0.025}^{+0.026}$ to $0.027_{-0.022}^{+0.023}$ when the lensing template is

<footnote n="6">⁶Note that this $\sigma(r)$ is computed in a six-dimensional parameter space as opposed to the eight-dimensional parameter space which is used when sampling. This is to maintain consistency with the BK14 paper. For a 8D search, $\sigma(r) = 0.026$ without the lensing template and $\sigma(r) = 0.023$ with. The relevant metric here is the fractional reduction in $\sigma(r)$ between the two simulation sets which is similar for the 6D and 8D searches.</footnote>

<footnote n="7">⁷The three-experiment Q/U combined LT gives $\sigma(r) = 0.0221$. We provide three significant figures for comparisons between the templates.</footnote>
<figure><img src="image_11.png" /><figcaption>FIG. 9. Posterior distributions of the baseline model parameters given the BK14 + LT dataset (black lines) compared with the BK14 dataset (red lines, which are the same as the black lines in Fig. 4 of the BK14 paper). The lensing template is constructed using combined Q/U maps from SPTpol, BICEP/Keck, and Planck (Sec. III A 3) and a CIB map as the $\phi$ tracer (Sec. III B 1). The 95% C.L. upper limit on the tensor-to-scalar ratio tightens from $r_{0.05} < 0.090$ to $r_{0.05} < 0.082$ with the addition of the lensing template. The parameters $A_d$ and $A_{sync}$ are the amplitudes of the dust and synchrotron B-mode spectra, where $\beta$ and $\alpha$ are the frequency and spatial spectral indices respectively. The dust-synchrotron correlation parameter is denoted by $\epsilon$. The up-turn of the 1D posterior distribution of $A_{sync}$ as it approaches zero comes from the increased volume allowed by the $\epsilon$ parameter as $\epsilon$ becomes ambiguous when $A_{sync} = 0$. In the 1D panels for the $\alpha, \beta$, and $\epsilon$ parameters, the blue dashed lines denote the priors for each parameter.</figcaption></figure>

included. The 95% C.L. upper limit on $r_{0.05}$ is reduced from 0.090 to 0.082. Some of the other constraints are $A_d = 4.2_{-0.9}^{+1.1} \mu\text{K}^2$ and $A_{sync} < 3.7 \mu\text{K}^2$ (95% C.L.).<fnref n="8" /> The maximum-likelihood model (including priors) in the 8D parameter space is: $r_{0.05} = 0.025$, $A_d = 4.0 \mu\text{K}^2$, $A_{sync} = 1.4 \mu\text{K}^2$, $\beta_d = 1.6$, $\beta_s = -3.1$, $\alpha_d = -0.17$, $\alpha_s = -0.95$, and $\epsilon = 0.00$. Against

<footnote n="8">⁸As noted, the model space is identical to BK14 to enable apples-to-apples comparison. However, we have since then made one model change in BK15 [14] and widened the prior range of the dust-synchrotron correlation parameter $\epsilon$ from $0 < \epsilon < 1$ to $-1 < \epsilon < 1$ (see Appendix E1 in BK15 for details). With this prior, the BK14 r peak and 68% credible regions reduce from $0.043_{-0.028}^{+0.031}$ to $0.038_{-0.024}^{+0.029}$ when a lensing template is included.</footnote>
<figure><img src="image_12.png" /><figcaption>FIG. 10. Posterior distributions on $r$ and $A_L$, a parameter used to scale the lensing BB power, from an alternative analysis in which the amplitude of lensing is a free parameter. With the addition of the lensing template, the probability of shuffling lensing power to other parameters is reduced, thus the degeneracy between $r$ and $A_L$ is reduced.</figcaption></figure>

this model, we compute $\chi^2 = (d-m)^{\dagger}\text{Cov} - 1(d-m) = 768$ for the 9 × 78 = 702 data bandpowers. We compare this number against the distribution in simulations finding a PTE of 0.15. We conclude that the model is a sufficient description of the data at present.

We perform a couple of variations to the baseline analysis to explore degeneracies amongst model parameters that are important to lensing and changes in $r$ with different input datasets. In the baseline analysis, the lensing BB spectrum is taken as the ΛCDM expectation in both normalization and shape. As an alternative we rescale this spectrum by the parameter $A_L$ and sample the posterior distribution in the ΛCDM + $A_L$ model space. Secondly, as is done in Sec. VI A, we form input lensing templates using $Q/U$ maps from one of the three experiments instead of combining them. We discuss the results of each variation in the following paragraphs.

When we allow $A_L$ to float, we note a $A_L - r$ degeneracy in the BK14 dataset, as shown in Fig. 10, and as was previously noted in an earlier BICEP/Keck analysis [15]. When the lensing template is added to the BK14 dataset, the degeneracy between $r$ and $A_L$ is reduced. In this model space, the peak and 68% credible regions of the marginalized $r$ distribution with and without the lensing template are $0.025_{-0.022}^{+0.023}$ and $0.009_{-0.009}^{+0.031}$, and the upper limits on $r$ are $r_{0.05} < 0.081$ and $r_{0.05} < 0.079$ respectively. The peak and 68% credible regions of $A_L$ with and without the lensing template are $1.03 \pm 0.10$ and $1.21 \pm 0.17$.

<figure><img src="image_13.png" /><figcaption>FIG. 11. The $r$ posterior curves from the baseline analysis, along with $r$ curves from analyses using lensing templates constructed from $Q/U$ maps from only one of the three experiments: BICEP/Keck, Planck, and SPTpol. The shifts in the curves are consistent with expectations from simulations.</figcaption></figure>

respectively.<fnref n="9" /> The shift in the peak $A_L$ is consistent with expectations from simulations, where 25% of the simulation realizations have $A_L$ shifts with absolute magnitude larger than that seen in data. We see that with the addition of the lensing template, we are able to better constrain the lensing power in the measured auto- and cross-spectra across the different frequencies and thereby reduce the probability of misassigning power to lensing.

We show in Fig. 11 the $r$ posterior distributions from analyses in the lensed-ΛCDM model space using lensing templates constructed from $Q/U$ maps coming from only one of the three experiments, SPTpol, BICEP/Keck and Planck. We see that the peaks of the $r$ posteriors from the BICEP/Keck-only and the Planck-only cases are close to the baseline case, while the width of the $r$ posterior from the Planck-only case is a bit larger than the baseline case. The larger $r$ posterior uncertainty is expected given the larger $\sigma(r)$ from the Planck-only simulation set in Sec. VI A. The peak of the $r$ posterior for the SPTpol-only case is shifted up slightly compared with the baseline case. This might seem slightly surprising given that the SPTpol $Q/U$ maps contribute most of the weight in the combined $Q/U$ maps over a broad range of angular scales. To quantify the probability of the observed shift between the baseline case and the SPTpol-only case, we extract the best-fit $r$ values from the baseline simulation set and the SPTpol-only case.

<footnote n="9">⁹We note that we have kept fixed a component of the noise bias in the LT autospectrum ($s_\phi * n_{QU}$ in Eq. (11)) which varies with $A_L$. It contributes <10% of the total noise bias and is only present in the LT × LT part of the data vector. Varying this noise component with $A_L$ would slightly tighten the constraint on $A_L$, but the qualitative conclusion would be changed.</footnote>
lensing template simulation set. Restricting to the subset
with positive best-fit *r* in the baseline setup, we count the
fraction of realizations that have larger best-fit *r* differences
between the SPTpol-only and the baseline set than is seen
in the data. We find 20% of the simulations fit this criterion,
and thus we conclude that what is observed in the data is
typical of the expected fluctuations.

VII. CONCLUSION

In this work, we build on the BICEP/Keck analysis
framework and demonstrate, for the first time, improve-
ments to constraints on the tensor-to-scalar ratio *r* with
delensing. With the addition of a lensing template, we
reduce the uncertainty of the *r* estimate by constraining
the lensing *B*-mode contribution to the observed *B*-modes.
We construct the lensing template using an undeflect-and-
difference approach, in which we undeflect the observed
Q/U maps by a *φ* tracer, and then difference the unde-
flected maps from the input maps. The Q/U maps we use
are a 150 GHz combination of SPTpol observations from
2013–2015, BICEP/Keck observations up to 2014, and the
Planck satellite full-mission observations. The *φ* tracer we
use is a CIB map constructed using the GNILC algorithm
from Planck data. The resulting lensing template is added
as a pseudofrequency band to the BK14 dataset, in which
BICEP/Keck WMAP and Planck maps are used to constrain
Galactic foregrounds and *r*.

We present two key results from this analysis. First, we
estimate σ(r) using our lensed-ACDM + dust + noise
simulation set. We find maximum likelihood values of
the baseline model parameters for each simulation reali-
zation and take the mean and standard deviation over the
499 realizations. We find that, with the addition of the
lensing template, σ(r) improves from 0.024 in BK14 to
0.022, a ~10% improvement. The second main result is the
posterior peak value, 68% credible region, and upper limit
on r when we add the lensing template to the BK14 dataset.
With delensing, the peak and 68% credible regions shift
from r = 0.028<sup>+0.026</sup><sub>-0.025</sub> to r = 0.027<sup>+0.023</sup><sub>-0.022</sub>, and the 95% C.L.
upper limit on r is reduced from 0.090 to 0.082.

We estimate the impact on *r* from potential biases in the
inputs used to construct the simulated lensing templates.
We find the biases to *r* from misestimating the cross-
spectrum of the CIB and *φ* to be small, and the biases to *r*
from biases in polarization efficiency of the CMB Q/U
maps to be negligible. We find negligible difference in *σ(r)*
due to modeling the non-Gaussian CIB field as Gaussian
for this dataset, and that the uncertainties in the CIB
autospectrum and the CIB × *φ* cross-spectrum contribute
subdominantly to *σ(r)*. We perform checks against poten-
tial unmodeled systematic contaminations to the lensing
template. This includes Galactic foregrounds leaking into
the lensing template through either the input Q/U maps or
the input CIB map. We show that the data lensing template
is sufficiently well-described by the simulations. Therefore

we conclude that the results are robust against these sources
of systematics given the current noise levels.

At the BK14 level of map noise and Galactic foreground
variance, simulations indicate that perfect delensing would
reduce σ(r) from 0.024 to 0.018. This implies that the
variance from lensing B-modes is not the dominant source
of uncertainty (<30%) when constraining r in this dataset.
However, with current and upcoming ground-based CMB
telescopes, e.g., BICEP Array [56], SPT-3G [57], AdvACT,
Simons Array, Simons Observatory [58], and CMB-S4
[59], the millimeter-wave sky will be mapped with ever
higher signal-to-noise. Lensing B-modes will become a
dominant source of uncertainty, and delensing will be
crucial to break the floor of σ(r) set by the lensing variance.
For example, while in the most recent BICEP/Keck r analysis
BK15 [14] lensing variance continues to be subdominant,
in the upcoming result BK18 lensing variance contributes
roughly half of the r uncertainty budget. Projecting further,
without delensing, the BICEP Array experiment σ(r) would
plateau at ~0.006. However, this σ(r) could be reduced by
a factor of about 2.5 with delensing using a φ field
reconstructed using CMB maps from the SPT-3G experi-
ment. This is a much more significant reduction in the
uncertainty on r than is achieved in this work.

To reach the target $\sigma(r)$ of $5 \times 10^{-4}$ for the next-
generation ground-based CMB experiment CMB-S4, more
than 90% of the lensing sample variance needs to be
removed [60]. Delensing to such low residual levels
requires high values of $\rho_L$, the correlation between the $\phi$
tracer and the underlying $\phi$ field. In addition to using $\phi$
maps reconstructed from low-noise, high-resolution CMB
observations [e.g., [24,61,62]], higher $\rho_L$ tracers could be
obtained by combining different tracers [e.g., [63,64]] and
using optimal methods [e.g., [65–68]]. We will be exploring
various approaches to delensing [e.g., [67]] in future joint
analyses of BICEP/Keck and SPT-3G data, confronting
delensing algorithms with real-world nonidealities and
developing techniques to mitigate systematics, readying
our analysis for the future of low-noise data and the
possibility of detecting PGWs.

ACKNOWLEDGMENTS

The authors thank Dominic Beck and Chang Feng
for useful comments on an early version of the draft.
The BICEP2/Keck Array projects have been made
possible through a series of grants from the National
Science Foundation including Grants No. 0742818,
No. 0742592, No. 1044978, No. 1110087, No. 1145172,
No. 1145143, No. 1145248, No. 1639040, No. 1638957,
No. 1638978, and No. 1638970, and by the Keck
Foundation. The development of antenna-coupled detector
technology was supported by the JPL Research and
Technology Development Fund, and by NASA Grants
No. 06-ARPA206-0040, No. 10-SAT10-0017, No. 12-
SAT12-0031, No. 14-SAT14-0009, and No. 16-SAT-16-
0002. The development and testing of focal planes were supported by the Gordon and Betty Moore Foundation at Caltech. Readout electronics were supported by a Canada Foundation for Innovation grant to UBC. Support for quasioptical filtering was provided by UK STFC Grant No. ST/N000706/1. Some of the computations in this paper were run on the Odyssey cluster supported by the FAS Science Division Research Computing Group at Harvard University. The analysis effort at Stanford and S. L. A. C. is partially supported by the U.S. DOE Office of Science. We thank the staff of the U.S. Antarctic Program and in particular the South Pole Station without whose help this research would not have been possible. Most special thanks go to our heroic winter-overs Robert Schwarz and Steffen Richter. We thank all those who have contributed past efforts to the BICEP-Keck Array series of experiments, including the BICEP1 team. S. P. T. is supported by the National Science Foundation through Grants No. PLR-1248097 and No. OPP-1852617. Partial support is also provided by the NSF Physics Frontier Center Grant No. PHY-1125897 to the Kavli Institute of Cosmological Physics at the University of Chicago, the Kavli Foundation and the Gordon and Betty Moore Foundation Grant No. GBMF 947. This research used resources of the National Energy Research Scientific Computing Center (NERSC), a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The Melbourne group acknowledges support from the University of Melbourne and an Australian Research Council's Future Fellowship (FT150100074). 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. Work at McGill is supported by the Natural Science and Engineering Research Council of Canada, the Canadian Institute for Advanced Research, and M. D. acknowledges a Killam research fellowship W. L. K. W. is supported in part by the Kavli Institute for Cosmological Physics at the University of Chicago through Grant No. NSF PHY-1125897, an endowment from the Kavli Foundation and its founder Fred Kavli, and by the Department of Energy, Laboratory Directed Research and Development program and as part of the Panofsky Fellowship program at SLAC National Accelerator Laboratory, under Contract No. DE-AC02-76SF00515. B. B. is supported by the Fermi Research Alliance LLC under Contract No. De-AC02-07CH11359 with the U.S. Department of Energy. We acknowledge the use of many Python packages: IPython [69], Matplotlib [70], SciPy [71], and HEALPY [72,73]. We also thank the Planck and WMAP teams for the use of their data. Some of the sky simulations used in this paper were developed by the WebSky Extragalactic CMB Mocks team, with the continuous support of the Canadian Institute for Theoretical Astrophysics (CITA), the Canadian Institute for Advanced Research (CIFAR), and the Natural Sciences and Engineering Research Council of Canada (NSERC), and were generated on the Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund—Research Excellence; and the University of Toronto.

## APPENDIX: LENSING TEMPLATE CONSTRUCTION METHODS

In this paper, we have used a map-space “undeflect-and-difference” method to construct the lensing template. Previous works have inferred the lensing B-modes $\hat{B}^{\text{lens}}(\ell)$ by lensing the observed E-modes to first order in $\phi$ given a $\phi$ tracer. Specifically,

$$ \hat{B}^{\text{lens}}(\ell) = \int \frac{d^2\ell'}{(2\pi)^2} W(\ell, \ell') \bar{\phi}(\ell - \ell') \bar{E}(\ell'), \quad (A1) $$

where $W(\ell, \ell') = \ell' \cdot (\ell - \ell') \sin(2\phi_{\ell, \ell'})$, and $\bar{E}$ and $\bar{\phi}$ are the Wiener filtered E-modes and $\phi$ tracer respectively [e.g., [16]]. An advantage to this formulation is that by acting on the E-modes of the observed sky only, noise in the lensing template is reduced versus the undeflect-and-difference method. This extra noise enters by undeflecting $Q/U$ maps which also contain B-modes. While the signal contribution from the B-modes is small, the level of noise fluctuations is similar to those in the E-modes, which contribute noise to the undeflect-and-difference templates. However, this is not a fundamental limitation to the map-space approach as implemented in this paper. One could Fourier transform the $Q/U$ maps to E/B-modes, null the B-modes, and then transform back to $Q/U$ maps before performing the undeflect-and-difference operation, which would remove this specific noise. In fact, we experimented with adding these steps and found that for the present case, the reduction in lensing template noise is fractionally very small for $\ell < 500$.

For future analyses, we will revisit the algorithm used to produce the lensing template to further improve its signal-to-noise. Besides removing the extra noise contribution, other possible improvements include filling in the region outside the SPTpol coverage (as seen in Fig. 2) using the information available from BICEP/Keck and Planck.
[1] M. Kamionkowski and E. D. Kovetz, The quest for B modes from inflationary gravitational waves, Annu. Rev. Astron. Astrophys. **54**, 227 (2016).

[2] Y. Akrami, F. Arroja et al. (Planck Collaboration), Planck 2018 results. X. Constraints on inflation, Astron. Astrophys. **641**, A10 (2020).

[3] U. Seljak and M. Zaldarriaga, Signature of Gravity Waves in the Polarization of the Microwave Background, Phys. Rev. Lett. **78**, 2054 (1997).

[4] M. Kamionkowski, A. Kosowsky, and A. Stebbins, A Probe of Primordial Gravity Waves and Vorticity, Phys. Rev. Lett. **78**, 2058 (1997).

[5] R. Adam, P. A. R. Ade et al. (Planck Collaboration), Planck intermediate results. XXX. The angular power spectrum of polarized dust emission at intermediate and high Galactic latitudes, Astron. Astrophys. **586**, A133 (2016).

[6] N. Krachmalnicoff, E. Carretti, C. Baccigalupi et al., S-PASS view of polarized Galactic synchrotron at 2.3 GHz as a contaminant to CMB observations, Astron. Astrophys. **618**, A166 (2018).

[7] A. Lewis and A. Challinor, Weak gravitational lensing of the CMB, Phys. Rep. **429**, 1 (2006).

[8] D. Hanson, S. Hoover, A. Crites et al., Detection of B-Mode Polarization in the Cosmic Microwave Background with Data from the South Pole Telescope, Phys. Rev. Lett. **111**, 141301 (2013).

[9] P. A. R. Ade, Y. Akiba et al. (Polarbear Collaboration), A measurement of the cosmic microwave background B-mode polarization power spectrum at sub-degree scales with POLARBEAR, Astrophys. J. **794**, 171 (2014).

[10] R. Keisler, S. Hoover, N. Harrington et al., Measurements of sub-degree B-mode polarization in the cosmic microwave background from 100 square degrees of SPTpol data, Astrophys. J. **807**, 151 (2015).

[11] T. Louis, E. Grace, M. Hasselfield et al., The Atacama Cosmology Telescope: Two-season ACTPol spectra and parameters, J. Cosmol. Astropart. Phys. 6 (2017) 031.

[12] P. A. R. Ade et al. (BICEP2 Collaboration and Keck Array Collaboration), Improved Constraints on Cosmology and Foregrounds from BICEP2 and Keck Array Cosmic Microwave Background Data with Inclusion of 95 GHz Band, Phys. Rev. Lett. **116**, 031302 (2016).

[13] P. A. R. Ade, M. Aguilar et al. (POLARBEAR Collaboration), A measurement of the cosmic microwave background B-mode polarization power spectrum at subdegree scales from two years of polarbear data, Astrophys. J. **848**, 121 (2017).

[14] P. A. R. Ade et al. (BICEP2 Collaboration and Keck Array Collaboration), Constraints on Primordial Gravitational Waves Using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season, Phys. Rev. Lett. **121**, 221301 (2018).

[15] P. A. R. Ade et al. (BICEP2/Keck Collaboration and Planck Collaboration), Joint Analysis of BICEP2/Keck Array and Planck Data, Phys. Rev. Lett. **114**, 101301 (2015).

[16] A. Manzotti, K. T. Story, W. L. K. Wu et al., CMB polarization B-mode delensing with SPTpol and Herschel, Astrophys. J. **846**, 45 (2017).

[17] J. Carron, A. Lewis, and A. Challinor, Internal delensing of Planck CMB temperature and polarization, J. Cosmol. Astropart. Phys. 5 (2017) 035.

[18] N. Aghanim, Y. Akrami et al. (Planck Collaboration), Planck 2018 results. VIII. Gravitational lensing, Astron. Astrophys. **641**, A8 (2020).

[19] S. Adachi, M. A. O. Aguilar Faúndez, Y. Akiba et al., Internal Delensing of Cosmic Microwave Background Polarization B-Modes with the POLARBEAR Experiment, Phys. Rev. Lett. **124**, 131301 (2020).

[20] D. Han, N. Sehgal, A. MacInnis et al., The Atacama Cosmology Telescope: Delensed power spectra and parameters, arXiv:2007.14405.

[21] W. Hu and T. Okamoto, Mass reconstruction with cosmic microwave background polarization, Astrophys. J. **574**, 566 (2002).

[22] G. Simard, D. Hanson, and G. Holder, Prospects for delensing the cosmic microwave background for studying inflation, Astrophys. J. **807**, 166 (2015).

[23] B. D. Sherwin and M. Schmittfull, Delensing the CMB with the cosmic infrared background, Phys. Rev. D **92**, 043005 (2015).

[24] W. L. K. Wu, L. M. Mocanu, P. A. R. Ade et al., A measurement of the cosmic microwave background lensing potential and power spectrum from 500 deg² of SPTpol temperature and polarization data, Astrophys. J. **884**, 70 (2019).

[25] N. Aghanim, M. Ashdown et al. (Planck Collaboration), Planck intermediate results. XLVIII. Disentangling Galactic dust emission and cosmic infrared background anisotropies, Astron. Astrophys. **596**, A109 (2016).

[26] P. A. R. Ade, N. Aghanim et al. (Planck Collaboration), Planck 2015 results. XV. Gravitational lensing, Astron. Astrophys. **594**, A15 (2016).

[27] W. Hu, Weak lensing of the CMB: A harmonic approach, Phys. Rev. D **62**, 043007 (2000).

[28] E. Anderes, B. D. Wandelt, and G. Lavaux, Bayesian inference of CMB gravitational lensing, Astrophys. J. **808**, 152 (2015).

[29] D. Green, J. Meyers, and A. van Engelen, CMB delensing beyond the B modes, J. Cosmol. Astropart. Phys. 12 (2017) 005.

[30] P. A. R. Ade, R. W. Aikin et al. (BICEP2 Collaboration), Detection of B-Mode Polarization at Degree Angular Scales by BICEP2, Phys. Rev. Lett. **112**, 241101 (2014).

[31] S. Hamimeche and A. Lewis, Likelihood analysis of CMB temperature and polarization power spectra, Phys. Rev. D **77**, 103013 (2008).

[32] J. E. Austermann, K. A. Aird, J. A. Beall et al., SPTpol: An instrument for CMB polarization measurements with the South Pole Telescope, in *Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy VI*, Vol. 8452 of Proc. SPIE (SPIE, Bellingham, 2012), p. 84521E.

[33] J. E. Carlstrom, P. A. R. Ade, K. A. Aird et al., The 10 meter South Pole Telescope, Publ. Astron. Soc. Pac. **123**, 568 (2011).

[34] J. W. Henning, J. T. Sayre, C. L. Reichardt et al., Measurements of the temperature and E-mode polarization of the
CMB from 500 square degrees of SPTpol data, Astrophys. J. **852**, 97 (2018).

[35] https://cosmologist.info/lenspix/

[36] K. M. Smith, D. Hanson, M. LoVerde, C. M. Hirata, and O. Zahn, Delensing CMB polarization with external datasets, J. Cosmol. Astropart. Phys. 12 (2012) 014.

[37] D. Lenz, O. Doré, and G. Lagache, Large-scale maps of the cosmic infrared background from Planck, Astrophys. J. **883**, 75 (2019).

[38] N. Aghanim, M. Ashdown et al. (Planck Collaboration), Planck intermediate results. XLVIII. Disentangling Galactic dust emission and cosmic infrared background anisotropies, Astron. Astrophys. **596**, A109 (2016).

[39] M. Remazeilles, J. Delabrouille, and J.-F. Cardoso, Foreground component separation with generalized Internal Linear Combination, Mon. Not. R. Astron. Soc. **418**, 467 (2011).

[40] A. Maniyar, G. Lagache, M. Béthermin, and S. Ilić, Constraining cosmology with the cosmic microwave and infrared backgrounds correlation, Astron. Astrophys. **621**, A32 (2019).

[41] G. Chon, A. Challinor, S. Prunet, E. Hivon, and I. Szapudi, Fast estimation of polarization power spectra using correlation functions, Mon. Not. R. Astron. Soc. **350**, 914 (2004).

[42] http://www2.iap.fr/users/hivon/software/PolSpice/README.html

[43] P. A. R. Ade, N. Aghanim et al. (Planck Collaboration), Planck 2015 results. XIII. Cosmological parameters, Astron. Astrophys. **594**, A13 (2016).

[44] Y. Omori, R. Chown, G. Simard et al., A 2500 deg² CMB lensing map from combined South Pole Telescope and Planck data, Astrophys. J. **849**, 124 (2017).

[45] N. Aghanim, Y. Akrami et al. (Planck Collaboration), Planck 2018 results. III. High frequency instrument data processing and frequency maps, Astron. Astrophys. **641**, A3 (2020).

[46] V. Böhm, M. Schmittfull, and B. D. Sherwin, Bias to CMB lensing measurements from the bispectrum of large-scale structure, Phys. Rev. D **94**, 043519 (2016).

[47] A. Lewis and G. Pratten, Effect of lensing non-Gaussianity on the CMB power spectra, J. Cosmol. Astropart. Phys. 16 (2016) 003.

[48] P. A. R. Ade, N. Aghanim et al. (Planck Collaboration), Planck 2013 results. XXX. Cosmic infrared background measurements and implications for star formation, Astron. Astrophys. **571**, A30 (2014).

[49] https://mocks.cita.utoronto.ca/index.php/WebSky_Extra-galactic_CMB_Mocks

[50] G. Stein, M. A. Alvarez, and J. R. Bond, The mass-Peak Patch algorithm for fast generation of deep all-sky dark matter halo catalogues and its N-body validation, Mon. Not. R. Astron. Soc. **483**, 2236 (2019).

[51] G. Stein, M. A. Alvarez, J. R. Bond, A. van Engelen, and N. Battaglia, The Websky extragalactic CMB simulations, J. Cosmol. Astropart. Phys. 10 (2020) 012.

[52] T. Namikawa and R. Takahashi, Impact of nonlinear growth of the large-scale structure on CMB B-mode delensing, Phys. Rev. D **99**, 023530 (2019).

[53] P. A. R. Ade, Z. Ahmed et al. (BICEP2 and Keck Array Collaborations), BICEP2/Keck Array V: Measurements of B-mode polarization at degree angular scales and 150 GHz by the Keck Array, Astrophys. J. **811**, 126 (2015).

[54] N. Aghanim, M. Ashdown et al. (Planck Collaboration), Planck intermediate results. XLVI. Reduction of large-scale systematic effects in HFI polarization maps and estimation of the reionization optical depth, Astron. Astrophys. **596**, A107 (2016).

[55] A. Lewis and S. Bridle, Cosmological parameters from CMB and other data: A Monte-Carlo approach, Phys. Rev. D **66**, 103511 (2002).

[56] A. Schillaci, P. A. R. Ade, Z. Ahmed et al., Design and performance of the first BICEP Array receiver, J. Low Temp. Phys. **199**, 976 (2020).

[57] A. N. Bender, P. A. R. Ade, Z. Ahmed et al., Year two instrument status of the SPT-3G cosmic microwave background receiver, in *Millimeter, Submillimeter, and Far-Infrared Detectors and Instrumentation for Astronomy IX*, Vol. 10708 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (SPIE, Bellingham, 2018), p. 1070803.

[58] P. Ade, J. Aguirre et al. (The Simons Observatory Collaboration), The simons observatory: Science goals and forecasts, J. Cosmol. Astropart. Phys. 02 (2019) 056.

[59] K. N. Abazajian, P. Adshead et al. (CMB-S4 Collaboration), CMB-S4 science book, first edition, arXiv:1610.02743.

[60] K. Abazajian, G. E. Addison et al. (The CMB-S4 Collaboration), CMB-S4: Forecasting constraints on primordial gravitational waves, arXiv:2008.12619.

[61] M. A. Faúndez, K. Arnold, C. Baccigalupi et al., Measurement of the cosmic microwave background polarization lensing power spectrum from two years of POLARBEAR data, Astrophys. J. **893**, 85 (2020).

[62] B. D. Sherwin, A. van Engelen, N. Sehgal et al., Two-season Atacama Cosmology Telescope polarimeter lensing power spectrum, Phys. Rev. D **95**, 123529 (2017).

[63] A. Manzotti, Future cosmic microwave background delensing with galaxy surveys, Phys. Rev. D **97**, 043527 (2018).

[64] B. Yu, J. C. Hill, and B. D. Sherwin, Multitracer CMB delensing maps from Planck and WISE data, Phys. Rev. D **96**, 123511 (2017).

[65] U. Seljak and C. M. Hirata, Gravitational lensing as a contaminant of the gravity wave signal in the CMB, Phys. Rev. D **69**, 043005 (2004).

[66] J. Carron and A. Lewis, Maximum a posteriori CMB lensing reconstruction, Phys. Rev. D **96**, 063510 (2017).

[67] M. Millea, E. Anderes, and B. D. Wandelt, Bayesian delensing delight: Sampling-based inference of the primordial CMB and gravitational lensing, Phys. Rev. D **102**, 123542 (2020).

[68] J. Caldeira, W. L. K. Wu, B. Nord et al., DeepCMB: Lensing reconstruction of the cosmic microwave background with deep neural networks, Astron. Comput. **28**, 100307 (2019).

[69] F. Pérez and B. E. Granger, iPython: A system for interactive scientific computing, Comput. Sci. Engg. **9**, 21 (2007).
[70] J. D. Hunter, Matplotlib: A 2d graphics environment, Comput. Sci. Engg. **9**, 90 (2007).

[71] P. Virtanen, R. Gommers, T. E. Oliphant et al., SciPy 1.0: Fundamental algorithms for scientific computing in Python, Nat. Methods **17**, 261 (2020).

[72] A. Zonca, L. Singer, D. Lenz et al., HEALPY: Equal area pixelization and spherical harmonics transforms for data on

the sphere in Python, J. Open Source Software **4**, 1298 (2019).

[73] K. M. Górski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann, HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere, Astrophys. J. **622**, 759 (2005).https://cosmologist.info/lenspix/