Modeling optical systematics for the Taurus CMB experiment

We simulate a variety of optical systematics for Taurus, a balloon-borne cosmic microwave background (CMB) polarisation experiment, to assess their impact on large-scale E-mode polarisation measurements and constraints of the optical depth to reionisation {\tau}. We model a one-month flight of Taurus from Wanaka, New Zealand aboard a super-pressure balloon (SPB). We simulate night-time scans of both the CMB and dust foregrounds in the 150GHz band, one of Taurus's four observing bands. We consider a variety of possible systematics that may affect Taurus's observations, including non-gaussian beams, pointing reconstruction error, and half-wave plate (HWP) non-idealities. For each of these, we evaluate the residual power in the difference between maps simulated with and without the systematic, and compare this to the expected signal level corresponding to Taurus's science goals. Our results indicate that most of the HWP-related systematics can be mitigated to be smaller than sample variance by calibrating with Planck's TT spectrum and using an achromatic HWP model, with a preference for five layers of sapphire to ensure good systematic control. However, additional beam characterization will be required to mitigate far-sidelobe pickup from dust on larger scales.


Introduction
Since the mid-1990s a concordance model of cosmology has emerged, known as "Lambda cold dark matter" (ΛCDM).The most basic ΛCDM model, defined with only six parameters, has been found to fit a wide range of observations from large-scale structure, supernovae and the cosmic microwave background (CMB) rather well [1].Each of these parameters tells us specific information about our Universe, and can also help constrain fundamental physics.Some of the strictest constraints on them come from observations of the CMB [2].
One of the six parameters is τ , the optical depth to reionisation.At face value, τ tells us the integrated density of free electrons along the line of sight since the end of reionisation, or equivalently the redshift z re at which an instantaneous reionisation would have taken place.During and after reionisation, CMB photons can scatter on free electrons.In the non-relativistic limit, this reduces to Thomson scattering that suppresses temperature anisotropies in the CMB by a factor of e −2τ at the power spectrum level.However, the differential Thomson scattering cross-section is polarisation dependent.An electron seeing a local quadrupolar anisotropy of source photons will re-scatter them with a net linear polarisation [3][4][5][6].Since our lines of sight as observers point to the last scatterings of the individual photons, the net linear polarisation from any particular line of sight will then depend on the CMB quadrupole moment as it appeared on the cosmological horizon at the location of that last scattering [7].The horizon size at reionisation subtends a large area on the sky, therefore we expect reionisation to induce polarisation correlations on large angles, which will lead an observer today to see a small excess in the polarisation power spectrum at small multipoles.
Knowledge of τ also informs attempts to estimate the sum of neutrino masses from the CMB.Primordial density perturbations seed structure growth, which is damped by relativistic neutrinos streaming out of areas of high gravitational potential [8].Since massive neutrinos become non-relativistic at late times, precision measurements of this damping as a function of scale (e.g., from gravitational lensing of the CMB) can be used to constrain the sum of neutrino masses, Σm ν [9].This requires a precise estimate of A s , the amplitude of the primordial density perturbations, but there is a degeneracy between A s and τ , as the amplitude of the CMB temperature power spectrum is a multiple of A s e −2τ [10].A precise estimate of Σm ν from the CMB thus relies on an estimate of τ .Today, the best constraints on A s and τ come from the Planck satellite [2,11,12].While A s is determined to percent level accuracy, Planck 's constraints on τ are an order of magnitude looser, with σ(τ ) = 0.006.This motivates an experiment dedicated to constraining τ independently from Planck.Furthermore, in the coming decade, new observations will change our understanding of galaxies and stars in the epoch of reionisation [13].This new history of reionisation can be compared with the CMB inferred τ in a valuable cross-check of ΛCDM.
In Figure 1 we can see how the amplitude of the EE spectrum at large angular scales depends on the value of τ even if A s e −2τ stays constant.Planck 's 1σ uncertainty on τ is larger than the EE spectrum's sample variance for 5 < ℓ < 15.Then an experiment that can measure the EE power spectrum around ℓ ∼ 10 could constrain τ .Beyond ℓ = 20, any variation in D EE ℓ due to the uncertainty on τ is compatible with cosmic variance.Depending on the observed sky fraction, the signal starts being dominated by sample variance for 15 ≤ ℓ ≤ 20.Therefore, to probe the angular power spectra on large scales, an experiment with a wide sky coverage is needed.Furthermore, to disentangle polarized foregrounds from the CMB, multi-frequency maps are required.Studies of our Galaxy will also benefit from foreground polarisation maps at higher-than-Planck sensitivity.Taurus, a balloon-borne CMB experiment flying at mid-latitudes, can assure both the wide spatial and spectral coverage.Taurus' current configuration is described in detail in [14], and its cryostat design in [15].To ensure the success of Taurus, the mission must be designed to minimize systematics.Modeling can inform us, however incompletely, on the amount of systematic error that the current Taurus design will face on deployment, and therefore influence the evolution of that design.In this paper, we simulate time-ordered data and sky maps for a Taurus mission profile, evaluating the impact that beam and HWP systematics might have.In Section 2, we introduce the Taurus mission.Section 3 describes how we simulate a Taurus flight.In Section 4 we describe the treatment of the simulated sky maps and introduce the various calibration parameters we apply.Results are presented in Section 5, and we discuss their meaning for Taurus in Section 6, outlining our future lines of inquiry.

Taurus: Experiment design
The baseline Taurus design has three cryogenically cooled receivers that will observe the sky at four frequencies using arrays of superconducting transition-edge sensor (TES) detectors coupled to feedhorns.All the feedhorn pixels are dichroic, meaning they feed two TESs sensitive to two distinct frequencies.This increases the amount of information that can be gathered from a given focal plane area.The four bands are grouped in pairs, with the bands  [2] (color bar), with the product A s e −2τ kept constant.Right: the relative variation from the Planck nominal value for 2 ≤ ℓ ≤ 24.The grey lines represent the sample variance associated with the EE spectrum for different observed sky fractions: the dotted line is for f sky = 0.4, the dashed line for f sky = 0.7 and the full line is the cosmic variance limit on σ(τ ).In our simulations, we estimate the power spectrum on roughly 40 % of the sky.centred at 150 and 220 GHz designated LF1 and LF2 respectively, while the bands centred at 280 and 350 GHz are HF1 and HF2.A batch of detectors, split between two receivers, observes the LF bands.Another batch, grouped in one focal plane, observes the HF bands.Each of the three receivers is a refractor with a 28 • field of view.There will be a polarisation modulator in front of each telescope's stop: due to the dichroic nature of the receivers, a broadband half-wave plate (HWP) is required.The HWP is stepped regularly to improve cross-linking of the observations and mitigate beam systematics.The three receivers are offset by 90 • in azimuthal pointing from one another.
Taurus will fly suspended from a super-pressure balloon (SPB) launched from New Zealand.We are planning for the flight duration to be of order one month, during which the balloon circumnavigates at mid-latitudes pushed by the prevailing stratospheric winds (as represented in Figure 2).During the day, Taurus will recharge its batteries by pointing its solar panels towards the Sun.Thanks to the long flight duration, limiting Taurus to night-time observations only modestly reduces the accessible sky coverage.The instrument rotates about its axis and observes at a fixed boresight elevation of 35 • ; because Taurus scans at night, it does not need to avoid the Sun and can rotate through the whole 360 • azimuthal range.The payload rotates at 30 °s−1 .The HWPs are stepped once per day by integer multiples of 22.5 • .With such a scan strategy, a one-month flight will get well-conditioned maps as shown in Figure 3: in our simulations, the condition number is smaller than 2.3 for 90 % of the observed pixels and smaller than 3 for 99 % of them.

Modeling the flight
There have been five SPB launches from Wanaka in New Zealand since the Columbia Scientific Ballooning Facility began its operations there in 2015.They are designed to minimize altitude fluctuations due to day/night temperature cycles.The SPBs had long hold times for three

COSI SuperBIT
Figure 2. Trajectories for the three longest NASA super-pressure balloon flights to date from Wanaka, NZ [16][17][18].Alternating shading along each track demarcates each day of flight.Our simulated trajectory follows the blue track (qualification) for thirty days.
out of five flights: the qualification flight lasted 32 days [16], the SuperBIT flight was 39 days long [17], and COSI stayed aloft for 46 days [18].Their trajectories are represented in Figure 2.For the purposes of this work, we have chosen to assume the trajectory of Flight 662 NT, the qualification flight, which had an average altitude of 33 500 m.Our simulation covers the first 30 days of that flight, from March 26 th to April 26 th .We let the simulated payload's altitude vary daily around an average value of 32 700 m, with a standard deviation of 200 m.The resulting sky coverage is visible in the left panel of Figure 3.At 45 • S, roughly 85 % of the sky is visible over a 24-hour period.Not all of this is visible for a balloon-borne telescope, as the atmosphere degrades performance at low elevation, and the balloon obscures the zenith.With Taurus' pointing and field of view, the maximum fraction of the sky visible is f sky ≈ 80 %.Depending on the season of the launch, parts of the sky will be only visible during the day, reducing the observed sky fraction.We follow the civil definitions of sunset/sunrise (center of the solar disk 5 • below the horizon), which we compute for a given position and altitude along the trajectory using PyEphem [19].A flight starting March 1 st will see f sky = 62 % while a June launch would get f sky = 70 %.The input CMB maps are generated from the Planck power spectrum with healpy1 [20,21], while the dust emission is based on the PySM [22] d9 model.In that model, the dust emission is given by a modified blackbody whose temperature T d and spectral index β vary from pixel to pixel.While more complex models exist in PySM, it is a good starting point to assess the impact of differences between the CMB and the dust on our instrument's performance.

Focal plane modeling
We model beams for 49 detector pairs of one of the LF receivers, covering a 20  1.In this study, we focus on simulations of the 150 GHz band.
For added realism, we use the TICRA-tools GRASP package2 to simulate Physical Optics (PO) beam models for each of the detector positions we study.The telescope model in GRASP comprises full-wave simulations of detector horns, two HDPE lenses, and a circular aperture representing the vacuum window.The aperture is 228 mm in diameter with an fnumber of f /1.6, and the lenses focus the light down to a 200 mm × 200 mm square focal plane.The fast optics result in a −10 dB edge taper averaged across the field of view.We neglect reflections on the walls of the optics tube.A Lambertian source at the aperture emits −20 dB of unpolarized power relative to the power passing through the aperture due to the input feedhorn beam.This is done to simulate scattering by the vacuum window or other optical elements near the aperture.The aperture also acts as the overall stop of the system.GRASP outputs both co-and cross-polar far-field beams, which can then be translated to Stokes (I, Q, U ) beams.We sample the beam finely within 3 • of the beam centroid (the main beam), and more loosely up to 30 • (to model beam sidelobes).The coarser resolution on sidelobes is reasonable as they are both much weaker and extend much further than the main beam.
Finally, those PO beams can also be approximated by a fit model of their co-polar component.We adopt the fitting algorithm described in [23] and [24] to create a symmetrized beam model of the central pixel's PO beam.An azimuthally averaged beam profile for the different beam models can be seen in Figure 4.In blue, a 33 ′ FWHM Gaussian beam.The green and red lines are the azimuthally averaged profiles for the PO simulations.In orange, an analytical fitting of the PO beam for the central detector.To evaluate the relative contribution of main lobe and side lobe for the PO beams, we can run simulations that ignore any beam power beyond 3 • of the beam centroid.This is represented by the gray outline.

Half-wave plate
A half-wave plate modulates polarisation by introducing a phase of π radians between the two linear polarisations.We often describe the incoming radiation using the transposed Stokes vector S t = (I, Q, U, V ), where I is total intensity, Q, U are the two linear polarisations, and V is circular polarisation.These four Stokes parameters are real-valued, with The action of the HWP on the Stokes vector is expressed by means of a Mueller matrix M HW P , a 4 × 4 square matrix that transforms S [25,26].For an ideal HWP, the Mueller matrix is diagonal, and When the ideal HWP rotates by an angle θ, the Stokes vector then gets modulated: A real HWP will have a more complex Mueller matrix: the diagonal elements will no longer be ±1, and the other elements can be non-zero [27].The off-diagonal elements will cause temperature to polarisation leakage, or complicate the mixing between the Q and U terms in Equation 3.2 [28,29].Furthermore, the elements of the Mueller matrix will be frequency dependent, so those non-idealities will be different for every band we observe, and can even be tricky to model within a band.It is therefore in our interest to design a HWP which has a stable behavior over a wide frequency band: an achromatic HWP (AHWP) [30].Each of the Taurus HWPs is designed to be achromatic, in order to cover the two LF or HF bands.AHWPs can be created either by stacking layers of birefringent materials with their ordinary axes at specific angles [30,31], or by using a wire mesh acting as a metamaterial [32].For Taurus, our simulations assume one, three or five layers of birefringent sapphire, which we will refer to as BR1, BR3, and BR5 respectively.Each HWP is covered in three layers of dielectric anti-reflection coating3 , whose dielectric constant gradually increases from that of free-space to that of sapphire (n 1 , n 2 , n 3 ) = (1.27,1.98, 2.86).The corresponding thicknesses of the AR coating layers are 0.34, 0.21, and 0.18 mm.BR3 is based on the POLARBEAR-2 HWP [29], BR5 on a potential HWP for LiteBIRD [34].At a fixed frequency, ν, the optimal thickness, d, of a sapphire layer is given by: with c the speed of light in a vacuum and ∆n the difference between the two indices of refraction of sapphire.In our simulations, ν is the mid-point between the two bands of LF or HF.All the sapphire layers in our HWP models have the same thickness for a given receiver: the values are given in Table 1.We calculate the Mueller matrices for stacks of dielectric and birefringent materials using software described in [35].

The beamconv library
Our modelling is built around beamconv [36], a lightweight algorithm/python library to simulate time-ordered data from beam-convolved CMB maps.In beamconv, both the beam and the sky map are decomposed into spin-weighted spherical harmonics.The coefficients can be linearly combined into spin maps, which are then sampled by beamconv following a sequence of pointings that represent the telescope's scan strategy, appropriately modulated by the HWP Mueller matrix.We therefore obtain time-ordered data (TOD) for each detector that gets projected into a map using a simple binning scheme.A real flight will discard a moderate fraction of this TOD due to fridge cycles, antenna noise, and other unavoidable instrumental noise, but ignoring this won't bias our simulations.To simulate the variation in either beam or sky behavior over a band, beamconv simulations for several sub-frequencies within a band can be co-added into a band-averaged map.When testing for frequency-dependent systematics, like HWP non-idealities or the influence of dust, we make nine simulations spaced every 5 GHz between 130 and 170 GHz.For each frequency, the behavior of the HWP and the sky will be different: we do not simulate variations of the beam within the band.beamconv has been used in the recent past to estimate the effect of HWP non-idealities on B-mode and cosmic birefringence searches [37,38].
4 Analysis of output maps

Power spectrum estimation
Power spectrum estimation is fundamentally limited by cosmic variance [39,40].Beyond cosmic variance, power spectrum estimators acting on cut-sky data are sensitive to several effects due to degeneracies among spherical harmonics when measured on only part of the sky.To choose a specific estimator, we quantify the effects of mode-mixing and EB leakage.We simulate 100 realisations of the CMB based on Planck's best fit power spectrum [41] and then mask areas outside of Taurus' scan range (see Figure 3).We also mask the 30 % of the sky most contaminated by the galactic plane, resulting in around 44 % of the sky remaining unmasked.We estimated the power spectrum of the masked maps, downgraded to a resolution of NSIDE=8, using three estimators: Polspice [42], xQML [43], and Namaster [44].We decide to use Namaster (without B-mode purification) for the rest of our analysis, as it recovers the input EE spectrum with relatively little bias over the multipole range 2 ≤ ℓ ≤ 22. xQML, while more consistent in its estimations, exhibits a bias.We note that after mode-coupling is corrected, the auto-spectra in Namaster are no longer necessarily positive on the very large angular scales.Since the reionisation bump in the EE spectrum is fairly narrow, we decide to divide the spectrum into bins of width ∆ℓ = 5.The reionisation signal will be contained within the first three bins, while other bins can be used for calibration.
To quantify the effect of non-idealities, we will look at the differences between maps simulated with a systematic (non-ideal maps) and other maps simulated without it (ideal maps).By estimating the power spectrum of the difference maps we can gauge the effect of the systematic at the C ℓ level.We can then compare that effect to the uncertainty in the EE power spectrum that we expect from sample variance for the observed sky fraction: We build our sample variance estimate following the Planck best-fit value τ = 0.0543 [2] and f sky = 0.44.

Half-wave plate efficiency
All microwave telescopes rely on some form of instrument calibration to make sense of the data they collect.In particular, they need to calibrate the beam of the instrument and the efficiency ϵ with which the instrument turns incoming radiation into measured Stokes I, Q, and U .This process involves the use of well-understood sources such as planets [24], thermal sources on drones and balloons [45], the CMB dipole [46], or lab-based holography [47] to provide estimates of the efficiency and construct beam maps.In the presence of a polarisation modulator like a HWP, the efficiency will depend on the non-idealities of that modulator, as represented by its Mueller matrix.The efficiency will be different for total intensity and polarisation: we can define an ϵ I for intensity and an ϵ P for the polarisation modulation efficiency, with In Figure 5 we plot ϵ P for the three HWP models that we study following the formula given in [34]: where the terms in the above equation all correspond to a frequency-dependent Mueller matrix element.Equation 4.1 describes the fraction of the incoming polarized light that will be modulated by 4θ, like in the idealized case of Eq. 3.2.The ϵ P of the one-layer HWP varies within each band and can be up to 40 % lower than for the three or five-layer AHWP.
The right panel of Figure 5 showcases a peculiarity of AHWPs.While it will have a near-ideal ϵ P over a large frequency domain, an AHWP will exhibit a frequency-dependent phase shift δα(ν), meaning that the plane of polarisation of light exiting the HWP has been rotated, changing some of the Q into U .This can be accounted for during map-making: we correct the HWP angle by that phase shift [48,49].Since it varies over the 150 GHz We have offset the 3BR phase by 32.67 • for ease of reading.The phase-shift is zero for 1BR and negligible for 5BR.
frequency band, the correction operation is accurate for the band-averaged value δα of δα(ν) where g(ν) is the spectral response of the telescope and detector, and I(ν) the intensity of the sky signal.In the above expression, we assume a spectral response function that defines a bandwidth ∆ν centered on ν 0 .The value of δα(ν) and g(ν) can be measured during pre-flight calibration.During observations, however, the sky signal, I(ν), represents the location-dependent sum of different components (including synchrotron emission, dust, CMB) in proportions that depend on the sky position.This means that the best-suited value of δα will depend on our assumptions about the spectral energy density of the sky we are observing.Taking for instance the BR3 HWP model, δα is 32.67 • for a pure CMB sky but 32.94 • for the dust emission described by PySM's d9 model.

Corrections at the power-spectrum level
The multiplicative gain due to ϵ I and ϵ P will propagate to maps and power spectra [50].
Ignoring the efficiency issue leads to smaller inferred amplitudes for the CMB anisotropies which obviously impacts attempts to measure τ or the other cosmological parameters.We therefore try to account for it by defining a HWP-related gain G XY : where G XY is the square root of the relative amplitude of the cross-or auto-spectrum XY of simulated maps due to the presence of the HWP non-idealities, and (ℓ 1 , ℓ 2 ) is the calibration range.Multiplying the simulated maps made with a non-ideal HWP by G XY compensates for the HWP efficiency losses.G T T is similar to ϵ −2 I , while ϵ P can be corrected by a good estimate of G EE and/or G BB .The difference maps m diff are then given by: (4.4) A real experiment will not know CXY ;ideal ℓ , but it would be reasonable to construct this quantity by convolving a known C XY ℓ from another experiment's spectra or maps with the effective beam transfer function B ℓ for the experiment: If we rely on another experiment, the calibration range has to be at multipoles where both Taurus and the other experiment can measure the T T or EE spectrum accurately.Beyond multipoles of ℓ ∼ 400, Taurus' angular resolution limits its sensitivity.At the other end, setting the calibration range to be the same as the multipoles we are trying to characterize (ℓ < 22) can intertwine our analysis with the systematics of previous experiments in a very direct way.It is difficult to make an independent measure of the EE spectrum on large angular scales if our calibration efforts depend on the measured Planck EE spectrum on those scales.In this paper, we settle on ℓ 1 = 50 and ℓ 2 = 100, where Planck was able to measure the polarisation power spectrum with high accuracy [41].
If the beam transfer function is correct, Eq. 4.3 remains unchanged.However, if the wrong B ℓ is assumed, then CXY ;ideal ℓ will have an ℓ-dependent bias that will propagate to G XY .This could happen, for instance, if the assumed B ℓ for the experiment ignores the beam sidelobes.Then their Cℓ s will have a different gain in the "science" bins at than in the calibration range, and our procedure to estimate the gain will introduce an error in the experiment's estimation of the CMB power spectrum.Let us try to illustrate this with a simple case: suppose our ideal model assumes a Gaussian beam, denoted with g, with FWHM f 1 .The transfer function for that beam is: Suppose, however, that the true beam had an extra Gaussian sidelobe aligned with the main beam, denoted with sl, with amplitude ε and FWHM f 2 .This can be seen as a reasonable approximation of the "tails" of a real beam beyond the main lobe.The transfer function of the true beam is: where B sl ℓ follows the form set in Eq. 4.6 but with the FWHM replaced by f 2 .Then, the bias per multipole is: The transfer functions for our fitted and PO beams, visible in Figure 6, can be approximated with ε = 0.03, f 1 = 33 ′ and f 2 = 5 • .The average bias over the calibration range (50 ≤ ℓ ≤ 100) is approximately −1.5 × 10 −3 .In our "science" multipole range, the bias in a 7 ≤ ℓ < 12 multipole bin for those two beams would however be −4.4×10−2 .Therefore, beam mismatch during gain calibration could lead to more than 4 % multiplicative bias in our estimate of the power spectrum.
A power spectrum prior can be used for more than just efficiency calibration.It should also be possible to calibrate the polarisation angles of the detectors by attempting to null the CMB EB power spectrum [51], which is supposed to be zero in ΛCDM.In the same vein, [38] discussed the degeneracy between non-zero C EB ℓ and an incorrect estimation of the HWP angle.We could therefore imagine determining δα by minimizing the EB power.That is particularly attractive as polarized point sources that can be used for calibration are very rare in the microwave sky.There are two downsides to that approach: first, it blinds us to detecting the possible non-zero CMB EB correlation caused by cosmic birefringence [52,53].Furthermore, the EB power spectrum could contain a non-zero contribution from Galactic dust emission [54], so EB minimization would not be a good option for the other bands of Taurus, which are dominated by dust emission.

Results
There is no shortage of potential instrument non-idealities that can complicate data analysis.With the model described in Section 3, we focus on beam and HWP non-idealities.Beam non-idealities include deviations from the assumed beam model, errors in pointing, and polarisation angle errors.HWP non-idealities cover both the effect of non-ideal elements in the Mueller matrix and the coupling between those elements and the beam.We present the power spectrum level residuals in the multipole bin 7 ≤ ℓ < 12, comparing them to value of the sample variance for C EE ℓ in that bin for f sky = 0.44, which is about 4.2 × 10 −5 µK 2 .Dust emission has a different SED than the CMB.As discussed in Section 4.2, this means that the optimal AHWP δα in a given frequency band will vary depending on the observed component.Furthermore, the dust is not distributed isotropically on the sky.Most of it is in or near the Galactic plane, which we mask before taking the power spectra.However, sidelobes could lead to emission from within the mask to leak into our observations.Therefore we perform the simulations of beam non-idealities, reflection ghosts, and non-ideal HWPs for the PySM d9 sky model.

Optical ghosts
Before it reaches a detector, a photon passing through a refracting telescope will interact with several optical elements, including the vacuum window of the cryostat, optical filters, the HWP, and the lenses.Each change of medium will lead to some reflection.This can lead to optical ghosts, weaker images of the main beam which in some cases are mirrored across the boresight.For a beam with pointing (x, y) with respect to the telescope's boresight, the associated ghost beam's centroid will be (−x, −y).While we sandwich lenses and HWP with anti-reflection coatings (ARCs) to suppress those reflections, it is challenging to totally avoid ghosting.With a wide field of view like ours, ghosting will introduce spurious correlations at large angular scales for a single detector.To estimate the effect of ghosting, we create simulations where each beam has a 1 % reflection ghost of the main beam.Because of our symmetric focal plane, the ghosts are symmetrized too.In a real experiment, dead detectors would break that symmetry and the ghosting would have a larger impact: therefore, for that part of the simulation work, we kill 7 randomly selected detector pairs out of our 49 (a roughly 14 % fraction).In our CMB simulations, the ghost beams have a negligible impact (at most a few 10 −4 of the sample variance target, 4.2×10 −5 µK 2 ) for all beam models or gain calibration options.We therefore will not discuss those simulations further.In Table 2 we show the impact of the ghost beams for dust scans is about 6 % of our self-imposed sample variance target but is mostly unaffected by the choice of gain parameter: we only notice that the fully co-polar beam models (Gaussian, fitted) gets slightly increased residuals when calibrating on EE, suggesting the pickup of the additional galactic dust by the sidelobes of the ghost beams can be suppressed by calibration at intermediate scales.

Pointing and polarisation angle errors
Small to intermediate-level pointing reconstruction errors, also known as beam centroid errors, are known to primarily impact power spectrum reconstruction on intermediate angular scales (500 ≤ ℓ ≤ 1000) [55].However, a recent reassessment in [56] did show that sub-pixel errors have an impact on the T T spectrum at scales ℓ ≤ 500.Polarisation angle errors, both random and systematic, have been studied extensively in [57].We group them with pointing errors as a polarisation-sensitive detector's pointing is defined by three angles on the sphere: (θ, ϕ, ψ), where θ and ϕ are the two normal spherical coordinates and ψ indicates the angle between the detector's orientation and the local meridian.Alternatively, we can express these angles relative to the orientation of the telescope's boresight: (θ, ϕ, ψ) → (∆az, ∆el, ξ).The topic of polarisation angle errors has been intensely discussed recently due to the interest in cosmic birefringence [52,53,58].We simulate two modes of pointing error: either a common offset by 1 ′ in azimuth and elevation for all detectors analogous to a bias in the boresight position or a random error for each detector drawn from a Gaussian distribution with 1 ′ standard deviation.In the same vein, we examine the impact of improperly characterizing the polarisation angle of the detectors during the flight.Here, the characteristic error is estimated to be as large as 1 • .
We find that both effects are much smaller than the sample variance in the 7 ≤ ℓ < 12 bin, with the random az-el offset being at most 2 % of the sample variance if the simulation includes PO beams with sidelobes and we attempt to calibrate on the EE spectrum.This is consistent with the simulations done for other balloon CMB experiments, like SPIDER [59,60], and LSPE [61].The polarisation angle offset also returns very small residuals.The negative sign is due to intrinsic variance from the Namaster mode decoupling algorithm, but the "true" value should be around the same amplitude.

Beam non-idealities
Perturbations to the Gaussian beam model are to be expected in realistic instrument designs.The main causes are non-uniform illumination of the aperture and extra reflections within the instrument.Ray-tracing or physical optics simulations can help us understand those extra components.While they are much weaker than the main beam, the large angular extension of sidelobes will cause a detector to pick up some stray radiation from other much brighter sources, like the Galaxy or the ground.Our physical optics (PO) simulations of the Taurus optics show an extended sidelobe with an amplitude between −40 dB and −60 dB, visible in Figure 4.The azimuthally averaged beam main beam profile is similar to a Gaussian beam with a 33 ′ FWHM, as can be seen in Figure 6.
When subtracting maps simulated with the Gaussian beam from maps simulated with the fitted beam model or the PO beams, this results in some residuals in the EE power spectrum.We can observe those in Figure 7, where we plot the power spectrum of the residual dust maps up to ℓ = 100.If we do not calibrate the HWP-related gain (i.e.G = 1), we see residuals that increase fairly monotonically until ℓ ∼ 50 and then plateau, more marked  3. Ratio of EE power spectrum residuals to sample variance in the 7 ≤ ℓ < 12 bin, as a fraction of sample variance for f sky = 0.44, due to substraction of a simulation with Gaussian beams.We bold the configurations where the residuals are larger than the sample variance.
for the PO+sidelobe case than the fitted or PO case.This is consistent with residuals from the difference maps getting deconvolved by the B ℓ associated with the main lobe (which is a more inadequate approximation in the PO+sidelobe case).In the calibrated cases, we can see the bias scenario we described in Section 4.3: residuals for both T T and EE calibration are suppressed between ℓ = 50 and ℓ = 100 compared to the un-calibrated case, but are enhanced at low multipoles.The value in the 7 ≤ ℓ < 12 bin can be seen in Table 3.The calibration mismatch leads to similar residuals for the T T and EE calibrations, which is to be expected since the HWP is kept ideal.For CMB simulations, the residuals are about a percent of the sample variance limit for the three beam models.We do not see a difference between T T and EE based calibration, as the ideal behavior of the HWP stops any gain difference between those two spectra.
For the dust map, the residuals are of the order of our sample variance target.On top of that, we see larger residuals for cases of the fitted beam and the physical optics beam with sidelobes, highlighting that sidelobes beyond a few FWHMs will significantly impact component characterization and separation.This stresses the importance of beam characterisation and calibration for CMB experiments like Taurus.

Half-wave plate non-idealities
We evaluate residuals due to HWP non-idealities for each of the beam models.As we are now probing a frequency-dependent effect, we apply the band-averaging method described in Section 3.4.
In Table 4 we see that, if the input beam is fully known, the performance of the HWP is only weakly dependent on the specific beam model chosen: for a given HWP, sky component, and gain determination method, the residuals are similar for all four beam models.Beyond that, we are able to see clear differences in performance between the three HWP models depending on the choice of calibration scheme.The residuals for the single layer model, in the top two rows, are about six to seven times as large as for the three-and five-layer AHWPs in the absence of gain calibration.When calibrating the gain on T T , the residuals decrease for the three HWP models: moderately for BR1 and fairly aggressively for BR3 and BR5.Calibrating on EE, the smallest residuals for the CMB simulation are achieved by the BR3 model, while the BR1 and BR5 models have percent-level residuals for the dust simulations.This suggests the BR5 model is the most suitable if Taurus plans to calibrate on the T T spectrum, and that the polarisation efficiency of both achromatic models can be calibrated effectively against the EE power spectrum.4. Residuals due to HWP non-idealities for all four beam configurations.Each set of three columns is a different beam configuration, each line corresponds to a HWP and sky model."No", "T T " and "EE" refer to different estimation methods of the gain as defined in Section 4.3.Results are expressed as a fraction of the sample variance in the 7 ≤ ℓ < 12 bin.We bold the configurations where the residuals are larger than the sample variance.

Half-wave plate angle error
The frequency-dependent phase shift of the AHWP models leads us to question another HWP-related non-ideality: what happens if the angle determination of the HWP is wrong?For an ideal HWP, an extra rotation of the HWP by an angle α is degenerate with a rotation of the polarisation-sensitive detectors by an angle 2α [57].In a real HWP with a nondiagonal Mueller matrix M , the M IQ , M IU , M QI and M U I terms will also be modulated by the rotation of the HWP, therefore an error on the HWP angle will propagate differently than an error on the detector polarisation angle.We evaluate the effect of a 0.5 • error in determining the angle of each HWP model for all four beam models.The CMB results, visible in Table 5, are nearly identical to those achieved with a detector polarisation angle offset in Section 5.2, suggesting that the impact of the HWP angle error on the modulation of the off-diagonal terms is minimal.However, we can observe that a half-degree error on the HWP's angle determination leads to residuals in the dust power spectrum that are more than half of the sample variance.This once again stresses the importance of calibration, this time of the HWP's angles.We note however that the error we chose is a very conservative number considering the performance achieved by SPIDER [62] where the HWP rotation mechanism's angle was known to ±0.1 • .

Beam and half-wave plate interplay
What if we have imperfect knowledge of both the beam and the HWP model?In the absence of gain calibration, the residuals are similar to those of Table 4, suggesting that the HWP non-ideality is the main systematic.That is sensible as the residuals for the uncalibrated cases in Table 3 are minor compared to the residuals caused by the HWP.When performing gain calibration, the residuals associated with the BR1 plate are divided by a factor of 3 when calibrating on T T and a factor of 10 when calibrating on EE, making the residuals for the dust map equivalent in size to the sample variance associated with f sky = 0.44.Interestingly, this makes the T T residuals smaller than in the case where we evaluate only the effect of HWP non-ideality.
For the BR3 HWP, the residuals are lowest when calibrating on the T T power spectrum, no matter the beam model chosen.In fact, it is the only HWP and calibration configuration with residuals smaller than the cosmic variance target for the "PO+Sidelobe" beam model.The "Fitted" and "PO+Sidelobe" beam models create higher residuals than in the case of a PO beam with no sidelobes, suggesting that characterizing sidelobes should be a more important focus for Taurus than getting an exact mapping of the main beam.The ∼ 30 % difference between the residuals of the PO+sidelobes model and the azimuthal-symmetric fitted beam do however concede that there is some importance to non-azimuthally symmetric beam modes, and that a calibration campaign should not ignore them.
For the BR5 HWP model, gain calibration has trouble reducing the power spectrum residuals due to the beam mismatch issue.Whereas, in Table 4, T T or EE calibration reduced the dust map residuals to less than 1 % of the sample variance, here the residuals are between 1 and 1.3 times the sample variance depending on the beam model, with slightly smaller residuals when calibrating on T T .The PO beam without sidelobes has once again smaller residuals, reinforcing the idea that it is the beam mismatch at low ℓ that makes calibration difficult.

Conclusion
On average, the residuals we see at power spectrum level for the CMB sky scans are smaller than the target we set ourselves at the start of this study: that the systematic errors due to optical effects be often smaller than the cosmic variance associated with the Planck EE spectrum at a 44 % sky fraction.In general, being able to reliably calibrate on the EE power spectrum at smaller angular scales reduces residuals for HWP-related effects.The nonideality with the largest effect, beam model mismatch, is unfortunately insensitive to that approach.As could be expected, replacing the single layer HWP with an AHWP increases the polarisation modulator's performance when scanning the CMB, although calibrating on the EE power spectrum brings the behavior of the one-layer model in line with the two AHWPs.
However, we should interrogate the possibility of such a calibration: it relies on having a simulation of Taurus' observations of the sky with a perfect instrument, but how will we know what the true sky looks like?If we use Planck's EE spectrum, then we will be implicitly relying on Planck 's model of its polarisation sensitivity for our analysis.This problem is more important for the three-layer AHWP, as its scans had smaller residuals when calibrating on EE than on T T .The five-layer AHWP therefore seems like the best option, as it performs very well when calibrating on the temperature power spectrum, with the residual power in the 7 ≤ ℓ ≤ 12 bin of the EE power spectrum at most 1 % of the sample variance in that bin for dust foregrounds.This is an easier feat to accomplish for us: we could conceivably develop a reliable temperature calibration from the CMB dipole, or use Planck 's T T spectrum but retain freedom for our polarisation analysis.
The other concerning effect appears to be the coupling of sidelobes to the polarized dust emission of our galaxy.We note that imperfect knowledge of the beam model, and in particular of its far-sidelobes, is one of the main sources of systematic error for LiteBIRD's B-mode effort [63].Taurus should have, at a minimum, an estimate of its average sidelobe amplitude within the instrument model akin to the fitted azimuthally symmetric beam we used in this paper.This is only made more urgent by the fact that T T or EE calibration increases the amplitude of the beam residuals, to the point where the residuals associated with the five-layer AHWP are larger than the sample variance limit.
On the other hand, we have only simulated the 150 GHz band and are not forecasting how successful Taurus' component separation efforts will be once the four bands are included.We will however note that other collaborations that have either conducted simulations [61,64] or gathered data [65] appear to state with confidence that beam systematics can be kept under control when studying large-scale CMB polarisation by using astrophysical sources to characterize the beam and the polarisation angle.Holographic measurements [47] could also provide pre-flight beam characterization.

General conclusion
In summary, we have simulated a number of systematics related to the beams or the HWP for a one-month flight of the Taurus balloon.Using both Gaussian and PO beam simulations of 49 detector pairs in the 150 GHz band, we generated timestreams and maps of simulated CMB and dust emission with beamconv and evaluated the impact of non-idealities in the EE power spectrum.We saw that the pointing and polarisation angle accuracy targets of Taurus led to negligible amounts of error in the EE spectrum.We have found that the three HWP models perform differently, with the five-layer AHWP performing the best if the beam is fully known while the three-layer AHWP seems to deal with beam mismatch better.We also observed that the coupling of beam sidelobes to polarized dust emission makes it difficult to self-calibrate the HWP polarisation efficiency if a Gaussian beam model is assumed.In future work, we would like to consider the impact of other systematics, like optical loading, and include more realistic timestreams with correlated detector noise.We will also include the other frequency bands of Taurus.derived using the healpy and HEALPix package.Taurus is supported in the USA by NASA award number 80NSSC21K1957.JEG acknowledges support from the Swedish Research Council (Reg.no.2019-03959) and the Swedish National Space Agency (SNSA/Rymdstyrelsen).This work is in part funded by the European Union (ERC, CMBeam, 101040169).

Figure 1 .
Figure 1.Left: the CMB EE power spectrum for 2 ≤ ℓ ≤ 200 for a variety of values of τ within Planck's 1σ uncertainty[2] (color bar), with the product A s e −2τ kept constant.Right: the relative variation from the Planck nominal value for 2 ≤ ℓ ≤ 24.The grey lines represent the sample variance associated with the EE spectrum for different observed sky fractions: the dotted line is for f sky = 0.4, the dashed line for f sky = 0.7 and the full line is the cosmic variance limit on σ(τ ).In our simulations, we estimate the power spectrum on roughly 40 % of the sky.

Figure 3 .
Figure 3. Left: Normalized hits map, in equatorial coordinates, for the scanning strategy described in Section 2 for 49 detector pairs sampling at 50.01 Hz.Grey pixels are not observed.The cyan area falls within our galactic mask and covers the 30 % of the sky most contaminated by the Galaxy according to Planck.Right: Condition number map.

Figure 4 .
Figure 4. Azimuthally-averaged Stokes I beam intensity for pixels of our simulated refractor.In blue, a 33 ′ FWHM Gaussian beam.The green and red lines are the azimuthally averaged profiles for the PO simulations.In orange, an analytical fitting of the PO beam for the central detector.To evaluate the relative contribution of main lobe and side lobe for the PO beams, we can run simulations that ignore any beam power beyond 3 • of the beam centroid.This is represented by the gray outline.

Figure 5 .
Figure 5. Figures of merit as a function of frequency from the Mueller matrices of the three HWP models studied, computed using[35].1BR, 3BR, 5BR refer to the number of birefringent layers in the HWP model.The shaded regions are the two LF bands of Taurus.Left: Polarisation modulation efficiency.Right: frequency-dependent phase shift.We have offset the 3BR phase by 32.67 • for ease of reading.The phase-shift is zero for 1BR and negligible for 5BR.

Figure 6 .
Figure 6.The azimuthal, m = 0, component of the beam transfer functions B ℓ used in our simulations, normalized to their average in the grey area, our calibration range.We can see how the B ℓ diverge from it in the colored areas from ℓ = 7 to ℓ = 22 that correspond to our science bins.

2 )Figure 7 .
Figure 7. EE power spectrum of residual dust maps due to beam mismatch.The light grey line indicates the binned sample variance for f sky = 0.7, the dark grey line for f sky = 0.4.Each subfigure corresponds to a different beam model, indicated at the top.Within a sub-figure, each color and marker style refers to different gain calibrations.No calibration is attempted for the blue round markers (we assume the HWP-related gain is 1), while T T and EE refer to the self-calibration G XY of Section 4.3.
• × 20 • square grid on the sky.The diagonal therefore covers the instrument's 28 • field of view.Each detector pair has two detectors that are aligned with ±Q or ±U (i.e. the detectors' polarisation angles are 90 • apart).Neighboring detector pairs are staggered by 45•to alternate between Q and U sensitivity.The Taurus optical design is still evolving, but the baseline refractor designs gives the beams a FWHM of at least 22.5 ′ (for the 220 GHz band) and at most 33 ′ (for the 150 GHz band).Values for all four frequency bands are given in Table

Table 1 .
Frequency-dependent properties of the instrument model.

Table 2 .
Ratio of power spectrum residuals to sample variance due to ghost beams for ideal HWP and PySM d9 sky model, 7 ≤ ℓ < 12 as a fraction of the sample variance in that bin.

Table 5 .
Residuals in the 7 ≤ ℓ < 12 bin, as a fraction of sample variance, due to 0.5 • HWP angle error (see Section 5.4.1).Lines and columns follow the definitions from Table4.The small negative values are due to cut-sky effects within Namaster and are consistent with zero power.

Table 6
we show the residuals versus simulations made with Gaussian beams and ideal HWPs.In

Table 6 .
7 ≤ ℓ < 12 multipole bin of the EE power spectrum, as a fraction of the sample variance, for residual maps of the HWP+beam systematic (see Section 5.5).