A&A 395, 417-421 (2002)
DOI: 10.1051/0004-6361:20021288
A. Balbi ^{1,2} - G. de Gasperis^{1} - P. Natoli^{1,2} - N. Vittorio^{1,2}
1 - Dipartimento di Fisica, Università di Roma ``Tor
Vergata'', via della Ricerca Scientifica 1, 00133, Roma, Italy
2 - INFN, Sezione di Roma II, via della Ricerca Scientifica 1, 00133, Roma, Italy
Received 12 March 2002 / Accepted 3 September 2002
Abstract
We use an iterative generalized least squares map-making
algorithm, in conjunction with Monte Carlo techniques, to obtain
estimates of the angular power spectrum from cosmic microwave
background (CMB) maps. This is achieved by characterizing and
removing the instrumental noise contribution in multipole space.
This technique produces unbiased estimates and can be applied to an
arbitrary experiment. In this paper, we use it on realistic
simulations of Planck Low Frequency Instrument (LFI) observations,
showing that it can lead to fast and reliable estimation of the CMB
angular power spectrum from megapixel maps.
Key words: cosmic microwave background - methods: data analysis
State of the art measurements of the cosmic microwave background (CMB) anisotropy are affected by non negligible and (in the case of one horned experiments) correlated instrumental noise. In order to minimize this contaminant one resorts to non trivial statistical techniques which almost invariably require knowledge of the correlation structure of underlying detector noise, to be measured from the data themselves. As a consequence, methods to estimate the noise correlation properties out of time ordered data (TOD) have been proposed (Doré et al. 2001; Stompor et al. 2002; Natoli et al. 2002).
Most cosmological information is encoded in the angular power spectrum of the CMB anisotropies. However, the size of modern CMB datasets and the presence of correlated instrumental noise in the observations make it unfeasible to extract the power spectrum from CMB maps by performing standard and robust matrix manipulations commonly used to solve linear systems. Computationally, obtaining a brute force maximum likelihood estimation of the power spectrum requires at least operations, where is the number of pixels in the map. Current CMB maps from balloon-borne experiments have -10^{5}. Brute force power spectrum estimation from these maps is already prohibitive, and will become totally unattainable for upcoming space missions such as NASA's MAP satellite^{} or ESA's Planck Surveyor^{}, whose maps will have .
A number of strategies have been proposed in the past to address the problem of power spectrum estimation in a computationally feasible way. The MADCAP package (Borrill 1999) uses a parallel implementation of a quadratic estimator technique (Bond et al. 1998) to obtain maximum likelihood estimates of the power spectrum. This method, however, will certainly be too time consuming for future satellite data sets. Doré et al. (2001) applied the same approach in a hierarchical fashion on subsets of large data sets, lowering somewhat the computational time requirements. Szapudi et al. (2001) adopted an entirely different strategy, extracting the power spectrum from the 2-point correlation function of the map. Although this approach is applicable to a generic data set and might in principle account for noise correlations, until now it has only been tested in the case of uniform white noise. Finally, other methods were based on simplifying assumptions tailored on specific experimental strategies (e.g., Oh et al. (1999) for the MAP satellite; Wandelt & Hansen (2001) for the Planck Surveyor).
In this work, we focus on characterizing and removing the instrumental noise contribution in multipole space, in order to obtain unbiased estimates of the CMB power spectrum. This approach is based on using map-making techniques to project estimates of the time stream noise on the sky, according to the experimental scanning strategy, and on Monte Carlo (MC) simulations to determine the statistical properties of the noise in multipole space. Rather than being based on a noise model or an approximation (such as, e.g., uncorrelated noise) the time stream noise properties are estimated directly from the data. A similar strategy was developed independently in the MASTER package by Hivon et al. (2002) and applied to the analysis of the BOOMERanG data (Netterfield et al. 2002). In this case, fast projection of the noise on the sky was achieved through non optimal map-making, and filtering of the data in time domain was required to reduce the correlated noise contribution. As such, MASTER does not attempt to create a useful map as an end-product. Instead, we use the iterative generalized least squares (IGLS) map-making algorithm described in Natoli et al. (2001) which is fast enough to allow for the generation of an appropriate number of simulations in a reasonable time. This algorithm has optimal statistical properties (i.e. produces unbiased and minimum variance estimates of the map) and hence it does not require any filtering of the time stream. As an application, we use realistic simulations of Planck Low Frequency Instrument (LFI) observations to show that we can successfully recover the CMB power spectrum even from megapixel maps. This is, to date, the first computationally feasible and unbiased pixel-based power spectrum estimator for Planck which uses information on the correct (i.e. estimated from the data themselves) instrumental noise covariance.
This paper is organized as follows. In Sect. 2 we discuss the method we use. In Sect. 3 we present the results of the application to Planck/LFI. Finally, in Sect. 4 we draw the main conclusions of this work.
We begin by reviewing some basic notions about the statistics of the
CMB on the sky sphere. The CMB temperature field as a function of the
direction of observation on the sky can be expanded in spherical
harmonics
as:
(1) |
(2) |
(4) |
(5) |
(6) |
(7) |
An approximate way to model the effect of partial sky coverage is to
change the number of DOF of the
distribution from
to
,
where
is the observed fraction of the sky.
The
estimates get biased by a
factor, and the rms
becomes (Scott et al. 1994; Hobson & Magueijo 1996):
Since the applications shown in this paper will regard estimation from nearly full-sky maps ( ), an exact treatment of the partial sky coverage effect is totally superfluous, as the approximation described above yields comparable accuracy. Indeed, we have verified that this simplified approach provides an excellent approximation to the exact treatment even when a considerable fraction of the sky remains unobserved. We have applied a rather severe Galactic cut ( symmetric around the equator, corresponding to ) to 100 Monte Carlo simulations of maps containing only CMB signal, and extracted the corresponding . Figure 1 shows the results of this test. When the estimates are corrected for the bias they are nearly undistinguishable from the true underlying theoretical model. From the same figure it is also apparent that Eq. (8) provides a very good approximation to the true dispersion of the estimated .
Figure 1: The effect of incomplete sky coverage on the power spectrum estimates. The blue dots are extracted from 100 Monte Carlo simulated maps containing only CMB signal, after removing a symmetric strip of around the Galactic equator ( ). The black dotted curves are the input theoretical used in the simulations and the upper and lower bounds obtained from Eq. (8). The blue continuous curves are the mean and upper and lower bounds estimated from the MC simulations. This shows that Eq. (8) is a good approximation even when a substantial fraction of the sky is unobserved. | |
Open with DEXTER |
The measurements performed by a given experiment can be thought of as a superposition of sky signal (S) and instrumental noise (N), which are independent statistical processes. When a map is produced from the observations, residual noise is left in the pixelized data. This noise contribution is minimal when a minimum variance map-making (such as the IGLS) is used to obtain the map.
We can write the observation in pixel p as:
(9) |
(10) |
(11) |
It can be easily shown that the rms uncertainty of this power
spectrum estimate is:
To determine the noise properties in multipole space, we first apply the IGLS map-making procedure on the TOD, ending up with a minimum variance map of the CMB and an estimate of the true noise properties in time domain. We then produce a number of MC realizations of time streams containing only instrumental noise, where the noise has the statistical properties measured from the data. These simulated time streams naturally incorporate all information about the experiment observational strategy. We apply again the IGLS map-making algorithm to these mock data sets, to construct a set of noise maps. These noise maps have the same statistical properties of the noise which contaminates the real map. They can then be used to characterize the noise statistical properties in multipole space, for example estimating the mean to be used in Eqs. (12) and (13).
Figure 2 shows the result of applying this procedure to
simulated observations of Planck/LFI 100 GHz channel. The simulations
were produced using the standard Planck scanning strategy (i.e. a
spinning frequency of 1 r.p.m. and
offset angle between the
pointing direction and spin axis of the telescope, with the latter
always on the ecliptic) as well as a realistic model of the receiver
noise (i.e. 1/f noise with
Hz and Planck
sensitivity goal, see e.g. Natoli et al. 2002). A comparison was
made with the case of non uniform uncorrelated noise. It is apparent
that such a simple model is not accurate enough to characterize the
noise properties in multipole space, particularly for small values of
(
). It is interesting to observe that, since we
only need to estimate an average from the Monte Carlo, i.e. the
to be used in Eqs. (12) and (13),
convergence is achieved after a fairly small number of simulations.
This is evident from Fig. 2.
Figure 2: Monte Carlo estimation of the noise properties in multipole space for Planck/LFI. The upper plot shows the average instrumental noise contribution to the power spectrum, . The blue curve was obtained from 22 realizations of maps containing only instrumental noise with realistic properties (including time correlations); the green curve is from 100 realizations of non uniform uncorrelated noise; the orange curve is from a subset of 22 realizations of the same non uniform uncorrelated noise maps. The simulations were performed using the nominal sensitivity of all combined 34 receivers at 100 GHz, for 7 months of observation. The lower plot shows the noise contribution to the error bars, from the same sets of simulations. | |
Open with DEXTER |
Figure 3 shows histograms of the noise angular
power spectrum
,
for four multipole values, obtained from
1000 Monte Carlo realizations Planck/LFI maps containing only
instrumental noise. In order to speed up the production of simulated
maps, we used an inhomogeneous white noise model in place of the
correlated noise used in the rest of the paper. This should not alter
the results of this test, since the deviation from the white noise
behaviour is relevant only at low multipoles (see Fig. 2),
where the noise contribution is sub-dominant with respect to the
cosmic variance. The comparison with the theoretical distribution,
i.e. a
with
degrees of freedom, shows that
use of Eq. (13) yields nearly optimal error bars for the
power spectrum estimates.
Figure 3: Distribution of the noise contribution in multipole space. The plots show the histogram of , normalized to its average value, extracted from 1000 Monte Carlo simulations Planck/LFI maps, containing only inhomogeneous uncorrelated instrumental noise, for four values of (from left to right and from top to bottom: 50, 200, 400, 1000). Also shown is the distribution with degrees of freedom (continuous line) and the corresponding 1 symmetric error bar, given by (blue lines). | |
Open with DEXTER |
Figures 4 and 5 show the results of applying our method to extract the CMB power spectrum from a simulation of Planck/LFI observations at 30 and 100 GHz. In simulating the observations we fully took into account the real Planck scanning strategy. The offset between the pointing direction and spin axis of the telescope leaves unobserved two areas of about 80 square degrees around the ecliptic poles. We modeled the optical response of the instrument as a symmetric Gaussian beam with FWHM of 33 arcmin for the 30 GHz detectors and of 10 arcmin for the 100 GHz detectors. The simulated TOD for the 30 GHz channel has time samples (14 months of observation), which are mapped into pixels. These numbers turn into (7 months of observation^{}) and for the 100 GHz channel (we remind that IGLS map-making scales approximately linearly with ). Going from the simulated data to the power spectrum estimates took about 12 h for the 30 GHz channel, and about 18 h for the 100 GHz channel, using a parallel implementation of the IGLS algorithm running on 16 processors of the Origin 3000 supercomputer at Cineca. These times are dominated by the production of the MC noise maps, 22 in our case. This number is enough to produce good estimates for the average noise contribution (see again Fig. 2). For a detailed discussion about memory requirements, CPU timing, and code scalability of our implementation of the IGLS map-making algorithm, see the paper by Natoli et al. (2001).
Comparison of our
estimates with the theoretical input model
yields a /d.o.f. of 1.06 and 1.10 for the 30 and 100 GHz case,
respectively. As mentioned in Sect. 2.1, we also estimated
the power
in bands b of
width
,
in order to reduce small spurious correlation
in multipole space. When we use a simple cubic spline algorithm to
interpolate these band-power estimates, the resulting curve is nearly
undistinguishable from the theoretical input model (see bottom panel
in Figs. 4 and 5).
Figure 4: CMB Power spectrum estimation for Planck/LFI 30 GHz channel. The simulated data were produced assuming the nominal sensitivity from all 4 combined radiometers and 14 months of observation (full mission). The top panel shows the power spectrum estimated directly from the map ( , green points), the MC estimation of the noise contribution ( , blue curve), and the recovered signal power spectrum ( , red points). The bottom panel shows the estimated band-powers and their 1 error bars (red boxes). The orange curve is the result of a cubic spline interpolation of the band-power estimates. In both panels the black curve is the input theoretical model used in the simulation (a COBE-normalized standard CDM). | |
Open with DEXTER |
Figure 5: Same as in Fig. 4, but for the 34 combined radiometers of Planck/LFI 100 GHz channel and 7 months of observation (half mission). | |
Open with DEXTER |
We have shown that IGLS map-making can be successfully applied, in conjunction with MC techniques, to the problem of estimating the angular power spectrum from CMB maps. The method discussed in this work is fast enough to be already applicable to megapixel maps such as those expected from the Planck Surveyor. No manipulation of the time stream (i.e. high-pass filtering) is required by this method. Furthermore, no unrealistic simplification of the instrumental noise behavior is needed, since the noise properties are estimated directly from the data. Our estimated noise angular power spectrum contains the right information on noise correlation, as well as on the details of scanning strategy, even if we never use explicitly the full pixel covariance matrix. Although we only presented in this paper an application to the case of nearly full-sky maps from the Planck Surveyor, the approach described here can be used for any arbitrary scanning strategy and sky coverage (for example to analyze balloon data), as long as the spectral information on the spatial window of the observation is taken into account.
Finally, we would like to mention a few possible developments of the work described in this paper. We considered a symmetric beam approximation in our analysis. This approximation may prove to be too simplistic when dealing with real experiments. Thus, we are currently generalizing our map-making algorithm to deal with asymmetric beam patterns. We also point out that this power spectrum estimation technique is easily generalizable to CMB polarization maps, since an implementation of the IGLS map-making algorithm for polarization observations exists (see Balbi et al. 2002). We shall discuss this topic in a forthcoming paper.
Acknowledgements
We thank F. K. Hansen and D. Marinucci for discussions and useful advice. The supercomputing resources used for this work have been provided by Cineca. We acknowledge use of the HEALPix package (http://www.eso.org/science/healpix/).