A Projection of maps on spherical harmonics

QuickPol: Fast calculation of effective beam matrices for CMB polarization

Key Words.:
Cosmology, polarization, systematic effects

Current 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 half-wave 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 electronics-related imperfections on the measured CMB temperature and polarization angular power spectra and their statistical isotropy. Systematic effects such as beam non-circularity, 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 WMAP2 (Smith et al. 2007; Hinshaw et al. 2007; Page et al. 2007) or Planck3 (Planck 2013-VII 2014; Planck 2013-XVII 2014; Planck 2015-XI 2016) satellite missions.
At the same time, several deconvolution algorithms and codes have been proposed to clean up the CMB maps from such beam-related effects prior to the computation of the power spectra, like PreBeam (Armitage-Caplan & 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 2013-IV 2014; Planck 2013-VII 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 Monte-Carlo 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 2013-VII (2014). It allows a quick and accurate computation of the leakage and cross-talk 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 time-independent 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 Monte-Carlo 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 half-wave 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 electro-magnetic radiation in the optical system. Let us consider a quasi monochromatic4 radiation propagating along the axis, and hitting the optical system at a position . The incoming electric field will be turned into , where is the 2x2 complex Jones matrix of the system.

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




and the sign representing the adjoint operation, which for a real rotation matrix, simply amounts to the matrix transpose. The measured signal is




We now introduce the Stokes parameters of the input signal (dropping the dependence on )


and of the (un-rotated) instrument response


to obtain


With the rotated instrument response:


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 :




for . The Stokes parameters of the instrument are then for .

If the beam is assumed to be perfectly co-polarized, that is, it does not alter at all the polarization of the incoming radiation, with and , then , , and , , , Eqs. (8, 9) become




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 co-polarized 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 co-polarized.

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 2013-VII (2014) and Planck 2015-VII (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):


although one usually prefers the scalar and fixed parity and components


such that for . In other terms




The sign convention used in Eq. (16) is consistent with Zaldarriaga & Seljak (1997) and the HEALPix5 library (Górski et al. 2005).

The response of a beam centered on the North pole can also be decomposed in SH coefficients


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


The elements of Wigner rotation matrices are related to the SH via (Challinor et al. 2000)


with .

If the beam is assumed to be co-polarized 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


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


while we show in Appendix G that Eq. (24) is true for arbitrarily shaped co-polarized 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 so-called Pxx coordinates in Planck parlance), reads


so that


We introduced to distinguish it from the value used in the map-making, 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


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 2013-VIII 2014), and we assumed the circular polarization to vanish. With the definitions introduced in Section 2.2, this becomes


The map-making formalism is set ignoring the beam effects, assuming a perfectly co-polarized detector and an instrumental noise (Tristram et al. 2011, and references therein), so that, for a detector , Eq. (12) becomes


where the leading prefactors are here again absorbed in the gain calibration. Let us rewrite it as


with (Shimon et al. 2008)


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


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 map-making process; Eq. (33) then becomes


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 so-called sub-pixel effects will be considered in Section 3.5.

2.4 Measured power spectra

To compute the cross-power spectrum of any two spin and maps, we first project each polarized component of on the appropriate spin weighted sets of spherical harmonics,


and average these terms according to


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 noise-weighted 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 non-circular beams of the pseudo-power spectra measured on a masked or weighted map (Hivon et al. 2002; Hansen & Górski 2003), and extends to polarization the Quickbeam non-circular beam formalism used in the data analysis conducted by Planck 2013-VII (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 full-fledged Planck-HFI 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.

Figure 1: Orientation of polarization measurements in Planck. The two left panels show, for an actual Planck detector, the maps of and respectively, where is the direction of the polarizer with respect to the local Galactic meridian, which contributes to the spin 2 term defined in Eq. (59). The right panel shows the power spectrum of , multiplied by .
Figure 2: Same as Fig. 1 for an hypothetical detector of a LiteBIRD-like mission, except for the right panel plot which has a different -range.
Figure 3: Effective beam window matrix introduced in Eq. (41) and detailed in Eq. (102a), for the cross-spectra of two simulated Planck maps discussed in Section 5. Left panel: raw elements of , showing for each how the measured map angular power spectrum is impacted by the input spectrum, because of the observation of the sky with the beams. Right panel: blown-up ratio of the non-diagonal elements to the diagonal ones:

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 LiteBIRD6 like detector (but without half-wave plate modulation) in which we assumed the detector to cover a circle of 45 in radius in one minute, with its spin axis precessing with a period of four days at 50 from the anti-sun direction. As expected for such a scanning strategy, the values of are pretty uniformly distributed over the range , which translates into a low amplitude of the and maps. Even if those maps do not look as smooth as those of Planck, their power spectra peak at fairly low multipole values.

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 sub-orbital 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




As derived in Appendix E.1, Eq. (39) reduces to a mixing equation relating the observed cross-power spectra to the true ones:


with .

In the smooth scanning case representative of past and forthcoming satellite missions, the effect of observing the sky with non-ideal 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


then Eqs. (39) and (40b) feature terms like , which when written in a matrix form, verify the equality


according to Eq. (79). The measured power spectra are then


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 cut-sky 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


with the normalization factors


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).

Figure 4: Computer simulated beam maps (, , and clockwise from top-left) for two of the Planck-HFI detectors (100-1a and 217-5a) used in the validation of QuickPol. Each panel is 1x1 in size, and the units are arbitrary.

3.5 Finite pixel size and sub-pixel effects

As shown in Planck 2013-VII (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


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 cross-correlating 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 :

  1. 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.

  2. 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 cross-spectra 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.

  3. Finally, using Eqs. (95-41) we note that , so that, for instance, for a given , the 3x3 matrix is computed by replacing in Eq. (95) its central term with its partial derivative, such as