Residual noise covariance for Planck low-resolution data analysis

R. Keskitalo; M. A. J. Ashdown; P. Cabella; T. Kisner; T. Poutanen; R. Stompor; J. G. Bartlett; J. Borrill; C. Cantalupo; G. de Gasperis; A. de Rosa; G. de Troia; H. K. Eriksen; F. Finelli; K. M. Górski; A. Gruppuso; E. Hivon; A. Jaffe; E. Keihänen; H. Kurki-Suonio; C. R. Lawrence; P. Natoli; F. Paci; G. Polenta; G. Rocha

doi:10.1051/0004-6361/200912606

Home

All issues

Volume 522 (November 2010)

A&A, 522 (2010) A94

Full HTML

Free Access

Issue		A&A Volume 522, November 2010


Article Number		A94
Number of page(s)		29
Section		Astronomical instrumentation
DOI		https://doi.org/10.1051/0004-6361/200912606
Published online		05 November 2010

A&A 522, A94 (2010)

Residual noise covariance for Planck low-resolution data analysis

R. Keskitalo¹^,2^,3, M. A. J. Ashdown⁴^,5, P. Cabella⁶^,7, T. Kisner⁸, T. Poutanen¹^,2^,9, R. Stompor¹⁰, J. G. Bartlett¹⁰, J. Borrill⁸^,11, C. Cantalupo⁸, G. de Gasperis⁷, A. de Rosa¹², G. de Troia⁷, H. K. Eriksen¹³^,14, F. Finelli¹²^,15^,16, K. M. Górski³^,17, A. Gruppuso¹²^,15, E. Hivon¹⁸, A. Jaffe¹⁹, E. Keihänen¹, H. Kurki-Suonio¹^,2, C. R. Lawrence³, P. Natoli⁷^,20, F. Paci²¹, G. Polenta⁶^,22^,23 and G. Rocha³

¹ University of Helsinki, Department of Physics, PO Box 64, 00014 Helsinki, Finland
e-mail: reijo.t.keskitalo@jpl.nasa.gov; reijo.keskitalo@gmail.com
² Helsinki Institute of Physics, PO Box 64, 00014 Helsinki, Finland
³ Jet Propulsion Laboratory, California Institute of Technology, 4800 Oak Grove Drive, Pasadena CA 91109, USA
⁴ Astrophysics Group, Cavendish Laboratory, J J Thomson Avenue, Cambridge CB3 0HE, UK
⁵ Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK
⁶ Dipartimento di Fisica, Universitá di Roma “La Sapienza”, p.le A. Moro 2, 00185 Roma, Italy
⁷ Dipartimento di Fisica, Universitá di Roma “Tor Vergata”, via della Ricerca Scientifica 1, 00133 Roma, Italy
⁸ Computational Cosmology Center, Lawrence Berkeley National Laboratory, Berkeley CA 94720, USA
⁹ Metsähovi Radio Observatory, Helsinki University of Technology, Metsähovintie 114, 02540 Kylmälä, Finland
¹⁰ Laboratoire Astroparticule & Cosmologie (APC) – UMR 7164, CNRS, Université Paris Diderot, 10 rue A. Domon et L. Duquet, 75205 Paris Cedex 13, France
¹¹ Space Sciences Laboratory, University of California Berkeley, Berkeley CA 94720, USA
¹² INAF-IASF Bologna, Istituto di Astrofisica Spaziale e Fisica Cosmica di Bologna Istituto Nazionale di Astrofisica, via Gobetti 101, 40129 Bologna, Italy
¹³ Institute of Theoretical Astrophysics, University of Oslo, PO Box 1029 Blindern, 0315 Oslo, Norway
¹⁴ Centre of Mathematics for Applications, University of Oslo, PO Box 1053 Blindern, 0316 Oslo, Norway
¹⁵ INFN, Sezione di Bologna, Via Irnerio 46, 40126 Bologna, Italy
¹⁶ INAF-OAB, Osservatorio Astronomico di Bologna Istituto Nazionale di Astrofisica, via Ranzani 1, 40127 Bologna, Italy
¹⁷ Warsaw University Observatory, Aleje Ujazdowskie 4, 00478 Warszawa, Poland
¹⁸ Institut d’Astrophysique de Paris, 98bis boulevard Arago, 75014 Paris, France
¹⁹ Department of Physics, Blackett Laboratory, Imperial College London, South Kensington campus, London, SW7 2AZ, UK
²⁰ INFN, Sezione di “Tor Vergata”, via della Ricerca Scientifica 1, 00133 Roma, Italy
²¹ Instituto de Fisica de Cantabria, CSIC-Universidad de Cantabria, Avda. Los Castros s/n, 39005 Santander, Spain
²² ASI Science Data Center, c/o ESRIN, via G. Galilei snc, 00044 Frascati, Italy
²³ INAF-Osservatorio Astronomico di Roma, via di Frascati 33, 00040 Monte Porzio Catone, Italy

Received: 31 May 2009
Accepted: 15 June 2010

Abstract

Aims. We develop and validate tools for estimating residual noise covariance in Planck frequency maps, we also quantify signal error effects and compare different techniques to produce low-resolution maps.

Methods. We derived analytical estimates of covariance of the residual noise contained in low-resolution maps produced using a number of mapmaking approaches. We tested these analytical predictions using both Monte Carlo simulations and by applying them to angular power spectrum estimation. We used simulations to quantify the level of signal errors incurred in the different resolution downgrading schemes considered in this work.

Results. We find excellent agreement between the optimal residual noise covariance matrices and Monte Carlo noise maps. For destriping mapmakers, the extent of agreement is dictated by the knee frequency of the correlated noise component and the chosen baseline offset length. Signal striping is shown to be insignificant when properly dealt with. In map resolution downgrading, we find that a carefully selected window function is required to reduce aliasing to the subpercent level at multipoles, ℓ > 2N_side, where N_side is the HEALPix resolution parameter. We show that, for a polarization measurement, reliable characterization of the residual noise is required to draw reliable constraints on large-scale anisotropy.

Conclusions. Methods presented and tested in this paper allow for production of low-resolution maps with both controlled sky signal error level and a reliable estimate of covariance of the residual noise. We have also presented a method for smoothing the residual noise covariance matrices to describe the noise correlations in smoothed, bandwidth-limited maps.

Key words: cosmic microwave background / cosmology: observations / methods: data analysis / methods: numerical

© ESO, 2010

1. Introduction

Over the past two decades observations of the cosmic microwave background (CMB) have led the way towards the high-precision cosmology of today – a process best emphasized recently by the high-quality data set delivered by the Wilkinson microwave anisotropy probe (wmap) satellite (Hinshaw et al. 2007). The next major step in a continuing exploitation of the CMB observable will be data analysis of data sets from another satellite mission, Planck. Planck will observe the entire sky in multiple frequency channels, promising to improve on the recent wmap constraints on many fronts. In particular, as a satellite, Planck will provide us with unique access to the largest angular scales, in which the total intensity has proven controversial and difficult for theoretical interpretation and is still poorly measured and exploited in the polarization. Planck will be the only CMB satellite deployed in the current decade. It is therefore particularly important that the constraints on a large angular scale derived from the anticipated Planck data not only are robust, but they also efficiently exploit the information contained in them. This will be certainly necessary if Planck is to set strong constraints on the CMB B-mode power spectrum (Efstathiou et al. 2009) – one of the most attractive potential science targets of the mission.

In this paper we develop tools necessary for the statistically sound analysis of constraints on large angular scales. These include robust approaches to producing low-resolution maps and techniques for estimating pixel-pixel correlations from their residual noise. The work is a natural extension of the extensive mapmaking studies conducted within the Planck Working Group 3 (CTP) (Poutanen et al. 2006; Ashdown et al. 2007a,b, 2009). The low-resolution maps are expected to compress nearly all the information relevant to the large angular scale in fewer pixels making them more manageable. Given our Gaussian noise assumption, the full statistical description of the map uncertainty is given by its pixel-pixel noise covariance matrix (NCM). This is defined as $N = ⟨ (\hat{s} - s) {(\hat{s} - s)}^{T} ⟩, where ⟨ (\hat{s} - s) ⟩ = 0,$ $\begin{equation} \label{eq:ncvm} \vec N =\left\langle \left(\hat{\vec s}-\vec s\right) \left(\hat{\vec s}-\vec s\right)^\mathrm T \right\rangle, \quad\textrm{where}\quad \left\langle\left(\hat{\vec s}-\vec s\right)\right\rangle = 0, \end{equation}$ (1)where s is the 3N_pix input sky map of Stokes I, Q, and U parameters, and $\hat{\vec s}$ is our estimate of s. In the absence of signal errors, the difference $\left(\hat{\vec s}-\vec s\right)$ contains only instrument noise. We note that N is symmetric and usually dense matrix, which in general will depend on the mapmaking method that produced the estimate. In the following, we consider a number of numerical approaches to calculating such a matrix for each of the studied low-resolution maps and then test their consistency with the actual noise in the derived maps.

The full noise covariance matrices have been commonly computed and used in the analysis of the small-scale balloon-borne, (e.g. Hanany et al. 2000; de Bernardis et al. 2000) and ground-based experiments, (e.g., Kuo et al. 2004). The cosmic background explorer differential microwave radiometer (cobe-DMR) team also used them to derive low-ℓ constraints, (e.g., Górski et al. 1996; Wright et al. 1996). In all those cases, however, no resolution downgrading has been required, unlike with Planck, as the calculations for those experiments could be done at full resolution. To this date, only the wmap team has encountered a similar problem. The instrument noise model employed by them is, in fact, similar to the one used for Planck. It is parametrized, however, in the time domain rather than in the frequency domain (Jarosik et al. 2003, 2007). Calculation of the wmap NCM is formulated in exactly the same manner as for our optimal mapmaking method and the wmap likelihood code¹ ships with an NCM very similar to what we present here.

Our analysis is made unique by the differences in the experiment design: wmap pseudo-correlation receivers are differencing assemblies (DA) with two mirrors, whereas Planck will use a single mirror design (HFI, the high-frequency instrument) or has a reference load in place of the second mirror (LFI, the low-frequency instrument) (Planck Collaboration 2005). Between these, the pixel-pixel correlations are different. In principle the balanced load systems of cobe and wmap should bring less correlated noise than the unbalanced Planck LFI. On the other hand, differencing experiments generate pixel-pixel noise correlations even from white noise, whereas in Planck they originate in the correlated noise alone. In addition, we also study here the so-called destriping algorithms, which have been proposed as a Planck-specific mapmaking approach (Delabrouille 1998; Maino et al. 2002; Keihänen et al. 2004, 2005).

Residual noise covariance for Planck-like scanning and instrument noise has been studied before, either via some simplified toy-models (Stompor & White 2004) or in more realistic circumstances (Efstathiou 2005, 2007). Those studies approached the problem in a semi-analytic way and thus needed to employ some simplifications, which we avoid in our work. They were also based solely on the destriping approximation, assuming that noise can be accurately modeled by relatively long (one hour) baseline offsets and white noise, and did not consider any other approaches. In this work we extend those analyses into cases where modeling the noise correlations requires shorter baselines and compare those with optimal solutions.

2. Algebraic background

2.1. Maps and their covariances

To formulate mapmaking as a maximum likelihood problem we start with a model of the timeline: $d = A s + B x + n,$ $\begin{equation} \label{eq:tod} \vec d = \vec A\vec s + \vec B\vec x + \vec n, \end{equation}$ (2)where the underlying microwave sky signal, s, is to be estimated. Here A is the pointing matrix, which encodes how the sky is scanned and n is a Gaussian, zero mean noise vector. x denotes some extra instrumental effects, usually taken hereafter to be constant baseline offsets, which we will use to model the correlated part of the instrumental noise and B – a “pointing” matrix for x describing how it is added to the time domain data. Convolution with an instrumental beam, assumed here to be axially symmetric, is already included in s.

The signal part of the uncalibrated data vector, d, is the detector response to the sky emission observed in the direction of pixel p. For a total power detectors, e.g., the ones on Planck, it is a linear combination of polarized and unpolarized contributions: $d_{t} = 𝒦 {(1 + ϵ) s_{p I} + (1 - ϵ) (s_{p Q} \cos (2 χ_{t}) + s_{p U} \sin (2 χ_{t})^{)}} + n_{t},$ $\begin{equation} \label{eq:det_response} d_t = \mathcal K\left\{ (1+\epsilon)s_{p\mathrm I}\! +\! (1-\epsilon)\left( s_{p\mathrm Q}\cos (2\chi_t) + s_{p\mathrm U}\sin (2\chi_t) \right) \right\} + n_t, \end{equation}$ (3)where it is implied that sample t is measured in pixel p and χ_t is a detector polarization angle with respect to the polarization basis and we have dropped the baseline term for simplicity. Equation (3) includes an overall calibration factor 𝒦, and a cross polar leakage factor ϵ, however, in what follows, we only consider the case of perfect calibration, setting 𝒦 = 1 with no loss of generality, and no leakage, ϵ = 0.

To simplify future considerations we introduce a generalized pointing matrix, A′, and a generalized map, s′. They are defined as $A' \equiv [A, B], and s' \equiv [\begin{matrix} s \\ x \end{matrix}] allowingustowrite d = A' s' + n .$ $\begin{equation} \vec A' \equiv \left[ \vec A, \vec B\right], \mathrm{and} \, \vec s' \equiv \left[ \begin{array}{c} \vec s\\ \vec x \end{array} \right] \;\textrm{allowing us to write}\; \vec d = \vec A'\vec s' + \vec n. \label{eq:todprime} \end{equation}$ (4)The detector noise has a time-domain covariance matrix 𝒩 = ⟨ nn^T ⟩ and the probability for the observed timeline, d, becomes the Gaussian probability of a noise realization n = d − A′s′: $P (d) = P (n) = {[(2 π)^{\dim n} \begin{matrix} \det \end{matrix} 𝒩^{]}}^{- 1 / 2} \exp (- \frac{1}{2} n^{T} 𝒩^{-1} n) P (x),$ $\begin{equation} \label{eq:map_prob} P(\vec d) = P(\vec n) = \left[(2\pi)^{\mathrm{dim}\,\vec n}\det\mathcal N\right]^{-1/2} \exp \left( -\frac{1}{2}{\vec n}^\mathrm T\mathcal N^{-1}\vec n \right) P\left(\vec x\right), \end{equation}$ (5)where the last factor is a prior constraint on the noise offsets, x, which hereafter we will take to be a Gaussian with a zero mean and some correlation matrix, 𝒫, i.e., $P (x) \propto \exp (- \frac{1}{2} x T 𝒫^{-1} x) .$ $\begin{equation} P\left(\vec x\right) \propto \exp\left( -{1\over2}\vec x^\mathrm T\mathcal P^{-1}\vec x \right). \end{equation}$ (6)By maximizing this likelihood with respect to the sky signal and baselines contained in s′, we find an expression for a maximum likelihood estimate which reads, $\hat{s}' = {(ℛ^{-1} + {A'}^{T} 𝒩^{-1} A'^{)}}^{-1} {A'}^{T} 𝒩^{-1} d, where ℛ^{-1} \equiv [\begin{matrix} 0 & 0 \\ 0 & 𝒫^{-1} \end{matrix}] .$ $\begin{equation} \label{eq:min_var_map} \hat{\vec s}'\! =\! \left( \mathcal R^{-1} + {\vec A'}^\mathrm T\mathcal N^{-1}\vec A' \right)^{-1}\!{\vec A'}^\mathrm T\mathcal N^{-1}\vec d, \, \textrm{where}\, \mathcal R^{-1} \equiv \left[ \begin{array}{c c} \vec 0 & \vec 0\\ \vec 0 & \mathcal P^{-1} \end{array} \right]. \end{equation}$ (7)The first part of the vector $\hat{\vec s}'$ is an estimate of the actual sky signal, $\hat{\vec s}$ , while all the rest is an estimate of the baseline offsets, $\hat{\vec x}$ . The map $\vec A'^{\mathrm T}\mathcal N^{-1}\vec d$ is called the noise-weighted map. In case of flat prior for x, expressions identical to Eq. (7) can be also derived from minimum variance or generalized least square considerations and we will refer to $\hat{\vec s}$ as either a minimum variance or optimal map in the future. We note that we have ignored here any pixelization effects that cause differences between d and A′s′ even for a noise-free experiment. This is usually true in the limit of the pixel size significantly smaller than the beam resolution of the instrument. If this condition is not fulfilled, the pixelization effects may be important and special methods may be needed to minimize them. We discuss specific proposals in Sect. 2.3. In the absence of such effects, the difference between a map estimate, $\hat{\vec s}'$ , and the input map, s′, is called residual pixel domain noise.

Let us now consider first the prefactor matrix in Eq. (7), $M' \equiv {(ℛ^{-1} + {A'}^{T} 𝒩^{-1} A'^{)}}^{-1},$ $\begin{equation} \vec M'\equiv \left({\mathcal R}^{-1} + {\vec A'}^\mathrm T\mathcal N^{-1}\vec A'\right)^{-1}, \end{equation}$ (8)a weight matrix combining both the baseline prior and the noise variance weights. It acts on the generalized noise-weighted map, producing estimates of the pixels and baselines.

Given that our generalized map is made of two parts: the actual sky signal and the baseline offsets, the matrix M′ has four blocks: two diagonal blocks, denoted M and M_x, and two off-diagonal blocks, each of which is a transposed version of the other and one of which is referred here to as M_o. Using inversion by partition we can write an explicit expression for each of these blocks. For example, for the sky-sky diagonal blocks we get, $\begin{matrix} M = {[A T 𝒩^{-1} A - (A T 𝒩^{-1} B) {(𝒫^{-1} + B T 𝒩^{-1} B)}^{-1} {(A T 𝒩^{-1} B)}^{T}]}^{-1} \\ \times M x {(A T 𝒩^{-1} B)}^{T} {(A T 𝒩^{-1} A)}^{-1}, \end{matrix}$ $\begin{eqnarray} \label{eq:MmatSS0} && \vec M \!=\! \left[\!\vec A^\mathrm T\mathcal N^{-1}\vec A \!-\!\left(\vec A^\mathrm T\mathcal N^{-1}\vec B\right)\! \left(\mathcal P^{-1}\!+\!\vec B^\mathrm T\mathcal N^{-1}\vec B\right)^{-1}\! \left(\vec A^\mathrm T\mathcal N^{-1}\vec B\right)^\mathrm T \right]^{-1}~~~~~~~~ \\ && \label{eq:MmatSS}\quad = \left(\vec A^\mathrm T\mathcal N^{-1}\vec A\right)^{-1} +\left(\vec A^\mathrm T\mathcal N^{-1}\vec A\right)^{-1} \left(\vec A^\mathrm T\mathcal N^{-1}\vec B\right)\nonumber\\ &&\qquad \qquad\times \vec M_ x \left(\vec A^\mathrm T\mathcal N^{-1}\vec B\right)^\mathrm T \left(\vec A^\mathrm T\mathcal N^{-1}\vec A\right)^{-1}, \end{eqnarray}$ while for the offset-offset part, $\begin{matrix} M x & = & [𝒫^{-1} + B T 𝒩^{-1} B - (B T 𝒩^{-1} A) \end{matrix}$ $\begin{eqnarray} \vec M_x &= & \bigg[ \mathcal P^{-1} + \vec B^\mathrm T\mathcal N^{-1}\vec B - \left(\vec B^\mathrm T\mathcal N^{-1}\vec A\right) \nonumber \\&& \times \left(\vec A^\mathrm T\mathcal N^{-1}\vec A\right)^{-1} \left(\vec B^\mathrm T\mathcal N^{-1}\vec A\right)^\mathrm T \bigg]^{-1}. \label{eq:MmatXX} \end{eqnarray}$ (11)With help of these equations we can now write explicit separate expressions for the estimated sky signal and offsets. The former is given by, $\hat{s} = (M A T + M o B^{T}) 𝒩^{-1} d$ $\begin{equation} \hat{\vec s} = \left( \vec M\vec A^\mathrm T + \vec M_\mathrm oB^\mathrm T \right)\mathcal N^{-1}\vec d \label{eq:destripeMap1} \end{equation}$ (12)while the latter, $\hat{x} = (M \begin{matrix} T \\ o \end{matrix} A T + M x B^{T}) 𝒩^{-1} d .$ $\begin{equation} \hat{\vec x} = \left( \vec M_\mathrm o^\mathrm T \vec A^\mathrm T + \vec M_ x B^\mathrm T \right)\mathcal N^{-1}\vec d. \end{equation}$ (13)We can also combine these two equations to derive an alternative expression for the sky signal estimate, which makes a direct use of the offsets assumed to be estimated earlier, $\hat{s} = {(A T 𝒩^{-1} A)}^{-1} (A T 𝒩^{-1} d - A T 𝒩^{-1} B \hat{x}) .$ $\begin{equation} \hat{\vec s} = \left(\vec A^\mathrm T\mathcal N^{-1}\vec A\right)^{-1} \left(\vec A^\mathrm T\mathcal N^{-1}\vec d - \vec A^\mathrm T\mathcal N^{-1}\vec B\hat{\vec x} \right). \label{eq:destripeMap2} \end{equation}$ (14)If the assumed data model, Eq. (2), and the time domain noise and baseline covariances are all correct, then the covariance of the residual pixel domain noise is $N' = ⟨ (\hat{s}' - s') (\hat{s}' - s')^{T} ⟩ = M' .$ $\begin{equation} \label{eq:residual_covariance} \vec N' = \left\langle (\hat{\vec s}'-\vec s')(\hat{\vec s}'-\vec s')^\mathrm T\right \rangle = \vec M'. \end{equation}$ (15)In particular, Eq. (1), the pixel-pixel residual noise covariance matrix, N, is equal to M and given by Eq. (10).

Fig. 1

Noise model (dashed lines) and the spectra of simulated noise (solid lines). The two sets of curves correspond to the two considered knee frequencies with f_min = 1.15 × 10^-5 and α = 1.7.

We note that a sufficiently high-quality estimate of the inverse time domain correlations, 𝒩^-1, is required in order to calculate both the minimum-variance map and its noise covariance. If it is misestimated the map estimate will still be unbiased, though not any more minimum variance or maximum likelihood, and its covariance will not be given any more by Eq. (15).

For example for computational reasons we will find later that using some other matrix, denoted here as ℳ^-1, rather than the actual inverse noise covariance, 𝒩^-1, in the calculation of the map estimates in Eq. (7) can be helpful. The corresponding noise correlation matrix for such a map is then given by (Stompor et al. 2002, ungeneralized case) $\begin{matrix} N' = ⟨ (\hat{s}' - s') (\hat{s}' - s')^{T} ⟩ = {(ℛ^{-1} + {A'}^{T} ℳ^{-1} A'^{)}}^{-1} \\ \times (ℛ^{-1} ℛ_{𝒩} ℛ^{-1} + {A'}^{T} ℳ^{-1} 𝒩 ℳ^{-1} A'^{)} {(ℛ^{-1} + {A'}^{T} ℳ^{-1} A'^{)}}^{-1}, \end{matrix}$ $\begin{eqnarray} \label{eq:generalized_residual_covariance} &&\vec N' = \left\langle (\hat{\vec s}'-\vec s')(\hat{\vec s}'-\vec s')^\mathrm T \right\rangle = \left( \mathcal R^{-1} + {\vec A'}^\mathrm T\mathcal M^{-1}\vec A' \right)^{-1}\nonumber\\ &&\;\;\times\left( \mathcal R^{-1}\mathcal R_\mathcal N\mathcal R^{-1} + {\vec A'}^\mathrm T\mathcal M^{-1}\mathcal N\mathcal M^{-1} \vec A' \right) \left( \mathcal R^{-1} + {\vec A'}^\mathrm T\mathcal M^{-1}\vec A' \right)^{-1}, \end{eqnarray}$ (16)where ℳ and ℛ define our mapmaking operator, whereas 𝒩 and ℛ_𝒩 are the true noise properties. This expression is significantly more complex and computationally involved than Eq. (15). Fortunately, as we discuss in the following, in many cases of interest, the latter expression turns out to be a sufficiently good approximation of the former with 𝒩^-1 replaced by ℳ^-1 at least for some of the potential applications.

The mapmaking methods considered for Planck fall into two general classes both of which are described by the equations derived above. The first class, called optimal methods, assumes the sufficient knowledge of the time domain noise correlations, does not introduce any extra degrees of freedom, x, (or alternately assumes a prior for them with 𝒫 vanishing). In this case one can derive from Eqs. (14) and (10), $\begin{matrix} \hat{s} & = & {(A^{T} 𝒩^{-1} A)}^{-1} A^{T} 𝒩^{-1} d, \\ N & = & {(A^{T} 𝒩^{-1} A)}^{-1}, \end{matrix}$ $\begin{eqnarray} \label{eq:optimal_map} \hat{\vec s} & = & \left( {\vec A}^\mathrm T\mathcal N^{-1}\vec A \right)^{-1}\, {\vec A}^\mathrm T\mathcal N^{-1}\vec d, \\ \label{eq:optimal_ncm} \vec N & = & \left( {\vec A}^\mathrm T\mathcal N^{-1}\vec A \right)^{-1}, \end{eqnarray}$ where the latter is a covariance of the residual noise of the former.

The second class of methods introduces a number of baseline offsets with a Gaussian correlated (so-called generalized destripers) or uncorrelated (standard destripers) prior on them and describe the noise as an uncorrelated Gaussian process. The destriper maps are evaluated via Eq. (12) or Eq. (14), with 𝒩 assumed to be diagonal. Clearly on the time-domain level the destriper model is just an approximation, therefore at least a priori we should use the full expression in Eq. (16) to estimate its covariance. The CTP papers (Poutanen et al. 2006; Ashdown et al. 2007a,b, 2009) have shown that for Planck the two methods produce maps which are very close to one another. Moreover, they have shown that using a generalized destriper, the derived maps eventually become nearly identical to those obtained with the optimal methods, if an appropriate length of the baseline and a number of the baseline offsets is adopted together with a consistently evaluated prior. This motivates using the simplified Eq. (10) as an approximation for the noise covariance of the destriped maps, i.e., N = M. We will investigate the quality and applicability of this approximation later in this paper.

2.2. Time domain noise

We assume that the time domain noise is a Gaussian process and for the simulations we take the noise power spectral density to have the form $P (f) = \frac{σ^{2}}{f_{sample}} \cdot \frac{f_{knee}^{α} + f^{α}}{f_{\min}^{α} + f^{α}},$ $\begin{equation} \label{eq:psd} P(f) = \frac{\sigma^2}{f_\mathrm{sample}}\cdot \frac{f_\mathrm{knee}^\alpha+f^\alpha}{f_\mathrm{min}^\alpha+f^\alpha}, \end{equation}$ (19)where the shape is defined by the slope, minimum and knee frequencies (α, f_min and f_knee respectively) and the scaling by the white-noise sample variance and sampling frequency (σ and f_sample). Two examples of the theoretical and simulated noise spectra can be seen in Fig. 1.

In the calculations of the maps using the optimal algorithms or generalized destripers we will assume that noise power spectrum is known precisely. As the noise simulated in the cases analyzed here is piece-wise stationary, with no correlations allowed between the data in the different pieces (see Sect. 5) the respective noise correlation matrix, 𝒩, is block Toeplitz with each of the blocks, describing the noise correlations of one of the stationary pieces, defined by the noise power spectrum. Given that we will approximate the inverse of 𝒩 as also a block Toeplitz matrix with each blocks given by an inverse noise power spectrum. Though this is just an approximation it has been demonstrated in the past that it performs exquisitely well in particularly in the cases with long continuous pieces of the stationary noise (Stompor et al. 2002), as it is a case in all simulations considered here.

2.3. Low-resolution maps

The high-resolution (down to ~5 arcmin) of the Planck data makes it impossible with current resources to evaluate and store the full pixel-pixel covariance of all map pixels. The solution is to produce a low-resolution dataset for the study of large scale structure of the microwave sky. Low-resolution mapmaking for a high-resolution mission is non-trivial for two reasons: first, the considered mapmaking methods assume that there exists no sky structure at the sub pixel scale. Secondly, cosmological information contained in the polarization signal has several orders of magnitude worse signal-to-noise ratio than intensity signal, causing completely different downgrading effects in the Stokes I, Q and U maps.

In this section we define three alternative methods of producing low-resolution maps from high-resolution observations. The first two, direct and (inverse) noise weighting methods, have already been used in the WMAP analysis (Jarosik et al. 2007). As the third option we consider smoothed (low-pass-filtered) maps and their noise covariance.

2.3.1. Direct method

The most straightforward method to produce a low-resolution map is to project the detector observations directly to the pixels of the final target resolution. Hereafter we will refer to it as the direct method.

Solving the low-resolution map directly involves less uncertainty. This means that the solved and subtracted correlated noise realization is more accurate. The absence of intermediate, high-resolution products makes noise in the direct maps easier to model. The noise covariance matrices for the direct maps can be computed using formalism presented in Sect. 2.1, for example, Eqs. (15) and (16).

On the other hand, if the sky signal varies within the pixels, non-uniform sampling of a pixel can cause a systematic error in the solved map (Ashdown et al. 2007a). This is an important aspect, as none of the considered algorithms deal with sub pixel structures by design.

In Sect. 6.3 we quantify the signal error for the different low-resolution maps and mapmakers.

Hereafter, we will use the direct method as a reference with respect to which we compare the other approaches.

2.3.2. Inverse noise weighting (INW)

In the case of nested pixelization schemes, such as Hierarchical Equal Area isoLatitude Pixelization (HEALPix, Górski et al. 2005) used in this paper, to downgrade a temperature only map, one may compute a weighted average of the sub pixel temperatures. A natural choice are the optimal, inverse noise weights that lead to minimum noise of the super pixel in the absence of correlated noise. We will refer to these maps as inverse noise weighted (INW) maps.

A similar weighting scheme exists for the polarized data as well. Algebraically, the entire procedure can be summarized on the map level as, $\hat{s}' = W ′′ X W \hat{s},$ $\begin{equation} \label{eq:noise_weighting} \hat{\vec s}' = \vec W''\,\vec X \,\vec W \,\hat{\vec s}, \end{equation}$ (20)where W and W′′ are weight matrices for the high and low resolution maps, respectively. They depend on A and A′′ – the pointing matrices at high and low resolution. X simply sums the pixels in resolution N_pix to resolution $N_\mathrm{pix}'$ , $X qp = {\begin{matrix} 1, & p subpixel of q, \\ 0, & otherwise \end{matrix}$ $\begin{equation} \label{eq:X} \vec X_{qp} = \left\{ \begin{array}{cc} 1,&p\;\mathrm{subpixel}\;\mathrm{of}\;q,\\ 0,&\mathrm{otherwise}\\ \end{array}\right. \end{equation}$ (21)weights are given by, $W ′′ = {({A ′′}^{T} 𝒩_{u}^{-1} A {′′}^{)}}^{-1} and W = {(A T 𝒩_{u}^{-1} A)}^{-1},$ $\begin{equation} \vec W'' = \left({\vec A''}^\mathrm T\mathcal N_\mathrm u^{-1}\vec A''\right)^{-1} \quad \textrm{and}\quad \vec W = \left(\vec A^\mathrm T\mathcal N_\mathrm u^{-1}\vec A\right)^{-1}, \end{equation}$ (22)and 𝒩_u is the time domain covariance matrix of the uncorrelated part of the noise, n.

The covariances for the maps obtained via such a procedure can be derived from Eq. (20) and the expressions described in Sect. 2.1.

Noise weighting reduces signal errors by first solving the map at a resolution where sub pixel structure is weak. In comparison to the direct method, the more accurate signal is gained at the cost of higher noise. Like the direct method, INW disregards any sub pixel structure but at high-resolution. In this case the instrument beam naturally smoothes out small scale structure causing the approximation to hold.

2.3.3. Harmonic smoothing

Both methods described above may result in the signal smoothing (or its band-width) being position-dependent as it is achieved via averaging of the observed high-resolution pixel amplitudes contained within each low-resolution pixel.

Applying a smoothing operator to the high-resolution map prior to resolution downgrading could alleviate such a problem. The smoothing operation needs to properly take care of the high frequency power contained in the maps, not to alias it as power at the scales of interest.

As the smoothing operation is usually performed in the harmonic domain, it requires that the high-resolution map is first expanded in spherical harmonics. If the sky is only partially covered, the expansion is not unique. Since Planck is a full sky experiment, we do not address the partial sky problem here but leave it for future work.

To suppress small angular scale power the expansion is convolved with an axially symmetric window function (e.g. a symmetric Gaussian window function Challinor et al. 2000), $\begin{matrix} {a_{ℓm}^{T}}_{˜} = W_{ℓ} a_{ℓm}^{T}, W_{ℓ} = e^{- \frac{1}{2} ℓ (ℓ + 1) σ^{2}} \\ {a_{ℓm}^{E}}_{˜} =_{2} W_{ℓ} a_{ℓm}^{E},_{2} W_{ℓ} = e^{- \frac{1}{2} [ℓ (ℓ + 1) - 4] σ^{2}}, \end{matrix}$ $\begin{eqnarray} \label{eq:beam_window} & & \tilde a^\mathrm T_{\ell m} = W_\ell a^\mathrm T_{\ell m}, W_\ell = {\rm e}^{-\frac{1}{2}\ell(\ell+1)\sigma^2} \\ & &\tilde a^\mathrm E_{\ell m} = {}_2W_\ell a^\mathrm E_{\ell m}, _2W_\ell = {\rm e}^{-\frac{1}{2}[\ell(\ell+1)-4]\sigma^2}, \end{eqnarray}$ chosen to leave only negligible power at angular scales that are not supported by the low target resolution. Finally the regularized expansion is synthesized into a low-resolution map by sampling the expansion values at pixel centers. We conduct most of our studies using a beam having a full width at half maximum (FWHM) of twice the average pixel side. For the N_side = 32 resolution² this is approximately 220′ ( $3\fdg 7$ ). Whenever transforms between harmonic and pixel space are conducted, it is important to consider the range of multipoles included in the transformation. We advocate using such an aggressive smoothing that the harmonic expansion has negligible power beyond ℓ = 3N_side and results are stable for any ℓ_max beyond this. For completeness we have set ℓ_max = 4N_side but stress that any residual power beyond ℓ = 3N_side will lead to aliasing.

The smoothing window does not need to be a Gaussian but it is preferable to avoid sharp cut-offs that may induce “ringing” phenomena. Benabed et al. (2009) suggest a window function that preserves the signal basically unchanged until a chosen threshold and then smoothly kills all power quickly above that angular resolution. Their window is $W_{ℓ} = {\begin{matrix} 1, & ℓ \leq ℓ_{1} \\ \frac{1}{2} [1 + \cos ((ℓ - ℓ_{1}) π / (ℓ_{2} - ℓ_{1}))], & ℓ_{1} < ℓ \leq ℓ_{2} \\ 0, & ℓ > ℓ_{2} \end{matrix}$ $\begin{equation} \label{eq:teasingbeam} W_\ell = \left\{ \begin{array}{cl} 1, & \ell \leq \ell_1 \\ \frac{1}{2}\left[ 1+\cos\left((\ell-\ell_1)\pi/(\ell_2-\ell_1)\right)\right], & \ell_1 < \ell \leq \ell_2 \\ 0, & \ell > \ell_2 \end{array} \right. \end{equation}$ (25)with the typical choice ℓ₁ = 5N_side / 2 and ℓ₂ = 3N_side.

In contrast to the two noise-driven methods, harmonic smoothing can be considered optimal from the (large-scale) signal viewpoint; however it may be suboptimal as far as the noise is concerned, in particular in the case of strongly inhomogeneous noise distribution.

The noise covariance matrices described in Sect. 2.1 need to be amended to accurately characterize the smoothed maps and thus we need to smooth the matrices as well. Smoothing of a map is a linear operation. For any linear operator, L, acting on a map, m, we can compute its covariance as $⟨ (L m) (L m)^{T} ⟩ = L ⟨ m m T ⟩ L T = L N L T = \sum_{i} λ_{i} L \hat{e} i \cdot (L \hat{e} i)^{T},$ $\begin{equation} \langle(\vec L\vec m)(\vec L\vec m)^\mathrm T\rangle = \vec L\langle\vec m \vec m^\mathrm T\rangle\vec L^\mathrm T = \vec L\vec N\vec L^\mathrm T = \sum_i \lambda_i\,\vec L\hat{\vec e}_i\cdot(\vec L\hat{\vec e}_i)^\mathrm T, \label{eq:covSmoothing} \end{equation}$ (26)where λ_i and e_i stand for eigenvalues and eigenvectors of the noise covariance, N, i.e., $\vec N =\sum_i\lambda_i\widehat{\vec e}_i\widehat{\vec e}_i^\mathrm T$ , and m is understood here to contain only noise. We note that in general one should replicate the same processing steps for the covariance matrix as are to be applied to maps and therefore the smoothing operation should be applied to the noise covariance of the high-resolution map and its resolution downgraded later.

In order to make the evaluation of Eq. (26) feasible, we must commute the order of the smoothing and downgrading operations contained in L. Though these two operations are clearly not exchangeable on the map level, due to the potential presence of the subpixel power, such an approach can be more justifiable for the noise covariances. In this case we can explicitly compute the low-resolution unsmoothed covariance matrices directly and subsequently smooth them with the signal smoothing kernel.

In the following we will apply the smoothing technique to both high and low resolution maps already downgraded using some of the other approaches. We will demonstrate that such a combined approach results in controllable properties of the residual noise on the one hand and well defined sky signal bandwidth on the other. Unlike both the direct and INW methods, the smoothed low-resolution maps are actually solved from a signal that lacks sub pixel structure.

3. Numerical calculations of residual noise covariance

This section presents numerical methods to compute the residual noise covariance matrix and describes briefly their implementations corresponding to three different mapmaking methods, the optimal method (MADping and ROMA implementations) and the generalized (Madam) and classical destriping (Springtide) methods.

3.1. Optimal map covariance

The noise covariance of optimal maps assuming exact noise model is given by Eq. (18). It is evaluated in two steps: first the inverse covariance matrix, A^T𝒩^-1A, is assembled and subsequently inverted. Given that the matrix can be singular the latter step may require computing the pseudo-inverse.

The pseudo-inverse is defined by the eigenvalue decomposition of the inverse noise matrix. Because the noise matrix is symmetric and, by definition, non-negative definite, its eigenvalues are real and non-negative (λ_i ≥ 0), and its eigenvectors form a complete orthogonal basis. This allows us to expand the matrix as $N -1 = U Λ U T = \begin{matrix} \sum \end{matrix} λ_{i} \hat{e} i \hat{e} \begin{matrix} T \\ i \end{matrix} .$ $\begin{equation} \label{eq:eigevaluedecomposition} \vec N^{-1} = \vec U \vec \Lambda \vec U^\mathrm T = \sum \lambda_i\,\hat{\vec e}_i \hat{\vec e}_i^\mathrm T. \end{equation}$ (27)Here λ_i are the eigenvalues and $\hat{\vec e}_i$ are the corresponding eigenvectors of the matrix N^-1. We can now invert N^-1 by using its eigenvalue decomposition, $N = \sum_{i} λ_{i}^{-1} \hat{e} i \hat{e} \begin{matrix} T \\ i \end{matrix},$ $\begin{equation} \label{eq:eigendecomp} \vec N = \sum_i \lambda_i^{-1}\hat{\vec e}_i\hat{\vec e}_i^\mathrm T, \end{equation}$ (28)and controlling the ill-conditioned eigenmodes. Any ill-conditioned eigenmode will have an eigenvalue several orders of magnitude smaller than the largest eigenvalue. By including in the sum only the well-conditioned eigenmodes we effectively project out the correlation patterns that our methods cannot discern.

We describe in Sect. A.3, that the covariance is indeed commonly expected to be singular or nearly so.

MADping is one of the codes of the MADCAP³ suite of CMB analysis tools (Borrill 1999). In this implementation, each processor scans through its sections of time-ordered data and accumulates its local piece of the inverse noise covariance. These pieces are then gathered and written to disk. The scaling of this technique is $N_{flops} ~ 𝒪 ({n_{samples} \cdot n_{correlation}}^{)},$ $\begin{equation} \label{eq:madping_scaling} N_\mathrm{flops} \sim \mathcal{O}\left( {n_\mathrm{samples}\,\cdot\,n_\mathrm{correlation}} \right), \end{equation}$ (29)where n_correlation is the filter length set by the noise autocorrelation length.

For a Planck-sized, full-sky dataset and using a reasonable pixel resolution (half a degree), the construction of the inverse noise covariance dominates over the computational cost of inverting this matrix. Nevertheless, inversion methods as well as the eigendecomposition scale as $\mathcal O(N_\mathrm{pix}^3)$ .

In order to correctly treat the signal component of the data map in Eq. (7), we must apply our low-resolution noise covariance to a noise-weighted map which has been downgraded from higher resolution. This downgrading process is equivalent to the technique discussed in Sect. 2.3.2, and ensures that signal variations inside a low-resolution pixel are accounted for. The high-resolution noise-weighted map consistent with the above formalism is constructed as the first step of the mapmaking carried out by the MADmap program (Cantalupo et al. 2010). The matrix eigenvalue decomposition is done using a ScaLAPACK⁴ interface that allows efficient parallel eigenvalue decomposition of large matrices using a divide-and-conquer algorithm.

The ROMA code (Natoli et al. 2001; de Gasperis et al. 2005) is an implementation of the optimal generalized least squares (GLS) iterative map making specifically designed for Planck, but also successfully used on suborbital experiments such as BOOMERanG (Masi et al. 2006). To estimate noise covariance as in Eq. (15) we start by calculating row i of its inverse, N^-1, by computing its action on the unit vector along axis i: $(N -1)_{i :} = (A T 𝒩^{-1} A) \cdot e i .$ $\begin{equation} (\vec N^{-1})_{i:} = (\vec A^\mathrm{T} \mathcal N^{-1}\vec A) \cdot \vec e_i. \end{equation}$ (30)This calculation is implemented as follows: i) projecting the unit vector into a TOD by applying A on $\hat{\vec e}_i$ ; ii) noise-filtering the TOD in Fourier space; iii) projecting the TOD into a map by applying A^T. By computing each column independently we reduce memory usage because we allocate memory only for one map instead of allocating memory for the full matrix. The computational cost of the full calculation is dominated by Fourier transforms that are repeated as many times as the number of columns, hence the scaling can be expressed as: $N_{flops} ~ 𝒪 (n_{samples} \cdot \log_{2} (n_{filter}) \cdot n_{pix}) .$ $\begin{equation} N_\mathrm{flops} \sim \mathcal O\left( n_\mathrm{samples} \cdot\log_2 (n_\mathrm{filter})\cdot n_\mathrm{pix} \right). \end{equation}$ (31)Alternately, we can compute the NCM column-by-column with help of multiple mapmaking-like operations, i.e., $y i : = (A T 𝒩^{-1} A)^{-1} \cdot e i,$ $\begin{equation} \vec{y}_{i:} = (\vec A^\mathrm{T} \mathcal N^{-1}\vec A)^{-1} \cdot \vec e_i, \end{equation}$ (32)where $\hat {\vec e}_i$ is a unit vector as defined above and y_i stands for a column of N. We rewrite the above equation as $(A T 𝒩^{-1} A) \cdot y_{i} = e i,$ $\begin{equation} (\vec A^\mathrm{T} \mathcal N^{-1}\vec A) \cdot {\vec y}_i = \vec e_i, \end{equation}$ (33)and solve it using the standard PCG mapmaking solver. In such case there is no need to store a full inverse noise covariance matrix in the memory of a computer at any single time as the operations on the lhs can be performed from right to left. As a result, this approach can be applied also for high-resolution cases for which the direct method described above would not be feasible.

We note that unlike in the previous approaches based on the direct matrix inversion in the latter case there is no special care taken of potential singularities. Though the presence of those does not hamper the PCG procedure (e.g., Cantalupo et al. 2010), nonetheless care must be taken while interpreting its results.

3.2. Destriped map covariance

In the destriping approach to mapmaking, we model all noise correlations by baseline offsets. Thus we write Eq. (2) as $d = A s + B x + n u,$ $\begin{equation} \label{eq:ds_tod} \vec d = \vec A\vec s + \vec B\vec x + \vec n_\mathrm u, \end{equation}$ (34)where n_u is a vector of uncorrelated white noise samples. Accordingly, we must replace the time-domain noise covariance matrix, 𝒩, by a diagonal matrix, 𝒩_u. All noise correlation is then included in the prior baseline offset covariance matrix, 𝒫.

If we now apply the destriping approximation to Eqs. (9) we find for the pixel-pixel residual noise covariance matrix: $M -1 = A T 𝒩_{u}^{-1} A - A T 𝒩_{u}^{-1} B {(𝒫^{-1} + B T 𝒩_{u}^{-1} B)}^{-1} B T 𝒩_{u}^{-1} A .$ $\begin{equation} \label{eq:madam_corr} \vec M^{-1} = \vec A^\mathrm T\mathcal N_\mathrm u^{-1}\vec A - \vec A^\mathrm T\mathcal N_\mathrm u^{-1} \vec B \left( \mathcal P^{-1}+\vec B^\mathrm T\mathcal N_\mathrm u^{-1}\vec B \right)^{-1} \vec B^\mathrm T\mathcal N_\mathrm u^{-1}\vec A. \end{equation}$ (35)The first term on the rhs is the binned white noise contribution (a diagonal or block-diagonal matrix for temperature-only and polarized cases respectively) and the second term describes the pixel-pixel correlations due to errors in solving for the baselines, i.e., the difference between the solved and actual baselines (Kurki-Suonio et al. 2009).

When making a map using destriping, one can use high resolution to solve for baselines and still bin the map at low resolution. Since this is equivalent to producing first a high-resolution map and then downgrading through inverse noise weighting, we will always assume the same pixel size for both of these steps.

3.2.1. Conventional destriping

Springtide (Ashdown et al. 2007b) is an implementation of the conventional destriping approach which solves for one baseline per pointing period. Since the baselines are so long, it allows for a number of optimizations in the handling of the data. During one pointing period, the same narrow strip of sky is observed many times, so the time-ordered data are compressed into rings before doing the destriping. Another effect of the long baselines is that the prior covariance matrix of the baselines, 𝒫, is strongly diagonal-dominant, so to a very good approximation can be assumed to be diagonal. As a consequence, the matrix $\left(\mathcal P^{-1}+\vec B^\mathrm T\mathcal N_\mathrm u^{-1}\vec B\right)^{-1}$ that appears in the expression for the inverse map covariance matrix (35) is also diagonal. Thus the number of operations taken to compute (35) is $N_{flops} ~ 𝒪 (n_{base} (n_{pix / base})^{2}) .$ $\begin{equation} N_\mathrm{flops} \sim \mathcal{O}\left(n_{\mathrm{base}}(n_{\mathrm{pix/base}})^2\right). \end{equation}$ (36)The number of baselines is small compared to the generalized destriping approach, so it is also possible to compute the inverse posterior covariance matrix of the baseline offsets explicitly, invert it and use it to compute the map covariance matrix. This method has the advantage that the resolution at which the destriping is performed is not constrained to be the same as the resolution of the final map covariance matrix. The destriping can instead be done at the natural resolution of the data to avoid subpixel striping effects. The inverse of the posterior baseline error covariance matrix can be calculated using Eq. (11) and the corresponding map covariance matrix is given by Eq. (10). However, the pointing matrices, A, need not be for the same resolution in both steps.

The structure of the inverse posterior baseline covariance matrix, Eq. (11), depends on the scanning strategy, but it is in general a dense matrix. Inverting the matrix and using it to compute Eq. (10) involves dense matrix operations, so this method is computationally more demanding than the other approach described above. However, the posterior baseline covariance matrix needs only to be computed once and stored, and then can be used many times to compute the map covariance matrix for any desired resolution. It is also possible to use Eq. (10) to compute the residual noise covariance matrix for a subset of the pixels in the map.

3.2.2. Generalized destriping

Madam (Keihänen et al. 2005, 2010) is an implementation of the generalized destriping principle. It is flexible in the choice of baseline length and makes use of prior information of baseline covariance (𝒫 is not approximated to be diagonal).

Even for a generalized destriper, all the matrices in Eq. (35) are extremely sparse. Most of the multiplications only require operations proportional to the number of pixels, baselines or samples. The matrices 𝒫^-1 and $\left(\mathcal P^{-1}+\vec B^\mathrm T\mathcal N_\mathrm u^{-1}\vec B\right)^{-1}$ are approximately circulant, band-diagonal matrices whose width is determined by the noise spectrum. For all cases studied in this paper, the latter matrix is limited to the order of 10³ non-negligible elements per row. We call this width the baseline correlation length, n_corr. n_corr corresponds to a few hours of samples and is inversely proportional to baseline length, n_bl. In addition to the correlated noise part, our formalism naturally includes pixel correlations due to white noise.

We evaluate the prior baseline offset matrix from the power spectral density (PSD) of the correlated noise, P(f), by Fourier transforming the baseline PSD, P_x(f), into the autocorrelation function. The baseline PSD is evaluated as (Keihänen et al. 2010): $\begin{matrix} P_{x} (f) & = & \frac{1}{t_{base}} \sum_{m = - \infty}^{\infty} P (f + m / t_{base}) g (f t_{base} + m), \\ where & g (x) = \frac{\sin^{2} πx}{(πx)^{2}} \cdot \end{matrix}$ $\begin{eqnarray} \label{eq:psd_x} P_ x(f) &= &\frac{1}{t_\mathrm{base}}\sum_{m=-\infty}^{\infty} P(f+m/t_\mathrm{base})g(ft_\mathrm{base}+m), \nonumber\\ \mathrm{where}&& \quad g(x) = \frac{\sin^2 \pi x}{(\pi x)^2}\cdot \end{eqnarray}$ (37)The sum converges after including only a few m around the origin. For stationary noise, any row of 𝒫^-1 can then be evaluated as a cyclic permutation of ℱ^-1 [ 1 / P_x(f) ] .

We evaluate (35) after computing (37) by approximating the inner matrix, $\left(\mathcal P^{-1}+\vec B^\mathrm T\mathcal N_\mathrm u^{-1}\vec B\right)^{-1}$ , as circulant. This allows us to invert the n_base × n_base band diagonal matrix by two short Fourier transforms. It turns out that the matrix multiplications are most conveniently performed by first evaluating the sparse $\vec A^\mathrm T\vec N_\mathrm u^{-1}\vec B$ matrix and then operating with it on the inner matrix from both sides.

Actually the inverse covariance matrix gains a contribution from all quadruplets (x_i,x_j,p,q), where baselines x_i and x_j are within baseline correlation length and hit pixels p and q. The number of operations required to complete the estimate is then proportional to $N_{flops} ~ 𝒪 (n_{base} \cdot n_{corr} \cdot (n_{pix / base})^{2}),$ $\begin{equation} \label{eq:madam_scaling} N_\mathrm{flops} \sim \mathcal{O}\left( n_\mathrm{base} \cdot n_\mathrm{corr} \cdot (n_\mathrm{pix/base})^2 \right), \end{equation}$ (38)where n_base and n_corr, the number of baselines per survey and correlation length respectively, are inversely proportional to the length of a baseline. In contrast, n_pix / base, pixels per baseline, is proportional to baseline length. For short baselines and low-resolution maps, the magnitude of n_pix / base is close to unity. Table 1 lists low-resolution N_side parameters and the baseline lengths that correspond to average pixel sizes. It can be used to estimate n_pix / base and shows that, for example, 1.25s baseline offset at N_side = 32 resolution covers approximately 4 pixels. The success of this approach is to replace the number of baselines squared, $n_\mathrm{bl}^2$ , in the computation complexity by $n_\mathrm{pix/base}^2$ .

Table 1

Pixel side to baseline length at 1 rpm spin rate.

3.3. Smoothed covariance matrices

We define the linear smoothing operator, L, of Eq. (26) as $L = Y † W Y,$ $\begin{equation} \vec{L}=\vec{Y}^{\dagger}\vec{W}\vec{Y}, \end{equation}$ (39)where each row of the matrix Y is a map of the relevant spin-0 and spin-2 spherical harmonic functions. Y acts on a map vector and returns its spherical harmonic expansion: a = Ym. W is the smoothing kernel matrix that suppresses power at the chosen scales according to Eq. (23).

The smoothing procedure on its own will often lead to singular eigenvectors of the smoothed covariance matrix, with the eigenvalues corresponding to those close to zero. Though at a first glance such eigenvectors may look like being strongly constrained by the data, their actual value in the analysis is negligible as the sky signal in those modes is also smoothed. The nearly vanishing variance of those modes will often spuriously exaggerate the smoothing and mapmaking artifacts likely present whenever any inverse noise weighting needs to be applied. To avoid such problems hereafter we apply the pseudo-inverse described in 3.1 The criterion for selecting the nearly singular modes will in general depend on the case at hand and we will discuss the matter with the χ² test in Sect. 6.2.

In some cases the eigenvalue decomposition of the smoothed matrix may not be readily available or its computation not desirable. We can then regularize the inversion of $\vec {\tilde N}$ by adding some low level of the pixel-independent uncorrelated noise. For consistency, a random realization of such noise should also be added to the corresponding maps. We note that both approaches are effectively equivalent and that the choice of the singularity threshold needed to select the singular eigenmodes corresponds roughly to the choice of the noise level to be added. We will commonly use the latter approach in some of the power spectrum tests discussed later.

4. Numerical tests and comparison metric

This paper has two main goals. On the one hand we propose and compare various methods devised to produce the low-resolution maps, searching for the mapmaking method which leads to low-resolution maps virtually free of artifacts such as those due to subpixel power aliasing. In parallel we develop the tools to estimate residual noise covariance for such maps, which properly describe the error of the estimated maps due to residual noise. The numerical results presented in the subsequent sections of the paper aim therefore at comparing and validating the algorithms which we have presented earlier. The discussed comparisons involve standard statistical tests, such Kolmogorov-Smirnov, χ², etc. and ones which are specifically devised in the light of the anticipated future applications of the maps and their covariances. We describe such tests in this section.

4.1. Quadratic maximum likelihood power spectrum estimation

One of the main applications for the low-resolution noise covariance matrices we envisage is to the estimation of power spectra, C_ℓ. This is often separated into the estimation of large and small angular scales, usually associated with high and low signal to noise regimes. The methods discussed here are relevant for large angular scales. The successful estimation of the underlying true power spectrum of the sky signal sets demanding requirements for the quality of the maps produced for such a purpose as well as the consistency of the estimated noise covariance and the actual noise contained in the map. For this reason we will use hereafter power spectrum estimation as one of the metrics with which to evaluate the quality of the proposed algorithms.

We will use the quadratic maximum likelihood (QML) method for the power spectrum estimation as introduced in (Tegmark 1997) and later extended to polarization in (Tegmark & de Oliveira-Costa 2001). Given a temperature and polarization map, m, the QML provides estimates of the power spectra, that read, $Ĉ \begin{matrix} X \\ ℓ \end{matrix} = \sum_{ℓ^{'}, X^{'}} (F -1)_{ℓ ℓ^{'}}^{X X^{'}} [m T E \begin{matrix} ℓ^{'} \\ X^{'} \end{matrix} m - tr (N E_{X^{'}}^{ℓ^{'}})^{]} .$ $\begin{equation} \hat{C}_\ell^X = \sum_{\ell' , X'} (\vec F^{-1})^{X \, X'}_{\ell\ell'} \left[ \vec m^\mathrm T\vec E^{\ell'}_{X'}\vec m -\mathrm{tr}(\vec N{\vec E}^{\ell'}_{X'}) \right]. \label{eq:QMLestim} \end{equation}$ (40)Here, $\hat {C}_\ell^X$ is an estimated power spectrum, X = TT,EE,TE,BB,TB, or EB, and $\vec F_{X X'}^{\ell \ell '}$ is the Fisher matrix. F and E are defined as $F \begin{matrix} ℓ ℓ^{'} \\ X X^{'} \end{matrix} = \frac{1}{2} tr [C^{-1} \frac{\partial C}{\partial C_{ℓ}^{X}} C^{-1} \frac{\partial C}{\partial C_{ℓ^{'}}^{X^{'}}}] and E \begin{matrix} ℓ \\ X \end{matrix} = \frac{1}{2} C^{-1} \frac{\partial C}{\partial C_{ℓ}^{X}} C -1,$ $\begin{equation} \vec F^{\ell\ell'}_{X X'} =\frac{1}{2}\mathrm{tr}\Big[{\vec C}^{-1}\frac{\partial {\vec C}}{\partial C_\ell^X}{\vec C}^{-1}\frac{\partial {\vec C}}{\partial C_{\ell'}^{X'}}\Big] \; \textrm{and} \; \vec E^\ell_X =\frac{1}{2}{\vec C}^{-1}\frac{\partial\vec C}{\partial C_\ell^X}\vec C^{-1}, \label{eq:QMLelle} \end{equation}$ (41)where C = S(C_ℓ) + N is the (signal plus noise) covariance matrix of the map, m. C_ℓ is a fiducial power spectrum needed for the calculation of the signal part of the covariance. In this paper we will take it to be given by the true power spectrum as used to produce the simulated skies. Though this would be an unfair assumption, while testing the performance of this power spectrum estimation technique, it is justified in our case, where the fact that it leads to the minimal estimation uncertainties increases the power of our test.

Gruppuso et al. (2009) describe a specific implementation of the QML method, nicknamed Bolpol, used in this work and discusses its performance in the application to the WMAP 5 year data. The QML estimator is known to be equivalent to a single iteration of a quasi-Newton-Raphson procedure to search for the true likelihood maximum (Bond et al. 1998).

Table 2

Frequently used symbols in this paper.

4.2. Noise bias

In the mapmaking methods considered in this paper residual noise in the maps is independent of the sky properties and completely defined by the time-domain noise properties and the scanning strategy. The noise present in the maps contributes to the power spectrum estimates of the map signal. We will therefore refer to this contribution as the noise bias and use it to quantify the noise level expected in the maps of our different methods in a manner more succinct and manageable than the full noise covariance.

The noise bias is defined to be the expectation value for angular power spectrum in the noise-only case $N_{ℓ}^{XY} = ⟨ {C_{ℓ}^{XY, noise}}^{⟩},$ $\begin{equation} \label{eq:noise_bias} N_\ell^{XY} = \left\langle C^{XY,\mathrm{noise}}_\ell\right\rangle, \end{equation}$ (42)where X and Y stand for T, E and B. In the general case of a quadratic estimator, such as QML, Eq. (40), the power spectrum estimates are given as a quadratic form of the input map, m, $C_{ℓ}^{XY} = m T Q \begin{matrix} XY \\ ℓ \end{matrix} m = tr Q \begin{matrix} XY \\ ℓ \end{matrix} m m T .$ $\begin{equation} \label{eq:qml_bias} C_\ell^{XY} = \vec m^\mathrm T \vec Q_\ell^{XY}\vec m = \mathrm{tr}\:\vec Q_\ell^{XY}\vec m\vec m^\mathrm T. \end{equation}$ (43)By taking m to be noise-only, the noise bias can be expressed as $N_{ℓ}^{XY} = ⟨ {C_{ℓ}^{XY, noise}}^{⟩} = tr Q \begin{matrix} XY \\ ℓ \end{matrix} N .$ $\begin{equation} \label{eq:noise_bias1} N_\ell^{XY} = \left\langle C^{XY,\mathrm{noise}}_\ell\right\rangle = \mathrm{tr}\:\vec Q_\ell^{XY}\vec N. \end{equation}$ (44)Given the eigenvalue decomposition of the noise covariance matrix, $\vec N = \sum_i\lambda_i\hat{\vec e}_i\hat{\vec e}_i^\mathrm T$ , we can evaluate the noise bias as $N_{ℓ}^{XY} = \sum_{i} λ_{i} \hat{e} \begin{matrix} T \\ i \end{matrix} Q \begin{matrix} XY \\ ℓ \end{matrix} \hat{e} i .$ $\begin{equation} \label{eq:noise_bias2} N_\ell^{XY} = \sum_i\lambda_i\hat{\vec e}_i^\mathrm T\vec Q_\ell^{XY}\hat{\vec e}_i. \end{equation}$ (45)In the following we will validate the noise covariance matrices, estimated for the recovered sky maps, with the help of their respective noise biases, which we will compare to results of Monte Carlo simulations. In this context it is particularly useful to consider a pseudo-C_ℓ estimator that assumes uniform pixel weights and full sky. For this estimator, the operator $\vec Q_\ell^{XY}$ is $Q \begin{matrix} XY \\ ℓ \end{matrix} = \frac{1}{2 ℓ + 1} {(Y {\begin{matrix} X \\ ℓ \end{matrix}}^{)}}^{†} Y \begin{matrix} Y \\ ℓ \end{matrix},$ $\begin{equation} \label{eq:pseudo_cl} \vec{Q}_\ell^{XY} = \frac{1}{2\ell+1} \left(\vec{Y}_\ell^X\right)^{\dagger}\vec{Y}_\ell^Y,\rm{} \end{equation}$ (46)where $\vec Y_\ell^X$ has 2ℓ + 1 rows that are maps of the appropriate spherical harmonics. It maps a map vector into a vector of spherical harmonic expansion coefficients $\{a_{\ell m}^X\}$ where m = − ℓ...ℓ.

The procedure we implement here involves two steps. First for every estimated noise covariance we compute analytically the noise bias using Eqs. (45) and (46). Second, we compute the bias using Monte Carlo realizations of the noise-only maps produced using the corresponding mapmaking procedure. The multiplications $\vec Y_\ell^X\vec u_i$ are conveniently implemented using the HEALPix (Górski et al. 2005) Fortran 90 subroutine map2alm. We note that as we consider hereafter only full-sky cases, the noise biases we compute as described above would be equal to those expected in the Maximum Likelihood estimates of the sky power spectrum, were it not for the imperfection of the sky quadrature due to pixelization effects (Górski et al. 2005).

5. Simulation

5.1. Scanning strategy

In this study the Planck satellite orbits around the second Lagrangian point (L2) of the Earth-Sun system (Dupac & Tauber 2005). The spin axis lies near the ecliptic plane, precessing with a six month period around the anti-Sun direction at an angle of $7\fdg 5$ . The focal plane line-of-sight forms an 85° angle with the spin axis. In addition to these modes, we include a nutation of the spin axis and slight variations to the 1rpm spin rate. Details of the scanning simulation can be found from Ashdown et al. (2009) where it was used in a mapmaking study.

5.2. Planck detectors

In this paper we study residual noise in the Planck 70 GHz frequency maps. Planck has twelve detectors at 70 GHz. In the focal plane they are located behind six horn antennas, a pair of detectors (“Side” and “Main” detectors⁵) sharing a horn. A pair of detectors measures two orthogonal linear polarizations. The horns are split in two groups (three horns in a group). The Side and Main polarization sensitive axes of a group are nearly aligned and the polarization directions of the second group differ from the first group nominally by 45°. Two horns from the different groups make a polarization pair that follows the same scan path in the sky (three pairs in total with slightly different scan paths). As a minimum the observations of a polarization pair are required to build a polarization map. Due to implementation restrictions the Side and Main polarization axes are not fully orthogonal and the polarization direction differences between the groups are not exactly 45°, but the deviations from these nominal values are small ( ≲ $0\fdg 2$ ). The Side polarization axes of the two groups differ by $+22\fdg 5$ and $-22\fdg 5$ from the scan direction.

The beams of the detectors were assumed circularly symmetric with a 14′FWHM beam width. The beams do not impact the residual noise maps or covariance matrices. None of the mapmaking methods studied here make an attempt to correct for beam effects in the maps.

5.3. Time-ordered data

We computed the NCM’s of our three mapmaking methods using our noise model spectrum and the one year pointing data of twelve 70GHz detectors. We produced NCM’s for both N_side = 8 and N_side = 32 pixel sizes. We wanted to compare these NCM’s to the noise maps made by the same mapmaking methods. For that purpose we simulated 50 noise-only timelines and made maps from them. Our correlated noise streams were simulated in six day chunks by inverse Fourier transforming realizations of the noise spectrum (Natoli et al. 2002). We assumed an independence between the chunks and between the detectors. Figure 1 contains a comparison between the power spectra of the generated noise streams and the model spectra.

Twenty-five of the surveys featured a relatively high 1 / f contribution having the knee frequency, f_knee, set to 50mHz. The other half was simulated to have a more favorable f_knee = 10mHz. It should be noted that these frequencies have been chosen above and below the satellite spin frequency, 1rpm ≈ 17mHz. The slope of the 1 / f noise power spectrum was α = − 1.7. The correlation timescale of the 1 / f noise was restricted to about one day. This made our noise spectrum flat at low frequencies (below a minimum frequency f_min = 1.15 × 10^-5Hzt1 / 24h). As we described earlier, it is the minimum frequency that determines the correlation length of the noise filter in the optimal mapmaking. In the noise covariance matrix of the generalized destriping, the baseline correlation length is, however, determined by the knee frequency.

We used a uniform white noise equivalent temperature (NET) of $204\,\mu\mathrm K\sqrt{\mathrm s}$ for all detectors⁶. We chose this NET because we wanted to produce noise maps and covariance matrices whose noise levels are compatible with another CTP study (Jaffe et al., in prep.). In all mapmaking and NCM computations we assumed a perfect knowledge of the detector noise spectrum.

The noise timelines were processed directly into both low-resolution (N_side = 32) and high-resolution (N_side = 1024) HEALPix maps using the discussed mapmaking codes. In the N_side = 1024 temperature and polarization maps the mean standard deviations of white noise per map pixel were 44 and 63μK (Rayleigh-Jeans μK). For N_side = 32 maps the corresponding values were 1.4 and 2.0μK.

The high-resolution maps were in turn downgraded to the low target resolution using the schemes detailed in Sect. 2.3.

For the signal error studies described in Sect. 6.3, we scanned simulated foreground maps into signal-only timelines. These we processed into low-resolution (N_side = 8) maps using the same methods as for the N_side = 32 (both directly and through high resolution). We then extracted the signal error part by subtracting a binned map from the destriped map. The foreground signal errors were summed with a CMB map to provide a worst case scenario of signal striping in otherwise perfectly separated CMB map.

5.4. Input maps

To study bandwidth limitation with respect to downgrading we simulated 117 high-resolution N_side = 1024 CMB skies corresponding to the same theoretical spectrum, C_ℓ. These maps were smoothed and downgraded to N_side = 8 using three different Gaussian beams of widths 5°, 10°, 20° and three apodized step functions with the choices of (ℓ₁,ℓ₂) being (20,24), (16,24) and (16,20). A seventh set of downgraded maps was produced by noise weighting according to the scanning strategy. To comply with this last case, the smoothing windows include also the N_side = 8 pixel window function from the HEALPix package.

For the signal error exercise we used the Planck sky model, PSM⁷ version 1.6.3, to simulate the full microwave sky at 70GHz. For diffuse galactic emissions we included thermal and spinning dust, free-free and synchrotron emissions. We then added a Sunyaev-Zeldovich map and finally completed the sky with radio and infrared point sources. The combined N_side = 2048 map was smoothed with a symmetric Gaussian beam and scanned into a timeline according to the scanning strategy.

For final validation the noise covariance matrices were tested in power spectrum estimation. Each noise map was added to a random CMB map drawn from the theoretical distribution defined by a fixed theoretical CMB spectrum. The theoretical spectrum is the WMAP first-year best fit spectrum and has zero BB mode.

All maps in this work are presented in the ecliptic coordinate system. This choice is useful for Planck analysis since the scanning circles and many mapmaking artifacts form circles that connect the polar regions of the map. In this coordinate system the galaxy is not positioned at the equator but rather forms a downward opening horse shoe shape around the center of the map.

6. Results

In this section we first focus on the noise covariance matrices computed using different mapmaking techniques. We discuss and compare the overall noise patterns implied by such matrices and test the quality of the destriper approximation as applied to the noise covariance predictions. In the second part of the section we discuss the low-resolution maps and evaluate their quality in the light of their future potential applications, such as those to the large angular scale power spectrum estimation. Due to the computational resource restrictions the low-resolution results presented here are obtained either with the HEALPix N_side = 8 – in particular tests in Sect. 6.3 and power spectrum estimation tests in Sect. 6.2.3 – or N_side = 32 – in most of the other sections.

6.1. Noise covariance matrices

First we discuss the noise covariance matrices computed for the low-resolution maps of the direct method. As explained in Sect. 2.3 we compute from these matrices the noise covariances of the other downgrading techniques. In the following we consider noise covariance calculated using 4 different ways. In the first way we compute the noise covariance using the optimal algorithm. For this purpose we have developed two codes MADping or ROMA, which are described in Sect. 3.1. However, as they are just two different implementations of the same algorithm, we derive most of the results presented in the following and involving the optimal covariance using MADping. We note that whenever results from the both codes are available they have turned out to be virtually identical within the numerical precision expected from this kind of calculations. The optimal noise covariance matrices are expected to provide an accurate description of the noise level found in the actual optimal maps. We will test this expectation in the following and use the optimal results as a reference with which to compare the destriping results.

Fig. 2

Eigenspectra of the inverse covariance matrices N^-1. MADping, ROMA and Madam results for f_knee = 50mHz overlap completely. Springtide results are for 10mHz.

The three remaining computations of the noise covariance are based on the destriping approach and correspond to different assumptions about the offset prior as well as baseline length. We consider the following specific cases: a classical destriper calculation with a baseline of 3600s (Springtide) and two generalized destriper computations with a baseline of 1.25 s and 60s (Madam). For each of these cases we will compare the covariance matrices with each other, with the optimal covariances and then test their consistency with the noise found in the simulated maps. We note again that this last property is not any more ensured given the approximate character of the destriper approach.

Figure 2 shows the eigenvalue spectra of some inverse NCMs. We note that all matrices possess a positive semi-definite eigenspectrum as is required for any covariance matrix, yet at the same time they all have one nearly ill-conditioned eigenmode⁸, which renders the condition number, i.e., the ratio of the largest and smallest eigenvalue, very large. This is in agreement with our expectations as described in Sect. A.3. Indeed the peculiar eigenmodes corresponding to the smallest eigenvalues of the inverse matrices are also found to be non-zero and constant for the I part of the vector and nearly zero for its polarized components, and thus close to the global offset vector discussed in Sect. A.3. The small deviations, on the order of 10^-3, exist, as expected, as none of the peculiar modes is in fact truly singular. We note that the MADping and ROMA results, both computed in this test, are seen in the figure to be indistinguishable. They also coincide very closely with the Madam results computed for the same f_knee = 50 mHz. The Springtide results, computed with f_knee = 10 mHz are close to, though not identical with, the Madam results for the very same value of f_knee when a longer (60s) baseline is used for Madam.

Figure 3 depicts the estimated Stokes I, Q and U pixel variances as well as the covariance between them. These quantities are dominated by the white noise contribution and all methods describe white noise in the same manner. The top right-most panel shows the reciprocal condition number (1/condition number) of the 3 × 3 blocks of the matrix $\vec A^\mathrm T\mathcal N^{-1}_\mathrm u\vec A$ for each of the sky pixels. These numbers define our ability to disentangle the three Stokes parameters for each of the pixels. Whenever they are equal to 1 / 2 not only can the parameters be determined but their uncertainties will not be correlated. If the reciprocal condition number for a selected pixel approaches 0, the Stokes parameters can not be constrained. In the cases considered here, the Stokes parameters can clearly be determined for all the pixels.

Figure 3 shows a strong asymmetry between the IQ and IU blocks. The fact that the polarization axes of the two detectors of a horn are not fully perpendicular makes the I noise of a pixel correlate with the Q and U noises of the same pixel. We can imagine an instrument basis, where one group of three horns measures the Q of this basis and the other group measures the U of the same basis. Because Side and Main polarization axis deviations from the orthogonality have similar magnitudes in the two groups, we expect the diagonals of IQ and IU blocks to be similar (symmetric in IQ and IU) in the instrument basis. In the map we use a different polarization basis⁹. Building the noise weighted map (in the mapmaking) rotates the Q and U from the instrument basis to the map basis. The IQ and IU blocks become asymmetric in this rotation.

Fig. 3

Top: MADping pixel variances for temperature and polarization and the reciprocal condition number of the pixel observation matrices. Bottom: correlation coefficient × 10³ between I − Q, I − U and Q − U pixels. This part of the noise covariance is dominated by white noise which is modeled equivalently in all three paradigms. Hence, MADmap, ROMA, Madam and Springtide results are nearly identical. Maps are rotated into galactic coordinates to show the structure near the ecliptic poles.

Fig. 4

A single column of the MADping covariance matrix corresponding to a pixel at the ecliptic equator. For both ROMA and Madam 1.25s baseline counterparts all visible characteristics remain unchanged. We plot the value of the correlation coefficient, R, multiplied by 10³. In order to enhance the features, we have halved the range of the color scale.

Figures 4, 5 show plots of a single column of the noise covariance matrix. The column corresponds to reference pixel number 0. In the HEALPix nested pixelization scheme for N_side = 32 resolution, pixel 0 is located at the equator. In the plots each pixel has the value ⟨ m_pm_q ⟩ normalized by $\sqrt{\langle m_p^2 \rangle \langle m_q^2 \rangle}$ . Thus the pixel values of the plots represent correlation coefficients. Due to this normalization, the reference pixel automatically gets unit value and is later set to zero in order to bring out smaller features of the other elements of the columns.

Planck scans the sky in nearly great circles that connect the ecliptic polar regions and center around the ecliptic plane. The displacement of the spin axis from the ecliptic plane depends on the phase of the satellite $7\fdg 5$ precession. The NCM column maps are characterized by bands of correlation along the scanning rings. Pixels near the equator, such as the reference pixel 0, are only observed during a few-hour window as the satellite scanning ring is rotated over the course of the survey. The two crossing bands of higher correlation correspond to two pointing periods half a year a part that observe the reference pixel 0.

For both generalized destriping and optimal mapmaking, there is a visible gradient in the correlation along the scanning ring. Pixels that are observed immediately before or after the reference pixel have the strongest correlation. The conventional destriping with its hour long baselines assumes constant correlated noise over the scanning ring and does not, therefore, show this feature.

Figures 4 and 5 show that the strongest cross-correlations (~1%) exist between Q and U noise maps. Side and Main polarization sensitive directions differ slightly from 90° and, as a result, small (~10^-4%) IQ and IU correlations remain.

Fig. 5

A single column of the Springtide covariance matrix corresponding to a pixel near the ecliptic equator. For description of the normalization, see text.

Figures 6, 7 show plots of the NCM columns of another reference pixel (number 2047). This pixel is located at the northern ecliptic pole region and it exhibits a very different correlation pattern compared to the previous case. Since the pole is visited frequently through the course of the survey, it becomes correlated with the rest of the map as a whole. Correlation amplitude is increased by a factor of 3 from the equator pixel case (the increase can be seen from the color bar ranges) and there is now a distinct asymmetry between northern and southern hemispheres. As one expects, the asymmetry only appears in optimal and generalized destriping estimates.

Fig. 6

A single column of the MADping covariance matrix corresponding to a pixel at the ecliptic pole. For both ROMA and Madam 1.25s baseline counterparts all visible characteristics remain unchanged. We plot the value of the correlation coefficient, R, multiplied by 10³. In order to enhance the features, we have halved the range of the color scale.

Fig. 7

A single column of the Springtide covariance matrix corresponding to a pixel near the ecliptic north pole. For description of the normalization, see text.

6.2. Noise covariance validation

In this section we report the results of three different validation tests. We performed a χ² test, compared the noise biases computed from the matrices and corresponding Monte Carlo maps and finally used the matrices and Monte Carlo maps as inputs to angular power spectrum estimation.

6.2.1. By χ²

Residual noise expected in the recovered maps is Gaussian due to the linear character of all the mapmaking methods considered here. Thus the noise can be completely described by its covariance matrix. More specifically, in the absence of any singular modes of the estimated residual covariance, N, the residual maps, $\vec m = \hat{\vec s} - \vec s$ , are drawn from a multivariate Gaussian distribution defined by N. Therefore, the χ² statistic, defined as χ² = m^TN^-1m, is drawn from a χ² distribution with 3N_pix degrees of freedom (d.o.f.).

If any singular mode is present we simply replace the matrix, N^-1, in the definition of χ², by a matrix, N^′ − 1, which is like N^-1 in all respects but has the eigenvalue corresponding to the singular vector set to 0 (see Sect. 3.1). We note that if the eigenvalue decomposition of the matrix N^-1 is not available or too costly to compute, we can achieve numerically the same effect by defining N′ = N + η²vv^T, where v is a singular vector we want to project out and η² is a large positive number for which however inversion of N′ is still stable (e.g. Bond et al. 2000). As we subtract one degree of freedom corresponding to the excluded, ill-conditioned eigenmode we expect that there are in total 3N_pix − 1 degrees of freedom left in our maps

The χ² test can be applied, with certain restrictions, to the smoothed maps and matrices. In general, a smoothed matrix will not have any longer the original number of degrees of freedom: some of the small scale eigenmodes have become ill-conditioned (see Sect. 3.3). In this section we regularize the smoothed matrices by omitting the singular eigenmodes.

Singular modes of the original matrix require special attention when smoothing and need to be treated explicitly if the smoothed version of the ill-conditioned eigenmode, v, does not belong to the null space of the inverse smoothed NCM, i.e., $\vec {\tilde N}^{-1}\left(\vec L\vec v\right) \nsimeq 0$ . To do so, we employ the same approach as before, replacing the regularized inverse of the smoothed covariance matrix, $\vec {\tilde N}^{-1}$ , by, $\begin{matrix} {\begin{matrix} ˜ \\ N \end{matrix}}^{-1} \to {[\begin{matrix} ˜ \\ N \end{matrix} + η^{2} L v {(L v)}^{T}]}^{-1} \overset{η^{2} \to \infty}{- \to} \begin{matrix} ˜ \\ N \end{matrix} -1 \\ - (\begin{matrix} ˜ \\ N \end{matrix} -1 L v) {[{(L v)}^{t} \begin{matrix} ˜ \\ N \end{matrix} -1 {(L v)}^{]}}^{-1} {(\begin{matrix} ˜ \\ N \end{matrix} -1 L v)}^{T} \end{matrix}$ $\begin{eqnarray} &&\tilde{\vec N}^{-1} \rightarrow \left[ \vec {\tilde N} +\eta^2\,\vec L \vec v \left(\vec L \vec v\right)^\mathrm T \right]^{-1} \stackrel{\eta^2 \rightarrow \infty}{\longrightarrow} \vec {\tilde N}^{-1} \nonumber\\&&\qquad \quad- \left( \vec {\tilde N}^{-1}\,\vec L \vec v \right) \left[ \left(\vec L \vec v\right)^t \vec {\tilde N}^{-1} \left(\vec L\vec v\right) \right]^{-1} \left( \vec {\tilde N}^{-1}\,\vec L \vec v \right)^\mathrm T \end{eqnarray}$ (47)\arraycolsep1.75ptwhere L is a smoothing operator, Eq. (26), and the last expression follows from the Sherman-Morrison-Woodbury formula (Woodbury 1950). This last operation additionally reduces the number of degrees of freedom by 1 (or whatever number of modes, v is to be projected out).

Fig. 8

Empirical χ² distribution functions from the 25 residual noise maps compared with the theoretical cumulative probability density. The black stair line is the empirical distribution function, the blue solid line is the theoretical χ² distribution for 3N_pix − 1 = 36863 degrees of freedom. The red dashed line is the least squares fit of the χ² distribution to the experimental distribution (d.o.f. being the fitting parameter). The horizontal axis is translated to the expected center of the distribution, ⟨ χ² ⟩ = d.o.f. = 36,863, and scaled by the expected deviation, $\sigma_{\chi^2}=\sqrt{2\mathrm{d.o.f.}}$

The Kolmogorov-Smirnov test can be used to test whether a set of samples conforms to some theoretical distribution. The test estimates the probability of the maximal difference between the empirical distribution function, $F_{n} (x) = \frac{1}{n} \sum_{i = 1}^{n} Θ (x_{i}),$ $\begin{equation} \label{eq:edf} F_n(x) = \frac{1}{n}\sum_{i=1}^n \Theta(x_i), \end{equation}$ (48)of the observations x_i (in our case the individual χ²) and the theoretical cumulative distribution function. Θ(x) is the Heaviside step function. We note that in this work we take an advantage of the fact that we can simulate the residual noise directly. Though this is clearly not the case when real data are considered, the tests described here can be applied to a difference of two sky maps produced by disjoint sets of detectors operating at the same frequency and can therefore be a useful test of the real life data processing integrity (e.g. Stompor et al. 2002).

Figure 8 shows the theoretical and Monte Carlo cumulative distribution functions for 25 noise maps in the case of the direct method. Reported p-values are the probabilities of observing this level of disagreement even if the noise description was exact. Conventionally, the level p < 0.05 is considered to be enough to reject the null hypothesis that the distributions match.

In the case of the direct method maps our results show that only optimal (noise) maps and their respective noise covariance are mutually consistent in the light of the χ² statistics. The good statistical agreement in this case does not depend on the time domain noise characteristics nor map resolution. This is expected given that the noise covariance estimator implemented in the optimal codes, Eq. (18), is an exact expression describing the noise properties in the pixel domain, and that we have assumed perfect knowledge of the time domain noise.

The level of consistency found in the destriping cases varies depending on the underlying time domain properties, i.e., f_knee, and on the assumed baseline length. In the case of f_knee = 50mHz, we have not found a satisfactory agreement in any of the considered destriping cases. For the lower f_knee = 10mHz the results obtained with the generalized destriper, Madam, are satisfactory for the short, 1.25s, baseline choice, and marginal for the long, 60s, choice.

Destriper performance is better using short baselines and at lower knee frequencies. Thus it is evident that failures of the destriping approach to model the residual noise are due to correlated noise that cannot be modeled by baselines of the chosen length. This unmodeled noise is a concern whenever the baseline frequency is lower or comparable to the knee frequency.

But where does the unmodeled noise show? Our three choices of the baseline length (3600s, 60s and 1.25s) correspond to baseline frequencies 0.28mHz, 17mHz and 800mHz. These baseline frequencies provide lower limits to the frequencies where unmodeled noise will manifest. If the lower limit is higher than the satellite spin frequency, 1 / min = 17mHz, it will begin to limit the angular scales at which the unmodeled noise can manifest. Translated to angular units, the 800mHz threshold is $7\fdg 5$ , multipole ℓ = 24 or size of a single N_side = 8 pixel. For short baselines or low knee frequencies the unmodeled noise manifests as a relatively small angular scale correlation and does not bias the power spectrum estimates at low ℓ.

Fig. 9

Madam empirical χ² distribution functions from 60 residual noise maps compared with the theoretical cumulative probability density. For a smoothed map, we count the degrees of freedom as the number of included eigenmodes in the inverse NCM. Top: three sets of low-resolutions maps using just one of the low-resolution methods at a time. The high-resolution maps for INW and smoothing methods had N_side = 1024. Bottom: since it is suboptimal to compute a low-resolution spherical harmonic expansion from a noisy high-resolution map, we test how well the smoothing approach works in conjunction with the two pixel-based downgrading methods.

Fig. 10

Analytical and mean Monte Carlo noise biases. Gray band is the 1-σ region for the average, computed by dividing the sample variance by $\sqrt{25}$ . Top row: Madam runs at f_knee = 50 mHz using a short 1.25s baseline. Middle row: Madam runs at f_knee = 50 mHz using a relatively long 60s baseline. Bottom row: springtide runs at f_knee = 10 mHz.

Fig. 11

Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 10mHz knee frequency.

We then proceed to study the agreement between downgraded noise maps and the noise covariance matrices. As a test case, we use the Madam NCM for 10mHz knee frequency, 1.25s baselines. We smoothed the covariance matrix using an apodized window function, setting the thresholds to 2N_side and 3N_side respectively. As expected, the smoothed matrix is extremely singular. We compute its inverse by including only the eigenvalues that are greater than 10^-2 times the largest eigenvalue, including 20882 of the 36864 available modes.

Figure 9 shows the empirical distribution functions of the χ². Even though the matrix is computed for the direct method, the inverse noise weighted (INW) maps conform well to it. However, when we apply the smoothing kernel to the high-resolution maps, there is a clear disagreement. This stems from the fact that in this downgrading the high-resolution pixels are not correctly (inverse noise variance) weighted when we compute the low-resolution map. If we first produce a low-resolution direct method or INW map and then smooth it, the agreement is much better. This is shown in the bottom row of Fig. 9.

We conclude that the noise covariance of the low-resolution maps produced by the destriping algorithm needs to be used with care. In the case of low knee frequency, the generalized destriper results are comparable with the optimal results, but we emphasize that, if the accuracy of the noise description in pixel domain is a major concern, then only the optimal techniques are suitable. In the next sections we will reconsider all these low-resolution mapmaking techniques in a context more specific to the large angular scale power spectrum estimation work, which is envisaged as the main application of the low-resolution maps and their covariances.

6.2.2. By noise bias

In this section we describe the calculation of the average angular power spectrum of noise maps, i.e., the noise bias – see Sect. 4.2, using a pseudo C_ℓ estimator for which both the NCM estimate and the Monte Carlo map averages are feasible to compute. Testing the noise covariance matrix by comparing estimated and measured noise biases can be viewed as complementary to the χ² tests described in the previous section. It can certainly provide more information than the plain χ² test, as instead of simple pass or fail indicator, the noise bias comparison will tell at which angular resolution the noise model agrees with the data. At the same time, the noise bias is less sensitive to the anisotropic features, in this case mainly ring-like structures, present in the residual noise. The noise bias test is more directly relevant for power spectrum estimation.

Figure 10 compares noise bias averages from 25 noise realizations of the maps of the noise residuals to the analytical estimates, Eq. (45), based on the estimated noise covariance matrix. Map spectra are computed using the HEALPix anafast utility and the estimated noise biases using the corresponding map2alm subroutine. We only show the auto-spectra, TT, EE and BB, as the scanning has decoupled the modes to large extent and only minimal coupling between the modes exist.

The error band around the averages is the standard deviation of the individual C_ℓ values divided by $\sqrt{25}$ . Each plot exhibits up to five curves: direct, inverse noise weighted and harmonic smoothed noise biases, and two analytical estimates.

Top row of Fig. 10 displays the noise bias comparison for the high 50mHz knee frequency with short 1.25s Madam baselines. The earlier χ² plot, Fig. 8, showed a significant disagreement between the maps and the matrix used to derive these biases. In harmonic space the correspondence is significantly better. Since the angular power spectrum does not contain information about the orientation of structures, it may be concluded that the disagreement averages out in the compression.

In contrast with the excellent short baseline results, middle row of Fig. 10 displays an apparent failure of the long Madam 60s baselines to model the high f_knee = 50mHz noise. To no surprise, the unmodeled noise above the low 17mHz baseline frequency has passed through the mapmaking process and projected onto the noise maps. As our baseline approximation of the noise is now worse than in Fig. 10, the measured noise bias is higher even though the estimated bias is lower.

Between these two extremes we find the low 10mHz knee frequency cases. For these time lines the level of correlated noise is so low that even extremely long hour baselines do an adequate job estimating the noise bias shown in bottom row of Fig. 10.

At 10mHz the Madam results for long 60s baselines are equivalent with Springtide results: both succeed to estimate noise bias at low knee frequency but should not be used for high knee frequencies as such.

We note that our analytic predictions are well in line with WMAP findings (Hinshaw et al. 2007; Page et al. 2007) and studies of the destriping framework (Efstathiou 2005, 2007).

Our noise bias estimates for EE and BB spectra are equal, but the averaged spectra for the Monte Carlo maps appear visually different in this respect. To ensure that this is only due to Monte Carlo noise we ran Madam in Monte Carlo mode, simulating noise on the fly and avoiding the costs associated with storing of the time-ordered data. After averaging over 117 N_side = 8 noise map spectra we found that the differences between EE and BB noise map spectra were at most 10%.

In the Appendix A.1 we re-plot some of these figures after subtracting the analytical bias and dividing by the Monte Carlo sample deviation to highlight the differences between the analytical and numerical results.

6.2.3. By power spectrum estimation

Our final validation procedure for the noise covariance matrices was to use them in C_ℓ estimation. Due to resource constraints this exercise was conducted at a lower N_side = 8 resolution. All three mapmaking codes produced maps from the 25 noise realizations. Each realization was paired with an independent realization of the CMB sky and the co-added maps were processed using the Bolpol code, an implementation of the QML estimator described in Sect. 4.1. The 25 power spectrum estimates for each multipole were then averaged over and the Bolpol-determined error bars were accordingly divided by $\sqrt{25}$ .

An example of the results is shown in Fig. 11. These estimates were obtained for the case with the low knee frequency, f_knee = 10mHz, using the conventional destriper, Springtide. We note they do not hint unambiguously at any problem with the estimated covariance, even if the χ² (strongly) and the noise bias (mildly) tests may indicate otherwise. This is likely in part due to a lower sensitivity of the power spectrum test on the one hand and on the other due to the fact that the lower resolution has been used in this last case.

Similar statistically good agreements can also be seen in the case of the higher value of f_knee, if the covariance is computed using either the optimal or generalized destriping technique with the short baselines of 1.25 s. If longer baselines are used, i.e., 60 s, the estimates of the polarized spectra, both E and B, are visibly discrepant with the assumed inputs. Similar disagreement can be seen if the off-diagonal elements of the covariance matrices are neglected. In both of these cases, no particular effect on the total intensity spectrum can be noticed in the range of investigated angular scales. We illustrate all these statements in Appendix A.2. These observations emphasize the importance of precise estimation of the noise covariance in particular for the polarized power spectra.

6.3. Low-resolution maps

In the previous sections we have discussed our ability to estimate correctly the properties of the noise present in the low-resolution maps. We have demonstrated that this is indeed the case for all considered resolution downgrading strategies and both optimal and destriper maps. Though in the latter case a baseline length needs to be carefully chosen depending on the time domain noise characteristic.

In this section we focus on the low-resolution maps themselves. We will look at them from three different perspectives, evaluating the level of the mapmaking artifacts left in the maps, the properties of the sky signal and the level of the noise.

In Fig. 12 we show the differences of the noise-free low-resolution maps computed using different downgrading approaches discussed earlier and the input map used for the simulations. The reference low-resolution version of the input map has been obtained via simple binning of the sky signal directly into low-resolution pixels.

In the case of the direct method we see clearly the extra power spread all over the sky in all three Stokes parameters. Any other proposed approach clearly fares much better leading to a substantial decrease of the level of the observed artifacts. This is quantified with the help of the pseudo-spectra in Fig. 13 and in Table 3, where we have collected the root mean square estimates for the signal residual maps.

Fig. 12

Examples of signal striping. We show the difference between a binned and destriped signal-only maps. Rows correspond to direct Madam results, 1.25s and 60s, and inverse noise weighted 1.25s baseline case respectively.

Table 3

Comparison of TT power (rms) in the stripe maps at N_side = 8.

The CMB part of the low-resolution maps also depends on the downgrading technique. These attempt to suppress the high-ℓ (sub pixel) power and therefore can potentially affect the CMB angular power spectrum even within the band of interest. Figure 15 shows the full-sky pseudo-C_ℓ spectra averaged over 117 CMB realizations, downgraded using inverse noise weighting and a number of smoothing kernels. In the QML method the situation is slightly more complex, as the quadratic map function is multiplied with the inverse Fisher matrix to produce the QML power spectrum estimate (see Sect. 4.1). The inverse Fisher matrix can correct some of the aliasing effects that cause bias in the power spectra. To simulate the effect of the inverse Fisher in QML we deconvolved our pseudo spectra with a mode coupling kernel that we computed from a map of ones.

Comparison of the spectra shows that it is hard to attain subpercent bias even at ℓ = 2N_side. If these were estimates from an actual power spectrum estimation code with correlated noise and a sky cut, the smoothed covariances would, however, be regularized by adding white noise effectively leading to a considerable uncertainty already at that multipole due to white noise and a sky cut. Methods that produce less than 5% bias for C_ℓ estimates at ℓ = 2.5N_side are the Gaussian 10° symmetric beam and the apodized step function for (ℓ₁,ℓ₂) = (2or2.5N_side,3N_side). The noise weighting and direct low-resolution mapmaking produce similar aliasing effects.

In Fig. 14 we show the actual sky signal spectra estimated for the low-resolution, noiseless maps. It can be seen that the estimated band powers are not drastically affected with respect to the estimates coming from the binned maps. Nevertheless, the case of estimates coming from the direct low-resolution maps show some deviation from the binned case. The bias is most prominent in the BB low multipole estimates. We note that since the signal covariance matrix is ill-conditioned it was regularized by adding a small white NCM (σ ≃ 1μK) and each signal map received a noise realization consistent with this white NCM.

Fig. 13

Comparison of CMB and stripe power spectra reveals that the striping can significantly bias the EE and BB power spectrum estimates. Stripe map spectra are computed from maps shown in Fig. 12.

Fig. 14

Band power estimates of CMB with CMB and foreground stripes.

Fig. 15

Bandwidth limiting the signal using various window functions. The three panels show the same curves first as raw estimates, then after deconvolving the smoothing and pixel windows and finally after subtracting and dividing by input model. The apodized step window functions (“cosine”) correspond to choices of the thresholds (ℓ₁,ℓ₂) as (20,24), (16,24) and (16,20) respectively. Here, solid lines are for Gaussian windows and dashed lines for the apodized step functions.

Fig. 16

TT noise bias computed from smoothed covariance matrices. Solid lines correspond to the Gaussian window functions and the dashed ones to apodized step functions. Left: linear vertical scale. Right: logarithmic vertical scale.

None of the proposed downgrading approaches can yield a noise level better than the direct method. This is because noise-weighted downgrading or smoothing both introduce departures from the optimal weighting of the noise present in the data. The expected level of the noise is therefore an important metric with which to compare the different downgraded maps. Figure 16 shows the analytical noise biases evaluated from an N_side = 8 Madam noise covariance matrix after smoothing using the 6 different beam window functions defined in Sect. 5.4 and the unsmoothed case that is a close match to the noise weighted maps. The narrowest Gaussian window function having an FWHM equal to 5°, slightly less than average pixel width, appears pathological due to aliasing effects. The rest of the test cases are more stable but feature a non-negligible amount of aliased power for multipole moments beyond ℓ = 3N_side (the C_ℓ are not normalized with the conventional ℓ(ℓ + 1) / 2π).

6.4. Resource requirements

Table 4 lists some CPU time costs for various N_side, baseline length and knee frequency combinations. Being a considerably expensive operation, we have not tested the scaling of the optimal calculation for this particular exercise. See Borrill (1999) for more discussion of scaling. The Madam resource cost scales roughly linearly with respect to the number of pixels and the so-called baseline correlation length.

Resource requirements of the three approaches vary. As with mapmaking, the destriping problem size is related to the chosen baseline offset length. The same consideration applies also for noise covariance estimation. Both the Madam generalized destriper and the MADping covariance calculations scale by the length of the noise filter. ROMA does the calculation in Fourier space and as a results scales as the logarithm of the noise filter length. The MADping algorithm has no dependency (ignoring communication and final file writing) on the number of pixels being used. Although the computation prefactors are not specified, we can see that at N_side = 32, MADping is already more resource efficient than ROMA (cf. rows two and three).

Table 4

CPU time costs for 1270 GHz detector years = 3 × 10¹⁰ samples.

7. Conclusion

We have presented the formalism and tools to compute the residual noise covariance matrix for three mapmaking paradigms studied for Planck (an optimal method and two destriping methods). The structure of these matrices follows from the scanning strategy but is modulated by the underlying noise model that defines the mapmaking method. The matrices were tested against Monte Carlo noise maps that were processed from correlated noise streams into maps using MADmap, Madam and Springtide mapmaking codes.

The most accurate correspondence between the covariance matrix and the noise maps is, as expected, between the optimal mapmakers, MADmap and ROMA, and their covariance matrices and the two codes produce nearly identical matrices. Both the generalized (Madam) and classical (Springtide) destripers are shown to disregard some medium frequency correlated noise that cannot be modeled by the chosen baseline offset length. It is shown that for a low knee frequency, 10mHz, the Springtide baseline length of 1 h is sufficient to model the correlated noise and compute the residual noise covariance. For a high knee frequency, 50mHz, even the Madam 60s baselines are too long to suffice. However, using a short 1.25s baselines (just 96 samples) the Madam results are extremely close to optimal results even for the high knee frequency.

As a concluding test we used the matrices in actual power spectrum estimation and verified that all methods model residual noise adequately when the noise approximation (baseline length) is short enough to model the correlated noise.

Resource costs of the methods vary greatly. Although both MADping and ROMA arrive at the same result, the implementations differ and the ROMA result scales with the resolution of the map. Both optimal implementations are extremely resource intensive. The Madam method can be used to produce good approximations of the optimal covariance matrices at a fraction of their cost. The covariance matrices for these tests were evaluated for two low resolutions, N_side = 8 and N_side = 32. It is possible to compute the matrices up to N_side = 64 (already 162 gigabytes) or even up to N_side = 128 (2.6TB) but the computational scaling of the methods using the matrices will likely set limits to the usefulness of such resolutions.

We studied two classes of downgrading strategies, those that make an attempt to limit the signal bandwidth and those that do not. The choice of the best downgrading approach depends on both the accuracy of the resulting noise and signal models. First measuring our ability to compute an accurate noise covariance matrix and the second describing our ability to control signal effects such as striping and aliasing.

All methods to produce low-resolution maps have their drawbacks. Direct mapmaking at low resolution produces an unacceptable level of signal striping that is caused by subpixel structure. Downgrading by noise weighting biases the power spectrum through aliasing effects. The frequently used Gaussian beam smoothing has a significant drawback of suppressing the signal at otherwise useful angular resolutions. We find that an apodized step function is able to retain a great deal of signal power up to ℓ_max = 2N_side but even then the power spectrum estimates will be biased beyond ℓ = 2.5N_side. However, to accurately evaluate the noise covariance matrix for a smoothed map, we would need to compute an unsmoothed covariance matrix at the high map resolution and then apply the same smoothing kernel to both the matrix and the map. Disregarding this requirement leads to disagreement between the map and the matrix that can be alleviated by combining two or more of the downgrading methods.

We note that, although not considered in this work, most of the complications associated with low resolution mapmaking can be directly addressed by constructing the map in harmonic domain instead of pixel domain (Challinor et al. 2002; Armitage & Wandelt 2004; Armitage-Caplan & Wandelt 2009). A study of residual noise covariance with such methods would make an interesting extension to that presented here.

Of the downgrading methods considered, we consider smoothing, with a suitable choice of the window function and possibly a intermediate downgrading step by inverse noise weighting, to produce the best possible low-resolution maps for power spectrum analysis.

We presented in this work a method to compute the residual noise covariance of a smoothed, bandwidth-limited map. The method was shown to produce an accurate description of the noise in the smoothed maps when both the map and the matrix agree prior to smoothing. Our method of smoothing the covariance matrix makes it possible to consider bandwidth limited low-resolution maps and produce subpercent level unbiased power spectrum estimates up to ℓ = 2.5N_side.

In this work we have assumed a single frequency channel, uncorrelated noise between detectors, noise that is white at high frequencies and full sky coverage. In further work any of these constraints can be lifted. The only application that we used the covariance matrix was in power spectrum estimation. The downgrading methods that suit power spectrum analyses best may not be optimal for different low-resolution analysis, e.g., study of large scale topology. A relevant future direction to explore is the use of the covariance matrices as inputs in a likelihood code for cosmological parameters.

¹

http://lambda.gsfc.nasa.gov/

²

HEALPix resolution is reported as either N_side or r, related by N_side = 2^r. A HEALPix map contains $12\times N_\mathrm{side}^2$ pixels.

³

http://crd.lbl.gov/~borrill/cmb/madcap/

⁴

http://www.netlib.org/scalapack/scalapack_home.html

⁵

Side and Main refer to two detector branches downstream from the orthomode transducer that separates the two orthogonal linear polarizations.

⁶

$\sigma = \mathrm{NET}\sqrt{f_\mathrm{sample}}$ , where σ and f_sample were defined in Eq. (19). The 70GHz detectors had f_sample = 76.8Hz.

⁷

See http://www.apc.univ-paris7.fr/APC_CS/Recherche/Adamis/PSM/psky-en.php

⁸

A double precision (64 bit) matrix is numerically ill-conditioned when the condition number exceeds approximately 10¹² (Press 1992). For our matrices the situation is not as dire but the first eigenmode still deserves special attention.

⁹

In a HEALPix map the Stokes parameters Q and U at a point in the sky are defined in a (x, y, z) reference coordinate, where the x-axis is along the meridian and points to south, the y-axis is along the latitude and points to east, and the z-axis points to the sky (Górski et al. 2005).

¹⁰

Actually, the global offset mode is only nearly singular. This is due to our assumed noise power spectrum, which is finite at zero frequency. In a more realistic case the offsets of the stationary segments will be however unknown, corresponding to an infinite amplitude at zero frequency. The noise covariance in such cases should be therefore considered to be strictly singular with the global offset being the singular eigenvector. In such cases the noise weighting on the right term of the map making equation, Eq. (17), will force the offset of each of the stationary data segments to be strictly zero. This may not be a sufficiently good approximation in particular for short stationary time segments. This could, however, be alleviated by introducing the segment offsets as extra degrees of freedom contained in the vector, x, (e.g. Stompor et al. 2002).

Acknowledgments

The work reported in this paper was done by the CTP Working Group of the Planck Consortia. Planck is a mission of the European Space Agency. We acknowledge the use of the CAMB (http://camb.info) code for generating theoretical CMB spectra. Some of the results in this paper have been derived using the HEALPix package (Górski et al. 2005). We acknowledge the use of version 1.6.3 of the Planck sky model, prepared by the members of Planck Working Group 2. This work has made use of the Planck satellite simulation package (Level-S), which is assembled by the Max Planck Institute for Astrophysics Planck Analysis Centre (MPAC). This research used resources of the National Energy Research Scientific Computing Center, which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. We thank CSC (Finland) for computational resources. RK wishes to thank the Jenny and Antti Wihuri Foundation and the Väisälä Foundation for financial support. This work was supported by the Academy of Finland grants 121703 and 121962. HKS thanks Waldemar von Frenckells stiftelse and Magnus Ehrnrooth Foundation for financial support. We acknowledge support from ASI, contract Planck LFI Activity of Phase E2. RS acknowledges partial support of the European Commission Marie Curie IR Grant, MIRG-CT-2006-036614.

References

Armitage, C., & Wandelt, B. D. 2004, Phys. Rev. D, 70, 123007 [NASA ADS] [CrossRef] [Google Scholar]
Armitage-Caplan, C., & Wandelt, B. D. 2009, ApJS, 181, 533 [NASA ADS] [CrossRef] [Google Scholar]
Ashdown, M. A. J., Baccigalupi, C., Balbi, A., et al. 2007a, A&A, 471, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ashdown, M. A. J., Baccigalupi, C., Balbi, A., et al. 2007b, A&A, 467, 761 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Ashdown, M. A. J., Baccigalupi, C., Bartlett, J. G., et al. 2009, A&A, 493, 753 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Benabed, K., Cardoso, J. F., Prunet, S., & Hivon, E. 2009, MNRAS, 400, 219 [NASA ADS] [CrossRef] [Google Scholar]
Bond, J. R., Jaffe, A. H., & Knox, L. 1998, Phys. Rev. D, 57, 2117 [NASA ADS] [CrossRef] [Google Scholar]
Bond, J. R., Jaffe, A. H., & Knox, L. E. 2000, ApJ, 533, 19 [NASA ADS] [CrossRef] [Google Scholar]
Borrill, J. 1999, in Proceedings of the 5th European SGI/Cray MPP Workshop, Bologna, Italy [Google Scholar]
Cantalupo, C. M., Borrill, J. D., Jaffe, A. H., Kisner, T. S., & Stompor, R. 2010, ApJS, 187, 212 [NASA ADS] [CrossRef] [Google Scholar]
Challinor, A., Fosalba, P., Mortlock, D., et al. 2000, Phys. Rev. D, 62, 123002 [NASA ADS] [CrossRef] [Google Scholar]
Challinor, A. D., Mortlock, D. J., van Leeuwen, F., et al. 2002, MNRAS, 331, 994 [NASA ADS] [CrossRef] [Google Scholar]
de Bernardis, P., Ade, P. A. R., Bock, J. J., et al. 2000, Nature, 404, 955 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
de Gasperis, G., Balbi, A., Cabella, P., Natoli, P., & Vittorio, N. 2005, A&A, 436, 1159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Delabrouille, J. 1998, A&AS, 127, 555 [Google Scholar]
Dupac, X., & Tauber, J. 2005, A&A, 430, 363 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Efstathiou, G. 2005, MNRAS, 356, 1549 [NASA ADS] [CrossRef] [Google Scholar]
Efstathiou, G. 2007, MNRAS, 380, 1621 [NASA ADS] [CrossRef] [Google Scholar]
Efstathiou, G., Gratton, S., & Paci, F. 2009, MNRAS, 397, 1355 [NASA ADS] [CrossRef] [Google Scholar]
Górski, K. M., Banday, A. J., Bennett, C. L., et al. 1996, ApJ, 464, L11 [NASA ADS] [CrossRef] [Google Scholar]
Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]
Gruppuso, A., de Rosa, A., Cabella, P., et al. 2009, MNRAS, 400, 463 [NASA ADS] [CrossRef] [Google Scholar]
Hanany, S., Ade, P., Balbi, A., et al. 2000, ApJ, 545, L5 [NASA ADS] [CrossRef] [Google Scholar]
Hinshaw, G., Nolta, M. R., Bennett, C. L., et al. 2007, ApJS, 170, 288 [NASA ADS] [CrossRef] [Google Scholar]
Jarosik, N., Barnes, C., Bennett, C. L., et al. 2003, ApJS, 148, 29 [NASA ADS] [CrossRef] [Google Scholar]
Jarosik, N., Barnes, C., Greason, M. R., et al. 2007, ApJS, 170, 263 [NASA ADS] [CrossRef] [Google Scholar]
Keihänen, E., Kurki-Suonio, H., Poutanen, T., Maino, D., & Burigana, C. 2004, A&A, 428, 287 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Keihänen, E., Kurki-Suonio, H., & Poutanen, T. 2005, MNRAS, 360, 390 [NASA ADS] [CrossRef] [Google Scholar]
Keihänen, E., Keskitalo, R., Kurki-Suonio, H., Poutanen, T., & Sirviö, A.-S. 2010, A&A, 510, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Kuo, C. L., Ade, P. A. R., Bock, J. J., et al. 2004, ApJ, 600, 32 [NASA ADS] [CrossRef] [Google Scholar]
Kurki-Suonio, H., Keihänen, E., Keskitalo, R., et al. 2009, A&A, 506, 1511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Maino, D., Burigana, C., Górski, K. M., Mandolesi, N., & Bersanelli, M. 2002, A&A, 387, 356 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Masi, S., Ade, P. A. R., Bock, J. J., et al. 2006, A&A, 458, 687 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Natoli, P., de Gasperis, G., Gheller, C., & Vittorio, N. 2001, A&A, 372, 346 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Natoli, P., Marinucci, D., Cabella, P., de Gasperis, G., & Vittorio, N. 2002, A&A, 383, 1100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335 [NASA ADS] [CrossRef] [Google Scholar]
Planck Collaboration 2005, ESA-SCI, 1, preprint [arXiv:astro-ph/0604069] [Google Scholar]
Poutanen, T., de Gasperis, G., Hivon, E., et al. 2006, A&A, 449, 1311 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Press, W. H. 1992, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed. (Cambridge University Press) [Google Scholar]
Stompor, R., & White, M. 2004, A&A, 419, 783 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
Stompor, R., Balbi, A., Borrill, J., et al. 2002, Phys. Rev. D, 65, 022003 [Google Scholar]
Tegmark, M. 1997, Phys. Rev. D, 55, 5895 [Google Scholar]
Tegmark, M., & de Oliveira-Costa, A. 2001, Phys. Rev. D, 64, 063001 [NASA ADS] [CrossRef] [Google Scholar]
Woodbury, M. A. 1950, Statistical Research Group, Memo. Rep., 42, 4 [Google Scholar]
Wright, E. L., Bennett, C. L., Górski, K. M., Hinshaw, G., & Smoot, G. F. 1996, ApJ, 464, L21 [NASA ADS] [CrossRef] [Google Scholar]

Appendix A: Additional material

A.1. Noise biases

For completeness, we present in Fig. A.1 the noise bias computed from the MADmap NCM and the 25 corresponding noise maps and the Madam bias for f_knee = 10 mHz and 60 s baseline.

Fig. A.1

Analytical and mean Monte Carlo noise biases. Gray band is the 1-σ region for the average, computed by dividing the sample variance by $\sqrt{25}$ . Like the TE, TB and EB biases are both consistent with zero and are not shown here. The analytical bias corresponds to the unsmoothed case presented in Sect. 6.3. Top row: MADmap runs with f_knee = 50 mHz. Bottom row: Madam runs with f_knee = 10mHz.

Fig. A.2

Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot contains the same curves as Fig. A.1.

Fig. A.3

Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot contains the same curves as the top row of Fig. 10.

Fig. A.4

Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot contains the same curves as the bottom row of Fig. 10.

To highlight differences between the estimates and simulated noise maps we also show the fractional differences in Figs. A.2–A.4. These plots complete the ones presented in Sect. 6.2.2.

A.2. Power spectra

Here we present another successful test of the covariance matrix used in power spectrum estimation, Fig. A.5. We also show how the power spectrum estimates can be used to pick out inaccurate residual noise covariances in Fig. A.6.

Fig. A.5

Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 50mHz knee frequency.

Fig. A.6

Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 50mHz knee frequency.

A.3. Singularities

As we have pointed out in Sect. 3.1, the inversion of the inverse noise covariance matrix, N^-1, Eqs. (18) and (9), often needs to be regularized due to the presence of singular or numerically singular modes. In this section we discuss the origin of such modes.

We first note that in all cases considered here the inverse covariance can be expressed as $N -1 = A T ℳ^{-1} A,$ $\appendix \setcounter{section}{1} \begin{equation} \vec N^{-1} = \vec A^\mathrm T \, \mathcal M^{-1} \, \vec A, \end{equation}$ (A.1)where ℳ^-1 is defined to be, $ℳ^{-1} \equiv {\begin{matrix} 𝒩^{-1}, fortheoptimalmaps; \\ 𝒩_{u}^{-1} - 𝒩_{u}^{-1} B \\ \times {(𝒫^{-1} + B T 𝒩_{u}^{-1} B)}^{-1} B T 𝒩_{u}^{-1}, forthedestripedmaps . \end{matrix}$ $\appendix \setcounter{section}{1} \begin{equation} \mathcal M^{-1} \equiv \left\{ \begin{array}{l l} \mathcal N^{-1},\quad\quad \; \qquad \qquad \qquad \qquad\hbox{\rm for the optimal maps;}\\ \mathcal N_\mathrm u^{-1} - \mathcal N_\mathrm u^{-1}\vec B\\ \times \left( \mathcal P^{-1} + \vec B^\mathrm T\mathcal N_\mathrm u^{-1}\vec B \right)^{-1}\vec B^\mathrm T\mathcal N_\mathrm u^{-1}, \hbox{\rm for the destriped maps.}\\ \end{array} \right. \end{equation}$ (A.2)}We assume that the pointing matrix, A, has full column rank, and thus Ax = 0 ⇒ x = 0. This is equivalent to an assumption that the sky signal can indeed be estimated from a given data set. Though this may not be always the case, in particular for the polarization sensitive experiments, it can usually be achieved if some of the ill-constrained pixels are removed from consideration. Given this assumption, the problem of the singular modes of N^-1 becomes that of the matrix, ℳ^-1, defined above.

Let us consider the optimal map case first. The matrix ℳ^-1 is equal to the inverse of the time-domain noise covariance, 𝒩^-1. The latter, Sect. 2.2, is a block Toeplitz matrix with each block defined by an inverse of the noise power spectrum, Eq. (19). Each of those blocks describes the noise properties of one of the stationary data segments assumed in the simulations. For each block the eigenmodes corresponding to the lowest frequencies as permitted by the length of the segment have eigenvalues vastly smaller than the high frequency modes. These modes can therefore lead to near singularities. This is specifically true for zero-frequency modes corresponding to an offset of each of the stationary data segments. The (near) null space of the full matrix will be therefore spanned by all such vectors corresponding to each of the segments.

Due to projection effects, not all of those modes result in singular modes of the final pixel domain noise covariance. However, if for a mode, t, from the null space of 𝒩^-1, there exists a pixel domain vector, x, such as t = Ax, the inverse noise covariance will be singular with eigenvector equal to x.

In the studied case, the scanning strategy is such that the sky areas observed in each of the stationary periods overlap, which efficiently removes most of the potential degenerate vectors. In fact only a single pixel-domain IQU vector of which the I part is one and all the others zero, called hereafter a global offset, can potentially be singular. We will indeed confirm these expectations via numerical results later¹⁰.

In the case of the maps produced with the destriper codes in the absence of any priors, i.e., 𝒫^-1 = 0 the ℳ^-1 matrix has as many singular vectors as the baseline offsets defined by the columns of the offset “pointing” matrix, B. However, as long as all of the offsets cross on the sky the only singular vector of the pixel-domain covariance will again correspond to the global offset vector as in the optimal map case. We note however that unlike in that case, this time this vector is exactly singular. If a prior is employed, as is the case in both the classical and generalized implementations of the destriper technique discussed here, the columns of the matrix B are no longer singular vectors of the matrix ℳ^-1, nor is the global offset vector a singular vector of N. Nevertheless, at least for some common choices of the prior the global offset vector remains nearly singular.

All Tables

Table 1

Pixel side to baseline length at 1 rpm spin rate.

In the text

Table 2

Frequently used symbols in this paper.

In the text

Table 3

Comparison of TT power (rms) in the stripe maps at N_side = 8.

In the text

Table 4

CPU time costs for 1270 GHz detector years = 3 × 10¹⁰ samples.

In the text

All Figures

	Fig. 1 Noise model (dashed lines) and the spectra of simulated noise (solid lines). The two sets of curves correspond to the two considered knee frequencies with f_min = 1.15 × 10^-5 and α = 1.7.
In the text

	Fig. 2 Eigenspectra of the inverse covariance matrices N^-1. MADping, ROMA and Madam results for f_knee = 50mHz overlap completely. Springtide results are for 10mHz.
In the text

Fig. 3

Top: MADping pixel variances for temperature and polarization and the reciprocal condition number of the pixel observation matrices. Bottom: correlation coefficient × 10³ between I − Q, I − U and Q − U pixels. This part of the noise covariance is dominated by white noise which is modeled equivalently in all three paradigms. Hence, MADmap, ROMA, Madam and Springtide results are nearly identical. Maps are rotated into galactic coordinates to show the structure near the ecliptic poles.

In the text

	Fig. 4 A single column of the MADping covariance matrix corresponding to a pixel at the ecliptic equator. For both ROMA and Madam 1.25s baseline counterparts all visible characteristics remain unchanged. We plot the value of the correlation coefficient, R, multiplied by 10³. In order to enhance the features, we have halved the range of the color scale.
In the text

	Fig. 5 A single column of the Springtide covariance matrix corresponding to a pixel near the ecliptic equator. For description of the normalization, see text.
In the text

	Fig. 6 A single column of the MADping covariance matrix corresponding to a pixel at the ecliptic pole. For both ROMA and Madam 1.25s baseline counterparts all visible characteristics remain unchanged. We plot the value of the correlation coefficient, R, multiplied by 10³. In order to enhance the features, we have halved the range of the color scale.
In the text

	Fig. 7 A single column of the Springtide covariance matrix corresponding to a pixel near the ecliptic north pole. For description of the normalization, see text.
In the text

Fig. 8

Empirical χ² distribution functions from the 25 residual noise maps compared with the theoretical cumulative probability density. The black stair line is the empirical distribution function, the blue solid line is the theoretical χ² distribution for 3N_pix − 1 = 36863 degrees of freedom. The red dashed line is the least squares fit of the χ² distribution to the experimental distribution (d.o.f. being the fitting parameter). The horizontal axis is translated to the expected center of the distribution, ⟨ χ² ⟩ = d.o.f. = 36,863, and scaled by the expected deviation, $\sigma_{\chi^2}=\sqrt{2\mathrm{d.o.f.}}$

In the text

Fig. 9

Madam empirical χ² distribution functions from 60 residual noise maps compared with the theoretical cumulative probability density. For a smoothed map, we count the degrees of freedom as the number of included eigenmodes in the inverse NCM. Top: three sets of low-resolutions maps using just one of the low-resolution methods at a time. The high-resolution maps for INW and smoothing methods had N_side = 1024. Bottom: since it is suboptimal to compute a low-resolution spherical harmonic expansion from a noisy high-resolution map, we test how well the smoothing approach works in conjunction with the two pixel-based downgrading methods.

In the text

	Fig. 10 Analytical and mean Monte Carlo noise biases. Gray band is the 1-σ region for the average, computed by dividing the sample variance by $\sqrt{25}$ . Top row: Madam runs at f_knee = 50 mHz using a short 1.25s baseline. Middle row: Madam runs at f_knee = 50 mHz using a relatively long 60s baseline. Bottom row: springtide runs at f_knee = 10 mHz.
In the text

	Fig. 11 Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 10mHz knee frequency.
In the text

	Fig. 12 Examples of signal striping. We show the difference between a binned and destriped signal-only maps. Rows correspond to direct Madam results, 1.25s and 60s, and inverse noise weighted 1.25s baseline case respectively.
In the text

	Fig. 13 Comparison of CMB and stripe power spectra reveals that the striping can significantly bias the EE and BB power spectrum estimates. Stripe map spectra are computed from maps shown in Fig. 12.
In the text

	Fig. 14 Band power estimates of CMB with CMB and foreground stripes.
In the text

Fig. 15

Bandwidth limiting the signal using various window functions. The three panels show the same curves first as raw estimates, then after deconvolving the smoothing and pixel windows and finally after subtracting and dividing by input model. The apodized step window functions (“cosine”) correspond to choices of the thresholds (ℓ₁,ℓ₂) as (20,24), (16,24) and (16,20) respectively. Here, solid lines are for Gaussian windows and dashed lines for the apodized step functions.

In the text

	Fig. 16 TT noise bias computed from smoothed covariance matrices. Solid lines correspond to the Gaussian window functions and the dashed ones to apodized step functions. Left: linear vertical scale. Right: logarithmic vertical scale.
In the text

Fig. A.1

Analytical and mean Monte Carlo noise biases. Gray band is the 1-σ region for the average, computed by dividing the sample variance by $\sqrt{25}$ . Like the TE, TB and EB biases are both consistent with zero and are not shown here. The analytical bias corresponds to the unsmoothed case presented in Sect. 6.3. Top row: MADmap runs with f_knee = 50 mHz. Bottom row: Madam runs with f_knee = 10mHz.

In the text

	Fig. A.2 Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot contains the same curves as Fig. A.1.
In the text

	Fig. A.3 Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot contains the same curves as the top row of Fig. 10.
In the text

	Fig. A.4 Averaged noise biases after subtracting the analytical estimate and normalizing with the standard deviation. This plot contains the same curves as the bottom row of Fig. 10.
In the text

	Fig. A.5 Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 50mHz knee frequency.
In the text

	Fig. A.6 Averaged power spectrum estimates over 25 noise and CMB realizations. The noise has a 50mHz knee frequency.
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.

[1] Armitage, C., & Wandelt, B. D. 2004, Phys. Rev. D, 70, 123007 [NASA ADS] [CrossRef] [Google Scholar]

[2] Armitage-Caplan, C., & Wandelt, B. D. 2009, ApJS, 181, 533 [NASA ADS] [CrossRef] [Google Scholar]

[3] Ashdown, M. A. J., Baccigalupi, C., Balbi, A., et al. 2007a, A&A, 471, 361 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[4] Ashdown, M. A. J., Baccigalupi, C., Balbi, A., et al. 2007b, A&A, 467, 761 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[5] Ashdown, M. A. J., Baccigalupi, C., Bartlett, J. G., et al. 2009, A&A, 493, 753 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[6] Benabed, K., Cardoso, J. F., Prunet, S., & Hivon, E. 2009, MNRAS, 400, 219 [NASA ADS] [CrossRef] [Google Scholar]

[7] Bond, J. R., Jaffe, A. H., & Knox, L. 1998, Phys. Rev. D, 57, 2117 [NASA ADS] [CrossRef] [Google Scholar]

[8] Bond, J. R., Jaffe, A. H., & Knox, L. E. 2000, ApJ, 533, 19 [NASA ADS] [CrossRef] [Google Scholar]

[9] Borrill, J. 1999, in Proceedings of the 5th European SGI/Cray MPP Workshop, Bologna, Italy [Google Scholar]

[10] Cantalupo, C. M., Borrill, J. D., Jaffe, A. H., Kisner, T. S., & Stompor, R. 2010, ApJS, 187, 212 [NASA ADS] [CrossRef] [Google Scholar]

[11] Challinor, A., Fosalba, P., Mortlock, D., et al. 2000, Phys. Rev. D, 62, 123002 [NASA ADS] [CrossRef] [Google Scholar]

[12] Challinor, A. D., Mortlock, D. J., van Leeuwen, F., et al. 2002, MNRAS, 331, 994 [NASA ADS] [CrossRef] [Google Scholar]

[13] de Bernardis, P., Ade, P. A. R., Bock, J. J., et al. 2000, Nature, 404, 955 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]

[14] de Gasperis, G., Balbi, A., Cabella, P., Natoli, P., & Vittorio, N. 2005, A&A, 436, 1159 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[15] Delabrouille, J. 1998, A&AS, 127, 555 [Google Scholar]

[16] Dupac, X., & Tauber, J. 2005, A&A, 430, 363 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[17] Efstathiou, G. 2005, MNRAS, 356, 1549 [NASA ADS] [CrossRef] [Google Scholar]

[18] Efstathiou, G. 2007, MNRAS, 380, 1621 [NASA ADS] [CrossRef] [Google Scholar]

[19] Efstathiou, G., Gratton, S., & Paci, F. 2009, MNRAS, 397, 1355 [NASA ADS] [CrossRef] [Google Scholar]

[20] Górski, K. M., Banday, A. J., Bennett, C. L., et al. 1996, ApJ, 464, L11 [NASA ADS] [CrossRef] [Google Scholar]

[21] Górski, K. M., Hivon, E., Banday, A. J., et al. 2005, ApJ, 622, 759 [NASA ADS] [CrossRef] [Google Scholar]

[22] Gruppuso, A., de Rosa, A., Cabella, P., et al. 2009, MNRAS, 400, 463 [NASA ADS] [CrossRef] [Google Scholar]

[23] Hanany, S., Ade, P., Balbi, A., et al. 2000, ApJ, 545, L5 [NASA ADS] [CrossRef] [Google Scholar]

[24] Hinshaw, G., Nolta, M. R., Bennett, C. L., et al. 2007, ApJS, 170, 288 [NASA ADS] [CrossRef] [Google Scholar]

[25] Jarosik, N., Barnes, C., Bennett, C. L., et al. 2003, ApJS, 148, 29 [NASA ADS] [CrossRef] [Google Scholar]

[26] Jarosik, N., Barnes, C., Greason, M. R., et al. 2007, ApJS, 170, 263 [NASA ADS] [CrossRef] [Google Scholar]

[27] Keihänen, E., Kurki-Suonio, H., Poutanen, T., Maino, D., & Burigana, C. 2004, A&A, 428, 287 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[28] Keihänen, E., Kurki-Suonio, H., & Poutanen, T. 2005, MNRAS, 360, 390 [NASA ADS] [CrossRef] [Google Scholar]

[29] Keihänen, E., Keskitalo, R., Kurki-Suonio, H., Poutanen, T., & Sirviö, A.-S. 2010, A&A, 510, A57 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[30] Kuo, C. L., Ade, P. A. R., Bock, J. J., et al. 2004, ApJ, 600, 32 [NASA ADS] [CrossRef] [Google Scholar]

[31] Kurki-Suonio, H., Keihänen, E., Keskitalo, R., et al. 2009, A&A, 506, 1511 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[32] Maino, D., Burigana, C., Górski, K. M., Mandolesi, N., & Bersanelli, M. 2002, A&A, 387, 356 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[33] Masi, S., Ade, P. A. R., Bock, J. J., et al. 2006, A&A, 458, 687 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[34] Natoli, P., de Gasperis, G., Gheller, C., & Vittorio, N. 2001, A&A, 372, 346 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[35] Natoli, P., Marinucci, D., Cabella, P., de Gasperis, G., & Vittorio, N. 2002, A&A, 383, 1100 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[36] Page, L., Hinshaw, G., Komatsu, E., et al. 2007, ApJS, 170, 335 [NASA ADS] [CrossRef] [Google Scholar]

[37] Planck Collaboration 2005, ESA-SCI, 1, preprint [arXiv:astro-ph/0604069] [Google Scholar]

[38] Poutanen, T., de Gasperis, G., Hivon, E., et al. 2006, A&A, 449, 1311 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[39] Press, W. H. 1992, Numerical Recipes in FORTRAN: The Art of Scientific Computing, 2nd ed. (Cambridge University Press) [Google Scholar]

[40] Stompor, R., & White, M. 2004, A&A, 419, 783 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]

[41] Stompor, R., Balbi, A., Borrill, J., et al. 2002, Phys. Rev. D, 65, 022003 [Google Scholar]

[42] Tegmark, M. 1997, Phys. Rev. D, 55, 5895 [Google Scholar]

[43] Tegmark, M., & de Oliveira-Costa, A. 2001, Phys. Rev. D, 64, 063001 [NASA ADS] [CrossRef] [Google Scholar]

[44] Woodbury, M. A. 1950, Statistical Research Group, Memo. Rep., 42, 4 [Google Scholar]

[45] Wright, E. L., Bennett, C. L., Górski, K. M., Hinshaw, G., & Smoot, G. F. 1996, ApJ, 464, L21 [NASA ADS] [CrossRef] [Google Scholar]

Residual noise covariance for Planck low-resolution data analysis

1. Introduction

2. Algebraic background

2.1. Maps and their covariances

2.2. Time domain noise

2.3. Low-resolution maps

2.3.1. Direct method

2.3.2. Inverse noise weighting (INW)

2.3.3. Harmonic smoothing

3. Numerical calculations of residual noise covariance

3.1. Optimal map covariance

3.2. Destriped map covariance

3.2.1. Conventional destriping

3.2.2. Generalized destriping

3.3. Smoothed covariance matrices

4. Numerical tests and comparison metric

4.1. Quadratic maximum likelihood power spectrum estimation

4.2. Noise bias

5. Simulation

5.1. Scanning strategy

5.2. Planck detectors

5.3. Time-ordered data

5.4. Input maps

6. Results

6.1. Noise covariance matrices

6.2. Noise covariance validation

6.2.1. By χ2

6.2.2. By noise bias

6.2.3. By power spectrum estimation

6.3. Low-resolution maps

6.4. Resource requirements

7. Conclusion

Acknowledgments

References

Appendix A: Additional material

A.1. Noise biases

A.2. Power spectra

A.3. Singularities

All Tables

All Figures

6.2.1. By χ²