Impact of half-wave plate systematics on the measurement of cosmic birefringence from CMB polarization

Polarization of the cosmic microwave background (CMB) can probe new parity-violating physics such as cosmic birefringence (CB), which requires exquisite control over instrumental systematics. The non-idealities of the half-wave plate (HWP) represent a source of systematics when used as a polarization modulator. We study their impact on the CMB angular power spectra, which is partially degenerate with CB and miscalibration of the polarization angle. We use full-sky beam convolution simulations including HWP to generate mock noiseless time-ordered data, process them through a bin averaging map-maker, and calculate the power spectra including $TB$ and $EB$ correlations. We also derive analytical formulae which accurately model the observed spectra. For our choice of HWP parameters, the HWP-induced angle amounts to a few degrees, which could be misinterpreted as CB. Accurate knowledge of the HWP is required to mitigate this. Our simulation and analytical formulae will be useful for deriving requirements for the accuracy of HWP calibration.

The CMB polarization can also probe new parity-violating physics [27]. For example, in the presence of a time-dependent parity-violating pseudoscalar field, the linear polarization plane of CMB photons would rotate while they travel towards us [28][29][30]. Because of its similarity with photon propagation through a birefringent material, this phenomenon is referred to as cosmic birefringence (CB). The so-called CB angle, β, denotes the overall rotation angle from last scattering to present times. Although the effect of β on the observed CMB angular power spectra is degenerate with an instrumental miscalibration of the polarization angle [31][32][33][34], the methodology proposed in [35][36][37], which relies on the polarized Galactic foreground emission to determine miscalibration angles, allowed to infer β = 0.35 ± 0.14 • at 68% C.L. [38] from nearly full-sky Planck polarization data [39]. Subsequent works [40][41][42] reported more precise measurements for β. The statistical significance of β is expected to improve with the next generation of CMB experiments, given the high precision at which they aim to calibrate the absolute position angle of linear polarization. This will make it unnecessary to rely on the Galactic foreground to calibrate angles and measure β [27], hence avoiding the potential complications highlighted in [43].
The unprecedented sensitivity goals of future surveys, aiming to detect faint primordial B modes, can only be achieved if systematics are kept under control. To this end, a promising strategy is to employ a rotating half-wave plate (HWP) as a polarization modulator. As shown by the previous analyses [44][45][46][47][48][49][50][51], a rotating HWP can both mitigate the 1/f noise component [44] and reduce a potential temperature-to-polarization (I → P ) leakage due to the pair differencing of orthogonal detectors [52,53]. Because of these advantages, HWPs are used in the design of some next-generation experiments, including LiteBIRD [16]. However, nonidealities in realistic HWPs induce additional systematics which should be well understood in order for future experiments to meet their sensitivity requirements. This necessity motivated a number of recent works, from descriptions of HWP non-idealities [46,[54][55][56] and their impact on measured angular power spectra [57] to mitigation strategies [58][59][60][61].
In this paper we study how HWP non-idealities can affect the estimated CMB angular power spectra if overlooked in the map-making step. We employ a modified version of the publicly available beam convolution code beamconv 1 [62,63] and simulate two sets of noiseless time-ordered data (TOD). The two simulations make different assumptions on the HWP behavior. In the first case the HWP is assumed to be ideal, while non-idealities are included in the second case. We then process the two TOD sets with a map-maker assuming the ideal HWP and compare the output power spectra. We also derive a set of analytic expressions for the estimated angular power spectra as functions of the input spectra and the elements of the HWP Mueller matrix. These formulae accurately model the output power spectra. Finally, we show that neglecting the non-idealities in the map-maker affects the observed spectra in a way that is partially degenerate with the CB and instrumental miscalibration of the polarization angle. This effect is evident in the simulations and the analytical formulae.
The rest of this paper is organized as follows. In section 2 we present a simple data model for the signal measured by a single detector; generalize it to a larger focal plane and a longer observation time; and introduce the bin averaging map-making method employed in the paper to convert the TOD to maps. In section 3 we discuss the instrument specifics we have implemented in the simulation and show the output angular power spectra. The interpretation of the result is the topic of section 4, where we derive some analytical formulae modeling it with good precision. In section 5 we show how the effect of the HWP non-idealities is partially degenerate with an instrumental miscalibration of the polarization angle, and can therefore be misinterpreted as CB. We quantify the HWP-induced miscalibration angle, which amounts to a few degrees for our choice of the HWP parameters. Conclusions and outlook are presented in section 6.

Mathematical framework: data model and map-maker
Data model for a single detector Polarized radiation can be described by the Stokes I, Q, U and V parameters or, more compactly, by a Stokes vector, S ≡ (I, Q, U, V ). In this paper we use the "CMB convention" for the sign of Stokes U [64] and define the Stokes parameters in right-handed coordinates with the z axis taken in the direction of the observer's line of sight (telescope boresight). The Stokes vector is transformed as S → S = R ϕ S by rotating the coordinates by an angle ϕ, where (2.1) Defining the position angle of the plane of linear polarization, θ, by Q ± iU = P e ±2iθ with P = Q 2 + U 2 and 2θ = arctan(U/Q), the rotation of coordinates shifts the position angle as θ → θ = θ − ϕ.
The action of any polarization-altering device on S can be encoded in a Mueller matrix M, so that the outgoing Stokes vector reads S = MS [65]. In our case of interest, S represents the incoming CMB radiation and M the Mueller matrix of a telescope that employs a rotating HWP as a polarization modulator, i.e.
where R ϕ is given in eq.  x y x sky y sky ψ sky coords.
x y x hwp y hwp φ HWP coords.
x y x det y det ξ detector coords. We can then model the signal d measured by one detector as where n represents an additional noise term.
Modeling the TOD In a realistic CMB experiment, n det detectors collect data by scanning the sky for an extended period of time, resulting in n obs observations for each detector. All together, these n det × n obs measurements constitute the TOD. We represent the TOD as a vector d given by where m denotes the {I, Q, U, V } pixelized sky maps, A the response matrix, and n the noise component. Eq. (2.5) generalizes eq. (2.4) to larger n obs and n det .
Bin averaging map-maker To extract physical information from the TOD, we convert them to the map domain via some map-making procedure. A simple method often employed in the CMB analysis is the bin averaging [66], that estimates the sky map as where A is the response matrix assumed by the map-maker. As long as the beam is axisymmetric and purely co-polarized, and the correlated component of the noise, such as 1/f , is negligible, the bin averaging can, in principle, recover the input {I, Q, U, V } maps. Whether the reconstructed maps actually reproduce the sky signal or not depends on how well the instrument specifics are encoded in the map-maker or, in other words, how close A is to A. When A = A and n is uncorrelated in time, m is the optimal (unbiased and minimum-variance) estimator of m.

Simulation setup and output
We generate statistically isotropic random Gaussian {I, Q, U } CMB maps with HEALPix 2 [67] resolution of n side = 512 (high enough to avoid aliasing effects) by feeding the best-fit 2018 one-year simulation Figure 2. Simulated boresight hit maps for one-month (left) and one-year (right) observations. The two panels share the same logarithmic color map, although the range shown is truncated for the one-month case. Unobserved pixels appear gray.
Planck power spectra [3] to the synfast function of healpy 3 (the Python implementation of HEALPix). We choose to neglect V here, since the circularly polarized component of CMB is expected to be negligible 4 . The observation of the input maps is simulated by a modified version of the publicly available library beamconv. This choice is motivated by beamconv's ability to simulate TOD with realistic HWPs, scanning strategies and beams, which makes it a promising framework to develop simulations for, among others, LiteBIRD-like experiments. The changes we have implemented to the library all aim to better tailor the simulations to LiteBIRD-like specifics. In particular: Scanning strategy We implement a LiteBIRD-like scanning strategy by mimicking the relevant functionalities of pyScan 5 . The values of the telescope boresight and precession angles, together with their rotation parameters, are specified in table 1. We simulate one year of observations to cover the full sky (see figure 2).
Instrument We work with 160 detectors from the 140 GHz channel of LiteBIRD's Medium Frequency Telescope (MFT) and read the relevant parameters from [16]: the HWP rotation rate, the full-width-at-half-maximum (FWHM) of the (Gaussian and co-polarized) beam, the instrument sampling frequency and the detectors' pointing offsets. See table  1 for their numerical values. HWP Mueller matrix The Mueller matrix elements for the MFT's HWP at 140 GHz are taken from [57], up to a coordinate change from International Astronomical Union (IAU) to CMB standards that flips the sign of the m iu , m qu , m ui and m uq elements: This is the HWP Mueller matrix we assume when including non-idealities 6 . Since the elements of M hwp are frequency-dependent, choosing a different frequency would result in slightly different output spectra.
We run two simulations for one-year observations. Noise is not included in either simulation to isolate the effect of HWP non-idealities in the signal; thus, using a different n det is almost free from consequences and our results do not change using fewer detectors. In the first simulation we assume the ideal HWP by setting M hwp = M ideal ≡ diag(1, 1, −1), while we account for non-idealities in the second one. We convert both TODs to {I, Q, U } maps by the bin averaging map-maker (see eq. (2.6)) whose response matrix A assumes the ideal HWP described by M ideal . We then calculate two sets of full-sky angular power spectra using the anafast function of healpy. We denote the first (second) set of output spectra ,hwp spectra are plotted in figure 3, together with the input spectra multiplied by the Gaussian beam transfer functions, D XY ,in . The simple map-maker recovers the input spectra with average deviations less than 0.1% in the plotted range when processing the TOD generated with M ideal , while important discrepancies arise for the non-ideal case.
We do not account for photometric calibration, although it represents a crucial step in any CMB analysis pipeline. Gain calibration, if perfect, would ensure intensity to be recovered exactly, hence compensating the lack of power in D T T visible in figure 3. The discrepancies in D T E and D T B would also be reduced, although not removed. The discussion and results presented in the following would however not change, reason why we omit the step.

Analytical estimate of the output spectra
To understand the simulation results, we derive approximate analytical formulae for the angular power spectra affected by HWP non-idealities. Since we are neglecting any circularly polarized component, the Stokes vector is given by S = (I, Q, U ). To obtain analytical formulae we apply the bin averaging map-maker of eq. (2.6) to a minimal TOD consisting to eq. (2.4), we obtain where m ss (s, s = i,q,u) are the elements of non-ideal M hwp and α denotes the sum of the HWP's (φ) and the telescope's (ψ) angles 8 : α ≡ φ + ψ. The quantities with the subscript "in" on the right hand side denote the sky signals, whereas S = ( I, Q, U ) on the left hand side are maps recovered by the map-maker. These formulae are applicable to our case as long as eq. (2.4) accurately models the TOD simulated by beamconv, which is the case for an axisymmetric and purely co-polarized beam. Eqs. (4.1) model S p reconstructed from four observations of the pixel p, one for each detector. If each of those 4 detectors were to observe that same pixel n p times, the change in eqs. (4.1) would amount to substituting for n = {2, 4}. If p is observed with a uniform enough sample of α t values and n p is large enough, these terms can be neglected, resulting in We expect this to be a good approximation, given the presence of a rapidly spinning HWP and the good coverage of the simulation (see figure 2). By expanding eq. (4.3) in spherical harmonics, we write the corresponding angular power spectra as a mixing of the input ones weighted by combinations of the non-ideal HWP's Mueller matrix elements: These analytical formulae explain quite well the non-ideal output spectra C XY ,hwp (see figure 4). They are especially accurate on large scales, 500, where average deviations between C XY ,hwp and C are less than 0.1%. Larger deviations on smaller scales are due to the approximate nature of eq. (4.3). Cosine and sine terms do not average out exactly, resulting in pixel-bypixel fluctuations on smaller scales. 8 For the simple 4-detector configuration we are considering, the response matrix can be expressed as A = BR ξ−φ MhwpR φ+ψ , where B accounts for the different ξ angles of the four detectors and happens to satisfy B T B = diag(1, 1/2, 1/2). As for the map-maker response matrix, A = BR ξ−φ M ideal R φ+ψ . All BR ξ−φ terms cancel out in eq. (2.6) and we are left with S = R T φ+ψ M ideal MhwpR φ+ψ Sin. The discrepancies between S and Sin can therefore only depend on φ + ψ.

Impact on cosmic birefringence
Next generation CMB experiments are expected to measure the CMB polarization with unprecedented sensitivity and improve the constraints on the CB angle, β, recently obtained from the Planck data [38,[40][41][42]. Here we discuss how HWP non-idealities can impact such constraints in the particular case of a LiteBIRD-like mission discussed so far.
First, we recall that the sign of β reported in the literature is also chosen to follow the CMB convention and a positive β corresponds to a clockwise rotation on the sky [27]. The isotropic CB angle, β, and a miscalibration of the instrument polarization angle, ∆α, affect the observed spectra identically, since both rotate the observed Stokes parameters in the same way. The observed spectra then satisfy the equations [77,78] C EB ,obs = where θ represents rotation in the position angle of the plane of linear polarization including β, ∆α, or their sum. Not accounting for the HWP non-idealities in the map-maker step is degenerate with θ, as it is evident from both our simulations and the analytical formulae given in eq. (4.4). We will refer to this additional rotation of the polarization plane as the HWP-induced miscalibration.
HWP-induced miscalibration from the simulated output spectra We separately fit the simulated C EB ,hwp and C T B ,hwp for the angles θ EB and θ T B , respectively, using the leastsquares method with variance given by for XY = {EB, T B}, respectively. The best-fit values, θ EB = 3.800 • ± 0.007 • and θ T B = 3.79 • ± 0.02 • , are compatible with each other in agreement with eqs. (5.1). The observed and best-fit spectra are plotted in figure 5 and are in good agreement.
HWP-induced miscalibration from the analytical formulae Using the fact that both C EB in and C T B in simply fluctuate around zero, eqs. (4.4) can be rearranged to express C EB and C T B similarly to the C XY ,obs of eqs. (5.1): in agreement with the best-fit values reported above.
If we were to relax all the underlying assumptions at once, we could not write θ this compactly. However, controlled generalizations do not necessarily spoil the simplicity of the analytical formulae. For instance, accounting for the frequency dependence of both the HWP Mueller matrix elements and the CMB signal, θ can be expressed as (see appendix A for the derivation): where S CMB (ν) denotes the CMB spectral energy distribution (SED).
Another assumption that can be relaxed without spoiling the simplicity of the analytical formulae is the absence of miscalibration angles in the map-maker. When the telescope, HWP, and detector angles are not exactly known, ψ = ψ + δψ, φ = φ + δφ, and ξ = ξ + δξ, where the hat denotes the values assumed by the map-maker. As long as we neglect the frequency dependence of δψ, δφ and δξ, we find (see appendix B for the derivation) The sign difference between the contributions from δξ and δψ + 2δφ is due to the presence of the HWP. Ideally, the HWP acts on a polarization vector by reflecting it over its fast axis. This causes counterclockwise rotations applied before the HWP to look clockwise after, meaning that δφ + δψ should be subtracted from δξ − δφ (see eq. (2.2)).

Conclusions and outlook
In this work, we studied how overlooking HWP non-idealities during map-making can affect the reconstructed angular power spectra of CMB temperature and polarization fields. We focused on the impact of non-idealities on the measurement of the CB angle, β.
As a concrete working case, we considered a single frequency channel (140 GHz) of a space CMB mission with LiteBIRD-like specifics: scanning strategy, sampling frequency, detectors' pointing offsets and their polarization sensitivity directions, FWHM of the Gaussian beam and HWP specifics (rotation frequency and Mueller matrix elements). We employed the publicly available beam-convolution code beamconv to simulate the noiseless TOD for the above instrument and scanning specifications. We ran two different simulations: the HWP has been assumed to be ideal in the first simulation, while a realistic Mueller matrix has been employed in the second. We then converted both TODs to maps by a bin averaging mapmaker that neglects the HWP non-idealities. As expected, the output spectra computed from the ideal simulation faithfully recovered the input spectra, while the spectra of the non-ideal maps showed a very different behavior ( figure 3). We also derived a set of analytical formulae (see eq. (4.4)) that accurately model the reconstructed angular power spectra as functions of the input spectra and the HWP Mueller matrix elements.
We studied the impact of the HWP non-idealities on β. We found that neglecting them in the map-making step induces an additional miscalibration of the polarization angle which might be erroneously interpreted as CB. For the concrete case we studied, the miscalibration angle induced by the HWP non-idealities amounts to θ 3.8 • . This value, obtained by fitting the output angular power spectra from the simulation, is compatible with the prediction from the analytical formulae (see eq. (5.4)).
Definitive confirmation of the current hint of CB [38,[40][41][42] requires the systematic uncertainty in the absolute position angle of linear polarization to be well below 0.1 • [27].
We must therefore acquire accurate knowledge of the Mueller matrix elements via calibration, so that the systematic uncertainty in θ due to HWP non-idealities is well below 0.1 • . With such knowledge, one can take into account HWP non-idealities either during the map-making step or when interpreting the angular power spectra. As one cannot know the Mueller matrix elements perfectly, any remaining mismatch between the true Mueller matrix and the matrix assumed by the map-maker still affects the power spectra. Our simulation and analytical formulae will be useful for deriving the required accuracy of HWP calibration to meet specific science goals.
The situation we considered in this paper is still simplistic: we simulated a single frequency channel in the absence of noise, and we used a Gaussian beam and a simple bin averaging map-maker. However, a similar analysis can be carried out for more complex cases. It is of utmost importance to make better predictions about how HWP non-idealities realistically affect the data collected by CMB experiments and, therefore, the cosmological information extracted from them. In this direction, we plan to carry on the following steps: i) drop the single frequency approximation, generalizing the results discussed here to a finite frequency bandwidth; ii) add a noise component to the TOD; iii) study the combined effect of beam asymmetries and HWP non-idealities; iv) include non-idealities in the map-maker and study how the uncertainties in our knowledge of non-idealities might propagate to the observed angular power spectra; and v) derive requirements for the accuracy of HWP calibration. We leave these topics for future work.

B Additional miscalibration angles
So far, we neglected any miscalibration angles in the map-maker, i.e. we assumed the response matrix A to encode the true values of the telescope, HWP, and detector angles: ψ ≡ ψ, φ ≡ φ, and ξ ≡ ξ, where the hat denotes the values assumed by the map-maker. We now consider a more general case by allowing for deviations: ψ = ψ + δψ, φ = φ + δφ, and ξ = ξ + δξ.
Single frequency Repeating the analysis presented in section 4 with miscalibration angles, eq. Therefore, the additional miscalibration angles simply shift θ, as expected.
Finite frequency bandwidth Taking into account a finite frequency bandwidth and miscalibration angles simultaneously is slightly more complicated, but does not spoil the analytic treatment as long as δθ is assumed to be frequency-independent. The generalization of eq.