Cosmology from CMB lensing and delensed <math display='inline'><mi>E</mi><mi>E</mi></math> power spectra using 2019

I. INTRODUCTION

Since the successful completion of the Planck satellite's sky surveys [1], the frontier in the study of the cosmic microwave background (CMB) anisotropies has been its polarization. Here, we present the strongest cosmological constraints to date derived from polarization data alone, constraints that are, by some measures, comparable to those from temperature data alone. Low-noise polarization measurements allow for new tests of the standard cosmological model, as the cosmological parameter inference relies on different signals, and they allow for additional robustness to many sources of potential systematic error.

The data we use come from the SPT-3G camera [2] installed on the South Pole Telescope (SPT) in early 2017. By October of 2024, SPT-3G had been used to survey 10;000 deg 2 , with the bulk of the Austral winter observing time spent on a 1500 deg 2 patch. The resulting data allow for powerful tests of the standard cosmological model, ΛCDM, via the primary and CMB lensing signals [3], improved constraints on primordial gravitational waves via delensing of the degree-scale observations of the BICEP/Keck Collaboration [4], and many other science applications [5][6][7][8][9][10][11].

This paper is the first of a series of papers on CMB power spectra, the CMB lensing power spectrum, and cosmological parameter estimates inferred from the 2019 and 2020 Austral winter observations of the 1500 deg 2 patch. The temperature and polarization maps used in these analyses are the result of 4 times as much observing time and twice as many detectors as were used for similar analyses that used a half season of data taken in 2018 with a partial SPT-3G focal plane [12][13][14][15].

This paper is distinct from the others in the coming series by its focus on polarization data only, and by its use of a novel algorithm for optimally and simultaneously estimating the CMB lensing power spectrum, the unlensedfoot_0 CMB E-mode polarization power spectrum (EE), and systematics parameters. This algorithm, the marginal unbiased score expansion (MUSE), was developed by Millea and Seljak [16] and Millea [17], with its precursors discussed in [18,19]. This is the first time the MUSE method has been applied to real data. We present the analysis in some detail and emphasize its advantages, including straightforward incorporation of systematic effects into a statistical model, a lensing potential map with lower noise than if it were produced by a quadratic estimator [e.g., [20] ], and an unlensed EE power spectrum with reduced sample variance relative to any inference of the lensed EE power spectrum [21].

Delensed power spectra from quadratic estimators have been estimated previously from data by Planck Collaboration et al. [22] and Han et al. [23]. The methodology for optimal lensing reconstruction has been principally developed in [16,[24][25][26][27][28][29]; and an application to SPTpol data was described in [30]. With the SPT-3G data we analyze here, we are crossing the threshold of polarization noise level below which optimal methods lead to significant improvements in lensing reconstruction.

We focus solely on polarization anisotropies in order to avoid modeling extragalactic foregrounds and their non-Gaussian properties. Such modeling would likely be necessary for a similar analysis which included temperature anisotropies, where the relative impact of foregrounds is larger. Due to the low noise and high angular resolution of the SPT-3G measurements, this initial conservative approach still leaves us with an inferred CMB lensing power spectrum with high signal-to-noise ratio. In the part of sky we observe, over the angular scales used in the analysis, galactic foregrounds are also negligibly small [12,31].

We compare with other CMB measurements in both temperature and polarization, including those from the Planck Collaboration [32] and the Atacama Cosmology Telescope (ACT) Collaboration [33][34][35][36]. We check for consistency of our results in the context of the ΛCDM model. We also determine cosmological parameters assuming the ΛCDM model and extensions from SPT data alone and in combination with these other datasets. Although the SPT observations we use here cover much less sky relative to Planck and ACT, the very low noise (4.9 μK-arcmin in coadded polarization maps) and arc-minute angular resolution enable measurements of the EE power spectrum and CMB lensing potential power spectrum (ϕϕ) that are the most precise measurements ever made at l > 2000 for EE and L > 350 for ϕϕ.

We discuss how the signals in our data allow for the determination of cosmological parameters assuming the ΛCDM model. The origins of these determinations have interesting differences compared to those from Planck CMB data. For example, the Hubble constant (H 0 ) inference depends heavily on the lensing information; without it, the error increases by more than 10 times. Consistency between the Planck and SPT inferences of H 0 places further constraints on attempts to solve the H 0 tension [37] with changes to cosmological models [38,39].

We pay particular attention to the implications of our data for the "negative neutrino mass" problem and its relation to excess lensing power pointed out by Craig et al. [40], and to the suggestion from Amon and Efstathiou [41] that the σ 8 tension is associated with nonlinear evolution, either due to mismodeling the ΛCDM predictions or by neglect of some new physics that modifies evolution on nonlinear scales.

The rest of this paper is constructed as follows. In Sec. II, we describe the process for making CMB maps from real and mock observations. We give an overview of MUSE and our data modeling in Sec. III. We present our pipeline and data validation processes, the blinding procedure, and results in Sec. IV. We show the resulting bandpowers and cosmological parameter constraints in Sec. V. We conclude in Sec. VI. Throughout this work, we use a series of abbreviations for external datasets tabulated in Table III.

II. DATA AND MAPMAKING

We use data taken with the SPT-3G camera during the 2019 and 2020 winter observing seasons. The 1500 deg 2 SPT-3G winter footprint extends from -42°to -70°in declination and -50°to 50°in right ascension. The footprint is divided into four subfields with declination centers of -44.75°; -52.25°; -59.75°, and -67.25°, and that allows for fitting different systematics parameters between patches that improve detector response consistency and linearity across the entire field. Observations are performed in twohour periods and map a single subfield at a time. During one observation, the telescope performs 72 constant-elevation (and thus constant-declination) scans with a scan speed of 1 deg =s in right ascension and an elevation step of 12.5 arcmin between scans. The subfields are observed 1034, 902, 772, and 578 times, respectively, across the two seasons, which produces approximately uniform map noise levels across subfields.

Each of the ∼16; 000 detectors in the SPT-3G camera samples the sky as the telescope scans, producing onedimensional "timestreams" that record sky brightness as a function of right ascension. The timestreams used in this analysis have on-sky sample rates between 0.3 arcmin=sample and 0.6 arcmin=sample depending on declination, corresponding to an effective Nyquist frequency well beyond our maximum multipole of interest.

The timestreams are low-pass filtered to reduce aliasing when later binning into maps, and high-pass filtered to remove large-scale instrumental and atmospheric noise. The high-pass filter is performed by deprojecting a set of Legendre polynomials up to 30th order and a set of sines and cosines up to an equivalent angular scale of l ≤ 300.

The filtering described above can produce undesired "ringing" features in the maps around the locations of point sources (a term we will loosely use here to describe both emissive galaxies and the Sunyaev-Zeldovich effect from galaxy clusters, the latter which can be spatially extended). To reduce these features, when fitting the deprojection coefficients, we ignore regions around point sources, specifically around all emissive sources with 150 GHz flux >6 mJy and all galaxy clusters with a signal-to-noise >10, for a total of 2655 objects. This procedure reduces ringing but does not mask the point sources, which instead happens at the map level later in the pipeline.

After filtering, we bin the timestream samples using inverse-variance weights into 0.5625 0 pixels given a Lambert projection centered on a right ascension of 0°a nd a declination of -59.55°. Weights are determined from the detector noise in the range 320 < l < 4000. The pixel size was chosen so that the 1500 deg 2 patch fits into maps with a number of pixels on each side that are powers of two, which maximizes the efficiency of FFTs. For this resolution, this yields 4096 × 8192 maps.

After maps have been made for each observation in the 2019 and 2020 winter observing seasons, we sum these maps into a full-depth "coadd" with a small number of data cuts. In addition to full-depth coadds, we also construct 500 "sign-flip" noise realizations by summing the observations with a random half of them multiplied by -1, canceling the signal but leaving the statistical properties of the noise unchanged. We use these to quantify the noise covariance, and to directly add to simulated signal maps to produce highly realistic simulations of the full dataset. We estimate that the noise realizations are correlated with each other at the percent level due to the finite number of observations in the dataset, but this is below the level which can significantly impact this analysis.

Next, we apply an isotropic antialiasing filter to the 0.5625 0 maps and rebin them to the final 2.25 0 resolution, 1024 × 2048 pixel analysis maps. We note that the reason for not binning the timestreams directly into the final 2.25 0 resolution is exactly to be able to apply this antialiasing filter, as otherwise nothing would prevent aliasing in the scan-perpendicular direction. The final Nyquist frequency for these maps corresponds to l ¼ 4800, high enough to accommodate the l max ¼ 4000 of our analysis while at the same time minimizing the pixel count and hence, the computational cost of MUSE.

Finally, we apply the pipeline masks described in Sec. III C. These processed and masked 1024 × 2048-pixel maps of Q and U stokes polarization at 95, 150, and 220 GHz are the basic data, d, which is input to our analysis pipeline.

III. METHODOLOGY: BAYESIAN INFERENCE WITH MUSE A. Bandpower and systematics inference

We perform a map-level simultaneous Bayesian inference of the gravitational lensing potential bandpowers, the unlensed CMB EE bandpowers and systematic parameters. This procedure is described in Millea and Seljak [16] and Millea [17], and summarized below with a few additions specific to this work. The central piece of the analysis is the joint posterior probability function, Pðf; ϕ; θjdÞ; ð1Þ where f refers to the unlensed CMB maps (here, just polarization), ϕ to the gravitational lensing potential, θ to the bandpower and systematics parameter vector, and d to the maps described in the previous section. The exact details of our modeling and construction of this function are given in the next section. Here, we first discuss how we use it to perform inference. The goal of our analysis is to marginalize over the "latent" f and ϕ variables to produce the marginal posterior on just the parameters of interest, θ,

The constraints on θ described by this function are formally optimal, satisfying our goal of an optimal analysis. This optimality is achieved because, implicitly, the marginalization sources information from moments of the data at all orders, unlike, e.g., quadratic estimators, which only use second-order moments. In particular, this means our results use information from all combinations of the data such as EE, EB, BB, EEE, EEB, etc… Various techniques exist for performing the difficult highdimensional integration in Eq. (2). A popular but computationally expensive choice is to use Hamiltonian Monte Carlo to yield the exact answer, up to sampling errors. MUSE instead performs an approximate but much faster marginalization, which yields a multivariate Gaussian approximation to the marginal posterior,

for posterior mean θ, covariance Σ MUSE , and unimportant constant, C. The Gaussian approximation is well-motivated due to the central limit theorem as long as many modes of f and ϕ contribute to the constraint on each θ, which is well-satisfied here.

The marginal posterior mean, θ, is given by solving the following equation for θ:

where s MAP i ðθ; dÞ is the gradient of the joint posterior with respect to parameter θ i , evaluated at the maximum a posteriori ("MAP") estimate of z ≡ ðf; ϕÞ, i.e.,

MAP i ðθ; dÞ ¼ d dθ i log Pðẑðd; θÞ; θjdÞj θ ; ð5Þ with ẑðd; θÞ ¼ arg max z Pðz; θjdÞ:

One way to understand this definition is to consider s MAP i as a data compression, which, for a given value of θ, takes as input our ∼ million-dimensional data vector and returns another vector of the same dimensionality as θ, which in our case is around a hundred bandpower and systematics parameters. The MUSE algorithm then proceeds by picking a starting guess for θ and compressing the data according to this, generating a suite of data simulations given the same θ and compressing them as well, and finally iteratively searching for the value of θ which make the data and simulations most similar, as compared via this data compression [Eq. ( 4)]. Intuitively, this yields the parameters θ, which are most likely to have actually generated the real data. More formally, Millea and Seljak [16] show that following this procedure leads to an asymptotically unbiased estimate of parameters, θ, for any Gaussian or non-Gaussian latent space, z. For a Gaussian latent space, the data compression is a sufficient statistic, and therefore, θ gives the exact marginal posterior mean. For non-Gaussian latent spaces such as the one here, the data compression can be considered "lossy," but this leads only to a possible loss of optimality, with the asymptotic bias remaining zero. For a lensing analysis similar to the one performed here, Millea and Seljak [16] bounds any loss of optimality to less than 10% of the statistical uncertainty.

The covariance, Σ MUSE , is computed from

where

and Σ prior is the covariance of any assumed prior. Note that this is not a pure Monte Carlo covariance of the MUSE estimate, which would be prohibitive computationally, but rather a cheaper Monte Carlo covariance of just s MAP , with analytic propagation to the resulting covariance in parameters. All derivatives in practice are computed with automatic differentiation, and Monte Carlo error stemming from J is reduced using a method described in Appendix C. One aforementioned aspect of MUSE, which is key for this analysis, is that it is unbiased even for "lossy" data compressions. This holds as long as the simulation-generating distribution, PðdjθÞ, is an accurate representation of the real data. There is not actually any requirement on the form of s MAP i for MUSE to remain asymptotically unbiased and for the covariance expressions to be valid. The default choice for s MAP i given in Eq. ( 5) is simply the one which makes MUSE exact for Gaussian problems.

Here, we keep this choice largely intact, but modify the joint posterior appearing in Eqs. ( 5) and (6), Pðf; ϕ; θjdÞ, to be slightly different from the likelihood function used to generate the simulations, PðdjθÞ. Normally, the two distributions would be related by

Instead, here we have one statistical model denoted as the posterior model, Pðf; ϕ; θjdÞ, and another called the simulation model, PðdjθÞ.

In the analysis, the posterior model is used in computationally costly high-dimensional maximizations, while the simulation model is only ever used to generate a few hundred samples from PðdjθÞ. This hints at the motivation for this choice of separate models, namely to allow defining a simulation model which fully meets our accuracy requirements, while keeping a simpler "surrogate" posterior model which is much faster to compute and makes the analysis more tractable. With an accurate simulation model, the COSMOLOGY FROM CMB LENSING AND DELENSED EE … PHYS. REV. D 111, 083534 (2025) MUSE analysis will remain unbiased, and we find little loss of optimality due to approximations in the posterior model. We note that our analysis can also be considered a form of "simulation-based inference" (SBI). This term is used to describe a class of algorithms which leverage the ability to generate samples from PðdjθÞ directly into the ability to perform inference and has been of increasing interest in various domains in cosmology [see, e.g., [42][43][44][45][46][47][48] ]. In typical high-dimensional SBI applications to-date, an important first step is to obtain a suitable data compression function to reduce the dimensionality of the problem, often by training a neural network on the generated samples themselves [49]. The compressed data is then used to infer parameters with something like approximate Bayesian computation (ABC) [50,51], or by training a neural network surrogate model for the likelihood which can be inverted to obtain the posterior [see [52] ].

MUSE is similar in that it also involves a data compression function, but rather than learned with a neural network, the generic and semianalytic s MAP i is used instead. Furthermore, rather than using something like ABC or training surrogate models for inference, we exploit our access to gradients and Jacobians through the posterior and simulation-generating function to derive semianalytic expressions for the posterior mean and covariance. In this way, MUSE can be considered one of the possible augmentations of SBI suggested by [52] which use additional tractable quantities from the simulator.

C. Map-level model summary

The posterior and simulation models are specified in Eqs. ( 12)- (19), with a summary given below and full details given in Appendix A.

Simulation model:

Posterior model:

where

Reading these models top-to-bottom describes how to produce a random sample from the model, with "∼" used to represent generation of a random sample from the distribution on the rhs and N ðμ; ΣÞ a multivariate Gaussian with mean and covariance μ and Σ.

To summarize, in the simulation model, we generate curved-sky unlensed CMB maps (f) and lensing potential maps (ϕ) from appropriate covariance operators (C), which are parametrized by the bandpower amplitudes which we will infer (A EE b and A ϕϕ b ). We apply the lensing operation (L) to produce lensed CMB maps and project to the flat-sky Lambert pixelization of the data (P). We apply an additional small deflection (G) to account for relativistic aberration due to our proper motion relative to the CMB rest frame. We apply the beam convolution operator [Bðβ n ; β ν pol Þ], which is parametrized by a set of beam eigenmode amplitudes (β n ) and the polarization fraction of the beam sidelobes at each frequency (β ν pol ). We rotate the global polarization angle of the maps [Rðψ pol Þ] by an angle ψ pol to model systematic errors in the overall polarization angle of our detectors. In each subfield, i, we scale the maps by A ν;i cal , to model the uncertain absolute calibration of the instrument. We apply the model transfer function (TF), and apply the pixel window function (PWF ). We add in temperature to polarization monopole leakage Q and U templates determined directly from our temperature observations (t Q and t U ) with subfield-dependent amplitudes (ϵ i Q and ϵ i U ), as well as add in a sample of the noise (n ν ) from one of the sign-flip noise realizations. Finally, we apply our chosen pixel, "trough," and Fourier masks (see Appendix A 9) to produce a simulated data vector, d ν .

The posterior model is nearly the same, with the only differences being that the covariances, C, assume flat-sky statistics, the noise is generated from a covariance and is unmasked, and aberration is ignored. For the posterior model, it is possible to calculate the posterior probability function, given in Eq. (20). Here, we also choose to include a "supersample" prior, PðϕÞ, in order to reduce the lensing mean-field (before its impact is automatically removed entirely in the process of obtaining the MUSE estimate).

The final inferred parameter vector is

containing 77 bandpower parameters and 47 systematics parameters. The MUSE estimate is a joint best-fit and a covariance between all of them. We note that the approximation of leaving the noise unmasked, while seemingly benign, is actually critical to obtaining a sufficiently fast posterior function since it avoids an additional large matrix inversion, while having little impact on the final results (see also Sec. 3.7 of [30]). We also note that the ability to directly use the sign-flip noise realizations in the simulation model is extremely advantageous. Without this, a Bayesian analysis would depend on an accurate statistical model of all moments of the noise (power spectrum, bispectrum, etc…) and would require both sampling and evaluating the likelihood of a noise map given such a model. Here, noise modeling is essentially a nonissue since MUSE requires only sampling the noise, and the sign-flip realizations provide these samples and empirically match the true noise distribution.

There are a few effects which we considered but ultimately chose not to include the model. The first effect is polarized foregrounds. The dominant expected foreground comes from polarized point sources. For the same point source flux cut used here, Dutcher et al. [12] construct a set of priors on the foreground amplitudes based on measurements of temperature foregrounds scaled to polarization with very conservative priors on the polarization fraction. Based on this, we find that the maximum possible contamination has a power 2 orders of magnitude lower than our noise spectra at all frequencies. That the foreground contribution is small is also consistent with [58], in which it is found that even for much deeper observations than ours, the bias to lensing from polarized foregrounds is fairly marginal. As mentioned, this finding that the foreground contribution is negligible in polarization for this analysis was key in our decision to focus only on polarization in this paper.

The second effect we ignore is the bias due to spatial correlation between the lensing potential and the location of our masked point source holes [59]. Using the Agora simulations [60], which include a model for these spatial correlations, we find that for our point source masking threshold, the maximum bias to C ϕϕ L would be ∼0.5%. However, Lembo et al. [59] also show that this maximum bias is highly suppressed in an optimally filtered ϕ map such as the MAP ϕ maps which enter our analysis. We thus conclude this effect is also negligible here.

End-to-end pipeline tests

To check that our pipeline is unbiased, we start with an end-to-end test of the pipeline on a set of 100 mock simulations. First, we compute MUSE estimates of bandpowers and systematic parameters on each mock simulation using all three 95 þ 150 þ 220 GHz polarization maps (hereafter, "all-frequency"). Then we repeat the same procedures on single-frequency polarization maps with the corresponding systematic parameters fixed to the all-frequency results (this is how we handle systematics in the single-frequency runs on real data as well).

In the left two panels of Fig. 1, we show the mean bandpowers over 100 mock simulations compared to the input truth. The PTE of the difference between the mean and input truth and its conversion to a number of σ is given in the legend. 5 No bias is detected at 3σ significance from the 100 mock simulations. The scatter of the mean bandpowers is at the level of 10% of the statistical uncertainty, as expected from a test with 100 realizations. In the right panel of Fig. 1, we also show the mean systematics parameters over 100 all-frequency mock simulations, here scaled by their 1σ uncertainty. These are also recovered with no significant bias. Single-frequency runs do not vary systematics parameters; hence, those estimates are not shown.

We next compare the empirical covariance of the 100 allfrequency estimates with the MUSE covariance given in Eq. (7). The left column of Fig. 2 compares the square root of the diagonal covariance entries. We see good consistency between the two except for A cal . This is because the mock simulations were made without varying systematic parameters within their prior, while the MUSE covariance takes the prior uncertainty of all parameters into account. This difference is particularly evident for prior-dominated parameters such as A cal .

In the right column of Fig. 2, we compare the first and second off-diagonal elements in the bandpower covariance matrices, this time fixing the systematics in Σ MUSE for a more like-to-like comparison. We see good agreement FIG. 1. Verification that our pipeline recovers unbiased bandpower and systematics parameters using 100 mock simulations. The colored lines are the average MUSE parameter estimates over simulations after subtracting the input simulation truth (in the final panel, we also divide by the error bars for one simulation due to the otherwise large dynamic range). The first two panels show ϕϕ and unlensed EE bandpower amplitude parameters (and also consider single-frequency runs), and the last panel groups together all the systematics parameters (see Sec. III C and Appendix A for description of the systematics considered). The PTEs of the colored lines relative to the expected scatter and 2-tailed numbers of σ are shown in the legends. To get the PTEs, we Monte Carlo sample the expected distribution of χ 2 values, accounting both for the look-elsewhere effect of four tests (three single-frequency runs and one all-frequency run) and the fact that we compute the expected scatter itself from 100 mocks. For reference, gray bands are the all-frequency 1 and 2σ errors (note that this is not the expected scatter of any of the lines, rather is a way to judge the size of any potential bias).

between the empirical covariance and the MUSE covariance. The anticorrelation in the first off-diagonal elements is caused by projection effects, masking, delensing, and the transfer function, where the projection effects explain about half of the anticorrelation. Here, the projection effects refers to our use of the flat-sky approximation in the posterior model. This effect has been found in the previous SPT-3G analysis [31].

In Fig. 3, we show the MUSE covariance with its diagonal entries rescaled to unity (i.e., the correlation matrix between all bandpowers and systematics parameters). We find no significant correlation between lensing bandpowers and any other parameter, except for A ϕϕ at the lowest L bin, which shows correlations with A EE at ≲10% level and with β pol at ∼12% level. This suggests a robust determination of the lensing spectrum with little uncertainty from systematic effects. The correlations within the unlensed EE block are due to masking, delensing, projection and correlation with systematics parameters. The EE bandpowers are significantly correlated with both calibration uncertainty and the sidelobe polarization fraction, both of which add appreciable systematic uncertainty to our final EE estimates.

We next obtain ΛCDM cosmological parameter estimates from the 100 all-frequency mock simulations. Unlike the main results which use Cobaya and CAMB to perform parameter estimation, here, for speed, we have trained and verified the accuracy of an emulator using Capse.jl [61], and sample using Hamiltonian Monte Carlo.

The one-dimensional cosmological parameter posterior distributions inferred from each mock simulation and the product over all of these posteriors is shown in Fig. 4. It is expected that the product over all posteriors should cover the input truth, which is found to be the case. The nondetection of any bias with 100 sims suggests that, at 68% confidence, we have bounded any possible pipelineinduced systematic error to be less than ≲0.1σ for each cosmological parameter.

To summarize, we have validated the pipeline using 100 highly realistic mock simulations. With 100 simulations, we are able to detect biases at the level of 10% of the statistical uncertainty. We do not detect any such biases to within 3σ of the expected scatter on the mean in either our estimates of ϕϕ lensing spectra, unlensed EE bandpowers, systematics, or in the cosmological parameters derived from them. We find good agreement of the bandpower covariance matrices between the one derived from mock simulations and the one calculated in MUSE using Eq. (7).

Intrafrequency null tests

Null tests involve splitting the full dataset into two parts and differencing them such that the signal cancels, then checking whether the remaining portion is consistent with FIG. 2. Verification of the accuracy of the all-frequency MUSE covariance, Σ MUSE . (Top left) Comparison of the square root of the diagonal elements between the MUSE covariance, σ MUSE ≡ ffiffiffiffiffiffiffiffi ffi diag p ðΣ MUSE Þ, and the empirical covariance from running the estimate on mock simulations. (Bottom left) The ratio of the two curves in the upper panel. The shaded band shows the 1 and 2σ standard error on σ empirical from the 100 simulations used. These are expected to match for fully likelihood-dominated parameters, and do so for all parameters except the partially priordominated A cal parameters. (Top and bottom right) Comparisons of the first and second off-diagonal elements of the MUSE and empirical correlation matrices. In this case, systematics are assumed fixed in the MUSE covariance, since they are fixed in the simulations.

the expected statistical properties of noise. Finding something outside of this expectation could indicate the presence of a systematic, which might then also be present in the full dataset coadd at some level. By canceling the signal and thus, the sample variance from the signal, these tests can probe systematics to levels well below even what could impact the final results, which are themselves limited by sample variance. Since temperature data does enter our analysis via leakage templates, we also present null TE spectra to more thoroughly vet the full temperature and polarization dataset. Further details of the null tests will be presented in a future work.

We consider six types of data splits: sun, moon, azimuth, year, scan, and wafer. The sun and moon tests are designed to check whether a portion of the dataset is significantly contaminated by the radiation from the sun or moon. In these tests, one part of the dataset comprises the data we took when the sun or moon was above the horizon, and the other part below the horizon. The azimuth test is designed to check whether a portion of the dataset is contaminated by thermal radiation from a building near the telescope. One part of the dataset comprises the data we took when the telescope's azimuth angle was within 90°of that of the building, and the other part azimuth angles outside that range. The year test is designed to check whether the data we took during the 2019 season are significantly different from the data we took in 2020. The scan test is designed to check whether the data we took when the telescope was scanning in one direction are significantly different from the data we took when the scanning was in the opposite direction. Lastly, the wafer test is designed to check whether the data from one half of the full detector array are significantly different from the data from the other half.

For the scan and wafer tests, we do not expect the null spectra to perfectly cancel the signal due to details of our map-making pipeline. For the wafer test, this is because we filter time-ordered data from different detector wafers differently. For the scan test, this is mainly because we do not correct for each individual detector's time constant effect when converting its time-ordered data into map pixel values. The time constant effects refers to the fact that each detector does not respond to changing radiation from the sky instantaneously. As a result, a map made from leftgoing scans contain a version of the CMB shifted to the left, and vice versa for right-going scans. For these two tests, we calculate the expected null spectra through simulations of the map-making pipeline and subtract it. This is computed assuming a fiducial cosmology, but because these expected null spectra are very small compared to the signal, it does not significantly imprint model dependence to these null tests or degrade their usefulness.

For each null test type, we divide each of the two halves of the data into 25 smaller parts and call them bundles. We create 25 null maps by differencing bundles from each half one-by-one, then compute all 300 cross spectra from these 25 null maps and bin them into the same bandpowers as the analysis. The average of these is what we call a null spectrum, and the standard error is our estimate of the expected scatter of the null spectrum.

The maps used for these tests are curved-sky projections of the same timestreams which enter the main analysis, with spectra computed with PolSpice [62]. For convenience, most of the tests have been computed using slightly different masks than the ones used to produce our main results. These tests, which we do not expect to depend significantly on the masking choice, use a pixel mask which includes slightly more sky, and a harmonic space mask, which zeros all a lm for l ∈ ½500; 680 and m ∈ ½350; 425. This harmonic mask removes contamination from a narrow-band signal around 1.1 Hz, which is similar, but not exactly the same as the "trough" mask used to remove the contamination in the MUSE pipeline. As the 1.1 Hz signal presents itself mainly in detector wafer and scan direction tests, we have recomputed these with the exact MUSE pixel and Fourier masking to verify the MUSE masking is sufficient to remove it. FIG. 4. Verification that our pipeline recovers unbiased cosmological parameters using 100 mock simulations. We infer ΛCDM parameters from each simulated all-frequency dataset and plot the one-dimensional marginalized posteriors above as different colored lines for each simulation. The shaded gray region in each panel is the product over all posteriors, which should cover the simulation truth, denoted as vertical black lines. Note that the smaller scatter for τ is because this parameter is mainly constrained by our prior. Across the other five parameters, the shaded region is consistent with the input truth with significance equivalent to 1.1σ, suggesting no detection of any bias at the level afforded by our 100 simulations.

After masking, computing null spectra, and subtracting expected contributions as necessary, we calculate a χ 2 with respect to the expected scatter, assuming no correlation between neighboring bandpowers. To compute a threshold probability-to-exceed (PTE) for the observed χ 2 , we Monte Carlo sample the expected distribution of χ 2 values, accounting both for the look-elsewhere effect of 54 tests and the fact that we compute the expected scatter itself from a finite 300 simulations. We consider the test passed if the PTE is >5% (corresponding to a <2σ deviation) with respect to this Monte Carlo distribution. This corresponds to a PTE >0.07% on any single individual test amongst our suite (this is roughly 5%=54 with an additional correction for the slightly noisy expected scatter).

A summary of all tests and their PTEs is given in Fig. 5. All tests pass at our threshold of 0.07%, and no obvious features are otherwise visible.

Agreement of systematics with alternative estimates

Our analysis infers systematic parameters jointly with bandpowers by incorporating systematics into a map-level model for the data. In this section, we compare our systematics estimates to independent determinations from alternate methods which have traditionally been used to estimate systematics in CMB analyses.

In this comparison, we consider the global polarization rotation angle, ψ pol , the amplitude of monopole temperature-to-polarization leakage, ϵ i Q and ϵ i U , and the calibration at each frequency, A ν;i cal . Details of how these are estimated with alternate methods are given in Appendix A. Briefly here, ψ pol can be estimated by assuming the sky EB spectrum is zero, and searching for the angle correction which achieves this in the data. ϵ i Q and ϵ i U can be estimated by assuming that the sky TQ and TU spectra are zero and finding the deprojection of leakage templates, which achieves this. A ν;i cal can be estimated taking ratios of power spectra to a known calibrated spectrum and averaging over very broad multipole ranges. Here, A 150 cal is fit via ratio to Planck, which performs its own absolute calibration, and A 95 cal and A 220 cal are fit from internal comparison to the Planck -calibrated SPT 150 GHz spectrum.

The source of MUSE estimates of these systematics is generally the same as described above, but happens implicitly and jointly inside the global MUSE fit, driven by our model definition (e.g., the model assumes the CMB EB spectrum is zero and allows a global rotation controlled by ψ pol to absorb any EB power observed in the data during the fit).

The best-fit systematic parameters inferred with MUSE from the SPT data are summarized in Table II. Broadly speaking, we find around 1% temperature to polarization leakage, a little less than half of a degree of global polarization angle miscalibration, and beam sidelobes (see Appendix A 4) which are polarized at around the 50% level.

In Fig. 7, we compare to systematic parameters derived from the alternative methods. Up to small differences expected due to slightly different choices of masking and filtering made in these alternative estimates, we find good agreement between our two sets of estimates. We also note that in general the MUSE systematic estimate is slightly tighter, likely due to simultaneously fitting for all systematics and increased optimality of the Bayesian approach.

C. Post-unblinding changes

With all of these validation tests passed, we unblinded results, fit cosmological models, and compared to data from other experiments. We found ΛCDM was a good fit to our own data and that our ϕϕ spectrum agreed very well with other experiments. However, the EE spectrum disagreed with the Planck theory prediction by a smooth slope (shown in the middle panel of Fig. 10), leading to an overall 3.5σ disagreement across five ΛCDM parameters (excluding τ) between Planck and SPT datasets. The slope was most significant at small scales and would have been undetectable by Planck, ACT, or previous SPT measurements, thus could have indicated a newly discovered departure from ΛCDM. We were thus motivated to exhaustively check the robustness of the slope. Doing so, however, we discovered that its likely origin was an overly strong assumption we had been making in our polarization beam modeling.

Prior to unblinding, we had assumed the polarized and temperature beams were the same at the scales used in the analysis. While the temperature beam sidelobes are mapped with high signal-to-noise via observations of bright point sources and planets, we lack sufficiently bright polarized sources to similarly map the sidelobe response of the polarized beams. To check the robustness of the EE slope to this assumption, we devised a model of the polarized beams to realistically account for how much these might differ from the temperature beams. Appendix A 4 describes the model in more detail, but in summary, we took the polarized beams to be the sum of two components. The first is a central main lobe which can be modeled physically using the basic optics of the telescope and is expected from physical modeling to be fully polarized and thus, the same

TABLE II. Posteriors on systematics parameters given our allfrequency MUSE run (top section) and priors which were input (bottom section). Blank entries or parameters not listed under priors have uniform priors on ð-∞; ∞Þ. Superscripts denote each of the four subfields. Note that no cosmological model is assumed here. Posterior 95 GHz 150 GHz 220 GHz A 1 cal 0.8977 AE 0.0074 0.9276 AE 0.0074 0.8612 AE 0.0084 A 2 cal 0.8928 AE 0.0072 0.9133 AE 0.0071 0.864 AE 0.0079 A 3 cal 0.8861 AE 0.0073 0.9399 AE 0.0075 0.849 AE 0.0083 A 4 cal 0.8775 AE 0.008 0.9272 AE 0.0081 0.8405 AE 0.0097 100ϵ 1 Q 0.291 AE 0.025 0.319 AE 0.021 0.492 AE 0.059 100ϵ 2 Q 0.402 AE 0.024 0.39 AE 0.022 0.418 AE 0.06 100ϵ 3 Q 0.555 AE 0.025 0.838 AE 0.021 2.096 AE 0.066 100ϵ 4 Q 0.603 AE 0.033 0.912 AE 0.03 2.221 AE 0.084 100ϵ 1 U 0.584 AE 0.027 0.74 AE 0.025 0.735 AE 0.064 100ϵ 2 U 0.648 AE 0.025 0.748 AE 0.023 0.642 AE 0.058 100ϵ 3 U 0.851 AE 0.027 1.238 AE 0.023 1.33 AE 0.063 100ϵ 4 U 0.83 AE 0.035 1.174 AE 0.03 1.121 AE 0.092 ψ pol ½° 0.393 AE 0.024 0.419 AE 0.021 -0.188 AE 0.079 β pol 0.44 AE 0.20 0.60 AE 0.28 0.51 AE 0.26 Prior 95 GHz 150 GHz 220 GHz A 1 cal 0.9275 AE 0.0094 A 2 cal 0.9095 AE 0.0089 A 3 cal 0.9386 AE 0.0096 A 4 cal 0.927 AE 0.011 β pol Uð0; 1Þ Uð0; 1Þ Uð0; 1Þ FIG. 7. Agreement of polarization angle rotation angle (ψ pol ), monopole T → P leakage (ϵ Q and ϵ U ) and calibration estimates (A cal ), between MUSE estimates and alternative methods which have been previously used. The groups of four points refer to the four different subfields in cases where those are estimated separately. Details of the alternative estimates are given in Appendix A. FIG. 8. No evidence for dipole, quadrupole, or octopole T → P leakage. These are MUSE posteriors from a run on all-frequency data that includes these extra leakage templates in addition to the other usual bandpower and systematics parameters. In our baseline run, we thus do not include such higher-order templates.

in polarization as in temperature. This model is fit to match temperature beams themselves at small scales where the model is a good description of the beam. The second component is a diffuse sidelobe stemming from effect such as scattering or reflection from of optical elements in the camera and scattering from panel gaps in the primary mirror, which is not necessarily expected to be fully polarized. We marginalize over the polarization fraction of the sidelobe component with a new parameter at each frequency, β ν pol . In this definition, our earlier assumption had corresponded to β ν pol ¼ 1, which makes the polarized beam identical to the temperature beam.

With these free parameters introduced, we reran the MUSE fit, finding β ν pol ¼ ½0.44 AE 0.20; 0.60 AE 0.28; 0.51 AE 0.26 at [95,150,220] GHz. That is, our CMB data alone are able to internally constrain the sidelobe polarization fractions, and it constrains them away from our previously assumed value of unity with mild significance. The ability to FIG. 9. Maximum a posteriori (MAP) maps from 95 þ 150 þ 220 GHz data at the MUSE estimate of theory and systematics parameters, θ. MAPs correspond to a filtering of the data which maximizes signal relative to noise (akin to a linear Wiener filter, but in the case of the MAP κ, a nonlinear filter). The statistical anisotropy visible particularly in the E map originates from this filtering. No other additional smoothing has been applied, although we have subtracted the simulation-mean of κ. Except for this subtraction, these MAPs correspond exactly to ẑ which enters the MUSE estimate in Eq. ( 5). We have chosen to show only unlensed E and lensed B (derived from the κ and unlensed E in the top two panels), as lensed and unlensed E look nearly identical at the scale of this figure, and unlensed B is effectively zero by assumption.

constrain these parameters at all arises from inter-frequency differences, as changing the β ν pol changes the beam at each frequency differently. Indeed, if we reexamine the inter-frequency agreement of our data with the β ν pol marginalized, we find slightly improved agreement. This is shown as the orange line in Fig. 6. Thus, even though our original blind inter-frequency tests did not show any statistically significant evidence for disagreement, there was indeed some small systematic disagreement hidden below the noise, which is alleviated by this improved beam modeling.

Regarding the slope in EE relative to Planck ΛCDM, the change from β ν pol ¼ 1 to the new marginalized posterior values largely removes it. It also increases our EE uncertainties at l > 2000 by as much as a factor of 2 for some bins, which, however, is highly correlated between bins, with the effect much like an additional relative calibration between high-l and low-l. The impact to our ϕϕ spectrum of this post-unblinding change is negligible, as optimal lensing estimation is largely insensitive to smooth changes to beams. Since much of our cosmological information is coming from lensing, it also means that our ΛCDM constraints from SPT data are not that strongly affected; for example, we originally found H 0 ¼ 67.3 AE 0.77 with the sidelobe polarization fraction fixed to unity and H 0 ¼ 66.81 AE 0.81 with it marginalized. A larger impact occurs in the ðn s ; Ω b h 2 Þ plane, with both parameters shifting just over 1σ and the posteriors broadening. These changes are all summarized in Fig. 10. Our final agreement with Planck across five ΛCDM parameters is 0.5σ, with the change roughly coming from equal parts central value shifts and posterior broadening in the ðn s ; Ω b h 2 Þ plane.

We consider this beam model sufficient to capture our polarized beam uncertainty because we find no significant shifts in results if we change the main lobe model to a simple Gaussian, suggesting the exact details of the main lobe modeling are not relevant. Additionally, allowing for a free linear dependence of the sidelobe polarization fraction as a function of angular distance from the peak response does not significantly shift results either, suggesting the model of a constant sidelobe polarization fraction is sufficient.

With this improvement to the polarized beam modeling, some additional small tweaks mentioned in Appenidx A which lead to insignificant changes, and all validation tests passing at the same or better levels than they were prior to unblinding, we adopted this as our final baseline model.

All results we will now present in the following section utilize this baseline model.

A. Bandpowers

In Fig. 11, we show our inference of ϕϕ and unlensed EE bandpowers. The ΛCDM model that best fits these data provides a good fit with χ 2 ¼ 75.8 given 71 degrees of freedom (76 bandpowers minus 5 cosmological parameters constrained beyond their prior), which corresponds to a PTE of 32.7%. By this test, the SPT data are consistent with the standard cosmological model. We also show as an orange band the 68% confidence level predictions of the ΛCDM model for the CMB lensing and unlensed EE power spectra given these SPT data. FIG. 10. Changes to our results made after unblinding, and their origin in changes to beams. In the left two panels, blue and orange data points show ϕϕ and unlensed EE bandpowers relative to Planck ΛCDM theory for our final result and our initial unblinded result, respectively. Colored lines in the EE panel show ratios of temperature (dashed lines) and polarization beams (solid lines) between final and initial unblinded versions. For the polarized beam changes, we also show the posterior uncertainty from the final polarized beam fit with shaded bands indicating 1σ. In the right panel, the same two colors show the corresponding changes to cosmological parameter results from the SPT dataset. TABLE III. Glossary of dataset abbreviations used in the text and in figure captions. When combining SPT with WMAP and Planck data, we sum the likelihood directly assuming no correlation, an assumption justified by the small sky area of the SPT main field (1500 deg 2 ). When combining ACT DR4 with WMAP and Planck data, we followed the procedure in [35]. When using both Planck NPIPE PR4 lensing and ACT DR6 lensing, we used the built-in likelihood implementation of act_dr6_lenslike package with variant=actplanck_baseline which accounts for correlations between ACT and Planck lensing due to overlapping sky coverage. a For the BAO dataset, we followed the proposal in Sec. III.3 of Adame et al. [65] to replace the DESI BGS and lowest redshift LRG measurements with SDSS MGS (z eff ∼ 0.15) and two SDSS DR12 LRG points (z eff ∼ 0.38, 0.51), as well as replacing the DESI Lyα point with joint eBOSS þ DESI Lyα BAO. With the exception of Sec. V I, any CMB dataset combinations we use includes a prior on τ of 0.051 AE 0.006 [64]. In Sec. V I, we instead drop the τ prior and PlanckT&E is to be understood as including the Planck low-lEE likelihood.

B. Origins of ΛCDM parameter constraints

The SPT bandpowers represent a significant advance in our capacity to constrain cosmological parameters with polarization data alone. In Fig. 13, we show constraints on four ΛCDM parameters from SPT, ACT, and Planck, assuming the ΛCDM model, and limiting to only EE spectra and polarization-based lensing reconstruction, or in the case of ACT, also including temperature in the lensing reconstruction. We see significantly tighter constraints on three of the parameters from SPT, while the higher dynamic range of the Planck EE data (which extends to lower l than shown in Fig. 12) leads to a tighter constraint on n s for Planck.

Constraints on ΛCDM parameters from SPT data and other dataset combinations, now including both temperature and polarization as detailed in Table III, are summarized in Table IV. The SPT constraints are comparable to full Planck results in the matter density, Ω m h 2 , and angular scale of the sound horizon, θ ⋆ . The constraining power on Ω m h 2 from SPT is mainly due to the CMB lensing bandpowers, which are sensitive to the amount of matter in the Universe. This is evidenced by the fact that adding just SPT lensing to WMAP data reduces the error on Ω m h 2 by over a factor of 3 from 0.0045 [72] to 0.0014, and that further adding SPT EE to this combination leaves the error unchanged.

That the sensitivity to Ω m h 2 primarily comes from the lensing reconstruction rather than primary power spectra (in this case unlensed EE) is a unique feature of these data compared to prior CMB datasets. The sensitivity of primary power spectra to matter density is weakened due to its weighting toward small scales that are sensitive to conditions at horizon crossing; at such early times the matter density is negligible compared to the radiation density. Conversely, the lensing power spectrum, if probed precisely enough, has high sensitivity to Ω m h 2 . This comes mainly (but not solely) from how CMB lensing power depends on the scale of matter-radiation equality. FIG. 13. Posterior distributions of ΛCDM cosmological parameters derived only from polarization EE spectra and polarization-based lensing reconstruction, or in the case of ACT [For the ACT DR4 EE likelihood, we adopted their default prior for yp2 ∼ ½0.5; 1.5], also including temperature in the lensing reconstruction. We see tighter constraints on Ω b h 2 , Ω c h 2 , and H 0 from the SPT data. We show the SH0ES constraint on H 0 for comparison. See also Table III for more details on dataset definitions. TABLE IV. Parameter estimates for various datasets assuming ΛCDM. The datasets are defined in Table III.

.045 AE 0.024 3.042 AE 0.011 3.046 AE 0.011 3.047 AE 0.012 3.042 AE 0.013 n s 0.944 AE 0.027 0.965 AE 0.004 0.9647 AE 0.0037 0.9637 AE 0.0076 0.961 AE 0.011 100θ MC 1.04162 AE 0.00066 1.04092 AE 0.00031 1.04104 AE 0.00028 1.04154 AE 0.00057 1.0394 AE 0.0022 Ω b h 2 0.02255 AE 0.00057 0.02237 AE 0.00014 0.02239 AE 0.00013 0.02248 AE 0.00027 0.02226 AE 0.00045 Ω c h 2 0.1227 AE 0.0018 0.1199 AE 0.0011 0.12039 AE 0.00094 0.1215 AE 0.0014 0.1210 AE 0.0014 τ reio 0.052 AE 0.006 0.0532 AE 0.0056 0.0547 AE 0.0058 0.052 AE 0.006 0.0515 AE 0.0059 H 0 [km=s=Mpc] 66.81 AE 0.81 67.4 AE 0.5 6 7 .28 AE 0.42 67.1 AE 0.6 6 6 .4 AE 1.1 100θ ⋆ 1.04181 AE 0.00069 1.0411 AE 0.0003 1.04125 AE 0.00027 1.04173 AE 0.00057 1.0396 AE 0.0022 Ω m h 2 0.1459 AE 0.0019 0.143 AE 0.001 0.14340 AE 0.00089 0.1446 AE 0.0015 0.1439 AE 0.0015 Ω m 0.327 AE 0.011 0.3147 AE 0.0068 0.3169 AE 0.0058 0.3211 AE 0.0085 0.326 AE 0.011 σ 8 0.8143 AE 0.0088 0.8104 AE 0.0051 0.8137 AE 0.0046 0.8177 AE 0.0065 0.8118 AE 0.0086 S 8 ¼ σ 8 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Ω m =0.3 p 0.850 AE 0.017 0.830 AE 0.012 0.8363 AE 0.0097 0.846 AE 0.016 0.847 AE 0.016

Additionally, the sensitivity increases at the small scales probed here, with C ϕϕ L ∝ ðΩ m h 2 Þ 2 at L ∼ 200 and asymptoting to ðΩ m h 2 Þ 3 by L ∼ 2000 [80].

The SPT data are also sensitive to the angular scale of the sound horizon, θ ⋆ . This sensitivity comes from the EE bandpowers, where the peak and trough locations of the EE spectrum are highly sensitive to θ ⋆ . The importance of the EE bandpowers for this constraint, as opposed to the lensing bandpowers, can be seen in the comparison of θ ⋆ constraints from SPT, WMAP þ SPT, and WMAP þ SPTϕϕ; the exclusion of SPT EE significantly degrades the θ ⋆ determination. We note that the constraint on θ ⋆ is known to be improved by delensing [81], a result which applies to our unlensed spectra inferences here as well, although we do not explicitly attempt to quantify exactly how much improvement we have achieved.

The constraints on τ from the SPT dataset are driven almost entirely by Planck low-l polarization data, included here as a τ prior. Although lensing can, in theory, help constrain τ via degeneracy breaking, the contribution at current noise levels is negligible. The constraint on logð10 10 A s Þ is similarly largely driven by Planck, since it is degenerate both with τ and the absolute calibration, the latter which also comes from Planck.

The SPT constraints on Ω b h 2 and n s are weaker than those inferred from Planck. Both n s and Ω b h 2 cause tilts of the CMB EE spectrum, since Ω b h 2 affects the photon diffusion damping at small angular scales and n s affects the shape of the primordial spectrum. Our observations lack the very largest scales probed by Planck, which, logarithmically, contain a significant range of scales useful for constraining a tilt. Additionally, at small scales our uncertainties are increased by polarized beam uncertainty, preventing what could otherwise be a slightly better constraint. Our final result is a 2.5% determination of the baryon density.

With the baryon density, matter density, and θ ⋆ determined, one has everything necessary, in ΛCDM, for calculating both H 0 and Ω m [e.g., [82] ]. Adding in A s and n s one can also calculate σ 8 and S 8 . We show constraints on these four derived parameters in Table IV.

G. Amplitude of nonlinear structure growth

One proposed phenomenological resolution to the S 8 tension is a suppression of the matter power spectrum at nonlinear scales, potentially due to baryonic feedback or dark matter structure growth which differs from the standard picture [41,95]. The CMB lensing reconstruction we present in this paper is the most precise measurement of the CMB lensing spectrum at the smaller angular scales where nonlinear effects are more important and can be used to further constrain the space of possible models which FIG. 16. Comparison of the ΛCDM constraints on Ω m inferred from different CMB observations (blue) with those from SNe observations (orange) and DESI BAO (green). The SNe observations include Pantheonþ [89], DES-Y5 [90], and Union 3 [91]. The biggest discrepancy is between DESI BAO and the DES-Y5 determination of Ω m , a 2.5σ difference. FIG. 17. Comparison of the ΛCDM constraints on S 8 inferred from different CMB observations (blue) with those from galaxy surveys of HSC Y3, KiDs-1000 and DES-Y3 (orange). The estimate from Planck þ ACT þ SPT differs with the one from HSC-Y3 at 2σ, the one from KiDs at 3.3σ, and with the one from DES-Y3 at 3σ.

resolve this tension. To do so, we introduce a new parameter, A CMB mod , equivalent to the A mod parameter of [41] under the Limber approximation, but here measured at CMB lensing scales and redshifts. We implement this by directly scaling the non-linear correction to the lineartheory CMB lensing spectrum,

where C ϕϕ L;lin and C ϕϕ L;nonlin are the linear and nonlinear CMB lensing spectra, the latter computed via our default Halofit. This modified lensing spectrum is then also propagated to changes in the lensed primary spectra.

As a small aside, we have verified that constraints on A CMB mod are insensitive to the choice of nonlinear and baryonic feedback models. Given that A CMB mod is the parameter most sensitive to the nonlinear prediction, our constraints on other cosmological parameters in this paper are insensitive as well.

Figure 18 shows results of fitting the ΛCDM þ A CMB mod model to various CMB datasets. We take the Planck primary CMB measurements as a baseline, which tightly limit the linear contribution to lensing, then explore the impact to constraints on the nonlinear contribution, A CMB mod , from adding lensing measurements from Planck, ACT, or the SPT results from this work. As expected from the spectra shown in Fig. 12, there are significant improvements in constraining the nonlinear contribution from Planck to ACT, and again from ACT to SPT. The Planck primary CMB together with SPT detect a non-zero A CMB mod at >3σ for the first time.

All measurements are consistent with the standard value A CMB mod ¼ 1, and the tightest combination of all lensing measurements yields A CMB mod ¼ 1.60 AE 0.39. The A mod inferred from DES-Y3 and KiDs is A mod ¼ 0.820 AE 0.042 [95]. This latter result measures the amplitude of nonlinear structure growth at redshift z ∼ 0.3 and scales of k ∼ 1 Mpc -1 , whereas our CMB lensing measurements are sensitive to higher redshifts and larger scales. Using a Fisher calculation based on our actual achieved error bars, in the top panel of Fig. 18, we have computed the range of ðk; zÞ which contribute most to our constraints on total lensing power, as well as to just the non-linear contribution to the lensing power, i.e., to A mod .

This plot, taken together with our A mod constraints, suggests (at ∼2σ) that if the S 8 tension is to be explained by some nonstandard physics suppressing the nonlinear matter power spectrum, these effects have not started at the higher redshift and larger scales of CMB lensing as compared to galaxy lensing, and evolve in the decade of separation in ðk; zÞ before suppressing the galaxy lensing power by ∼20% at later times and smaller scales. We note that with even better future measurements of small-scale CMB lensing, sensitivity to A CMB mod will push to smaller scales; e.g., significantly tightening up measurement of lensing power in a broad band centered on L ¼ 4000 will lead to sensitivity to galaxy weak lensing scales, although at z ≃ 1 to 2 rather than z ¼ 0.3.

H. Massive neutrinos

The ΛCDM model, as usually understood, includes a strong assumption about the spectrum of neutrino masses, namely that the sum is 0.058 eV. This particular value follows from assuming that the lightest neutrino mass is effectively zero, that normal ordering applies, and that the usual interpretation of atmospheric and solar neutrino oscillations is correct. Here, we relax this assumption and FIG. 18. Top panel: Redshifts and scales to which our lensing measurements are sensitive. Colored regions are iso-contours in the contribution to the square root of the Fisher information from logarithmic bins in k and z. The outermost contour corresponds to 10% of the peak. Red is for the Fisher information in the total lensing power, whereas blue is for the Fisher information just in A CMB mod (and is weighted more towards higher L). Under the Limber approximation, a given ðk; zÞ contributes to a single CMB lensing multipole, L, and dashed lines denote this mapping for scales of L ¼ 30 and L ¼ 4000. Bottom panel: Constraints on A CMB mod , defined as the amplitude of nonlinear contributions to the CMB lensing spectrum [Eq. ( 29)], from Planck primary CMB measurements (labeled PlanckT&E) in combination with lensing measurements from Planck, ACT DR6, and the SPT dataset presented here. The orange band shows the 68% and 95% constraints of A mod inferred from DES-Y3 and KiDs datasets in [95]. study constraints on the sum of neutrino masses, Σm ν . For a review of neutrino physics and its cosmological significance, see [96].

Relative to the case of massless neutrinos, neutrino rest masses increase the mean energy density and therefore, the expansion rate. To keep the angular size of the sound horizon fixed, the cosmological constant must then decrease, with a result that the expansion rate is boosted at z > 1 and decreased at z < 1 [97]. The increased expansion rate at z > 1 reduces growth on scales below the neutrino free-streaming length. Increasing mass also decreases the neutrino free-streaming length. On scales above the free-streaming length, neutrinos can cluster, approximately canceling the impact on growth of the increased expansion rate. For neutrino masses consistent with the data we consider, the free-streaming lengths are on very large scales so that the clustering of neutrinos only impacts the CMB lensing power spectrum at L ≲ 200 [98]. For this reason, nearly all of the constraining power of CMB data on neutrino masses comes from the changes to the expansion rate, rather than the changed free-streaming scale [99]. This is not to say that CMB lensing does not contribute to the constraints, as it is indeed sensitive to the reduction in growth that occurs as a result of the increased expansion rate.

In Table V, we show the constraints on Σm ν from various datasets. Combining the CMB observations of Planck, SPT, and ACT, we report the following constraint on neutrino mass:

Additionally including BAO measurements yields

which puts extreme pressure on the inverted hierarchy minimal mass sum of about 0.098 eV. These constraints on Σm ν should be interpreted with some caution. The need for that caution is well-illustrated by Fig. 3 in Loverde and Weiner [98]. Here, they show that, due to constraints from primary CMB data on θ s , ω bc , and r d (where ω bc ¼ ω b þ ω c ), and constraints from BAO data on Ω m and ω m r 2 d (where here ω m ¼ ω b þ ω c þ ω ν and likewise for Ω m ), the preference of the combined constraints is for Ω m < Ω bc ; i.e., for an (unphysical) negative neutrino mass. Further, the Pantheon þ SH0ES constraint in this plane is completely at odds with both the CMB and BAO (either DESI or SDSS) constraints. Loverde and Weiner [98] argue that the main cause of DESI BAO data leading to tighter constraints on Σm ν is the central value of their constraints in the Ω m -ω m r 2 d plane, rather than a decrease in the uncertainty in this plane.

Similar reason for caution is given by Craig et al. [40] and Green and Meyers [100]. They point out that there is some weak (2.5σ) evidence for an excess in lensing power present in CMB data, relative to ΛCDM predictions given CMB and BAO data. This is, again, opposite the expected influence of an increase to Σm ν from its minimal (ΛCDM) value of 0.058 eV. We turn to this question of excess lensing power in the next subsection.

I. Excess lensing power?

As mentioned, Craig et al. [40] pointed out that the lensing power inferred from the Planck lensing reconstruction and from the ACT DR6 lensing reconstruction are somewhat in excess of ΛCDM predictions given DESI BAO and primary CMB data. Here, we investigate this claim, first with L-independent lensing template parameters and then with the phenomenological "negative neutrino mass" parameter used in [40], Σ mν , which can be thought of as an L-dependent lensing template parameter. TABLE V. Summary of constraints on selected extensions to the ΛCDM model. The errors indicate the 68% confidence region, and the upper limits on the sum of neutrino masses are 95% upper limits.

þ0.014 -0.011 -0.0070 AE 0.0056 -0.0098 AE 0.0052 -0.0076 þ0.0053 -0.0047 0.0014 AE 0.0014 w 0 -1.42 þ0.29 -0.47 -1.48 þ0.20 -0.37 -1.52 þ0.18 -0.35 -1.52 þ0.19 -0.35 -1.075 AE 0.046 w 0 -1.10 þ0.50 -0.69 -1.18 þ0.46 -0.62 -1.16 þ0.48 -0.60 -1.18 þ0.48 -0.58 -0.65 AE 0.23 w a -1.1 þ0.6 -1.9 -1.06 þ0.61 -1.94 -1.18 þ0.52 -1.82 -1.21 þ0.57 -1.79 -1.15 þ0.71 -0.59 N eff 2.95 AE 0.33 2.89 AE 0.18 2.85 AE 0.16 2.66 AE 0.14 2.86 AE 0.13 Σm ν <0.38 eV <0.31 eV <0.17 eV <0.20 eV <0.075 eV N eff 2.95 AE 0.32 2.89 AE 0.19 2.84 AE 0.17 2.67 AE 0.14 2.83 AE 0.13 Σm ν <0.38 eV <0.31 eV <0.17 eV <0.21 eV <0.061 eV N eff 2.79 AE 0.52 2.8 AE 0.3 2 .77 AE 0.26 2.60 AE 0.23 2.89 AE 0.23 Y P 0.256 AE 0.028 0.248 AE 0.018 0.250 AE 0.016 0.246 AE 0.014 0.241 AE 0.015

To conduct our analysis, we define three L-independent lensing template parameters:

(1) A recon scales the ΛCDM model lensing power used to predict the reconstructed lensing power, (2) A 2pt scales the model lensing power used to compute any lensed power spectra (2-point correlation functions), and (3) A lens scales both of these model power spectra; i.e., it is what we call A recon and A 2pt when we force them to be equal to each other. Note that what we call A 2pt has often been called A L .

We consider four model spaces, each an extension of ΛCDM. They are ΛCDMþA recon , ΛCDM þ A 2pt , ΛCDMþA lens , and ΛCDM þ A recon þ A 2pt . Taking θ to be the ΛCDM parameters, the extensions to the model spaces that include these template parameters are done via

where X ¼ recon, 2pt, or lens and C ϕϕ L ðθÞ on the right-hand side is the ΛCDM model evaluated at θ. 7 To assist in comparison with [40,100], in this subsection only we take "PlanckT&E" to include the Planck low-lEE likelihood and remove the τ prior, and we restrict our use of ACT data to ACT DR6 lensing only. Our choice of BAO data can be found in Table III.

In the dash-dotted curves of the top panel of Fig. 19, we see the well-known "A lens" or "A L " anomaly, which in our language is an A 2pt anomaly driven by Planck PR3TT data [101]. This is a preference, in the lensed TT, TE, and EE spectra, for excess lensing power beyond what one gets from the ΛCDM model. In the lower left panel, we can see that if both A recon and A 2pt are allowed to vary freely, the 95% confidence regions for the four considered data combinations have central values greater than 1 in both dimensions and no overlap with the ΛCDM point of

Turning to the question of an excess of reconstructed lensing power, we see in the solid lines of the right panel that if we marginalize over A 2pt (solid lines) then there is some preference for A recon > 1. The significances of these preferences are shown as ranging from 1.8σ to 2.6σ. From these significance levels and from the blue, orange, and green curves we see that we get a very similar result whether we use Planck, ACT DR6, or SPT lensing reconstructions, with the ACT DR6 and SPT curves nearly overlapping. If we remove DESI BAO data 8 (a case not shown in the figure), the preferences for A recon > 1 persist, although with somewhat reduced significance ranging from 1.4σ to 2.2σ. The preference for A recon > 1 is weakened much more if we fix A 2pt to 1 rather than marginalize over it (dashed curves); any evidence for excess reconstructed lensing power is then quite weak.

This dependence of the A recon posteriors on the treatment of A 2pt (fixing to 1 or marginalizing) makes sense given the correlation we see in the lower left panel. The lower the assumed A 2pt , the lower the distribution of A recon values. Physically, this correlation emerges because as A 2pt is increased, the ΛCDM parameters used to fit the Planck T&E data adjust in a way that decreases lensing power (both Ω m h 2 and A s decrease). With the ΛCDM lensing power decreased, a given level of reconstructed CMB lensing power needs a larger A recon to fit it.

Readers who are familiar with [100] may be surprised at this sensitivity of A recon to choice of treatment of A 2pt (marginalizing it or fixing it to 1), given arguments made in that paper about the unimportance of the 2-point lensing anomaly to their conclusion of excess lensing power. But these two different results do not contradict each other. Green and Meyers [100] define a slightly L-dependent lensing template parameter, A L ðΣ mν Þ, and an L-independent template parameter, B lens , which are FIG. 19. Posteriors of the lensing power phenomenological scaling parameters A 2pt and A recon . Different combinations of CMB lensing data added to the Planck primary CMB and BAO are shown in different colors. The 2D contours and solid lines show the marginal posterior from the ΛCDM þ A 2pt þ A recon model; the dot-dashed and dashed lines show the results from ΛCDM þ A 2pt and ΛCDM þ A recon , respectively. The σ values give the level at which A recon (A 2pt ) values less than 1 are excluded, for the case of marginalization over A 2pt (A recon ); i.e., for the solid curves. 7 Recall that in this paper the ΛCDM model has Σm ν fixed to 0.058 eV. 8 When removing DESI, we also add back in BAO data we had removed (see Table III), namely the SDSS DR16 LRG, ELG, and QSO BAO [76] and the DR12 BAO bin with z eff ¼ 0.61 [75].

approximately related to our parameters by A recon ¼ A L ðΣ mν Þ and A 2pt ¼ A L ðΣ mν ÞB lens . 9 The resulting correlation between B lens and Σ mν is much smaller than the correlation we see in Fig. 19, which explains why they see the posterior for Σ mν not change much if they switch from fixing B lens ¼ 1 to marginalizing over it.

We have a strong theoretical prior that the lensing power one would infer from lensed CMB spectra (given the right cosmological model) should be equal, within the expected uncertainties, to what one would infer in a CMB lensing reconstruction. 10 So, even though the bulk of the probability lies to one side of the A recon ¼ A 2pt line, it is worth exploring what we find if we force these two parameters to be equal. This is the model space we explore with an even greater variety of dataset combinations in Fig. 20, and it is also very similar to the model space chiefly studied by Green and Meyers [100] (with their B lens ¼ 1).

One of the dataset variations we consider throughout Fig. 20 is the replacement of the Planck PR3 Plik T&E likelihood with the CamSpec PR4 T&E likelihood [102]. 11The Planck PR4 maps use more timestream data than does PR3 and the CamSpec PR4 T&E likelihood uses more sky than does the Plik PR3 T&E likelihood. If the A 2pt anomaly is driven by some fluctuations that drive the spectra away from the mean (in a manner similar to enhanced lensing) then we would expect the support for it to go down as we add in more data. For Plick PR3 we have, from TT; TE; EE þ lowE, A L ¼ 1.180 AE 0.065 [101]. With the additional PR4 data, we instead have A L ¼ 1.095 AE 0.056 [102]. Note that there is one other PR4 T&E likelihood [103], in which they find an even lower A L ¼ 1.039 AE 0.052.

The lower half of the entries in Fig. 20 do not include any BAO data. We see that without the BAO data any conclusion of A lens > 1 by more than 2σ is not very robust, occurring only for one of the dataset combinations. With the addition of BAO data (top half of the figure), such a conclusion becomes much more robust. The combinations that include Planck T&E, BAO, and all three lensing reconstructions become robust to the difference between PR3 and PR4. With PR3 specifically we find

and a 2.7σ exclusion of A lens < 1. Note that the constraints are weakened by removing DESI from the BAO dataset. Also, if one does so, the results become dependent again on the choice of Planck T&E. Removing DESI (in the same manner as done above) and switching from Planck PR3 T&E to Planck PR4 T&E, we find A lens ¼ 1.064 AE 0.032.

We have checked that switching to the L-dependent lensing template used in [40,100] in place of A lens delivers a similar result. Their template depends on their neutrinomass-like parameter

and a 2.7σ exclusion of Σ mν > 0.058 eV. We note that this evidence for excess lensing power arises from three types of constraints: 1) lensing power measurements, 2) inferences of ΛCDM parameters (other than τ) from the combination of CMB and BAO data, and 3) inference of τ from large angular-scale CMB polarization. We have seen that the preference for A lens > 1 is robust to selection of different datasets for (1). Of particular relevance for this paper, the preference is present when the Planck and ACT lensing measurements are replaced with the SPT one. The SPT lensing is inferred entirely from polarization, making it much less susceptible to possible biases from extragalactic foregrounds. Regarding (2), we see that a more than 2σ preference for A lens > 1 requires either BAO data (including DESI) or PR3 T&E. The next release of BAO data from DESI should be informative, as will the forthcoming lensed TT, TE, and EE spectra from SPT-3G. Note that τ (item 3 above) is inferred from Planck low-lEE data. The importance of τ comes from its implications for A s as the combination A s e -2τ is very well determined by CMB data: a higher τ would lead to higher A s FIG. 20. Constraints on the lensing power phenomenological scaling parameter A lens (68% C.L.) using Planck PR3 Plik T&E likelihood (in blue) and CamSpec PR4 T&E likelihood (in orange). and hence, a higher lensing power prediction. Independent and competitive constraints on τ may come via the kinematic SZ effect [11,104].

The conclusion that there is an excess of lensing power also relies on assumption of the ΛCDM model. Craig et al. [40] introduced a number of changes to the model that could possibly restore concordance. Lynch et al. [105] noted that nonstandard recombination can simultaneously reduce the H 0 tension and lead to a higher predicted lensing power. The DESI Collaboration points out that with free w 0 and w a the constraints on Σm ν are significantly relaxed [65].

J. Other ΛCDM extensions

We also explore ΛCDM extensions with curvature, the density of light relics, primordial abundance of helium, and a time-varying dark energy equation of state. Constraints on additional parameters are summarized in Table V.

For the ΛCDM extension with free curvature and for the ones with changing dark energy equation of state, we find that all CMB datasets are consistent within 2σ with the fiducial values of the ΛCDM model. The constraints from different CMB datasets are consistent. Adding the BAO data significantly tightens the constraints. This is due to the BAO sensitivity to the expansion rate and the distanceredshift relation over redshifts when the dark energy begins to dominate the density.

The energy density of additional relativistic particles other than photons is parameterized by the effective number of neutrino species, N eff , defined via

where ρ γ is the photon energy density. Assuming the standard cosmological model, N eff ¼ 3.044 [106][107][108].

The CMB power spectrum is sensitive to N eff [109][110][111][112][113].

By increasing the density of light relics, one increases the expansion rate during the radiation-dominated era. To keep the precisely measured angular scale of the sound horizon unchanged, the expansion rate after last scattering also increases [e.g., [114] ]. Therefore, there is a positive correlation between H 0 and N eff . The marginal posterior distribution of N eff and H 0 is shown in the left panel of Fig. 21. Combining Planck, SPT, ACT, and BAO, we get

where N eff is consistent with 3.044 from the standard model prediction, but H 0 is still in 5.4σ tension with the SH0ES [83] measurement. One can further extend the model to allow for a free Σm ν as well. The marginal posterior of N eff and Σm ν is shown in the middle panel of Fig. 21. There is nearly no degeneracy between N eff and Σm ν . Adding SPT data to Planck improves the constraints on N eff slightly and Σm ν by about a factor of 2 as discussed in Sec. V H. The lensing information in SPT data is useful for constraining the neutrino mass sum, as discussed in Sec. V H, and the improvement of the N eff constraint is mostly from the more precisely measured EE bandpowers.

The density of light relics is also partially degenerate with the primordial fraction of baryonic mass in helium, Y P . Both N eff and Y P affect the photon diffusion damping of the CMB power spectra at small angular scales. In the right FIG. 21. Marginal posterior distributions of the ΛCDM extension models involving N eff . The horizontal red dashed line shows N eff ¼ 3.044. The yellow vertical band in the left panel shows the constraint on H 0 from Breuval et al. [83]. The vertical brown dashed line shows Σm ν ¼ 0.058 eV. The purple vertical band shows the inference of the primordial helium fraction from Aver et al. [115]. panel of Fig. 21, we show the marginal posterior of N eff and Y P . Compared to the constraints using only Planck data, the uncertainty is reduced by ∼41% with the addition of SPT data. 12 When we fit jointly for N eff and Y P using Planck, ACT, SPT, and BAO measurements, we find values of N eff that are consistent with the standard model prediction of 3.044, N eff -Y P contours that are consistent with BBN predictions, and Y P constraints that are consistent with direct measurements of HII regions in metalpoor galaxies [115].

VI. CONCLUSIONS

In this paper, we presented an optimal joint inference of the CMB lensing potential power spectrum and the unlensed CMB EE power spectrum, derived entirely from CMB polarization maps. These maps were made from observations of 1500 deg 2 taken during the SPT-3G 2019 and 2020 Austral winter observing seasons. The inference was performed using a Bayesian framework called the marginal unbiased score expansion (MUSE) method, which also enables propagation of systematic uncertainties into the estimated bandpower uncertainties. The inferred bandpowers have smaller uncertainties than previously published results at l > 2000 for the EE spectrum and at L > 350 for the lensing spectrum, with a 38σ detection of lensing power.

After checking that the ΛCDM model provides a good fit to the SPT data (combined with a prior on τ), we reported constraints on the ΛCDM parameters. We found H 0 ¼ 66.81 AE 0.81 km=s=Mpc, a significantly tighter constraint than from prior polarization-only measurements. It is consistent with the Planck TT=TE=EE plus Planck lensing constraints assuming ΛCDM, and in 5.4σ tension with the most precise distance ladder measurement [83]. The SPT constraints on Ω c h 2 and Ω b h 2 are also tighter than from prior polarization-only or polarization plus CMB lensing measurements. We also found S 8 ¼ 0.850 AE 0.017, consistent with the result from Planck data and with comparable uncertainty.

We estimated ΛCDM parameters for several other CMB datasets and checked for consistency between these estimates. These consistency tests are stringent tests of the ΛCDM model. All differences investigated were within expectations given the measurement uncertainties and the ΛCDM model. The comparison between SPT and Planck provides a particularly interesting test of the ΛCDM model, because the SPT constraints are more heavily weighted toward small-scale polarization power, and especially lensing power, than is the case for Planck.

We estimated ΛCDM parameters from SPT data in combination with other CMB datasets. When combined with primary CMB and lensing from Planck and ACT, we find H 0 ¼ 67.33 AE 0.37 km=s=Mpc, a 6.2σ tension with the distance ladder [83], and S 8 ¼ 0.8380 AE 0.0084, in 3.3σ tension with S 8 from the 3 × 2 pt analysis of KiDs-1000 [93] and also from the 3 × 2 pt analysis of DES-Y3 [94].

Motivated by the work of Amon and Efstathiou [41] and Preston et al. [95] regarding the S 8 tension, we investigated the sensitivity of our lensing measurement to nonlinear corrections. By introducing the A CMB mod parameter that scales the amplitude of nonlinear corrections to the CMB lensing spectrum, we found A CMB mod ¼ 1.60 AE 0.39. This is the first >3σ "detection" of the nonlinear correction for CMB lensing, and is ∼2σ higher than the value needed for resolving the S 8 tension between galaxy lensing and CMB measurements. Our data clearly do not favor any suppression of the nonlinear contribution, though the uncertainties are still large. A suppressed S 8 could arise at lower redshifts and from slightly smaller scales than the ones we are probing.

The CMB lensing potential is also sensitive to the neutrino mass sum, Σm ν . We investigated constraints on Σm ν with the SPT lensing spectrum. In the ΛCDM þ Σm ν model, we find Σm ν < 0.20 eV ð95% C:L:Þ from the joint fit to Planck, ACT, and SPT data. Adding BAO data, the constraint is further tightened to Σm ν < 0.075 eV ð95% C:L:Þ, placing pressure on the inverted hierarchy minimal mass sum of 0.098 eV.

With the addition of BAO data to the CMB datasets, the strongest challenge to the standard model we see in our results comes from an investigation of lensing power amplitudes. We introduced A 2pt to scale the ΛCDM model CMB lensing spectrum used to calculate lensed model TT=TE=EE spectra and A recon to scale that same ΛCDM model CMB lensing spectrum, but for comparison to the reconstructed CMB lensing spectrum instead. With both of these parameters free, the ΛCDM point in the A recon ; A 2pt plane is outside the 95% confidence region for a variety of dataset combinations. After marginalization over A 2pt , Planck lensing, ACT DR6 lensing, and SPT lensing all favor A recon > 1 at ∼2σ.

We also set A recon ¼ A 2pt and called the single template parameter A lens . We investigated the posterior distribution of A lens for a variety of datasets. For our most comprehensive combination of datasets, we find A lens ¼ 1.083 AE 0.032 and a 2.7σ exclusion of A lens < 1. Following Craig et al. [40] and Green and Meyers [100], we found a similar constraint on their neutrino-mass-like L-dependent lensing template parameter, Σ mν . For the same dataset, we found a 2.7σ exclusion of Σ mν > 0.058 eV. We also saw that these constraints are weakened by dropping DESI BAO, or by switching from Planck PR3 T&E data to Planck PR4 T&E data. Doing both, we find A lens ¼ 1.064 AE 0.032.

Although prior work has shown that delensing can improve constraints on N eff and Y P [81,116], we suspect that the improvement we see here is largely unrelated to the unlensed nature of our bandpowers; i.e., it would be quite similar even if our bandpowers were lensed bandpowers.

We also discussed the cosmological inference with SPT CMB lensing potential bandpowers and the unlensed EE bandpowers on other one-and two-parameter extensions of the ΛCDM model. These model extensions include varying the amount of spatial curvature, the dark energy equation of state and its evolution, the density of light relics, and the primordial fraction of baryonic mass in helium. The constraints on the extended parameters are summarized in Table V. Combining the SPT results with Planck, ACT, and BAO, we saw no significant preference for any of these extensions.

Although the 2019-2020 SPT-3G CMB polarization data analyzed here primarily constrain H 0 and S 8 , SPT-3G constraints on other ΛCDM parameters and extensions will be improved significantly in the near future. The addition of temperature information to both the power spectrum and lensing analyses will significantly improve the constraints on the inferred cosmological parameters, in particular the effective number of neutrino species, Σm ν , and curvature [3]. With polarized beam uncertainties identified as an important contribution to our systematic error budget, we have planned additional deep observations of polarized point sources to better map out the polarized beam response and reduce this source of uncertainty in future analyses. The 1500 deg 2 field has since been observed for three more seasons. The noise level of the full 2019-2023 1500 deg 2 data is expected to be 1.7 times lower than the 2019-20 observations, and the figure of merit for ΛCDM parameters is expected to be 2 times better [3]. The addition of ∼8500 deg 2 more sky mapped with SPT-3G by the end of 2024 will strengthen constraints even further, eventually surpassing those from Planck [3]. The analysis of these data, particularly the ultradeep 1500 deg 2 data, will benefit even more from optimal analyses such as the MUSE method demonstrated here. This work has demonstrated the strength of joint Bayesian CMB lensing potential bandpower reconstruction and unlensed CMB EE bandpower estimation, in particular its ability to propagate systematic uncertainties, and paves the way for future analyses of SPT-3G observations and similarly deep datasets in the future.

these operations for the autodifferentiation (AD) system. The Jacobian of a spherical harmonic transformation (SHT), f ¼ Y a and the Jacobian with respect to f of the lensing operation, f ¼ LðϕÞf, are trivially just Y and L themselves since they are linear operators. The Jacobian of lensing with respect to ϕ is derived from

where x is a pixel location and ϵ is a small perturbation to the lensing potential. The second term demonstrates that the Jacobian is the operator ∇fðx þ ∇ϕÞ • ∇, and involves the spatial gradient of a lensed map, which itself is computed with Lenspyx. With these AD rules implemented, our curved-sky lensing model is fully compatible with the MUSE covariance calculation.

Posterior model.-In the posterior model, we assume the flat-sky approximation, so simulated maps are generated from white-noise maps FFT-convolved with appropriate kernels. We implement lensing with LenseFlow [120], which is GPU-compatible and automatically differentiable up to second-order as needed by the MUSE covariance calculation.

Bandpower parameters.-Both simulation and posterior models take as input a theory unlensed CMB spectrum and gravitational lensing potential spectrum. We introduce free parameters which we will infer, A EE b and A ϕϕ b , which control the amplitude of these spectra relative to a fiducial model within each bin. Our EE bandpowers have linearly spaced bins with Δl ¼ 50 between 350 and 3500, and our ϕϕ bandpowers have 12 logarithmically spaced bins between 20 and 3000. We take the fiducial model to be the Planck best-fit spectrum, but this introduces no dependence on Planck, rather is just a choice of definition of the bandpower amplitude parameters. Within a bin, the parameters scale the theory spectrum uniformly across l's. One can show that with this choice, our inferred A EE b and A ϕϕ b parameters correspond to a constraint on an inversevariance weighted average of the theory spectrum across the bin. Since this inverse-variance weighting arises only implicitly, we also construct explicit bandpower window functions (used in parameter estimation) from an interpolation of the diagonal of the MUSE covariance, Σ MUSE , to every l. The tests on simulations in Sec. IVA 2, which demonstrate no bias and correct scatter of cosmological parameters derived from simulations, provide verification of this procedure.

NFFT-based projections

Within our model, we perform interpolation between maps in different projections in four different places: 1) as mentioned, in the simulation model, the lensed CMB is projected from HEALpix N side ¼ 2048 to a Lambert pixelization, 2) in both simulation and posterior models, the transfer function model involves projecting from a Lambert to an Equirectangular pixelization, applying some weights, then projecting back to Lambert, 3) when computing the transfer function, a set of end-to-end pipeline simulations are compared to a direct projection of the HEALpix N side ¼ 8192 maps which are input to the endto-end simulations, and 4) to model deflections due to relativistic aberration.

To ensure easy differentiability, GPU compatibility, and the best possible accuracy, we perform these interpolations using nonuniform FFTs (NFFTs). NFFT interpolation implicitly uses information from all pixels to compute the interpolated value at a given location, as opposed to other typically-used interpolations such as bilinear interpolation, which uses only the nearest pixel neighbors. It is exact for a band-limited signal (i.e., when there is no power beyond the Nyquist frequency) and remains very accurate for nearly band-limited signals such as ours, where power beyond the Nyquist frequency is highly suppressed by the beam.

We use the CPU and GPU compatible NFFT library NFFT.jl [121]. An interpolation is performed by a pair of forward and adjoint NFFTs, F , with appropriate grid points, γ, which describe the relative location of the pixel centers in the map projections being interpolated between

The adjoint of the projection operator, which arises in automatic differentiation, is then simply

We verify these projections are accurate enough in Fig. 22. For use cases (1) and ( 2), in the left panel, we generate a simulated lensed CMB map in the Lambert pixelization. We then project to equirectangular and back to Lambert, as well as to HEALpix N side ¼ 2048 and back to Lambert. We then examine the change to the spectrum of the original map. Figure 22 demonstrates this change is ≲0.3% up to the maximum multipole used in the analysis, l < 3500. For use case (3), in the right panel, we compare the spectrum of the projected HEALpix N side ¼ 8192 maps to the input theory spectrum, again finding ≲0.5% changes up to the maximum multipole. These differences are well below our uncertainties; for example, typical EE bandpower errors are at the 5% to 10% level.

a. Temperature beams

The first step in determining our polarized beams is to determine the temperature beams. At low radii (≲2 arcmin), temperature beams are determined from observations of several bright AGN in our observing field. The AGN, however, do not have sufficient brightness to fully map out the beam at larger scales (referred to as the beam "sidelobes"), where instead we use observations of Saturn. In turn, observations of Saturn cannot be used at small scales, since its intensity saturates detectors and changes the response as compared to typical CMB observations. Therefore, we stitch the AGN and Saturn maps into projections of the HEALpix N side ¼ 8192 maps, which are input to mock observation (labeled "projskies") vs the input theory which generated those HEALpix maps. Blue bands represent 1 and 2σ Monte Carlo error from the finite number of simulations used to compute the blue curve. This confirms the "projskies" sufficiently recover the input theory so as not to bias the transfer function calculation, where they are used. For reference, this plot also shows the expected change to the power spectra if bilinear interpolation was instead used, which shows much worse performance than NFFT interpolation. a composite beam map, which uses the AGN measurements at low radii, and Saturn for the sidelobes.

To create the composite map, we must account for several systematic effects. First, Saturn is very slightly extended. Before stitching, we convolve the AGN maps with a disk of 8.55 arcsec (the average angular diameter of Saturn during our observations). We later deconvolve the same disk from the final beam. The SPT-3G detectors have a finite response time, which spreads signals along the scan direction. While the response time is constant, the on-sky scan speed varies with declination. Therefore, we deconvolve the map-space effect of the time constants from each map (both for individual AGN and Saturn), and then reconvolve with the mean effect of the time constant over the SPT-3G winter field footprint. We also attempt to remove background signals (primarily the CMB) from both the Saturn maps and AGN maps. Since Saturn moves across the sky, we are able to subtract direct observations of the backgrounds from the Saturn maps. To do this, we mock observe Planck maps for each of our Saturn observations, and subtract the mockobserved maps from our Saturn maps. For the AGN maps, we are unable to do this. However, we are only concerned with small scale (≲7 arcmin) information in the AGN maps. On these scales, the CMB varies very little, and can be modeled with only low-order modes. Therefore, we fit a simple background model to each AGN consisting of a slope in radius and an offset. We fit these parameters (using the Saturn maps as a reference for the beam) at 95 and 150 GHz, since the CMB is expected to be common between bands, and the 220 GHz maps are significantly noisier (both from instrumental noise, and higher astrophysical backgrounds). We subtract the background model from all three bands.

Finally, we apply a very basic model of intra-band frequency dependence to the beam. We model the beam at each band as a Gaussian made of two components, one which is fixed across the band, and one which scales with frequency. We model the frequency scaling as σðνÞ ∝ 1=ν, which is the expected behavior for a diffraction-limited beam. We use the same constant of proportionality across all three bands. In order to fit these parameters, we assume that the 95 GHz beam is diffraction limited (i.e., the frequency independent component is nonexistent). We approximate the "effective beam frequency" for each band using the SPT-3G mean bandpass functions, and the source spectrum. Then, we fit the three remaining parameters (the constant of proportionality in the frequency dependent portion, and the two frequency independent portions of the 150 and 220 GHz beams) using observations of the brightest AGN. This simple model allows us to approximate the change in the beam by calculating the effective beam frequency for different source spectra, and convolving the beam maps (or multiplying their spectra) with the appropriate Gaussian. We apply this model to transform the final beam from the mean AGN spectrum to the CMB spectrum, and to transform the Saturn maps to match the AGN spectrum before stitching.

With the systematic effects accounted for, the stitching process is relatively simple. For each pair of AGN and Saturn maps, we fit for an amplitude to account for the difference in intensity between the sources. As discussed above, we cannot use the innermost portion of the Saturn maps. The AGN maps also become noise-dominated at large distances from the source, so we perform this fit over an annulus with inner radii of (2.25, 1.5, 1.0) arcmin, and outer radii of (3.5, 3.0, 2.5) arcmin at (95,150,220) GHz. Finally, each composite map is constructed by linearly transitioning from the AGN map to the Saturn map over the same annulus.

We make composite maps from all pairs of AGN and Saturn maps. We also vary several systematic parameters: the choice of stitching radii, uncertainty in the time constants, the choice of the background fitting region. For each of these variations, we create the full set of composite maps. Then, we construct cross spectra of composite maps such that none of the component maps are the same. If we included the same component map more than once, we would incur a noise bias in the beam. To measure the final beam, we take the mean of these cross spectra and apply the two final effects described above (specifically, deconvolving the Saturn disk, and converting from the AGN spectrum to the CMB spectrum).

The beam uncertainty is measured from the cross spectra. By including both different Saturn and AGN maps, as well as the systematic variations, we capture both the map noise, systematic uncertainty from our analysis choices, and systematic differences between the AGN. Using the Monte Carlo method, this uncertainty is propagated to a covariance between temperature B ν l 's at each frequency, ν, and multipole, l. Finally, we add one additional uncertainty directly to the covariance matrix, which accounts for uncertainty in the detector bandpasses and the AGN spectral index. This covariance is decomposed into eigenmodes, and in the analysis the amplitude of five eigenmodes, β n is free to vary jointly with the bandpower and other systematics parameters. We note that these are eigenmodes of the joint covariance between frequencies and multipoles, so each β n changes the beam simultaneously at each frequency (in different ways at each frequency). The total uncertainty in the determination of the temperature beams is 0.2% at l ¼ 3000, and is subdominant to the uncertainty on polarized beams which we now describe.

b. Polarization beams

While the temperature beams are mapped out to very high signal-to-noise at angular scales relevant to this analysis via observations of bright point sources, we lack sufficiently bright polarized sources to similarly map the sidelobe response of the polarized beams. Prior to unblinding, we had assumed the polarized and temperature beams were identical at the scales used in the analysis. However, the temperature beam has significant diffuse sidelobes due to diffraction and scattering that may not be fully polarized. To account for this, we devised a model of the beams which allows the polarization of the beam sidelobes to vary. This model assumes the beams are made up of a fully polarized main central beam, and a diffuse beam sidelobe which may be partially unpolarized.

We consider four components in the model for the main beam 1) Gaussian illumination pattern of the primary by the detector lenslets, 2) truncation of that illumination pattern by the Lyot stop, 3) averaging the frequency dependent beam over the detector bandpass weighted by a CMB spectrum, and finally 4) Gaussian broadening of the diffracted beam due to a frequency independent geometric aberration.

We model the electric field illumination pattern on the primary mirror as a truncated Gaussian,

where σ 0 is the width of the Gaussian beam from the lenslets at some arbitrary fiducial frequency ν 0 , and ΠðRÞ is a disk of radius R,

For SPT-3G, R ≈ 3.75 m. The illumination on the sky can then be computed from the Fraunhofer diffraction equation,

Uðθ; νÞ ¼ 2π

where z is the distance to the image plane, ρ is the distance from the origin on the image plane, and J 0 is the zeroth order Bessel function. In the second line, we have taken the small angle approximation ρ=z ¼ θ to put this in terms of angle on the sky. The beam on the sky is the intensity which is the square of the electric field pattern, Bðθ; νÞ ideal ¼ Uðθ; νÞ 2 : ðA10Þ

The beam for our wide frequency bands is then integrated over the bandpass and source spectrum. We integrate over the beam intensity, since the electric field phases will decohere at any significant frequency separation. Adjacent frequencies where phase coherence matters have very nearly the same illumination pattern, because both the spectra of our sources and bandpasses vary slowly. Thus, we have

where tðνÞ is the bandpass function, and IðνÞ ¼ dB dT j T¼T CMB is the observed spectrum.

Finally, we need to account for the effect of geometric aberrations. This is simply a matter of convolving BðθÞ ideal with a Gaussian, a representation of geometric aberrations, so the final beam is

In principle, we expect σ geom to be the same for all three bands. This model for the main beam has five free parameters, σ 0 , σ geom , and an overall amplitude at each frequency, A ν , and predicts the main beam simultaneously in each of our bands. We fit this model to the observed temperature beams as determined in the previous section, limiting to an inner radius of θ < 0.75 arcmin where we observe a good fit, suggesting the main beam dominates and the fit is not biased by the presence of temperature beam sidelobes. For this fit, we construct a χ 2 of the form

2 , which ignores uncertainties in the temperature beams, since we expect the systematic uncertainty from our choice of physical model to dominate to the total uncertainties in fit.

With this fiducial main beam determined, we denote the difference between this and the temperature beam as the sidelobes, and model the polarized beams as the sum of the fixed main beam and the sidelobes, the latter scaled by some unknown scale-independent polarization fractions, β ν pol . That is, the polarized beams used in our analysis are

where the temperature beam at each frequency is

where B 0 is the best-fit beam and B n is the nth eigenmode of the beam covariance. The β n is the amplitude of the nth eigenmode, which is the free parameter. The B n is normalized such that β i has a unit normal prior. In Fig. 23, we show the beam used in this analysis.

Transfer function

The transfer function accounts for the filtering of the timestreams during mapmaking described in Sec. II. In our transfer function model TF , we must approximate the impact of this filtering directly to a Lambert-projected map, since our model is not at the level of timestreams. The nontrivial procedure of sequentially deprojecting templates while masking point sources still corresponds to a linear operator, but not one which is diagonal in any easily accessible basis.

Our model begins by first reprojecting the map to an equirectangular projection using NFFTs, such that individual scans are straight lines in the horizontal direction. We then create a single map-level deprojection template by multiplying the point source mask assumed in mapmaking and applying a 1D filtering in the horizontal direction using FFTs where the 1D filter is a free function which we will fit for using simulations. This template is subtracted from the entire map (including the masked regions), before projecting back to Lambert. The operation can be represented as

where P is the NFFT projection from Lambert to equirectangular, M ptsrc is the point source mask and TFðl x Þ is the fitted FFT kernel in the horizontal direction.

To fit for the appropriate function, TFðl x Þ, we compare 1) mock observations of the signal, m, which simulate the entire timestream filtering and mapmaking process, and 2) direct projections of the simulated signal maps which were fed into the mock observations, p. The relation between these maps are

The factor of the HEALpix N side ¼ 8192 pixel window function is due to the interpolation of p onto mock timestreams during the simulation process, and does not arise in the real data. The factor of Lambert 2.25 0 pixel window function, PWF 2.25 0 arises from the binning of these timestreams into 2.25 0 resolution maps. We then minimize the function km truem model k 2 by fitting TFðl x Þ at a set of spline points. We find that a single pair of m and p simulations is sufficient to obtain a stable estimate and FIG. 23. Beams used in this analysis, in real space (top row) and multipole space (bottom row), at each frequency (each column). Black lines show the temperature beams. Blue and orange lines show our best fits to the main beam, where the fitting region is denoted by the gray shaded region. We used the physical main beam for our final result, but we also verified that using the Gaussian main beam leads to no significant shifts in results (see Sec. IV C). Green lines and green bands represent the posterior mean and 1σ uncertainty on the polarized beams from the MUSE analysis, which includes marginalization over β n , β pol , and all other systematics and bandpower parameters.

have verified the result changes negligibly between a few different simulations.

The left panel of Fig. 24 shows the fitted TFðl x Þ. As expected, it is roughly zero until around l x ¼ 300, then transitions to unity at higher multipoles. We choose to allow some small negative values as it gives a better overall fit with no other observed impact. The right two panels of the figure then show the comparison of EE and BB spectra of the mock observations with those computed by applying the best-fit model transfer function to p. Across the multipole range used for the analysis, the model is accurate to ≲0.2% in EE and ≲2% in BB, both less than 10% of our uncertainties (which are denoted as gray bands).

APPENDIX B: CALIBRATION PRIORS AND ALTERNATIVE SYSTEMATICS ESTIMATES

Here, we describe the priors on A 150;i cal which enter the MUSE analysis, and the alternative estimates of A 95;i cal , A 220;i cal , ϵ Q , ϵ U , and ψ pol which are compared against the MUSE results.

In practice, calibrating to Planck is performed as two steps. First, we determine a calibration factor for temperature maps in each subfield, T i cal . Then, we determine a correction on top of this for polarization maps, P cal , which is assumed the same across all subfields. By our convention, these factors multiply the SPT data, rather than multiplying the theory model like the MUSE definition of A ν;i cal . The relation is thus

The temperature calibration is determined by computing the cross-spectrum of the full-depth 150 GHz SPT-3G data with the full-depth 143 GHz Planck map passed through the SPT-3G observing pipeline normalized by the crossspectrum of two half-depth SPT-3G maps,

150 cal;external ¼ SPT 150 full × Planck 143 full SPT 150 half1 × SPT 150 half2 :

The uncertainty on T cal;external is computed from the scatter in the quantity when it is computed using 100 different subsets of SPT observations entering the SPT maps. After unblinding, we also realized we had neglected to fold in a 0.25% uncertainty [125] on the Planck absolute calibration [126], which leads so insignificant changes, but is fixed in our final results.

Although not used in MUSE, we also obtain temperature calibrations of 95 and 220 GHz SPT-3G maps by crosscorrelating against the 150 GHz map (this yields tighter constraints than calibrating these against Planck). The calibrations are

95 cal;internal ¼ SPT 150 half1 × SPT 150 half2 SPT 95 half1 × SPT 150 half2 T 150 cal;external ; T 220 cal;internal ¼ SPT 150 half1 × SPT 150 half2 SPT 220 half1 × SPT 150 half2 T 150 cal;external :

We recalibrate the Q and U maps by this factor, then correct them for temperature-to-polarization leakage and polarization angle calibration as described below. Then, we determine P 150 cal by cross-correlating the SPT 150 GHz polarization maps to Planck 143 GHz polarization maps in a manner analogous to the absolute temperature calibration in Eq. (B3). The estimated P 150 cal is then corrected by an additional factor of ffiffiffiffiffiffiffiffiffiffiffi 0.966 p to match the final recalibration factor applied to the 143 GHz Planck polarization bandpowers [Eq. ( 45), [125] ], and additionally has a 0.5% uncertainty folded in corresponding to the uncertainty on the Planck polarization calibration. Similarly used only for comparison, we then determine P 95 cal and P 220 cal by internally cross-correlating the SPT maps.

Next, we describe the alternative estimate of temperature-to-polarization leakage and global polarization angle correction which is used above in determining the calibration factors and is compared against the MUSE results in Fig. 7. The main difference from the MUSE estimate is that these are performed at the power spectrum rather than at the map level, and are done with an iterative procedure which assumes a theory cosmological model rather than simultaneously and jointly with all systematics and bandpower parameters.

Independent estimates of the monopole leakage coefficients, ϵ Q;TT , ϵ Q;TQ , ϵ U;TT , are formed by fitting a template to the C TQ l and C TU l cross-spectra using cross-spectra of the half-depth data maps and cross-spectra from simulations that have no monopole leakage, C TQ l;template ¼ ϵ Q;TT C TT l;data þ ϵ Q;TQ C TQ l;sim ; ðB4Þ C TU l;template ¼ ϵ U;TT C TT l;data þ ϵ U;TU C TU l;sim :

In practice, C TU l;sim scatters around null, and therefore ϵ U;TU can take arbitrary values. For this reason, the second term is