© ESO, 2016

1. Introduction

Several methods have been proposed for deconvolution mapmaking for cosmic microwave background (CMB) measurements (Armitage & Wandelt 2004; Armitage-Caplan & Wandelt 2009; Harrison et al. 2011; Keihänen & Reinecke 2012). The aim of deconvolution mapmaking is to produce a map that is free from effects of beam asymmetry. For instance, the image of a point source in a deconvolved map appears symmetric, instead of being elongated according to the beam shape. More importantly from the point of view of cosmology, deconvolution mapmaking produces a map with a unique effective beam window that takes the same form in all map pixels. This happens at the cost of more complicated noise structure. A simple binning operation keeps white noise uncorrelated: white noise in the input time-ordered information (TOI) translates into uncorrelated noise in the binned map. This is not the case for deconvolved maps. The deconvolution operation creates correlations between pixels. Correlated noise present already in the input TOI complicates the noise structure further.

In this work we study the impact of beam deconvolution on noise properties in CMB maps. We build our treatment around one particular deconvolution code, artDeco. The artDeco code (Keihänen & Reinecke 2012) was developed for beam deconvolution of absolute CMB measurements. The code takes as input the time-ordered data stream along with pointing information and a beam description, and produces as primary output the harmonic coefficients a_Xℓm (X = T,E,B) of the sky. From these coefficients one can further construct a sky map that is free from beam asymmetry effects. The computation is carried out in harmonic space. In doing so the method breaks the assumption shared by pixel-based mapmaking methods that the sky signal is constant within a pixel.

We investigate two kinds of input noise: pure white noise and destriped 1 /f noise. The starting point in the latter is a situation where TOI is first destriped with a general destriper with noise prior. One source of such data is the Madam map-making code (Keihänen et al. 2005, 2010), which we use in this work. The code subtracts an estimate of the correlated noise component from the TOI and yields a cleaned TOI stream where residual noise is dominated by white noise. This cleaned TOI is then provided as input to the artDeco beam deconvolver.

The noise in the output map is correlated for two reasons:

• There is residual correlated noise in the TOI stream that cannot be entirely removed by destriping techniques.

• The deconvolution process itself changes the noise properties.

We aim at deriving a noise covariance matrix (NCVM) that captures both phenomena.

We validate our results using simulated Planck LFI data. However, we aim at keeping the method more general. The required conditions are as follows:

• The experiment records the absolute signal, not differential.

• The experiment offers full or nearly full sky coverage.

• Beam shapes are known.

• Sky sampling is dense compared to beam size.

• Noise is piecewise stationary and its spectrum is known.

• The data is cleaned of correlated noise through destriping or another method with equivalent results.

We do not require a particular scanning pattern, or a particular noise spectrum.

Our calculations are carried out in harmonic space, which is the natural domain for beam deconvolution. We construct a harmonic noise covariance matrix that describes the noise properties of the harmonic a_Xℓm coefficients, which are the primary output of artDeco. Because the problem is huge, it is not possible to construct the covariance matrix for all multipoles. In practice, the noise covariance matrix approach is limited to multipoles .

We further produce a pixel noise covariance matrix, which describes the noise properties of deconvolved sky maps, constructed from the a_Xlm coefficients through spherical harmonic transform. A noise covariance matrix is an essential ingredient in many CMB power spectrum estimation methods. A method for producing a noise covariance matrix for un-deconvolved maps has been presented in Keskitalo et al. (2010). Its application to Planck LFI data is presented in Planck Collaboration VI (2016).

The paper is structured as follows. In Sect. 2 we perform a series of full-resolution simulations, to gain a general view of the impact of deconvolution on residual noise. In Sect. 3 we focus on white noise. We derive a white noise covariance matrix, and compare its predictions with Monte Carlo simulations. The main result of this paper is the derivation of a full low-resolution noise covariance matrix, presented in Sect. 4. In Sect. 5 we perform a thorough validation of the harmonic covariance matrix through Monte Carlo simulations. We show results from χ² tests and from noise bias comparison. In Sect. 6 we perform analysis in pixel space. We summarise our results in Sect. 7.

2. Simulations

We start by performing a series of full-resolution noise simulations to assess the impact of deconvolution. Our simulations mimic the data from the LFI instrument of Planck experiment. We created a time-ordered data stream, which we filled with simulated noise, and fed it as input to the artDeco deconvolver. The simulations were run on Sisu, the Cray XC40 supercomputer of CSC, Finland.

2.1. Generating noise

We regenerated the LFI detector pointing with LevelS software (Reinecke et al. 2006). The simulated mission covers 4 yr of data. We used the internal noise generator of the Madam map-making code to generate a noise time stream. The noise properties mimic those of the real Planck data. We considered the three LFI channels: 30 GHz, 44 GHz, and 70 GHz. The noise parameters used in simulations were taken from Planck Collaboration II (2016). For convenience, the parameters are listed in Table 1.

We consider two types of noise. In the first case we generated pure white noise, with variances taken from Table 1. The variance is assumed constant in time, but varies between radiometers. In the second series of simulations we generated realistic 1 /f noise, again with parameters from Table 1. We destriped the noise stream with Madam, and stored the cleaned data on disk. The data was then fed as input to artDeco.

The destriping procedure follows the one applied to LFI mapmaking (Planck Collaboration VI 2016). We applied the same destriping mask, and used the same HEALPix map resolution N_side = 1024 (3.44′), and same baseline length (1 s for 70 GHz and 0.25 s for 30 GHz) as in real Planck mapmaking. We refer to these simulations as realistic noise simulations, in contrast to the white noise simulations. In both simulations we applied realistic Planck flags to the TOI stream. This excludes roughly 8% of the data. Madam has the capacity of storing the destriped TOI on disk in “4D map” format. The four dimensions in this context refer to the three angles θ,φ,ψ which define the detector pointing, and pointing period index as fourth dimension. Pointing period refers to a period of typical duration of 40 min, where the Planck scanning pattern follows a fixed ring on the sky.

Madam generates a three-dimensional grid by discretising the detector pointing angles. The pixelization obeys HEALPix pixelization in θ and φ, and the ψ angle is divided uniformly into N_psi bins. Another parameter N_side controls the resolution in θ,φ. The TOI is binned on the grid according to its pointing, and the grid is stored by pointing period. The data volume is compressed typically by a factor of 20 compared to the full TOI. Radiometer-specific parameters required in subsequent analysis steps are stored as keywords in the 4D object header. For the present study we used resolution N_side = 1024, N_psi = 4096. At this resolution the 4D maps take roughly 8 GB of disk space per detector.

In white noise simulations, we use Madam only to generate the noise, to apply the flags, and to store the data stream in 4D map format. No destriping is applied.

The 4D map is an intermediate data object that is useful for many purposes. In particular, it can further be compressed into the “3D map” format that serves as input to artDeco. This step is trivial, it consists of coadding all pointing periods on a single grid. Further, we can construct a binned map from the same 4D map input. The binning operation represents “normal” mapmaking without beam deconvolution. The combined operation of destriping, 4D map production, and map binning, is nearly equivalent to the operation of producing a destriped map directly by Madam. The only difference is the discretization of the ψ angle, which at the chosen resolution N_psi = 4096 has a negligible effect on the final map.

To reduce the scatter in results, we produced several noise realizations for each case studied. At 30 GHz, we generated 40 realizations for each case. At 44 GHz and 70 GHz, where the simulations are more demanding, we produced 10 realizations. The spectra shown are averaged over available realizations. A complete list of simulations is given in Table 5.

2.2. Destriping options

The full simulated TOI consists of four years of data for all LFI radiometers. In cases where we deconvolve only a subset of data from a frequency channel, we have two options for the destriping step. Either we can use the full data set for destriping, and then pass a subset of the cleaned TOI to deconvolution, or we can use the same data subset in both processing steps. We refer to the two options as full destriping and independent destriping, respectively. Full destriping leads to a smaller residual noise level in the final data products, as there is more information available for the destriping solution and removal of correlated noise. On the other hand, independent destriping has the benefit over full destriping that the residual noise in the final products is uncorrelated from one data set to another.

We compared the two destriping options for horn pair 18/23. This represents one third of the full 70 GHz data set. The horn pair data set was selected for closer examination on the basis of low-resolution simulations (Sect. 5), where it showed as the “worst case” subset.

2.3. Binned maps

We compare the deconvolved spectrum to the binned map spectrum, which we constructed as follows. We bin the same destriped TOI into a sky map as ${m}\mathrm{=}\mathrm{(}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}{\mathrm{)}}^{-1}\mathrm{P}{\mathrm{C}}_{\mathrm{n}}^{-1}\hat{{y}}\mathit{.}$(1)Here P is the pointing matrix, $\hat{{y}}$ is the input time-ordered data (TOI), and C_n is a diagonal white noise covariance. If several detectors are involved, their data are appended to a single TOI stream.

Binning represents the usual mapmaking process, which does not attempt to correct for the beam shapes. We then compute the pseudo spectrum of the map with the Anafast tool of the HEALPix package. We refer to the spectrum obtained this way as map spectrum. A comparison of the map spectrum and deconvolved spectrum gives insight to the impact of the deconvolution process on noise. We perform the map-binning operation starting from the same 4D map objects that serve as input to the deconvolution procedure. That way we can be sure that we have exactly the same inputs in both procedures.

We deviate from the mapmaking procedure applied in actual Planck analysis in one aspect. LFI mapmaking applies the horn-uniform radiometer weighting scheme, when combining several radiometers into one map (Planck Collaboration VI 2016). Horn-uniform weighting reduces leakage of temperature signal to polarization through beam shape mismatch, at the cost of slightly increased noise compared to the other option, noise weighting. In noise weighting, the radiometer weights are constructed as ${\mathit{w}}_{\mathit{j}}\mathrm{=}\mathrm{1}\mathit{/}{\mathit{\sigma }}_{\mathit{j}}^{\mathrm{2}}$ where σ_j is the white noise rms for radiometer j, from Table 1. In horn-uniform weighting the weights are made equal for the M and S radiometers of same horn. Because the beam shapes are explicitly accounted for in deconvolution, horn-uniform weighting provides no benefit over noise weighting, and artDeco internally applies noise weighting. Since our aim is to study the effects of beam deconvolution, and not those of different radiometer weighting schemes, we build also the binned maps with noise weighting. Horn-uniform weighting is, however, applied in the destriping phase. This way we have two mapmaking methods that take as input the same exact destriped TOI, and differ only in the way they deal (or do not) with the beam. The detector weighting scheme has a small effect on the residual white noise level. We return to this in Sect. 3.3.

2.4. Deconvolution

ArtDeco yields as output the harmonic coefficients a_Tℓm, a_Eℓm, a_Bℓm. with ℓ = 0...ℓ_max. Here ℓ_max is an input parameter that defines the multipole range considered. Another input parameter k_max sets the level of beam asymmetry taken into account. The harmonic beam expansion is constructed for b_sℓk, where 0 <ℓ ≤ ℓ_max and | k | ≤ k_max. The maximum ℓ_max that can be reached is dependent on the beam width. For the full-resolution simulations we chose ℓ_max = 700 for 30 GHz, ℓ_max = 1000 for 44 GHz, and ℓ_max = 1500 for 70 GHz. In all cases we used k_max = 6.

The artDeco deconvolution operation can formally be written as ${a}\mathrm{=}\mathrm{(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}{\mathrm{)}}^{-1}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\hat{{y}}\mathit{.}$(2)The solution is built on the assumption that the noise in the input TOI is white, or at least dominated by white noise. Matrix A is given by (Keihänen & Reinecke 2012), ${\mathit{A}}_{\mathit{s\ell m}}^{\mathit{j}}\mathrm{=}\sum _{\mathit{k}}{\mathit{b}}_{\mathit{s\ell k}}^{\mathrm{\ast }}{\mathit{D}}_{\mathit{mk}}^{\mathit{\ell }\mathrm{\ast }}\mathrm{\left(}{\mathit{\omega }}_{\mathit{j}}\mathrm{\right)}\mathit{,}$(3)where b is the harmonic representation of the beam, D is a Wigner function, which takes as argument the pointing angles θ, φ, and ψ, and ω_j is a short for the triplet of pointing angles for sample j.

The spin index s takes values 0 (temperature) and ± 2 (polarization). Index k takes values from 0 to ± k_max. Indices ℓ,m are the usual indices of spherical harmonics, ℓ being related to angular scale and m to orientation. In this notation, the output harmonic coefficients are written as a_sℓm and the input TOI stream as y_j. The E,B coefficients are further constructed as $\begin{array}{ccc}{\mathit{a}}_{\mathit{E\ell m}}& \mathrm{=}& \\ {\mathit{a}}_{\mathit{B\ell m}}& \mathrm{=}& \mathrm{-}\frac{\mathrm{1}}{\mathrm{2}\mathit{i}}\mathrm{(}{\mathit{a}}_{\mathrm{2}\mathit{\ell m}}\mathrm{-}{\mathit{a}}_{\mathrm{-}\mathrm{2}\mathit{\ell m}}\mathrm{)}\mathit{.}\end{array}$(4)In this work we are interested in the impact of deconvolution on noise. Since both destriping and deconvolution are linear processes, we can examine the signal and noise components in isolation. From now on we assume that the input TOI consists of pure noise (white and correlated). We compute the deconvolved spectrum from the output a_Xℓm coefficients as $\mathit{Ĉ}\mathit{XY\ell }\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}\mathrm{=}\mathrm{-}\mathit{\ell }}^{\mathit{\ell }}{\mathit{a}}_{\mathit{X\ell m}}^{\mathrm{\ast }}{\mathit{a}}_{\mathit{Y\ell m}}$(5)where X,Y stand for T,E,B. We further average the spectrum over the available noise realizations.

2.5. Beams

Throughout this paper we use Planck LFI beams (Planck Collaboration IV 2016). We included in the deconvolution process the main and intermediate beam components. The contribution of far sidelobes is assumed to be removed before the destriping step (Planck Collaboration II 2016).

Table 2 lists some characteristics of the beams per frequency. The FWHM and ellipticity values are based on effective FEBeCoP beams (Mitra et al. 2011), and include the main beam only. The values are taken from Planck Collaboration IV (2016), and are shown here to give an idea of typical beam shape at each frequency. Ellipticity is defined as the ratio of beam widths along the two main axes of the beam. A symmetric beam will have ellipticity equal to 1.

For simulations we used the combined main and intermediate radiometer scanning beams. The characteristics of these can be found in Planck Collaboration IV (2016). The last column of Table 2 gives the beam efficiency at each frequency. These are computed from of the harmonic beam expansion of the noise-weighted frequency beam. Beam efficiency is proportional to the monopole component of the expansion, and normalised so that the full 4π beam has efficiency 1. The squared value gives the effect in the spectral domain. The missing fraction reflects the power absorbed by far sidelobes.

For comparison purposes we constructed a symmetrized beam window as follows. First we computed the weighted sum individual radiometer beams, $\mathit{b̅}\mathit{\ell m}\mathrm{=}{\left[\sum _{\mathit{j}}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{j}}^{\mathrm{2}}}\right]}^{-1}\sum _{\mathit{j}}\frac{\mathrm{1}}{{\mathit{\sigma }}_{\mathit{j}}^{\mathrm{2}}}{\mathit{b}}_{\mathit{\ell m}}^{\mathit{j}}$(6)where j labels the radiometers, and σ_j are their respective white noise rms. The m = 0 component of the beam expansion, ${\mathit{W}}_{\mathit{\ell }\mathrm{[}\mathit{m}\mathrm{=}\mathrm{0}\mathrm{]}}\mathrm{=}\sqrt{\frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}}\mathit{b̅}\mathit{\ell }\mathrm{0}$(7)represents a beam that is symmetrized by averaging over the azimuthal angle around the beam centre. Factor 4π/ (2ℓ + 1) serves to normalize the window function to unity at ℓ = 0 for a beam with efficiency equal to 1. The power spectrum of the beam expansion provides an alternative way of constructing a symmetrized beam, as proposed by Wu et al. (2001), ${\mathit{W}}_{\mathit{\ell }\mathrm{\left[}\mathrm{ms}\mathrm{\right]}}^{\mathrm{2}}\mathrm{=}\frac{\mathrm{4}\mathit{\pi }}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}\sum _{\mathit{m}}\mathrm{|}\mathit{b̅}\mathit{\ell m}{\mathrm{|}}^{\mathrm{2}}\mathit{.}$(8)where “ms” stands for mean of squares. The beam window of Eq. (8) is always larger than that from Eq. (7), thus representing a narrower beam. For LFI beams the difference is small. The squared window functions ${\mathit{W}}_{\mathit{l}}^{\mathrm{2}}$ for both definitions are plotted in Fig. 1. The squared window can directly be compared with power spectra.

2.6. Results from frequency simulations

Figure 2 shows the 30 GHz TT spectrum for white noise and realistic noise inputs. All spectra have been averaged over 40 noise realizations. We show in the same plot the deconvolved spectrum and the binned map spectrum. The main effect from deconvolution is evident. While the white noise map spectrum is uniform, the deconvolved spectrum rises steeply towards higher multipoles. As a first approximation, the effect of deconvolution is to scale the binned map spectrum by the inverse of the symmetric beam window squared (shown in dashed line type). This is, however, only true as a first approximation. The real deconvolution spectrum rises even more steeply due to beam asymmetry. We investigate this further in Sect. 3.

The primary deconvolution output is the harmonic coefficients a_Xlm. For many practical purposes one may want to transform the result into an ordinary sky map. The procedure of constructing a sky map from the deconvolved a_Xlm coefficients involves a smoothing operation that brings the noise spectrum back down to the level of un-deconvolved noise. We discuss this procedure in Sect. 6. In this and the following three sections, where we do analysis in harmonic space, we always show the unsmoothed raw noise. In harmonic space the smoothing is a simple ℓ-dependent scaling operation which does not change the relative differences between spectra. The relative differences between spectra are thus more important than their absolute level.

Residual correlated noise dominates the lower multipoles, but decreases towards higher multipoles. Figure 3 shows the ratio of correlated residual noise spectrum and white noise spectrum. The quantity shown is computed as (C_ℓreal − C_ℓwn) /C_ℓwn, where C_ℓreal and C_ℓwn are the spectra from realistic and white noise simulations, respectively.

An interesting and perhaps surprising observation is that deconvolution leaves relatively less correlated residual noise on top of the white noise component at high multipoles, than does normal mapmaking. This has important consequences. The level of correlated residual noise is a measure of the error we make if we neglect the correlated component when estimating the residual noise. In Sect. 3 we derive a model for the deconvolved white noise spectrum, which we use as an approximation for the full noise at high multipoles. The accuracy of the approximation can be read from Fig. 3. Further, MC simulations involving only white noise spectrum are much cheaper than simulations with realistic noise, as the destriping step is avoided.

We now look closer at the white noise spectrum. The spectra vary strongly due to the limited number of realizations. However, when comparing spectra derived from the same simulation with different methods, we can make use of the fact that we have the same noise realizations in both cases. We take the white noise map spectrum as reference level, and divide the other spectra by it. This reduces the scatter in the spectra and shows the relative differences clearly. Spectral ratios constructed this way are plotted in Fig. 4. The white noise map spectrum is uniform by construction. The shape of the deconvolved white noise spectrum reflects the beam window function. The deconvolved spectrum rises above the binned map spectrum at the lowest multipoles by approximately 2%. This is the effect of the missing power lost into the far sidelobe component of the beam. (The y = 99% efficiency of Table 2 translates to an 1 − y² ≈ 2% effect in spectral domain.) Deconvolution scales the signal up to correct for the missing power, and the same increase is observed in the noise spectrum.

The purple line in Fig. 4 is the result of a simulation where we deconvolved for a perfect delta beam (“deconvolution without deconvolution”). Even if the beam is a delta peak, the deconvolution procedure is not identical to the combined map-binning and anafast procedure. Deconvolution performs a maximum-likelihood fit of the harmonic coefficients to the TOI, while anafast computes a harmonic expansion over the sky map, without taking into account the varying hit counts. The two procedures are identical only in the case where the distribution of observations is uniform over the sphere. The spectrum from deconvolution with a delta beam is slightly below the reference level of the white noise map spectrum, but the difference is small. Above ℓ = 600 the spectrum drops, which is a typical boundary effect in deconvolution. The drop occurs at a much higher multipole than the drop in the spectral ratio of Fig. 3. The exact mechanism behind the latter drop remains unclear.

2.7. Horn pair LFI18/23

Finally we compare the two destriping options discussed in Sect. 2.2. We performed two simulations, where we deconvolved the data from horn pair LFI18/23. In the first set of simulations we used for destriping the same 18/23 horn-pair data set, in the second one the full 70 GHz data set. The results are shown in in Fig. 5. We show the deconvolved spectra for white noise and for realistic noise, with both destriping options. The white noise map spectrum is shown in the same plot for reference. The scatter is large due to the small number of realizations (10), but the main conclusion is still clear. Full destriping leaves considerably less residual noise at low multipoles.

In the lower panel we replot the spectra the same way we did for full 30 GHz data in Fig. 4. We divided the spectra by the white noise map spectrum to reduce scatter, and to bring out the differences more clearly.

3. White noise covariance

3.1. White noise covariance

We proceed to study the properties of deconvolved noise in more detail. The deconvolution operation of Eq. (2), when operating on a noise time stream n, yields a vector of harmonic coefficients whose properties are described by the harmonic noise covariance matrix given by $\mathrm{C}\mathrm{=}\mathrm{(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}{\mathrm{)}}^{-1}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{⟨}{n}{{n}}^{\mathit{T}}\mathrm{⟩}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}\mathrm{(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}{\mathrm{)}}^{-1}\mathrm{.}$(9)It describes the noise correlation between elements a_slm.

In this section we study the simpler case where the input TOI consists of pure white noise, with diagonal variance ${\mathrm{C}}_{\mathrm{n}}\mathrm{=}\mathrm{⟨}{n}{{n}}^{\mathit{T}}\mathrm{⟩}\mathrm{.}$(10)The harmonic noise covariance matrix simplifies into $\mathrm{C}\mathrm{=}\mathrm{(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}{\mathrm{)}}^{-1}\mathrm{.}$(11)The noise covariance is thus given by the inverse of the deconvolution matrix $\mathrm{\left(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}\mathrm{\right)}$. This is a product we can extract from artDeco. The size of the matrix rapidly increases with increasing ℓ_max. The rank of the matrix with polarization is 3(l_max + 1)². The memory requirement increases as square of the rank, and the CPU time needed for inversion in the third power, or ${\mathit{l}}_{\max }^{\mathrm{6}}$. The full matrix inversion rapidly becomes prohibitive. In practice we can only construct it for a limited multipole range.

3.2. Noise bias

We want to compare the noise spectra from simulations to the covariance matrix. We derive from the covariance matrix a noise bias, a prediction for the expectation value of the noise spectrum of Eq. (5).

We define g_Xs as factors that convert the spin harmonics a_sℓm into a_Xℓm, X = T,E,B. We write Eq. (4) as ${\mathit{a}}_{\mathit{X\ell m}}\mathrm{=}\sum _{\mathit{s}}{\mathit{g}}_{\mathit{Xs}}{\mathit{a}}_{\mathit{s\ell m}}$(12)Now

$\begin{array}{ccc}& & {\mathit{g}}_{\mathit{T,}\mathrm{0}}\mathrm{=}\mathrm{1}\\ & & {\mathit{g}}_{\mathit{E,}\mathrm{2}}\mathrm{=}{\mathit{g}}_{\mathit{E,}\mathrm{-}\mathrm{2}}\mathrm{=}\mathrm{-}\mathrm{1}\mathit{/}\mathrm{2}\\ & & {\mathit{g}}_{\mathit{B}\mathrm{2}}\mathrm{=}\mathrm{-}{\mathit{g}}_{\mathit{B,}\mathrm{-}\mathrm{2}}\mathrm{=}\mathit{i}\mathit{/}\mathrm{2}\\ & & {\mathit{g}}_{\mathit{T,}\mathrm{2}}\mathrm{=}{\mathit{g}}_{\mathit{T,}\mathrm{-}\mathrm{2}}\mathrm{=}{\mathit{g}}_{\mathit{E,}\mathrm{0}}\mathrm{=}{\mathit{g}}_{\mathit{B,}\mathrm{0}}\mathrm{=}\mathrm{0.}\end{array}$(13)The noise bias estimate is obtained from the noise covariance C as ${\mathit{C}}_{\mathit{XY\ell }}\mathrm{=}\sum _{\mathit{s}{\mathit{s}}^{\mathrm{\prime }}\mathit{m}{\mathit{m}}^{\mathrm{\prime }}{\mathit{\ell }}^{\mathrm{\prime }}{\mathit{\ell }}^{\mathrm{\prime \prime }}}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}{\mathit{g}}_{\mathit{Xs}}{\mathit{g}}_{\mathit{Y}{\mathit{s}}^{\mathrm{\prime }}}^{\mathrm{\ast }}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime \prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}{\mathrm{C}}_{{\mathit{\ell }}^{\mathrm{\prime }}\mathit{ms,}{\mathit{\ell }}^{\mathrm{\prime \prime }}{\mathit{m}}^{\mathrm{\prime }}{\mathit{s}}^{\mathrm{\prime }}}\mathit{.}$(14)Formally we can write ${\mathit{C}}_{\mathit{XY\ell }}\mathrm{=}\mathrm{Tr}\mathrm{\left[}\mathrm{C}{\mathrm{G}}^{\mathit{XY\ell }}\mathrm{\right]}$(15)where ${\mathrm{G}}_{{\mathit{\ell }}^{\mathrm{\prime }}\mathit{ms,}{\mathit{\ell }}^{\mathrm{\prime \prime }}{\mathit{m}}^{\mathrm{\prime }}{\mathit{s}}^{\mathrm{\prime }}}^{\mathit{XY\ell }}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}\mathit{\ell }\mathrm{+}\mathrm{1}}{\mathit{g}}_{\mathit{Xs}}{{\mathit{g}}^{\mathrm{\ast }}}_{\mathit{Y}{\mathit{s}}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime }}}{\mathit{\delta }}_{\mathit{\ell }{\mathit{\ell }}^{\mathrm{\prime \prime }}}{\mathit{\delta }}_{\mathit{m}{\mathit{m}}^{\mathrm{\prime }}}\mathit{.}$(16)The noise bias is extracted from a block-diagonal part of the NCVM, where each 3 × 3 block corresponds to a pair of (ℓ,m) indeces, and s index takes all allowed values (0, ± 2).

Figure 6 shows the TT, EE, and TE, noise spectra for three LFI frequencies. The BB spectra are quite similar to EE, and are not shown. We show the deconvolved white noise spectrum from MC simulation, along with a noise bias estimate obtained from the white noise covariance. We show also the realistic noise spectra, which rise above the white noise spectra at low multipoles. The white noise component dominates the noise at high multipoles. Consequently, the white noise covariance provides a good model for the full noise at these high multipoles. We show also the raw estimate based on the symmetric beam window. It provides a reasonable approximation at intermediate multipoles, but underestimates the noise at high multipoles.

Because of the large size of the deconvolution matrix, we cannot compute the exact noise bias for the full multipole range. We constructed the 30 GHz bias in pieces. First we constructed and inverted the matrix for range ℓ = 0−110 to obtain the noise bias for multipoles ℓ = 0−100. We dropped the last 10 multipoles to eliminate boundary effects. Similarly, noise bias for multipoles ℓ = 100−200 was constructed from a matrix computed for ℓ = 90−210. Multipoles ℓ = 200−300 were constructed from two pieces in range ℓ = 190−260 and ℓ = 240−310. Because the number of m multipoles increases with increasing ℓ, we had to make the pieces narrower as we proceeded towards higher multipoles. Range ℓ = 300−400 was composed from pieces of 20 multipoles each, range ℓ = 400−600 from pieces of 10, and finally ℓ = 600−700 from pieces of 5 multipoles. In all cases above ℓ = 300 we included 5 extra multiples at both ends, which we then dropped when combining the pieces. The whole computation took 20 000 CPU hours. The estimate constructed this way is far more accurate than the one based on symmetric beam window, as can be seen from Fig. 6. However, the bias still underestimates the MC noise spectrum by 2−3% because it neglects correlations between distant multipoles. At 30 GHz this is not a major problem, because we make an error of similar magnitude anyway by neglecting the correlated noise component (see Fig. 3).

The computation becomes increasingly heavy towards higher multipoles. We therefore adopted a different procedure for 44 GHz and 70 GHz. We constructed the bias in one piece for multipoles ℓ = 0 − 200. At higher multipoles we concentrated the computational resources on a number of individual multipoles, for which we estimated the bias as accurately as possible. We optimized the matrix inversion routine in such a way that it only computes the diagonal elements actually needed for the bias. We constructed the white noise covariance for the multipoles from ℓ₀–Δℓ to ℓ₀ +Δℓ, and picked the central value l₀ for plotting. Values of ℓ₀ were taken to be multiples of 100.

At 44 GHz we used Δℓ = 30 for ℓ₀ = 300−500, and Δℓ = 20 for l₀ = 400−800. A fortunate coincidence comes to help here. The correlations between distant multipoles become weaker with increasing ℓ, so that even though we had to make the pieces narrower, the error in the bias rather stayed constant. The noise bias values obtained for ℓ₀ are shown by circles in Fig. 6. The differences between predicted noise bias and MC spectra are of the order of 1%, and in many cases the predicted value lies within the statistical scatter of the MC spectrum.

At 70 GHz the computation is more demanding. We applied a similar procedure as at 44 GHz. We constructed the full covariance for ℓ = 0−200, then evaluated the bias for multiples of 100 above that. We used Δℓ=30 for l₀ = 300−500, Δℓ = 20 for l₀ = 600−800, and Δℓ = 10 for l₀ = 900−1400. We reach an accuracy of 1−3%.

The importance of parameter Δℓ is illustrated in Fig. 7. We computed the noise bias for ℓ = 500 at 30 GHz, and for ℓ = 800 at 70 GHz, for increasing values of Δℓ. At 30 GHz we show the result for Δℓ = 5, 10, 20, 40. The values approach the Monte Carlo result as Δℓ increases, as expected. The computational burden rapidly increases with increasing Δℓ. While computing the Δℓ = 20 point took 94 min on 960 cores, the Δℓ = 40 point already took three hours on 1920 cores.

At 70 GHz we computed the bias for Δℓ = 5, 10, 20. The last point took 4 hours on 1920 cores. Again the results are approaching the MC result, but the three data points are not sufficient to conclude if there is a real convergence. Also the statistical scatter in the MC spectrum, based on 10 realizations, is too large for us to give an accurate estimate of the remaining error, other than saying that it is of the order of 1%.

The 1−3% error in bias estimation may still be too large for many applications. To compute the bias more accurately, we would need a more efficient matrix inversion routine, or one that constructs the diagonal of the inverse. Although more efficient methods surely exist, they are not likely to bring us very far, since the largest matrices we invert already have rank over 200 000. Should higher accuracy be required, the best approach is probably to produce a larger number of Monte Carlo simulations and to extract the bias from them. This will require significant optimization of the pipeline.

The realistic noise simulations at 70 GHz took 15 min of wall-clock time per MC cycle on 960 cores, corresponding to 250 CPU hours. Roughly half of that went into noise generation and destriping in Madam, another half into the deconvolution process. The cost of white noise simulations was 52 CPUh per realization. At 30 GHz, the cost was 52 CPUh for realistic and 41 CPUh for white noise simulations.

The limiting factor in our current pipeline is the disk space, rather than CPU time, since we are writing the intermediate 4D map products on disk. The pipeline could be made a lot more efficient by if the step of writing disk I/O could be avoided. In case of white noise simulations, this could be done by generating a noise 3D map directly in artDeco memory, the same way we generate a noise TOI in Madam. This remains a topic for future development.

At 44 GHz, white noise simulations give sub-percent accuracy above multipoles ℓ = 800, as can be seen from Fig. 3. At 70 GHz, this accuracy is reached already around ℓ = 400. At 30 GHz, the contribution from the correlated component is above 2% over the whole multipole range.

3.3. Resolution

The auto spectra (TT, EE, BB) of a binned white noise map become uniform, when averaged over a large number of realizations. The level of the spectrum is dependent on the hit count distribution. For a given total number of hits and fixed noise rms per sample, the lowest spectral level is obtained when the hits are uniformly distributed over the sky. Downgrading the resolution of a white noise map in general reduces the spectral level, as the distribution becomes more uniform.

The spectral levels can be predicted from a known hit distribution, as described in Planck Collaboration VI (2016). We tabulate the predicted white noise levels for our simulation parameters for different resolutions in Table 3. We show the results for two radiometer weighting schemes: noise weighting and horn-uniform weighting. Horn-uniform weighting was used in actual LFI mapmaking.

In deconvolution the situation is different. The resolution of the input 4D map does not have a significant effect on residual noise. Instead, resolution dependence enters through parameter ℓ_max. This is demonstrated in Fig. 8. We deconvolved the same 30 GHz data set with ℓ_max = 50, 100, 200, 400. We show also the corresponding noise bias estimates for ℓ_max = 50,100, 200, constructed from the noise covariance matrix. For ℓ_max = 400 we were unable to compute the exact noise bias, due to the huge size of the covariance matrix. For plotting purposes the spectra have been divided by the spectrum of an N_side = 1024 (3.44′) white noise binned map. The noise bias are scaled to the same dimensionless units by dividing them by the theoretical white noise level of 3.2180 × 10^-15 K², taken from Table 3. The noise level at low multipoles increases with increasing ℓ_max, as there are more multipoles to be solved from same amount of data.

The result has relevance for the following section, where we construct the full low-resolution noise covariance matrix for multipole range ℓ = 0−50. In order to have an exact match between the covariance matrix and the actual a_Xℓm coefficients, it is required that both have been computed with same ℓ_max.

4. Full noise covariance

We now consider the more complicated case where the input consists of 1 /f type noise, which has been destriped to remove the correlated component. We aim at deriving a full noise covariance matrix, which takes into account the correlated noise component. We assume a destriping process that involves a noise prior. One such implementation is the Madam code presented in Keihänen et al. (2010). The map binned from the destriped TOI approaches that of generalised least squared (GLS) mapmaking methods (Poutanen et al. 2006; Ashdown et al. 2007a,b, 2009) when the baseline length approaches one sample.

4.1. Destriping

In the destriping approach the initial noise stream is modelled as ${{n}}^{\mathrm{\prime }}\mathrm{=}\mathrm{F}{b}\mathrm{+}{{n}}_{\mathrm{wn}}\mathrm{,}$(17)where n_wn represents white noise, and the correlated component is modelled by a sequence of baselines b. Matrix F formally spreads the noise baselines into TOI.

The destriping solution gives an estimate for the baseline vector b as $\hat{{b}}\mathrm{=}\mathrm{(}{\mathrm{F}}^{\mathit{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{ZF}\mathrm{+}{\mathrm{C}}_{\mathrm{b}}^{-1}{\mathrm{)}}^{-1}{\mathrm{F}}^{\mathit{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{Z}{y}\mathrm{,}$(18)where $\mathrm{Z}\mathrm{=}\mathrm{I}\mathrm{-}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}\mathrm{(}{\mathrm{P}}^{\mathit{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}{\mathrm{)}}^{-1}{\mathrm{P}}^{\mathit{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{.}$(19)Here P is the pointing matrix, and C_b is the a priori covariance of the baselines, constructed from the known noise spectrum. The cleaned TOI is constructed as $\hat{{y}}\mathrm{=}{y}\mathrm{-}\mathrm{F}\hat{{b}}\mathrm{=}{y}\mathrm{-}\mathrm{F}\mathrm{(}{\mathrm{F}}^{\mathit{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{ZF}\mathrm{+}{\mathrm{C}}_{\mathrm{b}}^{-1}{\mathrm{)}}^{-1}{\mathrm{F}}^{\mathit{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{Z}{y}\mathrm{.}$(20)

Consider then the residual noise n in a TOI stream that has been processed according to Eq. (20). We aim at deriving a formula for ⟨ nn^T ⟩, and finally for the covariance of Eq. (9), where n is now the destriped noise stream. We assume that the original noise stream n′ obeys the model (17), so that $\mathrm{⟨}{{n}}^{\mathrm{\prime }}{{n}}^{\mathrm{\prime T}}\mathrm{⟩}\mathrm{=}{\mathrm{C}}_{\mathrm{n}}\mathrm{+}\mathrm{F}{\mathrm{C}}_{\mathrm{b}}{\mathrm{F}}^{\mathrm{T}}\mathrm{.}$(21)The baseline length, hidden inside matrix F, is a key parameter in the calculation. The baseline length plays a double role in the noise covariance. On the one hand, it limits the frequencies that are taken into account by the covariance matrix. The frequencies above the inverse of the baseline length are neglected. If we choose too long a baseline length, we will underestimate the residual noise. For optimal results, we should thus choose a baseline length as short as possible.

On the other hand, the baseline length in F represents the baseline length that was used in destriping. If we select a shorter baseline for covariance computations, we will not model the destriping process correctly. However, simulations with Planck LFI data (Planck Collaboration VI 2016) show that the destriping result converges already around a baseline length of one second. There is correlated noise at higher frequencies, but destriping is unable to remove it, and making the baseline shorter does not change the results. We can turn this around and argue that we do not make a significant error if we assume a shorter baseline length than was actually used for destriping, as long as the actual baseline length was short enough that the optimal solution was already reached.

All this put together, we argue that we get the optimal noise covariance when we make as short a baseline as possible. We bring this to the extreme and set the baseline length to one TOI sample. With this assumption, matrix F reduces into a unity matrix, F = I. The destriped TOI (20) simplifies into $\hat{{y}}\mathrm{=}{y}\mathrm{-}\mathrm{(}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{Z}\mathrm{+}{\mathrm{C}}_{\mathrm{b}}^{-1}{\mathrm{)}}^{-1}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{Z}{y}\mathrm{=}\mathrm{(}{\mathrm{C}}_{\mathrm{b}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{Z}\mathrm{+}\mathrm{1}{\mathrm{)}}^{-1}{y}\mathrm{.}$(22)

4.2. Covariance matrix

We can rewrite Eq. (22) for the destriped noise component ${n}\mathrm{=}\mathrm{(}{\mathrm{C}}_{\mathrm{b}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{Z}\mathrm{+}\mathrm{1}{\mathrm{)}}^{-1}{{n}}^{\mathrm{\prime }}\mathit{.}$(23)We reformulate this with the aim of constructing the time domain covariance ⟨ nn^T ⟩. Inserting Z we obtain ${n}\mathrm{=}{\mathrm{C}}_{\mathrm{n}}\mathrm{[}{\mathrm{C}}_{\mathrm{n}}\mathrm{+}{\mathrm{C}}_{\mathrm{b}}\mathrm{-}{\mathrm{C}}_{\mathrm{b}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}\mathrm{(}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}{\mathrm{)}}^{-1}{\mathrm{P}}^{\mathrm{T}}{\mathrm{]}}^{-1}{{n}}^{\mathrm{\prime }}\mathit{.}$(24)We apply the Woodbury formula to get $\begin{array}{ccc}{n}& \mathrm{=}& {\mathrm{C}}_{\mathrm{n}}[{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{+}{\mathrm{C}}_{\mathrm{t}}^{-1}{\mathrm{C}}_{\mathrm{b}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}\mathrm{[}\mathrm{I}\mathrm{-}\mathrm{(}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}{\mathrm{)}}^{-1}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}{\mathrm{C}}_{\mathrm{b}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}{\mathrm{]}}^{-1}\\ & & \end{array}$(25)where ${\mathrm{C}}_{\mathrm{t}}\mathrm{=}{\mathrm{C}}_{\mathrm{n}}\mathrm{+}\mathrm{C}b$(26)denotes the full (white+correlated) time domain noise covariance of un-destriped TOI. We bring $\mathrm{(}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}{\mathrm{)}}^{-1}$ inside the brackets and note that ${\mathrm{C}}_{\mathrm{t}}^{-1}{\mathrm{C}}_{\mathrm{b}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{=}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{-}{\mathrm{C}}_{\mathrm{t}}^{-1}$. Equation (25) simplifies to ${n}\mathrm{=}{\mathrm{C}}_{\mathrm{n}}{\mathrm{C}}_{\mathrm{t}}^{-1}{{n}}^{\mathrm{\prime }}\mathrm{+}{\mathrm{C}}_{\mathrm{b}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}\mathrm{[}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}{\mathrm{]}}^{-1}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}{{n}}^{\mathrm{\prime }}\mathit{.}$(27)Assuming that the true noise obeys the noise model, we have C_t = ⟨ n′n^′T ⟩, and we obtain for the time domain covariance of the destriped noise $\begin{array}{ccc}\mathrm{⟨}{n}{{n}}^{\mathit{T}}\mathrm{⟩}& \mathrm{=}& {\mathrm{C}}_{\mathrm{n}}{\mathrm{C}}_{\mathrm{t}}^{-1}{\mathrm{C}}_{\mathrm{n}}\mathrm{+}\mathrm{P}\mathrm{[}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}{\mathrm{]}}^{-1}{\mathrm{P}}^{\mathrm{T}}\\ & & \mathrm{-}{\mathrm{C}}_{\mathrm{n}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}\mathrm{[}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}{\mathrm{]}}^{-1}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}{\mathrm{C}}_{\mathrm{n}}\mathit{.}\end{array}$(28)Inserting this into Eq. (9) we finally obtain for the harmonic noise covariance matrix $\begin{array}{ccc}\mathrm{C}& \mathrm{=}& \mathrm{\left(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}{\mathrm{)}}^{-1}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{A}\mathrm{\right(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}{\mathrm{)}}^{-1}\\ & & \mathrm{+}\mathrm{\left(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}{\mathrm{)}}^{-1}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}\mathrm{\right(}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}{\mathrm{)}}^{-1}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}\mathrm{(}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}{\mathrm{)}}^{-1}\\ & & \end{array}$(29)Equation (29) is the main result of this section.

4.3. Cross-correlation

As discussed in Sect. 2.2, there are two destriping options available when deconvolving data subsets. Either we destripe the full data set (full destriping) or only the deconvolved data set (independent destriping). In the former case, the common destriping step generates correlation in residual noise between data sets. This effect too we can assess through the noise covariance formalism. We can construct a cross-covariance matrix ${\mathrm{C}}_{\mathrm{12}}\mathrm{=}\mathrm{⟨}{{a}}_{\mathrm{1}}{{a}}_{\mathrm{2}}^{\mathit{T}}\mathrm{⟩}$(30)where a₁ and a₂ are the harmonic coefficients obtained by deconvolution of the two data sets. The cross-covariance matrix is constructed in very much the same way as the covariance matrix of Eq. (29), only now we have on both sides two destriping operations A₁ and A₂. Formally both operate on the same full data stream, but yield zeros when operating to a part of TOI that does not belong to the data set. If the two data sets do not overlap, we have ${\mathrm{A}}_{\mathrm{1}}^{†}{\mathrm{C}}_{\mathrm{t}}^{-1}{\mathrm{A}}_{\mathrm{2}}\mathrm{=}\mathrm{0.}$(31)With this assumption, the cross-covariance matrix is given by $\begin{array}{ccc}{\mathrm{C}}_{\mathrm{12}}& \mathrm{=}& \mathrm{\left(}{\mathrm{A}}_{\mathrm{1}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}{\mathrm{A}}_{\mathrm{1}}{\mathrm{)}}^{-1}[{\mathrm{A}}_{\mathrm{1}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}\mathrm{\right(}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}{\mathrm{)}}^{-1}{\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}{\mathrm{A}}_{\mathrm{2}}\\ & & \end{array}$(32)The auto-covariance for the full destriping case is given by the same formula (29) as for independent destriping, only P is taken to be the pointing matrix of the full data set, while in case of independent destriping it is interpreted as the submatrix that contains the rows for the data set in question.

4.4. Implementation

Equations (29) and (32) appear complicated, but actually consist of only five ingredient matrices that appear in different combinations: ${\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}$, ${\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{A}$, ${\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}$, ${\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}$, and ${\mathrm{P}}^{†}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}$. We can further reduce the computation load as we describe in the following. As described in Keihänen & Reinecke (2012), artDeco performs deconvolution through a process where the TOI is first compressed into a “3D map”. The data is binned on a 3-dimensional grid, spanned by the pointing angles θ,φ,ψ. In θ and φ, artDeco are pixelized as in HEALPix, and ψ is divided uniformly into N_psi cells. We can formally write the deconvolution operation $\mathrm{A}{\mathrm{C}}_{\mathrm{n}}^{-1}$ as ${\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{=}\mathrm{G}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathit{,}$(33)where S^T, operating on a TOI, represents the operation of binning the TOI into a 3D map, and G the remaining part of the deconvolution operation.

Matrix P^T represents the operation of binning a TOI into a two-dimensional sky map. This too can be split into a two-step operation where one first constructs a 3D map, and from that a sky map. We write ${\mathrm{P}}^{\mathrm{T}}\mathrm{=}{\mathrm{X}}^{\mathrm{T}}{\mathrm{S}}^{\mathrm{T}}\mathit{.}$(34)An element of a 3D map is identified through a combined index qi, where q refers to sky pixel, and i = 0...N_psi−1 to the bin defined by beam rotation angle ψ_i. The elements of a sky map are referred to as Ip, Qp, or Up. In this notation, the X matrix is given by $\begin{array}{ccc}{\mathrm{X}}_{\mathit{qi,Ip}}& \mathrm{=}& {\mathit{\delta }}_{\mathit{q}\mathrm{\in }\mathit{p}}\\ {\mathrm{X}}_{\mathit{qi,Qp}}& \mathrm{=}& \cos \mathrm{\left(}\mathrm{2}\mathrm{\right(}{\mathit{\psi }}_{\mathit{i}}\mathrm{+}{\mathit{\psi }}_{\mathrm{pol}}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}{\mathit{\delta }}_{\mathit{q}\mathrm{\in }\mathit{p}}\\ {\mathrm{X}}_{\mathit{qi,Up}}& \mathrm{=}& \sin \mathrm{\left(}\mathrm{2}\mathrm{\right(}{\mathit{\psi }}_{\mathit{i}}\mathrm{+}{\mathit{\psi }}_{\mathrm{pol}}\mathrm{\left)}\mathrm{\right)}\hspace{0.17em}{\mathit{\delta }}_{\mathit{q}\mathrm{\in }\mathit{p}}\end{array}$(35)where ψ_pol is the detector-specific angle of polarization sensitivity. To add to the generality, we allow the 3D map pixel resolution to be higher than that of the sky map. In this notation δ_{q ∈ p} = if pixel q belongs into the larger pixel p. However, in our computations we always use same resolution for both, so that δ becomes the usual Kronecker delta. Operation X^T is thus a trivial coaddition operation involving weighting by the polarization angle.

The five ingredient matrices become $\begin{array}{ccc}{\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{A}& \mathrm{=}& {\mathrm{G}}^{†}\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{S}\mathrm{\right)}\mathrm{G}\\ {\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{A}& \mathrm{=}& {\mathrm{G}}^{†}\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{S}\mathrm{\right)}\mathrm{G}\\ {\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{P}& \mathrm{=}& {\mathrm{G}}^{†}\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{S}\mathrm{\right)}\mathrm{X}\\ {\mathrm{A}}^{†}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}& \mathrm{=}& {\mathrm{G}}^{†}\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{S}\mathrm{\right)}\mathrm{X}\\ {\mathrm{P}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{P}& \mathrm{=}& {\mathrm{X}}^{\mathrm{T}}\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{S}\mathrm{\right)}\mathrm{X}\mathit{.}\end{array}$(36)In the end, there remain three matrices to evaluate: $\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{S}\mathrm{\right)}$, $\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{S}\mathrm{\right)}$, and G. Operation X is carried out algorithmically where needed. Matrix G is part of the artDeco deconvolver, and we can extract it from the code. Matrix $\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{S}\mathrm{\right)}$ represents the white noise covariance between elements of a 3D map. It is diagonal, with elements given by σ²/N_pi, where σ² is the white noise variance in TOI domain, and N_pi is the number of hits to bin pi.

To proceed with matrix $\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{S}\mathrm{\right)}$, we first use the Woodbury formula to rewrite ${\mathrm{C}}_{\mathrm{t}}^{-1}$ as ${\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{=}\mathrm{(}{\mathrm{C}}_{\mathrm{n}}\mathrm{+}{\mathrm{C}}_{\mathrm{b}}{\mathrm{)}}^{-1}\mathrm{=}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{-}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{(}{\mathrm{C}}_{\mathrm{b}}^{-1}\mathrm{+}{\mathrm{C}}_{\mathrm{n}}^{-1}{\mathrm{)}}^{-1}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathit{.}$(37)This formulation has the advantage that C_n appears as inverse, which makes dealing with flagged data easier. The operation of Eq. (37) can be interpreted as a filtering operation in time domain. To construct one column $\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{S}\mathrm{\right)}$, we pick a TOI from a 3D map with one non-zero element, apply the filter of Eq. (37) to it, and compress the TOI back into a 3D map. The FFT technique is usually efficient in filtering, but since we are dealing with TOIs with very few non-zero elements, we find more efficient a procedure where we first use FFT to construct the filter, but then to perform the actual filtering in time domain. When dealing with Planck data we assume that noise is uncorrelated from one pointing period (typically one hour) to another, and do the filtering by pointing period.

We further benefit from the fact that matrices of Eq. (36) are additive. If we have two sections of TOI, with no noise correlation between, we can compute the contribution from each section to the ingredient matrices independently, and coadd the results. In particular, several detectors are combined simply by coadding their matrices.

We construct the ingredient matrices for each radiometer quadruplet and each half-year survey, write them on disk, and then coadd them in different combinations in a separate step. That way we avoid doing repeated computations when dealing with different data combinations.

4.5. Flagged data

To add to the generality of the method, we allow for a case where part of the data is flagged as unusable. Madam handles flags by setting formally ${\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{=}\mathrm{0}$ for the flagged samples. This is easy to take into account in $\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{S}\mathrm{\right)}$. Matrix $\mathrm{\left(}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{S}\mathrm{\right)}$ requires one more manipulation step. Equation (37) is expensive to construct, since the middle matrix ${\mathrm{C}}_{\mathrm{b}}^{-1}\mathrm{+}{\mathrm{C}}_{\mathrm{n}}^{-1}$ is no more stationary. We therefore make the approximation $\mathrm{(}{\mathrm{C}}_{\mathrm{b}}^{-1}\mathrm{+}{\mathrm{C}}_{\mathrm{n}}^{-1}{\mathrm{)}}^{-1}\mathrm{\approx }\mathrm{(}{\mathrm{C}}_{\mathrm{b}}^{-1}\mathrm{+}{\mathrm{C}}_{\mathrm{wn}}^{-1}{\mathrm{)}}^{-1}\hspace{0.17em}$(38)where C_wn is the original white noise variance, where flags have not been applied. This can be evaluated efficiently using Fourier techniques. We thus have for the last ingredient matrix ${\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{t}}^{-1}\mathrm{S}\mathrm{\approx }{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{S}\mathrm{-}{\mathrm{S}}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{(}{\mathrm{C}}_{\mathrm{b}}^{-1}\mathrm{+}{\mathrm{C}}_{\mathrm{wn}}^{-1}{\mathrm{)}}^{-1}{\mathrm{C}}_{\mathrm{n}}^{-1}\mathrm{S}\mathit{.}$(39)In absence of flags the formula is exact. The approximation misestimates the correlation between distant 3D map elements, but captures the main effect of flagging: the reduction in the amount of available data and the consequent increase in residual noise.

4.6. Destriping mask and detector weighting

The mapmaking procedure of Planck LFI involves a destriping mask (Planck Collaboration VI 2016). The mask is used for the purpose of reducing errors arising from high signal gradients and from bandpass mismatch between detectors.

The proper way of treating the mask in the NCVM would be to introduce another white noise covariance C_mask, whose inverse ${\mathrm{C}}_{\mathrm{mask}}^{-1}\mathrm{=}\mathrm{0}$ for samples that fall in the masked region. The masked covariance is then applied in the destriping phase, while the normal white noise covariance is applied in the deconvolution phase. Unfortunately, unlike in the case of flags, the different covariance terms appear in combinations that do not cancel out, with the result that the noise covariance matrix of Eq. (29) is no more valid, and not easily enhanced. Note that replacing the pointing matrix P by its masked version does not have the desired effect, since reducing P has the effect of reducing the residual noise, as there is a smaller number of pixels to solve, while the effect of destriping mask is to increase the residual noise level.

The mask thus affects the residual noise in a way that is not captured by our noise covariance matrix. This has to be accepted as a source of uncertainty in noise modelling.

In principle, a similar limitation arises from detector weighting. The destriping procedure of Planck LFI replaces the white noise covariance C_n by a slightly different modified covariance, which is made equal between a pair of detectors. This is referred to as horn-uniform weighting. The effect from this is expected to be small.

4.7. Inputs

The required inputs for NCVM computation are detector pointing, harmonic beam expansion, and noise spectrum. It is not required that noise parameters remain the same throughout the mission, although in our simulations we have assumed so. The TOI can be assembled of sections with different noise spectra, and each section’s contribution added to the ingredient matrices.

Four user parameters control the computation. Parameters ℓ_max and k_max control the deconvolution process and enter matrix G. Parameters ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ and ${\mathit{N}}_{\mathrm{psi}}^{\mathrm{3}\mathrm{D}}$ set the resolution of the 3D map dimension of matrices S and G. The meaning of ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ is twofold. On the one hand it sets the internal resolution of the deconvolution calculation. As a rule of thumb, 1.5${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ should exceed ℓ_max. Thus, for deconvolution with ℓ_max = 50 one should have at least ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ = 32, and for ℓ_max = 100 at least ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ = 64. On the other hand, ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ also sets the destriping resolution. The covariance matrix models a destriping procedure where the sky signal is assumed constant within a pixel of resolution N_side = ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$.

The required computational resources increase rapidly with increasing ℓ_max and ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$. The complexity of the computation increases as ℓ_max⁵, and the size of the matrix as ℓ_max⁴. In practice we can produce a full covariance matrix for only the low multipoles .

4.8. Constructing NCVM for Planck LFI frequencies

We have computed NCVM products for Planck LFI frequencies, for a number of radiometer and survey combinations. A complete list is given in Table 4. For each radiometer and survey combination we produced two sets of NCVM products: one with the white noise matrix, given by Eq. (11), and one with the full NCVM matrix given by Eq. (29). We used the same pointing, noise parameters, and beam, as in the simulations of Sect. 2.

We constructed the harmonic NCVM matrix for multipoles up to ℓ_max = 50. Beam elements with high k only affect the high ℓ multipoles (Keihänen & Reinecke 2012). We chose k_max = 2 for this low-multipole computation. Value ℓ_max = 50 requires a minimum pixel resolution of ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ = 32 (110′). To realistically model the mapmaking process of Planck LFI, we should have ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ = 1024 (3.44′) , which is out of reach. As a trade-off, we used ${\mathit{N}}_{\mathrm{side}}^{\mathrm{3}\mathrm{D}}$ = 64 (55′) and ${\mathit{N}}_{\mathrm{psi}}^{\mathrm{3}\mathrm{D}}$ = 256. The effect of this approximation on matrix quality is assessed in Sect. 5.

For all cases, the set of NCVM products includes: a harmonic NCVM matrix with ℓ_max = 50 (930 MB), a noise bias estimate computed according to Eq. (15), and a pixel-pixel covariance matrix at resolution N_side = 16 (720 MB). All products include polarization. The construction of the pixel-pixel matrix from the harmonic matrix is described in Sect. 6.

In addition to the full frequency matrices for 30, 44, and 70 GHz, we produced matrices for 70 GHz horn pairs 18/23, 19/22, 20/21, and for individual years (years 1−4) involving all 70 GHz radiometers. For each of these, we considered both full and independent destriping.

The construction of the ingredient matrices as described in Sect. 4.4 took 2000 CPU hours for each year and horn pair, adding up to 24 000 CPU hours for the whole 70 GHz data set. From those we constructed the NCVM for different horn pair and year combinations. The process of constructing the final NCVM from the ingredients was dominated by the inversion of the middle matrix, and took about 1500 CPU hours per matrix.

The predicted noise bias for LFI frequencies is shown in Fig. 9. We show three auto spectra (TT, EE, BB) and three cross spectra (TE, TB, EB). The cross spectra are small, and are multiplied by 10 in the plot to show the structure more clearly.

Figure 10 illustrates the matrix structure in ℓ,m space. We show the real part TT block of the 30 GHz NCVM for lowest multipoles (ℓ< 7). The matrix elements are arranged according to the combined index i = ℓ² + ℓ + m. The strongest correlations appear between elements with same ℓ, or same m and ℓ separated by an multiple of 2.

5. Validation of the covariance matrix

5.1. Low-resolution simulations

We validated our NCVM products through MC simulations. A complete list of simulations is given in Table 5.

The simulations followed the same two-step procedure as the high-resolution simulations reported in Sect. 2. We first generated and destriped a noise stream with the Madam destriping code, and then ran the artDeco deconvolution code, to yield an array a_Xlm of harmonic coefficients. The noise and beam parameters are the same as in Sect. 2. This time we ran the deconvolution step with parameters ℓ_max = 50, k_max = 2, in line with the settings used in NCVM construction. With this low an ℓ_max value the deconvolution step is fast, and the computation time is dominated by generation of noise rather than deconvolution. We produced 100 noise realizations for each case studied. As before, we considered two types of input noise: white and realistic 1 /f noise.

We cannot expect the NCVM to perfectly model the residual noise, due to simplifications we have to do when constructing the matrix, either due to limitations of the formalism, or due to limited computing resources. We can list four factors:

• 1.

  The covariance matrix assumes destriping resolution of N_side = 64, while actual LFI mapmaking uses N_side = 1024.

• 2.

  LFI mapmaking involves a destriping mask, neglected by the NCVM.

• 3.

  LFI mapmaking applies horn-uniform weighting in the destriping phase, NCVM computation intrinsically uses noise weighting.

• 4.

  NCVM assumes a baseline length of one sample, while the actual baseline length is 1.0 or 0.25 s.

We performed two series of simulations, which we refer to as “ideal” and “realistic”. The destriping parameters are given in Table 5. The realistic simulations mimick the LFI mapmaking procedure as closely as possible: we use a destriping mask, apply horn-uniform weighting, and destripe at a high resolution N_side = 1024.

In the ideal simulations, the destriping parameters follow the intrinsic assumptions of NCVM. Destriping is performed at low resolution N_side = 64, no destriping mask is applied, and horn-uniform weighting is replaced by noise weighting. To exactly match the assumptions of NCVM computation, we ought to use a baseline length of one sample. This becomes, however, computationally very demanding. As a trade-off we selected a baseline length of a quarter of the one used in actual mapmaking.

The aim of the ideal simulation set is to verify that the NCVM is constructed correctly. The realistic simulation set tests how accurately the matrix produced models the actual LFI data products, given the approximations involved.

5.2. χ² tests

We validate the noise covariance matrices through χ² tests. We generated 100 realizations of harmonic vectors a_Xℓm, as described above. For each vector and corresponding matrix we compute the value ${\mathit{\chi }}^{\mathrm{2}}\mathrm{=}\frac{{{a}}^{\mathrm{T}}{\mathrm{C}}^{-1}{a}}{{\mathit{N}}_{\mathrm{d}\mathrm{.}\mathrm{o}\mathrm{.}\mathrm{f}\mathrm{.}}}\mathrm{·}$(40)Here C is the covariance matrix, and N_d.o.f. is the number of degrees of freedom. Multipole ℓ contributes 2ℓ + 1 degrees of freedom, except for the monopole and dipole, which are identically zero for E and B. This yields in total N_d.o.f. = 3(ℓ_max+1)²−8 degrees of freedom, which for ℓ_max = 50 gives N_d.o.f. = 7795. For a matrix that perfectly describes the noise properties, the results obey the Chi distribution, For large values of N_d.o.f. the distribution becomes Gaussian, with expectation value 1 and standard deviation of $\sqrt{\mathrm{2}\mathit{/}{\mathit{N}}_{\mathrm{d}\mathrm{.}\mathrm{o}\mathrm{.}\mathrm{f}\mathrm{.}}}$, which for N_d.o.f. = 7795 is 0.016.

We compared the white noise covariance matrix of Eq. (11) against white noise simulations, and the full noise matrix of Eq. (29) against ideal and realistic simulations. Results from tests at the three Planck LFI frequencies, 30, 44, and 70 GHz are shown in Fig. 11. We show a histogram of the χ² values for 100 realizations, and quote their mean. The white noise matrix models the white noise residual nearly perfectly, χ² mean ranging from 0.9992 to 1.0004, Comparison between the full matrix and ideal simulation gives slightly larger values, from 1.0002 at 70 GHz to 1.0037 at 44 GHz, while the expected one-sigma deviation for 100 realizations is 0.0016. Deviations from one are larger in the case of realistic simulation, as expected, given the approximations involved. At 30 GHz we find a χ² mean of 1.0196, at 70 GHz 1.0037.

Results for 70 GHz subsets of are shown in Fig. 12. In the first column we compare the white noise matrix against a pure white noise simulation. Second and third column show results for the full matrix. Results are shown for two destriping options: full 70 GHz destriping, and independent destriping.

In all cases, we find very good agreement with the white noise matrix and corresponding simulations. The same can be said about the full matrix and the ideal simulations. The χ² mean is at maximum 1.0038. In the case of realistic simulation, deviations are again larger. The agreement is, however, still quite good for cases where deconvolution and destriping involve the same data (independent destriping), the maximal χ² value being 1.0079.

The case of full destriping, where the whole 70 GHz data set is destriped together, but only a subset of the cleaned data is provided to deconvolution, turns out to be more difficult to model with the covariance matrix. While the ideal simulation still shows good agreement, the mean χ² results for the realistic simulation range from 1.0146 to 1.0495. The discrepancy is largest for horn pair LFI18/23, which we now pick up for a closer examination.

We performed yet another series of simulations, where we applied full destriping, and deconvolved the LFI18/23 horn pair. We took as starting point the ideal simulation, and turned on, one at a time, each of the features that distinguish the realistic simulation from the ideal one. That is: destriping resolution N_side = 1024 instead of 64, horn-uniform weighting instead of noise weighting, baseline length 1 s instead of 0.25 s, and destriping mask. We generated 50 noise realizations for each case. The results of the χ² test for these simulations are shown in Fig. 13. We observe that detector weighting and destriping resolution have very little effect on the result. The finite baseline length and destriping mask both contribute significantly. Both effects alone raise the χ² mean to 1.03.

As discussed in Sect. 4.6, the mask affects the residual noise in a way that we cannot model in the covariance matrix. The effect of baseline length, instead, could be cured by running mapmaking with a shorter baseline.

5.3. Effect of ℓ_max

The noise covariance matrix models a deconvolution process where the harmonic coefficients are solved up to multipole ℓ_max. For a good match between the matrix and the harmonic coefficients, the matrix and the actual deconvolution must use the same value of ℓ_max. One might prefer running deconvolution with a high ℓ_max to avoid aliasing effects, and then pick the lowest multipoles for a low-resolution data set. This, however, degrades the match between the matrix and the data.

We illustrate this in Fig. 14. We deconvolved the 10 realizations of our 70 GHz full noise simulation with different ℓ_max, pick the lowest ℓ = 0 − 50 multipoles, and ran a χ² test with the ℓ_max = 50 NCVM. The matrix gives a good match with the coefficients obtained with ℓ_max = 50, but at ℓ_max = 100 there is already a significant deviation. Going from ℓ_max = 200 to ℓ_max = 1500 changes the result only little. This is in line with results obtained for white noise, shown in Fig. 8.

It is beyond the scope of this work to answer the question if the low ℓ_max can cause signal distortion. Should this be the case, we still have a few options left. We can improve the match by constructing the NCVM to higher ℓ_max. Computational cost increases rapidly with ℓ_max, but ℓ_max = 50 used in this work is by no means an upper limit. Luckily, Fig. 14 suggests that the resolution effect has already largely saturated at ℓ_max = 100, so increasing ℓ_max even a little may already solve the problem.

Another option worth exploring is to construct the full matrix as combination of two components, a white noise component, computed to higher resolution (as we have done in Sect. 3) and from a full noise component computed to lower resolution.

5.4. Validation through noise bias

Another way of validating the NCVM product is to compare the predicted noise bias of Eq. (15) against the MC spectrum of Eq. (5). Figure 15 shows the comparison for the full 70 GHz data. We show all six spectral components. The MC spectra shown are averaged over 100 noise realizations. We show both the idealised and the realistic simulation. For auto-spectra (TT, EE, BB) we show also the white noise simulation and the corresponding noise bias.

The predicted auto-spectra agree reasonably well with the simulation, except for multipole ℓ = 2, where the noise bias underestimates the true noise level. In the case of the much smaller cross spectra (TE, TB, EB), the 100 realizations do not provide enough statistics for a meaningful comparison. In the subsequent analysis we therefore concentrate on TT and EE spectra.

In Fig. 16 we show the TT and EE spectra for 30 and 44 GHz channels. The agreement is again good, and the deviation at ℓ = 2 is smaller than in case of 70 GHz. The difference between ideal and realistic simulations is small.

Figure 17 shows a comparison for horn pair LFI18/23. We consider again two destriping options. Simulations show that full destriping yields a much lower residual noise level, and this is captured by the NCVM as well.

Figure 18 shows the cross spectrum between 70 GHz horn pairs 18/23 and 19/22. When the horn pairs are destriped independently, there is no correlation in the noise. When the whole 70 GHz data set is destriped together, there is a significant cross-correlation at low multipoles. The NCVM models this correctly. We show in the same plot also the auto spectrum for horn pair 18/23 (same as in Fig. 17). The theoretical noise bias agrees well with the idealised simulation, but underestimates the realistic noise by a few per cent, in agreement with the χ² result for the same case.

6. Covariance in pixel domain

The primary output of artDeco is an array of harmonic coefficients. Consequently, the primary noise covariance matrix is constructed in harmonic space.

From the harmonic coefficients one can further construct a sky map via spherical harmonic transform. The harmonic coefficients must be multiplied by a smoothing kernel before performing the expansion. Smoothing serves two purposes. Firstly, it removes ringing artefacts that arise around point sources and other regions of high signal gradients. Secondly, smoothing compensates for the rise in noise level towards high multipoles.

6.1. High-resolution maps

We consider first high-resolution maps. We do not attempt to build a covariance matrix for high resolution, but study the noise properties through Monte Carlo simulations.

We transformed our high-resolution simulations to pixel space. Eliminating ringing artefacts in general requires a smoothing width larger than the average width of the actual beam. The same conclusion is supported by noise analysis. As shown in Sect. 3, the deconvolved noise spectrum rises more steeply than the inverse symmetrized beam window. Smoothing by this window will thus not be sufficient.

We applied a Gaussian smoothing kernel of width FWHM = 40′ (30 GHz), 30′ (44 GHz), or 20′ (70 GHz). The widths were chosen so as to remove ringing effects at each respective frequency. Figure 19 shows the TT noise spectra with smoothing applied. All spectra vanish at high multipoles.

Figure 20 illustrates the effect of deconvolution on noise in the pixel domain. We show a simulated 1 Jy point source at 30 GHz, located just above the origin of the Galactic coordinate system. In order to have the source in the centre of an N_side = 1024 pixel we place it not at the exact Galactic centre, but at (0.0373°, 0) (lat, lon), which is the centre of the pixel just above. We show also as example one realization of noise at the same location. The noise in the binned map is dominated by white noise. The smoothing operation associated with deconvolution suppresses the noise level in the deconvolved map.

Owing to the smoothing, the image of a point source on a deconvolved map is not point-like. It is, however, perfectly symmetric. At the same time the deconvolution procedure generates correlations in the noise component, in the direction perpendicular to the direction of beam main axis. To quantify this, we compute an estimate of the correlation of the central pixel with its surroundings, as ${\mathit{C}}_{\mathrm{map}\mathit{XYj}}\mathrm{=}\mathrm{⟨}{\mathit{n}}_{\mathit{Xi}}{\mathit{n}}_{\mathit{Yj}}\mathrm{⟩}$(41)where X,Y = I,Q,U are the Stokes components, n_j is the noise in pixel j, index i represents the reference pixel, and j labels all other pixels of the sky. We evaluate C_map as an average over the 40 available noise realizations. The II component map is plotted in Fig. 22. As expected, noise is correlated in a direction at right angles to the beam main axis.

6.2. Low-resolution maps and pixel NCVM

We proceed to produce low-resolution maps, and related covariance matrices. We transformed the low-resolution simulations of Sect. 5 to pixel space. We adopted the smoothing procedure that was applied to LFI 2015 release data (Planck Collaboration VI 2016). We applied Gaussian smoothing with FWHM = 440′ and constructed maps at resolution N_side = 16 (220′).

The same operations must be applied to the covariance matrix. The smoothing operation reduces the effective number of degrees of freedom, making the covariance matrix singular or badly scaled. This is cured by adding a small amount of white regularization noise within each pixel of the map. A corresponding constant variance was added on the diagonal of the covariance matrix. Following the LFI procedure, we added white regularization noise with rms = 2 μK (I) or rms = 0.02 μK (Q and U) on each pixel.

The operation on map can formally be written as ${m}\mathrm{=}\mathit{H}{a}\mathrm{+}{{n}}_{\mathrm{reg}}\mathit{.}$(42)Here H represents the harmonic expansion combined with smoothing, a is the vector of harmonic coefficients, as produced by artDeco, m is the beam-deconvolved map, and n_reg represents regularization noise.

The noise in map of Eq. (42) is characterized by the covariance matrix ${\mathrm{C}}_{\mathrm{pix}}\mathrm{=}\mathit{H}\mathrm{C}{\mathit{H}}^{\mathrm{T}}\mathrm{+}{\mathrm{N}}_{\mathrm{reg}}$(43)where C is the harmonic covariance matrix, and N_reg is diagonal matrix representing the regularization noise.

The operation HCH^T is computationally intensive and broken down into the following individual steps:

• 1.

  Each row of C is transformed via SHT. This is done using the efficient libsharp implementation (Reinecke & Seljebotn 2013). Since this library only supports real-valued maps, two SHTs (for the real and imaginary part, respectively) must be carried out for each row; depending on the parameter choice, these are either scalar or fully polarized. Since a matrix row is located entirely on a single MPI task, no communication is necessary for this step.

• 2.

  To perform the analogous operation on every column of the intermediate result CH^T, this matrix must be effectively transposed in memory, so that columns reside on a single MPI task, as the rows did before. This is a communication-intensive step.

• 3.

  SHTs are carried out over the rows of (CH^T)^T in a fashion analogous to step 1, resulting in ${\mathrm{C}}_{\mathrm{pix}}^{\mathit{T}}$.

• 4.

  In another communication step, this matrix is again transposed.

Figure 22 illustrates the structure of the pixel covariance matrix. We show as a map one column of the 70 GHz covariance matrix, which represents II correlation between the reference pixel and the rest of the sky. We pick the reference pixel just above the galactic centre. The reference pixel is centred at (2.388°, 0) and has pixel number 1440 in the ring pixeling scheme of Healpix at resolution N_side = 16.

The reference pixel is strongly correlated with the neighbouring pixels. A weaker correlation structure follows the scanning direction. The map has a significant positive minimum of roughly 0.012 μK², reflecting the difficulty of distinguishing the map monopole from a global noise offset. The correlation structure may vary a lot according to the reference pixel chosen. The plot shown here represents one example.

6.3. χ² tests in pixel space

We validate the pixel covariance matrices through χ² tests, similarly to the validation tests presented in Sect. 5. The number of degrees of freedom is now given by the number of pixels in the Healpix map, times 3 Stokes components, which for N_side = 16 gives N_d.o.f. = 9216. Consequently, the one-sigma statistical variation becomes 0.0147, and rms of mean for 100 realizations becomes 0.00147.

As in Sect. 5, we compared the white noise covariance matrix against white noise simulations, and full covariance matrix against ideal and realistic simulations. The χ² results for full four-year frequency maps and for 70 GHz subsets are shown in Figs. 23 and 24.

The mean reduced χ² values for pixel covariance are in general closer to one than in the case of harmonic covariance. This is the effect of regularization noise. We are adding to the maps a noise component that is accurately modelled by the extra term in covariance matrix, thus improving the agreement with the matrix and the maps.

The χ² mean is 1.01 for 30 GHz 4-yr data, 1.0076 for 44 GHz, and 1.0016 for 70 GHz. The latter is only slightly larger than the one-sigma statistical variation. For 70 GHz subsets the values range from 1.0017 to 1.0287.

6.4. Pre-whitened noise

The power of the χ² test is limited in the sense that it collapses a complicated noise structure into a single number, which is not constrained from below. It is easy to conceive a situation where two errors cancel to produce a χ² result close to one. For instance, a simple overestimation of the noise rms in TOI domain leads to an artificially low χ² result.

We complemented the validation tests by analysis of pre-whitened noise. The types of tests we performed here require a real-valued data set. For this reason we can perform the test in pixel space only.

Consider a pre-whitened noise map, computed as ${x}\mathrm{=}{\mathrm{M}}^{-1}{m}$(44)where m is the noise map, and M is the “square root” of the covariance matrix C, for which C = MM^T. Matrix M is can be constructed in multiple ways. We took it to be the Cholesky decomposition of C. If matrix C accurately describes the properties of the noise map m, x will be a vector of uncorrelated Gaussian noise with unit variance, ⟨ xx^T ⟩ = I.

We build histograms of values of x and compare them to theoretical Gaussian distribution, in a fashion similar to the analysis presented in Stompor et al. (2002). We combine the 100 noise realizations into the same plot to improve statistics. The variance of the distribution yields the usual χ² test. Comparison of the full distribution against the expected Gaussian distribution provides a more stringent test.

We show the distributions for two selected simulation cases in Fig. 25, both with realistic noise. The first example is one where the NCVM models the noise well, the 70 GHz frequency map. The second example is the worst case from earlier tests, the LFI18/23 horn pair map with full destriping. In both cases, noise before whitening is very non-uniform. We show for comparison the Gaussian distribution with variance equal to that of the data. The variance is dominated by the high tails of the distribution, leading to a broad distribution that does not fit the data anywhere. On the right we show the pre-whitened noise, together with the ideal Gaussian distribution with unit variance. In the case of 70 GHz frequency map, the distribution of whitened noise follows the Gaussian closely, indicating that the NCVM models the noise well. In the horn pair case there is a visible deviation. There are no high tails in the distribution, the deviation manifests rather as a overall broadening of the distribution. In fact, the distribution matches very well a Gaussian shape with standard deviation of 1.015.

To be more quantitative, we performed Kolmogorov-Smirnov test on the pre-whitened noise. We constructed the cumulative distribution of the values of x, and found the maximum difference with respect to the Gaussian distribution. The probability to exceed (PTE) values computed according to Kolmogorov-Smirnov statistics are given in Table 6. The number of degrees of freedom is N_e = 921 600 (100 maps of 9216 pixels each). In general, the K-S test yields the same conclusions as the χ² test. In white noise and ideal simulation cases the NCVM models the noise well. The worst match is found for 70 GHz subsets where we combine realistic simulation with full destriping. The 18/23 horn pair again stands out as the worst case. The full frequency maps and submaps with independent destriping are modelled well. We can identify a couple of case where the χ² test indicates a better match than the K-S test. In particular, for horn pair 19/22 with independent destriping and realistic simulation we get the excellent result χ² = 1.0003, but only 0.01 in the K-S test. This is likely an example of a case where two opposite errors compensate in the χ² test.

The distribution of pre-whitened noise for horn pair 18/23, shown in Fig. 25, indicates that the error made in NCVM computation when ignoring non-idealities of realistic data, results in an overall underestimation of the noise, and can to some degree be accounted for scaling the NCVM by an appropriate constant. To test this hypothesis, we performed another round of K-S tests where we rescaled the pre-whitened noise x by its stdev value, bringing the variance to unity. This is equivalent to scaling the NCV matrix by the variance of x. The required scaling factor is directly given by the χ² result, and the χ² result for this modified test is by construction exactly one. The K-S test can still reveal discrepancies in the shape of the noise distribution. The PTE values from the modified test are given in the last column in Table 6. We see that a simple rescaling of the covariance matrix leads to a much better agreement between the matrix and the maps.

7. Conclusions

We performed an analysis of noise properties of beam deconvolved maps, both in the harmonic and pixel domain. The analysis is built on the formalism of the artDeco deconvolution code. As the starting point we have a situation where a correlated 1 /f type noise stream is first destriped to remove the correlated component as accurately as possible and then fed as input to the deconvolver.

The main result of this work is the low-resolution noise covariance matrix; this matrix describes the residual noise in the harmonic coefficients which is the primary output of artDeco. The covariance matrix models the contribution of correlated residual noise that is left in the data after destriping and which contributes strongly to the residual at lowest multipoles. The covariance matrix also captures the effect of the deconvolution process itself on the noise. From the harmonic covariance matrix we further constructed a pixel-pixel covariance matrix, which describes the noise properties of deconvolved low-resolution maps.

As a practical application we computed covariance matrices for the Planck LFI four-year mission, using real LFI beams and noise spectra, and reconstructed detector pointings. We constructed the matrices up to multipole ℓ_max = 50. We constructed further pixel-pixel covariance matrices at resolution N_side = 16.

We compared the covariance products against Monte Carlo simulations. For validation purposes we performed a series of ideal simulations, which closely follow the intrinsic assumptions of the covariance matrix. The mean reduced χ² values in harmonic space range from 0.9991 to 1.0071, while the one-sigma statistic variation is 0.0016. In pixel space the agreement is even better, ranging from 0.9986 to 1.0039.

The covariance matrices are not expected to perfectly model realistic LFI noise due to the simplifications we have to do in matrix construction. We performed another series of realistic simulations to find out how accurately the matrices model real LFI data. The mean reduced χ² values in harmonic space are 1.0037 for 70 GHz, 1.0138 for 44 GHz, and 1.019 for 30 GHz. In pixel space the corresponding numbers are 1.0016 (70 GHz), 1.0076 (44 GHz), and 1.01 (30 GHz).

Various subsets of 70 GHz data gave χ² results ranging from 1.0003 to 1.0495. We studied more closely the case with the largest discrepancy, and found that the main contributors to the discrepancy are the use of a Galactic mask, and finite baseline length, both related to destriping. The mask is required in LFI analysis to control the signal error, but our NCVM does not model its effect correctly. The inaccuracy arising from baseline length could be removed by running the destriping step with a shorter baseline.

We considered two destriping options. In full destriping the full data of a frequency channel is destriped together, and a subset of it used for deconvolution. In independent destriping we use the same data subset in both processing steps. The NCVM is found to model the residual noise well in the latter case. The full destriping option yields a lower general noise level, with the side effect that the systematic error sources become more important. Consequently, the agreement between the NCVM and the data is poorer and the largest χ² values are obtained with this option. Fortunately, although the full destriping option yields the lowest residual noise for a particular submap, the independent destriping option is more important from the point of view of cosmological analysis. Noise covariance matrices are required input in many methods of cosmological parameter estimation. These methods typically input either the full frequency maps, or a set of submaps with independent noise, which we produced using the independent destriping option.

We complemented the validation tests with Kolmogorov-Smirnov tests in pixel space. We used the NCV matrix to pre-whiten the simulated noise maps, and test the hypothesis that the whitened noise is indeed white. We observe that in cases where the NCVM underestimates the noise, the discrepancy is in most cases, at least to a first approximation, well modelled by an overall rescaling of the matrix. The value of the optimal scaling factor remains to be determined from simulations.

It is instructive to compare the χ² results for low-resolution maps to similar tests performed on the official LFI data products, as presented in Planck Collaboration VI (2016). The official products include pixel-pixel covariance matrices at resolution N_side = 16. χ² results from for full mission frequency maps are 1.042 (30 GHz), 1.027 (44 GHz), and 1.010 (70 GHz). For the same data sets we obtain 1.01, 1.0076, and 1.0016, respectively. Although our tests indicate a better agreement, the results must be interpreted with care. Our simulation procedure mimics that of official LFI pipeline. We applied same smoothing width (440′), the same target resolution, and the same amount of regularization noise. There are, however, also important differences. The results of Planck Collaboration VI (2016) are based on 10 000 realizations with a time-dependent noise model, while in our simulations the noise properties remains the same throughout the mission. A reliable comparison would require tests performed on the same simulation data, and NCVM constructed with the same noise model. The FFP8 simulations used by Planck Collaboration VI (2016) are available in the form of high-resolution maps. For deconvolution, however, we need alternatively the TOI, or the 4D map objects constructed from it, which are not available

At higher multipoles the residual noise is dominated by the white noise component of TOI, deformed by the deconvolution process. Our simulations show that deconvolution leaves relatively less residual correlated noise on top of the white component at high multipoles than does the ordinary binning operation. From this it follows that (with the exception of 30 GHz where the correlated residual is still large) the full noise bias at high multipoles can be estimated to high accuracy from the white noise covariance matrix, which is cheaper to compute than the full covariance. Also, white noise simulations are much cheaper computationally than full noise simulations, which require the destriping step. A possible future development is to equip the deconvolution code with an internal noise generator, which fills the input data structure with white noise. This would avoid the bottleneck of writing the simulated TOI on disk, thus allowing extensive MC simulations.

We computed the white noise covariance to ℓ_max = 200. At higher multipoles we constructed partial matrices that cover a narrow multipole range around the region of interest. The noise bias obtained this way underestimates the noise spectrum obtained from Monte Carlo simulations by 1−3%. Should higher accuracy be required, the recommended approach is to estimate the bias through Monte Carlo simulations. That will require optimization of the MC pipeline, as described above.

In cases where a subset of available data is used for deconvolution, there are two options for destriping. One can either destripe the same data set, or the full data set. The second option significantly reduces the level of residual correlated noise, but also leaves the different subsets correlated. We derive a covariance matrix that describes the level of intercorrelation, and validate it on simulations involving 70 GHz horn-pair data sets.

Finally we examined the noise properties in deconvolved high-resolution maps at a qualitative level. We show that the residual noise is correlated in the direction perpendicular to the beam axis. The overall noise level is low owing to the smoothing operation involved in map construction.

Acknowledgments

Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package. M.R. is supported by the German Aeronautics Center and Space Agency (DLR), under pr.ogram 50-OP-0901, funded by the Federal Ministry of Economics and Technology. EK is supported by Academy of Finland grant 253204. V.L. and K.K. are supported by Academy of Finland grant 283497. A. S. has been supported by Väisälä (2014) and Wihuri (2015) foundations. We acknowledge that the results of this research have been achieved using the PRACE-3IP project (FP7 RI-312763) resource Sisu based in Finland at CSC. We thank CSC – the IT Center for Science Ltd. (Finland) – for computational resources. We thank H. Kurki-Suonio and J. Väliviita for careful reading of the manuscript and helpful comments. We thank E. Palmgren for help with plots.

References

All Tables

All Figures