© ESO, 2012

1. Introduction

Removing systematic effects plays an important role in the data analysis of modern high-sensitivity cosmic microwave background (CMB) experiments, such as the WMAP and Planck missions (Jarosik et al. 2011; Tauber et al. 2010). In this work we concentrate on one source of systematic effects, the beam asymmetry. We present an efficient beam deconvolution code, called artDeco. It can be applied to any absolute CMB experiment with sufficient sky coverage.

The required input consists of the time-ordered data (TOD), the corresponding detector pointings, and known beam shapes. The primary output of the code consists of the harmonic a_Tlm, a_Elm, and a_Blm coefficients of the sky. From these one can derive temperature and Q and U polarisation maps.

This is the usual deconvolution map-making problem (see, e.g., Armitage & Wandelt 2004; Armitage-Caplan & Wandelt 2009; Harrison et al. 2011). Our method differs from the earlier works by an efficient reformulation of the deconvolution problem through heavy use of Wigner functions. The formulation is general and makes no assumptions on the scanning pattern. Similar ideas have been studied by Prézeau, independently of us (priv. comm.).

In the CMB literature the concept of map-making often refers to a procedure that involves removal of correlated noise. Methods for doing this have been treated in several papers (Keihänen et al. 2005, 2010; Poutanen et al. 2006; Ashdown et al. 2007a,b; Kurki-Suonio et al. 2009; Sutton et al. 2009). The methods discussed in these works construct pixelised temperature and Q, U polarisation maps, but do not attempt to correct for an asymmetric beam shape.

Our method steps in at the point where the correlated noise has already been cleaned from the data. The input TOD is assumed to include signal and uncorrelated noise only. We are assuming here that the correlated noise component can be removed from the data to a sufficient degree prior to the beam deconvolution step. The noise removal step itself is affected by beam asymmetries because signal differences caused by beam mismatch can be falsely interpreted as noise. Fortunately, in destriping methods it is possible to strongly reduce this effect by masking out the regions with the strongest foreground contamination that generate most of the effect. Also, destriping methods produce as natural output the noise-cleaned TOD, which is the required input in deconvolution map-making. Accordingly, destriping methods are the natural counterparts of our deconvolution method.

Since the aim of this paper is to present a new deconvolution method, rather than assess the performance of any particular instrument, we tested the code with a somewhat idealised simulation. In particular, we ignored other systematic effects such as frequency response. On the other hand, we chose strongly asymmetric beam shapes to show the beam effects more clearly.

2. Beam deconvolution

2.1. Constructing the linear system

In this section we present the deconvolution problem using a formalism based on Wigner functions. We write the solution in a form that allows us to solve the problem numerically in an efficient way. Some details are left for Appendix A.

We start by writing the signal t_j received by a detector at a given time τ_j as ${\mathit{t}}_{\mathit{j}}\mathrm{=}\sum _{\mathit{slmk}}{\mathit{a}}_{\mathit{slm}}{\mathit{b}}_{\mathit{slk}}^{\mathrm{\ast }}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}\mathrm{+}{\mathit{n}}_{\mathit{j}}\mathrm{.}$(1)The a_slm (s = 0, ± 2) are harmonic coefficients that represent temperature and polarisation of the sky signal, b_slk are corresponding beam coefficients, n_j represents noise, and ${\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}$ is a Wigner function that defines a rotation of the beam from a fiducial orientation at the north pole to its actual position and orientation. For exact definitions, see Appendix A.

We have adopted a compact notation where ω_j denotes a combination of three angles  { φ_j,θ_j,ψ_j } , which define the detector’s orientation. Angles θ_j and φ_j define a point on the celestial sphere, and ψ_j defines the beam orientation. The angles can be interpreted as Euler rotation angles.

We perform a simple linear least-squares fit of a_slm to the data up to some multipole l_max.

In general, a linear least-squares fit leads to an equation of the form $\mathit{x}\mathrm{=}{\mathrm{(}{\mathit{A}}^{}{\mathit{C}}_{\mathrm{n}}^{-1}\mathit{A}\mathrm{+}{{\mathit{C}}_{\mathit{s}}^{-1}}^{\mathrm{)}}}^{-1}{\mathit{A}}^{}{\mathit{C}}_{\mathrm{n}}^{-1}\mathit{y,}$(2)where x is the vector of unknowns, y is the data vector, C_n is noise covariance, and C_s represents a priori information on the covariance of x (Wiener 1949). Matrix A relates the unknowns to the data through $\mathit{y}\mathrm{=}\mathit{Ax}\mathrm{+}\mathit{n}\mathit{.}$(3)In the present case matrix A is ${\mathit{A}}_{\mathit{slm}}^{\mathit{j}}\mathrm{=}\sum _{\mathit{k}}{\mathit{b}}_{\mathit{slk}}^{\mathrm{\ast }}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}$(4)and x represents a_slm.

We made two simplifying assumptions, which are not of fundamental nature, however. Firstly, we assume the noise is white with constant variance, and assign equal weights to all samples j. Secondly, we use no a priori information on the sky signal, i.e. the harmonic coefficients are assumed to be uncorrelated and can host a signal of infinite variance. Under these assumptions the covariance matrices drop out of the equation (C_n = const., ${\mathit{C}}_{\mathit{s}}^{-1}\mathrm{=}\mathrm{0}$). The least-squares solution for a_slm is obtained by solving the linear system of equations $\sum _{{\mathit{s}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\sum _{\mathit{j}}{\mathit{A}}_{\mathit{slm}}^{\mathit{j}\mathrm{\ast }}{\mathit{A}}_{{\mathit{s}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{\mathit{j}}{\mathit{a}}_{{\mathit{s}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{=}\sum _{\mathit{j}}{\mathit{A}}_{\mathit{slm}}^{\mathit{j}\mathrm{\ast }}{\mathit{t}}_{\mathit{j}}\mathit{.}$(5)In the following we process Eq. (5) further to bring it into a form where it can be solved in an efficient way. We treat the right- and left-hand sides separately.

2.1.1. Right-hand side

The right-hand side yields $\sum _{\mathit{j}}{\mathit{A}}_{\mathit{slm}}^{\mathit{j}\mathrm{\ast }}{\mathit{t}}_{\mathit{j}}\mathrm{=}\sum _{\mathit{kj}}{\mathit{b}}_{\mathit{slk}}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}{\mathit{t}}_{\mathit{j}}\mathrm{=}\sum _{\mathit{k}}{\mathit{b}}_{\mathit{slk}}{\mathit{T}}_{\mathit{mk}}^{\mathit{l}}\mathit{.}$(6)Here we define ${\mathit{T}}_{\mathit{mk}}^{\mathit{l}}\mathrm{=}\sum _{\mathit{j}}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}{\mathit{t}}_{\mathit{j}}\mathit{.}$(7)The triplets ω_j =  { φ_j,θ_j,ψ_j }  define the detector orientation and are provided as input to the code. To numerically evaluate the transform, we replace the sum over TOD samples by a sum over a three-dimensional (3D) grid as follows. We divide the space spanned by all possible triplets ω =  { φ,θ,ψ }  into bins of equal volume and denote by t(ω) the sum of values t_j falling into one bin. In θ,φ we use HEALPix¹ pixelisation (Górski et al. 2005), and in ψ we divide the space uniformly. Most elements of t(ω) are zero, since a particular HEALPix pixel is typically observed only in a few orientations of the beam.

We call quantity t(ω) the 3D signal map, and quantity ${\mathit{T}}_{\mathit{mk}}^{\mathit{l}}\mathrm{=}\sum _{\mathit{\omega }}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{t}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$(8)the Wigner transform of the signal map.

Here  ∑ _ω denotes a sum over all bins. Since most bins are empty, the sum needs to be carried out only over a subset of all ω values. Expansion (8)can be evaluated numerically quite efficiently using fast Fourier transforms (FFTs) and fast algorithms for the generation of Wigner functions.

Quantities t(ω) and ${\mathit{T}}_{\mathit{mk}}^{\mathit{l}}$ carry no memory of the order in which the sky was scanned. They just record which pixels of the sky were observed in which orientation of the beam.

2.1.2. Left-hand side

Consider then the matrix on the left-hand side of Eq. (5), $\sum _{\mathit{j}}{\mathit{A}}_{\mathit{slm}}^{\mathit{j}\mathrm{\ast }}{\mathit{A}}_{{\mathit{s}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{\mathit{j}}\mathrm{=}\sum _{\mathit{jk}{\mathit{k}}^{\mathrm{\prime }}}{\mathit{b}}_{\mathit{slk}}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}{\mathit{b}}_{{\mathit{s}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathit{.}$(9)The multiplier matrix is too large to be inverted. Instead, we use the conjugate gradient (CG) iteration method to solve Eq. (5).

We pixelise the 3D pointing space in the same way we did with the signal, and denote by n(ω) the number of hits into bin ω. In line with the definition of the 3D signal map, we call n(ω) the 3D hit map.

We replace the sum over samples j by a sum over bins and write $\sum _{\mathit{j}}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\omega }}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{n}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{.}$(10)As the next step, we write the Wigner functions as ${\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta ,\phi ,\psi }\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{m\phi }}{\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k\psi }}\mathit{,}$(11)where ${\mathit{d}}_{\mathit{mk}}^{\mathit{l}}$ are the reduced Wigner functions. Inserting this into (10), we obtain $\begin{array}{ccc}& & \sum _{\mathit{j}}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}\mathrm{=}\sum _{\mathit{\theta }}{\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{d}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\\ & & \mathrm{\times }\sum _{\mathit{\phi ,\psi }}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{m}\mathrm{-}{\mathit{m}}^{\mathrm{\prime }}\mathrm{)}\mathit{\phi }}{\mathrm{e}}^{\mathrm{i}\mathrm{(}\mathit{k}\mathrm{-}{\mathit{k}}^{\mathrm{\prime }}\mathrm{)}\mathit{\psi }}\mathit{n}\mathrm{\left(}\mathit{\theta ,\phi ,\psi }\mathrm{\right)}\mathrm{=}\sum _{\mathit{\theta }}{\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{d}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{N}}_{\mathit{m}\mathrm{-}{\mathit{m}}^{\mathrm{\prime }}\mathit{,k}\mathrm{-}{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathit{,}\end{array}$(12)where we have defined ${\mathit{N}}_{\mathit{mk}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\sum _{\mathit{\phi ,\psi }}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{m\phi }}\mathit{n}\mathrm{\left(}\mathit{\theta ,\phi ,\psi }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k\psi }}\mathit{.}$(13)The sums over θ,φ,ψ are carried out over 3D bins. Substituting this back into (9)and rearranging, we bring the left-hand side of the deconvolution equation into the form $\begin{array}{ccc}& & \\ & & \sum _{\mathit{k}}{\mathit{b}}_{\mathit{slk}}\sum _{\mathit{\theta }}{\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\sum _{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}{\mathit{N}}_{\mathit{m}\mathrm{-}{\mathit{m}}^{\mathrm{\prime }}\mathit{,k}\mathrm{-}{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\sum _{{\mathit{l}}^{\mathrm{\prime }}{\mathit{s}}^{\mathrm{\prime }}}{\mathit{d}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathit{b}}_{{\mathit{s}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{\mathrm{\ast }}{\mathit{a}}_{{\mathit{s}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\mathit{.}\end{array}$(14)This can be evaluated much more quickly than the original formula (9), since the summation over the long TOD vector has been replaced by an operation on the more compact N object. Also, the convolution operations in indices k and m can be carried out efficiently with FFTs.

The terms are written in the order in which they are applied in the code.

2.2. Formulation through 3j coefficients and preconditioner

Equation (9)can be written in a yet more compact form.

The product of two Wigner functions can be expressed in terms of 3j symbols (Appendix A.3) as $\begin{array}{ccc}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}& & \mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{m}}^{\mathrm{\prime }}\mathrm{+}{\mathit{k}}^{\mathrm{\prime }}}\sum _{{\mathit{l}}_{\mathrm{2}}}\mathrm{(}\mathrm{2}{\mathit{l}}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{\left)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}\mathrm{-}\mathit{m,}{\mathit{k}}^{\mathrm{\prime }}\mathrm{-}\mathit{k}}^{{\mathit{l}}_{\mathrm{2}}}\mathrm{\right(}\mathit{\omega }\mathrm{)}\\ & & \mathrm{\times }\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\mathit{.}\end{array}$(15)Substituting this into (10)and that into (9), we arrive at the formula $\begin{array}{ccc}\sum _{\mathit{j}}{\mathit{A}}_{\mathit{slm}}^{\mathit{j}\mathrm{\ast }}{\mathit{A}}_{{\mathit{s}}^{\mathrm{\prime }}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}^{\mathit{j}}& & \\ & & \mathrm{\times }\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\mathit{.}\end{array}$(16)Indices l and l′ run from 0 to l_max. The sum over l₂ covers the values for which |l − l′| ≤ l₂ ≤ l + l′. We have defined the Wigner transform of the hit map as ${\mathit{W}}_{\mathit{mk}}^{\mathit{l}}\mathrm{=}\sum _{\mathit{\omega }}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{n}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\mathit{.}$(17)We use the formula (16)to explicitly construct the diagonal of the matrix, which serves as a preconditioner to speed up the CG iteration.

The formulation presented here is general, and makes no assumptions on the scanning pattern or on the sky coverage.

3. Implementation

3.1. User interface

The artDeco code takes two lists of input data objects containing the beam coefficients b_slk and the so-called 3D maps for the respective detector for all detectors that are to be processed simultaneously. The beam coefficients can be read from FITS files in the format that is also being used by the HEALPix package for storing spherical harmonic coefficients.

The 3D maps contain the quantities n(ω) and t(ω) for a given detector; they are created from time streams of detector pointings and sky signal, respectively, by another stand-alone module called todtobin. This module achieves the angular discretisation of the detector pointings by sorting them into pixels of a HEALPix map with a user-defined N_side parameter, which are in turn subdivided into n_ψ bins for the discretisation in ψ. The choice of N_side and n_ψ has direct consequences on the run-time of the deconvolver and the accuracy of its results: if high values are chosen, the 3D map will become very large, and the deconvolver will require more memory and CPU time. On the other hand, too low values will result in noticeable discretisation errors and inaccurate output a_Xlm. As a guideline, N_side should be chosen in such a way that the smallest features of all beams are well resolved at the angular resolution $\mathit{\delta \alpha }\mathrm{\approx }\mathrm{3500}\prime \mathit{/}{\mathit{N}}_{\mathrm{side}}$(18)associated with this N_side. In a similar fashion, n_ψ should be large enough to resolve azimuthal features of the beams (assuming the beams are centred on a pole). For completely axisymmetric beams, n_ψ = 1 is sufficient. Generally, n_ψ = 10k_max (where k_max is the highest azimuthal multipole used to describe the beams) is a good choice.

The computational complexity of todtobin is approximately , where n_samp is the total number of TOD samples, since its runtime is dominated by the sorting of the input TOD according to HEALPix pixel number and ψ bin. Consequently, N_side and n_ψ only have a significant effect on the size of the output 3D maps. Typical execution times for todtobin with several billion samples (as used for the experiments shown in this paper) are of the order of a few minutes.

ArtDeco takes several user-supplied parameters that influence its behaviour. Most important among these are l_max, which specifies the maximum multipole up to which deconvolution is carried out, and k_max ≤ l_max, which limits the maximum azimuthal multipole used for the beams. The explicit introduction of k_max is advantageous since most beams can be expressed accurately even when using a k_max ≪ l_max (even k_max = 0 if the beam has rotational symmetry and is unpolarised), and the reduction in k_max reduces CPU and memory requirements dramatically.

ArtDeco’s output consists of a set of a_lm in HEALPix-compatible format. Depending on whether the user requested polarisation, the file contains only a_Tlm or additional a_Elm and a_Blm coefficients.

3.2. Algorithm overview

The deconvolution algorithm has been implemented in C++, making use of the MPI library for parallelisation. The quantities that depend on θ, such as d and N, are distributed to MPI tasks according to the θ angle. This allows a fairly good load balance. Summation over θ, when needed, is done efficiently through the collective MPI_Reduce() routine. The a_slm coefficients are distributed according to index m. The beam coefficients b_slk are fully present on all tasks.

The procedure of evaluating Eq. (14)can be summarised as follows. First we multiply the input a by the beam coefficients b. Then we perform a loop over index m′, inside which we broadcast the a_s′l′m′ coefficients for the current m′ to all processes, multiply by b, construct the reduced Wigner matrices for the local values of θ, and accumulate over l′. The multiplication by matrix N is a convolution operation in indices m,k and can be performed efficiently using 2D FFTs. Each MPI task performs the multiplication by N independently for the local θ. After that, we make another loop over index m, construct again the Wigner matrices, sum over θ, collect the result according to m to the owner process, and multiply by the beam coefficients.

The solution is assumed to have converged sufficiently as soon as the residual of the current a_slm estimate is less than 10^-12 times the residual of a vector containing all zeroes.

To generate the 3j symbols needed to construct the preconditioner, we use a C implementation of the recursive algorithm described by Schulten & Gordon (1975).

For a more detailed technical description of the algorithm, see Appendix B. CPU and memory requirements of various runs are presented and discussed in Appendix C.

4. Simulation setup

The deconvolution algorithm was tested using mock data produced by the Planck simulation package (Level-S, Reinecke et al. 2006).

4.1. Mission parameters

We assumed a scanning strategy that caused the telescope’s spin axis to perform two cycloidal excursions from the ecliptic plane per year with an amplitude of 7.5°; in azimuthal direction, the spin axis always pointed towards the Sun.

Time-ordered data were generated for a time span of 10 000 h, corresponding to slightly more than two full sky surveys; their length is of the order of 10⁹ samples.

4.2. Detector geometry

For the first set of experiments, data were simulated for two horns, i.e. two pairs of detectors, which were modelled roughly following Planck’s 30 GHz channel. In particular, the detectors had an average FWHM of 32$\stackrel{\mathrm{\prime }}{.}$ 2, the polarisation directions in each detector pair were rotated by about (but not exactly) 90° against each other, and the polarisation directions of the two horns were in turn shifted by 45°. Both horns followed the same scan path on the sky, although with a little time delay.

The beam shapes were assumed to be elliptical, with ellipticities of 1.7 and 1.8 for each detector in a pair, respectively. These values are considerably higher than in a typical experiment and were chosen to demonstrate the performance of the method even under unfavourable conditions.

In order to show the algorithm’s performance at higher resolutions, another data set was produced using two horns sensitive at 70 GHz and with an average FWHM of 13′.

Table 1 gives an overview over the detector properties.

4.3. Input data

To generate a theoretical power spectrum for the CMB emission, we used the CAMB² code with the cosmological parameters Ω_b = 0.045, Ω_cdm = 0.255 and Ω_λ = 0.7. Using the HEALPix tool syn_alm_cxx, we created a random, unconstrained realisation of polarised a_lm from this power spectrum.

Emission maps of foreground components were generated using the Planck pre-launch Sky model³ (PSM, Delabrouille et al. 2012), version 1.6.6. We chose to include the contributions from galactic infrared emission (galactic synchrotron radiation, free-free emission, as well as thermal and spinning dust).

The signal from strong point sources was taken into account as well; to achieve this, a PSM-generated catalogue of point sources was directly expanded into spherical harmonic coefficients, using a sufficiently high cutoff value l_max that the beam response was negligible at this scale.

Polarised emission was disabled for all contributions except for the CMB, which makes it very easy to detect leakage of the strong temperature signal in the galactic plane into the polarisation signal.

We neglected bandpass effects and assumed an identical delta frequency response for all detectors.

The power spectra of the signal components at 30 GHz are plotted in Fig. 1.

4.4. Noise

In some of our simulations, detector noise was added to the idealised sky signal. Since the deconvolution method requires correlated noise to be removed from the TOD beforehand, we only added Gaussian white noise with an appropriately chosen σ to the simulated samples (see Table 1).

4.5. 3D map generation

The mock TOD were processed with the todtobin tool to generate the 3D input maps required by the deconvolver. We set N_side = 1024 for the 30 GHz case, and N_side = 2048 for 70 GHz; n_ψ was always set to 256.

5. Results

5.1. Constructing a sky map

The primary output of the deconvolution code consists of the a_Xlm (X = T, E, B) coefficients of the sky, up to l_max. From these one can construct the usual temperature and Q, U polarisation maps through harmonic expansion of the coefficients. We did this with the alm2map_cxx tool of the HEALPix package.

For comparison, we also produced binned maps directly from the TOD. A binned map is constructed by assigning each TOD sample completely to the pixel containing the centre of the beam. The binning operation can be formally written as $\mathit{m}\mathrm{=}{\mathrm{(}{\mathit{P}}^{\mathit{T}}{\mathit{C}}_{\mathrm{n}}^{-1}{\mathit{P}}^{\mathrm{)}}}^{-1}{\mathit{P}}^{\mathit{T}}{\mathit{C}}_{\mathrm{n}}^{-1}\mathit{y}\mathit{.}$(19)Here P is a pointing matrix, C_n is white noise covariance, and m is a (I, Q, U) map triplet.

Binned maps are distorted by beam effects if the beam is asymmetric. In particular, the polarisation component of a binned map is contaminated by temperature leakage through beam shape mismatch. We aim at demonstrating that deconvolution can reduce these undesired effects.

Since our simulation includes point sources, the input a_Xlm coefficients are not band-limited. However, we can only solve the coefficients up to some limited l_max. This leads to ringing artefacts if a map is constructed directly from the raw a_Xlm coefficients. Ringing is particularly prominent around point sources and other sharp features.

Ringing is caused by the sharp cut-off of the spectrum at l_max, and it occurs even if the deconvolved coefficients are fully correct up to l_max. For this reason it is necessary to smooth the coefficients with a symmetric Gaussian beam prior to constructing a map. The required amount of smoothing depends on the spectrum of the underlying sky, and on l_max. The highest l_max that one can solve in turn depends on the beam width. In other words, it is not possible to recover structures significantly smaller than the beam itself. This gives the following rule of thumb: the required smoothing width is given by the smallest features of the original beam.

In Fig. 2 we demonstrate the ringing around a point source. We show an image of a point source, constructed from the input a_Xlm, but including only coefficients up to l_max = 800. Ringing artefacts are evident. On the right we show the same region of the sky after the map is smoothed with a FWHM = 32′ Gaussian symmetric beam. Smoothing removes the ringing, but naturally also smoothes out the smallest structures in the map.

5.2. Signal-only simulations with elliptic beam

We analyse first the noise-free 30 GHz simulation in detail. A detailed description of the simulation was given in Sect. 4. We ran the code with various combinations of parameters l_max, k_max, and analysed the results both in pixel space and in harmonic space. We included data from all four detectors and considered both temperature and polarisation.

5.2.1. Temperature map

In Fig. 3 we show a full-sky Mollweide projection of the binned temperature map together with that of the deconvolved map. We used deconvolution parameters l_max = 800, k_max = 6.

We chose the value FWHM = 32′ as the baseline smoothing width for the 30 GHz maps. This corresponds to the mean width of the input beam. In some cases it was necessary to apply even stronger smoothing.

The difference between the binned and the deconvolved map does not appear dramatic in the full-sky map image, but becomes evident when we zoom into a point source. Figure 4 shows a randomly selected strong point source below the galactic plane. While the source in the binned map is clearly elongated due to beam shape, it appears circular in the deconvolved map.

The deconvolved image can be compared to the right panel of Fig. 2, where we have shown the same patch of the sky constructed from the input a_Xlm coefficients. The maps are nearly identical.

5.2.2. Polarisation maps

The beam effects show up more dramatically in the polarisation maps (Fig. 5). The temperature signals recorded by two detectors differ due to the different beam shapes. The binning procedure interprets this erroneously as polarisation signal. This leads to a leakage of temperature signal into the polarisation maps. Since our simulation did not include polarised foregrounds, the whole galactic structure seen in the binned map is temperature leakage through beam mismatch.

We show deconvolved maps with two levels of smoothing below the binned map. Smoothing with a FWHM = 32′ beam is not sufficient to fully remove the galactic residual, but if we smooth the map further to FWHM = 45′, the residual disappears almost completely.

In Fig. 6 we show a zoom into the Q polarisation map, just below the galactic plane. In the binned map we see a strong galactic residual, and typical “four-leaf clover” structures that arise from temperature leakage at the location of a point source. We show the deconvolved map again for two levels of smoothing: FWHM = 32′ and FWHM = 45′. The additional smoothing removes the clover structure nearly completely. To show that this is not only an effect of more aggressive smoothing, we smoothed the binned map further with a Gaussian beam of width FWHM = 32′, which, combined with the 32′ of the original beam, gives a total smoothing of 45′. This map is shown below the unsmoothed binned map. The additional smoothing has little effect on the galactic residual.

5.2.3. Harmonic space

We proceed to analyse the deconvolution results in harmonic space. We computed the TT, EE, and TE spectra of the deconvolved a_Xlm coefficients and compared them to the corresponding input spectra. The spectrum is constructed as

${\mathit{C}}_{\mathit{l}}^{\mathit{XY}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{l}\mathrm{+}\mathrm{1}}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}{\mathit{a}}_{\mathit{Xlm}}^{\mathrm{\ast }}{\mathit{a}}_{\mathit{Ylm}}\mathrm{,}$where X and Y stand for T and E.

We ran the code for all combinations of k_max = 4, 6 and l_max = 400, 500, 600, 700, 800. Results for different combinations are shown in Fig. 7. We applied no smoothing to the coefficients.

We observe that the optimal value for k_max is coupled to l_max. Results for k_max = 4 and k_max = 6 begin to diverge only above l = 400. The results for different l_max with k_max = 6 are practically on top of each other, except for the highest multipoles near l_max, which are clearly distorted. When constructing a map, we cut out the last 20 multipoles before making the harmonic expansion.

For comparison we also show the spectrum obtained from a harmonic expansion of the binned map, which we computed using the anafast tool of the HEALPix package. This decreases rapidly as a function of l as a result of beam smoothing. To correct for the smoothing, we divided the spectrum by the window function of a symmetric Gaussian beam with FWHM = 32′. This corrects for the mean smoothing and allows one to recover the correct spectrum up to l_max = 300. Above that, the result begins to deviate from the input spectrum due to beam ellipticity, which cannot be corrected for by application of a window function only.

Deconvolution recovers the TT input spectrum well up to l = 800. At higher multipoles there is little information left in the signal, and the convergence properties of the code become poor. The weaker EE and TE spectra are recovered up to l = 400. Owing to temperature leakage, the binned map only gives the correct EE spectrum for the ten first multipoles in EE, and the TE spectrum is useless.

We did one run with parameters l_max = 600, k_max = 8 (not shown), to see if it would help to recover the EE spectrum to even higher multipoles. The results, however, did not differ significantly from the l_max = 600, k_max = 6 case.

5.2.4. Reconstruction error of harmonic coefficients

A comparison between the reconstructed spectrum and the input spectrum is not a fully reliable test of the correctness of the deconvolution result. It is possible to construct a set of harmonic coefficients that provides the correct spectrum, although the individual coefficients are incorrect.

To assess directly the accuracy of the a_Xlm coefficients, we constructed a quantity we call the error spectrum. It is computed as follows. We take the difference between the a_Xlm coefficients we obtain from deconvolution and the corresponding input coefficients ${\mathit{a}}_{\mathit{Xlm}}^{\mathrm{0}}$. As before, X stands for T or E. We compute the spectrum of this difference as

${\mathit{C}}_{\mathit{Xl}}^{\mathit{err}}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{l}\mathrm{+}\mathrm{1}}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}{\mathrm{(}{\mathit{a}}_{\mathit{Xlm}}\mathrm{-}{{\mathit{a}}_{\mathit{Xlm}}^{\mathrm{0}}}^{\mathrm{)}}}^{\mathrm{\ast }}\mathrm{(}{\mathit{a}}_{\mathit{Xlm}}\mathrm{-}{{\mathit{a}}_{\mathit{Xlm}}^{\mathrm{0}}}^{\mathrm{)}}\mathrm{.}$The resulting spectrum is a measure of the error in the reconstruction of the a_Xlm coefficients, as a function of scale.

We show the error spectra of the 30 GHz runs for T and E coefficients in Fig. 8. For comparison, we show the input TT or EE spectrum in the same figure. The a_Xlm coefficients are recovered with good accuracy when the error spectrum is well below the input spectrum. The meaning of the blue line is the same as in Fig. 7, that is, we compute the harmonic expansion of the binned map and divide them by the window function of a symmetric Gaussian beam (FWHM = 32′).

We made another comparison with an ideal window function, which we computed as the ratio of the uncorrected spectrum of the binned map, and the input spectrum (green and black lines in the top panel of Fig. 7). The error spectrum for this is shown in light blue in Fig. 8 (TT spectrum only). By construction, applying the ideal window function to the spectrum of the binned map returns exactly the correct input spectrum. The error spectrum, however, is not improved with respect to the result obtained with Gaussian window function, which shows that the failure of the binning procedure to recover the correct spectrum above l = 300 in Fig. 7 was not just a result of a poorly chosen window function.

The error spectrum reveals that the deconvolved a_Xlm coefficients are more accurate than those obtained by harmonic expansion of the binned map already at the lowest multipoles.

5.3. Simulations with noise

The first simulation was quite idealised, since the data were noise-free. We made another simulation, where we added white noise to the TOD. In other aspects the simulation was identical to the first one. No correlated noise was added. We are implicitly assuming that correlated noise components, if present, were removed to sufficient accuracy prior to deconvolution.

As expected, noise made the deconvolution process more difficult, and the code did not converge for the highest values of l_max. We reached full convergence for the combinations l_max = 600, k_max = 6, and l_max = 700, k_max = 4. Case l_max = 700, k_max = 6 did not converge anymore.

Figure 9 shows a zoom into a point source, the same as in Fig. 4, in presence of white noise. The deconvolution parameters were l_max = 600, k_max = 6. The deconvolved map smoothed with the fiducial FWHM = 32′ is contaminated by distorted noise. Also, we see some ringing around the point source due to the moderately low l_max. We therefore smoothed the map further to FWHM = 45′, which was enough to remove the distortion. For comparison we also smoothed the binned map to comparable resolution by applying a smoothing by another FWHM = 32′ to it. This, combined with the original beam, gives a total resolution of roughly FWHM = 45′. These two maps are shown in the bottom row. Even after the additional smoothing, the point source in the binned map is elongated. Deconvolution restores the circular shape of the source.

Figure 10 shows a zoom into the Q polarisation map. We show the same patch of the sky as in Fig. 6. The raw binned map (not shown) is fully noise-dominated. We show maps smoothed to FWHM = 45′. Again, deconvolution removes the galactic leakage that is apparent in the binned map.

Figure 11 shows the spectra constructed from the deconvolved a_Xlm coefficients for the parameter combinations for which the code converged fully. Again we also show the spectrum of the binned map, and the same corrected for a symmetric Gaussian beam with FWHM = 32′. Deconvolution recovers the true TT spectrum up to nearly l = 600, after which noise begins to dominate. The spectrum of the binned map begins to deviate already below l = 400. In the case of the EE spectrum, both deconvolution and binning are unable to recover the true spectra except for the very lowest multipoles. In the case of the TE spectrum, deconvolution performs considerably better than binning, and the result follows the input spectrum up to l = 300, although it is noisy.

Finally, in Fig. 12 we show the error spectra in presence of noise. Deconvolution recovers the a_Tlm coefficients well up to l = 600. In the case of a_Elm coefficients, the error is well above the input spectrum for all but the lowest multipoles.

It is a well-known and often quoted fact that deconvolution tends to amplify noise at high multipoles. We see this clearly in the upper deconvolved map in Fig. 9. The effects can, however, be suppressed by applying a sufficiently strong smoothing to the map, as we see from the lower panels of the same figure.

We can look at the amplification of noise also in harmonic space (Fig. 11). While the TT spectrum of the binned map (green line) continues to fall until it levels out at l = 900 (not seen in the plot), the deconvolution spectrum begins to rise above l = 600. Note, however, that for the whole multipole range where we have a deconvolution result, the extra power in the deconvolved map is still below that in the binned map corrected for a symmetric beam (blue line).

In the EE spectrum, the power of the deconvolved map exceeds that in the binned map corrected for a symmetric beam, around l = 300. This is in the domain where we have full noise domination in any case.

5.4. 70 GHz simulation

To test the performance of the code at higher multipoles, we generated one noise-free simulation mimicking the data voulme of the Planck 70 GHz channel. The used elliptic beams have an average FWHM of 13′; for more parameters see Table 1. We ran the deconvolution code with parameters l_max = 1700, k_max = 6.

Temperature and polarisation maps constructed from the a_Xlm with resolution FWHM = 13′ (T), or FWHM = 18′ (Q) are shown in Figs. 13 and 14 together with the corresponding binned maps. We show the same region of the sky as with 30 GHz simulations. Again, deconvolution has worked well, returning the circular shape of the point source, and removing the galactic leakage in the Q polarisation map.

In Fig. 15 we show the TT and EE error spectra for the 70 GHz simulation. The a_Tlm are recovered well nearly up to l = 1700, and the a_Elm coefficients well above l = 1000. Thus we have shown that, with sufficient computational resources, our deconvolution method can be applied to realistic data sets of volume and resolution of the Planck 70 GHz channel.

5.5. Convergence

The convergence criterion we had set in our code is rather strict. We require that the square of the residual vector has fallen below a fraction of 10^-12 of its initial value.

To see if a less strict criterion would help in reducing the computation time, we studied the convergence of the solution. We dumped the a_Xlm coefficients after every 50 steps, and computed the error spectrum with respect to the input coefficients.

We show the TT and EE spectra after 200 and 500 iteration steps in Fig. 15, along with the final results. We observe that the solution changes only very little after the first few hundred iterations, and 200 steps already give a significant improvement over the harmonic expansion of the binned map. In fact, after the first 100 iteration steps we can see no visible improvement in the reconstructed sky map. We could thus reduce the computation time by a large factor if we used a more relaxed convergence criterion and would not significantly lose accuracy.

5.6. Complicated beam shape

To demonstrate the power of the deconvolution method, we made yet another simulation with an unconventional D-shaped beam (D for “deconvolution”), which is shown in Fig. 16. The linear dimension of the beam was approximately 12° squared. For demonstration purposes we included only point sources and considered temperature only. For the binning process into the 3D map, we chose N_side = 1024 and n_ψ = 256. Deconvolution parameters were l_max = 500, k_max = 30, and only a single full-sky survey was simulated. The comparatively high value of k_max was required to properly take into account the complicated azimuthal structure of the beam.

A map binned directly from the TOD is shown in the left panel of Fig. 17. The beam transforms the image of a point source into a D-shaped pattern, resembling the beam. The deconvolved map is shown on the right; it was smoothed by a symmetric Gaussian beam with FWHM = 1° to reduce ringing. Even with this complicated beam shape we are able to recover the underlying point sources.

It is interesting that in this case it was sufficient to smooth the map with a beam of FWHM = 1° to obtain a very satisfactory map, although the full dimension of the beam is much larger.

6. Conclusions

We have presented an efficient beam deconvolution code, designed for absolute CMB experiments, and tested it on simulated data.

We looked at maps constructed from the recovered a_Tlm, a_Elm, and a_Blm coefficients, and examined the coefficients directly in harmonic space. These reflect different aspects of the solution. We compared the results to a harmonic expansion of a binned map. We have also shown that with sufficient computational resources, we can extend the method to l_max = 1700, which would be sufficient for the Planck 70 GHz channel.

In absence of noise we could recover the a_Tlm coefficients to a high accuracy, and remove the effects of beam asymmetry. When white noise was added, we were able to reach a lower value of l_max, but the advantage of deconvolution over binning was still clear.

In the case of polarisation, deconvolution worked well in absence of noise. We were able to almost completely remove the temperature leakage due to beam mismatch. When noise was added, results were less clear. Deconvolution removed the visible galactic residual that arises from beam mismatch, but in harmonic space deconvolution did not seem to bring a clear benefit over binning. The recovered a_Elm coefficients were dominated by noise at nearly all multipoles, but this was true for both deconvolved and binned maps.

Our simulation used beams that were highly asymmetrical and also had a strong mismatch between them, which leads to a significant polarisation leakage. We did this to show the beam effects more clearly, but at the same time we also made the deconvolution problem very challenging. The accuracy we can expect when the code is applied to a real experiment strongly depends on the beam shapes of the particular experiment.

Some topics for future study can be mentioned. We have set a very strict convergence criterion for the CG iteration. Our convergence study indicates that a much more relaxed criterion could be sufficient for practical purposes. A future improvement would be to define a more suitable convergence criterion. This could reduce the computation time by a significant factor.

We applied the method to a simulation data set that provides full sky coverage. Our code is built on a formalism that makes no assumptions on the sky coverage. In practice, good convergence requires nearly full coverage. This is intuitively evident, since we are solving the harmonic coefficients, which necessarily represent the complete celestial sphere. If the input data only cover a part of the sky, the coefficients cannot be well constrained. A straighforward extension of the method would be to insert a prior that constrains the solution in the region which is not covered by the data, but this is beyond the scope of this paper and is a subject for future study.

Yet another topic for further study is determining the noise bias present in the TT and EE spectra at high multipoles.

Though the simulations used in this work mimic some aspects of the Planck mission, the parameters used do not reflect the properties of the actual instrument. The results presented in this paper are thus not representative of the sensitivity of the Planck experiment, and the authors do not represent the Planck collaboration in this context.

---

Acknowledgments

Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package. We acknowledge the use of the Planck Sky Model, developed by the Component Separation Working Group (WG2) of the Planck Collaboration. M.R. is supported by the German Aeronautics Center and Space Agency (DLR), under program 50-OP-0901, funded by the Federal Ministry of Economics and Technology. E.K. is supported by Academy of Finland grant 253204. This work was granted access to the HPC resources of CSC made available within the Distributed European Computing Initiative by the PRACE-2IP, receiving funding from the European Community’s Seventh Framework Programme (FP7/2007-2013) under grant agreement RI-283493. We thank CSC – the IT Center for Science Ltd. (Finland) – for computational resources. We thank Torsten Enßlin for constructive discussions and feedback, Valtteri Lindholm for help in using the Planck Sky Model, and Ronnie James Dio for the D-beam.

References

Appendix A: Definitions

In this appendix we present definitions and conventions used in the paper. We follow the conventions of Varshalovich et al. (1988).

A.1. Harmonic expansion of the CMB radiation field

The CMB temperature and polarisation fields can be presented by spin-weighted harmonic functions with spin s = 0, ±2. Spin s = 0 represents the temperature field, while components s =  ±2 represent polarisation.

The temperature field on the celestial sphere is expanded as $\mathit{T}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{=}\sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}{\mathit{a}}_{\mathrm{0}\mathit{lm}}{\hspace{0.17em}}_{\mathrm{0}}{\mathit{Y}}_{\mathit{lm}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathit{.}$(A.1)The Q and U Stokes fields, which describe linear polarisation, can be written as $\mathrm{(}\mathit{Q}\mathrm{\pm }\mathrm{i}\mathit{U}\mathrm{)}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{=}\sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}{\mathit{a}}_{\mathrm{\pm }\mathrm{2}\mathit{lm}}{\hspace{0.17em}}_{\mathrm{\pm }\mathrm{2}}{\mathit{Y}}_{\mathit{lm}}\mathrm{.}$(A.2)Here ${\mathit{s}}_{}{\mathit{Y}}_{\mathit{lm}}$ are spin-weighted harmonic functions, and a_slm are the corresponding (complex) coefficients. The harmonic functions and coefficients obey the symmetry relation $\begin{array}{ccc}{}_{\mathit{s}}{\mathit{Y}}_{\mathit{lm}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}& & \mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{m}\mathrm{-}\mathit{s}}{}_{\mathrm{-}\mathit{s}}{\mathit{Y}}_{\mathit{l,}\mathrm{-}\mathit{m}}^{\mathrm{\ast }}\mathrm{(}\mathit{\theta ,\phi }\mathrm{)}\\ {\mathit{a}}_{\mathit{slm}}& & \mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{m}\mathrm{-}\mathit{s}}{\mathit{a}}_{\mathrm{-}\mathit{s,l,}\mathrm{-}\mathit{m}}^{\mathrm{\ast }}\mathrm{.}\end{array}$(A.3)The harmonic functions are related to the Wigner functions ${\mathit{D}}_{\mathit{sm}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta ,\phi ,\psi }\mathrm{\right)}$ through ${\mathit{s}}_{}{\mathit{Y}}_{\mathit{lm}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathit{l}}{\mathit{D}}_{\mathit{m,}\mathrm{-}\mathit{s}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}\mathit{\phi ,\theta ,}\mathrm{0}\mathrm{\right)}\mathit{,}$(A.4)where ${\mathit{F}}_{\mathit{l}}\mathrm{=}{\left(\frac{\mathrm{2}\mathit{l}\mathrm{+}\mathrm{1}}{\mathrm{4}\mathit{\pi }}\right)}^{\mathrm{1}\mathit{/}\mathrm{2}}\mathrm{·}$(A.5)It is convenient to define further $\begin{array}{ccc}{\mathit{a}}_{\mathit{Tlm}}\mathrm{=}& & \hspace{0.17em}{\mathit{a}}_{\mathrm{0}\mathit{lm}}\\ {\mathit{a}}_{\mathit{Elm}}\mathrm{=}& & \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}}\mathrm{(}{\mathit{a}}_{\mathrm{2}\mathit{lm}}\mathrm{+}{\mathit{a}}_{\mathrm{-}\mathrm{2}\mathit{lm}}\mathrm{)}\\ {\mathit{a}}_{\mathit{Blm}}\mathrm{=}& & \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{i}}\mathrm{(}{\mathit{a}}_{\mathrm{2}\mathit{lm}}\mathrm{-}{\mathit{a}}_{\mathrm{-}\mathrm{2}\mathit{lm}}\mathrm{)}\mathrm{.}\end{array}$(A.6)The CMB radiation field is expected to be a statistically isotropic Gaussian spin 0, 2 field. The statistical properties of the harmonic expansion coefficients are then fully characterised by four spectra ${\mathit{C}}_{\mathit{l}}^{\mathit{T}}\mathit{,}{\mathit{C}}_{\mathit{l}}^{\mathit{E}}\mathit{,}{\mathit{C}}_{\mathit{l}}^{\mathit{B}}\mathit{,}{\mathit{C}}_{\mathit{l}}^{\mathit{C}}$, such that $\begin{array}{ccc}\mathrm{⟨}{\mathit{a}}_{\mathit{Tlm}}^{\mathrm{\ast }}{\mathit{a}}_{\mathit{T}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{⟩}& & \mathrm{=}{\mathit{C}}_{\mathit{l}}^{\mathit{T}}{\mathit{\delta }}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\\ \mathrm{⟨}{\mathit{a}}_{\mathit{Elm}}^{\mathrm{\ast }}{\mathit{a}}_{\mathit{E}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{⟩}& & \mathrm{=}{\mathit{C}}_{\mathit{l}}^{\mathit{E}}{\mathit{\delta }}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\\ \mathrm{⟨}{\mathit{a}}_{\mathit{Blm}}^{\mathrm{\ast }}{\mathit{a}}_{\mathit{B}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{⟩}& & \mathrm{=}{\mathit{C}}_{\mathit{l}}^{\mathit{B}}{\mathit{\delta }}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\\ \mathrm{⟨}{\mathit{a}}_{\mathit{Tlm}}^{\mathrm{\ast }}{\mathit{a}}_{\mathit{E}{\mathit{l}}^{\mathrm{\prime }}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{⟩}& & \mathrm{=}{\mathit{C}}_{\mathit{l}}^{\mathit{C}}{\mathit{\delta }}_{\mathit{l}{\mathit{l}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\mathrm{.}\end{array}$(A.7)

A.2. Wigner functions and rotations on a sphere

The harmonic coefficients transform under rotations as follows.

Assume we are given the coefficients a_slm in one coordinate system. We then rotate the coordinate system, as described by Euler angles φ,θ,ψ:

• 1.

  Rotate by angle ψ around the z-axis (in positive direction).

• 2.

  Rotate by angle θ around the initial y-axis.

• 3.

  Rotate by angle φ around the initial z-axis.

The rotation is equivalent to rotating first by angle φ around the z-axis, then by angle θ around the new y-axis, and finally by angle ψ around the new z-axis. The harmonic coefficients in the rotated coordinate system are given by ${\mathit{a}}_{\mathit{slm}}^{\mathrm{\prime }}\mathrm{=}\sum _{{\mathit{m}}^{\mathrm{\prime }}}{\mathit{a}}_{\mathit{sl}{\mathit{m}}^{\mathrm{\prime }}}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}\mathit{m}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right)}\mathrm{,}$(A.8)where ${\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}\mathit{m}}^{\mathit{l}}$ are Wigner functions.

A Wigner function may be split into a product of three terms, each of which only depends on one of the angles: ${\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta ,\phi ,\psi }\mathrm{\right)}\mathrm{=}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{m\phi }}{\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k\psi }}\mathit{.}$(A.9)Functions ${\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$ are the reduced Wigner functions. They have the property ${\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}\mathrm{=}{\mathit{\delta }}_{\mathit{mk}}\mathit{.}$(A.10)They are real-valued and obey the symmetry relation ${\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{m}\mathrm{+}\mathit{k}}{\mathit{d}}_{\mathrm{-}\mathit{m,}\mathrm{-}\mathit{k}}^{\mathit{l}}\mathrm{(}\mathit{\theta }\mathrm{)}\mathit{.}$(A.11)The symmetry relation for the full Wigner functions is ${\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right)}\mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{\mathit{m}\mathrm{+}\mathit{k}}{\mathit{D}}_{\mathrm{-}\mathit{m,}\mathrm{-}\mathit{k}}^{\mathit{l}\mathrm{\ast }}\mathrm{(}\mathit{\phi ,\theta ,\psi }\mathrm{)}\mathit{.}$(A.12)

A.3. Wigner expansion

Wigner functions ${\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right)}$ for integer values of l,m,k are defined in domain V: $\mathit{V}\mathrm{:} \mathrm{0}\mathrm{\le }\mathit{\phi }\mathrm{\le }\mathrm{2}\mathit{\pi ,} \mathrm{0}\mathrm{\le }\mathit{\theta }\mathrm{\le }\mathit{\pi ,} \mathrm{0}\mathrm{\le }\mathit{\psi }\mathrm{\le }\mathrm{2}\mathit{\pi }\mathrm{.}$(A.13)Domain V can be interpreted as the space of Euler angles. The volume element in domain V is dV = sinθdθdφdψ, and V’s total volume is 8π².

Any finite function f(φ,θ,ψ) defined in domain V may be expanded in a series of Wigner functions as $\mathit{f}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right)}\mathrm{=}\sum _{\mathit{l}\mathrm{=}\mathrm{0}}^{\mathrm{\infty }}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}\sum _{\mathit{k}\mathrm{=}\mathrm{-}\mathit{l}}^{\mathit{l}}{\mathit{c}}_{\mathit{mk}}^{\mathit{l}}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right)}\mathrm{,}$(A.14)where the coefficients ${\mathit{c}}_{\mathit{mk}}^{\mathit{l}}$ are given by ${\mathit{c}}_{\mathit{mk}}^{\mathit{l}}\mathrm{=}\frac{\mathrm{2}\mathit{l}\mathrm{+}\mathrm{1}}{\mathrm{8}{\mathit{\pi }}^{\mathrm{2}}}{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}\mathit{\phi }{\mathrm{\int }}_{\mathrm{0}}^{\mathit{\pi }}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }{\mathrm{\int }}_{\mathrm{0}}^{\mathrm{2}\mathit{\pi }}\mathrm{d}\mathit{\psi f}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right)}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right)}\mathrm{.}$(A.15)Expansion (A.14)can be understood as the 3D equivalent of the spherical harmonics expansion.

In particular, a product of two Wigner functions ${\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}$ can be expanded as a series of Wigner functions as $\begin{array}{ccc}& & \\ & & \sum _{{\mathit{l}}_{\mathrm{2}}{\mathit{m}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{2}}}{\mathit{D}}_{{\mathit{m}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{2}}}^{{\mathit{l}}_{\mathrm{2}}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}\frac{\mathrm{2}{\mathit{l}}_{\mathrm{2}}\mathrm{+}\mathrm{1}}{\mathrm{8}{\mathit{\pi }}^{\mathrm{2}}}\mathrm{\int }{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{\omega }}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{2}}}^{{\mathit{l}}_{\mathrm{2}}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{\omega }}^{\mathrm{\prime }}\mathit{.}\end{array}$(A.16)An integral over a product of three Wigner functions is given by Wigner 3j symbols as $\begin{array}{ccc}\int {\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}^{\mathrm{\prime }}\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}{\mathit{\omega }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{=}{\mathit{D}}_{{\mathit{m}}_{\mathrm{2}}{\mathit{k}}_{\mathrm{2}}}^{{\mathit{l}}_{\mathrm{2}}\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}^{\mathrm{\prime }}\mathrm{\right)}\mathrm{d}{\mathit{\omega }}^{\mathrm{\prime }}& & \\ \mathrm{\times }\mathrm{8}{\mathit{\pi }}^{\mathrm{2}}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{m}}^{\mathrm{\prime }}\mathrm{+}{\mathit{k}}^{\mathrm{\prime }}}\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)& & \end{array}$(A.17)

(Varshalovich et al. 1988).

The following properties of the 3j symbols help in restricting the range of summation. A 3j symbol is non-vanishing $\left(\begin{array}{c}\\ {\mathit{l}}_{\mathrm{1}}& {\mathit{l}}_{\mathrm{2}}& {\mathit{l}}_{\mathrm{3}}\\ {\mathit{k}}_{\mathrm{1}}& {\mathit{k}}_{\mathrm{2}}& {\mathit{k}}_{\mathrm{3}}\end{array}\right)\ne \mathrm{0}$(A.18)only if the two conditions ${\mathit{k}}_{\mathrm{1}}\mathrm{+}{\mathit{k}}_{\mathrm{2}}\mathrm{+}{\mathit{k}}_{\mathrm{3}}\mathrm{=}\mathrm{0}$(A.19)and $\mathrm{|}{\mathit{l}}_{\mathrm{1}}\mathrm{-}{\mathit{l}}_{\mathrm{2}}\mathrm{|}\mathrm{\le }{\mathit{l}}_{\mathrm{3}}\mathrm{\le }{\mathit{l}}_{\mathrm{1}}\mathrm{+}{\mathit{l}}_{\mathrm{2}}$(A.20)are fulfilled. We can thus eliminate the summation over k₂ and m₂. Product (A.16)becomes $\begin{array}{ccc}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}{\mathit{k}}^{\mathrm{\prime }}}^{{\mathit{l}}^{\mathrm{\prime }}}\mathrm{\left(}\mathit{\omega }\mathrm{\right)}& & \mathrm{=}\mathrm{(}\mathrm{-}\mathrm{1}{\mathrm{)}}^{{\mathit{m}}^{\mathrm{\prime }}\mathrm{+}{\mathit{k}}^{\mathrm{\prime }}}\sum _{{\mathit{l}}_{\mathrm{2}}}\mathrm{(}\mathrm{2}{\mathit{l}}_{\mathrm{2}}\mathrm{+}\mathrm{1}\mathrm{\left)}{\mathit{D}}_{{\mathit{m}}^{\mathrm{\prime }}\mathrm{-}\mathit{m,}{\mathit{k}}^{\mathrm{\prime }}\mathrm{-}\mathit{k}}^{{\mathit{l}}_{\mathrm{2}}}\mathrm{\right(}\mathit{\omega }\mathrm{)}\\ & & \mathrm{\times }\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\left(\begin{array}{c}\\ & & \\ & & \end{array}\right)\mathit{.}\end{array}$(A.21)

A.4. Beam

A.4.1. Detector signal and definition of beam coefficients

Consider a detector measuring CMB temperature and linear polarisation.

We assume now, very generally, that the signal recorded by a detector, apart from noise, is proportional to the a_slm coefficients. We define the beam coefficients b_slk as follows. When the detector is at some (predefined) fiducial position and orientation (for instance at the north pole), the signal detected in the absence of noise is ${\mathit{t}}_{\mathrm{fid}}\mathrm{=}\begin{array}{}\sum \end{array}{\mathit{a}}_{\mathit{slm}}{\mathit{b}}_{\mathit{slm}}^{\mathrm{\ast }}\mathit{.}$(A.22)It is clear from the definition that the beam coefficients must contain all information not only of the beam shape, but also of the possible polarisation sensitivity and orientation.

We may then rotate the beam to another position and orientation, as described in Sect. A.2. The rotation brings the beam onto location θ,φ in spherical coordinates, in an orientation defined by angle ψ. The harmonic coefficients transform under rotations according to (A.8). The signal seen by the rotated detector is then $\mathit{t}\mathrm{=}\sum _{\mathit{slmk}}{\mathit{a}}_{\mathit{slm}}{\mathit{b}}_{\mathit{slk}}^{\mathrm{\ast }}{\mathit{D}}_{\mathit{mk}}^{\mathit{l}\mathrm{\ast }}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right)}\mathit{.}$(A.23)In the following we derive expressions for the beam coefficients in some special cases.

A.4.2. Non-polarised beam

Consider first a detector with an ideal beam, beam width equal to zero, and no sensitivity to polarisation. The signal detected by the beam is directly given by (A.1). We see readily, by comparing (A.1)and (A.22), that if the chosen fiducial beam position is at the north pole (θ = 0,φ = 0), the beam coefficients must be ${\mathit{b}}_{\mathit{slm}}\mathrm{=}{\mathit{s}}_{}{\mathit{Y}}_{\mathit{lm}}^{\mathrm{\ast }}\mathrm{\left(}\mathrm{0}\mathit{,}\mathrm{0}\mathrm{\right)}{\mathit{\delta }}_{\mathit{s}\mathrm{0}}\mathrm{=}{\mathit{F}}_{\mathit{l}}{\mathit{\delta }}_{\mathit{s}\mathrm{0}}{\mathit{\delta }}_{\mathit{m}\mathrm{0}}\mathit{.}$(A.24)A general non-polarised beam with sensitivity g(θ,φ) gives for the s = 0 components $\begin{array}{ccc}{\mathit{b}}_{\mathrm{0}\mathit{lm}}& & \mathrm{=}”\mathit{g}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{}_{\mathrm{0}}{\mathit{Y}}_{\mathit{lm}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }\\ & & \mathrm{=}{\mathit{F}}_{\mathit{l}}”\mathit{g}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathit{D}}_{\mathit{m}\mathrm{0}}^{\mathit{l}}\mathrm{\left(}\mathit{\phi ,\theta ,}\mathrm{0}\mathrm{\right)}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }\end{array}$(A.25)and b_slm = 0 for s =  ± 2.

A.4.3. Polarisation measurements

A polarisation-sensitive detector requires careful treatment, because the Stokes parameters Q and U are not defined at the pole. The beam coefficients can still be defined in a meaningful way, as we show in the following.

Consider again a beam at a fiducial position and orientation at the north pole. A general beam may have different polarisation sensitivity at different locations. We therefore define the beam shape g(θ,φ) separately for temperature and polarisation.

The signal detected may be written in a general form as $\begin{array}{ccc}\mathit{t}& & \mathrm{=}”\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }\{{\mathit{g}}_{\mathit{T}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathit{T}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\\ & & \mathrm{+}{\mathit{g}}_{\mathit{P}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\left[\mathit{Q}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\cos \mathrm{\left(}\mathrm{2}\mathit{\psi }\mathrm{\right(}\mathit{\theta ,\phi }\mathrm{\left)}\mathrm{\right)}\mathrm{+}\mathit{U}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\sin \mathrm{\left(}\mathrm{2}\mathit{\psi }\mathrm{\right(}\mathit{\theta ,\phi }\mathrm{\left)}\mathrm{\right)}\right]\}\mathrm{.}\end{array}$(A.26)Here ψ(θ,φ) is the local direction of polarisation sensitivity, measured as an angle from the local meridian. Functions g_T(θ,φ), g_P(θ,φ), and ψ(θ,φ) together define a general beam.

We proceed to calculate the beam coefficients. We write the trigonometric functions in exponential form, to obtain $\begin{array}{ccc}\mathit{t}& & \\ & & \mathrm{+}\frac{{\mathit{g}}_{\mathit{P}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}}{\mathrm{2}}\mathrm{[}\mathrm{(}\mathit{Q}\mathrm{-}\mathrm{i}\mathit{U}\mathrm{)}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathrm{e}}^{\mathit{i}\mathrm{2}\mathit{\psi }\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}}\mathrm{+}\mathrm{(}\mathit{Q}\mathrm{+}\mathrm{i}\mathit{U}\mathrm{)}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{2}\mathit{\psi }\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}}\mathrm{]}\}\mathrm{.}\end{array}$(A.27)We can now insert the harmonic expansions (A.1)and (A.2)into (A.27)and read the beam coefficients $\begin{array}{ccc}{\mathit{b}}_{\mathrm{0}\mathit{lk}}\mathrm{=}& & ”{\mathit{g}}_{\mathit{T}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\hspace{0.17em}}_{\mathrm{0}}{\mathit{Y}}_{\mathit{lk}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }\\ {\mathit{b}}_{\mathrm{\pm }\mathrm{2}\mathit{lk}}\mathrm{=}& & \frac{\mathrm{1}}{\mathrm{2}}”{\mathit{g}}_{\mathit{P}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\hspace{0.17em}}_{\mathrm{\pm }\mathrm{2}}{\mathit{Y}}_{\mathit{lk}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathrm{e}}^{\mathrm{\pm }\mathrm{i}\mathrm{2}\mathit{\psi }\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }\mathrm{.}\end{array}$(A.28)We then use definition (A.4)and write the harmonic functions in terms of Wigner functions. In terms of the reduced Wigner functions we have ${\mathit{s}}_{}{\mathit{Y}}_{\mathit{lk}}^{\mathrm{\ast }}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{=}{\mathit{F}}_{\mathit{l}}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k\phi }}{\mathit{d}}_{\mathit{k,}\mathrm{-}\mathit{s}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{.}$(A.29)The ψ-term can be transferred inside the Wigner function. We may thus write for a general polarised beam $\begin{array}{ccc}{\mathit{b}}_{\mathrm{0}\mathit{lk}}& & \mathrm{=}{\mathit{F}}_{\mathit{l}}”{\mathit{g}}_{\mathit{T}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathit{D}}_{\mathit{k,}\mathrm{0}}^{\mathit{l}}\mathrm{\left(}\mathit{\phi ,\theta ,}\mathrm{0}\mathrm{\right)}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }\\ {\mathit{b}}_{\mathrm{\pm }\mathrm{2}\mathit{lk}}& & \mathrm{=}\frac{{\mathit{F}}_{\mathit{l}}}{\mathrm{2}}”{\mathit{g}}_{\mathit{P}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathit{D}}_{\mathit{k,}\mathrm{\mp }\mathrm{2}}^{\mathit{l}}\mathrm{\left(}\mathit{\phi ,\theta ,\psi }\mathrm{\right(}\mathit{\theta ,\phi }\mathrm{\left)}\mathrm{\right)}\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\mathrm{d}\mathit{\phi }\mathrm{.}\end{array}$(A.30)The above formulas for the beam coefficients are completely general as long as the beam response is linear.

We now make the simplifying assumption that the polarisation sensitivity is “in the same direction” everywhere on the beam. By “same direction” we mean that the direction of polarisation sensitivity at an arbitrary location is obtained by rotating the polarisation direction vector at the north pole to the desired location along a meridian. The local polarisation orientation angle then becomes $\mathit{\psi }\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}\mathrm{=}{\mathit{\psi }}_{\mathrm{pol}}\mathrm{-}\mathit{\phi }\mathrm{,}$(A.31)where ψ_pol is the deviation of polarisation direction from φ = 0 at the north pole.

At this point it is convenient to write the Wigner functions in terms of the reduced Wigner functions (A.9). The coefficients for a beam with unique polarisation direction become $\begin{array}{ccc}{\mathit{b}}_{\mathrm{0}\mathit{lk}}& & \mathrm{=}{\mathit{F}}_{\mathit{l}}\mathrm{\int }\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\hspace{0.17em}{\mathit{d}}_{\mathit{k,}\mathrm{0}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{\int }{\mathit{g}}_{\mathit{T}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathit{k\phi }}\mathrm{d}\mathit{\phi }\\ {\mathit{b}}_{\mathrm{\pm }\mathrm{2}\mathit{lk}}& & \mathrm{=}\frac{{\mathit{F}}_{\mathit{l}}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{\pm }\mathrm{i}\mathrm{2}{\mathit{\psi }}_{\mathrm{pol}}}\mathrm{\int }\sin \mathit{\theta }\mathrm{d}\mathit{\theta }\hspace{0.17em}{\mathit{d}}_{\mathit{k,}\mathrm{\mp }\mathrm{2}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}\mathrm{\int }{\mathit{g}}_{\mathit{P}}\mathrm{\left(}\mathit{\theta ,\phi }\mathrm{\right)}{\mathrm{e}}^{\mathrm{-}\mathrm{i}\mathrm{(}\mathit{k}\mathrm{\pm }\mathrm{2}\mathrm{)}\mathit{\phi }}\mathrm{d}\mathit{\phi }\mathrm{.}\end{array}$(A.32)We readily see that for a symmetric beam the only non-vanishing coefficients are those with k + s = 0. Especially, for a perfect beam with zero width and perfect polarisation sensitivity (g_P = g_T) we obtain $\begin{array}{ccc}{\mathit{b}}_{\mathrm{0}\mathit{lk}}& & \mathrm{=}{\mathit{F}}_{\mathit{l}}{\mathit{d}}_{\mathit{k,}\mathrm{0}}^{\mathit{l}}\mathrm{\left(}\mathrm{0}\mathrm{\right)} \mathrm{=} {\mathit{F}}_{\mathit{l}}{\mathit{\delta }}_{\mathit{k,}\mathrm{0}}\\ {\mathit{b}}_{\mathrm{\pm }\mathrm{2}\mathit{lk}}& & \mathrm{=}\frac{{\mathit{F}}_{\mathit{l}}}{\mathrm{2}}{\mathit{d}}_{\mathit{k,}\mathrm{\mp }\mathrm{2}}^{\mathit{l}}\mathrm{\left(}\mathrm{0}\mathrm{\right)}{\mathrm{e}}^{\mathrm{\pm }\mathrm{i}\mathrm{2}{\mathit{\psi }}_{\mathrm{pol}}} \mathrm{=} \frac{{\mathit{F}}_{\mathit{l}}}{\mathrm{2}}{\mathrm{e}}^{\mathrm{\pm }\mathrm{i}\mathrm{2}{\mathit{\psi }}_{\mathrm{pol}}}{\mathit{\delta }}_{\mathit{k,}\mathrm{\mp }\mathrm{2}}\mathrm{.}\end{array}$(A.33)We have assumed normalisation ${}^{\mathrm{\int }}\mathit{g}_{\mathit{T}}\mathrm{\left(}\mathrm{\Omega }\mathrm{\right)}\mathrm{d\Omega }\mathrm{=}\mathrm{1}$.

In this work we have constructed the beam coefficients according to (A.32), but the assumption of unique polarization direction is not crucial for the deconvolution method. A more general beam presentation is given by (A.30).

Appendix B: Methods used to accelerate computation

B.1. Load balancing

In an MPI-parallel code, good load balancing among the individual tasks is crucial. As was mentioned in Sect. 3.2, some quantities are distributed across tasks depending on their colatitude θ or the index m. Since the computational cost of the algorithm is not independent of θ and m, but rather an unknown, smoothly varying function of these two numbers, it is advisable not to assign sequential ranges of those indices to individual MPI tasks, because choosing index ranges that distribute the load evenly is very hard in that scenario.

To avoid this complication, we adopted a “round robin” strategy for data distribution. Assuming n MPI tasks numbered 0 to n − 1, we assign to a particular task i the θ indices θ_i, θ_i + n, θ_i + 2n, ..., and deal with the index m in analogous fashion. This results in near-perfect load balancing.

B.2. Convolution optimisation

A significant part of total CPU time is spent in the FFTs necessary for the convolution with the quantity N (see Sect. 3.2). Fortunately, several measures can be taken to decrease the resource requirements. First of all, we are making use of the highly optimised FFTW⁴ library (Frigo & Johnson 1998), which provides close to optimal performance for FFTs of specific size.

Furthermore, it is important to notice that the execution time for an FFT of length n depends sensitively on the prime factorisation of n; in general, an n composed of only small prime factors is much preferable.

Since we are performing a convolution operation, we have some freedom in the choice of array dimensions; they can be increased and zero-padded if desired. In our case, the minimal array dimensions are 4l_max + 1 and 4k_max + 1, respectively; if necessary, we increase both of them to the next larger number, which is a product of the prime factors 2, 3, and 5 only. Depending on the prime factorisation of the original array dimensions, this can decrease the run-time for this part of the code by integral factors.

Finally, the quantities that are to be convolved exhibit the symmetry

${\mathit{V}}_{\mathit{i,k}}\mathrm{=}{\mathit{V}}_{{\mathit{n}}_{\mathit{i}}\mathrm{-}\mathit{i,}{\mathit{n}}_{\mathit{k}}\mathrm{-}\mathit{k}}^{\mathrm{\ast }} \mathrm{for}\mathit{i,k}\mathit{>}\mathrm{0}\mathrm{,}$which allows us to employ the specialised “real-valued” FFT, resulting in another substantial speedup and a size reduction of the precomputed Fourier transform of N (see Sect. B.5).

B.3. Efficient computation of the reduced Wigner functions

During code execution, the reduced Wigner functions ${\mathit{d}}_{\mathit{mk}}^{\mathit{l}}\mathrm{\left(}\mathit{\theta }\mathrm{\right)}$ are required twice in every conjugate gradient iteration step. Storing them in memory is prohibitively expensive, so it is important to regenerate them efficiently whenever needed. For this purpose we used code originally presented by Prézeau & Reinecke (2010), which was extended to make use of the SSE2 instruction set present in all modern Intel-compatible CPUs. We also used the symmetry properties of the d matrix to avoid redundant computations.

B.4. Shortcuts for symmetrical beams

If a beam used for deconvolution is symmetric with respect to a 180° rotation around its axis (which is true for elliptical beams, for example), all its b_slm coefficients for odd m vanish by definition. This fact is used to skip unnecessary computations as well as save space by not storing any array elements that are known to be zero. It also reduces the minimum second dimension of the convolution arrays from 4k_max + 1 to 2k_max + 1, saving even more time on the FFT operations and reducing memory consumption.

B.5. Optional precomputation

The convolution step mentioned in Sect. 3.2 requires the Fourier transform of array N, which is constant throughout a run. The code permits one to precompute and store this quantity, which reduces total CPU time (because only two FFT operations are then needed for each convolution step, instead of three). However, especially for runs with high l_max and k_max this increases the total memory consumption considerably. So for situations where the available main memory is insufficient to store this precomputed array, we provide the option of recomputing the Fourier-transformed N whenever necessary, approximately halving the memory consumption of the code, but also slowing it down by 5–15 percent. For the 70 GHz runs described in this paper, we had to make use of this space-saving feature.

B.6. Supplying an initial guess

While not strictly a performance optimisation, the code allows the user to supply an initial guess from which to start the CG iteration, which is then used instead of a vector of zeroes. This is particularly helpful if one tries to compute series of deconvolutions on the same input data set, but with different parameters. For example, if a solution for a deconvolution problem with a given l_max and k_max has already been obtained, it can be used as a good starting point for other calculations with higher l_max and/or k_max. In any case, the convergence criterion will still be evaluated with respect to an all-zero a_slm, not the supplied guess.

Appendix C: Resource usage

Most computations discussed in this paper were performed on a computer with 64 GB of RAM and 24 Intel Xeon E7450 cores. Since this computer is used interactively, we only ran our calculations on 16 of these cores and ensured that interactive usage did not exceed the capability of the remaining 8 cores. Nevertheless, the measured timings contain small uncertainties and should be interpreted as upper limits.

Table C.1 gives an overview of the parameters of the various runs and their impact on resource usage. Runs 30a–g are the 30 GHz calculations discussed in Sect. 5.2; runs 70a, 70x, and 70y are related to the 70 GHz high-resolution study (Sect. 5.4), and run “D” is the deconvolution with the D-shaped beam presented in Sect. 5.6. For all runs except the D-beam, the beam was assumed to be symmetrical with respect to a 180° rotation around its axis. Also, for all runs except the 70 GHz ones the FFT of array N was cached (see Sect. B.5).

For the 70 GHz run on the test machine described above (run 70a), each CG iteration step needed significantly more time than for the calculations using the 30 GHz data set; in combination with the high number of iterations, this leads to an inconveniently long computation time. As a consequence, deconvolution runs of this magnitude are better suited for execution on massively parallel computers. We therefore repeated the computation on the Louhi supercomputer of CSC (Cray XT4/XT5) on 128 and 512 processor cores, respectively. The corresponding performance figures are listed as runs 70x and 70y in Table C.1; memory consumption was not measured for these runs.

The computational power of a single CPU core is roughly the same on both computers used for our tests. Looking at the total wall clock time for the 70 GHz runs, it becomes obvious that the scaling from 16 to 128 cores is close to ideal, but becomes much worse when going from 128 to 512 cores. In this range, MPI communication and redundant computations begin to dominate. For deconvolution problems with larger l_max, we expect that this transition will occur at a higher number of MPI tasks.

The numbers of iterations required for the 70 GHz runs are not identical, as one might intuitively expect. This is because in an MPI-parallel code some computation steps – most notably summations and other reductions of quantities residing on all tasks – are not carried out in a strictly defined order, but depend on the relative timing of the tasks during such an operation, which is nondeterministic. Since floating-point algebra with finite precision is non-associative, the results of such operations can differ slightly from run to run, causing the conjugate gradient solver to take different paths towards the solution, but of course ultimately solving the same numerical problem.

Compared to earlier publications on deconvolution map-making, it is evident that both memory and CPU requirements of the presented algorithm are quite modest and allow deconvolution up to values of l_max that were prohibitively expensive up to now. Making use of massively parallel computing clusters, we expect that the algorithm can be applied to data sets with even higher resolution than those presented here.

All Tables

All Figures