Impact of half-wave plate systematics on the measurement of CMB B-mode polarization

Polarization of the cosmic microwave background (CMB) can help probe the fundamental physics behind cosmic inflation via the measurement of primordial B modes. As this requires exquisite control over instrumental systematics, some next-generation CMB experiments plan to use a rotating half-wave plate (HWP) as polarization modulator. However, the HWP non-idealities, if not properly treated in the analysis, can result in additional systematics. In this paper, we present a simple, semi-analytical end-to-end model to propagate the HWP non-idealities through the macro-steps that make up any CMB experiment (observation of multi-frequency maps, foreground cleaning, and power spectra estimation) and compute the HWP-induced bias on the estimated tensor-to-scalar ratio, r. We find that the effective polarization efficiency of the HWP suppresses the polarization signal, leading to an underestimation of r. Laboratory measurements of the properties of the HWP can be used to calibrate this effect, but we show how gain calibration of the CMB temperature can also be used to partially mitigate it. On the basis of our findings, we present a set of recommendations for the HWP design that can help maximize the benefits of gain calibration.

Inflation sources initial conditions for cosmological perturbations via primordial vacuum quantum fluctuations [23][24][25][26].The relative amplitude of the resulting scalar and tensor perturbations is quantified in terms of the tensor-to-scalar ratio, r.Since tensor perturbations [27,28] would leave a distinct B-mode signature on the CMB polarization [29][30][31][32], r can be inferred from the angular power spectrum of the primordial B modes.To date, CMB observations have only placed upper bounds on r, the tightest being r < 0.032 (95% CL) [33] (see also [11,34,35]).Future surveys aim for unprecedentedly low overall uncertainties, which, depending on the true value of r, would lead to a detection or a tightening of the upper bounds, both of which would allow us to place strong constraints on inflationary models.
In this paper, we present a simple framework to propagate the HWP non-idealities through the three macro-steps that characterize any CMB experiment: observation of multifrequency maps, foreground cleaning, and power spectra estimation.We exploit the simplicity of the harmonic internal linear combination (HILC) foreground cleaning method [56] to keep the treatment semi-analytical.This choice, along with our working assumptions, makes the analysis computationally inexpensive1 and reflects our intention to develop an intuitive understanding of how the HWP affects the observed CMB.
The remainder of this paper is organized as follows.In section 2 we generalize the arguments presented in [49] and provide a simple model for multi-frequency maps observed through a rapidly spinning HWP.We then introduce the HILC foreground cleaning method and present the procedure we will use to infer r.In section 3, we discuss the specific choices we make to model sky, noise, and beams, and present the results of the analysis in two cases.First, we assume that the HWP is ideal and verify that the pipeline recovers the input CMB signal.Second, we consider LiteBIRD-like instrument specifics and assume realistic HWPs.We find that, for our choice of HWPs and r true = 0.00461 in input, the HWP non-idealities introduce an effective polarization efficiency that suppresses the polarization signal, resulting in r = (4.30+0.56  −0.53 )×10 −3 .We also show how including gain calibration of the CMB temperature in the map model can partially mitigate this effect.In section 4, we derive a set of design recommendations that can help maximize the benefits of the gain calibration step.We also review the simplifying assumptions underlying the model and briefly discuss how they might be relaxed.Conclusions and perspectives are presented in section 5.

Mathematical framework
In this section we present a simple model for multi-frequency maps observed through a rapidly spinning HWP.We also introduce the HILC foreground cleaning method and derive an explicit expression for the B-mode angular power spectrum of its solution, C BB ℓ,hilc , given the modeled multi-frequency maps.Finally, we present the methodology we use to estimate the tensor-toscalar ratio parameter, r, from C BB ℓ,hilc .

Modeling the observed maps
We describe linearly polarized radiation2 by the Stokes I, Q and U parameters defined in right-handed coordinates with the z axis taken in the direction of the observer's line of sight (telescope boresight), according to the "CMB convention" [66].Given an incoming Stokes vector S ≡ (I, Q, U ), the effect of a polarization-altering device on S can be described by a Mueller matrix M, so that S ′ = MS [67].Assuming azimuthally symmetric and purely co-polarized beams, we can approximate the entire telescope's optical chain by means of a Mueller matrix acting on appropriately smoothed input Stokes parameters.This setup allows us to write the telescope response matrix 3 , A, analytically, and to obtain simple expressions for both time-ordered data (TOD), d, and binned maps, m [68]: where m denotes the pixelized {I, Q, U } sky maps smoothed to the resolution of the instrument, n the noise contribution to the TOD, and Â the response matrix assumed by the map-maker.
If the telescope's first optical element is a rapidly rotating HWP with Mueller matrix the maps reconstructed from the TOD of the i channel's detectors by an ideal binning mapmaker that assumes where the sum over λ spans different sky components (CMB, dust, and synchrotron emission), the integral represents a top-hat bandpass with a bandwidth of ∆ν i ≡ ν i min − ν i max , the superscript i in m i λ stresses that the input map is smoothed with the beam of the frequency channel i, and n i denotes the noise maps.
Eq. (2.3) approximates the observed maps well when the cross-linking is good, that is, when each sky pixel is observed with a variety of scan angles.This condition is ensured by the rapid HWP rotation and the good LiteBIRD sky coverage, which guarantee that the scan angles are sampled uniformly enough for each pixel [49].As a consequence, our model neglects intensity-to-polarization leakage, the effects of which have been shown to be correctable [55].
If we also make the simplifying assumption that the spectral energy distribution (SED) of each component is uniform throughout the sky, we can rewrite each sky map as m λ (ν) ≡ a λ (ν)m λ (ν * ), where ν * is some reference frequency.This is equivalent to using the s0d0 option in the Python Sky Model (PySM) package [69], which has often been used in the literature for the study of systematics (e.g., [70,71]).The reason for this assumption is twofold.First, it is often useful to separate the effects of systematics from the complexity of the foreground emission.Second, as shown in [70], the study of systematics is strongly influenced by the specific class of component separation methods, that is, whether it is a blind method, such as HILC [56], or a parametric method, such as FGbuster [72].In this paper, we use HILC and leave the study based on a parametric method for future work.
The factorization, m λ (ν) = a λ (ν)m λ (ν * ), allows us to rewrite eq.(2.3) as where we have dropped the ν * dependence for the sake of simplicity and defined ) The coefficients in these equations have a clear physical interpretation: g i λ is an effective gain for the temperature data, ρ i λ and η i λ are effective polarization gain (or polarization efficiency) and cross-polarization coupling, respectively, caused by the non-idealities of the HWP.
Including photometric calibration Photometric calibration is a crucial step in any CMB analysis pipeline that allows us to map the instrumental output to the incoming physical signal [73].Here, we assume that the CMB temperature dipole [74,75] is used as a calibrator, as is commonly done in CMB experiments, and we neglect any imperfections in calibration.In other words, we assume to know gi = g i CMB exactly after calibration.The photometrically calibrated counterpart of eq.(2.4) reads Spherical harmonics coefficients To apply the HILC method to the modeled maps, we expand eq.(2.6) in spin-0 and spin-2 spherical harmonics and write the corresponding B-mode spherical harmonics coefficients as where a E ℓm,λ and a B ℓm,λ are the Eand B-mode coefficients of the unsmoothed maps at some reference frequency ν * (implicit here), and B i ℓ is the beam transfer function of the channel i.

Harmonic internal linear combination
The internal linear combination (ILC) [76] is a blind foreground cleaning method.It can be implemented in both map and multipole space, the latter case being referred to as HILC [56].
Given the spherical harmonics coefficients, a X,i ℓm with X = (T, E, B) and i ∈ {1, . . ., n chan }, of the maps observed by each of the n chan frequency channels, the HILC solution is given by [56] a X ℓm,hilc = where e is a column vector with n chan elements all equal to one, and C ℓ is the n chan × n chan covariance matrix of the observed maps: C ij ℓ = ⟨a i * ℓm a j ℓm ⟩.By construction, the weights minimize the variance of the final map and add to unity, ∑ i w i ℓ = 1, preserving the frequency independence of the CMB black-body spectrum.However, the frequency dependence of g i CMB , ρ i CMB , and η i CMB can violate this sum rule.This is the main point we study in this paper.
Modeling the HILC solution To apply the HILC to the analytical predictions discussed in section 2.1, we could simply use eq.(2.7); however, since different channels are characterized by different beams, it is preferable to perform the HILC on unsmoothed spherical harmonic coefficients, a i ℓm ≡ âB,i ℓm /B i ℓ and write the covariance matrix as We use eq.(2.9) to compute the HILC weights, w ℓ , and the spherical harmonics coefficients of the HILC solution according to eq. (2.8).The corresponding angular power spectrum reads (2.10)This is the main equation from which we derive all of our results.
Even at this early stage, we can make some educated guesses about which terms will contribute the most to the final angular power spectrum.By construction, the HILC tries to select the component λ whose ρ i λ and/or η i λ are nearly constant across all frequency channels, i.e., a black-body spectrum.For example, if m qq (ν) − m uu (ν) or m qu (ν) + m uq (ν) depended on frequency as the inverse of the SED of the foreground emission, the foreground would leak into the HILC solution.However, the Mueller matrix elements of realistic HWPs do not exhibit such behavior.We therefore expect foreground-to-CMB leakage to be small in the final angular power spectrum.
Focusing on the CMB, eq.(2.10) tells us that there are two potential contaminations: E-to-B leakage, which can occur if the effective cross-polarization coupling, η i CMB , is nearly constant across the frequency channels, and suppression of the B modes, which is instead driven by the effective polarization efficiency, ρ i CMB .The relative importance of these effects depends on the specific design choice of the HWP.

Maximum likelihood estimate of the tensor-to-scalar ratio
The modeled angular power spectrum is where C GW ℓ is the primordial B-mode power spectrum with r = 1 [31,32], C lens ℓ is the lensed B-mode power spectrum [77], A lens is its amplitude with A lens = 1 being the fiducial value, and N BB ℓ is the HILC solution for the total noise power spectrum [the last term in eq.(2.10)].The probability density function (PDF) of the observed B-mode power spectrum for a given value of r and A lens , P (C BB l,obs | r, A lens ), is given by (e.g., [78]) where f sky is the sky fraction used to evaluate C BB ℓ,obs .We use f sky = 0.78, for which our sky model is defined (see table 1 for details).Given the PDF, the likelihood function is We use ℓ max = 200, which is the fiducial value for LiteBIRD [16].Using Bayes' theorem, the posterior PDF of r with A lens marginalized over a flat prior is The frequentist profile likelihood is given instead by maximizing the bidimensional likelihood with respect to A lens for a set of values {r 0 , . . ., r n } Regardless of whether L(r) ≡ L m (r) or L(r) ≡ L p (r) is chosen, we define r as the maximumlikelihood estimate (MLE), i.e., the value of r that maximizes L(r).We compute the corresponding uncertainty as [78] where L(r) is normalized as Eq. (2.15) defines the variance associated with a Gaussian random variable, which is characterized by a likelihood that is symmetric with respect to its maximum.More generally, however, L(r) may be asymmetric, and we estimate uncertainties as asymmetric 68% CL intervals.

Analysis
We apply the framework presented in section 2 to extract the bias on r caused by a particular choice of HWP design.Given M hwp , our code5 performs the following steps: 1. Compute the covariance matrix, C B,ij ℓ , as in eq.(2.9), 2. Invert C B,ij ℓ to obtain the HILC weights, w i ℓ , as in eq.(2.8), 3. Use the w i ℓ to compute the BB spectrum of the HILC solution, C BB ℓ,hilc , as in eq.(2.10), 4. Compute the two-dimensional likelihood L(r, A lens ) from C BB ℓ,hilc , according to eq. (2.13), 5. Obtain the one-dimensional posterior PDF, L m (r), by marginalizing over A lens , and the profile likelihood, L p (r), by maximization, 6.Return r and σ r , defined as in eq.(2.15), computed from L m (r) and L p (r).
To validate our end-to-end model and code, we first perform the analysis for an ideal HWP and then move on to more realistic cases.However, before presenting our results, we review the additional assumptions that go into the explicit computation of the HILC covariance matrix C B ℓ , with the exception of the HWP choice.CMB, dust and synchtrotron spectral responses For maps in thermodynamic units, the a λ (ν) functions entering in eqs.(2.5) read (see appendix A for a complete derivation)  [80].Right panel: The power-law parameters for the angular power spectra of synchrotron and thermal dust emission entering in eq.(3.2) as reported in [80] for the Commander [81] analysis with f sky = 0.78.
where B ν (T ) denotes a black-body spectrum at temperature T , x ≡ hν/(k B T 0 ) and T 0 = 2.725 K is the average temperature of the CMB [79].The values of the remaining parameters entering in eqs.(3.1) are specified in table 1.
CMB, dust and synchtrotron angular power spectra The CMB angular power spectrum is computed with CAMB [82] assuming the best-fit 2018 Planck values for the cosmological parameters [3], except for the tensor-to-scalar ratio, which is set to r true = 0.00461.This is the same fiducial value as assumed in [16], and corresponds to Starobinsky's R 2 inflationary model [83] with the e-folding value of N * = 51.
As for the polarized foreground emission, we parameterize their angular power spectra as a power law [80] Specific values of the parameters are reported in table 1 for both dust and synchrotron.Note that we neglect any intrinsic EB correlation in the input, which is inaccurate (polarized dust emission has been observed to have non-zero T B correlation [84,85], which implies the presence of a EB correlation [86,87], and cosmic birefringence [22] would also result in a non-zero EB).When presenting our results in section 3.2, we comment on this assumption and argue that allowing non-zero EB in input would not dramatically affect the analysis.
Instrument specifics To simulate LiteBIRD's design, we consider an instrument that mounts three different telescopes at low (LFT), medium (MFT), and high frequency (HFT).
The specific frequency ranges of each telescope and frequency channel are taken from [16].
Noise covariance matrix Using a rotating HWP as polarization modulator suppresses the polarized 1/f noise component [36].Being left with white noise only, we parameterize where n i p is the noise in Stokes parameters Q or U per pixel with solid angle Ω pix = 1 arcmin 2 .The specific values assumed for each n i p are taken from [16].
Beams Since we assume the beams to be Gaussian and perfectly co-polarized, the B i ℓ coefficients only depend on the beam's full width at half maximum (FWHM).Specific FWHM values for each channel are taken from [16].

Validation: ideal HWP
An ideal HWP is described by a frequency-independent Mueller matrix with elements In this case, the coefficients g i λ and ρ i λ reduce to the average of the correspondent a λ (ν) function over the band i [eq.(2.5)], which we will denote a i λ .The η i λ coefficients go instead to zero.According to eq. (2.6), the multi-frequency maps reduce to While the CMB component is not affected by the presence of the ideal HWP, the foreground emission suffers from a color correction, and the noise term is rescaled channel-by-channel.In this simple situation, the HILC should perform well and recover the CMB signal plus some noise bias given by We should therefore check that, for M ideal = diag(1, 1, −1), the HILC output is in good agreement with the input CMB angular power spectrum, once the noise bias is removed.
In figure 1, we show the angular B-mode power spectrum of the HILC solution, together with the input angular power spectra of CMB, dust, and synchrotron.For completeness, we also show the foreground residual and the noise bias.The noise bias has been removed from both the HILC solution and the foreground residual.The agreement between the HILC solution and the input CMB power spectrum is excellent up to ℓ ≃ 325, roughly corresponding to LiteBIRD's beam resolution.
In figure 2 we show the HILC weights for the three telescopes.All MFT channels have positive weights, consistent with them being CMB channels.On the other hand, some of LFT and HFT channels (at very low and very high frequencies, respectively) have negative weights, resulting in foreground subtraction.
The code returns the MLE r = (4.64+0.57−0.54 ) × 10 −3 , which is compatible with the fiducial value of r true = 0.00461, as the bias ∆r ≡ r − r true = 0.03 × 10 −3 is a small fraction of the uncertainty.Similarly, A lens is also unbiased: Âlens = 1.00 ± 0.01.This is what we expect, given the good agreement between the debiased HILC solution and the input CMB shown in figure 1.
To test that the good agreement between the estimated r and r true is not just due to the specific value chosen for r true , we repeat the analysis for a sample of the currently allowed values of r true .The bias remains a small fraction of the error bar for all the values considered.In particular, for r true = 0 the best fit is r = 0 with 68% C.L. upper bound 0.00016.

More realistic HWPs
For this analysis, we consider more realistic HWPs for each telescope.For LFT, we consider the Pancharatnam-type multi-layer sapphire symmetric stack design described in [88], provided with an anti-reflection coating (ARC) as presented in [89].For the metal-mesh HWPs of MFT and HFT, we use the same input simulations and working assumptions as in [50].For an ideal HWP, the rescaled angular power spectrum, D BB ℓ , of the HILC solution (dashed teal line) overlaps the input CMB spectrum (black solid line) for a wide range of multipoles.For large ℓ, the two spectra begin to diverge as we approach the instrumental resolution.This can be seen by looking at the dotted gray line, representing the residual noise, which intersects the input spectrum at ℓ ∼ 325.For completeness, we also plot the input dust and synchrotron D BB ℓ (orange and yellow, respectively) and the foreground residual (red dotted line).The noise bias has been removed from both the HILC solution and the foreground residual spectra.The w i ℓ weights corresponding to the HILC solution are shown in figure 2. HILC weights, w i ℓ , for each of the three telescopes with an ideal HWP.In each plot, different colored lines correspond to different frequency channels: from purple to yellow for lower to higher central frequencies (see [16] for details on the channels' specifics).The corresponding BB angular power spectrum is shown in figure 1 (dashed teal line).
We manipulate each set of Mueller matrices by performing a rotation of the angle θ t that minimizes the integral For LFT, we consider a symmetric stack design [88] provided with ARC [89], compute its Mueller matrix elements, and rotate them of 55.02 ○ , to express them in a reference frame with the x axis parallel to the HWP optic axis.Instead, the Mueller matrix elements for MFT and HFT are obtained by following the same procedure and input simulations as done in [50], and rotating them of 0.29 ○ .The dashed gray lines represent the ideal values of each element.
over the entire frequency band of each telescope, specified by t = {l, m, h}.This choice is ultimately motivated by the specific design we assume for LFT, since there is no unique way to determine the position of the HWP's optical axes for a symmetric stack.Rotating M hwp,l of θ l then amounts to calibrate the HWP Mueller matrix and express it in a coordinate system aligned with the optical axes.Instead, the HWPs of MFT and HFT employ meshfilter technology [90], for which optical axes can be more easily identified.However, for the sake of consistency, we choose to perform analogous rotations on the Mueller matrices of MFT and HFT metal-mesh HWPs.Rotation angles that minimize eq.(3.7) are 55.02 ○ for LFT and 0.29 ○ for M-HFT.The rotated Mueller matrix elements of each HWP are shown as a function of frequency in figure 3.
Given the elements of the Mueller matrix, we compute the coefficients ρ i λ and η i λ according to eq. (2.5) and repeat all the steps outlined at the beginning of section 3. The HILC solution, D BB ℓ,hilc , is shown in figure 4.Although the foreground residual (red dotted line) shows more features than in the ideal case of figure 1, its contribution to D BB ℓ,hilc is still subdominant.This confirms our intuition that reasonably optimized HWPs do not cause strong foreground leakage in the HILC solution [see the discussion below eq.(2.10)].Note that, given the negligible foreground leakage, taking C EB ℓ,dust = C EB ℓ,synch = 0 in input is not such a strong assumption.Even if we allowed non-zero EB correlations, they would not contribute Compared to the ideal HWP case shown in figure 1, the non-ideal HILC solution slightly differs from the input CMB at low multipoles.For comparison, we also show the residual noise bias (dotted gray line) and the foreground residual (red dotted line).They both show more features than their counterparts in figure 1.The w i ℓ weights corresponding to the HILC solution are shown in figure 5. significantly to the HILC solution.
In figure 5 we also show the HILC weights for the three telescopes.The weights look qualitatively similar to their ideal counterparts shown in figure 2.
To give more precise considerations, figure 6 shows the power spectra on large angular scales in more detail.We show the two independent terms that contribute to D BB ℓ,hilc component-by-component: ρ-only (polarization efficiency) and η-only (cross-polarization coupling).These were obtained using the full covariance matrix C ℓ given in eq.(2.9) to compute the HILC weights, while neglecting some of the terms entering in eq.(2.10).For instance, the ρ-only dust contribution reads Intuitively, it makes sense for the effective polarization efficiency component to dominate in the CMB contribution.While η i CMB can be both positive and negative, all ρ i CMB are constrained to be smaller than 1.This means that, while the average ⟨η i CMB ⟩ across all frequency channels can be close to zero, ⟨ρ i CMB ⟩ cannot be arbitrarily close to 1.The HILC, which looks for the solution that minimizes the variance, may then be able to get rid of all crosspolarization coupling, while it cannot undo the average suppression due to the polarization efficiency.As a consequence of the smallness of the cross-polarization coupling component relative to the polarization efficiency, we argue that relaxing the C EB ℓ,CMB = 0 assumption for the input spectra would not significantly change our results.
Interestingly, the HILC solution approximately satisfies with 10 −5 relative tolerance and 10 −8 absolute tolerance for a wide range of multipoles, 25 ≤ ℓ ≤ 372.The upper limit has a simple interpretation: it roughly corresponds to the instrumental resolution.
Bias on the tensor-to-scalar ratio We finally employ the methodology introduced in section 2.3 to propagate the small discrepancy between the input CMB and the HILC solution shown in figure 4 into a bias on r.We compare the marginalized posterior PDF, L m (r), with .Normalized profile likelihood, L(r) = L p (r), obtained from the HILC solution, ĈBB ℓ,hilc , given the HWP specifics presented in section 3.2 (teal solid line).The likelihood has a maximum at r = 0.00430.The shaded region identifies the 68% CL interval, and goes from r −0.00053 to r +0.00056 .The solid red line represents the input tensor-to-scalar ratio parameter, r true = 0.00461.The dotted light teal line shows the normalized profile likelihood obtained from the HILC solution when the gain calibration for the CMB temperature is not included.
the profile likelihood, L p (r) [as defined in eqs.(2.14a) and (2.14b), respectively], and find that they are identical up to relative discrepancies of ≲ 10 −3 .
We show L(r) = L p (r) in figure 7 (teal solid line), together with a red vertical line corresponding to the input value, r true = 0.00461.The MLE is r = (4.30+0.56  −0.53 ) × 10 −3 .The bias, ∆r = −0.31× 10 −3 , is comparable to the uncertainty.We find that this bias is caused by the HWP polarization efficiency being lower than one.The B-mode signal is suppressed and r is underestimated.Note that the suppression due to the HWP polarization efficiency also affects the observed lensing amplitude: Âlens = 0.9548 +0.0093 −0.0096 .We also find non-detectable bias in the r true = 0 case: the best fit is r = 0 with 68% C.L. upper bound 0.00017, similarly to the ideal HWP case (see Section 3.1).

The weight of gain calibration
The inclusion of the gain calibration for the CMB temperature in the modeling of multi-frequency maps may seem inconsequential, but it has strong implications.We repeat the analysis of section 3.2, except that we now skip the gain calibration, i.e., we model the mi as in eq.(2.4) instead of eq.(2.6).The corresponding spherical harmonic coefficients read where the w/o subscript stresses that we are not calibrating the maps.By retracing the same steps as presented in section 2.2, we end up with an expression for the BB angular power spectrum of the HILC solution that reads ) where the w i ℓ,w/o are the HILC weights corresponding to the spherical harmonic coefficients of eq.(3.10).The corresponding normalized profile likelihood is shown in figure 7 (dotted light teal line).We now find a much lower MLE of the tensor-to-scalar ratio, r = (3.94+0. 52 −0.50 ) × 10 −3 , which is incompatible with r true , as the bias ∆r = −0.67 × 10 −3 is larger than the uncertainty.Similarly, the bias on the lensing amplitude is also stronger than the case when photometric calibration is included: Âlens = 0.913 ± 0.009.

Discussion
Clearly, gain calibration can partially mitigate the suppression of primordial B modes caused by the HWP.Of course, one can characterize the non-idealities in laboratory measurements and correct for them in the data.However, if HWPs are properly designed, gain calibration for the CMB temperature allows us to mitigate the effects of non-idealities on polarization in-flight for space missions.The ability to perform in-flight calibration is always valuable.
To this end, we derive some realistic recommendations that can help maximize its benefits.In section 4.2, we also discuss the assumptions underlying our end-to-end model and comment on the possibility of relaxing some of them.

HWP design recommendations
We express the relevant combinations of Mueller matrix elements in terms of a set of 7 independent values that uniquely determine the components of M hwp : the HWP Jones parameters, h 1,2 , β, ζ 1,2 and ξ 1,2 (see appendix B for their definitions).The loss parameters h 1,2 describe the deviation from the unitary transmission of E x,y ; β parametrizes the deviation from π of the phase shift between E x and E y ; ζ 1,2 and ξ 1,2 describe the amplitude and phase of the cross-polarization coupling.We write g(ν) ≡ m ii (ν), ρ(ν) ≡ [m qq (ν) − m uu (ν)]/2, and η(ν) ≡ [m qu (ν) + m uq (ν)]/2 as [50] where any dependence on ν is kept implicit for the sake of compactness.Designing a perfectly ideal HWP with identically vanishing Jones parameters is technically impossible.However, some parameters are easier to minimize than others.For example, ζ 1,2 (ν) ∼ 10 −2 can be achieved for both metal-mesh and multi-layer HWPs.If that is the case, the Taylor expansion of the above expressions for small ζ 1,2 (ν) yields, up to first order, We can further simplify these expressions by requiring h 1,2 ∼ 10 −2 , which implies ρ(ν) = g(ν) cos 2 [β(ν)/2] up to relative corrections of O(10 −4 ).Alternatively, by keeping h 1,2 free while requiring |h 1 − h 2 | to be small, we ensure that ρ(ν) = g(ν) cos 2 [β(ν)/2] still holds up to relative corrections of O(|h 1 − h 2 |).On the other hand, we cannot require β(ν) to be arbitrarily small due to the limitation of current technology.Keeping β(ν) free, we have If at least one of h 1 (ν) + h 2 (ν) and cos 2 [β(ν)/2] = [1 + cos β(ν)]/2 is slowly varying within the band, we find that ρ i CMB ≃ A i g i CMB , where A i is an appropriate factor that depends on β.Then, if we know A i with good precision, its effect can be undone by multiplying each multi-frequency polarization map by 1/A i .In this way, the gain calibration for the CMB temperature can partially mitigate the impact of the HWP polarization efficiency.
Regarding cross-polarization coupling, we argue that there are two strategies to keep its effects under control.First, we could simply require η(ν) ≲ 10 −3 so that the E → B leakage is negligible.However, this might be technically challenging.Another strategy is to exploit the fact that the HILC weights minimize the variance.Even if η(ν) is not vanishing small, as long as the η i CMB fluctuate around zero, the HILC should be able to mitigate their effect.HWP angle miscalibration An imperfect calibration of the HWP angle can dramatically affect the considerations we have presented so far.If an HWP with g i CMB ≃ ρ i CMB and ⟨η i CMB ⟩ ≃ 0, is rotated by some angle θ, its effective gain, polarization efficiency, and crosspolarization coupling are transformed as On the one hand, this causes the cross-polarization coupling coefficients to fluctuate around some non-zero value, making it impossible for the HILC to filter them out.On the other hand, the polarization efficiency and gain coefficients might strongly deviate from each other, reducing the benefits of gain calibration.Therefore, a good calibration of the HWP position angle, θ, is crucial to ensure the validity of our considerations and recommendations.Derotating the polarization maps by θ prior to the foreground cleaning step, as suggested in [70], would allow us to account for potential differences in the miscalibration angles of the HWPs.

Reviewing the underlying assumptions
We derived the model for multi-frequency maps and their spherical harmonics coefficients [eqs.(2.6) and (2.7), respectively] under several assumptions.We list them in order of appearance: Assumptions 1 and 2 cannot be relaxed while maintaining the semi-analytical treatment, since more complex beams and more refined map-makers can only be included in numerical simulations.On the other hand, assumptions 3 and 5 can be straightforwardly relaxed within our simple analytical model (given our focus on the HWP non-idealities, however, we chose not to play around with the bandpass shape or imperfect temperature gain calibration).
Assumption 4 can also be relaxed easily, but allowed us to analytically model the foreground cleaning step.Indeed, as soon as the SED of the foreground emission becomes anisotropic, the simple implementation of the HILC presented in section 2.2 is no longer able to recover the CMB signal accurately, and more elaborate methods such as Needlet ILC [91] and its moment [92] and Multiclustering [93] extensions will be needed.Although our quantitative results may be affected, qualitative conclusions will remain valid as long as the method is still based on ILC.
It would be interesting to relax some of these assumptions and check whether the recommendations presented in section 4.1 still ensure that gain calibration for the CMB temperature can mitigate polarization systematics due to the HWP non-idealities.We leave this analysis for future work.

Conclusions and perspectives
In this work, we presented a simple framework to propagate the HWP non-idealities through the three macro-steps of any CMB experiment: observation of multi-frequency maps, foreground cleaning, and power spectra estimation.We focused on the impact of non-idealities on the tensor-to-scalar ratio parameter, r.
We generalized the formalism presented in [49] to include the polarized Galactic foreground emission (dust and synchrotron), foreground cleaning using a blind method (HILC), bandpass integration, noise, beam smoothing, and gain calibration for the CMB temperature.As a concrete working case, we considered a full-sky CMB mission with LiteBIRD-like specifics [16].
We validated the code against an ideal HWP and confirmed that the MLE r had negligible bias.Then, we employed more realistic Mueller matrix elements for each of the three telescopes of LiteBIRD and found r = (4.30+0.56  −0.53 ) × 10 −3 .We showed how the suppression is mostly due to the effective polarization efficiency of the HWP, which averages to a value lower than 1.The effective cross-polarization coupling and the foreground residual are found to be subdominant in our output B-mode power spectrum.
We found that the bias in r significantly worsens if gain calibration for the CMB temperature is not included in the modeled multi-frequency maps: r = (3.94+0.52  −0.50 ) × 10 −3 , which is incompatible with the input value.Gain calibration would perfectly remove the HWP effects if ρ i CMB = g i CMB and η i CMB = 0, which are, however, unrealistic requirements.Still, we showed that an effective mitigation can be achieved if we can factorize ρ i CMB ≃ A i g i CMB , we have good knowledge of the A i coefficients, and ⟨η i CMB ⟩ ≃ 0. These considerations helped us to formulate some recommendations on the HWP design in terms of the HWP Jones parameters: ▷ Cross-polarization coupling should be small, ζ 1,2 ≲ 10 −2 , which can be achieved for both metal-mesh and multi-layer HWPs; ▷ The loss parameters should also be small, h 1,2 ≲ 10 −2 , or, alternatively, |h 1 − h 2 | ≲ 10 −3 ; ▷ At least one of h 1 (ν) + h 2 (ν) and [1 + cos β(ν)]/2 should be slowly varying within the band, so that ρ i CMB ≃ A i g i CMB ; ▷ Cross-polarization coupling can be kept under control by requiring ζ 1,2 to be even smaller, or alternatively, by ensuring that η i CMB fluctuates around zero.
One can characterize the non-idealities of the HWP in laboratory measurements, and a requirement for the smallness of a bias in r gives a requirement for the accuracy of the calibration in the laboratory.However, if the above recommendations are implemented in the design of the HWP used for space missions, the in-flight gain calibration for the CMB temperature can also be used to check and correct for the effects of HWP non-idealities in the data, complementing the laboratory calibration.
Some of the recommendations above depend strongly on the class of foreground cleaning methods we used in our end-to-end model.We used a blind method (HILC), but if one were to use a parametric component separation method to derive design recommendations, they would likely be different from those listed above.This highlights the importance of developing analysis strategies together with hardware designs.This work represents a first generalization of the model presented in [49] towards a more realistic account of how the HWP non-idealities affect the observed CMB.However, being semi-analytical, this framework still relies on several simplifying assumptions (see section 4.2).One of the most crucial is the isotropy of the foreground SED.It would be interesting to relax this assumption and repeat the analysis carried out in this paper, using more elaborate ILCbased methods (e.g., [92,93]).This would help us test the robustness of our recommendations for the design of HWPs in a more realistic context.We leave this study for future work.
where h 1,2 are loss parameters describing the deviation from the unitary transmission of E x,y ; β parametrizes the deviation from π of the phase shift between E x and E y ; ζ 1,2 and ξ 1,2 describe the amplitude and phase of the cross-polarization coupling.All Jones parameters tend to zero in the ideal limit.

Figure 1 .
Figure1.For an ideal HWP, the rescaled angular power spectrum, D BB ℓ , of the HILC solution (dashed teal line) overlaps the input CMB spectrum (black solid line) for a wide range of multipoles.For large ℓ, the two spectra begin to diverge as we approach the instrumental resolution.This can be seen by looking at the dotted gray line, representing the residual noise, which intersects the input spectrum at ℓ ∼ 325.For completeness, we also plot the input dust and synchrotron D BB

Figure 2 .
Figure 2. HILC weights, w iℓ , for each of the three telescopes with an ideal HWP.In each plot, different colored lines correspond to different frequency channels: from purple to yellow for lower to higher central frequencies (see[16] for details on the channels' specifics).The corresponding BB angular power spectrum is shown in figure 1 (dashed teal line).

Figure 3 .
Figure 3. HWP Mueller matrix elements for LFT (purple), MFT (red) and HFT (orange) as function of frequency.For LFT, we consider a symmetric stack design[88] provided with ARC[89], compute its Mueller matrix elements, and rotate them of 55.02 ○ , to express them in a reference frame with the x axis parallel to the HWP optic axis.Instead, the Mueller matrix elements for MFT and HFT are obtained by following the same procedure and input simulations as done in[50], and rotating them of 0.29 ○ .The dashed gray lines represent the ideal values of each element.

Figure 4 .
Figure 4. Same as figure1but for the realistic HWP discussed in section 3.2 (dashed teal line).Compared to the ideal HWP case shown in figure1, the non-ideal HILC solution slightly differs from the input CMB at low multipoles.For comparison, we also show the residual noise bias (dotted gray line) and the foreground residual (red dotted line).They both show more features than their counterparts in figure1.The w i ℓ weights corresponding to the HILC solution are shown in figure5.

Figure 5 .
Figure 5. Same as figure 2 but for the Mueller matrix elements given in figure 3. The corresponding BB angular power spectrum is shown in figure 4 (dashed teal line).

2 synchrotronFigure 6 .
Figure 6.Different contributions to the B-mode power spectrum of the HILC solution (teal solid line).We focus on a different component (CMB, dust, and synchrotron) in each of the panels.The effective polarization efficiency and cross-polarization coupling components are shown in dashed and dotted, respectively.The largest contribution comes from the polarization efficiency component of the CMB.

Figure 7
Figure 7. Normalized profile likelihood, L(r) = L p (r), obtained from the HILC solution, ĈBB ℓ,hilc , given the HWP specifics presented in section 3.2 (teal solid line).The likelihood has a maximum at r = 0.00430.The shaded region identifies the 68% CL interval, and goes from r −0.00053 to r +0.00056 .The solid red line represents the input tensor-to-scalar ratio parameter, r true = 0.00461.The dotted light teal line shows the normalized profile likelihood obtained from the HILC solution when the gain calibration for the CMB temperature is not included.

Table 1 .
Left panel: SED parameters entering in eqs.(3.1) for each component as reported in