CLASS Data Pipeline and Maps for 40 GHz Observations through 2022

The Cosmology Large Angular Scale Surveyor (CLASS) is a telescope array that observes the cosmic microwave background over 75% of the sky from the Atacama Desert, Chile, at frequency bands centered near 40, 90, 150, and 220 GHz. This paper describes the CLASS data pipeline and maps for 40 GHz observations conducted from 2016 August to 2022 May. We demonstrate how well the CLASS survey strategy, with rapid (∼10 Hz) front-end modulation, recovers the large-scale Galactic polarization signal from the ground: the mapping transfer function recovers ∼67% (85%) of EE and BB (VV) power at ℓ = 20 and ∼35% (47%) at ℓ = 10. We present linear and circular polarization maps over 75% of the sky. Simulations based on the data imply the maps have a white noise level of 110μKarcmin and correlated noise component rising at low-ℓ as ℓ −2.4. The transfer-function-corrected low-ℓ component is comparable to the white noise at the angular knee frequencies of ℓ ≈ 18 (linear polarization) and ℓ ≈ 12 (circular polarization). Finally, we present simulations of the level at which expected sources of systematic error bias the measurements, finding subpercent bias for the Λ cold dark matter EE power spectra. Bias from E-to-B leakage due to the data reduction pipeline and polarization angle uncertainty approaches the expected level for an r = 0.01 BB power spectrum. Improvements to the instrument calibration and the data pipeline will decrease this bias.


INTRODUCTION
The quest for a precise understanding of cosmology has propelled the development of cosmic microwave background (CMB) observations.Satellite missions like COBE, (Boggess et al. 1992), WMAP (Bennett et al. 2003(Bennett et al. , 2013;;Hinshaw et al. 2013), andPlanck (Tauber et al. 2004;Planck Collaboration et al. 2020a) have made all-sky measurements of the CMB anisotropy in temperature and polarization, which are cornerstones supporting the standard model of cosmology.Related ground-based and balloon-borne efforts provide first-look results, cross-checks, and extended capabilities (e.g., higher resolution and/or higher sensitivity) to the satellites.Suborbital efforts also develop new technologies (e.g., cryogenic detectors and high-throughput optics) to space-readiness levels and train the next generation of experimental cosmologists.These experiments typically observe patches of the sky ranging from degree scale to ∼ 40% of the sky (e.g., Aiola et al. 2020).Over the past decade, suborbital observations have offered complementary views of the CMB and tightened the constraints on cosmological parameters through improved measurement of the damping tail of the temperature power spectrum (Choi et al. 2020;Dutcher et al. 2021), the linear polarization spectra at and below ∼ 5 • angular scales (Polarbear Collaboration et al. 2020;BICEP/Keck Collaboration et al. 2022;Spider Collaboration et al. 2022), the gravitational lensing potential (Polarbear Collaboration et al. 2014;BICEP2 Collaboration et al. 2016a;Bianchini et al. 2020;Madhavacheril et al. 2023), deep surveys of galaxy clusters (Bleem et al. 2020;Hilton et al. 2021), and characterization of the Galactic foregrounds (Ruud et al. 2015;Harper et al. 2022;Rubiño-Martín et al. 2023).
However, ground-based CMB polarimetry has been largely unexplored on the largest angular scales (ℓ ≲ 30, Ruud et al. 2015;BICEP2 Collaboration et al. 2016b;Kusaka et al. 2018) due to fluctuations in atmospheric emission and other sources of systematic error arising from the interaction of the telescope with its environment.This has become an impediment to the percent-level characterization of the reionization history of the universe (Zaldarriaga 1997;Pagano et al. 2020;Watts et al. 2020) and to the search for tensor perturbations on the largest angular scales (Guth 1981;Kamionkowski et al. 1997;Seljak & Zaldarriaga 1997;Tristram et al. 2021).While the search for tensor perturbations has progressed at ℓ ≳ 30 led by BICEP/Keck Collaboration et al. ( 2021), the largest-scale B modes would provide the distinctive "reionization peak" feature and would be most unambiguously separable from the late-time lensing effect (Zaldarriaga & Seljak 1998).The largest angular scales also provide access to beyond-the-standard-model physics (e.g., Muir et al. 2018;Hogan 2019;Hogan et al. 2023;Shi et al. 2023) and the physics of the interstellar medium (e.g., Caldwell et al. 2017).It is the goal of the Cosmology Large Angular Scale Surveyor (CLASS) project to develop the technology and techniques needed to explore the large-scale CMB polarization from the ground.
CLASS is a telescope array located in the Atacama Desert in Chile (Essinger-Hileman et al. 2014) that observes at frequencies around 40, 90, 150, and 220 GHz and surveys 75% of the sky every day.Access to the largest angular scales is enabled through rapid front-end modulation with a variable-delay polarization modulator (VPM; Chuss et al. 2012b;Harrington et al. 2018Harrington et al. , 2021)), which also suppresses instrumental polarization.Compared to other modulation technologies, such as the half-wave plate (Kusaka et al. 2014;Takakura et al. 2017), the employment of a VPM enables CLASS's unique sensitivity to circular polarization (Petroff et al. 2020b;Padilla et al. 2020).CLASS also employs boresight rotations, an optical design that prioritizes signal fidelity over maximizing throughput (Eimer et al. 2012;Xu et al. 2020;Datta et al. 2023), and a fully enclosed comoving ground shield, to map the largest angular scales.CLASS measurements complement large-scale data from future satellite (LiteBIRD Collaboration et al. 2023) and balloon-borne (Lazear et al. 2014;Benton 2020;Addamo et al. 2021) telescopes as well as other ground-based strategies (Lee et al. 2020;Addamo et al. 2021;Rubiño-Martín et al. 2023).Major upcoming international-scale groundbased surveys target scales ℓ ≳ 30 (Simons Observatory Collaboration et al. 2019;CMB-S4 Collaboration et al. 2022).
In this paper, we describe the data reduction pipeline and polarization maps of the CLASS 40 GHz data taken through 2022.Angular power spectra and other map-based results are presented in a companion paper (Eimer et al. 2023, hereafter E23).The paper is organized as follows.Section 2 overviews the design of the 40 GHz telescope and the survey.Section 3 explains the main data processing steps.Mapmaking and the survey maps are presented in Section 4. The impact of several systematic issues is reviewed in Section 5.The Appendix provides a description of the pointing model.

OVERVIEW 2.1. The CLASS 40 GHz Telescope
The CLASS telescope array is sited on Cerro Toco at 5140 m in the Atacama Desert of northern Chile (latitude −22.96 • , longitude −67.79 • ), a location long recognized for its microwave-transparent atmosphere due to the combination of high elevation and low precipitable water vapor (PWV; Cortés et al. 2020).During the observations presented in this paper, PWV quartiles were [0.63, 1.10, 1.98] mm. 1iven its proximity to the equator, the site also provides access to most of the sky, which is essential for large-angularscale measurement.
Figure 1 shows schematics for the 40 GHz telescope, both in its single-telescope configuration (2016-2018; left panel) and when it was paired with the 90 GHz telescope (2018present; right panel).The telescope is shown on its mount structure, which includes three axes of rotation: azimuth, elevation, and boresight.The boresight axis is aligned with the Lower-side cage panels are not shown at left to reveal the telescope.Light enters the enclosure by the forebaffle extension, first encountering the VPM, then the primary (obscured) and secondary mirrors, and finally the receiver.The enclosed multireflection design limits spurious signals from stray light.In 2018, the enclosure and forebaffle extension were extended to accommodate the 90 GHz telescope on the same mount.
center of the telescope's field of view and has a full range of motion from −45 • to 45 • with respect to a nominal central position.With the boresight axis so defined, the azimuth and elevation coordinates give the direction of the boresight axis.Since only one linear polarization state is modulated by the VPM at a time, execution of boresight rotation is needed for a full sampling of the linear polarization signal.The elevation axis allows for 26 • − 86 • rotation, but the VPM drive system restricts polarized observations to 40 • − 55 • .The azimuth has a 720 • range centered on the geographic south.Atop the mount, the cryogenic receiver is secured to a baseplate.Supporting instrumentation, including the helium compressor, gas handling system, and drive electronics sit on a platform that moves with the telescope in azimuth, simplifying the routing of cables and hoses for the large azimuth scans.An aluminum cage structure rises above the receiver and supports the telescope mirrors and the VPM.Aluminum honeycomb panels are bolted to the cage to enclose the telescope, blocking radiation paths from the ground.For the majority of observations analyzed here, the aluminum panels were coated on the inside by microwave absorbers (Eccosorb HR-10).Light enters the cage enclosure through an extension cone, which is rolled at the top to decrease diffraction from sources away from the telescope boresight.We refer to the whole upper part of the cage enclosure above the VPM as a forebaffle as it serves to reject incoming stray light through reflection or absorption.This consists of two parts: the forebaffle roof and the forebaffle extension as labeled in Figure 1.
To accommodate the 90 GHz telescope on the same mount in 2018, the cage structure and forebaffle were expanded as shown in the right side of Figure 1.Other notable changes were implemented at this time such that we designate the time before 2018-02-22 as Era 1 and after 2018-06-22 as Era 2. See Section 2.3 for further discussion.Era 1 has a total of 540 days (1.48 yr).Era 2, which was interrupted by the pandemic, has 1038 days (2.84 yr).
After light passes the forebaffle, it first encounters the VPM.Positioning the VPM as the first optical element modulates the incoming polarization such that it can later be separated from instrumental polarization, which is not modulated (e.g., Miller et al. 2016).After the polarization state is encoded by the modulator, the signal is guided by two ambient-temperature mirrors into a cryogenic receiver.The receiver uses a horizontally mounted Bluefors2 pulse-tube-backed helium dilution refrigerator that continuously cools the receiver optics and focal plane (Iuliano et al. 2018).It employs a combined strategy of multilayer foam, reflective metal mesh, and plastic absorptive filters to block infrared radiation.Plastic lenses produce a telecentric image of the sky on a feedhorncoupled detector array at the focal plane.The telescope's average beam FWHM is 1.54 • (E23), and its field of view spans 20 • in the azimuth direction and 15 • in elevation for zero boresight rotation.The absolute polarization angle calibration will be discussed in Section 5.1.1.
In the focal plane, smooth-walled feedhorns (Zeng et al. 2010) couple incoming radiation to 36 detector chips.Onchip ortho-mode transducers cleanly separate the ±45 • linear polarization states, the power levels of which are detected by separate transition-edge-sensor (TES) bolometers (Chuss et al. 2012a).Therefore, the telescope has 72 feedhorncoupled TES bolometers (an orthogonally polarized pair for each horn/chip).Bandpass filters on the detector chip define a frequency band centered on 38 GHz with a 12.3 GHz bandwidth (Dahal et al. 2022).A cryogenic superconducting quantum interference device (SQUID) based timedivision multiplexing architecture and ambient-temperature Multi-Channel Electronics (MCE) read out the detectors (Reintsema et al. 2003;Battistelli et al. 2008).Details regarding the array construction and performance can be found in Rostem et al. (2012), Appel et al. (2014Appel et al. ( , 2019)), and Dahal et al. (2022).The dilution refrigerator cools the entire focal plane to ∼ 40 mK, which serves as the base temperature for the bolometers that are voltage-biased at their 150 mK transition temperature.At these low temperatures, the primary source of bolometric noise is from the stochasticity of the incoming radiation load (1.2 pW median; Appel et al. 2019;Dahal et al. 2022).For Era 1, the nominal detector array sensitivity3 was 32 µK √ s.Changes to the instrument described in Dahal et al. (2022) increased the number of working detectors, improved optical efficiency, and reduced optical loading and resulted in a decrease to 30.5 µK √ s in Era 2.

Observations
This paper covers 40 GHz observations from their beginning on 2016-08-31 through 2022-05-19 (2089 days; 5.72 yr).Nominally, the CMB survey observations were conducted 24 hr per day with no restriction on the time of year.In practice, scheduled calibrations, maintenance, and instrument upgrades as well as unscheduled weather and other events interrupted the survey.Figure 2 shows the survey observation efficiency over time with 100% indicating nominal 24-hr-per-day operation.Two major stoppages totaling 502 days are shown in pink fill.These were due to the addition of the 90 GHz telescope on the same mount in 2018 and to a VPM grid failure in 2020 that was not repairable until COVID travel restrictions were eased a year later.Because these did not have to do with the regular operation of the 40 GHz telescope, we excluded the associated time when estimating the total possible data volume and observing efficiency.With this exclusion, the total time under consideration is 1587 days (4.32 yr), which is divided between Era 1 (1.48 yr) and Era 2 (2.84 yr) as described in Section 2.1.In addition to the two major stoppages shown in pink, a number of other longer periods (not demarcated in Figure 2) were excluded.This could be for sustained bad weather, such as the two roughly month-long periods in the austral fall and winter of 2017, or for systematic instrument malfunction such as after the 90 GHz deployment through early 2019 when radiofrequency interference (RFI) compromised the 40 GHz data.After these exclusions, 64.6% of the data remained (Table 1, "Timeline").Data quality cuts reduced the usable data volume further and will be discussed in Section 3.
During CMB survey observations, the telescope scanned in azimuth with the elevation and boresight angles fixed.With few exceptions, the elevation was held at 45 • .Each day the telescope boresight rotation was set at an angle between −45 • and 45 • in 15 • increments; the seven boresight angles were visited once a week during observations.The azimuth scan covered the full 720 • range, centered on the geographic south (180 • ).Therefore, the telescope traveled from −180 • to 540 • (a forward sweep) and then back from 540 • to −180 • (a backward sweep).This simple back-and-forth circle scan was repeated throughout the observations.Because the telescope is located at latitude −23 • , during each sweep the boresight traced out a large circle on the celestial sphere centered on decl.−23 • with a 45 • radius.Accounting for the footprint of the field of view about the boresight pointing, the CLASS data extend from decl.30 • in the north to −76 • in the south.As the Earth rotates, the circle scans swept through the full 24 hr of R.A. to cover 75% of the sky. Figure 3 shows the CLASS survey region in celestial coordinates overlaid on the Planck synchrotron temperature map (Planck Collaboration et al. 2020c).The footprints of other CMB surveys are also shown.
One exception to the simple circle scan described above occurred during daytime observations when the Sun passes over the scan path.The telescope maintained a separation greater than 20 • between the boresight and the Sun position.Therefore, the scans were truncated to less than 360 • azimuth as the Sun rose or set through the scan path.

Changes to the Instrument and Observations
Several adjustments in instrument configuration or observing strategy have taken place since the beginning of the survey, many of which were due to instrument updates, replacement of failed components, and remedies for systematics guided by analysis.In Section 2.1, we already noted the enlargement of the cage and forebaffle extension to accommodate the 90 GHz telescope-alterations that demarcate parts Cage Camera Interference -On 2017-08-09, webcams were installed inside the cage enclosure to monitor the telescope.These cameras were later found to introduce RFI around the CLASS modulation frequencies and were turned off for CMB observation from 2018-06-22 onwards.The effects in demodulated data are at frequencies below the scanning frequency and have little impact on the result; see Section 4.1.4.
Forebaffle Roof Blackening and Geometry Change -To improve the rejection of stray light, microwave absorbers (Eccosorb HR-10) were attached to the inside top and walls of the forebaffle roof first on 2016-10-25.When the forebaffle roof was replaced in Era 2, it was not only extended to accommodate the 90 GHz telescope; its geometry (e.g., angle of the roof) was also changed to improve the rejection of stray light.Eccosorb HR-10 was also applied to the inside of the new forebaffle roof on 2019-02-25.(The data taken while the new Era-2 forebaffle roof was unblackened were rejected due to RFI; see below, unrelated to the blackening.)The time during which the forebaffle roof was blackened is hatched in Figure 2.
Forebaffle Extension Blackening -To prevent fractionally (∼ 10 −3 ) polarizing reflections seen in detectors toward the edge of the field of view during moon observations (Figure 19 of Xu et al. 2020), the forebaffle extension interior was first blackened with a microwave absorber (Eccosorb HR-10) on 2017-07-20.The blackening was retained through the end of Era 1 and attached to the new baffle at the beginning of Era 2. After suffering damage, the blackening was removed on 2019-01-16 and not replaced until late in Era 2. The time of the survey when the forebaffle extension was blackened is hatched in Figure 2. See also Section 4.1.3.
Variable Speed Azimuth Scan Test -From 2017-04-01 to 2017-05-18, a variable scan speed strategy (dec scan) was explored to even the integration time across the decl.range.This was abandoned for constant velocity scans (az scan) when detector-warming vibrations were induced at certain velocities explored by the dec scans.
Focal Plane Fix -The initial deployment of the 40 GHz receiver had 64 of 72 optical detectors working, 56 of which were in pairs.The inoperable channels were due to multiplexer failures, likely from static discharge on address lines.The receiver updates for Era 2 recovered all optical detectors.
VPM Grid 2 -The Era-1 VPM wire grid had brown oxidation on its wire grid and imperfections toward its lower edge.Suspecting these may be responsible for larger-than-expected grid emissions, a new grid was installed for Era 2. However, no significant change in performance was observed.
Infrared Filter Changes, RFI, and Thin Grille Filter -Infrared filtering changes made between Era 1 and Era 2 increased the telescope's optical efficiency such that the noise-equivalent temperature (for both white and correlated components) dropped by 20%.However, additional RFI, either due to a new Era 2 component (e.g., new VPM drive electronics) or increased susceptibility of the receiver, required the introduction of a warm thin grille filter (TG-filter) at the receiver window on 2019-01-12.This canceled the efficiency gains from the increased in-band transmission of the new infrared filters (Dahal et al. 2022, Cleary et al. in prep.).
Azimuth Scan Speed Increase -The az-scan speed was increased from 1 to 2 deg s −1 on 2019-03-04.This was found to decrease the low-frequency noise in the demodulated data at 40 GHz (Cleary et al. in prep.).This indicates that the correlated noise component may be due to the residual ground signals after filtering (Section 4.1) and/or the temporal correlations in the turbulent atmosphere emission (Morris et al. 2022) leaking into polarization; in the latter case, the faster signal modulation would permit better separation of the correlated noise from the sky signal.
VPM Grid 3 -The second VPM grid failed on 2020-01-08, likely due to heating of the grid-securing epoxy during exposure to direct sunlight.Delayed by pandemic travel restrictions, a third VPM grid was installed, and observations resumed on 2021-02-11.
Closeout Removal -A thin (17.8 µm) plastic environmental seal, the "closeout" in Figure 1, was used where the light enters at the base of the forebaffle extension (diameter ∼1.3 m).
We found that the closeout produced polarization systematics when deformed by the wind (Section 4.1.2).Since 2021-09-27, the closeout has been removed for most CMB observations and has only been temporarily reinstalled during bad weather.

DATA PROCESSING AND DEMODULATION
The CLASS data reduction pipeline is designed to ingest the raw data and output well-characterized maps and spectra.The use of the VPM for signal modulation naturally divides the pipeline into two parts.First, raw detector time-ordered data (TOD) are vetted, calibrated, demodulated, and downsampled into linear and circular polarization TOD.The subsequent mapmaking pipeline identifies and removes polarization systematics and solves for sky maps given the demodulated TOD. Figure 4 provides an overview of this pipeline.In this section, we describe the procedures up to and including the demodulation and defer the discussion of the rest of the pipeline to Section 4.

Data Acquisition and Selection
The detector data and the encoder data from the telescope mount and the VPM are measured synchronously with a rate of 201 samples per second and grouped into clock-aligned 10 minute packages.These data, in combination with various asynchronously collected housekeeping data, are collated and saved on disk as data packages (DataPkg), which are the smallest units for data storage and transmission.(See Petroff et al. 2020a, for details of the data acquisition pipeline.) A total of 224, 781 DataPkgs (5 TB, in compressed volume) were collected for 40 GHz during the 4.32 yr of observations considered here.These are represented by the gray area in the timeline of Figure 2. As discussed in Section 2.2, several periods of sustained poor weather and instrument malfunction were discarded, leaving 64.6% of the total data (Table 1, "Timeline"; unmarked in Figure 2).After this initial exclusion, an additional downselection was performed based on the metadata of each DataPkg.Data acquired during unfavorable conditions (instrument maintenance, short-term poor weather, etc.) were dropped.Furthermore, only the Data-Pkgs collected when the mount elevation was above 40 • were kept for processing.As a final step, DataPkgs were discarded if the cloud cover as monitored by the site webcams (Li et al. 2023) was consistently too high to avoid highly variable data triggered by enhanced optical loading or polarization from clouds (Takakura et al. 2019).These data cuts left 102, 003 (45.4%)DataPkgs (Table 1, second row); their distribution is shown in yellow in Figure 2.
The selected contiguous DataPkgs were concatenated to form spans for data processing.A span typically has ∼ 130 DataPkgs (22 hr) of data, interrupted every day at around noon when the boresight angle was incremented by 15 • and detector gain calibrations were performed.However, spans may also be shorter if observations were interrupted (e.g., for planet observations).While the previous section described data cuts that removed entire DataPkgs, additional analysis was performed on the resulting spans to flag data samples that are of low quality or susceptible to systematic issues.The second half of Table 1 enumerates these flags along with the fraction of data retained at each step.After the TOD selections, there remained 28% of the total data volume, corresponding to 86.77 detector•years of data for mapmaking.These data are represented in green in the timeline of Figure 2. A description of each TOD selection step is given below.

TOD Selection
Survey Interruption -Over the course of a span, any incidental interruptions to the nominal CMB survey, e.g., for targeted source observations, were flagged.Similarly, interruptions to key instrument system performance were flagged.Monitored systems included the VPM, the full cryostat health, the telescope mount, and the detector readout system.Environmental factors were also monitored; in particular, data taken when the PWV exceeded 5 mm were flagged.This threshold was selected as a pragmatic limit by which strong atmospheric effects were avoided with only a modest impact on overall observing efficiency.NOTE-The quoted retained fraction at each step depends on the order in which selections are applied.
Transient Detector Cuts -The detector TOD were analyzed to flag periods of anomalous behavior.Data from detectors with constant output were flagged.Data for which the SQUID flux-locked closed-loop detector readout became unstable were flagged.To recover these unlocked detectors and reestablish the SQUID tuning state optimized at the start of each span, the first stage SQUID feedback for all detectors was occasionally relocked.The data at the time of relock were flagged.Additionally, the SQUID readout can experience sudden jumps that manifest as discontinuities in the detector data; such jumps were flagged.During any window of 100 samples (0.5 s), if the total array of detectors experienced more than 10 jumps, all detectors were flagged for that window.Finally, the 100 samples (0.5 s) surrounding the azimuth drive direction reversals were flagged.Further information regarding the operation of the SQUID multiplexing amplifiers can be found in Reintsema et al. (2003) and Doriese et al. (2016).VSS Amplitude Cuts -The strongest signal at the modulation frequency of 10 Hz is the VSS, which serves as a good indicator of a detector's optical response.An estimator of the strength of the VSS was computed every 10 minutes for each detector across a span.A detector's data over the entire span were flagged if the detector's median VSS strength estimator was less than five times the standard deviation of the detector's VSS strength estimator over the span (i.e., the VSS strength estimator had a signal-to-noise, S/N, ratio below five).
Conditioning Cuts -The final mapmaking operates on 10sweep segments of data for noise modeling (Section 4).The 10-sweep segments of data of each detector were dropped for analysis if the retained data fraction was below 52% to improve the stability of filtering and noise model estimation.

Pointing
During data acquisition, the pointing model was used to convert the commanded position to the encoder positions used by the servo system.As the DataPkgs were read-in to form the span, the telescope pointing was reconstructed from the recorded encoder data by inverting an updated pointing model through simple iterative successive substitution.In this way, the pointing model used for the analysis was generally not the pointing model used to acquire data: it was revised based on additional calibrations.It is shown in Xu et al. (2020) that the long-term stability of the 40 GHz telescope pointing model is better than 1.4 arcminute.For a detailed description of the pointing model, see Appendix A. The computation of detector pointings from the pointing model used an adapted version of qpoint.4

Calibration
The spans constructed from the DataPkgs were encoded in raw data acquisition units.All span data were first calibrated to physical units of power detected at each bolometer by applying TES responsivity gain factors extracted from the current-voltage characteristic (I-V load curve) of the detector.The uncertainty of each CLASS TES responsivity factor is approximately 0.5% (Appel et al. 2022).Previous CMB experiments have also derived their TES detector calibration from I-V data (e.g., Dünner et al. 2013;Kusaka et al. 2018).The power detected at the bolometers was converted to sky thermodynamic temperature taking into account the optical efficiency of each detector, the optical depth of the atmosphere at the detector elevation pointing, and the CMB power-to-temperature relationship evaluated at the detector frequency bandpass.The relative and absolute optical efficiencies of each detector were obtained from source observations such as the Moon or Jupiter (Appel et al. 2019;Xu et al. 2020;Dahal et al. 2022).

Detector Time Constants and Butterworth Filter
The sky signal was low-pass filtered twice before being recorded on the disk.The detectors have time constants that depend on their electro-thermal properties at any given time.The detector electro-thermal behavior varies in response to changes in optical loading and bias current.The detector filtering can be modeled as a single-pole low-pass filter with 3 − 4 ms time constant ( f 3dB = 40 − 50 Hz) applied to the data (Appel et al. 2019;Dahal et al. 2022).Furthermore, the MCE sampled the detectors at around 23 kHz (Dahal 2020) and applied a fourth-order low-pass Butterworth filter with 50 Hz cutoff frequency (3 dB) before the data were downsampled to 201 Hz and recorded.Both of these filters are causal and shift the detector response in time.They must be deconvolved before further analysis, which requires synchronizing the detector response with the VPM encoder data and detector pointing data.The first panel of Figure 6 illustrates these steps.The Butterworth filter is a known property of the MCE and was deconvolved as the first step after the data calibration (blue curve).Due to the presence of the VSS, the effect of the time-constant convolution shows up as a distinct hysteresis between the VSS and the VPM mirror motion (Appel et al. 2019).Taking ad-vantage of this, the time constants were measured by minimizing the hysteresis and were deconvolved from the data as shown by the orange curve.The detector time constants were found to be mostly stable; therefore, we chose the average value per observation era per detector for the analysis.We assess the potential systematic errors from this choice in Section 5.1.2.

VPM and Demodulation
Polarized sky signals as seen by the CLASS detectors are modulated by the front-end VPM; therefore, a demodulation step is needed to recover the sky signal before making maps.This section provides an overview of the VPM modulation physics and the demodulation processing.In-depth descriptions are presented in Chuss et al. (2006Chuss et al. ( , 2012b)), Harrington (2018), and Harrington et al. (2021); the summary here highlights the calibration of the VPM transfer function for cosmological analysis.

Monochromatic VPM Modulation
For monochromatic radiation at frequency ν, the intensity reaching each linearly polarized detector, I ν , as a function of time t and the grid-mirror distance z(t) is:5 where A z (z)i VPM includes the emission from the VPM grid and mirror that makes up the majority of the VSS at 40 GHz; S i,q,u,v terms are modulation functions, which we describe below; the lowercase i, q, u, v are used for Stokes I, Q,U,V -like signals in the VPM coordinate system, where u is the polarization state that is modulated by the VPM, and q is the orthogonal state.The linear polarization recorded by each detector after modulation by a common VPM is related to the sky polarizations [Q,U,V ] by a rotation angle ϕ P such that: where γ is the detector position angle on the sky, and ϕ P corresponds to the VPM wire direction as seen from each detector (Figure 7), both of which are functions of the detector pointing offset from the center of the focal plane.The throw z(t) is changing at about 10 Hz so that the polarization modulation functions S u/v move the signal in u(t) and v(t) to higher frequencies (around the harmonic series of the VPM fundamental frequency), away from the slowly varying noise component in the unpolarized term i(t) (Tatarskii 1961; Church 1995;Morris et al. 2022).The throw range is chosen to optimize the linear polarization observation so that S u/v is mainly nonzero around 10/20 Hz.For an ideal VPM with zero emissivity, A z = 0 and i, q are not modulated (S i = 1 and S q is constant).The microwave frequency dependence of the Stokes parameters and the modulation functions in Equation 1 should be implicitly assumed.We refer the reader to Harrington (2018) for detailed derivations of the modulation functions S u,v .

VPM Transfer Functions
The modulated power D detected by the TES bolometers is the bandpass and beam integrated version of Equation 1: where η is the optical efficiency that accounts for all loss terms in the optics, B ν (Ω) and A e (ν) are the beam profile and the effective telescope aperture area, f ν is the detector bandpass, and λ is the wavelength.The identity (Pawsey & Bracewell 1955) is assumed for a single-mode detector observing a beamfilling source (e.g., CMB and extended foreground emission).
For unpolarized radiation, the received power is calibrated to the thermodynamic temperature unit by assuming a blackbody spectrum I BB : Here, we use ∆ to denote the spatial fluctuation of the quantity; x ≡ hν/k B T CMB where h and k B are the Planck and Boltzmann constants, and T CMB is the CMB temperature (Mather et al. 1994).The coefficient ∆T CMB /∆D i is the intensity calibration factor obtained from the Moon observations (Appel et al. 2019;Dahal et al. 2022).
For polarization modulated by a VPM, Equation 3 can be formally cast into a matrix product where Here we have ignored the unmodulated linear polarization component q(t) and substituted the polarization terms u(t) and v(t) from Equation 1 into Equation 3. The polarization intensity u(t)/v(t) are separated as the product of their effective temperature T u/v,s evaluated at an arbitrary reference frequency, and the associated spectral shapes s u/v (ν) such that u(ν) ≡ 2k B T u,s s u (ν)λ −2 (likewise for v(ν)).Note that the effective temperature depends on the spectrum of the source, which may be different for linear and circular polarization and might vary across the sky.For cosmological analysis, the source effective temperature is calibrated to thermodynamic temperature T u/v6 through Therefore, the calibrated transfer function relates the modulated power to polarization intensity in thermodynamic units, where the first term corresponds to the power-to-thermodynamic temperature unit conversion factor for intensity, and the second term is the uncalibrated transfer function.

Demodulation
The modulated time stream calibrated in thermodynamic temperature units is7 where the bold vector symbol is used to emphasize that each of the quantities is a time stream of length n samp .The transfer function time streams M u/v are evaluated for the VPM gridmirror distance encoder data (z(t), synchronously sampled with the data).The goal of demodulation is to solve the sky polarization signal T u/v from the raw time stream D. The least-squares solution to Equation 12 is (see also, Harrington et al. 2021) where M = [M u M v ] is the matrix form of the transfer function.Prior to the least-squares solving, a bandpass filter was applied that includes the first five harmonics of the VPM frequency (∼ 10 − 50 Hz) with a margin of f c = 0.5 Hz for the 1 deg s −1 az scan and f c = 1.0 Hz for the 2 deg s −1 scan, as informed by the beam-crossing timescale of the 40 GHz telescope.The bandpass filter was applied to both M and D so that the effect of this filter does not bias the solution.The solution was then filtered with an antialiasing low-pass filter (L) with cutoff frequency f c before downsampling to data rates of 1.45 Hz for the 1 deg s −1 az scan and 2.42 Hz for the 2 deg s −1 scan.Unlike the bandpass filter, which is accounted for in the demodulation least-squares solution, the low-pass filter in principle can remove signal as it is not accounted for in the mapmaking step.In practice, however, the cutoff frequency f c was chosen to have minimal suppression of the signal beyond beam smoothing; its residual impact on the mapping transfer function is characterized in Section 4.6.The middle panel of Figure 6 visualizes this process.The least-squares fit finds the best solution that matches the amplitude of the u/v transfer functions, which are interpreted as the sky polarization intensity through Equation 2.
A preliminary gap-filling was performed after the demodulation to fix the discontinuity of the data due to data selection (Section 3.2) and splitting of data for consistency tests (Section 5.2).The demodulated data were separated into low-and high-passed components by a rolling top-hat kernel of 50 samples in width.The gaps in the demodulated data were filled with the linear interpolation of the low-pass filtered component and then injected with white noise sampled from the white noise estimation of the high-pass filtered component.This treatment ensures the basic continuity of the data for subsequent operations but does not necessarily preserve the low-frequency noise structure.A dedicated gapfilling for this purpose is introduced in Section 4.3.

VPM Parameters
Calibration of the VPM parameters is essential to the recovery of the polarization signals.The parameters used for cosmological analysis were determined through a minimization process of mixing between the linear and circular polarizations (polarization leakage).The model parameters included the incident angle of light onto the VPM (per VPM grid period), an overall offset in the grid-mirror distance (per grid period), an overall shift in the bandpass centroid (per Era), the spectral index of the linear polarization, and an additional power-law correction to the atmospheric circular polarization spectrum to account for model uncertainties (Petroff et al. 2020b).
For each set of parameters, the demodulated data u/v were swapped such that the linear polarization stream u was mapped into horizontal coordinates as an intensity-like signal (to probe for leakage from atmospheric circular polarization into u) and the circular polarization stream v was mapped into equatorial coordinates as linear polarization signals (to probe for leakage from Galactic linear polarization into v).We defined the following polarization-leakage statistic: where for each of the three VPM grids labeled by j, A u is the amplitude of the atmospheric circular polarization model (Petroff et al. 2020b, see their Figure 3) fit to the horizontalcoordinate maps created from linear polarization (u), and A v is the correlation of the WMAP Q-band linear polarization maps around the Galactic plane with the "linear polarization" equatorial-coordinate maps generated from circular polarization (v).For linear polarization, the (quasi-)azimuthsynchronous signals (Section 4.1) were not filtered and left systematics in the horizontal coordinates; therefore, the azimuth range −10 • to 110 • was avoided when computing the circular polarization model amplitude A u .For the circular polarization, maps were made separately for detectors on the top/bottom of the focal plane (see A v,top versus A v,bot ) to better break the degeneracy between the tilt of the VPM and the grid-mirror distance.The uncertainties of these leakage amplitudes σ were evaluated with noise-only simulations.
The χ 2 values over the entire VPM parameter space were first sparsely explored with 250 Latin hypercube samples (McKay et al. 2000).Demodulated data and maps were made for each of the samples, and the χ 2 values were computed according to Equation 14.Another 250 samples were drawn near the (approximate) minimum of χ 2 and evaluated in the same way.With a total of 500 evaluations, the rest of the parameter space was parameterized with Gaussian process regression.Finally, we used the Markov Chain Monte Carlo method to determine the best-fit parameters that globally minimize the χ 2 values.
For the cosmology maps, we adopted the instrumental parameters and the circular polarization spectrum correction from the minimization process above but used the linear polarization spectrum from Krachmalnicoff et al. (2018) with a spectral index of −0.7 because the minimization process focuses on the Galactic region where the linear polarization index could be different from the rest part of the sky due to the variation of the synchrotron index (predominately from synchrotron; Gold et al. 2009;Fuskeland et al. 2014;Krachmalnicoff et al. 2018;Planck Collaboration et al. 2020c;Rubiño-Martín et al. 2023) and the mixture with the CMB.The impact of the uncertainties of the linear polarization spectral index and the other parameters from the leakage minimization are further characterized in Section 5.1.5.

MAPMAKING
The raw CLASS data can be formally modeled as where M(t) is the modulation transfer function introduced in Equation 12, P(t) is the pointing matrix that transforms polarization Stokes parameters from the sky coordinates to the VPM coordinates (i.e., the matrix in Equation 2), m = [Q,U,V ] are the sky polarization maps, and n represents the raw data noise.The demodulation process described in the previous section partially solves this equation and yields intermediate demodulated data where n u and n v are noise in the demodulated data.This section describes the process of solving the polarization maps from d(t).We start with filtering the demodulated data to reduce systematics not accounted for by the data model.We then describe the noise model and gap-filling methods for the demodulated time streams and how they are applied for maximum-likelihood mapmaking.Finally, we present the maps and the associated transfer functions due to the filtering.The CLASS mapmaking algorithms are developed based upon the public code minkasi8 (Romero et al. 2020).

Filtering
The modulation technique and the demodulation process produce polarization time streams that are mostly immune to atmospheric fluctuations and intensity-like systematics from the sky and the ground.However, these data were found to have systematic signals that may be traced back to polarized environmental emission, T -to-P leakage, electronic pickup from the instrument, etc. Time-domain filters were designed for each of the cases and were jointly fit and removed from the polarization time streams as where F is a collection of column vectors that include all systematic signal models, which we describe in the rest of this section; M is a time-domain mask to prevent biasing the filter.The mask comprises the TOD selection mask (Section 3.2) and a linear polarization mask that vetoes the brightest 3.6% of the sky in synchrotron polarization (Planck Collaboration et al. 2020c) around the Galactic plane.Nevertheless, this approach filtered out sky modes that mimic systematics in the time domain, especially at large angular scales.
The impact is quantified by the mapping transfer functions in Section 4.6.

Azimuth Servo Motor Signal
A spurious signal was found to be synchronous with the az-servo motor current (thus, with the az-velocity of the telescope), having a peak-to-peak amplitude up to several times the VSS amplitude (∼ 5 fW).A set of harmonic components was filtered to remove this signal: Here, ϕ 8π is the azimuth of the telescope with the 8π period (which accounts for both the positive and negative azimuth velocity sweeps).This component was filtered from the linear and circular polarization data every 3 hr.

Wind-induced Signal
A quasi-azimuth-synchronous signal is present in the demodulated linear polarization data that correlates with the wind recorded by the WeatherHawk 10 weather station installed close to the CLASS telescope.The weather station provides wind speed information through a cup anemometer (starting threshold at 0.78 m s −1 ) and wind direction through a vane.
Figure 8 shows the demodulated linear polarization signal for each of the −45 • oriented detectors as a function of the bearing angle of the wind with respect to the telescope azimuth pointing and the wind speed (the radial axis).The quadrupole feature across the focal plane that peaks around 0 • and at high wind speed is due to the deformation of the plastic closeout film at the telescope's optical entrance when pressed by the wind.The blue/red features toward the south are due to arbitrary baseline adjustments for this plot.The 9 The bracket with sub/superscripts indicates a set of filter basis and its parameter range. 10http://www.weatherhawk.com/s232dcwind signal was found to be consistent at different boresight configurations and scales roughly linearly with the wind speed.As shown in Figure 9, the prevailing wind at the site came from the northwest during the austral winter and had significant contributions from the east in summer, with a slight shifting in direction throughout the night (not shown in the plot); therefore, this wind-induced signal left a systematic error that is covariant with azimuth pointing and thus with sky signals.
Instrumental mitigation of this issue started in September 2021 by removing the covering plastic during observations.For the time period with the closeout on, the filter components took the form: where γ w and ϕ w are the wind speed and wind direction measured by the weather station, ϕ(t) is the azimuth pointing of the telescope with 2π period.Wind data were up-sampled from the original rate at 0.5 Hz to align with the demodulated data.This filter was applied to the linear polarization time streams every 2 hr.

Azimuth-synchronous Signal
Several systematic issues were found to contribute to an azimuth-synchronous signal.The metal surface of the telescope cage reflected ground emission and produced signals in both intensity and linear polarization.This was mitigated sequentially by blackening the interior of the cage and the forebaffle with microwave absorbers.At the beginning of Era 2, the circular forebaffle extension was replaced by the double baffle, which has an asymmetric shape and a larger opening angle to accommodate the new 90 GHz telescope.The forebaffle extension blackening was removed at the beginning of 2019.Figure 10 presents an example of this signal in the linear polarization stream for a pair of detectors in Era 1 before and after the blackening. 11The linear polarization signal at each boresight angle is binned in the telescope azimuth coordinates as seven separate rings.The reflection picked up terrestrial emission that correlated with the Cerro Toco mountain toward the northeast and depended on the boresight rotation of the telescope (which changed both the polarization angle and the pointing elevation of the detector).The atmospheric circular polarization due to the Zeeman splitting of molecular oxygen magnetic dipole transitions is defined by the azimuth angle from the magnetic North and the pointing elevation (Petroff et al. 2020b).This is a smooth azsynchronous systematic for sky circular polarization and was also a potential bias to the linear polarization measurement through leakage from imperfect modeling of the VPM transfer function (Section 3.6.3).The aforementioned wind signal The panels are laid out by the detector pointing offsets, similar to Figure 7, with the detectors pointing at higher elevations at the top.Within each panel, the wind-related systematics are mostly confined in the quadrant when the telescope is facing toward the wind.The signal amplitude varies across the focal plane with a quadrupole pattern that increases from the center toward the edge of the focal plane and also rises with the wind speed.The whirlpool-shaped structure in some of the edge detectors is due to the ground pickup (Section 4.1.3)may also have left a residual in the azimuth due to inaccuracy in the wind data and modeling errors.
In addition to these relatively stable components, it was found that the electronic coupling to the detector caused an azimuth-synchronous signal that varies on time scales of a few hours.This is likely related to the wiring of the detectors since the detectors on the right-hand side of the focal plane, which connect to shorter cryogenic wires, show more stable signals at lower amplitudes.This electronic pickup was also notably improved since the deployment of Era 2.
A set of harmonic filter components was employed to mitigate these signals: where the u/v time streams were fit with 15/10 harmonic components, respectively.The filters for circular polarization were evaluated every 3 (4) hr for Era 1(2); for linear polarization, the timescale was chosen to be 2 hr for the left detectors, and 3 (6) hr for Era 1(2) for the right detectors.The amplitude of each of the components in F az was determined separately for the positive and negative az-velocity regions to further mitigate the difference in the az-synchronous signals that were correlated with the az-servo motor.Of all of the filters described in this section, these removed the most celestial signal.

Camera RFI
Beginning on 2017 August 9 through 2018 June, the 40 GHz detectors experienced RFI pickup from a camera installed inside the mount cage.The RFI (and its harmonics) around the VPM modulation frequencies produced a slowly varying harmonic structure in the demodulated polarization time streams.Since linear and circular polarization were predominantly modulated by the first and second harmonics of the VPM modulation frequency, the contamination in the demodulated data was mostly confined around 0.51 and 1.53 mHz for linear polarization, and around 1.02 mHz for circular polarization.
For the time period with cage cameras on during observations (Section 2.3), a set of harmonic lines at these frequencies was fit and removed from the demodulated data every 3 hr.This is a gentle filter to the sky signal as the frequencies are below the scanning frequency at 2.7 mHz (for 1 deg s −1 scan).

Noise Model
Like other ground-based CMB experiments, CLASS (demodulated) data are noise dominated; therefore, the noise term in Equation 16needs to be carefully modeled to achieve optimal sensitivity in the maps.
The demodulated noise n u and n v were similarly white at high frequencies and correlated over long time scales ( f knee , ≈ 5 − 20 mHz, Harrington et al. 2021;Cleary et al. 2022) but were also distinct due to the difference in the modulation functions and the nature of linear and circular polarization.In particular, atmospheric signals sourced long-timescale correlated noise across the focal plane (Tatarskii 1961;Church 1995;Wollack et al. 1997;Lay & Halverson 2000;Errard et al. 2015;Morris et al. 2022), some of which were linearly polarized (Takakura et al. 2019), or could have impacted linear polarization through T -to-P leakage.The emission from the VPM was predominantly covariant with the linear polarization signal, and its slow temporal variation was another potential source of long-term instability (Miller et al. 2016;Harrington et al. 2021).The electronics in the readout system (Reintsema et al. 2003;Doriese et al. 2016) could also have contributed to the correlated noise in the demodulated data, which was mainly manifest as common correlation features between pairs of detectors and, to a lesser extent, within each readout column (Dünner et al. 2013;Harrington et al. 2021).Due to the covariance of the modulation transfer functions, n u and n v were also expected to be correlated at all scales.
Formally, the noise model takes the following form in Fourier space: where N is the covariance matrix among all detectors and between linear and circular polarizations, N M contains the power spectra of the common modes that are projected to each detector through V, and N D is a collection of power spectra per feedhorn.The construction of the common modes N M is informed by the noise properties of the data.Following Dünner et al. (2013), a singular value decomposition (SVD) was performed on the low-frequency (below 0.1 Hz) part of the data to identify modes with dominant singular values (above 3.5 times the median); subsequently, a second SVD was applied on the full frequency range to further find common modes after the removal of the low-frequency modes.An SVD per readout column was then used to search for residual correlated modes unique to each readout column.
There are eight columns in the 40 GHz telescope focal plane, each containing eight or ten optical detectors.This hierarchical construction typically found ∼ 20 modes among the 144 detector-Stokes time streams in total.The per-feedhorn component N D captured the rest of the noise power as a blockdiagonal matrix with four-by-four blocks that describe the covariance between the two polarization states of the two paired detectors associated with the same feedhorn.This noise model was estimated for all the data over every 10 sweeps, which is about 2 hr for 1 deg s −1 scan.To facilitate fast evaluations of the noise model and its inversion, the spectra were logarithmically (linearly) binned above (below) 50 mHz.Despite the filtering described in Section 4.1, the azimuth-synchronous systematics were not completely removed as the shapes of the ground pickup, wind-induced signals, and electronic coupling gradually vary over time, and contribute to the noise model.The binned power spectrum has the advantage of capturing this residual power and allows for the down weighting of the data at around the scanning frequencies.The matrix construction above ensures a low-rank N M and trivially invertible N D , allowing for an efficient inversion using the Woodbury identity (Woodbury 1950).
The noise model was used to optimally weight the data for mapmaking and can be directly sampled to create noise simulations.However, the model estimated above can be biased due to (1) the missing data in the time streams from data selection, and (2) the direct estimation of the covariance matrix from data that contains both noise and signal.The second type of bias would further bias the sky signal estimation when used in Equation 16.These issues are addressed, respectively, with iterative methods in Sections 4.3 and 4.4.

Gap-filling
As mentioned in Section 3.6.3,the preliminary gap-filling does not preserve the low-frequency noise properties of the data within the gap, and the noise model directly estimated from this may be biased.To improve this, we express the mask-aware data likelihood as by introducing the better gap-filled data d.Here, N is the noise model estimated from the preliminary gap-filled data or from a previous iteration, and N ∞ is the noise model for the masked data that has an infinite variance for samples within the mask and zero elsewhere.The second row expresses the likelihood function as the marginalized conditional probability against an auxiliary variable t that has noise properties described by the white noise component of the noise model N w and its residual, d − t, that follows the red noise part N r = N − N w (Huffenberger & Naess 2018).This formalism permits the "regeneration" of the gapfilled data by Gibbs sampling of the conditional probability functions.After 10 steps of sampling, we took the resultant d as the updated version of the gap-filled data for another iteration of noise model estimation.Based on the Kullback-Leibler divergence between the result and simulation inputs with known noise properties, this iteration converged quickly most of the time when the mask fraction was low and when the masking was not correlated among the detectors.
For some of the null tests (E23) where the data were heavily masked (50%) by splits between, e.g., positive and negative azimuth velocity scans, longer iterations were needed to obtain an accurate recovery of the noise model.We chose to run this iteration five times for each noise model to accommodate the extreme cases.The gap-filled data from Gibbs sampling were only used for updating the noise model and were not projected into the maps.

Maximum-likelihood Mapmaking
The maximum-likelihood solution to Equation 16given the noise model N from Sections 4.2 and 4.3 is This was solved iteratively using the preconditionedconjugate gradient method (Shewchuk et al. 1994), where we used the inverse of the hits map as the preconditioner Here, M and P are the modulation transfer function and the projection matrix defined in Equation 15.This is a proxy for the covariance of the sky map assuming constant white noise in the raw data around the modulation frequencies.Using it as the preconditioner enabled fast convergence within 50 steps of conjugate gradient iteration.Note that m solved from Equation 25is biased due to the filtering applied in Section 4.1, and this bias will be characterized in Section 4.6.
The healpix pixelization at N side = 128 was used for all map products at 40 GHz, unless otherwise noted.
The noise covariance estimate N in Section 4.2 has ignored the signal term in Equation 16.While this is a good approximation, as the per-sweep data have very low S/N, it would slightly bias the signal map, especially on large angular scales, where the signal is degenerate with the correlated component in the noise model.To mitigate this bias, we performed multiple template iterations by projecting out the estimated sky signal from the data for an updated noise estimation (as illustrated in Figure 4, without additional Gibbs-sampling gap-filling).This technique has also been employed in experiments with similar mapmaking strategies (Dünner et al. 2013;Aiola et al. 2020;Romero et al. 2020).Figure 11 demonstrates this effect by showing the signal power spectra of the mapping pipeline, estimated from the cross-correlation of two independent splits at each template iteration.The input signals were Gaussian realizations from power-law spectra following the best-fit synchrotron model of Planck (Planck Collaboration et al. 2020c) with noise from the noise model estimated from the demodulated data.The bias on large angular scales was corrected for over a few template iterations; we found that using five iterations was sufficient for the 40 GHz maps with a remaining bias below 2.5% for ℓ ≤ 5.The remaining bias was verified to be caused by noise modeling since simulations using the input noise model (or other static noise models independent  The bottom panels show the residual spectra between the output spectra and the input spectra normalized by the input amplitude, in the same color scheme.The correction converged quickly at high ℓ where the corresponding time streams are noise dominated; at low-ℓ the iteration converged by the fifth iteration, with a bias below 2.5% for ℓ ≤ 5, much smaller compared to the sample variance. of the data) for weighting (i.e., Equation 25) did not show this deficit of power at low ℓ.The 40 GHz polarization maps did not suffer from the large-scale bias due to subpixel errors as pointed out in Naess & Louis (2022), since CLASS demodulated data have a small dynamic range between largeand small-scale noise (Harrington et al. 2021, Cleary et al. in prep.); this prevents the subpixel residuals from outweighing the noise model at low frequencies and causing large-scale bias.This was verified by simulations that use higher-pixelresolution sky maps as input.

Maps
The 40 GHz linear and circular polarization maps made from the CLASS observations through 2022 are presented in Figure 12.A battery of self-consistency null tests and a comparison with satellite missions will be presented in E23.In Figure 13, we show the hits count map defined in Equation 26.The diagonal components are the integration time of each Stokes map; the off-diagonal terms reflect the covariance between maps.The u-v covariance through the VPM modulation is integrated down in QV but not in the UV component due to the projection effect. 12The QU component has minimal covariance due to the design of the CLASS scanning strategy.In addition, the bottom-left corner of the figure is a cross-linking map from the CLASS scanning strategy, defined as where the sum is over the N TOD falling within a pixel, and Γ j is the angle between the scanning direction with respect to the local meridian for the j th sample (Aiola et al. 2020;McCallum et al. 2021).This value reflects the uniformity of the scanning direction coverage and is a good proxy for the (inverse) large-scale noise related to atmospheric emissions (Atkins et al. 2023).Together, the cross-linking and hits maps show complementary information about the CLASS sensitivity on the sky at different scales.
The noise power spectra estimated from simulations are shown in Figure 15, which reach white noise levels 110 µK arcmin in EE/BB/VV .After correcting for the mapping transfer function, the large-scale noise has a logarithmic slope close to −2.4 and knee angular scales (at which the spatially correlated noise equals the white noise) of ℓ = 12 and 18 for circular and linear polarization, respectively.

Harmonic-domain Transfer Matrix
The filtering performed on both raw data (Section 3.6) and demodulated data (Section 4.1) removed power from the sky signal.This effect can be modeled as where ℓ and ℓ ′ are extended multipole indices that run through multipole moments of {VV, EE, BB}, and F ℓℓ ′ is the mapping transfer matrix.Here we have assumed isotropic filtering on each mode, which is only a good approximation for statistically isotropic sources at small angular scales.Despite this, we found through simulations that the resulting harmonic transfer function results in unbiased spectra for ℓ > 4. Estimation of the filter transfer matrix was obtained by mapping signal simulations with known input spectra, performing the same filtering and noise weighting as the data, and comparing the resultant spectra.The off-diagonal components, i.e., mode mixing among adjacent multipoles and among Stokes parameters, were found to be mostly insignificant; only the covariance over 10 adjacent multipoles in the EE and BB blocks and five adjacent multipoles below ℓ = 30 in the EE −BB cross blocks were modeled, and the remaining elements were fixed to zero. Figure 14 shows the transfer function due to low-pass filtering in demodulation (Equation 13) and the diagonal com-  ponent of the mapping transfer matrix F ℓℓ .For linear polarization, about 35% of the power is retained at ℓ = 10, and the signal sensitivity peaks at ℓ ≈ 40 (after accounting for the beam window function and noise power spectrum).

Pixel-space Transfer Matrix
At large angular scales, the anisotropic effect of both the foreground signals and filtering is more prominent, and the harmonic-domain transfer matrix is a less-robust representation.A pixel-space transfer matrix can be introduced at low resolution for this situation: where m in/out are the input and filtered maps downgraded to N side = 16 resolution, F i j is the transfer matrix estimated from an ensemble of signal-only simulations, and the subscripts denote the map pixel index.Although the mapmaking pipeline is linear, this equation is only an approximation for the downgraded maps due to the noncommutativity of the mapping and the downgrade operation and showed a ∼ 10% discrepancy at the largest angular scales when compared to the reobservation (Section 4.7).This can be improved in the future by using a separate low-resolution mapping pipeline.This transfer matrix can be used for pixel-space analyses, and it can be integrated with quadratic estimators (e.g., Vanneste et al. 2018) to optimally correct for the bias in the power spectra.For the latter, we found that the transferfunction-corrected power spectra are unbiased for ℓ ≥ 4, comparable to the pseudo-C ℓ estimator with the harmonicdomain transfer function correction, but that they are statistically more optimal at low-ℓ than the pseudo-C ℓ approach.

Reobservation
To facilitate direct comparison with other experiments, in particular the all-sky maps from WMAP and Planck, we applied the CLASS filtering and weighting on these maps to forward-model the filtering effects.The reobservation started with convolving the input map to 1.5 • (FWHM) resolution.The linear polarization components of the input map were then projected to the CLASS demodulated time streams 13 using the CLASS pointing model, and the circular polarization time streams were set to zero.These data were then filtered in the same ways as the CLASS demodulated data and projected back to the maps using the fixed noise model from the last template iteration of the CLASS mapmaking procedure-no template iteration was performed for reobservation.Since all mapping operations are linear, this is an accurate description of the filtering that the CLASS data have undergone. 13In principle, the reobservation should go through the modulation and demodulation process as well; however, in practice, it is much more efficient to start the simulation from the demodulated stage.Although the filtering effect due to the low-pass filtering in the demodulation (the pink curve in Figure 14) is not captured, its effect is subdominant compared to the beam transfer function and is safe to be neglected.

DATA VALIDATION AND SYSTEMATICS 5.1. Systematics and Simulation
In this section, we characterize several types of systematic errors and assess their impact on the scientific result.These issues were studied through simulations with CMB realizations as input, and the resultant bias to the power spectrum was characterized by the difference between the systematicsincluded auto-power spectrum and the input power spectrum.The simulations were drawn from the Planck best-fit parameters (Planck Collaboration et al. 2020a) with the B-mode amplitudes set to zero.Since the input had no power in the B mode and in circular polarization, the effects in BB and VV were dominated by the auto-correlation of the systematics residuals, i.e., a second-order effect of the systematics, while the residual EE spectrum were dominated by the crosscorrelation between the residual E model power and the original signal, i.e., a first-order effect.Figure 15 summarizes the results, which we describe below.

Detector Polarization Angles
Calibration of the absolute polarization angle is critical for accurate separation of the E/B-mode signal (Hu et al. 2003) and the search for parity-violating physics (Finelli & Galaverni 2009).Systematic uncertainties associated with the alignment of individual detector pixels, offsets between the focal plane and the VPM wire, pointing errors in the telescope boresight rotation, and modeling errors in the optics can lead to bias in the polarization angle.
We used the bright polarized source Tau A as the main calibrator whose polarization angle was measured by CLASS to be −87.02± 0.2 • in Galactic coordinates 14 where the statistical error is derived from noise simulations.However, as shown in Figure 7, only the central bottom part of the focal plane covered Tau A at all boresight rotations; detectors on the sides observed Tau A only at certain boresight rotations, and part of the top detectors never saw Tau A. This boresightdependent partial coverage of Tau A limits its ability to characterize systematic errors of the polarization angle.Based on the optics of the telescope and the distribution of the wire direction ϕ P (Equation 2) across the focal plane, we used the discrepancy in the Tau-A polarization angle measured between splits of the data to assess the systematic errors.The "quadrupole" split defined by the sign of ϕ P (see details in E23) probes the systematic effects in the optics modeling.Similarly, a split between scans with positive and negative boresight rotation relies on either side of the blue detectors in Figure 7 to measure Tau A and is therefore sensitive to the optics model as well.The Tau A polarization angle differences in these two ways of splitting the data are 0.70 • (quadrupole) and 0.80 • (boresight).The Tau A angle measured from the data using each of the three VPM grids shows a maximum discrepancy of 0.37 • , which indicates the level of the error caused by an angular offset between the VPM grid and the 14 The polarization angles in this paper follow the IAU convention.focal plane.Combining these two factors, we assign 0.7 • to the systematic error in the polarization angle calibration.The final measurement of Tau A polarization angle from CLASS, −87.02 ± 0.20 (stats.)± 0.70 (sys.)• , is consistent with that from WMAP (−87.3 ± 0.2 ± 1.5  15 as the pink curve based on the analytical model (Keating et al. 2013).

Time Constants
The time constants used for data reduction are the medians of values estimated from each DataPkg per detector per observation era.So for each detector, there is a single time constant estimate used for all data in Era 1 and another for Era 2. The time constants among the detectors have a typical value of 3 ms (median), but are sensitive to the optical loading from the sky and show correlations with the air temperature, PWV, and the telescope boresight rotation.Most notably, air temperature accounts for a shift of 0.2 ms (comparable to the standard error of the time constant estimations) in the time constants of all detectors over its normal range (−10 • C, 6 • C).The time constants are also affected by the thermal history of the detectors; some detectors jump between states of time constants differing by approximately 0.2 ms when the focal plane temperature warms up above 0.1 K.In the following analysis, we take half the value of 0.1 ms to study the impact of a systematic bias.The statistical uncertainties from averaging all time constants per era are insignificant compared to this.
These variations are not considered in the pipeline, and deconvolution with biased detector time constants has dual impacts on the CLASS data.First, the actual pointings of the telescope would be offset from the calibrated encoder values, but this effect is negligible compared to the beam scale of the telescope and is further diminished by the forward and back-ward scanning of the telescope.More importantly, the phase delay from the VPM encoder due to the biased time constants would cause leakage between the linear and circular polarization.Figure 15 shows the effect of the time constants biased by 0.1 ms.This is an approximately 10 −3 effect in the mixing of the polarization states; therefore, for simulation with pure E mode input, the residuals in the EE/VV power spectra are at the 10 −3 /10 −6 level, respectively.No E-to-B leakage is detected above the level of the lensing B-mode at the largest angular scales.

Ghosting Beam
Beam ghosting caused by the internal reflection of the telescope is detected in Moon observations at a level 3 × 10 −3 of the main beam at the opposite position of the focal plane for each detector.To simulate this effect, a Gaussian beam centered on the opposite point of the focal plane was assigned to each detector, with the peak amplitude of the beam consistent with the reflection amplitude measured from the Moon maps.These ghosting beams were then convolved with the sky simulations for mapping.The residual power spectra of the ghosting beam, which are shown in red in Figure 15, have the greatest impact at angular scales greater than the field of view of the telescope (≳ 20 • ).

Temperature-to-polarization Leakage
The placement of VPM as the first element in the optical path is to prevent polarization due to oblique reflection from being modulated.Moon maps made from dedicated Moon scans (Xu et al. 2020) have shown that the monopole T -to-P leakage is at the 4 × 10 −5 level for pairs of detectors.The Moon maps also reveal a dipolar pattern that takes opposite signs in linear polarization for a pair of detectors with its amplitude and orientation independent of the telescope boresight rotation (top row of Figure 16).Similar patterns are also observed for circular polarization, but at lower magnitudes and with the orientations of the dipole offset by approximately 45 • from those in linear polarization.This effect is consistent with a misalignment of the VPM, where the tilt between the grid and the mirror creates differential pointing and leads to an additional term in Equation 1 proportional to the brightness temperature gradient along the tilt direction.This term is modulated at the VPM frequency and is covariant with the linear (primarily) and circular polarization modulation function and is picked up by demodulation (Section 4.1.2 of Harrington 2018).
The bottom-left panel of Figure 16 shows the pair-null linear polarization map made by differencing polarization maps made with the +45 • detectors from those made with the −45 • detectors.Because of polarization modulation, the linear polarization can be recovered with the +45 • and −45 • detectors separately.Thus, differencing removes the polarization signal, and enhances the oppositely signed dipolar T -to-P leakage (Figure 16, top row).Furthermore, the T -to-P effect can be simulated by convolving the sky temperature signal with the dipolar beam estimated from the Moon maps for each detector and each VPM grid.The convolution was performed in the pixel space using pisco (Fluxá et al. 2020).The bottomright of Figure 16 shows a simulation of the bottom-left panel made by convolving the WMAP Q-band temperature maps scaled to the CLASS bandpass by the dipolar leakage beams and differencing the simulated +45 • and −45 • leakage maps.The agreement between the data and simulation shows that the dipolar leakage measured from the moon is in agreement with that measured in the CLASS survey maps via Tau A.
Although the dipolar leakage is on average 0.3% compared to the main beam for a single detector, its impact on the final maps is further diminished when polarization data from pairs of +45 • /−45 • detectors are averaged, and the dipoles in the top row of Figure 16 cancel each other (instead of reinforce, as in the bottom row).Only eight of the 72 detectors in Era 1 were unpaired (due to readout failures), six (two) of which are +45 • (−45 • ) oriented, and this was reduced to a single unpaired detector in Era 2 (due to data selection).
To assess the impact of this leakage in the angular power spectrum, we convolved CMB temperature map simulations with the dipole beam in linear and circular polarization.We then took cross spectra between the output maps and the sum of the input and output (i.e., the main beam plus the dipole beam).The resultant EE power spectrum has a contribution from the auto-correlation of the dipole systematics and the cross-correlation with the CMB E mode due to the CMB T E correlation; the effects on the BB/VV power spectra are solely from the auto-correlation.These results are shown in Figure 15 as blue curves, and in all cases, the effects of the dipole T -to-P leakage are subdominant compared to the noise level and/or the cosmology signals of interest.

VPM Transfer Function Uncertainty
The best-fit VPM transfer function parameters were determined with a combination of instrument characterization and a polarization-leakage minimization process as outlined in Section 3.6.4.We assessed the impact of VPM parameter uncertainties by modulating a single realization of the sky signal with different VPM parameters drawn from the likelihood chain in Section 3.6.4that are within the 95% confidence region around the best fit.These simulations were then demodulated with the best-fit VPM parameters and mapped in the same way as the data.The maximum absolute differences between the output and input power spectra are plotted in Figure 15 as the navy blue curves with shades.The VPM parameter uncertainties typically translate to a 1% error in map amplitudes.
The demodulation pipeline assumed a single spectral index for the linear polarization, which is a simplification to the real-world case where both the spectral index of the dominant synchrotron emission and the mixing between different components vary across the sky.The PySM (Thorne et al. 2017) simulation with realistic input from synchrotron, CMB, and dust suggests that the aggregated effect corresponds to a standard deviation across the sky of 0.3 at 40 GHz. Figure 15 shows in brown a simulation with a uniform +0.3 bias in the linear polarization spectral index; the leakage effect manifests as a transfer of the E mode power into circular polarization.This result should be considered conservative since the variation of the index bias over the sky should partially cancel out and leave less residual in the power spectra.

Filtering Artifacts
The demodulated data filtering performed before mapmaking removes and redistributes the sky signal over large angular scales.Although its effect in the harmonic domain has been modeled by the transfer matrix (Section 4.6), the insufficiency in the modeling could still lead to bias in the corrected power spectra.The purple curve in the BB panel of Figure 15 shows the transfer function-corrected BB power spectra from filtered E-mode-only simulations.

Internal Consistency Test
The validation of the CLASS data product is checked through a series of internal consistency tests.The tests are executed by splitting the demodulated data into two similarsized subsets (denoted A and B) that are expected to expose certain types of systematic error.Two homogeneous temporal split maps are made for each A/B split through the same mapmaking pipeline.The cross spectra of the difference (null) maps A − B from the two temporal splits are computed and compared to an ensemble of simulations to check for consistency.Details of the design of the split and the systematic errors probed by each null test, as well as the final results, are presented in E23.

CONCLUSION
We have presented a detailed description of the CLASS data reduction pipeline for 40 GHz observations conducted from August 2016 to May 2022.When weather, instrument upgrades, and other interruptions permitted, observations were conducted continuously, regardless of time-of-day or season-of-year.After all data cuts, the analysis incorporated 86.77 detector-years of data, representing ∼ 20% of the possible data volume.These data cover 75% of the sky and extend from −76 • to 30 • in declination.
The sky polarization signal in the data was amplitudemodulated at the VPM frequency (10 Hz) and its harmonics.Therefore, the selected data were filtered and demodulated to remove the dominant time-correlated noise in the raw data below 5 Hz and retain the polarization signal in the 0.5−1 Hz wide side-bands around the modulation harmonics.Isolating the polarization signal from the correlated noise in this way is the most important aspect of the CLASS strategy for achieving the large-angular-scale measurement.The demodulated polarization data were then filtered to remove systematic effects, such as azimuth-synchronous and windinduced signals.The type and level of filtering were tuned to enable the data to pass internal consistency "null" tests (E23).The demodulated and filtered data were then input to an iterative preconditioned-conjugate-gradient-descent algorithm to jointly solve for the maximum-likelihood Stokes Q, U, and V maps.Due to the filtering, the maps are biased low on large angular scales.We used simulations to show that the bias to the angular power spectra is ∼ 67% (∼ 85%) at ℓ = 20 and ∼ 35% (∼ 47%) at ℓ = 10 for linear (circular) polarization.After correcting this bias, the noise level in the angular power spectra was found in data-based simulations to be 110 µK arcmin [1 + (ℓ knee /ℓ) 2.4 ] 1/2 , with ℓ knee ≈ 12 for circular polarization and ℓ knee ≈ 18 for linear polarization.With these maps, CLASS is pushing the limits of what has been achieved from a suborbital platform at the largest angular scales.
Multiple sources of systematic error were quantitatively studied with simulations.The bias induced in the ΛCDM EE angular power spectrum was found to be subpercent.Leakage from the EE to BB spectra was found to be comparable to the predicted B mode spectrum with r = 0.01.Improvements in calibration and the data pipeline will reduce the leakage.For CLASS, the 40 GHz data are intended to measure the synchrotron foreground and not to constrain ΛCDM.Therefore, this study of the impact of systematic errors on the ΛCDM spectra is provided as an initial benchmark of the analysis on the way to analyzing the multifrequency dataset.Additionally, the subpercent bias found for the ΛCDM EE spectrum should be similar to the expected bias level for the EE and BB spectra of the diffuse synchrotron emission.This is the first demonstration of the full data pipeline for CLASS.At the time of writing, the methods developed here for demodulation and mapping were being applied, adapted, and improved on data from the other CLASS frequency bands.Several hardware improvements were made (guided by data), and both software and further hardware improvements were desired and planned.These results are therefore preliminary.Together with previous, ongoing, and planned improvements to the instrument and measurement strategy, future analyses will provide an independent view of the CMB polarization at the largest angular scales.

APPENDIX
A. POINTING MODEL The CLASS pointing model is a 34-parameter model for determining encoder positions to achieve a desired pointing in azimuth, elevation, and boresight angle.We present the model here, as it is unique in how it handles errors inherent to a three-axis mount.For a mount with one telescope, there is one pointing model, and the mount is positioned such that the commanded position corresponds to array center in azimuth and elevation at a boresight angle with respect to the zero boresight angle of the array.For a mount with two telescopes, each telescope has its own pointing model.While the individual models are used in the data analysis, the mount uses the average of the two for positioning such that the commanded position corresponds to a point on the sky halfway between the two array centers at a boresight angle with respect to a zero angle halfway between the zero boresight angle of the two arrays.Henceforth, when referring to the telescope mount, this average position will be denoted as the array center.

A.1. Boresight Pointing Model
The telescope mount boresight is defined as the axis of rotation of the boresight platform that houses the telescopes.As the boresight azimuth, elevation, and rotation angle, as read from the encoders, are not perfectly aligned with the sky, we need a pointing model to correct this misalignment.We use a boresight pointing model that contains 21 terms: 11 in azimuth and 10 in elevation shown in Table 2.Each term cor-responds to a physical effect that affects the alignment of the boresight.Of these 21 terms, 12 are currently in use; seven in azimuth and five in elevation.The four tilt terms are reduced to two coefficients in the pointing model data reduction, as they are not independent and are used in linear combination to describe the tilt of the mount as a rotation: ∆ az = α sin(el) cos(az) + β sin(el) sin(az) (A1) ∆ el = β cos(az) − α sin(az), (A2) where α represents a tilt of the mount to the West and β represents a tilt to the North.Here az and el are the commanded position of the mount boresight and ∆ az and ∆ el are the tiltrelated pointing corrections in azimuth and elevation.All azimuth pointing corrections are in units of true arc on the sky and are subsequently multiplied by sec(el) to yield the azimuth coordinate offsets that are applied to the azimuth axis encoder.In addition to the fixed tilt given by these coefficients, the mount employs a two-axis tilt meter.The signals from this tilt meter are passed through a 1 Hz low-pass filter and then used as additional corrections, which are applied to the encoders and recorded to disk for use in pointing reconstruction for data analysis.The tilt-meter pointing corrections are a combination of the residual tilt left over from leveling the mount, temporal tilts, and any zero offset of the meter itself.The most significant temporal tilt is an approximately 10 millidegree tilt away from the direction of the Sun during the day caused by the expansion of the sunlit side of the pedestal of the mount.

A.2. Boresight Angle Pointing Model
Since rotating the boresight platform by a given angle as read by the boresight axis encoder does not correspond to the true angle of rotation on the sky, we need a model to correct the boresight angle pointing.The model we use consists of parabolas in commanded boresight angle (bo) whose coefficients are themselves parabolas in elevation.This is best visualized as a 3 × 3 matrix as shown in Table 3.

A.3. Array Pointing Model
The array pointing model is described in the instrument frame in what we call the receiver coordinate system.This is a spherical coordinate system centered on the equator of the sphere at zero longitude.The principle axes are X and Y , where X corresponds to azimuth and Y corresponds to elevation at a boresight angle of zero.Each detector has an offset ∆x and ∆y with respect to the origin at array center.Since the fields of view of our telescopes are large, these offsets are computed using spherical trigonometry.The array center itself is not aligned with the boresight of the telescope mount, so we need a pointing model to describe the offset of each telescope's array center with respect to the boresight.While the array center offsets are small, we use spherical trigonometry to describe them for consistency with the detector offsets from the array center.The model uses two coefficients Azimuth Axis Warp a a These terms, while included in the model, were found to be insignificant and have their coefficients set to zero.b This term, while included in the model, was found to be degenerate with the gravity terms over the elevation range of the observations and has its coefficient set to zero. in each of X and Y as shown in Table 4.The elevation dependency is due to gravitational deflections in the optics that change with elevation.

A.4. Pointing Model Data Reduction
Pointing model data reduction begins with the analysis of a set of drift scans of the Moon during which the telescopes are scanned back and forth at a constant elevation while the Moon rises or sets through their fields of view.This analysis

Figure 1 .
Figure1.Schematic of the CLASS 40 GHz Telescope.The telescope sits atop the three-axis telescope mount in an enclosed cage structure.Lower-side cage panels are not shown at left to reveal the telescope.Light enters the enclosure by the forebaffle extension, first encountering the VPM, then the primary (obscured) and secondary mirrors, and finally the receiver.The enclosed multireflection design limits spurious signals from stray light.In 2018, the enclosure and forebaffle extension were extended to accommodate the 90 GHz telescope on the same mount.

Figure 2 .Figure 3 .
Figure 2. CLASS daily observation efficiencies from 2016-08-31 to 2022-05-19.The gray region is the total detector time recorded when the VPM is working; the yellow region shows the amount of data initially selected at the DataPkg level (Section 3.1), and the green region shows the fraction after the TOD-level selection that is used for mapping (Section 3.2).The time spans for Era 1 and Era 2 are shown in orange at the top; three periods with different VPM grids are shown in blue.The gray hatched regions indicate times when different parts of the forebaffle are blackened.Critical changes to the instrument configuration are marked by the vertical bars, with the text pointing at the direction where the annotation applies.The closeout was first taken off on 2021-09-27 but was occasionally reinstalled to guard against bad weather.
Source Avoidance -Although the scanning schedule had a 20 • boresight avoidance from the Sun, no such strategy was applied for the Moon and planets.With the 1.5 • beamwidth at 40 GHz, we flagged all detector data when the Moon is within 3 • of any detector's pointing direction to prevent the impact of the moonlight through detector crosstalk.Despite the Sun avoidance incorporated into the telescope scan, we were motivated to remove additional data when the Sun was above the horizon due to pickup observed in all detectors at the −70 dB level (100 µK for the Sun) when the Sun encroached on the telescope boresight.The spurious signal in every detector appeared in the same place relative to the telescope boresight position, independent of the pointing of the individual detectors across the 20 • field of view.This suggests that the issue was not due to direct pickup of sunlight (as is the case with far-sidelobes;Xu et al. 2020), but another indirect systematic effect, such as a change in the VPM-synchronous signal (VSS;Harrington et al. 2021) due to the exposure of the VPM to sunlight.The spurious signal and the corresponding data cuts are shown in Figure5.In Era 1, we flagged all data when the boresight of the telescope was within 60 • of the Sun.With the redesigned Era-2 forebaffle extension and roof, we were able to decrease this zone of solar exclusion to 40 • .However, a fan-shaped solar keepout region extending to 60 • in the direction of the 90 GHz telescope forebaffle opening was still required.The lower plots of Figure5show the undetectable impact of the spurious signal when the polarization measurements of +45 • and −45 • oriented detectors are combined.In this case, the spurious signal modulated at around 10 Hz decreases by the same amount (100 µK) in both +45 • and −45 • detectors.For +45 • (−45 • ) detectors, this produces negative (positive) spurious "polarization" signals that cancel one another upon combination.Because our survey maps incorporate all of the detectors and cancel the spurious solar signal as in the bottom half of Figure5, this avoidance cut represents a conservative measure to ensure data quality.

Figure 5 .
Figure 5.The CLASS Sun-centered linear polarization maps and the Sun avoidance cut.The maps are made by coadding all +45 • oriented detectors (top row) and all detectors (bottom row) in each era in the telescope boresight coordinates (independent of the boresight rotation) and are presented under orthographic projection.Some diffuse features around the compact structure and with an opposite sign in amplitude are artifacts due to the baseline adjustment in the mapmaking.The location of the lobe on the right in Era 2 corresponds to the location of the opening in the cage of the 90 GHz telescope.The features in Sun-centered maps resemble temperatureto-polarization (T -to-P) leakage that has a canceling effect between +45 • and −45 • oriented detectors.When all detectors are combined, the residual effect is subdominant to the noise.The inner black circles show the commanded 20 • Sun avoidance during the survey; the outer black circles/wedges are the extended boundaries for data flagging based on the relative position of the Sun with respect to the telescope boresight (60 • for Era 1 and 40 • − 60 • for Era 2) when the Sun is above the horizon.The position of the horizon changes with the boresight rotation in the telescope boresight coordinates, making regions below the dashed line also accessible from certain boresight angles.

Figure 6 .
Figure6.Demonstration of pre-processing and demodulation of a single detector's data over a 150 s window.Top: the raw data calibrated into power units (black) and processed data with the MCE readout Butterworth filtering (blue) and detector time constant (orange) deconvolved.The inset plot zooms into a shorter ∼ 0.45 s window (in red) to see the impact over a few VPM modulation cycles.The dominant 10 Hz component is due to the VSS.Middle: the deconvolved data were further bandpass filtered around the modulation harmonics lines (black).The VPM transfer functions M u/v are shown as green and purple dashed lines, with the amplitudes set to the best-fit values from the data segment, i.e., the demodulated u/v amplitudes.Here the calibration factor ∆TCMB/∆Di defined in Equation 6 is divided out to keep the demodulated data in power units.Bottom: demodulated u/v time streams.The corresponding range to the zoomed-in region is marked in red.The azimuth angle scanned over during this time period is marked on the top axis.The flatness of the demodulated curves indicates the stability of the polarization measurement.

Figure 7 .
Figure 7. CLASS 40 GHz focal plane layout and the VPM wire grid angle ϕP as seen by each detector.Each circle represents a pair of orthogonal detectors and is plotted by its pointing offsets on the sky.The orientation of the line within each circle denotes the VPM wire direction as seen by that detector, ϕP(×10), from the geometric modeling; the polarization sensitivity is 45 • from the wire direction.ϕP varies by 1.6 • across the focal plane.The detectors on the bottom of the plot, (shown in green) see the main polarization calibrator Tau A at the nominal 45 • elevation scan at all boresight rotations; the blue detector pairs on either side only have partial Tau A coverage depending on the boresight rotation angle, and the unfilled pairs have no Tau A coverage.

Figure 8 .
Figure 8.The wind-induced signal across the focal plane when the plastic closeout was on and the forebaffle extension was blackened.Each panel shows the demodulated linear polarization signal u binned in the wind bearing angle and wind speed coordinates for the −45 • oriented detector in each feedhorn; the other detector in the pair shows a similar signal.As indicated by the legend, the azimuth angle of the polar plots corresponds to the wind bearing angle, and the radial axis marks the wind speed in units of meters per second.The panels are laid out by the detector pointing offsets, similar to Figure7, with the detectors pointing at higher elevations at the top.Within each panel, the wind-related systematics are mostly confined in the quadrant when the telescope is facing toward the wind.The signal amplitude varies across the focal plane with a quadrupole pattern that increases from the center toward the edge of the focal plane and also rises with the wind speed.The whirlpool-shaped structure in some of the edge detectors is due to the ground pickup(Section  4.1.3)

]Figure 9 .
Figure 9. Wind direction distribution at the Atacama site as measured by the CLASS weather station over the 6 yr survey.The wind direction distribution is shown as the histogram in the outer and inner ring for the Austral winter (June to August) and summer (December to February), respectively.The histograms are color-coded by the mean wind speed in each bin.

Figure 10 .
Figure 10.Ground pickup signals and the landscape around the CLASS site.The concentric annuli show the demodulated linear polarization signal (u) for a pair of detectors at the bottom edge of the focal plane (linear polarizations after demodulation have the same sign between the pair and are averaged to enhance the signal) binned in the telescope azimuth angles for the seven boresight rotations.Left: data taken prior to 2017-06 when the circular baffle was not blackened.Of note is the Cerro Toco mountain peak (15 • elevation) at around 45 • in azimuth, which aligns with the peak of the u signal.For this edge detector pair, the signal from the mountain appears earlier or later in the telescope azimuth coordinates depending on the boresight rotations.Right: data from 2017-07 to 2018-01, where the interior of the baffle extension was blackened, and before the replacement with the asymmetric double baffle extension and new baffle roof.The Era 2 configuration gave a reduction in ground pickup comparable to the right panel in all but the outermost detectors, regardless of the blackening state.

Figure 11 .
Figure11.The effect of template iteration in correcting the large angular scale power bias.The top panels show the signal power spectra for EE and BB. Green curves are the input synchrotron power spectrum model from PlanckCollaboration et al. (2020c), and the associated sample variance for a 75% sky coverage.The colored curves show the estimated signal spectra, averaged over 2000 simulations, at each template iteration.The error bars are the standard error of the mean and are only shown for the last iteration for visual clarity.The bottom panels show the residual spectra between the output spectra and the input spectra normalized by the input amplitude, in the same color scheme.The correction converged quickly at high ℓ where the corresponding time streams are noise dominated; at low-ℓ the iteration converged by the fifth iteration, with a bias below 2.5% for ℓ ≤ 5, much smaller compared to the sample variance.

Figure 12 .
Figure 12.CLASS 40 GHz polarization maps in equatorial coordinates under Mollweide projection.The linear (Stokes Q/U) and circular polarization (Stokes V ) are the final products of the data pipeline.The maps are smoothed to 2 • resolution to enhance the large-scale features.The gray-shaded regions are not surveyed.The color scale is linear below 5 µK and logarithmic above to show the structure in the map where bright synchrotron radiation dominates.Due to the designed VPM throw, the noise levels are similar in the Q, U, and V maps; therefore, the fluctuations in the V map approximately represent the noise in the Q/U maps.The apparent signals in the Q/U maps are explored in E23.

Figure 13 .Figure 14 .
Figure13.CLASS 40 GHz hits maps and cross-linking map.The upper triangle of the plot shows the hits maps in units of the integration time per Nside = 128 pixel (0.46 deg 2 ).This is proportional to the inverse-variance map assuming constant white noise in the raw data around the VPM modulation frequencies.The bottom-left corner shows the cross-linking map in arbitrary units.Higher values reflect even coverage of the scanning direction and therefore suppression of the scanning-related low-frequency noise.

Figure 15 .
Figure15.Summary of the effect of multiple systematic errors on the power spectra.The black curve shows the map noise estimated from an ensemble of simulations with the mapping transfer matrix corrected.For reference, the signal spectra are plotted as dashed curves for the diffuse synchrotron signal at 40 GHz(Planck Collaboration et al. 2020c, dark green )  and the CMB(Planck Collaboration et al. 2020b, light green , a range of tensor-to-scalar ratio 0 < r < 0.01 is represented by the shaded green area in the BB panel), respectively.The signal spectra are convolved with the 1.5 • FWHM beam window function.The solid/dotted curves are the measured positive/negative systematic bias; the curves with downward arrows represent the 2σ confidence level upper limit for systematics that are not detected given the sample variance in the simulations.Pink : systematic error from the polarization angle calibration uncertainty (Section 5.1.1).Orange : effect of a 0.1 ms bias in the detector time constants (Section 5.1.2).Red : effect of a 3 × 10 −3 level beam ghosting across the focal plane (Section 5.1.3).Light blue : effect of the dipolar T -to-P leakage (Section 5.1.4).Navy blue : uncertainty in the VPM transfer function parameters.The curves indicate the maximum variation in the residual power spectra for VPM parameters drawn from the 2σ confidence interval of the VPM parameter optimization process (Section 5.1.5).The shaded regions in EE and VV highlight that these are the result of an ensemble of parameters, but are not quantitative depictions of the spread.Brown : effect of a +0.3 bias in the linear polarization spectral index assumed for demodulation (Section 5.1.5).Purple : residual E-to-B mixing due to the mapping filters after the transfer matrix correction (Section 5.1.6).

Figure 16 .
Figure 16.Effect of dipolar T -to-P leakage.Top: the demodulated linear polarization (Stokes u) Moon maps in the VPM coordinates for a pair of detectors.In this coordinate system, the VPM wire grid is horizontal.The colors are scaled to the peak amplitude of the Moon temperature maps.The 1.54 • FWHM main beam is marked by the white dashed circle.Bottom left: the differential linear polarization maps of Tau A between all +45 • /−45 • detectors.The maps are rotated by the Tau A polarization angle at 40 GHz (Weiland et al. 2011) so that any residual Tau-A polarization signal would appear entirely as a point source in the Stokes-Q map (denoted Qrot).Because the +45 • /−45 • signals are differenced, their opposite-signed dipoles (top row) average constructively.Bottom right: simulation of the dipolar leakage effect from convolving scaled WMAP Q-band temperature maps with the dipolar leakage beam from all detectors.

Table 1 .
Data Selection and Processing Flags