# CLASS Angular Power Spectra and Map-Component Analysis for 40 GHz
Observations through 2022

Joseph R. Eimer Joseph R. Eimer<EMAIL_ADDRESS>William H. Miller III
Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD
21218, USA  Yunyang Li (李云炀) William H. Miller III Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  Michael K.
Brewer  William H. Miller III Department of Physics and Astronomy, Johns
Hopkins University, Baltimore, MD 21218, USA  Rui Shi (时瑞) William H. Miller
III Department of Physics and Astronomy, Johns Hopkins University, Baltimore,
MD 21218, USA  Aamir Ali William H. Miller III Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  John W. Appel
William H. Miller III Department of Physics and Astronomy, Johns Hopkins
University, Baltimore, MD 21218, USA  Charles L. Bennett William H. Miller III
Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD
21218, USA  Ricardo Bustos Departamento de Ingeniería Eléctrica, Universidad
Católica de la Santísima Concepción, Alonso de Ribera 2850, Concepción, Chile
David T. Chuss Department of Physics, Villanova University, 800 Lancaster
Avenue, Villanova, PA 19085, USA  Joseph Cleary William H. Miller III
Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD
21218, USA  Sumit Dahal NASA Goddard Space Flight Center, 8800 Greenbelt Road,
Greenbelt, MD 20771, USA  William H. Miller III Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  Rahul Datta
Department of Astronomy and Astrophysics, University of Chicago, 5640 South
Ellis Avenue, Chicago, IL 60637, USA  William H. Miller III Department of
Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
Jullianna Denes Couto William H. Miller III Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  Kevin L. Denis
NASA Goddard Space Flight Center, 8800 Greenbelt Road, Greenbelt, MD 20771,
USA  Rolando Dünner Instituto de Astrofísica, Facultad de Física, Pontificia
Universidad Católica de Chile, Avenida Vicuña Mackenna 4860, 7820436, Chile
Centro de Astro-Ingeniería, Facultad de Física, Pontificia Universidad
Católica de Chile, Avenida Vicuña Mackenna 4860, 7820436, Chile  Thomas
Essinger-Hileman NASA Goddard Space Flight Center, 8800 Greenbelt Road,
Greenbelt, MD 20771, USA  Pedro Fluxá  Centro de Astro-Ingeniería, Facultad de
Física, Pontificia Universidad Católica de Chile, Avenida Vicuña Mackenna
4860, 7820436, Chile  Johannes Hubmayer National Institute of Standards and
Technology, Boulder, CO 80305 USA  Kathleen Harrington High Energy Physics
Division, Argonne National Laboratory, 9700 S. Cass Avenue, Lemont, IL 60439,
USA  Department of Astronomy and Astrophysics, University of Chicago, 5640
South Ellis Avenue, Chicago, IL 60637, USA  Jeffrey Iuliano Department of
Physics and Astronomy, University of Pennsylvania, 209 South 33rd Street,
Philadelphia, PA 19104, USA  William H. Miller III Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  John Karakla
William H. Miller III Department of Physics and Astronomy, Johns Hopkins
University, Baltimore, MD 21218, USA  Tobias A. Marriage William H. Miller III
Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD
21218, USA  Carolina Núñez William H. Miller III Department of Physics and
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  Lucas Parker
Space Remote Sensing and Data Science, Los Alamos National Lab, Los Alamos, NM
87545  Matthew A. Petroff Center for Astrophysics, Harvard & Smithsonian, 60
Garden Street, Cambridge, MA 02138, USA  Rodrigo A. Reeves CePIA, Astronomy
Department, Universidad de Concepción, Casilla 160-C, Concepción, Chile
Karwan Rostem NASA Goddard Space Flight Center, 8800 Greenbelt Road,
Greenbelt, MD 20771, USA  Deniz A. N. Valle William H. Miller III Department
of Physics and Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA
Duncan J. Watts Institute of Theoretical Astrophysics, University of Oslo,
P.O. Box 1029 Blindern, N-0315 Oslo, Norway  Janet L. Weiland William H.
Miller III Department of Physics and Astronomy, Johns Hopkins University,
Baltimore, MD 21218, USA  Edward J. Wollack NASA Goddard Space Flight Center,
8800 Greenbelt Road, Greenbelt, MD 20771, USA  Zhilei Xu (徐智磊) MIT Kavli
Institute, Massachusetts Institute of Technology, 77 Massachusetts Avenue,
Cambridge, MA 02139, USA  Lingzhen Zeng Center for Astrophysics, Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA

###### Abstract

Measurement of the largest angular scale ($\ell<30$) features of the cosmic
microwave background (CMB) polarization is a powerful way to constrain the
optical depth to reionization, $\tau$, and search for the signature of
inflation through the detection of primordial $B$-modes. We present an
analysis of maps covering nearly 75% of the sky made from the ground-based
$40\,\mathrm{GHz}$ channel of the Cosmology Large Angular Scale Surveyor
(CLASS) from August 2016 to May 2022. Using fast front-end polarization
modulation from the Atacama Desert in Chile, we show this channel achieves
higher sensitivity than the analogous frequencies from satellite measurements
in the range $10<\ell<100$. After a final calibration adjustment, noise
simulations show the CLASS linear (circular) polarization maps have a white
noise level of $125\,(130)\,\mathrm{\mu K\,arcmin}$. We measure the Galaxy-
masked $EE$ and $BB$ spectra of diffuse synchrotron radiation and compare to
space-based measurements at similar frequencies. In combination with external
data, we expand measurements of the spatial variations of the synchrotron
spectral energy density (SED) to include new regions of the sky and measure
the faint diffuse SED in the harmonic domain. We place a new upper limit on a
background of circular polarization in the range $5<\ell<125$ with the first
bin showing $D_{\ell}<0.023$ $\mathrm{\mu K^{2}_{CMB}}$ at 95% confidence.
These results establish a new standard for recovery of the largest-scale CMB
polarization from the ground and signal exciting possibilities when the higher
sensitivity and higher frequency CLASS channels are included in the analysis.

Early Universe (435); Cosmic microwave background radiation (322);
Observational Cosmology (1146); Polarimeters (1127); Astronomy Data Analysis
(1858)

††software:  numpy (Harris et al., 2020), scipy (Virtanen et al., 2020),
matplotlib (Hunter, 2007), astropy (Astropy Collaboration et al., 2022),
HEALPix (Górski et al., 2005), camb (Lewis et al., 2000), fastcc (Peel et al.,
2022), pysm (Thorne et al., 2017), PolSpice (Chon et al., 2004), NaMaster
(Alonso et al., 2019)

##  1 Introduction

Measurements of the cosmic microwave background temperature and polarization
angular power spectra (e.g., Hinshaw et al., 2013; Planck Collaboration et
al., 2020a; Aiola et al., 2020; Balkenhol et al., 2021), in combination with
surveys of large-scale structure (e.g., Bocquet et al., 2019; Asgari et al.,
2021; Amon et al., 2022; Dalal et al., 2023; Madhavacheril et al., 2023),
baryon acoustic oscillations (e.g., Alam et al., 2021; Abbott et al., 2022),
and cosmic expansion (e.g., Brout et al., 2022), have led to a well-
constrained cosmological model, which includes a cold dark matter (CDM)
component and parameterizes dark energy as a cosmological constant, $\Lambda$.
The resultant $\Lambda$CDM model describes the matter and energy content of
the universe, as well as its age, ionization history, and a statistical
description of its initial conditions.

A hypothesized early period of accelerating expansion, termed inflation,
provides an explanation for the nearly flat, homogeneous, and isotropic
universe we observe, as well as a mechanism for producing the nearly scale-
invariant spectrum of density perturbations needed to seed large-scale
structure within $\Lambda$CDM cosmology (Starobinskiǐ, 1979; Guth, 1981;
Linde, 1982). If such an inflationary period occurred, a background of
primordial gravitational waves (tensor modes) is expected to have imprinted an
odd-parity, “$B$-mode,” component in the CMB polarization, along with even-
parity, “$E$-mode,” polarization and temperature anisotropy (Seljak &
Zaldarriaga, 1997; Kamionkowski et al., 1997). While anisotropy in the
temperature and $E$-mode polarization are expected and observed from density
(scalar) perturbations, primordial $B$-mode polarization is only expected from
tensor modes. The $BB$ angular power spectrum amplitude can be related to the
energy scale of the processes driving the inflationary expansion through the
tensor-to-scalar ratio, $r$. The tightest constraint ($r<0.036$) comes from
the BICEP/Keck $BB$ measurement at $\ell>30$ supported by Planck and WMAP data
(BICEP/Keck Collaboration et al., 2021). Other recent $B$-mode measurements at
$\ell>30$ include Kusaka et al. (2018); Polarbear Collaboration et al. (2020);
Spider Collaboration et al. (2022). $BB$ measurements on the largest angular
scales ($\ell<30$) probe a distinctive CMB polarization enhancement due to
Thomson scattering since recombination by the ionized intergalactic medium
during and after reionization. The tightest constraint on the tensor-sourced
$BB$ spectrum at $\ell<30$ is an upper limit from Planck ($r<0.274$, de
Belsunce et al., 2023).

Tensor modes from inflation also source low-$\ell$ $E$-mode polarization, but,
given current limits on $r$, the tensor contribution to $E$-mode signal is
dominated by contributions from scalar perturbations. Like the $BB$ spectrum,
the $EE$ spectrum at $\ell<30$ is enhanced by Thomson scattering of CMB
photons by the ionized intergalactic medium, and therefore contains
information regarding the processes by which the first stars converted neutral
hydrogen into its current ionized state (Hu & Holder, 2003; Heinrich et al.,
2017; Watts et al., 2020). This so-called “reionization peak” has an amplitude
proportional to the square of the optical depth to reionization ($\tau^{2}$,
Zaldarriaga, 1997). Of the six standard $\Lambda$CDM parameters, the least
well constrained is $\tau$: while other parameters are measured with sub-
percent precision, the reionization optical depth uncertainty is 10% (Pagano
et al., 2020). Improvement in the measurement precision of $\tau$ will break a
critical degeneracy with the amplitude of the initial scalar curvature
fluctuations $A_{s}$. Improved knowledge of $\tau$, $A_{s}$, and upcoming
measurements of large-scale structure will provide constraints on the sum of
the neutrino masses at a level sufficient to address the fundamental question
of the mass eigenstate hierarchy of neutrinos (Allison et al., 2015; Liu et
al., 2016).

The Cosmology Large Angular Scale Surveyor (CLASS) is a project to measure the
largest angular scale CMB polarization from the ground (Essinger-Hileman et
al., 2014; Harrington et al., 2016). Located at 5140 m elevation in the
Atacama Desert of Northern Chile, the CLASS telescopes observe in frequency
bands near 40, 90, 150, and 220 $\mathrm{GHz}$ (Dahal et al., 2022). CLASS was
designed to address sources of systematic error, such as instrumental
polarization and ground pick-up, that traditionally limit measurements at the
largest angular scales. Design elements include front-end polarization
modulation (Chuss et al., 2012; Harrington et al., 2018, 2021), a low-spill
optical design with entrance pupil at the modulator (Eimer et al., 2012; Xu et
al., 2020), and high-efficiency detectors with well-defined frequency/mode
acceptance (Rostem et al., 2012; Appel et al., 2014, 2019). A unique aspect of
CLASS is that its polarization modulator measures both linear and circular
polarization (Padilla et al., 2020; Petroff et al., 2020b). Other projects
targeting the largest angular scale microwave polarization include ground-
based [GroundBIRD, (Lee et al., 2020), LSPE-strip, (Addamo et al., 2021),
C-BASS (Harper et al., 2022), S-PASS (Carretti et al., 2019), and QUIJOTE
(Rubiño-Martín et al., 2023)], balloon-borne [LSPE-SWIPE (Addamo et al., 2021)
and Taurus (Benton, 2020)], and space-based [LiteBIRD (LiteBIRD Collaboration
et al., 2023)] efforts.

In this paper, we present polarization maps, angular power spectra, and
analyses of Galactic and cosmological components based on 40 $\mathrm{GHz}$
observations from August 2016 to May 2022. In a companion paper, we describe
the 40 $\mathrm{GHz}$ instrument, observations, and data reduction to maps (Li
et al., 2023, hereafter L23), and we include an overview of these topics in
Sections 2 and 3. In Section 4, we describe the method of angular power
spectrum estimation for data and simulations. Internal consistency tests of
the maps, comprising a battery of jackknife “null” tests, are presented in
Section 5. In Section 6, we present comparisons to other data sets, including
a final calibration from comparison to the bright synchrotron signal near the
galactic plane established by WMAP. Galactic synchrotron and cosmological
component analyses are given in Sections 7 and 8. We conclude in Section 9.

##  2 CLASS overview

###  2.1 Instrument Description

The CLASS telescope array resides in the Atacama Desert of northern Chile and
surveys the sky with single-frequency telescopes at 40 $\mathrm{GHz}$ and 90
$\mathrm{GHz}$, as well as a dual-frequency 150/220 $\mathrm{GHz}$ telescope.
With its observing site at $23^{\circ}$ S latitude, large instantaneous field
of view, and ability to perform constant-elevation scans over $720^{\circ}$ in
azimuth, the array can efficiently survey nearly 75% of the sky. This full
area is mapped daily, apart from a Sun avoidance region that changes
throughout the year, with boresight angle changed each day to improve
polarization angle coverage and check for systematic errors, e.g., due to
unwanted ground pickup. The CLASS frequency coverage enables characterization
of both the CMB and polarized Galactic emission from synchrotron and dust with
sensitivity centered near the polarized foreground minimum.

The telescopes share a common architecture. Key to enabling the stability
required to recover large angular scale features on the sky, the first optical
element in each telescope is a variable-delay polarization modulator (VPM).
The VPM consists of a wire grid polarizer in front of a moving mirror which
modulates one linear Stokes parameter into circular polarization and vice
versa, which provides CLASS sensitivity to both linear and circular
polarization (Chuss et al., 2012; Harrington et al., 2018). Each VPM has a
$60\,\mathrm{cm}$ diameter clear aperture. The fast (10 Hz) front-end
polarization modulation elevates the signal band out of the dominant noise at
low frequencies, e.g., $1/f$ noise from bright atmospheric emission partially
polarized after the VPM, while mitigating systematic errors associated with
polarization induced by other optical elements in the system (Harrington et
al., 2021,  Cleary et al. in prep.).

Following the VPM, the signal is guided via two ambient-temperature mirrors
into a cryogenic receiver where cold ($<4$ K) dielectric lenses focus the
signal onto a focal plane array at $\sim 50\,\mathrm{mK}$ (Iuliano et al.,
2018). The focal planes use custom profiled, smooth-walled feedhorns (Zeng et
al., 2010) to guide the signal onto planar microwave circuits containing band-
defining and out-of-band rejection filters before terminating the signal on a
pair of superconducting transition edge sensor (TES) bolometers—one TES for
each linear polarization (Rostem et al., 2012; Appel et al., 2014, 2019; Dahal
et al., 2018, 2022). The detectors are read out through a series of cryogenic
time-division multiplexing amplifiers with ambient-temperature multichannel
electronics (MCE) (Reintsema et al., 2003; Battistelli et al., 2008; Doriese
et al., 2016).

The telescopes are housed within a comoving ground screen with a flared
forebaffle extension opening toward the sky. For most of the survey described
in this paper, a thin polypropylene sheet was held taut at the base of the
forebaffle extension to environmentally shield the telescope from dust and
potentially damaging weather events. Identified as a source of wind-induced
systematic error, the closeout film was removed for periods of good weather
after September 2021.

Pointing is achieved through a three-axis mount—boresight over elevation over
azimuth. During normal operation, the boresight and elevation pointing are
held fixed while the azimuth is scanned.

This paper uses data from the 40 $\mathrm{GHz}$ telescope observations
beginning 31 August 2016, after initial commissioning, and extending to 19 May
2022, when the telescopes were shut down for maintenance and instrument
upgrades. We define the period before 22 February 2018 as _Era 1_ and after 22
June 2018 as _Era 2_ —various instrument changes and upgrades were performed
in the intervening time. More details can be found in L23.

###  2.2 Observing Strategy

Observations were organized and scheduled on a 24-hour cycle using custom
software (Petroff et al., 2020a). Usually, the daily observing schedule began
at midday when the boresight is set to one of seven values, one for each day
of the week, from $-45^{\circ}$ to $+45^{\circ}$ relative to vertical in
$15^{\circ}$ increments. After the cryogenic system stabilized following the
boresight adjustment, detector loading was measured by performing a current-
voltage (_IV_) sweep, and optimal detector bias values were selected (Appel et
al., 2022). The set of data collected after the IV measurement until the next
IV or until the schedule terminated is called a _span_. At nighttime, the
telescope nominally scanned $\pm 360^{\circ}$ in azimuth at $45^{\circ}$
elevation and constant boresight. The azimuth scanning direction typically
reversed at $180^{\circ}$ azimuth (due South) during the night. In the
daytime, the azimuth range was reduced to keep the array boresight center
$>20^{\circ}$ from the Sun. (No other celestial objects were avoided by the
telescope.) For a schedule dedicated to the CMB survey, observations spanned
nearly 24 hours, ultimately terminating in preparation for the next day’s
boresight adjustment. At times, however, the CMB survey was interrupted by
weather, instrument malfunction, and dedicated observations of the Moon or
planets conducted for calibration purposes (Xu et al., 2020; Dahal et al.,
2021). Further details regarding the survey and instrumental configuration
changes can be found in L23.

##  3 Data processing and maps

The data collected by CLASS are unique in the field of CMB surveys. These are
the first and only data collected by telescopes employing a front-end VPM to
stabilize the polarization measurement and guard against certain systematic
errors. In this section we summarize the collection and processing of these
data. A more complete description is given in L23.

Figure 1: The $Q$/$U$/$V$ Stokes parameter maps for the CLASS 40
$\mathrm{GHz}$ channel are shown in Galactic coordinates using a Mollweide
projection. Gray regions indicate portions of the sky not fully sampled by
CLASS. To show both bright synchrotron features within the Galactic plane and
the faint regions elsewhere in the map, the color scale is linear within $\pm
5\,\mu$K and logarithmic beyond this range. The dipolar atmospheric $V$ signal
has been removed through the mapmaking filters, and no remaining Stokes $V$
signal is detected—the map is consistent with noise. The survey coverage is
non-uniform; see L23 for hit maps. In addition to anisotropy in the signal,
the noise also changes over the sky.

###  3.1 Data Acquisition and Pre-processing

The total data volume includes 4.32 years, excluding downtime for the 90
$\mathrm{GHz}$ telescope installation and the COVID-19 pandemic. After data
taken during periods of poor weather or systematic instrument malfunction were
removed, 45.4% of the total volume remained. At this point, contiguous data
were grouped into spans. Most spans comprise 22 hours of data, interrupted at
midday for changing boresight and IV calibration. However, spans may also be
shorter if the observations were interrupted. We additionally cut data taken
for calibration purposes or that were of lower quality, leaving 28% of the
total volume, equivalent to 86.77 detector$\cdot$years of data.

The raw data were affected by the finite electro-thermal response of the TES
bolometers, which is sufficiently modeled as a single-pole filter with
$3$–$4\,\mathrm{ms}$ time constant. The MCE readout electronics also applied
an anti-aliasing filter to the data. Both of these filters were deconvolved
from the selected data. The data were then calibrated from binary MCE-units to
power using the IV-based calibration (Appel et al., 2022).

###  3.2 Demodulation

As described in more detail in L23, the effect of the VPM is to modulate the
amplitude of the polarization signal at a stable frequency near
$10\,\mathrm{Hz}$ and its harmonics. This moves the signal in the raw data to
sidebands of the modulation frequency components. _Demodulation_ is the
process by which the modulated signal is recovered. For this analysis, the
first step was to filter the data with a passband of approximately
$9$–$51\,\mathrm{Hz}$, which preserved the encoded sky signal in the modulated
data and band limited the noise beyond the instrument signal band. The next
step was to use knowledge of the VPM modulation (Chuss et al., 2006), the
source spectra (Petroff et al., 2020b; Thorne et al., 2017), the telescope
bandpass (Dahal et al., 2022), and the geometric orientation of the VPM within
the telescope (Eimer et al., 2012) to solve for the polarization signal
incident on the VPM. In the process of solving for the signal, the data were
further low-pass filtered with a cutoff frequency between 0.5 and 1
$\mathrm{Hz}$—with the value chosen depending on the telescope scan speed. The
data were then downsampled, and gaps corresponding to cut data were filled by
linear interpolation with white noise to limit the introduction of spurious
artifacts when treating the data as a complete, uniformly sampled record.

###  3.3 Filtering and Mapmaking

Before mapping, the demodulated data were filtered to remove systematic errors
that were not modeled by the mapmaking procedure. The first filter eliminated
pick-up from the mount azimuth motor current by removing the first five
harmonics of each full $1440^{\circ}$ (distinguishing between the positive and
negative scan velocity) azimuth scan. The second filter targeted spurious
polarization due to the wind-induced deformation of the thin polypropylene
environmental seal at the base of the forebaffle extension. Beginning in
September 2021, we removed the environmental seal during good weather, and
data taken thereafter did not require this filter. A third filter attempted to
remove any azimuthally synchronous signal (e.g., ground pick-up, electric
pickup, or, for circular polarization, Zeeman-splitting of atmospheric oxygen
in the Earth’s magnetic field). For linear (circular) polarization, this
filter removed fifteen (ten) azimuthal harmonics. These three filters were
estimated piecewise over timescales ranging from 2 to 6 hours depending on the
detector set, the era of the survey, and the polarization states. In other
words, although the harmonic functions removed are defined relative to the
telescope azimuth scan, they are estimated and held fixed over many scan
periods, decreasing their impact on the celestial signal (Section 3.4).
Finally, from August 2017 to the end of Era 1, the data were notch-filtered at
frequencies below the scan frequency due to radio frequency interference from
web cameras installed in the telescope. During Era 2, the cameras were kept
off during observations. Among these filters, the third filter targeting the
azimuthally synchronous pick-up removed the most celestial signal on large
scales. A more complete description of the filters and the spurious signals
they mitigate can be found in L23. Prior to mapping, the covariance of the
filtered demodulated data was estimated by identifying common-mode covariance
through singular value decompositions and modeling the remaining noise power
as covariant between the detector pairs and Stokes states of each feedhorn
(Dünner et al., 2013, L23). In this way, covariance was estimated as a
function of frequency (bins), between detectors, and between linear and
circular polarization. The noise estimation bias due to missing data was
iteratively corrected by filling the gaps with mock data that preserved the
noise covariance through Gibbs sampling (Huffenberger & Næss, 2018, L23). With
this noise covariance, linear and circular polarization maps were made via a
preconditioned conjugate gradient solver with an adapted version of
minkasi111https://github.com/sievers/minkasi (Romero et al., 2020). Once a map
solution was achieved, it was subtracted from the demodulated data, the noise
covariance was re-estimated, and a new map solution was found. This “template
iteration” was used to make the noise covariance estimation independent from
the signal (Dünner et al., 2013; Aiola et al., 2020, L23)—a requirement for
recovering the signal in the maps, especially on the largest angular scales.
L23 found that five template iterations were sufficient to recover $\geq
97.5\%$ of the signal in the filtered demodulated data. The final linear and
circular polarization maps are shown in Figure 1. These are represented in
HEALPix format with a resolution parameter of $N_{\rm side}=128$ (Górski et
al., 2005). We have made maps and associated data products available to NASA
LAMBDA for public release.222https://lambda.gsfc.nasa.gov/product/class/

For map-based comparison to previous measurements, we have developed
simulations that apply the effect of the CLASS mapping pipeline on other maps.
To generate such a simulation, the $Q$/$U$ polarization maps of another
experiment were first smoothed with a beam defined by the ratio of the CLASS
beam, discussed below, and the beam of the other experiment, thus deconvolving
the intrinsic beam from the map and leaving the map smoothed as if observed by
CLASS. Using the CLASS pointing model, the maps were then projected into the
time domain as though they were CLASS demodulated data. The data selection,
filtering, and mapping processes were then applied to form so-called
_reobserved_ maps. For this work, we have used reobserved WMAP K-, Ka-, and
Q-band (Bennett et al., 2013) and the Planck 353 $\mathrm{GHz}$ map (Planck
Collaboration et al., 2020c). When using a reobserved map, we prepend the name
with “r”, e.g., rWMAP K-band or simply rK.

###  3.4 Mapping and Beam Transfer Functions

To evaluate the impact of time-stream filtering on the maps, L23 used
simulations to estimate the harmonic domain mapping transfer function. This
has the form

$C_{\mathrm{\ell,\mathrm{out}}}=F_{\ell\ell^{\prime}}C_{\ell^{\prime},\mathrm{in}},$
(1)

where $C_{\ell^{\prime},\mathrm{in}}$ refers to the $EE$, $BB$, and $VV$
angular power spectra of the unfiltered sky. The transfer function
$F_{\ell\ell^{\prime}}$ relates these to the angular power spectrum as
observed by CLASS, $C_{\mathrm{\ell,\mathrm{out}}}$, accounting for
correlations introduced between multipole moments and between spectra.
Simulations showed that when this transfer function is used with a
pseudo-$C_{\ell}$ estimator (Chon et al., 2004), that the recovered estimate
of the input spectra are unbiased for $\ell>4$. The effect of the mapping
transfer function is to bias signal power low on the largest angular scales.
For linear (circular) polarization, the transfer function was found to leave
approximately 85% (93%), 67% (85%), 34% (47%), 6% (17%) of the signal in the
angular power spectrum at $\ell=30,\,20,\,10,\,\mathrm{and}\,5$, respectively.
See also Figure 14 of L23.

The mapping transfer function shown in L23 assumed a uniform weight within the
CLASS survey region. That form of $F_{\ell\ell^{\prime}}$ is useful when
simulating the impact of the CLASS filtering strategy on the power spectrum of
other experiments, as we have done when simulating the noise of other
experiments in our calibration (Section 6.3) and when studying the spatial
variation in the synchrotron spectral energy density (SED) (Section 7.3). When
correcting for the transfer function bias on CLASS cross spectra, however, a
modified transfer function correction that assumes the non-uniform CLASS
weight is used, cf. spectra in Sections 7.1 and 7.2. The uniform and CLASS
survey weighted mapping transfer functions are both included in this data
release.

The beam window function for Era 1 was estimated from Moon observations as in
Xu et al. (2020). However, for this paper we have updated that estimate to
include data from Era 2 as well as new analysis steps. The updated analysis is
described in Datta et al. (2023), and summarized in Appendix A. The aggregate
beam for the survey combines beam maps from individual detectors weighted by
their relative contribution to the survey map. The aggregate beam is further
symmetrized by combining beams from across the field of view with different
elliptical orientations and by including the effects of boresight rotations
and scan cross-linking in the survey map. The resulting beam has a FWHM of
$1.559^{\circ}\pm 0.006^{\circ}$, a solid angle of $799.4\pm 0.5\,\mathrm{\mu
sr}$, and an eccentricity of $0.07\pm 0.01$, where the errors are formal
statistical errors from the beam modeling. Furthermore, the beam estimates
made separately from Era 1 and Era 2 agree well within uncertainties. The beam
transfer function was estimated as the square of the harmonic transform of the
beam (a.k.a. the beam window function). This function indicates the measured
angular power spectrum will be reduced to 86%, 74%, 49%, 9% of its intrinsic
value at $\ell=20,\,30,\,50,\,\mathrm{and}\,100$, respectively. For more
details on the beam and beam window function estimate, see Appendix A. Finally
we note that when comparing the Moon-derived beam profile to that of
unresolved sources in the survey maps, one must account for the mapping
transfer function, which attenuates the larger scale modes.

##  4 Power spectrum estimation

Angular power spectra for synchrotron signal analysis and all null tests were
estimated by using the pseudo-$C_{\ell}$ method as implemented in PolSpice
(Chon et al., 2004). This method produces non-biased $EE$ and $BB$ power
spectra accounting for varying survey weight across the sky and the effect of
beam convolution. For any given frequency, or _channel_ , $i$, the
polarization weight $w_{i}$ was modeled as

$w_{i}=\left(h_{\mathrm{QQ}}h_{\mathrm{UU}}-h_{\mathrm{QU}}^{2}\right)^{1/2}.$
(2)

Here, $h_{XY}$ are the hits maps; proportional to the inverse per-pixel
variance of the map. The CLASS hits maps were described in L23. Having used
cross spectra exclusively, instrumental noise bias was avoided and systematic
errors not shared between the crossed sources are reduced. When estimating the
power spectra at a single frequency and single experiment, the cross spectrum
was estimated using a data split. For the CLASS data, the base split is
defined in Section 5.1. When evaluating cross spectra for other experiments,
split definitions depend upon the experiment. For WMAP channels, covariance
weighted coadded even and odd-numbered years define the split. For the Planck
30 $\mathrm{GHz}$ channel, no splits were considered; only cross spectra with
other experiments were used.

In the following discussion, it is convenient to adopt shorthand expressions
for the cross spectra. We use the band letter designation K, Ka and Q to
indicate a WMAP channel. The numbers 30, 100, 143, and 353 are used to
indicate the Planck channels. We use $C$ to indicate the CLASS 40
$\mathrm{GHz}$ data. Thus $\mathit{K}\times\mathit{C}$ indicates the cross
spectrum between the full $K$-band maps and the full CLASS maps, and
$\mathit{Ka}\times\mathit{Ka}$ denotes the cross spectrum between the even and
odd years of the WMAP Ka-band.

When estimating the CLASS signal cross power spectra, the filtering described
in Section 3.3 causes power to be removed from the map and thus biases the
resulting spectrum. The biased spectrum is then corrected by

$\hat{C}_{\mathrm{\ell}}=F^{-1}_{\ell\ell^{\prime}}\hat{C}^{\mathrm{polspice}}_{\mathrm{\ell^{\prime}}},$
(3)

where $F_{\ell\ell^{\prime}}$ is the harmonic domain transfer function
described in Equation 1. To correct for bias when estimating the cross
spectrum between CLASS and another experiment, the matrix square root of
$F_{\ell\ell^{\prime}}$ is used to model the required effective combined
transfer function; any imaginary portion of the matrix is very small and
dropped. In cases when portions of the sky were masked, e.g. to block bright
sources or the galactic plane, the initial binary mask was apodized using the
NaMaster333https://github.com/LSSTDESC/NaMaster (Alonso et al., 2019) C2
apodization tool with a $3^{\circ}$ scale. The map weight used for spectrum
estimation was then the product of the weight from Equation 2 with the
apodized mask.

Power spectra were evaluated in celestial coordinates, as these are the native
coordinates of the CLASS maps. The publicly available WMAP and Planck maps
used in this work were first rotated at their full resolution from Galactic
coordinates to celestial and smoothed with an anti-aliasing cosine-like filter
(Benabed et al., 2009)

$f_{\ell}=\begin{cases}1&\text{ for }\ell\leq\ell_{1},\\\
\frac{1}{2}\left[1+\cos\pi\frac{\ell-\ell_{1}}{\ell_{2}-\ell_{1}}\right]&\text{
for }\ell_{1}<\ell\leq\ell_{2}\\\ 0&\text{ for }\ell>\ell_{2},\end{cases},$
(4)

where $\ell_{1}=128$, and $\ell_{2}=3\times\ell_{1}$. The maps were then
downgraded to $\text{N}_{\mathrm{side}}=128$. In cases when it was necessary
to have different channels expressed in a common resolution, the map was
convolved with the beam ratio of the intrinsic beam and the target beam. Note
for synchrotron-dominated sky regions studied in this paper, the $B$-mode
power is a significant fraction of the $E$-mode power. In this case, $E/B$
purification methods are not required and were not applied in the spectra
presented here.

Finally, all binned power spectra follow the same binning method; bins are
defined between $\ell=[5,125)$ with fixed width $\Delta\ell=10$. Binned
spectra use the typical $\ell(\ell+1)$ weight within each bin.

In an effort to quantify the uncertainty in the derived spectra, all estimates
were compiled with an ensemble of simulations. Each simulation was then
processed through the same power spectra estimation method as the data, and
the resulting ensemble was used to directly estimate the covariance of the
spectra. For visualization, the error bars for power spectra are the square
root of the diagonal of the covariance matrix. For any model fit, however, the
full covariance matrix is used. The details of the simulations used for the
covariance estimation vary depending on the quantity being estimated, but they
can broadly be broken down into two categories: instrument noise only and sky
signal plus instrument noise simulations.

The CLASS noise simulations were generated using the Fourier domain noise
model inferred during the mapmaking procedure; see L23. Random time domain
realizations constrained to be consistent with the noise spectrum were then
mapped. Depending on the application, the noise simulation could model the
full CLASS data set or half of the CLASS data, i.e., splits. For Planck, the
300 FFP10 end-to-end channel simulations were used (Planck Collaboration et
al., 2016c). For WMAP simulations, large-scale noise realizations consistent
with the published low-$\ell$ pixel-pixel noise covariance matrix were
combined with high-$\ell$ noise realizations consistent with the polarization
hits (Larson et al., 2011). The WMAP and Planck noise simulation maps were
smoothed with an anti-aliasing filter and downgraded in resolution to match
the CLASS maps. The noise estimated from the CLASS noise simulations is shown
in comparison with the WMAP Q-band and the Planck 44 $\mathrm{GHz}$ full sky
simulated noise spectra in Figure 2. After accounting for filtering, the CLASS
noise was found to improve upon the comparable space-based measurements in the
$10<\ell<100$ range.

Figure 2: Comparison of the map noise between the full sky regions observed by
WMAP and Planck and that of CLASS as measured in the $EE$ cross spectrum. The
CLASS noise spectra are determined via noise-only simulations, for WMAP
through sampling the WMAP covariance matrix, and for Planck through the
published noise simulations. At each frequency, these simulations have been
shown to be representative of the noise in the respective total maps. The
dashed line illustrates the CLASS noise corrected for the transfer function
(TF). The upward curve at higher $\ell$ is a manifestation of the resolution
capability of the CLASS 40 $\mathrm{GHz}$ telescope. CLASS has the lowest
noise in the range $10<\ell<100$.

##  5 Internal Consistency Checks

The internal consistency of the survey and stability of the calibration
procedure were verified by performing a series of null tests. Broadly
speaking, a “null map” is the difference between two CLASS maps made from
nominally independent sets of data. The null test procedure checks the
statistical consistency between (1) the cross power spectra of two independent
CLASS null maps and (2) the cross spectra of null maps created from
simulations of the data. Simulations include a sky model and the statistical
noise simulations described in Section 4. An inconsistency would indicate the
presence of systematic errors that may bias the CLASS cross spectra results.
No such inconsistency was found.

###  5.1 Null Test Method

The time-ordered data spans (Section 3.1) were first distributed into two base
sets, $m_{1}$ and $m_{2}$, which we call the _base split_. The span
distribution was performed controlling for the boresight rotation of the
telescope, azimuth scanning style, presence of the environmental closeout
film, blackening state of the telescope and baffle structures, and observation
Era. Each base split consists of approximately half of the full mapped data
and includes each distinct observational condition in nearly the same
proportion as in the full survey map. The null maps of the base split itself
are shown in Figure 3.

Figure 3: The base split, $m_{1}-m_{2}$, null maps for the $Q$/$U$ data are
shown. For clear visualization, the maps have been smoothed with a $3^{\circ}$
FWHM Gaussian beam. The amplitude of the noise variations visually track the
variations in the survey depth, see the hit maps in L23. While not used for
the suite of null-tests, these null maps demonstrate the overall overall
calibration stability for the base split — the null maps are dominated by
noise.

The base splits were further subdivided into sub-splits, $A$ and $B$,
according to other criteria, e.g., distinct detector sets, defining four maps,
$m_{i\rho}$ with $i\in\\{1,2\\}$ and $\rho\in\\{A,B\\}$—each map being formed
through the same full mapping procedure of L23.444The data splits were
performed after the filtering of the demodulated data and before the noise
modeling and map-making. Therefore, temporal splits like East/West (see
definitions below), were filtered in the same way as the combined data. The
differences (null maps), $\Delta m_{i}=m_{iA}-m_{iB}$, null the sky signal
leaving a biased noise map. The full list of splits is enumerated below. An
example of a typical set of null maps is shown in Figure 4.

Figure 4: Example null map constructions are shown. The top row shows the
Stokes $Q$/$U$/$V$ base null maps $m_{1A}-m_{1B}$ for the 8 hour in/out
split—maps were smoothed with a $3^{\circ}$ FWHM Gaussian beam to make large
scale features more apparent. The middle rows show the equivalent split from a
reobserved WMAP K band map scaled to the CLASS bandpass — the features show
the expected map difference resulting from imbalanced filtering in the splits.
The bottom row is the residual difference between the data and the template,
which is statistically consistent with the noise.

Differences in filtering and sky coverage between the $A$/$B$ splits will bias
the null cross-spectra, and we use simulations to account for these effects.
For each split, a realization of the CLASS noise model for the split in
question was added to a polarization sky signal model projected into the time
domain. The sky model assumed an atmospheric $V$ dipole (Petroff et al.,
2020b) and the Stokes $Q$/$U$ maps from WMAP 9-year Q-band smoothed with a
$1.5^{\circ}$ Gaussian beam (Bennett et al., 2013). The synthetic time streams
were then mapped as described in L23. In principle, one could then compare the
cross power spectra of the null maps $\Delta m_{1}$ and $\Delta m_{2}$ to
analogous spectra from simulations. Such a comparison would fail to compare
possible spatial differences between the null maps based on the data versus
those from simulations. To account for this effect, a null bias template
$\Delta m^{b}_{i}$ was removed from each respective null map to produce a
bias-corrected null map:555While the bias template corrects for the signal
imbalance between the split, the mapping transfer function bias is still
present in the null map and the simulations.

$\Delta m^{c}_{i}=\Delta m_{i}-\Delta m^{b}_{i}.$ (5)

For the simulated null maps, the $\Delta m^{b}_{i}$ templates are the average
of the simulated null maps $\langle\Delta m_{i}\rangle$ in which the noise is
greatly diminished. For the null maps made from data, the $\Delta m_{i}^{b}$
template is made from a signal-only simulation created by smoothing the
synchrotron-dominated WMAP 9-year K-band map (Bennett et al., 2013) with a
$1.5^{\circ}$ Gaussian beam and scaling it by the approximate ratio ($0.17$)
in synchrotron intensity between the WMAP K-band and the CLASS
$40\,\mathrm{GHz}$ band.

Finally, the bias-corrected null spectra from CLASS were compared to the
corresponding ensemble of simulations via the $\chi^{2}$ test

$\chi^{2}=\left(C_{\ell}^{1\times 2}-\langle
C_{\ell}^{\text{sim}}\rangle\right)^{T}\text{Cov}_{\ell\ell^{\prime}}^{-1}\left(C_{\ell^{\prime}}^{1\times
2}-\langle C_{\ell^{\prime}}^{\text{sim}}\rangle\right),$ (6)

where $C_{\ell}^{1\times 2}$ is the binned CLASS null cross-spectra between
$\Delta m_{1}^{c}$ and $\Delta m_{2}^{c}$, $\langle
C_{\ell}^{\text{sim}}\rangle$ is the average of the corresponding simulated
binned null power spectra,666Seventy simulations were used for each data
split. Using additional simulations did not significantly change the results
of the null tests. $\text{Cov}_{\ell\ell^{\prime}}$ is the bin-bin covariance
matrix determined from the simulations, and summation is implied over repeated
indices. We compute the probability-to-exceed (PTE) for each null test as the
fraction of corresponding simulations having a $\chi^{2}$ greater than that of
the data.

One important difference between the null tests described here and the angular
power spectrum analysis of the CLASS data described in Section 7 is the choice
of masking. The null tests are evaluated _with no Galactic mask_. In the
angular power spectrum analysis of the data, we mask regions of brightest
polarized synchrotron emission (Section 6.1). Therefore, the null tests hold
the data and analysis pipeline to the same stringent self-consistency
standards in the bright Galactic center as they do in the synchrotron-diffuse
regions measured by the angular power spectra.

###  5.2 Split Definitions

Here we describe the A/B data splits used in the null test suite and the
potential systematic errors that motivate the test. Two categories of splits
were used. The first category, detector-based splits, divide the $A/B$-split
by choosing subsets of half the available detectors; the full collection of
detector splits are summarized in Figure 5. In addition to testing for
sensitivity to possible instrument based systematic effects, some detector
splits also test for possible errors in VPM modeling. The second category of
splits, temporal splits, divides the $A/B$-split based upon a temporary
condition that evolves over time and verifies the total maps are insensitive
to that evolving state.

Figure 5:  Detector-splits definition of the CLASS 40 $\mathrm{GHz}$
telescope. The detector pairs coupled to the same feedhorn are shown as
semicircles with the $+45^{\circ}$/$-45^{\circ}$ oriented (viewed from the
focal plane to the sky) detectors shown on the left/right and are positioned
by their pointing offset on the sky. In the “readout rows” panels, the
boundaries of the readout columns are delineated by the gray lines and shades.
During the Era 2 deployment, some of the detector readouts are re-mapped to
optimize detector yields; thus the split based on the odd/even readout row is
different for Era 1 and Era 2. The panel above shows the current configuration
(Era 2) and has the detectors with different assignments in Era 1 outlined in
black.

In the following list, we summarize the $A/B$ splits used in the null tests
and describe the systematic errors tested by the specific splits.

#### Focal plane position

The splits labeled top/bot, left/right, radial, horizontal, and vertical each
break the symmetries of the optical system in different ways and therefore
test for beam and VPM based systematic errors that could also follow these
subsets. The top/bot split also tests for sensitivity to the atmospheric
loading, circular polarization emission from the atmosphere, and residual
ground pickup.

The quadrupole split is designed to probe the impact from the wind signal as
discussed in L23. This split is also sensitive to the modeling of the VPM
transfer function.

#### Detector addressing

Two splits were defined based on their multiplexing (MUX) address location.
The odd/even row split is sensitive to the electric crosstalk between the
adjacent detectors in the same column. The half rows split divides the first
and second halves of the detectors within each column based on the row index
and is another check on possible readout effects.

#### Polarization orientation

The plus/minus split divides the detectors according to their polarization
orientation relative to the telescope. Within a pair, this split shows
opposite signs in the demodulated data for T-to-P signal leakage and VPM
synchronous modulated signals. While the total map is less susceptible to
these issues, the null test from this split can examine this effect at the
extreme limit. Single detector Moon observation revealed a faint dipolar
residual T-to-P signal. This effect is shown to be negligible for cosmological
analysis, but its impact from leaking bright foreground temperature signal
around the Galactic plane into polarization is significant for the plus/minus
null test. Following the prescription in L23 Section 5.1.4, templates of this
leakage are built by convolving the WMAP Q-band temperature maps with the
dipole leakage beam inferred from the CLASS Moon-centered maps for the plus
and minus splits. The template is subtracted from the CLASS null splits before
comparing them to the noise simulations.

The VPM syn. split is the product of the left/right split with the plus/minus
split. This split divides the detectors into groups that systematically
measure more/less synchronous VPM emission and test for sensitivity to
residual emission in the maps.

#### Mount orientation

The bs in/out split divides the spans into the boresight set
$\\{-15^{\circ},0^{\circ},15^{\circ}\\}$ and
$\\{-45^{\circ},-30^{\circ},30^{\circ},45^{\circ}\\}$—i.e., the inner versus
outer boresight angles. The split tests for map sensitivity between the less
tilted versus the more tilted instrument configurations. For example, the
cryogenic system could behave differently or there could be systematic
pointing variations. In a similar way, the bs pos/neg split separates spans
based on the sign of the boresight rotation—boresight 0 is omitted. This split
also serves as a diagnostic for changes in the optical alignment of the
telescopes as any drift in the array center position with respect to the mount
boresight will always show up as a separation in either RA or Dec, depending
on whether the drift is in elevation or azimuth. Due to the projection effect,
the leakage from the sky linear polarization to the circular polarization
through errors in the VPM transfer function would have the largest contrast
when comparing the positive and negative boresight maps.

#### Azimuth scanning

The east/west split divides the data according to the orientation of the
telescope boresight with the east half of the split corresponding to the
$0-180^{\circ}$ azimuth range. This split tests sensitivity to residual ground
pickup since the horizon is different between the two directions. For example,
most of the nearby mountain features are in the east, but most population
centers are in the west. The split also tests sensitivity to observations of
the rising versus setting sky, especially in regards to sunrise and sunset.
The harmonic domain, azimuth-based filters (Section 3.3 and L23) will
generally include data from both splits. For this test, we elect to use the
same azimuth filter as defined for the full data and therefore allow the
filter definition to include the constraining power from both splits. The
approach is duplicated in the simulations to which the split is compared. The
az velocity split divides the data by the sign of the azimuth scanning
direction. This split tests for residual systematic errors related to the
azimuth servo motors. Similarly, the half sweep split divides each
$720^{\circ}$ scan of the mount into the first and second $360^{\circ}$ scans
and would be sensitive to the impact of the azimuth servo in a different way.
For example, the first (second) $360^{\circ}$ scan always includes the mount
accelerating after (decelerating before) a turn-around.

#### Sun and Moon splits

The midnight split divides each span at the solar midnight. Several
environmental factors evolve during the night; for example, wind direction
often shifts by $\sim 40^{\circ}$ toward the North before midnight and the
secular warming/cooling cycles of the ground are different when the sun is
setting versus rising. The 8h in/out split divides each span into the middle 8
hours close to the solar midnight and the remaining time close to the
twilight. This is designed to test the effect of the Sun illumination and the
associated environmental change. The diurnal cycle of air temperature, wind
speed, and PWV follows well with local time and peaks in the afternoon at
around 14:00, 16:00, and 18:00 (UT1-4), respectively. Therefore this split
approximately divides the data into times when the temperature, wind, and PWV
are low and high. It is also sensitive to the clouds, as the cloud cover at
the site is observed to decrease during the night and rise again in the
morning. Finally, the Moon split divides each span depending on whether the
Moon elevation is below/above the horizon. This is sensitive to the Moon’s
illumination through the far sidelobes.

#### Half survey

The final split, survey, divides the data before/after 2019 May 21, which is
approximately the half point of the survey. This probes a potential secular
change in the instrument or sky during the six-year survey.

###  5.3 Null Test Result Summary

Using the two bias-corrected null maps derived from Equation 5, the cross
power spectra were computed for each of the splits described above. The
resulting null spectra from the detector-based splits are collected in Figure
6, and the temporal splits are in Figure 7. In each spectrum, the points shown
are the ratio of the null spectra from the maps to the diagonal of the noise
simulation-based covariance matrix. The reduced $\chi^{2}$ value for each
spectrum is consistent with unity and is shown in the lower right corner of
each subplot; this is computed using the full covariance matrix as in Equation
6. The $\chi^{2}$ PTEs for all spectra are collected in Table 5.3. As the full
survey maps combine all of these data, the final results benefit from reducing
any residual systematic by restoring the symmetry that individual tests are
designed to break. In fact, the symmetries, e.g., polarization angle, focal
plane side, etc., and degeneracies, e.g., boresight angle, scan direction,
etc., are _designed_ in the instrument and survey specifically to reduce
possible systematic errors in the combined data set.

The presented null tests are not necessarily independent of one another, and
therefore the collected PTEs are not sampled from a uniform distribution.
Rather than a strict comparison to a uniform distribution, these results
should be viewed as a demonstration of the overall robust rejection of
systematic errors even when the data are split in ways that maximize the
impact of the potential systematic errors probed by each split.

Table 1:  Null test PTE table. The PTE values are tabulated for every null
test split and each of the six polarization (cross) spectra. The cell color
indicates the PTE value—red (purple) for PTE values less (greater) than 0.5,
meaning that the simulation shows less (more) scatter than the data. The
saturation of the color scales with the deviation of the PTE from $0.5$. The
detector splits and the temporal splits are separated by the white space in
the middle. A subscript third significant digit is included in cells otherwise
rounding to 1 or 0.

Split | $EE$ | $BB$ | $EB$ | $VV$ | $VE$ | $VB$
---|---|---|---|---|---|---
top/bot | 0.03 | 0.52 | 0.41 | 0.07 | 0.47 | 0.24
left/right | 0.33 | 0.10 | 0.29 | 0.23 | 0.40 | 0.84
radial | 0.31 | 0.42 | 0.54 | 0.25 | 0.51 | 0.63
horizontal | 0.51 | 0.22 | 0.81 | 0.63 | 0.70 | 0.69
vertical | 0.87 | 0.46 | 0.79 | 0.08 | 0.63 | 0.50
quadrupole | 0.99 | 0.40 | 0.86 | $0.99_{\footnotesize 9}$ | 0.80 | 0.60
MUX halves | 0.26 | 0.96 | 0.78 | $0.00_{\footnotesize 2}$ | 0.45 | 0.74
MUX parity | 0.10 | 0.79 | 0.51 | 0.35 | 0.74 | 0.31
plus/minus | 0.30 | 0.19 | 0.05 | $0.00_{\footnotesize 4}$ | 0.19 | 0.47
VPM syn. | 0.03 | 0.45 | 0.42 | 0.10 | 0.85 | 0.87
bs in/out | 0.01 | 0.81 | 0.43 | 0.58 | 0.11 | 0.23
bs pos/neg | 0.22 | 0.02 | 0.30 | 0.24 | 0.09 | 0.65
east/west | 0.02 | 0.90 | 0.48 | 0.22 | 0.43 | 0.12
az velocity | 0.22 | 0.95 | 0.47 | 0.03 | 0.43 | 0.78
half sweep | 0.24 | 0.50 | 0.54 | 0.83 | 0.48 | 0.03
midnight | 0.50 | 0.34 | 0.10 | 0.59 | 0.01 | 0.28
8h in/out | 0.31 | 0.39 | 0.39 | 0.14 | 0.26 | 0.32
Moon | 0.13 | 0.46 | 0.67 | 0.58 | 0.30 | 0.37
survey | 0.43 | 0.69 | 0.04 | 0.13 | 0.20 | 0.55
Figure 6:  Null test spectra for detector splits. Each row corresponds to the
six null spectra of a null test split. The null spectra (colored points) are
shown as the $\chi$ value defined by the ratio of the null spectra band power
and the diagonal of the simulation covariance matrix. For visualization
purposes, the error bars of all the data points are set to unity. In the null
test, the $\chi^{2}$ statistics are based on the full covariance matrix
(Equation 6), and the reduced-$\chi^{2}$ values are displayed in the lower
right corner of each panel. The corresponding PTE values are collected in
Table 5.3. Figure 7:  Same as Figure 6, but for the temporal splits.

##  6 Comparison and Calibration to Previous Measurements

###  6.1 Masks

For our analysis, we have defined a series of sky masks. The sky masks
considered were defined to mask sequentially larger regions with the brightest
polarized synchrotron (sX) or dust (dX) emission, where X is an integer $0-9$.
For example, s0 masks the smallest region with the most intense synchrotron
polarization. (Details of the Galactic mask construction are provided below.)
Additionally, all masks exclude data outside of the CLASS survey boundary at
declinations $\delta=-76^{\circ}$ and $\delta=30^{\circ}$. There is also a
common source mask, described below.

The sX and dX masks exclude regions of polarized intensity based on the Planck
Commander synchrotron and dust component maps, respectively (Planck
Collaboration et al., 2020e). For both sX and dX masks, we used a fiducial
angular power spectrum,
$C_{\ell}^{\mathrm{fid}}=0.64\,(\ell/10)^{-1}\,\mathrm{\mu K^{2}}$, the form
of which was fine-tuned to produce well defined emission regions as follows.
The Commander maps and $C_{\ell}^{\mathrm{fid}}$ were smoothed with ten
logarithmically spaced scales, from $0.5^{\circ}$ for s0/d0 to $12^{\circ}$
for s9/d9, and pixels in the smoothed map with values greater than the RMS of
the smoothed spectrum were flagged for masking. Flagged (unflagged) pixels
were set to unity (zero), and the resulting binary mask was then smoothed
again with a $5^{\circ}$ Gaussian beam to remove sharp features. This smoothed
mask was then reset to a binary state with a threshold at 0.25. Built in the
native Commander resolution of $N_{\rm side}=2048$, the Galactic masks were
then downgraded to the CLASS resolution of $N_{\rm side}=128$, masking any
lower-resolution pixel if it contains any flagged $N_{\rm side}=2048$ data.
The Galactic masks for s3, s7, and s9 boundaries are shown in purple, teal,
and yellow, respectively, in Figure 8. Because different Commander maps were
used in constructing the maps for each component, the contour index is not
directly relatable between foreground components. For instance, s1 and d1 do
not share common properties such as masked sky fraction.

As seen in Figure 8, the morphologies of the sX masks follow the Galactic
plane with excursions to mid-to-high latitudes in the regions of the polar
spurs. However, the synchrotron intensity distribution has a marked low
intensity region in the western Galactic plane region. The dX masks tend to be
more confined to the plane with less interruption in longitude direction.

Unless otherwise stated, all masks also account for bright sources based on
the Planck 30 and 44 $\mathrm{GHz}$ polarized compact source catalog that fall
within the CLASS survey region (Planck Collaboration et al., 2016b). In cases
where sources have angular separation less than $1^{\circ}$, only the
brightest source of the collection was kept. Of the remaining catalog, the
brightest 70 sources were masked. The source mask has a $3^{\circ}$ diameter
exclusion region centered at each of the selected points, and custom regions
are tailored to surround a few extended sources. Finally, on account of its
exceptionally high intensity, a $5.4^{\circ}$ diameter region is excluded
around Tau A. The source mask is shown in red in Figure 8.

Figure 8: The CLASS 40 $\mathrm{GHz}$ analysis masks s3 , s7 , and s9 are
shown. The masks share the CLASS declination limitation (gray with black
border) and bright-source mask . They differ in terms of the masked Galactic
synchrotron emission region. The background grayscale image shows the
polarization intensity of the Planck component-separated synchrotron map using
the Commander algorithm (Planck Collaboration et al., 2020e).

###  6.2 Color correction

To facilitate the comparison of CLASS maps with other experiments for
astrophysical sources with non-CMB spectra, we provide the color correction
factors based on the bandpass presented in Dahal et al. (2022). Following
Bennett et al. (2013), the measured source temperature in thermodynamic units
is

$T_{\mathrm{CMB}}=\frac{\Delta T_{\mathrm{CMB}}}{\Delta
T_{\mathrm{A}}}\sum_{i=1}^{3}\omega_{i}T_{\mathrm{A}}(\nu_{i}),$ (7)

where $\Delta T_{\mathrm{CMB}}/\Delta T_{\mathrm{A}}=1.038$ is the antenna-to-
thermodynamic unit conversion factor at 38 $\mathrm{GHz}$ (Dahal et al.,
2022)–the approximate band center of the CLASS 40 $\mathrm{GHz}$ channel;
$T_{\mathrm{A}}(\nu_{i})$ are the source brightness temperatures evaluated at
fixed frequencies, and $\omega_{i}$ are the weights that encode the shape of
the bandpass function. The weighting here is the Gaussian quadrature
approximation to the actual bandpass integral and is accurate for sources with
a smooth spectral shape over the bandpass. The frequencies $\nu_{i}$ are
somewhat arbitrary and were chosen to even out the weights $w_{i}$. Table 2
summarizes these quantities for the CLASS 40 $\mathrm{GHz}$ maps. Since the
CLASS effective bandpass has changed in Era 2 due to the deployment of an
optical low-pass filter (Dahal et al., 2022), we use the average bandpass
function from Era 1 and Era 2 weighted by their relative statistical weights
in the final maps, see Figure 9.

Figure 9: The era-weighted effective bandpass model of CLASS $40\,\mathrm{GHz}$ is shown. For comparison, the detector average of the WMAP Ka band and WMAP Q band measurements (Jarosik et al., 2003) are also included. Each average bandpass has been normalized to unity as the maximum value. Table 2: 40 $\mathrm{GHz}$ Interpolation Data for Color Correction using Equation 7aaAs explained in the text, these values result from a detailed consideration of the detector bandpass, but they should not be interpreted as bandpass values themselves.. $\nu_{1}$bbFrequencies are given in $\mathrm{GHz}$. The choice of frequencies is arbitrary; therefore, the value reported here should not be taken as the central frequencies or the bandwidth of the detector. |  $\nu_{2}$bbFrequencies are given in $\mathrm{GHz}$. The choice of frequencies is arbitrary; therefore, the value reported here should not be taken as the central frequencies or the bandwidth of the detector. |  $\nu_{3}$bbFrequencies are given in $\mathrm{GHz}$. The choice of frequencies is arbitrary; therefore, the value reported here should not be taken as the central frequencies or the bandwidth of the detector. | $\omega_{1}$ | $\omega_{2}$ | $\omega_{3}$
---|---|---|---|---|---
33.9 | 38.0 | 42.1 | 0.3060 | 0.3448 | 0.3491

###  6.3 Calibration

The initial calibration of the maps, as described in Li et al. (2023), was
based on measurements of the Moon and Jupiter (Appel et al., 2019; Xu et al.,
2020; Dahal et al., 2022). Following this, a final adjustment to the absolute
polarization calibration was made by comparing the CLASS data to the bright
diffuse synchrotron signal in reobserved WMAP Ka and Q bands.

The CLASS $40\,\mathrm{GHz}$, and rWMAP Ka/Q-band maps were smoothed to a
common resolution of $2^{\circ}$ and downgraded to $N_{\rm side}=32$. The
bright source mask was applied prior to smoothing to diminish the impact from
bright sources, and only $N_{\rm side}=32$ pixels derived from more than 5 of
the possible 16 unmasked $N_{\rm side}=128$ sub-pixels were kept for further
usage. To select regions with bright diffuse synchrotron signal but avoid
regions with the largest impact from thermal dust emission, we kept pixels
inside s1 but outside d1, and we denote this mask as s1-d1, see the inset in
Figure 10. Even with the smallest dust mask, d0, applied, a template
subtraction with rPlanck $353\,\mathrm{GHz}$ scaled to $38\,\mathrm{GHz}$
following a modified blackbody spectrum with $\beta_{d}=1.53$ and
$T_{\mathrm{dust}}=19.6\mathrm{K}$ only resulted in a $0.001$ level change in
the calibration value.

The uncertainty per $N_{\rm side}=32$ map pixel for each dataset is the
standard deviation across 200 different combinations of the corresponding
noise simulations and the CMB simulations. For rWMAP, full-sky noise
simulations were generated at $N_{\rm side}=16$ by sampling according to the
full noise covariance matrix.777wmap_band_quninv_r4_9yr_Ka_v5.fits
wmap_band_quninv_r4_9yr_Q_v5.fits White noise realizations were created at
$N_{\rm side}=128$ by sampling according to the per-pixel Q/U covariance
matrices. The simulations were converted to spherical harmonic coefficients
$a_{\ell m}$, and composite spectra were made by concatenating the full-
covariance $a_{\ell m}$’s at $\ell\leq 48$ with the white noise simulations at
$\ell>48$. The combined spectra were then converted back to $N_{\rm
side}=128$, smoothed to $2^{\circ}$ FWHM, and downgraded to $N_{\rm side}=32$.
For CLASS, details of the noise simulations can be found in L23. These are
masked, smoothed, and downgraded in the same way as the CLASS data. Gaussian
random CMB simulations were generated at $N_{\rm side}=128$ according to the
Planck best-fit $\Lambda$CDM CMB power spectra888COM_PowerSpect_CMB-base-
plikHM-TTTEEE-lowl-lowE- lensing-minimum-theory_R3.01.txt and then masked and
smoothed to $2^{\circ}$ FWHM, and downgraded to $N_{\rm side}=32$. To create
the effect of reobservation, we applied the CLASS harmonic domain mapping
transfer function (Equation 1) to all simulations except the CLASS noise
simulations (already filtered).

The calibration factor was obtained via a two-step fitting. We first estimated
the spectral index $\beta_{s}$ within the region selected by s1-d1 as:

$\hat{\beta}_{s}=\underset{\beta_{s}}{\mathrm{argmin}}\sum_{p}\frac{\left[\tilde{d}_{Y}(p,\nu,\beta_{s})-\tilde{d}_{X}(p,\nu,\beta_{s})\right]^{2}}{\tilde{\sigma}_{Y}^{2}(p,\nu,\beta_{s})+\tilde{\sigma}_{X}^{2}(p,\nu,\beta_{s})},$
(8)

where $d_{X}$, $\sigma_{X}$ are the amplitude and uncertainty of Stokes
$Q$/$U$ measured in band $X$, and the sum is over all pixels $p$ in s1-d1. The
notation $\tilde{x}(\nu,\beta_{s})$ means that the quantity $x$ has been
converted into antenna-temperature units and color corrected to frequency
$\nu$ assuming a power-law spectrum with index $\beta_{s}$. For this fitting,
we used the rWMAP Ka and rWMAP Q bands data color corrected to the same
frequency $\nu=39\,\mathrm{GHz}$ using parameters from Table 20 of Bennett et
al. (2013). (The frequency $\nu$ was chosen arbitrarily.) The best-fit value
is $\hat{\beta}_{s}=-3.23$. In the second step, we began by computing the
effective frequency for CLASS $40\,\mathrm{GHz}$ ($\nu_{\mathrm{eff}}$) using
Equation 7, by setting
$\sum_{i=1}^{3}\omega_{i}(\nu_{i}/\nu_{\mathrm{eff}})^{\hat{\beta}_{s}}=1$.
The effective frequency is $\nu_{\mathrm{eff}}=37.56~{}\mathrm{GHz}$. We
converted CLASS $40\,\mathrm{GHz}$, rWMAP Ka, and Q data into antenna
temperature units and then color corrected the rWMAP Ka and Q band data to
$\nu_{\mathrm{eff}}$. We fit the model $y=kx$ with $y$ as either rWMAP Ka or Q
band and $x$ as CLASS $40\,\mathrm{GHz}$, taking into consideration the
uncertainties on both axes. In these fits, the slope $k$ is the calibration
factor. As a check of the first step, we also fit rWMAP Q against rWMAP Ka
with the expectation that the slope be consistent with unity. The correlations
between pixels and between the Stokes $Q$/$U$ signal within a pixel were not
considered. The slope, the Pearson correlation coefficient $r$, the $\chi^{2}$
of the fit, and the PTE values are shown in Figure 10. In all the fits, a
strong linear relationship is present with $r\geq 0.94$. As expected, the
slope between WMAP Ka and Q bands is consistent with unity. The calibration
factors (slopes) derived from the fits between rWMAP Ka/Q bands and CLASS
$40\,\mathrm{GHz}$ are $\eta_{\mathit{Ka}}=1.177\pm 0.012$ and
$\eta_{\mathit{Q}}=1.179\pm 0.020$, where the smaller error on from the rWMAP
Ka band fit is a result of its higher synchrotron signal-to-noise ratio (S/N).

We express the absolute calibration factor as
$\eta=\eta_{0}\pm\sigma_{\eta}^{\mathrm{(stat)}}\pm\sigma_{\eta}^{\mathrm{(sys)}}$.
The calibration value $\eta_{0}=1.18$ is the inverse-variance weighted average
of $\eta_{Ka}$ and $\eta_{Q}$. The statistical uncertainty
$\sigma_{\eta}^{\mathrm{stat}}=0.01$ is obtained by propagating the
uncertainties on each of the two slopes to the weighted average. We note that
this treatment neglected the correlation caused by the fact that CLASS
$40~{}\mathrm{GHz}$ data was used in both, though the correlation would be
minimal as the CLASS $40~{}\mathrm{GHz}$ data has the smallest uncertainty.
The PTE values of the fitting are slightly low because the intrinsic
scattering of the synchrotron signal starts to matter at this relatively high
S/N region. The uncertainty estimated by bootstrapping has the same value. The
systematic uncertainty $\sigma_{\eta}^{\mathrm{sys}}=0.05$ results from the
$0.5~{}\mathrm{GHz}$ uncertainty on the band center of CLASS 40 $\mathrm{GHz}$
(Dahal et al., 2022), which is the difference in the calibration value when
converting rWMAP Ka and Q bands to antenna temperature at
$\nu_{\mathrm{eff}}\pm 0.5~{}\mathrm{GHz}$. As a final check, we tried all
combinations of sX-dY with X,Y ranging from 0 to 9, and for all cases with
$r>0.9$ and PTE of slopes within $(0.01,0.99)$ (9 out of 100 in total), the
calibration values are consistent with $\eta_{0}=1.18$.

Figure 10: Correlation of Stokes $Q$ (blue) and $U$ (orange) between different
bands. All data are converted to antenna temperature at
$\nu_{\mathrm{eff}}=37.56~{}\mathrm{GHz}$. Left: rWMAP Ka to CLASS
$40~{}\mathrm{GHz}$. Middle: rWMAP Q to CLASS $40~{}\mathrm{GHz}$. Right:
rWMAP Q to rWMAP Ka. The pink solid line is the best fit for the data, and the
black dashed line has a slope equal to 1 for reference. The top left of each
panel shows, from top to bottom, the best fit slope, the Pearson correlation
coefficient $r$, the $\chi^{2}$, and the PTE. The number of degrees of freedom
of the fitting is 417. The slightly low PTE values result from intrinsic
scattering of the synchrotron signal within the selected region — a
bootstrapping analysis confirmed the same uncertainties on the slope. The fact
that the fitted slopes from the left two panels match each other well and that
from the right panel is consistent with 1 means that the frequency spectrum of
the data we selected can be well explained by a single power law with index
$\beta_{s}=-3.23$. The inset within the right panel shows the region selected
by s1-d1 mask (black) at $N_{\rm side}=32$, where the CLASS declination limit
is in gray.

The final calibration factor is $\eta=1.18\pm 0.01\mathrm{(stat)}\pm
0.05\mathrm{(sys)}$. At this time, it is unclear why this polarization
calibration differs from the Moon and Jupiter temperature-based calibrations.
A portion of this final correction includes accounting for the stable, $\sim
3\%$, atmospheric opacity that is not included in the baseline temperature
calibration. The uncertainty of the initial Moon and Jupiter based calibration
is $\sim 2\%$. However the remaining discrepancy is $\sim 14\pm 5\%$.

As a check on the survey-based calibration, we confirmed that the measured
polarized flux of Tau A, using aperture photometry, is consistent with
previous measurements; see Figure 11. We measured the Tau A flux by
integrating the polarization $P=\sqrt{Q^{2}+U^{2}}$ in antenna temperature
within a disk of radius 1.85∘ centered at RA = 83.633∘ and DEC = 22.014∘. The
CLASS beam profile suppresses the Tau A $P$ signal to 0.6% of the peak value
at a radius of 1.85∘. At this radius, the Tau A signal amplitude matches the
background noise level of the CLASS $P$ map. An estimate of the nearby diffuse
Galactic emission was subtracted from the $Q$ and $U$ disk amplitudes to
account for the Galactic signal that the large CLASS beam integrates together
with Tau A (see Figure 13). This baseline correction increased the measured
flux by $\sim 8$%. An additional $\sim 4$% flux was added to account for the
beam solid angle outside the disk radius (see Appendix A). The integrated Tau
A antenna temperature (in units of K-sr) was divided by
$\Gamma\Omega=c^{2}/(2k\nu_{\mathrm{eff}}^{2})$ to obtain the Tau A flux in
units of Janskys (Jy). Finally, we note that the survey maps were produced
assuming a VPM transfer function based on the flux density of the sky having a
spectral index $\alpha=-1$ as opposed to the known Tau A index $\alpha=-0.36$.
To correct for this difference, based on simulations testing the impact of
this effect, the map-based CLASS Tau A polarization flux was divided by 1.04.
We applied the same aperture photometry method to re-observed WMAP and Planck
maps and found Tau A flux values consistent with the published values.
Simulations showed that the CLASS mapping transfer function (Equation 1)
diminishes flux density measurements of unresolved sources (like Tau A) by
$\leq 1\%$. For simplicity, a correction of this bias was not applied to this
Tau A measurement.

Figure 11: The measured Tau A polarization flux of ${21.6\pm
1.2\,\mathrm{Jy}}$ at an effective frequency of 38.1 $\mathrm{GHz}$ is
consistent with previous measurements. The error on the CLASS Tau A flux
measurement is driven by the uncertainty on the map absolute calibration
factor. The polarization flux densities measured by WMAP (Weiland et al.,
2011) and Planck (Planck Collaboration et al., 2016b) are referenced to the
year 2019 by correcting for the expected secular decrease in Tau A flux
(Weiland et al., 2011). The blue line shows the best fit WMAP Tau A
polarization flux density model (referenced to 2019) as a function of
frequency. The shaded region is the 1$\sigma$ contour of the model’s flux
prediction including spectral and time evolution uncertainty.  Figure 12: A
comparison of the CLASS 40 $\mathrm{GHz}$ survey maps (middle column) to the
rWMAP K-band (left column) and rWMAP Q-band (right column) channels. Stokes
$Q$ ($U$) maps are in the top (bottom) row. All bands are smoothed with a 1.5∘
Gaussian beam for visualization purposes. The synchrotron foreground is much
brighter at the lower frequency (23 $\mathrm{GHz}$) K-band, so the maps were
scaled (by $0.2$) to roughly match the CLASS amplitude. The agreement between
the high signal rWMAP K-band map and the CLASS map is striking—as is the
sensitivity improvement relative to rWMAP Q-band. The CLASS color scale is the
same as Figure 1.

###  6.4 Relation to Previous Measurements

The CLASS 40 $\mathrm{GHz}$ survey complements the collection of available
surveys targeting the synchrotron foreground component of the polarized
millimeter sky—with the WMAP Q-band ($41\,\mathrm{GHz}$, Bennett et al., 2013)
being the closest in frequency and sky footprint to the maps presented in this
work. To the extent that the emission is dominated by synchrotron, the SED is
well described by a power law in antenna temperature units with index
$\beta\approx-3.1$; though this value is known to vary over the sky (Weiland
et al., 2022), a point we consider in Section 7.3. Taking advantage of this
spectral “lever arm,” we will also use WMAP K-band ($23\,\mathrm{GHz}$,
Bennett et al., 2013) data with their high S/N polarized synchrotron emission
to qualitatively validate the CLASS maps in this section.

Figure 12 shows full-sky comparisons between CLASS polarization maps and
reobserved WMAP K\- and Q-band maps. The rWMAP K-band maps have been scaled by
a factor of 0.2 to approximately reproduce the expected signal level in the
CLASS maps. No scaling has been applied to rWMAP Q-band. Note that for the
WMAP data, the reobservation process already smooths the map with a
$1.5^{\circ}$ beam—no additional smoothing was performed to those maps. For
visualization purposes, the CLASS maps have been smoothed with a $1.5^{\circ}$
beam. Smoothing the CLASS map has the effect making the noise comparable
between the three maps, but the signal in the CLASS maps are smoothed to a
slightly lower resolution — this effect is too small to be visible in the
figures. There are two main qualitative conclusions to draw from Figure 12.
First, the CLASS maps have significantly less noise than the rWMAP Q-band
maps. This is consistent with angular noise spectra of CLASS being the lowest
in the range $10<\ell<100$ in Figure 2. The second is that the synchrotron
emission qualitatively matches the rWMAP K-band signal on the largest angular
scales.

Figures 13 and 14 provide a zoomed-in comparison between the CLASS
polarization maps to the rWMAP K-band data - no further image processing is
performed. Figure 13 gives a view of the Galactic plane in the range
$-25^{\circ}<b<25^{\circ}$. Overall there is excellent visual agreement
between the two data sets. Some prominent Galactic features, such as Tau A at
$(\ell,b)\approx(185^{\circ},-5^{\circ})$ and the Galactic Center, appear
discrepant due to frequency spectra that differ significantly from the overall
spectra implicit in the applied scaling. Otherwise, the overwhelming majority
of features, including familiar sources such as the Gum Nebula at
$(\ell,b)\approx(260^{\circ},0^{\circ})$ and the Centaurus A radio galaxy at
$(\ell,b)\approx(310^{\circ},20^{\circ})$, match in morphological detail.
Figure 14 is a more detailed look at a subset of diffuse polarized features in
the maps. The amplitude of the features is approximately $5-10\,\mathrm{\mu
K}$. At this brightness contrast, the apparent discrepancies between the maps
on the smallest scales is due to the rWMAP K-band data having lower noise
after being rescaled by the 0.2 factor to match the CLASS synchrotron
brightness. On larger scales, we see agreement between CLASS and rWMAP for
extended lower-surface brightness features extending $40^{\circ}$ across the
maps.

Figure 13: A comparison of the CLASS 40 $\mathrm{GHz}$ survey maps to the re-
observed WMAP K-band along the Galactic plane. Stokes $Q$ ($U$) maps are in
the top (bottom) two panels. The maps have been smoothed with a 1.5∘ Gaussian
beam for visualization purposes, and the rK-band map has been scaled by a
factor of 0.21 to enable direct visual comparison with the same color scale.
Given the diversity of astrophysical processes contained within this map and
the substantial difference in frequencies at which the signals are measured,
the differences near compact sources and in some diffuse regions can be
expected. Overall the CLASS result shows remarkable similarity to rWMAP,
despite being from a ground based observation. Figure 14: A comparison of the
CLASS (bottom row) and rWMAP K-band (top row) maps in an area of relatively
diffuse emission centered on $(\ell,b)=(22^{\circ},-50^{\circ})$. The columns,
from left to right, show the Stokes $Q$, $U$, and polarization intensity.
CLASS has replicated the rWMAP space measurement and with improved sensitivity
using a ground-band observatory. This comparison includes lower surface
brightness ($5-10\,\mathrm{\mu K}$ at 40 $\mathrm{GHz}$) features extending
$40^{\circ}$ across the map.

##  7 Galactic synchrotron analysis

Polarized diffuse Galactic emission predominantly arises from two physical
mechanisms: synchrotron emission from relativistic electrons spiraling about
Galactic magnetic field lines, and thermal emission from dust grains that are
preferentially aligned within that field. Emission from these components
significantly surpasses the polarized CMB signal on large angular scales, and
thus accurate removal of the Galactic emission is essential to studies of the
CMB polarization. In general, polarized foreground emission from synchrotron
dominates over that of thermal dust at frequencies $\nu\lesssim 60$
$\mathrm{GHz}$, although the exact mixture is dependent on sky position and
angular scale (Planck Collaboration et al., 2020e; Weiland et al., 2022). As
the primary component for CLASS 40 $\mathrm{GHz}$, we examine behavior of
synchrotron emission in both harmonic and pixel space.

In units of antenna temperature, synchrotron emission is characterized as a
power law, with $T\propto\nu^{\beta_{s}}$, where $\beta_{s}$ is the spectral
index. Previous studies using data from, e.g., WMAP, Planck, S-PASS, and
QUIJOTE have firmly established spatial variations of $\beta_{s}$ with a mean
value near $-3.1$ (e.g., Planck Collaboration et al., 2020e; Krachmalnicoff et
al., 2018; Rubiño-Martín et al., 2023; de la Hoz et al., 2023; Fuskeland et
al., 2014; Weiland et al., 2022), though the degree to which $\beta_{s}$
varies has not yet been well characterized over the full sky.

In this section, we consider how the CLASS measurement contributes to
furthering the characterization of the synchrotron foreground component.

###  7.1 Synchrotron Power Spectra

The distribution of synchrotron radiation is highly anisotropic—emission is
dominated by structures within the Galactic plane—and non-Gaussian, and
therefore suffers information loss when compressed to a power spectrum
representation. Nevertheless, the power spectrum can be computed and formally
interpreted as a representative measurement of how the angular power varies
within the region of interest, which, for CMB studies, generally excludes the
brightest emission from the plane and other Galactic features that are the
most anisotropic and non-Gaussian.

When evaluating properties of the synchrotron signal, this analysis used only
the pseudo-$C_{\mathrm{\ell}}$ estimator (PolSpice, Chon et al., 2004) to
estimate the cross power spectrum over the $\ell$-range (5, 125). While less
optimal than a quadratic estimator (e.g., xQML, Vanneste et al., 2018) for
$\ell<20$, there are fundamental difficulties in combining the two estimators
in a single spectrum. For example, a low-resolution mask appropriate for xQML
evaluation that is exactly upgraded to the full map resolution has sharp edges
which cause excess mode mixing in the PolSpice spectrum estimation. Normally
this pseudo-$C_{\ell}$ issue is ameliorated by apodizing the mask
edges—fundamentally decreasing the power relative to the low-resolution map
version. The estimator could adjust for the mask distinction only if the
underlying field were isotropic and purely Gaussian. Rather than compromising
the consistency of the measured power between the estimators, we simplify this
analysis by using only the single PolSpice estimator at the expense of less
optimal uncertainties for the first two band powers.

The statistical properties of the CLASS 40 $\mathrm{GHz}$ maps were summarized
by estimating the angular power spectrum with the s9 mask applied, as
described in Section 6.1. Following the cross spectrum estimation approach
described in Section 4, the $EE$ and $BB$ power spectra are shown in Figure
15. For comparison, the the cross spectrum for the WMAP Q\- and Ka-bands
within the same sky region is also shown. The spectra are in antenna
temperature units assuming a constant spectral index $\beta_{s}=-3.1$ and
color corrected to a common frequency of 38 $\mathrm{GHz}$ using the factors
in Table 2 for CLASS and the analogous table in Bennett et al. (2013),
Appendix E, for WMAP. To isolate the synchrotron component, the Planck best-
fit CMB angular power spectrum (Planck Collaboration et al., 2020a) was
subtracted from the $EE$ spectra prior to color correction. The resulting
uncertainty is estimated using the respective noise simulations described in
Section 4. For CLASS, noise simulations for the base data split were randomly
paired to form an ensemble of 10000 CLASS noise cross-spectra simulations —
individual split simulations were allowed to be reused, but no pair was
duplicated. The WMAP Q/Ka-band noise cross spectra was estimated by generating
500 WMAP cross-spectra simulations with independent noise realization for the
even/odd year split as described in Section 4. The uncertainties on the
spectrum are the square root of the the diagonal of the inferred binned
covariance matrix. The measured CLASS spectra are in good agreement with the
WMAP Q-band result within the same sky region. Differences between the spectra
might result from the bandpass differences between the two channels–recall the
WMAP Q-band approximately includes 38–46 $\mathrm{GHz}$ (8 $\mathrm{GHz}$
bandwidth, Page et al., 2003), the WMAP Ka-band is 28–37 $\mathrm{GHz}$ (6.9
$\mathrm{GHz}$ bandwidth, Page et al., 2003), while the CLASS 40
$\mathrm{GHz}$ channel uses a somewhat larger 32.8-43.6 $\mathrm{GHz}$ (10.9
$\mathrm{GHz}$ bandwidth) range (Dahal et al., 2022). See also Figure 9. Given
the anisotropic nature of the foreground signal, the differences between the
map weights of each experiment could manifest as somewhat different resulting
spectra as well.

The measured spectra follow the expected overall trend found in the Commander
synchrotron model (Planck Collaboration et al., 2020e) which is shown as the
shaded band in Figure 15. Note that these spectra, outside the s9 mask, are
very diffuse and yet polarized synchrotron emission still dominates the
polarized CMB signal. The new CLASS measurement is seen to improve upon the
measured spectra from the nearest frequency space-based measurement.

Figure 15: _Top panels._ The angular cross power spectra of the dominant
synchrotron component of the CLASS 40 $\mathrm{GHz}$ data are shown in units
of $\mathrm{\mu K}^{2}$ antenna temperature within the s9 mask. For
comparison, the synchrotron component of the WMAP Q-band and WMAP Ka-band
(Bennett et al., 2013) cross-spectra between the even and odd years are also
shown. Each spectrum has been color-corrected to 38 $\mathrm{GHz}$ assuming a
fixed synchrotron spectral index of $\beta_{s}=-3.1$. The shaded swaths are
bounds of simulations of synchrotron $EE$ and $BB$ models from the Planck
Commander component analysis (Planck Collaboration et al., 2020e) which has
been referenced to $38\,\mathrm{GHz}$ also assuming $\beta_{s}=-3.1$. No
additional correction was applied for the differing sky regions used for the
Planck model and that used for these spectra. Since the s9 mask was designed
to block a large fraction of the anisotropic synchrotron power that was
included in the region described by the Planck model, it is expected to find
the CLASS and WMAP spectra to fall somewhat lower than the model. The Planck
best-fit $EE$ angular power spectrum (Planck Collaboration et al., 2020a) was
subtracted from the $EE$ data sets as plotted to isolate the synchrotron
components. _Bottom panels._ The uncertainty-normalized difference between the
CLASS 40 and WMAP-Ka/Q band.

###  7.2 Synchrotron SED

The SED of the diffuse synchrotron radiation has been measured to follow a
power law with index $\beta_{s}\approx-3$ with small differences depending
upon the sky region and angular scale measured, cf. Krachmalnicoff et al.
(2018); Martire et al. (2022); Planck Collaboration et al. (2020e); Choi &
Page (2015). In this section we explore the inclusion of the CLASS data to
those publicly available data dominated by synchrotron, but not significantly
impacted by polarization decorrelation and rotation from Faraday rotation
effects. To this end, we considered data from WMAP K-, Ka-, and Q-bands,
Planck 30 $\mathrm{GHz}$, and the CLASS 40 $\mathrm{GHz}$ channel.

In the sky regions where polarized synchrotron radiation is faint, the
harmonic domain description of the data can be used to quantify the brightness
for a given channel. Following the conventions from Krachmalnicoff et al.
(2018) and Rubiño-Martín et al. (2023), we modeled the diffuse synchrotron
emission as

$C^{XY}_{\ell}=A^{XY}_{40}\left(\frac{\ell}{40}\right)^{\alpha_{XY}},$ (9)

where $\\{X,Y\\}$ indicated the data sets being crossed. Provided the index
$\alpha$ are consistent between different channels, the amplitude, $A_{40}$,
measures the synchrotron SED at the effective frequency
$\nu_{\mathrm{eff}}=\sqrt{\nu_{X}\nu_{Y}}$. The frequencies $\nu_{X(Y)}$ are
the effective frequency for each individual channel accounting for the finite
bandwidth of the detectors and the shape of the input spectrum.

The SED for each of the $E$ and $B$ components are assumed to follow a simple
model

$D^{\nu_{1}\times\nu_{2}}=A_{0}\left(\frac{\nu_{1}\nu_{2}}{\nu_{0}^{2}}\right)^{\beta_{s}},$
(10)

where we have adopted the pivot frequency $\nu_{0}=30$ $\mathrm{GHz}$ as it is
near the geometric mean of the range of effective frequencies considered.

Measurement of the SED followed two methods. In the first method, the power
law in Equation 9 was fit to each cross spectra $C^{XY}$ for each $EE$ and
$BB$ pair. As explained in Section 4, power spectra covariance are estimated
using simulations. For spectra not including CLASS data, each map simulation
comprised of the noise sampled as described in Section 4, CMB realizations
following the Planck best fit $\Lambda$CDM model (Planck Collaboration et al.,
2020a), and Gaussian realizations of a fiducial synchrotron model.999The exact
parameters of this model would affect the fitting result through its
contribution to the sample variance. We have iterated this procedure with the
best-fit results from the data and find it converged within five steps Power
spectra were estimated for these map simulations in the same way as the data;
therefore, their covariance naturally incorporates the statistical variance
from the map noise, the signal sample variance, and the inter-bin covariance
from the estimator due to mode coupling.

For CLASS data, the noise simulations provided by the pipeline (L23) represent
the map noise after filtering. To correct for the mapping transfer function,
the simulated cross power spectra were computed separately for the signal,
noise, and cross terms; taking advantage of the associative property of the
two-point statistic. The harmonic transfer function correction was only
applied to terms where the CLASS noise simulations were involved — there is no
need for signal term corrections as the filtering was not included in the
signal simulations. For cross spectrum terms that included map pairs which
both required transfer function correction, the spectrum was corrected with
$F^{-1}_{\mathrm{\ell\ell^{\prime}}}$, and for terms with only one CLASS noise
map, the correction was done with $F^{-1/2}_{\mathrm{\ell\ell^{\prime}}}$.

The resulting $A_{40}$ amplitudes and uncertainties for the s3, s7, and s9
masks are collected in Table 3 and shown in Figure 16. The cross spectra
including CLASS data are indicated with the open circled, relatively larger,
and bolder colored points in the plot.

Table 3: $A_{40}$ values for cross-spectra synchrotron SED  |  $EE$ [$10^{2}\times\mathrm{\mu K}^{2}$] |  $BB$ [$10^{2}\times\mathrm{\mu K}^{2}$]
---|---|---
ChannelaaThe cross spectrum channel is indicated using an abbreviation of the channels involved. The abbreviations $K$, $Ka$, $Q$ are used for the eponymous WMAP bands, the Planck 30 GHz channel is denoted 30, and the CLASS 40 $\mathrm{GHz}$ channel is denoted $C$. |  $\nu_{\mathrm{eff}}$ [$\mathrm{GHz}$] | s3 | s7 | s9 | s3 | s7 | s9
K$\times$K | 22.4 | 10.36 $\pm$ 0.76 | 5.52 $\pm$ 0.56 | 4.12 $\pm$ 0.45 | 4.64 $\pm$ 0.65 | 2.34 $\pm$ 0.49 | 1.47 $\pm$ 0.39
K$\times$30 | 25.2 | 5.47 $\pm$ 0.31 | 3.14 $\pm$ 0.21 | 2.34 $\pm$ 0.17 | 2.43 $\pm$ 0.23 | 1.20 $\pm$ 0.17 | 0.71 $\pm$ 0.13
K$\times$Ka | 27.1 | 3.23 $\pm$ 0.17 | 1.84 $\pm$ 0.15 | 1.44 $\pm$ 0.15 | 1.75 $\pm$ 0.16 | 1.08 $\pm$ 0.13 | 0.70 $\pm$ 0.13
K$\times$C | 29.2 | 2.07 $\pm$ 0.12 | 1.23 $\pm$ 0.06 | 0.91 $\pm$ 0.05 | 0.93 $\pm$ 0.09 | 0.41 $\pm$ 0.05 | 0.23 $\pm$ 0.04
K$\times$Q | 30.1 | 1.89 $\pm$ 0.14 | 1.12 $\pm$ 0.10 | 0.91 $\pm$ 0.12 | 0.72 $\pm$ 0.12 | 0.37 $\pm$ 0.09 | 0.20 $\pm$ 0.10
30$\times$Ka | 30.5 | 1.71 $\pm$ 0.09 | 1.03 $\pm$ 0.08 | 0.80 $\pm$ 0.07 | 0.86 $\pm$ 0.08 | 0.50 $\pm$ 0.07 | 0.25 $\pm$ 0.06
Ka$\times$Ka | 32.7 | 0.94 $\pm$ 0.12 | 0.52 $\pm$ 0.15 | 0.49 $\pm$ 0.15 | 0.55 $\pm$ 0.13 | 0.39 $\pm$ 0.13 | 0.18 $\pm$ 0.12
30$\times$C | 32.9 | 1.06 $\pm$ 0.06 | 0.61 $\pm$ 0.03 | 0.45 $\pm$ 0.03 | 0.50 $\pm$ 0.04 | 0.20 $\pm$ 0.02 | 0.14 $\pm$ 0.02
30$\times$Q | 33.8 | 1.02 $\pm$ 0.06 | 0.66 $\pm$ 0.05 | 0.54 $\pm$ 0.07 | 0.46 $\pm$ 0.05 | 0.24 $\pm$ 0.04 | 0.17 $\pm$ 0.05
Ka$\times$C | 35.2 | 0.66 $\pm$ 0.05 | 0.42 $\pm$ 0.03 | 0.32 $\pm$ 0.03 | 0.31 $\pm$ 0.04 | 0.16 $\pm$ 0.03 | 0.10 $\pm$ 0.02
Ka$\times$Q | 36\. 3 | 0.61 $\pm$ 0.08 | 0.36 $\pm$ 0.06 | 0.34 $\pm$ 0.06 | 0.19 $\pm$ 0.06 | 0.12 $\pm$ 0.05 | 0.05 $\pm$ 0.04
C$\times$C | 38.0 | 0.56 $\pm$ 0.05 | 0.32 $\pm$ 0.04 | 0.23 $\pm$ 0.03 | 0.24 $\pm$ 0.04 | 0.14 $\pm$ 0.03 | 0.08 $\pm$ 0.03
C$\times$Q | 39.1 | 0.39 $\pm$ 0.03 | 0.24 $\pm$ 0.03 | 0.19 $\pm$ 0.03 | 0.16 $\pm$ 0.02 | 0.09 $\pm$ 0.02 | 0.06 $\pm$ 0.02
Q$\times$Q | 40.3 | 0.44 $\pm$ 0.10 | 0.40 $\pm$ 0.10 | 0.35 $\pm$ 0.10 | 0.15 $\pm$ 0.08 | 0.08 $\pm$ 0.08 | 0.07 $\pm$ 0.09
Figure 16: Measurement of the synchrotron SED for the $EE$ ($BB$) modes are
shown on the left (right). Colors indicate the mask applied when evaluating
the signal with s3 , s7 , and s9 masks. The points and error bars are measured
by fitting a power law model to each individual angular cross spectrum–open
circle, heavier lined points indicate cross spectra that include the CLASS
data. The channel, best fit amplitude, and uncertainties are collected in
Table 3. The curve, and legend values, are from the joint fit of all spectra
together assuming the amplitudes follow the power law SED in Equation 10.

The second method was to fit jointly Equation 10 to the full set of cross
power spectra. In this case, the power law from Equation 9 is still fit to the
angular power spectrum of each frequency, but now the amplitudes were
constrained to follow Equation 10. The resulting best-fit power laws for the
s3, s7, and s9 masks are shown in Figure 16 and show good agreement with the
single frequency measurements. For this fit, the $EE$ and $BB$ spectra were
fit separately. Furthermore, the $168\times 168$ covariance matrix (12 bins
$\times$ 14 cross spectra) was assumed to be block diagonal with only the
internal cross spectra $12\times 12$ sub-block considered to be non-zero. The
best fit power law index, $\beta_{s}$, is shown in the Figure legend. The
ratio of the best fit $BB/EE$ amplitudes is $0.47\pm 0.015$, $0.39\pm 0.02$,
and $0.03\pm 0.02$ for the s3, s7, and s9 mask. Since the shapes of the $EE$
and $BB$ SED profiles differ, this amplitude ratio is frequency dependent. It
still summarizes the relative contributing power in this frequency range
critical for CMB component separation efforts.

In this study, we have not investigated the possible $\beta_{s}$ variation
with angular scale as has been hinted at by other investigations (Choi & Page,
2015; Krachmalnicoff et al., 2018), but that topic will be interesting to
pursue with future work. Here we are focusing on the overall consistency of
the CLASS data with previous satellite measurements and the improved
characterization of the large scale synchrotron properties near 40
$\mathrm{GHz}$.

###  7.3 Pixel-based Spectral Index Estimation

Beyond improving knowledge of the galactic processes involved in generating
the polarized synchrotron emission, improved measurements of the spatial
variation on the spectral index may be essential for high precision component
separation–especially for the detection of the primordial $B$-modes. For
example, Errard & Stompor (2019) found that foreground cleaning with a single
spectral index could result in spurious detection of the tensor-to-scalar
ratio of order 0.01.

In this Section, we present a polarized synchrotron spectral index
($\beta_{s}$) map made by analyzing the CLASS 40 $\mathrm{GHz}$ data together
with the rWMAP $K$-band data. The $K$-band data were chosen due to their high
S/N measurement of synchrotron emission and the fact that they are not
significantly modified by Faraday rotation and depolarization effects in the
plane of the Galaxy. Using the _reobserved_ WMAP data mitigates the bias on
$\beta_{s}$ that would result from comparing the raw (unbiased) WMAP data to
the filtered CLASS data. (The filtering both decreases and redistributes the
polarization signal in the CLASS data.) We investigate potential residual bias
from the reobservation processing below.

In addition to the CLASS declination limit, the bright-source mask described
in Section 6.1 was applied to both maps at $N_{\rm side}=128$ to diminish the
bias on $\beta_{s}$ from small bright regions with different SEDs. Similar to
the process used in Section 6.3, the masked maps were smoothed to $2^{\circ}$
FWHM resolution and downgraded to $N_{\rm side}=32$, where only $N_{\rm
side}=32$ pixels derived from more than five unmasked $N_{\rm side}=128$ sub-
pixels were kept for the analysis. To mitigate the impact of polarized thermal
dust emission in the CLASS ${40\,\mathrm{GHz}}$ maps, we scaled and subtracted
the rPlanck 353 $\mathrm{GHz}$ data from the CLASS ${40\,\mathrm{GHz}}$ data
following two strategies. In the first, we model the dust as a modified black
body spectrum with a uniform $\beta_{d}=1.53$ and $T_{d}=19.6~{}\mathrm{K}$
across the sky (Planck Collaboration et al., 2020e), and in the second we
consider a model with a spatially varying frequency-scaling model. The dust
contribution to the polarized signal in the rWMAP K-band is expected to be
negligible.

Simulations were used to understand the uncertainty on the fitted $\beta_{s}$.
The rWMAP K-band noise simulations were generated and processed in the same
way as was described in Section 6.3 for rWMAP Ka-/Q-bands.101010With
wmap_band_quninv_r4_9yr_K_v5.fits We used the same CMB simulations (filtered
by the CLASS harmonic domain mapping transfer function) and CLASS
${40\,\mathrm{GHz}}$ noise simulations as in Section 6.3. All the simulations
were converted to antenna temperature units according to Table 4.

Table 4: Frequencies, color correction, and thermodynamic-to-antenna unit
conversion ($\Delta T_{\mathrm{A}}/\Delta T_{\mathrm{CMB}}$) factors.

|  $\nu_{K}^{s}$ |  $\nu_{\mathrm{40}}^{s}$ |  $\nu_{\mathrm{353}}^{d}$
---|---|---|---
Center Freq. [$\mathrm{GHz}$] | 22.8 | 40 | 353
Color correction | 0.9496 | 0.8243 | 0.9036
$\Delta T_{\mathrm{A}}/\Delta T_{\mathrm{CMB}}$ | 0.9867 | 0.9634 | 0.0773

Note. — From left to right: frequency for synchrotron at WMAP K-band
$\nu_{K}^{s}$, synchrotron at CLASS 40 $\mathrm{GHz}$ $\nu_{\mathrm{40}}^{s}$
and thermal dust at Planck 353 $\mathrm{GHz}$ $\nu_{353}^{d}$. We used Table
20 in Bennett et al. (2013) to compute the factors for WMAP, Table 2 for CLASS
and the public code fastcc (Peel et al., 2022) to compute the color correction
factor for Planck. When computing the color correction factors in the table,
$\beta_{s}=-3.1$ is assumed for synchrotron, and $\beta_{d}=1.53$,
$T_{d}=19.6~{}\mathrm{K}$ is assumed for thermal dust.

The $\beta_{s}$ values were fit in the Galactic coordinate system at
resolution $N_{\rm side}=8$ ($7.3^{\circ}$ pixels) following the linear
correlation method that has been adopted in several previous works (Fuskeland
et al., 2014, 2021; Weiland et al., 2022). Within each $N_{\rm side}=8$
“super-pixel”, an $x-y$ scatter plot was formed using the contained unmasked
$N_{\rm side}=32$ Stokes $Q$ and $U$ data from rWMAP K-band ($\bm{d}_{x}$) and
CLASS ${40\,\mathrm{GHz}}$ ($\bm{d}_{y}$). The plot was fit using a line with
slope $k$ and zero intercept. The estimation of the slope in each $N_{\rm
side}=8$ super-pixel was obtained by total least squares fitting as

$\hat{k}=\underset{k}{\mathrm{argmin}}\left[\bm{d}(k)^{T}\mathsf{\Sigma}^{-1}\bm{d}(k)\right],$
(11)

where $\bm{d}(k)\equiv\bm{d}_{y}-k\bm{d}_{x}$, and $\mathsf{\Sigma}\equiv
k^{2}\mathsf{N}_{x}+\mathsf{N}_{y}$ is the covariance matrix. The
$\mathsf{N}_{x}$ and $\mathsf{N}_{y}$ are the full Q/U covariance matrices
obtained from rWMAP K-band noise and CMB + CLASS $\mathrm{40~{}GHz}$ noise
simulations, respectively. The uncertainty on the slope $\Delta k$ was
obtained by taking half of the difference between 0.16 and 0.84 percentiles of
the likelihood function
$\mathcal{L}\propto\exp\left[-\frac{1}{2}\bm{d}(k)^{T}\mathsf{\Sigma}^{-1}\bm{d}(k)\right]$.

We converted the fitted slope $\hat{k}$ to $\hat{\beta}_{s}$ as

$\hat{\beta}_{s}=\frac{\log{\hat{k}}}{\log(\nu_{\mathrm{40}}^{s}/\nu_{K}^{s})},$
(12)

where $\nu_{\mathrm{40}}^{s}$ and $\nu_{K}^{s}$ can be found in Table 4. We
began with the color correction factors computed at $\beta_{s}=-3.1$ across
the sky and then performed the fitting iteratively with the updated
$\hat{\beta}_{s}$ values for color correction per $N_{\rm side}=8$ pixel until
the fractional difference between two iterations is smaller than 0.01%.
Convergence was achieved within 5 steps for all the values shown in Figure 17.
The scaled rPlanck 353 $\mathrm{GHz}$ maps were removed from the corrected
CLASS ${40\,\mathrm{GHz}}$ data in each iteration. Because the spectral slope
of the dust is not changed during the fitting, the rPlanck 353 $\mathrm{GHz}$
data are converted to antenna temperature according to Table 4 only once.

Figure 17:  Top: the $\hat{\beta}_{s}$ map at $N_{\rm side}=8$ resolution,
fitted with rWMAP K-band and CLASS $40~{}\mathrm{GHz}$ data. Middle: the total
uncertainty on the spectral index $\Delta\beta_{s}^{\mathrm{tot}}$. Bottom:
the PTE values. Regions with $\Delta\beta_{s}^{\mathrm{tot}}>0.2$ and those
masked by the declination limit (light gray) and point-source mask are
excluded. In the middle and bottom panels, pixels with $\mathrm{PTE}<0.01$ are
highlighted by green borders.

Figure 17 shows the $\beta_{s}$ map, its total uncertainty
$\Delta\beta_{s}^{\mathrm{tot}}$ map, and a map of the PTE associated with the
linear fit. Regions with $\Delta\beta_{s}^{\mathrm{tot}}>0.2$ and those masked
by the declination limit (light gray) and the bright-source mask are excluded.
In the $\Delta\beta_{s}^{\mathrm{tot}}$ and $\mathrm{PTE}$ maps, pixels with
$\mathrm{PTE}<0.01$ are highlighted by green borders. The total uncertainty is
defined as

$\Delta\beta_{s}^{\mathrm{tot}}\equiv\sqrt{\left(\Delta\beta_{s}^{\mathrm{stat}}\right)^{2}+\left(\Delta\beta_{s}^{\mathrm{sys}}\right)^{2}},$
(13)

where $\Delta\beta_{s}^{\mathrm{stat}}$ is the statistical uncertainty. We
adopt two different ways to estimate the $\Delta\beta_{s}^{\mathrm{stat}}$.
For pixels with $\mathrm{PTE}<0.01$, the uncertainty was obtained by
bootstrapping. When estimating the bootstrapping uncertainty
($\Delta\beta_{s}^{\mathrm{BS}}$) at each $N_{\rm side}=8$ pixel, we fit for
slopes from 5000 different data resamplings using Equation 11 with only the
diagonal components of the covariance matrix (resampled in the same way), and
convert them to the spectral indexes. Finally, the
$\Delta\beta_{s}^{\mathrm{BS}}$ was estimated as the standard deviation of the
spectral indexes fitted from different samples. As shown in the bottom panel
of Figure 17 and Figure 18, regions with high synchrotron S/N ratio tend to
have extremely low PTE values, which is at least in part because the intrinsic
scatter of the synchrotron signal is not considered in the covariance matrix
used in Equation 11. For other pixels, the statistical uncertainty was
obtained by propagating the $\Delta k$ estimated with the full covariance
matrices as $\Delta\beta_{s}^{\mathrm{COV}}=\Delta
k/[\hat{k}\log(\nu_{\mathrm{40}}^{s}/\nu_{K}^{s})]$. The noise correlation
between pixels is not negligible in the relatively low S/N ratio region. In
our case, neglecting the off-diagonal components in the covariance matrix
tends to give smaller estimations of the uncertainty, which is reflected by
the ratio in Figure 18 with larger PTE values.

Figure 18: The ratio of the statistical uncertainty on $\beta_{s}$ using the
bootstrapping method ($\Delta\beta_{s}^{\mathrm{BS}}$) to that obtained by
propagating the $\Delta\hat{k}$ estimated with the full covariance matrices
($\Delta\beta_{s}^{\mathrm{COV}}$), as a function of the $\mathrm{PTE}$
values. The $x$-axis goes to log-scale at $\mathrm{PTE}<0.01$. The vertical
dashed line ($\mathrm{PTE}=0.01$) and the horizontal line
($\Delta\beta_{s}^{\mathrm{BS}}/\Delta\beta_{s}^{\mathrm{COV}}=1$) are
included for reference. All data points with $\mathrm{PTE}<10^{-3}$ are
shifted to $\mathrm{PTE}=10^{-3}$ and marked with horizontal left arrows to
show that they have smaller $\mathrm{PTE}$ values. For pixels with
$\mathrm{PTE}<0.01$ (orange triangles), we adopt the
$\Delta\beta_{s}^{\mathrm{BS}}$ as the statistical uncertainty. For the rest
pixels (blue circles), the $\Delta\beta_{s}^{\mathrm{COV}}$ is chosen as the
statistical uncertainty.

We used simulations to verify that the $\hat{k}$ in Equation 11 is an unbiased
estimation and that $\Delta k$ is near optimal. In the fiducial simulations,
the synchrotron amplitude template was chosen to be the PySM (Thorne et al.,
2017) s1 model (unrelated to s1 masks in Section 6.1) at WMAP K band, and we
scaled the amplitude template to the CLASS $40\,\mathrm{GHz}$ band following
the spectral index model of s1 at $N_{\rm side}=8$. The simulated synchrotron
amplitude maps were created at $N_{\rm side}=128$ with a $2^{\circ}$ FWHM beam
and were re-observed through the CLASS $40\,\mathrm{GHz}$ mapping pipeline. We
combined the synchrotron templates with 100 different realizations of rWMAP
K-band noise, CLASS ${40\,\mathrm{GHz}}$ noise, and the CMB simulations, all
created as described above. The CMB is only included in the CLASS
${40\,\mathrm{GHz}}$ simulation under the assumption that rWMAP K band data
are synchrotron-dominated. We fit for $\hat{k}$ and $\Delta k$ from each pair
of simulations. The per-pixel difference between the simulation ensemble mean
of $\hat{k}$ and the s1 model, normalized by the ensemble standard deviation
on the mean, follows a standard normal distribution with values mostly within
$\pm 3$ across the 100 $N_{\rm side}=8$ super-pixels in the map, which means
that the estimation of $\hat{k}$ is unbiased. The fractional difference
between the ensemble average of $\Delta k$ to the ensemble standard deviation
of the $\hat{k}$ is largely within $\pm 20\%$, which means that the $\Delta k$
is near optimal. To understand the impact of the CMB, we compared the ensemble
standard deviation caused by only CMB simulations to that by only noise
simulations, and the former is found to be at the $10-25\%$ level of the
latter.

Since the CLASS mapping transfer function mixes the signal from pixels in
different directions and since $\beta_{s}$ is known to be non-uniform on the
sky (Krachmalnicoff et al., 2018; Fuskeland et al., 2021; Weiland et al.,
2022), the method of fitting $\beta_{s}$ from the CLASS maps and rWMAP K band
data induces systematic errors that depend on the morphology of both the
synchrotron amplitude and the spectral index. We perturb both the fiducial
synchrotron amplitude and spectral index template (one at a time) based on the
current constraints on them to estimate the systematic error caused by the
method. For the amplitude template, we perturbed the rWMAP K-band data with
the rWMAP K-band noise simulations, and the perturbed template was then scaled
to CLASS ${40\,\mathrm{GHz}}$ according to the s1 $\beta_{s}$ model and re-
observed before the fitting. For the spectral index template, we perturbed the
s1 model with Gaussian random realizations that have a uniform standard
deviation equal to 0.1 per $N_{\rm side}=8$ pixel over the sky (Planck
Collaboration et al., 2020e). Each template was perturbed 100 times, and we
obtained the systematic uncertainty as the standard deviation of the
difference between the fitted and the input $\beta_{s}$ template. No noise and
CMB were included in the simulations. We found that the systematic errors on
$\beta_{s}$ are negligible for the variations in the amplitude morphology but
not for the variations in the spectral index, and we denote the latter the
systematic uncertainty $\Delta\beta_{s}^{\mathrm{sys}}$ on the $\beta_{s}$
estimation. The distribution of the $\Delta\beta_{s}^{\mathrm{sys}}$ has a
median and 0.16/0.84 percentiles of $0.033^{+0.035}_{-0.015}$.

While the tests described in the previous paragraph measure the _standard
deviation_ in recovered $\beta_{s}$ relative to an input template when either
the amplitude or spectral index are perturbed, they do not test the _mean_
offset between recovered $\beta_{s}$ and the input due to the mapping transfer
function. To estimate this bias, we scaled the synchrotron amplitude template
at WMAP K band to the CLASS $40~{}\mathrm{GHz}$ band following a power law
with the $\hat{\beta}_{s}$ map as shown in Figure 17, where the missing values
were filled by $-3.1$. We then re-observed the new template at CLASS
$40~{}\mathrm{GHz}$ and fit for the $\beta_{s}$ from the re-observed noiseless
simulations using the same covariance matrices as for the main result. The
difference between the fitted $\beta_{s}$ and the input model is consistent
with the $\Delta\beta_{s}^{\mathrm{sys}}$ template. Therefore the bias in
$\beta_{s}$ introduced by the mapping transfer function is captured within the
systematic uncertainty estimate.

The distribution of spectral indices in Figure 17 has a median and 0.16/0.84
percentiles of $\hat{\beta}_{s}=-3.10^{+0.19}_{-0.10}$. This is consistent
with previous pixel-based studies that computed $\beta_{s}$ along the plane
(Kogut et al., 2007; Ruud et al., 2015; Weiland et al., 2022). To understand
how the uncertainty on the thermal dust frequency scaling model impacts the
$\beta_{s}$ estimation, in addition to the uniform $\beta_{d}=1.53$ and
$T_{d}=19.6~{}\mathrm{K}$ templates used for the main results, we also used
the $\beta_{d}$ and $T_{d}$ templates with spatial variation, derived with the
Commander component separation algorithm using the Planck 2015 data (Planck
Collaboration et al., 2016a). The difference of the $\hat{\beta}_{s}$ fitted
with different dust frequency scaling models
($\delta\hat{\beta}_{s}(\mathrm{dust})$) scaled by the
$\Delta\beta_{s}^{\mathrm{tot}}$ is shown in Figure 19. The
$\delta\hat{\beta}_{s}(\mathrm{dust})$ has a median and 0.16/0.84 percentiles
of $-0.004^{+0.01}_{-0.01}$, and we found that
$|\delta\hat{\beta}_{s}(\mathrm{dust})|<0.2\times\Delta\beta_{s}^{\mathrm{tot}}$
is achieved in most directions, meaning that the uncertainty in the dust
frequency scaling model does not contribute much to the uncertainty on
$\hat{\beta}_{s}$ in most regions. The $\delta\hat{\beta}_{s}(\mathrm{dust})$
is significant in regions closer to the Galactic plane, as expected from the
complexity of astrophysical processes contributing to the dust signal within
that region. This shows that a better characterization of the thermal dust
frequency scaling is becoming more and more important for precise $\beta_{s}$
estimation, and vice versa.

Figure 19: The difference of the $\hat{\beta}_{s}$ fitted with different dust
frequency scaling models ($\delta\hat{\beta}_{s}(\mathrm{dust})$) scaled by
$\Delta\beta_{s}^{\mathrm{tot}}$. The $\hat{\beta}_{s}$ fitted using the
$\beta_{d}$ and $T_{d}$ templates with spatial variation is subtracted from
that fitted with uniform $\beta_{d}=1.53$ and $T_{d}=19.6~{}\mathrm{K}$.
Regions with $\Delta\beta_{s}^{\mathrm{tot}}>0.2$ and those masked by the
declination limit (light gray) and point-source mask are excluded.

Finally, the PTE of the hypothesis that all $\hat{\beta}_{s}$ are consistent
with their mean value is below $10^{-6}$ for both the uniform dust model and
the one with spatial variation, meaning that we detected the spatial variation
on the $\hat{\beta}_{s}$ at $N_{\rm side}=8$ resolution, using rWMAP K and
CLASS $40~{}\mathrm{GHz}$ bands data. We have yet to consider the impact from
CLASS ${40\,\mathrm{GHz}}$ 5% calibration uncertainty on the fitted
$\beta_{s}$. Changing the CLASS data by $\pm 5\%$ is equivalent to shifting
the $\beta_{s}$ measurement in all $N_{\rm side}=8$ pixels of Figure 17 by a
common $\pm 0.1$. The variation in the spectral index about its median value
is essentially independent of a 5% change in the normalization of the
${40\,\mathrm{GHz}}$ data.

##  8 Cosmology

The CLASS ${40\,\mathrm{GHz}}$ channel was designed specifically to target
Galactic synchrotron emission with sufficient sensitivity to enable high-
precision cleaning of regions on the sky with relatively low Galactic
foreground contamination. Using external data sets, however, these foreground-
dominated maps have shown sufficient sensitivity to probe CMB power.

###  8.1 CMB $EE$ Spectrum

On large scales ($\ell<30$), the $EE$ angular power spectrum contains key
information for understanding the reionization history of the Universe
(Zaldarriaga, 1997). The difficulties in measuring this signal are underscored
by the fact that the parameter directly tied to the $EE$ amplitude on these
scales, the optical depth to reionization, $\tau$, remains the least well
measured $\Lambda$CDM parameter (Pagano et al., 2020). To demonstrate the
progress made with these new maps collected by a ground based platform, we
compute cross-spectra between CLASS and external data to verify the large
angular scale CMB power has been preserved in the CLASS ${40\,\mathrm{GHz}}$
maps. We use a pseudo-$C_{\ell}$ estimator (Chon et al., 2004) for this
purpose, while recognizing it is non-optimal for the largest-scale bins. A
more optimal estimator, such as xQML (Vanneste et al., 2018), will be used for
future work—especially as the remaining higher-sensitivity CLASS channels are
incorporated in our analysis framework.

The measurement of the CMB $EE$ spectrum using CLASS ${40\,\mathrm{GHz}}$
$\times$ Planck ${100\,\mathrm{GHz}}$ and ${40\,\mathrm{GHz}}$ $\times$ Planck
${143\,\mathrm{GHz}}$ is shown in Figure 20—the Planck PR3 maps (Planck
Collaboration et al., 2020c) were used to avoid the extra computational
complication needed to estimate and correct for the mapping transfer functions
associated with the later PR4 release (Planck Collaboration et al., 2020d).
Evaluation of the PR4 mapping transfer function within our sky region could be
used in future work. The method to estimate these spectra is described below.

Figure 20: _Top._ CMB $EE$ spectrum computed over the range $5<\ell<125$ from
foreground reduced CLASS 40 $\mathrm{GHz}$ $\times$ Planck PR3 100
$\mathrm{GHz}$ or Planck PR3 143 $\mathrm{GHz}$. Uncertainties are estimated
using simulations and are dominated by CLASS noise. As the synchrotron
foreground channel, it is expected that the CLASS 40 $\mathrm{GHz}$ noise
would dominate here — it is the CLASS 90 $\mathrm{GHz}$ channel that is
optimized for sensitivity to the CMB. Horizontal bars indicate the bins. As
described in the text, template subtraction has been used to remove an
estimate of the foreground emission within the s9 mask. The solid line is the
Planck best-fit $\Lambda$CDM spectrum (Planck Collaboration et al., 2020a).
_Bottom._ The uncertainty normalized difference between each spectra and the
binned $\Lambda$CDM expectation is shown.

To target the CMB, synchrotron and dust templates are subtracted from each
channel. The rWMAP K-band (WMAP Ka-band) was used as the synchrotron template
for CLASS (Planck channels). The Planck 353 $\mathrm{GHz}$ channel was used as
a dust template in the Planck 100 and $143\,\mathrm{GHz}$ maps — within the
mask region considered and the level of sensitivity achieved, dust is expected
to be a negligible component of the CLASS $40\,\mathrm{GHz}$ maps. The
foreground-cleaned maps are then defined as

$\hat{m}_{x}=\frac{{m}_{x}-c_{1}{m}_{\mathrm{sync}}-c_{2}{m}_{353}}{1-c_{1}-c_{2}},$
(14)

where $m_{x}$ is the $Q$ or $U$ map of the channel to be cleaned,
$m_{\mathrm{sync}}$ is rWMAP $K$-band (WMAP $Ka$-band) for cleaning
synchrotron from CLASS (Planck), and $m_{353}$ is the Planck 353
$\mathrm{GHz}$ channel used for dust removal. The final template coefficients
are collected in Table 5. The dust coefficients are directly adopted from
Pagano et al. (2020). For the synchrotron coefficients, we assume the very
simple model of $\beta_{s}=-3.1$, and scale the rWMAP $K$-band (WMAP
$Ka$-band) to the CLASS (Planck) channel. A more careful analysis minimizing
residual foregrounds and considering more complicated foreground models will
be needed as the remaining higher-sensitivity CLASS channels are included in
future analysis. For our purpose of demonstrating recovery of CMB power
without attempting to place constraints on particular cosmological parameters,
the current simple approach is sufficient.

Table 5: Polarized foreground removal template coefficients.

Channel |  $c_{1}$ (synchrotron) |  $c_{2}$ (dust)
---|---|---
CLASS 40 | 0.207 | 0
Planck 100 | 0.039 | 0.0186
Planck 143 | 0.0165 | 0.0394

Note. — Template scaling coefficients apply to maps in thermodynamic
temperature. The synchrotron template is WMAP K-band for CLASS 40
$\mathrm{GHz}$ and WMAP Ka-band for Planck 100 $\mathrm{GHz}$ or 143
$\mathrm{GHz}$. The dust template is the Planck PR3 353 $\mathrm{GHz}$ map.

The maps are weighted using an estimate for the inverse variance maps with the
foreground-reduction procedure applied. These weight maps are estimated using
simulations. The CLASS noise was modeled by coadding random pairs of the 200
split noise simulations — repeated use of a single split was allowed, but no
pairing was duplicated. The FFP10 noise simulations were used for the Planck
channels (Planck Collaboration et al., 2016c); each simulation was anti-alias
filtered, rotated to celestial coordinates, and downsampled to
$N_{\mathrm{side}}=128$. The WMAP noise simulations were performed following
the same procedure described in Section 4. Finally, since the reobserved
$K$-band was used, the $K$-band noise simulations were reobserved as well.
Since $Ka$-band was used as the synchrotron cleaning template for both Planck
channels, all data and simulations used for the Planck channels were smoothed
to match the $Ka$-band resolution. The weight map for the cleaned CLASS map
was estimated by applying the cleaning procedure in Equation 14 to the
combined CLASS noise and r$K$-band noise simulations. The pixel covariance was
then directly estimated from an ensemble of 5000 noise simulations (allowing
repeat instances, but no repeat of any particular combination), from which the
weight map was computed following Equation 2. An analogous procedure was used
to estimate the weight in the Planck channels. Finally, the CLASS weights were
multiplied by the apodized s9 mask and the Planck channels were multiplied by
the respective Plik frequency masks combined with the total Planck point-
source mask (Planck Collaboration et al., 2020b).

Uncertainties were estimated using the same suite of simulations described
above, except now including a reobserved CMB realization as part of the CLASS
channel and the same CMB instance (but not reobserved) for each simulation in
the Planck channel. Cross spectra were estimated assuming the map weights
described above. The square root of the diagonal of the resulting estimated
covariance matrix is shown as the uncertainties in the $40\times 100$ and
$40\times 143$ cross spectra in Figure 20. Both cross spectra are found to be
consistent with each other and with the best fit $\Lambda$CDM model (Planck
Collaboration et al., 2020a). Efforts for the inclusion of the more sensitive
higher-frequency CLASS channels are well underway.

###  8.2 Circular Polarization Limits

A possible circular polarization background has been predicted by multiple
theories, e.g. Lorentz violating physics (Alexander et al., 2009; Caloni et
al., 2023) and Faraday conversion (De & Tashiro, 2015), among others. The
unique capacity for CLASS to measure circular polarization allows us to search
for such a signal. In the current work, we improve on the results already
demonstrated (Padilla et al., 2020). The new mapping procedure, L23, and null
test framework, Section 5.1, have been applied to the extended data described
here. The $V$ maps shown in Section 3 appear to be free from any significant
signal above the noise. Similar to the transfer function correction described
in Section 4, the $V$ spectra are corrected to create a non-biased estimate of
the signal power. Note the filtering applied to the $V$ channel removed less
power than the filtering applied to the linear polarization maps, and
therefore the $VV$ spectra require less of a correction. The $VV$ cross-power
spectra, evaluated on the base data split and using the respective $w_{VV}$
pixel weights over the full survey area is shown in Figure 21. The band power
signal covariance is estimated using 10000 noise-only CLASS simulations—where
the mapping transfer function is corrected in each simulation as well. The
binned $VV$ spectrum has $\chi^{2}=10.6$ with 12 degrees of freedom and is
therefore consistent with noise. The resulting $95\%$ quantiles of the band
power distribution were then interpreted as upper limits in Figure 22. After
accounting for binning differences, the new spectrum reduces the upper limit
on large-scale circular polarization by over two orders of magnitude over
results from others (Nagy et al., 2017; Mainini et al., 2013), and
approximately an order of magnitude over our previous CLASS result (Padilla et
al., 2020). The improvement in this work comes from the increase in the data
volume and the optimal map-making that suppresses the correlated noise at
large angular scales (L23).

Figure 21: The CLASS ${40\,\mathrm{GHz}}$ $VV$ spectrum is shown in
$\mathrm{\mu K^{2}}$ thermodynamic units. The spectrum is evaluated over the
full survey area accounting for the mapping and beam window functions. The
result is consistent with no circular polarization signal in our maps.  Figure
22: Upper limits ($95\%$ confidence level) on cosmological circular
polarization spectrum. Former measurements from MIPOL (Mainini et al., 2013)
at $33\,\mathrm{GHz}$ and Spider (Nagy et al., 2017) at $95/150\,\mathrm{GHz}$
are shown in comparison with the CLASS measurement at $40\,\mathrm{GHz}$. The
constraints from Padilla et al. (2020) used only the Era 1 nighttime data, and
the maps were inverse-variance weighted.

##  9 Conclusions

In this paper, we have presented the initial results collected from the 40
$\mathrm{GHz}$ channel of CLASS—the first of four frequency channels. Data
from the survey operating from 31 August 2016 to 19 May 2022 have been reduced
to maps and have been demonstrated to have improved noise performance over
previous measurements at nearby frequencies in the critical $10<\ell<100$
range. The low-$\ell$ performance is set by the current filtering strategy
used in the mapmaking process to reduce large angular scale systematic errors,
and the resultant maps have been shown to pass a broad array of null tests.
The maps and associated data products have been made available to NASA LAMBDA
for public release.

The significant increase in sensitivity over previous measurements surveying
large fractions of the sky near 40 $\mathrm{GHz}$ is striking. Furthermore, we
have found remarkable agreement between reobserved higher signal-to-noise
measurements of synchrotron emission from lower frequencies once scaled to be
comparable to the CLASS 40 $\mathrm{GHz}$ maps; this includes bright regions
along the Galactic plane and lower brightness regions well removed from the
plane. CLASS has replicated and improved upon measurements of the same sky
signal as compared to space missions.

As a synchrotron-dominated channel, we have characterized the power law
behavior of the angular power spectrum and found general agreement with
previous measurements. With the inclusion of external data sets, we have
measured the synchrotron SED off the Galactic plane by comparing the amplitude
of the signal in the harmonic domain. Furthermore, we have obtained a new
$\beta_{s}$ map at $N_{\rm side}=8$ resolution, and measured the variations of
the SED in expanded region along the Galactic plane relative to previous
efforts.

The unique capability of the CLASS array to measure circular polarization
continues to allow for the search for a cosmic circular polarization signal.
We have presented new upper limits on the existence of such a signal,
improving on the CLASS-set state of the art by a significant margin.

We have shown, through cross spectra between foreground-reduced versions of
the CLASS and Planck 100 and 143 $\mathrm{GHz}$ maps, that the $EE$ power of
the CMB has been measured even in the CLASS 40 $\mathrm{GHz}$ foreground
channel. This result foreshadows CMB-focused analyses when the other CLASS
frequency channel data are included — an effort well underway. The low-$\ell$
noise performance of the 40 $\mathrm{GHz}$ survey already matches the
optimistic forecast of upcoming major observatories (e.g., Simons Observatory
Collaboration et al., 2019). Through our self-consistency tests, we have made
significant progress in understanding the nature of the systematic errors that
must be overcome to make progress at the lowest $\ell$. Therefore, the CLASS
project is at the forefront of pushing ground-based polarization measurements
to the lowest $\ell$’s to provide access to the reionization optical depth and
eventually primordial $B$-modes through the reionization signal.

##  10 Acknowledgments

We acknowledge the National Science Foundation Division of Astronomical
Sciences for their support of CLASS under Grant Numbers 0959349, 1429236,
1636634, 1654494, 2034400, and 2109311. We thank Johns Hopkins University
President R. Daniels and the Deans of the Kreiger School of Arts and Sciences
for their steadfast support of CLASS. We further acknowledge the very generous
support of Jim and Heather Murren (JHU A&S ’88), Matthew Polk (JHU A&S Physics
BS ’71), David Nicholson, and Michael Bloomberg (JHU Engineering ’64). The
CLASS project employs detector technology developed in collaboration between
JHU and Goddard Space Flight Center under several NASA grants. Detector
development work at JHU was funded by NASA cooperative agreement
80NSSC19M0005. We acknowledge the use of the Legacy Archive for Microwave
Background Data Analysis (LAMBDA), part of the High Energy Astrophysics
Science Archive Center (HEASARC). HEASARC/LAMBDA is a service of the
Astrophysics Science Division at the NASA Goddard Space Flight Center. CLASS
is located in the Parque Astronómico Atacama in northern Chile under the
auspices of the Agencia Nacional de Investigación y Desarrollo (ANID).

We acknowledge scientific and engineering contributions from Max Abitbol,
Fletcher Boone, David Carcamo, Manwei Chan, Benjamín Edwing Fernández Ríos,
Joey Golec, Dominik Gothe, Francisco Espinoza, Ted Grunberg, Mark Halpern,
Saianeesh Haridas, Kyle Helson, Gene Hilton, Connor Henley, Lindsay Lowry,
Jeffrey John McMahon, Nick Mehrle, Carolina Morales Perez, Ivan L. Padilla,
Gonzalo Palma, Lucas Parker, Sasha Novack, Bastian Pradenas, Isu Ravi, Carl D.
Reintsema, Gary Rhoades, Daniel Swartz, Bingjie Wang, Qinan Wang, Tiffany Wei,
Ziáng Yan and Zhuo Zhang. For essential logistical support, we thank Jill
Hanson, William Deysher, Joseph Zolenas, LaVera Jackson, Miguel Angel Díaz,
María José Amaral, and Chantal Boisvert. We acknowledge productive
collaboration with Dean Carpenter and the JHU Physical Sciences Machine Shop
team.

S.D. is supported by an appointment to the NASA Postdoctoral Program at the
NASA Goddard Space Flight Center, administered by Oak Ridge Associated
Universities under contract with NASA. K.H. was supported by NASA under award
number 80GSFC21M0002. R.R. acknowledges support from ANID BASAL project
FB210003. Z.X. is supported by the Gordon and Betty Moore Foundation through
grant GBMF5215 to the Massachusetts Institute of Technology. R.D. thanks ANID
for grants BASAL CATA FB210003 and FONDEF ID21I10236.

This paper uses data products derived from observations obtained with Planck
(http://www.esa.int/Planck), an ESA science mission with instruments and
contributions directly funded by ESA Member States, NASA, and Canada.

## References

  * Abbott et al. (2022) Abbott, T. M. C., Aguena, M., Allam, S., et al. 2022, Phys. Rev. D, 105, 043512
  * Addamo et al. (2021) Addamo, G., Ade, P. A. R., Baccigalupi, C., et al. 2021, J. Cosmology Astropart. Phys, 2021, 008
  * Aiola et al. (2020) Aiola, S., Calabrese, E., Maurin, L., et al. 2020, J. Cosmology Astropart. Phys, 2020, 047
  * Alam et al. (2021) Alam, S., Aubert, M., Avila, S., et al. 2021, Phys. Rev. D, 103, 083533