CLASS Observations of Atmospheric Cloud Polarization at Millimeter Wavelengths

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, sub-classified 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 area in liquid/solid phases. As one of the exceptionally arid sites for mm-astronomy, the Chajnantor area experiences a median $1~\mathrm{mm}$ PWV (Cortés et al., 2020). Water in the hydrometeor form is long-tail distributed, with LWP and IWP below $10^{-3}~\mathrm{kg\,m^{-2}}$ 80% of the time but could rise above $10^{-2}~\mathrm{kg\,m^{-2}}$ for $10\%$ of the time (Kuo, 2017, austral summer months excluded). This is also consistent with the typical $10^{-2}~\mathrm{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^{\circ}\times 0.25^{\circ}$ pixel centered at $67.8^{\circ}\,\mathrm{W},23.0^{\circ}\,\mathrm{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 $\lesssim 2~\mathrm{km}$ above the ground ($\lesssim 7~\mathrm{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.

In the Rayleigh limit where the scatterer size is small compared to the wavelength, the electric field component of the scattered radiation $\mathbf{E}_{o}$ at direction $\mathbf{r}$ from the scatterer follows¹¹1In most of the cases (Figure 2), the cloud distance is well within the far-field region of the telescopes ($\gg 200~\mathrm{m}$ for the $220~\mathrm{GHz}$ telescope).

where $V$ is the scatterer volume; $\nu$ is the frequency of the incoming radiation $\mathbf{E}_{i}$; $c$ is the speed of light. The polarizability tensor, $\boldsymbol{\alpha}$, 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 notations in Takakura et al. (2019), the polarizability tensor is parameterized as

where $\epsilon$ is the relative permittivity; $\mathbf{\Delta}=\mathrm{diag}\{\frac{1-\Delta_{z}}{2},\frac{1-\Delta_{z}}{2},\Delta_{z}\}$ is the depolarization matrix for spheroidal scatterers (with $\Delta_{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: $\mathbf{S}_{i}=[T(\theta),0,0,0]^{\mathrm{T}}$, with intensity profile

where $\theta$ is the zenith angle; $\theta_{\mathrm{h}}\approx\sqrt{2h/R_{\mathrm{e}}}$ is the astronomical horizon for clouds at altitude $h$ above the ground ($R_{\mathrm{e}}$ is the radius of the Earth); $T_{\mathrm{g}}$ is the ground temperature, and $T_{a}(\theta,\nu)$ is the brightness temperature profile of the atmosphere. This setup is illustrated in Figure 1. For spherical particles with $\boldsymbol{\alpha}\propto(\epsilon-1)/(\epsilon+2)\mathbf{I}$, the Stokes parameters from cloud reflection evaluated in the telescope frame can be expressed as

where $\mathbf{S}_{o}(\nu)$ is the radiation received by the telescope at $(-\delta,0)$, with $\delta$ being the telescope pointing elevation; Matrix $\mathbf{R}$ rotates Stokes vectors into the telescope frame from the scattering frame in which the phase matrix $\mathbf{M}$ is evaluated. The column of $\mathbf{M}$ for unpolarized incoming radiation is $(1/2)[(\cos^{2}\Theta+1),(\cos^{2}\Theta-1),0,0]^{\mathrm{T}}$, where $\Theta(\theta,\phi,\delta)$ is the scattering angle. The Rayleigh scattering optical depth is

where $\sigma_{\mathrm{sca}}$ is the Rayleigh scattering cross section; $n(\mathbf{r})$ is the number density of cloud particles with diameter $D$. The typical crystal size in cirrus is $100~\mathrm{\mu m}$, although the distribution is widespread (Lawson et al., 2019; Austin et al., 2009). At mm-wavelengths the dielectric permittivity for ice is, $\epsilon\simeq 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 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 $\delta=45^{\circ}$ elevation, from horizontally aligned spheroidal cloud particles with random azimuthal orientations as a function of the depolarization factor $\Delta_{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~\mathrm{g\,m^{-2}}$ LWP and a particle size $100~\mathrm{\mu m}$ (geometric mean diameter). The atmospheric brightness profile $T_{a}(\theta,\nu)$ 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 non-spherical scatterers with extreme aspect ratios and peaks at $\sim 30\%$ for oblate crystals at $45^{\circ}$ elevation. 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 $\lesssim 1\%$ for $\delta=45^{\circ}$. 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).

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 $\propto\nu^{4}$; (2) the forward scattering is enhanced compared to other directions; (3) the scattering phase matrix contains non-zero 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.²²2https://github.com/scottprahl/miepython 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~\mathrm{GHz}$ is suppressed in the Mie solution in a way that depends sensitively on particle size in the $300$ – $400~\mathrm{\mu m}$ range. While the majority of ice crystals are below this threshold, the size distribution is long-tailed (McFarquhar & Heymsfield, 1997; Heymsfield et al., 2002; Austin et al., 2009), and the larger particles can contribute significantly to the polarization signal (roughly $\propto D^{6}$). Assuming the size distribution of ice crystals with log-normal mean size $46~\mathrm{\mu m}$ (corresponding to the best-fit parameters at $-70~\mathrm{{}^{\circ}C}$, Austin et al., 2009) and standard deviation $0.226$, the black curve shows the population-averaged spectrum ratio. While Rayleigh scattering is a good approximation of cloud polarization at frequencies below $150~\mathrm{GHz}$, the Mie scattering effect cannot be ignored for higher frequencies. For realistic high cloud composition, the cloud polarization signal at the $220~\mathrm{GHz}$ band is a factor of a few smaller compared to the Rayleigh solution and is comparable to the signal at $150~\mathrm{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.

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^{\circ}$ scan 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 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 A. The CNN classifier labels each of the $240\times 240$ pixel sub-regions 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 sub-regions 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^{\circ}$ elevation ($\mathrm{CM}_{10}$) as a proxy of the cloud condition at the site.

The high-cadence all-day monitoring during the six-year survey of CLASS provides an assessment of the cloud conditions at the site on Cerro Toco. Figure 6 shows the $\mathrm{CM}_{10}$ statistics from 2016 through 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.³³3https://www.apex-telescope.org/ns/weather-data 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 around 50%, consistent with the cloud detection rate of daytime images taken from the neighboring POLARBEAR site (Takakura et al., 2019).

CLASS observes at a constant $45^{\circ}$ elevation and scans continuously in azimuth by $720^{\circ}$ before turning around at a rate of $1$–$2^{\circ}\,\mathrm{s}^{-1}$. Polarization is measured through the variable-delay polarization modulator (VPM, Chuss et al., 2012; Harrington et al., 2018, 2021), which modulates between the linear polarization (Stokes $U$ in the telescope coordinates) and circular polarization at around $10~\mathrm{Hz}$. For cosmological observations, the linear polarization angle coverage is achieved by the rotation of the Earth and the telescope boresight rotation from $-45^{\circ}$ to $45^{\circ}$ at $15^{\circ}$ increments every day. For a cloud transient event, however, the VPM is only sensitive to one polarization direction — e.g., $\pm$ Stokes $Q$ for $\mp 45^{\circ}$ boresight and Stokes $\pm U$ for $0^{\circ}$ boresight. Nevertheless, polarimetry is still possible since the $\sim 30^{\circ}$ azimuthal differential pointing of the detectors on the focal plane translates into a $\sim 15^{\circ}$ spread in position angle on the sky at $45^{\circ}$ 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).

The radiometer data analyzed in this work are taken from the CLASS survey from 2016 August to 2022 May. As the 90 GHz 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 mid-day. 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.

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 $\mathbf{D}$ is the raw time stream modulated by the VPM, $\mathbf{M}$ is the VPM modulation transfer function, and $\mathfrak{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~\mathrm{Hz}$ and so contributes to the single-detector demodulated linear polarization as a baseline offset⁴⁴4This synchronous signal drifts slowly as the temperature difference between the VPM and the air evolves (Harrington, 2018). On shorter time scales, this can be thought of as baseline offset.. 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 pair-differencing for the unmodulated experiment, Harrington et al., 2021, Cleary et al., in prep.) for the rest of the analysis. We keep the notations $u/v$ for the pair-add data for simplicity. The demodulation process depends on 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 of 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⁵⁵5Throughout this work, we distinguish the polarization signals in the detector frame and on the sky by small and capitalized letters. $[Q,U,V]$ through

where $\gamma$ is the detector position angle on the sky and $\phi_{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 $\mathbf{P}$ is the pointing matrix in Equation 9 and $\mathbf{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-time-scale 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^{\circ}$, and temporal resolution was chosen such that each continuous $360^{\circ}$ scan was binned into the same temporal bin. Depending on the scanning speed, the temporal resolution is around $180$–$360~\mathrm{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$ az-maps are shown in Figure 7.

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~\mathrm{GHz}$ $u$-az-map and the cloud detection flags for the span between 2022-01-25 and 2022-01-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.

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, non-zero 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 $\mathrm{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 $\hat{B}_{vu}$. For $90~\mathrm{GHz}$ we found distinct trends before and after 2019-05-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 $\hat{B}_{vu}$ measurements can be used for $B_{uu}$ estimates in the absence of precise VPM parameter determination. These $\hat{B}_{uu}$ 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.

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 ($\mathrm{CM}_{10}$) associated with CLASS radiometer detections. The average $\mathrm{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 single-sided Mann–Whitney U test suggests that the $\mathrm{CM}_{10}$ during the time period with radiometer polarization detection is significantly higher than the global distribution to the at least $22.7\sigma$ level (for the lowest $40~\mathrm{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.

The polarization angle for each candidate was estimated as

where $\bar{U}$ and $\bar{Q}$ are the background-subtracted averaged Stokes amplitudes in the flagged region (Section 4.4). For boresight angles $\pm 45^{\circ}$ and $0^{\circ}$, 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 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 $\pm 45^{\circ}$ boresight, 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\pm 1.5^{\circ}$, $91.1\pm 0.5^{\circ}$, $89.6\pm 0.5^{\circ}$, and $86.4\pm 0.6^{\circ}$ for the $40$, $90$, $150$, and $220~\mathrm{GHz}$ samples, respectively. Despite the unfavorable conditions for the polarization angle measurement, the distribution is consistent with the predicted $90^{\circ}$ 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^{\circ}$ ($180^{\circ}$). 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).

CLASS’s multifrequency cloud observations provide a unique opportunity to study the spectrum of the cloud 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~\mathrm{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, $\hat{B}_{uu}$ 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 four. The spectral index $\alpha$ 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~\mathrm{GHz}$ with $\mathrm{S/N}>4$, we find the index $\alpha=3.90\pm 0.06$, consistent with the Rayleigh scattering model of ice particles as outlined in Section 2. However, when the $40$ or $220~\mathrm{GHz}$ samples are included, the best-fit $\alpha$ parameters significantly deviate from the nominal value of $4$. These results suggest excess/deficit of power at $40/200~\mathrm{GHz}$, which we discuss below.