© ESO 2022

1. Introduction

In past decades, the cosmic microwave background (CMB) has played a fundamental role in helping to improve our knowledge of the Universe. Measurements of CMB anisotropies in temperature and polarization (E-modes and lensing-induced B-modes) have been decisive in shaping the current cosmological model, from the quantum mechanical origin of the Universe, to its current energy composition (Hinshaw et al. 2013; Planck Collaboration I 2020; Ade et al. 2018, 2014; Aiola et al. 2020; Bianchini et al. 2020). Major advances in the observation of the polarized CMB signal are expected from the forthcoming generation of CMB experiments such as the ground-based Simons Observatory (SO; Ade et al. 2019) and CMB-S4 (The CMB-S4 Collaboration 2020) as well as the LiteBIRD satellite mission (Hazumi et al. 2020). The most ambitious target is the measurement of the primordial B-mode signal. A high-significance detection of the latter would constrain the amplitude of primordial gravitational waves, parameterized in terms of the tensor-to-scalar ratio r. A combination of state-of-the-art cosmological data (Ade et al. 2018) provides the upper bound, r < 0.06 at 95% C.L., updated to r < 0.044 (95% CL) based on a recent re-analysis of Planck data (Tristram et al. 2021). The detection of r would strongly support the validity of the inflation paradigm. On the other hand, a tighter upper bound for r < 0.001 would rule out a large proportion of consistently viable early Universe models (single-field models with typical inflaton excursion that is much larger than the Planck mass scale).

The ambitious sensitivity goals of future surveys, namely: σ(r)≃0.002 from SO, r < 0.001 at 95% CL from CMB-S4, and σ(r)≃0.001 from LiteBIRD, require extraordinary control over systematic effects and noise contamination (LiteBIRD collaboration, in prep.; Ade et al. 2019).

To that end, the use of a polarization modulator like a half-wave plate (HWP) has been included in the design of future surveys, from the Small Aperture Telescopes (SAT) of SO (Ade et al. 2019) to LiteBIRD (Hazumi et al. 2020) and LSPE (Aiola et al. 2012). The HWPs have been already deployed in many polarization-sensitive experiments: MAXIPOL (Johnson et al. 2007), SPIDER (Rahlin et al. 2014), ABS (Kusaka et al. 2014), POLARBEAR (Hill et al. 2016), PILOT (Misawa et al. 2014), BLAST (Galitzki et al. 2016), and EBEX (Reichborn-Kjennerud et al. 2010). These experiments have shown that the use of a HWP can reduce both the 1/f noise (in the case of countinuous spinning; Johnson et al. 2007) and systematic uncertainties related to the pair differencing of orthogonal detectors (Bryan et al. 2016; Essinger-Hileman et al. 2016). However, pernicious systematic effects induced by non-idealities in manufactured HWPs can propagate through the analysis pipeline and bias the final estimation of cosmological parameters, including r. Therefore, a study of the impact of HWP non-idealities on high-level science products is required.

The aim of this work is to provide an exhaustive summary of the mathematical formalism that fully characterizes the behavior of a non-ideal HWP in the context of CMB measurements. We applied this formalism to simulate the effect of HWP non-idealities on the observation of the full sky, using a LiteBIRD-like strategy (Hazumi et al. 2020). A simple analysis is performed for single-frequency observations to showcase the effect of the different HWP systematic parameters at the power-spectrum level. We also conducted a more realistic, multi-frequency analysis where the impact of HWP non-idealities is quantified in terms of a bias in the determination of the tensor-to-scalar ratio, r.

This paper is organized as follows. In Sect. 2, we lay down the mathematical formalism employed in our analysis. Particular care is devoted to clarify a common misunderstanding when dealing with the choice of the matrix formulation (Jones and Mueller) to describe propagation of light through optical systems. In Sect. 3, we describe the scanning strategy and map-making procedure adopted in our simulations. In Sect. 4, we present a simple monochromatic analysis. In Sect. 5, we present the multi-frequency study and the requirements we set on each non-ideal parameter to keep the bias on r under a pre-defined threshold (Δr ≤ 10⁻⁵). Our conclusions are presented in Sect. 6.

2. Matrix representation of HWP optical effects

In this section, we review the two main mathematical formalisms employed to characterize the optical effect of a HWP on incident radiation, namely, on the Jones and Mueller matrix formalisms.

First, we begin with some basic assumptions. Supposing that a quasi-monochromatic wave propagates along a direction orthogonal to the surface of an optical device, we define a coordinate system x − y on the surface of the optical device, so that the incoming wave can be decomposed into an x-component, E_x and a y-component E_y. A wave plate (or retarder) is a phase-shifter, that is, a non-depolarizing linear optical device that modifies the phase of the incident wave. An ideal HWP induces a phase shift of π between the two orthogonal components, E_x, y, of the incident wave. The phase-shift is due to the fact that the components of the incident wave propagate through the HWP with a different index of refraction. The physical properties (e.g., thickness of the plate in case of HWP made of birefringent crystal, design of the stack of mesh filters in case of mesh-HWPs; Pisano et al. 2012) of the HWP can be tuned at the manufacturing stage in such a way that the difference between the optical paths of the two components of the incident wave result in a phase shift of π once the signal emerges from the HWP. The optical axis of the HWP with the highest (or lowest) index of refraction is called the “slow” (or “fast”) axis.

The linear response of the HWP to the incoming signal allows us to represent the output signal emerging from it via linear algebra, that is, via a simple matrix transformation of the input signal. The non-depolarizing property means that the HWP does not decorrelate or randomize the amplitude and phase of the orthogonal components of the incident wave. The non-depolarizing nature of the HWP allows use of the Jones matrix formalism as the matrix representation of the HWP. We will see later that a Mueller matrix approach is also allowed and entirely equivalent to the Jones formalism in this case (i.e., a non-depolarizing device). In the following, we make use of the Jones formalism to provide a much clearer description of the physical effects of HWP non-idealities. The Mueller formalism will be handy for the application of our analysis to future CMB missions. We would like to stress that the choice of the matrix representation of the optical element is independent from the polarization state of the incoming signal. Whether or not an optical element can be represented in terms of a Jones matrix does only depend on the nature of the optical system. In particular, the choice of the matrix representation stems from the non-depolarizing nature of the device. It can be proven (Azzam & Bashara 1977, Sect. 2.10) that for a non-depolarizing device, the degree of polarization of the outgoing signal is always greater than or equal to the degree of polarization of the incoming signal. In contrast, a depolarizing device transfers power out of polarized states into unpolarized states. As such, the degree of polarization of the signal coming out from a depolarizer can be lower than the degree of polarization of the incoming signal. This is the only effect that the Jones formalism is unable to capture. When dealing with such devices, it is better to rely on alternative formalisms. All other non-depolarizing optical systems, including the HWP, can be adequately represented with Jones matrices (Azzam & Bashara 1977).

2.1. Jones matrix formalism

The Jones matrix of an optical system, including that of a HWP, is a 2 × 2 complex matrix applicable to the (E_x, E_y) Jones vector. It is fully characterized by seven real parameters: the real and imaginary part of each matrix element, minus a global phase that is not measurable. The Jones matrix of an ideal HWP with fast axis either along the x-axis or y-axis is expressed as:

$$\begin{array}{c}\hfill J_{\mathrm{HWP},\mathrm{id}}\equiv \left(\begin{array}{cc}1& 0\\ 0& e^{ı\pi }\\ \end{array}\right)=\left(\begin{array}{cc}1& 0\\ 0& -1\\ \end{array}\right).\end{array}$$(1)

Equation (1) has a straightforward interpretation: the field along x is left unchanged by the optical element while the phase of the y-component is shifted by π. However, the behavior of a real HWP can deviate from the ideal case. An expression that also accounts for small deviations of the HWP matrix elements from the ideal case (systematic effects) is as follows (O’Dea et al. 2007):

$$\begin{array}{c}\hfill J_{\mathrm{HWP}}=\left(\begin{array}{cc}1+h_1& {\zeta }_1{e}^{i{\chi }_1}\\ {\zeta }_2{e}^{i{\chi }_2}& -(1+h_2)e^{i\beta }\\ \end{array}\right)\equiv \left(\begin{array}{cc}{A}_1& B_1\\ B_2& A_2\\ \end{array}\right),\end{array}$$(2)

where A₁ is real and A₂, B_1, 2 are complex numbers. The meaning of these non-ideal parameters is as follows:

h₁ and h₂: loss parameters describing the deviation from the unitary transmission of light components E_x, E_y. They are negatively defined parameters, with a range of [−1,0] (light absorption + reflection). In the ideal case, h₁ = h₂ = 0;

β = ϕ − π: where ϕ is the phase shift between the two directions. It accounts for variations of the phase difference between E_x and E_y with respect to the nominal value of π for an ideal HWP. In the ideal case, β = 0;

ζ_1, 2 and χ_1, 2: amplitudes and phases of the off-diagonal terms, coupling E_x and E_y. In practice, if the incoming wave is fully polarized along x (y), a spurious y (x) component would show up in the outgoing wave. Hereafter, we refer to this effect as “cross-polarization”, In the ideal case, ζ_1, 2 = χ_1, 2 = 0.

So far, we have omitted the dependence of the HWP Jones matrix elements on the frequency of the incident wave. The manufacturing of a HWP is always tuned such that a phase shift of π between orthogonal components is realized at a given frequency. Therefore, the matrix elements in both Eqs. (1) and (2) are function of the incident frequency. We assess in Sect. 5 the relevance of this aspect in the context of CMB observations. We have also omitted the dependance on the incident angle, which we neglect in this study (see Sect. 2.3).

In our analysis, we are interested in the possibility that a rotating HWP is employed to modulate the polarization signal. When the HWP is rotated by θ(t)≡ωt, where ωt is the (time-dependent) angle between the HWP fast axis and the x-axis and ω is the angular velocity of the HWP, the Jones matrix is transformed as follows:

$$\begin{array}{c}\hfill J_{\mathrm{RHWP}}(\theta )=R^T(\theta )J_{\mathrm{HWP}}R(\theta )=\left(\begin{array}{cc}{J}_{11}(\theta )& J_{12}(\theta )\\ J_{21}(\theta )& J_{22}(\theta )\\ \end{array}\right),\\ \hfill \mathrm{with}R(\theta )=\left(\begin{array}{cc}\cos \theta & \sin \theta \\ -\sin \theta & \cos \theta \\ \end{array}\right),\end{array}$$(3)

where the time dependence is understood. The explicit expressions of the matrix elements of J_RHWP are:

$$\begin{array}{cc}\hfill J_{11}(\theta )& =(1+h_1){cos}^2\theta -(1+h_2)e^{i\beta }{sin}^2\theta \hfill \\ \hfill & -({\zeta }_1{e}^{i{\chi }_1}+{\zeta }_2{e}^{i{\chi }_2})cos\theta sin\theta \equiv A_1{cos}^2\theta +A_2{sin}^2\theta \hfill \\ \hfill & -(B_1+B_2)cos\theta sin\theta ,\hfill \\ \hfill J_{12}(\theta )& =[(1+h_1)+(1+h_2)e^{i\beta }]cos\theta sin\theta +{\zeta }_1{e}^{i{\chi }_1}{cos}^2\theta \hfill \\ \hfill & -{\zeta }_2{e}^{i{\chi }_2}{sin}^2\theta \equiv [A_1-A_2]cos\theta sin\theta \hfill \\ \hfill & +B_1{cos}^2\theta -B_2{sin}^2\theta ,\hfill \\ \hfill J_{21}(\theta )& =[(1+h_1)+(1+h_2)e^{i\beta }]cos\theta sin\theta \hfill \\ \hfill & +{\zeta }_2{e}^{i{\chi }_2}{cos}^2\theta -{\zeta }_1{e}^{i{\chi }_1}{sin}^2\theta \equiv [A_1-A_2]cos\theta sin\theta \hfill \\ \hfill & +B_2{cos}^2\theta -B_1{sin}^2\theta ,\hfill \\ \hfill J_{22}(\theta )& =(1+h_1){sin}^2\theta -(1+h_2)e^{i\beta }{cos}^2\theta \hfill \\ \hfill & +({\zeta }_1{e}^{i{\chi }_1}+{\zeta }_2{e}^{i{\chi }_2})cos\theta sin\theta \equiv A_1{sin}^2\theta +A_2{cos}^2\theta \hfill \\ \hfill & +(B_1+B_2)cos\theta sin\theta .\hfill \end{array}$$(4)

In the ideal case, Eq. (4) is reduced to

$$\begin{array}{c}\hfill J_{\mathrm{RHWP}}^{\mathrm{ideal}}(\theta )=\left(\begin{array}{cc}cos2\theta & sin2\theta \\ sin2\theta & -cos2\theta \\ \end{array}\right).\end{array}$$(5)

The matrix J_RHWP(θ) is in the reference frame of the telescope. If we refer instead to a fixed reference frame on the sky, we have to take into account also the instrument orientation angle, ψ, such that the expression for the rotated matrix becomes (Bryan et al. 2010):

$$\begin{array}{c}\hfill J_{\mathrm{RHWP}}(\theta )\to J_{\mathrm{RHWP}}(\theta )R(\psi ),\end{array}$$(6)

which is equivalent to Eq. (3) but including the substitution:

$$\begin{array}{c}\hfill \theta (t)=\omega t+\frac{\psi (t)}{2}·\end{array}$$(7)

In this work, we need to take into account that the signal modulated by a rotating HWP is then collected by a polarization-sensitive detector. We consider pairs of polarization-sensitive detectors with orthogonal orientations, as those usually employed in CMB experiments in order to reconstruct the input sky signal more efficiently. The full optical chain traversed by incoming light that is perpendicular to the HWP is described via:

$$\begin{array}{cc}\hfill J_{\mathrm{tot},(x,y)}(\theta )& =J_{\mathrm{pol},(x,y)}{J}_{\mathrm{RHWP}}(\theta ),\mathrm{where}{J}_{\mathrm{pol},x}=\left(\begin{array}{cc}1& 0\\ 0& 0\\ \end{array}\right)\hfill \\ \hfill & \mathrm{or}{J}_{\mathrm{pol},y}=\left(\begin{array}{cc}0& 0\\ 0& 1\\ \end{array}\right).\hfill \end{array}$$(8)

Here, and in the following, the subscript _tot indicates the Jones matrix of the complete optical chain, while the subscripts _out and _in refer to the fields that are, respectively, at the output and input of the optical chain.

The result for both polarizations is:

$$\begin{array}{c}\hfill J_{\mathrm{tot},x}=\left(\begin{array}{cc}{J}_{11}& J_{12}\\ 0& 0\\ \end{array}\right);J_{\mathrm{tot},y}=\left(\begin{array}{cc}0& 0\\ J_{21}& J_{22}\\ \end{array}\right),\end{array}$$(9)

where the θ dependence is understood.

2.2. Coherency matrix

So far, we present the case of a quasi-mono-chromatic fully polarized wave. The CMB signal is only partly polarized and it cannot be easily represented in terms of a Jones vector. The stochastic nature of the quasi-polarized incoming signal requires a statistical description that goes beyond that introduced in the previous section; in other words, the quantity that we can really measure is the time-averaged intensity

$$\begin{array}{c}\hfill P\equiv ⟨\mathbf{E}{\mathbf{E}}^{†}⟩=\left(\begin{array}{cc}T+Q& U-iV\\ U+iV& T-Q\end{array}\right),\end{array}$$(10)

where T, Q, U are the Stokes parameters that describe the polarization state of the wave. They are defined through the time average of the electromagnetic field:

$$\begin{array}{cc}\hfill T& =⟨|E_x{|}^2⟩+⟨|E_y{|}^2⟩,Q=⟨|E_x{|}^2⟩-⟨|E_y{|}^2⟩,U=2\mathrm{Re}[⟨E_x^{\ast }{E}_y⟩],\hfill \\ \hfill V& =2\mathrm{Im}[⟨E_x^{\ast }{E}_y⟩].\hfill \end{array}$$(11)

So, for the observed polarized signal:

$$\begin{array}{cc}\hfill P_{\mathrm{out}}& ={\left(\begin{array}{cc}T+Q& U-iV\\ U+iV& T-Q\end{array}\right)}_{\mathrm{out}}=⟨{\mathbf{E}}_{\mathbf{out}}{\mathbf{E}}_{\mathbf{out}}^{†}⟩\hfill \\ & =⟨J_{\mathrm{tot}}{\mathbf{E}}_{\mathbf{in}}{\mathbf{E}}_{\mathbf{in}}^{†}{J}_{\mathrm{tot}}^{†}⟩=J_{\mathrm{tot}}{\left(\begin{array}{cc}T+Q& U-iV\\ U+iV& T-Q\end{array}\right)}_{\mathrm{in}}{J}_{\mathrm{tot}}^{†}.\hfill \end{array}$$(12)

The signal collected by a total power detector is proportional to the Stokes parameter, T, which can be obtained as half the trace of P_out. Plugging in Eq. (12) each of the expressions for J_tot, (x/y) given in Eq. (9) and taking [P_out]/2, we obtain the expression for the total power collected by the x/y oriented detectors, d_obs, (x/y), as follows:

$$\begin{array}{cc}\hfill d_{\mathrm{obs},x}& =\frac{1}{2}({\left|J_{11}\right|}^2+{\left|J_{12}\right|}^2)T+\frac{1}{2}({\left|J_{11}\right|}^2-{\left|J_{12}\right|}^2)Q\hfill \\ \hfill & +\mathfrak{R}\left(J_{11}{J}_{12}^{\ast }\right)U+\mathfrak{I}\left(J_{11}{J}_{12}^{\ast }\right)V,\hfill \end{array}$$(13a)

$$\begin{array}{cc}\hfill d_{\mathrm{obs},y}& =\frac{1}{2}({\left|J_{22}\right|}^2+{\left|J_{21}\right|}^2)T+\frac{1}{2}({\left|J_{21}\right|}^2-{\left|J_{22}\right|}^2)Q\hfill \\ \hfill & +\mathfrak{R}\left(J_{22}^{\ast }{J}_{21}\right)U-\mathfrak{I}\left(J_{21}^{\ast }{J}_{22}\right)V.\hfill \end{array}$$(13b)

In the case of an ideal HWP, Eqs. (13a)–(13b) become:

$$\begin{array}{cc}\hfill d_{\mathrm{obs},x}^{\mathrm{ideal}}& =\frac{1}{2}[T+cos(4\theta )Q+sin(4\theta )U],\hfill \end{array}$$(14a)

$$\begin{array}{cc}\hfill d_{\mathrm{obs},\mathit{y}}^{\mathrm{ideal}}& =\frac{1}{2}[T-cos(4\theta )Q-sin(4\theta )U].\hfill \end{array}$$(14b)

From Eq. (14), it is clear that the effect of a rotating HWP is to modulate the detected signal from an input linear polarization four times per rotation of the plate.

A useful decomposition of the coherency matrix that is subsequently shown to be useful is given in terms of the Pauli matrices:

$$\begin{array}{c}\hfill P=⟨\mathbf{E}{\mathbf{E}}^{†}⟩=T{\sigma }_T+Q{\sigma }_Q+U{\sigma }_U+V{\sigma }_V,\end{array}$$(15)

where:

$$\begin{array}{c}\hfill {\sigma }_T=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right),{\sigma }_Q=\left(\begin{array}{cc}1& 0\\ 0& -1\end{array}\right),{\sigma }_U=\left(\begin{array}{cc}0& 1\\ 1& 0\end{array}\right),{\sigma }_V=\left(\begin{array}{cc}0& -i\\ i& 0\end{array}\right).\end{array}$$(16)

It is also useful to express the Stokes vector s = (T, Q, U, V) in terms of the elements of the coherency matrix:

$$\begin{array}{c}\hfill \mathit{s}=A\mathbf{P},A=\left(\begin{array}{cccc}1& 0& 0& 1\\ 1& 0& 0& -1\\ 0& 1& 1& 0\\ 0& ı& -ı& 0\\ \end{array}\right),\mathbf{P}=\left(\begin{array}{c}{P}_{11}\\ P_{12}\\ P_{21}\\ P_{22}\\ \end{array}\right),\end{array}$$(17)

where P = ⟨E × E^†⟩ is the Kronecker product of the incoming signal with itself. Although the incoming signal is no longer represented as a Jones vector as it was for the fully polarized wave, we note that the Jones formalism still allows a full description of the effects of the train of optical elements on the incoming signal.

2.3. Mueller matrix formalism

To express directly how the Stokes parameters, s = (T, Q, U, V), get transformed by the observation, the Mueller formalism can be adopted. Analogously to the Jones formalism, the observed parameters become s_obs = Ms, where M is the Mueller matrix of the whole optical element.

In moving on from the Jones to the Mueller matrix, we can easily see that

$$\begin{array}{cc}& {\mathbf{P}}_{\mathbf{out}}=(J_{\mathrm{tot}}\times J_{\mathrm{tot}}^{†}){\mathbf{P}}_{\mathbf{in}}\to (A^{-1}{s}_{\mathrm{out}})=(J_{\mathrm{tot}}\times J_{\mathrm{tot}}^{†})(A^{-1}{s}_{\mathrm{in}})\hfill \\ \hfill & s_{\mathrm{out}}=A(J_{\mathrm{tot}}\times J_{\mathrm{tot}}^{†})A^{-1}{s}_{\mathrm{in}}=M{s}_{\mathrm{in}},\mathrm{where}M\equiv A(J_{\mathrm{tot}}\times J_{\mathrm{tot}}^{†})A^{-1}.\hfill \end{array}$$(18)

In a similar fashion, using the decomposition of the coherency matrix in terms of the Pauli matrices, we can show that:

$$\begin{array}{c}\hfill M_{\mathit{ij}}=\frac{1}{2}\mathrm{Tr}({\sigma }_i{J}_{\mathrm{tot}}{\sigma }_j{J}_{\mathrm{tot}}^{†}),\end{array}$$(19)

where i, j = {T, Q, U, V}.

It is easy to show that the Mueller matrix of an ideal HWP, for a spinning angle θ = 0 is simply:

$$\begin{array}{c}\hfill M_{\mathrm{HWP}}^{\mathrm{ideal}}=\left(\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& -1& 0\\ 0& 0& 0& -1\\ \end{array}\right).\end{array}$$(20)

In the case of a non-ideal HWP, the elements along the diagonal will deviate from unity, and the off-diagonal elements could be also populated. The most general expression of the Mueller matrix of a realistic HWP is expressed as:

$$\begin{array}{c}\hfill M_{\mathrm{HWP}}=\left(\begin{array}{cccc}{T}_1& {\rho }_1& a_1& b_1\\ {\rho }_2& T_2& a_2& b_2\\ a_3& a_4& c_1& -s_1\\ b_3& b_4& s_2& c_2\\ \end{array}\right).\end{array}$$(21)

By transforming the most general Jones matrix in Eq. (2) according to Eq. (18) or (19), we can find the relation between the Mueller matrix elements and the parameters of the HWP non-idealities introduced in the Jones formalism. The complete expression of the Mueller matrix elements can be found in Appendix A. Here, we would like to make note of the following: In the same framework as that of Eq. (20), if no cross-polarization is present (i.e., ζ₁ = ζ₂ = 0), then the Mueller matrix in Eq. (21) becomes block-diagonal, with a_i = b_i = 0 for i = 1, 2, 3, 4. In addition, we obtain T₁ = T₂ ≡ T, ρ₁ = ρ₂ ≡ ρ, c₁ = c₂ ≡ c, and s₁ = s₂ ≡ s. This is the expression that can be commonly found in the literature (compare e.g., Bryan et al. 2010).

Similarly to the Jones formalism, the Mueller matrix for the complete optical system is simply given by the multiplication of the individual matrices¹.

If we set the matrix elements of such a Mueller matrix to:

$$\begin{array}{c}\hfill M_{x/y}=\left(\begin{array}{cccc}{M}^{\mathit{TT}}& M^{\mathit{TQ}}& M^{\mathit{TU}}& M^{\mathit{TV}}\\ M^{\mathit{QT}}& M^{\mathit{QQ}}& M^{\mathit{QU}}& M^{\mathit{QV}}\\ M^{\mathit{UT}}& M^{\mathit{UQ}}& M^{\mathit{UU}}& M^{\mathit{UV}}\\ M^{\mathit{VT}}& M^{\mathit{VQ}}& M^{\mathit{VU}}& M^{\mathit{VV}}\\ \end{array}\right),\end{array}$$

where the subscript x/y implies that the optical train ends with a polarizer along x/y, and the same subscript is understood in each of the matrix elements.

The total power collected by a single detector – being proportional to the Stokes parameter T – corresponds to taking the first row of the Mueller matrix and multiplying it by the Stokes vector of the input signal. The general expression of the signal obtained by one detector is:

$$\begin{array}{c}\hfill d_{\mathrm{obs},(x/y)}=M_{x/y}^{\mathit{TT}}T+M_{x/y}^{\mathit{TQ}}Q+M_{x/y}^{\mathit{TU}}U+M_{x/y}^{\mathit{TV}}V.\end{array}$$(22)

From Eq. (18), we can see the equivalence between each of the coefficients $M_{x/\mathit{y}}^{\mathit{TT}},M_{x/\mathit{y}}^{\mathit{TQ}},M_{x/\mathit{y}}^{\mathit{TU}},M_{x/\mathit{y}}^{\mathit{TV}}$ and the terms appearing in Eqs. (13a)–(13b), derived from the Jones formalism.

2.4. Study in context

In this work, we only consider the specific case-study of a detector at the boresight collecting signal coming from light hitting the HWP perpendicularly. In the absence of beam convolution, the light rays convolved by the optical system on the detector at boresight are the orthogonal ones (Lamagna et al. 2020). That is why our Mueller matrix elements depend only on the spinning angle of the HWP and does not depend on the incidence angle, which would be the case in general. To take into account non-orthogonal incidence, we should include the dependence on the angle of incidence in the computation of the Mueller matrix elements (Salatino et al. 2018; Essinger-Hileman et al. 2016) and convolve the total matrix by the beam (Duivenvoorden et al. 2021). In this work, we have chosen a simplified approach neglecting the coupling between beam convolution and HWP non-idealities. This coupling can be a source of additional systematic effects, such as temperature-to-polarization leakage, which are not included in this study and instead deferred to a future work (Patanchon et al., in prep.). We refer to Patanchon et al. (in prep.); Salatino et al. (2018), D’Alessandro et al. (2019), Duivenvoorden et al. (2021), Essinger-Hileman et al. (2016) for studies that include the effect of slant incidence.

By considering orthogonal incidence only, we are implicitly assuming a symmetric beam. However, we note that even if the beam is asymmetric in its shape, it is symmetrized to some extent due to the scanning strategy. This symmetrization effect further motivates our choice of restricting the study to orthogonal incidence.

Based on Eq. (22) combined with Eq. (21), we may notice that in the ideal case, we have M^TV = 0. An optical system employing a realistic (non-ideal) HWP allows detection of V-mode signal (Nagy et al. 2017). However, in this work, we will restrict to the case V = 0 as expected in the standard cosmological model. We note that several mechanisms have been proposed to generate a certain amount of CMB circular polarization (Lembo et al. 2021; Zarei et al. 2010; Alexander et al. 2009, 2020; Sadegh et al. 2018; Inomata & Kamionkowski 2019; Vahedi et al. 2019; Bartolo et al. 2019). Nevertheless, the predicted signal is very faint, and therefore justifies our choice of assuming V = 0 in the next section.

A final note before we move on to discuss other aspects of our analysis. In principle, for the simple case we are studying here (i.e., a non-depolarizing optical system or a normal incidence) we could have worked with the Jones formalism, provided that the input and output signals were described in terms of the coherency matrix. Nevertheless, we decided to switch to the Mueller formalism since it provides a more direct handle to the Stokes parameters.

3. Application to future CMB missions: Scanning strategy and mapmaking

In the previous section, we laid down the mathematical formalism to describe the effects of a continuously rotating, non-ideal HWP. In this section, we present the experimental setup we want to investigate, in which such a HWP is employed. We are interested in quantifying the impact of HWP non-idealities in the context of future CMB observations. In particular, we focus on simulating the performance of a HWP on board of a LiteBIRD-like satellite experiment. LiteBIRD (Hazumi et al. 2020) is a satellite mission expected to be launched in the late 2020s, whose main scientific target is the detection of an inflationary signal with a precision on r of σ_r ≲ 10⁻³. LiteBIRD will perform a full-sky survey over three years at the Sun-Earth Lagrangian point L2, using three telescopes (Low-Frequency Telescope LFT, Medium-Frequency Telescope MFT, and High-Frequency Telescope HFT) that observe in 15 frequency bands between 34 and 448 GHz. Each telescope will use a HWP as a polarization modulator. The parameters defining the LiteBIRD scanning strategy are listed in Table 1 and a map of the expected satellite footprint is shown in Fig. 1.

Fig. 1. Map displaying the number of samples collected in each pixel, for a pair detectors at boresight. Healpix resolution: Nside = 512.

We developed a software package that simulates a realistic scanning strategy of a satellite mission, and subsequently reconstructs maps of the T, Q, and U Stokes parameters from the simulated observations. As noted earlier in this paper, we considered a single pair of polarization-sensitive detectors located at boresight (i.e., perfectly centered on the instrument focal plane). Both detectors share a view of the same sky patch, but they are oriented at 90° with respect to each other, so as to remain sensitive to orthogonal polarization directions. The full optical chain as viewed by the sky signal entering the telescope is then composed by a continuously rotating HWP followed by a pair of orthogonal polarization sensitive detectors. We also account for the relative orientation of the telescope with respect to the local coordinate system that locally identifies Q and U. The Mueller matrix of the full optical chain is obtained as the product of the individual Mueller matrices:

$$\begin{array}{c}\hfill M_{\mathrm{full},\mathrm{j}}(t)=M_{\mathrm{pol},j}{M}_{\mathrm{rot}}^T(t)M_{\mathrm{HWP}}{M}_{\mathrm{rot}}(t)M_{\psi }(t),j=x,y,\end{array}$$(23)

where the time dependence has been made explicit where relevant. As explained in the previous section, M_pol, j is the Mueller matrix for the polarizers along the j = x, y direction, M_rot is the Mueller version of R(θ), M_HWP is the Mueller matrix of a non-ideal HWP (defined as in Eq. (21)). Finally, M_ψ is the Mueller version of the rotation matrix R(ψ), that takes into account the angle ψ between the instrument and the local coordinate system on the sky. The total power collected by a single detector at a given time, t, also known as time-ordered data (TOD), is then

$$\begin{array}{c}\hfill d_{\mathrm{obs}}(t)={\displaystyle \sum _j}{M}_{\mathrm{full},j}^{\mathit{TT}}T+M_{\mathrm{full},j}^{\mathit{TQ}}Q+M_{\mathrm{full},j}^{\mathit{TU}}U.\end{array}$$(24)

Since the instrument has finite angular resolution, the sky is discretized in small patches (pixels). At a given time, t_i, one pixel, p, is observed. We follow the HEALPix pixelization scheme (Gorski et al. 2005). The sky is divided into $12\times N_{\text{side}}^2$ pixels, where N_side = 512 is chosen in such a way that the size of each pixel is smaller than the angular resolution of the experiment (0.5 deg at 100 GHz; Hazumi et al. 2020). Ignoring beam convolution², the total power collected by a single detector at a given time t_i is then

$$\begin{array}{c}\hfill d_{\mathrm{obs}}(t_i)=M_{\mathrm{full},p}^{\mathit{TT}}(t_i)T(p)+M_{\mathrm{full},p}^{\mathit{TQ}}(t_i)Q(p)+M_{\mathrm{full},p}^{\mathit{TU}}(t_i)U(p)+n_i,\end{array}$$(25)

where the sum of j = x, y is now understood. We also allow for the possibility of instrumental noise, n_i, to be added to the i-th time sample. We note that the Mueller matrix of the optical system also depends on the observed pixel through the relative orientation with respect to local coordinate system. The TOD equation can be arranged in a matrix notation:

$$\begin{array}{cc}& {\mathbf{d}}_{\mathrm{obs}}(t)=\left(\begin{array}{c}⋮\\ d_{\mathrm{obs}}(t_i)\\ ⋮\end{array}\right)\hfill \\ & =\left(\begin{array}{ccccccc}\ddots & ⋮& ⋮& ⋮& ⋮& ⋮& \\ \dots & 0& M_{\mathrm{full},p}^{\mathit{TT}}(t_i)& M_{\mathrm{full},p}^{\mathit{TQ}}(t_i)& M_{\mathrm{full},p}^{\mathit{TU}}(t_i)& 0& \dots \\ & ⋮& ⋮& ⋮& ⋮& ⋮& \ddots \end{array}\right)\left(\begin{array}{c}⋮\\ T(p)\\ Q(p)\\ U(p)\\ ⋮\end{array}\right)+\left(\begin{array}{c}⋮\\ n_i\\ ⋮\end{array}\right)\hfill \\ & \equiv \mathbf{A}(t){\mathbf{m}}_{\mathrm{in}}+\mathbf{n}(t),\hfill \end{array}$$(26)

where A is the “pointing matrix” with the following dimension:

N_samples × (3 × N_pixels), m_in is the vector of Stokes parameters with the dimension: (3 × N_pixels), and n is the vector of instrumental noise contributions with dimension (N_samples). The rows of the pointing matrix have non-zero elements only for the samples i in which the pixel p is observed. Finally, d_obs(t) is the full TOD vector, with the dimension (N_samples). The number of samples is easily computed from the mission duration and the data sampling rate (see Table 1). Equation (26) can be inverted to reconstruct the Stokes maps from the TOD with the mapmaking procedure, as, for example, in Tegmark (1997), Natoli et al. (2001), Keihanen et al. (2005). In the case of non-correlated noise (i.e., ⟨nn^T⟩=σ²𝕀, where σ is the uniform noise standard deviation), it can be shown that the reconstructed sky signal in matrix form is

$$\begin{array}{c}\hfill {\mathbf{m}}_{\mathrm{out}}={({\mathbf{B}}^T\mathbf{B})}^{-1}{\mathbf{B}}^T{\mathbf{d}}_{\mathrm{obs}}={({\mathbf{B}}^T\mathbf{B})}^{-1}{\mathbf{B}}^T\mathbf{A}{\mathbf{m}}_{\mathrm{in}}+{({\mathbf{B}}^T\mathbf{B})}^{-1}{\mathbf{B}}^T\mathbf{n},\end{array}$$(27)

where B is the estimated pointing matrix, usually constructed from the actual attitude of the telescope and from a pre-launch measurement of the instrument optical elements.

In the next sections, we will use the following formalism: A is the “real” pointing matrix and includes all the HWP systematic effects that might affect our data (we will also refer to it as “TOD HWP”); B is the estimated pointing matrix, which we refer to as the “solver” matrix or as “map-making HWP”. Here, B is to be used in the mapmaking process and we go on to construct it either as identical to A (to correctly recover the input sky signal at the right hand side of Eq. (27)) or to be different from A (to propagate the effect of unaccounted systematics).

Armed with this basic formalism, we are now ready to follow the steps implemented in the code: (1) at each time step t_i, the observed pixel p is identified, the signal d_obs(t_i) as given in Eq. (25) is computed (using the matrix A) for both detectors observing that pixel; (2) from Eq. (27), the two quantities (B^TB)_p and (B^Td_obs)_p are computed for the pixel p (it should be noted that the algorithm does not require storage of the TOD vector); (3) every time a pixel p is observed, the two quantities above are summed to those already computed in previous steps for the same pixel. The number of samples falling in each pixel is also stored to produce a coverage map (Fig. 1) at the end; (4) we cycle over the first three points for all the time samples collected by the instruments; (5) at the end of the mission time, the Stokes maps are estimated using Eq. (27).

In the next few sections, we demonstrate how we applied the algorithm above to two case studies to quantify the effects of HWP systematics in the context of future satellite missions. In Sect. 4, we consider the simple case of single-frequency observations of a CMB-only sky with white noise. This allows us to easily understand the impact of each class of HWP non-idealities on the output signal. In Sect. 5, we consider a more realistic scenario that consists in multi-frequency observations of a more complex sky with the CMB signal contaminated by the presence of (frequency-dependent) foreground emission. Since the systematic effects treated in this paper do not cause temperature-to-polarization leakage, we only focus on polarized emission. The analysis setup and results for both the mono-chromatic and multi-frequency studies are described in detail in the following sections.

4. Single-frequency case

In this section, we describe how we applied the algorithm described in Sect. 3 to a simple case study. We simulated observations with a future CMB satellite. We assumed perfectly monochromatic detectors, that is, we observed the sky at one frequency. The input sky is given by the CMB signal only, both in temperature and polarization. We also assumed a simple model for the instrumental noise, properly rescaled to take into account that our simulated observations only employ two detectors. Finally, we assumed that no systematics but those related to HWP non-idealities are present. To quantify the impact of HWP non-idealities, we compared the BB power spectrum residuals to the ideal BB power spectrum that would be observed in the absence of HWP systematics.

4.1. Input sky and experimental setup

The input sky is composed of the CMB signal only. We computed 100 T,  Q,  U map realizations from the same fiducial set of TT,  TE,  EE,  BB CMB power spectra. The latter were computed with the default values of the Boltzmann solver CAMB (Lewis et al. 2000). We assumed a flat ΛCDM cosmology with three families of active neutrinos with total mass 0.06 eV. We note, however, that the choice of a different cosmology would have had negligible impact on the results presented in this section. We allowed for a non-zero value of the tensor-to-scalar ratio, which we set to r = 0.003. We use the lensed version of the spectra as generated by CAMB, to account for the extra variance from lensing particularly relevant for the BB signal. Maps were generated with the HEALPix routine synfast as implemented in the python package healpy. The generation of a large number of CMB realizations is needed to account for cosmic variance: in the absence of systematic effects, the ensemble average of the observations is expected to reproduce the fiducial input spectra (Gerbino et al. 2020).

We adopted the publicly available instrumental specifications of the future LiteBIRD satellite (summarized in Table 1). We set N_side = 512 and smooth the signal with a FWHM = 30.8 arcmin gaussian beam to simulate a LiteBIRD-like angular resolution at 150 GHz (Table 2; Hazumi et al. 2020). When generating maps, we also take into account the effect of the pixel window function. Each of the 100 sky realizations is used as the input sky for the scanning and mapmaking algorithm described in the previous section. We include experimental noise to highlight differences in the noise bias due to different choices of the map-making matrix B. Since the main focus of this work is to study the effects of HWP non-idealities, we only consider a simple noise contribution that is isotropic and uncorrelated. These properties translate to a white noise spectrum in harmonic space: N_ℓ = Σ_mσ²/(2ℓ+1) = σ². In practice, n_i in Eq. (25) is drawn from a zero-mean Gaussian distribution with variance ${\stackrel{\sim }{\sigma }}^2$. The noise variance is determined from the intrinsic detector sensitivity NET and the number of samples t_p in which each pixel is observed. Assuming the specifications in Table 2, we obtain the following: $\stackrel{\sim }{\sigma }=3.16\mathrm{\mu }\mathrm{K}$.

Each sky realization was observed twice. First, we ran the code setting the pointing matrix A to be equal to the solver matrix B, both with an ideal HWP. In a second run, we instead imposed B to be different from A; thus, A takes into account a non-ideal HWP while B can either consider an ideal or a non-ideal HWP. The exact values of the HWP parameters entering A and B are given in Sect. 4.2. In both cases, the output of the code consists of a reconstructed map that is the sum of the observed CMB signal and instrumental noise. In total, we have two sets of 100 output maps: one set of 100 ideal output maps and another set of 100 realistic output maps. We computed the output TT,  TE,  EE,  BB spectra from each output map employing the HEALPix routine anafast from the healpy python package. When generating the output spectra, we assume full-sky observations and we correct for the beam smoothing effect and for the pixel window function. As stated at the beginning of this section, we focus on BB residuals to quantify the impact of HWP non-idealities. Therefore, going forward, we only focus on BB spectra $C_{\mathcal{l}}^{\text{BB}}$. In total, we have two sets of 100 output spectra: a set of i = 1, 2, …, 100 ideal spectra $C_{\ell ,i}^{\mathrm{BB}}$(ideal), and another set of i = 1, 2, …, 100 realistic spectra $C_{\ell ,i}^{\mathrm{BB}}$(realistic) $=C_{\ell ,i}^{\mathrm{BB}}$(w/ systematics)$+N_{\ell ,i}^{\mathit{BB}}$³. A schematic picture of the procedure described above is depicted in Fig. 2.

Fig. 2. Scheme of the procedure for the monochromatic analysis. From a set of n = 100 input maps, we obtain two sets of ideal and realistic output spectra, depending on whether we allow for the TOD matrix A to be equal to or different from the mapmaking matrix B, respectively.

To get the BB residuals due to systematics, we first need to noise-debias the observed spectra. The noise bias is obtained as the average over noise spectra computed from 100 noise maps drawn from the noise covariance matrix 𝒩 ≡ σ²(B^TB)⁻¹. Clearly, different choices of the B matrix would lead to different noise on the maps.

We considered two different scenarios. In both scenarios, the map-making matrix B is kept fixed while we consider different choices for the pointing matrix A, so to illustrate the effects of unaccounted systematics.

First, we took B to be the map-making matrix of an ideal optical system, that is, the solver as computed from Eq. (23) when using Eq. (20) for the ideal HWP. We reconstructed the output map for different choices of the pointing matrix A = A(h₁, h₂, ζ₁, ζ₂, β, χ₁, χ₂) to highlight the impact of each class of HWP non-idealities represented by the parameters h₁, h₂, ζ₁, ζ₂, β, χ₁, χ₂. In detail, we considered the following three classes: non-ideal transmittance with A ≡ A(h₁, h₂ = const, ζ₁ = ζ₂ = 0, β = 0, χ₁ = χ₂ = 0), non-vanishing cross-polarization with A ≡ A(h₁ = h₂ = 0, ζ₁, ζ₂ = const, β = 0, χ₁, χ₂ = const), and non-ideal phase-shift with A ≡ A(h₁ = h₂ = 0, ζ₁ = ζ₂ = 0, β = const, χ₁ = χ₂ = 0). We extended the second class by perturbing also the phases χ_1, 2 of the cross-polarization terms, even though in general they can be reabsorbed in a redefinition of ζ. Within each class, we explored different values of the non-vanishing non-ideal parameters. A summary of these values is reported in Table 3.

A second setup was then considered. We took B to be the solver matrix of a more realistic optical system, that is, we allowed for the non-ideal parameters to be non-vanishing one at a time. We again computed the output maps for different choices of the pointing matrix A, similarly to what was done in the previous case. Values employed to build the pointing matrix are reported in Table 4, where we highlight in boldface the values used to build the map-making matrix B (we use the subscript s for them to indicate the “solver”). The first setup allows for a characterization of the residuals when the HWP systematics are not accounted for in the map-making; the second one is when they are accounted for, but with a mismatch between their estimate in the solver and their actual value in the pointing matrix. In this second setup, we did not consider perturbations of the phases χ_1, 2, as the residuals are mainly driven by the value of ζ.

We adopted exaggerated values for the non-ideal parameters in order to make their effect on the power spectra well visible. In Sect. 5, we consider more realistic values in order to propagate their effects to r.

4.2. Results of the single-frequency analysis

Here, we present and discuss the results obtained for the single-frequency analysis. The BB power spectra for the two cases of ideal mapmaking matrix B and non-ideal B discussed in the previous section are presented in Figs. 3 and 4.

Fig. 3. $C_{\ell }^{\mathit{BB}}$ power spectra with CMB only (no foregrounds) and noise simulations. This figure summarizes the results in the case of ideal map-making matrix B (ideal HWP). We report in blue the CMB (with r = 0.003) + noise spectra in the ideal case (TOD matrix equal to the mapmaking matrix, A = B). The standard deviation of the 100 ideal CMB realizations is shown as a shaded blue region. Spectra from the output maps obtained with the choice A ≠ B→ (systematic+noise) are shown in dashed-dotted and dashed lines. In the case of h (or ζ), the dashed lines have the same h1 + h2 (or ζ1 + ζ2), whereas this sum is different in the case shown in the dashed-dotted lines: we see that the dashed lines are overlapping, as the 2θ terms in Eq. (28) and (30) are canceled out in the mapmaking process. The noise bias is shown in dotted. In the case of ideal B, the noise bias is the same regardless from the parameter perturbed. We refer to the main text for a more detailed discussion.

Fig. 4. Percent difference of the noise de-biased $C_{\ell }^{\mathit{BB}}$(CMB + noise + systematics) with respect to the ideal $C_{\ell }^{\mathit{BB}}$(CMB only), when assuming a non-ideal mapmaking matrix B. The shaded gray region shows the standard deviation of the ideal CMB realizations, normalized to their mean. The parameters with and without the subscript s enter the mapmaking matrix B/TOD matrix A. In each panel, we fix the value of each class of solver parameters (indicated in the corresponding title), except for the solid blue curve, showing (for reference) the case with non-ideal pointing matrix A and ideal mapmaking matrix B. In orange dashed lines, we report the results in the case A = B. In the dashed-dotted lines, we report the results in the case A ≠ B (different colors correspond to different values of the HWP parameters, see legend). In the right-most panel, the case with ideal B (blue solid) is divided by a factor of 10 with respect to the actual signal to ease the comparison with the other curves. The lines are wiggly because of the noise de-biasing. See the text for discussion.

4.2.1. Ideal map-making matrix B

The main findings in the case of an ideal solver matrix B are summarised in Fig. 3, where we show: in blue solid line, the average spectrum obtained from 100 realizations of ideal CMB maps, to which we sum the noise bias $⟨C_{\ell }^{\mathrm{BB}}$(ideal)$⟩+⟨N_{\ell }^{\mathit{BB}}⟩$ (this is our reference spectrum); in colored dashed or dashed-dotted line, we show the average spectrum over 100 CMB realizations affected by noise and one kind of systematics at a time $⟨C_{\ell }^{\mathrm{BB}}$(realistic)⟩; in purple dotted line, we show the noise spectrum $⟨N_{\ell }^{\mathit{BB}}⟩$.

Since the map-making matrix B is fixed to the ideal case, the systematic effects due to non-ideal parameters are not taken into account in the map-making stage. Because of this, the noise spectrum shown in Fig. 3 is the same for all the panels (see Sect. 4.1).

Next, we discuss the impact of each class of non-idealities. The parameters h₁, h₂ have the effect of shifting the spectra to lower amplitudes. This can clearly be seen by expanding $M_x^{TX}$ (Eq. (A.3)) for small h, with all the parameters but h set to zero:

$$\begin{array}{cc}\hfill M_x^{\mathit{TT}}& \simeq \frac{1}{2}(1+h_1+h_2)+\frac{1}{2}(h_1-h_2)cos(2\theta ),\hfill \\ \hfill M_x^{\mathit{TQ}}& \simeq \frac{1}{2}(h_1-h_2)cos(2\theta )+\frac{1}{2}(1+h_1+h_2)cos(4\theta ),\hfill \\ \hfill M_x^{\mathit{TU}}& \simeq \frac{1}{2}(h_1-h_2)sin(2\theta )+\frac{1}{2}(1+h_1+h_2)sin(4\theta ).\hfill \end{array}$$(28)

The top left panel of Fig. 3 shows the effects of h_1, 2 on the BB spectrum.

In the dashed lines, we report the results for two combinations of h_1, 2 which share the same value of h₁ + h₂ = −0.2 but different h₁ − h₂. The two dashed curves overlap almost perfectly. This can be explained with the fact that the 2θ terms in Eq. (28), scaled by h₁ − h₂, are canceled out in the map-making procedure by considering orthogonally polarized detectors⁴.

The case of β ≠ 0 is shown in the top right panel of Fig. 3. Expanding Eq. (A.3) with respect to β and setting the other parameters to zero, we obtain:

$$\begin{array}{cc}\hfill M_x^{\mathit{TT}}& =\frac{1}{2},\hfill \\ \hfill M_x^{\mathit{TQ}}& =\frac{1}{4}(1-cos\beta )+\frac{1}{4}(1+cos\beta )cos(4\theta ),\hfill \\ \hfill M_x^{\mathit{TU}}& =\frac{1}{4}(1+cos\beta )sin(4\theta ).\hfill \end{array}$$(29)

The 4θ terms are now scaled by (1 + cos β), which acts to reduce the output signal when we fix β_s = 0. Indeed, in the ideal case of β = 0, we should measure a signal with an ideal phase-shift of exactly 180°. However, the actual signal is detected with a slightly different phase-shift and a fraction of the input power is not transfered. This can be observed in the plot, where the shift toward lower amplitude of the BB spectrum is more enhanced for higher values of β (always smaller than 90°).

The cases with ζ ≠ 0 are shown in the lower panels of Fig. 3. On the left, we set the phases χ = 0 and show: in dashed, the cases with fixed ζ₁ + ζ₂ = 0.2; in dashed-dotted, a case with a different ζ₁ + ζ₂. In the right panel, we set χ ≠ 0 and ζ₁ = ζ₂ = 0.1. The expanded expressions of Eq. (A.3) with only the ζ and χ different from zero are the following:

$$\begin{array}{cc}\hfill M_x^{\mathit{TT}}& \simeq \frac{1}{2}+\frac{1}{2}({\zeta }_{1}cos{\chi }_1-{\zeta }_{2}cos{\chi }_2)sin(2\theta ),\hfill \\ \hfill M_x^{\mathit{TQ}}& \simeq \frac{1}{2}cos(4\theta )-\frac{1}{2}({\zeta }_{1}cos{\chi }_1-{\zeta }_{2}cos{\chi }_2)sin(2\theta )\hfill \\ \hfill & -\frac{1}{2}({\zeta }_{1}cos{\chi }_1+{\zeta }_{2}cos{\chi }_2)sin(4\theta ),\hfill \\ \hfill M_x^{\mathit{TU}}& \simeq \frac{1}{2}sin(4\theta )+\frac{1}{2}({\zeta }_{1}cos{\chi }_1-{\zeta }_{2}cos{\chi }_2)cos(2\theta )\hfill \\ \hfill & +\frac{1}{2}({\zeta }_{1}cos{\chi }_1+{\zeta }_{2}cos{\chi }_2)cos(4\theta ).\hfill \end{array}$$(30)

In both bottom panels of Fig. 3, we can appreciate the effect of cross-polarization: the shape of the BB spectrum is modified by the E → B leakage, enhancing the final spectrum. This should be compared with the effect of h and β, which instead act to rescale the input BB spectrum. For this reason, a value of ζ of the same order of magnitude of h causes a more prominent effect on the final spectra. Of course, the effect is stronger in combination with ζ_1, 2.

In the left panel, we can see that the two dashed lines, corresponding to two different choices of ζ₁ − ζ₂, overlap: the 2θ terms, canceled out by the map-making procedure, don’t impact on the BB spectrum. In the right panel, the three curves are similar but not perfectly overlapping, as we are varying the sum (ζ₁ cos χ₁+ζ₂ cos χ₂) by changing the phases of χ. We note that the difference between the curves is mainly driven by the high value of ζ_1, 2. At first order, the χ phases act as a small real multiplicative factor (see Eq. (30)), so perturbing them has the same effect of slightly changing the module of ζ. For that reason, in the following, we act only on ζ and keep χ = χ_s = 0.

4.2.2. Non-ideal mapmaking matrix B

In Fig. 4, we report the results for the study with a non-ideal mapmaking matrix. We plot the percent difference of the average over 100 BB (noise de-biased) spectra affected by systematics with respect to the average over the spectra from the same CMB realizations, not affected by systematics.

In this case, we use the following values for the parameters in the solver matrix B (indicated with the subscript s): h_1, s = h_2, s = −0.1,  β_s = ζ_s = 0 (left panel), β_s = 15° , h_s = ζ_s = 0 (middle panel), ζ_1, s = ζ_2, s = 0.1,  β_s = h_s = 0 (right panel). Results corresponding to this choice are shown in Fig. 4, using dashed and dashed-dotted lines. For reference, we also include the residual spectra obtained when assuming the ideal B (blue solid). We want to stress that since the map-making matrix B is different between the panels, also the noise bias is slightly different (see Sect. 4.1).

Different cases are shown in Fig. 4:

• (i)

  A = B, that is, TOD HWP equal to map-making HWP (orange dashed line in all panels), which would perfectly correct for systematic effects in a noiseless case. The correction is less visible in our case due to noise;

• (ii)

  for x ≡ h, ζ, the green (red) dashed-dotted line in the leftmost panel corresponds to x₁ + x₂ = x_1, s + x_2, s (x₁ + x₂ ≠ x_1, s + x_2, s). We note that the green line overlaps with the orange line, as expected in the case A = B;

• (iii)

  for β, the green (red) dashed-dotted line in the middle panel has β < β_s (β > β_s), which gives a slightly positive (negative) shift. In fact, we would expect a correction of order cos β_s, while a higher (smaller) cos β enters in the TOD matrix;

• (iv)

  in the case of ζ (rightmost panel), having a mismatch of the type: ζ₁ + ζ₂ ≠ ζ_1, s + ζ_2, s always causes a positive bias because it provides E → B leakage.

It is interesting to observe that the red dashed-dotted line for h, which refers to |h₁ + h₂|< |h_1, s + h_2, s|, corresponds to an overall shift of the spectrum toward higher values, as we are over-correcting for h. The differences between the ideal and realistic $C_{\ell }^{\mathit{BB}}$ can be quantified by estimating the corresponding bias on the tensor-to-scalar ratio r (see Sect. 5.4 for detail). However, we defer this detailed discussion to the more realistic multi-frequency analysis, in the following Sections. Nevertheless, as a qualitative prediction in this simple setting, we would expect, in general, a negative Δr for |h|> |h_s|, a positive (negative) Δr for cos β > ( < ) cos β_s, and a positive Δr for ζ ≠ ζ_s.

5. Multi-frequency case

In this section, we consider a more realistic scenario. In particular, we move from the Dirac-delta response of the detectors in frequency as implicitly assumed in the previous section to a top-hat response. This allows us to take into account two main effects that were previously neglected. First of all, we allowed for HWP matrix elements to be frequency-dependent: we assumed a specific LiteBIRD MHWP design with certain frequency profiles within each band. The MHWP performance was computed using realistic models developed for previous waveplate applications (Pisano et al. 2022).

Secondly, we included a frequency-dependent foreground component in the input sky maps. We describe below the details regarding the inclusion of the two new effects. The analysis follows the same steps detailed in Sect. 3.

5.1. Setup for the multi-frequency analysis

We considered four frequency bands corresponding to the MFT channels of the proposed LiteBIRD satellite (Hazumi et al. 2020) and one frequency band each for LFT and HFT, the closest ones to the CMB channels. The central frequency, band-width, and FWHM of the Gaussian beam for each channel are summarized in Table 2.

The input sky is different from that used in Sect. 4. We only considered one CMB realization from the fiducial spectra chosen by LiteBIRD collaboration (in prep.), with r = 0 and τ = 0.0544. We go on to explain the reasoning behind this later on in this work. We decided to work in μK_CMB units, so that the CMB signal is independent from the frequency in the frequency range that we consider in this work⁵.

We added frequency-dependent maps of foreground emissions to the CMB maps. To do so, we generated a set of foreground maps for the range of frequency used in this work. Foreground maps are generated with the PySM software package (Thorne et al. 2017). We adopted the [d1,s1,a1,f1] model available in PySM, with spatially varying spectral indices of dust and synchrotron. Both the CMB and foreground maps have a resolution of N_side = 512 and are smoothed with a Gaussian beam with FWHM given by the LiteBIRD resolution at the given frequency channel (see Table 2). To keep things simple and to set the focus on the possible chromaticity of systematic effects, in this multi-frequency analysis, we neglected the contribution of instrumental noise. The signal, d_obs(t_i, p, ν), observed at a given time sample, t_i, in a certain pixel, p, in a given frequency channel, ν, is simply the weighted average of the sky signal in that frequency band. Weights are given by the frequency-dependent Mueller matrix elements, M_full(ν), of the optical system.

In general, HWP parameters show a non-trivial dependence on the frequency of the incident signal due to fabrication details. Finite-element modeling and a laboratory characterization of the HWP devices allows reconstruction of the expected profiles of the HWP matrix elements (Pisano et al. 2012, 2014, 2022).

In this analysis, we adopted a model of the HWP derived for the MFT channels (see Figs. 5 and 6; Montier et al. 2020; Lamagna et al. 2020). In a simulated Mesh-HWP the cross polarization parameters ζ are exactly null because of the supposed exact symmetry of the system. Because of this, we are not able to access the frequency profiles of ζ₁, ζ₂ and so, for simplicity, we assumed them to be constant and equal to 10⁻². In reality, the symmetry of the configuration could be spoiled and ζ could be as large as the value we consider. Also for the L/HFT bands we take constant values for all the systematic parameters (for LFT: h = −0.015, β = 7.19° ,ζ = 0.01, for HFT: h = −0.01, β = 15° ,ζ = 0.01). We checked that this does not affect significantly the final result (see Appendix C). The solver matrix B(ν) is built on this model. The pointing matrix A(ν) is a perturbed version of the same model, as we explain later in the text. The Mueller matrix elements of the realistic HWP as given by our model are shown in Fig. A.1 for each frequency channel.

Fig. 5. Simulated profiles of the MFT HWP transmissions h1, h2, for the four selected MFT frequency bands. The subscript s indicates that they are used in the solver matrix B.

Fig. 6. Simulated profiles of the HWP phase-shift β, for the four selected MFT frequency bands of LiteBIRD. The subscript s indicates that they are used in the solver matrix B.

We are thus ready to express the TOD sample d_obs(t_i, p, ν) as a generalized version of Eq. (25), splitting the contribution from CMB (in CMB units) and from the foregrounds (FG, in Rayleigh-Jeans units):

$$\begin{array}{cc}& d_{\mathrm{obs}}(t_i,p,\nu )=\hfill \\ \hfill & \frac{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )[M_i^{\mathit{TT}}(\nu )T_{\mathrm{CMB}}(p)+M_i^{\mathit{TQ}}(\nu )Q_{\mathrm{CMB}}(p)+M_i^{\mathit{TU}}(\nu )U_{\mathrm{CMB}}(p)]}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )}\hfill \\ \hfill & +\frac{\int \mathrm{d}\nu F_{\mathrm{FG}}(\nu )[M_i^{\mathit{TT}}(\nu )T_{\mathrm{FG}}(\nu ,p)+M_i^{\mathit{TQ}}(\nu )Q_{\mathrm{FG}}(\nu ,p)+M_i^{\mathit{TU}}(\nu )U_{\mathrm{FG}}(\nu ,p)]}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )}.\hfill \end{array}$$(31)

The terms F_CMB, FG in Eq. (31) are given by:

$$\begin{array}{cc}\hfill F_{\mathrm{CMB}}(\nu )& =\frac{\partial BB(\nu ,T)}{\partial T_{\mathrm{CMB}}}{\tau }_c(\nu ),\hfill \end{array}$$(32a)

$$\begin{array}{cc}\hfill F_{\mathrm{FG}}(\nu )& =\frac{\partial B{B}_{\mathit{RJ}}(\nu ,T)}{\partial T_{\mathit{RJ}}}{\tau }_c(\nu ).\hfill \end{array}$$(32b)

where τ_c(ν) is the top-hat bandpass, including the throughput factor, ν², and the black-body derivatives take into account the conversion from CMB and RJ units, respectively (see footnote 5). The common denominator in Eq. (31) sets the final units as μK_CMB. We note that in Eq. (31) T_CMB, Q_CMB, U_CMB are independent of frequency and may, therefore, be taken out from the first integral. In matrix form, we can express the following

$$\begin{array}{c}\hfill d_{\mathrm{obs}}(\nu )=A_{\mathrm{CMB}}{m}_{\mathrm{CMB}}+{\displaystyle \int }{A}_{\mathrm{FG}}(\nu )m_{\mathrm{FG}}(\nu )\mathrm{d}\nu ,\end{array}$$(33)

where m_CMB is the (3N_pixels) vector of CMB Stokes parameters and m_FG(ν) is the (3N_pixels) vector of foreground Stokes parameters as a function of frequency. The pointing matrix A_CMB is a frequency-independent (3N_pixels × N_samples) matrix with elements given by the first integral in Eq. (31) over the Mueller elements⁶

$$\begin{array}{cc}& A_{\mathrm{CMB},i}=\hfill \\ \hfill & (\frac{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )M_i^{\mathit{TT}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )},\frac{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )M_i^{\mathit{TQ}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )},\frac{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )M_i^{\mathit{TU}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )}).\hfill \end{array}$$(34)

The subscript i indicates the time sample, as before. The pointing matrix A_FG(ν) is a frequency-dependent (3N_pixels × N_samples) matrix with elements given by the Mueller elements in the second integrand in Eq. (31):

$$\begin{array}{c}\hfill A_{\mathrm{FG},i}(\nu )=(\frac{F_{\mathrm{FG}}(\nu )M_i^{\mathit{TT}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )},\frac{F_{\mathrm{FG}}(\nu )M_i^{\mathit{TQ}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )},\frac{F_{\mathrm{FG}}(\nu )M_i^{\mathit{TU}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )}).\end{array}$$(35)

The mapmaking procedure then consists of inverting Eq. (33) along the same lines as described in Sect. 4. The estimated map is given by:

$$\begin{array}{cc}\hfill m_{\mathrm{out}}& =N^{-1}{\displaystyle \left(\sum _i{B}_i^T{d}_{\mathrm{obs}}(t_i)\right)}\hfill \\ \hfill & =N^{-1}{\displaystyle (\sum _i{B}_i^T{A}_{\mathrm{CMB}},i{m}_{\mathrm{CMB}})}+N^{-1}{\displaystyle (\sum _i{B}_i^T\int A_{\mathrm{FG},i}(\nu )m_{\mathrm{FG}}(\nu )\mathrm{d}\nu )},\hfill \end{array}$$(36)

where

$$\begin{array}{c}\hfill N^{-1}={{\displaystyle \left(\sum _i{B}_i^T{B}_i\right)-1}}^{}\end{array}$$

and again B is the solver matrix integrated over frequency. To build the solver matrix in the multi-frequency case, we adopt the following procedure. We define B as

$$\begin{array}{cc}& B_i=\hfill \\ \hfill & (\frac{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )M_{s,i}^{\mathit{TT}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )},\frac{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )M_{s,i}^{\mathit{TQ}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )},\frac{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )M_{s,i}^{\mathit{TU}}(\nu )}{\int \mathrm{d}\nu F_{\mathrm{CMB}}(\nu )}),\hfill \end{array}$$(37)

where $M_{s,i}^{\mathit{TT}}(\nu ),M_{s,i}^{\mathit{TQ}}(\nu ),M_{s,i}^{\mathit{TU}}(\nu )$ contains our model profiles of the frequency-dependent Mueller matrix elements in Eqs. (34)–(35). As before, the subscript i refers to the time samples. The subscript s labels the matrices entering the map-making matrix B, in parallel with the solver parameters defined in Sect. 4.1. The profiles for h and β are shown in Figs. 5 and 6. In practice, the frequency profile x(ν) of the HWP parameter x = h, ζ, β is used to build the matrix B in Eq. (37), while the perturbed profile x(ν)+Δx(ν) is used to build the matrices A_CMB and A_FG in Eqs. (34)–(35). The perturbation Δx(ν) is treated as a Gaussian fluctuation around the “true” value of x(ν) with a variance of σ_Δx, and can therefore be either positive or negative. An example of the perturbed profiles is shown in Fig. 7.

Fig. 7. Example of perturbations of the HWP profiles: Dashed-dotted lines show one realization of perturbed profiles for h1, h2 with σΔh = 0.001, while the dotted lines show the same realization with higher σΔh = 0.002.

We note that for each parameter, we considered uncorrelated perturbations in the frequency ($\langle \mathrm{\Delta }x(\nu )\mathrm{\Delta }x({\nu }^{\prime })\rangle ={\delta }_{\nu {\nu }^{\prime }}{\sigma }_{\mathrm{\Delta }x}^2$). In addition, perturbations in one parameter are also uncorrelated with perturbations in a different parameter ($\langle \mathrm{\Delta }x(\nu )\mathrm{\Delta }y(\nu )\rangle ={\delta }_{xy}{\sigma }_{\mathrm{\Delta }x}^2$). With this treatment, we aim to simulate the realistic, albeit simplified case of having a mismatch between the model for the profile of each systematic to be used in the solver matrix B as well as in the profiles actually entering in the TOD matrix A. This mismatch is, for example, produced by calibration errors on the parameters. Modeling more realistic error distributions would require a clear knowledge of the optical chain and calibration setup used to perform laboratory measurements. Since this is not available at the current stage, we prefer to defer this topic to a detailed study to future publications. Notwithstanding, the simplified modeling employed in this work remains a valid approach for the pedagogical purposes of this analysis. It is worth noting that even in the ideal case of B = A_CMB, that is, in a perfect calibration of the HWP profiles, the HWP non-idealities coupled with the mapmaking procedure intrinsically produce a distortion of the foreground field that needs to be deprojected (see Sect. 5.2; Vergès et al. 2021).

For each non-ideal parameter h, ζ, β, we considered different choices of σ_Δx, as summarized in Table 5. To avoid oversensitivity to the specific realization of the pointing matrix A, we ran ten simulations for each systematic parameter and each band, resulting in ten realizations of the pointing matrix A. This way, each value of σ_Δx can share the same realizations and the comparison between cases with different σ_Δx thus depends just on the amplitude of the error (within the same band and systematic parameter). To keep matters simple, we used one CMB and FG realization, that is, we kept a fixed input sky. In fact, the residuals coming from systematic effects are mostly dominated by the much brighter foregrounds and not particularly affected by the actual CMB realization. The choice of keeping the same foreground model in the analysis is considered further in Sect. 5.2.

As done in Sect. 4, we vary one parameter at a time. In addition, we also consider the case of joint variation of multiple parameters to test for possible correlations between systematic effects induced by HWP non-idealities.

With regard to the units of Δσ_x in Table 5, we note that in this work, we assume that we are able to reconstruct the HWP profiles with a resolution of Δν = 1 GHz, so that Δσ_x effectively refers to the accuracy on the x-parameter per unit resolution. Had a different frequency resolution been chosen, the corresponding accuracy would be scaled as $1/\sqrt{\mathrm{\Delta }\nu }$.

5.2. Deprojection template

In Sect. 4, residuals were shown with respect to the ideal case of perfect knowledge of the optical system. The reason was that we wanted to focus on the physical effects induced by HWP non-idealities to the observed quantities. In this section, we follow a closer approach to that can be applied to realistic observations. Residuals are shown with respect to a template of our best estimate of the input sky. To build this template, we assume an input CMB sky m_CMB and an input FG model m_FG(ν). We then simulate observations of CMB+FG following the same algorithm described above for multi-frequency observations, with the only difference that we employ the matrix B both as the pointing matrix and the solver matrix. At the end of the mapmaking procedure, we are left with the template:

$$\begin{array}{cc}\hfill m_{\mathrm{templ}}& =N^{-1}{\displaystyle \left(\sum _i{B}_i^T{B}_i{m}_{\mathrm{CMB}}\right)+N^{-1}(\sum _i{B}_i^T\int B_{\mathrm{FG},i}(\nu )m_{\mathrm{FG}}(\nu )\mathrm{d}\nu )}\hfill \\ \hfill & =m_{\mathrm{CMB}}+N^{-1}{\displaystyle (\sum _i{B}_i^T\int B_{\mathrm{FG},i}(\nu )m_{\mathrm{FG}}(\nu )\mathrm{d}\nu ),}\hfill \end{array}$$(38)

where B_FG, i(ν) is the matrix in Eq. (35) with $M_i^{\mathit{TX}}(\nu )=M_{s,i}^{\mathit{TX}}(\nu )$.

The residuals ℛ are then computed as the difference between the output maps (Eq. (36)) and the template, ℛ ≡ m_out − m_templ. The template m_templ is built with the same CMB map and the same FG model used as an input to generate the output maps m_out. In principle, we could have made a different choice, as, indeed, in the real case, it may well be that the estimated sky used to generate a template does not perfectly match the “true” observed sky. This would of course be the cause of differences between the output maps and the template map, even if instrumental systematics were perfectly corrected for. However, it has been proved that an iterative approach can be employed (Planck Collaboration III 2020; Delouis et al. 2019) that converges quickly to the best template estimate, even when the initial sky models are much different from the observed sky. Therefore, in this work, we use the same sky model to generate both the deprojection template and the output maps. In doing so, we rather focus on possible differences between the output maps and the template that are only due to unaccounted HWP systematics and neglect the coupling of component separation with systematics. Here, we include the uncertainties from component separation individually in the noise term that we add to the power spectrum (see Sect. 5.4). By construction, the residual due to unaccounted-for systematics is null if A_CMB = B (TOD HWP = map-making HWP). This requires measuring the systematic parameters with sufficient precision.

We want to stress that the deprojection procedure does not only reduce the residual systematics in case of imperfect knowledge of the HWP profiles, but it also corrects for the intrinsic foreground distortion due to the frequency-dependent non-idealities. Indeed, as is evident in making the comparison between Eqs. (34), (35), and (37), the mapmaking procedure allows for the recovery of unbiased estimates of the CMB component. However, the same is not true for the foregrounds, even if A_CMB = B. The frequency dependence of the HWP parameters introduces a band integration for each sample that not only depends on the direction of observation (as in a usual band-pass integration), but also on the rotation angle of the HWP. This peculiarity of band integration in the presence of a frequency-dependent HWP is particularly dangerous, requiring not only an accurate optical characterization but also HWPs with as flat as possible in-band properties (Vergès et al. 2021).

Finally, we would like to note that the template-subtraction procedure is not intended to serve as a component separation technique. Rather, this procedure has to be intended as a deprojection algorithm that allows for the isolation of the propagation of systematics in the output (or, in real scenarios, observed) maps. As such, this technique has been already proven to be efficient and has been employed in data analysis pipelines for recent CMB experiments (Planck Collaboration III 2020; Ade et al. 2015; Delouis et al. 2019). From the definition of the residuals, ℛ, we can easily see that in our case, a non-vanishing residual is clearly due to the mismatch between the pointing and the solver matrices.

5.3. Results of the multi-frequency analysis

The residual power spectra caused by perturbations in each systematic parameter are presented here. The residual maps, ℛ = m_out − m_templ, computed in each frequency band are masked with a f_sky = 70% galactic mask M₇₀, previously apodized with 5° apodization scale. Then the power spectra are simply computed with the healpy routine anafast, correcting for the sky fraction, f_sky, and the beam window function, b_ℓ:

$$\begin{array}{c}\hfill C_{\ell }^{\mathrm{res}}=C_{\ell }(\mathcal{R}\times M_{70})/(f_{\mathrm{sky}}{b}_{\ell }^2).\end{array}$$(39)

We checked that the use of more sophisticated power spectrum estimators (Alonso et al. 2019) does not make a significant difference in our analysis.

Indeed, the residuals are dominated by the foregrounds. The foreground emission in EE and BB is of the same order of magnitude. The E-to-B leakage due to partial sky coverage (not corrected for when using the anafast estimator) is therefore much weaker than what would be expected in the case of a CMB-only signal. As a result, we can safely neglect the EE − BB mixing when using the anafast estimator.

Furthermore, when applying E − B purification techniques, the estimator may still be biased if the mask apodization procedure is not optimized individually in each multipole bin (Ferté et al. 2013, 2015). Performing this optimization is beyond the scope of this paper, so we decided to employ the simplest estimator available.

In Fig. 8, we report an example of the BB residual power spectra for the 140 GHz frequency band. The residuals clearly follow a power-law behavior. This is expected, since most of the residual power due to the mismatch between the pointing matrix and the mapmaking matrix comes from foregrounds. As already explained in Sect. 5.1, we computed ten realizations of the perturbed profiles for each systematic parameter in each frequency band with the values of σ_Δx in Table 5. We finally computed the residual power spectrum for each σ_Δx as the average of the ten spectra corresponding to each error realization.

Fig. 8. Residual BB power spectra for the frequency band centered at 140 GHz. The thick blue lines are the fiducial $C_{\ell }^{\mathit{BB}}$(CMB) (lensing BB, r = 0). The dashed-dotted lines show the residual spectra $C_{\ell }^{BB,\mathrm{res}}$ obtained as the average over ten Gaussian error realizations of the perturbations applied to each systematic parameter. Different colors correspond to different values of the variance σΔx of the error realizations applied to the parameter x. The dotted line is the noise bias $C_{\ell }^{BB,\mathrm{noise}}$, as described in Sect. 5.4.

We also compute the coaddition of residual maps from different channels, resembling a rough component separation procedure. We implemented it as the weighted average of the maps from each channel:

$$\begin{array}{c}\hfill m_{\mathrm{res},\mathrm{tot}}^x=\frac{{\sum }_i{m}_{\mathrm{res},i}^x{w}_i}{{\sum }_i{w}_i},\end{array}$$(40)

where x = {h, β, ζ}, the sum i runs over the frequency channels, and w_i is the corresponding map of weights obtained from the component separation procedure for the foreground model adopted in Ref. (Poletti, Errard et al., in prep.). The average of value of w_i over the pixels is reported in Table 6: we note that the weights are all positive, as they correspond to CMB channels. The maps $m_{\text{res},i}^x$ have been obtained with values of σ_Δx listed in Table 6. These values correspond to the highest standard deviations among the ones in Table 5 that also satisfy the requirements summarised in Table 7.

In Fig. 9, we show the residual spectra of $m_{\text{res},\text{tot}}^x$. These residuals are lower than those obtained individually from each $m_{\text{res},i}^x$. A similar result is obtained in Sect. 5.4. Here, we anticipate that there is a compensation between residual maps of different channels, even though they share the same error realizations on the frequencies where the bands overlap. From this comparison, we argue that the component separation procedure could reduce the residual caused by mismatches in the HWP systematic parameters, at least when the perturbations to the HWP parameters are uncorrelated in frequency (as we assume in this work).

Fig. 9. Residual BB power spectra from the coadded residual maps. The thick blue line is the fiducial $C_{\ell }^{\mathit{BB}}$(CMB) with r = 0. The dashed-dotted lines are the residual $C_{\ell }^{BB,\mathrm{res}}$ of the coadded maps, see Eq. (40). Different colors correspond to residual spectra due to a different class of HWP systematic parameter. The dotted line is the noise spectrum $C_{\ell }^{BB,\mathrm{noise}}$, as described in Sect. 5.4.

5.4. Propagation to the tensor-to-scalar ratio

In the multi-frequency case, we quantify the impact of HWP non-idealities in terms of a possible bias on the estimate of the tensor-to-scalar ratio r. To evaluate the bias,

we compute the likelihood of the output spectra:

$$\begin{array}{c}\hfill {\stackrel{\sim }{C}}_{\ell }^{\mathrm{BB}}=C_{\ell }^{BB,\mathrm{fid}}+C_{\ell }^{BB,\mathrm{res}}+C_{\ell }^{BB,\mathrm{noise}},\end{array}$$(41)

where $C_{\ell }^{BB,\mathrm{fid}}\equiv C_{\ell }^{BB,\mathrm{lensed}}$ is the fiducial BB power spectrum (lensing only, r = 0), $C_{\ell }^{BB,\mathrm{res}}$ is the BB power spectrum of the residual map ℛ as computed in Sect. 5.2 and $C_{\ell }^{BB,\mathrm{noise}}$ is the noise spectrum due to foreground residual from component separation and instrumental noise for LiteBIRD (see Poletti, Errard et al., in prep., for details on how this spectrum was obtained). Hereafter, we drop the superscript BB for simplicity. The residual $C_{\ell }^{\mathrm{res}}$ is treated as if it were a spurious cosmological signal leading to a bias in the estimate of r.

Next, we adopt the exact likelihood distribution (Gerbino et al. 2020; Hamimeche & Lewis 2008):

$$\begin{array}{cc}\hfill -2\ln {\stackrel{\sim }{\mathcal{L}}}_i(r)& =-2\ln \mathcal{L}({\stackrel{\sim }{C}}_{\ell ,i}|C_{\ell }(r)+C_{\ell }^{\mathrm{noise}})\hfill \\ \hfill & =f_{\mathrm{sky}}\sum _{\ell }(2\ell +1)[\frac{{\stackrel{\sim }{C}}_{\ell ,i}}{C_{\ell }(r)+C_{\ell }^{\mathrm{noise}}}-\ln \left(\frac{{\stackrel{\sim }{C}}_{\ell ,i}}{C_{\ell }(r)+C_{\ell }^{\mathrm{noise}}}\right)],\hfill \end{array}$$(42)

where the i = 1, …, n index stands for one specific realization of the B_i solver matrix, ${\stackrel{\sim }{C}}_{\ell }$ is the observed power spectrum, and $C_{\ell }(r)=C_{\ell }^{\mathrm{lensed}}+C_{\ell }^{\mathrm{tens}}(r)$ is the theoretical BB power spectrum for a given value of r. The likelihood analysis is restricted to the multipole range 2 ≤ ℓ ≤ 200 of interest for LiteBIRD (Hazumi et al. 2020). The log-likelihoods are then averaged over all the i realizations and renormalized to the peak of the distribution. We use a flat prior on r.

We define ${\stackrel{\sim }{p}}^{\mathrm{res}}(r)$ as the posterior distribution corresponding to the likelihood averaged over all the error realizations. The bias on r due to HWP systematics, Δr, is quantified as the maximum probability value⁷ of the posterior distribution: Δr = r_{peak(Δsyst.)}.

As a final remark, we note that the residual map is computed as the difference with respect to the template map, containing the input CMB map m_CMB. The latter leaves a distorted signal in the residual map that is in principle dependent on the actual CMB realization. However, this contribution is completely negligible with respect to the foreground-induced residual. For this reason, we can safely neglect the scatter due to the cosmic variance and add the residual power spectrum directly to the fiducial one.

5.5. Requirements on the sensitivity for the systematics from the bias Δr

The bias Δr is estimated for each value of σ_Δx for each systematic parameter and frequency band. We expect a quadratic relation between the bias on r and the variance of the error realization since $\mathrm{\Delta }r\propto C_{\mathcal{l}}\propto {\sigma }_{\mathrm{\Delta }x}^2$. Indeed, if we plot Δr vs. σ_Δx, we find that a quadratic fit works well, especially for smaller error variance. In Fig. 10, the bias on r due to the perturbation of β is larger than the bias coming from the other two parameters (h, ζ). This is due to the fact that we are considering larger σ_Δβ because of a wider dynamical range of β in the HWP model profiles. In the same figure, the error on Δr is reported as ${\sigma }_{\mathrm{\Delta }r}/\sqrt{10}$, that is the standard deviation of Δr_i from each i = 1, …, 10 error realization divided by the square root of the number of realizations. We derive the accuracy requirements on each systematic parameter in each band so that Δr ≲ 10⁻⁵. This threshold is set as 1% of the expected sensitivity on r from LiteBIRD, that is, σ_r ∼ 10⁻³. The accuracy requirements are quoted in Table 7.

Fig. 10. Quadratic fit of the relation Δr  −  σΔx, for the band centered at 140 GHz. The green solid line marks the threshold Δr = 10−5, which we have set to derive a requirement on the highest tolerable σΔx when we perturb only one systematic at a time.

We also checked the effects of letting two systematic parameters at a time to be perturbed in the analysis. In this case, we find that the bias on r is approximated by the sum of the bias induced individually in the case when one parameter is perturbed, that is, Δr_x, y ≃ Δr_x + Δr_y, for x, y = {h, β, ζ}. To allow for a more robust check, we increased the number of error realizations per parameter to 200 for a single band. We first generated maps in which one kind of parameter was perturbed at a time and computed the corresponding averaged Δr_x. Then we used the same error realizations to generate maps in which two classes of parameters x, y are jointly perturbed. We derived the corresponding Δr_x, y and compared it against the sum Δr_x + Δr_y, finding no significant difference. This amounts to having no clear correlation between different parameters within the error $\frac{{\sigma }_{\mathrm{\Delta }r}}{\sqrt{200}}$. This could be due to the fact that we are perturbing the model profiles in a way that is independent both on frequency and on the systematic parameters. A more complex and realistic modeling of the perturbations could change this result. For example, we could allow for some degree of correlation between errors on different systematics, which might be the case if measurements of the HWP parameters were performed at the same time. However, this choice would require a realistic modeling of how measurements are performed, which goes beyond the scope of this paper.

Finally, we estimate the bias on r caused by the weighted average of residual maps from each channel (Eq. (40)). We obtain a $\mathrm{\Delta }{r}_{\text{tot}}^x$ which is always smaller than the weighted average $\mathrm{\Delta }{r}_{\mathrm{weight}}^x=\frac{{\sum }_{j=1}^{N_{\mathrm{pixels}}}({\sum }_i{\mathit{w}}_i^2\mathrm{\Delta }{r}_i^x/{\sum }_i{\mathit{w}}_i^2)}{N_{\mathrm{pixels}}}$ of the $\mathrm{\Delta }{r}_i^x$ corresponding to each $m_{\text{res},i}^x$ (in Table 6):

$$\begin{array}{cc}\hfill \mathrm{\Delta }{r}_{\mathrm{tot}}^h=1.3\times {10}^{-6}& <\mathrm{\Delta }{r}_{\mathrm{weight}}^h=3.9\times {10}^{-6}\hfill \\ \hfill \mathrm{\Delta }{r}_{\mathrm{tot}}^{\beta }=1.2\times {10}^{-6}& <\mathrm{\Delta }{r}_{\mathrm{weight}}^{\beta }=8.3\times {10}^{-6}\hfill \\ \hfill \mathrm{\Delta }{r}_{\mathrm{tot}}^{\zeta }=2.6\times {10}^{-6}& <\mathrm{\Delta }{r}_{\mathrm{weight}}^{\zeta }=5.9\times {10}^{-6}.\hfill \end{array}$$(43)

We weight the biases with $w_i^2$, as Δr ∝ C_ℓ, which is quadratic in the map. It is possible that a non-parametric component separation procedure would relax the requirements shown in Table 7.

6. Conclusions

In this work, we study the impact of non idealities of the half-wave plate (HWP) in the context of future cosmic microwave background (CMB) observations, focusing on the case of a LiteBIRD-like satellite mission. We consider the following classes of non-idealities: departure from unitary transmission (h), spurious phase shift (β), and cross-polarization (mixing of orthogonal polarization components, ζ). Any mismatch between the measured properties of the HWP and the actual properties that enter in the construction of the time-ordered data (TOD) during observations can propagate throughout the analysis pipeline and bias the final science products down to the tensor-to-scalar ratio, r, estimate. We have first presented at length the formalism describing how light propagation is affected by a non-ideal rotating HWP. We have developed an agile simulation suite to quickly reproduce the LiteBIRD scanning strategy and find the relative on-the-fly mapmaking solution, with little computational cost. To do so, we considered a simplified scenario where light propagates with normal incidence through the optical system (including the HWP) and is collected by a single pair of polarization-sensitive detectors at boresight. Because of this simplified setting, we are not able to capture some additional systematic effects, such as the HWP synchronous signal, which has been observed in different experiments (Ritacco et al. 2017; Johnson et al. 2007; Kusaka et al. 2014).

The optical system has been described in the Jones formalism and we have also shown the conversion to the Mueller formalism. The full expressions of the Mueller matrix elements of the non-ideal (not rotating) HWP are presented in Appendix A.

We first focused on the case of an input CMB-only sky observed at a single frequency. This is motivated by the fact that we wanted to single out the effects of HWP non-idealities on the reconstructed CMB spectra while neglecting any other source of contamination (e.g., color-correction due to bandpass integration). We have shown results obtained in two scenarios: (a) in the case in which a mismatch persists between the HWP parameters entering the TOD and the ones used in the map-making solution; (b) in the case in which the two set of parameters are identical, albeit they are still non-ideal. As expected, our results show that scenario (b) minimizes the propagation of HWP-induced systematic effects to CMB spectra.

We then moved on to a more realistic study with a frequency-dependent input signal (including also foregrounds) modulated by a frequency-dependent HWP profile. We have considered the four MFT frequency bands of LiteBIRD centered at [100, 119, 140, 166] GHz, and the closest LFT/HFT bands to the CMB channels (centered at 100/195 GHz; see Table 2). We assumed a top-hat bandpass profile, for simplicity. In this multi-frequency study, we have only focused on the case in which the profile of the TOD HWP does not match the profile of the mapmaking HWP. We adopted simulated frequency profiles for the departure from unitary transmission, h, and the non-ideal phase shift β in the MFT frequency bands, provided by finite-element simulations of the MHFT LiteBIRD MHWPs. The waveplate designs are based on previous developments and realizations (Pisano et al. 2022).

The profiles for ζ are always fixed to a realistic (Pisano et al. 2012) level of 0.01 in all the bands, as this parameter was not included in the suite of simulations at our disposal. In the LFT (also not included in the simulation suite) and HFT bands, we used constant profiles also for h and β. To simulate a mismatch between the TOD HWP and the map-making HWP, all the profiles were perturbed with Gaussian-distributed errors, uncorrelated both in frequency and among the different parameters. We noted that this simple procedure allow our results to be basically independent from the initial shape in frequency. This justifies our choice a posteriori of fixing the parameters for the LFT/HFT bands, as well as the value of ζ in each channel, to a constant value.

In this multi-frequency study, a template map obtained from the observation of the same input sky (CMB and foregrounds) with an ideal HWP was deprojected from the realistic output maps to obtain maps of residuals. The template has been generated with exactly the same foreground model adopted for the input maps. This is equivalent to assuming a perfect knowledge of the foreground sky. Of course, this may not be the case with actual observations – however, it is justified by the need to not include uncertainty on the foreground modeling on top of the effect of the systematics in the residual maps. Our assumption guarantees that by construction, the residual maps vanish when the TOD HWP parameters perfectly match those associated with the mapmaking.

In the multi-frequency case, the CMB spectra have been extracted from the residual maps after applying a galactic mask (f_sky = 70%). These residual spectra, in addition to foreground residual from component separation and instrumental noise for LiteBIRD, have been fed to an exact likelihood to quantify their induced bias Δr on r with respect to the fiducial estimate obtained in absence of systematics residuals. The bias has been quantified for each class of systematic effects individually and in each individual frequency band. By imposing that Δr ≤ 10⁻⁵ (1% of the expected sensitivity on r from LiteBIRD), we set a requirement on the accuracy needed on each HWP parameter in each band.

We repeated the analysis by allowing for pairs of non-ideal parameters to be perturbed simultaneously to check for correlated effects between classes of non-idealities. We found that the bias Δr from a joint variation is consistent with the sum of the biases corresponding to perturbing each of the two parameters at a time, with the same error. This is enough to exclude significant correlations between non-ideal HWP parameters given our experimental setup. In fact, our assumption that Gaussian perturbations fully capture the error in the measurement of HWP parameters is likely to be unrealistic. For example, it is possible to have errors that are correlated within the same frequency band. In addition, errors on different parameters might be correlated if their measurements are simultaneous. To implement this kind of perturbation scheme, we would need a realistic model of how measurements of HWP properties are performed. We defer this study to a future work.

We also provide the results of the coaddition of residual maps from the different frequency channels. We find a general reduction of Δr for the final coadded map. This could point to the fact that a component separation procedure might mitigate the impact of HWP non-idealities thanks to a cancellation among the frequency channels – at least in our setup.

Nonetheless, allowing for a correction of these systematic effects in the mapmaking process remains key to mitigating their impact on science products. We showed that considering an ideal HWP in the map-making procedure could lead to Δr ≈ 𝒪(10⁻³ − 10⁻²), depending on the amplitude of each systematic parameter. Furthermore, some calibration procedures could be attempted for the parameters h and β, which mediate the polarization efficiency (see Eqs. (28), (29)). In Appendix B, we showed that ζ behaves similarly to a rotation of the polarization angle (see Appendix B), and could be thus reabsorbed in the calibration of the latter. Some complications could arise from the frequency dependence of those parameters, however. When this work was in preparation, a study by Duivenvoorden et al. (2021) was published on the same topic, and we would like to point out that our analysis nicely complements the findings of those authors. Here, we offer a pedagogical approach to the use of a HWP in CMB experiments. We also provide a thorough discussion of the complementarity between the Jones and Mueller formalisms in the context of CMB polarimetry. Finally, we highlight a significant, and possibly problematic, effect: the fact that the in-band variation of the properties of a non-ideal HWP can affect the observed signal and the reconstructed sky maps by introducing an effective band integration that depends also on the HWP rotation angle (see Sect. 5.2). This effect can potentially lead to a direction-dependent bandpass mismatch.