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

Measurement of the largest angular scale (ℓ < 30) features of the cosmic microwave background (CMB) polarization is a powerful way to constrain the optical depth to reionization and search for the signature of inflation through the detection of primordial B-modes. We present an analysis of maps covering 73.6% of the sky made from the 40 GHz channel of the Cosmology Large Angular Scale Surveyor (CLASS) from 2016 August to 2022 May. Taking advantage of the measurement stability enabled by front-end polarization modulation and excellent conditions from the Atacama Desert, we show this channel achieves higher sensitivity than the analogous frequencies from satellite measurements in the range 10 < ℓ < 100. Simulations show the CLASS linear (circular) polarization maps have a white noise level of 125(130)μKarcmin . 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 sky regions and measure the diffuse SED in the harmonic domain. We place a new upper limit on a background of circular polarization in the range 5 < ℓ < 125 with the first bin showing D ℓ < 0.023 μKCMB2 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.


INTRODUCTION
Measurements of the cosmic microwave background temperature and polarization angular power spectra (e.g., Hinshaw et al. 2013;Aiola et al. 2020;Planck Collaboration et al. 2020a;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, Λ.The resultant Λ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 Λ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 (Kamionkowski et al. 1997;Seljak & Zaldarriaga 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 ℓ > 30 supported by Planck and WMAP data (BICEP/Keck Collaboration et al. 2021).Other recent B-mode measurements at ℓ > 30 include Kusaka et al. (2018), Polarbear Collaboration et al. (2020), and Spider Collaboration et al. (2022).BB measurements on the largest angular scales (ℓ < 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 ℓ < 30 is an upper limit from Planck (r < 0.274; de Belsunce et al. 2023).
Tensor modes from inflation also source low-ℓ 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 ℓ < 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 (τ 2 ; Zaldarriaga 1997).Of the six standard ΛCDM parameters, the least well constrained is τ .While other parameters are measured with subpercent precision, the reionization optical depth uncertainty is 10% (Pagano et al. 2020).Improvement in the measurement precision of τ will break a critical degeneracy with the amplitude of the initial scalar curvature fluctuations A s .Improved knowledge of τ , 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).
In this paper, we present polarization maps, angular power spectra, and analyses of Galactic and cosmological components based on 40 GHz1 observations from 2016 August to 2022 May.In a companion paper, we describe the 40 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.

Instrument Description
The CLASS telescope array resides in the Atacama Desert of northern Chile and surveys the sky with singlefrequency telescopes at 40 and 90 GHz, as well as a dualfrequency 150/220 GHz telescope.With its observing site at latitude 23 • S, large instantaneous field of view (FOV), and ability to perform constant-elevation scans over 720 • in azimuth, the array can efficiently survey 73.6% 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 foreground minimum depends upon the exact sky region considered, and tends to be near 80-90 GHz when large fractions of the sky are considered (Planck Collaboration et al. 2020e).
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 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 preparation).
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 ∼ 50 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-ofband 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. 2014Appel et al. , 2019;;Dahal et al. 2018Dahal et al. , 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 2021 September.
Pointing is achieved through a three-axis mountboresight 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 GHz telescope observations beginning 2016 August 31, after initial commissioning, and extending to 2022 May 19, when the telescopes were shut down for maintenance and instrument upgrades.We define the period before 2018 February 22 as Era 1 and after 2018 June 22 as Era 2 -various instrument changes and upgrades were performed in the intervening time.
The 40 GHz telescope includes 36 feedhorns (72 TES bolometers) with array average 1.56 • FWHM beam (using a Gaussian beam approximation).The FOV is 20 • wide in azimuth and 15 • wide in elevation at zero boresight (Eimer et al. 2012;Xu et al. 2020).More details can be found in L23.

Observing Strategy
Observations were organized and scheduled on a 24 hr 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 • to +45 • relative to vertical in 15 • 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 ±360 • in azimuth at 45 • elevation and constant boresight.Over the course of the scan, the FOV of the telescope would sweep out a swath roughly 15 • -20 • in elevation, depending on the boresight angle, centered on the 45 • boresight direction.The azimuth scanning direction typically reversed at 180 • azimuth (due south) during the night.In the daytime, the azimuth range was reduced to keep the array boresight center > 20 • from the Sun to avoid overheating and damaging the instrument.No other celestial objects were avoided by the telescope, and bright objects, e.g., the Moon, are removed in data cuts.More details regarding data selection are discussed in L23.For a schedule dedicated to the CMB survey, observations spanned nearly 24 hr, 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).A summary of key properties of the survey and instrument are presented in Table 1.Further details regarding the survey and instrumental configuration changes can be found in L23.

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.The total survey dates included in these results span a period of 4.32 yr; this range excludes downtime for the 90 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 hr 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•years of data.

Data Acquisition and Pre-processing
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 ms time constant.The MCE readout electronics also applied an antialiasing 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 /emphIVbased calibration (Appel et al. 2022).

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 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 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 (Thorne et al.Gaussian beam for visualization purposes.To show both bright synchrotron features within the Galactic plane and the faint regions elsewhere in the map, the color scale is linear within ±5 µK and logarithmic beyond this range.The dipolar atmospheric V signal has been removed from the V map through the mapmaking filters, and no remaining Stokes V signal is detected-the map is consistent with noise.The survey coverage is nonuniform; see L23 for hit maps.In addition to anisotropy in the signal, the noise also changes over the sky. 2017; Petroff et al. 2020b), 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 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.

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 pickup from the mount azimuth motor current by removing the first five harmonics of each full 1440 • (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 2021 September, 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 pickup, electric pickup, or, for circular polarization, Zeeman-splitting of atmospheric oxygen in the Earth's magnetic field).For linear (circular) polarization, this filter removed 15 (10) azimuthal harmonics.These three filters were estimated piecewise over timescales ranging from 2-6 hr 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 2017 August 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 pickup 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 & Naess 2018, L23).With this noise covariance, linear and circular polarization maps were made via a preconditioned conjugate gradient solver with an adapted version of minkasi2 (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 ≥ 97.5% of the signal in the filtered demodulated data.
An initial calibration of the maps in temperature units was based on measurements of the Moon and Jupiter (Appel et al. 2019;Xu et al. 2020;Dahal et al. 2022) and is described in Li et al. (2023).Each detector's calibration was found to be stable to the subpercent level over the course of the survey, a consequence of the stable atmospheric emission from the CLASS site at 40 GHz (Appel et al. 2022).A final calibration correction was performed by comparing the bright regions in the CLASS maps directly to those same regions in WMAP Ka and Q band maps.This procedure is described in Section 6.3.The stability of the calibration is further verified by a series of null tests described in Section 5.1.
The final linear and circular polarization maps are shown in Figure 1.These are represented in HEALPix format with a resolution parameter of N side = 128 (Górski et al. 2005) (∼ 27.5 ′ pixels over-sampling the telescope beam).We have made maps and associated data products available to NASA LAMBDA for public release. 3or 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 re-observed maps.For this work, we have used re-observed WMAP K -, Ka-, and Q-band (Bennett et al. 2013) and the Planck 353 GHz map (Planck Collaboration et al. 2020c).When using a re-observed map, we prepend the name with "r," e.g., rWMAP K -band or simply rK.

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 where C ℓ ′ ,in refers to the EE, BB, and V V angular power spectra of the unfiltered sky.The transfer function F ℓℓ ′ relates these to the angular power spectrum as observed by CLASS, C ℓ,out , accounting for correlations introduced between multipole moments and between spectra.Simulations showed that when this transfer function is used with a pseudo-C ℓ estimator (Chon et al. 2004), the recovered estimate of the input spectra is unbiased for ℓ > 4. The effect of the mapping transfer function is to bias signal power to be low on the largest angular scales.For linear (circular) polarization, the transfer function was found to leave approximately 85% (93%), 67% (85%), 34% (47%), and 6% (17%) of the signal in the angular power spectrum at ℓ = 30, 20, 10, 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 ℓℓ ′ 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 nonuniform CLASS weight is used, for example see 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 the Appendix.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 FOV with different elliptical orientations and by including the effects of boresight rotations and scan cross-linking in the survey map.The resulting beam has an FWHM of 1.559 • ± 0.006 • , a solid angle of 799.4 ± 0.5 µsr, and an eccentricity of 0.07 ± 0.01, where the uncertainties are formal statistical uncertainties 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%, and 9% of its intrinsic value at ℓ = 20, 30, 50, and 100, respectively.For more details on the beam and beam window function estimate, see the Appendix.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.

POWER SPECTRUM ESTIMATION
Angular power spectra for synchrotron signal analysis and all null tests were estimated by using the pseudo-C ℓ method as implemented in PolSpice (Chon et al. 2004) 4 .This method produces nonbiased 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 the scalar weight (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.In this work, we have used crossspectra exclusively.As a result, 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 evenand odd-numbered years define the split.For the Planck 30 GHz channel, no splits were considered; only crossspectra 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 GHz data.Thus, K ×C indicates the cross-spectrum between the full K-band maps and the full CLASS maps, and Ka × Ka denotes the crossspectrum 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 where F ℓℓ ′ 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 ℓℓ ′ 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 NaMaster5 (Alonso et al. 2019) C2 apodization tool with a 3 • scale.The map weight used for spectrum estimation was then the product of the weight from Equation 2 with the apodized mask.Simulations assuming the angular power spectrum followed a power-law model or a nominal fixed synchrotron sky model plus instrument noise were used to confirm the estimators were unbiased.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 antialiasing cosine-like filter (Benabed et al. 2009) where ℓ 1 = 128, and ℓ 2 = 3 × ℓ 1 .The maps were then downgraded to N 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 ℓ = [5, 125) with fixed width ∆ℓ = 10.Binned spectra use the typical ℓ(ℓ + 1) weight within each bin.Note, that we adopted this binning to help visualize the results of the null test (see Section 5.1) and to help mitigate mode coupling.To be consistent, we adopted these bins for all angular power spectrum results.In all cases, simulations were used to determine the full EE/EB/BB bin-bin covariance matrix.Off-diagonal elements are dominated by the nearest neighboring bins; these elements are ≲ 0.1 for the full region observed by CLASS and can be slightly larger but still ≲ 0.15 when considering the largest foreground masks.
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)-note these simulations include noise and residual systematic errors.For WMAP simulations, large-scale noise realizations consistent with the published low-ℓ pixel-pixel noise covariance matrix were combined with high-ℓ noise realizations consistent with the polarization hits (Larson et al. 2011).The WMAP and Planck simulation maps were smoothed with an antialiasing 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 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 ≲ ℓ ≲ 100 range.

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.
The filtering outlined in Section 3.3 was designed to mitigate systematic effects, some of which were initially discovered through early null test investigations.While the current filtering was found to be sufficient to remove systematic errors to the level probed by these consistency checks, the detailed filtering parameters were not pre-cisely tuned to minimize the signal loss.More optimal filtering strategies will be considered in future work.

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.Each span was placed into one of the two base sets while 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.
The base splits were further subdivided into subsplits, A and B, according to other criteria, e.g., distinct detector sets, defining four maps, m iρ with i ∈ {1, 2} and ρ ∈ {A, B}.The A/B maps were solved simultaneously similar to the procedure described in L23.This approach ensured that the noise modeling of data from each split was maintained regardless of the detector or temporal division and that the data were combined in the same way as in the total map.The differences (null maps), ∆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.
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-yr Q-band smoothed with a 1.5 • 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 ∆m 1 and ∆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 ∆m b i was removed from each respective null map to produce a bias-corrected null map: For the simulated null maps, the ∆m b i templates are the average of the simulated null maps ⟨∆m i ⟩ in which the noise is greatly diminished.For the null maps made from data, the ∆m b i template is made from a signal-only simulation created by smoothing the synchrotron-dominated WMAP 9 yr K -band map (Bennett et al. 2013) with a 1.5 • Gaussian beam and scaling it by the approximate ratio (0.17) in synchrotron intensity between the WMAP K band and the CLASS 40 GHz band.
Finally, the bias-corrected null spectra from CLASS were compared to the corresponding ensemble of simulations via the χ 2 test where C 1×2 ℓ is the binned CLASS null cross-spectra between ∆m c 1 and ∆m c 2 , ⟨C sim ℓ ⟩ is the average of the corresponding simulated binned null power spectra,7 Cov ℓℓ ′ 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 χ 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 selfconsistency standards in the bright Galactic center as they do in the synchrotron-diffuse regions measured by the angular power spectra.

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/Bsplit by choosing subsets of half the available detectors; the full collection of detector splits is 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.
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, 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  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 rotationboresight 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 R.A. or decl., 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 • 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 • scan of the mount into the first and second 360 • scans and would be sensitive to the impact of the azimuth servo in a different way.For example, the first (second) 360 • 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 ∼ 40 • 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 hr 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 6 yr survey.

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 (and summarized χ 2 PTEs in Figure 6) from the detector-based splits are collected in Figure 7, and the temporal splits are in Figure 8.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 χ 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 fullcovariance matrix as in Equation 6.The χ 2 PTEs for all spectra are collected in Figure 6.The split names indicated in the Figures follow the definitions outlined in the text.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 lowest PTE value, 0.01, occurs in the midnight split of the V E cross-spectrum.Careful study of these null maps and spectrum do not indicate a clear problem.The low PTE value is driven by the prevalence of slightly negative band powers in the null spectrum, seen in Figure 8.It is also possible there is some very small signal or unknown systematic effect remaining in this particular split; it could also just be noise.The number of PTE values near 1 indicate that the noise might be slightly overestimated in some cases.All uncertainties, however, retain this noise model.
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.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.Additionally, all masks exclude data outside of the CLASS survey boundary at declinations δ = −76 • and δ = 30 • and a bright-source mask.
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).Each mask was defined following a five-step process: 1.Each Commander map, e.g., the synchrotron polarized intensity map, was smoothed with a Gaussian  0.78 (0.67) Figure 7. 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 χ 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 of the data points are set to unity.In the null test, the χ 2 statistics are based on the full-covariance matrix (Equation 6); the reduced-χ 2 and PTE (in parentheses) values are displayed in the lower-right corner of each panel.The corresponding PTE values are collected in Figure 6.beam with FWHM selected from 0.5 • -12 • in 10 logarithmically spaced intervals.
2. A threshold, t X , was chosen at each smoothing scale.Pixels exceeding the threshold were masked to form a series of 10 masks for each component map.
3. The resulting binary mask was smoothed again with a 5 • Gaussian beam to remove sharp features and then reset to a binary state with a threshold at 0.25. 5. Any pixel beyond the CLASS survey region or within the bright-source mask, described below, was also masked.
The thresholds, t X , used in step two, are defined as where b X,ℓ is the smoothing window function for each respective Gaussian beam used in step one.This construction corresponds to the pixel variance for a Gaussian random field with a power-law spectrum.The only free parameter is the overall scaling, which was fine-tuned to visually reject the Galactic features at each scale.
The Galactic masks for s3, s7, and s9 boundaries are shown in purple, teal, and yellow, respectively, in Figure 9.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 9, 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 approximate l ∼ (180 • , 300 • ) range in Galactic latitude.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 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 • , 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 • diameter exclusion region centered at each of the selected points (twice the CLASS beamwidth), and custom regions are tailored to surround a few extended sources.Finally, on account of its exceptionally high intensity, a 5.4 • diameter region is excluded around Tau A. The source mask is shown in red in Figure 9. Nearly all of the masked sources reside within the Galactic plane and would already be covered by the foreground mask.The individual sources within the plane are still useful when considering the bright synchrotron near the plane, but not associated with particular compact regions.All masks were included in the public release hosted on the NASA LAMBDA site.

Color correction
To collect as much signal from the sky as possible, incoherent detectors used for surveys observing near the CMB foreground minimum are typically designed to have fractional bandwidths in the range of 0.2-0.3.Such a broadband detector integrates the incoming signal over the SED of all sources along a line of sight-different SED shapes can produce different detector responses.One method to account for this effect, called color correction, is described in Appendix E of Bennett et al. (2013).To facilitate the comparison of CLASS maps with other experiments for astrophysical sources with nonCMB 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 where ∆T CMB /∆T A = 1.038 is the antenna-tothermodynamic unit conversion factor at 38 GHz (Dahal et al. 2022)-the approximate band center of the CLASS 40 GHz channel; T A (ν i ) are the source brightness temperatures evaluated at fixed frequencies, and ω 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 ν i are somewhat arbitrary and were chosen to even out the weights w i .Table 2 summarizes these quantities for the CLASS 40 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 10.

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 re-observed WMAP Ka and Q bands.Thus, the CLASS calibration depends upon the low calibration uncertainty from the nearby frequency channels of WMAP and removes any dependence on the initial Moon and Jupiter calibration.
The CLASS 40 GHz, and rWMAP Ka/Q-band maps were smoothed to a common resolution of 2 • and downgraded to N side = 32 (1.8 • pixels).The bright source mask was applied prior to smoothing to diminish the impact from bright-sources, and only N side = 32 pixels derived from more than 5 of the possible 16 unmasked N side = 128 subpixels 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 11.
The uncertainty per N 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 side = 16 by sampling according to the full noise covariance matrix. 8White noise realizations were created at N side = 128 by sampling according to the per-pixel Q/U covariance matrices.The simulations were converted to spherical harmonic coefficients a ℓm , and composite spectra were made by concatenating the full-covariance a ℓm 's at ℓ ≤ 48 with the white noise simulations at ℓ > 48.The combined spectra were then converted back to N side = 128, smoothed to 2 • FWHM, and downgraded to N side = 32.For CLASS, details of the noise simulations can be found in L23.These were masked, smoothed, and downgraded in the same way as the CLASS data.Gaussian random CMB simulations were generated at N side = 128 according to the Planck best-fit ΛCDM CMB power spectra 9 and then masked and smoothed to 2 • FWHM, and downgraded to N side = 32.To create the effect of reobservation, we applied the CLASS harmonic domain mapping trans-fer function (Equation 1) to all simulations except the CLASS noise simulations (already filtered).
The calibration factor was obtained via a two-step process.First, the spectral index β s was estimated within the region selected by s1-d1 as: where d X , σ 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 x(ν, β s ) means that the quantity x has been converted into antenna-temperature units and color corrected to frequency ν assuming a power-law spectrum with index β s .In this step, we used the rWMAP Ka and rWMAP Q bands data color corrected to the same frequency ν = 38 GHz using parameters from Table 20 of Bennett et al. (2013).The best-fit value is βs = −3.23.
In the second step, we began by computing the effective frequency for CLASS 40 GHz (ν eff ) where the color correction factor is unity for the best-fit βs .The ν eff can be solved using Equation 8, by setting 3 i=1 ω i (ν i /ν eff ) βs = 1.The effective frequency is ν eff = 37.56 GHz.
We then converted CLASS 40 GHz, rWMAP Ka, and Q data into antenna-temperature units and then color corrected the rWMAP Ka-and Q-band data to ν eff .Provided the signal in the selected sky region scales like the assumed power law, all three maps should have the same signal up to an overall calibration factor, η.To find this factor, the model y = ηx was fit with y as either rWMAP Ka-or Q-band and x as CLASS 40 GHz, taking into consideration the uncertainties on both axes.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 reduced χ 2 of the fit, and the PTE values are shown in Figure 11.In all of the fits, a strong linear relationship is present with r ≥ 0.94, and the reduced χ 2 is close to unity.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 GHz are η Ka = 1.18 ± 0.01 and η Q = 1.18 ± 0.02, where the smaller uncertainty from the rWMAP Ka band fit is a result of its higher synchrotron signal-to-noise ratio (S/N).The fact that the fitted slopes using the CLASS data match each other well and that from the rWMAP fit 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 β s = −3.23.
We express the absolute calibration factor as η = η 0 ± σ . The calibration value η 0 = 1.18 is the inverse-variance weighted average of η Ka and η Q .The statistical uncertainty σ 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 GHz data were used in both, though the correlation would be minimal as the CLASS 40 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 is about 10% larger than the propagated uncertainty, which rounds up to the same value of 0.01.The systematic uncertainty σ sys η = 0.05 results from the 0.5 GHz uncertainty on the band center of CLASS 40 GHz (Dahal et al. 2022), which is the difference in the calibration value when converting rWMAP Ka and Q bands to antenna temperature at ν eff ± 0.5 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 η 0 = 1.18.
The final calibration factor is η = 1.18 ± 0.01(stat) ± 0.05(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, ∼ 3%, atmospheric opacity that is not included in the baseline temperature calibration.The uncertainty of the initial Moon and Jupiter based calibration is ∼ 2%.However the remaining discrepancy is still ∼ 14 ± 5% and cannot be fully accounted for by the known calibration uncertainties and the atmospheric opacity correction.
To check that the d1 mask provides robust protection against calibration bias from dust contamination, we performed the same procedure using the smaller d0 mask and subtracted a dust template formed by scaling rPlanck 353 GHz-38 GHz following a modified blackbody spectrum with β d = 1.53 and T dust = 19.6K.The derived calibration differed by only 0.001 indicating dust negligibly impacts our calibration.
As a check on the survey-based calibration, we confirmed that the measured polarized flux density of Tau A, using aperture photometry, is consistent with previous measurements; see Figure 12.We measured the Tau A flux density by integrating the polarization P = Q 2 + U 2 in antenna temperature within a disk of radius 1.85 • centered at R.A. = 83.633• and decl.= 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 14).This baseline correction increased the measured flux density by ∼ 8%.
An additional ∼ 4% flux density was added to account for the beam solid angle outside the disk radius (see the Appendix).The integrated Tau A antenna temperature (in units of K-sr) was divided by ΓΩ = c 2 /(2kν 2 eff ) to obtain the Tau A flux density in units of janskys (Jy).We subtracted a 0.2 Jy noise bias from the CLASS Tau A flux measurement.The bias was estimated by applying the same algorithm to compute the flux density across the 200 noise only CLASS map simulations.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 α = −1 as opposed to the known Tau A index α = −0.35(Weiland et al. 2011).To correct for this difference, based on simulations testing the impact of this effect, the map-based CLASS Tau A polarization flux density measurement was divided by 1.04.We applied the same aperture photometry method to re-observed WMAP and Planck maps and found Tau A flux density values consistent with the published values.Simulations showed that the CLASS mapping transfer function (Equation 1) diminishes flux density measure- ments of unresolved sources (like Tau A) by ≤ 1%.For simplicity, a correction of this bias was not applied to this Tau A measurement.Gaussian beam for visualization purposes.This causes the signal features to be smoothed to 2.25 • in the CLASS maps.The noise between the maps can still be visually compared.The synchrotron foreground is much brighter at the lower-frequency ( 23GHz) 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.

Comparison with Previous Measurements
The CLASS 40 GHz survey complements the collection of available surveys targeting the synchrotron foreground component of the polarized millimeter sky-with the WMAP Q-band (41 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 β ≈ −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 GHz, Bennett et al. 2013) data with their high-S/N polarized synchrotron emission to qualitatively validate the CLASS maps in this Section.
Figure 13 shows full-sky comparisons between CLASS polarization maps and re-observed WMAP K -and Qband 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 • beam-no additional smoothing was performed to those maps.For visualization purposes, the CLASS maps have been smoothed with a 1.5 • beam.Smoothing the CLASS map has the effect making the noise comparable between the three maps, but the signal in the CLASS maps is smoothed to a slightly lower resolution, 2.25 • -this effect is too small to be visible in the figures.There are two main qualitative conclusions to draw from Figure 13.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 < ℓ < 100 in Figure 2. The second is that the synchrotron emission qualitatively matches the rWMAP K -band signal on the largest angular scales.
Figures 14 and 15 provide a zoomed-in comparison between the CLASS polarization maps to the rWMAP K -band data -no further image processing is performed.Figure 14 gives a view of the Galactic plane in the range −25 • < b < 25 • .Overall there is excellent visual agreement between the two data sets.Some prominent Galactic features, such as Tau A at (ℓ, b) ≈ (185 • , −5 • ) 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 (ℓ, b) ≈ (260 • , 0 • ) and the Centaurus A radio galaxy at (ℓ, b) ≈ (310 • , 20 • ), match in morphological detail.Figure 15 is a more detailed look at a subset of diffuse polarized features in the maps.The amplitude of the features is approximately 5 − 10 µK.At this brightness contrast, the apparent discrepancies between the maps on the smallest scales are 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 • across the maps.

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.Polarized free-free emission is expected to be negligible and is ignored in this work.Spinning dust, a hypothesized explanation for anomalous microwave emission (AME), is also expected to be negligible (Draine & Hensley 2016) in line with current upper limits (Herman et al. 2023).Polarized emission from magnetic dust is also expected to be very faint (Hoang & Lazarian 2016).Thus, we do not consider AME and magnetic dust in this work.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 ν ≲ 60 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 GHz, we examine behavior of synchrotron emission in both harmonic and pixel space.
In units of antenna-temperature, synchrotron emission is often approximated as a power-law over a limited range of frequencies, with T ∝ ν βs , where β s is the spectral index.Previous studies using data from, e.g., WMAP, Planck S-PASS, and QUIJOTE have firmly established spatial variations of β 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 β 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.

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 ℓ estimator (PolSpice, Chon et al. 2004) to estimate the cross power spectrum over the ℓ-range (5, 125).While less optimal than a quadratic estimator (e.g., xQML, Vanneste et al. 2018) for ℓ < 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 that cause excess mode mixing in the PolSpice spectrum estimation.Normally this pseudo-C ℓ 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 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 16.For comparison, 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 β s = −3.1 and color corrected to a common frequency of 38 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 binned and directly 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 10,000 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 diagonal of the inferred binned covariance matrix.The measured CLASS spectra are in reasonable agreement A comparison of the CLASS 40 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.This causes the signal features to be smoothed to 2.25 • in the CLASS maps.The noise between the maps can still be visually compared.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.
with the WMAP results 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 GHz (8 GHz bandwidth, Page et al. 2003), the WMAP Ka-band is 28-37 GHz (6.9 GHz bandwidth, Page et al. 2003), while the CLASS 40 GHz channel uses a somewhat larger 32.8-43.6GHz (10.9 GHz bandwidth) range (Dahal et al. 2022).See also Figure 10.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 16.These spectra indicate that diffuse polarized synchrotron emission still dominates the polarized CMB signal outside the s9 mask.The new CLASS measurement is seen to improve upon the measured spectra from the nearest frequency space-based measurement.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.Visually, there is a suggestion that the angular power spectrum is flattening, especially for ℓ ≥ 50 in the WMAP Q-band spectra.While this is consistent with the WMAP Ka and CLASS measurements, we do not attempt to constrain a flattened component using the low-S/N measurements in this region.As a check on this trend, we extended the WMAP cross-spectra to higher multipoles, but the signal is too faint for a clear indication of a flattened region.Since point sources can contribute a flat spectral component (e.g., Wright et al. 2009;Krachmalnicoff et al. 2018), we also tried masking the remaining Planck identified 30 and 44 GHz polarized point sources.This had a negligible impact on the spectra.Future measurements with higher sensitivity will help improve our understanding of this portion of the spectra.

Synchrotron SED
Within the frequency range from a few gigahertz up to 100 GHz, the SED of the diffuse synchrotron radiation has been measured to follow a power law with index β s ≈ −3 with small differences depending upon the sky region and angular scale measured, see for example Krachmalnicoff et al. (2018)  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 (Vidal et al. 2015).To this end, we considered data from WMAP K -, Ka-, and Q-bands, Planck 30 GHz, and the CLASS 40 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 where {X, Y } indicated the data sets being crossed.Provided the index α is consistent between different channels, the amplitude, A 40 , measures the synchrotron SED at the effective frequency ν eff = √ ν X ν Y .The frequencies ν 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 SEDs for each of the E and B components are assumed to follow a simple model where we have adopted the pivot frequency ν 0 = 30 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 10 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 the noise sampled as described in Section 4, the CMB realizations following the Planck best-fit ΛCDM model (Planck Collaboration et al. 2020a), and the Gaussian realizations of a fiducial synchrotron model.10Power 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 that both required transfer function correction, the spectrum was corrected with F −1 ℓℓ ′ , and for terms with only one CLASS noise map, the correction was done with F −1/2 ℓℓ ′ .The resulting A 40 amplitudes and uncertainties for the s3, s7, and s9 masks are collected in Table 3 and shown in Figure 17.The cross-spectra including CLASS data are indicated with the open circled, relatively larger, and bolder colored points in the plot.The power-law index from Equation 10 stays roughly constant for all cross-spectra and masks.The index for the s9 mask, for example, is α EE = −2.42± 0.14 and α BB = −2.84± 0.25 where we report the mean value and median error taken over the 14 cross power spectra.
The second method was to fit jointly Equation 11 to the full set of cross power spectra.In this case, the power law from Equation 10 is still fit to the angular power spectrum of each frequency, but now the amplitudes were constrained to follow Equation 11.The resulting best-fit power laws for the s3, s7, and s9 masks are shown in Figure 17 and show good agreement with the single-frequency measurements.For this fit, the EE and BB spectra were fit separately.Furthermore, the 168 × 168 covariance matrix (12 bins × 14 cross-spectra) was assumed to be block diagonal with only the internal cross-spectra 12 × 12 sub-block considered to be nonzero.The best-fit power-law index, β s , is shown in the Figure legend.The ratio of the best-fit BB/EE amplitudes is 0.42±0.02,0.39±0.02,and 0.33±0.02for the s3, s7, and s9 masks.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 β 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 GHz.

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 of 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 r of the order of 0.01.
In this Section, we present a polarized synchrotron spectral index (β s ) map made by analyzing the CLASS 40 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 lack of significant Faraday rotation and depolarization effects.Using the re-observed WMAP data mitigates the bias on β s that would result from comparing the raw (unbiased) WMAP data to the filtered CLASS data.We investigate potential residual bias from the reobservation processing below.
The bright-source mask, described in Section 6.1, is applied at N side = 128 to avoid bias on β s from the different physical conditions within those regions.To mitigate the impact of polarized thermal dust emission in the CLASS 40 GHz maps, we scaled and subtracted the rPlanck 353 GHz data assuming a modified blackbody spectrum with a uniform β d = 1.53 and T d = 19.6K across the sky (Planck Collaboration et al. 2020e).Later in this Section, we test the impact of this dust model assumption.Following the process used in Section 6.3, the masked and dust subtracted maps were smoothed and downgraded to N side = 32 (1.8 • pixels).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 11. 353 .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)  The maps of β s were computed in the Galactic coordinate system at resolution N side = 8 (7.3 • pixels) following the linear correlation method (Fuskeland et al. 2014(Fuskeland et al. , 2021;;Weiland et al. 2022).Within each N side = 8 "superpixel", an x-y scatter plot was formed using the unmasked where ν s 40 and ν s K can be found in Table 4.As d(k) depends on the color correction, and thus β s , the values were computed iteratively, with initial value β s = −3.1.Convergence to within 0.01% was achieved in five iterations.
Figure 18 shows the β s map, its total uncertainty ∆β tot s map, and a map of the PTE associated with the linear fit.The distribution of spectral indices has median βs = −3.09+0.23  −0.11 (16th/84th percentiles)-consistent with previous studies that computed β s along the plane (Kogut et al. 2007;Ruud et al. 2015;Weiland et al. 2022).
The total uncertainty is defined as where ∆β stat s (∆β sys s ) is the statistical (systematic) uncertainty.Regions with ∆β tot s > 0.2 and those masked by the decl.limit (light gray) and the bright-source mask are excluded.In the ∆β tot s and PTE maps, pixels with PTE < 0.01 are highlighted by green borders.As shown in the bottom panels of Figure 18, regions with high synchrotron S/N tend to have these 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 12.In these pixels, it could be possible to resolve finer-scale spectral index variations.For this investigation, however, we capture this scatter in the statistical uncertainty using bootstrapping.We determined the bootstrapping uncertainty, ∆β BS s , by fitting k for 5000 different data resamplings using only the diagonal components of the covariance matrix (resampled in the same way).The ∆β BS s was estimated as the standard deviation of the spectral indexes inferred from different samples.
For all other pixels, the statistical uncertainty was obtained by propagating the ∆k estimated with the fullcovariance matrices as ∆β COV s = ∆k/[ k log(ν s 40 /ν s K )].A comparison of the two statistical uncertainty estimates for all mapped pixels is shown in Figure 19.The noise correlation between subpixels is most significant in the relatively low-S/N region.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 19 with larger PTE values.
We estimated the systematic uncertainty, ∆β sys s , by considering the impact of variations in the synchrotron morphological and spectral properties.As a basis for comparison, we adopted the PySM (Thorne et al. 2017) synchrotron amplitude model at WMAP K -band.First, to test for a dependence on the morphology of the synchrotron signal, we generated an ensemble of perturbations of this template by adding to it 100 K -band noise simulations.Each instance was then scaled to CLASS 40 GHz according to the PySM s1 spectral index variation model (unrelated to s1 masks in Section 6.1)11 .The ensembles of K -band and scaled maps were then re-observed with the CLASS pipeline and combined to estimate β s at N side = 8 following the procedure described above.The differences between the recovered β s map and the PySM model were negligible, and therefore the uncertainty on the synchrotron morphology does not contribute to the systematic uncertainty.Second, to check for a dependence of the assumed spectral index variation model, we fixed the K -band template and perturbed the s1 model.Gaussian random values with standard deviation equal to 0.1 were added to each pixel in the s1 template.The single K -band model was then scaled to CLASS 40 GHz using 100 different perturbed s1 models.As before, the K -band map and scaled maps were re-observed with the CLASS pipeline and combined to estimate β s .In this case, the variations in the recovered spectral index had a median deviation 0.033 +0.035 −0.015 (16th/84th percentiles), which we identify as the systematic uncertainty on the β s estimation, ∆β sys s .As a check on a possible mean offset between the recovered β s and the input model due to the mapping transfer function, we scaled the K -band template to the CLASS 40 GHz band using the estimated values shown in Figure 18 with missing values filled by −3.1.After both maps were re-observed, we fit for β s using the same covariance matrices as used for the main result.The difference between the fitted β s and the input model is consistent with the ∆β sys s estimate.Therefore, the bias in β s introduced by the mapping transfer function is captured within this systematic uncertainty estimate.
The impact from a multiplicative bias can be bounded by the CLASS calibration uncertainty.A shift in calibration by ±5% is equivalent to shifting the β s measurement in all N side = 8 pixels of Figure 18 by a common ±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 GHz data.The impact of this uncertainty is not included in the ∆β tot s estimate, but it should be kept in mind when interpreting the map of β s .
To verify k in Equation 12 is unbiased and ∆k is near optimal, we performed an additional series of simulations.We again used the PySM model at WMAP K -band and a scaled map at CLASS 40 GHz using the s1 index model.An ensemble of K -band maps was made by adding 100 K -band noise simulations after which each instance was re-observed with the CLASS pipeline.Similarly, 100 CLASS noise simulations and CMB simulations are added to the synchrotron model at the CLASS 40 GHz band and re-observed with the CLASS pipeline.The pairs of simulated maps were then used to generate β s maps following the procedure used for the main result.The per-pixel difference between the ensemble mean of β s and the s1 model, normalized by the ensemble standard deviation on the mean, follows a standard normal distribution with values mostly within ±3 across the map, indicating k is unbiased.The fractional difference between the ensemble average of ∆k to the ensemble standard deviation of the k is largely within ±20%, indicating ∆k is near optimal.We note, using an analogous set of simulations with only CMB or only noise being added to the CLASS 40 GHz band synchrotron model, that we found the CMB contributed 10%-25% of the error relative to the CLASS noise alone.
Next, we tested how these results depend upon the assumed dust model.Rather than using a constant spectral index for the polarized dust emission, we used the β 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 βs fitted with different dust frequency scaling models (δ βs (dust)) scaled by the ∆β tot s is shown in Figure 20.The δ βs (dust) has a median of −0.005 +0.01 −0.02 (16th/84th percentiles), and we found that |δ βs (dust)| < 0.2 × ∆β tot s is achieved in most directions, meaning that the uncertainty in the dust frequency scaling model does not contribute much to the uncertainty on βs in most regions.The δ βs (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.Better characterization of the thermal dust frequency scaling is required for more precise β s estimation in these regions.

COSMOLOGY
The CLASS 40 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.

CMB EE Spectrum
On large scales (ℓ < 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, τ , remains the least well-measured Λ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 GHz maps.We use a pseudo-C ℓ estimator (Chon et al. 2004) for this purpose, while recognizing it is nonoptimal for the largestscale bins.A more optimal estimator, such as xQML (Vanneste et al. 2018), will be used for future workespecially as the remaining higher-sensitivity CLASS channels are incorporated in our analysis framework.
The measurement of the CMB EE spectrum using CLASS 40 GHz × Planck 100 GHz and 40 GHz × Planck 143 GHz is shown in Figure 21-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.
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 GHz channel was used as a dust template in the Planck 100 and 143 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 GHz maps.The foreground-cleaned maps are then defined as where m x is the Q or U map of the channel to be cleaned, m sync is rWMAP K-band (WMAP Ka-band) for cleaning synchrotron from CLASS (Planck), and m 353 is the Planck 353 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 β 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 when 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.
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 antialias filtered, rotated to celestial coordinates, and downsampled to N side = 128.The WMAP noise simulations were performed following the same procedure described in Section 4. Finally, since the re-observed K-band was used, the K-band noise simulations were re-observed 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 15 to the combined CLASS noise and rK-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 re-observed) 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 × 100 and 40 × 143 cross-spectra in Figure 21.Both cross-spectra are found to be consistent with each other and with the best-fit ΛCDM model (Planck Collaboration et al. 2020a).Efforts for the inclusion of the more sensitive higher-frequency CLASS channels are well underway.

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 nonbiased estimate of the signal power.Note that the filtering applied to the V channel removed less power than the filtering applied to the linear polarization maps, and therefore the V V spectra require less of a correction.The V V cross power spectra, evaluated on the base data split and using the respective w V V pixel weights over the full survey area is shown in Figure 22.The band power signal covariance is estimated using 10,000 noise-only CLASS simulationswhere the mapping transfer function is corrected in each simulation as well.The binned V V spectrum has χ 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 23, e.g., the first bin places the upper limit D ℓ < 0.023 µK 2 CMB .After accounting for binning differences, the new spectrum reduces the upper limit on large-scale circular polarization by over 2 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 mapmaking that suppresses the correlated noise at large angular scales (L23).

CONCLUSIONS
In this paper, we have presented the initial results collected from the 40 GHz channel of CLASS-the first of four frequency channels.Data from the survey operating from 2016 August 31 to 2022 May 19 have been reduced to maps and have been demonstrated to have improved noise performance over previous measurements at nearby frequencies in the critical 10 < ℓ < 100 range.The low-ℓ 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 systematic effects motivating the filtering described in L23 are responsible for the decreased sensitivity as one moves to the lowest ℓ values.As explained in L23, the cause of some of these systematic effects has already been fixed.The remaining systematic effects are the object of ongoing instrument improvements from which we anticipate further improvements on the largest scales.
The increase in sensitivity over previous measurements surveying large fractions of the sky near 40 GHz is significant-map sensitivity is up to a factor of 2.1 (1.9) times higher than WMAP Q (Planck 44).Furthermore, we have found agreement between re-observed higher S/N measurements of synchrotron emission from lower frequencies once scaled to be comparable to the CLASS 40 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.In the most restrictive synchrotron mask, s9, we find the diffuse synchrotron is characterized by a power law with β s = −2.95± 0.09 for EE and β s = −3.28 ± 0.24 for BB.The ratio of the amplitude BB/EE for this region is 0.33 ± 0.02.Furthermore, we have obtained a new β s map at N side = 8 resolution, and measured the variations of the SED of the total polarized signal in an 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 an order of magnitude-D ℓ < 0.023 µK 2 CMB at 95% confidence for the lowest ℓ bin.
We have shown, through cross-spectra between foreground-reduced versions of the CLASS and Planck 100 and 143 GHz maps, that the EE power of the CMB has been measured even in the CLASS 40 GHz foreground channel.This result foreshadows CMB-focused analyses when the other CLASS frequency channel data are included-an effort well underway.The low-ℓ noise performance of the 40 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 ℓ.Therefore, the CLASS project is at the forefront of pushing ground-based polarization measurements to the lowest ℓ's to provide access to the reionization optical depth and eventually primordial B-modes through the reionization signal.

Figure 1 .
Figure 1.The Q/U /V Stokes parameter maps for the CLASS 40 GHz channel are shown in Galactic coordinates using a Mollweide projection and (l, b) = (0, 0) placed at the center of each plot.Gray regions indicate portions of the sky not observed by CLASS or included too few observations to be well constrained.The maps have been smoothed with a 1.5• Gaussian beam for visualization purposes.To show both bright synchrotron features within the Galactic plane and the faint regions elsewhere in the map, the color scale is linear within ±5 µK and logarithmic beyond this range.The dipolar atmospheric V signal has been removed from the V map through the mapmaking filters, and no remaining Stokes V signal is detected-the map is consistent with noise.The survey coverage is nonuniform; see L23 for hit maps.In addition to anisotropy in the signal, the noise also changes over the sky.

Figure 2 .
Figure2.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+systematic error 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 ℓ is due to the CLASS 1.56 • beam.CLASS has the lowest noise in the range 10 ≲ ℓ ≲ 100.

Figure 3 .Figure 4 .
Figure 3.The base split, m1 − m2, null maps for the Q/U data are shown.For clear visualization, the maps have been smoothed with a 3 • 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 calibration stability for the base split-the null maps are dominated by noise.

Figure 5 .
Figure5.Detector-splits definition of the CLASS 40 GHz telescope.The detector pairs coupled to the same feedhorn are shown as semicircles with the +45 • /−45 • 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.

Figure 8 .
Figure8.Same as Figure7, but for the temporal splits.

4.
The binary masks, at the native Commander resolution of N side = 2048, were downgraded to N side = 128, masking any lower-resolution pixel if it contained any flagged N side = 2048 pixels.

Figure 9 .
Figure 9.The CLASS 40 GHz analysis masks s3 , s7 , and s9 are shown.The masks share the CLASS decl.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).

Figure 10 .a
Figure10.The era-weighted effective bandpass model of CLASS 40 GHz is shown.For comparison, the detector average of the WMAP Ka-band and WMAP Q-band measurements(Jarosik et al. 2003) is also included.The bandpass model for CLASS 40 GHz era1 (era2) is shown by solid (dashed) gray line.

Figure 11 .
Figure 11.Correlation of Stokes Q (blue) and (orange) between different bands.All data are converted to antenna temperature at ν eff = 37.56 GHz.Left: rWMAP Ka to CLASS 40 GHz.Middle: rWMAP Q to CLASS 40 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 the best-fit slope, the Pearson correlation coefficient r, the reduced χ 2 , and the PTE.The number of degrees of freedom of the fitting is 417.The slightly low PTE values result from the intrinsic scattering of the synchrotron signal within the selected region-bootstrapped uncertainties are ∼ 10% larger than the uncertainties shown here, which rounds up to the same values.The inset within the right panel shows the region selected by s1-d1 mask (black) at N side = 32, where the CLASS decl.limit is in gray.

Figure 12 .
Figure 12.The measured Tau A polarization flux density of 21.4 ± 1.2 Jy at an effective frequency of 38 GHz is consistent with previous measurements.The uncertainty on the CLASS Tau A 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 density(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σ contour of the model's prediction including spectral and time evolution uncertainty.

Figure 13 .
Figure 13.A comparison of the CLASS 40 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.This causes the signal features to be smoothed to 2.25 • in the CLASS maps.The noise between the maps can still be visually compared.The synchrotron foreground is much brighter at the lower-frequency (23 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 Figure1.

Figure 15 .
Figure 15.A comparison of the CLASS (bottom row) and rWMAP K -band (top row) maps in an area of relatively diffuse emission centered on (ℓ, b) = (22 • , −50 • ).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 µK at 40 GHz) features extending 40 • across the map.

Figure 16 .
Figure 16.Top panels: the angular cross power spectra of the diffuse synchrotron component of the CLASS 40 GHz data are shown in units of µK 2 antenna temperature outside of 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 years are also shown.Each spectrum has been color corrected to 38 GHz assuming a fixed synchrotron spectral index of βs = −3.1.Error bars are the binned noise of each measurement and do not include sample variance.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 GHz also assuming βs = −3.1.The Planck best-fit CMB 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.

Figure 17 .
Figure 17.Measurement of the synchrotron SED for the EE (BB) modes is 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 Table3.The curve, and legend values, are from the joint fit of all spectra together assuming the amplitudes follow the power-law SED in Equation11.
to compute the color correction factor for Planck.When computing the color correction factors, βs = −3.1 is assumed for synchrotron, and β d = 1.53,T d = 19.6K for thermal dust.
N side = 32 Stokes Q and U subpixel data from rWMAP K -band (d x ) and CLASS 40 GHz (d y ).The data were fit using a line with slope k and zero intercept.The estimation of the slope in each N side = 8 superpixel was obtained by total least-squares fitting as k = argmin k d(k) T Σ −1 d(k) , (12) where d(k) ≡ d y − kd x , and Σ ≡ k 2 N x + N y is the covariance matrix.The N x and N y are the Q/U covariance matrices obtained from rWMAP K -band noise and CMB + CLASS 40 GHz noise simulations, as described in Section 6.3.The uncertainty on the slope ∆k was obtained by taking half of the difference between the 16th and 84th percentiles of the likelihood function L ∝ exp − 1 2 d(k) T Σ −1 d(k) .We converted the fitted slope k to βs

PTEFigure 18 .
Figure 18.Top image: the βs map at N side = 8 using rWMAP K -band and CLASS 40 GHz data.The vertical white line in the color bar marks the median βs = −3.09.Middle image: the total uncertainty on the spectral index ∆β tot s .Bottom image: the PTE values.Regions with ∆β tot s > 0.2 and those masked by the decl.limit (light gray) and point-source mask are excluded.In the middle and bottom panels, pixels with PTE < 0.01 are highlighted by green borders.

Figure 19 .
Figure 19.The ratio of the statistical uncertainty on βs using the bootstrapping method (∆β BS s ) to that obtained by propagating the ∆ k estimated with the full-covariance matrices (∆β COV s ), as a function of the PTE values.The xaxis goes to log-scale at PTE < 0.01.The vertical dashed line (PTE = 0.01) and the horizontal line (∆β BS s /∆β COV s = 1) are included for reference.All data points with PTE < 10 −3 are shifted to PTE = 10 −3 and marked with horizontal left arrows to show that they have smaller PTE values.For pixels with PTE < 0.01 (orange triangles), we adopt the ∆β BS s as the statistical uncertainty.For the remaining pixels (blue circles), the ∆β COV s

Figure 20 .
Figure 20.The difference of the βs fitted with different dust frequency scaling models (δ βs(dust)) scaled by ∆β tot s .The βs fitted using the β d and T d templates with spatial variation is subtracted from that fitted with uniform β d = 1.53 and T d = 19.6K. Regions with ∆β tot s > 0.2 and those masked by the declination limit (light gray) and point-source mask are excluded.

χ
Figure 21.Top: CMB EE spectrum computed over the range 5 < ℓ < 125 from foreground-reduced CLASS 40 GHz × Planck PR3 100 GHz or Planck PR3 143 GHz.The spectrum is plotted in D ℓ = ℓ(ℓ + 1)C ℓ /2π.Uncertainties are estimated using simulations and are dominated by CLASS noise.As the synchrotron foreground channel, it is expected that the CLASS 40 GHz noise would dominate here-it is the CLASS 90 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 ΛCDM spectrum (Planck Collaboration et al. 2020a).Bottom: The uncertainty-normalized difference between each spectra and the binned ΛCDM expectation is shown.

Figure 22 .Figure 23 .
Figure22.The CLASS 40 GHz V V spectrum is shown in µ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 maps.

Table 1 .
CLASS 40 GHz instrument and survey summary.

Table 1 .
1 Figure 6.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.

Table 3 .
Synchrotron power spectra amplitude values, A40 values for cross-spectra synchrotron SED.EE (10 2 × µK 2 ) BB (10 2 × µK 2 ) The cross-spectrum channel is indicated using an abbreviation of the channels involved.The abbreviations K, Ka, Q are used for the WMAP bands, the Planck 30 GHz channel is denoted 30, and the CLASS 40 GHz channel is denoted C.

Table 4 .
Frequencies, color correction, and thermodynamic-toantenna unit conversion (∆TA/∆TCMB) factors.Note-From left to right: frequency for synchrotron at WMAP K -band ν s K , synchrotron at CLASS 40 GHz ν s 40 and thermal dust at Planck 353 GHz ν d

Table 5 .
Polarized foreground removal template coefficients.Note-Template scaling coefficients apply to maps in thermodynamic temperature.The synchrotron template is WMAP K -band for CLASS 40 GHz and WMAP Ka-band for Planck 100 GHz or 143 GHz.The dust template is the Planck PR3 353 GHz map.