A S T U D Y O F N U M E R I C A L E R R O R S I N S I M U L A T I O N S O F T H E C O S M I C M I C R O W A V E B A C K G R O U N D P O L A R I Z A T I O N By Bell L L . Chen B . Sc. (Physics & Astronomy), University of British Columbia A THESIS SUBMITTED IN PARTIAL F U L F I L L M E N T OF T H E REQUIREMENTS FOR T H E D E G R E E OF M A S T E R OF SCIENCE in T H E F A C U L T Y OF G R A D U A T E STUDIES PHYSICS & A S T R O N O M Y We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH COLUMBIA October 2000 ©Bell L L . Chen, 2000 In presenting this thesis in partial fulfillment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. Physics & Astronomy The University of British Columbia 129-2219 Main Mall Vancouver, Canada V6T 1Z4 Date: Abstract There have been numerous studies of data analysis issues involving temperature anisotropics on the microwave sky, but far less attention paid to the polarization sig-nals. The production of maps and their reduction to power spectra proceeds by choosing a particular way of dividing the sphere into pixels. The Equidistant Cylindrical Projection (ECP) is the geometrically simplest pixelization scheme and might have been sufficient in the early days of Cosmology, when the uncertainties from numerical calculations were ususally overshadowed by experimental errors. But with future satellite missions (MAP and Planck) on the horizon, we need to make sure that the pixelization scheme we choose does not add to the small experimental errors, in order to determine the cosmological parameters as accurately as possible. Numerical errors are a small part of the whole po-larization data analysis process but one that is easily dealt with in comparison to other more complicated analyses of polarization. In this thesis, we begin by stating a consistent set of equations for calculating the C M B polarization sky and power spectrum. This will be implemented into an E C P scheme, which allows us to study the numerical errors intro-duced by various effetcs related to pixelization. We will show that all these errors can be optimally reduced by numerical techniques. And finally, we will apply these techniques to a HEALPix pixelization scheme to obtain the most accurate polarization power spectra. n Table of Contents Abstract ii List of Figures iv 1 Introduction 1 2 Simulation of Temperature and Polarization Maps 8 2.1 Spin-Two Spherical Harmonics Transform 10 2.2 Recurrence Relation of Legendre Polynomials 15 2.3 Pixelization 17 3 Inverse Spin-Two Spherical Harmonics Transform 21 4 Numerical Errors 23 4.1 Stability of the Recurrence Relation of Legendre Polynomials 23 4.2 Orthogonality Issues Related to Pixelization 26 5 Other Relevant Numerical Issues 34 5.1 Numerical Error as a Function of Map Resolution 34 5.2 Approximation of Inverse Transform Integration 35 5.3 E C P with Equal-Area Pixels 37 6 Conclusion 38 References 39 in List of Figures 1-1 Southern C M B sky observed by B O O M E R A N G 4 1-2 C M B data with a fitted Cosmological model 4 1- 3 Some C M B power spectra produced by C M B F A S T 6 2- 1 A map of C M B temperature anisotropy 10 2-2 A map of C M B Polarization: parameter Q 13 2-3 A map of C M B Polarization: parameter U 13 2-4 Spin-two spherical harmonics at removable singularities 17 2-5 HEALPix and E C P of the same resolution in Aitoff projection 20 4-1 Associated Legendre Polynomials of even degree 24 4-2 Errors in recurrence relation: Legendre polynomials 25 4-3 Errors in recurrence relation: spin-two spherical harmonics 25 4-4 Leakage in E and B polarization spectra 27 4-5 First-order approximation results 29 4-6 First-order approximation residual 30 4-7 First-order approximation results for small B polarization 31 4- 8 Second-order approximation results 33 5- 1 Effect of sampling on overall error 35 5-2 Diagram for a four-corner average method 36 iv 1. Introduction The most fundamental aspect of the theoretical background for modern Cosmology was set in the early nineteenth century by Einstein's General Relativity. Together with facts like Olbers' Paradox and Hubble's Law, this theory led to the postulation of the Big Bang model. In simple terms, it states that our current Universe evolved from a primordial sea of photons kept in thermal equilibrium with the surrounding particles through electromag-netic interactions. This allowed the photons to remain in thermal equilibrium until the Universe was cool enough to become neutral (epoch of recombination). After that, matter and radiation evolved separately from each other (see e.g. Peebles, 1993). The Big Bang model was little more than a hypothesis before the discovery of the Cosmic Microwave Background (CMB, Penzias & Wilson, 1965), which is the remnant of the primordial pho-tons. Further measurements showed that the C M B is a nearly-isotropic radiation field with a Planckian spectrum at 2 . 7 2 5 ± 0 . 0 0 1 K (COBE's FIRAS, Mather et al. 1990, 1998). There are physical mechanisms that redshift or redistribute the C M B photons after the recombination epoch. This generally creates some small features in the overall isotropy. For example, two photons at some angular separation at the epoch of recombination will be subject to gravitational lensing by the matter distribution along their path, resulting in a different angular separation (e.g. Seljak, 1996). Or the photons might scatter through hot electron gas in a galaxy cluster, resulting in a slightly distorted planckian spectrum (Sunyaev & Zel'dovich, 1980). And evolving gravitational potential wells can redshift the C M B (Rees & Sciama, 1968). Underlying these effects is the primordial anisotropy that carries the information of the quantum fluctuations from the Big Bang -— the initial conditions that create the structures of the Universe — together with the physics of acoustic oscillations in the primordial plasma (White, Scott, &; Silk, 1994 and references therein). The anisotropy is around 10~5 K, and it was first mapped by COBE with an angular 1 resolution of around 10° (Smoot et al, 1992). Present-day technology allows us to probe the anisotropy with higher angular resolution so that we have maps of parts of the sky with a resolution near 10' (e.g. M A X I M A and B O O M E R A N G ) . In the future, MAP and Planck will extend these maps to the whole sky. In addition to the anisotropy, such technology also allows us to measure the C M B polarization, which is generally at the fiK level, or 5 to 10% of the temperature anisotropy. Polarization is generated by C M B photons Thomson-scattering off the electrons at the last scattering surface just before recombination (Rees, 1968). The Thomson scattering differential cross-section is proportional to |e • e'|, the polarization directions before and after scattering (Chandrasekhar, 1960). If a radiation field scatters off an electron, the resulting polarization of the radiation will be largely parallel to the incident one. For an isotropic radiation field, there is no net polarization because the contributions from all directions add to zero. An incident dipole radiation field similarly gives the same result because of symmetry arguments (see e.g. Seljak & Zaldarriaga, 1998, Hu & White, 1996). A quadrupolar incident radiation, on the other hand, produces a net polarization pattern because of the excess radiation intensity in the polar lobes. Since polarization is produced by Thomson scattering, most of the polarization was created just before the epoch of recombination. This is known as the last scattering surface. After that free eletrons are too scarce to have much effect. Secondary polarization could also have been produced once the Universe was reionized, but this polarization signal will be quite weak and only on the largest scales (scales corresponding approximately to the horizon at the reionization epoch). The primary constituent of the C M B polarization will be that produced during last scattering. In this way, polarization is a more direct probe than the anisotropics of the physics at last scattering. Polarization is also important because it is an unavoidable prediction of the models for structure formation— hence its detection will provide firm confirmation of the Big Bang paradigm. 2 The signals returned by a satellite to measure the anisotropy and polarization are in the form of a time series, which can be processed into a full-sky map according to the satellite's rotation pattern, or the scan strategy employed for that particular mission. We expect to obtain most of the cosmological information by turning the map into a power spectrum. In general terms, the map is a linear combination of orthogonal polynomials called spin-weighted spherical harmonics. A spin-weighted spherical harmonic of degree I characterizes features with an angular scale of around j . If a map has many features with this angular size, its power spectrum as a function of £ will be strongly peaked near that I. The shape of the power spectrum, Ce, encodes all the essential cosmological information. For example, one of the most notable features of the temperature anisotropy power spectrum is the first acoustic peak produced by the oscillations of the coupled photon-baryon fluid (see Figure 1-1). While gravity tried to compress the fluid, photon pressure resisted the compression, setting up acoustic oscillations. Once the baryons and photons decoupled, those acoustic signatures remained in the C M B as hot and cold spots (e.g. Hu, Sugiyama, & Silk, 1996). 3 -300 fjK ^m^^^^mm 300 /iK -.'loo -aoo -loo o IOO 2oo ;>oo Figure 1-1: An example of the southern microwave sky for an 1800 square degree patch taken by BOOMERANG. Here we can clearly see cold and hot spots of 1° size, which indicates that the power spectrum will be strongly peaked at the corresponding multipole I. This picture is taken from the BOOMERANG homepage (see also de Bernardis et al, 2000). o o o u o a I •a a <0 1° I -.7 10' ° Other data * Maxima-1 - BOOM 98 H-f— i * i rnr _L 200 400 600 Multipole I 800 Figure 1-2: This plot shows some power spectrum estimates from the map in Figure 1-1 together with data from the Maxima-1 experiment (Hanany et al., 2000) and a compilation of previous data by Pierpaoli et al., 2000. 4 Since various cosmological parameters define the shape of the Ct spectrum, we would ex-pect to be able to determine their values once we have completely mapped the temperature anisotropy. Things are, unfortunately, not so simple. Some combinations of parameters can produce the same features in the power spectrum that another set of combinations produces. This is commonly referred to as "parameter degeneracy". This degeneracy is breakable if we have some additional information, independent of Ct- We can use the C M B polarization power spectra just for this task. Except for its much lower amplitude, the polarization spectra contain just as much information as the anisotropy power spec-trum. In fact, polarization provides a much cleaner snapshot of the last scattering than the anisotropy because fewer physical mechanisms can alter polarization after z ~ 1000 (redshift, z, at last scattering), thus constraining the parameters more precisely. In the C M B polarization convention, polarization is separated according to its parity, ( —1)*: that with even parity (+1) is referred to as E-type and that with odd parity (—1) as the B-type (analogous with the E and B fields of electromagnetism, Seljak &; Zaldarriaga, 1998). In additionl to these two, there are also cross-correlations of E and B with the temperature, which are usually referred to as T - E and T-B cross-correlations. The T-B cross-correlation is zero because the functions that generate T and B have opposite parities (Kamionkowski et al, 1996). Thus there are four non-zero power spectra we can obtain: temperature anisotropy (T), E - and B-type polarization, and the T - E cross-correlation. Standard theories also suggest that the amplitude of the B-type spectrum may be quite small, in agreement with the current limits of tensor perturbations which could produce B-type polarization (Zibin et al., 1999). Of the three polarization-related spectra, the T - E cross-correlation will be the most easily measurable because of its higher overall amplitude. Figure 1-3 shows an example of these four power spectra (produced by C M B F A S T ) . 5 Figure 1-3: Examples of the T, E , B , and T-E power spectra produced by C M B F A S T . In a more realistic scenario, the map will not just contain the C M B alone. There will be the Galactic and extragalactic foregrounds, perhaps Zodiacal light, and of course some instrumental noise. In the case of C M B polarization, we need to worry about Galactic synchrotron radiation (Smoot, 1999), polarized dust emission (Lazarian Sz Prunet, 1999), and Faraday rotation induced by the Galactic magnetic field (Giovannini, 1995), among others. It is not an easy task to remove all these to obtain the underlying C M B . Before the age of space missions, most C M B measurements only focussed on small patches of the sky. There are some advantages to this, of course. For example, we can choose a field of view away from dusty regions of the Galaxy, and the subsequent data analysis to obtain the power spectrum will only involve a simple Fast Fourier Transform (FFT) . Staying away from the Galactic plane means we only have to deal with instrumental (and usually atmospheric) noise, which can be modelled independently, and the F F T is reasonably fast and accurate (for estimating the power spectrum). The most obvious drawback to only focussing on a small piece of the sky is that we are Hmited to probing the part of the power spectrum with upper I set by angular resolution and lower £ set by the map's dimension. Future space missions are full-sky surveys, so we can probe the whole C M B spectrum under the constraint of how well we can remove the foregrounds and how well we can calculate the C / s from the spin-weighted harmonics. The goal of this thesis is to study the numerical effects associated with the latter issue. In section 2, we will briefly summarize the mathematics for simulating the anisotropy sky and focus on the polarization and T -E cross-correlation. We will also describe how we can implement a simple H E A L P i x algorithm from an Equidistant Cylindrical Projection based on their symmetry properties. In section 3, we will discuss how we obtain the C M B spectra from the full-sky maps. In section 4, we will describe in detail the method we derived to remove the numerical errors caused by pixelization, and how well it can do so. And finally, section 5 lists some of the methods that further reduce numerical errors. o 7 2. Simulation of Temperature and Polarization Maps The goal of this section is to establish a numerical routine to calculate the C M B tem-perature anisotropy and polarization maps from some given power spectra. We will also summarize a few geometrical characteristics of the two pixelization methods we will con-sider. This routine, together with the inversion we will establish in the next section, will form a basis for our discussion of the numerical uncertainties which arise from pixelization. Making temperature or polarization maps is based on the same principle, known as or-thogonal polynomial expansion, and the polynomials in question are strongly related to each other. We will begin with the simpler of the two: making temperature anisotropy maps. The anisotropics in C M B ( A T or sometimes are scalar quantities, in other words, they remain invariant on a sphere (the sky) under rotation. Thus, naturally, we can calculate with an expansion in spherical harmonics for a given Cf. —{9A) = E E «ft»*WM). (i) 1=2 m=—l The i^rn's are the spherical harmonics defined in terms of the Legendre polynomials, Ptm-y^A) = ^ ( { ^ P U ^ e ) e ^ . (2) The aim coefficients are related to Ci (variance of a^m) by m= — t where the angular brackets denote an expectation value. In equation (1), the lower sum-mation limit for £ is 2 because we do not include the C M B blackbody temperature, £=0 and its dipole term, £=1 for the anisotropy map. The upper limit is some ^ m a x set by computational fimitations. Generally this makes the expansion an approximation, since we lose the information contained in modes above £m&x- In our case, however, the primor-dial Ci spectrum is usually damped to zero around £ ~ 3000 (except for some eccentric 8 cosmological models). Thus if our £ m a . x is high enough, the expansion will be essentially exact. A direct implementation of equation (1) makes the expansion an extremely time-consuming process. There have been a few suggestions for speeding up the computation, the most notable being the separation of variables (Muciaccia, Natoli, and Vittorio, 1997, hereafter M N V ; also in Mohlenkamp, 1997). In summary, the summations r^aax ^ ^max ^max SI - £ £ cover the same part of the £m plane; therefore, equation (1) is equivalent to A rji ^max ^max m= — lma.*l=\m\ Together with equation (2), we can factor out the only (£-dependent part, etrn^, so that A rp ^max — (M) = YI « i m **«W. (4) 771— ^max where , ,„x l2l + l(l-m)\ „ , £=\m\ V v 1 ; M N V suggested using an F F T on the eim part since 6m(#)'s are essentially coefficients of etm^. In this way, the computing time for equation (4) is reduced from JV| to iV^logTV^, and the whole routine requires computing time proportional to JV| log (usually stated as i V | since log N) = (Wlm(0)±iXlm(9))e±im^ • (6) with Wem(6) and Xim(0) defined below (Kosowsky, 1999; Kamionkowski, 1996; Zaldar-riaga, 1997): (f£—\ra\2 £(£ — l)\ cost9 \ Wtm = -\lm . 2 J + ^ 0 ' Pem(cOS0) - (£ + \ m \ ) - — P t _ 1 (cOsO) | \ \ sin 0 1 ) sin 9 J Xim = -Xim ((£ ~ l)cos.0P/m(cos0) - (£ + |m|)P*_i,m(cos d)j ; (7) _ r l2£ + l(£-\m\)\(£-2)l m " V 8 V 4TT (£+\m\)\(£ + 2)V Equations (6) and (7) say that e±im^ is the only part of the function that depends on <£, so the separation of variables that applies to the spherical harmonics also applies here. In other words, we calculate Q±iU= Y E ± 2 a ^ m ± 2 ^ m - (8) m=—tms.x(.= \m\ Since the Q and U parameters must necessarily be real (Stokes parameters are time aver-ages of physical quantities), we must be able to write equation (8) in terms of real functions only. We do this for computational simplicity and possibly efficiency on some computers. So let us define another set of coefficients in place of the ± 2 ^ m (Zaldarriaga, 1997): aE,lm — — -^{ + 2aim + - 2 ^ m ) ; 0>B,lm — — X r ( + 2 ^ m ~ - 2 < ^ m ) -Both of these coefficients can be calculated for a given cosmological model, and they are more often used in the literature because of the parity properties associated with the E 11 and B mode polarization. And although in principle we can fully expand equation (8), i.e., Q±iU= £ £ (aEflm ±aB,em){Wem(0) ± iXtm{9))e±im*, m=-lm!lyrl=\m\ we can simphfy this a little by assigning cosmcf) to m>0 and sinra^ to m<0. This is a valid operation because (a# ± aB)(Wtm ± iXimj is always real. This is the same process as changing the complex Fourier series into a series of sines and cosines. Therefore, we can write (^max 1 \ ^max "Yl cosmcf) + ^ s inm^ J ^ {aEWem ~ aBXim) m=0 m = - £ m a x / t=\m\ ( e —1 \ / ^ ' *-max ± \ *-max ^ cosrn(f) + E s inm^ J ^ (a#Xe m + aBWem), m=0 m = - £ m a x / e=\m\ where every term is explicitly real. We also tried implementing the fully expanded Q and U of equation (8). The result gives the same map but takes approximately twice the computing time of equation (9). This is because the latter contains twice as many arithmetic operations. 12 Figure 2-2: Amplitude of Q for the standard cold dark matter model in Figure 1-3. We have reduced ^ m a x to 200 for this map so that some features become more apparent. Figure 2-3: Amplitude of U produced by the same simulation for the same model at ^ m a x = 200. 13 Figures 2-2 and 2-3 show polarization in the form of the two Stokes parameters. A n alternative way is to show maps of polarization amplitude, P, and angle, a, which are defined as (e.g. Kosowsky, 1996) P= ^/Qi + U2, and 1 -iU a = - tan —. 2 Q We do not show these maps here because both contain the same amount of information, and our eyes will not see any significant difference between the two sets. The last of the four maps relevant in the analysis of C M B data is the T - E cross-correlation map. In the case of Gaussian statistics (what we have assumed so far for all the temperature and polarization maps), making realizations of temperature and polarization is easy since and for any Gaussian random number r\ and r2 with unit variance. However, the presence of a T - E cross-correlation means that r\ and r2 are constrained such that / T* E \ \atmalm) — W where " X " denotes cross-correlation. Thus, independently-chosen T and E realizations such as those in Figures 2-1, 2-2, and 2-3 will not give a T - E map predicted by the same cosmological model. Since cross-correlation is a scalar quantity, making a T - E map requires only the spherical harmonics used for the anisotropy map. We have ignored the cross-correlation here, since there is no effect on numerical error analysis. 14 Up to this point we have only discussed the separation of variables method, also known as the summation exchange, into the routine. There are also other symmetries in the Legendre polynomials that we can exploit to save additional computing time. One of them is a symmetry in m (Abramowitz & Stegun, 1972): Pe,+m = P^-m- (10) This saves quite a significant amount of time (close to a factor of two), since generating the Legendre polynomials is the slowest part of the algorithm. Another symmetry is ( M N V , 1997; Turok & Crittenden, 1998) P (n™ d\ - I +Pim(cos(Tr - 9)), if £ + m is even; n ^, ^m{COStf)~\-Pemlcos(7r-9)), if £ + m is odd. This means we only have to calculate the Legendre polynomials above the equator. Those below just differ by a minus sign for odd £ + m. These symmetries can be directly imple-mented into equation (2). But because of the more complicated 6 and P ^ m dependence of Wim and Xim (equation (7)), this gives us another exercise in keeping track of minus signs. 2.2 Recurrence Relation of Legendre Polynomials Now we come to the question of calculating the Legendre polynomials. This is usually done using recurrence relations. There are several of these recurrence relations available, but most of them involving m contain a factor of (1—a:2) -^ with x= cos 9, which is numerically unfavourable because of the singularity (Press et al, Chapter 6.8). For our calculations, we adopt the relation Ym-(2£-l)xYm / ( l + ™ - l ) ( l - ^ - l ) r , n i Y^ = V2£^rlxY^1. 15 This recurrence relation is the one mentioned in Zaldarriaga (1997) and seems to be the most accurate and fastest subroutine available (downloadable at the C M B F A S T home-page). In section 4.1, we will show that the recurrence relation is stable when we discuss numerical uncertainties. For the spin-weighted harmonics, we simply evaluate equation (7) for each £ and m. How-ever, there are places where sin 2 0 in the denominator is near zero. We can avoid this by assigning the functions' values at these points. Specifically, H m Ptm{cos0) = H m - s m f l P p o s t f ) = _1 0-+o sin 2 0 «-o 2 sin 0 cos 0 2 £ m V ; P'lrrSX) * s O I u y n o n - z e r o f ° r m — 2 because / j \ \m\ / j \ |ra| + l and the (1 — x2) that makes P'im(x) zero for x = 1 disappears from the first term when \m\ = 2. So U m Pim(cos9) = f - l ^ - i ) ^ + i)(^ + 2), |m| = 2; o-+o sin 2 0 \ 0, all others. If we apply this to equations (7), we find, for all £, the values of Wn{9) and X£2(9) where sin0 is zero: WE2(9) sin 0=0 V 47T So how important is this? In Figure 2-4, we plotted A (sin 0== 10 1 0 ) which is defined as V M s i n 0 = l ( n 1 0 ) - W £ a 2 n a l y t i c (s in0=O). We also plotted A W « at two other values of sin 9 for reasons that will become apparent shortly. Since 1 0 — 1 0 is quite small, naively one would explain the difference as the numerical 16 inaccuracy as the result of dividing by something close to zero. This is why we also plotted the other two curves in Figure 2-4: as sin# approaches zero, AW 12 also approaches zero asymptotically. So what we see in A W ^ s i n # = 1 0 - 1 0 ) is not a numerical error, but a real difference. We can further justify this statement by the hmits of double-precision floating points — 1.7E ± 308. So 10~ 1 0 is well within this range. The analytic formula is useful when we need to calculate the value of the pixels precisely at the poles (sin 6 = 0 exactly). And there are times when this is necessary. 0 .010 0 .005 CM 0.000 - 0 . 0 0 5 - 0 . 0 1 0 500 1000 1500 2000 2 5 0 0 3 0 0 0 M u l t i p o l e I Figure 2-4: A plot showing Wt2 at various values of sin#. This also demonstrates the recurrence relation's ability to accurately generate Wtm and Xlm at very small angles 2.3 Pixelization Throughout this section we have always been using the separation of variables method. This is because for large simulations, the direct implementation becomes too slow to be useful. We have also seen that for this summation exchange to work properly, the pix-elization scheme has to contain equal-latitude rings — rings at constant (j) (Crittenden & Turok, 1998). This is an essential feature for fast spherical harmonic transforms. There 17 are many pixelization schemes that have such a feature, and in this thesis we will focus on two of them: the Equidistant Cylindrical Projection (ECP), because it is simple and obviously flawed, and the HEALPix pixelization, because it has been commonly adapted as the standard spherical harmonics analysis package for its numerous advantages. Other pixelization schemes that have been proposed include the quadrilateralized spherical cube (Chan & O'Neill, 1975), which was used in the COBE data analysis; the icosahedron-based method (Tegmark, 1996), which has a certain degree of spherical symmetry (not just azimuthal); and IGLOO pixelization (Crittenden & Turok, 1998), which is rather sim-ilar to HEALPix in terms of geometrical properties. These share some, but not all of the advantages of HEALPix . To implement this routine into an E C P pixelization is a trivial matter because ^jF-(6,(j>) will simply form a rectangular array. The 0 coordinate is incremented evenly from 0 to 7r and the (j) coordinate evenly from 0 to 2TT. When we project this array onto a sphere, we see that the pixel concentration is much higher near polar areas than near the equator. This introduces a few problems: pixel shape distortion, unequal pixel area, and unnecessarily-wasted computating time. In section 4, we will focus on these problems directly and find methods to remedy them as much as possible. Of the three problems we stated for the E C P pixelization, the HEALPix pixelization avoids the latter two while having only minimal pixel distortion (Gorski, 2000). At base resolution, it has 12 equal-area pixels: 4 around the north pole, 4 along the equator, and 4 around the south pole. Each step up in resolution divides each pixel into 4 equal-area pixels. The hierarchy allows for an efficient way to locate neighboring pixels, which is essential if a quick analysis demands only a lower resolution map. A more detailed description of H E A L P i x can be found at the HEALPix home page at http://www.eso.org/~kgorski/healpix/. The webpage provides a numerical implementation of the HEALPix pixelization scheme in Fortran 90. Here we will show that we can modify the E C P code only sightly so that it 18 can produce maps in the H E A L P i x scheme (i.e. we use the H E A L P i x algorithm in our code). H E A L P i x is a curvilinear partition of a sphere into equal-area quadrilaterals with varying shapes. The E C P is a simple-minded pixelization that happens to have a few nice geomet-rical properties which are essential for fast spherical harmonic transform and for exploiting the various symmetries of the spherical harmonics. On first observation, the pixel shapes look very different for the two. However, we must note that unless we perform an exact integration in the inversion (see section 3), all that the program "sees" are the pixel cen-ters, whose values are each taken as an approximation for the entire pixel. So in both pixelizations, we have points that He on rings of constant latitude. The only difference is that H E A L P i x tries to conserve pixel area by having fewer and fewer points near polar caps (see Figure 2-5). We can integrate this fact into the E C P by calculating only the necessary ezm^ around the poles. For example, Figure 2-5 shows the two pixelizations at the same resolution: 22.5°, with H E A L P i x at its first partition (each base pixel divided into four) and the E C P being an 8 x 16 rectangular array (8 divisions along equator). In the context of modifying the E C P , all we need to do is evaluate eim^ at the four locations in the first H E A L P i x ring, the eight in the second ring, and so on, instead of the same number of eim<^ for every latitude (ring). This is only a matter of coding. In addition, the e , m*'s need not be calculated whenever we start with a different latitude. The angles are the same for the northern and southern hemispheres, and we can pre-compute all the e im^"s as soon as we decide what resolution we choose. So all the symmetry properties are applicable to this pixelization. This modification results in a slightly faster and more accurate algorithm because of the fewer number of polar pixels. The issue of speed depends very much on the computing facility the routines are performed, but generally the H E A L P i x routine requires about 10% less time than the E C P routine. The speed is not so dramatically different because we 19 still have to calculate the same number of ^-intervals, which are the most time-consuming. The most notable improvement is in the accuracy of inversion. We performed an inversion for some maps using the resolution seen in Figure 2-5, the resulting residual (difference between obtained spectrum and input) for the HEALPix routine was 100 times smaller then that for the ECP. Thus, if we can reduce the numerical error found in the E C P routine substantially, the numerical error in the HEALPix routine will become even smaller. Figure 2-5: This figure shows that all we need to do is modify the -summation a little in order to implement the HEALPix pixelization. The one to the left is HEALPix, and to the right is ECP. 20 3. Inverse Spin-Two Spherical Harmonics Transform Now we turn our attention to an equally important issue — the inversion to the power spectrum. C M B anisotropy data have been available for almost a decade, so not sur-prisingly, the spherical harmonics inversion has been investigated several times, the most notable being the Quadrilateralized Spherical Cube of COBE, the fast E C P inversion uti-lizing F F T , and H E A L P i x . In this section we will mainly focus on the inversion of C M B polarization data, in other words, inversion involving the spin-2 spherical harmonics. This is an otherwise trivial matter if not for the various sign conventions seen throughout the literature. In section 2 we have written down an explicitly real expression for calculating the Q and U Stokes parameters by defining a set of coefficients from ± 2 & £ m that can be predicted by various cosmological models. Our goal now is to do the same for the CLE and ag coefficients. The functions ±2 i ^ m ' s are mutually orthonormal on a spherical surface, ±2Y(*m ±lYvm' Sm0 dO d(f) = Sa1 &mm'• (13) We apply this to equation (8) and solve for the coefficients ±2dtm'-±2aim = J J(Q±iU)(W ±iX)* (e±im(t>)* sin0 d6 d. For numerical calculations, we must convert the integral into a summation with a certain "grid size" (resolution) and consider 1 up to some limit of ^ m a x - In general, if we are inverting a numerically simulated map, the £ m a x will be set by the ^ m a x of the map since the simulated map itself contains no information above that. For the inversion of an actual data set, ^ m a x is limited only by the angular resolution of the experiment. In complex notation, we perform ±2) AO. (14) j=i ' 21 The (W ±iX)*j and eim<^ are usually approximated by their values at the pixel center, but at times we might need to integrate over the entire pixel or perform a better approximation instead of just using its central value. This is the only approximation necessary in the mapping and inversion processes. Equations (1) and (8) are still exact whether ^ m a x is finite or infinite. After some algebra we obtain aE and aB from equation (13). And as before, we choose to assign cosm to ra>0 and smrruj) to m<0 instead of fully expanding equation (14): i = l i=l Ne Nf aB=Y, sin 9j ^{QijXj + UHXi) { tinrnt m < u! i=i i=i The expressions above are written in variable-separated form in which the (^-summation is outside of the ^-summation. So again we can implement this with the E C P and HEALPix pixelizations. 22 4. Numerical Errors The purpose of all these expansion and inversion routines is that we need to know how much numerical error is contained in the obtained spectra, T, E , B , and T - E , in the absence of experimental (instrumental) errors. In the hypothetical case where the inversion routine leaves absolutely no numerical error, the only part we need to worry about is the recurrence relation, which gives errors independent of the pixelization scheme. So for all purposes, if we can reduce the overall numerical error to a magnitude comparable to that made by the recurrence relation, our work is done. That is why we test the recurrence relation first in this section. 4.1 Stability of the Recurrence Relation of Legendre Polynomials We can test the stability of the P^ m ' s recurrence relation by looking at their values at certain 6 where analytic expressions exist. For example, when cosf? is 0 (Abramowitz & Stegun, 1972), om _ TV L _ l _ m 4 . I") ^ ( 0 ) = ^ c o 8 ( ( / + m ) | ) r ^ ^ ; ) J . • (16) Since we start the recurrence relation (equation (12)) at Pu: the error will be the largest for PtQ. Figure 4-1 is a plot of Peo(0) for even £ up to 3000 which we generated, overlaid with the values given by equation (16). Those for odd £ are zero due to the cosine term. We intentionally add 0.01 to those produced by equation (16) because the two differ only by a small amount. (We do not want the reader to have the impression that we forgot to plot one of them.) Figure 4-2 is the difference of the two. The increasing error is expected because each Pio(0) is the result of 2£+l recursions. Other than this fact, the errors are small and randomly distributed near 10~ 7. The spin-2 harmonics are computed from the Legendre polynomials, so the uncertainty is of the same order of magnitude. For example, the Wn and Xn generated recursively and through analytical expressions (equations (7) and (16)) behave in the same way as the Pio we have discussed. We do not want to fill the 23 page with identical figures, so they are shrunk and put into Figure 4-3. We chose m=l because Xio=0 for all m, which cannot tell us anything about errors. The error in Xn is much smaller because the non-zero term in Xt\ is a factor of £ less than the non-zero term in Wu (refer back to equation (7)). The routine Gammln(x) of Numerical Recipes (Press et al, Chapter 6.1) we used to calculate Gamma functions is based on an approximation derived by Lanczos in 1964, which looks like r ( z + 1 ) = + T + 0 * + 2 c-<*+^*> (co + J2 7 ^ ) for certain integer 7 and N and coefficients C{. The routine we chose includes up to C§ and gives errors less than 2 X 1 0 - 1 0 , which is sufficient for our test where errors are near I O - 7 . In Figures 4-2 and 4-3, there are some features around £ ~ 1400, 2100, 2700...etc. We believe this pattern appears in the recurrence relation since Press et al. claims that the accuracy of Gammln(x) is within 2 x 10~ 1 0 . This pattern does not affect the overall accuracy in any way because the magnitude is so low. cu 0.30 0.20 g-0.10 ^ - 0 . 0 0 P-- 0 . 1 0 - 0 . 2 0 F; - 0 . 3 0 0 500 1000 1500 2000 2500 3 0 0 0 M u l t i p o l e I Figure 4-1: P^o(0)'s for even I up to 3000. The blue, dashed, lines are from equation (16) shifted vertically by 0.01, and the black lines are from the recurrence relation. 24 o X OH -a . i—i m CD 1.0 0.5 0.0 -0.5 h 1.0 •• : A > ' : : $ : *•-••V." 500 1000 1500 2000 2500 3000 Multipole I (even) Figure 4-2: The difference between equation (15) and the P^o(0)'s generated by the recurrence relation (without that constant we added). ca a "tn 0.10 0.05 t" 0.00 -0.05 b -0.10 0 500 1000 1500 2000 M u l t i p o l e I (odd) 2500 3000 o x CO T 3 V) CP 500 1000 1500 2000 M u l t i p o l e I (even) 2500 3000 Figure 4-3: The same plot as Figure 4-2 but for the functions Wtl and Xtl. 25 4.2 Orthogonality Issues Related to Pixelization Since the errors from the recurrence relation define a limit on how accurate the inver-sion can be, we can set our goal to reducing the other numerical errors to a comparable level. One of them comes from the fact that we only consider up to a finite £ m a , x of the orthogonal polynomials. This issue is described in detail in Crittenden & Turok (1998) for spherical harmonics. In summary, by imposing a pixelization scheme on the sky, we break the orthogonality of spherical harmonics. A degree-^ spherical harmonic becomes an approximate linear combination of the rest of the spherical harmonics of the set {0,^m a x}. And in the process of inversion, the coefficient of each mode picks up something extra from this "linear combination". This effect on spherical harmonics can be efficiently reduced by fine resolution, since generating the spherical harmonics and summing over all ^-modes have become a relatively fast procedure. After this all we see are the errors from the recurrence relation. Thus, at high enough resolution, the orthogonality issue for spherical harmonics is of little concern. To correct for this effect at low resolution, we can apply the unbiased Ct spectrum estimator shown in Crittenden & Turok (1998), which depends only on the pixelization scheme: estimated Ct = Mtt'Cf, where M i i ' = 27TlS (jlAPWp,tmWp,tm) E ^ ^ m ^ m ' E A Q W Q ^ ^ Q ^ ' m ' . mm' V P J P Q Letter A denotes pixel area, P and Q are pixel indices, and the Wpixeit£m here is the window function of the pixel (-^ within the pixel and 0 elsewhere) integrated with Ytm over the entire pixel. This removes the majority of the leakage between multipoles, but it ignores the off-diagonal correlations between estimated atm- However, because of this fact, calculating Mtt> is relatively fast (Crittenden &: Turok, 1998). 26 The above issue could have been called something along the lines of "orthogonaUty-breaking as a result of pixelization", but the word "leakage" has become the description of such an effect. We can expect to see this same kind of error in the Q and U maps. Equations (9) and (15) say that and as together determine Q and U, which in turn are required to calculate ag and a#, so there are effectively twice as many multipoles to leak into each other. Combined with the fact that both Wim and Xim depend on Ptm and J R g _ i ) m , it is not surprising that the leakage is as much as that seen in Figure 4-4. In reality, the amplitudes of the B-mode spectrum is likely to be small (or zero), which greatly reduces the leakage in the E-mode. We exaggerated Figure 4-4 so that we can obtain a worst-case estimate of the errors. This is an important consideration particularly for trying to measure any weak B-mode signal in the presence of leakage from the E-mode. 10 8 6 u 4 2 0 10 100 1000 Multipole I 10 8 »- 6 10 100 1000 Multipole I Figure 4-4: A plot on the leakage problem for the polarization spectra. The blue lines are the spectra as the direct result of inversion. The red, dashed lines are the spectra we used to calculate polarization. The normalization for both spectra is arbitrary. Now we will present a method to remove this leakage. Let us suppose we invert some 27 polarization map (data) which gives us overestimated E and B spectra. Since 1 e 9 m=-l the 's are more useful when we want to know how much extra contribution they get from other modes of ±2Ytrn. And from here on we will use their complex notations, i.e., equation (8) for the expansion, Q±iU = ^ Y2 ± 2 ^ m ± 2 ^ m , (17) e—o m=-t and a simpler form of equation (14) for the inversion, ±2aim= £ (Q±iU)±2Ye*m. (18) a l l p i x e l s We will also assume that the expansion is exact based on the reasons discussed earlier in section 2 (in short, due to that fact that the spectra are zero after some £). Thus the only significant numerical error here other than that made by recurrence relations is due to pixelization. For a particular pixelization, what equation (18) does is introduce an extra additive constant for each £ and m. Therefore, ±2imi which contains the true power spectrum and a little offset, and calculate Q ± iU using equation (17): ^max (Q±lU)1 = ^ ^ ± 2 ^ m ± 2 ^ m 1=0 m=-t ^max ^ / \ = E E Ua?-e + A ± 2 a £ r ) ± 2 ^ m (20) = 0 m=-e 28 Then we find that we can estimate A±2<^m by another inversion: True E A(Q±iU)±2Y* all pixels (21) where A(Q ± ill) is (Q ± iU)1—(Q ± iU). We say approximately equal in equation (21) because A ± 2 ^ m is also subject to the same pixelization effects as ±2«£m- This gives us a first-order approximation to ±2aJmae- Examples of this are in Figures 4-5 and 4-6. 10 100 Mult ipole Multipole I 1000 1000 Figure 4-5: This shows the first-order approximation to ±2aLu e in blue. The purple line is the ±2=0, equation (14)), the summation becomes an integral. But how much can we increase the resolution before the increased computing time outweighs the gain in accuracy? In reality we can not increase the resolution indefinitely since each experiment has a beam size, which experimentally supresses information. And there are also computational limitations. But on a purely theoretical basis, we can ask what is a good compromise between accuracy and computing time, which is directly related to resolution. To answer this question, we performed a test with a low lmeLX map to see how the accuracy improves with resolution. For the E C P at a fixed £ m a . x , the computing time generally increases linearly with the number of pixels. In Figure 4-4, we plotted the variance of the residual spectrum (difference of input and inversion output) versus sampling factor, where a sampling factor of 1 means we sample at the Nyquist frequency given by ^m a x> and higher means we oversample the map by that factor times the Nyquist frequency. We see that sampling twice as often as the Nyquist factor gives a reasonably small numerical error while keeping the computational time short. This is precisely what is expected for MAP and Planck, which are planned to over-sample the beam size by about a factor of three. 34 CD O C 1.10 cti £ 1.08 1 1.06 -a * i—i W v 1.04 u 1.02 N • I 1 "cd fi i - o o u o ~i—I—i—i—i—i—r i i — i — i — _ i i i i _ 1.0 1.5 2 .0 2 .5 3 .0 3 .5 4 . 0 S a m p l i n g F a c t o r ( x N y q u i s t f r e q u e n c y ) Figure 5-1: A plot showing how the numerical error decreases with increasing number of pixels. The vertical scale is normalized to the point on the far right. Computing time scales with the square of the Nyquist frequency if ^ m a x is fixed. 5.2 Approximation of Inverse Transform Integration Basically, everything we described in section 4 can only be done after we have obtained the power spetra. We assumed that the pixel-center approximation for the integral in equation (14) was sufficiently accurate, but in most cases, it is not. This is especially true for low-resolution maps where the value of ±2Ytm,s can vary greatly within a pixel (longitudinally). In this case, the pixel-center values are no longer sufficient, and we must look for better ways to calculate ±2 a i m = J J{Q±iU)(W ±iX)* (e±im)* sind d$ d, which, in its variable-separated form, looks like = / s\nQ(W±iXy f (Q±iU) (e±irn)* d d8. (28) (29) 35 One option is to consider averaging the ±$Yim from the four corners of each pixel (Figure 5-1). This is advantageous because it is a better representation of ±^Yim for each pixel, and it defines a pixel "boundary" which the computer does not know otherwise. However, to do this, we must resort to an inversion without separation of variables (equation (28)), since summing over the four e ± i m < ^ values and then multiplying by the four (W ± iX)* will produce twelve unwanted and unremovable cross terms. This fact alone prohibits the use of four-corner averaging. Any other similar n-point averaging technique will have the same problem as long as it has points lying on different latitudes. Figure 5-2: This illustrates the four-corner averaging scheme. Unfortunately, its incompatibility with separation of variables renders it useless. Another option is to only average points on the same latitude (points with subscripts 2 and 4 in Figure 5-1). This works in unison with equation (29)'s variable separation, and since it does not require any additional calculations of (W ± iX), the algorithm does not suffer any significant loss in speed. We performed a low-resolution simulation with £ = 100, and it showed that this averaging can reduce the overall error (sum of errors in each multipole) by 25%. This is not much, considering it is only a reduction in overall error (the improvement in each individual multipole is approximately 2.5 x 10 - 3 ) , but it 36 is still better than nothing. The time for the small simulation was about 20% longer than a simulation without averaging. 5.3 E C P with Equal-Area Pixels We know the E C P pixelization loses its accuracy partly because of varying pixel areas. Once the rectangular pixels on a plane are projected onto a spherical surface, each pixel's area is altered by a sin# factor. In an attempt to correct for this, we tried to calculate 0's at intervals such that each pixel has the same area, regardless of shape. Apparently, this increases pixel shape distortion because now we have very-long polar pixels. We found that this neither improves the inversion accuracy nor computing time. Although this seems like a worthless study, we were able to conclude that both equal shape and area are necessary constituents of a good pixelization scheme, as previous studies have shown (e.g. Crittenden & Turok, 1998). The majority of the pixels in the original E C P pixelization have similar shapes. The better result given means that pixel shape is a more important effect than pixel area. This is reasonable because pixel area effects can be removed by proper normalization, while the normalization for shape is more difficult to quantify. In practice we can improve results of E C P (for the same total number of pixels) by re-ducing the number of polar pixels and increasing the number around the equator. This can be done without any unnecessary computing time because the number of operations involving ± 2 ^ m remains the same. Continuing this process, we eventually would arrive at a pixelization scheme very similar to HEALPix or IGLOO pixelization. 37 6. Conclusion The main numerical error in the calculation of C M B temperature and polarization power spectra comes from pixelization schemes which cannot fully describe spin-zero and spin-two spherical harmonics. We found that such numerical error can be removed to within 1 0 ~ 6 of the original spectra by an iterative method that calculates how much leakage one multipole receives due to a particular pixelization. We also investigated some other numerical issues that might be relevant for future polarization studies. Together with previous studies on some similar issues for temperature anisotropy and polarization, our results ensure that the numerical part of the polarization data analysis should not contribute any significant errors to the determination of the cosmological parameters. 38 References Abramowitz, M . and Stegun, I., 1972, Handbook of Mathematical Functions, Dover New York. Chan, F. K . and O'Neill, E. M . , 1975, Feasibility Study of a Quadrilateralized Spherical Cube Earth Database. Computer Sciences Corporation. Chandrasekhar, S., 1960, Radiative Transfer, Chapter 1. Dover, New York. Crittenden, R., astro-ph/9811273 Crittenden, R. and Turok, N . G., astro-ph/9806374 de Bernardis, P. et al., 1999, American Astronomical Society Meeting 195. Doroshkevic, A. G. et al, astro-ph/9901399 Giovannini, M . et al., 1995, Physical Review Letters, 75, 3796. Gorski, K. , 2000, The HEALPix Primer, HEALpix Homepage. Gradshteyn, I. and Ryznik, I., 1994, Table of Integrals, Series, and Products. Academic Press. Hanany, S., 2000, in Press. Hu, W., Sugiyama, N . , and Silk, J., Nature, 1996. Hu, W. and White, M . , 1997a, Astronomy and astro-physics, 321, 8. Hu, W. and White, M . , 1997b, New Astronomy, 2, 323. Hu, W. and White, M . , 1997c, Physical Review, D, 56, 596. JafFe, A. H., Kamionkowski, M . , and Wang, L. , astro-ph/9909281. Kamionkowski, M . , et al, 1997, Physical Review D, 55, 7368-7388. Kosowsky, A. , 1996, New Astronomy Reviews, 43, 157-168. Lanczos, C , 1964, SIAM Journal on Numerical Analysis, ser.B, 1, 86-96. Lazarian, A. and Draine, B., 1997, Astrophysical Journal, 512, 740-754. Lazarian, A . and Prunet S., 1999, Microwave Foregrounds, ASP Conference Series #181, 113, San Francisco. Mather, J . et al., 1990, Astrophysical Journal Letters, 354, L37. Mather, J. et al., 1998, Astrophysical Journal Letters, 508, L25. Mohlenkamp, M . J., 1997, A Fast Transform for Spherical Harmonics, Yale Ph.D. Thesis. 39 Muciaccia, P. Natoli, P. and Vittorio, N . , 1997, astro-physical Journal Letters, 488, L63. Newman, E. T., et al, 1967, Spin-s Spherical Harmonics and d. Journal of Mathematical Physics, 8, 11. Newman, E. T., Penrose, R. , 1966, Note on the Bondi-Metzner-Sachs Group, Journal of Mathe-matical Physics. 7, 5. Peebles, P. J . E., 1993, Physical Cosmology, Princeton University Press, Princeton. Penzias, A . and Wilson, R., 1965, Astrophysical Journal, 142, 419. Pierpaoli, E., Scott, D. and White, M . , 2000, Science. Press, et al., 1993, Numerical Recipes, Chapter 6, 2nd edition, Cambridge University Press, London. Rees, M . J . and Sciama, D. N . , 1968, Nature, 217, 511. Seljak, U . , 1996, astro-physical Journal, 436, 509-516. Seljak, U . , 1997, astro-physical Journal, 482, 6. Seljak, U . and Zaldarriaga, M . , 1998, 19th Texas Symposium on Relativistic Astrophysics and Cosmology, 14, Paris. Smoot, G. F. et al, 1992, Astrophysical Journal L, 396, L13. Smoot, G. F., 1999, Microwave Foregrounds, ASP Conference Series #181, 61, San Francisco. Sunyaev, R. and Zel'dovich, Y . B. , 1980, Annual Review of Astronomy and astro-physics, 18, 537. Tegmark, M . , 1996, Astrophysical Journal L, 470, L81. White, M . , Scott, D., and Silk, J., 1994, Annual Review of Astronomy and astro-physics, 32, 319-370. White, M . , 1998, Evolution of Large Scale Structure. Proceedings of the MPA-ESO Cosmology Conference, Garching, Germany, 25. Zaldarriaga, M . , 1998, Astrophysical Journal, 503, 1. Zibin, J . Scott, D, and White, M . , 1999, Physical Review, D, 60. BOOMERANG: http://www.physics.ucsb.edu/~boomerang/ CMBFAST: http://www.sns.ias.edu/~matiasz/CMBFAST/cmbfast.html HEALPix: http://www.eso.org/~kgorski/healpix/ M A X I M A : http://cfpa.berkeley.edu/group/cmb/ 40