FAST PIXEL SPACE CONVOLUTION FOR COSMIC MICROWAVE BACKGROUND SURVEYS WITH ASYMMETRIC BEAMS AND COMPLEX SCAN STRATEGIES: FEBeCoP

Cosmic microwave background (CMB) anisotropy measurements have played a leading role in precision cosmology. Detection of the CMB monopole by Penzias & Wilson (1965) firmly established that the big bang model is consistent with observation. The Cosmic Background Explorer (Smoot et al. 1992) detected the tiny ∼10⁻⁵ fluctuations, which justified how structures could form in a statistically isotropic universe starting from seed density perturbations. The Wilkinson Microwave Anisotropy Probe (Larson et al. 2011) measured cosmological parameters—e.g., the flatness of the universe—with unprecedented accuracy. There have been many suborbital experiments (Torbet et al. 1999; Miller et al. 1999; Padin et al. 2001; Pearson et al. 2003; Kuo et al. 2004; de Bernardis et al. 2000; Hanany et al. 2000) that have probed specific aspects of CMB anisotropy. Yet only a small fraction of the information embedded in the CMB temperature anisotropy has so far been measured, and measurements of the CMB polarization anisotropy have just begun. Precise measurement of the CMB temperature and polarization anisotropy is the most important next step in building a rigorous model of the content and evolution of the universe. The Planck mission is a substantial step forward in CMB anisotropy measurement—it aims to measure almost all the information available in temperature anisotropy and will make a full-sky, high-resolution, low-noise measurement of CMB polarization anisotropy. It will also image thousands of point sources of astrophysical and cosmological interest.

Planck (The Planck Collaboration 2006; Tauber et al. 2010a) is a full-sky scanning experiment launched in 2009 May, taking data since 2009 August. It will produce high-resolution maps with tens of millions of pixels. The low-frequency instrument (LFI) covers 30, 44, and 70 GHz; the high-frequency instrument (HFI) covers 100, 143, 216, 353, 545, and 857 GHz. From the second Lagrangian point of the Earth–Sun system (L₂), Planck scans nearly great circles on the sky, covering the full sky twice in one year (Dupac & Tauber 2005). The satellite spins at 1 rpm around an axis that is re-pointed roughly 30 times per day along a cycloidal path, with the spin axis moving in a 75 circle around the anti-Sun direction with a period of 6 months. This ensures that all feeds completely cover the sky, including the ecliptic polar regions.

High-resolution, low-noise measurement of the CMB anisotropy is challenging. Cosmic-variance-limited measurement at high resolution demands not only precision, but very high accuracy—all the systematic effects must be accounted for; otherwise systematic errors can become significant or dominant compared to statistical errors. The effect of asymmetric beams is one of the most important sources of systematic error in current and future CMB missions, which invest a significant effort in measuring the polarization. Problems with asymmetric beam systematics, which could transfer power from large scales to small, could be acute for measurements of the rapidly falling power spectra at high resolution. For polarization (Ashdown et al. 2009), where the signal is an order of magnitude below the temperature anisotropy and the inherent directionality of the polarization field gets coupled to the asymmetry of the beam, the effect could lead to large deviations in the polarization anisotropy. In general, precision measurement of the polarization power spectrum at high resolution requires a correction for asymmetric beams.

Planck beams range from ∼30' at 30 GHz to ∼5' at 217 GHz and higher frequencies. The lowest frequency channels are farthest from the center of the focal plane and exhibit the most asymmetric beams. For example, the 30 GHz channel has beams with ellipticities of the order 1.37. In Figure 1, we show realistic Planck beams, for one 30 GHz and one 143 GHz detector, simulated by a full diffraction analysis of the telescope using GRASP9 (Sandri et al. 2002, 2010; Maffei et al. 2010; Yurchenko et al. 2004), while in Figure 2 we show the focal plane layout (Dupac & Tauber 2005). The related problems of measurement of beam shapes (Huffenberger et al. 2010) and the incorporation of their uncertainties in the analysis (Rocha et al. 2010) are also active areas of research. Our focus in this paper is to study the effect of asymmetries in the beam shapes on the sky maps.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Broadly speaking, the problem of asymmetric beams at high resolution has been addressed in the literature in two ways: estimating the bias of the "pseudo-C_l" estimator semi-analytically in the spherical harmonic basis (Mitra et al. 2009; Hinshaw et al. 2007; Mitra et al. 2004; Fosalba et al. 2002; Souradeep & Ratra 2001); or deconvolution of observed maps (Armitage-Caplan & Wandelt 2009; Burigana & Saez 2003). However, to handle a complicated instrument like Planck, with non-trivial scan strategy and noise characteristics, we aim to use a simple forward approach, i.e., to estimate the systematic effect by studying Monte Carlo simulated maps convolved with detector beams using the real scan path. This requires a fast convolution algorithm for scanning CMB experiments like Planck. Existing simulation methods (Reinecke et al. 2006), involving spherical harmonic based convolution algorithms, are relatively slow, and even computationally prohibitive for a large number of simulations (typically a few thousand in Monte Carlo studies). This is because they generate time-ordered data (TOD) prior to mapmaking. Improving simulation capabilities was the main motivation for developing a fast, easy-to-use map convolution method which would include the effect of asymmetric beams and scan.

Convolution in pixel space is fast and flexible because in scanning CMB experiments the main beam is highly localized in pixel space. The standard choice of pixel number N_pix for a given beam size ensures that only N_beam ∼ 200 pixels are required to accurately describe the beam. This is much smaller than the number of spherical harmonic moments needed to describe the beam with adequate accuracy. The total computation cost for pixel space convolution (after effective beam computation) scales as N_pix N_beam. This is very small compared to any other existing method. For example, using convolution in the harmonic basis for an idealized scan (where each direction has a fixed beam orientation angle), performed by the fast multipole based method "total convolution" (Wandelt & Górski 2001; Challinor et al. 2000), scales as ∼ℓ³_max m_max, where ℓ_max is the highest multipole moment that contains non-negligible information about the anisotropies in the map and m_max is the highest azimuthal multipole needed to adequately describe beam asymmetry. This is more demanding than the pixel space convolution. There has been recent progress in making this step faster using recursive algebra (Prezeau & Reinecke 2010). In practice, though, convolution is much more complicated than the ideal case, as the scan can take the beam to any arbitrary orientation and one has to go through the whole scan path and interpolate between orientation angles obtained from the fixed Total Convolution grid, which is a tedious process. Considering all these issues, pixel space is often a more convenient choice.

We have developed a HEALPix⁹ (Górski et al. 2005) based pixel space convolution method by pre-computing the "effective beams" for each pixel for both temperature and polarization maps, incorporating asymmetric beams of the instruments and the real scan strategy. The effective beam for a given pixel is essentially the average of (discretized) beam functions for every pointing sample that fall in that pixel (a formal mathematical definition is provided later in the paper). The method has been successfully implemented for the Planck mission as a computer code FEBeCoP, acronym of "Fast Effective Beam Convolution in Pixel space." FEBeCoP efficiently utilizes large distributed computing facilities, hence it is fast, but provides a very easy interface, which enables external applications to use this facility in a convenient way. Different analyses, which require realistic Monte Carlo simulations of observed sky maps, but not TODs, can benefit from this method.

In scanning experiments which target static objects (not transient sources), the data products are better understood by studying the effective beam rather than the detector intrinsic beam function. For instance, consider an elliptic beam function that moves on the sky in a scanning pattern similar to that of Planck. If the scan is along the minor axis of the ellipse, repeated hits in a pixel will tend to make the effective beam more symmetric than the true beam, hence the final map will have less severe high-resolution systematics. On the other hand, if the scan is along the major axis, the elongation will be enhanced in the effective beam and the final map will have more high-resolution systematics. Effective beams thus provide a more intuitive picture of the systematic effects in the observed map as compared to the detector intrinsic beam function. This feature would be useful not only to understand the non-trivial effects of beams on observed CMB anisotropy, but also for precise measurement and extraction of point sources and Sunyaev–Zel'dovich (SZ) clusters from Planck data (The Planck Collaboration 2006) using the accurate knowledge of the point-spread function (PSF) obtained from the effective beams.

The paper is organized as follows. In Section 2, we introduce effective beams. Corresponding mathematics for temperature-only and polarized effective beams are given in Section 3. Section 4 analyzes the accuracy of the method comparing with theoretically exact results. Section 5 describes the details of computational implementation, which is followed by two types of applications of the effective beams. Simulation of maps for Planck and comparison with existing simulations are shown in Section 6. Study of collective effect of the beams via elliptical Gaussian fits and estimation of approximate scalar transfer function for pseudo-C_l estimator by direct spherical harmonic decomposition of beams are covered in Section 7. We conclude with summary and discussions in Section 8.

The finite angular response of a detector is described by the beam function, which is closely related to the PSF in image space. The effective beam is the average of all beam functions pointing at a certain direction for a given scan strategy. Since, in practice, a scan strategy may not pass through exactly the same point twice, the effective beam concept makes sense only when the sky is discretized, then one can define an effective beam for each pixel on the sky.

The mapmaking algorithms (Ashdown et al. 2007; Hinshaw et al. 2007) commonly used in CMB analyses do not take into account the finite size of the instrumental beams, so the resulting map carries the signature of the instrumental beam function and the satellite scan strategy. In fact, the observed map is a convolution of the true CMB sky with the effective beams. Effective beams can therefore be treated as the primary object for studying the combined systematic effects of beams and scan strategy. To appreciate this fact consider the hypothetical case of symmetric detector intrinsic beam. Since the beam center (pointing direction) need not coincide with the pixel center, if the distribution of pointing directions in a pixel is elongated (say only one set of scan rings passes though that pixel), effective beams can become asymmetric.

If the effective beams were symmetric, high-resolution CMB data analysis would be computationally trivial. However, in practice such an advantage is absent. The effective beams in CMB experiments are (significantly) asymmetric due to several unavoidable reasons, which include

• 1.  
  intrinsic optics of the instruments,
• 2.  
  not all detectors can be placed close to the principal optic axis, so there are aberrations,
• 3.  
  a sample is an integral along the scan path, and
• 4.  
  pointing directions are not uniformly distributed in a pixel.

Moreover, the shape of effective beams varies across the sky. Our aim is to estimate the effect of asymmetric effective beams and account for it in data analysis.

One of the immediate advantages of effective beams is that the PSFs, the transpose of effective beams, become easily derivable at every pixel, which is otherwise non-trivial. For instance, in order to compute PSF at neighboring pixels using point-source convolution methods, convolutions must be performed separately for each pixel to prevent overlap between PSFs. Limitation of the effective beams' method is that the PSFs are centered at the pixel centers, though it may not be an issue in (pixelized) map-based applications.

As mentioned before, Monte Carlo simulations are extremely important for studying systematic effects. Consider a sky map observed with certain detector response and scanning. To study the systematics in this map, one would simulate a number of model input maps and convolve them with the same detector response and scanning to estimate the average effect on these simulated maps. This can then be used to draw inferences on the map observed by the instrument. Since at each simulation scanning and beam are the same, huge amount of computation can be saved by pre-computing the effective beams, which include the effect of beam and scan once and for all for the whole set of simulations.

Once the effective beam is computed, it is rather straightforward and notably fast to make simulated "observed" maps. This is indeed the main goal of this paper—to show that the effective beam convolution is fast and accurate. In the end, we also provide some insights on the observed systematic effects directly by studying the properties of the effective beams.

The cost of effective beam computation scales as N_tod × N_beam, where N_tod is the total number of samples, and N_beam is the number of pixels required to adequately describe the beam, and cost of convolution scales as N_pix × N_beam, where N_pix is the total number of equal area pixels in the map. In general, N_beam = N_pix, as the far side lobe of the detector response can span the whole sky. In practice, however, the beam is negligible at most of the pixels; in fact, for CMB experiments of interest, non-negligible part of a main beam is very small and localized,¹⁰ leading to a big saving in computation and input–output (I/O) handling.

The finite size of the instrumental beam response causes exponential damping of angular power spectra at high multipoles, where noise becomes dominant. As a result, a CMB map observed with a certain full width at half-maximum (FWHM), FWHM_beam, has usable information content typically up to a multipole of ℓ_max ∼ 4/FWHM_beam. Also, a CMB map is pixelated in such a way that the number of pixels is minimized while retaining all non-negligible information. According to the HEALPix guidelines, this can be achieved by choosing N_pix ∼ 200/FWHM²_beam. HEALPix divides the sky into N_pix = 12 n²_side pixels, where n_side is an even power of 2 (Górski et al. 2005). Thus the above choice of pixelization translates to n_side ∼ ℓ_max = 4/FWHM_beam. As an example, a map observed with ∼30' beam is usually tessellated with N_pix = 3, 145, 728 HEALPix pixels (i.e., n_side = 512). Then, if we assume that the beam is significant only up to an axial distance of r times FWHM, the number of pixels required to describe the beam would be

$$\begin{eqnarray} && \! \! \! N_{\rm beam}\ \approx \ \frac{N_{\rm pix}}{4\pi } \times {\rm area} \end{eqnarray} \tag{ 2.1 }$$

$$\begin{eqnarray} &&\qquad \qquad \qquad \qquad \ \, \sim \frac{200}{4\pi \, {\rm FWHM}_{\rm beam}^2} \pi \, (r \, {\rm FWHM}_{\rm beam})^2 \nonumber \\ && \qquad \qquad \qquad \qquad \ \, \sim \ 50 \, r^2, \end{eqnarray} \tag{ 2.2 }$$

which means that, if the beam is significant up to a diameter of four times FWHM, the number of pixels needed to describe it with adequate accuracy is ∼200. Note that, a symmetric Gaussian beam falls to 2⁻¹⁶ ≈ 0.0015% (−48 dB) of its peak value at this diameter, and, in Section 4.1, we explicitly show that the net power (effective area) outside this diameter is ∼0.001% of the total beam power.

We compare the computational efficiency of pixel space convolution with the fast harmonic space convolution algorithm "Total Convolution" (Wandelt & Górski 2001; Challinor et al. 2000). Total Convolution uses a fast Fourier transform on the sphere and convolution in the conjugate space to generate a "Total Convolution Data Cube" (TCDC)—sets of convolved sky maps for an equispaced grid of position and orientation coordinates of the beam. Total Convolution is theoretically exact at those grid points for band limited beam functions. Even though Total Convolution itself does not incorporate satellite scan strategy, the algorithm sets the limit on computation cost and accuracy that can be achieved by harmonic space simulation methods. Below we provide approximate theoretically estimated computation cost scaling of effective beam and total convolution; later we provide numerical comparison of accuracy. Note, however, that effective beam method also includes scan strategy (in the precomputation stage), to do the same in harmonic space one goes through the whole scan strategy and interpolates a TCDC, adding a computation cost which scales poorly due to machine I/O.

Computation cost of Total Convolution scales as ∼ℓ³_max m_max (Wandelt & Górski 2001; Challinor et al. 2000; Prezeau & Reinecke 2010), where ℓ_max is the highest multipole moment that contains non-negligible information about the anisotropies in the map and m_max is the highest azimuthal multipole needed to adequately describe beam asymmetry. In order to maintain minimum accuracy in convolution, one needs to use $\ell _{\rm max} \sim \sqrt{\vphantom{A^R}\smash{\hbox{${N_{\rm pix}/3}$}}} \ (= 2 \, n_{\rm side})$ and m_max ≳ 10, which means that the computational cost would scale as ${\gtrsim} N_{\rm pix}^{3/2}/\sqrt{3}$. Since convolution in pixel space scales as ∼N_pix N_beam, pixel based convolution would be faster than total convolution as long as N_pix > 3 × N²_beam ∼ 10⁵, as N_beam was shown to be few hundred for a reasonable beam cutoff radius.

As effective beam computation and convolution are the focus of the paper, the rest of the paper provides details of computation and accuracy analysis considering the main beam only. In addition, though our method and implementation can handle any arbitrary beam shape, in order to reliably compare our results with existing Planck simulations we consider only elliptic Gaussian beam for numerical results.

Computation of effective beams involves large number of mathematical operations. Hence, it is important to formulate it in a computationally minimalist manner. In this section, we list the formulae we use in our code. Before we proceed further we define some quantities which will be used throughout the paper:

• 1.  
  Effective beam for pointing pixel i and beam pixel j: $\mathcal {B}_{ij} \equiv \bm {\mathcal {B}}$.
• 2.  
  Pointing matrix for pixel p and sample (or time) t:

  $$\begin{eqnarray} \mathbf {A} \equiv \ A_{tp} = \left\lbrace \begin{array}{@{}cl}1, & \hbox{if sample lies in pixel}\\ 0, & \hbox{otherwise} \end{array}\right.\!. \end{eqnarray} \tag{ 3.1 }$$
• 3.  
  Measured signal or TOD at sample t: d_t ≡ d.
• 4.  
  Pointing direction at a sample t: $\hat{\mathbf {r}}(t) \equiv [\theta (t), \phi (t)]$.
• 5.  
  Polarization angle, the angle between the polarization axis of the detector and the local meridian: ψ_t ≡ ψ(t).(For unpolarized detectors, the polarization axis is the x-axis of the fiducial detector frame.)
• 6.  
  The set of three pointing angles that define the instantaneous geometry of each detector $\hat{\mathbf {p}}_t \equiv [\theta (t),\phi (t), \psi (t)].$
• 7.  
  Angle between polarization axis and the x-axis of the fiducial detector frame: ψ_pol.
• 8.  
  Angle between the major axis of an elliptical beam (or any preferred axis relative to which the beam profile has been defined) and the x-axis of the fiducial detector frame: ψ_ell.
• 9.  
  Detector intrinsic beam function at a direction $\hat{\mathbf {r}}_b \equiv [\theta _b, \phi _b]$ with the pointing angles $\hat{\mathbf {p}}$: $b(\hat{\mathbf {r}}_b, \hat{\mathbf {p}})$.
• 10.  
  Direction vector corresponding to pixel number p: $\hat{\mathbf {n}}_p$.

Note that the polarization angle is being used to describe the orientation of the detector without any loss of generality. Since the polarization axis of the detector is fixed to the fiducial detector frame, any other orientation angle in this frame, e.g., the orientation angle of the major axis of an elliptic beam, ψ_ell, can be derived by adjusting the appropriate offset angles. For an unpolarized detector, the angle between the polarization axis and the fiducial detector frame is taken as zero.

The detector intrinsic beam function is generally given in a two-dimensional coordinate system on the beam plane, with the origin at the beam center and the x-axis conveniently chosen—e.g., for an elliptic Gaussian beam the major axis is chosen as the x-axis. Here we consider the beam to be defined on the tangent plane, in general the beam could as well be specified in the "curved" polar coordinates at the north pole, then we could take the ϕ = 0 direction as the x-axis.

3.1. Temperature-only Effective Beam

By definition, effective beams are the objects whose convolution with the true CMB sky should produce the observed sky map. So we write the expression for the observed CMB map (without noise) in a convolution form and then define effective beams.

Each observed sample is a convolution of beam and the CMB temperature anisotropy field $T(\hat{\mathbf {r}})$,

$$\begin{eqnarray} d_t &=& \int \mathrm{d}\Omega _{\hat{\mathbf {r}}_b} \, b(\hat{\mathbf {r}}_b,\hat{\mathbf {p}}_t) \, T(\hat{\mathbf {r}}_b) \end{eqnarray} \tag{ 3.2 }$$

$$\begin{eqnarray} &\approx & \ \Delta \Omega \sum _j b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t) \, T_j, \end{eqnarray} \tag{ 3.3 }$$

where T ≡ T_p is the temperature map in pixel basis. Typical mapmakers, however, ignore the effect of extended beam, as it is computationally prohibitive, and it "bins" (averages) the observed data at each pixel (in the absence of correlated noise). So the binned map is made assuming an infinitely narrow (Dirac δ-function) beam, in the discrete pixel space, which is equivalent of setting the beam response to 1 in the pointing pixel and 0 outside (that is, a Kronecker δ). Then the observed binned map $\widetilde{\mathbf {T}} \equiv \widetilde{T}_p$ can be written as

$$\begin{equation} \widetilde{T}_i \ :=\ \frac{\sum _{t} A_{ti} \, d_t}{\sum _{t} A_{ti}}. \end{equation} \tag{ 3.4 }$$

Here the numerator adds all the detected signal if the pointing direction falls in pixel i, and the denominator counts the number of observations in that pixel. However, since the beam is extended, the observed map is connected to the true underlying anisotropy map through the following relation (see Equation (3.3)):

$$\begin{eqnarray} \widetilde{T}_i &:=& \frac{\sum _{t} A_{ti} \, d_t}{\sum _{t} A_{ti}} = \frac{\sum _{t} A_{ti} \, \Delta \Omega \sum _j b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t) \, T_j}{\sum _{t} A_{it}} \end{eqnarray} \tag{ 3.5 }$$

$$\begin{eqnarray} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\! &=& \Delta \Omega \sum _j \left[ \frac{\sum _{t} A_{ti} \, b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t)}{\sum _{t} A_{ti}} \right] \, T_j. \end{eqnarray} \tag{ 3.6 }$$

Then, in order to satisfy the definition of effective beam,

$$\begin{equation} \widetilde{\mathbf {T}} \ =: \ \Delta \Omega \, \bm {\mathcal {B}} \cdot \mathbf {T}, \end{equation} \tag{ 3.7 }$$

one must choose

$$\begin{equation} \mathcal {B}_{ij} \ = \ \frac{\sum _t A_{ti} \, b(\hat{\mathbf {r}}_j, \hat{\mathbf {p}}_t)}{\sum _t A_{ti}}. \end{equation} \tag{ 3.8 }$$

This is indeed the crude definition of effective beam stated in Section 2, which directly applies to the temperature-only case. In fact, the reason for arriving at the expression above in a formal way is to aid our understanding for the case of polarized detectors, where the definition is not straightforward.

Equation (3.8) is essentially what we compute for temperature only beams, the detailed method of how it is done in an efficient way is discussed in this section and Section 5.

It is clear from Equation (3.8) that for each sample t, the beam function b must be evaluated at all the beam pixels j by rotating the detector intrinsic beam function. We do this by transforming the pixel center coordinates to the coordinates in which the detector intrinsic beam function is specified. Note that these operations are the most expensive mathematical operations in the effective beam computation, as they are used about trillions times (N_tod × N_beam), so we take good care in writing them in efficient forms.

• 1.  
  First, we obtain the coordinates of the point in the (locally flat) orthogonal coordinates defined by the polar and azimuthal unit vectors ($\hat{\bm {\theta }},\hat{\bm {\phi }}$) with the beam center ($\hat{\mathbf {r}}$) as the origin.The radial, polar, and azimuthal unit vectors in spherical polar coordinates are, respectively, given by

  $$\begin{eqnarray} \hat{\mathbf {r}}&=& \sin \theta \cos \phi \, \hat{\mathbf {\vphantom{a^a}i}}\, + \, \sin \theta \sin \phi \, \hat{\mathbf {j}}\, + \, \cos \theta \hat{\mathbf {k}}\nonumber \\ \hat{\bm {\theta }}&=& \cos \theta \cos \phi \, \hat{\mathbf {\vphantom{a^a}i}}\, + \, \cos \theta \sin \phi \, \hat{\mathbf {j}}\, - \, \sin \theta \hat{\mathbf {k}}\nonumber \\ \hat{\bm {\phi }}&=& -\!\sin \phi \hat{\mathbf {\vphantom{a^a}i}}\, + \, \cos \phi \hat{\mathbf {j}}. \end{eqnarray} \tag{ 3.9 }$$

  Then, with respect to the beam center, the coordinates (Δx, Δy) of a pixel at $\hat{\mathbf {r}}_b$ can be written as

  $$\begin{eqnarray} \hspace{-12pc}\Delta x &:=& (\hat{\mathbf {r}}_b - \hat{\mathbf {r}}) \cdot \hat{\bm {\theta }} \end{eqnarray} \tag{ 3.10 }$$

  $$\begin{eqnarray} \!\!\!\!\! \!\!&=& \sin \theta _b \, \cos \theta \, \cos (\phi _b-\phi) - \ \sin \theta \, \cos \theta _b \end{eqnarray} \tag{ 3.11 }$$

  $$\begin{eqnarray} \hspace{-12pc}\Delta y &:=& (\hat{\mathbf {r}}_b - \hat{\mathbf {r}}) \cdot \hat{\bm {\phi }} \end{eqnarray} \tag{ 3.12 }$$

  $$\begin{eqnarray} \hspace{-8pc}&=& \sin \theta _b \, \sin (\phi _b - \phi). \end{eqnarray} \tag{ 3.13 }$$

  We also provide a less accurate but faster (by ∼10% of the whole computation) set of formulae to get these coordinates. We have tested that these formulae provide enough accuracy for practical purposes:

  $${\fontsize{8.5}{11}\selectfont{\begin{eqnarray} \hspace{-4.4pc}\Delta x &=& - \!\theta + \left\lbrace \begin{array}{@{}l}\theta _b \, \cos (\phi _b - \phi), \quad {\rm if} \ \theta \le \pi /2 \\ \pi - (\pi -\theta _b) \, \cos (\phi _b - \phi), {\rm otherwise} \end{array} \right. \end{eqnarray}}} \tag{ 3.14 }$$

  $${\fontsize{8}{11}\selectfont{\begin{eqnarray} \Delta y &=& \sin \theta \times \left\lbrace \begin{array}{@{}l}\sin (\phi _b - \phi), {\rm if} \ |\phi _b - \phi | > 1^\circ \ \ {\rm and}\ |\Delta x| < {\rm FWHM}_{\rm beam}\\ \phi _b - \phi, \quad {\rm otherwise}. \end{array}\right.\nonumber\\ \end{eqnarray}}} \tag{ 3.15 }$$
• 2.  
  We then rotate the coordinates by the angle ψ_b = ψ(t) − ψ_pol + ψ_ell to transform to the beam coordinate system:

  $$\begin{eqnarray} &&\left(\begin{array}{@{}c@{}}\Delta x_{\rm beam}\\ \Delta y_{\rm beam} \end{array}\right) = \left[\begin{array}{@{}rr@{}}\cos \psi _b & \sin \psi _b\\ -\sin \psi _b & \cos \psi _b \end{array}\right] \left(\begin{array}{@{}c@{}}\Delta x\\ \Delta y \end{array}\right). \end{eqnarray} \tag{ 3.16 }$$
• 3.  
  The beam function is a known input in the new Δx_beam, Δy_beam coordinates. For example, if the intrinsic beam is elliptical Gaussian with major and minor axes σ_a, σ_b, respectively, the value of the rotated beam function at $\hat{\mathbf {r}}_b$ will be

  $$\begin{equation*} \frac{1}{2 \pi \sigma _a \sigma _b} \exp \left[-\frac{\Delta x_{\rm beam}^2}{2 \, \sigma _a^2} - \frac{\Delta y_{\rm beam}^2}{2 \, \sigma _b^2}\right]. \end{equation*}$$

In this paper, all the numerical exercises have been done assuming elliptical Gaussian beam shape with parameters read from a database of Planck focal plane parameters and the accurate formulae for beam coordinates, Equations (3.11) and (3.13), have been used.

3.2. Polarized Effective Beams

To define effective beams for polarized detectors we again formally express the output of a mapmaker in the noiseless case. This expression is then generalized to the extended beam case and the polarized effective beam is defined. In case the beam is significantly extended, we also need to consider the coordinate changes within the beam (as the local meridians—the polarization reference axes—at different pixels may become significantly different). We supply the correction formulae in the end of this subsection, though for Planck main beams we achieved adequate accuracy without using this correction.

3.2.1. Polarized Binned Map

The binning process is independent of the beam size. However, if the beam has infinite resolution (δ-function) the binned map should be identical to the true map in the noiseless case. So we derive the formula for making binned polarized maps starting with such an idealized case, with the view that the procedure will be valid for any beam.

A detector in a CMB experiment records only one number for each sample, a linear combination of three Stokes parameters $I(\hat{\mathbf {r}}), Q(\hat{\mathbf {r}}), U(\hat{\mathbf {r}})$. For the idealized detector, one can write

$$\begin{eqnarray} d_t &=& w_{1}(\psi _t) \, I(\hat{\mathbf {r}}(t)) + w_{2}(\psi _t) \, Q(\hat{\mathbf {r}}(t)) + w_{3}(\psi _t) \, U(\hat{\mathbf {r}}(t)),\nonumber \\ \end{eqnarray} \tag{ 3.17 }$$

where w_t ≡ [w₁(ψ_t), w₂(ψ_t), w₃(ψ_t)] is the weight vector at sample t. The above equation can be written in a discrete form using the pointing matrix A_tp as

$$\begin{equation} d_t = \sum _{p} A_{tp} \mathbf {w}^T_t \left(\begin{array}{@{}c@{}}I\\ Q\\ U\end{array}\right)_{p}. \end{equation} \tag{ 3.18 }$$

For the polarization sensitive instruments used in the Planck satellite, it is a good approximation to write the weight factors as

$$\begin{equation} \mathbf {w}_t := \left(\begin{array}{@{}c@{}}1 \\ \gamma \, \cos (2 \psi _t) \\ \gamma \, \sin (2 \psi _t) \end{array}\right), \end{equation} \tag{ 3.19 }$$

where γ := (1 − )/(1 + ) is the polarization efficiency, with being the cross-polar leakage.

Since the detector provides one number for each sample, but one needs to estimate three numbers (I, Q, U) for each pixel, there must be at least three observations per pixel and the set of linear equations must be solved. Ignoring the issues of degeneracy (not enough independent detector configurations), below we describe how one can solve the equations to make a binned map.

The goal of making a binned map is that, if the beam was indeed very narrow and there was no noise, one should get the true sky map as the result. We repeat the procedure presented in Jones et al. (2007). Multiply Equation (3.18) through by $\mathbf w_t$ and sum over the samples in a pixel (enumerated by the pointing matrix A), to obtain

$$\begin{equation} \sum _t A_{tp} \mathbf {w}_t d_t = \sum _t \sum _{p^{\prime }} A_{tp} A_{tp^{\prime }} \mathbf {w}_t \mathbf {w}^T_t \left(\begin{array}{@{}c@{}}I\\ Q\\ U\end{array}\right)_{p^{\prime }}. \end{equation} \tag{ 3.20 }$$

Because the pointing matrix has only one unit entry per sample, the sum over p' collapses to

$$\begin{equation} \sum _t A_{tp} \mathbf {w}_t d_t = \mathbf {H}_p \left(\begin{array}{@{}c@{}}I\\ Q\\ U\end{array}\right)_{p}. \end{equation} \tag{ 3.21 }$$

On the left-hand side we have binned weighted TOD and on the right we have the true sky times the 3 × 3 polarized "number of observation matrix" at pixel p

$$\begin{equation} \mathbf {H}_p \ := \ \sum _t A_{tp} \mathbf {w}_t \mathbf {w}^T_t. \end{equation} \tag{ 3.22 }$$

To recover the sky, multiply the binned weighted TOD by the inverse of number of observation matrix. So the polarized binned map is given by

$$\begin{eqnarray} &&\left(\widetilde{I} \quad \widetilde{Q} \quad \widetilde{U} \right)_p^T \ = \ \mathbf {H}_p^{-1} \, \sum _t A_{tp} \mathbf {w}_t d_t \nonumber \\ && \qquad = \ \left[ \sum _t A_{tp} \mathbf {w}_t \mathbf {w}^T_t \right]^{-1} \sum _t A_{tp} \mathbf {w}_t d_t. \end{eqnarray} \tag{ 3.23 }$$

As mentioned above, this is the general binning formula for any beam size. Also note that, not surprisingly, it is similar to the formula for the temperature only case (Equation (3.5)).

3.2.2. Polarized Effective Beams

We can now calculate the effect of binning on a data set observed with a realistic beam of finite size and thereby arrive at the formula for polarized effective beams. We follow the same path as we did for temperature only beams.

Following Jones et al. (2007), the observed signal at sample t is simply the integral of the contraction:

$$\begin{eqnarray} d_t &=& \int \mathrm{d}\Omega _{\hat{\mathbf {r}}_b} \, b(\hat{\mathbf {r}}_b,\hat{\mathbf {p}}_t) \, \mathbf {W}^T(\hat{\mathbf {r}}_b,\hat{\mathbf {p}}_t) \cdot (I \quad Q \quad U)^T_p \qquad \end{eqnarray} \tag{ 3.24 }$$

$$\begin{eqnarray} \ \ \ \ \ &\approx & \Delta \Omega \sum _{p} b(\hat{\mathbf {n}}_p,\hat{\mathbf {p}}_t) \, \mathbf {W}^T(\hat{\mathbf {n}}_p, \hat{\mathbf {p}}_t) \cdot (I \quad Q \quad U)_{p}^T,\qquad \end{eqnarray} \tag{ 3.25 }$$

where $\mathbf {W}(\hat{\mathbf {r}}_b,\hat{\mathbf {p}}_t)$ is the polarization weight vector for general extended beam,

$$\begin{equation} \mathbf {W}(\hat{\mathbf {r}}_b, \hat{\mathbf {p}}_t) := \left(\begin{array}{@{}c@{}}1 \\ \gamma \, \mathcal {P}(\hat{\mathbf {r}}_b,\hat{\mathbf {p}}_t) \cos (2 \psi _t) \\ \gamma \, \mathcal {P}(\hat{\mathbf {r}}_b,\hat{\mathbf {p}}_t) \sin (2 \psi _t) \end{array}\right), \end{equation} \tag{ 3.26 }$$

where $\mathcal {P}(\hat{\mathbf {r}}_b,\hat{\mathbf {p}})$ is the normalized beam response function defined as

$$\begin{equation} \mathcal {P}(\hat{\mathbf {r}}_b,\mathbf {0}) \ := \ \frac{P_\parallel (\hat{\mathbf {r}}_b) - P_\perp (\hat{\mathbf {r}}_b)}{P_\parallel (\hat{\mathbf {r}}_b) + P_\perp (\hat{\mathbf {r}}_b)}, \end{equation} \tag{ 3.27 }$$

expressed in terms of the co- and cross-polar beam response functions $P_\parallel (\hat{\mathbf {r}}_b)$ and $P_\perp (\hat{\mathbf {r}}_b),$ respectively. $\mathcal {P}(\hat{\mathbf {r}}_b,\hat{\mathbf {p}})$ can be obtained for any given pointing $\hat{\mathbf {p}}$ following the same procedure as that for the beam profile $b(\hat{\mathbf {r}}_b,\hat{\mathbf {p}})$ described explicitly for temperature only beams.

On substitution into the mapmaking equation (Equation (3.23)), the observed (convolved) sky map becomes

$$\begin{eqnarray} &&\hspace{-3pc}\left(\begin{array}{@{}c@{}}\widetilde{I}\\[3pt] \widetilde{Q}\\[3pt] \widetilde{U}\end{array}\right)_i \ = \ \mathbf {H}_i^{-1} \, \sum _t A_{ti} \, \mathbf {w}_t \nonumber\\ && \hspace{-3pc}\ \ \ \quad \qquad \times \Delta \Omega \sum _j b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t) \, \mathbf {W}^T(\hat{\mathbf {n}}_j,\hat{\mathbf {p}}_t) \cdot \left(\begin{array}{@{}c@{}}I\\ Q\\ U\end{array}\right)_j. \end{eqnarray} \tag{ 3.28 }$$

Thus, in order to satisfy the definition of polarized effective beam,

$$\begin{equation} (\widetilde{I} \quad \widetilde{Q} \quad \widetilde{U})_i^T \ =: \ \Delta \Omega \sum _j \bm {\mathcal {B}}_{ij} \cdot (I \quad Q \quad U)_j^T, \end{equation} \tag{ 3.29 }$$

one must choose

$$\begin{eqnarray} \bm {\mathcal {B}}_{ij} &=& \left[ \sum _t A_{tp} \mathbf {w}_t \mathbf {w}^T_t \right]^{-1} \ \sum _t A_{ti} \, b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t) \, \mathbf {w}_t \mathbf {W}^T(\hat{\mathbf {n}}_j,\hat{\mathbf {p}}_t). \nonumber \\ \end{eqnarray} \tag{ 3.30 }$$

The effective beam for each pointing pixel i now has two parts. There is a 3 × 3 symmetric matrix H_i, which is a generalized version of number of observation at that pixel. Then, for each neighboring pixel j there is a 3 × 3 matrix $\sum _t A_{ti} \, b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t) \, \mathbf {w}_t \mathbf {W}^T(\hat{\mathbf {n}}_j,\hat{\mathbf {p}}_t)$, which encapsulates the beam profile and need not be a symmetric matrix in general. Computation of this quantity can be done in a similar way as in the temperature only case, first computing the coordinates of neighboring pixels with respect to the beam center and then evaluating the polarized beam profiles $b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t) \mathbf {W}^T(\hat{\mathbf {n}}_j,\hat{\mathbf {p}}_t)$ using the measured/modeled intrinsic detector response functions.

The cross polar beam is normally much smaller than the co-polar beam and, in practice, the co-polar beam might be the only response that could be accurately measured. These facts lead to a reasonable assumption that $\mathcal {P}(\hat{\mathbf {r}}_b) \approx 1$, so that $\mathbf {W}(\hat{\mathbf {r}}_b,\hat{\mathbf {p}}_t) = \mathbf {w}_t$, which, in turn, simplifies the computation of effective beams,

$$\begin{equation} \bm {\mathcal {B}}_{ij} = \left[ \sum _t A_{tp} \mathbf {w}_t \mathbf {w}^T_t \right]^{-1} \sum _t A_{ti} \, b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t) \, \mathbf {w}_t \mathbf {w}^T_t, \end{equation} \tag{ 3.31 }$$

as the beam profile for all the polarization becomes the same, except for weight factors (not dependent on $\hat{\mathbf {r}}_j$), and makes all the component matrices symmetric. Note that, if the effective beam matrices are not symmetric,¹¹ computational storage requirement increases by 50% (nine instead of six components for each matrix). Computation time is not much affected by this modification.

FEBeCoP does include the ability to compute effective beam matrices of both kinds, symmetric or otherwise. However, the numerical results and comparisons presented in this paper are for the possibly more practically achievable scenario of $\mathcal {P}(\hat{\mathbf {r}}_b) = 1$. So, unless otherwise specified, the matrices representing polarized effective beams computed for accuracy and performance analysis of this method would be symmetric.

3.2.3. Transformation from Beam Frame to Sky Frame

If the beam is significantly extended, the pixels that are far (that is, on very different longitudes) would suffer from an extra polarization rotation effect, this is because the far-away longitudes are no longer close to parallel (except when they are very close to the equator), so the frame of polarization also rotates. The mathematics that account for that extra rotation is presented below.

If the detector orientation at sample t is represented by Euler rotation matrix D(ϕ_t, θ_t, ψ_t), then the polarization field for the beam in pixel j is evaluated at a direction $\mathbf {D}^{-1}(\phi _t,\theta _t,\psi _t) \cdot \hat{\mathbf {n}}_j$ in the beam frame. The orthonormal polarization bases differ between frames by a rotation. Stoke I is a scalar, unchanged by the rotation. To transform (Q, U), we need to know the angle between the polarization bases in the sky and beam frame (Figure 3). Imagining the spherical triangle between the north pole, the position of sample i, and pixel p, we may deduce the angle between the orthonormal bases at the pixel location, denoted by ζ_ip, with two applications of the law of haversines,

$$\begin{eqnarray} \hspace{-3.6pc}{\rm haversin} (\Delta \Phi _{ip}) & = & {\rm haversin} (\theta _i -\theta _p) + \sin (\theta _i) \sin (\theta _p) \nonumber \\ \hspace{-3.6pc}&& \times{\rm haversin} (\phi _i -\phi _p) \end{eqnarray} \tag{ 3.32 }$$

$$\begin{eqnarray} \hspace{-2.5pc}{\rm haversin} (\theta _i) & =& {\rm haversin} (\theta _p - \Delta \Phi _{ip}) + \sin (\theta _p) \sin (\Delta \Phi _{ip}) \nonumber \\ \hspace{-2.5pc}&& \times {\rm haversin} (\zeta _{ip}), \end{eqnarray} \tag{ 3.33 }$$

where haversin(θ) = sin²(θ/2). The first equation solves for ΔΦ_ip, representing the great circle distance from the pixel to the sample position, which is required in the second equations to solve for ζ_ip.

Zoom In Zoom Out Reset image size

This extra rotation would lead to a transformation of the observed sky in the following way:

$$\begin{equation} \left(\begin{array}{@{}c@{}}Q^{\prime } \\ U^{\prime }\end{array}\right) \ = \ \left[\begin{array}{@{}rr@{}}\cos 2\zeta _{tp} & \sin 2\zeta _{tp} \\ -\!\sin 2\zeta _{tp} & \cos 2\zeta _{tp} \end{array}\right] \left(\begin{array}{@{}c@{}}Q \\ U\end{array}\right). \end{equation} \tag{ 3.34 }$$

This would change the beam weight factor in the expression for effective beams to $\mathbf {W}(\hat{\mathbf {n}}_p,\hat{\mathbf {p}}_t) \rightarrow \mathbf {W}^{\prime }(\hat{\mathbf {n}}_p, \hat{\mathbf {p}}_t) = \mathbf {W}(\psi _t + \zeta _{tp},\hat{\mathbf {p}}_t)$ and the expression would become

$$\begin{eqnarray} \bm {\mathcal {B}}_{ij} &=& \left[ \sum _t A_{tp} \mathbf {w}_t \mathbf {w}^T_t \right]^{-1} \ \sum _t A_{ti} \, b(\hat{\mathbf {r}}_j,\hat{\mathbf {p}}_t) \,\mathbf {w}_t \mathbf {W}^{\prime }(\hat{\mathbf {n}}_j, \hat{\mathbf {p}}_t). \nonumber\\ \end{eqnarray} \tag{ 3.35 }$$

It is easy to see that the effective beam matrix at each pixel is no longer symmetric in this case, even if the cross-polar beam is ignored.

Effective beams defined over a grid could have quadrature accuracy issues if the pixelization was not fine enough. Also, the radial cutoff we impose on the beam leads to some information loss. In this section we study the truncation error in integrating a two-dimensional Gaussian over a HEALPix grid to determine the minimum value of the cutoff radius and also show that quadrature error is sub-dominant here. Then, by comparing effective beam convolved CMB maps using Total Convolution for a toy scan model, we show that the default map pixel resolution (discussed in Section 2) provides sufficient accuracy.

4.1. Effective Area Loss of a Truncated Gaussian

By definition, the surface integral of a two-dimensional Gaussian over an infinite plane is unity. However, if we cut off a zero mean Gaussian with variance σ² at $r_{\rm max} = n_{\rm FWHM} \, \times \, {\rm FWHM} \ = \ 2 \, n_{\rm FWHM} \,\sqrt{2\ln {2}} \, \sigma$, the "area loss" turns out to be

$$\begin{eqnarray} &&\epsilon := 1 - \frac{1}{\sigma ^2} \int _0^{r_{\rm max}} d r \ r \, e^{-\frac{r^2}{2\sigma ^2}} \ = \ e^{-\frac{r_{\rm max}^2}{2\sigma ^2}} = 2^{- 4 \, n_{\rm FWHM}^2}.\qquad \end{eqnarray} \tag{ 4.1 }$$

The order of magnitude error is thus

$$\begin{equation} \log _{10}{\epsilon } \ = \ 4 \, n_{\rm FWHM}^2 \, \log _{10}{2} \ \approx \ 1.2 \, n_{\rm FWHM}^2. \end{equation} \tag{ 4.2 }$$

That is, if we choose the cutoff at n_FWHM = 2, the area loss will be ∼10⁻⁵, which should provide enough accuracy for convolution. It is worth noting that in the context of CMB maps the error is actually the difference between convolved values with and without cutoff, which is dependent on the underlying CMB map. Since the maps have zero expectation, the expected error would also be zero, but the fractional rms error would be the same as the effective area loss.

In practice convolution is done in discrete space. To confirm that the above analytical estimate still holds and the usual choice of pixel resolution does not add any significant quadrature error, we do the exercise numerically. We compute the area loss for a Gaussian with FWHM ∼32' (Planck LFI 30 GHz beam size) using two sets of beam cutoff radius, n_FWHM = 2 and n_FWHM = 4, and pixelization, N_pix = 12 × 512² and N_pix = 12 × 2048².

Figure 4 shows plots of area loss for the four cases as functions of latitude (the evaluation was done along the meridian without any loss of generality). Clearly, area loss depends on the cutoff radius and cannot be improved by increasing N_pix to higher than the usual value for the given beam size. There is some quadrature error near the poles, but we can ignore that as the solid angle which is affected by this is very small (only few pixels). There is a sharp features near each HEALPix polar circle for higher cutoff radius, only if the target accuracy is better than 10⁻⁵. Our default choice would be n_FWHM = 2, which provides adequate accuracy without bringing in any significant quadrature issues.

Zoom In Zoom Out Reset image size

4.2. Comparison with Total Convolution

Next we compare effective beam convolved maps for toy models of beam orientations with TCDCs for different choices of pixel resolution. TCDCs are theoretically exact for band limited beam functions on an equispaced grid of beam location and orientation. We used those values of orientation in the effective beam computation and convolved with the same map (in pixel basis) from which the TCDCs were derived. TCDCs are computed at equispaced intervals of R.A. and decl., so the beam centers do not coincide with the HEALPix pixel centers. This is similar to the real scans, which need not pass right through the pixel centers. One small approximation is unavoidable here, a Gaussian beam is not purely band limited, specially for asymmetric beam there is a cutoff on the azimuthal index "m." This error is not severe unless the accuracy requirement is very stringent. In fact, this error adds to true differences making the test even more strict for the pixel space convolution.

Figure 5 shows that for a ∼32' beam, the natural choice of pixel resolution n_side = 512 provides very good accuracy and further increase in resolution does not improve it by any significant amount. At low n_side the beam is very poorly sampled, which leads to orders of magnitude better accuracy for each step of n_side improvement. However, beyond n_side = 512 the amount of extra information added is low, so the accuracy improvement is marginal and is also limited by single-precision floating point accuracy. So, we may conclude that the usual choice of pixelization is the optimum choice.

Zoom In Zoom Out Reset image size

In the previous sections, we have presented the theoretical statement of the problem and error analysis. The practical implementation of the method is a bigger challenge due to the organizational complexities of efficient distributed computation. This section outlines the possible routes of computation, their relative merits and demerits, and some final outcomes.

5.1. Quantities to Compute

5.2. Algorithm

5.3. Implementation for Planck

We have fully implemented the effective beam method for the Planck mission. The mission specific inputs required for the effective beam computation are the detector characteristics and the pointing information. Planck's detector characteristics—beam parameters, polarization angles, cross polar leakage, etc.—are read from a focal plane database. The detector pointing information is extracted from compressed satellite pointing using the Generalized Compressed Pointing library through the MADCAP3 interface.¹²

The beams are stored in binary formats which are accessible to both C and Fortran standard I/O routines. Each process writes its own binary file, effective beams and other data for each pixel are stored as one "record" in one of those files. To preserve compatibility between different computing languages, the record size is kept the same for each pixel, which is determined by maximum N_list. The address of a pixel is stored in a combined index file, which is actually a sky map—each map value encodes the file number and location for that pixel as a single precision integer. A software interface to conveniently load partial or full effective beams has also been developed.

The cost of effective beam computation scales as N_tod × N_beam, where N_tod is the total number of samples (sum of number of observations at all pixels). For 15 months of Planck HFI 143 GHz channel (simulated) data, sampled at 172 Hz (total 6.8 × 10⁹ samples), with N_pix = 12 × 2048² pixels, the computation cost is ∼130 CPU hr per detector. For Planck LFI 30 GHz channel, the cost time is ∼30 CPU hours. These costs could be reduced by ∼10% by using the little less accurate formulae for beam computation provided in Section 3.1.

The total size of the effective beam can be estimated as follows: for each pixel with N_list beam pixels, total volume of data is essentially the number of neighboring pixels plus number of beam matrix elements. So for polarized beams, if we use single precision numbers, the size of each record is approximately 4(1 + 6)N_list = 28N_list bytes, so the total size of the polarized effective beam will be 28N_listN_pix bytes. For Planck HFI 143 GHz with n_side = 2048 and N_list ∼ 225 the total volume of effective beam turns out to be ∼290 GB. This is approximately the disk space required to store the beams, the total amount of (distributed) memory requirement is greater than this, but of the same order. Note that, if we chose to use different beams for unpolarized and polarized components, the polarized effective beam matrices will not be symmetric, and 50% more storage would be required. The other factor which may affect the size is the cutoff radius. Increasing this may improve the accuracy slightly for polarized detectors, but the total beam object size scales as the square of the radius, so the increase in this data volume can be significant.

We illustrate the features of effective beams and PSFs derived from them by showing a few figures for Planck 30 GHz (simulated data). The distribution in the number of observations due to the scan pattern is shown in Figure 6. In Figure 7, we show logarithmic plots of effective beams (left) and PSFs (right) at different locations on the sky, where the locations have been marked in Figure 6 (in Galactic coordinates), respectively, by "cross" (ecliptic north pole), "square" (a cusp near ecliptic north pole), "triangle" (a mid-declination point), and "asterisk" (ecliptic equatorial point). We overlay the isocontours at 50%, 10%, 1%, and 0.1% levels for the numerical maps and also elliptical Gaussian fits to the beams and PSFs. The fitted values, FWHM, ellipticity, orientation angle with respect to the local meridian (vertical axis), and the normalized rms fit error have been mentioned in the subtitle. The center of the central pixel (effective beam/PSF center) and the center of the fitted Gaussian have been marked with "+" and "×," respectively, when they coincide the markers form an "eight arm asterisk." These figures tell us that the effective beam/PSF has significant variations across the sky, implying that simplifying assumptions for asymmetric beam studies may not be very accurate. Specially, near the ecliptic poles and the cusps of the scans, the asymmetry can be large.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The main purpose of the effective beams is very fast convolution and the FEBeCoP implementation handles it efficiently. Below first we provide some approximate computation cost estimates for the effective beam convolutions; next we test the accuracy by comparing with the results derived using existing Planck simulation software "level-S" and analyzed by HEALPix utility anafast (Górski et al. 2005).

6.1. Efficiency

6.2. Accuracy: Comparison with Existing Planck Simulations

Next, we assess the accuracy of the convolution by comparing with the simulated maps generated using the standard Planck Level-S TOD simulation pipeline in combination with the Springtide or MADAM map maker (Ashdown et al. 2007). Starting from the focal plane database, the satellite scan strategy and (the spherical harmonic transforms of) an input map, we generate TCDC with m_max = 2 for symmetric beams and m_max = 14 for asymmetric beams. The TCDC are then interpolated by multimod (Reinecke et al. 2006), to generate TOD at each scan sample. These (noiseless) TODs are then binned by a mapmaker to generate maps. The effective beams are computed using the same pointing and focal plane database and convolved with (the pixelized version of) the input map. The comparison between the maps and the power spectra for two Planck channels, 30 GHz (all detectors: 27a, 27b, 28a, 28b) and 143 GHz (four detectors: 1a, 1b, 4a, 4b), are shown below. Among all the Planck beams, the 30 GHz ones have the highest asymmetry (ellipticity ∼1.35–1.40) and the 143 GHz ones have the lowest asymmetry (ellipticity ∼1.02–1.11).

We first show simulated full sky temperature and polarization maps for Planck 30 GHz channel using Level-S and FEBeCoP and their differences in left and right panels of Figure 8. The maps from the two methods are clearly similar and the difference is very small—the fractional rms difference for temperature maps is few × 0.01% and for polarization amplitude the difference is few × 0.1% (note that the contrast in the difference maps have been amplified by using a histogram equalized color map). We have also plotted zoomed versions of the above maps for north polar and equatorial regions in Figure 9. The agreement between FEBeCoP and levelS is good in all the places, expect for only few pixels (among millions) right at the poles, which is not significant. (We omit these plots for 143 GHz, where they are uninformative.) The power spectra of the maps from two methods for both 30 GHz and 143 GHz¹³ channels are shown in Figures 10 and 12, the spectra are almost indistinguishable (the lines nearly coincide). Their fractional differences reveal more details (Figures 11 and 13). What is important to note in these comparisons is that the difference between the effective beam convolved maps and the Level-S plus mapmaker maps, and their power spectra, is much smaller than the corresponding difference between a Level-S and symmetric Gaussian smoothed map, obtained using the HEALPix utility synfast (Górski et al. 2005). Especially for polarization the symmetric beam assumption systematic is large (∼few 10%), while the effective beam error is ≲0.1%. Even for symmetric detector intrinsic beams, small asymmetry is introduced to the effective beams by scan and that too is faithfully reproduced by FEBeCoP. This can be seen in the top panel of Figure 10. Here the spectrum of theoretical symmetric Gaussian convolved map computed using synfast differs from the spectrum of levelS map with symmetric beam and simulated Planck scan, but the spectrum of FEBeCoP map for the same inputs lies exactly on top of the LevelS spectrum. These exercises demonstrate that effective beam convolution provides very good accuracy in mimicking the systematic effects in the detectors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Effective beams are not merely useful for fast convolution, they can provide a more realistic way of understanding and interpreting systematic effects. To demonstrate the potential applications, we provide two examples in this section.

7.1. Statistics of the Effective Beams

Complex combination of scanning, beam and pixel shapes lead to different shapes of effective beams and PSFs, few examples were shown in Figure 7. Software tools developed for FEBeCoP readily allow one to obtain useful statistics of the beams and PSFs. In Figure 14, we show histograms characterizing the variation of elliptical Gaussian fit parameters to the 30 GHz simulated effective beams, sampled at HEALPix pixel centers of n_side = 16 (this method ensures uniform sampling, though it has a very high probability of missing the "cusp" regions in the scan pattern, where the beams/PSFs can be drastically different). Typically the effective beams are few percent bigger than the intrinsic detector beams as their centers are spread over a pixel, adding extra smoothing.

Zoom In Zoom Out Reset image size

7.2. Understanding Observed Power Spectra

In the previous section, we presented comparisons between symmetric and asymmetric beam convolved maps. The observed power spectra of these maps were different by few percent. There is, however, no trivial way of interpreting this difference—the reason is that the effective beams have varying shapes across the sky. Previous attempts to explain this systematic deviation using mathematical models assume constant beam shape (with possibly rotation) across the sky, though the real scan strategy does not have a simple analytical model. This problem can be partially or completely tackled by studying the effective beams. Here we use a leading-order approximation to show that the ratio of temperature power spectra observed using asymmetric and symmetric beams can be modeled using the rotated spherical harmonic transforms of the effective beams.

For a symmetric beam of spatially non-varying size, it is trivial to show that the observed power spectrum $\widetilde{C}_l$ is related to the true power spectra through the relation

$$\begin{equation} \widetilde{C}_l \ = \ B_l^2 \, C_l, \end{equation} \tag{ 7.1 }$$

where B_l are the Legendre transforms of the beam. In general, though, the effective beam shapes are not constant and they vary across the sky, as shown in Figure 7. So the transfer function is no longer as simple as B²_l. The easiest way to estimate the transfer function is by performing Monte Carlo simulations and estimate the ratio $\widetilde{C}_l/C_l$ for each multipole, which is, in fact, one of the motivations for the effective beam approach. However, to extract more insights directly form the effective beams, we crudely estimate average B²_l of effective beams by taking some sample effective beams from different latitude, computing the rotated harmonic transform and incorporating the leading order correction following the algebra in Mitra et al. (2004).

Even though these approximations are crude, the results look quite satisfactory. In Figure 15, we plotted the ratio of power spectra of simulated maps for Planck 143-5 channel with asymmetric and symmetric beams. For the same two cases, the ratio of B²_l estimated from effective beams is overlaid on the plot. Clearly the match is quite good for demonstration purposes. There are accuracy issues at high multipoles due to many reasons—for example, the formula we used is valid for beams with reflection symmetry, though the effective beams do not have that symmetry—the effective beams are constructed from elliptic or circular beams which have that symmetry, but since the beam centers and orientations are not distributed uniformly in a pixel, the reflection symmetry is slightly broken, which may introduce extra first-order correction terms.

Zoom In Zoom Out Reset image size

The ratio of B²_l of effective beams at certain latitudes (not averaged) is also overlaid to show that the characteristics of the effective beams vary significantly at different sky locations, which may break simplified assumptions often used in simulations. The variation of B²_l with latitude need not represent any trend, as the effective beams can vary significantly from pixels to pixels and also Planck scanning is not completely azimuthally symmetric.

Realistic simulations of high-precision experiments play an extremely important role in studying systematic effects of the instrument. There is no exception for scanning CMB experiments, like Planck. In the last several years, there has been significant success in developing an end-to-end (E2E) simulation pipeline to simulate very realistic Planck TOD going through the full scan strategy. However, such simulations involve a very detailed computational procedure and are usually computationally prohibitive in applications to adequate Monte Carlo studies (that is, repetitive many thousand runs, or more). Many CMB analyses do not require full TOD for the study of systematic effects, as simulated realistic observed maps provide all the ingredients required for these studies. In fact, these studies may require many more simulations than it is feasible to obtain by performing the E2E simulations and then going through mapmaking. Beam asymmetry is one of the most important potential instrumental systematic effects in Planck. Pseudo-C_l and other power spectrum estimators could make use of large number of realistic sky map simulations to estimate the correction needed to estimate the true power spectra. This was the motivation for developing a pixel-based convolution method using effective beams, which we have implemented using the computer code FEBeCoP. The method precomputes the effective beams and makes it possible to perform very fast convolutions to simulate realistic sky maps. FEBeCoP also provides PSFs at every pixel on the sky, which is otherwise non-trivial to obtain. We have successfully constrained the computation cost to a moderate value as well as managed to properly utilize the available distributed computing power. The computation requirement for FEBeCoP is modest compared with the modern scientific computing facilities available for Planck analyses. Due to the simplicity of FEBeCoP interface and its efficiency, we believe that it would be very easy and useful to incorporate effective beam convolution in many systematic effect studies.

We thank B. Crill, M. Seiffert, and F. Bouchet for useful discussions on Planck beams. We acknowledge C. Cantalupo and T. Kisner for support with MADcap3 interface, I. O'Dwyer for assistance in devising data formats for storage of effective beams, and R. Keskitalo for aid in usage of MADAM mapmaker. Part of the research described in this paper was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.