CLASS Data Pipeline and Maps for 40 GHz Observations through 2022

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), and Planck (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;Keck/Keck Collaboration et al. 2022;Spider Collaboration et al. 2022), the gravitational lensing potential (Polarbear Collaboration et al. 2014;Keck2 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;Keck2 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 Keck/ 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 variabledelay 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 ground-based surveys target scales ℓ  30 (Simons Observatory Collaboration et al. 2019;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, which includes updated calibration that impacts some results presented here). 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.

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. 15 Given its proximity to the equator, the site also provides access to most of the sky, which is essential for large-angular-scale 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 (2018-present; 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 center of the telescope's field of view and has a full range of motion from -45°-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 February 22 as Era 1 and after 2018 June 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 Bluefors 16 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 feedhorn-coupled 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. On-chip orthomode 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 feedhorn-coupled 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 time-division 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 sensitivity 17 was m 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 m 30.5 K s in Era 2.

Observations

This paper covers 40 GHz observations from their beginning on 2016 August 31 through 2022 May 19 (2089 days; 5.72 yr). 16 bluefors.com 17 The referenced numbers correspond to the noise-equivalent temperature of the raw data in the white noise regime. With the VPM modulation, the polarization sensitivities are degraded by a factor corresponding to the (inverse of) modulation efficiency, which is around 0.72/0.48 for linear/circular polarizations, respectively (see Equation (15)).

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 radio-frequency 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°a nd 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 of the survey are labeled Era 1 (before 2018 February 22) and Era 2 (after 2018 June 22). Important changes are noted in Figure 2 and explained below in approximate time order.

Cage Camera Interference-On 2017 August 9, 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 June 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 October 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 February 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 July 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 January 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 April 1 to 2017 May 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-thanexpected 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 noiseequivalent 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 January 12. This canceled the efficiency gains from the increased in-band transmission of the new infrared filters (Dahal et al. 2022;J.C. Cleary et al. 2023, in preparation).

Azimuth Scan Speed Increase-The az-scan speed was increased from 1-2 deg s -1 on 2019 March 4. This was found to decrease the low-frequency noise in the demodulated data at 40 GHz (J.C. Cleary et al. 2023, in preparation). 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 January 8, 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 February 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 September 27, the closeout has been removed for most CMB observations and has only been temporarily reinstalled during bad weather.

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 DataPkgs collected when the mount elevation was above 40°w ere kept for processing. As a final step, DataPkgs were discarded if the cloud cover as monitored by the site webcams (Y.L. Li et al. 2023, in preparation) 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).

TOD Selection

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. 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.

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).

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 Figure 5. 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 keep-out region extending to 60°i n the direction of the 90 GHz telescope forebaffle opening was still required. The lower plots of Figure 5 show 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°a nd -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 Figure 5, this avoidance cut represents a conservative measure to ensure data quality.

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 ratio, S/N, below five).

Conditioning Cuts-The final mapmaking operates on 10-sweep 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.

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 advantage 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.

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:foot_3

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 f P such that:

sin 2 0 sin 2 cos 2 0 0 0 1 , 2

where γ is the detector position angle on the sky, and f 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 beam-filling 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

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/v 20 Figure 7. CLASS 40 GHz focal plane layout and the VPM wire grid angle f 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, f P ( × 10), from the geometric modeling; the polarization sensitivity is 45°f rom the wire direction. f 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. through ( )

x e e ,s 1

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 isfoot_5

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 leastsquares solution to Equation ( 12) is (see also, Harrington et al. 2021)

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 lowand 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 gap-filling 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-)azimuth-synchronous 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 minkasifoot_6 (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;  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.

Wind-induced Signal

A quasi-azimuth-synchronous signal is present in the demodulated linear polarization data that correlates with the wind recorded by the WeatherHawkfoot_8 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 wind 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 2021 September by removing the covering plastic during observations. For the time period with the closeout on, the filter components took the form:

where γ w and f w are the wind speed and wind direction measured by the weather station, f(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. 25 The 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 az-synchronous 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 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 timescales 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 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.

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.

Noise Model

Like other ground-based CMB experiments, CLASS (demodulated) data are noise dominated; therefore, the noise term in Equation ( 16) needs 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 timescales ( 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 lowfrequency (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 block-diagonal 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 of 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 azimuthsynchronous 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 maskaware data likelihood as

t  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 gap-filled 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 ( 16) given the noise model N from Sections 4.2 and 4.3 is

This was solved iteratively using the preconditioned-conjugate gradient method (Shewchuk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 ( 25) is 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 crosscorrelation 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 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 large-and small-scale noise (Harrington et al. 2021;J.C. Cleary et al. 2023 in preparation); 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-pixel-resolution 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. 26 The 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 jth 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) largescale 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 m 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 component 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 (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. resolution for this situation:

where m in/out are the input and filtered maps downgraded to N side = 16 resolution, F ij 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.

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 28 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 f 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 f 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 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°; Jarosik et al. 2007;Weiland et al. 2011) and Planck (-88.65 ± 0.79 ± 0.50°; Planck Collaboration et al. 2016) at similar frequencies. At around 40 GHz, the CLASS beam at FWHM 1.54°is wider than that of WMAP (0.49°) and Planck (0.47°), and the resolution confusion can contribute a subdegree discrepancy in the polarization angle determination. Since no significant discrepancy in the polarization angle is found within internal comparison or externally with other experiments (Aumont et al. 2020), we do not apply any correction to the detector angle.

The E-B mixing from a 0°. 7 polarization angle error is shown in Figure 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 (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 backward 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.

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°f rom 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 (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°d etectors. Because of polarization modulation, the linear polarization can be recovered with the +45°and -45°d etectors 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 bottom-right panel 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°l eakage 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 crosscorrelation with the CMB E mode due to the CMB TE 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.4 that 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 realworld 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.

Conclusion

We have presented a detailed description of the CLASS data reduction pipeline for 40 GHz observations conducted from 2016 to 2022 May. When weather, instrument upgrades, and other interruptions permitted, observations were conducted continuously, regardless of time-of-day or seasonof-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 decl.

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 sidebands 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 largeangular-scale measurement. The demodulated polarization data were then filtered to remove systematic effects, such as azimuthsynchronous and wind-induced 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-conjugategradient-descent algorithm to jointly solve for the maximumlikelihood 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 [ ( ) ] m + ℓ ℓ 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 Λ cold dark matter (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 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 data set. 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. 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 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 corresponds 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:

( ) ( ) ( ) ( ) ( ) a b D = + az el az el az sin cos sin sin A1 ( ) ( ) ( ) b a D = el az az cos sin , A2

where α represents a tilt of the mount to the west, and β represents a tilt to the north. Here, az and el are the 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 ( ) el sec 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.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 scanned back and forth at a constant elevation while the Moon rises or sets through their fields of view. This analysis yields the position of each detector on the sky at the point when the Moon is at the beam center along with the mount position at that time and an average boresight angle offset from the boresight angle determined through a minimization of the detector positions with respect to their positions as given by the current array pointing model for each telescope projected onto the sky. These data are then used in an iterative nonlinear least-squares fitting procedure to determine the pointing model coefficients for each telescope. The residual rms in azimuth, elevation, and boresight angle is typically less than 1′ for a given set of Moon scans. The standard errors of the individual detector offsets from array center are typically less than 10″.