© ESO 2019

1. Introduction

The BICEP2 experiment, designed to measure the cosmic microwave background (CMB) polarization, claimed the first detection of primordial B-modes in 2014 (BICEP2 Collaboration 2014) by measuring the tensor-to-scalar ratio as $r=0.2_{-0.05}^{+0.07}$. One year later, a joint effort involving the Planck and BICEP2 collaborations (BICEP2/Keck & Planck Collaborations 2015) revised this bound by publishing an upper limit of r <  0.12, which was obtained by removing the residual dust contamination from the BICEP2 maps. More recently, BICEP2 and Keck Array experiments (Keck Array & BICEP2 Collaborations 2016) further reduced this constraint to r <  0.07, which represents the strongest constraint to date on inflationary gravitational waves.

Further improving the constraint on the tensor-to-scalar ratio represents a hard challenge for current and future CMB experiments, which must accurately take into account gravitational lensing and foreground removal (Errard et al. 2016) in addition to having an excellent control of systematic effects (Wallis et al. 2017). Concerning the control of systematics, some experiments are designed with the capability of self-calibrating and thus of removing some systematic effects (Piat et al. 2012). For other experiments that lack this feature, an accurate instrumental systematics prevision and laboratory calibrations are mandatory in order to separate systematic errors from scientific data (Natoli et al. 2018; Inoue et al. 2016; D’Alessandro et al. 2015; Johnson et al. 2017).

A standard device for polarization analysis is a Stokes polarimeter, composed of a half-wave plate (HWP) and a polarizer. The HWP (Pisano et al. 2008) induces a phase shift on the linear polarization by birefringence, and the polarizer selects one component, so that by rotating the HWP it is possible to modulate the linearly polarized fraction of the incoming radiation and extract the Stokes components (T, Q, U), given a reference frame. The HWP can spin fast (>1 Hz), see, for example, LSPE Collaboration (2012), Columbro et al. (2019), EBEX Collaboration (2018) and Thornton et al. (2016), in order to modulate the signal at high frequency, and it can also rotate step by step (<1 Hz) (Piat et al. 2012; Bryan et al. 2016) depending on the scanning strategy of the experiment. Systematic effects such as temperature-to-polarization leakage, chromatic HWP behavior, and scanning strategy, have been evaluated in Essinger-Hileman et al. (2016), Salatino et al. (2011, 2018), and Takakura et al. (2017) and were also measured by Kusaka et al. (2014).

We here analyze the systematic effect that is induced in observation of the CMB polarization by the precession of the HWP, in the specific case of a transmissive plate. We assume that the HWP free from other systematic effects.

In Sect. 2 we review the dynamic of a spinning cylindrical plate and define the precession rate. In Sect. 3 we provide an analytic study of the effect that is induced by the spinning and preceding HWP using the 3D Jones formalism for fully polarized radiation and the 3D generalization of the Müller formalism for partially polarized radiation. In Sect. 4 we use the equations derived before to show some results on electric field components produced with numerical simulations. Finally, we describe in Sect. 5 the effect that is induced on full-sky observation of the CMB, assuming a satellite mission with different scanning strategies.

2. Precessing half-wave plate theory

In this section we present equations to describe a precessing body that rotates along one symmetry axis. We introduce the main variables and show their evolution with time. We report here only the main equations that are essential for the results shown in the following sections. All the details of the computation are provided in Appendix A.

We approximate the HWP as a cylindrically symmetric rigid body, like a coin, and we define a reference system x̂–ŷ–ẑ, with the x̂–ŷ plane coincident with the base of the cylinder, see Fig. 1. We hypothesize that the HWP has a large spin angular momentum L_s = I_sω_s along the symmetry axis, where I_s and ω_s are the moment of inertia and the angular velocity, respectively. In the unperturbed case, L_s coincides with the ẑ axis.

Fig. 1. Illustration of rigid-body precession where we show the angular momentum dislocations with respect to the x̂ and ŷ axes.

We now introduce a small angular perturbation (θ₀) of the spin momentum L_s, expanding it in its components along the x̂ and ŷ axes, η and ξ. By introducing the precession rate γ, defined as the rotation velocity of L_s around the ẑ axis, η and ξ are simply defined as

$$\eta ={\theta }_0\cos (\gamma t),$$(1)

$$\xi ={\theta }_0\sin (\gamma t).$$(2)

These equations describe the torque-free precession of a spin axis that rotates in space with fixed angle θ₀.

The precession rate is given by γ  =  ω_sI_s/I_⊥ (Kleppner & Kolenkow 1973), after I_s and I_⊥ are defined as the inertia momenta parallel and perpendicular to the ẑ axis. Assuming that the cylindrical plate is a thin disk, we can approximate I_s = 2I_⊥ and γ = 2ω_s, which means that the disk wobbles twice as fast as it spins.

The apparent rate of a thin-disk precession for an observer on the rigid body is given by

$$\begin{array}{c}\hfill {\gamma }^{\prime }=\gamma -{\omega }_{\mathrm{s}}={\omega }_{\mathrm{s}}\left(\frac{I_{\mathrm{s}}-I_{\perp }}{I_{\perp }}\right)\sim {\omega }_{\mathrm{s}}.\end{array}$$(3)

In a torque-free precession we can identify two different rotations according to the reference frame we consider. In the fixed laboratory frame, the angular velocity vector rotates around the fixed ẑ axis (where the angular momentum vector lies), tracing the so-called space cone. In the reference frame integral with the rotating body, the angular momentum and angular velocity vector describe a circle around the symmetry axis of the cylinder, tracing the so-called body cone with a precession rate γ′.

3. 3D generalization of a Stokes polarimeter

A description of Stokes polarimetry by using Müller and Jones matrices is available in the literature, see, for instance, Bryan et al. (2010a), O’Dea et al. (2007), and Chuss et al. (2012). In this section and in Appendices B and C, we derive the analytic equation for a traditional Stokes polarimeter and for a Stokes polarimeter where the HWP has a precessional motion.

3.1. Jones formalism

We now use an extended Jones formalism (Jones 1941) to retrieve a formula for electromagnetic (EM) field intensity after the Stokes polarimeter. The Jones matrices used are defined in Appendix B followed by explicit calculations. The traditional Jones formalism can describe the state of fully polarized light with a two-dimensional vector i = (E_x, E_y) and optical elements that change the state of input radiation with 2 × 2 matrices. We extend the Jones vector in three dimensions, i = (E_x, E_y, 0), where the ẑ-axis is the propagation direction of the EM field, and the x̂–ŷ plane is perpendicular to ẑ, thus any optical element is described by 3 × 3 matrices (Sheppard 2011, 2014; Ortega-Quijano & Arce-Diego 2013). Combining these matrices, we obtain a general formula for the Jones vector for the on-axis detector of a precessing polarimeter:

$$\begin{array}{c}\hfill i_{\mathrm{out}}=J_{\mathrm{pol}}·J_{\mathrm{RotY}}^{-1}(\xi )·J_{\mathrm{RotX}}^{-1}(\eta )·J_{\mathrm{HWP}}(\theta )·J_{\mathrm{RotX}}(\eta )·J_{\mathrm{RotY}}(\xi )·i_{\mathrm{in}}.\end{array}$$(4)

For the ideal case with η = ξ = 0°, when no precession occurs, we obtain the outgoing intensity from an ideal polarimeter:

$$\begin{array}{c}\hfill I=\frac{1}{2}[T+Qcos(4\theta )+Usin(4\theta )].\end{array}$$(5)

For the general case with η, ξ ≠ 0° we find the intensity at the detector:

$$\begin{array}{c}\hfill I=i_{{\mathrm{out}}_x}^2+i_{{\mathrm{out}}_y}^2+i_{{\mathrm{out}}_z}^2={(E_{x}g-E_{y}f)}^2,\end{array}$$(6)

where we define the modulating functions g(ξ, θ) and f(η, ξ, θ):

$$\begin{array}{cc}\hfill g& ={sin}^2\left(\xi \right)+{cos}^2\left(\xi \right)cos\left(2\theta \right),\hfill \\ \hfill f& =2sin\left(\theta \right)cos\left(\xi \right)(sin\left(\eta \right)sin\left(\xi \right)sin\left(\theta \right)-cos\left(\eta \right)cos\left(\theta \right)).\hfill \end{array}$$(7)

We can write the outgoing intensity through the Stokes parameters as follows:

$$\begin{array}{c}\hfill I=\frac{1}{2}(g^2+f^2)T+\frac{1}{2}(g^2-f^2)Q+gfU,\end{array}$$(8)

which is valid only if T² = Q² + U². The last equation gives the intensity for the on-axis detector of a Stokes polarimeter whose modulating element describes a torque-free precession. It is not merely a function of the HWP orientation about the ẑ-axis, but also depends on the displacements from the x̂- and ŷ-axes due to the precession.

3.2. Müller formalism

We now use a 3D extended Müller formalism (Sheppard et al. 2016; Samim et al. 2016) to retrieve a formula for an EM field intensity outgoing from the Stokes polarimeter. The Müller formalism (fully defined in Appendix C) is required for the case of CMB experiments because we wish to propagate partially polarized radiation through a polarimeter in order to extract the information about its polarization state. Starting from 3D extended Jones matrices (Appendix B), the 9 × 9 Müller matrix corresponding to each optical element can be easily obtained from Eq. (C.1):

$$\begin{array}{c}\hfill M_{\mathit{ij}}=tr({\sigma }_i·J·{\sigma }_j·J^{†}),\end{array}$$(9)

where σ_n (n = [0, …, 8]) are the trace-normalized Gell-Mann matrices (Eq. (C.2)). The Stokes polarimeter Müller matrix becomes

$$\begin{array}{cc}\hfill M_{{\text{SP}}_{\mathrm{wob}}}(\eta ,\xi ,\theta )=& M_{\mathrm{pol}}·M_{\mathrm{RotY}}^{-1}(\xi )·M_{\mathrm{RotX}}^{-1}(\eta )·M_{\mathrm{HWP}}(\theta )\hfill \\ \hfill & \times M_{\mathrm{RotX}}(\eta )·M_{\mathrm{RotY}}(\xi ).\hfill \end{array}$$(10)

The output Stokes vector therefore is

$$\begin{array}{c}\hfill s_{\mathrm{out}}=M_{{\text{SP}}_{\mathrm{wob}}}(\eta ,\xi ,\theta )·s_{\mathrm{in}}.\end{array}$$(11)

Each component of the Stokes vector is defined in Eq. (C.4), and assuming the field entering the polarimeter has E_z = 0,

$$\begin{array}{c}\hfill s_{\mathrm{in}}^T=({\mathrm{\Delta }}_0,{\mathrm{\Delta }}_1,0,0,{\mathrm{\Delta }}_4,0,0,0,\frac{1}{\sqrt{2}}{\mathrm{\Delta }}_0).\end{array}$$

Based on Eq. (C.5), which describes the relation between the 3D Stokes components and the classic definition of Stokes parameters, we can determine the outgoing intensity:

$$\begin{array}{cc}\hfill I& =\frac{1}{\sqrt{6}}({\mathrm{\Delta }}_0+3(m_{01}{\mathrm{\Delta }}_1+m_{04}{\mathrm{\Delta }}_4+m_{08}{\mathrm{\Delta }}_8))\hfill \\ \hfill & =(\frac{1}{3}+\frac{\sqrt{2}}{2}{m}_{08})T+\frac{\sqrt{6}}{2}{m}_{01}U+\frac{\sqrt{6}}{2}{m}_{04}Q,\hfill \end{array}$$(12)

where m_ij are the components of the Stokes polarimeter Müller matrix and are a function of the wobbling angles η, ξ. When η = ξ = 0, the Eq. (12) gives the well-known formula of the Stokes polarimeter (Eq. (5)).

4. Phenomenology

Starting from Eq. (B.12), we built a simulation to show the effects of the precession for a specific case. We set the spinning frequency at 1 Hz and the ratio I_⊥/I_s = 0.502. The Jones vector used as input field in the simulation was i_in = (1, 0, 0), corresponding to a (1, 1, 0, 0) Stokes vector. Typical outputs of the simulations are reported in Fig. 2. We show the fractional residuals, defined as the difference between the intensity from a precessing HWP and from the ideal non-precessing case, normalized to the maximum intensity of the ideal polarimeter.

Fig. 2. (a) Fractional residuals between the total intensities in a Stokes polarimeter with and without HWP precession. We fix the ratio I⊥/Is = 0.502 and consider different HWP precession angles θ0 = [1°,0.5°,0.25°], normalizing the difference with respect to the maximum intensity of an ideal Stokes polarimeter. (b) Fractional residuals, as before, but in a shorter time interval to highlight the beats induced by precession. (c) Power spectra with peaks at multiples and submultiples of the HWP spinning frequency with an amplitude related to θ0.

The time lines reported in Fig. 2a show the fractional residual for different amplitudes of the precession angle θ₀. These time lines clearly show several beats with an amplitude depending on the precession angle θ₀ in a nonlinear way.

In Fig. 2c we report the power spectra of time lines for the ideal case (red) and for the precessing cases. The ideal case presents a single line at 4 Hz, as expected from a wave plate that spins at 1 Hz. The beats from a precessing HWP produce spurious peaks at multiples and submultiples of the spin frequency. The most noteworthy peak in the spectra is the one at lower frequencies because it is the result of the very slow modulation imposed by the precession, which in Fig. 2a has a period of about 60 s, implying a peak at 0.016 Hz in the spectrum. In the next section we illustrate how these beats depend on the geometrical parameters of the plate.

4.1. Spinning speed and ratio I_⊥/I_s effects

For a cylindrical HWP, including its support, with mass m, thickness h, and radius R, the components of the moment of inertia with respect to the principal axes are

$$I_{\text{s}}=\frac{1}{2}m{R}^2;$$(13)

$$I_{\perp }=\frac{1}{12}m\left(3R^2+h^2\right),$$(14)

where we assume a diagonal inertia tensor:

$$\begin{array}{c}\hfill \mathbf{I}=\left[\begin{array}{ccc}{I}_{\perp }& 0& 0\\ 0& I_{\perp }& 0\\ 0& 0& I_{\mathrm{s}}\\ \end{array}\right].\end{array}$$(15)

The frequency for the precessional motion is directly linked to the HWP regular spin frequency f_s and to the ratio I_⊥/I_s (Eq. (A.5)).

We can note that this ratio only depends on the cylinder height and radius as

$$\begin{array}{c}\hfill I_{\perp }/I_{\mathrm{s}}=0.5+\frac{1}{6}{\left(\frac{h}{R}\right)}^2.\end{array}$$(16)

We therefore explore different configurations as shown in Fig. 3, where we show the fractional residual with respect to the ideal case. We consider input radiation that is horizontally polarized, a precession angle θ₀  =  1°, spin frequencies f_s  =  [0.1, 0.5, 1.0, 2.0] Hz and I_⊥/I_s = [0.501, 0.502, 0.506, 0.513]. As an example, these values for the I_⊥/I_s ratio correspond to an HWP with mass m = 1 kg, radius R  =  16 cm, and thickness h  ≃  [1.2, 1.7, 3.0, 4.5] cm; then I_s is fixed to 0.0128 kg m². In practice, the ratio of the components of the moment of inertia does not only depend on the plate thickness, but also on the structure of the HWP support.

Fig. 3. Fractional residual for different ratios I⊥/Is and different HWP spinning frequencies. Each column shows a different fs value and each row a different I⊥/Is ratio. A thinner HWP, I⊥/Is → 0.500, has beats in the intensity over an extended period and a thicker HWP has shorter beats. By looking at the figure from left to right, it is clear that the spinning speed compresses the beats. The maximum amplitude remains the same because it only depends on θ0, which is fixed to 1° in this particular simulation.

Figure 3 shows the dependence of the effect from the HWP spinning speed and the inertia tensor: the simulation shows that a thinner HWP, I_⊥/I_s → 0.500, has beats in the intensity over an extended period and a thicker HWP has shorter beats. This is clear by looking at Fig. 3 from top to bottom. The value 0.5 is not possible because it corresponds to a null thickness.

The effect of different spinning speeds is to compress the beats. This is clear by looking at Fig. 3 from left to right. The maximum amplitude remains the same because it depends only on θ₀, which is fixed to 1° in this particular simulation.

In Fig. 4 we report the power spectra of the time lines for different ratios I_⊥/I_s. The spectra exhibit the effect discussed above, showing that the beat frequency moves to lower values as the ratio I_⊥/I_s is reduced toward the minimum value of 0.5 (for I_⊥/I_s = 0.513, the period of the beats is 10 s, corresponding to a peak in the spectrum at 0.1 Hz).

Fig. 4. Power spectra extracted from time lines showing the shifting of the spurious peak at low frequencies. The simulation considers input light that is horizontally polarized, a precession angle θ0  =  1°, and I⊥/Is  =  [0.501, 0.502, 0.506, 0.513].

4.2. Minor effect

As reported for example in Gu et al. (2018), a tilted HWP also changes its transmission properties because path of radiation inside the plate is different from the case of normal incidence. The resulting effect is a variation of retardance, which can be as high as 40% in the case discussed in Gu et al. (2018) for a source at 347 nm with a tilt of 5°. This effect is multiplicative with the ratio wavelength to thickness, which is much higher in the case of plates used in the optical bandwidth than in the case of the millimeter wavelength. In instrumentation devoted to observations of the CMB polarization, the plate thickness is on the same order of magnitude as the wavelength: 3.1 mm POLARBEAR-2 (Hill et al. 2016), 3.05 mm SPIDER (Bryan et al. 2010b), 1.62 mm EBEX (EBEX Collaboration 2018), 3.2 mm QUBIC (Aumont et al. 2016), 3 mm LSPE-SWIPE (LSPE Collaboration 2012), and 3.05 mm ABS (Kusaka et al. 2014). The impact of the effect for millimeter astronomy can be examined through the variation between the input and output Stokes vector, and it is estimated to be orders of magnitude smaller than the effect due to electromagnetic field projection that is analyzed in the rest of this paper.

5. Full-sky simulations

In order to test the effect of the HWP precessional motion on CMB observations, we built an algorithm that can generate a realistic satellite scanning strategy in presence of spinning HWP and produces data time-lines. We completed this software with a map-making algorithm that collapses data time-lines into maps. All simulations are noise-free to better capture the impact of systematic effects.

5.1. Simulation pipeline

The scan simulator takes as inputs the details of a satellite scanning strategy, three spin rates, and two angles (see Das & Souradeep (2015) for a detailed description of the geometrical configuration):

• Earth revolution velocity ω₁,

• precession velocity ω₂,

• satellite spin ω₃,

• precession angle α, that is, the angle between the satellite spin axis and the Sun-Earth direction,

• boresight angle β, that is, the angle between the focal plane direction and the spin axis.

We only simulated a single detector placed at the center of the focal plane by illuminating a spinning HWP with f_s frequency. The systematic affecting the HWP was included in the data at the time-line level, and a simple rebinning map-making was used to average all the samples in T,Q,U Stokes parameters maps (Tegmark 1997). In this paper we consider Planck (Tauber et al. 2006), WMAP (Bennett et al. 2003), COrE (Natoli et al. 2018), and LiteBIRD (Sekimoto et al. 2018) scanning strategies. The input parameters we employed for these scanning strategies are listed in Table 1. As sampling rate we used 60 Hz.

5.2. Maps and results

5.2.1. Input map

The input sky map we used for full-sky simulations contained the solar dipole and Galactic diffuse foregrounds in temperature and a CMB realization, both in temperature and polarization. We decided to include foregrounds only in temperature in order to highlight the temperature-to-polarization leakage. The input 𝒞_ℓ used for the CMB realization is compatible with the best fit of the Planck 2015 release (Ade et al. 2016) with no tensor perturbations.

The foreground field was generated from the Commander solution delivered with the Planck 2015 release (Adam et al. 2016). It includes the primary temperature emissions (Planck Collaboration ES 2015): synchrotron, free-free, spinning dust, CO, and thermal dust emission, without considering their polarization contribution.

This map was modeled in order to highlight the temperature-to-polarization leakage induced by the HWP precession during the observations. We set the resolution parameter of the input map at HEALPix (Gorski et al. 2005) resolution N_side  =  256, comfortable enough with respect to the Gaussian beam with FWHM  =  60 arcmin. In order to evaluate the effect of the parameters chosen for the simulation, we ran a case with N_side  =  128, finding the same results in terms of angular power spectra residual, except at smaller scales, where the pixel size matters independently of the presence of systematic effects.

5.2.2. Output maps

We performed several simulations with different configurations for the HWP. We varied the spin frequency, precession angle, and I_s/I_⊥ ratio. For each simulation we compared input and output maps and computed the B-mode power spectrum. As a visual example, we show in Fig. 5 (top panel) the output maps for a mission for which we adopted a LiteBIRD-like scanning strategy and solved the Stokes parameters through an HWP with a spin frequency of 1 Hz, a precession angle of θ₀ = 1°, and I_⊥/I_s = 0.514.

Fig. 5. Top panel: T, Q, and U maps reconstructed through a Stokes polarimeter where the HWP wobbles. Bottom panel: difference between input and output T, Q, and U maps showing the effect of the wave-plate precession on observed maps.

The residual maps (i.e., difference between output and input maps) in Fig. 5 (bottom panel) show the effect of the HWP wobble that induces variations of a few percent with respect to the input map. The effect is noticeable close to the Galactic plane and close to the maximum and minimum of the solar dipole, where the intensity emission is higher.

The effect on the maps is generated by the coupling between the satellite spin and the polarization modulation, which is affected by the precession, therefore we decided to test several conditions. In particular, when we slowed the HWP spin speed down to 0.1 Hz, the effect of the precession is emphasized. We show this in Fig. 6, where the systematic effect induced in the T, Q, and U maps, reported in histogram-equalized color scale, is at the same level as in the input map.

Fig. 6. Temperature and polarization difference maps showing the effect of the wave-plate precession for a slow-spinning modulator. Modulation parameters: fs = 0.1 Hz, θ0 = 1°, I⊥/Is = 0.514.

5.2.3. Results

The angular power spectra from the output maps shown in Fig. 5 are reported up to l ∼ 200, given the limit imposed by the beam. The relative variations for both E- and B-modes (Fig. 7) show the effect of the precession, combined with the satellite spin, which dominates at small angular scales (l >  150).

Fig. 7. Input minus output relative difference of EE and BB power spectra computed from the maps shown in Fig. 5.

When a 0.1 Hz spinning HWP is used, the synchronism with the satellite is slightly different and spurious B-mode polarization is detected at different angular scales (l >  100). These two cases show different matching between the systematic effect and the satellite spin (Figs. 8a and 8b).

Fig. 8. Panels a and c: BB and EE power spectra extracted from the input and output maps of Fig. 6; in panel c the two lines overlap. Panels b and d: normalized difference between input and output power spectra.

5.3. Scan strategy comparison

The few simulations presented so far, assuming a LiteBIRD-like scanning strategy, show the effect of the HWP precession on full-sky maps and angular power spectra. Because the scan strategy can have a role in mitigating this systematic effect that couples temperature and B-mode polarization (Wallis et al. 2017), we implemented simulations, as described in Sect. 5.1, that can reproduce different satellite observational strategies. The results obtained from analyzing these simulations are reported as residual maps (i.e., difference between output and input maps). They are shown in Fig. 9 as root mean square (hereafter RMS) of the residual maps, shown in Table 2, and as B-modes angular power spectra, shown in Fig. 10.

Fig. 9. Difference between polarization Q-maps reconstructed through a Stokes polarimeter showing the direct effect of the wave-plate precession with different scan strategies. Modulation parameters: fs = 0.1 Hz, θ0 = 1°, I⊥/Is = 0.508 for Planck, I⊥/Is = 0.514 for COrE and LiteBIRD, and I⊥/Is = 0.510 for WMAP. Scan strategy from the left to the right: Planck, WMAP, COrE, and LiteBIRD.

Fig. 10. B-mode power spectra (a) recovered from simulations including the systematic and the difference output minus input of BB power spectra (b) normalized with respect to the input spectrum for the Planck, WMAP, COrE, and LiteBIRD scan strategies with different HWP spinning speeds (fs = [0.1, 1] Hz). The dashed line represents the primordial B-mode angular power spectrum in the case r = 0.001, divided by the lensing only B-modes; the other cosmological parameters are assumed from Planck Collaboration I (2018).

In Fig. 9 we report the Q residual maps, in histogram-equalized color scale, for the case f_s = 0.1 Hz, θ₀ = 1° (the U maps show variations with a similar pattern and similar dynamic range). We made the simulations with several HWP physical parameters. Here we report the results for the following values of the I_⊥/I_s ratio as representative cases: 0.508 for Planck 0.514 for COrE and LiteBIRD, and 0.510 for WMAP.

The result of this analysis shows that the residual does not only depend on the scanning strategy, but mostly on the combination of scanning strategy, HWP rotation speed, and I_⊥/I_s ratio.

For example, all scanning strategy simulations show the largest effect in the case of an HWP spinning at 0.1 Hz, and they show the lowest residual in the case of 1 Hz spin frequency. This is also evident from the RMS value, reported in Table 2, and from the angular power spectrum in Fig. 10. Some strategies produce a spurious peak in the angular power spectrum that might be induced by a resonance between satellite spin and HWP wobbling, that is, COrE 0.1 Hz at ℓ ∼ 40, WMAP 0.1 Hz at ℓ ∼ 85, or LiteBIRD 0.1 Hz at ℓ ∼ 120.

The Planck-like scanning strategy (Tauber et al. 2006) does not create particular patterns or structures on larger angular scales, as can be seen in the BB power spectra that we plot in orange and blue in Fig. 10.

On the other hand, the COrE-like simulation (with slowly spinning HWP) shows the worst coupling between the satellite spin and wave-plate precession in terms of angular structures at large scales, as is visible in terms of spurious B-modes (plotted in brown in Fig. 10). These simple cases show that the large-scale patterns in the residual map are not related to the whole quality of the map, which is better described by the angular power spectrum.

The power spectra and the maps we recovered show the contamination that is generated by the half-wave plate precession systematic for a specified scan strategy. Repeating the analysis with different precession angles, we conclude as expected that the larger the precession angle, the stronger the spurious B-mode signal, and the higher the HWP spin frequency, the greater the mitigation of the systematic effect.

In terms of research for primordial B-modes, the faster rotation of the wave plate helps to mitigate the systematic effect that is induced by the precession of the modulating element in a Stokes polarimeter by moving the contamination at high ℓ. Fig. 10b illustrates the fractional residual B-modes due to observation with a wobbling HWP in the case of no-tensor perturbations, but only lensing-induced B-modes. The fractional residual power is a good figure of merit of the contamination because next-generation CMB polarization experiments are designed to reach a sensitivity that is usually quantified as a fraction of the lensing-induced B-modes level.

5.4. Temperature-to-polarization leakage

When we include a map with only temperature foregrounds in the simulation input, we can highlight the temperature-to-polarization leakage effect that is induced by the systematic for various scenarios. We verified that polarization foregrounds, removed with the ideal component separation method, leave a residual that is one order of magnitude lower in terms of polarization-to-polarization leakage. The HWP wobble induces B-modes whose amplitude is proportional to θ₀, shown in Fig. 11, where we report the BB power spectra for different precession angles [0.5°, 1°, 1.5°, 2°] extracted from the maps scanned by a COrE-like satellite. The induced B-mode signal exceeds the gravitational lensing contribution already for θ₀ = 1°.

Fig. 11. B-mode power spectra recovered from simulations including the systematic for the COrE scan strategy with fixed HWP spinning speed (fs = 0.1 Hz) and different precession angles. The input map includes only temperature contributions; the polarization contribution in the output map arises from the HWP precession.

The output polarization components Q and U are shown in Table 3.

5.5. Neighboring detector to mitigate the systematic

The systematic effect induced by the wobbling can be mitigated by redundancy. Observing the same sky pixel with different phases of the wobbling plate averages the contamination out. This can also be obtained by combining multiple detectors and observing the same pixel at different times. In order to determine this mitigation, we simulated the observation with two different detectors, pointing to different boresight angles β, shifted by 1°, for all the proposed scanning strategies. The modulation parameters were f_s  =  1 Hz, θ₀  =  1°, I_⊥/I_s  =  0.508 for Planck, I_⊥/I_s  =  0.514 for COrE and LiteBIRD, and I_⊥/I_s  =  0.510 for WMAP. The combination of the data from the two detectors resulted in a single map with a reduced contamination with respect to the single-detector maps, as reported in Table 4. In the table we report the RMS of the difference between the map with and without the induced systematic effect. This RMS of the residual is very similar for the two single-detector maps, and is reduced in the map we produced when they were combined. The RMS of the residuals scales with a factor $\sim \sqrt{2}$, indicating that the contamination is rather uncorrelated among the two detectors.

6. Conclusions

We described the systematic error induced by the precession of the half-wave plate modulator in a Stokes polarimeter and the effects on CMB polarization measurements. In Sect. 3 we derived the analytical equation (Eq. (12)) for the power on the central detector of the polarimeter when the HWP precesses with an angle θ₀ and with a frequency γ imposed by the physical properties of the device (the spin frequency f_s and the I_⊥/I_s ratio, including its support). Using Eq. (12), we performed several simulations to assess the level of the systematic effect that is induced by the precession. We found a strong dependence on θ₀ and γ both for the fractional residual of the signal and for the power spectrum.

In Sect. 5 we developed the simulation of full-sky observation by a satellite mission, including the systematic effect, and quantified its effect on B-modes that were retrieved from the output map. The HWP precession produces spurious effects on the maps at different angular scales depending on the strategy used to scan the sky. In particular, we implemented four scan strategies, WMAP-, Planck-, COrE-, and LiteBIRD-like. With the specific mechanical properties of the implemented HWP, our simulations show a B-mode induced by leakage from the intensity signal that dominates at different scales for the strategies selected: at ℓ ∼ 40 for COrE-like satellite, at ℓ ≃ 85 for a WMAP-like, and ℓ ≃ 120 in the LiteBIRD-like case. The new-era CMB experiment will gain some order of magnitude in sensitivity (Abitbol et al. 2017), a few 10⁻⁴μK, compared to today. The analytical formula of the systematic effect that is induced by the HWP wobbling enables removing or at least to predicting its effect on data.

In general, the precession effect induces beats in the time lines. In amplitude and frequency, these beats are related to the physical properties of the moving parts. Their effect on the maps depends on specific scanning parameters and on possible coupling with the beat frequency. We recommend that this effect be taken into account in the design of an observation strategy by modeling and measuring the inertia tensor of the moving parts. When the tensor is measured (or modeled), it can be inserted into Eq. (12) to simulate the effect on the timelines. The scanning strategy must avoid any synchrony with the beat frequency. In this case, redundancy helps to cancel most of the contamination, but considering the targeted sensitivity of future CMB experiments (Sekimoto et al. 2018), the precession must still be considered as a possible source of systematic effect.

Acknowledgments

We acknowledge the support of the ASI-COSMOS programme, Elia Battistelli, for CPU time on the SPINDUST server. The work has been supported by University of Rome “Sapienza” research funds. LP acknowledges the support of the CNES postdoctoral program.

References

Appendix A: Precession theory

We report in this appendix the derivation of the main equations that describe a precessing disk. In particular, we detail the derivation of Eq. (3).

The contribution to L_x by the rotation about the x-axis is $L_x=\frac{\mathrm{d}(I_{\mathit{xx}}\eta )}{\mathrm{d}t}=I_{\mathit{xx}}\frac{\mathrm{d}\eta }{\mathrm{d}t}$. We can treat I_xx as a constant because moments of inertia about the principal axes are constant for small angular displacements. In addition, the rotation about the y-direction contributes to L_x by giving a component L_s sin (ξ) on the x-axis. Combining these contributions, we obtain

$$\begin{array}{c}\hfill L_x=I_{\mathit{xx}}\frac{\mathrm{d}\eta }{\mathrm{d}t}+L_{\mathrm{s}}sin(\xi ),\\ \hfill L_y=I_{\mathit{yy}}\frac{\mathrm{d}\xi }{\mathrm{d}t}-L_{\mathrm{s}}sin(\eta ).\end{array}$$(A.1)

Because I_xx = I_yy = I_⊥, and by exploring small-angle first-order approximation, Eq. (A.1) read

$$\begin{array}{c}\hfill L_x=I_{\perp }\frac{\mathrm{d}\eta }{\mathrm{d}t}+L_{\mathrm{s}}\xi ,\\ \hfill L_y=I_{\perp }\frac{\mathrm{d}\xi }{\mathrm{d}t}-L_{\mathrm{s}}\eta .\end{array}$$(A.2)

Furthermore, based on the same approximation, L_z = I_sω_s. Because we consider a torque-free system (dL/dt = 0), both L_s and ω_s are constant, leading to

$$\begin{array}{cc}\hfill I_{\perp }\frac{{\mathrm{d}}^2\eta }{\mathrm{d}{t}^2}+L_{\mathrm{s}}\frac{\mathrm{d}\xi }{\mathrm{d}t}& =0,\hfill \\ \hfill I_{\perp }\frac{{\mathrm{d}}^2\xi }{\mathrm{d}{t}^2}-L_{\mathrm{s}}\frac{\mathrm{d}\eta }{\mathrm{d}t}& =0.\hfill \end{array}$$(A.3)

By introducing ${\omega }_x=\frac{\mathrm{d}\eta }{\mathrm{d}t}$ and ${\omega }_y=\frac{\mathrm{d}\xi }{\mathrm{d}t}$, Eq. (A.3) becomes

$$\begin{array}{cc}\hfill I_{\perp }\frac{\mathrm{d}{\omega }_x}{\mathrm{d}t}+L_{\mathrm{s}}{\omega }_y& =0,\hfill \\ \hfill I_{\perp }\frac{\mathrm{d}{\omega }_y}{\mathrm{d}t}-L_{\mathrm{s}}{\omega }_x& =0.\hfill \end{array}$$(A.4)

In order to solve this coupled system of differential equations, we can differentiate one and substitute the other:

$$\begin{array}{c}\hfill \frac{{\mathrm{d}}^2{\omega }_x}{\mathrm{d}{t}^2}+{\gamma }^2{\omega }_x=0\mathrm{with}\gamma =\frac{L_{\mathrm{s}}}{I_{\perp }}={\omega }_{\mathrm{s}}\frac{I_{\mathrm{s}}}{I_{\perp }}·\end{array}$$(A.5)

The solution for the harmonic motion is (with A and ϕ arbitrary constants)

$$\begin{array}{c}\hfill {\omega }_x=Asin(\gamma t+\varphi ),\end{array}$$(A.6)

while for ω_y we obtain

$$\begin{array}{c}\hfill {\omega }_y=-\frac{I_{\perp }}{L_{\mathrm{s}}}\frac{\mathrm{d}{\omega }_x}{\mathrm{d}t}=-\frac{I_{\perp }}{I_{\mathrm{s}}{\omega }_{\mathrm{s}}}A\gamma cos(\gamma t+\varphi )=-Acos(\gamma t+\varphi ).\end{array}$$(A.7)

Integrating ω_x and ω_y, we obtain

$$\eta =-\frac{A}{\gamma }\cos (\gamma t+\varphi )+{\theta }_{x_0},$$(A.8)

$$\xi =-\frac{A}{\gamma }\sin (\gamma t+\varphi )+{\theta }_{y_0}.$$(A.9)

In the small-angle approximation we impose that A/γ ≪ 1. These equations reveal that the spin axis rotates around a fixed direction in space. If that direction is along the z-axis, then θ_x0 = θ_y0 = 0. Assuming the initial conditions η(t = 0) = θ₀ and ξ(t = 0)=0, and assuming that A/γ = θ₀, we obtain

$$\eta ={\theta }_0\cos (\gamma t),$$(A.10)

$$\xi ={\theta }_0\sin (\gamma t).$$(A.11)

The last equations describe the torque-free precession of the spin axis that rotates in space at a fixed angle θ₀ with respect to the z-axis with a frequency of the precession motion given by γ = ω_sI_s/I_⊥ (Kleppner & Kolenkow 1973).

Considering a thin disk, we obtain I_s = 2I_⊥ and so γ = 2ω_s, which means that the disk wobbles twice as fast as it spins.

Finally, the apparent rate of a thin-disk precession to an observer on the rigid body is

$$\begin{array}{c}\hfill {\gamma }^{\prime }=\gamma -{\omega }_{\mathrm{s}}={\omega }_{\mathrm{s}}\left(\frac{I_{\mathrm{s}}-I_{\perp }}{I_{\perp }}\right)\sim {\omega }_{\mathrm{s}}.\end{array}$$(A.12)

Appendix B: 3D Jones matrix calculation

We detail in this appendix the 3D Jones matrix description of a rotating HWP. The Jones matrices used here are those that describe the rotation about the two axes that are orthogonal to the propagation direction (with their inverse transformations):

$$\begin{array}{cc}\hfill J_{\mathrm{rotX}}(\eta )& =\left[\begin{array}{ccc}1& 0& 0\\ 0& cos(\eta )& -sin(\eta )\\ 0& sin(\eta )& cos(\eta )\\ \end{array}\right],\hfill \\ \hfill J_{\mathrm{rotY}}(\xi )& =\left[\begin{array}{ccc}cos(\xi )& 0& sin(\xi )\\ 0& 1& 0\\ -sin(\xi )& 0& cos(\xi )\\ \end{array}\right].\hfill \end{array}$$(B.1)

The angles η and ξ are the precession angles η and ξ derived in Sect. 2 (Eqs. (1) and (2)).

The Jones matrix for an HWP with the fast axis at angle θ with respect to the horizontal axis is

$$\begin{array}{c}\hfill J_{\mathrm{HWP}}(\theta )=\left[\begin{array}{ccc}cos(2\theta )& sin(2\theta )& 0\\ sin(2\theta )& -cos(2\theta )& 0\\ 0& 0& 1\\ \end{array}\right].\end{array}$$(B.2)

In the end, the matrix for a linear polarizer that transmits the horizontal component of a light beam is

$$\begin{array}{c}\hfill J_{\mathrm{pol}}=\left[\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ \end{array}\right].\end{array}$$(B.3)

B.1. No-wobble case

By combining these matrices, we obtain a general formula for the Jones vector at the on-axis detector of a polarimeter with modulating element precesses:

$$\begin{array}{c}\hfill i_{\mathrm{out}}=J_{\mathrm{pol}}·J_{\mathrm{RotY}}^{-1}(\xi )·J_{\mathrm{RotX}}^{-1}(\eta )·J_{\mathrm{HWP}}(\theta )·J_{\mathrm{RotX}}(\eta )·J_{\mathrm{RotY}}(\xi )·i_{\mathrm{in}}.\end{array}$$(B.4)

These calculations were realized with the Python package SymPy¹, which allows symbolic computations with matrices and vectors.

For the ideal case with η = ξ = 0°, when no precession occurs, we obtain the outgoing Jones vector for an ideal polarimeter:

$$\begin{array}{c}\hfill i_{\mathrm{out}}=J_{\mathrm{pol}}·J_{\mathrm{HWP}}(\theta )·i_{\mathrm{in}}=\left[\begin{array}{c}{E}_{x}cos\left(2\theta \right)+E_{y}sin\left(2\theta \right)\\ 0\\ 0\end{array}\right].\end{array}$$(B.5)

So that the intensity is

$$\begin{array}{c}\hfill I=i_{{\mathrm{out}}_x}^2+i_{{\mathrm{out}}_y}^2+i_{{\mathrm{out}}_z}^2={(E_{x}cos\left(2\theta \right)+E_{y}sin\left(2\theta \right))}^2,\end{array}$$(B.6)

where we recognize the modulation terms, at four times the frequency of the HWP rotation, within the expressions cos²(2θ), sin²(2θ) and cos(2θ)sin(2θ). Passing through the Stokes parameters (T, Q, U), we retrieve the intensity at the same detector for an ideal polarimeter:

$$\begin{array}{cc}\hfill I=& (E_{x}cos\left(2\theta \right)+E_{y}sin\left(2\theta \right))(E_x^{\ast }cos\left(2\theta \right)+E_y^{\ast }sin\left(2\theta \right))\hfill \\ \hfill =& E_x^2{cos}^2(2\theta )+E_y^2{sin}^2(2\theta )\hfill \\ \hfill & +sin(2\theta )cos(2\theta )(E_x{E}_y^{\ast }+E_y{E}_x^{\ast })\hfill \\ \hfill =& \frac{T+Q}{2}{cos}^2(2\theta )+\frac{T-Q}{2}{sin}^2(2\theta )+sin(2\theta )cos(2\theta )\frac{U}{2}\hfill \\ \hfill =& \frac{1}{2}[T+Qcos(4\theta )+Usin(4\theta )],\hfill \end{array}$$(B.7)

where the Stokes components T, Q, U are then defined as usual:

$$\begin{array}{cc}\hfill T& =⟨E_x{E}_x^{\ast }⟩+⟨E_y{E}_y^{\ast }⟩Q=⟨E_x{E}_x^{\ast }⟩-⟨E_y{E}_y^{\ast }⟩\hfill \\ \hfill U& =⟨E_x{E}_y^{\ast }⟩+⟨E_y{E}_x^{\ast }⟩.\hfill \end{array}$$(B.8)

B.2. Wobbled case

For a general polarimeter with a wobbling HWP, we obtain a general formula for the Jones vector at the on-axis detector:

$$\begin{array}{c}\hfill i_{\mathrm{out}}=J_{\mathrm{pol}}·J_{\mathrm{RotY}}^{-1}(\xi )·J_{\mathrm{RotX}}^{-1}(\eta )·J_{\mathrm{HWP}}(\theta )·J_{\mathrm{RotX}}(\eta )·J_{\mathrm{RotY}}(\xi )·i_{\mathrm{in}}.\end{array}$$(B.9)

For the general case with η, ξ ≠ 0°, we find the intensity at the detector:

$$\begin{array}{c}\hfill I=i_{{\mathrm{out}}_x}^2+i_{{\mathrm{out}}_y}^2+i_{{\mathrm{out}}_z}^2={(E_{x}g-E_{y}f)}^2,\end{array}$$(B.10)

where we define the modulating functions g(ξ, θ) and f(η, ξ, θ):

$$\begin{array}{cc}\hfill g=& {sin}^2\left(\xi \right)+{cos}^2\left(\xi \right)cos\left(2\theta \right),\hfill \\ \hfill f=& 2sin\left(\theta \right)cos\left(\xi \right)(sin\left(\eta \right)sin\left(\xi \right)sin\left(\theta \right)\hfill \\ \hfill & -cos\left(\eta \right)cos\left(\theta \right)).\hfill \end{array}$$(B.11)

We can write the outgoing intensity through the Stokes parameters as follows:

$$\begin{array}{cc}\hfill I& =(E_{x}g-E_{y}f)(E_x^{\ast }g-E_y^{\ast }f)\hfill \\ \hfill & =E_x^2{g}^2+E_y^2{f}^2-E_x{E}_y^{\ast }gf-E_y{E}_x^{\ast }gf\hfill \\ \hfill & =\frac{T+Q}{2}{g}^2+\frac{T-Q}{2}{f}^2-Ugf\hfill \\ \hfill & =\frac{1}{2}(g^2+f^2)T+\frac{1}{2}(g^2-f^2)Q+gfU.\hfill \end{array}$$(B.12)

The last equation gives the intensity at the on-axis detector of a Stokes polarimeter with its modulating element describing a torque-free precession. It is not merely a function of the HWP orientation about the ẑ axis, but also depends on the displacements from the x̂- and ŷ-axes due to the precession.

Appendix C: 3D Müller matrix calculation

The Müller formalism is required to propagate partially polarized light (like that of the CBM), therefore we need to calculate Eqs. (B.7) and (B.12) using the Müller formalism. 3D Müller matrices M_ij are related to Jones matrices (Samim et al. 2016) by

$$\begin{array}{c}\hfill M_{\mathit{ij}}=tr({\sigma }_i·J·{\sigma }_j·J^{†}),\end{array}$$(C.1)

where J is the associated Jones matrix, and σ_n (n = [0, …, 8]) are the Gell-Mann matrices. Eq. (C.1) is valid if and only if we use trace-normalized Gell-Mann matrices (Gell-Mann 1962; Sheppard et al. 2016), which are defined as follows:

$$\begin{array}{cc}\hfill {\sigma }_0& =\frac{1}{\sqrt{3}}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],{\sigma }_1=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0& 1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right],{\sigma }_2=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0& 0& 1\\ 0& 0& 0\\ 1& 0& 0\end{array}\right],\hfill \\ \hfill {\sigma }_3& =\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0& -i& 0\\ i& 0& 0\\ 0& 0& 0\end{array}\right],{\sigma }_4=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}1& 0& 0\\ 0& -1& 0\\ 0& 0& 0\end{array}\right],{\sigma }_5=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 1\\ 0& 1& 0\end{array}\right],\hfill \\ \hfill {\sigma }_6& =\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0& 0& i\\ 0& 0& 0\\ -i& 0& 0\end{array}\right],{\sigma }_7=\frac{1}{\sqrt{2}}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& -i\\ 0& i& 0\end{array}\right],{\sigma }_8=\frac{1}{\sqrt{6}}\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& -2\end{array}\right].\hfill \end{array}$$(C.2)

Using the common polarization matrix

$$\begin{array}{c}\hfill \mathbf{P}=⟨\mathbf{E}{\mathbf{E}}^{\ast }⟩=\left[\begin{array}{ccc}⟨E_x{E}_x^{\ast }⟩& ⟨E_x{E}_y^{\ast }⟩& ⟨E_x{E}_z^{\ast }⟩\\ ⟨E_y{E}_x^{\ast }⟩& ⟨E_y{E}_y^{\ast }⟩& ⟨E_y{E}_z^{\ast }⟩\\ ⟨E_z{E}_x^{\ast }⟩& ⟨E_z{E}_y^{\ast }⟩& ⟨E_z{E}_z^{\ast }⟩\end{array}\right],\end{array}$$(C.3)

we can define the Stokes vector in the 3D formalism:

$$\begin{array}{c}\hfill s=\left[\begin{array}{c}{\mathrm{\Delta }}_0\\ {\mathrm{\Delta }}_1\\ {\mathrm{\Delta }}_2\\ {\mathrm{\Delta }}_3\\ {\mathrm{\Delta }}_4\\ {\mathrm{\Delta }}_5\\ {\mathrm{\Delta }}_6\\ {\mathrm{\Delta }}_7\\ {\mathrm{\Delta }}_8\end{array}\right]=\left[\begin{array}{c}\frac{1}{\sqrt{3}}(⟨E_x{E}_x^{\ast }⟩+⟨E_y{E}_y^{\ast }⟩+⟨E_z{E}_z^{\ast }⟩)\\ \frac{1}{\sqrt{2}}(⟨E_x{E}_y^{\ast }⟩+⟨E_y{E}_x^{\ast }⟩)\\ \frac{1}{\sqrt{2}}(⟨E_z{E}_x^{\ast }⟩+⟨E_x{E}_z^{\ast }⟩)\\ \frac{i}{\sqrt{2}}(⟨E_x{E}_y^{\ast }⟩-⟨E_y{E}_x^{\ast }⟩)\\ \frac{1}{\sqrt{2}}(⟨E_x{E}_x^{\ast }⟩-⟨E_y{E}_y^{\ast }⟩)\\ \frac{1}{\sqrt{2}}(⟨E_y{E}_z^{\ast }⟩+⟨E_z{E}_y^{\ast }⟩)\\ \frac{i}{\sqrt{2}}(⟨E_z{E}_x^{\ast }⟩-⟨E_x{E}_z^{\ast }⟩)\\ \frac{i}{\sqrt{2}}(⟨E_y{E}_z^{\ast }⟩-⟨E_z{E}_y^{\ast }⟩)\\ \frac{1}{\sqrt{6}}(⟨E_x{E}_x^{\ast }⟩+⟨E_y{E}_y^{\ast }⟩-2⟨E_z{E}_z^{\ast }⟩)\end{array}\right].\end{array}$$(C.4)

The conventional 2D Stokes parameters are related to the 3D Stokes parameters (optical ordering) by

$$\begin{array}{c}\hfill T=\sqrt{\frac{2}{3}}({\mathrm{\Delta }}_0+\frac{1}{\sqrt{2}}{\mathrm{\Delta }}_8),Q={\mathrm{\Delta }}_4,U={\mathrm{\Delta }}_1,V={\mathrm{\Delta }}_3.\end{array}$$(C.5)

C.1. No-wobble case

From Eq. (C.1) we can find the analogous 3D Müller matrices for the Jones HWP matrix (Eq. (B.2)) and for the Jones polarizer matrix (Eq. (B.3)),

$$\begin{array}{cc}& M_{\mathrm{polx}}=\left[\begin{array}{ccccccccc}\frac{1}{3}& 0& 0& 0& \frac{\sqrt{6}}{6}& 0& 0& 0& \frac{\sqrt{2}}{6}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{\sqrt{6}}{6}& 0& 0& 0& \frac{1}{2}& 0& 0& 0& \frac{\sqrt{3}}{6}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{\sqrt{2}}{6}& 0& 0& 0& \frac{\sqrt{3}}{6}& 0& 0& 0& \frac{1}{6}\end{array}\right]\hfill \\ \hfill \end{array}$$(C.6)

$$\begin{array}{}{M}_{\mathrm{HWP}}(\theta )=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& -cos\left(4\theta \right)& 0& 0& sin\left(4\theta \right)& 0& 0& 0& 0\\ 0& 0& cos\left(2\theta \right)& 0& 0& sin\left(2\theta \right)& 0& 0& 0\\ 0& 0& 0& -1& 0& 0& 0& 0& 0\\ 0& sin\left(4\theta \right)& 0& 0& cos\left(4\theta \right)& 0& 0& 0& 0\\ 0& 0& sin\left(2\theta \right)& 0& 0& -cos\left(2\theta \right)& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& cos\left(2\theta \right)& -sin\left(2\theta \right)& 0\\ 0& 0& 0& 0& 0& 0& -sin\left(2\theta \right)& -cos\left(2\theta \right)& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1\end{array}\right]\hfill \end{array}$$(C.7)

and by combining the previous matrices, we can determine the 3D polarimeter Stokes vector (η = ξ = 0°):

$$\begin{array}{cc}\hfill s_{\mathrm{out}}& =M_{\mathrm{polx}}·M_{\mathrm{HWP}}(\theta )·s_{\mathrm{in}}=M_{\mathrm{polx}}·M_{\mathrm{HWP}}(\theta )·\left[\begin{array}{c}{\mathrm{\Delta }}_0\\ {\mathrm{\Delta }}_1\\ {\mathrm{\Delta }}_2\\ {\mathrm{\Delta }}_3\\ {\mathrm{\Delta }}_4\\ {\mathrm{\Delta }}_5\\ {\mathrm{\Delta }}_6\\ {\mathrm{\Delta }}_7\\ {\mathrm{\Delta }}_8\end{array}\right]\hfill \\ \hfill & =\left[\begin{array}{c}\frac{1}{3}{\mathrm{\Delta }}_0+\frac{\sqrt{6}}{6}sin\left(4\theta \right){\mathrm{\Delta }}_1+\frac{\sqrt{6}}{6}cos\left(4\theta \right){\mathrm{\Delta }}_4+\frac{\sqrt{2}}{6}{\mathrm{\Delta }}_8\\ 0\\ 0\\ 0\\ \frac{\sqrt{6}}{6}{\mathrm{\Delta }}_0+\frac{1}{2}sin\left(4\theta \right){\mathrm{\Delta }}_1+\frac{1}{2}cos\left(4\theta \right){\mathrm{\Delta }}_4+\frac{\sqrt{3}}{6}{\mathrm{\Delta }}_8\\ 0\\ 0\\ 0\\ \frac{\sqrt{2}}{6}{\mathrm{\Delta }}_0+\frac{\sqrt{3}}{6}sin\left(4\theta \right){\mathrm{\Delta }}_1+\frac{\sqrt{3}}{6}cos\left(4\theta \right){\mathrm{\Delta }}_4+\frac{1}{6}{\mathrm{\Delta }}_8\end{array}\right].\hfill \end{array}$$(C.8)

Through Eq. (C.5), we can derive the intensity:

$$\begin{array}{c}\hfill I=\frac{T}{2}+\frac{Q}{2}cos\left(4\theta \right)+\frac{U}{2}sin\left(4\theta \right).\end{array}$$(C.9)

Equations (C.9) and (B.7) give the same result.

C.2. Wobbled case

For a wobbling HWP we need to calculate the Müller rotation matrices (Eq. (B.1)) from Eq. (C.1):

$$\begin{array}{c}\hfill M_{\mathrm{RotX}}(\eta )=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& cos\left(\eta \right)& 0& 0& 0& sin\left(\eta \right)& 0& 0& 0\\ 0& 0& cos\left(2\eta \right)& 0& -\frac{1}{2}sin\left(2\eta \right)& 0& 0& 0& -\frac{\sqrt{3}}{2}sin\left(2\eta \right)\\ 0& 0& 0& cos\left(\eta \right)& 0& 0& 0& -sin\left(\eta \right)& 0\\ 0& 0& \frac{1}{2}sin\left(2\eta \right)& 0& \frac{1}{2}{cos}^2\left(\eta \right)+\frac{1}{2}& 0& 0& 0& -\frac{\sqrt{3}}{2}{sin}^2\left(\eta \right)\\ 0& -sin\left(\eta \right)& 0& 0& 0& cos\left(\eta \right)& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& sin\left(\eta \right)& 0& 0& 0& cos\left(\eta \right)& 0\\ 0& 0& \frac{\sqrt{3}}{2}sin\left(2\eta \right)& 0& -\frac{\sqrt{3}}{2}{sin}^2\left(\eta \right)& 0& 0& 0& -\frac{3}{2}{sin}^2\left(\eta \right)+1\end{array}\right]\end{array}$$(C.10)

$$\begin{array}{c}\hfill M_{\mathrm{RotY}}(\xi )=\left[\begin{array}{ccccccccc}1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& cos\left(\xi \right)& -sin\left(\xi \right)& 0& 0& 0& 0& 0& 0\\ 0& sin\left(\xi \right)& cos\left(\xi \right)& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& cos\left(\xi \right)& 0& 0& -sin\left(\xi \right)& 0& 0\\ 0& 0& 0& 0& \frac{1}{2}{cos}^2\left(\xi \right)+\frac{1}{2}& \frac{1}{2}sin\left(2\xi \right)& 0& 0& \frac{\sqrt{3}}{2}{sin}^2\left(\xi \right)\\ 0& 0& 0& 0& -\frac{1}{2}sin\left(2\xi \right)& cos\left(2\xi \right)& 0& 0& \frac{\sqrt{3}}{2}sin\left(2\xi \right)\\ 0& 0& 0& sin\left(\xi \right)& 0& 0& cos\left(\xi \right)& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& \frac{\sqrt{3}}{2}{sin}^2\left(\xi \right)& -\frac{\sqrt{3}}{2}sin\left(2\xi \right)& 0& 0& -\frac{3}{2}{sin}^2\left(\xi \right)+1.\end{array}\right].\end{array}$$(C.11)

Finally, we can calculate the 3D Müller matrix for a Stokes polarimeter with a wobbling HWP:

$$\begin{array}{c}\hfill M_{{\text{SP}}_{\mathrm{wob}}}=M_{\mathrm{polx}}·M_{\mathrm{RotY}}^{-1}(\xi )·M_{\mathrm{RotX}}^{-1}(\eta )·M_{\mathrm{HWP}}(\theta )·M_{\mathrm{RotX}}(\eta )·M_{\mathrm{RotY}}(\xi )=\left[\begin{array}{ccccccccc}\frac{1}{3}& m_{01}& m_{02}& 0& m_{04}& m_{05}& 0& 0& m_{08}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{\sqrt{6}}{6}& \frac{\sqrt{6}}{2}{m}_{01}& \frac{\sqrt{6}}{2}{m}_{02}& 0& \frac{\sqrt{6}}{2}{m}_{04}& \frac{\sqrt{6}}{2}{m}_{05}& 0& 0& \frac{\sqrt{6}}{2}{m}_{08}\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0\\ \frac{\sqrt{2}}{6}& \frac{\sqrt{2}}{2}{m}_{01}& \frac{\sqrt{2}}{2}{m}_{02}& 0& \frac{\sqrt{2}}{2}{m}_{04}& \frac{\sqrt{2}}{2}{m}_{05}& 0& 0& \frac{\sqrt{2}}{2}{m}_{08}\end{array}\right]\end{array}$$(C.12)

For the sake of clarity, we define the following equations that we use from now on:

$$A={\sin }^2\left(\eta \right){\sin }^2\left(\xi \right)+{\sin }^2\left(\xi \right)-1$$(C.13)

$$B=\sin \left(\eta \right)\sin \left(2\xi \right)$$(C.14)

$$C=\sin \left(2\eta \right){\sin }^2\left(\xi \right)$$(C.15)

$$D=\sin \left(2\xi \right)\cos \left(\eta \right)$$(C.16)

$$E=3{\sin }^2\left(\eta \right){\sin }^2\left(\xi \right)-3{\sin }^2\left(\xi \right)+1$$(C.17)

$$\begin{array}{cc}\hfill m_{01}=& \frac{\sqrt{6}}{12}(Esin\left(2\eta \right)sin\left(\xi \right)-2(Asin\left(4\theta \right)\hfill \\ \hfill & +Bcos\left(4\theta \right))cos\left(\eta \right)cos\left(\xi \right)\hfill \\ \hfill & -(-Acos\left(4\theta \right)+Bsin\left(4\theta \right))sin\left(2\eta \right)sin\left(\xi \right)\hfill \\ \hfill & -2(Csin\left(2\theta \right)+Dcos\left(2\theta \right))sin\left(\eta \right)cos\left(\xi \right)\hfill \\ \hfill & -2(Ccos\left(2\theta \right)-Dsin\left(2\theta \right))sin\left(\xi \right)cos\left(2\eta \right))\hfill \end{array}$$(C.18)

$$\begin{array}{cc}\hfill m_{02}=& \frac{\sqrt{6}}{12}(Esin\left(2\xi \right){cos}^2\left(\eta \right)+2(Asin\left(4\theta \right)\hfill \\ \hfill & +Bcos\left(4\theta \right))sin\left(\eta \right)cos\left(2\xi \right)\hfill \\ \hfill & -\frac{1}{2}(Acos\left(4\theta \right)-Bsin\left(4\theta \right))(3-cos\left(2\eta \right))sin\left(2\xi \right)\hfill \\ \hfill & -2(Csin\left(2\theta \right)+Dcos\left(2\theta \right))cos\left(\eta \right)cos\left(2\xi \right)\hfill \\ \hfill & +(Ccos\left(2\theta \right)-Dsin\left(2\theta \right))sin\left(2\eta \right)sin\left(2\xi \right))\hfill \end{array}$$(C.19)

$$\begin{array}{cc}\hfill m_{04}=& \frac{\sqrt{6}}{12}(E({cos}^2\left(\eta \right){cos}^2\left(\xi \right)-2{cos}^2\left(\eta \right)+1)-(Asin\left(4\theta \right)\hfill \\ \hfill & +Bcos\left(4\theta \right))sin\left(\eta \right)sin\left(2\xi \right)-(Acos\left(4\theta \right)\hfill \\ \hfill & -Bsin\left(4\theta \right))(2{cos}^2\left(\eta \right)+2{cos}^2\left(\xi \right)-{cos}^2\left(\eta \right){cos}^2\left(\xi \right)-1)\hfill \\ \hfill & +(Csin\left(2\theta \right)+Dcos\left(2\theta \right))sin\left(2\xi \right)cos\left(\eta \right)+\frac{1}{2}(Ccos\left(2\theta \right)\hfill \\ \hfill & -Dsin\left(2\theta \right))(cos\left(2\xi \right)-3)sin\left(2\eta \right))\hfill \end{array}$$(C.20)

$$\begin{array}{cc}\hfill m_{05}=& \frac{\sqrt{6}}{12}(-Esin\left(2\eta \right)cos\left(\xi \right)-2(Asin\left(4\theta \right)\hfill \\ \hfill & +Bcos\left(4\theta \right))sin\left(\xi \right)cos\left(\eta \right)+(-Acos\left(4\theta \right)\hfill \\ \hfill & +Bsin\left(4\theta \right))sin\left(2\eta \right)cos\left(\xi \right)-2(Csin\left(2\theta \right)\hfill \\ \hfill & +Dcos\left(2\theta \right))sin\left(\eta \right)sin\left(\xi \right)+2(Ccos\left(2\theta \right)\hfill \\ \hfill & -Dsin\left(2\theta \right))cos\left(2\eta \right)cos\left(\xi \right))\hfill \end{array}$$(C.21)