Introduction

Recent advances in cosmology from cosmic microwave background (CMB) observations greatly benefit from the use of ground-based telescopes (Aiola et al. 2020;Dutcher et al. 2021;BICEP/Keck Collaboration et al. 2022) that complement space missions with high sensitivity polarization and small-scale intensity measurements. CMB observations from the ground are limited by the atmosphere in many ways. Strong emission lines, mostly from oxygen and water vapor, define the transmission windows. The turbulent and bulk motion of the atmosphere causes fluctuations in its brightness temperature on large time/angular scales, especially for inhomogeneously mixed species like water vapor, which leads to long-term instabilities of CMB measurements (Errard et al. 2015;Morris et al. 2022). Polarization measurements taking advantage of rapid modulation (Johnson et al. 2007;Chuss et al. 2012) are less sensitive to the unpolarized component of atmospheric turbulence, but the fluctuations could still degrade polarization sensitivity through instrument polarization.

Despite being mostly unpolarized (Battistelli et al. 2012), the atmosphere is known to be circularly polarized due to the Zeeman effect from the magnetic dipole transition of the oxygen molecules in the Earth's magnetic field (Lenoir 1967(Lenoir , 1968;;Petroff et al. 2020). This is a systematic issue for experiments with polarization modulation (Chuss et al. 2012;Nagy et al. 2017) and presents a challenge to ground-based pursuits of astrophysical/cosmological circular polarization signatures (King & Lubin 2016;Padilla et al. 2020;Eimer et al. 2023).

At millimeter wavelengths, linear polarization can be generated from unpolarized radiation via scattering by a hydrometeor in its liquid or ice phase. Polarization from the cloudy atmosphere was first detected at 100 GHz by Troitsky & Osharin (2000) and modeled as the scattering of the upwelling ground emission and atmospheric emission by the cloud particles. Similar observations have been made across millimeter wavelengths from both ground (Troitsky et al. 2003(Troitsky et al. , 2005) ) and remote-sensing satellites (Chepfer 1999;Davis et al. 2005;Prigent et al. 2005;Defer et al. 2014;Gong & Wu 2017). The cloud polarization from groundbased observations is mostly horizontally aligned (i.e., negative Stokes Q parameter. See Figure 1 for a schematic illustration) albeit with some exceptions of +Q, which are attributed to spherical cloud particles (Troitsky et al. 2003) or mixed-phase cloud (e.g., melting ice crystals, Troitsky & Osharin 2000;Troitsky et al. 2005). For nonspherical cloud particles (e.g., ice crystals), several mechanisms are at play to align the particle orientations and thus enhance polarization. Vonnegut (1965) points out that ice crystals could be vertically aligned by a discharge event in the cloud (Hendry & McCormick 1976;Prigent et al. 2005). More commonly, ice crystals are found to be horizontally aligned (Magono 1953;Ono 1969) since a hydrodynamic force acts on ice particles as they settle, tending to orient them with their broadside to the fall direction to maintain maximum air resistance (Klett 1995). Due to the counter-effect from turbulence (Bréon & Dubrulle 2004;Hashino et al. 2014;Zeng et al. 2023), the concentration of horizontal alignment is often low, and the alignment is most frequently found in warm clouds ( -20°C), with a typical occurrence of about 20%-80% depending on the temperature and the ice crystal habits (Noel & Chepfer 2010;Westbrook et al. 2010;Stillwell et al. 2019).

The impact on CMB observations from polarized cloud emission is found by Pietranera et al. (2007) to be twofold, namely, the backward scattering of the ground thermal emission and the forward scattering of the sky signal. While the latter distorts the CMB polarization signal, the former is a more substantial issue due to the large difference between the ground temperature (which sources the polarized cloud reflection) and the CMB polarization amplitude. The POLARBEAR team (Takakura et al. 2019) made the first connection between the cloud events seen in camera images during the daytime and burst signals in their 150 GHz detectors. The polarization of the bursts is consistent with the scattering of cirrus clouds composed of horizontally aligned ice crystals.

The Cosmology Large Angular Scale Surveyor (CLASS; Essinger-Hileman et al. 2014) is a CMB polarization experiment located on Cerro Toco in the Atacama Desert that aims to measure polarization on large angular scales from the ground. Here we provide an assessment of the cloud signal in multifrequency CLASS observations at 40, 90, 150, and 220 GHz between 2016 August and 2022 May. Table 1 summarizes the effective frequencies of the sources (assuming a Rayleigh scattering spectrum T RJ ∝ ν 4 , Dahal et al. 2022) and the start date of each detector array. The polarization modulation technique employed by CLASS enables both linear and circular polarization sensitivity. The fast surveying of the sky (scan through 36% of the half-sky twice every 10 minutes), accompanied by all-day optical image monitoring, offers a holistic view of the cloud conditions on Cerro Toco.

This paper begins with a review of the radiative properties of clouds in Section 2, and describes optical cloud identification from camera images in Section 3 and microwave polarization observations in Section 4. Section 5 summarizes the population analyses on cloud polarization, and we conclude in Section 6.

Cloud Physics

Cloud Characterization

Altitude is the main factor in determining the cloud composition and its classification. At 5200 m on Cerro Toco, the predominant types are high-level clouds (above 6 km in the tropical region, e.g., Cirrus, Cirrostratus, Cirrocumulus, subclassified based on the morphology) that are almost entirely composed of ice crystals and mid-level clouds (up to 8 km, e.g., Altostratus, Altocumulus) that may contain both ice and (supercooled) water droplets.

The water content in the atmosphere is often quantified by the precipitable water vapor (PWV) and liquid/ice water path (L/IWP), which measure the total mass of the water per unit The two orthogonal linear polarization states are labeled; they correspond to the Q (0°, in orange) and -Q (90°, i.e., horizontal, in blue) Stokes parameter evaluated in the horizontal frame of the telescope. Clouds are depicted as cumulus clouds instead of cirrus for their aesthetic appeal. Objects are not to scale.

Table 1 CLASS Multifrequency Information Band 40 90 150 220 ν eff (GHz) 38.5 94.2 148.2 220.0 Start date 2016 Aug 2018 Jun 2019 Oct 2019 Oct Note. The numbers are calculated assuming the fiducial Rayleigh scattering spectrum.

area in liquid/solid phases. As one of the exceptionally arid sites for millimeter astronomy, the Chajnantor area experiences a median of 1 mm PWV (Cortés et al. 2020). Water in the hydrometeor form is long tail distributed, with LWP and IWP below 10 -3 kg m -2 80% of the time but could rise above 10 -2 kg m -2 for 10% of the time (Kuo 2017, austral summer months excluded). This is also consistent with the typical 10 -2 kg m -2 IWP found for cirrus clouds (Davis et al. 2007;Austin et al. 2009).

Figure 2 shows the vertical profile of the cloud fraction and liquid and ice water content within the atmosphere above the CLASS site. The data near the CLASS site, within a 0°. 25 × 0°. 25 pixel centered at 67°. 8 W, 23°. 0 S, were queried from the global reanalyses ERA5 (Hersbach et al. 2020). The annual averaged cloud profile shows a bimodal distribution; mid-level clouds with liquid water constituents are 2 km above the ground (7 km above sea level), and ice clouds are higher in altitude, at the boundary of the troposphere. We note that the water content and cloud fraction distributions are highly skewed and that the records are zero most of the time.

Rayleigh Scattering

In the Rayleigh limit where the scatterer size is small compared to the wavelength, the electric field component of the scattered radiation E o at direction r from the scatterer followsfoot_1

where V is the scatterer volume; ν is the frequency of the incoming radiation E i ; and c is the speed of light. The polarizability tensor, α, depends on the dielectric properties and shape of the particle. Water droplets are often spherical, whereas ice crystals have a variety of habits ranging from regular shapes including quasi-spheroids, elongated columns, thin plates, and branched dendrites (Figure 1) depending on the temperature and vapor supersaturation (Libbrecht 2017) to highly asymmetric aggregates (Lawson et al. 2019). For simplicity, we model all particle types as spheroidal; following the notations in Takakura et al. (2019), the polarizability tensor is parameterized as

where ò is the relative permittivity and

z is the depolarization matrix for spheroidal scatterers (with Δ z less or greater than 1/3 for prolate or oblate crystals, respectively).

For cloud reflection, we define the coordinates centered on the scatterer and evaluate the radiation in terms of Stokes parameters. The incoming radiation is assumed to be unpolarized: S i = [T(θ), 0, 0, 0] T , with the intensity profile of

where θ is the zenith angle; q » h R 2 h e is the astronomical horizon for clouds at altitude h above the ground (R e is the radius of Earth), T g is the ground temperature, and T a (θ, ν) is the brightness temperature profile of the atmosphere. This setup is illustrated in Figure 1. For spherical particles with α ∝ (ò -1)/(ò + 2)I, the Stokes parameters from cloud reflection evaluated in the telescope frame can be expressed as

where S o (ν) is the radiation received by the telescope at (-δ, 0), with δ being the telescope pointing elevation; matrix R rotates Stokes vectors into the telescope frame from the scattering frame in which the phase matrix M is evaluated. The column of M for unpolarized incoming radiation is

, where Θ(θ, f, δ) is the scattering angle. The Rayleigh scattering optical depth is In both panels, the solid curves show the mean value at given pressure levels; the shaded regions fill between the 50th and 95th percentiles of the distribution. Due to aridness at the queried region, the 50th percentiles are essentially zero. The vertical dotted line marks the altitude where the median air temperature is -40°C, above which (altitude) homogeneous nucleation occurs.

2 10 100 GHz 100 m IWP 10 g m , 6 5 4 3 2

where σ sca is the Rayleigh scattering cross section; n(r) is the number density of cloud particles with diameter D. The typical crystal size in cirrus is 100 μm, although the distribution is widespread (Austin et al. 2009;Lawson et al. 2019). At millimeter wavelengths, the dielectric permittivity for ice is, ò ; 3.2 (Warren & Brandt 2008), and being primarily of interest in this setting (see Figure 3) is adopted in evaluating Equation (6). For liquid water, the permittivity has a larger magnitude and can be approximated by a Debye dielectric function (Ellison 2007;Turner et al. 2016, also see the discussion in Section 5.3).

Figure 3 shows the total linear polarization (Stokes Q, since U polarization in the horizontal frame is zero under parity symmetry) and the polarization fraction, observed at δ = 45°e levation, from horizontally aligned spheroidal cloud particles with random azimuthal orientations as a function of the depolarization factor Δ z (also converted to the aspect ratio of the spheroids assuming uniform permittivity). The absolute polarized emission is evaluated assuming an optical depth corresponding to a 10 g m -2 LWP and a particle size of 100 μm (geometric mean diameter). The atmospheric brightness profile T a (θ, ν) is computed with am (Paine 2022), using the annually averaged atmosphere profile at the Chajnantor site compiled from the reanalysis of MERRA-2 (Molod et al. 2015) and the CLASS bandpass (Dahal et al. 2022). The linear polarization fraction is highest for nonspherical scatterers with extreme aspect ratios and peaks at ∼30% for oblate crystals at 45°e levation. For spherical particles, the Q polarization could be positive due to the thermal emission from the atmosphere, which is more prominent at higher frequencies, and for clouds at lower altitudes where the air temperature and water vapor above the cloud are higher. This is the main modeling difference from Takakura et al. (2019), which did not consider the brightness temperature contribution from the atmosphere. Ignoring the forward scattering of atmospheric emission in Equation (3) and solving for the linear polarization component with Equation (4), we find negative Q polarization for all particle shapes, with the lowest amplitude from spherical scatterers:

and the polarization fraction is 1% for δ = 45°. This agrees with the derivation in Takakura et al. (2019).

The derivation here only considers single scattering and a full radiative transfer calculation is needed to self-consistently account for the Rayleigh scattering and contributions from absorption/emission of the atmosphere; nevertheless, the result is consistent with the full radiative transfer model (Troitsky & Osharin 2000).

Mie Scattering

At wavelengths comparable to the size of the cloud particles, the Rayleigh scattering approximation breaks down and the full treatment using the Mie theory is needed. The Mie solution differs from the Rayleigh approximation in three major aspects:

(1) the resonance effect related to the size and refractive index of the scatterer is considered, which yields complex frequency dependency and is generally shallower than the Rayleigh scattering spectrum ∝ν 4 ; (2) the forward scattering is enhanced compared to other directions; (3) the scattering phase matrix contains nonzero linear-to-circular polarization coupling terms.

The first two points alter the cloud reflection signal. As a simplified model, we consider the cloud composition as homogeneous spherical particles and ignore the atmospheric emission in Equation (3). The Rayleigh solution for polarization is given in Equation (7), and we compute the corresponding result for Mie scattering using miepython. 16 The results are shown in Figure 4 as the ratio of the linear polarization spectrum from the Mie solution and the Rayleigh scattering approximation. The polarization signal at around 220 GHz is suppressed in the Mie solution in a way that depends sensitively on particle size in the 300-400 μm range. While the majority of ice crystals are below this threshold, the size distribution is long-tailed (McFarquhar & Heymsfield 1997;Figure 3. Polarized signal from Rayleigh scattering of thermal emission from the ground and the atmosphere for different cloud particle shapes. Assuming homogeneous spheroidal cloud particles with uniform dielectric properties, the particle shape is parameterized by its aspect ratio, or equivalently, the depolarization factor Δ z . The linear polarization from an ensemble of randomly oriented cloud particles with the short axis aligned vertically is computed for each of the four CLASS frequency bands (color coded), for different cloud heights above the ground (5200 m above sea level), with a total 2 mm PWV. The cloud height and PWV determine the optical loading above the cloud (Equation (3)), and impact the polarization from the forward scattering of atmospheric emission, especially at higher frequencies. Top: the Stokes Q polarization signal in units of brightness temperature for a reference optical depth τ ν equivalent to IWP = 10 g m -2 and particle diameter D = 100 μm (geometric mean of the three axes). Bottom: the corresponding polarization fraction. Signed values are used to distinguish the sign of Stokes Q polarizations. Heymsfield et al. 2002;Austin et al. 2009), and the larger particles can contribute significantly to the polarization signal (roughly ∝D 6 ). Assuming the size distribution of ice crystals with log-normal mean size 46 μm (corresponding to the best-fit parameters at -70°C, Austin et al. 2009) and a standard deviation of 0.226, the black curve shows the populationaveraged spectrum ratio. While Rayleigh scattering is a good approximation of cloud polarization at frequencies below 150 GHz, the Mie scattering effect cannot be ignored for higher frequencies. For realistic high cloud composition, the cloud polarization signal at the 220 GHz band is a factor of a few smaller compared to the Rayleigh solution and is comparable to the signal at 150 GHz.

In the Rayleigh scattering regime, circular polarization from directly scattering atmospheric circular emission (Petroff et al. 2020) is negligible due to parity. For Mie scattering, the coupling between linear and circular polarization (Bohren & Huffman 1983) can convert the linear polarization induced by preceding scatterings into circular polarization (Kawata 1978;Slonaker et al. 2005). However, the required multiple scattering is improbable given the low scattering optical depth.

Cloud Measure

Camera Images

CLASS has deployed several cameras at its observing site to monitor telescope operations. Some of these cameras have a large fraction of sky coverage, which can be used for cloud detection. These images are taken at a 10 minute cadence, approximately the time for the telescope to complete a 720°s can in azimuth. Figure 5 shows the positioning and pointing directions of the seven cameras at the site and the landscape around the CLASS telescopes. The three primary cameras (labeled "1-3") are functional during daytime (colored) and nighttime (grayscale), whereas the rest only have daytime data due to the limited shutter speed at night, and were installed in the middle of the survey. For each stable period of each camera pointing, we built the World Coordinate System (WCS) for the images. For cameras 1-3, this is accomplished by solving the star positions with astrometry.net (Lang et al. 2010), whereas for daytime-only cameras, the camera models were solved from the terrain shapes, with geographical information from Google Earth. The fields of view (FoVs) of the cameras are shown as black wedges (solid and dashed) in Figure 5, which have a minimal intersection with the sky swept by the telescopes (green shades). Therefore, the connections we make in the following analyses have to assume similar cloud conditions seen by the telescopes and the cameras.

The cloud cover in each of the site images was estimated using a convolutional neural network (CNN) classifier. The details of the implementation and the validation of the classifier are described in Appendix. The CNN classifier labels each of the 240 × 240 pixel subregions in an image as one of the three classes: cloud, clear, and obscured (by the ground, the telescope, or celestial objects that prohibit a definitive cloud identification), and the cloud measure is defined as the ratio of subregions classified as cloud versus the total number that is not obscured (cloud + clear). Combined with the WCS of each site image, the cloud measure can be reported with finite angular resolution. In the rest of this work, we use the spatial average (over the sky covered by all cameras) of the cloud measure above 10°elevation (CM 10 ) as a proxy of the cloud condition at the site.

Cloud Conditions at Cerro Toco

The high-cadence all-day monitoring during the 6 yr survey of CLASS provides an assessment of the cloud conditions at the site on Cerro Toco. Figure 6 shows the CM 10 statistics from 2016-2022. Over the course of a day, the cloudiness begins increasing after midday and peaks around dusk; this agrees with the diurnal cycle of cirrus and cumulonimbus clouds over land (Eastman & Warren 2014;Feofilov & Stubenrauch 2019), which develop during the afternoon due to the mixing of the boundary layer and the increase in thermal convection. This variation is steeper at summer times when the convection is most efficient and peaks around January, consistent with the record from POLARBEAR (Takakura et al. 2019) at the same site. The annual variation is also evident and tracks well with the PWV measurement from the APEX weather station. 17 The apparent yearly modulation of the cloudiness around twilight reflects the synchronous shift of the diurnal peak of cloud cover with solar motion; however, part of this is the systematic uncertainty due to the higher false-positive cloud detection caused by aerosol scattering around the Sun.

The average cloud cover at Cerro Toco is around 32% during daytime and 20% during nighttime, modulated by a seasonal variation of around 50%, consistent with the cloud detection rate of daytime images taken from the neighboring POLAR-BEAR site (Takakura et al. 2019).

Cloud Signal in Microwave

CLASS Instruments and Polarimetry

CLASS observes at a constant 45°elevation and scans continuously in azimuth by 720°before turning around at a rate of 1°-2°s -1 . Polarization is measured through the variabledelay polarization modulator (VPM; Chuss et al. 2012;Harrington et al. 2018Harrington et al. , 2021)), which modulates between the linear polarization (Stokes U in the telescope coordinates) and circular polarization at around 10 Hz. For cosmological observations, the linear polarization angle coverage is achieved by the rotation of Earth and the telescope boresight rotation from -45°to 45°at 15°increments every day. For a cloud transient event, however, the VPM is only sensitive to one polarization direction-e.g., ± Stokes Q for m45°boresight and Stokes ±U for 0°boresight. Nevertheless, polarimetry is still possible since the ∼30°azimuthal differential pointing of the detectors on the focal plane translates into a ∼15°spread in position angle on the sky at 45°elevation. This finite position angle coverage permits the determination of the polarization angle from a single scan, albeit with lower precision compared to that from rapid modulation between the two linear polarization states through a continuously rotating half-wave plate (e.g., Takakura et al. 2019).

Data Selection

The radiometer data analyzed in this work are taken from the CLASS survey from 2016 August to 2022 May. As the 90 and 150/220 GHz receivers were deployed for less time, their data covered a smaller date range. The beginning date of each frequency is summarized in Table 1. For the purpose of investigating the clouds, the data selection is less stringent compared to that for cosmological analysis. The data were chunked into spans defined by a segment of uninterrupted observations. Typically, a span is a day long and is only separated by the detector calibration (Appel et al. 2022) and boresight rotation change at midday. Similar to the data selection described in Li et al. (2023), data with instrument failure (VPM, cryostats, mount, readout system) or abnormal detector response (constant detector output or excessive flux jumps) were discarded. Data were masked when the telescope scans close to the Sun or the Moon. No other source avoidance was made. Other environmental conditions such as the high PWV, and high cloud fraction were not used as selection criteria to avoid rejecting the cloud signals.

Azimuthal Maps

To find cloud signals among other burst-like events in the data, and to obtain better polarization angle measurements, we made polarization maps combining multiple detectors from the focal plane for each of the Stokes parameters.

The linear (u) and circular (v) polarization signals were obtained by demodulating the time stream with the VPM transfer function (Harrington et al. 2021;Li et al. 2023):

where D is the raw time stream modulated by the VPM, M is the VPM modulation transfer function, and L is the low-pass filter defined by the signal bandwidth given the beam size and scanning speed. As the VPM modulates, its synchronous signal is also modulated at around 10 Hz and so contributes to the single-detector demodulated linear polarization as a baseline offset. 18 To minimize its impact on the single scan data, we leveraged its asymmetry between the orthogonal detector pairs and used the pair-add demodulated data (equivalent to pairdifferencing for the unmodulated experiment, Harrington et al. (2021), J. Cleary et al. 2023, in preparation) for the rest of the analysis. We keep the notations u/v for the pair-add data for simplicity. The demodulation process depends on the modeling of the VPM transfer function, and its error could result in a bias in the polarization amplitude or cause leakage between linear and circular polarization. One of the main factors in the model is the assumed sky spectrum (the cloud linear polarization and atmospheric circular polarization) due to the frequency dependence the VPM modulation efficiency within the detector bandpass. CLASS data were demodulated assuming a sky spectrum consisting of CMB and Galactic foreground contribution in linear polarization and atmospheric emission in circular polarization (Petroff et al. 2020). The assumed linear polarization spectrum is shallower than the fiducial spectrum from Rayleigh scattering from clouds, which inevitably leads to polarization bias. We denote the bias to linear polarization and linear-to-circular polarization leakage from VPM transfer function error by B uu and B vu , respectively. The modeling of these quantities requires accurate calibration of the VPM parameters; however, in the absence of such calibration, we elect to empirically determine these bias factors from the data in Section 4.5.

The demodulated data can be related to the sky polarization signal 19

where γ is the detector position angle on the sky and f P is a small angle correction accounting for the VPM wire direction as viewed by each detector (Li et al. 2023). The sky maps were solved by inverse-variance weighting the demodulated data from all detectors:

where P is the pointing matrix in Equation (9) and N is the diagonal noise matrix. To guard against the bright cloud signal biasing the variance estimates, each detector time stream was first outlier-removed using an iterative median absolute deviation (MAD) algorithm where for each iteration samples with absolute deviation from the median greater than five times the median of the absolute deviation from the median were flagged as outliers. A singular-value decomposition was then performed on the outlier-removed data to isolate the common modes. The per span variance of each detector was estimated after projecting out modes with singular values four times greater than the median. Prior to the map-making, the median of the outlier-removed data was subtracted from each detector to prevent the residual VPM synchronous emission or other long-timescale drifts from biasing the result.

Limited by the instantaneous detector coverage of the transient cloud events, we binned the data into azimuthal maps (hereafter az-maps), where close-by data in the time-azimuth coordinates are averaged over regardless of the elevation pointing. For all frequencies, the azimuth resolution was set to 1°, and temporal resolution was chosen such that each continuous 360°scan was binned into the same temporal bin. Depending on the scanning speed, the temporal resolution is around 180-360 s. The az-maps trade off spatial resolution in the elevation direction for polarization sensitivity; therefore, elevation structures of the cloud smaller than the focal plane size are lost. Examples of the multifrequency Stokes Q/U azmaps are shown in Figure 7.

Cloud Detection

The search for cloud signals in CLASS data aims to find any burst-like linear polarization signal that has temporal and azimuthal contiguity without assuming a priori knowledge of the cloud polarization properties. Instead of operating on the sky az-maps introduced in Equation (10) for cloud detection, we made separate az-maps for the detector frame u signal by mapping it as an intensity signal, i.e., disregarding the position angle dependency in Equation (10). After the map-making as outlined in Section 4.1, we filtered out from the u az-map a low-order polynomial component in the temporal direction with outlier rejection based on MAD. The polynomial order was chosen to be the integer number of hours of each span.

Depending on the boresight angle of the telescope, the -Q polarization from the cloud can appear as bursts in the u az-map with positive (positive boresight) or negative (negative boresight) signs. The search was conducted for both cases, where pixel values exceeding seven times the map white noise level (propagated from the u time stream variances) were flagged. Assuming the smoothness of the cloud signal, the boundaries of flagged pixels were further expanded and tailored using consecutive morphological closing and opening algorithms. These morphological operations take into account the periodic boundary condition in the azimuth direction. Each isolated contiguous flagged region was identified as a cloud candidate, and an additional morphological dilation operation was applied to expand the region for background estimations. The sky polarization signals were extracted from the sky az-maps for each Stokes parameter at the corresponding region with background subtraction. The variance in the background region was assigned as the measurement uncertainty.

As an example, the last two panels of Figure 7 show the 90 GHz u az-map and the cloud detection flags for the span between 2022 January 25 and 2022 January 26. For the population analysis below, the same extraction procedure was performed independently for each of the four CLASS bands for every span. In total, we have identified 1481, 3677, 4110, and 4178 cloud candidates at 40, 90, 150, and 220 GHz, respectively.

Circular Polarization and Demodulation Bias

Clouds are not expected to produce circular polarization, and the reflection mechanisms outlined in Section 2 should not yield significant circular polarization due to the low Rayleigh scattering optical depth. However, nonzero circular polarization signals were found in az-maps at all frequencies. Figure 8 shows the measured circular polarization as a function of the linear polarization amplitude measured in the detector frame (u) for samples with a signal-to-noise ratio of S/N > 5. Since circular polarization is correlated with the detector frame quantity instead of the sky signal (Q/U), the majority of the signal must be the result of a polarization leakage presumably due to inaccurate VPM transfer function modeling. Considering the uncertainties of the measurement in both directions, we use orthogonal distance regression to fit a scaling factor for each frequency. The results are empirical estimates of the aforementioned demodulation bias B vu for cloud polarization. These numbers are summarized in Table 2 and denoted as Bvu . For 90 GHz we found distinct trends before and after 2019 May 20, which is likely from a perturbation to the VPM system during telescope maintenance. By simulating the modulation and demodulation process with various VPM parameters and sky spectra, we found tight linear correlations between B vu and B uu ; therefore, the Bvu measurements can be used for B uu estimates in the absence of precise VPM parameter determination. These Buu estimates are also summarized in Table 2, and used to correct for the linear polarization bias in the following analysis.

The significant linear polarization and its unique spectral shape make clouds promising targets (among other sources with different spectral dependencies) for VPM transfer function calibration (Li et al. 2023). This possibility will be explored and applied to the cosmological analysis for future CLASS data.

Cloud Population Analysis

Confirmation with Camera Images

In addition to the visual confirmation from individual examples (e.g., Figure 7), we compare the radiometer detection to the cloud measure from the site cameras. Figure 9 shows the distribution of the cloud measure (CM 10 ) associated with CLASS radiometer detections. The average CM 10 values within a 15 minute window around each of the cloud candidates were extracted. For detections at all four frequencies, the associated higher-than-global cloud measures indicate a significant correspondence between the radiometer detection and the clouds optically identified by the cameras. The singlesided Mann-Whitney U test suggests that the CM 10 during the time period with radiometer polarization detection is significantly higher than the global distribution to at least the 22.7σ level (for the lowest 40 GHz detections). Due to the aforementioned FoV mismatch between the site camera and the telescopes and the finite cadence of the camera images, the significance shown here should be considered as a lower limit.

Polarization Angle

The polarization angle for each candidate was estimated as

where Ū and Q are the background-subtracted averaged Stokes amplitudes in the flagged region (Section 4.4). For boresight angles ±45°and 0°, either Q or U may have large uncertainties due to the degeneracy between the cloud polarization direction and the detector angle, and the simple estimate here would have large uncertainty and/or bias. Therefore, we avoid the interpretation of polarization angles from individual measurements and only study the cloud population as a whole.

Figure 10 shows the polarization angle distribution of each Table 2 CLASS Demodulation Bias Band 40 a 90 b 150 220 Buu 1.080/1.078 0.983/0.935 0.966 0.943 Bvu -0.215/-0.186 -0.124/-0.255 -0.200 -0.257 Notes. a The estimated bias terms are given separately for data before and after 2020 January 20. b The estimated bias terms are given separately for data before and after 2019 May 20. cloud candidate weighted by the event observation time. It is worth noting that the seven boresight rotations of the telescopes have different sensitivities to cloud polarization. At ±45°b oresight, the VPM and detectors have the most sensitivity to the Q polarization, but the uncertainty on the Stokes U component is the greatest; therefore, the polarization angle cannot be reliably determined. We find the average polarization angle to be 85°.7 ± 1°.5, 91°.1 ± 0°.5, 89°.6 ± 0°.5, and 86°.4 ± 0°.6 for the 40, 90, 150, and 220 GHz samples, respectively. Despite the unfavorable conditions for the polarization angle measurement, the distribution is consistent with the predicted 90°polarization angle from the scattering of a population of horizontally aligned cloud particles. The distributions in Figure 10 also show a slight preference for angles around 0°(180°). By visually inspecting the camera images, we found that they are associated with clouds at very low altitudes close to the ground (but not exclusively), and the polarization angle may be explained by +Q polarization from the scattering of the atmosphere radiation by spherical cloud particles like water droplets (Section 2).

Polarization Spectrum

CLASS's multifrequency cloud observations provide a unique opportunity to study the spectrum of the cloud's polarization and to verify the polarization mechanisms as outlined in Section 2.

In order to compare the multifrequency data, we chose the candidates detected from the 90 GHz data and performed forced photometry on the other three bands. The forced photometry finds corresponding pixels in other bands by matching the temporal and azimuthal bins, which might not always be possible due to the different pointing directions of the telescope mounts and/or the missing data. To account for the bias resulting from the VPM transfer function error, Buu was divided out from the polarization measurements from each band. The power measured by the CLASS bolometers also needs to be corrected for the absorption by the water vapor. This is especially relevant for cloud observations that typically occur during high PWV conditions. However, we lack dedicated PWV measurement to accurately determine the real-time PWV, and resort to the PWV values reported by the nearby APEX weather station at 1 minute cadence. For each cloud candidate, we queried the maximum PWV reported by APEX within a 10 minute window and calculated the atmospheric transmission based on the detector bandpass (Appel et al. 2022).

Figure 11 summarizes the polarization amplitudes (-Q) from the forced photometry for all the candidates with S/N greater than 4. The spectral index α can be obtained by jointly fitting the polarization amplitudes across two or more frequency bands. Using the 375 samples with detections at 90 and 150 GHz with an S/N > 4, we find the index α = 3.90 ± 0.06, consistent with the Rayleigh scattering model of ice particles as outlined in Section 2. However, when the 40 or 220 GHz samples are included, the best-fit α parameters significantly deviate from the nominal value of 4. These results suggest excess/deficit of power at 40/200 GHz, which we discuss below.

Effects from Liquid Water

At temperatures above -40 °C, liquid water may coexist with ice crystals in its supercooled state (Korolev & Milbrandt 2022). The dielectric permittivity of liquid water within the microwave range has a strong frequency dependency; therefore, the Rayleigh scattering spectrum of liquid water has an additional term |(ò ν -1)/(ò ν + 2)| 2 (Equation ( 5)). However, the liquid water fraction is expected to be low in high clouds; moreover, the spherical shape and smallness of the water droplets (about 20 μm, Lasher-Trapp et al. 2005) also make liquid water scattering a subdominant contributor to the observed polarization signals. Their absorption effect, however, could be important since the absorption cross section, in the Rayleigh regime, is proportional to D 3 (whereas the scattering is proportional to D 6 , where D is the particle diameter) and the liquid water absorption efficiency is a few orders-of-magnitude higher than that of ice due to its large imaginary part of the dielectric permittivity (Turner et al. 2016). Considering a toy model with a layer of liquid water at the base of the ice cloud, the liquid water absorption can be incorporated by inserting into the integral of Equation (4):

where k m is the mass absorption coefficient computed from the complex permittivity that also depends on the frequency and temperature T; LWP is the column mass density of the liquid water layer. The two secant terms correspond to the absorption of the incident and scattered light, respectively. The nadir angle (θ) dependency of the absorption adjusts the scattering contribution from different directions; therefore, the final result also depends on the particle shape/orientation in the ice cloud.

We compared two models of the polarization spectrum below: R-α is a power-law model that has the spectral index α as a free parameter; R-abs includes the liquid water layer absorption and assumes the Rayleigh scattering spectrum (α = 4). We jointly fit these models to 131 samples with joint detection at all four frequency bands with an S/N > 4, and the residuals are displayed in Figure 12. As mentioned in the section above, a purely Rayleigh scattering spectrum is not a good fit across all four frequencies. Specifically, the R-α model yields a best-fit α = 3.17 ± 0.05. The excess of power at 40 GHz and the favor of a lower spectral index are evident in Figure 12, and the mean χ 2 per data degree of freedom is 3.7. By including additional liquid absorption, model R-abs brings more consistency between 40 and 90 GHz, as indicated by the improved normality of the distributions and the lower mean χ 2 = 2.3. We note that the absorption model depends on the cloud temperature, LWP, and particle shape (through Δ z in Equation (2)), and the regression problem is underdetermined with the four-frequency photometry. Therefore, the result here only demonstrates the possibility that the presence of liquid water could play a role in the spectrum of cloud polarization, but the interpretation of the best-fit parameters should be conducted with caution.

Spectra at 220 GHz

To assess the potential deviation from Rayleigh scattering at around 220 GHz, we selected 375 samples with joint detection at 150/220 GHz with an S/N greater than 4. Figure 13 plots in black the ratio of the Stokes Q signal measured at 220 GHz compared to the 150 GHz measurement scaled by (220/148.2) 4 , assuming the Rayleigh scattering spectrum. The left panel shows the deficit ratio against PWV, and the right panel is the histogram weighted by the cloud observation time, centering around the mean value 0.61 ± 0.02. A weak correlation with PWV is found with Pearson coefficient -0.33 ± 0.03; however, it is unlikely that the deficit is due to the unaccounted-for water vapor absorption. The additional 40% absorption would require PWV values greater than 20 mm, which is improbable at the CLASS site.

Nonetheless, as discussed in the previous section, it is possible that supercooled liquid water is contributing to the absorption. We jointly fit the 90, 150, and 220 GHz data (315 samples with an S/N > 4) with the model R-abs, but fixed the cloud temperature to -15 °C to reduce the number of free parameters. The result of the 220 GHz deficit ratio after correcting for the absorption with the best-fit R-abs parameters is shown in Figure 13 in blue. The best-fit model can explain most of the deficit with a wide spread of LWP values-the 50, 84, and 95% quantiles are 0.4, 1.2, and 2.0 kg m -2 , respectively. These values are much higher than the typical value in the Chajnantor area (below 0.01 kg m -2 for 90% of the time, Kuo 2017) based on MERRA-2 or the maximum value of 0.6 kg m -2 found in the ERA5 database over the survey periods. However, these records are based on reanalyses and are averaged values over a wide area of the atmosphere, and might not be applicable to individual cloud measurements. Lacking targeted LWP measurement for each sample, we cannot preclude the possibility that the clouds are in a mixed phase and the liquid water absorption contributes (at least partially) to the deviation of the Rayleigh scattering spectrum at 220 GHz. Alternatively, the apparent deficit can be interpreted as evidence of a population of large ice particles that contribute to polarization through resonant scattering. As is shown in Figure 4, the polarization amplitude is sharply suppressed at around 220 GHz for a particle diameter 350 μm. Considering a population of cloud particles, depending on size distribution (which primarily depends on the temperature, Austin et al. 2009), the deficit ratios range between 20 ∼ 80% for temperature -50 ∼ -80 °C. We note that the calculation above assumed Mie scattering from spherical particles, and the exact deficit ratio might be sensitive to the size, shape, and orientation of the cloud particles, which is beyond the scope of this work.

Conclusion

We presented multifrequency observations of cloud polarization from CLASS at frequencies centered around 40, 90, 150, and 220 GHz. The az-mapping technique identifies transient sky polarization signals and is less susceptible to time-domain glitches than searches in time-ordered data. A blind search for localized polarization transient signals in daytime and nighttime throughout the 2016-2022 observation period yielded significant detection of cloud events that concur with the cloud images captured by the site cameras and correlate with the cloud measure extracted from the cloud images based on a CNN classifier.

Leveraging the polarization modulator, finite on-sky position angle offsets of the detectors at the 45°elevation, and the azmaps, we were able to combine the data from multiple detectors to perform polarization measurement of the transient signal from clouds on a single scan. At the threshold chosen for cloud detection, we identified 1481, 3677, 4110, and 4178 cloud candidates at 40, 90, 150, and 220 GHz, respectively. The polarization angle distribution of the majority of cloud events is consistent with the predicted 90°polarization angle from the scattering of a population of horizontally aligned ice crystals. At a lower significance, we also found a group of events with vertical polarizations, which may be attributed to the polarization effect of spherical cloud particles (e.g., liquid water).

Focusing on the cloud candidates selected by the 90 GHz data, we compared the -Q polarization signal at all four frequencies. The joint spectrum fit from 90/150 GHz data shows a mean spectral index of α = 3.90 ± 0.06 consistent with the ice crystal Rayleigh scattering. The polarization amplitudes at 40/220 GHz show an excess/deficit compared to the α = 4 spectrum. Using a toy model with a layer of supercooled liquid water at the base of the ice cloud, we found better agreement between the data and Rayleigh spectra with liquid water absorption. Although the fitting result is underdetermined due to the finite frequency sampling of CLASS data and lack of targeted LWP measurement, the result hints that the presence of water might be important to explain the spectral shape of cloud polarization, but the exact solution would depend upon the modeling of cloud particle size/shape distribution, their orientation, the ice/liquid water mixing, and cloud temperature. Alternatively, the 60% deficit ratio observed at 220 GHz measurement compared to the scaled 150 GHz data assuming α = 4 can be explained by a population of cloud particles with diameters greater than 350 μm that contribute to the polarization through the Mie scattering. This picture is also consistent with the long tail distribution of the ice crystal size in high clouds (Austin et al. 2009) within the temperature range of -50 ∼ -80 °C.

No circular polarization signal was detected from the clouds. All apparent circular polarization signals are consistent with a linear-to-circular polarization leakage from VPM modeling uncertainty. Leveraging the large amplitude of the cloud signal, this leakage signal will be used in future CLASS analyses to constrain the VPM parameters for the cosmological data set.

As a final remark, the cloud events analyzed in this work were deliberately selected from the brightest end of the population. For the cosmological analysis, the time-ordered data around these events are likely flagged in the pipeline either due to the anomalously large signal or coincident bad weather and/or excessive detector flux jumps. The fainter end of the cloud population, on the other hand, is more difficult to avoid with data selection and could impact the long-term stability of the instrument as a polarization noise (Harrington et al. 2021, Cleary et al. 2023, in preparation) Since 2022, we have installed an all-sky camera to better monitor the cloud conditions and to guide the survey. The methods developed and lessons learned here will be used to improve the CMB data reduction for CLASS, and are available for future experiments (Abazajian et al. 2016;Ade et al. 2019) that aim to push the limit of CMB observations further from the ground.