BICEP/Keck. XVII. Line-of-sight Distortion Analysis: Estimates of Gravitational Lensing, Anisotropic Cosmic Birefringence, Patchy Reionization, and Systematic Errors

We present estimates of line-of-sight distortion fields derived from the 95 and 150 GHz data taken by BICEP2, BICEP3, and the Keck Array up to the 2018 observing season, leading to cosmological constraints and a study of instrumental and astrophysical systematics. Cosmological constraints are derived from three of the distortion fields concerning gravitational lensing from large-scale structure, polarization rotation from magnetic fields or an axion-like field, and the screening effect of patchy reionization. We measure an amplitude of the lensing power spectrum ALϕϕ=0.95±0.20 . We constrain polarization rotation, expressed as the coupling constant of a Chern–Simons electromagnetic term g a γ ≤ 2.6 × 10−2/H I , where H I is the inflationary Hubble parameter, and an amplitude of primordial magnetic fields smoothed over 1 Mpc B 1Mpc ≤ 6.6 nG at 95 GHz. We constrain the rms of optical depth fluctuations in a simple “crinkly surface” model of patchy reionization, finding A τ < 0.19 (2σ) for the coherence scale of L c = 100. We show that all of the distortion fields of the 95 and 150 GHz polarization maps are consistent with simulations including lensed ΛCDM, dust, and noise, with no evidence for instrumental systematics. In some cases, the EB and TB quadratic estimators presented here are more sensitive than our previous map-based null tests at identifying and rejecting spurious B-modes that might arise from instrumental effects. Finally, we verify that the standard deprojection filtering in the BICEP/Keck data processing is effective at removing temperature to polarization leakage.


INTRODUCTION
Even with many orders of magnitude improvement in the precision of measurements, primordial CMB fluctuations remain statistically isotropic, such that their statistics are well described by angular power spectra.On the other hand, multiple secondary effects after recombination distort the primary CMB fluctuations inducing new correlations among observed CMB fluctuations.Examples include gravitational lensing by large-scale structure (Zaldarriaga & Seljak 1998), patchy re-ionization that modulates the amplitude of the CMB fields (Hu 2000;Dvorkin et al. 2009), and cosmic birefringence that rotates the CMB polarization angle (Marsh 2016;Yadav et al. 2012).There are also various instrumental systematics that can generate spurious B-modes by distorting the incoming T, Q, U fields, most notably the temperature to polarization (T to P ) leakage caused by beam and gain mismatches (BICEP2 Collaboration III 2015; Keck Array and BICEP2 Collaborations XI 2019), and E to B leakage from errors in polarization angle calibration.A comprehensive investigation of the statistical properties of the temperature and polarization maps can be used as a powerful tool to distinguish the sources of the observed B-modes, deciding whether they are cosmological or instrumental.
The secondary and instrumental effects listed above are similar in that they can be described as distortion effects that mix the Stokes T , Q, and U fields along or around each line-of-sight direction n.Yadav et al. (2010) characterize distortions of the primordial CMB fluctuations by introducing 11 distortion fields which depend on the line-of-sight direction.The B-modes generated by these map distortions would have correlations with E or T that do not exist in the primordial signal in standard ΛCDM cosmology.Thus EB and TB correlations can be used to reconstruct the distortion fields and study the physical processes and instrumental systematic effects that are associated with specific types of distortions.
In this paper, we reconstruct the 11 distortion fields by applying the minimum variance EB and TB quadratic estimators derived in Yadav et al. (2010) and Hu & Okamoto (2002) to our observed B-mode signal and use their power spectra ĈDD L to constrain cosmological models and systematics.We will be referencing the previous publications from the BICEP/Keck (BK) experiments, (BICEP2 Collaboration I (2014), hereafter BK-I; BICEP2 Collaboration II (2014), hereafter BK-II; BICEP2 Collaboration III (2015), hereafter BK-III; Keck Array and BICEP2 Collaborations VII (2016), hereafter BK-VII; Keck Array and BICEP2 Collaborations VIII (2016), hereafter BK-VIII; Keck Array and BICEP2 Collaborations IX (2017), hereafter BK-IX; Keck Array and BICEP2 Collaborations X (2018), hereafter BK-X; BICEP/Keck Collaboration XIII (2021), hereafter BK-XIII).This paper is organized as follows: In Section 2, we give an overview of the different distortion fields and some background on the cosmological effects that correspond to some of the distortions.In Section 3, we describe our data and simulations used for the distortion field analysis.In Section 4, we outline the analysis method, including how to go from Q and U maps to an unbiased distortion field power spectrum and how to combine the distortion power spectra from two data sets.In Section 5, we use the power spectra of three of the reconstructed distortion fields to set constraints on gravitational lensing, patchy reionization, and cosmic birefringence.In Section 6, we discuss the instrumental effects that could produce distortion effects in our data and test for residual systematic effects in the BICEP/Keck data with the distortion field spectra.

INTRODUCTION TO THE DISTORTION FIELDS
In Hu et al. (2003), systematic effects in CMB polarization maps are described as modifications to the Stokes Q and U maps by distortions along the line-ofsight n.Following Yadav et al. (2010), we model these distortions with 11 distortion fields as where T , Q and Ũ stand for the un-distorted primordial CMB intensity and polarization fields.There are 11 terms; p(n) is a two-dimensional vector.Quantities in the top line correspond to distortions along a unique line of sight while the second line shows field mixing in the neighborhood of a single line n.The quantity σ denotes a chosen length scale for these terms and makes the distortion fields unitless.The operators ∂ 1 and ∂ 2 represent the covariant derivatives along the RA and Dec. directions, and Yadav et al. (2010) further show that these distortion fields can be estimated directly using quadratic combinations of the data.The filter weights f EB ℓ1,ℓ2 and f T B ℓ1,ℓ2 used in the construction of each of the eleven distortions from power spectra are shown in Table 2.Note that this table differs from a similar table in Yadav et al. (2010) in that we use a different notation for pixel-space and harmonic-space quantities, denoting the latter with alphabetical instead of numerical subscripts.Also, the weights to construct perturbations to E and B from these distortions have been omitted, since we do not use them.See Appendix A and Section 4 for more details.
Each of the distortion fields can be matched with a specific source, offering a rich phenomenology.The τ (n) field describes a modulation of the amplitude of the polarization maps, α(n) describes rotation of the polarization angle, f 1 (n) and f 2 (n) describe the coupling between the two polarization field spin states, the components of p(n), p 1 (n) and p 2 (n) describe the change in photon direction, γ 1 (n) and γ 2 (n) describe monopole T to P leakage, d 1 (n) and d 2 (n) describe dipole T to P leakage, and q(n) describes quadrupole T to P leakage.All 11 distortion fields can correspond to specific potential instrumental systematic effect.We will discuss them in depth in that context in Section 6.
Among the distortion fields in Eq. 1, there are three that correspond to known or conjectured cosmological signals.These are p(n): change of direction of the CMB photons, τ (n): amplitude modulation, and α(n): rotation of the plane of linear polarization.
CMB photons traveling from the last scattering surface are deflected by the intervening matter along the line of sight (Zaldarriaga & Seljak 1998).The change of photon direction, p(n), is referred to as the weak gravitational lensing of the CMB.
The lensing potential is commonly decomposed into gradient and curl lensing potentials, Φ and Ω (Hirata & Seljak 2003;Cooray et al. 2005;Namikawa et al. 2012), such that the lensed Q and U maps can be described as where the gradient ∇Φ has components ∂ i Φ and the curl ∇ × Ω has components ϵ ij ∂ j Ω, where ϵ ij is the antisymmetric symbol.To leading order, we obtain the map distortions which allows us to identify ∇ × Ω and ∇Φ with the curl and gradient mode of p(n), respectively.The gradient component of CMB lensing, Φ, is generated by the linear order density perturbations, while the curl component, Ω, is only generated by second order effects in scalar density perturbations or lensing by for example gravitational waves or cosmic strings (Dodelson et al. 2003;Cooray et al. 2005;Yamauchi et al. 2012).We expect these cosmological signals to be negligible (Hirata & Seljak 2003;Pratten & Lewis 2016;Fabbian et al. 2018).
The (gradient) CMB lensing potential power spectrum has been measured to high precision by many experiments using temperature, polarization, or both (Sherwin et al. 2017;Planck Collaboration et al. 2020a;Wu et al. 2019;Faúndez et al. 2020;Carron et al. 2022).
The distortion field τ (n) (amplitude modulation) can be generated by inhomogeneities in the reionization pro-cess, also referred to as patchy reionization.In addition to the kinematic Sunyaev-Zeldovich (kSZ) signal generated by the peculiar motion of ionized gas (Sunyaev & Zeldovich 1970, 1980), patchy reionization causes an uneven screening effect of photons (Dvorkin et al. 2009).The screening effect is described as where τ 0 (n) is the optical depth to recombination that varies for different line-of-sight directions n.Taylor expanding Eq. 4, the screening effect from patchy reionization generates the distortion field τ (n).
The details of the patchy reionization process are still largely unknown.Recent searches for the redshifted 21cm signal from neutral hydrogen by EDGES put a lower bound on the duration of reionization as ∆z ≳ 0.4 (Monsalve et al. 2017).The kSZ power obtained from the South Pole Telescope prefers ∆z ≲ 4.1 (Reichardt et al. 2021;Gorce et al. 2022).The constraints from Planck CMB temperature and polarization power spectra suggest that reionization occured at z re ≈ 8 with a duration of ∆z ≲ 2.8 (Planck Collaboration et al. 2016a).Previous work that studies patchy reionization through τ (n) reconstructions with CMB temperature and polarization include Gluscevic et al. (2013) and Namikawa (2018).We constrain the same crinkly-surface model of patchy reionization where the power spectrum of the optical-depth is given by with the amplitude A τ and the coherence length L c (Gluscevic et al. 2013).
The distortion field α(n) can be generated by a cosmic birefringence field that rotates the primordial Q and Ũ according to Two potential physical processes that can cause a rotation field α(n) are the coupling of CMB photons with pseudo-scalar fields through the Chern-Simons term, also described as parity-violating physics, and Faraday rotation of the CMB photons due to interactions with background magnetic fields.A massless axion-like pseudo-scalar field a that couples to the standard electromagnetic term has the Lagrangian density (Carroll et al. 1990): where g aγ is the coupling constant between the axionlike particles and photons, and F µν is the electromagnetic field tensor.The amount of rotation is given by: When the pseudo-scalar field fluctuates in space and time, the change of the field integrated over the photon trajectory, ∆a, varies across the sky and generates an anisotropic cosmic rotation field α(n).For a massless scalar field, the large-scale limit (L ≲ 100) of the expected cosmic rotation power spectra is described by (Caldwell et al. 2011) where H I is the inflationary Hubble parameter.
A second physical process that could generate a cosmic rotation field α(n) is Faraday rotation of the CMB photons by primordial magnetic fields (PMFs).In the large-scale limit (L ≲ 100), the cosmic rotation power spectra generated by a nearly scale-invariant PMF is (De et al. (2013); Yadav et al. (2012)): ) where ν is the observed CMB frequency, and B 1Mpc is the strength of the PMFs smoothed over 1 Mpc.The ν −2 frequency scaling of the Faraday rotation angle implies that a lower frequency offers better leverage for PMF measurement.

DATA AND SIMULATIONS
In this paper, we use BICEP/Keck maps that use data up to and including the 2018 observing season, referred to as the BK18 maps.In particular, we will focus on the two deepest maps: the 150 GHz map from BICEP2 and Keck Array data which achieves 2.8 µK-arcmin over an effective area of around 400 square degrees, and the 95 GHz map from BICEP3 which achieves 2.8 µKarcmin over an effective area of around 600 square degrees.These two data sets with the lowest noise levels are the most interesting for studying both the cosmological and instrumental effects related to the distortion fields.
We construct an apodization mask that down-weights the noisier regions of the T , Q, and U maps.For the polarized Q and U map, we use a smoothed inverse variance apodization mask similar to that used in BK-XIII.For the T map, we add constant power of 10 µK 2 to the smoothed noise variance and invert it to construct the apodization mask.The mask is similar to a Wiener filter with flat weights in the central region dominated by sample variance and an inverse variance weight at the edges of the map.Additionally, for the analysis of T to P distortions, we mask the 20 point sources with the largest polarized fluxes from a preliminary SPT-3G catalog by applying a 0.5 • wide Gaussian divot at the location of each point source in the apodization mask.The effects of point sources are discussed in Section 6.3 and Appendix E.
We reuse the standard sets of simulations described in BK-XIII and previous papers: lensed ΛCDM signalonly simulations constrained to the Planck T map (denoted by lensed-ΛCDM), sign-flip noise realizations, and Gaussian dust foreground simulations, each having 499 realizations.The details of the CMB signal and noise simulations are described in Section V of BK-I and the dust simulations are described in Section IV.A of BICEP2/ Keck and Planck Collaborations (2015) and Appendix E of Keck Array and BICEP2 Collaborations VI (2016).For estimating the noise bias of the distortion spectra constructed with TB estimators (described in Section 4.3), an additional set of lensed CMB signalonly simulations with unconstrained temperature are generated.
In addition to the standard simulation sets, we also generate simulations of random Gaussian realizations of the distortion fields, D(n), that are characterized by certain power spectra.For pipeline verification and calibration of the normalization factors (Eq.22), we use simulations described by a scale-invariant distortion spectrum, with fiducial amplitudes, A D fid. , and their specific values for each distortion field type given in Table 1.
For the amplitude modulation field τ (n), we generate Gaussian simulations of τ (n) according to the 10 −3 10 −4 10 −6 10 −4 10 −10 10 −14 10 2 Table 1.The fiducial amplitude for the scale-invariant power spectrum used as input for the Gaussian simulations for the calibration of the normalization.
power spectrum in Eq. 6.For comparing the sensitivity between quadratic estimators and BB power spectra for detecting distortion systematics, we generate Gaussian realizations of distortion fields with a scaleinvariant spectrum within a narrow range of multipoles (∆L = 50).
The distortion field simulations and the unconstrained temperature simulations are generated with the observation matrix R described in BK-VII.This matrix captures the entire map-making process including the observing strategy, the timestream filtering, and the deprojection of leading order beam systematics.Simulations are rapidly generated with matrix multiplications: where Q in and U in are input signal maps, Q noise & U noise are sign-flip noise realizations, and Q obs and U obs are "as observed" output maps.Because of T E correlation in ΛCDM such simulations are not fully accurate when the input T sky is not the same as that assumed in the construction of the deprojection operation which is built into the observing matrix.

Quadratic Estimator Construction
Since the BK-observed patch is relatively small (1-2% of the total sky), we work in the flat-sky limit using Fourier transforms.A complex field, D 1 (n) ± iD 2 (n), of spin s can be represented by its Fourier transform where ϕ L = cos −1 (n • L).In particular, we note that τ , ω and q are spin-0 fields, p 1 ± ip 2 and d 1 ± id 2 are spin-1 fields, γ 1 ± iγ 2 are spin-2 fields and f 1 ± if 2 are spin-4 fields.We transform between even-parity modes D a and odd-parity modes D b , and modes aligned with the RA/Dec coordinate system of the underlying maps, D 1 and D 2 , with a rotation for fields with even-valued spin or for fields with odd-valued spin.
In Appendix A, we show that one can construct unbiased minimum variance TB and EB quadratic estimators for each distortion field given by DXB where L = l 1 +l 2 , X may be T or E and C XX l1 , C BB l2 are the total observed power spectra including contributions from the noise and lensing.These estimators directly reconstruct the Fourier transform of the map distortions introduced in Eq. 1, which are denoted by alphabetical subscripts a, b.The specific filter functions f D,XB l1,l2 for each distortion field D and estimator XB are listed in Table 2. Eq. 20 shows the correction for the mean-field bias, which is estimated from simulations (Namikawa & Takahashi 2014a).
for the different distortion field estimators as introduced in Eq. 19.Here X can either be T or B in order to obtain the filter functions for the T B or EB estimator, respectively.CT X l 1 and CXE l 1 are lensed CMB power spectra corresponding to our fiducial model.Note that for the distortion fields pa and p b we use the notation prevailing in CMB lensing, Ω and Φ.
The analytical normalization factor is given by In practice, we obtain the normalization factor empirically by running Monte Carlo simulations: where D in L are the input distortion field Fourier modes, and D sim L are the un-normalized, reconstructed distortion modes.To obtain the input distortion Fourier modes D in L , the same apodization mask for the T , Q, and U maps are applied prior to the Fourier transform.We use the scale-invariant distortion input simulations (Eq.13) to calibrate the normalization factor for all the distortion fields except for lensing (p(n) in Eq. 1), where the standard lensed-ΛCDM simulations are used.
We can construct TB estimators sensitive to the distortion fields only concerning polarization (τ , α, f 1 , f 2 , p) due to the non-zero T E correlation in the CMB.Therefore all 11 distortion fields can be probed by both the EB and TB estimators with the weights listed in Table 2.However, the polarization-only distortion fields are better measured with the EB estimators, while the TB estimators have higher sensitivity to the distortion fields involving T to P leakage (γ 1/2 , d 1/2 , q).

Input E and B-modes for reconstruction
As described in BK-VII, the mixing of E and B-modes due to map filtering and apodization are taken care of by the matrix-based purification method with purification matrices Π B and Π E .The purified E and B-mode-only maps are: where Q obs and U obs can be either a simulation or the real map.The Fourier transform of the purified QE/B and Û E/B are used to construct the purified Ê and B modes which are the input to the quadratic estimators.
Eq. 19 minimizes the variance in the ideal case, ignoring beam smoothing, filtering, and apodization.In practice, transfer functions due to these effects must be compensated for in addition to E and B-mode purification.The observed BK maps, the Ê and B Fourier modes are corrected by: where C XX,in l and C XX,out l are the mean input and output spectra of the lensed-ΛCDM signal-only simulations, and t X l is the transfer function.
In practice, we find that using Fourier modes up to multipole of 600 yields the best signal-to-noise, whereas the l = 600 − 700 modes are noisy and can worsen the signal-to-noise of the reconstruction due to a misestimation of t X .Therefore, we use l max = 600 as our baseline reconstruction parameter.To avoid potential contamination by dust, we mask out the lowest multipoles and use l B min = 100 for 95 GHz and l B min = 150 for 150 GHz.The l B min cutoff is chosen such that the observed dust Bmode power spectrum in the BK patch of sky is lower than the lensing B-mode power spectrum at l > l B min using BK18 and Planck ΛCDM best-fit B-mode power spectra.

Estimating the Distortion Field Power Spectra
The power spectrum of a distortion field can be estimated by squaring the estimator D(L) from Eq. 19: where ĈDD L is the observed distortion field power spectrum and N DD L is the noise bias.When there is no distortion field present, the main contribution for N DD L is the disconnected N 0 bias, which can be estimated by the realization-dependent method described in Namikawa et al. (2013) and BK-VIII: where Ê and B are the real E and B modes, or a given simulation realization.The 499 simulation realizations are divided into two sets of roughly equal size, and the subscripts 1 and 2 stand for the first and second sets of simulations.The first term is averaged over the first set of simulations, while the second term is averaged over the first and second set of simulations.
For the TB estimators that we use for systematics checks in Sec.6, the realization-dependent bias is estimated in a similar manner to Eq. 28 but with T instead of E. Since the standard lensed-ΛCDM simulations are generated with the temperature sky fixed to the Planck T map (see BK-I), an additional set of simulations with unconstrained T are used as the simulation sets 1 and 2 to be averaged over.The realization-dependent bias is evaluated for the observed | D T B L | 2 and for each of the reconstructed distortion bandpowers of the 499 standard constrained-T simulations.See Sec.6.4 for more details.
When there is a distortion field signal, apart from the disconnected N 0 bias, there is an additional bias term that is proportional to the amplitude of the distortion field signal, referred to as the N 1 bias (Kesden et al. 2003).The N 1 bias can be estimated with two sets of simulations sharing the same distortion field realization (Story et al. 2015): where the subscripts 1,2 stand for the two sets of simulations with different CMB/noise realizations that share the same set of distortion field inputs, and ⟨ N 0 L ⟩ is the ensemble average of Eq. 28.
Higher-order bias terms are either mitigated by our choice of weights (Hanson et al. 2011) or are found to be small for our sensitivity levels (Beck et al. 2018;Böhm et al. 2018).Likewise, we do not expect a significant bias from galactic foregrounds in our polarization-based estimators (Beck et al. 2020) or from masking extragalactic sources in our temperature maps (Lembo et al. 2022).
Since the CMB signal contains gravitational lensing, the correlation between the lensing distortions and the various quadratic estimators can create a lensing bias N Lens .This is estimated with the mean reconstructed distortion field spectrum ⟨C DD L ⟩ of the lensed but otherwise un-distorted ΛCDM simulations.
We verified that after the normalization from Eq. 22 and accounting for the N 0 , N 1 , and N Lens , the input distortion spectra are recovered.Fig. 1 shows as an example the polarization rotation Ĉαα L spectra for a scale-invariant distortion input (Eq.13) and its N 0 , N 1 and N Lens biases.Since the N 1 bias is proportional to the distortion field spectra, it is included as part of the signal when we constrain the amplitudes of cosmological models with the distortion field power spectra.
We measure the distortion field power spectrum in multipole bins with widths of ∆L = 70, and the binned power spectrum values are referred to as bandpowers.The horizontal lines show the binned theory input for an ACB = 1 spectrum.Ĉαα L , the mean simulation bandpowers, matches the input spectrum after accounting for N 0 , N 1 , and the lensing bias N Lens .The error bars show the standard deviation of the simulation realizations.
The cosmological results in Section 5 are derived from L ∈ [1, 350) since the constraining power only comes from the low multipole modes, whereas in Section 6 we use L ∈ [1, 700) to perform systematics checks.In certain applications where the lowest multipole modes are important, i.e. constraining cosmic birefringence models and performing distortion systematics tests, an additional L ∈ [1, 20) bin is separated out from the L ∈ [1, 70) bin.

Joint analysis of two sets of maps
When using the distortion fields as systematics checks, the two frequency maps are examined independently since the BICEP3 (95 GHz) and BICEP2/Keck (150 GHz) maps may have different instrumental systematics.However, for studying cosmological signals, it is desirable to combine the results from the two frequencies into a single more powerful measurement.For the inference of cosmological information we will only consider the most sensitive EB estimators in the combination of the two frequency maps.
Our approach is to form distortion field estimators with all possible combinations of the E and B modes: DE1,B1 , DE1,B2 , DE2,B1 , and DE2,B2 , where 1,2 stands for 95 GHz and 150 GHz respectively.In our analysis for the tensor-to-scalar ratio r (BK-I, BK-X, BK-XIII), we examine all possible auto and cross spectra from the multiple frequencies and experiments without forming a combined map.In this distortion field analysis, we follow a similar approach where we construct all combinations of DEi,Bj , combine their auto and cross spectra, and derive a joint cosmological constraint.While the cross spectra approach might not necessarily yield the highest signal-to-noise compared to an analysis on the combined map, the different combinations of spectra can provide consistency checks between data sets.
In Eq. 27, the squares of the distortion field estimators are used to derive the auto power spectra.Similarly, we take the four estimators DE1,B1 , DE1,B2 , DE2,B1 , and DE2,B2 and compute all the possible auto-and cross-spectra.With 4 individual distortion field estimates, we get a total of 10 spectra: 4 auto-spectra and 6 cross-spectra.Similarly to Eq. 28, we compute the realization-dependent bias for the general scenario including the case of cross spectra.The more general form of the realization-dependent N 0 bias with the two E maps W/Y and the two B maps X/Z is (Eq.A17 of Namikawa & Takahashi (2014a)): where the subscripts 1,2 represent the two sets of simulations with different CMB/noise realizations, and D * is the complex conjugate of D. Ŵ , Ŷ can be either the observed 95 GHz or 150 GHz E modes, and X, Ẑ can be either the observed 95 GHz or 150 GHz B modes.It can be verified that Eq. 30 reduces to Eq. 28 when all four maps W, X, Y, Z come from the same frequency map.
The 10 possible auto and cross-spectra are combined linearly with appropriate weights so that the variance of the combined bandpowers is minimized, where i stands for the 10 spectra indices, b stands for the bins of the bandpowers, and the weights w b,i are a function of both the bins and the spectra indices.The weights are: where Cb,i is the mean power from the 499 simulations of spectrum i and bin b, and Cov −1 b,ij is the covariance ma-trix of the bandpowers of bin b from the 10 spectra.The minimum variance bandpowers C b from Eq. 31 combine the statistical power from the two frequency maps and are used to constraint the corresponding cosmological processes.

RESULTS: COSMOLOGY FROM DISTORTION FIELDS
In the following subsections, we present the observed distortion field spectra and the derived cosmological constraints from Φ(n), τ (n), and α(n) corresponding to gravitational lensing, patchy reionization and cosmic birefringence.

Gravitational Lensing
In Table 2, the weights for reconstructing the gradient, Φ, and the curl component, Ω, are listed.We use the gradient part to constrain the amplitude of the lensing signal parametrized as A ϕϕ L , while using the curl part as a systematics check in Section 6.
It is often more convenient to work with the lensingmass (convergence) field κ since the lensing potential has a red spectrum while the lensing-mass field has a nearly flat spectrum (Planck Collaboration et al. 2016b).The lensing convergence κ is related to the lensing potential Φ as: For the Fourier transform, we have: where L = |L|.Similarly, we also define the analogous quantity for the lensing rotation ω as: In Eq. 22 where we empirically calibrate the normalization factor, we correlate the input distortion field with the reconstruction.Before performing the Fourier transform and cross correlation, the inverse variance apodization masks are applied to the input distortion fields.Because the κ spectrum is relatively flat compared to the red Φ spectrum, it is better to apply the apodization mask to the κ map instead of to the Φ map to avoid mode mixing (BK-VIII): The 10 possible auto and cross-spectra for the lensing reconstruction are shown in Fig. 2, where the lensing convergence spectrum is plotted.The 4 diagonal subplots are the auto spectra from the 4 possible Φ Ei,Bj , while the other 6 are from cross correlating the different Φ Ei,Bj .The top left subplot is derived from only 95 GHz, and the bottom right subplot is derived from only 150 GHz.The other 8 subplots combine some information from both the 95 GHz and 150 GHz.
In Fig. 3 where Ĉb = Ĉκκ L is the observed lensing convergence bandpowers, C f b is the mean bandpowers from the lensed-ΛCDM simulations, and Cov bb ′ is the bandpower covariance matrix from the same lensed-ΛCDM simulations.
The best fit A ϕϕ L are: A ϕϕ L = 0.89 ± 0.23, for 95 GHz only, (38) A ϕϕ L = 0.95 ± 0.20, for 10 spectra combined.( 40) Since the lensing reconstruction is close to sample variance limited in the central parts of the map, the larger map coverage from BICEP3 95 GHz produces a tighter σ(A ϕϕ L ) = 0.23 compared to the 150 GHz σ(A ϕϕ L ) = 0.33.When all the cross spectra between the two frequencies are combined, we achieve σ(A ϕϕ L ) = 0.20, an ≈15%  (Planck Collaboration et al. 2014).The top left subplot is the auto spectrum from 95 GHz, while the bottom right is the auto spectrum from 150 GHz.The other subplots contain information from both 95 GHz and 150 GHz.We examine the 10 individual spectra, all 10 spectra combined, and some other data combinations, and we find that they are all consistent with the lensed-ΛCDM predictions.reduction compared to 95 GHz only.This is around a factor of 2 improvement from the previous BK-VIII result of A ϕϕ L = 1.15 ± 0.36.However, note that the lensing amplitude is better constrained by the B-mode power spectrum with A BB L = 1.03 +0.08 −0.09 in BK-XIII.We compile these constraints on the lensing amplitude in Fig. 4.  (Carron et al. 2022), the BK18 B-mode power spectrum measurement (BK-XIII) using the same data set and the previous measurement from the lensing potential reconstruction using BK14 data (BK-VIII).Bullet points denote constraints from the lensing potential auto-power spectrum, squares are used for measurements using CMB twopoint functions.

Patchy Reionization
Following Gluscevic et al. (2013) and Namikawa (2018), we use the τ (n) reconstruction from the EB estimator to constrain a simple crinkly-surface model.The model describes a scenario in which the universe goes suddenly from neutral to ionized but with a reionization surface that is crinkled on a comoving scale of R c ≈ 200Mpc (L c /150) −1 .The predicted power spectrum in Eq. 6 consists of white noise smoothed on an angular scale of θ C = π/L c .Fiducial model spectra are shown in Fig. 5 for L c = 100, 200, 400 and 800.The use of this parameter space is only valid in the assumption of this simplified model as it assumes an instantaneous reionization.As such the parameter L c has no physical meaning and limitations can be evaded with a more realistic model of the reionization history.The amplitude A τ is constrained with a log likelihood based on Hamimeche & Lewis (2008): where C f b are the mean bandpowers from the simulations of the fiducial model, g(x) = sign(x − 1) 2(x − ln x − 1), and Rb is the per bin ratio of the observed bandpowers over the fiducial bandpowers including the N 0 , N 1 , and lensing bias N Lens : Using the method outlined in Section 4.4, Fig. 5 shows the reconstructed Ĉττ L for 150 GHz auto-spectra, 95 GHz auto-spectra, and for all 10 auto-and cross-spectra combined.We see that the BK data are consistent with zero, offering no evidence for a patchy reionization signal-consistent with earlier limits derived from WMAP and Planck temperature maps (Gluscevic et al. 2013;Namikawa 2018).
We proceed to set upper limits on A τ in Eq. 6 using the log likelihood of Eq. 41.Fiducial simulations of L c at 100, 200, 400, and 800 are used to derive the constraints.In Table 3, the 2σ (95% C.L.) upper limits on A τ for 95 GHz only, 150 GHz only, and with all 10 spectra combined are listed.The sensitivity primarily comes from the 95 GHz map.Table 3.The 2σ upper limits for A τ (Lc) derived from the 95 GHz auto-spectrum, 150 GHz auto-spectrum, and the all 10 spectra combined.The corresponding Ĉττ L is shown in Fig. 5, and the fiducial model spectrum is Eq. 6.
Following Gluscevic et al. (2013) and Namikawa (2018), we plot the constraints of Table 3 in the A τ vs. L c parameter space.In Fig. 6, the constraint derived from our data is seen to be between the constraints from WMAP TT and Planck TT.The noise level of the τ reconstruction in the BICEP patch is roughly the same as the reconstruction from the Planck TT estimator.However, Planck's wider sky coverage significantly reduces the overall sample variance.According to Gluscevic et al. (2013), the lower limit on the duration of reionization ∆z ≳ 0.4 obtained by EDGES (Monsalve et al. 2017) can be translated to A τ ≳ 0.1, so only a narrow allowed band remains.

Cosmic Birefringence
The two physical processes that can lead to anisotropic cosmic birefringence, parity-violating physics and primordial magnetic fields, produce the predicted power spectra given in Eq. 11 and Eq.12.These are both of the form of L(L + 1)C αα L = constant.Following previous conventions (BK-IX, Namikawa et al. Planck TT WMAP TT WMAP TB EDGES BK 95GHz + 150GHz remaining parameter space Figure 6.Constraints on the amplitude A τ (Lc) of Eq. 6 obtained from this work (green) and in the previous work (Gluscevic et al. 2013;Namikawa 2018).The colored regions are excluded.
(2020); Bianchini et al. ( 2020)), we parametrize the power spectra with A CB , In standard BICEP/Keck analysis the overall polarization angle is adjusted to minimize the observed TB and EB power spectra.After this self-calibration, the polarization maps lose sensitivity to a uniform polarization rotation but are still sensitive to anisotropic rotations.

Constraints on parity violating physics
The best constraints from our data set on the coupling constant g aγ between axion-like particles and photons (Eq.11) are derived using the combined minimum variance Ĉαα L of the two frequency maps.Following the method outlined in Section 4.4, the reconstructed Ĉαα L for 95 GHz, 150 GHz, and all 10 spectra combined are shown in Fig. 7.

The Ĉαα
L in Fig. 7 are consistent with the un-rotated lensed-ΛCDM+dust+noise simulations.In a similar approach to Namikawa et al. (2020) and Bianchini et al. (2020), we use a log likelihood based on Hamimeche & Lewis (2008) to evaluate the 95% 2σ upper limit for A CB .This likelihood is the same as Eq.41 and 42 but with A CB in place of A τ .We obtain a 95% confidence upper limit of A CB ≤ 0.044.Using Eq. 11, this corresponds to an upper limit on the coupling constant g aγ , This is a factor of 3 improvement over our previous results using the BK14 maps in BK-IX: g aγ ≤ 7.2 × 10 −2 /H I .It is also somewhat better than the constraints from ACT: g aγ ≤ 4 × 10 −2 /H I (Namikawa et al. 2020) and SPT: g aγ ≤ 4 × 10 −2 /H I (Bianchini et al. 2020).

Constraints on Primordial Magnetic Fields
To derive constraints on PMFs, we study the two frequency maps separately since the polarization angle rotation from Faraday rotation scales with frequency as ν −2 (Eq.12).With the 95 GHz only and 150 GHz only spectra shown in Fig. 7, we again use the log likelihood Eq. 41 to derive 95% upper limits on A CB , where we obtain A CB ≤ 0.097 for 95 GHz and A CB ≤ 0.17 for 150 GHz.With Eq. 12, we convert these constraints on A CB to the following constraints on PMFs B 1Mpc , We show the previously published constraints from CMB four-point function measurements in Fig. 8, which are derived from 150 GHz maps, are SPT: B 1Mpc ≤ 17 nG (Bianchini et al. 2020) and BK-IX: B 1Mpc ≤ 30 nG.The leading constraint on this parameter is B 1Mpc ≤ 1.2 nG for a nearly scale-invariant PMF, derived from a combination of Planck and SPT two-point power spectra (Zucca et al. 2017).Through the effect of PMFs on the post-recombination ionization history, Paoletti et al. (2022)  The A CB constraint derived from BICEP3 95 GHz alone (A CB ≤ 0.097) is comparable to the upper limits from SPT and ACT, but the resulting constraint on B 1Mpc is considerably better because of the advantage of the lower frequency leverage with the ν −2 scaling.

Consistency checks and null tests
In this subsection, we discuss consistency checks and jackknife null tests for the three distortion fields τ (n), α(n), κ(n) that have been used to derive science constraints.We want to emphasize that the BK18 data set has passed a comprehensive set of data validations (BK-I, BK-III, BK-XIII), most importantly the jackknife null tests on the EE/BB power spectra.In the next section (Section 6), we study all the distortion effects in Eq. 1 as systematics checks, providing further evidence that the BK18 data set has systematic effects controlled below the level of statistical uncertainty.In this subsection we focus on demonstrating the robustness of the reconstructed distortion field spectra with different analysis choices, and present some additional null tests for the three distortion fields being used to derive science results.

Consistency checks
As consistency checks of the κ(n), τ (n) and α(n) reconstructions, the distortion bandpowers are constructed while altering some choices of the analysis.We summarize the conclusions here and provide detailed PTE values in Appendix C.
• Input E/B-mode multipoles: Similarly to BK-VIII and BK-IX, we lower the maximum multipole ℓ max from 600 to 400, raise the minimum multipole ℓ min to 200, or lower the B-mode maximum multipole ℓ B max from 600 to 350.The results from these three alternate choices are all consistent with the lensed-ΛCDM+dust+noise simulations.Additionally, the shift of the observed bandpowers from the three alternate choices vs. the baseline are also consistent with the shift in the simulation bandpowers for both 95 and 150 GHz and all three fields α, τ and κ.
• Differential beam ellipticity: In BK analysis the T to P leakage from differential gain and differential pointing are filtered out with the technique we call deprojection (BK-III).However, the T to P leakage from differential beam ellipticity cannot be treated with a direct filtering operation because the CMB TE correlation would cause a bias.Instead, we subtract a leakage template derived from the measured differential beam map ellipticity.We repeat the analysis without the subtraction and find very small changes in the reconstructed spectra, similar to what has been seen by Mirmelstein et al. (2021).
• Alternate foreground models: In Appendix D, the different foreground models explored in the main line BK18 analysis (BK-XIII) are used instead of the Gaussian dust simulations.We find that with the realization-dependent method and the baseline choice of ℓ B min = 100/150 for 95 GHz and 150 GHz, the shifts in the bandpowers when switching to the alternate foreground models, or to no foreground, are negligible.

Effects of absolute calibration error
Although the distortion fields κ(n), τ (n) and α(n) are dimension-less quantities, the EB quadratic estimator construction will lead to a distortion spectra C DD L that depends on the overall amplitude of the polarization map.An absolute calibration uncertainty of δ on the polarization map will translate to a systematic uncertainty of 4δ on either the A CB /A τ upper limits or the amplitude of the lensing potential A ϕϕ L (BK-VIII).The absolute calibration procedure which correlates the observed T with the Planck T map (BK-I) is estimated to have an uncertainty of 0.3%.The polarization efficiency is high (≈ 99%) with an uncertainty of ≲ 0.5% (BK-I, BK-II).Therefore, we estimate that the systematic uncertainty on A ϕϕ L from the absolute calibration of the polarized map is around 4δ ≲ 3%.

Jackknife null tests
We perform distortion field reconstruction on the 14 flavors of differenced (jackknife) maps that are designed to target different systematics in the main line analysis (BK-I, BK-III, BK-X, BK-XIII).The distortion field reconstruction is done with the full E-modes and the jackknife B-modes, written DE full ,B jack .We are interested in probing systematic effects that can potentially bias the distortion field reconstructions.For our B-mode search in the main analysis we are most worried about E-to-B leakage terms.Further, B-mode systematics that can be interpreted as a distortion field coupled with the full E-modes would be the most concerning contamination in terms of biasing the distortion field science results.Hence we focus on the particular combination of DE full ,B jack as opposed to DE jack ,B full .While the latter would be more sensitive if Eq. 1 would be a perfect model of our instrumental systematic contamination and systematic effects act symmetrically on E and B-modes, we decide to perform a more focused search of E-to-B leakage terms with the DE full ,B jack estimator.
With the method outlined in Section 4.3, we reconstruct the observed ĈDD where Cov bb ′ is the bandpower covariance matrix from 499 simulations, Ĉb is the observed distortion field bandpowers, and ⟨C b ⟩ is the mean bandpowers from the simulations.We also compute the χ and χ 2 values for each of the 499 simulation realizations and evaluate the probability-to-exceed (PTE) or p-value by counting the percentage of simulations that have larger χ or χ 2 .
There are 14 (jackknives) × 3 (fields) × 2 (frequency maps) = 84 PTE values for both χ and χ 2 statistics.These values are histogramed in Fig. 9.The value for χ PTE closest to zero or unity is 0.006 and the lowest χ 2 PTE is 0.008.Taking into account the look-elsewhere effect, we can construct a global statistical test that compares these real data values to the simulations.The specific procedure is as follows: χ extreme PTE: overall extreme PTE: p all = min(p χ , p χ 2 ) , where p in Eq. 49 are the 84 χ 2 PTE values and p in Eq.50 are the 84 χ PTE values for the real data, or a given simulation realization.The quantities p χ 2 and p χ are the most extreme χ 2 and χ PTEs, and the overall most extreme PTE p all is the smaller one of p χ and p χ 2 .We find that the most extreme value for the real bandpowers is p obs all = 0.006.Comparing p all between the observation and simulations, the probability to get a value smaller than the observed value is 0.59.Therefore, we conclude that there is no evidence of spurious B-modes in the κ(n), τ (n), and α(n) reconstructions from the jackknife maps.

RESULTS: DISTORTION FIELDS AS SYSTEMATICS TESTS
In this section, we comprehensively investigate the distortion fields potentially caused by systematics and consider different types of instrumental effects that could produce these distortions.In the main line analysis for tensor-to-scalar ratio r (BK-I, BK-VI, BK-X, BK-XIII), the most fundamental guard against systematics are the map jackknife (or null) tests.Maps are made splitting the data into (approximate) halves according to criteria which would be expected to result in nearly equal signal, but potentially different systematic contamination.The split maps are then differenced and the EE, BB and EB spectra of the result compared to simulations of signal plus noise.Well chosen jackknife splits can amplify systematics which cancel in the full coadd map (BK-III).It must be emphasized again that the published BK measurements of the tensor-to-scalar ratio, including the latest BK18 release, have passed all these null tests.
However, systematics detection and mitigation based on the distortion fields may complement and enhance the standard map jackknife tests in several ways.When using line-of-sight distortion fields for systematics we are checking for spurious B-modes in the full Q/U maps.Therefore, any spurious B-modes detected will indeed be present in the data set used to derive the science results.There are systematics that naturally cancel out with many detectors, or with the instrument boresight rotation of the observation strategy (BK-III Section 2.3 and Section 4).In this case failing a jackknife test might not necessarily mean that there was significant contamination in the full coadd maps.Conversely, there could hypothetically be systematic contamination in the full coadd map that somehow cancels in all considered jackknife splits.
In addition, going beyond the two-point statistics offers more information about the observed maps.Each of the distortion field corresponds to a certain type of systematics.Therefore, failing the systematics check for a certain distortion field can offer hints as to where to investigate.Furthermore, we will show that the quadratic estimators for distortion fields are usually more sensitive compared to the BB spectrum at detecting the corresponding distortions.Any spurious B-modes from distortion fields would be detected by its quadratic estimator before it significantly effects the BB spectrum.
In the map-making process, the Q and U modes that are potentially contaminated through beam systematics, in particular differential gain and differential pointing, are filtered out by the deprojection procedure (BK-I).The choice of deprojection time scale of around 10 hours is a compromise-a shorter deprojection time scale guards against systematics that vary over short periods, but at the same time removes more modes and reduces the overall statistical power (BK-I, BK-III).The distortion fields γ 1 , γ 2 are sensitive to the modes corresponding to the differential gain, while d 1 , d 2 are sensitive to the modes corresponding to differential pointing.Therefore, the distortion field estimators as systematics checks can guard against beam systematics that vary faster than the 10 hour deprojection time scale, eluding deprojection.However, there are many classes of systematic contamination that do not correspond to any of the distortion fields.Therefore, the distortion field systematics tests should be treated as a useful complementary check for the standard jackknife tests rather than a replacement.
For different experiments with different ways of measuring CMB polarization (e.g.pair differencing vs. rotating half wave plate), the mapping between detector systematics to the final line-of-sight distortion field can vary.We will discuss the case for experiments similar to BK which take the pair difference signal from pairs of detectors with orthogonal polarization directions, and then use boresight angle rotation to get a distribution of polarization angles to be able to solve for Q/U (BK-II).Hu et al. (2003) offers a more general discussion of the connection between instrumental systematics and distortion fields.

Instrumental Systematics and Distortion Fields
Since the telescope is constantly scanning the sky, a time-varying spurious systematic effect will translate into a position-dependent error in the map, which can be associated with different distortions (Hu et al. 2003;Yadav et al. 2010).On the other hand, a spurious systematic effect that stays constant in time but varies from detector to detector will also cause a positiondependent error in the map, since different detectors cover different regions of the map.
Miscalibration of detector gains (pair sum timestream signal) would be captured in the amplitude modulation field τ (n).The miscalibration can be time-varying, due to the uncertainties in the elevation-nod gain calibration (BK-II) between each hour of observation.There are also gain variations that stay constant in time but vary among detector pairs.When making the full season map, the observed T map is correlated with the Planck T map to derive one overall normalization factor to calibrate the amplitude of the map, referred to as the absolute calibration (BK-II).
Variations of the actual absolute calibration values between detectors can translate into spatial amplitude modulation of the coadded maps.This amplitude variation could also be introduced by bandpass mismatches (BK-II) between pairs of detectors.Since the detector gain is calibrated with the atmospheric response, and the atmospheric emission and CMB have different spectra, a mismatch in detector bandpasses will lead to a gain mismatch in the observed CMB signal.The gain mismatches discussed above can also happen between intra pair detectors.In this case, instead of an amplitude modulation distortion τ (n), we will get monopole T to P leakage (or differential gain leakage) which corresponds to the γ 1 (n) and γ 2 (n) fields.
Miscalibration of the orientation of the detectors will translate to the rotation of the plane of polarization field α(n).The overall rotation of the map is calibrated out by minimizing the EB and TB spectra (BK-I).However, variations of the orientation from detector to detector can translate to an anisotropic rotation distortion field whose amplitude can be limited by the quadratic reconstructions.
The f 1 (n) and f 2 (n) fields can be generated from a coupling between a gain miscalibration and the boresight angle rotation (for a telescope with such capability).For example, if at the boresight angles that contribute more to the Q map, the detectors consistently exhibit a higher gain than the boresight angles that contribute more to the U map, we will effectively get a higher amplitude Q map compared to U , which corresponds to an f 1 (n) distortion.With the observation strategy of BK, different pairs of detectors cover different RA and Dec. ranges on the sky.Therefore, gain variations among detector pairs can stochastically lead to f 1/2 (n) distortion fields.Although there is no physical mechanism known to us that can produce this type of systematic coupling between the detector gains and the boresight angles in the BK experiments, we constrain f 1 /f 2 for completeness.
The p(n) (κ(n)/ω(n)) fields capture changes in the CMB photon directions.The corresponding instrumental systematic is miscalibration of the beam center locations.In the main line analysis, the beam centers are derived from cross correlating the observed T maps with the Planck T map (Section 11.9 of BK-II).Any miscalibration or uncertainty that varies from detector pair to detector pair will produce the p(n) distortion fields.
The second line in Eq. 1 involves T to P leakage.These distortions arise from a mismatch of the beams of pairs of orthogonal detectors A and B. Consider a Gaussian beam: ] , where b is the beam offset, σ is the mean beam width, and e is the ellipticity in the direction of detector polar-ization (plus ellipticity).In the BK beam map measurements, the differential cross ellipticites are subdominant compared to the differential plus ellipticites (Keck Array and BICEP2 Collaborations XI 2019).Therefore, the cross ellipticity and its corresponding distortion field are ignored in this paper.A mismatch of the beam parameters between A and B will translate to the distortion fields as follows: In Table 4 we summarize the correspondence of instrumental systematics to the distortion fields.To demonstrate the connection between the systematic effects and the distortion fields, we generate simulations that contain some of the systematic effects and reconstruct their distortion field spectra.The systematic effects are added to the pair maps (BK-I) before the map coaddition step to reduce the cost of computation.The complete analysis is presented in Appendix B.Here we present two representative cases, one where the systematics are constant in time but vary over detectors (a 10 • random Gaussian detector polarization angle rotation), and another where the systematics vary over time (10% random Gaussian differential gain fluctuation varying from hour to hour).
The shift in the α, ω, EE, and BB spectra caused by the randomized detector rotation angles are shown in Fig. 10.For the extreme 10 • case simulated C αα L shows a very strong signal at low L while the BB spectrum remains unaffected.We stress that this is not anymore true if we didn't calibrate the overall rotation of the maps and there would be a non-zero mean angle calibration error, causing a scale dependent signal in the distortion field power spectrum Mirmelstein et al. (2021).The curl component of the lensing field ω also detects the rotation field due to the correlation between the α and ω estimators.
In Fig. 11, we show the spectra for γ 1 , γ 2 , EE, and BB for the 10% random differential gain fluctuations.Compared to the detector angle rotation where the instrumental effect is constant in time but varies over detector pairs, the time-varying gain mismatch simulations generate distortion power that are distributed over a wider range of multipoles.This systematic T to P leakage shows up strongly in both γ and BB.
In Appendix B, we observe that other kinds of systematics that are constant in time but vary over detectors  also generate distortions at large scales (low L).This is because each detector pair covers a significant portion of the map, and variations among detector pairs therefore primarily create large scale distortions.On the other hand, time-varying instrumental systematics can generate distortion power over a much wider range of multipoles.Depending on the type of systematic effects being studied, we can design systematics tests that focus on different multipole ranges of the distortion field spectra.

Quadratic Estimators vs. BB power spectra for detecting distortion fields
In the main line analysis (BK-I; BK-V; BK-VI; BK-X; BK-XIII) EE, BB and EB spectra of map difference splits are used to test for instrumental systematics.In this section, comparisons are made between BB power spectra and quadratic estimators in their ability to detect various systematics.We highlight cases in which the latter are more sensitive at detecting the spu-  4. A summary of the instrumental systematics that correspond to each distortion field in Eq. 1.
rious B-mode produced by the distortion fields.To that end, we generate simulations with Gaussian realizations of distortion fields within a narrow range of multipoles (∆L = 50).Any Gaussian distortion field with a smooth spectrum can be considered as a combination of multiple ∆L distortions, To make sensitivity comparisons between quadratic EB/T B estimators vs. BB power spectra, we use distortion field simulations with different L range inputs as the fiducial model and see how well the amplitude of that fiducial distortion spectra can be constrained by the standard simulations (un-distorted lensed-ΛCDM+dust+noise).We define the sensitivity ratio as: b stands for the mean bandpower from the distortion simulations characterized by Eq. 56, which we take as the "signal" of the particular systematic effect that we want to measure or constrain.σ(A D ) is the standard deviation of the best fit A D amplitude for the level of systematics from the 499 undistorted lensed-ΛCDM+dust+noise simulations.The estimator that is more sensitive in detecting the systematics will have a larger signal-to-noise in measuring the amplitude ÂD , and therefore a smaller σ(A D ) value.When the sensitivity ratio defined in Eq. 58 is greater than 1, the quadratic estimator is more sensitive than the BB power spectrum at detecting the particular distortion at that angular scale.
In Fig. 12, we demonstrate that the quadratic estimators are more sensitive than the BB power spectrum at detecting the distortion fields between L = 1 − 400.For all distortion fields, the quadratic estimators perform better when the distortion power is at larger scale (lower L).Among the polarization-only distortions, α, f 1 , and f 2 in particular are detected by the EB quadratic estimators with high sensitivity relative to the BB spectra.The distortions involving CMB temperature, d 1 , d 2 , and q are also very sensitively measured by the TB quadratic estimators.The above is as we would like it to be-we can detect systematics using the distortion fields before they significantly bias the BB spectrum.
Figure 12. sensitivity ratio as defined in Eq. 58.A sensitivity ratio σ(A EB D )/σ(A BB D ) above 1 means that the quadratic estimator is more sensitive than the BB spectrum at detecting the distortion field at that angular scale.The left plot shows the sensitivity ratio for the fields involving only the polarization, while the right plot shows the fields involving T .
Since we are reconstructing all the distortion field spectra C DD L simultaneously and the different distor-tion fields are not necessarily orthogonal to each other, it is important to study whether any spurious signal detected by a particular D 1 estimator can be reliably pointed to as an actual distortion signal from that field.To this end, we use the same set of simulations described by Eq. 56 to test for the cross sensitivity or correlation between the different distortion fields.We define the sensitivity ratio the same way as in Eq. 58, but in this analysis a different quadratic estimator D 2 is applied to try to detect the D 1 distortion input.
In Fig. 13, we show the correlation between different distortion fields using the L = 1 − 50 distortion simulations with the diagonal normalized to 1.The fact that the diagonal terms are much larger than the off-diagonal terms means the distortions would be much more strongly detected with the corresponding quadratic estimator before they are detected by another estimator.One exception is the correlation of q with d 1 /d 2 at low L. The existence of a large T to P dipole leakage can swamp the q estimator as we will see in Section 6.5.We see that q and d1/d2 have the strongest correlations, but in general all the distortion signals are best measured with their own estimators.

Effects of Point Source Contamination
The brightest point sources at the frequencies relevant for CMB observations are flat-spectrum radio sources (Battye et al. 2011).Since these point sources are brighter at lower frequencies relative to the CMB spectrum, the 95 GHz data set is much more affected by point sources in the distortion field analysis.
We estimate the effect of point source contamination in the distortion field analysis by injecting simulated point sources from a preliminary catalog obtained by private communication with the SPT-3G collaboration.The fluxes are taken from preliminary SPT-3G 95 GHz data, and are on average 2.5% polarized with an approximately exponential distribution.The T , Q, and U fluxes are converted to equivalent CMB temperature and added to the pixels closest to the location of the sources in the input map.The BICEP3beam and observation matrix R is then applied on to generate a point source simulation map (Eq.14).We find that the mean shift of the EB distortion spectra caused by these point sources is negligible.However, the 95 GHz TB estimators strongly detect the point sources, with the brightest few accounting for most of the contribution.
From this simulation we determine that the point source contribution becomes negligible after masking the 20 sources with the highest polarized fluxes.These are then added to the apodization mask by injecting Gaussian divots with 0.5 • width at the 20 locations.In Table 5, the χ/χ 2 PTEs for the real data, with and without the point source mask, for the TB estimators are listed.We find that the point source mask is necessary for the BICEP3 95 GHz data to pass the distortion field systematics tests, but does not affect the BICEP2/Keck 150 GHz data much.
In BK-XIII Appendix F it is estimated that the polarized flux from point sources may produce a bias on r at a level of ≈ 1 − 3 × 10 −3 .While small compared to our present uncertainties, point source contamination and its mitigation will become more important in future analysis, and the TB quadratic estimators can be a powerful diagnostic tool.In Appendix E, we discuss the reasons why TB quadratic estimators are sensitive to polarized point sources, derive estimators that are even more powerful for detecting point sources, and compare the performance of the different estimators at point source detection.

Distortion Field Systematics Tests on BK Real Data
With the connections between the various systematics and distortion fields established, in this section we present the results of the distortion field systematics tests for the two real data maps at 95 and 150 GHz.The distortion field systematics tests are performed with the  5.The χ/χ 2 PTEs with and without point source masks (psm) for the TB-reconstructed distortion fields.The bolded value is derived from the theoretical χ distribution since the observed value is outside of the 499 simulation distribution.The point source mask removes the brightest 20 sources in polarization according to a preliminary SPT-3G catalog.Without the point source mask, 95 GHz would fail the γ1 systematics test and in general have lower PTEs.For 150 GHz, there is no significant change to the PTEs.
same method as in Section 5.4.3 and Eq.47-48.The only difference is that the bandpowers Ĉb here correspond to the reconstructions from the full E and B map instead of the jackknife B map.In Fig. 14 we plot the difference of the real data reconstructed distortion field bandpowers and the mean of simulations, divided by the standard deviation of the simulations to show the significance of detection.
The C DD L spectra use the realization-dependent bias estimation outlined in Section 4.3.For the fields reconstructed with TB estimators (right side of Fig. 14), we substitute in the same observed T map for all the simulations.This is because the standard lensed-ΛCDMsimulations are constrained to the real CMB T map (BK-I), and T is measured with such high signalto-noise that the noise contribution is negligible.Since we are much more interested in testing for systematics in our B map rather than in the T map, we elect to fix to the same extremely well measured observed BK T map for both observation and simulations.
With the ĈDD L in Fig. 14, we evaluate the χ and χ 2 PTEs.We list these values in Table 6, and show histograms in Fig. 15.All the χ and χ 2 values lie within the 499 simulation distributions.There is one low PTE at 0.002 for the χ 2 of f 2 from 95 GHz.Examining the spectra in Fig. 14, the low χ 2 PTE can be traced to the second band power that fluctuates high.With the same method as Eq.49-51, we take into account the lookelsewhere effect and evaluate the global PTE statistic that compares the most extreme value among the 44 numbers in Table 6 to the simulations.We find that the probability to get a global value less than 0.002 is 0.08, offering no evidence for contamination in the data.6. χ and χ 2 PTEs derived from the distortion field spectra shown in Fig. 14.The results of TB estimators are derived with the point source mask applied, therefore the PTEs from d1 to q are identical to the "with psm" case of Table 5.The global PTE for the most extreme χ/χ 2 PTE (0.002 here) is 0.08.

Effectiveness of Deprojection
As discussed in Section 6.1, many of the distortion fields correspond to specific forms of beam mismatch.Beam systematics have been very important and well studied in the BK experiments (BK-III).Starting from BICEP2, the deprojection method has been developed to filter out potentially spurious signals that correspond to T to P leakage modes (BK-I, BK-VII).We have 95 GHz 150 GHz  L ).Plots on the left are reconstructed with the EB quadratic estimators, while the ones on the right are reconstructed with the TB estimators, and have the point source mask applied.The corresponding χ and χ 2 PTEs are listed in Table 6.also carried out extensive far field beam measurement campaigns every year as well as published a beam systematics paper to model and quantify how the beam systematics can affect the measurement of the tensor-toscalar ratio r (Keck Array and BICEP2 Collaborations XI 2019).
In Section 6.4, with the standard differential gain and differential pointing deprojection, the distortion field TB systematics tests pass with no evidence of any residual T to P leakage.We will now investigate whether the distortion field quadratic estimators can detect any spurious signal if we do not perform the differential gain and differential pointing deprojections.With the data products available on disk, it is simple to add the components that are filtered out by the deprojections back in, and construct maps without deprojection.We then perform the same χ/χ 2 distortion field systematics tests by comparing the observed ĈDD L with the simulations.
In Fig. 16, the BICEP2/Keck 150 GHz maps fail the d 1 , d 2 , and q systematics tests spectacularly without differential pointing deprojection.On the other hand, the BICEP3 95 GHz maps have much lower differential pointing and do not see much of a change in the reconstructed spectra when differential pointing depro-  6.
jection is turned off.Without differential gain however, we would detect a strong large scale γ 2 distortion in 95 GHz.We note that the differential pointing systematic in 150 GHz is also strongly detected by the q TB estimator.This is consistent with the results in Fig. 13, where we showed that q has significant correlation with d 1 and d 2 at large scale.
In Table 7, we show the χ and χ 2 PTEs for different deprojection options.In general, the χ and χ 2 PTEs decrease when either differential gain or differential pointing deprojection is disabled.This offers strong evidence that the deprojections are indeed successful in filtering out the monopole and dipole T to P leakage in the real data when compared to simulations with no such systematics.We note that the suite of jackknife tests for the power spectrum analysis includes targeted tests to detect these types of systematic contamination, which would lead to the same conclusion.

CONCLUSIONS
The line-of-sight distortion effects of the CMB can be characterized to first order with 11 fields.Three of these correspond to known or conjectured cosmological signals: gravitational lensing and κ(n), patchy reionization and τ (n), and cosmic birefringence and α(n).Combining the sensitivity from our two deepest maps: 150 GHz from BICEP2/Keck, and 95 GHz from BICEP3, we constrained physical models that can generate these distortion fields.For gravitational lensing, we measured the lensing amplitude to be A ϕϕ L = 0.97 ± 0.19, which is a factor of two improvement from our previous results in BK-VIII.For cosmic birefringence, we constrained the amplitude of cosmic birefringence and the related cosmological parameters to A CB ≤ 0.044, g aγ ≤ 2.6 × 10 −2 /H I , and B 1Mpc ≤ 6.6nG.This is a factor of three improvement on g aγ and a factor of four improvement for B 1Mpc compared to our previous BK-IX (BK14) analysis, resulting in the tightest constraint to date from the CMB four-point function.For patchy reionization, while not competitive compared to the Planck TT results (Namikawa 2018), we achieve the best constraint with CMB polarization on the τ (n) amplitude.
Treating the distortion fields as systematics tests, we demonstrated with simulations the connections between the distortion fields and the various instrumental effects in experiments with similar designs to BK.We further show that the EB/TB distortion field estimators are more sensitive than the BB spectrum at detecting random Gaussian realizations of distortions especially at larger scales.Additionally, we find that the TB estimators are very sensitive to contamination from polarized point sources.
We perform instrumental systematics tests on the 95 GHz and 150 GHz maps by comparing the distortion field bandpowers of the real data to lensed-ΛCDM+dust+noise simulations, and confirm that the 11 observed distortion spectra are consistent with the simulations.We also verify that the differential gain and differential pointing deprojections in our standard map-making pipeline are effective at filtering out the T to P leakage.Without the differential gain deprojection, we would detect an excess γ 2 power in the 95 GHz map, while without differential pointing deprojection, we would detect a very strong excess in the q, d 1 and d 2 fields of the 150 GHz map.This also confirms that the quadratic estimators are powerful tools to guard against T to P leakage systematics in the absence of differential gain and pointing deprojections.
With this first demonstration of quadratic estimators as instrumental systematics diagnostics on real data, we pave the way towards their future application as tools to self-calibrate upcoming data sets (Yadav et al. 2010;Williams et al. 2021).
The BICEP/Keck projects have been made possible through a series of grants from the National Science Foundation including 0742818, 0742592, 1044978, 1110087, 1145172, 1145143, 1145248, 1639040, 1638957, 1638978, & 1638970, and Yadav et al. (2010) and construct a minimum variance quadratic estimator of the distortion field from EB and TB correlations.See Eq. 1 for the definitions of the distortions.A scalar field such as CMB temperature T can be expanded in the Fourier basis as: A complex field (S 1 ± iS 2 )(n) of spin ±s can be expanded in the Fourier harmonics basis as: where ϕ l = cos −1 (n • l).We can also directly Fourier transform the S 1 ± iS 2 (n) fields as: One well known example is the transformation from Q, U in map space to the Q l , U l , E l , B l Fourier modes.In this case (S 1 ± iS 2 )(n) = (Q ± iU )(n) is a spin ±2 field.Its Fourier transform is (Q l ± iU l ) and its Fourier harmonics are (E l ± iB l ).The Fourier harmonics are not dependent on the coordinates, whereas the direct Fourier transform of the individual spin fields will transform into each other with a rotation of the coordinates.
For f 1 /f 2 , d 1 /d 2 , γ 1 /γ 2 , we reconstruct Fourier transform quantities which are coordinate dependent.They are used as systematics checks, and it is convenient to be able to connect them directly to the Q and U maps.For the lensing deflection p, we reconstruct the curl (Ω) and gradient (Φ) components, which are independent of the coordinates.
Assuming zero primordial B-mode, we write down to leading order the E L and B L with a distortion field D, where l 2 = L − l 1 and W B , W E are weights that can be derived for the individual distortion fields.Similarly, the E L and B L generated by the distortions that can be sourced by T to P leakage (γ 1/2 , d 1/2 , q) are written as: where the weights W B l1,l2 for generating B-modes are listed in Table 8.
Table 8.Weights and filters for the different distortion fields, where ϕ l = cos −1 (n • l) (Yadav et al. 2010).f D,XB l 1 ,l 2 are filter functions in Eq. 19, and the weight functions W B l 1 ,l 2 describe the B-modes generated from the distortions in Eq.A4 and A6.See Eq. 15 and Eq.A3 for the relation between the Fourier harmonic basis with subscript a/b and the Fourier transform of the distortion fields with subscript 1/2.
For the primordial un-distorted CMB fields, the power spectra are: where X, X ′ = T, E, B, and ⟨ Ẽl1 Bl2 ⟩ = 0, ⟨ Tl1 Bl2 ⟩ = 0. From Eq. A4 and Eq.A6, we can calculate the ensemble ⟨X l1 B l2 ⟩ correlation averaged over CMB realizations, where the filter functions f D,XB are listed in Table 2 for X = T /E, and L = l 1 + l 2 .There is only one universe and one CMB realization available for observation.However, given a L, there are many different combinations of l 1 and l 2 that satisfies l 1 + l 2 = L. Therefore, we can write down a linear combination of X l1 B l2 with some weight factor F D,XB l1,l2 that would minimize the variance of the DL ∝ estimator.The weight F D,XB l1,l2 can be derived from: where the bracket stands for the average over both CMB and distortion fields realizations.
For X = E or T , CXB l = 0.In this case, where f D,XB is exactly the factor in Eq.A9, and the C XX l1 , C BB l2 are the total observed power including contributions from the noise and the distortion fields (usually just lensing).Up to a normalization factor A D,XB L , the quadratic estimator for the distortion field can be written as: where L = l 1 + l 2 , and the analytical form for the normalization factor A D,XB L is: The mean-field bias ⟨ DXB L ⟩ is estimated from the simulations.After applying the correction for the mean-field bias, we have: DXB

SIMULATIONS OF SYSTEMATIC EFFECTS THAT GENERATE LINE-OF-SIGHT DISTORTIONS
In Section 6.1, we presented results for systematics simulations of random polarization and differential gain fluctuations.In this Appendix, we show the results from two more systematics simulations and offer a more in-depth discussion of each of the systematic effects.The four systematics simulations are: 1. 10 • random detector polarization angle rotation (Fig. 17   4. Dipole component of the T to P leakage from beam map simulations (Fig. 17(d)).
The well-designed BK observation strategy leads to a very high degree of cancellation of systematics with an increasing number of detectors and observing time.Since the purpose of these simulations are to clearly establish the connections between the different instrument systematics and the distortion fields, we inject large systematic errors in the polarization rotation, gain fluctuation, and differential gain.The simulations are for BICEP3 95 GHz except for the beam map simulations where we have simulations for both 95 GHz and 150 GHz.In each figure, the error bars represent the scatter on the mean spectra over 49 realizations of the systematics simulations.For clarity, we only show the distortion field spectra that are expected to detect the injected systematics.The distortion spectra that are not plotted do not show an elevated signal.
For the detector polarization angle errors, we generate misestimated angles by drawing from a Gaussian distribution of mean zero and standard deviation 10 • for each detector pair for the entire observing season.We rotate the detector polarization angle assumed in the map-making by these angles to produce a set of simulations including polarization angle systematics.The shift in the α, ω, EE, and BB spectra caused by the random detector rotation are shown in Fig. 17(a).Unsurprisingly, the random detector angle rotation creates a washout effect that reduces C EE ℓ .However,  the reduction caused by the washout effect and the distortion generated BB power roughly cancel, leaving the total C BB ℓ largely unchanged.For the distortion reconstruction spectra of the random rotation simulations, C αα L shows a strong signal at low L as expected.However, ω also detects the rotation field due to the strong correlation between α and ω estimators.
For the pair-averaged gain systematics, we simulate a 20% random Gaussian fluctuation on the pair gain.The injected relative gain error is constant over time and only varies from detector pair to detector pair.Due to the observation strategy, most of the detector pairs do not cover the same sky area at multiple boresight rotation angles.This means that a gain fluctuation over detector pairs also sporadically generates f 1 and f 2 distortions in addition to the amplitude modulation τ field.We show the f 1 , f 2 , τ , EE, and BB spectra generated by the gain fluctuation in Fig. 17(b).We observe a clear signal in f 1 and f 2 that corresponds to the injected systematics.However, we do not detect excess power in τ as one naively expects.One reason is that the τ EB estimator is not sensitive enough to detect the distortion effect at the level of 20% gain fluctuation with only 49 realizations of the systematics simulation.Another reason is the smooth apodization mask combined with the purification matrix which degrades the sensitivity to τ at the lowest multipole range where most of the τ distortion power is expected.
For the gain mismatches (differential gain) between detector pairs, we simulate a 10% random Gaussian fluctuation in (g A − g B )/2 for every detector pair and for every hour of observation.One possible contribution to gain mismatch is the uncertainties in the elevation-nod derived calibration factors (BK-I).With the differential gain deprojection that removes T to P leakage over 10 hour time scales, residual systematics can still arise from a differential gain that varies over shorter periods.Differential gain systematics lead to a monopole T to P leakage that corresponds to γ 1 and γ 2 .In Fig. 17(c), we show the spectra for γ 1 , γ 2 , EE, and BB.Compared to the detector angle rotation and gain variation simulations (Fig. 17(a) and (b)) where the instrumental effects are constant in time but vary over detector pairs, time-varying gain mismatches generate distortion power that is distributed over a larger range of multipoles.
For the instrumental systematics caused by the beam mismatch between orthogonal pairs of detectors, we make use of the T to P leakage template from the beam map simulations described in BK-III and BK-XI.With the high signal-to-noise far-field beam map measurements, the expected T to P leakage signal from the measured beam mismatch is simulated for both 95 GHz and 150 GHz.Since the beam map simulations are constant in time, the standard differential pointing deprojection completely removes the dipole leakage in the template and no signal is detected with d 1 and d 2 .As a sanity check, and to demonstrate the power of the quadratic estimators, in Fig. 17(d), we show the d 1 , d 2 , q, BB, and EE spectra generated by the dipole component of the leakage signal without deploying the differential pointing deprojection filter.Without deprojection, d 1 and d 2 spectra can detect the leakage signal at the lowest multipole with much higher signal-to-noise compared to the EE and BB spectra.In addition, q also detects the dipole leakage signal strongly because of its high correlation with d 1 and d 2 field.
In Table 9, we quantify the impact on the BB power from the four systematic simulations with an estimator ρ that represents the equivalent tensor-to-scalar ratio r level of the contamination (Keck Array and BICEP2 Collaborations XI 2019).The estimator is constructed in a similar way to the estimator for the lensing amplitude in Eq. 37, where Ĉb are the BB bandpowers from the systematics simulations, C r=1 b is the mean BB bandpowers for an r = 1 signal, and Cov bb ′ is the BB bandpower covariance matrix of the lensed-ΛCDM+dust+noise simulations.Even with the conservatively high levels of systematics in our simulations, the ρ estimates are relatively low at ≲ 1 × 10 −3 .Note that the large ρ value for T to P dipole leakage are for the case without differential pointing deprojection shown just for demonstration.In the main line analysis with deprojection enabled, the dipole T to P leakage will be entirely filtered out.

Systematics
sensitive distortion fields equivalent level of r (ρ) sensitivity ratio B3 10 • random polarization angle rotation α, ω < 5.4 × 10 −4 20 B3 20% pair averaged gain fluctuation f1, f2 < 4.7 × 10 −4 3.9 B3 10% time-varying differential gain γ1, γ2 < 1.2 × 10 −3 1.2 B3 95 GHz dipole T to P leakage (no deproj.)d1, d2, q 6.4 × 10 −3 2.6 150 GHz dipole T to P leakage (no deproj.)d1, d2, q 8.6 × 10 −2 4.2 Table 9.A summary of the four systematics simulations.The level of BB power from the systematics simulations are characterized by ρ, the equivalent level of r contamination.The sensitivity ratio shows the detection significance of the distortion field quadratic estimators vs. BB spectrum at measuring the systematics.A ratio larger than 1 means that the quadratic estimators are more sensitive to the systematic effect.Note that the dipole T to P leakage shown here is the case without the differential pointing deprojection.
From the distortion and BB bandpowers in Fig. 17, it is evident the relevant distortion spectra are able to detect the systematics with higher significance compared to the BB spectrum.Applying the same formalism described by Eq. 58 and using the mean systematics bandpowers as the fiduical C f b , we again evaluate the sensitivity ratio to compare the performance of quadratic estimators vs. BB spectra.When the sensitivity ratio is greater than 1, the quadratic estimator is more sensitive than the BB power spectra in detecting the systematics.In Table 9, we show the sensitivity ratio of the combined sensitivity of all relevant distortion spectra vs. BB spectra.The sensitivity ratio for the four systematics considered here are all larger than 1, which means that the quadratic estimators for distortion fields are more sensitive than BB at detecting the spurious B-modes from these systematics.The ratio for random parity profile.However, with the filter functions in Table 2 applied, the contribution from a polarized point source to γ 1 and γ 2 will be non zero.The power spectra of the quadratic TB estimators for γ 1 and γ 2 are effectively measuring the 4-point ⟨T QT Q⟩ and ⟨T U T U ⟩ correlations only using the B-mode and not the E-mode component of the CMB polarization signal.Compared to a direct correlation of the full ⟨T QT Q⟩ and ⟨T U T U ⟩, the γ 1/2 TB estimators will be more sensitive to the point sources because the sample variance is much lower without the contributions from the bright ΛCDME-modes.
The TB quadratic estimators with filter functions in Table 2 are designed to measure the distortion fields and not point sources.It is possible to design better estimators of the same quadratic form as Eq.19 that specifically target point sources.The point source estimators from the CMB temperature signal are described in Osborne et al. (2014), andNamikawa &Takahashi (2014b) extend the formalism to include polarization-only point source estimators.
Here we extend the point source estimators described in Section 3.1.2in Namikawa & Takahashi (2014b) to include temperature and polarization correlation.We will only consider the 1 source terms (see Section III of Osborne et al. (2014)), ignoring any contribution from clustering of the sources.
Let us consider a point source model with sky signal [T p (n), U p (n), Q p (n)] that is uncorrelated between pixels.Assuming the points sources are partially polarized with random orientations, we have: where where the bracket stands for mean over point source realizations given our point source model, and ⟨σ 1 (n)⟩ src = ⟨σ 2 n⟩ src .
Recall that a correlation of Eq.A9 would lead to a minimum variance quadratic estimator of Eq. 19.Therefore, we evaluate ⟨X l1 X ′ l2 ⟩ CMB with XX ′ = EE, EB, BB, T B. Note that we are taking the ensemble average over the CMB realizations here.With [E p l ± iB p l ] = d 2 ne −in•l [Q p ± iU p ](n)e ∓2iϕ l , we have: ).For the TB estimator, we have f T B,1 l1,l2 = − sin(2ϕ l2 ) and f T B,2 l1,l2 = cos(2ϕ l2 ) that probe ⟨T QT Q⟩ and ⟨T U T U ⟩ respectively.It is also possible to compute the cross spectrum of the two TB estimators to probe ⟨T QT U ⟩.Note that the two TB point source estimators for Sg 1 and Sg 2 have the same geometric terms (sin(2ϕ l2 ) and cos(2ϕ l2 )) as the distortion field estimators for γ 1 (n) and γ 2 (n).The only difference is that there is no longer the C T T ℓ factor in the filter function f .This change in the weights result in a factor of ∼ 2 sensitivity improvement at detecting point sources in the BICEP3 95 GHz maps.A sensitivity ratio of larger than 1 on the y-axis means the 4-point estimator is more sensitive than the BB power spectrum at detecting that population (polarized flux and fraction) of point sources.⟨EEEE⟩ has too low sensitivity to be probed by our simulations and is therefore not shown.
The polarized point sources also produce excess C BB ℓ power especially at higher ℓ.With the same method as Section 6.2, we compare the signal-to-noise of detecting the polarized point sources with the 4-point point source estimators vs. the 2-point C BB ℓ by running simulations.The parameter space explored are the polarized point source flux for individual sources (from 2 to 16 mJy) and the fraction of polarization (from 0.5% to 8%).For each set of parameters, we generated 99 realizations of 40 point sources with the same flux at random locations in the map and with random polarization orientations.In Fig. 19, we find that all the point source estimators (4-point functions) have sensitivity that scales as flux 2 relative to the 2-point function C BB ℓ , since the signal in a 4-point functions is ∝ flux 4 while the signal in a 2-point function is ∝ flux 2 .The ⟨BBBB⟩ estimator is much more sensitive than the ⟨EBEB⟩ and ⟨EEEE⟩ estimators, since ⟨BBBB⟩ does not have the large sample variance contribution from ΛCDM E-modes.Additionally, the ⟨T BT B⟩ estimator is more sensitive than the polarization only point source estimators when the polarization fraction is small.When the point sources are less than ≈ 3% polarized, the ⟨T BT B⟩ estimator is more sensitive compared to the ⟨BBBB⟩ estimator.

Figure 1 .
Figure 1.Example of the distortion field power spectra pipeline verification for the polarization rotation field α(n).The horizontal lines show the binned theory input for an ACB = 1 spectrum.Ĉαα L , the mean simulation bandpowers, matches the input spectrum after accounting for N 0 , N 1 , and the lensing bias N Lens .The error bars show the standard deviation of the simulation realizations.
we show the reconstructed Ĉκκ L of 95 GHz only, 150 GHz only, and with all 10 spectra combined.With the bandpowers in Fig. 3, we fit for the amplitude of the lensing potential power spectrum by taking a weighted mean of the real bandpowers over the fiducial simulation bandpowers (BK-VIII).With a linear model Ĉb = A ϕϕ L C f b of the noise-debiased power spectrum, where C f b is the fiducial model corresponding to the Planck ΛCDM prediction from their 2013 release 1 (Planck Collaboration et al. 2014), the least squares fit for A ϕϕ L is:

Figure 2 .
Figure2.The 10 different ways to combine two sets of E and B maps to form the lensing convergence spectrum estimator Ĉκκ L .The red line is the theoretical lensing convergence spectrum corresponding to our fiducial model assuming Planck 2013 cosmological parameters(Planck Collaboration et al. 2014).The top left subplot is the auto spectrum from 95 GHz, while the bottom right is the auto spectrum from 150 GHz.The other subplots contain information from both 95 GHz and 150 GHz.We examine the 10 individual spectra, all 10 spectra combined, and some other data combinations, and we find that they are all consistent with the lensed-ΛCDM predictions.

Figure 3 .
Figure3.The lensing convergence power spectrum C κκ L for 95 GHz auto-spectrum, 150 GHz auto-spectrum, and for all 10 auto-and cross-spectra combined.

Figure 5 .
Figure 5. Data and model power spectra of patchy reionization.The solid lines show fiducial model spectra from Eq. 6 for Lc = 100, 200, 400 and 800.The data points show Ĉττ L for 150 GHz only, 95 GHz only, and all 10 spectra combined.
2σ upper limit on A τ Lc 150 GHz 95 GHz all 10

Figure 7 .
Figure 7. Cosmic birefringence power spectra Ĉαα L for 95 GHz only, 150 GHz only, and all 10 spectra combined.The black line is the fiducial spectra from Eq. 43 with ACB = 1.Compared to Fig. 5 and Fig. 3, one additional bin of L ∈ [1, 20) is separated out from the L ∈ [1, 70) bin since the lowest multipoles are important for constraining ACB.

L
from DE full ,B jack and compare it to C DD L from the 499 lensed-ΛCDM+dust+noise simulations by evaluating χ and χ 2 values,

Figure 10 .
Figure10.Mean shift in α, ω, EE, and BB spectra in simulations with a 10 • random detector polarization angle scatter divided by the standard deviation of the non-rotated simulations.α and ω pick up a strong signal at low L, while C EE ℓ is suppressed at higher ℓ because the random per-detector polarization rotations average out near the center of the map and reduces the overall amplitude.

Figure 11 .
Figure 11.Mean shift in γ1, γ2, EE and BB spectra in simulations with a 10% Gaussian fluctuation of differential gain that varies from detector pair to detector pair and from hour to hour divided by the standard deviation of the nonfluctuated simulations.
where XX' represents EB, TB, or BB.C EB/T B b are the distortion field reconstruction bandpowers (4-point) from quadratic EB/TB estimators, and C BB b are the 2-point BB bandpowers.C f

Figure 13 .
Figure13.The correlation matrix of the sensitivity ratio as defined in Eq. 58 for all combinations of input distortion and quadratic estimators.The results plotted use distortion input from L = 1 − 50.The horizontal axis shows the input distortion injected in the simulations, and the vertical axis shows the quadratic estimators used to detect the signal.We see that q and d1/d2 have the strongest correlations, but in general all the distortion signals are best measured with their own estimators.

Figure 14 .
Figure 14.The fractional deviations of the 11 real data distortion field bandpowers ( ĈDD L − CDD L )/σ( ĈDDL).Plots on the left are reconstructed with the EB quadratic estimators, while the ones on the right are reconstructed with the TB estimators, and have the point source mask applied.The corresponding χ and χ 2 PTEs are listed in Table6.

Figure 16 .
Figure 16.The fractional deviations of the 11 real data distortion field bandpowers ( ĈDD L − CDD L )/σ( ĈDD L ) with and without the differential gain and differential pointing deprojections.95 GHz fails the γ2 systematics test without differential gain deprojection.150 GHz fails the d1, d2, and q systematics tests without differential pointing deprojection.
(a)10 % random polarization angle rotation (b) 20% random pair averaged gain fluctuation (c)10 % hour to hour differential gain fluctuation (d) Dipole T to P leakage from beam map simulations with differential pointing deprojection turned off.

Figure 17 .
Figure 17.The relevant distortion field, EE and BB difference spectra over the errorbar ∆C XX /σ(C XX ) from the different systematics simulations.Panels (a) and (c) are identical to Figs. 10 and 11 of Sec.6.1.

Figure 19 .
Figure19.The signal-to-noise for detecting point sources of the different 4-point estimators compared to 2-point C BB ℓ power in the BK18 95 GHz simulations.A sensitivity ratio of larger than 1 on the y-axis means the 4-point estimator is more sensitive than the BB power spectrum at detecting that population (polarized flux and fraction) of point sources.⟨EEEE⟩ has too low sensitivity to be probed by our simulations and is therefore not shown.
(Zucca et al. 20172020)he amplitude of the magnetic fields to √ B 2 < 0.69 nG.Comparing the constraint on the strength of primordial magnetic fields smoothed over 1 Mpc, B1Mpc, of this work with constraints from SPTpol anisotropic birefringence reconstruction(Bianchini et al. 2020)and Planck and SPTpol CMB power spectra(Zucca et al. 2017).

Table 7 .
χ and χ 2 PTEs derived from the distortion field spectra shown in Fig.16.The bolded numbers are derived from theoretical χ and χ 2 distributions when the real values are outside of the 499 sim distributions.
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 06-ARPA206-0040, 10-SAT10-0017, 12-SAT12-0031, 14-SAT14-0009 & 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 quasi-optical filtering was provided by UK STFC grant ST/N000706/1.The computations in this paper were run on the Odyssey/Cannon cluster supported by the FAS Science Division Research Computing Group at Harvard University.The analysis effort at Stanford and SLAC 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, Steffen Richter, Sam Harrison, Grantland Hall and Hans Boenish.We thank all those who have contributed past efforts to the BICEP/Keck series of experiments, including the BICEP1 team.We also thank the Planck and WMAP teams for the use of their data, and are grateful to the Planck team for helpful discussions.