Systematic effects induced by Half Wave Plate precession into Cosmic Microwave Background polarization measurements

TL;DR

This paper investigates the systematic effects caused by HWP precession in CMB polarization measurements, providing analytical formulas, simulations, and implications for B-mode detection in satellite experiments.

Contribution

It introduces an analytical model for HWP precession effects on polarization data and assesses their impact on B-mode measurements in future satellite missions.

Findings

HWP precession induces non-negligible systematic signals in polarization data.

The analytical formula relates precession parameters to Stokes parameters.

Systematic effects must be carefully controlled in future B-mode experiments.

Abstract

The most accessible method to measure polarization features of the CMB radiation is by means of a Stokes Polarimeter based on the rotation of an Half Wave Plate. The current observational cosmology is starting to be limited by the presence of systematic effects. The Stokes polarimeter with a rotating Half Wave Plate (HWP) has the advantage of mitigating a long list of potential systematics, by modulation of the linearly polarized component of the radiation, but the presence of the rotating HWP can by itself introduce new systematic effects, which must be under control, representing one of the most critical part in the design of a B-Modes experiment. In this paper we present, simulate and analyse the spurious signal arising from the precession of a rotating HWP. We first find an analytical formula for the impact of the systematic effect induced by the HWP precession on the propagating radiation, using the 3D generalization of the Muller formalism. We then perform several numerical simulations, showing the effect induced on the Stokes parameters by this systematic. We also derive and discuss the impact into B-modes measured by a satellite experiment. We find the analytical formula for the Stokes parameters from a Stokes polarimeter where the HWP follows a precessional motion with an angle $\theta_0$. We show the result depending on the HWP inertia tensor, spinning speed and on $\theta_0$. The result of numerical simulations is reported as a simple timeline of the electric fields. Finally, assuming to observe all the sky with a satellite mission, we analyze the effect on B-modes measurements. The effect is not negligible giving the current B-modes experiments sensitivity, therefore it is a systematic which needs to be carefully considered for future experiments.

Tables

Table 1. Table 1: Parameters for the scanning strategies adopted in the simulation pipeline.

Scanning	$\alpha \left[deg\right]$	$\beta \left[deg\right]$	${\omega }_2\left[deg/min\right]$	${\omega }_3\left[deg/min\right]$
Planck like	7.5	85.0	0.00139	360.0
WMAP like	22.5	70.5	6.0	167.0
COrE like	30.0	65.0	0.0625	180.0
LiteBIRD like	45.0	50.0	3.8709	36.0

Table 2. Table 2: RMS values in   µ ​ K times absent microkelvin \text{\,}\mathrm{\SIUnitSymbolMicro K} for intensity and polarization for different scan strategies. Maps correspond on the input/output difference maps for each case.

	Planck			WMAP			COrE			LiteBIRD
Frequency [Hz]	T[$\mu K$]	Q[$\mu K$]	U[$\mu K$]	T[$\mu K$]	Q[$\mu K$]	U[$\mu K$]	T[$\mu K$]	Q[$\mu K$]	U[$\mu K$]	T[$\mu K$]	Q[$\mu K$]	U[$\mu K$]
1.0	0.600	0.032	0.032	0.601	0.071	0.071	0.600	0.071	0.071	0.600	0.071	0.071
0.1	0.600	0.050	0.050	0.600	0.072	0.072	0.600	0.088	0.087	0.600	0.076	0.076

Table 3. Table 3: RMS values in   µ ​ K times absent microkelvin \text{\,}\mathrm{\SIUnitSymbolMicro K} for intensity and polarization residual maps considering a COrE -like scan strategy with the HWP spin speed at 0.1Hz and with different precession angles.

${\theta }_0$ [∘]	T[$\mu K$]	Q[$\mu K$]	U[$\mu K$]
$0.5$	0.150	0.533	0.531
$1$	0.600	0.540	0.537
$1.5$	1.350	0.569	0.569
$2$	2.401	0.640	0.640

Table 4. Table 4: RMS values, in   µ ​ K times absent microkelvin \text{\,}\mathrm{\SIUnitSymbolMicro K} , of the difference between the maps with and without the effect of the HWP wobbling. The RMS is calculated for the maps from the two detectors, and for the map from the two detectors combined. The RMS of the residual is reduced by a factor ∼ 2 similar-to absent 2 \sim\sqrt{2} combining two detectors. Modulation parameters: f s = 1 ​ H ​ z subscript 𝑓 𝑠 1 𝐻 𝑧 f_{s}=1Hz , θ 0 = 1 ∘ subscript 𝜃 0 superscript 1 \theta_{0}=1^{\circ} , I ⟂ / I s = 0.508 subscript 𝐼 perpendicular-to subscript 𝐼 𝑠 0.508 I_{\perp}/I_{s}=0.508 for Planck  and COrE , I ⟂ / I s = 0.514 subscript 𝐼 perpendicular-to subscript 𝐼 𝑠 0.514 I_{\perp}/I_{s}=0.514 for LiteBIRD and I ⟂ / I s = 0.510 subscript 𝐼 perpendicular-to subscript 𝐼 𝑠 0.510 I_{\perp}/I_{s}=0.510 for WMAP.

	$\beta$ [∘]	Q[$\mathrm{nK}$]	U[$\mathrm{nK}$]
Planck	$85.0$	31.9	32.0
	$86.0$	31.8	31.8
	combination	22.6	22.6
WMAP	$70.5$	70.7	70.8
	$71.5$	71.8	71.8
	combination	49.8	49.9
COrE	$65.0$	71.0	70.7
	$66.0$	72.6	72.6
	combination	50.9	50.5
LiteBIRD	$50.0$	70.9	70.9
	$51.0$	70.1	70.1
	combination	50.1	50.1

Equations

$$\eta$$

$$\xi$$

$${\gamma }^{\prime }=\gamma -{\omega }_s={\omega }_s\left(\frac{I_s-I_{\perp }}{I_{\perp }}\right)\sim {\omega }_s$$

$$i_{out}=J_{pol}\cdot J_{RotY}^{-1}(\xi )\cdot J_{RotX}^{-1}(\eta )\cdot J_{HWP}(\theta )\cdot J_{RotX}(\eta )\cdot J_{RotY}(\xi )\cdot i_{in}$$

$$I=\frac{1}{2}\left[T+Q\cos (4\theta )+U\sin (4\theta )\right].$$

$$I=i_{ou{t}_x}^2+i_{ou{t}_y}^2+i_{ou{t}_z}^2={\left(E_{x}g-E_{y}f\right)}^2$$

$$g$$

$$f$$

$$I=\frac{1}{2}\left(g^2+f^2\right)T+\frac{1}{2}\left(g^2-f^2\right)Q+gfU$$

$$M_{ij}=tr({\sigma }_i\cdot J\cdot {\sigma }_j\cdot J^{†})$$

$$M_{{\text{SP}}_{\text{wob}}}(\eta ,\xi ,\theta )=M_{pol}\cdot M_{RotY}^{-1}(\xi )\cdot M_{RotX}^{-1}(\eta )\cdot M_{HWP}(\theta )\cdot M_{RotX}(\eta )\cdot M_{RotY}(\xi )$$

$$s_{out}=M_{{\text{SP}}_{\text{wob}}}(\eta ,\xi ,\theta )\cdot s_{in}$$

$$s_{in}^T=({\Delta }_0,{\Delta }_1,0,0,{\Delta }_4,0,0,0,\frac{1}{\sqrt{2}}{\Delta }_0)$$

$$I=\frac{1}{\sqrt{6}}\left({\Delta }_0+3\left(m_{01}{\Delta }_1+m_{04}{\Delta }_4+m_{08}{\Delta }_8\right)\right)=\left(\frac{1}{3}+\frac{\sqrt{2}}{2}{m}_{08}\right)T+\frac{\sqrt{6}}{2}{m}_{01}U+\frac{\sqrt{6}}{2}{m}_{04}Q$$

$$I_s$$

$$I_{\perp }$$

$$\mathbf{I}=\left[\begin{array}{ccc}{I}_{\perp }& 0& 0\\ 0& I_{\perp }& 0\\ 0& 0& I_s\end{array}\right].$$

$$I_{\perp }/I_s=0.5+\frac{1}{6}{\left(\frac{h}{R}\right)}^2.$$

$$L_x$$

$$L_y$$

$$L_x$$

$$L_y$$

$$I_{\perp }\frac{d^2\eta }{d{t}^2}+L_s\frac{d\xi }{dt}$$

$$I_{\perp }\frac{d^2\xi }{d{t}^2}-L_s\frac{d\eta }{dt}$$

$$I_{\perp }\frac{d{\omega }_x}{dt}+L_s{\omega }_y$$

$$I_{\perp }\frac{d{\omega }_y}{dt}-L_s{\omega }_x$$

$$\frac{d^2{\omega }_x}{d{t}^2}+{\gamma }^2{\omega }_x=0\mathrm{with}\gamma =\frac{L_s}{I_{\perp }}={\omega }_s\frac{I_s}{I_{\perp }}.$$

$${\omega }_x=A\sin (\gamma t+\varphi ).$$

$${\omega }_y=-\frac{I_{\perp }}{L_s}\frac{d{\omega }_x}{dt}=-\frac{I_{\perp }}{I_s{\omega }_s}A\gamma \cos (\gamma t+\varphi )=-A\cos (\gamma t+\varphi ).$$

$$\eta$$

$$\xi$$

$$\eta$$

$$\xi$$

$${\gamma }^{\prime }=\gamma -{\omega }_s={\omega }_s(\frac{I_s-I_{\perp }}{I_{\perp }})\sim {\omega }_s$$

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

Full text

*Aims. *The current observational cosmology is starting to be limited by the presence of systematic effects. The Stokes polarimeter with a rotating Half Wave Plate (HWP) has the advantage of mitigating a long list of potential systematics, by modulation of the linearly polarized component of the radiation, but the presence of the rotating HWP can by itself introduce new systematic effects, which must be under control, representing one of the most critical part in the design of a B-Modes experiment. It is therefore mandatory to take into account all the systematic effects the instrumentation can induce. In this paper we present, simulate and analyse the spurious signal arising from the precession of a rotating HWP.

*Methods. *We first find an analytical formula for the impact of the systematic effect induced by the HWP precession on the propagating radiation, using the 3D generalization of the Müller formalism. We then perform several numerical simulations, showing the effect induced on the Stokes parameters by this systematic. We also derive and discuss the impact into B-modes measured by a satellite experiment.

*Results. *We find the analytical formula for the Stokes parameters from a Stokes polarimeter where the HWP follows a precessional motion with an angle $\theta_{0}$. We show the result depending on the HWP inertia tensor, spinning speed and on $\theta_{0}$. The result of numerical simulations is reported as a simple timeline of the electric fields. Finally, assuming to observe all the sky with a satellite mission, we analyze the effect on B-modes measurements.

*Conclusions. *The effect is not negligible giving the current B-modes experiments sensitivity, therefore it is a systematic which needs to be carefully considered for future experiments.

1 Introduction

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

Further improve the constraint on the tensor-to-scalar ratio represents a hard challenge for the current and future CMB experiments which must take into account, accurately, gravitational lensing and foreground removal (Errard et al. 2016) in addition to an excellent control of systematic effects (Wallis et al. 2017). Concerning the control over systematics, some experiments are designed with the capability of self-calibrating and thus of removing some systematic effects (Piat et al. 2012). Others experiments, which do not have 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 the Stokes Polarimeter, composed by an 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 by rotating the HWP it is possible to modulate the linearly polarized fraction of the incoming radiation, and to extract the Stokes components (T, Q, U), given a reference frame. The HWP can spin fast ($>1Hz$), see e.g. The LSPE collaboration et al. (2012); Columbro et al. (2019); The EBEX collaboration et al. (2018); Thornton et al. (2016) in order to modulate the signal at high frequency and it can also rotate step by step ($<1Hz$) (Piat et al. 2012; Bryan et al. 2016) depending on the experiment scanning strategy. Systematic effects like T-P leakage, chromatic HWP behaviour, scanning strategy, are already evaluated in Essinger-Hileman et al. (2016); Salatino et al. (2011); Takakura et al. (2017); Salatino et al. (2018), and also measured by Kusaka et al. (2014a).

In this paper we analyse the systematic effect induced in observation of the CMB polarization by the precession of the HWP, in the specific case of a transmissive plate. We assume the HWP free from other systematic effects.

In Sections 2 we review the dynamic of a spinning cylindrical plate, and define the precession rate. In section 3 we provide an analytic study of the effect induced by the HWP which spins and precesses, by using the 3D Jones formalism, for fully polarized radiation, and the 3D generalization of the Müller formalism, for partially polarized radiation. In section 4 by using the equations derived before, we show some results on electric field components produced with numerical simulations. Finally, in section 5 we describe the effect 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 describing a precessing body rotating along one symmetry axis. We introduce the main variables showing their evolution with time. We report here only the main equations, essential for the results shown in the following sections. All the details of the computation are provided in Appendix B.

We approximate the HWP as a cylindrically symmetric rigid body, like a coin, and we define a reference system $\hat{x}$-$\hat{y}$-$\hat{z}$ with the $\hat{x}$-$\hat{y}$ plane coincident with the base of the cylinder, see Figure1. We hypothesize that the HWP has a large spin angular momentum $L_{s}=I_{s}\omega_{s}$ along the symmetry axis, where $I_{s}$ and $\omega_{s}$ are respectively the moment of inertia and the angular velocity. In the unperturbed case $L_{s}$ coincides with the $\hat{z}$ axis.

We now introduce a small angular perturbation ($\theta_{0}$) of the spin momentum $L_{s}$, expanding it in its components along the $\hat{x}$ and $\hat{y}$ axes, $\eta$ and $\xi$. By introducing the precession rate $\gamma$, defined as the rotation velocity of $L_{s}$ around the $\hat{z}$ axis, $\eta$ and $\xi$ are simply defined as:

[TABLE]

Those equations describe the torque-free precession of the spin axis that rotates in space with fixed angle $\theta_{0}$.

The precession rate is given by $\gamma=\omega_{s}I_{s}/I_{\perp}$ (D. Kleppner & R. Kolenkow 1973), having defined $I_{s}$ and $I_{\perp}$ as the inertia momenta respectively parallel and perpendicular to the $\hat{z}$ axis. Assuming the cylindrical plate as a thin disc we can approximate $I_{s}=2I_{\perp}$ and $\gamma=2\omega_{s}$, thus the disc wobbles twice as fast as it spins.

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

[TABLE]

In a torque-free precession we can identify two different rotations accordingly to the reference frame we consider. In the fixed laboratory frame, the angular velocity vector rotates around the fixed $\hat{z}$ axis (where the angular momentum vector lies), tracing the so-called ”space cone”. In the reference frame integral with the rotating body, we can see both angular momentum and angular velocity vector describing a circle around the symmetry axis of the cylinder, tracing the so-called ”body cone” with a precession rate $\gamma^{\prime}$.

3 3D Generalization of Stokes Polarimeter

In literature a description of Stokes Polarimetry by using Müller and Jones matrices is already present, see e.g. Bryan et al. (2010), O’Dea et al. (2007); Chuss et al. (2012). In this section, and in Appendix C and D, we derive the analytic equation both 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 C 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 2x2 matrices. We extend the Jones vector in three dimensions, $i=(E_{x},E_{y},0$), where the $\hat{z}$ axis is the propagation direction of the EM field, and the $\hat{x}$-$\hat{y}$ plane is perpendicular to $\hat{z}$, thus any optical element is described by 3x3 matrices (Sheppard 2011, 2014; Ortega-Quijano & Arce-Diego 2013). Combining such matrices we get a general formula for the Jones vector for the on-axis detector of a precessing polarimeter:

[TABLE]

For the ideal case with $\eta=\xi=0^{\circ}$, when no precession occurs, we get the outgoing intensity from an ideal polarimeter:

[TABLE]

For the general case with $\eta$, $\xi\neq 0^{\circ}$ we find the intensity at the detector:

[TABLE]

where we define the modulating functions $g(\xi,\theta)$ and $f(\eta,\xi,\theta)$:

[TABLE]

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

[TABLE]

valid only if $T^{2}=Q^{2}+U^{2}$. The last equation gives the intensity for 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 $\hat{z}$ axis, it also depends on the displacements from the $\hat{x}$ and $\hat{y}$ 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 EM field intensity outgoing from the Stokes polarimeter. The Müller formalism (fully defined in Appendix D) is required for the case of CMB experiments since we want 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 C), the $9x9$ Müller matrix corresponding to each optical element can be easily obtained from Eq. 41:

[TABLE]

where $\sigma_{n}$ ($n=[0,...,8]$) are the trace-normalized Gell-Mann matrices 42. The Stokes polarimeter Müller matrix becomes:

[TABLE]

So, the output Stokes vector is:

[TABLE]

Each component of the Stokes vector is defined in Eq. 44 and assuming the field entering the polarimeter has $E_{z}=0$,

[TABLE]

Thanks to Eq. 45 which describes the relation between the 3D Stokes components and the classic definition of Stokes parameters we can find the outgoing intensity:

[TABLE]

where $m_{ij}$ are the components of the Stokes polarimeter Müller matrix and are function of to the wobbling angles $\eta,\xi$. Note that if $\eta=\xi=0$ the Eq. 12 gives the well-known formula of the Stokes polarimeter (Eq. 5).

4 Phenomenology

Starting from Eq. 40, we built a simulation to show the effects of the precession for a specific case. We set the spinning frequency at 1Hz and the ratio $I_{\perp}/I_{s}=0.502$. The Jones vector used as input field in the simulation is $i_{in}=(1,0,0)$, corresponding to (1,1,0,0) Stokes vector. Typical outputs of the simulations are reported in Figure 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.

The timelines reported in Figure 2(a) show the fractional residual for different amplitudes of the precession angle $\theta_{0}$. These timelines clearly show several beats with an amplitude depending on the precession angle $\theta_{0}$ in a non-linear way.

In Figure 2(c) we report the power spectra of timelines, for the ideal case (red), and for the precessing cases. The ideal case presents a single line at 4Hz, as expected from a wave-plate spinning at 1Hz. The beats from a precessing HWP produce spurious peaks at multiple and sub-multiple of the spin frequency. The most noteworthy peak in the spectra is the one at lower frequencies as it is the result of the very slow modulation imposed by the precession, that in Figure 2(a) has a period of about 60s, implying a peak at 0.016Hz 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_{\perp}/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 respect to the principal axes are:

[TABLE]

where we are assuming a diagonal inertia tensor:

[TABLE]

The frequency for the precessional motion is directly linked to the HWP regular spin frequency $f_{s}$ and to the ratio $I_{\perp}/I_{s}$ (Eq. 21).

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

[TABLE]

We therefore explore different configurations as shown in Figure 3 where we show the fractional residual with respect to the ideal case. We consider input radiation horizontally polarized, a precession angle $\theta_{0}=1^{\circ}$, spin frequencies $f_{s}=[0.1,0.5,1.0,2.0]Hz$ and $I_{\perp}/I_{s}=[0.501,0.502,0.506,0.513]$. As an example, these values for the $I_{\perp}/I_{s}$ ratio correspond to an HWP with mass $m=1kg$, radius $R=16cm$, and thickness $h\simeq[1.2,1.7,3.0,4.5]\leavevmode\nobreak\ cm$; then $I_{s}$ is fixed to $0.0128$ $kg\cdot m^{2}$. In practice the ratio of the components of the moment of inertia does not depend only on the plate thickness, but also on the structure of the HWP support.

The 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_{\perp}/I_{s}\rightarrow 0.500$, has beats in the Intensity over an extended period while a thicker one has shorter beats. It is clear by looking the Figure 3 from top to bottom. Anyway the value $0.5$ is not possible since it corresponds to a null thickness.

The effect of different spinning speeds is to compress the beats. This is clear by looking the Figure 3 from left to right. The maximum amplitude remains the same because it depends only on $\theta_{0}$ that is fixed to $1^{\circ}$ in this particular simulation.

In Fig. 4 we report the power spectra of the timelines for different ratios $I_{\perp}/I_{s}$. The spectra exhibit the effect discussed above, showing that the beats frequency moves to lower values as the ratio $I_{\perp}/I_{s}$ is reduced towards the minimum value of 0.5 (for $I_{\perp}/I_{s}=0.513$ the period of the beats is 10s, corresponding to a peak in the spectrum at 0.1Hz).

4.2 Minor effect

As reported e.g. in Honggang Gu & Liu (2018), a tilted HWP also changes its transmission properties due to the different path of radiation inside the plate with respect to 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 Honggang Gu & Liu (2018), for a source at 347 nm, with 5 degrees tilt. This effect is multiplicative with the ratio wavelength/thickness, which is much higher in the case of plates used in the optical bandwidth with respect to the case of the millimeter wavelength. In instrumentation devoted to observations of the CMB polarisation, the plate thickness is of the same order of magnitude as the wavelength: 3.1mm POLARBEAR-2 (Charles A. Hill 2016), 3.05mm SPIDER (Sean A. Bryan 2010), 1.62mm EBEX (The EBEX collaboration et al. 2018), 3.2mm QUBIC (Aumont et al. 2016), 3mm LSPE-SWIPE (The LSPE collaboration et al. 2012), 3.05mm ABS (Kusaka et al. 2014b). The impact of the effect for millimeter astronomy can be examined through the variation between the input Stokes vector and the output one, and it is estimated to be orders of magnitude smaller than the effect due to electromagnetic field projection analyzed in the rest of this paper.

5 Full-Sky Simulations

In order to test the impact of the HWP precessional motion on CMB observations we build an algorithm able to generate a realistic satellite scanning strategy in presence of spinning HWP, producing data timelines. We complete this software with a map-making algorithm which collapses data timelines 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), namely:

• •

Earth revolution velocity $\omega_{1}$,

• •

precession velocity $\omega_{2}$,

• •

satellite spin $\omega_{3}$,

• •

precession angle $\alpha$, i.e. the angle between the satellite spin axis and the sun-earth direction,

• •

boresight angle $\beta$, i.e. the angle between the focal plane direction and the spin axis.

We simulate only a single detector placed at the centre of the focal plane illuminating a spinning HWP with $f_{s}$ frequency. The systematic affecting the HWP is included in the data at the timeline level and a simple re-binning map-making is 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 employ for those scanning strategies are listed in Tab. 1. As sampling rate we use $60$ Hz.

5.1.1 Input map

The input sky map, used for full-sky simulations, contains: Solar Dipole and Galactic diffuse foregrounds in temperature and a CMB realization, both in temperature and polarization. We decide to include foregrounds only in temperature in order to highlight the temperature to polarization leakage. The input $\mathcal{C}_{\ell}$ used for the CMB realization is compatible with the best fit of Planck 2015 release (Ade et al. 2016b) with no tensor perturbations.

The foreground field is 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.

Such map is modelled 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_{\mathrm{side}}=256$, comfortable enough respect to the Gaussian beam with FWHM=60 arcminutes. In order to evaluate the impact of parameters chosen for the simulation, we have run a case with $N_{\mathrm{side}}=128$, finding the same results in terms of angular power spectra residual, except in the smaller scales, where the pixel size matters independently of the presence of systematic effects.

5.2 Maps and results

5.2.1 Output maps

We perform several simulations with different configurations for the HWP. We vary the spin frequency, precession angle and $I_{s}/I_{\perp}$ ratio. For each simulation we compare input and output maps and compute the B-modes power spectrum. As visual example, we show in Figure 5 (top panel) the output maps for a mission adopting a LiteBIRD-like scanning strategy and solving the Stokes parameters through an HWP with a spin frequency of $1Hz$, a precession angle of $\theta_{0}=1^{\circ}$ and $I_{\perp}/I_{s}={0.514}$.

The residual maps (i.e. difference between output and input maps) in Figure 5 (bottom panel) show the effect of the HWP wobble that induces variations of 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 larger.

Since the effect on the maps is generated by the coupling between the satellite spin and the polarization modulation, affected by the precession, we decided to test several conditions. In particular, slowing the HWP spin speed down to 0.1Hz the effect of the precession is emphasized as you can see in Figure 6, where the systematic effect induced in the T, Q and U maps, reported in histogram equalized color scale, is at the same level of the input map.

5.2.2 Results

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

If a $0.1$Hz spinning HWP is used, the synchronism with the satellite is slightly different and spurious B-mode polarization shows up at different angular scales ($l>100$). What differs in these two cases is the matching between the systematic effect and the satellite spin (Figure 8-(a)(b)).

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. Since 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 sec. 5.1, able to reproduce different satellite observational strategies. The results obtained analyzing those simulations are reported as residual maps (i.e. difference between output and input maps), shown in Figure 9, as root mean square (hereafter RMS) of the residual maps, shown in Table 2, and as B-modes angular power spectra, shown in Figure 10.

In Figure 9 we report the Q residual maps, in histogram equalized color scale, for the case $f_{s}=0.1Hz$, $\theta_{0}=1^{\circ}$ (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_{\perp}/I_{s}$ ratio as representative cases: $0.508$ for Planck, $0.514$ for COrEand LiteBIRD and $0.510$ for WMAP.

The result of this analysis shows that the residual does not depend only on the scanning strategy, but mostly on the combination of scanning strategy, HWP rotation speed and $I_{\perp}/I_{s}$ ratio.

For example, all scanning strategy simulations show the largest effect in the case of an HWP spinning at 0.1 Hz, while 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 Figure 10. Anyway, some strategies produce a spurious peak in the angular power spectrum possibly induced by a resonance between satellite spin and HWP wobbling, i.e. COrE $0.1Hz$ at $\ell\sim 40$, WMAP $0.1Hz$ at $\ell\sim 85$ or LiteBIRD $0.1Hz$ at $\ell\sim 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 orange and blue colored in Figure 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 visible in terms of spurious B-modes (brown colored in Figure 10). These simple cases show that the large scale patterns arisen in the residual map are not related to the whole quality of the map better described by the angular power spectrum.

The power spectra and the maps recovered show the contaminations 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 is the precession angle, the larger is the spurious B-mode signal; the higher is the HWP spin frequency, the greater is 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 induced by the precession of the modulating element in a Stokes polarimeter by moving the contamination at high $\ell$. Fig.10-(b) 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, given that next-generation CMB polarization experiment are designed to reach a sensitivity which is usually quantified as a fraction of the lensing-induced B-modes level.

5.4 Temperature to Polarization leakage

Including in the input for the simulations a map with only temperature foregrounds, we can highlight the temperature-to-polarization leakage effect induced by the systematic for various scenarios. We verified that polarization foregrounds, removed with ideal component separation method, leave one order of magnitude lower residual in terms of P-P leakage.

The HWP wobble induces B-modes which amplitude is proportional to $\theta_{0}$ as you can see in Figure 11 where we report the BB power spectra for different precession angles [$0.5^{\circ}$, $1^{\circ}$, $1.5^{\circ}$, $2^{\circ}$] extracted from the maps scanned by a COrE-like satellite. The induced B-mode signal exceeds the gravitational lensing contribution already for $\theta_{0}=1^{\circ}$.

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

5.5 Neighbouring 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 out the contamination. This can be also obtained by the combination of multiple detectors, observing the same pixel at different times. In order to check this mitigation, we simulate the observation with two different detectors, pointing to different boresight angles $\beta$, shifted by 1∘, for all the proposed scanning strategies. The modulation parameters used are:$f_{s}=1Hz$, $\theta_{0}=1^{\circ}$, $I_{\perp}/I_{s}=0.508$ for Planck, $I_{\perp}/I_{s}=0.514$ for COrEand LiteBIRD and $I_{\perp}/I_{s}=0.510$ for WMAP. The combination of the data from the two detectors results 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 produced with their combination. 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 section 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 $\theta_{0}$ and with a frequency $\gamma$ imposed by the physical properties of the device (the spin frequency $f_{s}$ and the $I_{\perp}/I_{s}$ ratio, including its support). By using Eq. 12 we performed several simulations to assess the level of the systematic effect induced by the precession. We found a strong dependence on $\theta_{0}$ and $\gamma$, both for the fractional residual of the signal and for the power spectrum.

In section 5 we developed the simulation of full-sky observation by a satellite mission, including the systematic effect, and quantified its impact on B-modes 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 intensity signal that dominates at different scales for the strategies selected: at $\ell\sim 40$ for COrE-like satellite, at $\ell\simeq 85$ for a WMAP-like and $\ell\simeq 120$ in the LiteBIRD-like case. The new-era CMB experiment will gain some order of magnitude in sensitivity (Abitbol et al. 2017), few $10^{-4}\mu K$, compared nowadays. By having the analytical formula of the systematic effect induced by the HWP wobbling, it is possible to remove it or, at least, to forecast its impact on data.

In general, the effect of the precession is to induce beats into the timelines. Those beats, in amplitude and frequency, are related to the physical properties of the moving parts. Their effect into the maps depends on specific scanning parameters, and on possible coupling with the beats frequency. We recommend to take this effect into account in the design of an observation strategy, by modelling and measuring the inertia tensor of the moving parts. Once the tensor is measured (or modelled) it can be inserted into Eq. 12 to simulate the impact into the timelines. The scanning strategy must avoid any synchrony with the beats frequency. In this case, redundancy helps cancelling 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.

Acknowledgements.

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

Appendix A Variable definitions

We list in this short appendix the main variables defined in the paper:

$E_{k}$ = Electric field components k=(x,y,z);

$T,Q,U$ = 2D Stokes Parameters;

$\Delta_{i}$ = Component of 3D Stokes vector;

$s$ = 3D generic Stokes vector;

$s_{in}$ = 3D input Stokes vector;

$s_{out}$ = 3D output Stokes vector;

$i_{in}$ = Input Jones vector;

$i_{out}$ = Output Jones vector;

$I$ = Intensity at detector;

$\sigma_{i}$ = Gell-Mann matrix;

$m$ = HWP support mass;

$h$ = HWP support thickness;

$R$ = HWP support radius;

$I_{s}$ = Moment of inertia component along z-axis;

$I_{\perp}$ = Moment of inertia component Along x(y)-axis;

$L_{k}$ = Angular momentum components k=(x,y,z);

$f_{s}=\omega_{s}/2\pi$ = HWP spin frequency respect to z-axis;

$\gamma$ = HWP precession frequency;

$\theta$ = HWP orientation angle respect to z-axis;

$\theta_{0}$ = HWP precession angle;

$\eta$ = HWP rotation angle respect to x-axis composing the precession;

$\xi$ = HWP rotation angle respect to y-axis composing the precession;

$r$ = Distance from the optical axis on the detector plane;

$\omega_{1}$ = Earth revolution velocity;

$\omega_{2}$ = Precession velocity in the satellite scanning strategy;

$\omega_{3}$ = Satellite spin;

$\alpha$ = Satellite precession angle, i.e. the angle between the satellite spin axis and the sun-earth direction;

$\beta$ = Boresight angle, i.e. the angle between the focal plane direction and the spin axis.

Appendix B Precession Theory

We report in this appendix the derivation of the main equations describing a precessing disc, in particular we detail the derivation of Eq. 3.

The contribution to $L_{x}$ due to the rotation about the x-axis is $L_{x}=\frac{d(I_{xx}\eta)}{dt}=I_{xx}\frac{d\eta}{dt}$. We can treat $I_{xx}$ as a constant since 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(\xi)$ on the x-axis. Combining such contributions we get:

[TABLE]

Since $I_{xx}=I_{yy}=I_{\perp}$ and exploring small angle first order approximation Equations 17 read:

[TABLE]

Furthermore, thanks to the same approximation $L_{z}=I_{s}\omega_{s}$. Since we are considering a torque-free system ($dL/dt=0$), both $L_{s}$ and $\omega_{s}$ are constant leading to:

[TABLE]

By introducing $\omega_{x}=\frac{d\eta}{dt}$ and $\omega_{y}=\frac{d\xi}{dt}$ Eq. 19 become:

[TABLE]

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

[TABLE]

The solution for the harmonic motion is (with $A$ and $\phi$ arbitrary constants):

[TABLE]

While for $\omega_{y}$ we get:

[TABLE]

Integrating $\omega_{x}$ and $\omega_{y}$ we obtain:

[TABLE]

In the small angle approximation we impose that $A/\gamma\ll 1$. Such equations reveal that the spin axis rotates around a fixed direction in space. If that direction is along the z-axis then $\theta_{x0}=\theta_{y0}=0$. Assuming the initial conditions $\eta(t=0)=\theta_{0}$ and $\xi(t=0)=0$, and assuming that $A/\gamma=\theta_{0}$ we get:

[TABLE]

The last equations describe the torque-free precession of the spin axis that rotates in space at a fixed angle $\theta_{0}$ respect to the z-axis with a frequency of the precession motion given by $\gamma=\omega_{s}I_{s}/I_{\perp}$ (D. Kleppner & R. Kolenkow 1973).

Considering a thin disc we get $I_{s}=2I_{\perp}$ and so $\gamma=2\omega_{s}$, thus the disc wobbles twice as fast as it spins.

Finally the apparent rate of a thin disc precession to an observer on the rigid body is:

[TABLE]

Appendix C 3D Jones Matrices calculation

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

[TABLE]

The angles $\eta$ and $\xi$ are the precession angles $\eta$ and $\xi$ derived in Section 2 (Eq. 1-2).

The Jones matrix for an half-wave plate with the fast axis at angle $\theta$ with respect to the horizontal axis is:

[TABLE]

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

[TABLE]

C.1 No wobbled case

Combining such matrices we get a general formula for the Jones vector at the on-axis detector of a polarimeter which modulating element precesses:

[TABLE]

These calculations have been realized with the Python package SymPy111http://www.sympy.org/en/index.html which allows symbolic computations with matrices and vectors.

For the ideal case with $\eta=\xi=0^{\circ}$, when no precession occurs, we get the outgoing Jones vector for an ideal polarimeter:

[TABLE]

So the Intensity:

[TABLE]

Where we recognize the modulation terms, at 4 times the frequency of the HWP rotation, within the expressions $\cos^{2}(2\theta)$, $\sin^{2}(2\theta)$ and $\cos(2\theta)\sin(2\theta)$. Passing through the Stokes parameters ($T$,$Q$,$U$) we get back the intensity at the same detector for an ideal polarimeter:

[TABLE]

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

[TABLE]

C.2 No wobbled case

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

[TABLE]

For the general case with $\eta$, $\xi\neq 0^{\circ}$ we find the intensity at the detector:

[TABLE]

where we define the modulating functions $g(\xi,\theta)$ and $f(\eta,\xi,\theta)$:

[TABLE]

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

[TABLE]

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 $\hat{z}$ axis, it also depends on the displacements from the $\hat{x}$ and $\hat{y}$ axes due to the precession.

Appendix D 3D Müller Matrices calculation

Since the Müller formalism is required to propagate partially polarized light (like che CBM one), we need to calculate the Eq. 35 and Eq. 40 by using Müller formalism. 3D Müller matrices $\textsc{M}_{ij}$ are related to Jones matrices (Samim et al. 2016) by:

[TABLE]

where $J$ is the associated Jones matrix, and $\sigma_{n}$ ($n=[0,...,8]$) are the Gell-Mann matrices. The Eq. 41 is valid if and only if we use trace-normalized Gell-Mann matrices (Gell-Mann 1962; Sheppard et al. 2016) defined as follows:

[TABLE]

By using the common polarization matrix:

[TABLE]

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

[TABLE]

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

[TABLE]

D.1 No wobbled case

From Eq. 41 we can find the analogous 3D Müller matrices for the Jones HWP matrix (Eq. 30) and for the Jones polarizer matrix (Eq. 31):

[TABLE]

and by combining the previous matrices we can find the 3D polarimeter Stokes vector ($\eta=\xi=0^{\circ}$):

[TABLE]

Thanks to Eq. 45 we can derive the intensity:

[TABLE]

Note that the Eq. 49 and Eq. 35 give the same result.

D.2 Wobbled case

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

[TABLE]

[TABLE]

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

[TABLE]

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

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

By assuming the field entering the polarimeter has $E_{z}=0$, this is true only for the on-axis detectors and for all the focal plane assuming a telecentric optic system, we can find the general Stokes vector for a wobbling HWP:

[TABLE]

where we used $\Delta_{8}=\frac{1}{\sqrt{2}}\Delta_{0}$, and thanks to Eq. 45 the equivalent intensity:

[TABLE]

Note that the components of the output Stokes vector which are not null (Eq. 63) are $\Delta_{0}^{\prime}$, $\Delta_{4}^{\prime}$, $\Delta_{8}^{\prime}$, so from the definition of the Stokes vector (Eq. 44) it is clearly $E_{z}^{out}=0$. This happens because the polarizing grid doesn’t permit $E_{z}^{out}\neq 0$ for on-axis rays. If we lost this assumption (i.e. for an off-axis detector) it is easy to verify that the output Stokes vector became a function of $\Delta_{0}$, $\Delta_{1}$, $\Delta_{2}$, $\Delta_{4}$, $\Delta_{5}$, $\Delta_{8}$, but the components of the output Stokes vector which are not null are always $\Delta_{0}^{\prime}$, $\Delta_{4}^{\prime}$, $\Delta_{8}^{\prime}$ (Eq. 65).

[TABLE]