NEW APPROACHES TO MAP-MAKING FOR COSMIC MICROWAVE BACKGROUND OBSERVATIONS BY CHARMAINE ROSE ARMITAGE-CAPLAN B.Sc, University of British Columbia, 2003 M.S., University of Illinois at Urbana-Champaign, 2005 M.S., University of Illinois at Urbana-Champaign, 2007 DISSERTATION Submitted in partial fulfillment of the requirements for the degree of Doctor of Philosophy in Physics in the Graduate College of the University of Illinois at Urbana-Champaign, 2009 Urbana, Illinois Doctoral Committee: Associate Professor Susan Lamb, Chair Associate Professor Benjamin Wandelt, Director of Research Associate Professor Robert Brunner Professor Aida El-Khadra UMI Number: 3391874 All rights reserved INFORMATION TO ALL USERS The quality of this reproduction is dependent upon the quality of the copy submitted. In the unlikely event that the author did not send a complete manuscript and there are missing pages, these will be noted. Also, if material had to be removed, a note will indicate the deletion. UMT Dissertation Publishing UMI 3391874 Copyright 2010 by ProQuest LLC. All rights reserved. This edition of the work is protected against unauthorized copying under Title 17, United States Code. ProQuest LLC 789 East Eisenhower Parkway P.O. Box 1346 Ann Arbor, Ml 48106-1346 © 2009 by Charmaine Rose Armitage-Caplan. All rights reserved. Abstract In this thesis, I examine systematic effects due to beams and foregrounds that will effect precision temperature and polarization measurements of the Cosmic Microwave Background (CMB) by satellite instruments. I have developed a unique, pixel-free solution to the general CMB map-making problem, which addresses the impact of beam asymmetries through beam deconvolution. The resulting algorithm, called PReBeaM, solves for the map directly in spherical harmonic space, avoiding pixelization artifacts. The most promising experimental means for testing the widely popular inflationary model of the early Universe is through the detection of primordial gravitational waves (or tensor perturbations) and the resulting imprint of B-mode polarization. As we enter the PLANCK era and become more sensitive to the CMB polarization signal on both large and small scales, the possible detection of inflationary B-modes seems achievable. Detection of this putative, minute signal will most certainly depend on our understanding of, and control over, sources of uncertainty due to systematic effects and foregrounds. My PhD research has been focused on developing a novel CMB map-making technique for removing systematics as a result of asymmetries in the beam. I demonstrate that my deconvolution method produces maps and spectra with less contamination when compared to methods that do not account for beam asymmetries in both temperature and polarization. I show that the true sky is recovered with greatly enhanced accuracy via the deconvolution method for temperature measurements in realistic experiments with asymmetric beams and far-side lobe Galactic foreground contamination. I also illustrate that standard map-making is unable to remove artifacts due to optical systematics in these cases. I then extend my method to operate on polarization data and apply my technique to the actual PLANCK 30 GHz simulated CMB data. In an actual CMB experiment, we are faced not only with challenges due to beam asymmetries, but also noise and astrophysical foregrounds. I study the noise properties of CMB maps from PReBeaM for PLANCK and include important foreground signals from dust, synchrotron and free-free emission. ii Extragalactic point sources are an additional foreground contaminant in CMB maps; proper identification and removal of these sources is a crucial step in the analysis pipeline. Prom the standpoint of Extragalactic Astrophysics, point source data from the WMAP instrument can provide important information about source variability on a range of different time scales. I establish a program to search the public time-ordered WMAP datasets for under-exploited multi-band data on the variability and evolutionary properties of point sources catalogued by the WMAP team. hi For David IV Acknowledgments There are many people deserving of acknowledgement for helping to bring me to this point in my academic career. I owe my advisor, Ben Wandelt, many thanks and my sincerest appreciation for guiding me so well through my PhD research. Many times, I have been inspired by Ben's excitement about CMB research and benefited from his breadth of knowledge. I credit my two summers spent doing undergraduate research at the Canadian Institute for Theoretical Astrophysics for starting me on the path of Cosmology and CMB research. It was there I was supervised by Dick Bond and Carlo Contaldi, both of whom were excellent mentors for someone who knew nothing about CMB power spectra or Fortran programming. During my time at Illinois, I have had the pleasure to work with many friendly and bright colleagues and I have received help from countless people. I thank David Larson, Chad Fendt, Amit Yadav, Rahul Biswas, and Rishi Khatri for helpful discussions and for sharing Ben's time and attention with me. I am grateful to many people at the NASA Jet Propulsion Lab who worked with me on the PLANCK collaboration. They include, but are not limited to, Charles Lawrence, Kris Gorski, Gary Prezeau, Jeff Jewell, and Graca Rocha. I am also thankful to Mark Ashdown at Cambridge for early conversations on map-making. I thank the team at the Lawrence Berkley Lab — Julian Borrill, Chris Canteloupo, and Ted Kisner — for being experts at computing and helping me with many of the issues I had with analysing the simulated PLANCK data. I was very fortunate to have been able to spend three weeks in June of 2007 at the Max Planck Institute for Astrophysics. I thank MPA for their generous hospitality. While I was there, and in many instances before and after, I received invaluable help and advice from Martin Reinecke on all things convolution and interpolation related. Without this help, it surely would have taken twice as long to debug my code. I would like to acknowledge the following sources of funding for this research: the National Center for Supercomputing Applications (NCSA), the Center for Advanced v Studies, NASA contract JPL1236748 by the National Computational Science Alliance under AST300029N and the University of Illinois, NASA contract JPL1371158, and NSF grant AST-05-07676. Additionally, I used resources at the National Energy Research Scientific Computing Center (NERSC), which is supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. Some of the results from Chapter 3 were computed on the NCSA computer Copper. Many of the results contained in my thesis, and particularly those in Chapter 4, were calculated on the NERSC computers Bassi and Seaborg. I thank the NERSC help consultants for being so prompt and thorough with all of my help requests. Additionally, I am indebted to the entire team who worked to produce the PLANCK simulated data on which I ran my code to produce the results in Chapter 4. Some of the work that I did in Chapter 5 was completed in collaboration with Torsti Poutanen at the University of Helsinki, who was always willing to answer my questions with prompt replies. None of my work on the WMAP variable point sources in Chapter 6 would have been possible without vital help from some members of the WMAP team: Gary Hinshaw, Nils Odegard, and Bob Hill. I also acknowledge extensive use of the Legacy Archive for Microwave Background Data Analysis (LAMBDA) site for WMAP timeordered data files, beam profiles, masks, and the point source catalogue. I thank my committee members: Susan Lamb, Robert Brunner and Aida ElKhadra for their interesting questions and thoughtful comments. I am indebted to Wendy Wimmer for answering dozens of my questions about dissertating and defending, for quelling my anxiety in the last weeks of my studenthood, and for striving to make my experience at Illinois as smooth and painless as possible. I thank Xiaoyue Guan for being an excellent office-mate, and for many funny conversations during long work days. I thank fellow physicists, Keiko Kircher and Madalina Colci-O'Hara for being good friends and for sharing the unique experience of entering motherhood during graduate school. It is Caplan. frequent much to without a doubt that I owe my deepest gratitude to my husband, David I can't imagine how I would have done this without his love, support, and pep talks. Finally, I thank my son, Emmanuel Caplan, for giving me so smile and laugh about every single day for this last year and a half. VI Table of Contents List of Tables x List of Figures xi Chapter 1 Introduction 1 Chapter 2 CMB Theory 2.1 Introduction 2.2 The Standard Cosmological Model 2.3 Models for the formation of perturbations 2.3.1 Fourier analysis of density fluctuations 2.3.2 Inflationary models 2.3.3 The cyclic model 2.4 Inflationary gravitational waves 2.5 Formation of Cosmic Microwave Background Anisotropics 2.6 Theory of CMB polarization 2.7 Review of Stokes parameters 2.8 Harmonic analysis for temperature and polarization anisotropics . . . 2.9 Conclusions 5 5 5 7 7 8 10 10 11 12 15 15 16 Chapter 3 Deconvolution Map-Making 3.1 Summary of chapter 3.2 Introduction 3.3 The map-making problem 3.3.1 Maximum-likelihood estimator of the true sky 3.4 Standard map-making methods 3.5 Deconvolution Map-Making 3.5.1 Regularization technique 3.5.2 Fast convolution on the sphere 3.5.3 Binning and simulation 3.6 Solving the Deconvolution Equations 3.6.1 Preconditioning 3.7 Standard map-maker for comparison 3.8 Test Cases 3.9 Results and Discussion 18 18 18 19 21 22 23 24 25 26 27 28 29 30 31 vii 3.10 Point source tests with elliptical beams 3.11 Conclusions 37 37 Chapter 4 Polarized map-making for Planck 4.1 Summary of Chapter 4.2 Introduction 41 41 41 4.3 PLANCK 42 4.3.1 The map-making challenge for PLANCK PReBeaM Method Comparison to standard map-making and other map-making methods Fast all-sky convolution for polarimetry measurements PReBeaM Implementation 4.7.1 Polynomial Interpolation and Zero-Padding 4.7.2 Parallelization Description 4.8 Simulations and Beams 4.9 Results and Discussion 4.9.1 Computational Considerations 4.10 Conclusion 4.4 4.5 4.6 4.7 43 43 46 48 49 49 50 51 54 62 65 Chapter 5 Noise and Diffuse Foregrounds 5.1 Noise model for PLANCK 5.2 Testing PReBeaM on CMB + white noise 5.3 Knowledge of polarized foregrounds 5.3.1 Synchrotron emission 5.3.2 Dust emission 5.4 Testing PReBeaM on polarized foregrounds 5.5 Conclusions 66 66 67 67 67 69 70 72 Chapter 6 Variable Point Sources 6.1 WMAP findings on point sources 6.2 Description of data extraction and processing 6.2.1 Loss imbalance parameter 6.2.2 Differencing assembly parameters 6.2.3 Beam normalization 6.2.4 Dipole corrections 6.3 Chi-square analysis 6.3.1 Null test 6.4 Results and discussion 6.5 Conclusions and further work 76 76 78 78 79 79 80 81 82 82 87 Chapter 7 Conclusions 7.1 Summary of Progress 7.2 Future work with PReBeaM 7.2.1 Improving the computational requirements 7.2.2 Interface with destriper for noise treatment 93 93 93 93 94 viii 7.3 7.2.3 Re-analysis of the WMAP data 7.2.4 Analysis of PLANCK data with far sidelobes Future work with the convolution routines References 94 95 95 97 Vita 102 IX List of Tables 3.1 Fractional RMS error for each of the six main test cases 32 4.1 The four detectors of the 30GHz PLANCK channel and their FWHM, ellipticity, angle of polarization measurement ippoi, and orientation of beam major axis tpeii 54 6.1 6.2 Beamsize (9b), loss imbalance parameter (xjm), U>A,B, and noise term (ao) for each of the differencing assemblies (DA) Point sources that we have identified as being variable. Sources are listed, following the WMAP convention, by their 5 GHz ID. Variability observed in either the A or B side (or both) is indicated by a checkmark for the K and Ka band separately. The peak K and Ka band flux is listed in columns 6 and 7, respectively. The type (G for galaxy, and QSO for quasar) is listed in column 8. v WMAP-identified as probable variable: x 2 > 36.7. v WMAP-identified as variable: \2 > 100 x 80 86 List of Figures 2.1 This plot shows a model for the dynamics of the scalar field during inflation and the transformation of fluctuations in the scalar field into density fluctuations. Accelerated expansion of the universe occurs during the slow-roll phase. Inflation ends when the field reaches the minimum of the potential, oscillates and decays into ordinary radiation and matter. Points on the potential perturbed by a fluctuation S(j) finish inflation at different times 5t, inducing a density fluctuation 5 = HSt. Reprinted with permission from [1] 2.2 Current measurements of the polarized CMB signal. The TE measurements (grey) are from the first-year WMAP data. The E-mode measurements (coloured) are from CAPMAP, CBI, DASI and Boomerang. The B-mode spectra is shown in red and blue for two cases, r = 0.30 and r = 0.01, respectively. B-mode due to lensing is shown in green. The black curves are the predicted one-sigma sensitivity for WMAP and PLANCK (neglecting correction for foreground emission). Figure reproduced from [2] On the left is a figure showing the orbit of the PLANCK satellite in the Earth-Sun system. The figure on the right shows the main beam of the telescope sweeping out large circles on the sky, while Galactic straylight is entering the detectors through the far sidelobe, seen at 90° from the main beam. Figure modified from [3] 3.2 A map of the spatial structure of the far sidelobes for one 30 GHz channel of PLANCK. The main beam would sit in the region at the pole. 3.3 Measured WMAP focal plane co- and cross-polar beams, shown on the left and right, respectively. The contours are spaced by 3 dB and the maximum value of the gain in dB is given next to the individual beams. Figure reproduced from [4] by permission of the A AS 3.4 Ratios of the spectra of the residual map to the spectrum of the input map for each of the beam models and both scanning strategies. The BSP results are plotted in the left column and the WSP results are plotted in the right column. Results for the sidelobe, elliptical and two-beam beam are shown in the top, middle, and bottom panels, respectively. The solid lines correspond to deconvolved spectra and the dashed lines correspond to the standard spectra 9 14 3.1 xi 20 20 21 33 3.5 Residuals after the first 100 iterations are shown in the first three rows. The figures on the left are for the basic scan path, and the ones on the right are for the WMAP-like scan path. First, second, and third rows correspond to the sidelobe, elliptical and two-beam beams, respectively. The true sky is shown in the fourth row. Note that solutions of the two-beam BSP and, to a lesser extent, the elliptical beam BSP test cases have not converged to sufficient accuracy. We chose to present the results for all cases after a fixed number of iterations to show the impact of scanning strategy and beam pattern on the condition number of the map-making equations 3.6 Convergence rates of the preconditioned conjugate gradient solver for each test case. The left panel refers to the basic scan path and the right panel to the WMAP-like scan path. The solid lines correspond to the sidelobe beam, dotted lines to the elliptical beam, and dot-dashed to the two-beam model 3.7 Deconvolving the effects of a large sidelobe in simulated observations of the WMAP Ka band map, using the coarsened WMAP scanning strategy described in the text. The top map is the input sky map, the middle map is the standard map-making result, and the bottom map is the deconvolved result 3.8 The top plot shows the distribution of point sources in the upper left quadrant of the CMB sky (note that the upper end of the color scale has been cut-off at 93.6^K; the original scale terminates at 374fJ,K). The middle and bottom plots show the difference maps between the input and output skies for m — 0 (middle) and m = 4 (bottom). It can be seen that the overall residual in the m — 4 case is roughly five times smaller than the m = 0 case. We also find that point sources in regions of the sky that are observed by the beam at a wide range of orientations are reconstructed with smaller errors 3.9 Examing the effect of beam ellipticity and subsequent deconvolution on point sources in the CMB. The actual point source (top) is scanned with an elliptical beam and reconstucted using deconvolution mapmaking with the asymmetry parameter set t o m = 0 (bottom left) and m — 4 (bottom right) 4.1 4.2 Comparison of one segment of the sky from a binned noiseless map (left) and the PReBeaM temperature map (right), both at a Healpix resolution of nside 1024. At this resolution, the binned map contains a number of unobserved pixels, of which three are visible in this frame. The PReBeaM map contains no unobserved pixels Forward interpolation from ring set to TOD element and transpose interpolation from TOD element to ring set xii 34 35 36 38 39 47 50 4.3 4.4 4.5 4.6 4.7 4.8 4.9 4.10 Depiction of hybrid parallelization scheme used in PReBeaM. Rectangles represent work done on a node, ellipses represent data products, and arrows represent transfer of data. The work done within a node (convolution, interpolation and their transpose operations) is parallelized using OpenMP. This figure shows a slice of two head nodes, however, the algorithm may operate on many more nodes Plot showing the footprint of the LFI main beams on PLANCK as seen along the optical axis looking towards the satellite. The pair of horns comprising the 30 GHz channel are labeled 27 and 28 in the plot. The scan direction and polarization orientations for the co-polar beams are also depicted. Figure reproduced from [5] with permission Contour plots in the wu-plane of the PLANCK 30 GHz main beam detectors. This coordinate system permits the beam to be mapped from a spherical surface to a plane. From left to right, the beams are: LFI-27a, LFI-27b, LFI-28a, and LFI-28b. Figure reproduced from [5] with permission EE, TE, and BB spectra of smoothed input map (black curve), binned map (blue curve) and PReBeaM (red curve). Note that the red and black curves are nearly on top of each other. The EE and TE spectra show the effect of temperature-to-polarization cross-coupling seen in the binned map spectra as shifts in the peaks and valleys and absent from the PReBeaM spectra. The input BB spectra is absent from the BB plot since the input B-modes were zero. For comparison, we show a theoretical B-mode spectrum (dashed curve) from a cosmological model with a 10% tensor-to-scalar ratio. TT spectra are omitted since differences in the three spectra are not apparent in this representation. Fractional difference in power spectrum for PReBeaM (red, blue, and cyan curves) and the binned map (black curve) for TT, EE, and TE. Spectra for PReBeaM are shown as a function of number of iterations to demonstrate convergence Power spectra of the difference map for the binned map (black) and for PReBeam (red). For comparison, we calculate the power spectra from the difference in a^m's (aemout — aemin) for the PReBeaM method (blue). In the case of the BB spectra, the blue curve lies exactly on top of the red curve since the input afm were zero The absolute value of the beam aem coefficients for m = 0 (solid), m = 2 (dotted), m = 4 (dashed), m = 6 (dash-dot), and m = 8 (dash-dot-dot-dot). We plot \ajm\ (left) and \afm\ (right) (we omit the B component plot since \afm\ = \afm\) BB power spectra as a function of asymmetry parameter mmax for m max = 2 (blue curve) , 4 (cyan curve), and 6 (red curve). The input BB spectra was zero so the smallest output BB spectra is most desirable. In this run, the PReBeaM input parameters interpolation order and zero-padding were set to 1 and 2, respectively xin 52 53 54 56 58 59 60 61 4.11 Worst case bias in estimation of cosmological parameters due to errors in the power spectra of PReBeaM (dashed curve) and due to the errors in the power spectra of the binned map (solid curve). These plots estimate the bias due to the 30 GHz channel only, and should not be taken as representative of PLANCK as a whole 4.12 The residuals between the input reference sky and PReBeaM output (left column) and the residual between the input reference sky and the binned map (right column) for Temperature (T), and the Stokes Q, and U parameters 5.1 5.2 5.3 5.4 5.5 TE spectrum of CMB and white noise for Springtide (blue curve) and PReBeaM (red curve). The smoothed input map is shown in black. Following the example in [5], we reduce £ to £ variation by filtering the spectra by a sliding average (A£ = 20). In this run, the PReBeaM input parameters interpolation-order and zero-padding were set to 1 and 4, respectively. While PReBeaM performs at least as well as Springtide in the TT, EE, and BB spectra, we omit these spectra since the detailed differences are difficult to assess without an in-depth Monte Carlo study Simulated map for synchrotron polarization amplitude at 30 GHz. . . Simulated map for dust polarization amplitude at 30 GHz Fractional difference in power spectrum for PReBeaM (red curve) and Madam (black curve) for TT, EE, BB, and TE Output spectra from PReBeaM (red curve) and Madam (black curve) for free-free emission only. Since free-free emission is temperatureonly, all of the resulting power in the polarization modes is a result of leakage from temperature to polarization through beam effects. Note that the PReBeaM TT spectra nearly exactly overlays the binned map TT spectra Map of the 390 point sources detected by the WMAP team. The shaded region shows the mask used to exclude extended foreground emission. The size of the plotted points indicates the flux of the source: the area of the dot scales like the maximum flux over the 5 WMAP bands plus 4 Jy. Image courtesy of the WMAP science team 6.2 Innermost parts of the normalized profiles for all 10 channels of the WMAP instrument 6.3 The 24 most variable point sources found on the A side of the K band. The distribution of the x2 probabilities for the null points plotted as a histogram for the A side (blue) and B side (green). The x2 probability, p, for the actual point sources are overplotted for the A side (dotted) and B side (dashed). This figure shows the point sources ordered in increasing p for side A. Note that the dotted and dashed lines for sources J0423-0120, J1229+0202, J1256-0547, J2258-2758, and J05384405 are difficult to see as they sit nearly exactly at zero 63 64 68 71 71 73 74 6.1 xiv 77 88 89 4 The 24 most variable point sources found on the B side of the K band. The distribution of the x 2 probabilities for the null points plotted as a histogram for the A side (blue) and B side (green). The x 2 probability, p, for the actual point sources are overplotted for the A side (dotted) and B side (dashed). This figure shows the point sources ordered in increasing p for side B. Note that the dotted and dashed lines for sources J0423-0120, J1229+0202, J1256-0547, J2258-2758, and J05384405 are difficult to see as they sit nearly exactly at zero 5 The 24 most variable point sources found on the A side of the Ka band. The distribution of the x 2 probabilities for the null points plotted as a histogram for the A side (blue) and B side (green). The x 2 probability, p, for the actual point sources are overplotted for the A side (dotted) and B side (dashed). This figure shows the point sources ordered in increasing p for side A. Note that the dotted and dashed lines for sources J1258-0547, J1229+0202, J2258-2758, J0423-0120, J0538-4405, and J1014-4508 are difficult to see as they sit nearly exactly at zero. . 6 The 24 most variable point sources found on the B side of the Ka band. The distribution of the x 2 probabilities for the null points plotted as a histogram for the A side (blue) and B side (green). The x 2 probability, p, for the actual point sources are overplotted for the A side (dotted) and B side (dashed). This figure shows the point sources ordered in increasing p for side B. Note that the dotted and dashed lines for sources J1258-0547, J0423-0120, J0538-4405, J1229+0202, J2258-2758, and J1014-4508 are difficult to see as they sit nearly exactly at zero. . xv 90 91 92 Chapter 1 Introduction Ever since its discovery by Penzias and Wilson in 1965 [6], the CMB has yielded a wealth of cosmological information. Immediately following the serendipitous discovery of the CMB, researchers at Princeton University led by Robert Dicke interpreted this excess antenna temperature to be the relic radiation of the Big Bang [7]. Twentyseven years later, in 1992, the Cosmic Background Explorer (COBE) mission first reported a measurement of the tiny fluctuations in the temperature of the CMB [8] as predicted by General Relativity and the theory of inflation. Most recently, five years of measurements of both the temperature and polarization anisotropies by the Wilkinson Microwave Anisotropy Probe (WMAP) 1 with 45 times the sensitivity of COBE and 33 times the angular resolution have constrained the model of our universe to fantastic precision[9]. By measuring the statistical properties of the fluctuations in the CMB, cosmologists have been able to refine a standard model of cosmology. This standard model consists of a flat universe in which roughly 74% of the total energy density is dark energy and 22% is dark matter, with the remainder in the baryon density [9]. The standard model also predicts nearly scale-invariant adiabatic Gaussian fluctuations. Cosmologists are now interested in answering new questions about our universe, including: Did inflation seed the primordial fluctuations? And, if so, what is the correct model for inflation? The inflationary model of the universe was first described in 1981 by Alan Guth as a means of solving two major puzzles in the standard model of the early universe [10]. The first puzzle is known as the horizon problem, in which causally disconnected regions of the universe are found to be homogeneous [10]. The second puzzle, referred to as the flatness problem, indicates that the critical value of the energy density of the universe today (which gives us a flat universe) requires fine-tuning of initial conditions to incredible precision. While current models of inflation (see [11] and [12] for recent reviews) are incredi1 http://map. gsfc.nasa.gov 1 bly successful at describing the observed universe, a distinct class of models has been developed which also predicts these same conditions. This new class of theories is called the ekpyrotic [13] or cyclic model [14], and is based on a multi-dimensional framework motivated by string theory. Fortunately, there exists a very promising means of experimentally verifying the standard inflationary cosmology model through the detection of the presence of a gravitational wave background (GWB). If present, a GWB will leave a unique signature, known as B-modes, in the polarization anisotropics of the CMB. Thus, a detection of B-mode polarization in the CMB would be a direct probe of inflationary gravitational waves. One of the most exciting prospects for current and up-coming CMB satellite missions, such as PLANCK 2 and CMBPol, is their potential to measure the B-mode polarization in the CMB. Measurement of this signal will depend on two things. First, it will depend on the actual amplitude of the B-mode power spectrum, which is uncertain and related to the energy scale of inflation. Second, it will depend on a successful treatment of the data in which the overwhelming signal from intervening foregrounds is removed and the systematics are corrected. PLANCK is a satellite mission of the European Space Agency that is due to be launched in May 2009. PLANCK is imaging the sky in nine frequency channels from 30 GHz to 857 GHz with angular resolution down to 5'. While PLANCK is not specifically designed to detect primordial B-modes, it does have the potential to do so. A better chance for B-mode detection will come with the future CMBPol mission, [15] designed specifically to measure this elusive signal. The process by which CMB data are collected suffers from both instrumental and astrophysical systematic effects. Instrumental effects include asymmetries of the beam in the main focal plane, and the pick-up of large foreground signals far from the focal plane in the side lobes. We use the term beam throughout this dissertation to describe the response of the detector instruments to the photons being measured. Astrophysical systematic effects include all Galactic and extra-galactic signals appearing in the frequencies where the CMB is measured. In an optimal treatment of CMB data, these effects must be corrected during the "map-making" step, the method by which time-ordered data (TOD) are turned into a sky map. Power spectra are generated from maps, which are then used to estimate cosmological parameters, such as the density of baryons fl^h2, the density of cold dark matter Qch2, and the amplitude of curvature perturbations A27£ . Systematic effects that are left untreated in the 2 http://astro.estec.esa.nl/Planck 2 maps may be propogated through the data analysis pipeline and may lead to errors in estimated parameters. In this dissertation, we describe the development of a novel technique for making maps of the CMB from satellite-generated instruments that fully addresses issues due to beams, optics, and foregrounds. The chapters are organized both by chronological development of our technique, and also by the increasing complexity of the problem being examined. In order to provide motivation for the quest for precision CMB polarization measurements, we review the main principles of CMB physics and inflationary theory in Chapter 2. We briefly compare inflationary models for the formation of metric perturbations in the universe with cyclic models. This is followed by a review of the theory, mathematics and formalism for the formation of CMB anisotropies and polarization. Chapter 3 describes in detail the map-making problem, some standard solutions, and our unique approach to the solution, known as deconvolution map-making. We compare the results from our approach to those from a standard map-making method in a variety of test cases that examine effects of scan strategies, beams, and foregrounds. Results from this chapter are restricted to temperature-only measurements of the CMB. In Chapter 4 we extend our description of deconvolution map-making from the previous chapter to include polarization measurements of the CMB. We describe the development of our technique under the framework of analysis for the PLANCK mission and we adopt the name PReBeaM for our tool. We introduce the objectives and parameters of the PLANCK mission and discuss the beams and simulated data. We describe interpolation and parallelization features that we implemented in order to meet the challenges that are faced when manipulating such large data sets and when attempting to make accurate polarization maps. Following the successful performance of PReBeaM in the cases of temperatureonly and temperature-plus-polarization-only CMB data, we introduce white noise and foregrounds into the data in Chapter 5. We compare the performance of our code to a map-making method currently employed by the PLANCK Algorithm Development group in the US. Point sources are an important foreground in CMB maps and some extraction methods rely on a complete physical characterization. In Chapter 6, we describe a program to identify point source variability in the public WMAP time-ordered data. An implicit assumption in our work with the convolution method is that neither 3 the sky nor the beam changes over the course of the survey. Variable sources are one contributor to time variation in the sky that have received little attention in the literature to date. This study of point source variability will help to provide important information about the structure of extragalactic radio sources over a wide frequency range. We conclude the thesis in Chapter 7 with a summary of the results and a plan for future work with the PReBeaM code. 4 Chapter 2 CMB Theory 2.1 Introduction In this chapter we introduce the standard model of cosmology and describe the formation of primordial density perturbations. The physics of the inflationary model for the beginning of the Universe is examined in detail. Following that, we discuss similarities and differences between the inflationary model and the cyclic class of models in terms of their specific predictions about the formation of scalar perturbations and gravitational waves. Next, we discuss the formation of temperature and polarization anisotropies in the CMB and the relationship between CMB polarization and gravitational waves. We conclude this chapter by reviewing the Stokes parameters, standard description for a polarized field and the application of harmonic analysis to temperature and polarization anisotropies. Over 40 years of study of the CMB has led to a very sound understanding of the underlying physical principles. There are many texts that present in great detail descriptions of the formation of initial perturbations, and temperature and polarization anisotropies. For a more detailed look at CMB theory than is presented in this chapter, we recommend the following texts: [16, 17, 1, 18]. 2.2 The Standard Cosmological Model Standard hot big bang cosmology describes a universe expanded from a hot dense state, dominated by thermal blackbody radiation, and evolved under the laws of Einstein's general theory of relativity. The Friedmann-Robertson-Walker (FRW) model [19, 20, 21, 22, 23] was developed as the simplest possible model obeying the basic assumptions that (i) the universe is the same at all points (homogeneous) and (ii) all spatial directions at a point are equivalent (isotropic). In this work, we assume the 5 FRW framework for the universe, the metric for which is given by ds2 = dt2 - R2(t)[, dr2 , „ + r2(d92 + sin2(e)d(f,2}. 1 — kr2 (2.1) Note that the speed of light c has been set to 1 throughout. This equation for the line element ds is characterized by the cosmological scale factor R(t) which determines proper distances in terms of the comoving coordinates, and the constant A; which describes the curvature of the universe. By a rescaling of the radial coordinate, we can restrict k to be either 0, -1 or 1 corresponding to a flat, open or closed universe respectively. R(t) is commonly normalized by R0 = R(t0), the value of R at the present epoch, defmining the dimensionless scale factor a(t) = R(t)/R0. The rate of change of the scale factor describes the background evolution of the universe and is known as the Hubble parameter H(t) = R(t)/R(t). Prom Einstein's field equations, we can define two cosmological equations of motion that govern the expansion and acceleration of the universe: and R R A 3 4irG 3 * + %>)• (2 3 ' > In the above equations, G is the gravitational constant, p is the energy density, A is the cosmological constant, and p is the isotropic pressure. In addition to the Hubble parameter, it is useful to define a handful of other cosmological parameters. A critical density can be derived from Equation 2.2 by requiring k = 0 when A = 0: AW - 8 ^ - (2-4) The cosmological density parameter Qtot is defined as the ratio of the energy density to the critical density Oot(t) = P(t)/Pc(t) (2.5) and we often consider the individual contributions of densities as Clx(t) — px(t)/pc(t). Today, our best estimates for the composition of our universe are: cold dark matter, Qc = 0.228±0.013; dark energy, 0 A = 0.726±0.015; and baryons, Qb = 0.0456±0.0015 [24]. In this section, we have described the background evolution of the universe un- 6 der the assumptions that it is homogeneous and isotropic on large scales. However, observations from telescopes show structures on a wide range of scales, from single galaxies to superclusters of galaxies. The existance of these structures tells us something about the initial conditions of the Big Bang. We review these special initial conditions for the formation of density fluctuations in the next section. 2.3 Models for the formation of perturbations In this section, we review a Fourier analysis of cosmological density fields. Then we describe the inflationary scenario and its predictions for the formation of the initial perturbations. We compare and contrast these predictions with those of the cyclic class of models for the early universe. 2.3.1 Fourier analysis of density fluctuations Gravitational instability is understood to leave an imprint of inhomogeneity and structure in the universe as a result of a primordial field of density fluctuations, <5(x). We briefly summarize the statistical measures used to study <5(x). The dimensionless density perturbation field is defined as *(x) ^ * * ^ M . (2.6) An important feature of 8 is that it occupies a universe that is isotropic and homogenous in its large scale properties. This suggests that the statistical properties of 5 should be similar to those of a stationary random process. We often expand 5 as a Fourier superposition and assume periodic boundary conditions in a cube of some large volume: 5(x) = ^ 4 e - t k x (2.7) where k is the wavenumber of the fluctuations. Then power spectrum is defined as ensemble-average of 5^ V(k) = (|4| 2 )- 7 (2.8) 2.3.2 Inflationary models Inflationary models [10, 25, 26] are based on the notion that there was a beginning of the universe in both space and time, with nearly infinite temperature and density. At a time of order 10~35 seconds (corresponding to an energy scale of 1016GeV), the universe experienced a period of accelerated expansion. In the simplest models, inflation is driven by a scalar field and excites fluctuations in all available modes. This period of accelerated expansion inflates quantum fluctuations in the scalar field to sizes much larger than the Hubble radius, after which gravitational instability leads to fluctuations in the density of the universe. Quantum fluctuations in the metric itself give rise to a spectrum of gravitational waves or tensor perturbations. We now look more closely at the inflationary scenario. Considering Equation 2.2, we find that each term on the right hand side has a different time-dependence. At high redshift, space curvature and the cosmological constant are negligible compared to the mass density so the cosmological equation for the expansion rate becomes h /8 \1/2 H{t) = -a= (JTTGPJ • (2.9) Local energy conservation implies that p=-3-(p + p). (2.10) a Inflation resolves the puzzle that the observable universe is homogeneous by assuming that there was a time when the net pressure was negative, V < -P/3 (2.11) and the universe underwent a period of accelerated expansion. The negative pressure in the inflation epoch is driven by a scalar field 4>, sometimes called the inflaton, with Langrangian density £ = \^4>jj ~ V{<j>). (2.12) In the limit that the field is nearly spatially homogenous, the energy density and 8 Figure 2.1: This plot shows a model for the dynamics of the scalar field during inflation and the transformation of fluctuations in the scalar field into density fluctuations. Accelerated expansion of the universe occurs during the slow-roll phase. Inflation ends when the field reaches the minimum of the potential, oscillates and decays into ordinary radiation and matter. Points on the potential perturbed by a fluctuation 5(f) finish inflation at different times St, inducing a density fluctuation S = H5t. Reprinted with permission from [1]. effective pressure in the field are P4> = <j>2/2 + V, P(j) = <j>2/2-V. (2.13) Thus, the condition for accelerated expansion dominated by the scalar field 0 is (d(/)/dt)2 < V(<p). This corresponds to the case when the kinetic energy of the field is much smaller than the potential energy, a condition which is referred to as the slow-roll of the scalar field. We depict a simple inflation scenario in Figure 2.1. Inflation ends as (p approaches a minimum of the potential and oscillates rapidly. Couplings between the scalar field and matter fields will produce particles, leading to what is called reheating. The key inflationary observables are the level of anisotropy due to scalar perturbations and the level of anisotropy due to tensor perturbations. The exact relationship between the level of anisotropy due to scalar and tensor perturbations and the inflaton potential depends upon the cosmological model (see, for example, [27]). If these can be measured, along with the power-law index, ns, that characterizes the density perturbations, then the scalar-field potential that drove inflation can be reconstructed. These parameters are described in Section 2.4. 9 2.3.3 The cyclic model The cyclic model [13, 14] predicts that the universe is infinite in space and time, undergoing endless cycles of "big bangs" and "big crunches". Each cycle begins with a bang, which is a transition from a contracting phase to an expanding phase, in which the density and temperature of the universe do not diverge. This transition results from a collision of two three-dimensional domain walls known as 'branes', and results in a period of radiation-domination and a subsequent period of matter-domination. Following this is a dark energy dominated epoch, where the dark energy is due to a scalar field (f> rolling down a shallow potential V(<p). Acting as dark energy, <fi causes a period of slow acceleration, converts the acceleration into deceleration, contracts into a "big crunch" and begins the cycle again. The dark energy phase is responsible for making the universe flat and homogeneous and provides the source for the presently observed dark energy and cosmic acceleration. Expansion slows down and reverses to contraction when <j> rolls to values where V{4>) is negative. Contraction occurs very slowly at first, and during this period quantum fluctuations have time to cause spatial variations in the rate of contraction. It has been shown that the resulting spectrum is scale-invariant [28], and evolves to become a spectrum of temperature and density perturbations after the bang. Current calculations [29] estimate that string-theory motivated models do not produce measurable gravitational waves. 2.4 Inflationary gravitational waves Inflationary gravitational waves are tensor perturbations in the background spacetime metric and were first proposed by [30]. Quantum fluctuations in the scalar field driving inflation lead to scalar metric perturbations, while quantum fluctuations in the metric itself lead to tensor metric perturbations. In a spatially-flat, FriedmannRobertson-Walker universe, the perturbations are described by the line element, ds2 = a2[-dr2 + (8i:j + hij)dxidxj] (2.14) where r is the conformal time, x% are comoving spatial coordinates, and hij is the tensor metric perturbation. The scalar and tensor power spectra predicted by most inflationary models ap- 10 proximately follow power laws A (/c) * = Sy{llZkl2) and 2k3 Al(k)^j^y2(\h+k\* K kns+l + \hxk\2)<xkn> (2 15) " (2.16) where A^ and A\ are the variance due to the scalar and tensor modes respectively. Here, 1Z is the curvature perturbation, and h+ and hx are the two polarization states of the primordial tensor perturbation. They are defined by <TC2> = / f A^(fc) (2.17) T*l(k). (2.18) rlk / The scalar spectral index ns and the tensor spectral index nt are nearly scale invariant (i.e., ns = 1 and nt = 0). The current best estimate on ns from WMAP with constraints from Type la Supernovae and Baryon Acoustic Oscillations (BAO) is ns = 0.960 ± 0.013 [9]. The tensor contribution, quantified by the tensor-to-scalar ratio r and evaluated at A;0 = 0.002Mpc_1, is The combination of the most recent WMAP results + Type la supernovae events + BAO puts an upper limit on the tensor-to-scalar ratio, r < 0.22 (95% CL) [24]. Since gravitational waves are small in amplitude they evolve according to linear theory. Thus, we can predict their evolution with high precision. Gravitational waves are a useful probe of the early universe as they imprint essentially unperturbed information on the CMB polarization anisotropy. 2.5 Formation of Cosmic Microwave Background Anisotropies In both the inflationary and cyclic model, our universe began as a hot, dense plasma of photons, baryons and electrons. This ionized matter underwent Thomson scattering, which effectively coupled radiation to matter. As the universe expanded, it cooled 11 adiabatically. At a temperature of roughly 3000K, it had cooled sufficiently so that protons and electrons could combine to form hydrogen. As a result, the opacity of the universe from Thomson scattering was reduced and radiation decoupled from matter. This event occurred about 380,000 years after the big bang and is known as last scattering. The CMB photons that we detect today have been redshifted by the changing scale-factor of the universe, thus the CMB now has a mean temperature of about 2.725K. The level at which the CMB is anisotropic is about one part in 100,000, or fluctuations on the microKelvin scale. These anisotropics are generated through well understood physics of the evolution of linear perturbations with a background FRW cosmology. The three mechanisms for the formation of primary anisotropies can be roughly separated by the range of scales that they each operate on: the Sachs-Wolfe effect (£ < 100); the acoustic oscillations (100 < £ < 1000); and the damping tail (*>1000). Photons which last-scattered from an overdense region (or potential well) lose energy from climbing out of the well and are redshifted, becoming cooler. Photons which last-scattered from an underdense region, gain energy and are blue-shifted and become hotter. This is the Sachs-Wolfe effect [31, 32], which contributes mainly to anisotropies on scales larger than about one degree, or £ < 100. On sub-degree scales, anisotropy is due to acoustic oscillations in the baryonradiation fluid at decoupling [33]. Perturbations in the gravitational field drive oscillations in the photon-baryon fluid. Regions of compressed fluid contain hot photons, while regions of rarefied fluid contain cooler photons. A damping of the anisotropies occurs at the highest £s because recombination is not instantaneous. The photon-baryon fluid exhibits imperfect couplings, and the amplitude of the oscillations decreases with time. This damping tail, also known as Silk damping [34], corresponds to scales smaller than that subtended by the thickness of the last scattering surface. 2.6 Theory of CMB polarization Polarization of the CMB is generated when an anisotropic radiation field Thomson scatters off free electrons. For incident radiation with a quadrupole moment, linear polarization results from the scattering [35, 36]. These polarization fluctuations will be an order of magnitude smaller than temperature fluctuations as the scattering will occur near the end of recombination, when there are fewer scattering sites. 12 Polarization is also generated by scattering at later times when the universe reionizes. Therefore, CMB polarization can also provide information about the reionization history of the universe and may shed light on the nature of polarized foregrounds, such as our Galaxy and thermal emission from dust. The contribution to anisotropy which leads to polarization may either be due to primordial density or gravitational wave perturbations. The primordial density fluctuations are scalar and can only generate E-modes, a curl-free pattern classified by polarization directions that are either purely radial or purely tangential to the hot and cold temperature spots. On the other hand, the spin-2 nature of tensor perturbations (or inflationary gravitational waves) permits both the generation of E-modes and B-modes [37, 38, 39], a divergence-free polarization pattern. The amplitude of B-modes depends linearly upon the amplitude of gravitational waves during inflation and is parameterized by the tensor-to-scalar power ratio (r), as shown in Figure 2.2. In 2002, CMB polarization was first detected at sub-degree angular scales by the DASI team [40, 41]. These results were confirmed and extended by results from CBI [42], CAPMAP [43], and Boomerang [44]. More recently, with the release of the third-year WMAP results, the first large-scale detection of E-mode polarization has been made [45]. On the largest angular scales, CMB polarization should also reveal information about the first stars that formed in our universe. In addition, gravitational lensing by clustered matter at low redshifts should be detected in polarization and temperature anisotropics. Though a measurement of primordial B-mode polarization in the CMB would be direct evidence of gravitational waves from inflation, a measurement of B-mode due to lensing may be more likely. E- and B-modes can mix through a weak lensing shear. The effect on B-modes due to this mixing is much more pronounced than for E-modes since the amplitude of E is much greater. Figure 2.2 demonstrates that primordial B-modes dominate on large scales, but B-modes due to lensing dominates on small scales. While the measurement of the power spectrum of B-modes due to weak-lensing would reveal important information about the intervening large-scale structure of the universe, only a primordial gravitational wave detection would confirm an inflationary epoch and set a corresponding energy scale. 13 OH < S 'd CD T 3 CP o « .2 cp - d O >> a <p CO r^H £P CJ> cP CP S-l cp [ i . pi .y CO o B CP cp co * CQ ho a « 2 2 •aH x) 03 r ^ l >-> . ^s § s-l co cp cp O s ° § PQ co a a ?i o S-l CO CO CO ' " S=! cp d O n o cp s-i S-l S CO CP ^ H-r OQ fr> O O O •s if Is CQ faO cp m CO pL, m <;••3 CP *• s^ ua g <• co J™ <p CL, CP 4 3 ^ H . 2 CP O Pi cc3 O CU O ,—. CP ^ £ m i.n„m.„,i„„,A. 9 o •JLumm&m pi AnAA. o o CO co Pi o CP II !B '-u S-l £- CO Pi co" CP CO CO [>tn] i v >> +j' CP a o •a* cp PI a3 cp CP CO Pi /T^ o a hO CP s-l S-l O S-l S a .2 ce cci pi CM cc3 14 Pi O The E-: red nd blue redicted from UJTT cp S-l -^ o ° o •e 2.2: C C <£ s-i .2 cp CB o CO bO cp Cu CP 2.7 Review of Stokes parameters In this section, I will review the standard description of CMB polarization following the conventions in [38, 46]. A monochromatic electromagnetic wave propagating in the z-direction with frequency ui0 has electric field components Ex = ex{t) cos[ui0t - 9x(t)}, (2.20) Ey = ey(t) cos[v0t - 6y(t)} (2.21) Linear polarization can be decomposed into two components, the Q and U Stokes parameters, denned as Q = (el)-(e2y), (2.22) U = {2exeycos(9x-ey)). (2.23) It can be shown [39, 38], that, for a coordinate system rotation of a, the polarization field is now described by Q' = Qcos(2a) + f/sin(2a), (2.24) U' = -Qsin(2a) + f/cos(2a) (2.25) We will show the connection between the Stokes Q and U parameters and the E and B modes in the next section. 2.8 Harmonic analysis for temperature and polarization anisotropies The CMB is described as a field of temperature fluctuations on the sky, AT/T(n), which can be expanded as a sum over the spherical harmonic coefficients ajm: AT — (n) = ^ a j m y , m ( n ) . (2.26) We say that the a^m are band-limited if they are zero for some £ > £max- The temperature power spectrum is given by C7 = (|aLl2>. 15 (2-27) Similarly, Q and U can be expanded in the spin-2 basis: (Q + U)(n) = J2a^rn2Yem{n), (2.28) (Q~U)(h) = (2.29) J]a_2,,m_2^m(n). Now it is convenient to introduce the linear combinations of a 2 / m and a_2,em'afm = -(a2,em + a-2,em)/2 (2.30) afm = i(a2,em - a_ 2 / m )/2 (2.31) It follows that afm and afm can be constructed from the Q and U Stokes parameters, afm ± iaBm = -J dh±2Y;m(h)(Q ± iU). (2.32) For real space calculations it is convenient to define two scalar quantities E(h) and B(h) corresponding to the divergence and curl parts of linear polarization a,EmY,m(n), B(n) = £ af m Y, m (n). E(n) = £ (2.33) The existence of linear polarization leads to six possible cross power spectra. The general form {almat'm<) = Ce Su'Smm' (2.34) for X, X' e {T, E, B}, which gives: CjT = (\aJJ), B 2 Cf Q CfE = (KEm|2) E = (\afm\ ), = iaemaem)i (2.35) CJ = (ajma^) (2.36) Ce (2.37) = {aernae£) Two of the spectra, CjB and CfB, are eliminated by parity. 2.9 Conclusions I have reviewed some of the physics of the formation of CMB temperature and polarization anisotropics and have shown the connection between inflationary induced primoridal gravitation waves and the unique B-mode polarization signal. In the next 16 chapter I will present the map-making problem and show how beam asymmetries and foregrounds may induce systematics that will distort maps. 17 Chapter 3 Deconvolution Map-Making 3.1 Summary of chapter We describe a new map-making code for cosmic microwave background (CMB) observations. It implements fast algorithms for convolution and transpose convolution of two functions on the sphere [47]. Our code can account for arbitrary beam asymmetries and can be applied to any scanning strategy. We demonstrate the method using simulated time-ordered data for three beam models and two scanning patterns, including a coarsened version of the WMAP strategy. We quantitatively compare our results with a standard map-making method and demonstrate that the true sky is recovered with high accuracy using deconvolution map-making.1 3.2 Introduction In Chapter 2, we reviewed the formation of the CMB and provided motivation for seeking out precision measurements of the temperature and polarization anisotropies in the CMB. We focused, in particular, on precison polarization measurements, and the significance of a B-mode detection for inflationary physics. We will revisit polarization of the CMB and B-modes in Chapter 4. In this chapter, we will focus on some of the challenges faced when taking measurements of the temperature anisotropies of the CMB. We will discuss the problem of going from time-ordered data to maps of the CMB and some standard procedures for this data-processing step. Finally, we will present our deconvolution solution to this problem. We describe the general map-making problem for satellite instruments measuring the temperature anisotropies of the CMB in Section 3.3. A review of some standard map-making approaches is given in Section 3.4. Section 3.5 details our deconvolution approach to the map-making problem. Some of the algorithmic details of the deconvolution method are outlined in Section 3.6. In Section 3.7, we create a sim1 This chapter contains material previously published in [48]. 18 pie binning map-maker to use as our standard of comparison. The test cases that we examine are listed in Section 3.8, followed by the results in Section 3.9 We show the results in Section 3.10 from an addtional test using point sources observed with elliptical beams. We summarize the chapter in Section 3.11. 3.3 The map-making problem Real microwave telescopes collect distorted information about the CMB anisotropics due to asymmetries in the beam shape [49] and stray light from sources such as the Galaxy [50, 51]. To correct for these systematic errors we must be able to remove the detector response at all orientations of the telescope over the whole sky. In an optimal treatment, this correction must be applied during the map-making step of the CMB data analysis pipeline, before the angular power spectrum can be reconstructed. The problem becomes increasingly important as new generations of CMB observations probe for ever fainter signals in the CMB sky, and especially as we are preparing to measure the polarization of the CMB with high sensitivity. A microwave telescope scans the CMB sky according to some scanning strategy, recording time-ordered scan data. The scanning of the satellite and the line-of-sight of the detector aboard the satellite is depicted in Figure 3.1 for the PLANCK satellite (here, we make brief mentions of the PLANCK satellite; we will be speaking in much more detail about the PLANCK instrument in Chapter 4). The process by which the data recorded by the detector is "wrapped" back on to the sphere to create an image of the CMB sky is known as map-making. The map-making process is made extremely difficult by the combination of two major effects: distorted beam shapes and large foreground signals. In an ideal experiment, a telescope would observe the CMB sky away from undesired microwave sources, with a perfect delta-function beam. In reality, while the gain of the optical system does drop very rapidly away from the main beam line-of-sight, the signal from the Galactic foreground can be 1000 times brighter than the fluctuations in the CMB. This Galactic straylight enters the detectors through the far sidelobes, as depicted in Figure 3.1 and also in Figure 3.2. It is therefore important to characterize the beam, and use that information during map-making to deconvolve beam effects. An ideal linearly polarized detector would have only a co-polar response and no cross-polar response to signals from the orthogonal linear polarization. In reality, optics will induce a cross-polar response in the detector. This has the effect of mixing the Q and U signal, and therefore the mixing of the much larger E-mode into the 19 Figure 3.1: On the left is a figure showing the orbit of the PLANCK satellite in the Earth-Sun system. The figure on the right shows the main beam of the telescope sweeping out large circles on the sky, while Galactic straylight is entering the detectors through the far sidelobc, seen at 90° from the main beam. Figure modified from [3]. Figure 3.2: A map of the spatial structure of the far sidelobes for one 30 GHz channel of PLANCK. The main beam would sit in the region at the pole. 20 A side * por* measured I A l i l L_j -2 1 1 1 0 1 i i I 2 Cress-polar response i i__j I 4t l__i -4 i i I -2 l l i I C- i i i I 2 t i i I 4 Cross-e'evation (deg) Figure 3.3: Measured WMAP focal plane co- and cross-polar beams, shown on the left and right, respectively. The contours are spaced by 3 dB and the maximum value of the gain in dB is given next to the individual beams. Figure reproduced from [4] by permission of the AAS. B-mode. The co- and cross-polarization beam profiles measured for the WMAP experiment are shown in Figure 3.3, clearly depicting a non-Gaussian beam response and a nonzero cross-polar component. 3.3.1 Maximum-likelihood estimator of t h e t r u e sky Now we will review the path from observations to maps. The TOD vector is effectively a convolution of the true sky with a beam function. Mathematically, we represent this as the matrix multiplication of the observation matrix A with the npj xe rvector containing the true sky s returning a vector d, containing the UTOD samples of the time-ordered data, A s = d, (3.1) The matrix A can encode information about both the scanning strategy and the optics of the scanning instrument, in which case it is called a convolution operator, or it can simply be a pointing matrix in the case of non-deconvolution methods. Each row of the pointing matrix has only one nonzero element corresponding to the direction of the observed TOD. The least-squares estimate of the true sky, s, is given by the normal equation ATAs = ATd 21 (3.2) where A T is the transpose convolution operator or transpose pointing matrix. A realistic model of the data should also include noise As + n = d, (3.3) where n is a time-ordered vector containing the noise. For the case where the noise is stationary and uncorrelated in the time-domain, then Equation 3.2 is exact. For more general noise cases, the normal equation is modified as ATN-1As = ATN-M (3.4) where N = (nnT) is the noise correlation matrix. 3.4 Standard map-making methods Standard map-making methods generally employ one of two basic algorithms for making maps from TOD. The first is an implementation of the maximum-likelihood solution (as described in the previous section) and the second is an approximate solver known as the destriping method. This second method is motivated by the presence of low-frequency noise in the data, which, if left untreated, leads to stripes in the final map. Destriping requires no prior knowledge of the characteristics of the instrument noise and is able to provide an estimate of the low-frequency part of the instrument noise. Although the map-making process is but one part of the data-analysis pipeline for CMB experiments, it has now been recognized as a critical step which requires efficient and exact computational methods [52, 53]. A search of the literature reveals at least six unique algorithms developed and tested for application to PLANCK data by international groups: Madam [54], MapCUMBA [55], ROMA [56], Polar [53], MADmap [53], and Springtide [53]. MapCUMBA, ROMA, and MADmap use optimal algorithms, computing the minimum-variance map in the case of Gaussiandistributed, stationary detector noise. All three codes solve the maximum-likelihood map-making equations using the iterative conjugate gradient method and FFT techniques. In order to produce accurate results, these codes require a robust estimate of the noise power spectrum. Springtide and Madam are destripers and produce equivalent results to the optimal algorithms in the limit of short baseline length (baseline parameters are a feature of destripers that allow for tuning the accuracy of the output map and computational resources required). Springtide achieves reduced memory 22 requirements by compressing the data into binned ring maps Early work on the map-making problem have relied on the brute force method of direct matrix inversion. However, current and future CMB experiments, like WMAP and the PLANCK satellite, return enormous data sets that render the brute force method intractable. More recent advancements include map-making methods applicable to the latest experiments; however, many treat the beam like a perfect deltafunction (e.g. [55, 57]) or assume a symmetric beam profile (e.g. [58]), and thereby relegate the problem of treating a non-Gaussian radial response of the beam to subsequent stages in the data analysis [59, 60, 55, 61]. In this class, special techniques exist to deal with differential measurements like that of the Differential Microwave Radiometer (DMR) on the Cosmic Background Explorer (COBE) satellite [62] or WMAP [63, 64]. A Fourier method has been developed [65] to perform deconvolution but only for nonrotating asymmetric beams, which renders it inapplicable to the PLANCK data. Lastly, [66] present a method to remove the main beam distortion over patches of the sky for asymmetric, rotating beams but operate in pixel-space which is computationally more expensive than spherical-harmonic-space algorithms (for the same level of accuracy [47, 67]) and too time-consuming for full-sky high resolution maps such as those from PLANCK. 3.5 Deconvolution M a p - M a k i n g We call our approach deconvolution map-making, a generalization of existing CMB map-making techniques to solve the maximum likelihood map-making problem for arbitrary beam shapes. For sufficiently high signal-to-noise data this technique allows super-resolution imaging of the CMB from time-ordered scans. We implement our method using the exact algorithms for the convolution and transpose convolution of two arbitrary function on the sphere — in this case the sky and the beam — as detailed by Wandelt and Gorski in [47]. These fast methods for convolution and transpose convolution are efficient because they make use of the Fast Fourier Transform (FFT) algorithm. They are guaranteed to work to numerical precision for band-limited functions on the sphere. A microwave telescope scans the CMB sky according to some scanning strategy, effectively convolving the true sky with a beam function, and returns a vector d, containing the UTOD samples of the TOD as described in Equation 3.1. We call our matrix A the convolution operator; A encodes both the scanning strategy and some level of information about the optics of the CMB instrument. 23 Each sample of the TOD is modeled as the scalar product of a row of the matrix A with the sky s. Each of the UTOD rows of A contains a rotated map of the beam. In a given row the beam rotation corresponds to the orientation of the antenna at the point in time when the sample is taken. We will assume the beam shape and pointing of the satellite to be known. The observation matrix A generalizes the notion of the pointing matrix, which is often used in expositions of map-making algorithms by including both optics and scanning strategy. This generalization is necessary for any map-making method that accounts for beam functions with azimuthal structure. Equation 3.2 is exact if the noise is stationary and uncorrelated in the time-ordered domain. For a more general noise covariance matrix in the time-ordered domain N, the normal equation is modified to Equation 3.4. We proceed, considering only white noise in this chapter; however, it is straightforward to generalize to non-white noise as indicated in equation 3.4. Indeed, the matrix-vector operations required for this generalization have already been implemented in publically available map-making codes (e.g., MADmap 2 ). 3.5.1 Regularization technique The least-squares estimate of the true sky, s, is given by Equation 3.2. The coefficient matrix in this system of equations, A T A , is a smoothing matrix and hence ill-conditioned. Inverting it to solve Equation 3.2 therefore poses a problem. We describe here a regularization technique for dealing with this problem. We split off the ill-conditioned part of A by factoring the convolution operator into A = BG where G is a simple Gaussian smoothing matrix, represented in harmonic-space by Ge = exp I 1 , (3.5) where a = FWHM/V81n2. Substituting the factorization into Equation 3.2, we get G T B T BGs = GTBTd (3.6) BTBx = BTd (3.7) where we are solving for x = Gs so as not to reconstruct the sky at higher resolution 2 http://crd.lbl.gov/~cmc/MADmap/doc/ 24 than that of the instrument. The method that we will employ to solve for x (discussed below) requires a symmetric coefficient matrix. We have constructed the symmetric coefficient matrix by removing the left-most factor, G T , on both sides of Equation 3.7. 3.5.2 Fast convolution on the sphere We now review the relevant formalism for general convolutions on the sphere and refer the reader to [47] for the full details of the fast convolution and transpose convolution of two functions in the spherical harmonic basis. The convolution of a band-limited beam function 6(7) with the sky 3(7) is given by the following integral over all solid angles T ( $ 2 , e , $ ! ) = y'dfy[jD($ 2 ,e,$ 1 )&](7)M7) (3-8) where D is the operator of finite rotations 3 , Db is the rotated beam, and the asterisk denotes complex conjucation. T($ 2 , @, $1) represents the beam-convolved signal for each beam orientation (<3>i,0,$2)Following a geometrically-motivated argument, the evaluation of Equation 3.8 is simplified by factorizing the rotation into two parts Z>(<E>2, e , $1) = D(4>E, eE, o)£>(0, e, u) (3.9) where the various angles (4>E,9E,<fi,9,u) are defined in [47]. Then, Equation 3.8 is rewritten as T(4>E,<I>,U) = f dn^\b{(f>E,eE^)b{<i>,d,u))b]^Ys{fi (3.10) for 9 and 9E fixed. The three-dimensional Fourier transform of T(<fiE, 0, u>) is defined as 1 /*27r Tmm'm" = J^y J d<f>Ed<f>du}T {<j>E, 0 , w )e^B-im'*-im''W) ( 3 > n ) and we find that Tmm'm" , s l m ^ m m ' (^g)^m'm" (^)^m"- = / (3.12) e The D-function, also known as the Wigner D-matrix, is related to the function dl 3 , Our Euler angle convention is denned as active right-hand rotations about the z, y, and z axes by $ 2 , 0 , $ i , respectively 25 by D^ m ,( a ,/?, 7 ) = e - i m a O ( ^ ) e - i m ' 7 . (3-13) The explicit expression for the d-function is given in [68] as dmmW)-^ l ) (e + m-t)l(l-m'-t)\(t + m'-m)\ x(cosp/2)2i+m-m'-2t(sinP/2)2t+m'-m { j (3.15) where the sum is taken over all values of t which lead to non-negative factorials. In practice, the d-matrices are more conveniently computed by making use of their recursion properties, as we do in our algorithm [69]. Following an analogous derivation as for the convolution, the transpose convolution of T($ 2 , Q, $i) is given by y*(7) = y , d<Med$ 1 [D($ 2) e,$ 1 )&]( 7 y:r($ 2 ,e,$ 1 ). (3.16) In spherical harmonic space, this becomes ylm = 2 ^ dmml(0E)dm,mll(6)b}m„Tmm>m». (3.17) m'm" An important feature of our approach is that it economizes the computational effort if the beam is nearly azimuthally symmetric. The parameter of the method that sets the degree to which asymmetries of the beam are taken into account is "i.max, the maximum m" in Equations 3.17 and 3.12. For m max = 0 we recover the computational cost of simple spherical harmonics transforms, (^(^ax)- Since m max is bounded from above by £max, the computational cost of the method never scales worse than C(^max)- For a mildy elliptical beam, we anticipate that just including the m max = 0 and m max = 2 terms will suffice, since the mmax = 1 term vanishes by symmetry. 3.5.3 Binning and simulation Now we will show how the fast convolution techniques of the previous section are implemented in the deconvolution map-making algorithm. For clarity, we rewrite Equation 3.7 in the compact spherical-harmonic basis (summing over repeated in- 26 dices) BVM>mm'm"^'rnmlm"LM^LM = B L ; M'mm'm" -^mrn'm"> (3.18) where B T acting on Tmm>m» is given by Equation 3.17 and B acting on xLM is given by Equation 3.12. To make matters even more concrete, we now explicitly describe the steps required to simulate time-ordered data d from a map — our "simulation" step. First, we convolve the beam bem with the map aem to obtain Tmm>m» following Equation 3.12. Then we inverse Fourier transform the Tmm/m» to get T($ 2 ) 0 , $i). Next, we must account for the scan path ($2(2), ©(*)> $i(0)> where $ 2 and G specify the position on the sphere and <&! specifies the orientation of the beam. This is achieved by extracting those values in T($ 2 , ©, $1) which fall on the scan path of the instrument. As a second example, we describe how to compute the right hand side of Equation 3.7. Start with the TOD d. For each sample in d, the scanning strategy specifies the orientation ( $ 2 , 0 , $1). The sampled temperature is added into the element of an initially empty array which is identical in size and shape to the array which stored T($ 2 , ©, $i)- We have effectively binned the TOD d, according to the position and orientation of the beam on the sky. Let us therefore refer to this operation as "binning". In order to minimize discreteness effects due to the gridded representation of T ( $ 2 , 0 , $ ! ) , more sophisticated interpolation techniques could be implemented. Additionally, the resolution of the grid into which the data is binned may be increased. These additional features will be revisited when we add in polarization measurements in Chapter 4. 3.6 Solving the Deconvolution Equations To obtain the optimal map estimate, we numerically solve the linear system of equations defined in Equation 3.7 for x^m. We have a choice between direct and iterative solution methods. An iterative method is advantageous compared to a direct method (such as Cholesky inversion) if the cost per iteration times the number of iterations required to converge to sufficient accuracy is less than the cost of the direct method. For the problem sizes of current and upcoming CMB missions, where the map contains a number of pixels npix ~ 106-107, direct solution methods would be prohibitive for two reasons. First, the required number of floating point operations scales as r?pix. Second, the amount of space required to store the coefficient matrix and its inverse scales as riiix. Therefore, the direct solution exceeds the capabilities of modern super- 27 computers by several orders of magnitude. For the PLANCK mission, direct solution would require of order 1021 floating point operations and hundreds of Terabytes of random access memory. We therefore advocate using an iterative technique, the Conjugate Gradient (CG) method [70, 71, 72]. The CG method is well suited to this problem, as it solves linear systems with a symmetric positive definite coefficient matrix and has advantageous convergence properties compared to other iterative methods such as the Jacobi method [70, 71, 72]. In order to apply the CG method, we must be able to apply the coefficient matrix on the left hand side of Equation 3.7 to our current guess for the solution x. In order to do so, we simply perform the two operations of "simulation" and "binning" in succession. The fast convolution and transpose convolution algorithms allow computing the action of the coefficient matrix on a map without ever having to store the matrix coefficient in memory. 3.6.1 Preconditioning It is desirable to minimize the number of iterations the CG method requires to converge to a given level of accuracy. This can be done by preconditioning the system of equations. Preconditioning amounts to multiplying on both sides with an approximation of the inverse of the coefficient matrix and solving this modified system. As long as the preconditioner is nonsingular the solution will be the same for the original and the preconditioned systems, but for a well-chosen preconditioner the number of iterations can be reduced significantly. A natural choice of the preconditioning matrix that we used to obtain the results in this paper is the diagonal matrix [diag(A T A)] _1 , which, for a delta-function beam, is just the inverse of the number of observations per pixel. At every iteration we have an approximate solution x for Equation 3.7. We assess convergence by computing the ratio of L2 norms L2[BTd] where Z^fx] = \ / | x • x|. 28 lJ Uj 3.7 Standard map-maker for comparison In order to compare our results from deconvolution map-making to more traditional techniques we also implemented a map-making code that solves the normal equation (Equation 3.2) assuming an azimuthally symmetric beam. In this implementation the observation matrix A becomes the pointing matrix, containing only a single entry on each row corresponding to the direction in which the main beam lobe is pointing at the time of sampling. Standard map-making therefore reconstructs a map which is smoothed by an effective beam whose shape varies as function of position on the map. This variation depends on the scanning strategy. More precisely, at any given position on the estimated map the effective beam shape depends on the various orientations of the beam as it passed through this position during the scan. For uncorrelated noise and an azimuthally symmetric beam the solution of the normal equation is simple to compute: bin the TOD into discrete sky pixels, summing over repeated hits, and dividing through by the number of hits per pixel. Numerical implementations of this algorithm and its generalization to correlated noise have been described in the literature [59, 61, 60, 55]. However, all of these treatments assume azimuthally symmetric beams. For experiments with highly asymmetric beams, and where contamination from the Galaxy is picked up in the sidelobes, we expect that this method will not fare well against our deconvolution method that also removes artifacts due to these optical systematics. In Section 3.9 we compare our deconvolution method to a standard method. We use the same TOD and scan path as for our deconvolution method. For the standard method, the data is binned into pixelized maps, rather than into the T($2, 0, $i) grid. Unless otherwise stated we use the HEALPix pixelization scheme [73] with resolution parameter nside = 64. The angular scale of a pixel is therefore just under 1°. Recall that our regularization method returns a smoothed map with an effective, azimuthally symmetric Gaussian beam. Thus, in order to compare the two methods we must make a similar modification to our standard map-making. We read out the resulting a^m of our standard map (using the HEALPix anaf a s t routine), after the binning step, and modify them in the following way where B( is the beam power spectrum and Ge is given in Equation 3.5. 29 3.8 Test Cases In this section we detail our tests and comparisons of the deconvolution and standard map-making methods. For the purposes of testing our method, we create several mock beam models, b^m. We test three possible beam shapes that break azimuthal symmetry at progressively stronger levels, two scanning patterns, and skies with and without Galactic emission. We test our algorithms on a simulated foreground- and Galaxy-free sky using a standard ACDM power spectrum and simulated spherical harmonic multipoles a(m up to £ = 128. We also use the first-year WMAP Ka-band (33 GHz) temperature map as our true sky containing Galactic emission. The first beam is a simple model of a sidelobe; it is composed of a Gaussian beam of FWHM = 1800' rotated at 90° to another Gaussian beam of FWHM = 180'. Both the main beam and the sidelobe are normalized such that they integrate to one. The second beam models a (somewhat exaggerated) elliptical shape, composed of two identical Gaussian beams with FWHM = 180' whose centers are on both sides of the optical axis, separated by 180'. The third beam is composed of two identical Gaussian beams (FWHM = 180') rotated at 140° from each other; we refer to this as the twobeam model. This case is motivated by the design of the WMAP satellite, which makes a differential measurement from two horns separated by 140° [74]. We set the asymmetry parameter m max for our three cases (sidelobe, elliptical, and two-beam) to 8, 38, and 128, respectively. Following Wandelt & Gorski (2001) [47], we first considered a basic scan path (BSP) in which the beam scans the full sky on rings of constant longitude with no rotation about its outward axis. To be clear, for the case of the sidelobe beam, the smaller beam follows this ringed-scan while the larger beam remains fixed at the equatorial longitude. Similarly, in the two-beam model, one beam follows the ringscan while the other rotates in smaller circles 140° away. The central lobe therefore covers the whole sky, while the offset beam remains within a band of ±50° centered on the Ecliptic plane. The elliptical beam simply follows the ring-scan, and is oriented such that its long axis remains perpendicular to the lines of longitude. A more realistic observational strategy has a beam that revisits locations on the sky in different orientations. Therefore, we model the one-year WMAP scan path followed by one horn. The WMAP scan strategy also covers the full sky and includes a spin modulation of 0.464 revolutions per minute and a spin precession of one revolution per hour [74]. We used a scaled-down model of the WMAP scan in which 30 the spin modulation is 0.00232 revolutions per minute and a step size of about 46 seconds (roughly 562 samples per period). This produces a pattern very similar to the spirograph-type pattern shown in Figure 4 of [74]. We refer to this as the WMAP-like scan path (WSP). The WSP has about six times as many samples as the BSP. For this strategy, the spirograph pattern is followed by the small beam of the sidelobe, the elliptical beam, and both beams of the two-beam model. In the two-beam case, both beams are offset from the spin axis of the satellite to mimic the WMAP scanning geometry. However, it is not differential in nature, since both beams have positive weight. We test each beam (sidelobe, elliptical, and two-beam) with both scanning patterns (BSP and WSP) on a sky without Galactic emission. We refer to these as the six main test cases. In reality, CMB experiments will also pick up a signal from the Galaxy. We use the first-year WMAP Ka-band temperature map, degraded to an nside of 64 and smoothed with a Gaussian beam of FWHM = 180', as our model of the true sky with Galactic emission. For our last test case, we convolve this sky with the sidelobe beam. For each test case, we assume that the beamshapes of the instrument are known and we use the deconvolution method to deconvolve the map with the same beam that was originally used to convolve our true sky. We attempt to recover features in the map corresponding to the smallest scale features of our test beams. We therefore set the width of our regularization kernel, represented by the matrix G in Equations 3.5 and 3.6, to FWHM = 180' in every case. We compare our map estimates to the true sky, smoothed by the regularization kernel. When we refer to the "true" sky in the following we mean this kernel-smoothed input sky. 3.9 Results and Discussion We present the results of the deconvolution algorithm for the six main test cases in the form of residual maps. We compare these results to the results from standard mapmaking by examining their power spectra and by calculating the root mean square (RMS) difference between the estimated and true sky. For the tests that include the Galactic signal we show the actual map estimates. In Figure 3.4 we plot ratios of the power spectra (ACe/Ce) of the residual maps (both standard and deconvolved) and the power spectrum of the input map. The BSP (WSP) results are plotted in the left (right) column. The solid (dashed) lines 31 represent the relative difference in Ci between the deconvolved (standard) map and true sky map. The standard map-making algorithm failed to give meaningful results for the two-beam test. We therefore excluded this case from the plot. For all cases we chose to present the results after a fixed number of iterations to show the impact of scanning strategy and beam pattern on the condition number of the map-making equations. We find that the deconvolution algorithm outperforms standard map-making by orders of magnitude in accuracy. For a fixed number of iterations, the BSP tests performed less well than the WSP tests. The two-beam BSP and, to a lesser extent, the elliptical beam BSP test cases have not converged to sufficient accuracy. There are several possible causes for this behaviour. The BSP leads to an extremely non-uniform sky coverage. Also, the BSP visits each pixel in a narrow range of beam orientations. Further, the number of sky samples is smaller for the BSP case than for the WSP case (as noted in Section 3.8). All of these aspects can contribute to increasing the condition number of the normal equation, which in turn leads to smaller error decay per iteration of the preconditioned CG solver. In Table 1 we summarize the RMS difference between the reconstructed and true sky. The RMS values are computed using the standard deviations the residual and true maps: stdev(residual map) (3.21) RMS = stdev(true map) where residual map = estimated map — true map. (3.22) The RMS values reflect the trends seen in the spectra in Figure 3.4. Beam Sidelobe Elliptical Two-beam BSP Standard 0.257467 0.186262 N/A Deconvolved 0.000216367 0.0213657 0.102778 WSP Standard 0.178828 0.129715 N/A Deconvolved 3.13741e-07 2.92714e-05 1.08579e-06 Table 3.1: Fractional RMS error for each of the six main test cases. The residual maps are shown in Figure 3.5. In order that the scale of the axes on the residual maps are meaningful, we also show the true sky map. Achieving a stably converging iterative solution method for the deconvolution problem is a success of our regularization technique. The convergence of our iterative solver as function of iteration number is plotted in Figure 3.6. 32 r-j CO 2. CO ^ -O C3 o •a 1 CO 2 o r^H CD CD H •-H f-i O CD s£ a £ a <H-H CO 0) - r t o T3 -EH +3 S-H a o « U a r-H o CO CD +3 o a CO r-H a n .a co -EH ^ CD -EJ a ^ co On*" O ~d 3 o - a OH CD 3 o O ^ CD 4 3 OH -^. ^3 "•^^^^ X. • "s X "***• \^x-* < > < -^ CD -E5 a .a 'co CD S-H CD -3 > o CD " EH O OH CO CD .a a a 8 & CO O -EH CD 3 co ^ o 3 ^ E5 ^^> , ^5* > ^ > pH of £ += t-H \1 CO <H-H O -u —i—r-ir—i—i—i—TT—iv"' l — T T — n — r - £ "3 O a a OH CD o3 03 XJ I co <-— ^"" """^ ^^ \y ^ J :} HJ J-H g OH OH CO CQ OT co CO _o '5b 'o/ov CO CD T3 EH - CD 0> X) •s J-H x) EH o CO CD -EH 5P «3 33 co E3 O 05 o +^ M-l S-. CD CD OH CO S-H O O H CO CD !-. o O Sidelabe WSP Sidclobe UHP '»+. *•' -.. A V - 3 . 1 e - 0 3 ^^K^^.---.—~ B ^ 2.7C-05 fiK KHlpUcol WSP KthpUciil BSP <i<MMp< Two-beam BSP ^Hii^Witii^iiliii •M * * Tvjo-beom I73P /.;. Figure 3.5: Residuals after the first 100 iterations are shown in the first three rows. The figures on the left are for the basic scan path, and the ones on the right are for the WMAP-like scan path. First, second, and third rows correspond to the sidelobe, elliptical and two-beam beams, respectively. The true sky is shown in the fourth row. Note that solutions of the two-beam BSP and, to a lesser extent, the elliptical beam BS.P test cases have not converged to sufficient accuracy. We chose to present the results for all cases after a fixed uumber of iterations to show the impact of scanning strategy and beam pattern on the condition number of the map-making equations. 34 BSP 1CT 10 - 6 1 "I 1 1 1 TTTT -I I I L_L I I i I WSP I ~l 1 I 1 10 number of iterations 1 I' I I irr ~1 1—I—I I I I I I 1 1—I—I 1 I I I ' I ' ' 10 number of iterations 100 100 Figure 3.6: Convergence rates of the preconditioned conjugate gradient solver for each test case. The left panel refers to the basic scan path and the right panel to the WMAP-like scan path. The solid lines correspond to the sidelobe beam, dotted lines to the elliptical beam, and dot-dashed to the two-beam model. In order to be able to compare the performance of our method for different beam patterns and scanning strategies we make the deliberate choice of limiting the number of iterations to 100 and comparing the best results obtained up to this point. Since our error estimate continues to drop stably (except in two cases, where we reach the single precision numerical accuracy floor after ~ 25 and ~ 70 iterations) it is clear that the accuracy of the reconstruction can be improved by allowing the system to iterate further, or by choosing a more sophisticated pre-conditioner. Our final test case consists of a model sky with Galaxy emission convolved with the sidelobe beam over the WSP. We present the output maps of both the standard and deconvolution methods in Figure 3.7, where the maps are shown in Ecliptic coordinates. One can see that the standard map contains a distorted image of the Galaxy and that the deconvolved map is virtually identical to the true map. 35 True sky with Galaxy Standard map Dcconvolution method Figure 3.7: Deconvolving the effects of a large sidelobe in simulated observations of the WMAP Ka band niap. using the coarsened WMAP scanning strategy described in the text. The top map is the input sky map, the middle map is the standard map-making result, and the bottom map is the deconvolved result. 36 3.10 Point source tests with elliptical beams At small angular scales, extragalactic radio sources are an important contamination in CMB maps. The characterization and removal of all foregrounds, including point sources, is a critical step in CMB data processing. In this section we will show that point sources can be used as tools to see how the effect of the beam manifests itself in the map. The reconstructed image of a point source in a map essentially shows an image of the effective beam at that location in the sky. We will revisit point sources in Chapter 6 when we examine the variability of identified radio sources in the WMAP data. In this test, we placed 13 artificial point sources on a grid in one quadrant of the sky as shown in Figure 3.8. We then scanned this simulated sky using the WSP and elliptical beam, described in Section 3.8, to generate a set of TOD. Deconvolution map-making was then run on the simulated TOD to generate skies with reconstructed point sources. We were interested to see how the reconstruction of the point sources depended on their positions in the sky. For scans like the WSP, some regions of the sky are revisited by the scan strategy with the beam having a wider range of orientations than other regions. This has the effect of symmetrizing the effective beam profile over those region. We were also interested to examine the effect of reconstruction on the asymmetry parameter m. The results from our first, low-resolution run at nside = 64 are shown in Figure 3.8. It can be seen that the overall residual in the m = 4 case is roughly five times smaller than the m = 0 case. We also find that point sources in regions of the sky which are observed by the beam at a wide range of orientations are reconstructed with smaller errors. We then increased our resolution to nside = 1024 and examined the reconstruction of an individual point source for m = 0 and m = 4. The results are shown in Figure 3.9. 3.11 Conclusions In this chapter, we have presented a deconvolution map-making method for temperature data from scanning CMB telescopes. Our method removes artifacts due to beam asymmetries and far sidelobes. We compare our technique with the standard map-making method and demonstrate that the true sky is recovered with greatly enhanced accuracy via our deconvolution method. Deconvolution map-making recovers 37 Origin**] map r • • • •* ,1' • •0 •'. \ Cs Diffrrence map, m - 0 Difference map, jn=4 Figure 3.8: The top plot shows the distribution of point sources in the upper left quadrant of the CMB sky (note that the upper end of the color scale has been cut-off at 93.6/i/v'; the original scale terminates at 374/xK). The middle and bottom plots show the difference maps between the input and output skies for m = 0 (middle) and m — 4 (bottom). It can be seen that the overall residual in the m = 4 case is roughly live times smaller than the ?n = 0 case. We also find that point sources in regions of the sky that are observed by the beam at a wide range of orientations are reconstructed with smaller errors. 38 fT Figure 3.9: Examing the effect of beam ellipticity and subsequent deconvolution on point sources in the CMB. The actual point source (top) is scanned with an elliptical beam and reconstucted using deconvolution map-making with the asymmetry parameter set to m = 0 (bottom left) and m = 4 (bottom right). 39 features of the CMB sky on the smallest scale of the beam, thereby achieving a form of super-resolution imaging. This extracts more of the information content in CMB data sets. One of the key difficulties encountered in deconvolution problems is that the systems of linear equations we need to solve are very nearly singular. We solve this problem by introducing a regularization method that allows us to solve the systems stably and recover maps at a uniform resolution and with an effective beam that is azimuthally symmetric and has a Gaussian profile. We tested the convergence speed of two particular scanning strategies and found that the WMAP-like scan is superior to the basic scan in both rate of convergence and true-sky recovery. We hypothesize that this is due to the nature of the BSP, where the poles are the location of the only beam crossings and thus receive many more hits than the rest of the sky. In addition, our implementation of the BSP had a smaller number of samples overall than our implementation of the WSP. We have also shown the relevance of this algorithm to the WMAP mission by demonstrating its operation using a WMAP-like scanning strategy and a two-beam model, which, while not differential, resembles the telescope orientations of the WMAP spacecraft. Our results underline the qualities of the WMAP scanning strategy compared to a BSP strategy for deconvolution map-making. Low-resolution and high-resolution tests with point sources and elliptical beams were performed. The results from these tests are another important way to examine the effect of beams on the maps and the effective symmetrization of asymmetric beams due to the scanning strategy. In order to decouple from issues that are not directly related to the optical performance of CMB instruments, we did not consider the effects of noise in our simulations. For a realistic assessment of the performance of our methods on real data this needs to be added. In particular, the choice of scale for the regularization kernel will depend on weighing the benefits of increased resolution against increased high-frequency. Noise issues will be revisited in Chapter 5. In the next chapter, we will extend our description of the deconvolution mapmaking method to include both temperature and polarization measurements. We will compare our results from polarized map-making with those from the PLANCK algorithm development group, using actual simulated PLANCK data. 40 Chapter 4 Polarized map-making for Planck 4.1 Summary of Chapter We describe a maximum likelihood, regularized beam deconvolution map-making algorithm for data from high resolution, polarization sensitive instruments, such as the PLANCK satellite. The resulting algorithm, which we call PReBeaM, is pixel-free and solves for the map directly in spherical harmonic space, avoiding pixelization artifacts. While Fourier methods like ours are expected to work best when applied to smooth, large-scale asymmetric beam systematics (such as far-side lobe effects), we show that our m-truncated spherical harmonic representation of the beam results in negligible reconstruction error — even for m as small as 4 for a polarized elliptically asymmetric beam. We describe a hybrid OpenMP/MPI parallelization scheme which allows us to store and manipulate the time-ordered data from satellite instruments with a typical full-sky scanning strategy. Finally, we apply our technique to noisy data and show that it succeeds in removing visible power spectrum artifacts without generating excess noise on small scales.1 4.2 Introduction As discussed in Chapter 2, the inflationary model of the universe can be tested by searching for B-mode polarization in the CMB radiation. The detectors on PLANCK will scan the full sky measuring the polarization of the CMB. The raw data will be processed into maps and then into power spectra. If present, the tiny primordial B-mode signal will be polluted by foreground contamination that leaks into the detector through the far sidelobes, and it will be mixed with the much larger E-mode signal by asymmetric beams with cross-polar responses. We continue work on the crucial step of processing data into maps, where we will use the deconvolution tool to remove systematic effects due to beam asymmetries. We will model the full PLANCK 1 This chapter contains material previously published in [75]. 41 simulation for a single channel with the largest beam asymmetries and investigate the effect on map-reconstruction due to beams. In the next chapter, we will examine the effect of noise and diffuse foregrounds on map reconstruction. In Section 4.3 we review the PLANCK satellite mission and science goals. In Section 4.4 we describe the deconvolution map-making algorithm for PReBeaM. Some standard methods for polarized map-making are described in Section 4.5, and we also detail the routine that we use as a baseline for comparison for the tests done in this chapter. Fast convolution for polarimetry measurements are described in Section 4.6. Special features that were implemented in our code in order to boost accuracy are explained in Section 4.7. The simulated data and beams are detailed in Section 4.8. We present results in Section 4.9 showing the effectiveness of PReBeaM in removing systematic effects due to beam asymmetry and we discuss computational considerations. We finish the chapter with our conclusions from this study in Section 4.10. 4.3 Planck The PLANCK 2 satellite will be launched in May 2009 and will measure the CMB anisotropics over the full sky in nine frequency bands from 30 GHz to 857 GHz. The resolution of the instrument in CMB-dominated bands will be 5 — 15' and the sensitivity will be A T / T ~ 2xl0~ 6 [76]. The PLANCK satellite is designed to extract essentially all of the information in the primordial temperature anisotropies, and to measure the polarization anisotropies to high accuracy for 2 < £ < 2500. The major deliverables of PLANCK include all-sky CMB maps at each frequency and all-sky foreground component maps of Galactic synchrotron, free-free, and dust emission. In addition, PLANCK aims to determine the parameters that describe the geometry and composition of our universe to about one percent. The scientific performance of PLANCK depends, in part, on the behavior of systematic effects which may distort the signal. One of the most exciting prospects for the upcoming PLANCK satellite is its capability to measure the polarization anisotropies of the CMB over the entire sky in nine frequency channels. The potential rewards from these measurements are many, and include tighter constraints on cosmological parameters, determination of the reionization history of the universe, and detection of signatures left by primordial gravitational waves generated during inflation [76]. 2 http://astro.estec.esa.nl/Planck 42 4.3.1 The map-making challenge for Planck Measurement of the CMB polarization signal presents a great experimental challenge as it is an order of magnitude smaller than the temperature signal, and it is especially susceptible to distortions due to optical systematics and foreground contaminants. If left untreated, leakage from the much stronger temperature signal will contaminate the polarization maps. Maps and spectra will also suffer from leakage from E-mode polarization to B-mode polarization, jeopardizing the potential detection of inflationary B-modes. At the resolution and sensitivity of the next generation of experiments, including the PLANCK mission, studies of primordial non-Gaussianity may also be sensitive to beam-induced systematics. A primary objective of PLANCK is to produce all-sky CMB maps at each frequency. The process by which the satellite's time-ordered data (TOD) is wrapped back on to the sphere to create an image is known as map-making. The map-making process becomes difficult due to a number of challenges: distortions in the beam, foreground contamination through far-side lobes, size of the data, and correlated noise effects. It is of critical importance to fully characterize the beam, and use this information during map-making to deconvolve beam effects. We will be most interested in studying the map-making problem for the 30 GHz channel of PLANCK, as this will be the most asymmetric beam. Within the PLANCK collaboration, the CTP working group has developed five map-making methods and compared their results using the simulated 30 GHz data in what is known as the Trieste paper [5]. The Trieste paper assessed the impact of beam asymmetries on the PLANCK spectra without attempting to treat the problem of beam asymmetry at the map-making level (an angular power spectrum correction method was developed based on simplifying assumptions). We will revisit some standard map-making methods in Section 4.5. 4.4 P R e B e a M Method We have previously described, in Chapter 3, a powerful map-making algorithm that implements the beam deconvolution technique for temperature measurements. In this section, we will extend that description to include polarization measurements. We refer to this new technique as PReBeaM: Polarized Regularized Beam deconvolution Map-making. While we focus on reconstructing the map with a uniform effective beam and realize corrections to the power spectrum as a consequence, other works 43 by [77] and [78] have focused on deriving corrections to the power spectrum due to asymmetric (noncircular) beam effects. First, we review the standard setup to the map-making problem for a solution of the least-squares (or maximum-likelihood) type. The TOD generated by a detector is effectively a convolution of the true CMB sky with a beam function. If we consider the sky as a pixelized vector, it will have length nPiXei x npoi where npoi = 3 for the / (total intensity), Q, and U Stokes components. The nToD-length TOD vector d is the result of a matrix multiplication of the observation matrix A with the sky s plus the noise n d = As + n. (4.1) In our implementation of the maximum-likelihood solution, we refer to A as the convolution operator. A encodes information about both the scanning strategy and the optics of the scanning instrument. The least-squares estimate of the true sky, s, is given by the normal equation ATAs = ATd (4.2) where A T is the transpose convolution operator. Equation 4.2 holds if the noise is stationary and uncorrelated in the time-ordered domain. The generalization to nonwhite noise is as follows ATN"1As = A ^ - M (4.3) where N is a noise covariance matrix. In this chapter we consider CMB only and we will examine CMB plus white noise in chapter 5. We modify the normal equation by introducing a regularization technique in order cope with the ill-conditioned nature of the coefficient matrix A T A. We split off the ill-conditioned part of A by factoring it into two parts: A = BG. The factor G is what we refer to as the regularizer in PReBeaM. In our study, we choose G to be a Gaussian smoothing matrix, defined in harmonic space as O- = exp ( ^ | ± I > ) or = exp (-w+D-4)) 44 m where a = FWHM/\/81n2. The superscript / refers to the intensity, and G and C refer to the gradient and curl components in the typical linear polarization decomposition. In general, the width of the regularizer can be set so as to reconstruct the sky at any target resolution. In practice, one would not want to choose a regularizer that is smaller than or close to the size of the sample spacing as this would likely introduce sampling effects into the map. We use the width of the angle-averaged detector beam and suggest this as a rule-of-thumb for choosing the regularizer's width. The choice of a regularizer will certainly affect the noise properties of the reconstructed map and power spectrum as it acts to smooth away noise on small scales. PreBeaM estimates the beam-corrected and regularized atm and the map is simply a visualization of this anm solution. If the sampling points were distributed in a way that aliased two different a^m, then the reconstruction would be singular. Our regularization scheme ensures that we are robust to the details of sampling as long as the sample space is finer than the scale set by the width of the regularizer. Our modified normal equation becomes BTBx = BTd (4.5) where we are solving for x = Gs. In this way, we are attempting to construct an image of the sky at the angular resolution of the chosen regularizer. A complete characterization of the beam includes both the main beam and the farside lobes. Sidelobes are located as far away as 90° from the main focal plane beam, and therefore require a large mmax parameter for a complete harmonic description. In [48] we demonstrated the full potential of our method using far-side lobes and maps with foreground signals. Here, we show the usefulness of PReBeaM for deconvolving main-beam distortions. In fact, we find that it makes sense to use PReBeaM for main beam effects since only a small mmax parameter is needed to capture the azimuthal structure of the main beam. In this way, we profit from the computational advantage of our method in the case of small mmax, allowing for the unified treatment of main beam and side lobe effects. 45 4.5 Comparison to standard map-making and other map-making methods In standard pixel-based optimal mapmaking (see, for example, [63, 79, 55, 58, 56]), in which one assumes that the observing beam is spherically symmetric, A is a sparsely filled pointing matrix. In the case of a single dish experiment, each row of A contains only three nonzero elements. The deconvolution map-making approach does not assume spherically symmetric beams, instead allowing for arbitrary beam shapes. We achieve this added complexity primarily by solving the normal equation in spherical harmonic space in order to make use of fast and exact algorithms for the convolution and transpose convolution of two arbitrary functions on the sphere [47, 80]. These algorithms are described in abbreviated form in Section 4.6. In addition to PReBeaM, another deconvolution map-making technique for PLANCK has been established by Harrison et al. [81]. Both methods allow for arbitrary beam shapes, and, in both cases, the asymmetry of the beam is parameterized by an asymmetry parameter mmax that can vary between 0 and £max- Our method scales computationally as 0 ( ^ a a . m m a x ) ; this is advantageous when large gains in accuracy can be achieved with small increases in mmax. In contrast, the Harrison approximate method scales as 0 ( ^ a 2 . ) thereby incurring a fixed computational expense for arbitrarily large mrnax and effectively setting a limit to the maximum £ at which the analysis can be done. Excluding TOD operations we achieve a speed-up of imax/f^max, which is of O(l02) for the case we examine herein. The Harrison method takes advantage of the PLANCK scanning strategy to condense the full TOD into phase-binned rings, thereby achieving a significant reduction in processing time. A secondary advantage of operating entirely in harmonic space is that artifacts due to pixelization (such as uneven sampling of the pixel) are completely avoided. Issues with undersampled pixels can be problematic for a pixel-based method — as in the case of pixels with less than three non-degenerate observations — and can result in pixels which are excluded from the map-making completely. Our approach is less sensitive to degenerate pixels resulting from poor sampling and does not result in excised pixels. The assumption of a band-limited signal, justified by finite resolution, regularizes the solution for degenerate pixels based on neighboring observations in the time stream. We demonstrate the difference between a pixel and harmonic-based method by comparing a PReBeaM temperature map at the Healpix [73] resolution of nside 1024 with a binned noiseless temperature map at the same resolution. A binning of the TOD from the 30GHz channel into an nside 1024 pixel map results 46 Binned Temperature Map PReBeaM Temperature Map *H # (8.0, 23.0) Ecliptic (6.0, 23.0) Ecliptic Figure 4.1: Comparison of one segment of the sky from a binned noiseless map (left) and the PReBeaM temperature map (right), both at a Healpix resolution of nside 1024. At this resolution, the binned map contains a number of unobserved pixels, of which three are visible in this frame. The PReBeaM map contains no unobserved pixels. in a number of unobserved pixels. The analogous PReBeaM map at nside 1024 contains no unobserved pixels. Figure 4.1 shows identical regions of the sky from the binned noiseless map and from PReBeaM. The difference between pixel and harmonic methods is even more distinct if we consider polarization maps. We quote results from [82] in which polarization maps were made from the 30GHz data at nside 1024. In section 3.7 of [82] it is noted that the nside 1024 binned polarization maps have 912,968 unobserved pixels (where they reject not only missed pixels, but also those with a condition number less than 0.01 — which indicates a degenerate pixel). In contrast, a PReBeaM polarization map at the same resolution has no missing pixels. While it is true that our approach would suffer in the case of gross undersampling or huge chunks of missing sky, we are not expecting to deal with these sorts of issues with CMB satellite missions such as PLANCK. Our method does have its own unique sampling issues that can be understood by examining the problem in harmonic space rather than in pixel space. The beam a^m used to simulate the TOD set a natural band limit where they approach zero. We choose an £max cutoff at which to reconstruct the sky aem. If we choose £max much less than the beam band limit then we may introduce small-scale features such as ringing 47 (of particular importance when we include foregrounds such as point sources). If we choose £max much larger than the £ where the beam effectively bandlimits the signal, then all agm for £ higher than that are noise dominated, and the regularized solution will therefore be driven to 0 for those large £. 4.6 Fast all-sky convolution for polarimetry measurements For a full presentation of the formalism for convolution of an instrument beam with a sky signal, the reader is referred to [80]. In compact spherical harmonic basis, Equation (4.2) is written as •^•L'M'mm'm"J^mm'm"LMsLM = AL,M'mm'm"^mm'm" (4-6) where SLM is the spherical harmonic representation of the sky and Tmmim" is defined as the result of a convolution of a band-limited function b with the sky s. The PLANCK Level-S software [83] nomenclature refers to T mm ' m " as a ring set. This is written in harmonic space as T -X^t^J h1* -L e G hG* S •l-mm'm" — / , l'2p f a % ' "•" btm°e.M' t +s?mb?M>)dernM(9E)deMM,(8) (4.7) where (6E, 0) are fixed parameters that define the scanning geometry. In Equation 4.7, demM(6E) and deMM,(8) are related to the Wigner D-matrices by BL'mfa 0, VO = e - ^ m ^ e - * " * - (4-8) Analogously, the transpose convolution in harmonic space is given by Vem = 2 ^ dmm/(0£)dm,m,,(0)6£7£,,Tmm/m// m'm" where P = I,G, C. 48 (4.9) 4.7 PReBeaM Implementation Now we outline the algorithmic steps taken to make a map from a TOD vector by PReBeaM. First we construct the right-hand side of Equation 4.6 in two steps: converting TOD to an Tmm'm" array and applying A T . Tmm'm" is constructed by transpose interpolating the TOD vector d. The transpose interpolation of the TOD vector onto the Tmm'm>i grid is akin to a binning step, where each element of the TOD is mapped, via interpolation weights, to several elements of the Tmm'm" cube according to the orientation and position of that data point in the scanning strategy. The interpolation scheme is described in greater detail in Section 4.7.1. Next, we transpose convolve the beam coefficients b(m with Tmm'm" according to Equation 4.9. Once the right-hand side has been computed, we use the conjugate gradient iterative method to solve equation (4.2). With each iteration, the coefficient matrix A T A is applied using the following procedure: 1. Apply the convolution operator, A, to project the sky a^m onto the convolution grid Tmm'm". 2. Inverse Fourier transform over the first two indices of Tmm'm" to get T<s>2,Q,m» (we omit the transform over m" as it is incorporated in the interpolation scheme). 3. Forward interpolate from T$2)eim» to a TOD vector. 4. Transpose interpolate from the TOD vector to a new ring set T$ 5. Fourier transform over the first two indices of T$2 e r n „ to get e „. Tmm,m„. 6. Apply the transpose convolution operator, A T , to project the ring set Tmrn,m„ back into a new sky a^m vector. 4.7.1 Polynomial Interpolation and Zero-Padding PReBeaM uses the same polynomial interpolation as implemented in the Level-S software used to generate the simulation TODs as described in [83]. The objective of forward interpolation is to construct a TOD element at a particular co-latitude, longitude, and beam orientation using several elements of the ring set T and their corresponding weights. Transpose interpolation operates in exactly the opposite manner as the forward interpolation: distributing a single element in the TOD to multiple 49 elements of the ring set. This is done using the same weights calculated for the forward interpolation. The entire operation of interpolation and transpose interpolation from ring set to TOD and back again is depicted in Figure 4.2 Ring Set • / / / _ ^ A_' / V / / , / / / // / / / / / / / , / / / / Ring Set 7 Transpose Interpolation •1 •/ ;j~~ \ - y mr > \ \ \ \ vs \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Forward J9 \ s :i•• k, \ \ k / i >i i Figure 4.2: Forward interpolation from ring set to TOD element and transpose interpolation from TOD element to ring set. PReBeaM also includes the option to zero-pad during the Fast Fourier transform (FFT) and inverse FFT steps. This means that the working array (either Tmm'm» or T$2ie,m») is enlarged and padded with zeros out to H.max,pad > (-max- This has the effect of decreasing the sampling interval. We found that the combined effects of small-order polynomial interpolation (order 1 or 3) and zero-padding of 2 x £max or 4 x £max dramatically reduced the residuals in our maps. 4.7.2 Parallelization Description PReBeaM employs a hierarchical parallelization scheme using both shared-memory (OpenMP) and distributed-memory (MPI) types of parallelization. The map-making was performed on the NERSC computer Bassi3. Bassi processors are distributed among compute nodes, with eight processors per node. OpenMP tasks occur within a node and MPI tasks occur between nodes. We show a diagram of our hybrid parallelization scheme in Figure 4.3. The full TOD and pointings are divided equally between the nodes for input and storage of pointings. Within an iteration loop, four head nodes are designated to perform the convolutions. All other tasks are shared by all nodes, including the head nodes. Since convolution is minimally time consuming, we are not concerned with the small amount of time the nonhead nodes spend idling while convolution is performed. The number of head nodes needed corresponds to the number of detectors; for the low-frequency 3 https://www.nersc.gov/nusers/systems/bassi/ 50 instrument (LFI) 30GHz channel, this is four. Each of these four nodes performs the convolution of the sky with one of the four detectors. The resulting arrays are then distributed to all nodes for interpolation over the segment of data stored locally and then gathered back onto the designated nodes for the transpose convolution. Finally, the aim are summed, using the MPI task mpi_reduce, into a single aem on a single node; this is the new estimate for the sky vector. This particular scheme was devised so that the four distinct convolutions that must occur (one sky with four different beams) can take place simultaneously, while the pointings are distributed among as many nodes as possible for maximum speed in interpolation. Both convolution and interpolation and their transpose operations make use of all processors on a node by using OpenMP directives. 4.8 Simulations and Beams The simulated PLANCK data on which PReBeaM was run were generated by the PLANCK CTP working group for the study of the performance and accuracy of five map-making codes summarized in the Trieste paper [5]. The 30 GHz channel is composed of two horns, labelled "27" and "28", each of which is comprised of two detectors tuned to orthogonal linear polarizations as shown in Figure 4.4. PLANCK will spin at a rate of approximately 1 rpm, with an angle between the spin axis and the optical axis of ~ 85°. We used the cycloidal scan strategy in which the spin axis follows a circular path with a period of six months, and the angle between the spin axis and the anti-Sun direction is 7 ? 5. TODs were generated for 366 days for the four 30 GHz LFI detectors. At a sampling frequency of 32.5 Hz, this corresponds to 1.028 x 109 samples per detector, for a total of over 65 Gb of data and pointings. The simulated data also included the effects of variable spin velocity and nutation (the option to include the effects of a finite sampling period was not included). The data were simulated with elliptical beams having a geometric mean FWHM of 32' 2352 and 32' 1377, and an ellipticity (maximum FWHM divided by minimum FWHM) of 1.3562 and 1.3929 for each pair of horns. The full details of the four detectors which make up the 30 GHz channel are given in Table 4.8 and contour plots of the main beam co-polar components are shown in Figure 4.5. The contour beams are plotted in a uv-plane (described in [5]) where u = sin(9)cos(p), v = sin(6)sin((p) and 9 and ip are the polar and azimuth angles of the spherical polar coordinates of the beam coordinate system. The quantity ifjpoi is the angle of the polarization measurement, defined with the beam in a fiducial orientation. The orientation of the 51 TOD convolution with beam 1 convolution with beam 2 mpi broadcast of ring array interpolation and transpose interpolation on TOD segment interpolation and transpose interpolation on TOD segment mpi reduction of ring array transpose convolution with beam 2 transpose convolution with beam 1 mpi reduction of aim Figure 4.3: Depiction of hybrid parallelization scheme used in PReBeaM. Rectangles represent work done on a node, ellipses represent data products, and arrows represent transfer of data. The work done within a node (convolution, interpolation and their transpose operations) is parallelized using OpenMP. This figure shows a slice of two head nodes, however, the algorithm may operate on many more nodes. 52 LFI Main B e a m s , , , . ! , , , 1 25 0.05 i i *"-"- J 22 j S C A N -24 D I R E C T I O N V ^ ^ 2 0 18 1 » 1 iK "I U •—| 0.00 ^} i - i i 1 _28 1 0.05 1' 26 i i i 1 i -0.05 i i i i 0.00 i J ... i i , ,L 1 0.05 Figure 4.4: Plot showing the footprint of the LFI main beams on PLANCK as seen along the optical axis looking towards the satellite. The pair of horns comprising the 30 GHz channel are labeled 27 and 28 in the plot. The scan direction and polarization orientations for the co-polar beams are also depicted. Figure reproduced from [5] with permission. 53 Figure 4.5: Contour plots in the uu-plane of the PLANCK 30 GHz main beam detectors. This coordinate system permits the beam to be mapped from a spherical surface to a plane. From left to right, the beams are: LFI-27a, LFI-27b, LFI-28a, and LFI-28b. Figure reproduced from [5] with permission. beam major axis is given by ipeu. The widths and orientations of the beams were different; this was referred to as beam mismatch in the Trieste paper. In spherical harmonic space, the simulation beams were described up to a beam mmax of 14. The same beams were used in PReBeaM to solve for the map, although we allowed the beam asymmetry parameter mmax to vary. Detector 27a 27b 28a 28b FWHM 32f2352 32^1377 32^2352 32:1377 Ellipticity 1.3562 1.3929 1.3562 1.3929 fj'pol tpell 0?2 89?9 -0?2 -89?9 101?68 100?89 78? 32 79? 11 Table 4.1: The four detectors of the 30GHz PLANCK channel and their FWHM, ellipticity, angle of polarization measurement Vv> a n d orientation of beam major axis ijjeii. 4.9 Results and Discussion For our tests, we make temperature and polarization maps from simulated one-year observations of the four 30 GHz detectors of the PLANCK LFI. We examine the case of a CMB signal only and relegate the case of CMB plus uncorrected (white) noise and foregrounds to Chapter 5. The data from the 30 GHz channel was an optimal choice for this analysis. These receivers sit farthest away from the center of the focal plane and therefore have the strongest main beam ellipticities of any PLANCK channel, while the low sampling rate and resolution minimize the data volume. PReBeaM operates entirely in harmonic space, solving for and producing as output, a set of a(m. For visualization purposes, maps were made from a^m's out to 54 Pmax 512 and at the Healpix resolution of nside 512 (~ 7" pixel size). Most of the results presented here were attained with an FFT zero-padding of factor 4, an interpolation order of 3, and an asymmetry parameter of m max = 4 (we will note where the parameters differ from this). To compare with the input signal, a reference map representing the true sky was created by smoothing the input a^m with a Gaussian beam of FWHM = 32/1865. Similarly, our regularizer G (in Equation 4.4) was set to have a FWHM of 32/1865 to match this smoothing. As noted in Section 4.8, the same data we use here have been processed by five map-making codes in [5]. We have chosen to compare our results with the analogous results from Springtide, one of the codes in this study. Springtide was chosen, out of the five codes, because it is the map-making code installed in and used by the PLANCK Data Processing Centers for the HFI and LFI instruments. It is sufficient to compare with Springtide only as no significant differences in accuracy were found between codes (with similar baselines and in the absence of noise; [5]). Springtide is a destriping algorithm designed to remove low-frequency correlated noise from the TOD by fitting for, and subtracting, offsets. While Springtide does fit for offsets even if the noise is absent or uncorrelated, the results in this case are not significantly different from that which one would derive from a simple binning algorithm. Because we are using Springtide to represent all non-beam-deconvolution methods we will refer to the Springtide maps as the binned maps. We begin by examining the spectra in the binned map, PReBeaM map and the smoothed input map shown in Figure 4.6. The effect of the beam mismatch is clearly seen where the peaks and valleys of the binned map spectra have been shifted toward higher multipoles. The detectors measure different Stokes / parameters, which translates to artifacts in the polarization map. Deconvolution suppresses leakage from temperature to polarization as evidenced by the PReBeaM spectra that overlays the input spectra. This shift is expected to remain apparent in the TE spectra of nonbeam-deconvolved maps even in the presence of noise, because of a larger temperature signal and the temperature-to-polarization cross-coupling. The fractional difference in the angular power spectrum (denned as (Ceout — Cein)/Cein) of the input and output maps is shown in Figure 4.7. We show the fractional difference spectra for the TT, EE, and cross-correlation TE signals, omitting the BB spectra since CfB is zero in the simulation of the CMB map. The results for PReBeaM are shown at three intervals: the 25th, 50th and 75th iterations. This shows the behavior of the power spectra as PReBeaM converges on the solution. The beam mismatch effect is also seen in Figure 4.7, where the fractional difference in the 55 G3 <3 e-S P O CO QJ z3 a; ss u^i b ^ "'" ' ' ' i""'' ' V'yr'jr f : Q. •£ • £ / a £ 5o "o o c oa E = . Cfl ID I ffl S )m a Q. o « -3, •- _ - a3 : \ f : - ( O co [i» r f ] 3 3 kx*] 'o CO C a> < o tc m 8 • « * 56 PReBeaM spectra lies closer to zero than the binned map spectra over the full range of multipole moments for EE and TE. Next, we examine the residuals by looking at the power spectra of the difference maps and comparing that with the power spectra of the difference in a^m's. The difference in aem can be calculated for PReBeaM as we have the option to output either a^m or maps, while the comparison method produces only maps. Figure 4.8 shows the power spectra of the difference map for the binned map and PReBeaM map. We overplot the power spectra calculated from the difference in a^s {flimout ~~ atmin) for PReBeaM. It can be seen that the spectra from the difference in a^m is smaller on small scales than the spectra of the difference maps. This is likely due to inaccuracies introduced during the transformation from pixel-map to afm-space to Q's. By-passing this extra step is another advantage of a harmonic-based method. As described earlier, PReBeaM allows for variation in the asymmetry parameter iTT'max- We examined the performance of PReBeaM as a function of mmax, setting it to 0, 2, 4, 6, and 8. In conjunction with this, we plot the beam a^m coefficients as a function of the first few m in Figure 4.9. A remarkable improvement in the power spectra was found by increasing mmax from 2 to 4, while an increase from 4 to 6 only resulted in marginal improvements. This effect is best seen in the BB power spectra as shown in Figure 4.10. There is a simple explanation for the fact that the PReBeaM a^m become very close to the input a^m at mmax 4, and thus do not change much beyond mmax 4: there is an order of magnitude drop in the beam aem coefficients from mmax = 4 to mmax = 6 for £ = 0 to £ — 512 as seen in Figure 4.9. Thus, while the input TOD was simulated with a beam having an mmax cutoff of 14, PReBeaM operates optimally at an mmax of just 4, thereby allowing us to capitalize on the computational property that PReBeaM scales as mmax. We define a quantity called na that models the expected bias due to beam asymmetry systematics to the x 2 statistic nai = £ \ £/=2 ^ . ,4,0, Planck^/ We use na to quantify the maximum, or worst-case bias beam systematics could induce in a cosmological parameter that happened to be degenerate with that parameter. The quantity apianck is the expected la errors for the LFI 30 GHz channel, computed as the diagonal elements of the covariance matrix for the simulated input spectra, assuming a sky fraction of 0.65. The na values are plotted in Figure 4.11 and show 57 lack nee. ,0 "•—*s P< CD M s—* £ a > £ o •d CD cj a> (3 r* (-1 -»-2 ,0 CO 0) X! O rH c -^» CaD -d -d o ves ^—^ _o S-( 3 o cya a CO -tj cS S-i cu _-t-3 o S-I T3 S3 a CD^ CD -3 a3 ?3 ^ip c S 3 o S3 - d " .2 <u : 4^> s-< o ^ S3 'o/'cw Cu CD PQ CD ff} OH f»4 £ CO 03 £ !g O ^3 CO CD a 53 CJ CD eg CD CO ^ % cu f ^ •S 2 r CD CJ 1%- CJ CD a a eg a -"! ^ eg s3 o 'o/ov Va .. H ^ .2 CD S-i y—, CD Go £ PH 58 CJ > CQ o -S CO jM'"1] 20V * ^ , » " "OV > o OH CD 1—( fa 59 a 3 0 OO 03 13 xi +^ CD S late ect ra, [ z »tf] ?3V [j»^] OV 0 3 IS 0 Cu CO o -a i co ^-^ e <o — 1 (''/"/ / / / ft < / /' / /' ; i/ i coef ficien .7 1 i i i / 1 / 1 i i 1 i i. r* i i i CD CQ = m +mJ - O ^ a, -v _ • II S +J ,—^ o T ) OH CD •.00 4-J \ \ \ \ \ • • • - » \ i <M ' O —i— i •7 ' •7 1 i F o 1 i .7 d ' nen / 1/) ' shec ince 1 -' O 3CM \ \ X \\ \\ \ .\ \\ \\ II '"v. B <") W OQ a> -^ •-• .^ • ~ i r--, o £ i 1 V: , P 1 : /' / / / 1: i i •' • K|* I CO Pi i 1 coefficients /• / 1 1: ' i » ! : I i u s E o o/ *M: 1 : \ « • \ \ \ \ 8 ^ ! i i ID i CD CD O \ -T-J .2 a £ ! i ! // ' / I I I I I I i 1 i • i ' : 1 \ CO \ '. \ S \ \ CD r3 \ c3 V. ,-—. O T3 > J. CD O CO ^ H 3 l>l en ° ° ^ II CD fa 60 e CO BB Spectra of CMB Map PReBeaM, mmax=2 ^A^^^^V-tyK,, PReBeaM, mmax=4 10"6b PReBeaM, mmax=6 JL 10"7 /WvV i A- u w V/ ur v/ *3fl io- B b |H %iMvllh^AvM^^^^'^ 10" 100 200 300 400 500 Figure 4.10: BB power spectra as a function of asymmetry parameter mmax for mm(tx = 2 (blue curve) , 4 (cyan curve), and 6 (red curve). The input BB spectra was zero so the smallest output BB spectra is most desirable. In this run, the PReBeaM input parameters interpolation order and zero-padding were set to 1 and 2, respectively. 61 that PReBeaM reduces the worst case bias due to untreated beam systematics by 1 or 2 orders of magnitude. We examine the resulting temperature and polarization (Q and U) maps. The output map for both PReBeaM and the binned map was subtracted from the smoothed input map at the same resolution to make the residual maps shown in Figure 4.12. PReBeaM residuals were plotted on the same color scale as the binned map, showing that PReBeaM attained smaller residuals for both temperature and polarization. 4.9.1 Computational Considerations The computational costs and advantages of our method should be noted. To perform a convolution up to £max requires 0(£^ o:r ra max ) for the general case. Since mmax is bounded by (,max, the cost never scales worse than C?(^ aa .) and is only 0(1^^) for the symmetric beam case. By comparison, a brute force computation in pixel space would require 0 ( ^ a x ) . In this study, data were simulated with beams having an asymmetry parameter of mmax = 14, but maps were made using a cutoff value of mmax — 4 in PReBeaM. We have demonstrated that computational cost can be minimized while still achieving the benefits of beam deconvolution It was found that an increase in the zero-padding factor from two to four produced superior results over an increase in the interpolation order from one to three. An optimal run of PReBeaM will therefore include the largest zero-padding possible, given machine memory constraints, in conjunction with a polynomial interpolation of order one or three. This is advantageous since the time spent in an FFT is nearly negligible and affected only minimally with an increase in zero-padding. In contrast, time for interpolation scales as the interpolation-order-squared, and, as this is a TODhandling step, it dominates over any cost incurred by convolutions. In the case of the results shown here, interpolation steps consume more than 90% of the wall-clock time per iteration. The results produced here were generated using 12 nodes on the NERSC computer Bassi (making use of all eight processors per node) and was complete in about 29 wallclock hours, for a total of 2797-CPU hours. The maximum task memory was 20 GB on a single node. 62 T3 v -^ o a3 1 1 ' i ' • i", 1 ' • ' ' n ~i—i : 1Y11 i i i—i—i 1111 i i i i , i 1¥11 I i-r-i—' O CQ CD 3 CO a3 OH 0) Xi \ Spectra D o a> O, L. y <D o. in : 1 c UJ h- c en * " « ; - \ \ a> v> o (J : \v \ o 5 • ^s^^ N . 1 ^ — ' • 1 • ' ' CD \ v N. " ^"-*«^^ J. • . S3 t *o3 i CO CD T2 v i rTT^trrpi^-r-: CD a JIIIII i i i—-* o ••p-n—i—i- i „ T -pn1! r- a i i i i U 2 ft < 03 J a &• CD d ^ cp 03 ft 1-° ^ t o3 'bO X I O -^ i Spectra -JMIII n—i cc3 to CO (-1 CD D . iiim <D . -o •. o inn 11 i—i co 0) 2 ^v l i i i 1 Nsir—*_ "ft \ ^~X \ • • ' en 1. " ' I _o Xl \ (/> V ~ \ \ \ o CD : \l \ \ m c Worst Case Bias Worst Case Bias a3 1 l CO O o * \ ' \ \ - '• ' ' \ \ '• \ • 2£ O CO CP o3 co C ~ is * - § .23 ^ g CO CD CD CO S3 i s-i CO 03 * , * \ O ft S i ^ \ CD CO <D 03 '' \ o •4-C i • Sis S-l CO <P r* ° ^ 3 +* ^ .. T3 ~cp ^^ 1 CD -—^ 2 o ? N fa 8 o 63 PKeifruM Residual T Sprinylide Kesidua! T PKI'IIPPM K M d n a l Q Springtide Residue! Q I'Hrijt'dM Residuwl U Springtide Residual U Figure 4.12: The residuals between the input reference sky and PReBeaM output (left column) and the residual between the input reference sky and the binned map (right column) for Temperature (T), and the Stokes Q, and U parameters. 64 4.10 Conclusion We have found that PReBeaM has outperformed the standard binned noiseless map using two measures: spectra and residual maps. We examined the fractional differences in the spectra and found markedly smaller differences in the PReBeaM spectra versus the binned map spectra across a range of multipole moments. We find that map-making codes that do not deconvolve beam asymmetries lead to significant systematics in the polarization power spectra measurements. The temperature-topolarization cross-coupling due to beam asymmetries is manifested as shifts in the peaks and valleys of the spectra. These shifts are absent from the PReBeaM spectra. We translated the errors found in the power spectra to an estimate of the statistical significance of the errors in a parameter estimation resulting from these spectra, which we call na. This analysis showed that the worst case parameter bias due to beam-induced power spectrum systematics could be tens of sigma while PReBeaM reduces the risk of parameter bias due to beam systematics to much less than la. We also found the /, Q, and U component residual maps to be smaller for PReBeaM than for the binned map, implying smaller map-making errors. We have discussed some advantages of using a harmonic-based map-making routine such as the avoidance of artifacts due to pixelization. An additional feature of our method is the direct generation of sky a^m. Analysis of CMB data often requires calculation of the power spectrum or bispectrum from the sky aem. A pixel-based map-making method generates a map which must first be transformed into a set of aim before computing the spectra. This extra step may introduce inaccuracies which will carry forward in the analysis done on the spectra. We have presented here the first results from PReBeaM for the straightforward test case of CMB only, and including only the effects of beams in the main focal plane. However, there is great potential for using PReBeaM to remove or assess systematics due to the combination of foregrounds and beam side lobes. Systematics introduced by side lobes will appear on the largest scales, potentially impeding the detection of primordial B-modes on the scales where they are most likely to be measured. We have already shown for temperature measurements [48] that our deconvolution technique can be used to remove effects due to side lobes. In the next chapter, we will examine the noise properties of PReBeaM maps and will include foregrounds from extragalactic sources and diffuse Galactic emission. 65 Chapter 5 Noise and Diffuse Foregrounds In the previous chapter, using the ideal case of CMB signal only, we showed the results from our polarized map-making algorithm. In this chapter we will examine the noise properties of CMB maps, and we will include foregrounds from diffuse Galactic emission. This chapter is organized as follows. In §5.1 we review some of the general noise properties of PLANCK. In §5.2 we test PReBeaM on CMB plus white noise and contrast these results with Springtide, our standard of comparison. Some of the diffuse Galactic foreground emissions that are expected to contaminate CMB measurements are described in §5.3. Results from a PReBeaM run on a map with only polarized foregrounds are presented in §5.4, where we also look at a temperature-to-polarization leakage case. We summarize our conclusions for this chapter on noise and foreground issues in §5.5. 5.1 Noise model for Planck Instrumental effects include not only beam effects but also noise, which can be white (uncorrelated) noise or correlated (1//) noise. White noise refers to random, uncorrected instrument noise. The expected white noise contribution to the RMS total error is proportional to \J(l/N0bs). We use the white noise from the Trieste [5] round of simulations. Its nominal white noise standard deviation per time sample integration time is a = 1350/J,K (CMB scale). Correlated 1 / / noise is characterized by a power spectrum that rises at low temporal frequencies. We test PReBeaM on white noise only, omitting correlated noise and noise due to temperature fluctuations in the instrument cooling system. We refer the reader to [5] for more details about 1 / / and cooler noise and also for tests that were done on other map-making routines with the full PLANCK noise model. 66 5.2 Testing PReBeaM on CMB + white noise We run PReBeaM on TOD containing CMB signal and white noise and compare the results to the smoothed input CMB spectrum and the analogous results from Springtide (in this case we refer to Springtide directly since this is not simply a binned map). The level of the uncorrelated noise is specified in the detector database, and has a nominal standard deviation per sample time of a — 1350/LJK [5]. PReBeaM achieves a noticeably superior fit to the input spectrum compared with Springtide from t ~ 150 to ~ 250. We acknowledge that the reconstructed power spectrum is sensitive to the choice of regularizer. Without exploring the full parameter space for this variable, we found that our choice for the regularizer fortuitously led to a recovered power spectrum free of excess noise on small scales. Assessing the relative performance of PReBeaM and Springtide in more detail would require performing Monte Carlo averages. We focus on the TE spectrum, shown in Figure 5.1, since the improvement is visible even for a single simulation. For the other spectra, PReBeaM performs as least as well as Springtide, but the detailed difference are more difficult to assess without a Monte Carlo study. 5.3 Knowledge of polarized foregrounds Foreground emission fluctuations are a limiting factor for precise CMB temperature and polarization anisotropy measurements. Foregrounds include any emissions that obscure the primordial CMB signal after last scattering. These may include extragalactic discrete point sources, synchrotron emission, dust emission, free-free (or Bremsstrahlung) emission, spinning dust emission, and atmospheric emission. Among these sources, synchrotron and dust emission are expected to be significantly polarized. The level at which these foregrounds contaminate the CMB signal depends on both the observing frequency and angular scale. In this section, we review what is currently known about synchrotron and polarized dust emission, and also highlight areas of limited understanding. At 30 GHz, a dominant source of extragalactic emission is radio point sources. We will study some properties of these sources in the next chapter. 5.3.1 Synchrotron emission Highly polarized synchrotron emission is the result of cosmic-ray electrons accelerated by the magnetic field of the Galaxy (e.g. [84]). Synchrotron emission is the main 67 TE S p e c t r a 60 of CMB + White Noise -T-r-r-r-r-i-pr Smoothed input 40 Springtide PReBeaM -40 -60 100 300 200 400 Figure 5.1: TE spectrum of CMB and white noise for Springtide (blue curve) and PReBeaM (red curve). The smoothed input map is shown in black. Following the example in [5], we reduce £ to £ variation by filtering the spectra by a sliding average (A£ = 20). In this run, the PReBeaM input parameters interpolation-order and zeropadding were set to 1 and 4, respectively. While PReBeaM performs at least as well as Springtide in the TT, EE, and BB spectra, we omit these spectra since the detailed differences are difficult to assess without an in-depth Monte Carlo study. 68 contaminant of the polarized CMB at frequencies below about 80 GHz. The intensity, frequency scaling, and intrinsic polarization fraction of the synchrotron emission depends on the properties of the cosmic rays. The density of the electrons follow a power law of index p, N(E) oc E~p and the synchrotron frequency dependence also follows a power law, S(u) oc vPs (5.1) where (3s = — (p + 3)/2. The spectral index of Galactic synchrotron emission, (3s, varies with frequency and position on the sky. The polarization fraction, /. = ^ + E ( , 2) / s = 3(p + l)/(3p + 7). (5.3) goes as For a typical value p = 3, the spectral index (3s = — 3 and fs = 0.75. A polarization fraction this large is not observed in practice as variations in both the polarization angle along the line-of-sight and beam averaging effects tends to depolarize the signal. Observational constraints from WMAP find an index averaged over regions of high signal-to-noise of (3s = —3.02±0.04 [85]. However, there is still significant uncertainty in the exact variability of the synchrotron spectral index. 5.3.2 Dust emission Thermal emission from warm (10-100K) interstellar dust grains dominates Galactic emission in the 100-6000 GHz frequency range. The thermal emissions originate from thermal fluctuations in the electric dipole moment of the grains. Observations of the polarization of starlight by dust grains suggests that elongated grains are partially aligned with the Galactic magnetic field, although the exact alignment mechanisms are not well understood (see, e.g., [86]). Interstellar dust grains are composed of many different collections of matter particles, which prohibits the construction of a single theoretical emission law for dust. However, an emission law can be fit to observational data. For a frequency range typical of CMB observations, it has been shown that dust emission intensity can be modelled by emission from a dual-component mixture of silicate and carbon grains [87]. The thermal emissions spectrum is modelled as a modified grey-body emission 69 with intensity Iv as Bu(T)ua temperature T. 5.4 where B is the Planck function at frequency v and Testing PReBeaM on polarized foregrounds In this section we describe a test in which we make noiseless foreground maps using the simulated data for the Trieste work [5]. Our motivation for separating out the foreground signal from the CMB signal in this test is to examine how well PReBeaM can recover the polarization amplitude and angle of these foregrounds. We will compare to Madam [54], one of the map-making routines developed by the PLANCK algorithm group for pre-launch data analysis tests. Madam is a destriping algorithm as described in Chapters 3 and 4. It is expected that Madam, which assumes an effective symmetric beam that varies from pixel to pixel, will recover a map and power spectra that distort the polarization parameters of the foreground signals. In this test, PReBeaM and Madam are run on TOD of the synchrotron, dust, and free-free emission (temperature only). There is no CMB signal in this test. The input map of the polarization amplitude of the synchrotron and dust emission are shown in Figures 5.2 and 5.3, respectively. These simulated maps were developed under the PLANCK Sky Model (PSM). The PSM is a simulation tool built by the PLANCK collaboration in preparation for data analysis [88]. Our input TOD also included free-free (bremsstrahlung) emission due to electronelectron scattering in the warm ionized gas of the interstellar medium. The spectrum of free-free emission follows a power-law with a spectral index of -2.1. Scattering directions are random, therefore free-free emission is intrinsically unpolarized. However, secondary polarization can arise at the edges of bright free-free areas such as HII regions (clouds of dust and plasma, undergoing star formation and rich in ionised atomic hydrogen) from Thomson scattering, causing substantial scattering in the Galactic plane. Far from the Galactic plane, at high latitudes, the expected residual polarization is < 1%. Thus, free-free emission is not expected to be a major foreground contaminant for polarization measurements. The fractional difference in the angular power spectrum (defined as {Ctout — Cein)/Cein) of the input and output maps is shown in Figure 5.4. We show the fractional difference spectra for the TT, EE, BB, and cross-correlation TE signals. The beam-mismatch effect described in Section 4.9 is clearly seen in these spectra where the fractional difference is measured to be as high as 30% for Madam at small scales in the EE spectra. On the same plot, the fractional difference in the PReBeaM 70 Synchrotron P (30 GHz) D - 2 . 2 Log (K) Figure 5.2: Simulated map for synchrotron polarization amplitude at 30 GHz. Dust P (30 GHz) ] - 3 . 9 Log (K) Figure 5.3: Simulated map for dust polarization amplitude at 30 GHz. 71 spectra remain closer to zero at all scales and remain under 10% for small scales. Similary, the BB plot shows that the PReBeaM residual remains under 1% for the full range of scales that we examined, while the Madam residual reaches 6% at high £. As a further test, we constructed PReBeaM maps of the free-free emission only. Because free-free emission contains no polarization signal, an examination of the Q and U maps produced by this test will show the effect of leakage from temperature to polarization. We compare to the Madam maps for this case in Figure 5.5. We refer to the Madam results as the binned maps because in the absense of noise, the destriping algorithm is akin to a TOD binning algorithm. We find that the PReBeaM algorithm iterates to a map containing stripe features in the Q component. This translates to some noise in the EE and TE spectra. Excluding a small region at t < 20 in the BB spectra, the PReBeaM polarization spectra have smaller residuals than the binned map spectra by up to three orders of magnitude. This test of the free-free emission only case demonstrates that PReBeaM's beam deconvolution aids in reducing temperature-to-polarization leakage due to the beam mistmatch effect. 5.5 Conclusions In Chapter 4, we examined the results of our PReBeaM algorithm for the case of temperature and polarization signal only for CMB measurements. In this chapter, we separated out noise and foregrounds in order to look at their effects in our maps and spectra. First, we described the noise model for PLANCK and showed that PReBeaM achieved a superior fit to the input TE spectra compared to Springtide, for the case of CMB plus white noise. Next, we briefly described the current state of the knowledge about polarized foregrounds including synchrotron and dust emission. We showed that PReBeaM reconstructed spectra with smaller residuals than Madam for TOD containing synchrotron, dust, and free-free emission. We then separated out the free-free emission, which contains no polarization component, in order to examine the leakage from temperature to polarization. We showed that PReBeaM polarization spectra contain less leakage from temperature than Madam polarization spectra over a wide range of £'s. A class of map-making techniques called destriping methods has been developed to deal with the low-frequency noise component [53]. This approach calculates a series of offsets which are subtracted from the TOD. A binned map can then be made from 72 -a o o _' ' ' i ' ' ' i .., ,i o ' ' ' ' ' ' ' 'o 'LO -! o o -<*• £ o o o i • > *j -i - o o -<* o o K> K> r O O CM O O CN O O O O 03 -d 03 "d , , . i , . . i a , • S-l OO O O O O O O O 0 0 CO ^ t CN CN - ^ o CD 2, o3 CD o o m o o CQ IT) oi o o -<* o o ^ o o o o o o CD PH (-< o to a 03 o o O CM CD CM O o O o CO ^ t CN O CN ^J- UD O O O O O O O o o o o d o o d o o c> ! I I CD .^ fc 03 O 03 ^ 3 . fa H 73 CD CD O -u> CD CD 3 •+J .S-< 1 s. CD IT) o o m O O -<t o o ^- o o K) o o K) O O <N o o <N o o <D P^ c/5 §H >> a g "3 «ttp T«3 o CD .2 'co « 3 CO CD CD O CD S-H -4-3 CD - CO £>> CD CD 3••& 3S3 <-d ^ £ 1 « ? O O O 2 bo S 3 3 "a CD tso O O oo <4=i CO • •—f <*> -—! 4-3 'QT > "2 * o o _ 3 a ^ it r^ I oo I a) i I o o o o tN ^ -* m to l l l l l r^ oo l o o o o o oO •2 g I 25^ oo [Z>W] 2o 3 [^] o o o m . CO o o o o 1 1 • C CQ J ' i' / / 1CD / If T> CD CO C CD o o t-n i" QL Si j*( a: Q- f y ; ] m m TJ 03 F=5 1 O O in o - o ^ o o to o o CM CM o o O O arizat ectra 10 I "o £* a ™ -1 cu r +J 3 •-*a S c6 > ss o S pq * CD 5 a aa. bfl CD -c 3 eBeaM (re of the r . Note t m I J3 '-3 +^ 5 S CD r-Q , CO OH a S3 -u o $ 1.. 00 LM^] h Uri] h _, o CD S3 1 lb 3 o LO lO CD >-< bO E 74 r-1 sA a* +J E <» 3 -Q CD p-J a o a J3 CD - u .22 o s '-S •2 a CO ' C .23 n3 a "3(X CD the resulting TOD. Future tests of the PReBeaM method could include joint work with a destriping algorithm in which the noisy TOD is pre-treated with a destriper before being processed by PReBeaM to remove beam effects. 75 Chapter 6 Variable Point Sources In their most recent release of data, the WMAP team has identified 390 bright point sources in their maps [89]. It is important to study these sources from the CMB standpoint in which extragalactic flat-spectrum radio sources are the most important foreground at small scales. If these sources are not removed properly, they can affect estimations of power spectra and tests of Gaussianity from CMB data [90]. The aim of extraction techniques is to identify and remove these sources as accurately as possible; a greater understanding of these sources and their variability will aid in that goal (see, for example, [91] and [92] for reviews of point source removal methods). Moreover, it is interesting to investigate these sources as astrophysical objects. Sources have been identified to be interesting astrophysical objects such as quasars, Active Galactic Nuclei (AGN), and BL Lac-type objects. The WMAP data allows for a study of the variability of extragalactic sources over a new, largely unstudied, frequency regime from 23-100 GHz. Variability measurements will provide important information on the structure of these sources. In this chapter, we establish a program to search the public WMAP TOD datasets for under-exploited data on the variability and evolutionary properties of point sources cataloged by the WMAP team. The chapter is organized as follows. We review work by the WMAP team on point source variability in Section 6.1. In Section 6.2, the procedure for data extraction and processing is described. Our method for determining excess variability through a x2 analysis is discussed in Section 6.3. Finally, we show our current results and discuss the implications in Section 6.4 and conclude with plans for future work in Section 6.5. 6.1 W M A P findings on point sources In this section, we briefly review the techniques used by the WMAP team to identify point sources in their maps and their findings on point source variability. Further details may be found in the five-year point source catalog paper [89]. 76 Figure 6.1: Map of the 390 point sources detected by the WMAP team. The shaded region shows the mask used to exclude extended foreground emission. The size of the plotted points indicates the flux of the source: the area of the dot scales like the maximum flux over the 5 WMAP bands plus 4 Jy. Image courtesy of the WMAP science team. The WMAP team finds point sources by searching their maps for bright spots that approximate the beam profile. They perform this search in order to generate masks that will be used to discard regions around the brightest sources before proceeding with data analysis. Figure 6.1 is a map showing the location of the 390 point sources that the WMAP team found by searching the individual band maps. The WMAP team cross-correlate detected sources with external surveys. The GB6 [93], PMN [94], and Kiihr et al. [95] catalogs are used to identify 5 GHz counterparts; we will use the same 5 GHz IDs to distinguish point sources. Of the 390 sources, 141 have been optically identified as quasars, 29 as galaxies, 19 AGNs, 19 as BL Lac-type objects, and one as a planetary nebulae [96]. In [89], the WMAP team looks at the variability of sources by forming fluxes from individual year maps. They measure the variability of a source by subtracting the five-year average map from each individual year. They then fit their data points to an arbitrary spectrum that is constant in time and perform a chi-square goodness-of-fit test. They find that 137 of the 390 sources have a x 2 > 37.6 and are therefore variable at greater than 99% confidence. It is found that 54 sources have x 2 > 100 and 15 sources have x 2 > 450. 77 Our analysis of source varability will be more in depth, as we will be looking not at the average year-to-year variation in sources, but rather at all the data from five years worth of observations for each point source. In some cases, this amounts to tens of thousands of data points. We will cross-check our results with the sources that the WMAP team found to be variable. 6.2 Description of data extraction and processing Fortunately, WMAP has made their data available for download on a public web page1. All of the data that we use, as well as the beam profiles, source catalogs, and masks, can be found on their site. We also made extensive use of the WMAP Interactive Data Language (IDL) Library 2 . The library contains a suite of procedures and routines to assist with reading and manipulating the data products. The WMAP instrument is composed of 10 differencing assemblies (DAs), covering 5 frequency bands from 23-94 GHz. The DAs are identified as follows: 1 DA at 23 GHz (Kl), 1 DA at 33 GHz (Kal), 2 DAs at 41 GHz (Q1,Q2), 2 DAs at 61 GHz (V1,V2), and 4 DAs at 94 GHz (W1-W4). Each DA is formed from two differential radiometers. We use the calibrated TOD files in which gain and baseline factors have been removed. Each TOD fits file contains 1875 records making up 24 hours of measuring time. One record has 30 major science frames of length 1.536 seconds. One major science frame consists from 12 measurements (for the K band) up to 30 measurements (for the W band). The calibrated archives have been coadded into two radiometer channels for each measurement so we read two vectors of data from each TOD fits file. The general flags and differencing assembly flags for each record are also read. Flags may indicate quality issues within the science frame, such as the sun, moon or a planet is visible in one of the beams. All records with any nonzero flags are rejected. 6.2.1 Loss imbalance parameter WMAP is a differential instrument consisting of two back-to-back Gregorian telescopes. The instrument records a TOD measurement as the difference between the 1 2 http://lambda.gsfc.nasa.gov/product/map/dr3/m_products.cfm http://lambda. gsfc.nasa.gov/product/map/current/m_sw.cfm 78 A-side and B-side sky signals. We will be interested in separating out the differential measurement into that signal which is seen by the A-side and that seen by the B-side. However, we will need to account for the fact that the output of each WMAP radiometer does not produce a perfect differential response. Power transmission coefficients, a A and aB, of the optics and waveguide components couple the sky signals from the A and B side beams to the radiometer. In reality, these coefficients are not equal. The output of a radiometer, S, when the beams observe regions of the sky with temperatures TA and TB becomes S = G{aATA-aBTB\ (6.1) where G is the gain of the radiometer. Jarosik et al [97] define the loss imbalance parameter as xwhere (3C = G(aA — 6.2.2 OLB)/2 (6.2) =&- and (3d = G(aA + «B)/2. Differencing assembly parameters For each DA, we define the following quantities: • 0b = beamsize (we choose beamsize to be the point at which the normalized beam profile falls to 75%) • Xim,i,2 — l° s s imbalance parameters in radiometer 1 and 2 • UJA = 1 + Xim • U>B = 1 - Xim • (To = noise per observation to ~ 0.1% uncertainty (in mK) The values for each of these parameters for each DA are given in Table 6.2.2. 6.2.3 Beam normalization In the five-year analysis, an azimuthally symmetrized radial beam profile is calculated by binning individual A- and B-side hybridized Jupiter observations. As such, the beam profiles are presented in units of mK with a constant bin size of 0.25'. For 79 DA Kl Kal Ql Q2 VI V2 Wl W2 W3 W4 8b (deg) •Eim 0.2521 0.0028 0.1979 0.00245 0.1521 0.0057 0.0089 0.1563 0.1063 0.00685 0.0094 0.1063 0.0688 0.00575 0.0646 0.00935 0.0646 0.0059 0.0688 0.02135 U1A UJB 1.0028 1.00245 1.0057 1.0089 1.00685 1.0094 1.00575 1.00935 1.0059 1.02135 0.9972 0.99755 0.9943 0.9911 0.99315 0.9906 0.99425 0.99065 0.9941 0.97865 Table 6.1: Beamsize (8b), loss imbalance parameter (xim), for each of the differencing assemblies (DA). a0 (mK) 1.436 1.470 2.254 2.141 3.314 2.953 5.899 6.565 6.926 6.761 U>A,B, and noise term (<70) the purposes of our analysis, we require a normalized beam profile. We calculate the normalization factor, Tpeak, following the description in Section 4.1 of the five-year beam paper [98], •'-peak = 1 Jupiter TZ. ^ ^beam l^-"/ The relevant values for Tjupiter, flfid, and f2&eam can be found in Tables 4 and 5 of [98]. We plot the inner parts of the normalized profiles in Figure 6.2. We show here only the region of the beams that we use for our analysis. Data points which fall at a distance from the point source where the beam profile has fallen to less than 75% are omitted. The peak is noisy because the inner region of the beam is directly from Jupiter data in the five year analysis. Some of the differencing assemblies have a peak in the innermost bin and some a dip. 6.2.4 Dipole corrections For unbalanced differential data, the TOD at time t is given by d(t) = UJAKPA) - UBKPB) + 0(5xim) (6.4) where \(PA) is the Stokes parameter I at pixel pA and \(PB) is the Stokes parameter I at pixel PBWe accept only TOD for which the following two conditions hold: (i) the distance between the pointing for the TOD element and the point source is less than 8b; (ii) the opposite horn is outside of a region that masks a point sources with a disc of radius 0.6°. 80 If the above conditions hold, then we calculate the variation for a point source observed in a single horn as _ dit) - {uAIA{pA) - LUBIB(PB)) x— - - dTf - dTv (o.oj Jbeam where IA is the Stokes parameter I from the uncleaned WMAP map at pixel pA, dTf is the fixed dipole component, dTv is the velocity dipole term, and fbeam is the value of the normalized beam profile at the pointing distance from the point source. The fixed dipole modulation term, dTf, due to the solar system barycenter's (SSB) motion with respect to the CMB is dTf = uAdA — u>BdB- (6.6) Here, dAtB are the values of the fixed dipole component at the A- and B-side, found by dotting the fixed dipole vector (given as d = [-0.233, -2.222,2.504] mK in [24]) with the position vectors for the A- and B-side. The dipole modulation term, dTv, due to WMAP's velocity with respect to the solar system barycenter (SSB) is calculated as T0 dTv = —v • (uAnA - w B n B ) (6.7) c where T0 is the CMB mean temperature, c is the speed of light, v is the WMAP velocity with respect to the SSB, and UA, fiB are the line-of-sight vectors for the Aand B-sides, respectively [24]. The error on x is calculated as (6.8) <yx = -pJ beam where 0.75 < fbeam < 1-0. 6.3 Chi-square analysis Given xA, XB, PX,A, and ax<B from the previous section, we calculate the chi-square minimum and chi-square probability following [99]. For n0bs repeated measurements (XJ, <7j) of the same quantitiy x, the x2 IS defined as: n obs / 81 ~\2 X2 is minimized by E(iM2) (6.10) Therefore, the chi-square minimum is given by Xmin = x V ) - (6.11) We compute xmin a n d compare with x2(n) t o determine the goodness of fit where n is the number of independent degrees of freedom. We find p, the probability that, in a x2 distribution with n degrees of freedom, a random variable is greater than Xmin'p = Prob{ X 2 (n) > xLn> (6-12) where, in our case, n = n0bs — 1, with one degree of freedom used to find x. 6.3.1 Null test Next, we do a null test as a means of determining a detection of excess variability in the actual point sources. For the null test we choose 1000 points outside of the WMAP mask, which are expected to contain CMB-only signal, and perform same analysis described above on the null points as we do for the actual point sources. This gives a set of x2 probabilities for the null points. We then histogram the probabilities and identify those actual point sources having a p that fall below some cutoff value on the null histogram. In this way, we may still be confident in our detection of excess variability even in the event that we have under or overestimated the noise term. 6.4 Results and discussion The x2 analysis described above was performed on the 390 actual point sources and 1000 null points using the K (23 GHz) and Ka (33 GHz) band data as a first test. The distribution of the x2 probabilities for the null points is plotted as a histogram, seen in Figures 6.3, 6.4, 6.5, and 6.6. The value of the x2 probability for those actual point sources whose value falls below the 5% cutoff on the null histogram is then overplotted as a vertical line. The results are summarized in Table 6.2, in which points are listed as being observed to be variable in either the A- side or B-side or both. Point sources are identified, as in the WMAP analysis, according to their 5 GHz ID. We also list the K and Ka band flux in Janskies (Jy), the type of object 82 as identified by [96] (if available), and a note which indicates the level of variability found in the WMAP analysis. K 5 GHz ID J0003-4752 PMN J0006-0623 PMN J0025-2602 GB6 J0029+0554B A B PMN J0204-1701 GB6 J0221+3556 PMN J0403-3605 Uy 0403+76 PMN J0423-0120 GB6 J0424+0226 PMN J0442-0017 PMN J0457-2324 PMN J0538-4405 PMN J0540-5418 GB6 J0542+4951 GB6 J0555+3948 PMN J0609-1542 PMN J0633-2223 PMN J0634-2335 GB6 J0646+4451 K Flux (Jy) / / 2.3 0.9 1.1 / 0.7 1.0 1.9 1.4 1.1 1.8 0.8 3.8 1.2 1.3 0.7 1.2 3.4 1.0 8.2 1.2 0.9 2.4 5.6 1.4 1.7 3.0 3.7 0.5 0.6 2.9 / / PMN J0038-2459 PMN J0050-0928 GB6 J0108+0135 GB6 J0108+1319 PMN J0125-0005 PMN J0132-1654 PMN J0133+5159 GB6 J0136+4751 GB6 J0152+2206 GB6 J0204+1514 Ka A B / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / Continued on next page 83 Ka Flux (Jy) 0.7 2.3 0.7 1.3 0.8 1.0 1.9 1.1 1.2 1.8 1.1 3.8 1.3 1.3 1.2 3.8 0.7 8.3 1.0 0.8 2.5 5.9 1.4 1.3 2.4 3.3 0.6 0.7 2.4 Type Note G V QSO V QSO V V QSO V QSO V QSO QSO V V V QSO V V V Table 6.2 - continued from previous page K Ka K Ka 5 GHz I D GB6 J0720+0404 GB6 J0753+5353 GB6 J0825+0309 GB6 J0841+7053 PMN J0847-0703 GB6 J0909+0121 GB6 J0914+0245 GB6 J0920+4441 GB6 J0927+3902 GB6 J0957+5522 PMN J1014-4508 PMN J1018-3123 GB6 J1022+4004 GB6 J1043+2408 J1053+8109 GB6 J1058+0133 PMN J1118-1232 PMN J l 147-3812 GB6 J1159+2914 PMN J1215-1731 GB6 J1219+0549 GB6 J1229+0202 GB6 J1230+1223 PMN J1246-2547 PMN J1256-0547 GB6 J1332+0200 GB6 J1333+2725 PMN J1337-1257 PMN J1427-3306 GB6 J1436+6336 GB6 J1443+5201 GB6 J1549+5038 A B A B Flux (Jy) Flux (Jy) 0.9 1.0 1.6 1.7 0.9 2.1 1.4 1.3 6.7 0.9 1.1 0.9 0.9 0.8 0.9 4.6 1.0 2.1 0.7 1.0 1.9 1.7 1.0 2.0 1.6 1.4 / / / / / / / / / / / / / / / / / / / 2.0 1.4 2.7 20.0 19.7 1.3 17.1 1.4 0.8 5.8 1.0 0.5 0.8 0.9 / / / / / / / / / / / / / / / / / / / / / / / / Continued on next page 84 5.8 0.9 0.8 0.9 0.9 0.8 0.9 4.4 0.9 2.3 2.2 1.2 2.2 18.4 15.5 1.4 17.9 1.4 0.9 6.0 1.4 0.9 0.8 Type Note V QSO V QSO V V QSO V QSO V V QSO V QSO G QSO QSO V V QSO V V V Table 6.2 - continued from previous page K 5 GHz I D A GB6 J1608+1029 J1633+8226 GB6 J1635+3808 GB6 J1642+3948 GB6 J1657+5705 GB6 J1658+0741 / GB6 J1727+4530 Uy 1803+78 PMN J1819-6345 PMN J1834-5856 GB6 J1848+3219 PMN J1923-2104 PMN J1924-2914 PMN J1939-6342 GB6 J1952+0230 PMN J2056-4714 GB6 J2123+0535 PMN J2131-1207 PMN J2134-0153 GB6 J2148+0657 GB6 J2202+4216 GB6 J2212+2355 PMN J2218-0335 PMN J2225-0457 GB6 J2253+1608 PMN J2258-2758 GB6 J2322+5057 GB6 J2330+1100 GB6 J2349+3849 PMN J2358-6054 / B Ka A B K Flux (Jy) Ka Flux (Jy) 2.0 2.0 1.5 4.3 6.0 0.6 1.5 1.0 1.7 1.5 1.1 1.3 3.9 6.5 0.5 1.4 0.9 1.8 1.7 1.1 0.7 2.3 12.3 0.9 0.8 2.2 2.2 2.7 2.0 8.0 3.4 1.3 2.3 5.2 7.4 5.2 0.9 1.0 0.8 1.9 / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / / Continued on next page 85 0.8 2.5 12.0 0.7 0.6 2.5 1.8 2.4 1.9 7.7 3.5 1.5 2.0 4.9 7.5 5.2 0.8 1.0 0.7 1.4 Type Nol V QSO QSO V V V V V V QSO QSO V V V QSO QSO V V V QSO QSO QSO QSO V V V V Table 6.2 - continued from previous p a g e K Ka^ K Ka^ 5 GHz ID A B A B Flux (Jy) Flux (Jy) Type Note Table 6.2: Point sources that we have identified as being variable. Sources are listed, following the WMAP convention, by their 5 GHz ID. Variability observed in either the A or B side (or both) is indicated by a checkmark for the K and Ka band separately. The peak K and Ka band flux is listed in columns 6 and 7, respectively. The type (G for galaxy, and QSO for quasar) is listed in column 8. v WMAP-identified as probable variable: x2 > 36.7. v WMAP-identified as variable: x2 > 100. Several observations can be drawn from Figures 6.3, 6.4, 6.5, 6.6 and Table 6.2. First, it was expected that the null points should have a uniform distribution of probabilities between 0 and 1. Rather, we find that they are peaked near p = 1 indicating that the null points are, on average, less variable than we expected. This suggests the possibility that the noise term has been overestimated. We find that reducing the noise term by 1.5% does indeed produce a uniform x2 probability distribution for the null points. This finding is in direct conflict with the results from a recent Bayesian analysis of the white noise levels in the WMAP data by Groeneboom et al. [100], who find an indication that the noise levels in the maps are underestimated by 0.5-1.0%. A second feature that one observes in Table 6.2 is the overall discrepancy between the x2 probabilities measured in the A and B side for the same point source. Of the 63 sources that we identify as being variable in either the A or B side for the K band, only eight are observed to be variable in both horns. Similarly, of the 43 sources identified as variable in the Ka band data, only seven are found to be variable in both horns. These sources that are seen to be variable in both horns also tend to be the brightest sources as indicated by the fluxes in Table 6.2. As we have only examined the data from the K and Ka band, it is possible that we have not reached the threshold of detectability for many of these sources. The test that we have performed to identify variability in point sources in each of the bands and horns individually is much more stringent than the WMAP test which combines all of the data. 86 6.5 Conclusions and further work Using the raw, calibrated TOD from the K and Ka band DAs of the WMAP instrument, we identify 10 point sources that display excess variability in both the A- and B-side, and in either the K band, Ka band, or both. These sources, and their optical identifications, are: PMN J0403-3605 (quasar), PMN J0423-0120 (quasar), PMN J0538-4405 (quasar), PMN J1014-4508, GB6 J1229+0202 (quasar), GB6 J1230+1223 (radio galaxy) PMN J1246-2547 (quasar), PMN J1256-0547 (quasar), PMN J13371257 (quasar), and PMN J2258-2758 (quasar). All of these sources (except PMN J1014-4508) have also been flagged as variable or probable-variable by the WMAP team. Subsequent work on the variability of point sources in the WMAP data will include an analysis of the data from the remaining differencing assemblies (Ql, Ql, VI, V2, Wl, W2, W3, and W4). We will be interested to see if the differing noise properities in these channels will produce results which are distinct from those that we found here for the K and Ka band analysis. 87 profile Ko1 E o j 1.0 "^-^V^_. " ^^v. •a 0.9 | o z 0.00 0.05 0.10 0.15 0.20 0.25 0.30 Distance from beam center (deg) Q1 1.00 0.95 0.8 i ^v i 0.7 0.00 0.05 0. 10 0.15 0.20 Distance from beam center (deg) Q2 1.00 y s ^ A •A^^ : ^^^v. \_ 0.90 0.80 \~, 0.90 V. 0.85 - 0.95 ^^\ \ 0.85 0.75 n 70 ^\ 0.80 X . 0.75 0.70 v. ' 0.00 0.05 0.10 0.15 Distance from beam center (deg) : ^\ 0.00 0.05 0.10 0.15 Distance from beam center (deg) V2 1.00 " 0.95 V 0.90 x^ 0.85 \ N. 0.80 \^ : 0.75 0 70 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Distance from beam center (deg) W2 0.00 0.02 0.04 0.06 0.08 0.10 0.12 Distance from beam center (deg) 1.00 0.95 0.90 0.85 0.80 0.75 0 70 1.00 "^~~\ 0.95 ^ v 0.90 \v : - ^X\ . 0.85 \/^~\ \ . >y \ 0.80 \ 0.75 0 70 \ \ 0.00 0.02 0.04 0.06 0.08 Distonce from beam center (deg) \ ' 0.04 0.06 0.08 0.00 0.02 Distance from beam center (deg) W4 -. 1.00 0.95 0.90 1 N^ : 0.85 0.80 0.75 0 70 0.00 0.02 0.04 0.06 0.08 Distance from beam center (deg) \ . \^ - 0.00 0.02 0.04 0.06 0.08 Distance from beam center (deg) Figure 6.2: Innermost parts of the normalized profiles for all 10 channels of the WMAP instrument. 88 . I i i i I 111 I i i i I i • H J1256-0547Q . . . i . . . i . . 0 0.2 0.4 0.8 0.8 I |."' j'64d3-3665 p F 0 .1 0.2 i i i 1 0.4 0.6 0.8 1 L J2131-12P7 . . . I • i . I • . • I • P t ''iJyu403'+7'6 50 tQ ct 0 • . i . • 0.2 1 i . • • i • • • i . 0.4 0.6 . . n 0.6 1 P rj,ie<l84iib29 i J2225-|0457[: I i r... I . . . i . . . I i f 11 I 1 1 1 1 1 1 i I 11 i'l 1 1 L " ' J ' ' ' I "T T T T T T T T J0424+02|26 100 j J1Q43+24C ;i i SOP I | I 1 I I I I I I 1 I I HI I I II " | n i | i i i | i n | k ; J2349+q849(j ,oo|| J1246-2547 Figure 6.3: The 24 most variable point sources found on the A side of the K band. The distribution of the x 2 probabilities for the null points plotted as a histogram for the A side (blue) and B side (green). The x2 probability, p, for the actual point sources are overplotted for the A side (dotted) and B side (dashed). This figure shows the point sources ordered in increasing p for side A. Note that the dotted and dashed lines for sources J0423-0120, J1229+0202, J1256-0547, J2258-2758, and J0538-4405 are difficult to see as they sit nearly exactly at zero. 89 rJ1266-6547 rj'22l2+2355Q 0 * * ' * * * * * * * * * '•' * * * * 0 0.2 0.4 0.6 0.B 1 P L"j'6669:-i542B t jb72b+64&4' 0 1 . . . 1.... 1 . . . 0 0 " ' " 0 0.2 0.4 0.6 1... 0.6 Ml I I I I I I I I I I I I I I I, I I r w r p » 1 1 1 1 1 j 111111 & J2218-0335 h J0646+4451 1 • ' • • • ' • • • ' • • • ' • • 0.2 0.4 0.6 0.8 1 P • • • 1 1 1 1 1 • 1 1 1,1 1 1 1 1 1 tj 71250+1223' 100 : jj J0633-2223 ; so 0 0 0.2 0.4 0.6 0.6 I P III Jtl 100 •4 J 0 4 5 7 - 2 3 2 4 i I I I I I I I I I I T I Fj J1147-381213 J'i635+'386'8 50 1 . , i . . i i . . • i •'. i i n i 0 0.2 0.6 04 0.8 1 P Figure 6.4: The 24 most variable point sources found on the B side of the K band. The distribution of the x2 probabilities for the null points plotted as a histogram for the A side (blue) and B side (green). The x2 probability, p, for the actual point sources are overplotted for the A side (dotted) and B side (dashed). This figure shows the point sources ordered in increasing p for side B. Note that the dotted and dashed lines for sources J0423-0120, J1229+0202, J1256-0547, J2258-2758, and J0538-4405 are difficult to see as they sit nearly exactly at zero. 90 0 i I r i i I i r » i i i i 11 i i I I i'i • i i I I I l i • i i l i i I i j J1256-0547, t J1229+0202J 0.2 0.4 0.6 0.8 1 0 li i i I i i i I i i i I i i i I i 0 0.2 0.4 0 . 6 0 . 8 t" j'6423-6'l2b j » * » • • ' • ' • 0.2 0 1 0.4 0.6 0.6 ! P P 11111 • 11 • 1 1 1 1 1 1 1 1 1 1 1111 0 PI 1 1 1 1 1 111 111 • i L"j'6i'62+?26'6 : Jl|939-6342j 100 r . . • i J• • i • . . i • . . 0.2 0 0.4 0.0 i...' OB 1 0 p . J0841+70Q3 St 0.2 0.4 0.6 0.8 1 P Mil I E j • • i • • • i • • • t • . . n. • • 0 0.2 0.4 O.S 0.8 P r . . . i .i.. i . . . i . . . i . . . i | I M L J195^+0230j 100 E 60 SO 0 .IM|llljim • • • I •• • ! • • • I • • i r i i ."j'C|i^6+475iJ I 100 . . • I . I . I . . . I I . . I . M 1 0.2 0.4 0.6 O.C P . . . . I . . • I .L . I . . • I . . . 0.2 0 0.4 0.6 0.8 I P J I I | I I I P I I I | I I I | I I I 1,1 i i I 111 I i i i I i i i i 11 i i J18l|9-6345| sot r . 0 • . I • . . 0.2 I • • • i • • L i . 0.4 0.8 • • O.E P II I I | I H J I II | I I l | l l l R J09,20+4441 F F; 0 I 0.2 •••••••• 0.4 0.6 0.8 P TTTTTTTTTTTTTTTTT J0006-0623 F ^ii[iiiiiiT|M r i m loot-! J1058+0133^ 60 -; Q I*-. . I . . . I . . . I . . . I . . I 0 0.2 0.4 0.6 0.8 1 P Figure 6.5: The 24 most variable point sources found on the A side of the Ka band. The distribution of the x2 probabilities for the null points plotted as a histogram for the A side (blue) and B side (green). The \2 probability, p, for the actual, point sources are overplotted for the A side (dotted) and B side (dashed). This figure shows the point sources ordered in increasing p for side A. Note that the dotted and dashed lines for sources J1258-0547, J1229+0202, J2258-2758, J0423-0120, J0538-4405, and J1014-4508 are difficult to see as they sit nearly exactly at zero. 91 TTTTTTTTTTTTTTTTTm 111 ' l i p I M 1 1 1 I I J1256-0547, 100 so(0 0 • • • 111 i i i i i . J0423-0120, 100 , i n |IIi n n|I IMIIIJ 1 1 1 i i i i i j : J0538-4405j 100 . J1229+0202j 100 50 0.2 0.4 0.8 0.8 1 P 0 o Q I . . . * 0 0.2 0.4 o.e o.e p 0.2 0.4 0.8 0.8 P 0 • | I I I I I I I I I I I I I I I J0403-3605, 100 50 0 i • i I i t | I • • . I 0 f • f I • • . 0.2 0.4 0.8 0.8 P P I 1,1 I I I M i l 1 I I I I I I I I L| J0108+0135 6 I i o-rf^ 0 0.2 0.4 0.8 0.8 P 0.2 0 4 0 6 0.8 P 1 I PI I I I ITT'M'I'Tn I 1.1 I I I I L| J1219+0549d J . ..I . . . I . . . I . . . I . . • 0 I . 32253+1608* 100 60h I... 0.2 0.4 0.8 0.8 P 0 iluli I i l i i i l i i i l l n i l 0.2 0.4 0.8 0.8 I P Figure 6.6: The 24 most variable point sources found on the B side of the Ka band. The distribution of the x2 probabilities for the null points plotted as a histogram for the A side (blue) and B side (green). The x2 probability, p, for the actual point sources are overplotted for the A side (dotted) and B side (dashed). This figure shows the point sources ordered in increasing p for side B. Note that the dotted and dashed lines for sources J1258-0547, J0423-0120, .10538-4405, J1229+0202, J2258-2758, and J1014-4508 are difficult to see as they sit nearly exactly at zero. 92 Chapter 7 Conclusions 7.1 Summary of Progress In this chapter, I summarize the achievements that were completed during the thesis work described in this dissertation. I developed deconvolution map-making for temperature observations; this included tests of several beam models and scanning strategies, and included Galactic foregrounds. This work was done at a relatively low resolution with beams of 3° or larger. I showed that deconvolution map-making was successful at removing systematic effects in CMB maps due to beam asymmetries. I then extended the deconvolution map-making technique to include polarization observations (PReBeaM) and tested the method on actual PLANCK simulated data. This study was done at a resolution of about 30', or about 0.5°. Realistic PLANCK beams and scan strategy were used in the tests and we showed that PReBeaM corrected shifts in the CMB power spectra due to beam asymmetry effects. Final tests of PReBeaM included white noise and polarized foreground signals. Finally, I developed a program to search public WMAP data archives for point source variability on a wide range of time scales. 7.2 Future work with PReBeaM There is still much that can be done to improve the PReBeaM algorithm and much that can be done to research the effect of beam systematics in CMB maps with the aid of our technique. Here I list a few areas of possible future work. 7.2.1 Improving the computational requirements The work described in this dissertation has resulted in the development of a mapmaking code which has been shown to be able to remove a significant level of sys- 93 tematic beam effects from maps and power spectra for the PLANCK satellite. Up to this point, the main focus of my PhD research has been to develop this code and to demonstrate its effectiveness at beam deconvolution, rather than to optimize the computational requirements of the algorithm. While efforts have been made to streamline the code, there is certainly room for improving both the memory requirements and the computational runtime. One such improvement that we have identified is to turn the entire operation of interpolation and transpose interpolation of all the TOD elements (described in Chapter 4) into a matrix operation. This has the dual advantage of omitting the need to store a TOD-sized object (less memory required) and removing the steps in which we run through each TOD element for every iteration (faster computational runtime). 7.2.2 Interface with destriper for noise treatment In Chapter 5, we presented the results from a PReBeaM run on PLANCK simulated data containing CMB and white noise. In reality, PLANCK and other satellite instrument data will contain correlated types of noise. A class of map-making techniques called destriping methods has been developed to deal with the low-frequency noise component [53]. This approach calculates a series of offsets which are subtracted from the TOD. A binned map can then be made from the resulting TOD. There is great potential for PReBeaM to be interfaced with a destriping algorithm in which the noisy TOD is pre-treated by the destriper before being processed by PReBeaM to remove beam effects. 7.2.3 Re-analysis of the W M A P data Through our involvement in the PLANCK collaboration, we have developed PReBeaM to run on simulated PLANCK data. However, the algorithm is easily adaptable to run on data from other CMB satellite instruments. A recent paper by Wehus et al. [101] has assessed the impact of asymmetric beams on cosmological parameters estimated from the WMAP temperature data. They find shifts of 0.4a in the amplitude of scalar perturbations and the physical density of cold dark matter, and a shift of 0.3<r in the spectral index of scalar perturbations. Wehus et al. point out that shifts this size should not be neglected as they are of the same order of magnitude or larger than marginalization over the Sunyaev-Zeldovich effect or unresolved point sources. It would be of interest to re-analyze the WMAP 94 data, which is freely available, using PReBeaM and the actual measured asymmetric WMAP beams in order to examine the level at which beam deconvolution can correct these observed shifts. It would be of great interest to repeat the full analysis on the polarization data which also suffers from beam effects and confusion due to leakage from the much larger temperature signal. 7.2.4 Analysis of Planck data with far sidelobes For the PLANCK instrument, a complete characterization of the beam includes both the main beam and the far sidelobes. Sidelobes are located as far away as 90° from the main focal plane and may induce systematics on the largest scales, potentially inhibiting the detection of B-modes. A major strength of the deconvolution technique is the ability to account for significant beam asymmetries such as sidelobes with relative ease using the harmonic method. In Chapter 3, we tested deconvolution mapmaking for temperature observations at low resolution with a simulated far-side lobe beam. We would like to repeat the PLANCK 30 GHz temperature and polarization analysis using the full beam model, that is main beam plus sidelobe. This gives rise to the unresolved question of how to best combine the benefits of a local, pixel-based approach for the innermost part of the main beam, with the advantages of harmonic space convolution for the remaining beam pattern. 7.3 Future work with the convolution routines At its most fundamental level, the standard model of cosmology assumes a homogeneous and isotropic universe. However, there have been numerous studies of the WMAP data resulting in claims of evidence for a preferred direction in the universe (e.g., most recently [102, 103, 104]). A preferred direction is defined as a direction in which local features in the CMB tend to be oriented towards. Recalling that the temperature field on the sphere can be expanded in harmonic coefficients aem and assuming an isotropic Gaussian process, each a^m is an independent Gaussian random variable and is fully characterized by the power spectrum CeOne finds no evident shape to the multipole if Ce is randomly and uniformly distributed among all ms. Therefore, the identification of an orientation for which an excess amount of power is concentrated in a particular m-mode suggests deviation from isotropy. Some of the numerical techniques that we develop for deconvolution map-making 95 can be naturally applied to this problem, giving an immediately science return when applied to the WMAP data. The convolution algorithm that we use was developed to compute the detector response at every pointing of the telescope, in other words, to perform the convolution of the beam with the sky for all orientations of the beam. We can use our fast convolution algorithm to search for correlations in the CMB anisotropies. These techniques will allow us to quickly and efficiently do this search out to high £ and for all orientations. In addition to searching for non-Gaussian signatures such as filaments and preferred directions, we can use our tool to expose systematic effects left in the data. For example, scan-synchronous effects may leave correlated noise stripes in the data products. Given full sky maps of Galactic dust emission and other astrophysical emissions, we can use the same formalism to search for anomalous structure in these signals. 96 References [1] J. Peacock, Cosmological Physics (Cambridge University Press, Cambridge, 1999). [2] J. Bock et al, Task Force on Cosmic Microwave Background Research, 2006, arXiv:astro-ph/0604101vl. [3] G. Efstathiou, C. Lawrence, and J. T. (eds), ESA-SCI(2005)-1, ESA Publications (2005). [4] N. Jarosik et al, The Astrophysical Journal 170, 263 (2007). [5] M. A. J. Ashdown et al., Astron.Astrophys. 493, 753 (2009). [6] A. Penzias and R. Wilson, ApJ 142, 419 (1965). [7] R. Dicke, P. Peebles, P. Roll, and D. Wilkinson, ApJ 142, 414 (1965). [8] G. Smoot et al, ApJL 396, LI (1992). [9] E. Komatsu et al, ApJS 180, 330 (2009). [10] A. H. Guth, Physical Review D 23, 347 (1981). [11] A. Linde, New Astronomy Review 49, 35 (2005). [12] D. Lyth and A. Riotto, Phys. Rept. 314, 1 (1999). [13] J. Khoury, B. A. Ovrut, P. J. Steinhardt, and N. Turok, Physical Review D 64, 123522 (2001). [14] P. J. Steinhardt and N. Turok, Science 296, 1436 (2002). [15] D. Baumann et al, CMBPol Mission Concept Study: A Mission to Map our Origins, 2008. [16] P. Peebles, Principles of Physical Cosmology (Princeton University Press, Princeton, New Jersey, 1993). [17] E. Kolb and M. Turner, The Early Universe (Addison-Wesley, USA, 1994). [18] W. Hu and S. Dodelson, ARA&A 40, 171 (2002). 97 [19] A. Friedman, Z.Phys 10, 377 (1922). [20] A. Friedman, Z.Phys 21, 367 (1924). [21] H. Robertson, ApJ 82, 284 (1935). [22] H. Robertson, ApJ 83, 187 (1936). [23] A. Walker, Proc.Lond.Math.Soc (2) 42, 90 (1936). [24] G. Hinshaw et ai, ApJS 180, 225 (2009). [25] A. Linde, Phys. Lett. B 108, 389 (1982). [26] A. Albrecht and P. Steinhardt, Phys. Rev. Lett. 48, 1220 (1982). [27] M. S. Turner and M. White, Phys. Rev. D 53, 6822 (1996). [28] J. Khoury et ai, Physical Review D 65, 086007 (2002). [29] L. A. Boyle, P. J. Steinhardt, and N. Turok, Physical Review D 69, 127302 (2004). [30] A. Starobinsky, JETP Lett. 40, 409 (1979). [31] R. Sachs and A. Wolfe, Astrophys.J. 147, 73 (1967). [32] M. White and W. Hu, Astron.Astrophys. 321, 8 (1997). [33] W. Hu and N. Sugiyama, Physical Review D 51, 2599 (1995). [34] J. Silk, Astrophys. J. 151, 459 (1968). [35] R. Crittenden, Polarization-Temperature Correlations in the Microwave Background, 1994. [36] A. Kosowsky, Annals Phys. 246, 49 (1996). [37] U. Seljak, Astrophys. J. 482, 6 (1997). [38] M. Zaldarriaga and U. Seljak, Physical Review D 55, 1830 (1997). [39] M. Kamionkowski, A. Kosowsky, and A. Stebbins, Physical Review D 55, 7368 (1997). [40] J. Kovac et ai, Nature 420, 772 (2002). [41] E. M. Leitch et al, Nature 420, 763 (2002). [42] A. C. S. Readhead et ai, Science v 306, 836 (2004). [43] D. Barkats et al, The Astrophysical Journal 619, L127 (2005). 98 [44] G. Polenta et al, Advances in Space Research 36, 1064 (2005). [45] L. Page et al, The Astrophysical Journal 170, 335 (2007). [46] Y.-T. Lin and B. D. Wandelt, Astroparticle Physics 25, 151 (2006). [47] B. Wandelt and K. Gorski, Phys.Rev.D 63, 123002 (2001). [48] C. Armitage and B. D. Wandelt, Physical Review D 70, 123007 (2004). [49] J. Arnau, A. Aliaga, and D. Saez, A&A 382, 1138 (2002). [50] C. Burigana et al, Astron.Astrophys. 373, 345 (2001). [51] C. Barnes et al, The Astrophysical Journal 148, 51 (2003). [52] T. Poutanen et al, Astron.Astrophys. 449, 1311 (2006). [53] M. A. J. Ashdown et al, Astron.Astrophys. 467, 761 (2007). [54] E. Keihanen, H. Kurki-Suonio, and T. Poutanen, Mon.Not.Roy.Astron.Soc. 360, 390 (2005). [55] O. Dore et al, A&A 374, 358 (2001). [56] G. De Gasperis et al, A&A 436, 1159 (2005). [57] D. Yvon and F. Mayet, Astron.Astrophys. 436, 729 (2004). [58] P. Natoli, G. De Gasperis, C. Gheller, and N. Vittorio, A&A 372, 346 (2001). [59] M. Tegmark, Astrophys.J. 480, L87 (1997). [60] R. Stompor et al, in From the time-ordered data to the Maximum-Likelihood temperature maps of the Cosmic Microwave Backgorund anisotropy, Proceedings of the MPA/ESO/MPE Workshop, Mining the Sky (PUBLISHER, Garching, 2000), astro-ph/0012418. [61] X. Dupac and M. Giard, MNRAS 330, 497 (2002). [62] C. H. Lineweaver et al, Astrophys.J. 436, 452 (1994). [63] E. L. Wright, Scanning and Mapping Strategies for CMB Experiments, 1996, astro-ph/9612006. [64] G. Hinshaw et al, The Astrophysical Journal 148, 63 (2003). [65] J. V. Arnau and D. Saez, New Astronomy 5, 121 (2000). [66] C. Burigana and D. Saez, A&A 409, 423 (2003). [67] A. D. Challinor et al, Mon.Not.Roy.Astron.Soc. 331, 994 (2002). 99 [68] D. M. Brink and G. R. Satchler, Angular Momentum (Clarendon Press, Oxford, 1962). [69] T. Risbo, Journal of Geodesy 70, 383 (1996). [70] J. Reid, Large Sparse Sets of Linear Equations (Academic Press, London, 1971), pp. 231-254. [71] J. R. Shewchuk, 1994, conjugate-gradient. ps. http://www.cs.cmu.edu/~quake-papers/painless- [72] W. Press et al., Numerical Recipes: the Art of Scientific Computing (Cambridge University Press, Cambridge, 1992). [73] K. M. Gorski et al, The Astrophysical Journal 622, 759 (2005). [74] C. L. Bennett et al, Astrophys.J. 583, 1 (2003). [75] C. Armitage-Caplan and B. D. Wandelt, ApJS 181, 533 (2009). [76] T. P. Collaboration, The Scientific Programme of Planck, 2005, astro-ph 0604069vl. [77] T. Souradeep et al., New Astronomy Review 50, 1030 (2006). [78] S. Mitra et al, Mon.Not.Roy.Astron.Soc. 394, 1419 (2009). [79] J. Borrill, MADCAP - The Microwave Anisotropy Dataset Computational Analysis Package, 1999, arXiv.org:astro-ph/9911389. [80] A. Challinor et al, Physical Review D 62, 123002 (2000). [81] D. Harrison, F. van Leeuwen, and M. Ashdown, 2009, in preparation. [82] M. A. J. Ashdown et al., Making Maps from Planck LFI 30GHz Data, 2007, submitted to A&A. [83] M. Reinecke et al., A&A 445, 373 (2006). [84] G. Rybicki and A. Lightman, Radiative Processes in Astrophysics (Wiley, New York, 1979). [85] J. Dunkley et al, arXiv:0811.4280vl (2008), submitted to ApJ. [86] A. Lazarian, Journal of Quantitative Spectroscopy and Radiative Transfer 106, 225 (2007). [87] D. P. Finkbeiner, M. Davis, and D. J. Schlegel, The Astrophysical Journal 524, 867 (1999). [88] J. Delabrouille et al., 2009, in preparation. 100 [89] E. L. Wright et al, ApJS 180, 283 (2009). [90] D. Babich and E. Pierpaoli, Physical Review D 77, 123011 (2008). [91] D. Herranz and J. L. Sanz, Matrix Filters for the Detection of Extragalactic Point Sources in Cosmic Microwave Background Images, 2008, arXiv.org:0808.0300. [92] R. Vio and P. Andreani. [93] P. C. Gregory, W. K. Scott, K. Douglas, and J. J. Condon, VizieR Online Data Catalog 103, 427 (1997). [94] M. R. Griffith, A. E. Wright, B. F. Burke, and R. D. Ekers, The Astrophysical Journal Supplement 90, 179 (1994). [95] H. Kuehr, A. Witzel, I. I..K. Pauliny-Toth, and U. Nauber, AA&AS 45, 367 (1981). [96] S. A. Trushkin, Bull.Spec.Astrophys.Obs.N.Caucasus 55, 90 (2003). [97] N. Jarosik et al, The Astrophysical Journal 148, 29 (2003). [98] R. S. Hill et al, The Astrophysical Journal 180, 246 (2009). [99] S. Meyer, Data Analysis for Scientists and Engineers (John Wiley & Sons, New York, 1975). [100] N. E. Groeneboom et al, Bayesian analysis of white noise levels in the 5-year WMAP data, 2009, arXiv.org:0904.2554. [101] I. K. Wehus, L. Ackerman, H. K. Eriksen, and N. E. Groeneboom, arXiv:0904.3998vl (2009), submitted to ApJ. [102] J. Hoftuft et al, Increasing evidence for hemispherical power asymmetry in the five-year WMAP data, 2009, arXiv.org:0903.1229. [103] N. E. Groeneboom and H. K. Eriksen, Bayesian analysis of sparse anisotropic universe models and application to the 5-yr WMAP data, 2008, arXiv.org:0807.2242. [104] F. K. Hansen et al, Power Asymmetry in Cosmic Microwave Background Fluctuations from Full Sky to Sub-degree Scales: Is the Universe Isotropic?, 2008, arXiv.org:0812.3795. 101 Vita CHARMAINE ARMITAGE-CAPLAN Contact Information Department of Physics University of Illinois at Urbana-Champaign 1110 West Green Street Urbana, IL 61801 USA Office: (217) 333-2807 E-mail: carmitag@illinois.edu www. astro. uiuc. edu/ ~ car mit ag Education University of Illinois at Urbana-Champaign, Urbana, Illinois USA • Ph.D., Physics, expected October 2009 • M.Sc, Finance, December 2007 • M.Sc, Physics, August 2005 University of British Columbia, Vancouver, British Columbia Canada • B.Sc, Combined Honours in Physics and Astronomy, May, 2003 Research Experience • Postdoctoral researcher in theoretical cosmology at the University of Oxford, Oxford, UK. Begins October 2009. • Graduate research assistant, Department of Physics, UIUC. Advised by Benjamin Wandelt. Summer 2003 - Summer 2009. • Undergraduate research for honours thesis, UBC, Vancouver, Canada. Advised by Douglas Scott. 2002-2003. 102 • Summer research assistant at the Canadian Institute for Theoretical Astrophysics (CITA), Toronto, Canada. Supervised by Richard Bond and Carlo Contaldi. Summer 2002. • Summer research assistant at CITA, Toronto, Canada. Supervised by Richard Bond and Carlo Contaldi. Summer 2001. Honors and Awards • UIUC: Ranked as Excellent Teacher, 2003 • CITA: Natural Sciences and Engineering Research Council (NSERC) Award, 2002 • CITA: NSERC Award, 2001 • UBC: Outstanding Student Initiative Award, 1997-2001 Teaching Experience • Teaching Assistant, How Things Work, UIUC, Spring 2005 • Teaching Assistant, Thermal and Quantum Physics, UIUC, Spring 2004 & Fall 2003 • Teaching Assistant, Physics & Astronomy Department, UBC, Spring 2002 Talks and Conferences • UIUC Center for Theoretical Astrophysics Seminar. Urbana, IL, 2006. Talk title: Probing the Early Universe through CMB Polarization Maps. • Preliminary Examination. Urbana, IL, April 10, 2006. Talk title: New Approaches to CMB Map-Making and Astrophysical Anomaly Detection. • American Physical Society April Meeting, Dallas, TX, 2006. Talk title: In Search of B-modes. • SF05 Cosmology Summer Workshop, July 5-22, 2005. Saint John's College, Santa Fe, New Mexico. 103 • American Physical Society April Meeting, Tampa, FL, 2005. Talk title: Deconvolution Map-Making for CMB Observations. • Planck HFI/LFI Joint Consortium Meeting, January 2005. Poster title: Deconvolution Map-Making for CMB Observations. • US Planck Algorithm Development Group Meeting, Pasadena, CA, October 27,2004. Talk title: Deconvolution Map-Making for CMB Observations. Publications • Armitage-Caplan, C. and Wandelt, B.D. "PReBeaM for Planck: A Polarized Regularized Beam Deconvolution Map-Making Method", ApJS 181, No.2, 533542 (2009). • Armitage, C , and Wandelt, B.D. "Generalized Beam Deconvolution MapMaking for Cosmic Microwave Background Observations", Physical Review D 70, 123007 (2004). Appearances in the Media • Interviewed on Channel 15 Sunrise This Morning News (Champaign, IL, 01/24/2007) to advertise the Saturday Astrophysics Honors Program. Computer Skills • Languages: Fortran, C++, Perl, OpenMP, MPI, IDL. • Operating Systems: Unix/Linux, Mac OS X, Windows. • Five years of experience in Linux-based high-performance parallel computing environment Professional Societies • American Physical Society • Union for Concerned Scientists 104 • Federation of American Scientists 105

1/--страниц