QuickPol: Fast calculation of effective beam matrices for CMB polarization
Key Words.:
Cosmology, polarization, systematic effectsCurrent and planned observations of the cosmic microwave background (CMB) polarization anisotropies, with their ever increasing number of detectors, have reached a potential accuracy that requires a very demanding control of systematic effects. While some of these systematics can be reduced in the design of the instruments, others will have to be modeled and hopefully accounted for or corrected a posteriori. We propose QuickPol, a quick and accurate calculation of the full effective beam transfer function and of temperature to polarization leakage at the power spectra level, as induced by beam imperfections and mismatches between detector optical and electronic responses. All the observation details such as exact scanning strategy, imperfect polarization measurements, and flagged samples are accounted for. Our results are validated on Planck high frequency instrument (HFI) simulations. We show how the pipeline can be used to propagate instrumental uncertainties up to the final science products, and could be applied to experiments with rotating halfwave plates.
1 Introduction
We are now entering an era of precise measurements of the cosmic microwave background (CMB) polarization, with
potentially enough sensitivity to detect or even characterize the primordial
tensorial modes, the smoking gun of inflation
(e.g., Zaldarriaga & Seljak (1997) and references therein). This raises
expectations about the control and the correction of contaminations by
astrophysical foregrounds, observational features, and instrumental
imperfections. As it has in the past, progress
will come from the synergy between instrumentation and data
analysis. Improvements in instrumentation call for improved precision in
final results, which are made possible by improved algorithms and the
ability to deal with more and more massive data sets. In turn,
expertise gained in data processing allows for better simulations that lead to new
instrument designs and better suited observations. An example of such joint
developments is the study of the impact of optics and electronicsrelated
imperfections on the measured CMB temperature and polarization angular power
spectra and their statistical isotropy. Systematic effects such as beam noncircularity,
response mismatch in dual polarization measurements and scanning strategy imperfections,
as well as how they can be mitigated, have been extensively studied
in the preparation of forthcoming instruments (including, but not limited to Souradeep & Ratra 2001; Fosalba et al. 2002; Hu et al. 2003; Mitra et al. 2004, 2009; O’Dea et al. 2007; Rosset et al. 2007; Shimon et al. 2008; Miller et al. 2009a, b; Hanson et al. 2010; Leahy et al. 2010; Rosset et al. 2010; Ramamonjisoa et al. 2013; Rathaus & Kovetz 2014; Wallis et al. 2014; Pant et al. 2016),
and during the analysis of data collected by
WMAP
At the same time, several deconvolution algorithms and codes have been proposed
to clean up the CMB maps from such beamrelated effects prior to the computation
of the power spectra, like PreBeam (ArmitageCaplan & Wandelt 2009),
ArtDeco (Keihänen & Reinecke 2012), and in Bennett et al. (2013) and
Wallis et al. (2015).
Finally, in a related effort, the FEBeCoP pipeline,
described in Mitra et al. (2011) and used in
Planck data analysis (Planck 2013IV 2014; Planck 2013VII 2014),
can be seen as a convolution facility, by providing,
at arbitrary locations on the sky, the effective beam maps and
point spread functions of a detector set, which,
in turn, can be used for a MonteCarlo based description
of the effective beam window functions for a given sky model.
In this paper, we introduce the QuickPol pipeline, an extension to polarization of the Quickbeam algorithm used in Planck 2013VII (2014). It allows a quick and accurate computation of the leakage and crosstalk between the various temperature and polarization power spectra (, , , , etc.) taking into account the exact scanning, sample flags, relative weights, and scanning beams of the considered set(s) of detectors. The end results are effective beam matrices describing, for each multipole , the mixing of the various spectra, independently of the actual value of the spectra. As we shall see, the impact of changing any timeindependent feature of the instrument, such as its beam maps, relative gain calibrations, detector orientations, and polarization efficiencies can be propagated within seconds to the final beam matrices products, allowing extremely fast MonteCarlo exploration of the experimental features. QuickPol is thus a powerful tool for both real data analysis and forthcoming experiments, simulations and design.
The paper is organized as follows. The mathematical formalism is exposed in Section 2 and analytical results are given in Section 3. The numerical implementation is detailed in Section 4 and compared to the results of Planck simulations in Section 5. Section 6 shows the propagation of instrumental uncertainties. We discuss briefly the case of rotating halfwave plates in Section 7 and conclude in Section 8.
2 Formalism
2.1 Data stream of a polarized detector
As usual in the study of polarization measurement, we will use
Jones’ formalism to study the evolution of the electric
component of an electromagnetic radiation in the optical system. Let us
consider a quasi monochromatic
A rotation of the optical system by around the axis can be seen as a rotation of both the orientation and location of the incoming radiation by in the detector reference frame, and the same input radiation is now received as
(1) 
with
(2)  
(3) 
and the sign representing the adjoint operation, which for a real rotation matrix, simply amounts to the matrix transpose. The measured signal is
(4) 
with
(5) 
We now introduce the Stokes parameters of the input signal (dropping the dependence on )
(6) 
and of the (unrotated) instrument response
(7) 
to obtain
(8) 
With the rotated instrument response:
(9a)  
(9b)  
(9c)  
(9d) 
Following Rosset et al. (2010), we can specify the instrument as being a beam forming optics, followed by an imperfect polarimeter in the direction , with , and having an overall optical efficiency :
(10) 
with
(11) 
for . The Stokes parameters of the instrument are then for .
If the beam is assumed to be perfectly copolarized, that is, it does not alter at all the polarization of the incoming radiation, with and , then , , and , , , Eqs. (8, 9) become
(12) 
where
(13) 
is the polar efficiency, such that with for a perfect polarimeter and for a detector only sensitive to intensity. In the case of Planck high frequency instrument (HFI), Rosset et al. (2010) showed the measured polarization efficiencies to differ by 1% to 16% from their ideal values, with an absolute statistical uncertainty generally below 1%. The particular case of copolarized beams is important because in most experimental setups, such as Planck, the beam response calibration is done on astronomical or artificial far field sources. Well known, compact, and polarized sources are generally not available to measure and and only the intensity beam response is measured. In the absence of reliable physical optics modeling of the beam response, one therefore has to assume and to be perfectly copolarized.
So far, we have only considered the optical beam response. We should also take into account the scanning beam, which is the convolution of the optical beam with the finite time response of the instrument (or its imperfect correction) as it moves around the sky, as described in Planck 2013VII (2014) and Planck 2015VII (2016). These time related effects can be a major source of elongation of the scanning beams, and can increase the beam mismatch among sibling detectors. If one assumes the motion of the detectors on the sky to be nearly uniform, as was the case for Planck, then optical beams can readily be replaced by scanning beams in the QuickPol formalism.
2.2 Spherical harmonics analysis
We now define the tools that are required to extend the above results to the full celestial sphere. The temperature is a scalar quantity, while the linear polarization is of spin , and the circular polarization is generally assumed to vanish. They can be written as linear combinations of spherical harmonics (SH):
(14)  
(15) 
although one usually prefers the scalar and fixed parity and components
(16) 
such that for . In other terms
(17) 
with
(18) 
The sign convention used in Eq. (16)
is consistent with Zaldarriaga & Seljak (1997) and the HEALPix
The response of a beam centered on the North pole can also be decomposed in SH coefficients
(19)  
(20) 
while the coefficients of a rotated beam can be computed by noting that under a rotation of angle around the direction , the SH of spin transform as
(21) 
The elements of Wigner rotation matrices are related to the SH via (Challinor et al. 2000)
(22) 
with .
If the beam is assumed to be copolarized and coupled with a perfect polarimeter rotated by an angle , such that in cartesian coordinates (or in polar coordinates), simple relations between and can be established. For a Gaussian circular beam of full width half maximum (FWHM) and of throughput Challinor et al. (2000) found
(23a)  
(23b) 
The factor in Eq. (23b) is such that for and for , and will be assumed to be from now on. For a slightly elliptical Gaussian beam, Fosalba et al. (2002) found
(24) 
while we show in Appendix G that Eq. (24) is true for arbitrarily shaped copolarized beams. This result can also be obtained by noting that an arbitrary beam is the sum of Gaussian circular beams with different FWHM and center (Tristram et al. 2004), each of them obeying Eq. (23b).
The detector associated to a beam is an imperfect polarimeter with a polarization efficiency and the overall polarized response of the detector, in a referential aligned with its direction of polarization (the socalled Pxx coordinates in Planck parlance), reads
(25) 
so that
(26) 
We introduced to distinguish it from the value used in the mapmaking, as described below.
2.3 Map making equation
A polarized detector pointing, at time , in the direction on the sky, and being sensitive to the polarization with angle with respect to the local meridian, measures
(27) 
The factor present in Eq. (8) is assumed to be absorbed in the gain calibration, performed on large scale temperature fluctuations, such as the CMB solar dipole (Planck 2013VIII 2014), and we assumed the circular polarization to vanish. With the definitions introduced in Section 2.2, this becomes
(28) 
The mapmaking formalism is set ignoring the beam effects, assuming a perfectly copolarized detector and an instrumental noise (Tristram et al. 2011, and references therein), so that, for a detector , Eq. (12) becomes
(29) 
where the leading prefactors are here again absorbed in the gain calibration. Let us rewrite it as
(30) 
with (Shimon et al. 2008)
(31)  
(32) 
and . Assuming the noise to be uncorrelated between detectors, with covariance matrix for detector , the generalized least square solution of Eq. (29) for a set of detectors is
(33) 
Let us now replace the ideal data stream (Eq. 29) with the one obtained for arbitrary beams (Eq. 27) and further assume that the noise is white and stationary with variance , so that . Let us also introduce the binary flag used to reject individual time samples from the mapmaking process; Eq. (33) then becomes
(34)  
(35) 
We have assumed here the pixels to be infinitely small, so that, starting with Eq. (28), the location of all samples in a pixel coincides with the pixel center. The effect of the pixel’s finite size and the socalled subpixel effects will be considered in Section 3.5.
2.4 Measured power spectra
To compute the crosspower spectrum of any two spin and maps, we first project each polarized component of on the appropriate spin weighted sets of spherical harmonics,
(36) 
and average these terms according to
(37)  
(38) 
where Eq. (85) was used. The detailed derivation of this relation and its associated terms is given in Appendix A. Suffice it to say here that terms are either or , terms are inverse noiseweighted beam multipoles, and terms are effective weights describing the scanning and depending on the direction of polarization, hit redundancy (both from sky coverage and flagged samples), and noise level of detector .
Equation (38) is therefore a generalization to noncircular beams of the pseudopower spectra measured on a masked or weighted map (Hivon et al. 2002; Hansen & Górski 2003), and extends to polarization the Quickbeam noncircular beam formalism used in the data analysis conducted by Planck 2013VII (2014). It also formally agrees with Hu et al. (2003)’s results on the impact of systematic effects on the polarization power spectra, with the functions absorbing the systematic effect parameters relative to detector . In the next sections, we present the numerical results implied by this result and compare them on fullfledged PlanckHFI simulations.
3 Results
We now apply the QuickPol formalism to configurations representative of current or forthcoming CMB experiments, and to a couple of idealized test cases for which the expected result is already known, as a sanity check. The effect of the finite pixel size is also studied.
3.1 A note about scanning strategies
To begin with, let us
consider the scanning strategy of Planck and of another satellite mission
optimized for the measurement of CMB polarization.
Figure 1 illustrates the orientation of the polarization
measurements achieved in Planck. It shows, for an actual Planck detector,
the maps of and respectively, where
is the direction of the polarizer with respect to the local Galactic
meridian. These quantities contribute to the spin 2 term
defined in Eq. (59). The large amplitude of these two maps is
consistent with the fact that for a given detector, the orientation of the
polarization measurements is mostly and , as expected when
detectors move on almost great circles with very little precession. Another
striking feature is the relative smoothness of the maps, which translate into
the power spectrum of
peaking at low values.
Figure 2 shows the same information for an hypothetical
LiteBIRD
3.2 Arbitrary beams, smooth scanning case
If one assumes that and vary slowly across the sky, as we just saw in the case of Planck and LiteBIRD  and probably a wider class of orbital and suborbital missions  then is dominated by low values and one expects because of the triangle relation imposed by the 3J symbols (see Appendix C). If one further assumes and to vary slowly in , then Eqs. (84) and (88) can be used to impose in Eq. (38) and provide
(39) 
with
(40a)  
(40b)  
(40c) 
As derived in Appendix E.1, Eq. (39) reduces to a mixing equation relating the observed crosspower spectra to the true ones:
(41) 
with .
In the smooth scanning case representative of past and forthcoming satellite missions, the effect of observing the sky with nonideal beams is therefore to couple the temperature and polarization power spectra at the same multipole through an extended beam window matrix , as illustrated on Fig. 3.
3.3 Arbitrary scanning, circular identical beams
If the scanning beams are now assumed to all be circular and identical, the measured will not depend on the details of the scanning strategy, orientation of the detectors, or relative weights of the detectors. We are indeed exactly in the ideal hypotheses of the map making formalism (Eq. 29) and get the well known and simple result that the effect of the beam can be factored out.
If one considers detectors with identical circular copolarized beams, and whose actual polarization efficiency was used during the map making: , such that
(42) 
then Eqs. (39) and (40b) feature terms like , which when written in a matrix form, verify the equality
(43) 
according to Eq. (79). The measured power spectra are then
(44) 
and for the Gaussian circular beam introduced in Eq. (23a). Obviously, these very simple results assume that the whole sky is observed. If not, the cutsky induced and coupling effects mentioned at the end of Section 2.4 have to be accounted for, as described, for example, in Chon et al. (2004), Mitra et al. (2009), Grain et al. (2009), and references therein.
3.4 Arbitrary beams, ideal scanning
Let us now consider the case of an ideal scanning of the sky, for which in any pixel , the number of valid (unflagged) samples is the same for all detectors , and each detector covers uniformly all possible orientations within that pixel along the duration of the mission. This constitutes the ideal limit aimed at by the scanning strategy illustrated in Fig. 2. The assumption of smooth scanning is then perfectly valid, and details of the calculations can be found in Appendix E.2. We find for instance that the matrix describing how the measured temperature and polarization power spectra are affected by the input spectrum reads
(45) 
with the normalization factors
(46) 
This confirms that in this ideal case, as expected and discussed previously (e.g., Wallis et al. 2014, and references therein), the leakage from temperature to polarization (Eq. 45) is driven by the beam ellipticity ( terms) which has the same spin as polarization. One also sees that the contamination of the and spectra by are swapped (e.g., ) when the beams are rotated with respect to the polarimeter direction by (), as shown in Shimon et al. (2008).
3.5 Finite pixel size and subpixel effects
As shown in Planck 2013VII (2014), in the case of temperature fluctuations, the effect of the finite pixel size is twofold. First, in each pixel, the distance between the nominal pixel center and the center of mass of the observations couples to the local gradient of the Stokes parameters to induce noise terms. Second, there is a smearing effect due to the integration of the signal over the surface of the pixel. Equation (41) then becomes
(47) 
with , and the squared displacement averaged over the hits in the pixels and over the set of considered pixels. As shown in Appendix F, the additive noise term, sourced by the temperature gradient within the pixel, affects both temperature and polarization measurements, with and , while the other spectra are much less affected, that is, . The sign of this noise term is arbitrary and can be negative when crosscorrelating maps with a different sampling of the pixels.
4 Numerical implementation
Numerical implementations of this formalism are performed in three steps, assuming that the individual beam is already computed for and :

For each involved detector , and for , one computes the th complex moment of its direction of polarization in pixel : defined in Eq. (59). Since this requires processing the whole scanning data stream, this step can be time consuming. However it has to be computed only once for all cases, independently of the choices made elsewhere on the beam models, calibrations, noise weighting, and other factors. As we shall see below, it may not even be necessary to compute it, or store it, for every sky pixel.

The computed above are weighted with the assumed inverse noise variance weights and polar efficiencies to build the hit matrix in each pixel, which is then inverted to compute the , defined in Eq. (62). Those are then multiplied together to build the scanning information matrix using its pixel space definition (Eq. 40b). The resulting complex matrix contains elements, where and are the number of detectors in each of the two detector assemblies whose crossspectra are considered. This step can be parallelized to a large extent, and can be dramatically sped up by building this matrix out of a representative subset of pixels. In our comparison to simulations, described in Section 5, and performed on HEALPix map with and pixels, we checked that using only pixels evenly spread on the sky gave final results almost identical to those of the full calculations.