# Multifrequency analysis of cosmic microwave background radiation and radiation transport in simulations of reionization

код для вставкиСкачатьMultifrequency analysis of cosmic microwave background radiation and radiation transport in simulations of reionization Kevin M. Huffenberger A Dissertation Presented to the Faculty of Princeton University in Candidacy for the Degree of Doctor of Philosophy Recommended for Acceptance by the Department of Physics January, 2006 UMI Number: 3198041 Copyright 2006 by Huffenberger, Kevin Michael All rights reserved. UMI Microform 3198041 Copyright 2006 by ProQuest Information and Learning Company. All rights reserved. This microform edition is protected against unauthorized copying under Title 17, United States Code. ProQuest Information and Learning Company 300 North Zeeb Road P.O. Box 1346 Ann Arbor, MI 48106-1346 c Copyright 2006 by Kevin M. Huffenberger. All rights reserved. Abstract We explore two means for probing cosmology, through multifrequency microwave background measurements and through future observations of the epoch of reionization. We use multi-frequency information in first year Wilkinson Microwave Anisotropy Probe (WMAP) data to search for the Sunyaev-Zel’dovich (SZ) effect. We derive an optimal combination of WMAP cross-spectra to extract SZ, limiting the SZ contribution to less than 2% (95% c.l.) at the first acoustic peak in W band. Under the assumption that the removed radio point sources are not correlated with SZ, this limit implies σ8 < 1.07 at 95% c.l. The next generation of microwave telescopes will study the sky at high resolution, scales where both primary and secondary anisotropies are important. We focus on the Atacama Cosmology Telescope (ACT), simulating observations in three channels, and extracting power spectra in a multifrequency analysis. We find that both radio and infrared extragalactic point sources are important contaminants, but can be effectively removed given three (or more) channels and a good understanding of their frequency dependence. However, improper treatment of the scatter in the point source frequency dependence introduces a large systematic bias. The kinetic SZ effect corrupts measurements of the primordial slope and amplitude on small scales. We discuss the nonGaussianity of the one-point probability distribution function as a way to constrain the kinetic SZ effect, developing a method for distinguishing this effect. We explore the simulation of maps for ACT, their application to the ACT survey geometry, and filtering techniques to recover signals. iii Recent work suggests that cosmological fluctuations in reionization develop on scales of tens or hundreds of comoving megaparsecs. We build models of ionizing sources from simulations, concluding that a large-scale simulation will require radiation transport from a large fraction of the grid cells. Simulations at a reasonable resolution will have optically thick cells. We extend the work of Cen (2002) to build a parallel radiation transport algorithm which conserves photon flux to a few percent, works well with optically thick cells, and is insensitive to the number of sources in the simulation box. iv Acknowledgements People My advisor, Uroŝ Seljak, has always recommended challenging projects with a forwardlooking scope. For this I thank him. He always found time for me when I needed, but also allowed me space and time to work things out for myself. Most of all I thank him for providing, with his own work and the projects he’s helmed, examples of solid science-in-progress, which I could observe close up. I thank Lyman Page and Joe Fowler for guidance, and for encouraging my participation in the ACT collaboration. Working in the collaboration helped to motivate me. In addition, I thank Joe for being the least intimidating faculty member in the history of physics departments everywhere. Discussion with fellow students helps to dislodge my thinking. In this regard, I am in immeasurable intellectual debt to my office-mate and friend, Nikhil Padmanabhan. I also thank Chris Hirata, who typically knows the answer I need before I fully understand my question. I have enjoyed fruitful collaborations with Alexey Makarov, Hy Trac, and Neelima Sehgal, and I thank Renyue Cen, Eiichiro Komatsu, and Pengjie Zhang for providing me with tools and data to help my work. Graduate studies have been more challenging than I anticipated. I expected, though underestimated, the academic challenges. What I did not expect were the mental challenges to self-confidence and motivation which accompany them. My meditation instructor, Fay Elliot Moore, taught me strategies for addressing these challenges which v have proved invaluable. I thank my fianc ée Angie Knapp for walking with me during this time, during happy times and unhappy. Finally, I thank my Mom and Dad, who first encouraged my scientific curiousity. Support and Resources During my first three years of my studies I was supported by a graduate research fellowship from the NSF. I acknowledge the Legacy Archive for Microwave Background Data Analysis (LAMBDA), used for a portion of this research. Support for LAMBDA is provided by the NASA Office of Space Science. Parts of this research used a Beowulf cluster at Princeton University, supported in part by NSF grant AST-0216105. Software Three publicly available software libraries proved to be immensely useful to my work. Although I originally used routines from Numerical Recipes (Press et al., 1992), I have since transitioned almost entirely to the high-quality GNU Scientific Library (GSL) 1 . In my opinion, it offers equal or better functionality, a better interface for C programming, and a much more reasonable licensing scheme for scientific and academic work. I continue to recommend the Numerical Recipes book for its useful introduction and discussion of algorithms and numerical issues. The number of fast Fourier transforms I have performed with the Fastest Fourier Transform in the West (FFTW) 2 library is beyond my estimation. Finally, NASA’s High Energy Astrophysics Science Archive Research Center produces a very useful input/output subroutine library for the Flexible Image Transport System (CFITSIO)3 . 1 http://www.gnu.org/software/gsl/ http://www.fftw.org/ 3 http://heasarc.gsfc.nasa.gov/docs/software/fitsio/fitsio.html 2 vi Contents Abstract iii Acknowledgements v Contents vii 1 Introduction 1 2 Sunyaev-Zel’dovich effect in WMAP and its effect on cosmological parameters 7 2.1 Original abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 2.3 Cross-spectra combination for SZ . . . . . . . . . . . . . . . . . . . . . . . 11 2.3.1 Estimator and data . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 2.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.4 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 ACT simulations and power spectrum analysis 22 25 3.1 Original Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26 3.3 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.3.1 Hydrodynamical simulations . . . . . . . . . . . . . . . . . . . . . 29 3.3.2 Lensed CMB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 vii 3.3.3 Sunyaev-Zel’dovich effect . . . . . . . . . . . . . . . . . . . . . . . . 32 3.3.4 Kinetic Sunyaev-Zel’dovich effect . . . . . . . . . . . . . . . . . . . 35 3.3.5 Point sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35 3.3.6 Detector noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 3.3.7 Caveats and other signals not included . . . . . . . . . . . . . . . 41 3.3.8 Summary of simulations . . . . . . . . . . . . . . . . . . . . . . . . 43 3.4 Power spectrum analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 3.4.1 Multifrequency estimator . . . . . . . . . . . . . . . . . . . . . . . 46 3.4.2 Relation to Wiener filter . . . . . . . . . . . . . . . . . . . . . . . . 50 3.4.3 Frequency dependence and scatter . . . . . . . . . . . . . . . . . . 51 3.4.4 Measured signal power spectrum . . . . . . . . . . . . . . . . . . . 52 3.4.5 Measuring primordial amplitude and slope . . . . . . . . . . . . . 61 3.5 One-point distribution function analysis . . . . . . . . . . . . . . . . . . . 62 3.5.1 Signal one-point distribution functions . . . . . . . . . . . . . . . . 64 3.5.2 Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3.5.3 Separation by distribution function . . . . . . . . . . . . . . . . . . 67 3.6 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69 4 Progress towards simulation of reionization with radiation transport 71 4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 4.2 Source model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2.1 Halo population model . . . . . . . . . . . . . . . . . . . . . . . . . 73 4.2.2 Source luminosity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76 4.3 Intergalactic Medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4 Radiation transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80 4.4.1 Cen’s Flux Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.4.2 Photon Absorption . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.5 Tests and future directions . . . . . . . . . . . . . . . . . . . . . . . . . . . 87 viii A Semi-analytic simulations of the Sunyaev-Zel’dovich effect 89 A.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89 A.2 Analytic halo profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90 A.3 Halo catalog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92 A.4 Semi-analytic Sunyaev-Zel’dovich maps . . . . . . . . . . . . . . . . . . . 93 B Importing flat-sky simulations into the ACT geometry 95 C Signal filtering in ACT maps 98 C.1 Linear (Wiener) filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98 C.2 An example: point source identification . . . . . . . . . . . . . . . . . . . 101 D Cross-spectrum contaminant estimator 103 References 109 ix Chapter 1 Introduction In the modern understanding of the beginnings of our universe, inflationary expansion stretches microscopic quantum fluctuations to beyond the horizon scale of the universe. Thus cosmology simultaneously examines physical laws on the largest and smallest accessible scales, providing one of the few avenues which approach certain phenomena too large, too small, or too energetic to study otherwise. Cosmology is also peculiar among the disciplines of physics due to the lack of the laboratory. We cannot turn the dial and run the universe with different conditions to see what happens. We have only a single trial of the “history of the universe” experiment. Attempts to explain some apparently peculiar features of the universe, such as low or coincidently aligned multipoles at low-l in COBE and WMAP (Spergel et al., 2003; Efstathiou, 2003; de Oliveira-Costa et al., 2004; Copi et al., 2004, among others), may be continually frustrated for this reason. Other cases are not so bleak: with the assumption of homogeneity (observationally valid on the scales of a few hundred Mpc) and an appeal to the ergodic hypothesis, we treat small patches of the universe as independent trials. This dissertation explores a few probes of cosmology, which we will come to shortly. Careful study of cosmological probes is important precisely because we do not have a “universe lab” to fiddle with. We get at fundamental physics by proxy, through a cosmological probe. Now is a rapidly advancing time in the field, driven by advances in 1 2 the quality and amount of data from many types of observations. The result is that the accounting of the universe’s contents has improved greatly, with some bizarre revelations: we know now to a few percent the total density of the universe, and know that only five percent is in normal baryonic matter, protons and neutrons. The contents of the universe remain ninety-five percent unknown. 1 This ninety-five percent, composed of the dark matter and the dark energy, lays out one of the striking puzzles and challenges for cosmology in the near term. A variety of complementary measurements support this accounting. We compare the known physical scale of acoustic oscillations in the epoch of recombination to their apparent scale observed in the cosmic microwave background (CMB). Combined with a weak prior on the current expansion rate, this constrains the geometry of the universe to be nearly curvature-free. This in turn forces the total energy density of the universe to be close to the critical density for eventual collapse, which depends only on the current expansion rate, or Hubble constant: ρtot ∼ ρcrit = 3H 2 , 8πG (1.1) or Ωtot = ρtot /ρcrit ∼ 1 (Spergel et al., 2003). We can apportion this total density into several constituents. The microwave back- ground makes up the bulk of the radiation energy density in the universe. However, this is currently only a small fraction of the total, Ω r ∼ 5×10−5 . Models of Big Bang nucleosynthesis and comparisons to the purest stores of primordial material (untainted by nuclear reprocessing by stars) suggest that Ω b ∼ 0.05 (Olive et al., 2000). The ratio of the peaks in the microwave background power spectrum confirms this measurement (Page et al., 2003). Thus we reach the conclusion that ninety-five percent of the universe is non-baryonic. The combination of CMB power spectrum with measurement of the distance-redshift relation with type Ia supernova tells yet a stranger story. Of this ninety-five percent, only a fraction is in matter which is gravitationally attractive. This dark matter accounts for ΩDM ∼ 0.25. The other portion, the dark energy, is grav1 Upon hearing this fact, an acquaintance recently asked, “How can you sleep at night?” 3 itationally repulsive, and accounts for the balance, Ω DE ∼ 0.7 (Perlmutter et al., 1999; Riess et al., 1998). This repulsive gravity drives the expansion in the present-day universe to accelerate. More probes confirm this basic story. Dynamical measurements of galaxy rotation curves and galaxy clusters have long indicated the presence of dark matter (Sofue & Rubin, 2001). Galaxy and Lyman-α forest measures of matter clustering are compatible with the primary CMB power spectrum. Combined with it, they limit the slope and the running of the primordial spectrum (e.g. Spergel et al., 2003; Seljak et al., 2005), constraining models of the universe’s inflationary epoch at early times. The dark energy, cause of the accelerated expansion, remains one of the most outstanding unknowns in the current cosmological model. A cosmological constant in Friedman’s equations would have this effect, perhaps caused by a vacuum energy density. Other proposals exist, such as scalar fields (quintessence, k-essence, phantom energy, etc.) or more elaborate notions (like brane collisions). We can parameterize pressure-density relation—the equation of state—of a perfect fluid by p = wρ. Normal gas has w > 0 and pressure-less dust has w = 0. The equation of state of a cosmological constant is constant with w = −1, while the other types of dark energy this relation- ship varies in time; the dark energy is dynamical. For constant w, observations imply w < −0.8 (Spergel et al., 2003). Considerable effort will be expended over the next few years to try to measure the dark energy equation of state, and any change in it over time. These measurements rely of effects on scale parameter evolution of the universe, measured directly (for example, with baryon wiggles in the matter density) or through the impact on the growth of structure. Useful probes of cosmology may be qualitatively placed on a continuum based on how they are calibrated. On this continuum, the most straightforward probes depend only on well known physics, calibrated by experiments in a terrestrial lab. An example of a probe of this type is the cosmic microwave background perturbation, which depends on atomic physics, Compton scattering, and linear-order general relativity, 4 all of which are well tested. At the opposite extreme are phenomena which are predictable, but not well understood physically. An example is the application of Type Ia supernovae as standard candles. The physics of the supernova detonation is not fully understood, yet calibration of the luminosity based the light curve enables measurement of the luminosity distance redshift relation. In some sense, it is good fortune, rather than deep understanding, which has led to their great success as cosmological probes. There is a tendency for future probes to depend on complicated physics (for example, utilizing gravitationally collapsing or collapsed objects) or to be hidden behind complicated foregrounds. The probes must be carefully characterized to understand their properties. We can hope that features of these probes turn out as well as for the supernovae. The following chapters focus on a few promising probes which are currently under study. Some of these probes will be observed with multifrequency observations of the cosmic microwave background. One of these is the Sunyaev-Zel’dovich (SZ) effect, which is a well established scattering effect of CMB light from gas in high-density structures. For this reason, it is a potentially powerful probe of the growth of structure, capable of constraining both the amplitude of matter fluctuations (parameterized by σ8 ) and the dark energy (through the growth’s dependence on the scale factor). Chapter 2 places a limit on the SZ contribution to the first year data from the Wilkinson Microwave Anisotropy Probe, countering some claims that this contribution is very large. Future observations of the CMB temperature (and later, polarization) at scales of a few arcminutes can probe fundamental physics, for example by detecting a non-zero neutrino mass (Kaplinghat et al., 2003). Three multifrequency experiments, APEXSZ,2 SPT 3 (Ruhl et al., 2004), and ACT4 (Kosowsky, 2003), will survey hundreds of square degrees at arcminute resolution in the range 100–300 GHz, beginning in 2005– 2 http://bolo.berkeley.edu/apexsz/ http://astro.uchicago.edu/spt/ 4 http://www.hep.upenn.edu/act/ 3 5 2007. These will constrain the CMB and provide detections of secondary anisotropies. On small scales, the CMB is dominated by non-noise contaminants, including SZ clusters and extragalactic point sources. Cleaning these out is essential, but tricky since the CMB fluctuations are Gaussian and the contaminants are not. No simply defined optimal cleaning method exists. In chapter 3, we build simulations of ACT observations, and discuss measurement of the power spectrum and certain non-Gaussianities. The conclusions are unsurprising: controlling the systematics caused by the contaminants crucially determines the reliability of the power spectrum measurement. The other probes we will examine relate to the epoch of reionization, the process by which light from early galaxies ionized the previously neutral intergalactic medium. This is an opening frontier in cosmology. Preliminary observations of this era, of the total optical depth to Thomson scattering since recombination (Kogut et al., 2003) and of the Gunn-Peterson (1965) trough in z ∼ 6 quasar spectra (e.g. Fan et al., 2004), indicate that the epoch of reionization must have been extended. A number of explanations have been offered, positing mechanisms of early partial (or even complete) ionization (e.g. Wyithe & Loeb, 2003; Chiu et al., 2003; Cen, 2003). These ideas lend hope that the epoch of reionization holds a rich store of information about diverse topics such as large scale fluctuations, gravitational collapse to early halos, chemistry and radiative processes in the intergalactic medium (IGM), and star formation in low metallicity environments. If so, the reionization epoch will provide fertile ground for study for many years to come. In chapter 4 we explore the process of reionization, and we develop a parallel radiation transport algorithm with features appropriate for its study in numerical simulations. Several appendices describe smaller projects, each related to multifrequency microwave measurement. Appendix A describes a prescription for simulating SZ based on halo catalogs from dark-matter-only numerical simulations. Appendix B briefly describes a method for importing numerical simulations from a cubical box into the ACT survey geometry. Appendix C describes a method to filter multifrequency maps, 6 like those from ACT, to find objects such as point sources and clusters. Appendix D describes a way to estimate the contamination in multifrequency CMB cross-power spectra. Chapter 2 Sunyaev-Zel’dovich effect in WMAP and its effect on cosmological parameters This chapter originally appeared, in slightly modified form, in Huffenberger, Seljak, & Makarov (2004). 2.1 Original abstract We use multi-frequency information in first year WMAP data to search for the SunyaevZel’dovich (SZ) effect. WMAP has sufficiently broad frequency coverage to constrain SZ without the addition of higher frequency data: the SZ power spectrum amplitude is expected to increase 50% from W to Q frequency band. This, in combination with the low noise in WMAP, allows us to strongly constrain the SZ contribution. We derive an optimal frequency combination of WMAP cross-spectra to extract SZ in the presence of noise, CMB, and radio point sources, which are marginalized over. We find that the SZ contribution is less than 2% (95% c.l.) at the first acoustic peak in W band. Under the assumption that the removed radio point sources are not correlated with SZ this limit 7 8 implies σ8 < 1.07 at 95% c.l. We investigate the effect on the cosmological parameters of allowing an SZ component. We run Monte Carlo Markov Chains with and without an SZ component and find that the addition of SZ does not affect any of the cosmological conclusions. We conclude that SZ does not contaminate the WMAP CMB or change cosmological parameters, refuting the recent claims that they may be corrupted. 2.2 Introduction Wilkinson Microwave Anisotropy Probe (WMAP, Bennett et al., 2003a) observations of the Cosmic Microwave Background (CMB) have ushered in a new era of high precision observational cosmology. Such a tremendous increase in data quality requires a corresponding increase in the care that goes into the data analysis and interpretation. One of the lingering concerns surrounding the analysis is the residual effect from additional sources of anisotropies. A very prominent candidate among these is the Sunyaev-Zel’dovich effect (SZ). Electrons in hot gas scatter photons and distort the blackbody spectrum of the CMB. Galaxy clusters, where gas is the hottest, contribute the bulk of the effect, shock-heating the gas in their potential wells. This effect is now well established observationally: radio and microwave band observations pointed at known clusters routinely yield SZ detections. The scattering preferentially raises the energy of the CMB photons, but the number of scatterings is low, so the process never achieves thermal equilibrium. Therefore SZ appears as a CMB temperature decrement at low frequencies and as an increment at high frequencies. This frequency dependence is well known: it is constant in the Rayleigh-Jeans (RJ) regime of the blackbody spectrum, and is universal (independent of gas temperature or density) for non-relativistic electrons. In the channels of WMAP, SZ is a temperature decrement. Table 2.1 provides WMAP’s frequency bands, and the SZ frequency dependence in those bands. The K and Ka bands are the most heavily polluted by galactic contamination, so the best prospects for identifying SZ in WMAP are in the differencing assemblies of the upper 9 Band K Ka Q V W ν (GHz) 23 33 41 61 94 ∆T SZ /(T CMB y) −1.97 −1.94 −1.91 −1.81 −1.56 Cl /ClRJ 0.972 0.945 0.915 0.820 0.611 Table 2.1: SZ contribution in the bands of WMAP. We note the band name, frequency, temperature perturbation relative to the comptonization parameter y, which does not depend on frequency, and power spectrum relative to the power in the Rayleigh-Jeans regime. In the WMAP bands, SZ is a temperature decrement. In the Rayleigh-Jeans regime, ∆T SZ /(T CMB y) = −2. three bands: two assemblies in Q, two in V, and four in W. While it is usually argued that SZ is indistinguishable from the CMB in the WMAP channels this is actually not so: in CMB temperature units the SZ power spectrum increases by 50% from W to Q channel and in W it is only 61% of the RJ power. We will show that this frequency dependence suffices to place strong constraints on SZ using WMAP data alone, a consequence of the remarkably low noise in WMAP data. In the literature, several groups have attempted to identify SZ in WMAP data using cross-correlations with other tracers of large scale structure (LSS). In Bennett et al. (2003b), they build an SZ template from the XBACs catalog of x-ray clusters (Ebeling et al., 1996), and fit for the amplitude of SZ, arguing for a 2.5σ detection. In Afshordi et al. (2003), they compute the cross-power spectra of the 2MASS extended source catalog (Jarrett et al., 2000) with WMAP’s Q, V, and W bands, then fit for the amplitude of the SZ signal, arguing for a 3.1–3.7σ detection of SZ, depending on their mask. In Diego et al. (2003), they compute the cross-power spectrum of the ROSAT diffuse x-ray background maps (Snowden et al., 1997) and a WMAP (Q−W) difference map, finding no detection. In Fosalba et al. (2003) they compute the cross-correlation function of SDSS DR1 (Abazajian et al., 2003) galaxy survey data with the V band and fit for the amplitude of SZ, finding a 2.7σ detection. In Fosalba & Gazta ñaga (2004) they compute the the cross correlation of the APM galaxy survey (Maddox et al., 1990) with 10 V band, W band, and the map of Tegmark et al. (2003), then fit for the amplitude of SZ, arguing for a 2.6–3.1σ detection, depending on scale. In Hern ández-Monteagudo & Rubiño-Martı́n (2004) they build templates for SZ from several optical and x-ray cluster surveys, then fit for the amplitude of these templates using maps from the W band and Tegmark et al. (2003), finding no detection for optical clusters and 2–5σ detections for x-ray clusters. Finally, in Myers et al. (2004) they compute the mean temperature in concentric rings about groups and clusters. They employ a friends-of-friends algorithm to identify galaxy concentrations in the APM galaxy survey (Maddox et al., 1990) and 2MASS extended source catalog (Jarrett et al., 2000) and use the ACO cluster catalog (Abell et al., 1989). They note a temperature decrement which they attribute to SZ. They interpret this decrement to extend to large scales, ∼ 1 degree, although the covariance of their bins is unclear. Redshift z < 0.2 clusters dominate their sample. Their most ex- treme model assumes that extended SZ emission is representative of the temperature and spatial clustering of gas to z = 0.5. In this model the SZ power is 30% of the first acoustic peak of the CMB. Thus they conclude the integrity of the WMAP cosmological parameter fits may be compromised. In this paper we seek to place limits on SZ from WMAP data alone, thus avoiding any uncertainties in connecting the results of WMAP-LSS cross-correlation to those of WMAP auto-correlation. The cross-correlation methods used before require a model for relating SZ temperature perturbations to some proxies for SZ, which are nonlinear structures such as clusters, and this connection can be quite uncertain. Our method is less model dependent. The downside is that WMAP-only methods sacrifice signal to noise because they do not focus on matter concentrations, where SZ should be strongest, so our method is not optimized for obtaining an SZ detection. However, our main goal is to investigate the amount to which the CMB analysis may be compromised by the residual SZ component, for which we need analysis as model independent as possible. 11 Our principal method, described in section 2.3, constructs a linear combination of the WMAP Q, V, and W band cross-power spectra which maximizes the contribution of SZ, while at the same time minimizing the radio point source and CMB contribution. We use this linear combination to fit for the SZ power spectrum amplitude, using a spectrum shape from halo model calculations (Komatsu & Seljak, 2002). In Huffenberger et al. (2004) we supplement this method by fitting for cosmological parameters from the WMAP temperature power spectrum using the Markov chain Monte Carlo method, allowing for an additional SZ component in the power spectrum with a free amplitude to be determined by the data. The two methods give consistent results, with the latter giving somewhat weaker constraints. Conclusions are presented in section 2.4. 2.3 Cross-spectra combination for SZ We begin our discussion with a simple example. In the absence of noise and point sources, the exact SZ power spectrum could be computed from the difference of any two cross-spectra from different bands. Take the maps from some of the WMAP differencing assemblies, which provide independent measurements of the sky, and compute their cross-power spectra. For example, take Q1 and Q2 from Q band, and W1 and W2 from W band. A cross-spectrum can be denoted e.g. (Q1Q2) l . One may find combinations of cross spectra which eliminate CMB and return the SZ power spectrum, using the coefficients from table 2.1. One combination which gives the SZ power in this idealized case is the difference (Q1Q2)l − (W1W2)l = 0.30 ClSZ,RJ (ideal). (2.1) Using other bands, other combinations of cross spectra would also yield the SZ spectrum, with different coefficients. The important things to note are that the coefficient is not very small and the CMB cancels exactly in this combination: any CMB cosmic variation in the Q band is the same as in the W band, and is subtracted out. In the 12 presence of instrument noise, this difference does not give the exact SZ power spectrum, but a noisy estimate for it. Given a power spectrum shape, we may estimate its amplitude by summing together different l bins, each weighted to give the proper normalization and to emphasize bins with high signal to noise. Our SZ estimator works in this fashion, except that in the presence of point sources, QQ−WW is a biased estimate: this combination will eliminate CMB, but the additional point source component with a different spectral dependence will remain. By including other cross-spectra, we can form a linear combination of the cross-spectra bins to eliminate both CMB and point sources, yielding an unbiased estimate for the SZ amplitude. In the following we find such a combination, and apply it to WMAP data. 2.3.1 Estimator and data In this section, we introduce our estimator and list the data we use. Our estimator, derived in appendix D, is a generalized version of the point source estimator of Hinshaw et al. (2003). We want to generate an estimate for SZ from the WMAP cross-spectra. We postulate that the data is the sum of the CMB and two contaminants, radio point sources and SZ. We marginalize point sources and estimate SZ in our main case, which we illustrate in detail below. We consider other cases also, but for brevity omit their details, which are similar. We can write the cross-spectra as hCli i = Cli,CMB + Cli,src + Cli,SZ , (2.2) showing explicitly the contribution from each part of the signal. Here the multipole bin is denoted by l and the pair of differencing assemblies in the cross-correlation is denoted by i = i1 i2 = W1W2, Q1V1, etc. No auto-power spectra are included, so we need not worry about noise subtraction. The scales we are most interested have little foreground contamination, and data we use has had a foreground model subtracted out, so we include no foreground contribution. Here and throughout, we assume point sources are uncorrelated with SZ, so we have included no cross term between point 13 sources and SZ. We discuss the implications of this assumption in this section below and again in the conclusions. We will marginalize over the CMB spectrum, which we denote by C lCMB . The window functions for each differencing assembly pair are w = {w lli 0 }. We boldface w because later we will think of it as a matrix. Therefore the contribution to the crossspectrum from the CMB is Cli,CMB = X wlli 0 ClCMB . 0 (2.3) l The spectra we use are in temperature units, and have already been beam-deconvolved, so the window functions wlli 0 = δll0 are trivial. However, it is necessary to keep the window functions explicit in the manipulations that follow. We denote the amplitude of the point source power spectrum by A. This amplitude relates to the cross-spectra via the frequency and shape dependence S = {S li }. Later we will think of S as a vector. We take radio sources to have a white noise spectrum with power law frequency dependence, given by: Cli,src = Sli A β β ν i2 ν i1 i , Sl = ν0 ν0 (2.4) where cross-spectrum i has channels at ν i1 and νi2 . The units of A are temperature squared. Well-resolved point sources have already been masked from the maps before the evaluation of the cross-spectra, so A represents unresolved sources only. Bennett et al. (2003b) found β = −2.0 and ν0 = 45 GHz for the resolved sources. Following Hinshaw et al. (2003), we take the same for the unresolved sources. We describe amplitude of SZ with B, the ratio of the SZ power spectrum and the predicted spectrum, assuming they have the same shape. Thus the SZ amplitude is dimensionless, and has a theoretically predicted value B = 1 for σ 8 = 0.9 using the halo models of Komatsu & Seljak (2002). At the scales relevant here, Komatsu & Seljak (2002) found that the halo model predictions for the SZ power show little dependence on halo concentration and are insensitive to the pressure profile of the cluster core. 14 The effect of radiative cooling is found to be modest. Preheating of cluster gas is found unlikely to substantially enhance the SZ power without compromising limits on the mean comptonization. The power is sensitive to the mass function, but not enough to make a significant change in the derived limits on σ 8 . In our notation, the amplitude B relates to the power spectra by the frequency dependence and shape for SZ, which we label Z = {Zli }. The frequency dependence of a temperature perturbation due to SZ is ∆T SZ ∝ −2f (x), T CMB (2.5) where the frequency dependence relative to the RJ regime is given by f (x) = 2 − x/2 , tanh(x/2) (2.6) where x = hν/kB TCMB . Note f → 1 in the RJ limit x → 0. Thus the SZ contribution to the cross-spectrum is: Cli,SZ = Zli B (2.7) Zli = ClSZ,RJ f (νi1 )f (νi2 ), where the spectrum shape ClSZ,RJ is from the halo model prediction in the RJ regime using σ8 = 0.9, as shown in Figure 2.3. The shape of SZ is roughly C lSZ,RJ ∝ l−1 , although it becomes slightly steeper for l greater than a few hundred. We organize the binned cross-spectra C li into a data vector D = {Cli }. We use a Gaussian model for the likelihood L of the power spectrum: 1 † −1 L ∝ exp − [D − hDi] Σ [D − hDi]) , 2 (2.8) where the covariance Σ = h(D − hDi)(D − hDi) † i can be written as Σ = {Σii ll0 }. We 0 derive B̄, an unbiased estimator for B, and its covariance Σ B in appendix D. The main result is: B̄ ≡ (Z† FZ)−1 Z† FD ΣB ≡ (Z† FZ)−1 , (2.9) 15 where we have defined the auxiliary matrices † −1 F ≡ E − ES S ES S† E −1 w † Σ−1 E ≡ Σ−1 − Σ−1 w w † Σ−1 w (2.10) In this notation, we consider D, S, and Z as column vectors with single index il, and w as matrix with indices il and l 0 . Σ, E, and F are matrices with indices il and i 0 l0 . This estimator marginalizes CMB and point sources, which is the most conservative treatment, since it assumes we know nothing about these components. A more aggressive treatment would be to estimate all 3 components simultaneously, but we defer this approach to a future analysis. Note that B̄ is a linear combination of the cross-spectra D, with weights (Z† FZ)−1 Z† F. The data we use consist of the 28 cross-power spectra from the eight differencing assemblies in the WMAP Q, V, and W bands. These spectra are provided at the Legacy Archive for Microwave Background Data Analysis. 1 A galactic foreground model (Bennett et al., 2003b) has already been subtracted out. WMAP’s temperature power spectrum is a linear combination of these 28 cross-spectra (Hinshaw et al., 2003). We bin the spectra in l, accounting for the number of modes at each l. The Legacy Archive does not provide the cross-spectrum covariance, which we need for Σ in equation (2.10), so we estimate the covariance from the data. We assume the covariance is diagonal in l, which is a good approximation (Hinshaw et al., 2003). Then in a single bin, we use the dispersion of the cross-spectra about WMAP’s combined temperature power spectrum to estimate Σ for that bin. (For this purpose we must first un-correct the combined spectrum for the point source contribution.) This procedure to obtain the covariance works best if the power spectrum variance does not change much within a bin and if the bin contains enough l’s to get a low-noise estimate of the variance. Our bin width of ∆l = 40 is fairly narrow, and gives a covariance which is numerically stable in our subsequent calculations. This technique does not account for the cosmic variance of the CMB, but as we note in the derivation in the appendix, our estimator is 1 http://lambda.gsfc.nasa.gov/product/map/ang power spec.cfm 16 insensitive to this source of variance. This makes sense because the CMB is completely projected out, as in our QQ−WW example at the beginning of section 2.3. As a test of our estimator, we repeat the estimate of Hinshaw et al. (2003) of the power spectrum amplitude of unresolved point sources. In this case we neither marginalize nor estimate SZ, but assume B = 0 in equation (2.7). We find a point source amplitude of A = 0.016 ± 0.001 µK 2 , which is roughly consistent with their quoted value of A = 0.015 µK2 . This gives us confidence in our estimator, despite our covariance matrix constructed from the data. 2.3.2 Results We show our estimates for SZ on a bin-by-bin basis, before combining different multipole bins to improve statistics. We marginalize over the CMB and the point source amplitude. Our best estimate for the binned SZ power spectrum is shown in Figure 2.1. Immediately we can see that the data do not tolerate a large SZ contribution. Next we turn to our main case, where we combine all the bins together and estimate the SZ amplitude. Figure 2.2 shows the weight (Z † FZ)−1 Z† F for each cross-spectrum in the amplitude estimator. The weights have been divided into groups based on the bands. For plotting, we sum the weights in each group. The bulk of the weight comes from l = 100 to l = 400. At lower l there are fewer modes and the SZ contribution is low. At higher l the statistics are limited by detector noise. The weights are difficult to interpret heuristically. From the frequency dependence, the strength of SZ decreases in order of QQ, QV, VV, QW, VW, and WW. We would expect our combination of cross-spectra would have the SZ-stronger bands minus the SZ-weaker bands. So it is easy to understand the positive weights of QV and QW and the negative weights of VW and WW. Point sources are strongest in Q band, so a negative weight for QQ makes sense in terms of a point source correction. QW and VW are positive except for bins which are strong in QV. These negative bins may also represent a point source correction 17 6000 5000 CMB 2 l(l+1) Cl / 2π (µK ) 4000 3000 2000 1000 0 SZ -1000 -2000 -3000 -4000 0 100 200 300 400 500 600 l Figure 2.1: Reconstruction of the SZ power spectrum in the RJ regime, in bins of ∆l = 40. There is no visible SZ detection, which would have been positive for this combination. The SZ error bars are correlated at less than the 1%–15% level. For comparison, we have included the WMAP combined temperature power spectrum, binned the same way, the WMAP best-fit theoretical CMB power spectrum, and a theoretical prediction of SZ in the RJ regime for σ 8 = 0.9 calculated from a halo model (thick dotted), which is barely distinguishable from zero (thin dotted) on this plot. 18 600 400 200 0 -200 -400 QQ -600 600 400 200 0 -200 -400 QV -600 600 400 200 0 -200 -400 QW -600 600 400 200 0 -200 -400 VV -600 600 400 200 0 -200 -400 VW -600 600 400 200 0 -200 -400 WW -600 0 100 200 300 400 500 600 700 800 900 l Figure 2.2: The contribution to the SZ amplitude estimator B̄ from the cross-spectra. The horizontal axis is the multipole l. The vertical axis show the weight of each spectrum in µK−2 . We have grouped the cross-spectra by band as noted. Within each group, we have summed the weights of the individual spectra. In the top frame, we have included the CMB power spectrum (arbitrary units) to give scale for l. 19 For the amplitude of SZ relative to the expected RJ amplitude from the halo model with WMAP parameters, our estimator gives us B̄ = −0.042 ± 1.685. The error is large, and is consistent with both no SZ and the expected B = 1. Including only the physical region B > 0, we integrate the likelihood until we include 95% of the probability. From this we quote an upper limit of B < 3.3 (95%). (2.11) We plot this limit, along with some models of the SZ power spectrum in Figure 2.3. At the first acoustic peak, SZ is less than 3% of the CMB in the RJ regime and less than 2% of CMB in W band, on which cosmological parameter estimation is heavily based. We find no evidence for a large SZ contribution. Assuming C lSZ ∝ (σ8 )7 , this gives a limit of σ8 < 1.07 at 95% confidence. To this one should add an additional modeling uncertainty at the level of 10% based on the comparison of predictions from different simulations (Komatsu & Seljak, 2002). There is a caveat in the upper limit derived here in that we are working with power spectra based on masked maps in WMAP, with more than 200 radio point sources removed. If the SZ signal is correlated with these point sources, which could happen if these radio sources sit in massive clusters (Holder, 2002), then more SZ may have been removed than expected based on the sky fraction of the mask. This would only affect the σ 8 limits and not the SZ contamination on the primary CMB in the WMAP power spectra. Since we are mostly concerned with the latter we do not explore this issue further. Note that the upper limit on σ 8 is already comparable to the predicted value based on detections from CBI and BIMA, which gives σ8 ∼ 0.95 − 1.05 (Readhead et al., 2004; Komatsu & Seljak, 2002). It is remarkable that WMAP first year data have sufficient sensitivity to place constraints on the SZ amplitude comparable to other small scale surveys. In the next application we jointly estimate the point source and SZ amplitude, rather than directly marginalizing over the point sources. We find the two parameters to be somewhat correlated (Figure 2.4), but not to the point of allowing very large SZ amplitude. An extremely strong SZ power is not allowed. 20 104 2 l(l+1) Cl / 2π (µK ) CMB 103 102 k Mar 5 in 9 ha ov C C % a 95 tr pec s-s ros σ8=1.1 t mi % li t limi SZ σ8=1.0 σ8=0.9 100 l Figure 2.3: Limit on SZ power spectrum. We present a 95% upper limit on the amplitude of the SZ power spectrum from a combination of WMAP band cross-spectra. We show the temperature fluctuation in the Rayleigh-Jeans regime. For comparison we display the WMAP best fit power spectrum and SZ power spectra from a halo model calculation using σ8 = 0.9, 1, 1.1. B ωb × 102 Ωm ωcdm τ σ8 h ns no SZ 0 +0.28 2.35+0.14 −0.13 −0.26 +0.07 +0.15 0.245−0.06 −0.10 +0.033 0.111+0.016 −0.015 −0.029 +0.07 +0.10 0.19−0.08 −0.14 +0.25 0.88+0.12 −0.11 −0.20 +0.06 +0.12 0.74−0.05 −0.09 +0.07 0.99+0.04 −0.04 −0.07 with SZ, no SZ prior < 7.1 (95%) +0.47 2.51+0.21 −0.18 −0.33 +0.07 +0.15 0.234−0.06 −0.10 +0.033 0.111+0.016 −0.016 −0.030 +0.07 +0.10 0.20−0.08 −0.14 +0.24 0.86+0.12 −0.11 −0.22 +0.07 +0.14 0.76−0.06 −0.11 +0.07 0.99+0.04 −0.04 −0.07 with SZ prior < 2.9 (95%) +0.30 2.43+0.15 −0.15 −0.29 +0.07 +0.15 0.243−0.06 −0.10 +0.033 0.111+0.016 −0.016 −0.031 +0.07 +0.10 0.19−0.08 −0.14 +0.23 0.87+0.12 −0.11 −0.22 +0.06 +0.13 0.75−0.05 −0.10 +0.07 0.99+0.04 −0.04 −0.07 Table 2.2: The first two columns contain the median value and 1- and 2σ constraints on cosmological parameters for two MCMCs without and with a Sunyaev-Zel’dovich component in the Cl power spectrum from WMAP data alone. For both chains there was an imposed prior of τ < 0.3. The third column shows the constraints when the limit on SZ from our cross-spectrum estimator is applied as a prior. 21 0.0175 Expected SZ 68% 2 Point source amplitude A (µ K ) 0.018 0.017 0.0165 Best fit 0.016 0.0155 WMAP est. unresolved src. 0.015 0.0145 0.014 -3 -2 -1 0 SZ amplitude B 1 2 3 Figure 2.4: A joint estimation of the SZ power spectrum amplitude and the point source spectrum amplitude. The ellipse encloses 68% of the likelihood. The best fit point is shown by a “+”. The vertical line at B=1 shows the theoretically predicted SZ amplitude. The horizontal line shows WMAP’s value for point sources, which does not take into account SZ. Because the power spectrum is a positive quantity, note that portions of the plane where either axis is negative are unphysical. 22 2.4 Conclusions We have analyzed the WMAP power spectra in two ways and found them to be free of serious contamination from SZ. The more powerful and less model dependent of the two methods is to use multi-frequency information: by combining the WMAP crossspectra, we limit at 95% confidence the amplitude of SZ to be below 2% of the CMB at the position of first peak in W band. By searching for an SZ-shaped component in the WMAP combined spectrum with a Markov chain Monte Carlo we also do not find any evidence of a signal, but we can only set a weaker limit. Combining the analyses, we show that the cosmological parameters are not affected by the SZ within the range allowed by the multi-frequency analysis. Our analysis uses WMAP data alone, so is less model dependent than cross-correlations with external data sets, at the cost of lower signal to noise. The mean redshift which contributes to the SZ power spectrum tends to be higher than for the galaxies in crosscorrelation methods, which additionally may have an unknown bias. Translation of the results of the cross-correlation methods into a measure of the SZ power spectrum is difficult. There are several possible improvements in the analysis that could further tighten the limits and may be applied to the second year data once they become available. In addition to having lower noise, the two year data are also expected to have better control of systematics such as beam uncertainties and noise correlations. First, we would like to improve our method for obtaining the cross-spectrum covariance matrix. The estimate for the covariance is noisy, which makes our weights for combining the power spectra noisier than they should be. If we double the size of the bins we average over to obtain the covariance, our limit changes slightly (B < 4.3 opposed to B < 3.3 at 95%) and the point source estimate is virtually unchanged. However, it is impossible to tell if this effect is due to a covariance matrix with less noise (which is better) or due to weights whose wide bins combine high- and low-signal-tonoise multipoles less optimally (which is worse). This can be improved if the WMAP 23 second year data release includes the cross-spectrum covariance matrix. The second possible improvement of our analysis is to improve the frequency dependence of point sources or the WMAP beam deconvolution, both of which can bias our estimate. Our estimator is only unbiased if the models one applies to it are faithful to the data. The models we have used are the best available, and this information will improve in the second year data. It seems unlikely that these improvements will strongly modify the conclusions in this paper. The third improvement is to perform the analysis on the maps which do not have point sources excised. As discussed above, masking of resolved point sources from the WMAP data may reduce the SZ signal if the two are correlated. This would weaken our limit on σ8 from the absence of SZ signal, but does not change our conclusions about the effect of SZ on the CMB and cosmological parameters. Still, given that the current limits are reaching the levels of interest it is worth exploring this possibility further. The fourth improvement is to verify the assumed shape for the SZ power spectrum. This is currently based on halo models, since hydrodynamic simulations cannot simulate sufficiently large volumes to make predictions reliable. If the power spectrum has a radically different shape, our limit is hard to interpret, but the level of contamination on primary CMB is likely to remain unchanged, as is clear from Figure 2.1. As an example we estimated the SZ in bins of ∆l = 120. In the bin containing the first acoustic peak, we limit the SZ power to 4.9 times the halo model power predictions at 95% confidence or 4% contamination in W. It is clear that as long as the SZ spectrum is smooth the limits on the contamination are not going to change significantly over the range with best sensitivity at 100 < l < 300. Finally, our procedure to marginalize over CMB and point sources is conservative and the errors could be further reduced if a joint estimation is used instead. We know something more about the CMB power spectrum than just its frequency dependence, and in this case we can use this information. We also note that a relatively trivial 24 modification to WMAP’s method for estimating the power spectrum could eliminate any contaminant with an SZ-like shape and frequency dependence, much in the same way that they eliminate point sources. However, the fact that we find no residual SZ effects implies that this procedure is not really required and justifies the treatment of the WMAP team in ignoring SZ. Even with future data detecting the SZ power spectrum from WMAP alone may be difficult. If the bulk of the covariance is from noise, it would take WMAP roughly 3 years to get a 1σ detection of SZ using the cross-spectrum combination method of this paper, if σ8 = 0.9. Elimination of non-noise sources of covariance, such as uncertainties remaining in the beams, may help. If σ 8 = 1 then a 2σ detection should be possible with 4 year data and one should start seeing a hint of the signal already with 2 year data. Theoretical predictions remain somewhat uncertain, so it is possible that these predictions have a systematic uncertainty of around 10%, but it is clear that it will be worth looking for SZ in future WMAP releases. Chapter 3 ACT simulations and power spectrum analysis This chapter originally appeared as Huffenberger & Seljak (2005). 3.1 Original Abstract A new generation of instruments will reveal the microwave sky at high resolution. We focus on one of these, the Atacama Cosmology Telescope, which probes scales 1000 < l < 10000, where both primary and secondary anisotropies are important. Including gravitational lensing, thermal and kinetic Sunyaev-Zel’dovich (SZ) effects, and extragalactic point sources, we simulate the telescope’s observations of the CMB in three channels, then extract the power spectra of these components in a multifrequency analysis. We present results for various cases, differing in assumed knowledge of the contaminating point sources. We find that both radio and infrared point sources are important, but can be effectively eliminated from the power spectrum given three (or more) channels and a good understanding of their frequency dependence. However, improper treatment of the scatter in the point source frequency dependence relation may introduce a large systematic bias. Even if all thermal SZ and point source effects 25 26 are eliminated, the kinetic SZ effect remains and corrupts measurements of the primordial slope and amplitude on small scales. We discuss the non-Gaussianity of the one-point probability distribution function as a way to constrain the kinetic SZ effect, and we develop a method for distinguishing this effect from the CMB in a window where they overlap. This method provides an independent constraint on the variance of the CMB in that window and is complementary to the power spectrum analysis. 3.2 Introduction Several survey telescopes (SPT1 , APEX2 , ACT3 ) will open a poorly probed regime of the microwave background: arcminute scales at frequencies up to a few hundred GHz. The sensitivities of these instruments will be a few micro-Kelvin. These efforts will provide a large amount of high-quality data, but face a different set of challenges than lower resolution surveys, such as WMAP 4 . For lower resolution experiments the primary anisotropies of the CMB are the principal signal. The major contaminant is the galaxy, whose emission is diffuse and drops off rapidly away from the galactic center. By contrast, for higher resolution experiments, secondary anisotropies dominate at small scales. Assuming the observations are away from the galaxy at high frequencies, the major contaminant is extragalactic point source emission. For a wide range of interesting scales, point sources contribute substantially more power than the instrument noise. The statistical properties of the primary CMB differ from those of these signals. Primary anisotropies, from the surface of last scattering, dominate the CMB on degree scales. The fluctuations are Gaussian, so the power spectrum sufficiently describes the anisotropy. Secondary anisotropies, such as lensing and the Sunyaev-Zel’dovich (SZ) effect, depend on the late-time, nonlinear evolution of the universe. Their imprints 1 http://astro.uchicago.edu/spt/ http://bolo.berkeley.edu/apexsz/ 3 http://www.hep.upenn.edu/act/act.html 4 http://map.gsfc.nasa.gov/ 2 27 need not be Gaussian. In this work, we investigate the prospects of these experiments to extract the power spectrum and other statistical information for the components that contribute at these angular scales and frequencies. We are particularly interested in determining the primordial power spectrum from the primary CMB. Such a determination would enable one to determine the amplitude and slope of the fluctuations on small scales and, in combination with large scale observations, would allow one to place powerful constraints on models of structure formation. We model our investigations on the Atacama Cosmology Telescope (ACT), a 6 meter off-axis telescope to be placed in the Atacama desert in the mountains of northern Chile (Kosowsky, 2003). ACT’s bolometer array is anticipated to measure in three bands from 145–265 GHz, with planned beams of 1–2 arcminutes. (See Table 3.1.) ACT plans to survey about 100 square degrees of the sky with high signal to noise. In section 3.3, we simulate observations based on the specifications of ACT. The first problem in extracting the primary CMB information from observations is separating the components on the sky with distinct frequency dependences. Although many signals at these scales are not Gaussian, the power spectrum is still the most important of the statistics, and in section 3.4 we extract the power spectrum with a multifrequency analysis. The power spectrum analysis requires knowledge of the frequency dependence of the signals and of the scatter in that relation. Our knowledge is presently incomplete, and any error biases the estimated power spectra. Our discussion explores the impact of estimating some of this frequency information from the data itself. The second problem in extracting the primary CMB information is that one of the secondary anisotropies, the kinetic Sunyaev-Zel’dovich (kSZ) effect (and its analogs like the Ostriker-Vishniac effect and patchy reionization effect), has the same frequency dependence as the primary CMB. This complicates their separation, requiring an explicit template for the kSZ power spectrum. Again our knowledge is incomplete, 28 Band (GHz) 145 210 265 Beam FWHM (arcmin.) 1.7 1.1 0.93 Sensitivity per pixel (µK) 2 3.3 4.7 Table 3.1: Frequency channels, beam full-width at half maximum, and detector noise (thermodynamic temperature) of the mock survey, based on ACT design goals. and any error biases the primary CMB power spectrum. We address in detail the impact of secondary anisotropies on the estimation of the amplitude and slope of primordial fluctuations. Non-Gaussian information can help distinguish between the components which cannot be separated using frequency information. It is difficult to make non-Gaussian analysis fully exhaustive, so one must choose some simple statistics. Here we explore the one-point probability distribution function, in section 3.5. This can provide a consistency check on a kSZ template and can assist in signal separation. Finally, in section 3.6 we present our conclusions. 3.3 Simulations We simulate the sky as five astrophysical components (CMB, kSZ, SZ, and two types of point sources) plus noise. The CMB is a Gaussian random field. The density field (used for lensing of CMB), kSZ, and SZ are produced from hydrodynamical simulations (Zhang et al., 2002). Point sources are Poisson distributed from source count models. We add the appropriate detector noise and convolve with the telescope beam (assumed Gaussian). We briefly examine the simulations in general, before addressing each signal in turn. The simulation contains forty fields, each 1.19 ◦ × 1.19◦ , totaling ∼ 57 square degrees, somewhat smaller than the ACT survey. The results will therefore be ap- propriate for a smaller, more conservative mock survey. Throughout, we have used the following cosmological parameters in a flat, ΛCDM model: Ω c = 0.32, Ωb = 0.05, 29 ΩΛ = 0.63, σ8 = 1.0, with H0 = 70 km/s/Mpc. We express our fields in temperature units, or as a temperature perturbation Θ( n̂) = ∆T (n̂)/TCMB for line of sight n̂. Sample temperature perturbation maps of each signal are shown in Figure 3.1. We will refer to each image as we discuss the signals in greater detail. Signals which have a blackbody frequency dependence (CMB, kSZ) have the same temperature perturbation at all observed frequencies. For all other signals, the temperature perturbation will depend on the observation frequency. In Fourier space, we use the flat sky approximation, Θ(n̂) = 1 (2π)2 Z d2 l Θl exp(in̂ · l). (3.1) The symbol l is a two dimensional wavevector. The power spectrum of a signal or channel is defined by h|Θl |2 i = Cl . For each signal, figure 3.2 displays the power spectra and cross-correlations between the bands. Below we include a more detailed look at these spectra. ACT will have channels near 145, 217, and 265 GHz, which in simulations we take to have delta function frequency responses. These bands are typical of ground based experiments because they lie in windows between atmospheric lines. The band at 217 GHz also takes advantage of the thermal SZ “null,” described below in section 3.3.3. Figure 3.3 show a simulated patch of sky observed by ACT in these three bands. We discuss this image further below. In the next section, we discuss the hydrodynamical simulations which were the source of several of our signals. Then we study the signals and noise in more detail. Section 3.3.7 lists caveats to our simulations and section 3.3.8 summarizes the discussion of the simulations. 3.3.1 Hydrodynamical simulations Hydrodynamical simulations provide the density field for lensing, the SZ and kSZ effects. These simulations are described in more detail elsewhere (Pen, 1998; Seljak 30 CMB kSZ x 10 Rad. pt. src. x 10 SZ 4.0 2.0 0.0 IR pt. src. ∆Τ (100 µΚ) −2.0 Figure 3.1: These simulated signals contribute to one of forty 1.2 ◦ × 1.2◦ fields. All fields are in temperature units with a common gray scale, although kSZ and radio source signals have been multiplied by 10 for visibility. SZ is shown at −T RJ = 2y. The CMB, SZ, and kSZ are not convolved with any beam. Point sources are shown in the channel that they are brightest and for visibility are convolved with the appropriate beam: radio at 145 GHz (1.7 arcminute FWHM) and IR at 265 GHz (0.9 arcminutes). 31 et al., 2001; Zhang et al., 2002, 2004). The simulation included gas and dark matter in a 100 h−1 Mpc periodic box at 5123 resolution with a moving mesh. At several redshift slices, the contents are projected onto a plane parallel to one of the box faces. These projections are stacked with random horizontal and vertical offsets to form maps which look back through a cone of the sky. The benefit to this method is that only one box must be simulated to create a much larger set of maps. In our case, we have forty maps from one simulation box. The obvious downside to this technique is that a single structure in the simulation may be seen several times, at different orientations and stages of development, in different fields. Thus the family of rare objects in the maps may not be truly representative. The authors of (Zhang et al., 2002) provided us with convergence maps for CMB lensing, comptonization parameter maps for thermal SZ, and ∆T /T maps for kinetic SZ. These are discussed below for each effect. 3.3.2 Lensed CMB We employed CMBFAST (Seljak & Zaldarriaga, 1996) to generate the unlensed CMB power spectrum to l = 20 000. This power spectrum was used to create a Gaussian random field, which is the CMB map before any gravitational lensing occurs. Lensing bends light from one part of the sky to another, so that the light observed from lineof-sight n̂ in reality came from n̂ + δ n̂. The real space displacement vector δ n̂ may be calculated via the Fourier space relation, δ n̂(l) = 2ilκ(l)/l 2 . The convergence map for CMB lensing κ(n̂) is defined as a line-of-sight projection of the density perturbation δ ≡ ρ/ρ̄ − 1: 3H02 Ωm κ(n̂) = 2c2 Z χLS 0 r(χ)r(χLS − χ) δ(n) dχ. r(χLS ) a (3.2) The integral is over the comoving radial coordinate χ. Here χ LS is the coordinate of the CMB’s surface of last scattering at recombination, the line element radial function r(χ) is simply χ for a flat universe, Ωm is the matter density parameter, n describes position (n̂, χ), and a is the expansion factor of the universe. Software from Zaldarriaga & Seljak (1999) computes the lensed maps from unlensed maps and the convergence (See 32 also Hirata & Seljak (2003)). Hydrodynamical simulations provided the convergence calculated from the gas density, δg . assuming gas traces matter. The primary CMB has Gaussian fluctuations, but lensing introduces non-Gaussianities by correlating the CMB with the nonlinear matter field. A sample lensed CMB map is shown in Figure 3.1. The lensing has negligible visible effect on the maps. The CMB power spectrum is shown in Figure 3.2. Because we examine the temperature perturbation power spectrum, any signal with a blackbody frequency dependence will contribute equally in the three bands. This is the case for the CMB. Since the maps are only 1.2◦ on a side the power spectrum bins have a width ∆l ≈ 300, so the acoustic peaks are not resolved. The CMB is the dominant signal on scales larger than l ∼ 1000. The major features are the damping tail, where the power falls exponentially at 1000 < ∼ 4000, and the contribution from lensing, which ∼l < dominates the power at l > ∼ 4000. We see below that lensing is not a major effect in our analysis. 3.3.3 Sunyaev-Zel’dovich effect The SZ (or thermal SZ) effect is a change in the spectrum of the CMB due to scattering off of hot gas in the potential wells of large structures (see Birkinshaw (1999); Carlstrom et al. (2002) for comprehensive reviews). Scattered CMB photons are preferentially raised in energy. This distorts the energy spectrum: in the direction of the hot gas, the CMB has fewer low-energy photons and more high-energy photons. The number of scatterings for a photon passing through a cluster is low, so the CMB photons do not achieve thermal equilibrium with the electrons. Therefore, the SZ effect does not follow a blackbody spectrum. If the electrons are not relativistic then the frequency dependence is (equation 2.6) x/2 Θ (n̂) = −2y(n̂) 2 − , tanh(x/2) SZ 33 1e-09 Combined l(l+1) Cl / 2π 1e-10 SZ CMB 1e-11 CMB+kSZ 1e-12 kSZ 1e-13 1e-09 l(l+1) Cl / 2π 1e-10 145-145 145-217 145-265 217-217 217-265 265-265 IR Sources 1e-11 Radio Sources Noise 1e-12 1e-13 1000 10000 1000 l 10000 1000 l 10000 l Figure 3.2: Auto- and cross-power spectra from the ACT bands, in bins of width ∆l ≈ 300. Cl measures the power in temperature perturbation. Line type describes the frequency dependence of the signals. “145-145” denotes the power spectrum in the 145 GHz channel. “145-217” denotes the cross correlation between the 145 GHz and 217 GHz channels. Others correlations are labeled similarly. Each signal is depicted in a separate panel. ACT measures the sum of all signals and noise, shown in the panel marked “combined.” At each frequency in the combined signal, the CMB dominates at large scales. Small scales are dominated by experimental noise and infrared point sources. Effects with a Planck spectrum (CMB/kSZ) do not vary with frequency. The SZ 145 GHz-265 GHz cross correlation is negative, and is shown in absolute value. 34 where x = hν/kB TCMB and the comptonization parameter y is given by the integral of the electron pressure along the line of sight: y(n̂) = Z σT ne (n) kB Te (n) adχ. me c2 (3.3) Here σT is the Thomson scattering cross section, n e is the electron density, and Te is the electron temperature. This expression is derived in the single scattering, nonrelativistic limit. In the low frequency (Rayleigh-Jeans) regime, x 1 and Θ SZ = −2y. At roughly 217 GHz, the SZ effect is absent. This is known as the SZ null. At frequen- cies lower than the null, the spectrum is underpopulated relative to a blackbody, and the effect is seen as a temperature decrement. At higher frequencies, the spectrum is overpopulated, and the effect is a temperature increment. This very specific deviation from a thermal spectrum is the SZ signature, and is vital for extracting the effect. A sample SZ map, showing negative Rayleigh-Jeans temperature, is depicted in Figure 3.1. The SZ effect is most striking in clusters. The resulting signal is nonGaussian. The SZ effect has been studied extensively with simulations (e.g. (Zhang et al., 2002; da Silva et al., 2000; Springel et al., 2001; Seljak et al., 2001; White et al., 2002; Refregier & Teyssier, 2002)) and semi-analytic methods (Komatsu & Seljak, 2002; Cooray, 2000). There are some differences in the results between different simulations and halo model predictions (see Komatsu & Seljak (2002) for a discussion), but these are most prominent on very small scales, while in the range around l ∼ 3000 the agreement in the power spectrum is quite good. The SZ power spectrum from simulations used here (Zhang et al., 2002) is shown in Figure 3.2. The effect of observation frequency is apparent. The SZ effect is most prominent, and negative, in the 145 GHz channel, absent at 217 GHz, and is substantial at 265 GHz. At very large scales, the discrete nature of clusters is apparent, and the power spectrum is like shot noise, with ClSZ constant. At smaller scales, the power goes as C lSZ ∝ l−2 for scales of 2000 < ∼l < ∼ 10000. At yet higher l, the power falls faster. 35 3.3.4 Kinetic Sunyaev-Zel’dovich effect The kinetic SZ effect is due to the scattering of photons off gas with a radial peculiar velocity. It may be expressed as an integral over the density weighted gas velocity projected in the line of sight direction. Θ kSZ (n̂) = Z σT Xe (n)ne (n) v(n) · n̂ exp[τ (χ)]adχ. c (3.4) Here Xe is the ionization fraction of the gas and τ (χ) is the optical depth to comoving radius χ. To a good approximation, this fluctuation affects only the radiation temperature, so the spectrum of this effect is still blackbody. Photons scattering from electrons with radial velocity away from us show a temperature decrement and vice versa. In the simulations, reionization was instantaneous and spatially homogeneous. We discuss this assumption in section 3.3.7. Figure 3.1 shows an example kSZ map. Note the association between the kSZ effect and SZ effect. The same gas causes both the kSZ and the SZ effect, but the cross-spectrum vanishes because, on average, half of clusters have a peculiar velocity towards us, and half away from us. The equivalent of kSZ at higher redshifts is called the Ostriker-Vishniac effect. The second order term dominates over the linear effect, which is strongly suppressed. This effect is thus also expected to be non-Gaussian. Figure 3.2 shows the kSZ power spectrum. Like the CMB, kSZ is equally represented in all three channels. At higher l, we have roughly ClkSZ ∝ l−2 for 4000 < ∼l < ∼ 15000 (Zhang et al., 2004) before falling off at high l (see (da Silva et al., 2001; Springel et al., 2001) for other simulations of the kSZ). Modeling the kSZ remains quite uncertain. It is sensitive to the details of reionization which are not yet well understood and the power spectrum is probably uncertain to an order of magnitude (Santos et al., 2003). 3.3.5 Point sources Two populations of extragalactic point sources are the main source of astrophysical contamination for this type of survey (White & Majumdar, 2004). These are radio 36 galaxies, brightest in the lowest frequency channel, and dusty galaxies, brightest in the highest frequency channel. Dusty infrared point sources are the more serious contaminant. Point sources dominate the secondary anisotropies in the two higher frequency bands, and obey Poisson statistics (assuming no correlations). Thus the statistics are non-Gaussian. In the following, we discuss each family of sources, point source removal by masking, and the frequency coherence of the source population as a whole. Radio point sources Synchrotron emitting point sources feature prominently in the WMAP maps. They are less prominent here because of the frequency dependence, but still make a contribution. We model radio sources using a fit to the source count model of De Zotti et al. (2002) based on Toffolatti et al. (1998). For each flux level, we compute the number of sources per field according to a Poisson distribution, where the average density is given by the source counts. We place the point sources without spatial correlation, as this is simplest, and any correlations, particularly for faint sources, are diminished by projection. The flux of these radio sources is thought to decreases with increasing frequency, perhaps scaling as ν −0.5 to ν −1 De Zotti et al. (2002). For our simulations, we use a uniform random deviate for the exponent of the frequency dependence for each source, and compute the flux in each of the ACT bands, converting to temperature units. This frequency dependence remains uncertain in the frequency range appropriate for ACT. Measurements of WMAP sources at lower frequencies yield a flatter frequency dependence (Bennett et al., 2003b; Trushkin, 2003), but it is unclear if these sources will maintain this dependence up to ACT frequencies, or begin to drop off. In section 3.4.3, we discuss further the impact of frequency dependence on the power spectrum measurement. The power in uncorrelated point sources is given by Z Scut dN 2 dB −2 1 ps S dS Cl = 2 dS TCMB dT TCMB 0 (3.5) 37 where the derivative of the Planck spectrum dB/dT is used to convert into temperature units, and dN/dS is the differential source counts. After masking of bright sources, described below in section 3.3.5, the radio point source power spectrum is plotted in Figure 3.2. The radio power decreases with frequency. Since the sources are assumed to be spatially uncorrelated, the power spectrum is white: C lrad is constant. Radio point sources do not dominate in any of the bands at any scale, but proper treatment of them is important to get unbiased power spectrum results. Infrared point sources Dusty galaxies reprocess starlight to emit strongly in the infrared. To model these sources, we use source counts at 353 GHz. These source counts are based on a model of the optical emission of galaxies combined with a model for interstellar dust re-emission in the far infrared of absorbed optical light (Totani & Takeuchi, 2002). This model’s counts are lower than the SCUBA counts (Barger et al., 1999, 1998; Eales et al., 1999; Holland et al., 1999; Hughes et al., 1998; Smail et al., 1997; Borys et al., 2003) at the same frequency, but within a factor of ∼ 2 . We place them in the fields in the same manner as the radio sources, with no spatial correlations. The flux of these sources depends on frequency roughly as ν 3.5 , with some scatter. Thus infrared sources become more prominent at higher frequency. For our simulation, we chose frequency dependences of ν 3 –ν 4 . After masking out the bright sources (section 3.3.5) the infrared point source power spectrum is plotted in Figure 3.2. The IR power increase with frequency. Again the power spectrum is white, since the sources are assumed to be spatially uncorrelated. Infrared point sources are more contaminating than the radio sources in the higher frequency bands, and are the major contaminant in these bands until the instrument noise dominates at high l. At 145 GHz, infrared sources are dimmer than the radio sources, but there are many more of them, so the power in each is about the same. The power in infrared sources increases with frequency. 38 Source Masking To some extent, the brightest point sources can be removed individually by thresholding, frequency modeling, and masking. We explore this procedure. At what level may we reasonably (and conservatively) detect and remove point sources? What is the impact on a power spectrum analysis? Is information from surveys outside the ACT survey helpful in constraining the point source contribution? For simplicity in our simulations, we do not attempt to extract or mask point sources individually. Rather we impose a flux limit in our simulation, imagining that the most offending point sources have been already excluded. Our strategy for identifying point sources focuses on the 145 and 265 ACT bands, where radio and infrared sources are most prominent, respectively. We spatially filter the maps to emphasize point sources and de-emphasize the CMB and detector noise. As a side effect, this also enhances SZ clusters in the maps. The filter we use is the reciprocal of the channel’s power spectrum (Figure 3.2), but any similar broad filter would work. We set a flux thresholds, and identify sources. A +5σ threshold at 145 GHz identified radio point sources down to 4 mJy. At 145 GHz SZ clusters are negative, so there is no chance to confuse SZ sources as point sources. At 265 GHz, a +5σ threshold identifies infrared point sources to 5.5 mJy. This threshold is sufficiently high that only the very brightest SZ clusters could be misidentified. Source shape and frequency information should eliminate any confusion. Thus we adopt these flux cuts in our simulations. As point sources are masked out, we reduce the survey area. The masked area can be calculated from the source counts. We assume the largest beam (1.7 arcmin) as the rough scale of a mask. The fraction of the survey remaining unmasked after N sources are excised is roughly (1 − (mask)/4πf sky )N , assuming masks can overlap. The radio and point source cuts described above will remove about 2 percent of the survey area. Might we make a more aggressive flux cut to reduce the power in power sources, using multifrequency information or a different point source detection scheme? Yes, 39 but to do so is not always useful. In a power spectrum analysis, a measured signal has variance like var(Clsignal ) ∝ (Clsignal + Clps + Clnoise )2 . fsky (3.6) The purpose of masking sources is to reduce C lps , thus reducing the variance of our measured signal. Note, however, that masking also reduces f sky , which can serve to increase the variance if we mask too many sources. We shall see that for infrared sources, even simple methods for identifying point sources identifies so many sources that masking them all would be counterproductive. This can be understood visually in Figure 3.3. A simple multifrequency combination can readily identify infrared point sources below the 5 mJy flux cut at 265 GHz applied that image. For example, compare the 265 and 145 GHz channels and imagine masking every visible point source. The bulk of the survey would be consumed by the mask. As a more concrete example, lowering the flux cut on infrared sources from 5.5 mJy to 1.6 mJy at 265 GHz removes ∼ 35 percent of the survey area. Such a flux cut reduces the power in such sources by ∼ 20 percent. For a signal whose variance is dominated by faint point sources, this aggressive flux cut improves the variance only slightly. For a signal comparable in amplitude to the point source noise, the variance actually worsens. At these flux levels the power in point sources is not dominated by a few bright sources, but by a mass of background sources. Thus it may be sufficient, or even advantageous, to use a more conservative point source removal. Whether or not to aggressively pursue point source removal depends on the particular application. Given a model for the signal and the source counts, one could in principle minimize the variance using expressions 3.5 and 3.6. Since we are examining signals over a wide range in l, and the relative power in point sources and other signals varies greatly, we will simply use the more conservative cut for infrared sources. The differential counts for radio sources at the level of our flux cut are shallower, 40 so masking is useful. It would also be helpful to compare the ACT survey at 145 GHz to radio surveys covering the same area, to assist in identification and masking of radio sources. This could reduce the variance caused by radio sources, if such information were available. However, we will proceed as if the fields from ACT are the only measurements of those sources available to us. Information from outside surveys to identify and mask offending radio point sources is helpful, but not necessary for a power spectrum analysis with ACT. Frequency Coherence For some purposes, it is useful to treat all sources as a single population. In this case, the details of the multiple source populations, frequency dependence, and scatter, may all be subsumed under a set of correlation coefficients which describe the frequency coherence across bands. If the correlation coefficient for point sources between two bands is unity, then one band is a scaled version of the other. Note that since significant power can come from just below the flux cut, the correlation coefficients depend on the cut. In our point source model, the correlation coefficient between the 217 and 265 GHz bands is 0.99, between 145 and 217 GHz it is 0.83, and between 145 and 265 GHz it is 0.76. IR sources dominate the point source contribution in the higher two bands, explaining their close correlation. The 145 GHz band mixes in a significant fraction of radio sources, which decreases the frequency coherence. This frequency incoherence has serious implications for the power spectrum, as we note below. 3.3.6 Detector noise We add Gaussian distributed detector noise to these simulated fields. In reality, the cosmic signals will be convolved with the beam of the telescope, and the detector will add white noise. Here we find it more convenient to take the mathematically equivalent approach of de-convolving the beams from the three channels. Thus the cosmic signals are untouched and the noise is correlated, with noise power increasing at small 41 scales. Assuming the noise and Gaussian beams in Table 3.1, the power spectrum of noise is given by: noise Clα 2 = (σα bα ) exp (lb )2 √ α 8 log 2 (3.7) Where α stands for the channel, σα is the RMS noise per beam-sized pixel, and the beam full width at half maximum bα must be converted to radians. Figure 3.2 shows the noise power spectrum in the three channels. At large scales the noise is similar in the three channels. At smaller scales, the noise dominates in the 145 GHz channel first, then at 217 GHz, then 265 GHz. 3.3.7 Caveats and other signals not included Here we list several caveats to our simulations. These issues will need to be addressed in further studies before the results of ACT may be applied to precision cosmology. The hydrodynamical simulations we used have σ 8 = 1 normalization. Except for the excess in the CBI power spectrum at high-l (Mason et al., 2003), this is at the high end of currently favored values. The SZ power varies strongly with C lSZ ∝ σ87 (Komatsu & Seljak, 2002), so this could be a substantial effect. At high electron temperatures, the formula for SZ frequency dependence needs relativistic corrections, which are a function of electron temperature. We have ignored these corrections in this analysis. These corrections displace the SZ null, but this is likely minor compared to the bandwidth of the actual experiment. The frequency dependence is a constant function of scale in the non-relativistic model, but relativistic corrections change the frequency dependence (because hotter clusters are more massive, larger, and tend to appear at lower redshift). We have computed this effect using halo model power spectra from Komatsu (private communication), and find the effective frequency dependence for SZ varies 15 percent when comparing 145 and 265 GHz for 100 < l < 10000. Another (cruder) model, where we have made SZ maps from combining cluster profiles and mock halo catalogs, puts this variation at about 7 percent. Despite this variation, the signal remains quite coherent across bands, in that 42 the cross correlation coefficient on all scales (in Fourier space) is very close to unity in absolute value. If necessary, we could fold this correction to the SZ frequency dependence into our power spectrum analysis. Another approach would be to identify and correct the pixels in the hottest clusters. Either approach will complicate and degrade our power spectrum recovery. The hydrodynamical simulations we incorporate into our mock survey assume a spatially uniform transition at reionization. In reality, the ionized fraction will depend on where ionizing photons are produced, their ability to propagate, and their efficiency to ionize. Since these quantities probably vary from place to place, reionization is not likely to be uniform. Modeling of this so-called patchy reionization is difficult, and the impact remains somewhat uncertain. Thus we have not included it in this analysis. However, its impact may be significant: in recent modeling by Santos et al. (2003), patchy reionization typically increased the power in the kSZ effect by an order of magnitude (see also Miralda-Escud é et al. (2000)). Overall, the modeling of point sources is still quite uncertain. Source number counts, clustering, frequency dependence, and correlation with other signals all require improvements as better information becomes available. For example, as a simplification, we ignored the likely correlation between sources and SZ clusters, which will have several effects. Correlations may lead to the unintended elimination of some of the SZ clusters during source masking, reducing the SZ signal. At low frequencies, radio source flux may fill in some of the SZ decrement, leading to a decrease in the SZ power (up to 20–30% at 30 GHz (Holder, 2002; Zhou & Wu, 2004)), although ACT’s multiple bands make it less susceptible to this error. If we take the above 30 GHz result and apply it to the worst-case scenario: all radio sources are flat spectrum and are exactly correlated with SZ, the correction to the SZ power spectrum in the ACT bands is at worst 10%, comparable to the current modeling uncertainty of the SZ signal (Komatsu & Seljak, 2002). We need to uncover the details of this correlation in order to use the SZ effect for precision cosmology. 43 Additionally, the infrared point sources may be correlated with CMB lensing, with this cross-correlation possibly detectable by Planck (Song et al., 2003) using the 545 GHz band. We have not included this effect and do not consider it. The peak scale of the effect at l ∼ 150, frequency scaling of point sources, and f sky of the experiments favor Planck, not ACT, to detect this effect. We have not included any galactic contaminants in these simulations. Because the survey covers a small fraction of the sky, we have tacitly assumed that a clean enough window could be found to look out of the galaxy with minimum contamination or that a slowly varying galactic component can be projected out of the analysis. Unfortunately, current dust templates ((Schlegel et al., 1998; Finkbeiner et al., 1999; Bennett et al., 2003b)) have poorer spatial resolution than ACT, making the second assumption difficult to test. Galactic dust has its own frequency signature, and a more detailed analysis should make an allowance for a dust component, which could bias a multifrequency power spectrum reconstruction if not accounted for. Finally, the model of the instrument we include is simplistic. The actual survey will have a nontrivial geometry, non-uniform noise with a 1/f component, and nonGaussian beams. 3.3.8 Summary of simulations We have constructed a 57 square degree mock CMB survey according to the specifications of ACT. It includes five astrophysical components. The CMB is a Gaussian random field, lensed by a convergence map from hydrodynamical simulations. The same simulations produced the gas pressure (for SZ) and momentum (for kSZ). Both SZ and kSZ have well defined frequency dependences. SZ is negative in the 145 GHz channel, zero in the 217 GHz channel, and positive in the 265 GHz channel. Kinetic SZ scales with frequency exactly as the primary CMB. The simulations include two types of point sources. Radio sources dim with frequency and have their greatest contribution at 145 GHz and little contribution at higher frequencies. Infrared point sources 44 145 GHz 217 GHz 4.0 2.0 265 GHz 0.0 ∆Τ (100 µΚ) −2.0 Figure 3.3: A simulated 1.2◦ × 1.2◦ field in the bands of the ACT telescope. Pixel value indicates the difference in temperature from T CMB . The mean flux increases with frequency because of the increasing background of galaxies shining in the far infrared. Small scale temperature decrements at 145 GHz indicate the SZ effect. 45 brighten with frequency, and are prominent in the 217 GHz and 265 GHz bands. We have conservatively removed the brightest radio and infrared sources. A more aggressive flux cut on radio sources can help a power spectrum analysis, if more information about the sources is available from other surveys, but we assume it is not. A more aggressive flux cut on infrared sources begins to degrade the survey area, and depending on the application may not improve the variance in a power spectrum analysis. We must include infrared sources in this analysis because they are too numerous to excise. We include instrument noise and beam from Table 3.1. The instrument has highest resolution at 265 GHz and lowest at 145 GHz. The assembled signals, convolved with the telescope beam and added to the noise, are depicted in the ACT channels in Figure 3.3. The striking features are the primary fluctuations on large scales, bright SZ clusters (as a decrement) at 145 GHz, and the substantial brightening of infrared point sources at higher frequencies. The power spectra of these channels are shown in Figure 3.2. The power spectrum in each of the three channels is dominated by the CMB on the larger scales and instrumental noise on the smaller scales. At the intermediate scales of this survey, SZ is the most significant component in the 145 GHz band. In the bands at 217 GHz and 265 GHz, infrared point sources are the major contribution to the point sources at intermediate scales. At this point we may draw several conclusions about these simulations. The CMB is sub-dominant on a wide range of scales accessible to this survey. The dominant cosmological features are secondary anisotropies, but these suffer from serious pollution from extragalactic point sources, particularly infrared sources in the high frequency bands. These astrophysical components exceed the noise over a wide range of interest. In section 3.4 we proceed with a power spectrum analysis and return to non-Gaussian features in 3.5. 46 3.4 Power spectrum analysis From our simulated maps, we seek to measure the power spectra of our various cosmic signals. We use a minimum-variance quadratic estimator, which can also be viewed as an iterative procedure to find maximum-likelihood solution (Bond et al., 1998). The power spectrum estimates are unbiased. For components that cannot be separated using the frequency information (formally in this case the Fisher matrix becomes singular and the errors infinite) it is better to treat them as a single component and then try to separate them using additional considerations such as a prior knowledge of the power spectrum shape for the signals. For some applications, power-spectrum estimators that involve the intermediate Wiener filtered maps are minimum variance (Seljak, 1998). However for this application, they are not necessarily minimum variance. Our approach here is a generalization of those treatments and should improve upon those estimators. As an added bonus, the effects of imperfect cleaning of one component on the mean and errors of other components are automatically included, for the estimator is unbiased and minimum variance (considering the imperfect cleaning, see Tegmark et al. (2000)). For non-Gaussian signals, the estimator is unbiased, but not minimum variance. In this section, we discuss this estimator and its relation to the Wiener filter. We discuss the impact of scatter in the frequency dependence of signals. We measure the power spectra in several cases, making different assumptions about our knowledge of the point source frequency dependence. Finally, we discuss the importance of cleaning out kSZ in determining the amplitude and slope of primary fluctuations. 3.4.1 Multifrequency estimator In this section we describe the power spectrum estimator we employ. We will reference a particular Fourier mode in the flat sky approximation by its wavevector l. We use l to denote the magnitude of this wavevector. We organize our data into a vector e = {e lα }. The subscript α refers to which band the data comes from. We will reserve Greek in- dices to refer to the channels of the detector, αβγ. . . . We say that the data is given by 47 the sum of the response of the instrument to the signals and the noise. Dropping the mode index, we write the data in channel α as e α , the signal as sα , and the noise as nα . For a particular multipole, we characterize the data, signals, and noise by their covariances heα e∗β i = Cαβ , hsα s∗β i = Sαβ and hnα n∗β i = Nαβ . Note that it is not necessary to assume the noise between channels is uncorrelated. As always, understanding the noise correlation matrix is vital, because an incorrect assumption about the noise immediately biases the result: unaccounted for noise correlations between channels will masquerade as an additional frequency-dependent component. We assume that the signal and noise are uncorrelated, and that the signal depends on a set of parameters Z = Zp . Thus, restoring the multipole dependence, the covariance matrix is l C = hee† i = Cαβ δll0 (3.8) l l l . (Z) + Nαβ = Sαβ Cαβ (3.9) where The covariance matrix C has both frequency channel and and multipole indices, and we have assumed that different multipoles are independent, hence the δ ll0 . In a Gaussian model, the likelihood function for the data is 1 L(e) = (2π)−N/2 det(C)−1/2 exp(− e† C−1 e). 2 (3.10) A second order expansion in Zp gives ln L(Z + δZ) = ln L(Z) + X ∂ ln L(Z) p + ∂Zp δZp 1 X ∂ 2 ln L(Z) δZp δZp0 , 2 0 ∂Zp ∂Zp0 (3.11) pp where ∂ ln L(Z) ∂Zp 2 ∂ ln L(Z) − ∂Zp ∂Zp0 −2 ∂C −1 ∂C C e − C−1 ) ∂Zp ∂Zp ∂C −1 ∂C −1 C = e† C−1 C e ∂Zp ∂Zp0 ∂C −1 ∂C 1 C ). − tr (C−1 2 ∂Zp ∂Zp0 = tr (e† C−1 (3.12) 48 The ensemble average of the second expression is the Fisher matrix 1 ∂C −1 ∂C Fpp0 = tr (C−1 C ). 2 ∂Zp ∂Zp0 (3.13) Note that for a given value of l, the trace accounts for the sum over (2l+1)f sky modes. At the maximum likelihood value the first derivative of the likelihood function vanishes, so we use the Newton-Raphson method to find the zero of the derivative. This leads to the minimum variance quadratic estimator for the parameters Z (Hamilton, 1997; Tegmark, 1997; Seljak, 1998). Ẑp = 1 X −1 † −1 ∂C −1 F 0 [e C C e − b p0 ] 2 0 pp ∂Zp0 p bp = tr [C−1 ∂C −1 C N]. ∂Zp (3.14) We use Ẑp to denote the estimate of Zp . In this formalism, the covariance matrix of the estimated parameters, written C pp0 , is given by the inverse Fisher matrix: (3.15) −1 Cpp0 = hẐp Ẑp0 i = Fpp 0 This technique requires knowledge of the covariance C. We may estimate this by taking a smoothed version of the data. If the covariance is incorrect, the estimator is still unbiased, but the error bars will be incorrect. The parameters we desire to estimate are the power spectra of the signals, defined in a way which does not depend on the frequency channels of our instrument. We denote our signal with s = slm . Here m indicates which signal we are referring to. l We will also use the (cross) power spectrum of signals, denoted S mn = hslm sln iδll0 . We 0 introduce the frequency response R αβmn to relate the spectra of the signals in the l , however we have chosen to define them. channels to the spectra of the signals S mn l Note that Smn is symmetric in m and n, so only one ordering of each pair is required. Then the cross power spectrum from the combination of channel α and channel β will be given as: l Sαβ = X mn l Rαβmn Smn . (3.16) 49 We discuss the frequency response in more detail below. To estimate the power specl . Note that in this case the index p enumerates trum of the signals we set Zp = Smn each combination of lmn. In this case the appropriate derivative of the covariance is given by 0 l ∂Cαβ ∂Zp 0 = l ∂Cαβ l ∂Smn X = m0 n0 Rαβm0 n0 δll0 (3.17) × [δmm0 δnn0 + δmn0 δnm0 − δmm0 δnn0 δmn0 δnm0 ] The final expression in the brackets is shorthand for: 0 0 0 1 if (m = m and n = n ) or (m = n and n = [. . .] = m0 ), 0 otherwise. (3.18) If we know the shape of the cross spectra and wish to estimate the amplitude, we can write the signal power in the bands as l Sαβ = X mn l Rαβmn Amn Smn (3.19) l where Amn encodes the information about the amplitude, and S mn is the shape. If we want to estimate this amplitude, then Z p = Amn . Thus the derivative of the covariance is given by l ∂Cαβ ∂Zp = l ∂Cαβ ∂Amn = X m0 n0 l Rαβm0 n0 Smn (3.20) × [δmm0 δnn0 + δmn0 δnm0 − δmm0 δnn0 δmn0 δnm0 ] If we know the shape of some signals, but not others, we may use the two expressions (3.17) and (3.20) in combination. In these expressions we have assumed that the power spectra of the components are unknown, but that their frequency dependence is known (and hidden in Rαβmn coefficients). If the parameters which determine the frequency scaling are not known then one can generalize this further, add these parameters to the estimate, and maximize the likelihood over them as well, provided the Fisher matrix does not become singular. The procedure of taking derivatives of the covariance matrix is similar, so we do not give such expressions here. 50 3.4.2 Relation to Wiener filter If all signals have a spatially uniform frequency dependence, then s lα = P m Rαm slm for all l (this defines Rαm ). The correlation coefficient of modes between frequency channels will be unity (neglecting noise). In this case we can factor the frequency response into a simpler form: Rαβmn = Rαm Rβn (3.21) We may put the components into a vector R m = Rαm . This allows us to compactly express the Wiener-filtered estimates of the signal maps, which minimize the variance in reconstruction errors. The same expression holds for all modes, so we drop the mode index l. The Wiener estimate is ŝm = ym = X Smm0 ym0 m0 −1 RT m C e. (3.22) (3.23) The ŝm are the Wiener filtered estimates of all signal maps. From these estimated maps of the signals, we can compute the power spectrum, taking into account the contamination from the other signals to the Wiener filtered maps. In the case of spatially uniform frequency dependence, this procedure is identical to the multifrequency estimate in the previous section. The Wiener estimate requires the power spectrum S im , or we could use a power spectrum estimate and iterate. For signals with non-uniform frequency dependence, it may not be possible to write a vector Rm which satisfies (3.21). A spatially-averaged frequency dependence may still be used to produce a Wiener-type estimate, but it is not guaranteed to be an unbiased or a minimum variance reproduction of the signal. In addition, the power spectrum obtained from these maps can be substantially biased. This is because a power spectrum estimate which uses a Wiener filtered map as an intermediate step does not account for the scatter in the frequency dependence. If this scatter is substantial, the Wiener filter may be inappropriate. (See Cooray et al. (2000) for making minimum variance maps of signals with good frequency coherence in the presence of incoherent 51 foregrounds. To make a map of an incoherent signal, like a point source template, this approach may not be minimum variance. See Tegmark et al. (2000) for a generic model of frequency incoherence.) 3.4.3 Frequency dependence and scatter To apply the multifrequency filter, one needs the noise power spectra and signal frequency dependences. We assume the noise power spectrum in the three channels is well known. The frequency scalings of CMB and kSZ signals are constant in temperature units, so the relevant components of R αβmn are unity. For SZ, the frequency dependence is from equation 2.6. The frequency dependence of point sources are poorly known. Significant modeling or outside observations are needed to fill in the remaining components of Rαβmn . These cannot be completely determined from ACT alone in an unbiased way. If they are added to the parameter vector Z p without any additional assumptions, the Fisher matrix becomes singular. This means these parameters are completely degenerate with the power spectra: there exists a family of power spectra and frequency dependences which reproduces the covariance, and it is impossible to distinguish between individual members of this family. It is possible, however, that for given models of the point source emission, the parameters of the model may be included in the estimate without forcing the Fisher matrix to be singular. This models could also absorb constraints from other measurements, such as the SCUBA counts. Additionally, we can use the high l part of the spectrum, which is dominated by point source emission, to help constraint frequency dependence, at the expense of biasing our estimator. We use this in the following section. In our simulations, the frequency scaling differs substantially from point source to point source. The correlation coefficient between modes at different frequencies is less than unity, so the cross spectrum is less than predicted from the auto-powers. It may be senseless to produce a Wiener filtered signal estimate based on a spatially averaged 52 frequency dependence. We do not use a spatially averaged frequency dependence in our power spectrum estimations, for we find it produces greater than 2σ biases in the power measurement. The conclusion is that it is crucial to address the scatter in the frequency dependence of point sources. ACT and similar experiments will identify a population of bright sources, which may be used to study their frequency scaling. One should apply this information to a power spectrum analysis with some caution, because the sources which contaminate the power spectrum analysis are dimmer than the population which is convenient to study. This may mean that contaminating sources are a different population, more distant, younger and less evolved, or some combination. Dust emission is complicated, and the emission spectra of the point sources likely depends on the relative abundances of the dust species present. If the populations are different, it will be difficult to build an accurate model of the emission from these sources. In addition, bright sources may result from the confusion of two dimmer sources. Measurements of the scatter of source frequency dependence is then difficult to interpret. High-resolution interferometric measurements would help in studying the source population. In particular, ALMA5 has frequency coverage in the ACT bands. It should have sufficient sensitivity to pick out sources which are confusing ACT, and sufficient angular resolution to resolve them. ALMA would be useful as a tool to explore the frequency dependence and scatter of sources in the larger, lower resolution ACT survey. 3.4.4 Measured signal power spectrum We measure the power spectrum in several cases assuming differing prior knowledge. Our list of cases is not exhaustive, but illustrates techniques one might explore. We always assume the frequency dependence of CMB, kSZ, and SZ. We enumerate our assumptions, and the parameters to estimate in each case. The first two cases assume we have perfect knowledge of the point source frequency 5 http://www.alma.nrao.edu/index.html 53 scaling and scatter. Case 1. We assume we know the frequency dependence of the point sources. The parameters to estimate are the power spectra, in bins, of CMB+kSZ, SZ, IR and radio point sources. Case 2. We assume we know the power spectrum shape and frequency dependence of the point sources. Since we know their shape, we only need to estimate their amplitude. The parameters to estimate are the binned power spectra of CMB+kSZ and SZ, and 2 amplitudes, for the IR and radio power spectra. The next two cases assume nothing about the point source frequency dependence. Case 3. We assume we know the shape of the point sources, and that the CMB+kSZ makes no contribution at l > 15000. The final condition ensures a non-singular Fisher matrix. We treat the point sources as a single population. Parameters to estimate are the binned power spectra of CMB+kSZ and SZ, and 6 amplitudes, corresponding to the point source amplitude in the cross-correlations between the bands. This is equivalent to estimating the frequency dependence of the sources. Case 4. Like case 3, but without assuming the shape of the point source power spectrum. We treat point sources as a single population, assuming that the scatter in the frequency dependence is constant as a function of scale, but we do not assume what it is. Parameters to estimate are the binned power spectra of CMB+kSZ and SZ, the 3 binned auto-power spectra of point sources in the three bands, and 3 correlation coefficients, which describe the scatter in the frequency relation. The final two cases separate kSZ from the CMB using a template for the kSZ power spectrum. Case 5. Like case 2, including point source frequency and power spectrum shape information, but with two more assumptions. First, we assume the shape of kSZ power spectrum. Second, we assume the CMB makes negligible contributions to 54 the power for l > 6000. The second additional assumption prevents a singular Fisher matrix. The parameters to estimate are the binned power spectra of CMB and SZ, and 3 amplitudes, for the kSZ, IR, and radio power spectra. Case 6. Like case 3, including only point source power spectrum shape information, but with the additional assumptions of case 5. Parameters to estimate are the binned power spectra of CMB+kSZ and SZ, 6 point source amplitudes, and the amplitude of the kSZ spectrum. If the parameters we choose to estimate can faithfully represent the real world, our estimate is unbiased. However, in any of our cases, an incorrect assumption can bias the estimate. Note that case 4 requires the least outside information. Where we use a frequency dependence, we use 265 GHz as our reference frequency for SZ and infrared sources, and 145 GHz for radio sources. This choice is a convention, and has no bearing on the quality of the estimate. We estimate the covariance on the measured power spectra in two ways. First we use the inverse Fisher matrix, so the error on the estimate of power bin p is (F −1 )pp . 1/2 Second, we computed the covariance from the estimates on our 40 realizations. In Figure 3.4, we show the errors on these power estimates of our simulations as a function of l. In the regime where instrumental noise dominates (at sufficiently high l), the inverse Fisher matrix provides a good estimate of the errors. Elsewhere, we find the inverse Fisher matrix underestimates the error. We tested the error estimate with Gaussian random fields (where the Fisher covariance gives the true covariance) substituting for the non-Gaussian effects, finding the estimated covariance much more reliably reproduces the Fisher covariance, compared to when the real signals are used. If the Fisher covariance is a good estimate of the true covariance, we can estimate analytically the probability that the errors from the realizations are so much higher by chance. Consider a single bin, and compute the ratio of the power difference from the mean and the Fisher error: z ≡ (C − C̄)/σFisher . Then z has zero mean and unit variance. Approximating the distribution of z as Gaussian, we have the sum over 55 realizations χ2 = P z 2 is χ2 -distributed with 39 degrees of freedom for 40 realizations (one degree of freedom is used up computing the bin’s mean). From Figure 3.4, which shows both RMS ∆C = (C − C̄) (with points) and σFisher (with lines), we can see that the ratio of computed error to Fisher error routinely exceeds ∼ 5. This ratio is (χ 239 /39)1/2 , so we have χ239 /39 ≥ 25, which has vanishing probability. Here we only considered a single bin, so the bound on the chance that the deviation is by random fluctuation is even smaller. We conclude that the Fisher estimate is a poor estimate for the true covariance. The increase in error over the Fisher matrix prediction is probably due to the assumption in equation 3.8 that the covariance is diagonal in l, implying mode crosscoupling, and due to the non-Gaussianities inherent in the SZ and point source signals. Experiments of this type will likely need to rely on Monte Carlo simulations for their errors. We now discuss each case in turn. The assumption in case 1 and 2 of point source frequency knowledge is probably dangerous if we do not have such knowledge. A poor model will bias the power spectrum in a way which may be hard to predict. On the other hand, if we do have a good understanding of the source population (from ALMA perhaps), cases 1 and 2 represent ideal situations. For case 1, the measured power spectrum is shown in Figure 3.5. In this plot we use the errors evaluated over our independent realizations. Except for radio sources, we recover the power spectra well. This provides information about the shape of the power spectrum, which we have not assumed. Our result in case 1 is unbiased. We compared the inverse Fisher matrix and the covariance matrix from the realizations. In case 1, the inverse Fisher matrix is diagonal in l. However, we find that the actual covariance matrix is not. The covariance between the CMB+kSZ power bins is nearly diagonal in this case, but to l of a few thousand, the SZ power bins are strongly correlated. For case 2, assuming the shape of the point source power spectrum improved the 56 l(l+1) ∆Cl / 2π 1e-11 1e-12 1e-12 1e-13 1e-13 1e-11 l(l+1) ∆Cl / 2π 1e-11 CMB+kSZ IR 1e-11 Monte Carlo Case 1 Case 2 Case 3 Case 4 1e-12 SZ Rad. Fisher Case 1 Case 2 Case 3 Case 4 1e-12 1e-13 1e-13 1000 10000 l 1000 10000 l Figure 3.4: The errors on the estimated power spectrum. Lines indicate errors derived −1 )1/2 . Symbols indicate the measured standard deviation from the Fisher matrix: (Fpp from 40 realizations. The errors in case 5 are similar to case 2, and the errors in case 6 are similar to case 3. For cases which assume the shape of the point source power spectrum, plotting the errors as a function of l is not possible, so the bottom panels show only case 1. 57 1e-09 1e-11 CMB+kSZ SZ l(l+1) Cl / 2π 1e-10 1e-11 1e-12 1e-13 l(l+1) Cl / 2π 1e-09 1e-12 1e-10 IR Rad. 1e-11 1e-10 1e-12 1e-11 1e-13 1e-12 1000 10000 l 1e-14 1000 10000 l Figure 3.5: Estimates of the power spectra in case 1, where we assume frequency information, but no information about the shape of any power spectrum. The solid line is the true value of the power spectrum. The × symbols show the estimated value of the power spectrum. The dashed lines show 1 standard deviation above and below the estimate, evaluated over our 40 realizations. This estimate is unbiased. The width of the bins is ∆l ≈ 300. 58 1e-09 1e-11 CMB+kSZ SZ l(l+1) Cl / 2π 1e-10 1e-11 1e-12 1e-13 1e-14 1000 10000 1e-12 1000 l 10000 l Figure 3.6: Estimates of the power spectra in case 3, where we assume the shape of the point source power spectrum, but assume nothing about frequency dependence. errors on the CMB+kSZ power spectrum by a factor of ∼ 4 relative to case 1. The errors on SZ are marginally improved. The power spectrum bins are somewhat correlated. In case 2, the error in the power spectrum amplitude in IR and radio point sources is ∆Cl /Cl = 5.0 × 10−3 and ∆Cl /Cl = 0.20 respectively. Recall that this case combines all the information on the point source spectra into two amplitudes. The inverse Fisher matrix underestimates this error by a factor of 2. We do not plot case 2 because it is unbiased and similar to case 1, with smaller error bars, although the bins are more correlated. Case 3 does not depend on prior frequency information, but has its own problems. Case 3 depends on an incorrect assumption that only point sources contribute at high l, while kSZ continues to have some small contribution at these l. Thus the power spectrum we deduce for point sources is slightly high, and this biases our measurement somewhat. This is relatively independent of the l cutoff. We plot the measured power spectrum in Figure 3.6, where the bias is apparent. The power spectrum bins are also somewhat correlated. However, we needed no model for the point source frequency dependence. Our model is biased by as much as 1σ at high l, but the bias is probably under better control than it would be for a model of frequency dependence, because at least we know in which direction our estimate is biased. Case 4 assumes less about the point sources than case 3, because it does not assume a shape for the power spectrum. We plot the power spectrum of CMB+kSZ and SZ in 59 1e-09 1e-10 CMB+kSZ SZ l(l+1) Cl / 2π 1e-10 1e-11 1e-11 1e-12 1e-13 1000 10000 1e-12 1000 l 10000 l Figure 3.7: Estimates for case 4. Assuming no information about the power spectrum shape or frequency dependence of the point sources, we lose the ability to reliably detect kSZ. The estimate for SZ is biased, but the statistical errors are large. The bins are highly correlated. Figure 3.7. This case does not successfully recover kSZ. There is a systematic bias in the SZ, but statistical errors dominate. The bins are strongly correlated. Systematic biases dominate statistical errors in the recovery of the point source power spectrum. The relative bias is greatest at 145 GHz, where the point source contribution is the least, and decreases at higher frequencies. Note that in case 4 especially, the Fisher matrix characterizes the errors poorly. In cases 5 and 6, the kSZ is separated from the CMB. This is necessary to allow parameter estimation based on the CMB. Because the CMB and kSZ are degenerate in frequency, this requires a model for the kSZ power spectrum. The estimator in this case uses the high l portion of the kSZ spectrum, where CMB is not significant, to calibrate the model of kSZ. kSZ is then subtracted at low l. We use a 1-parameter model for the kSZ power spectrum where the shape is known but the amplitude is not. The approach generalizes to a multiparameter model. We do not plot case 5. For case 6, we plot the recovered power spectrum (Figure 3.8). We do not detect a bias on the CMB in either case, although the SZ power spectrum is biased in case 6 as it was in case 3. The CMB power bins for l > 3000 are correlated by the kSZ template subtraction. SZ bins are also strongly correlated. We now summarize our power spectrum analysis results. Given a good under- 60 1e-09 1e-11 CMB SZ l(l+1) Cl / 2π 1e-10 1e-11 1e-12 1e-13 1e-14 1000 10000 l 1e-12 1000 10000 l Figure 3.8: Estimates of the power spectra in case 6, where we assume the shape of the kSZ and point source power spectra, but assume nothing about the point source frequency dependence. standing of the frequency dependence of point sources, we may successfully estimate the power spectrum of each component, provided we treat the CMB and kSZ as a single component. With no prior knowledge of the frequency dependence of point sources, we may make a somewhat biased estimate of the power spectra if we know the shape of the point source power spectrum. Knowing neither the frequency dependence nor the shape of sources, we detect SZ, but cannot positively identify kSZ. For us to have confidence in our estimates, we must trust our knowledge of the point source frequency dependence. This prospect is unclear. If we do not have confidence in our frequency knowledge, we may still make an estimate of the power spectrum, but it will be somewhat biased. We may speculate on the impact of the addition of a fourth band to this type of experiment. For deriving cluster parameters from the relativistic SZ effect, (Knox et al., 2004) indicates that a 30 GHz band is important for breaking the degeneracy between the comptonization parameter y and the cluster gas temperature. However, for a power spectrum measurement, it may be more important to constrain the frequency dependence of IR sources, which would be better served by a higher frequency band (∼ 350 GHz). In our more optimistic cases, if we want to estimate CMB separate from kSZ, we require a parameterized model for kSZ. The power spectrum template will likely come from simulation. However, any parameterized fit to the power spectrum of the signal components may introduce a bias. 61 3.4.5 Measuring primordial amplitude and slope In our model, the kSZ power is already comparable to the size of the CMB error bars at l ∼ 2000, so if we do not know the kSZ power spectrum at all then we should not use the CMB power spectrum beyond there. With a model for the kSZ power spectrum we may use bins with higher l, as in our cases 5 and 6. To illustrate this difference, we compute the covariance matrix for primordial amplitude A and slope n about the fiducial values A = 1 (COBE normalized) and n = 1. To convert from errors on binned power to errors on parameters we use C−1 qq 0 = X dCl ll0 dq C−1 ll0 dCl dq 0 (3.24) where Cqq0 is the covariance matrix for parameters q, C ll0 is the covariance matrix for the binned power, consisting of the covariances computed for CMB. For this we use the CMB block of the covariance matrix for case 6. (The result using case 5 is similar.) We compute the derivatives of the power spectrum using CMBFAST. Then for primordial amplitude and slope, we have the following covariance matrix if we do not clean out kSZ. C(l < 2000) = A n A 2 × 10−2 −4 × 10−3 n 1 × 10−3 However, if we can extract the CMB, the covariance improves. A n −3 −4 . C= A 4 × 10 −9 × 10 2 × 10−4 n . (3.25) (3.26) Thus the covariance increases by roughly a factor of 5 for A and n if the kSZ is not extracted. If inhomogeneous reionization raises the kSZ power, shortening further the lever arm of the CMB, it becomes more important to clean out kSZ. To gain this additional information from the CMB on scales l > 2000 we need a template for kSZ. This template will come from simulations, and may be calibrated 62 on the kSZ measured at higher l. However, there will be uncertainties in the template because of missing physics, such as patchy reionization. Therefore the power spectrum analysis should be checked in as many ways as possible against the data. One way we explore in the next section is using the non-Gaussian information. 3.5 One-point distribution function analysis CMB and kSZ cannot be distinguished in a multifrequency analysis. We have seen this in our power spectrum analysis in the previous section, in which we could only estimate CMB separately from kSZ by knowing the shape of the kSZ power spectrum beforehand. However, if the data do not represent a Gaussian random field, and our data do not, then the power spectrum does not extract all the information available. Statistics which highlight non-Gaussianity are also particularly interesting because in the standard picture, primary CMB fluctuations are Gaussian. Any nonGaussianities indicate deviations from the standard picture for primary fluctuations, or in our case, indicate the presence of secondary sources of anisotropy. In principle then, it is possible to use non-Gaussianities to tease apart the CMB and kSZ. In the following, we show this is possible, and explore one method for doing so. The statistic we choose is the histogram of pixel temperatures, that is, the 1-point probability distribution function (pdf). In the context of weak lensing, the pdf provides stronger constraints than simple statistics such as skewness (Jain et al., 2000). This makes sense because the pdf incorporates all the moments of the field evaluated at a single point, such as variance, kurtosis, and so on. For a Gaussian field the pdf is Gaussian. For a highly non-Gaussian field, such as kSZ, the pdf function may show strong tails. The information in the pdf is complementary to the information in the power spectrum. The pdf contains information about non-Gaussianities, but not about spatial correlations, because it is evaluated at a single point. The power spectrum contains information about spatial correlations, but not about non-Gaussianities. Below, we devise a method for examining the pdf of our maps to extract parameters 63 3e-11 Total @ 217 GHz IR Sources CMB, kSZ, CMB+kSZ Radio Sources Noise Filter l(l+1) Cl / 2π 2.5e-11 2e-11 1.5e-11 1e-11 5e-12 1000 2000 3000 4000 5000 6000 l Figure 3.9: The spatial filter used for the non-Gaussian analysis and the components at 217 GHz. The vertical scale of the filter in this plot is arbitrary. of the component signals, given a set of templates. Then we apply this method. We want to separate CMB from kSZ, so we consider the sky at 217 GHz. This is simplest because it minimizes the signal from thermal SZ. We want to focus on scales where neither the CMB nor kSZ completely dominates the other. Thus, before evaluating the pdf, we spatially filter the maps. In Fourier space, we multiply the maps by a Gaussian centered at l ≈ 3600 with standard deviation σ l ≈ 750. In figure 3.9, this filter is shown with the the power spectra of the various components of the sky at 217 GHz. In the next sections, we describe the pdfs of the major signals at 217 GHz in this window, describe our method for separating kSZ from CMB, and apply it to our simulated maps. We make no claim that this is an optimum method for separating CMB and kSZ. At this point, this method is intended as a demonstration of the power of using simple statistics which are sensitive to non-Gaussianities. 64 3.5.1 Signal one-point distribution functions Here we examine the pdfs of the signals with the goal of identifying features that distinguish CMB and kSZ. In addition, the infrared sources will be an unavoidable contaminant at 217 GHz. The radio sources are not very significant, and as we have mentioned, may be helped by additional masking. The noise after spatial filtering is very low. After spatial filtering, we show the histograms for the major signals at 217 GHz in Figure 3.10. The pdf of the unlensed CMB is Gaussian. After lensing, the distribution is still Gaussian because lensing simply remaps locations on the sky. After spatial filtering, the lensed CMB need not have a Gaussian distribution function, but we find the deviations in this case to be insignificant. The distribution of kSZ should be symmetric, because there should be roughly as many clusters moving towards us as away. We expect the pdf of kSZ, or any signal sourced by discrete objects, to have a non-Gaussian point distribution function. This is because the signal strength is highly localized: the values of pixels are extreme in the direction of the object, and zero elsewhere. The kSZ histogram shows tails well beyond that expected from a Gaussian pdf of the same variance. A pdf with power law tails, such as Student’s t-distribution, is a much better fit, although we find the t-distribution to be too poor of a model to apply our method. We might expect the point sources to have a non-Gaussian pdf, since they are positive definite, which would seem to argue for strong skewness. This is counteracted by two effects. For a high density of sources, several sources are included in each beam, and the non-Gaussianity is diminished due to central limit theorem considerations. In addition, our spatial filtering removes long wavelength modes, causing ringing about each source which can wash out skewness. We show the histogram for infrared point sources, after spatial filtering, in figure 3.10; it is nearly Gaussian. A wider window function reduces ringing, and produces a point source pdf which is more positively skewed. 65 2.5x105 2.0x105 CMB kSZ IR sources Combined 5 1.5x10 5 1.0x10 5.0x104 5 10 4 10 103 102 1 10 100 -2 -1 0 1 2 -2 -1 0 1 2 -2 5 -1 0 1 2 -3 -2 -1 0 1 2 3 ∆T x10 (K) Figure 3.10: The histograms of the major signals at 217 GHz, after spatial filtering. Each pair of plots shows the same data, except the vertical axis above has a linear scale and below has a log scale. The distributions have been normalized to unity. The CMB is well fit by a Gaussian. kSZ is poorly fit by a Gaussian with the same variance. Infrared point sources are slightly skewed in the positive direction, compared to a Gaussian. At right we show the histogram of the sum of these maps, along with our best fit pdf from fitting and convolving the template pdfs, as described in the text. We see that the tails of the kSZ distribution are a unique feature with which we can try to distinguish it from the other signals. 3.5.2 Likelihood Given a set of pixelized data, we want to examine the likelihood of a trial distribution function. Let Pa (∆T ) be a trial pdf, which depends on a set of parameters a. For now, think of it as an abstract parameterized model for our data histogram. We construct our trial pdf in the next section. For a set of independent pixels ∆T i , the likelihood of parameters a is given by L(a) = Y Pa (∆Ti ), (3.27) i which is just the probability of the first pixel, times that of the second, and so on. Computationally it is more convenient to work with the logarithm of the likelihood. We also work with a binned pdf, so we combine all the pixels for each bin, and sum 66 over bins: log L(a) = X Nj log Pa (∆Tj ), (3.28) j where Nj is the number of pixels in histogram bin j, and ∆T j is the value at the center of bin j. This gives us the likelihood of parameters a for model pdf P a (∆T ) in terms of the histogram of pixel values Nj . If the pixels are not all independent, this is not the correct likelihood function. The true likelihood would take into account the covariance between pixels (or alternatively, between bins). However, we find in practice that the peak of the given likelihood function in (3.28) is a good enough estimate for this application. The correlation length of our maps after spatial filtering is much less than the size of the maps. Thus our maps contain a large number of uncorrelated patches. If these patches are thought of as larger, independent pixels, then the likelihood we have written has the right peak, but may have a wrong normalization and width. In our method, we use only the location of the likelihood peak, while the errors are assessed using Monte Carlo simulations. Our goal is to find the set of parameters a which maximize the likelihood. We use Powell’s direction set method (Press et al., 1986). Because the data (the histogram of pixel temperatures) is noisy, the likelihood function is plagued by local maxima. We repeat our maximum likelihood search with randomized initial positions and search directions in parameter space, to ensure that the maximum we find is reliable. We compute errors on our parameters by the bootstrap method. We sample our sky map realizations with replacement. From this set of sky maps (which may include some duplicates) we generate a synthetic data histogram. From it we find the best-fit parameters for that particular set of sky maps. We repeat the process with a new set of sky maps many times, and use the distribution of the best-fit parameters to describe the probability distribution. 67 3.5.3 Separation by distribution function To evaluated the likelihood function (3.28), we need a parameterized model pdf for the data histogram at 217 GHz. To get unbiased estimates of the parameters, our model must faithfully reproduce the data. We proceed by making models of each signal in turn. The pdf of the sum of independent random variables is given by the convolution of the pdfs of the individual random variables. Therefore, by modeling the pdfs of the individual signals, we can produce a trial pdf for the total data by convolution. We use binned functions for our models, and evaluate the convolution with fast Fourier transforms. Our procedure is not very sensitive to the width of the bins unless the model pdfs contain sharp features. In this application, our pdfs are smooth functions, and we find that roughly one thousand bins covering the pdf suffices. We make a 2 parameter model for the pdf. Our parameters are the quantities we seek to estimate, in this case the standard deviation of the CMB, which we label σ CMB , and the standard deviation of the kSZ, which we label σ kSZ . For the model pdf of the 2 . For the kSZ pdf, we use the CMB we simply use a Gaussian with variance σ CMB actual histogram of the kSZ as a template. In the real experiment, this would not be available, but could be based on simulations. To generate a distribution with a given standard deviation we rescale the template, preserving the distribution’s unit integral. For the model pdf of the point sources, we again use the actual point source histogram as a template. If we understand the point source population well enough to make a good power spectrum analysis, such that we wish to attempt a separation of CMB and kSZ, the we should be able to model the point source pdf. Note that our non-Gaussian analysis will provide a constraint on σ CMB and σkSZ individually. Contrast this to the 2 , albeit 2 + σkSZ power spectrum analysis, which constrains only the total variance σ CMB in a narrower window of multipoles. We apply our maximum likelihood parameter estimate to this model. The result is a > 10σ detection of the CMB and a > 7σ detection of the kSZ in our window (see Figure 3.11). The estimate is somewhat degenerate in the direction which preserves 68 2.8 2.6 6 σkSZx10 (K) 2.4 2.2 2 1.8 1.6 1.4 1.2 2 2.2 2.4 6 2.6 σCMBx10 (K) 2.8 3 Figure 3.11: Separating CMB from kSZ using the pdf. This plot shows the result of our fitting procedure for the standard deviation of the kSZ and CMB. The “+” shows the best fit point. The “×” marks the true value. The dots show the distribution of the fits for 1000 bootstrap realizations. The contour denotes the 68% region computed from the covariance of the bootstrap points (χ 2 with 2 degrees of freedom). The fine dotted line shows the circle of constant total variance which passes through the true value. 2 . 2 + σkSZ the total variance σCMB The power of this method arises from the template for kSZ, which fixes the relation between the kSZ pdf’s variance and tails. If this relationship is not fixed, then a narrow kSZ distribution with wide tails could mimic a wide kSZ distribution with narrow tails, by adjusting the variance of the CMB pdf to compensate. Also, it is crucial to have the point source contribution constrained. Otherwise it would be degenerate with the CMB. We must depend on simulations for the relationship between the variance and tails of kSZ. This non-Gaussian method may also be viewed as a consistency check for kSZ template simulations. A simulation of kSZ must both satisfy the power spectrum at high l, and provide the correct non-Gaussianity where kSZ and CMB overlap. If the 69 measurements differ, one should not have confidence in the template power spectrum for kSZ provided by the simulation. 3.6 Conclusions We have developed and examined sky simulations using the frequency channels, beams, and detector noise anticipated for ACT. Any high-resolution CMB working at the same three frequencies should produce similar results. The simulations include CMB, kSZ, SZ, and radio and infrared point sources. Secondary anisotropies dominate the CMB over a wide range of interest. These in turn are heavily polluted by point sources. Infrared point sources must be dealt with in the power spectrum, and cannot be completely dealt with by masking. While we included the latest constraints on the secondary anisotropies, significant uncertainties on the amplitude and correlation properties remain and could affect the conclusions reached in this paper. Particularly uncertain is the contribution of patchy reionization to kSZ. An optimal multifrequency filter can extract the power spectrum, but the results are subject to assumptions put into the analysis. We considered several cases, making different assumptions about our knowledge of contaminating extragalactic point sources. For optimal extraction knowing both the mean and scatter of the frequency dependence of point sources is crucial. If we do not use a prior model for the frequency dependence of the point sources we may still extract the power spectrum assuming spatial correlations are known, but the results may still be biased. Furthermore, using only frequency information in a power spectrum analysis, kSZ and CMB cannot be separated using data alone. This degrades the ability to measure the CMB at l > 2000 and to determine parameters from it. We also explored the ability to use the non-Gaussian information to constrain further the individual components. Non-Gaussianities in the histogram of pixel temperatures at 217 GHz should be observable by ACT. Non-Gaussianities created by kSZ can allow one to distinguish it from Gaussian CMB on scales where the two cannot 70 be separated with the power spectrum analysis, assuming that simulations are able to provide an accurate template of kSZ. It is unclear how well a kSZ template from simulations will work when applied to, for example, kSZ created by patchy reionization. Therefore, our results are preliminary and more investigation of these techniques is needed. Combining the power spectrum template analysis with the non-Gaussian template analysis allows one to perform consistency checks among these methods. Our results suggest that extracting the primordial power spectrum information at high precision from small scale primary CMB will be challenging. Scatter in the frequency scaling of point sources, as well as their possible correlations, make the point source separation from CMB difficult. Assuming this is accomplished, removing kSZ is even more difficult and can only be done with reliable templates for the power spectrum and for non-Gaussian signals. The final precision is likely to be limited by the modeling uncertainties and not by the statistical precision of the observations. It is hard to prognosticate the final accuracy given how uncertain the kSZ amplitude from patchy reionization is. In this context it is worth mentioning that CMB polarization may have less severe secondary anisotropies than the CMB temperature anisotropy and so may ultimately be more promising as a tool to study the primordial power spectrum fluctuations on small scales. Chapter 4 Progress towards simulation of reionization with radiation transport 4.1 Introduction Two probes, the optical depth measurement from the CMB and the observation of Lyman-α Gunn-Peterson trough in z > 6 quasars, are the first of many that cosmologists hope will shed light on the processes of the reionization epoch (Loeb & Barkana, 2001; Barkana & Loeb, 2001). The James Webb Space Telescope 1 (launch estimated 2011) will study several of these effects. Bright sources should have Lyman-α halos, photons from the source caught in the opaque IGM until the Hubble expansion took them out of the transition resonance. Photon heating at reionization should suppress the formation of low mass galaxies: the light will unbind gas from small halos. As a result, galaxy counts and the low end of the luminosity function should drop after reionization, in turn prompting a drop in the total star formation rate. Other phenomena are probed by different experiments. Redshifted 21 cm emission from the neutral 1 http://www.stsci.edu/jwst/ 71 72 hydrogen hyperfine transition may allow tomography of ionizing regions. Radio telescopes such as PAST (Peterson et al., 2004), LOFAR 2 , and SKA3 may observe this effect, if the challenges of observing at ∼ 100 MHz can be met. As a final observational signal, reionization could cause a CMB Doppler anisotropy, commonly called patchy reionization. This second order effect, requiring both a velocity and electron density perturbation, is small if the reionization of the universe is uniform. Estimates of this anisotropy are still uncertain (Santos et al., 2003; McQuinn et al., 2005), though it may be observed by ACT (Kosowsky, 2003), SPT (Ruhl et al., 2004), and Planck. 4 ACT and SPT have beams of 1–2 arcmins, and Planck has beams from 5–33 arcmin. For scale, 1 arcmin corresponds to about 250 kpc (physical) at z = 10 in a flat ΛCDM cosmology (or 2.8 Mpc comoving). Recent work suggests that spatial fluctuations in the reionization process may be significant, indicating the scales where ionized regions overlap may be many tens of comoving Mpc (Barkana & Loeb, 2004; Furlanetto et al., 2004), so that a box of 100 Mpc or more may be necessary to capture the large scale fluctuations. Lacking large scale cosmic variance, a simulation in a smaller box will naturally be more uniform; reionization will therefore proceed more quickly once it begins. The box size requirement poses a problem. The necessary dynamical range, from sources on kpc scales to ionization fluctuations on 100 Mpc scales, is beyond current computational resources, and prevents a direct simulation of these phenomena at all relevant scales. If we want to explore the large scale features of reionization, we are forced to abandon our desire to simultaneously simulate the small scale physics. In this work we begin to develop a framework to allow semi-analytic models of small scale reionization processes (halo formation, source intensity, clumping, etc.) to incorporate long range radiation transport. We use numerical simulations as a basis to model the population of ionizing sources, developing criteria necessary for an algo2 http://www.lofar.org/ http://www.skatelescope.org/ 4 http://www.rssd.esa.int/index.php?project=PLANCK 3 73 rithm to transport the radiation (section 4.2). We examine the intergalactic medium (section 4.3). We then extend an existing transport algorithm so that it satisfies these criteria (sections 4.4 and 4.5). Our work is similar in spirit to the work of Zahn et al. (2005) in that numerical simulations give the large scale fluctuations and we fold the details of the sources and IGM into semi-analytic models. Although we use different semi-analytic models, the chief conceptual difference between the thrust of this work and theirs is that we seek to transport explicitly the radiation from source to absorber, while this function is part of their modeling. 4.2 Source model The population of ionizing sources needs to reflect the fluctuations in density in our simulation box. We assume that these sources reside in dark matter halos, and that the ionizing luminosity of these sources depends only on the mass, the dynamical time of the halo, and the time since halo formation. This recipe has two ingredients. First we describe the halo population of each cell, which depends on the cell’s density, specifying the number of halos, their masses, and their formation times. Then we assign luminosities to these sources. 4.2.1 Halo population model We begin with a 2563 simulation in a 100 h−1 Mpc comoving box provided by Iliev. This simulation contains dark matter only. We have 32 simulation outputs from z ∼ 15 to z ∼ 6. At yet higher redshift, we extrapolate the overdensity using linear theory. The mass function describes our halo population. We use a prescription by Mo & White (1996), which allows one to write the mass function as a function of density in terms of the conditional mass function. The conditional mass function specifies the probability for a halo to be later incorporated into another halo. In this case, the prescription imagines the matter in a simulation cell will, sometime in the future, 74 collapse to form a single large halo. Based on this assumption, the conditional mass function gives the mass function at the present time. We use the extended Press-Schecter (EPS) formalism (e.g. Lacey & Cole, 1993), which gives the conditional mass function as f (s, δ1 |S, δ0 )ds = δ1 − δ 0 ds (δ1 − δ0 )2 p , exp − 2(s − S) (s − S) 2π(s − S) (4.1) where s and S are variances in density in a top-hat window with a scale given by a mass (i.e. S = σ 2 (M )). Symbols δ0 and δ1 represent the density critical for spherical collapse at two epochs, extrapolated to the present day using linear theory. (In an Einstein-de Sitter universe, δ0 ≈ 1.686(1 + z0 ), and similarly for δ1 .) Stating the meaning of the conditional mass function precisely, it gives the probability that a fluctuation with variance (on mass scale m) in ds about s at z 1 has collapsed, given that the same fluctuation will have variance S at z 0 (on scale M ). Using n(m, z1 |M, z0 ) dm = M dm f (s, δ1 |S, δ0 ) ds m ds (4.2) casts the conditional mass function in more familiar form, as a number per unit mass. The key step in the Mo & White (1996) prescription specifies the effective threshold for collapse for a region with non-linear overdensity δ. Sheth & Tormen (2002) modify the original formula to make it applicable to ΛCDM cosmologies of interest: δ̃0 (δ, z0 ) = 1.35 δ0 1.68647 − 1.68647 (1 + δ)2/3 1.12431 0.78785 − + . (1 + δ)0.58661 (1 + δ)1/2 (4.3) Substitution of δ̃0 into equation 4.2 yields the conditional mass function, and may also be interpreted as giving the redshift for collapse for a cell with mass M . Sheth & Tormen (2002) shows that the EPS mass function conditioned on density over-predicts the abundance of massive halos by a factor of two at redshift zero. However, in the 75 densest regions, this formalism appears to be accurate: by tracking the densest region in at simulation box back to z ∼ 50, Gao et al. (2005) show that in the highest density regions (1.7 < δ < 4.3), the EPS conditional mass function and Mo and White prescription predict the abundance of halos nearly to the Poisson accuracy of the simulation. After we compute the mass function for each cell, we next need to compute halo formation as a function of mass. We mimic the approach of Chiu & Ostriker (2000), basing the halo formation on the time rate of change of the mass function and a computed halo destruction rate. The change in the mass function is straightforward. In a mass bin ∆m wide, h∆N i = Z m+∆m m [nt+∆t (m) − nt (m)] dm (4.4) is the change in the mean number of halos after time ∆t, where n t (m) is the mass function in the cell computed at time t with our previous formalism. The mean number of halos changes because new halos form and because old halos are destroyed by mergers into larger systems. For the destruction rate we use the prescription of Sasaki (1994). This assumes that the process of halo destruction is independent of scale. Between two epochs, the probability of a halo surviving is given by the ratio of the linear growth factors: psurvive (t|t0 ) = D(t0 ) . D(t) (4.5) Thus the mean number of halos destroyed (that is, incorporated) during a time step is simply, D(t) h∆Ndest i = hN i × 1 − D(t + ∆t) , (4.6) regardless of mass. From the destruction rate, we compute the halo formation rate necessary to maintain the halo mass function as it changes in time. The mean number of halos formed in this time step is therefore h∆Nform i = h∆N i + h∆Ndest i. (4.7) Using the mean value, we sample from a Poisson distribution to compute the halo formation ∆Nform (m, tf ) in mass bins as a function of formation time t f . Changing 76 the random seed for the Poisson sampling ultimately has negligible effect on the mean source luminosity. We find that 500 mass bins are sufficient. As expected, the mass function of the semi-analytic halo population aggregated over the whole simulation box agrees well with the Press-Schecter mass function. 4.2.2 Source luminosity For halos above a minimum source mass m min = 106 M , we spread the UV luminosity of the halo over a dynamical time (again following Chiu & Ostriker, 2000, and references therein) mc2 `(m, t, tf ) = ε tdyn t − tf tdyn t − tf exp − , tdyn (4.8) where ε is an efficiency factor and the dynamical time is based on the virial density at the formation epoch: tdyn = (3π/32Gρvir (tf ))−1/2 . We leave the efficiency for photons to be created, escape their source halo, and proceed into the IGM as a free parameter in this model. The total luminosity for the simulation cell is integrated over the halo formation and mass: L(t) = Z t dtf 0 " X mi # ∆Nform (mi , tf ) `(mi , t, tf ) . (4.9) Figure 4.1 shows the total source luminosity accumulated over the whole box as a function of redshift. The luminosity of sources consistently rises, with a rapid increase between the redshifts of 18 and 10. The wiggles in the luminosity at z < 9 are not due to the Poisson sampling of sources, for they prove independent of the random seed. In figure 4.2, we plot an accumulated distribution of the source luminosity. This helps us gauge the importance of the ray-tracing algorithm’s scaling with regard the number of sources. On the same figure we also show two hypothetical cases for contrast: first where sources are uniformly distributed and second where a few sources dominate the total light. From this plot we can see that in our application (100 h −1 Mpc comoving box at 323 resolution), we are not dominated by a few sources. The brightest cells at z ∼ 6 account for less than twenty percent of the total luminosity. To 77 1.2 3 10 × L/εVc (solar masses/second/Mpc ) 1 0.8 2 0.6 6 0.4 0.2 0 6 8 10 12 14 16 18 20 22 24 26 z Figure 4.1: Source luminosity density, scaled by the efficiency, as a function of redshift, computed in our semi-analytic source model. get 90 percent of the flux, we need to include the brightest 77 percent of cells at z ∼ 6. This works out to about 25,000 sources, showing that the assumption (necessary, as we will see, for efficient ray-tracing) that only a few cells produce the bulk of flux is not valid in this context. As we increase the number of cells, the number of sources will increase (although the fraction of cells containing sources may decrease). 4.3 Intergalactic Medium For simplicity, we consider only hydrogen reionization. We ignore scattering and secondary photons from recombination. The absorption function, also called the photoionization opacity or optical depth per unit length, κ = σ PI nH 0 has units of inverse length. The cross section for photo-ionization of neutral hydrogen is σPI = 6.3 × 10−18 (ν/ν0 )−3 cm2 (4.10) 78 100 h-1 Mpc box, 323 resolution 1 z=6 z = 15 uniform few sources Fraction of light 0.8 0.6 0.4 0.2 0 0 0.2 0.4 0.6 0.8 1 Fraction of cells Figure 4.2: The accumulated distribution of light in the box. We bin and rank the cells by luminosity and plot parametrically the fraction of dimmer cells versus the fraction of light coming from dimmer cells. Left-to-right and bottom-to-top, the contribution comes from brighter sources. A horizontal line marks ten percent of the light. To include ninety percent of the luminosity, one must include all cells to the right of (brighter than) where the distribution crosses this mark. We show our source model at z = 15 and z = 6, and also two hypothetical cases. A uniform distribution of light is a straight line y = x on this graph: the fraction of cells always equals the fraction of light. When a few bright sources dominate the light, the majority of cells contribute little, and the graph hugs the horizontal axis until the fraction of cells gets close to one. 79 where ν0 = 3.29 × 1015 Hz is the Lyman limit for ionization. This and other process rate coefficients are compiled in Cen (1992) and Maselli et al. (2003). Assuming gas traces matter, the number density of neutral hydrogen is given by nH 0 = (1 − x)(1 + δ)n̄H = (1 − x)(1 + δ)(1 − Y )Ωb ρc µH (4.11) where x is the ionization fraction, δ is the local overdensity, n̄ H is the mean hydrogen number density (neutral and ionized), Y is the mass fraction in helium, and µ H is hydrogen’s molecular mass. For a flat cosmology with h = 0.7 and Y ≈ 0.24 from big bang nucleosynthesis, the comoving opacity (optical depth per unit length) is κ = 3.62(1 − x)(1 + δ)(ν/ν0 )−3 (1 + z)2 Mpc−1 (comoving) (4.12) In physical units, the opacity grows faster by one power of (1 + z). For completeness we include now the equations by which the ionization in the IGM evolve. We break our ionization equation into three pieces for convenience: dx = P (x) − R(x). dt (4.13) P (x) the photo-ionization (determined by the radiation transport) and R(x) is the recombination term. Each has units of inverse time. The recombination parameter is (4.14) R(x) = C a(T )(1 + δ)nH x2 where C is the dimensionless sub-grid clumping factor (which enhances recombinations) and the temperature-dependent reaction coefficient is a(T ) = 8.40 × 10 −11 T 1K −0.5 T × 1+ 106 K T 103 K 0.7 !−1 −0.2 cm3 s−1 (4.15) The magnitude of the opacity (equation 4.12) is relevant for the radiation transport algorithm. For many radiation transfer applications it is necessary, or convenient, to 80 divide the simulation box into optically thin cells. In this application, this is an impossibility. For an average density, neutral cell to be optically thin at z = 10, the dimension of the cell l must be small, l 230 kpc comoving. A 100 Mpc comoving box then re- quires many more than (5 × 104 )3 ≈ 1014 cells, well beyond what is computationally feasible. Therefore we conclude that any radiation transfer algorithm appropriate for reionization on 100 Mpc comoving scales must be able to handle optically thick cells. 4.4 Radiation transport The energy of radiation passing through area dA perpendicular to n̂, headed in a direction within dΩ about n̂, in a frequency interval dν about ν, during a time interval dt at t, is given by I(~x, n̂, ν, t) dA dΩ dν dt. This define the specific intensity I. Two features of the radiation transport problem complicate the basic goal of computing the specific intensity. First, the computational domain is large (position × 3 + direction × 2 + fre- quency + time). Second, the radiation field cannot be computed locally. Ray-tracing is a brute-force solution. On a cubical grid with N cells, conventional ray-tracing scales as N 2 (every cells transmits to every other). This is expensive, although several schemes which reduce the number of rays cast can help (Abel & Wandelt, 2002; Maselli et al., 2003; Rijkhorst et al., 2005), particularly if only a small fraction of the cells contain sources, in effect reducing the pre-factor for the scaling. Another approach is the “optically thin variable method Eddington tensor” (OTVET) method, wherein some moments of the radiation field are approximated. The OTVET method scales better as N log N , and conserves the flux of transmitted photons, but not the direction, which allows flux to propagate behind optically thick barriers. We are particularly interested in large scales (∼ 100 Mpc comoving) so even at 128 3 resolution our cells are quite large. We have seen that a successful radiation transport algorithm will handle a large number of sources, and be accurate when the opacity is high. Neither conventional ray-tracing nor the OTVET technique provides an ideal tool. 81 4.4.1 Cen’s Flux Method Here we explore another appealing approach to the radiation transport problem, suggested by Cen (2002). The scaling, roughly as N (log N ) 2 , improves on conventional ray-tracing, and more importantly, is independent of the number of sources. We omit many of the details here, but will explain the main contrast between this and more conventional ray tracing schemes. Conventional ray-tracing schemes broadcast radiation: a source emits light along rays which propagate throughout the space. As distance increases, the number of rays necessary to intersect every cell increases rapidly. By contrast, the method of Cen (2002) focuses on each cell as a radiation receiver. Every cell gathers light from a discretized set of solid angles, and there is no reason to refine this set at large distance. About each cubical grid cell, the sources and the intervening absorbers are mapped to a discrete spherical coordinate system. Cen’s method computes the flux, F (~x, n̂, ν, t) = Z I(~x, n̂0 , ν, t) dn̂0 , (4.16) ∆Ω for some point in each grid cell, where ∆Ω is the solid angle element in direction n̂ we are considering. The origin of the spherical coordinates defines the point in each cell where we compute the flux. A convolution permits us compute simultaneously the flux and optical depth contribution to each grid cell from a particular spherical element or “window function” (a solid angle element and radial range, see Fig. 4.3); this is done efficiently on the grid with a fast Fourier transform. Working outward, we then integrate over radial direction to compute the flux (as a function of direction) at the same point in each cell. The light from a spherical element is attenuated by all the material in elements interior to it. The flux contribution falls quickly with distance from the source, as r −2 at least, and exponentially faster if there is absorption. Therefore this form can allow considerable savings: an increase in the number of grid cells does not require an increase in angular resolution, and requires only a slow increase in the radial resolution if the radial coordinate coarsens as the 82 Figure 4.3: The key step in Cen’s flux method for radiation transport, shown in a twodimensional schematic. For the solid angle under consideration, example cells A, B, and C (and the others) collect flux from sources in spherical volume elements. Two such elements are drawn for each example cell. Using the spherical element’s shape and position as a window function, the emissivity j and opacity κ of the analogous spherical elements for each cell are computed simultaneously by convolution. distance from the receiver increases. This method parallelizes trivially over the solid angle elements; all are computed independently. (Many ray tracing codes also parallelize straightforwardly.) We have implemented Cen’s flux method using the Message Passing Interface 5 (MPI) for use on a Beowulf cluster. At each radiation transport computation, each node ideally computes the contribution to flux from one solid angle element. (In practice on the cluster we use, fewer nodes are available, requiring each to compute more.) This implementation may have many uses, not limited to the reionization context, but in any situation which requires very large numbers of sources. Although here we consider the light travel time to be zero for simplicity, this framework admits the possibility of a non-zero travel time using retarded source and absorption functions. (Wyithe & Loeb (2004) argue that this effect is vital, and determines both the epoch and duration of reionization.) One objection to this method is that at large distance shadows are cast imperfectly, 5 http://www.mpi-forum.org/ 83 due to the fixed angular resolution. This is not a fundamental flaw in the method, but merely a product of finite computational resources. The angular resolution can be chosen as fine as desired, if sufficient resources are available. For the example calculations shown here, we break the sphere into roughly 200 equal area solid angle elements, which for convenience, have constant lines of latitude and longitude for their borders. 4.4.2 Photon Absorption As we have noted, Cen’s method is a method for computing the flux F for some point in each cell. For reionization simulations, this is not the quantity of immediate interest. Rather we are concerned with the number of photoionization in each cell, or equivalently with the amount of power absorbed. We denote this quantity by A, which has units of power per unit frequency. Computing the absorption from the flux is generically non-trivial. This is in contrast to more conventional ray-tracing methods, which may explicitly balance the number of photons absorbed and emitted: when a ray crosses a cell, the number of photons deposited equals the number of photons removed from the ray. The simplest case for computing absorption from flux occurs when the optical depth per cell is low, and the emissivity and opacity are very uniform on small scales. In this case, the flux is roughly constant throughout the cell, and the absorbed power per unit frequency is well approximated by A ≈ κF l 3 , where the flux F is evaluated anywhere inside the cell (say, at the center) and l is the dimension of the cell. This approxima- tion fails spectacularly when optical depth per cell is high or when the optical depth or emissivity show rapid spatial variations. This is chiefly because the flux is not constant throughout the cell; the radiation field can have large fluctuations on scales smaller than the grid size. When the optical depth per cell is high, the flux inside the cell falls exponentially along the photon path. If the emissivity or optical depth of adjacent cells are very different, the flux incident on different portions of the cell can also 84 be different. The results of employing the faulty “optically thin” approximation often leads to a severe suppression of the computed absorption. Simply boosting the absorption to account for attenuation in the cell can lead to computationally unstable results, problems with floating point overflows, and machine precision rounding errors. Unfortunately, reionization exhibits these challenging behaviors. A common scenario involves a source concentration ionizing a single cell, leaving the cell optically thin but surrounded by optically thick cells. If the absorption is suppressed by a faulty conversion of flux to absorption, the ionization front will stall at this stage, spoiling the computation. In the remainder of this part, we examine radiation absorption and design an algorithm which handles these difficult circumstances to a few percent accuracy. The specific intensity a distance s along the photon path evolves according to dI(s) = −κI + j, ds (4.17) where j is the emissivity and the quantities I, κ, and j are evaluated at the point ~x +sn̂. To simplify the notation, we have suppress the time and frequency dependence of these quantities, as well as the n̂-dependence of I. This is a linear first-order equation, readily solved with an integrating factor. Here we consider only the simplified case of a homogeneous medium (κ, j constant), appropriate in the interior of a grid cell. Then this equation yields the solution I(I0 , s) = (I0 − j/κ) exp(−κs) + j/κ, (4.18) where I0 is the specific intensity at s = 0. This in turn may be integrated to yield the absorption per solid angle in a uniform medium, Z s I(I0 , s0 ) κ ds0 a(I0 , s) = 0 = (I0 − j/κ)[1 − exp(−κs)] + js. (4.19) Quantity a has the same units as I, and represents the energy deposited along the photon path, from 0 to s, per time per area per frequency per solid angle. Written in this form, a suffers no numerical problems in computations with large κ. 85 Figure 4.4: A schematic view of the absorption calculation for a grid cell for light traveling in direction n̂. The specific intensity I 0 is computed at a point on the surface of the cell where the flux is incident. The thickness of the cell at that point is labeled s. These quantities are sufficient to compute the absorption when integrated over the cell’s projected area, labeled S. The specific intensity is computed by Cen’s flux method, symbolized by the cones on the incident surfaces of the cell. To compute the absorption in a cell for a specific solid angle element ∆Ω in direction n̂, we take the specific intensity at the incident surface of the cubic cell, evaluate equation (4.19) using for s the distance through the cell, and integrate over solid angle and the projected area of the cube: Z Z dA A∆Ω = S dΩ a(I0 (~x), s(~x)), (4.20) ∆Ω where ~x refers to a point in S, which represents the two-dimensional projected area of the cube. This geometry is shown schematically in Figure 4.4. For convenience, we in practice integrate over the three faces of the cube where the light is incident, Z 3 Z X A∆Ω = dA dΩ a (I0 (~x), s(~x)) |êi · n̂| , i=1 face i (4.21) ∆Ω where êi indicates the surface normal vector for face i in the final term, which converts the cell surface area to the projected area. Now we can make the connection to flux and absorption. For a given solid angle ∆Ω and point on the surface of the cube, we compute the flux F using Cen’s method. Then we compute the absorption along the line of sight with equation (4.21), using 86 Figure 4.5: Grids for refining the absorption integral over the cell surface, seen from the point of view of incoming light. The first level of refinement evaluates only at the filled square. The next level reuses these computations, and evaluates additionally at the empty squares. Next sequentially are the circles and triangles. I0 = F/∆Ω. To compute the total absorption in a cell, we compute the integral over the cube face as outlined below, and sum over all solid angle elements, A= X k A∆Ωk . (4.22) The absorption calculation is also implemented in parallel through MPI, again with different nodes handling different solid angle elements. We now describe the integration over the cell face. We evaluate the integrand of equation (4.21) on a refinable square grid on each incident cell face. This grid is designed to make a two-dimensional analog of the one-dimensional trapezoidal rule. (Compare to chapter 4 Press et al. (1992) for the 1-d version.) The mth level of refinement divides the cube face into 4 m squares, and evaluates the integrand (2 m + 1)2 times. However, as in the 1-d case, the cube face grid is refined without wasting earlier computations. The evaluation points for the first several grid refinements are show in Figure 4.5. We use this trapezoidal rule as the basis for a Romberg integrator, employing Richardson’s deferred approach to the limit: we extrapolate successive refinements 87 Figure 4.6: Radiation absorption from a single source at the middle of the box, in a transparent cell. The emissivity of the source is j = 1 cell −1 , and the opacity of all other cells is κ = 1 cell−1 . The mismatch between absorbed and emitting radiation is 6 percent. The central cell is dark because transparent cells absorb no radiation. of the integral to the continuum limit. For this we use a polynomial extrapolation, computing the grid to at least the m = 2 level of refinement. We estimate the error by the difference of the extrapolated limit and the best refinement computed. 4.5 Tests and future directions The algorithm preforms well when the emissivity and opacity are uniform, yielding 2 percent error regardless of the optical depth per cell. We examine the difficult case of the isolated source in an optically thin cell surrounded by optically thick cells. We find that the m = 2 refinement (61 evaluations on the cell surface) is sufficient to get 6 percent accuracy, regardless of optical depth in the surrounding cells. The source cell is given an emissivity j = 1 cell−1 and no opacity, and the opacity elsewhere κ = 1 cell−1 in our first case. The result is shown in Figure 4.6. In this case the error is 6 percent too few photons. If w increase the opacity in the surrounding cells to κ = 100 or κ = 1000 we find 5 percent too many photons. A few percent accuracy in the radiation absorption is probably sufficient to address reionization at the current time. The models of the sources, for example, do not exceed this accuracy. We can achieve 88 higher accuracy by going to higher order in the integral approximation, but at a steep computational cost. The next application of this algorithm will be to apply it to a radiation transfer scenario, using the source model laid out in section 4.2, and evolving the box from a neutral initial conditions at high redshift. This requires some additional considerations. One is the size of the timestep for the reionization code. Photoionizations and recombinations affect the opacity, which affects the radiation transport, which affects the photoionization. Mellema et al. (2005) expresses the essential “Courant condition” for radiation transfer as follows: an ionization front should not cross a cell in a single timestep, requiring ∆t < l/vfront . They summarize several other conditions used in the literature as radiation transport timesteps. These include the cell light-crossing time ∆t = l/c ∼ 1014 (1 + z)−1 (l/(1 Mpc)) s, which is a conservative estimate, several times the ionization rate ∆t = Γ−1 , and a fraction of the time for the ionized fraction to change ∆t = x/(dx/dt). (Mellema et al. (2005) also present a ray-tracing method which relaxes these conditions.) The epoch from z ∼ 25 to z ∼ 6 lasts approx 2.5 × 10 16 seconds. Fewer than, but probably on the order of, 2 × 10 3 ((1 Mpc)/l) timesteps are required, based on the cell light crossing time. With the present algorithm, the radiation transport computation (at a single frequency) in a 323 box takes 1.3 total cpu-hours on cluster of 2.4 GHz Xeon processors. A 643 box takes 11.0 cpu-hours, roughly corresponding to the increase in the number of cells. Since this is the time for a single timestep, this algorithm will need to be sped up to be practical for reionization simulations. The repeated flux calculations on the cell surface offer one obvious way to increase the efficiency. At a distance from the cell, the window functions for a set of points on the cell surface coincide almost exactly. Thus the flux gathered from a distance is nearly the same for all the points on the cell surface. This opens the possibility of using a common origin to compute the flux from distant cells, while only recomputing the flux contributions from nearby windows. This should lead to a substantial improvement in speed. Appendix A Semi-analytic simulations of the Sunyaev-Zel’dovich effect A.1 Motivation Hydrodynamical simulations are much more costly than N -body simulations. A Eulerian hydrodynamical code for an SZ simulation is 10 slower than a Particle-Mesh (PM) N -body. ACT will generate a catalog of some hundreds of clusters selected by their SZ signatures. We wish to probe the SZ signal using many realizations of large volume boxes (several hundred Mpc per side) without the expense of hydrodynamical simulations. With this goal, in this chapter we marry together an N -body simulation and an analytic gas model to generate maps of the SZ signal. Recall the expression for the SZ signal, equation (2.6). If the electrons are not relativistic then the frequency dependence is x/2 , Θ (n̂) = −2y(n̂) 2 − tanh(x/2) SZ where x = hν/kB TCMB and the comptonization parameter y is given by the integral of the electron pressure along the line of sight: y(n̂) = Z σT ne (n) 89 kB Te (n) dr. me c2 90 To compute the SZ signal, then we need to know the the electron density n e , and temperature Te . One approach to SZ using an N -body only simulation is to consider these electrons to be part of matter halos. In an N -body simulation where the dark matter halos have been identified, we can apply an analytic model for the profiles of these halos, and thus compute the SZ signal. This is the approach we use in this section. A.2 Analytic halo profile Komatsu & Seljak (2001, 2002) developed an analytic model for the gas density and temperature profiles of dark matter halos. In those works, they predict SZ profiles and the SZ power spectrum. Here we supply the necessary details of their model and apply it to a catalog of dark matter halos. This model depends on three assumptions. First, the gas is in hydrostatic equilibrium with the dark matter gravitational potential. Second, the gas density traces the dark matter density in the outer portions of the halo. Third and finally, the gas in the halo follows a polytropic equation of state P gas ∝ ργgas , where Pgas and ρgas are the pressure and density of the gas, and γ is the polytropic index. Under these assump- tions, the gas density, temperature and pressure profiles of the dark matter halo can be computed from the dark matter density. Calling y gas the gas density relative to the central density, we have: ρgas (r) = ρgas (0) ygas (r) γ−1 Tgas (r) = Tgas (0) ygas (r) (A.1) In recent years, computational N -body simulations have provided important insights into the non-linear formation of structure by the dark matter. Navarro et al. (1997) (NFW) observed a ‘universal’ density profile for dark matter halos, suggesting that the profiles depend only on the mass of the halo, with relatively little scatter between them. They also noticed that the halo concentration parameter, describing the de- 91 gree of profile steepness (defined more precisely below) decreases with increasing halo mass, though the exact relation depends on the cosmology at hand. Several formulas have been suggested for this relation (Navarro et al., 1997; Bullock et al., 2001; Eke et al., 2001) which rely on fitting some mechanism to halos culled from N -body simulations. These studies found that the concentration dependence on the power spectrum is caused primarily by the slope of the linear power spectrum. Both this and the mass dependence of concentration can be explained by the epoch of the halo formation, since less massive halos form at higher redshift when the universe was denser and this leads to a more concentrated profile. N -body simulations suggest a family of dark matter density profiles: ρdm (r) = (r/rs ρs , + r/rs )3+α )−α (1 (A.2) where ρs is a characteristic density, rs is the radius where the profile has an effective power law index of −2, and −1.5 < α < −1 (Navarro et al., 1997; Moore et al., 1998). Here we use α = −1 (or NFW profile). The quantities we wish to compute are not sensitive to the inner parts of the halo. As is conventional, we re-parametrize ρ s and rs in equation (A.2) in terms of a halo mass and concentration. We define the virial radius r vir of a halo to be the radius of a sphere with some characteristic mean density, discussed below. Insisting that the mass contained inside rvir is M fixes ρs for a given rs . We introduce the concentration parameter, the ratio c ≡ rvir /rs . Often we reparameterize the radius with x = r/r s , which means that x(r = rvir ) = c. Using the NFW profile for the dark matter, Komatsu & Seljak (2001) found ygas (x) ≡ B ≡ 1 γ−1 ln(1 + x) 1−B× 1− x 1 −1 3 γ − 1 ln(1 + c) − η(0) γ c 1+c (A.3) The polytropic index γ and the normalization η(0) have yet to be fixed. These are worked out analytically in Komatsu & Seljak (2001), and we simply cite the fitting 92 formulas they provide: γ = 1.137 + 8.94 × 10−2 ln(c/5) − 3.68 × 10−3 (c − 5) η(0) = 2.235 + 0.202 (c − 5) − 1.16 × 10−3 (c − 5)2 (A.4) The temperature normalization is given by Tgas (0) = 8.80 keV × η(0) M/1015 h−1 M rvir /1 h−1 Mpc (A.5) and the gas is normalized at the virial radial to have the universal density compare to the dark matter: ρgas (c) = ρgas (0) ygas (c) = ρdm (c) × Ωb /Ωc . With this model, we can now specify the gas density and temperature profile of a halo with known mass, redshift and concentration. Integrating these quantities along the line of sight (equation A.1), we obtain the comptonization parameter, and so can compute the SZ signal per halo. We note the frequency dependence of SZ is only an approximation. We include the approximation of Shimon & Rephaeli (2004) for the relativistic effect. In this case we compute for each halo the profile of SZ temperature perturbation per unit length, then integrate the profile along the line of sight. A.3 Halo catalog Trac produced a series of N -body simulations with corresponding halo catalogs. The simulations were in 200 Mpc/h and 400 Mpc/h boxes with periodic boundary conditions, each with 9603 particles on a (4 × 960)3 grid. The redshift outputs correspond to every 100 Mpc/h radial distance from the observer. The information in the catalogs is quite extensive: they provide mass estimates based the mass within both ellipsoidal and spherical overdensities. The length and semi-major axes of the enclosing ellipsoid, the radius of the enclosing sphere, the halo’s position within the simulation box, and the halo’s peculiar velocity. For this application, only the mass (we use the spherical estimate) and the position are useful. 93 It would be additionally useful to know the radius of the enclosing sphere at two different overdensities. This would give us a measure of the halo concentration, which is necessary to fully specify the gas profile model. This has been suggested as an improvement in future halo catalogs. For the mass and redshift dependence of concentration we choose a power law parametrization: 10 c(M, z) = 1+z M M∗ (z = 0) −0.2 (A.6) where the nonlinear mass M∗ is the mass scale where top-hat filtered fluctuations in the dark matter overdensity field have RMS δ c , which is the threshold for spherical collapse. As an option in our code, we can add a scatter in ln c, as suggested by N -body studies Bullock et al. (2001), where they found that the 1σ scatter was ∆ ln c ≈ 0.18 regardless of mass up to M = 6 × 1014 h−1 M . A.4 Semi-analytic Sunyaev-Zel’dovich maps Map producing code proceeds by looping over halo catalog redshift slices with a defined field of view, projecting the SZ signal of each slice in turn. Each redshift slice has several characteristics chosen at random: the offset of the box in 3 dimension, the axes corresponding to the line of sight, and the orientation of the other 2 axes. Then the code examines the halo catalog, disregarding any halos which fall outside the field of view or below the mass cut. The 3 dimensional profile for the halo is computed, then integrated along the line of sight to give the 2 dimensional projected SZ. This is done on a small map, with fewer pixels, but at higher resolution than the final map. Finally, the small map is added to the final map. Figure A.1 shows the result of this procedure for mass cuts from 10 12 to 1014 h−1 M . The foremost feature is the remarkable amount of confusion in the map from the smaller halos. This is to the extent that the halos in the final map no longer look spherical, although by definition in our model here, they are spherical. As a future application, I wish to compare these maps to hydro simulations with a 94 Figure A.1: This is a series of SZ maps with varying mass cut. The images are 1.75 ◦ on a side. From left to right and top to bottom, the mass cut increases, excluding more halos: 1012 h−1 M , 1013 h−1 M , 3 × 1013 h−1 M , and 1014 h−1 M . The brightest cluster in this field is has central temperature 2 ×10 −4 K at 265 GHZ. halo catalog. This will allow a direct evaluation of the quality of the this model, and allow an evaluation of its usefulness in place of full hydro simulations. Appendix B Importing flat-sky simulations into the ACT geometry On several occasions in this work (see sections 3.3.1 and A.4), we have used maps from cubical simulation boxes output at different redshifts. The technique involves offsetting, flipping, stacking, and projecting simulation boxes down to a 2-dimensional map. This is fine for a rectangular map of a few square degrees, because the effect of the celestial sphere’s curved surface is minimal. However, the ACT survey strip is a ring around the entire celestial sphere at δ ≈ −55 ◦ , so curvature cannot be ignored. The simulation boxes typically have periodic boundary condition, so we can produce tiled maps of large area, the usual projections are flat, and there is no way to map the surface to the sphere without discontinuities. Fortunately, there is a method to map the simulations to the sphere with little distortion and only a single discontinuity, which then can be hidden in the galactic plane. We use a conical equal area projection (Calabretta & Greisen, 2002), which maps the surface of the sphere to a cone which intersects the sphere at two lines of latitude, called “standard parallels.” The length of these parallels is preserved by the projection, as is the area between them. We use spherical coordinates (φ, θ) on the celestial sphere, and Cartesian coordi95 96 Celestial Sphere C on e δ = −55 deg P ro je ct io n ACT Strip Figure B.1: A schematic of the ACT survey strip at δ ≈ −55 ◦ (not to scale), and the tangent cone, onto which the simulation images are projected. nates (x, y) on the cone. Using the standard parallels θ 1 and θ2 , the equations of the projection are x = Rθ sin(Cφ) y = −Rθ sin(Cφ) + Y0 (B.1) The offset for y so that the origin in the cone coordinate system occurs at (0, (θ 1 +θ2 )/2). The parameters of the projection are C = γ/2 Rθ = Y0 = 180◦ 2 p 1 + sin θ1 sin θ2 − γ sin θ π γ 180◦ 2 p 1 + sin θ1 sin θ2 − γ sin((θ1 + θ2 )/2) π γ γ = sin θ1 + sin θ2 (B.2) Each redshift box is then tiled onto the cone surface with a random offset and projection direction. 97 Simulation Box ACT Survey Strip Figure B.2: Mapping tiled simulation boxes onto a segment of the ACT survey in the conical equal area projection (distortion is not to scale). Appendix C Signal filtering in ACT maps This appendix discusses general signal filtering to produce maps in the presence of noisy backgrounds. Signal filtering will be vital for an experiment like ACT because of the series of scientifically interesting signals layered on top of one another. Example applications include identification of point sources, which will allow study of their frequency dependence and number counts, and finding SZ clusters, which may allow study about the growth of structure in the universe, and ultimately dark energy. Filters can be limitless in their complexity. However, complicated filters can make errors difficult to evaluate. In the next sections, we keep it simple, and discuss an optimal linear filter. In last section, we work an example, and filter for point sources in simulated ACT data. C.1 Linear (Wiener) filter Define the vector of data as d = Fs + n (C.1) where s is the vector of signals, F is the instrument response matrix to the signals, and n is the uncorrelated noise. For simplicity, we assume all of these quantities have zero mean. We define a linear estimate of the signal, s̃ = Wd, by minimizing the matrix of 98 99 reconstruction errors: (s − s̃)(s − s̃) † . The computation of this minimization (omitted here) proceeds by differentiating the errors by W and setting to zero. The estimate optimized in this way is called the Wiener filter: W = SF† [FSF† + N]−1 , (C.2) where S = hss† i and N = hnn† i are the signal and noise covariances. Note also that the term which is inverted is the data covariance matrix: C = hdd† i = [FSF† + N]. (C.3) The covariance of the estimate is also easy to write down: hs̃s̃† i = WCW† . (C.4) If the signals and noise are spatially uniform and isotropic, then the covariance matrices are diagonal in Fourier space (or spherical harmonic space, depending on the geometry) and depend only on |~k|; the covariances are simply the power spectra. Then filtering may be computed efficiently mode by mode: s̃~k = Wk d~k , (C.5) where the signal and data are still vectors to account for multiple signals and detector channels. These vectors are small, d~k has 3 components for ACT, compared to millions of components for d, so the matrix inversion necessary in the computation in W k is trivial. In a single frequency experiment with a single signal (statistically spatially uniform), the Wiener filter has a particularly simple form: s~k = F 2P Ps (k) d~k s (k) + Pn (k) (C.6) If we have no notion of frequency dependence of signals, which appears in the response matrix F, we can of course treat a multifrequency experiment like several single frequency ones. This amounts to an incorrect assumption about the form of the data 100 covariance matrix, but is sometimes a necessary staring point, as in our example signal filtering below. For ACT, the signals (except for the galaxy, which we ignore in this discussion) will be spatially uniform, but the noise variance will vary from place to place. On the ACT strip, the center of the strip will be observed longer than the edges, and the placement of the bolometer arrays inside the cameras will cause portions of the survey strip to have varying coverage in the different bands. When signals or noise are not uniform, then the matrix inversion in W = SF † C −1 becomes much less trivial, because the matrix is not automatically sparse in the Fourier basis. For filtering, we need to compute C−1 d, and in this case a better approach is to use a conjugate gradient method to solve the equivalent Cy = d (C.7) for y, by minimizing f (y) = 1 † y Cy − d† y, 2 (C.8) which has the preceding relation as its gradient. The multiplication Cy proceeds quickly only if the matrix C is sparse. We have already noted that the signal covariance S is sparse (in fact diagonal) in Fourier space. The instrument noise covariance N , on the other hand, is sparse in real space, with perhaps only a few neighboring pixels correlated. It therefore is more efficient to split the covariance into a number of pieces C = CFourier + CReal space + . . ., each of which is sparse in a different basis. Then the trial y can be distributed to each of these terms in the covariance, transformed into the appropriate basis, and multiplied in a sparse manner. The most difficult step is the transform to a new basis, which scales as N log N for fast Fourier and fast spherical harmonic transforms. 101 0.007 0.035 145 GHz 0.006 0.03 0.005 0.025 0.004 0.02 0.003 0.015 0.002 0.01 0.001 0.005 0 0 -0.001 0 2 4 6 8 10 -0.005 217 GHz 0 2 4 arcmin 0.05 8 1 265 GHz 0.045 10 145 GHz 217 GHz 265 GHz 0.9 0.04 0.8 0.035 0.7 0.03 0.6 0.025 0.5 0.02 0.4 0.015 0.3 0.01 0.2 0.005 0.1 0 -0.005 6 arcmin 0 2 4 6 8 10 0 100 1000 arcmin 10000 l Figure C.1: The first three plots show the point source filters in real space for the three ACT channels. At bottom right, the same filters, in Fourier space. C.2 An example: point source identification As an example application, here we present filters for finding point sources in the ACT maps. This will be one of the first steps in the ACT analysis, and will be necessary to establish the point source contamination and frequency dependence. We use three single-band filters. We make the simplifying assumption of uniform noise. The filtering is performed in Fourier space using (C.6), with filters built from the power spectrum in the three bands (compare Figure 3.2). Figure C.1 show the filters in real and Fourier space. The filters become narrower at higher frequency because the beam is smaller. 102 Figure C.2: Above, maps of ACT channels before filtering. Left to right are 145, 217, and 265 GHz. Below, the same channels after filtering to emphasize point sources. Appendix D Cross-spectrum contaminant estimator In this appendix we present a generalized version of the point source estimator of Hinshaw et al. (2003), using a vector notation. We discuss contaminants in general and do not mention point sources or SZ specifically. We assume we know the power spectrum shape and frequency dependence of all contaminants. Our data are the cross-power spectra from the experiment. Let vector D = {C li } be these spectra. The multipole (or multipole bin) is denoted by l and the cross correlation is denoted by i = W1W2, W1V1, etc. We use a Gaussian model for the likelihood L(D) of the set of cross-spectra: (D.1) −2 log L + constant = [D − hDi] † Σ−1 [D − hDi] where the covariance Σ = h(D − hDi)(D − hDi) † i can be written as Σ = {Σii ll0 }. The 0 constant ensures the likelihood is properly normalized. We postulate that the data D is the sum of the CMB and a number of contaminants. For each of these contributions to the data we need a model, and parameters for that model. Some of these parameters we will estimate, and some we will marginalize. First we consider the CMB. The parameters which describe the CMB are the band 103 104 powers: Cl . We organize these band powers into a vector C = {C l }. We will later marginalize over these parameters. To relate the CMB parameters to the data we use the window function w = {wlli 0 } for each cross spectrum channel pair. We may think of the window function as the response of the instrument to a given set of CMB Cl ’s. In the absence of noise and contaminants, the data would be given by the matrix multiplication D = wC. We use as similar notation for the contaminants. We describe each contaminant by a power spectrum shape, a frequency dependence, and an amplitude. We take the shapes and frequency dependences as given, and the amplitudes as parameters. We divide the amplitude parameters into two groups, those to marginalize and those to estimate. The amplitude parameters to marginalize we denote by vector A = {A α }, where α runs over the components to marginalize. The amplitude parameters to esti- mate we denote by vector B = {Bβ }, where β runs over the components to estimate. We may think of the shape and frequency dependences as the response of the in- strument to a given set of contaminant amplitudes. Thus the shape dependence already includes the influence of the window function. Similarly to the amplitudes, we divide the shape and frequency dependences into two groups. The shape and frequency i }. The shape and dependence of the contaminants to marginalize we write as S = {S lα i }. In the frequency dependence of the contaminants to estimate we write as Z = {Z lβ absence of CMB and noise, the data would be given by D = SA + ZB. If we include all components and uncorrelated noise, then we can write the expected data as hDi = wC + SA + ZB. (D.2) We are considering cross-spectra only, and have included no noise term. We substitute into the likelihood expression: −2 log L + const. = [D − (wC + SA + ZB)] † Σ−1 [D − (wC + SA + ZB)] . (D.3) 105 We seek to estimate B. This we accomplish by repeatedly completing the square in the likelihood expression, as follows. This allows us to easily integrate out the contaminants to our signal. Note that completing the square to marginalize out a component is equivalent to letting the covariance for that component tend to infinity (Rybicki & Press, 1992). Thus we disregard all information about that component. The likelihood may be rewritten as a quadratic equation in the CMB power spectrum, 1 −2 log L + const. = C† aC + (b† C + C† b) + c† Σ−1 c, 2 (D.4) where we have introduced the shorthand a ≡ w † Σ−1 w b ≡ −2w † Σ−1 (D − SA − ZB) c ≡ D − SA − ZB. (D.5) We complete the square for variable C: 1 −1 † 1 −1 −2 log L + const. = C + a b a C + a b 2 2 1 +c† Σ−1 c − b† a−1 b 4 Note that the final term may be rewritten in terms of our auxiliary variable c, −1 1 † −1 −1 † −1 † † −1 b a b=c Σ w w Σ w c. w Σ 4 (D.6) (D.7) Thus if we define the symmetric matrix −1 E ≡ Σ−1 − Σ−1 w w † Σ−1 w w † Σ−1 , (D.8) we may compactly express the log likelihood as a term which depends on the CMB C, and a term which does not. 1 −1 1 −1 † (D.9) −2 log L + const. = C + a b a C + a b + c† Ec. 2 2 R Let us define a marginalized likelihood, L C ≡ dC L. Completing the square has made it trivial to integrate the likelihood over all C. When we exponentiate and integrate, 106 the exp(−c† Ec/2) factors out of the integral, which integrates to some (unimportant) normalization factor. We find −2 log LC + const. = c† Ec = [D − SA − ZB]† E [D − SA − ZB], (D.10) where the normalization constant is different than before. Note the similarity to equation (D.3). The matrix E is the new inverse covariance, once the CMB is marginalized out. We wish to repeat this sequence to marginalize the variable A. Thus we write 1 −2 log LC + const. = A† fA + (g† A + A† g) + h† Eh 2 (D.11) where we have introduced f ≡ S† ES g ≡ −2S† E(D − ZB) h ≡ D − ZB. (D.12) Again we complete the square, now for variable A: −2 log LC + const. = † A + f −1 g f A + f −1 g 1 +h† Eh − g† f −1 g. 4 (D.13) Noting the last term may be re-written in terms of h, −1 1 † −1 † † † g f g = h ES S ES S E h, 4 (D.14) −1 F ≡ E − ES S† ES S† E. (D.15) we write with analogy to equation (D.8), 107 Define LC,A ≡ R dCdA L = R dA LC , and as before we integrate out our contaminants: −2 log LC,A + const. = h† Fh = [D − ZB]† F [D − ZB]. (D.16) To estimate the amplitudes B, we express the marginalized likelihood as 1 −2 log LC,A + const. = B† pB + (q† B + B† q) + r 2 (D.17) where we have introduced p ≡ Z† FZ q ≡ −2Z† FD r ≡ D† FD. (D.18) We complete the square one final time for variable B: 1 1 −2 log LC,A + const. = [B − p−1 q]† p[B − p−1 q] 2 2 1 † −1 +r − q p q. 4 (D.19) Now LC,A (B) ∝ exp[− 21 (B − hBi)† (ΣB )−1 (B − hBi)]. If we estimate B̄ ≡ (Z† FZ)−1 Z† FD, (D.20) we have hB̄i = hBi, which indicates our estimator is unbiased. Moreover, the estimator is unbiased even if the original covariance Σ from equation (D.1) is incorrect. This is shown by integrating the (flawed) estimator with the likelihood to take the ensemble average. We also note the covariance on the B parameter estimates, ΣB = (Z† FZ)−1 . (D.21) We note that neither the estimator B̄ nor its variance ΣB are sensitive to the cosmic variance from the CMB. This makes sense because the CMB is completely projected 108 out. The estimator does not care about the value of any CMB multipole, so the cosmic variance of the multipoles is also immaterial. This independence from cosmic variance may be shown algebraically. The estimate appears to depend on the cosmic variance contribution to E (equation D.8) through F (equation D.15). However, if we explicitly write the cosmic variance Σ C = h(C − hCi)(C − hCi)† i contribution to the variance, we see that it projects out of the estimate. Let us define Σ0 as the balance of the variance, the part not due to cosmic variance of the CMB. Write the covariance: Σ = wΣCw † + Σ0 , (D.22) then substitute into E, E ≡ −1 wΣCw † + Σ0 −1 − wΣCw † + Σ0 w −1 −1 † C † 0 w w wΣ w + Σ −1 w † wΣC w † + Σ0 . (D.23) −1 −1 † 0 −1 . w w Σ w w † Σ0 (D.24) If we were to expand each of the matrix inverses in geometric series, we would find that E does not depend on ΣC at all: E = Σ0 −1 − Σ0 −1 We see that the cosmic variance of the CMB has been projected out of the estimator. Thus it cannot impact the estimate B̄, or its variance ΣB . References Abazajian K. et al., 2003, AJ, 126, 2081 Abel T., Wandelt B. D., 2002, MNRAS, 330, L53 Abell G. O., Corwin H. G., Olowin R. P., 1989, ApJS, 70, 1 Afshordi N., Loh Y., Strauss M. A., 2003 Barger A. J., Cowie L. L., Sanders D. B., 1999, ApJ, 518, L5 Barger A. J., Cowie L. L., Sanders D. B., Fulton E., Taniguchi Y., Sato Y., Kawara K., Okuda H., 1998, Nature, 394, 248 Barkana R., Loeb A., 2001, Phys. Rep., 349, 125 Barkana R., Loeb A., 2004, ApJ, 609, 474 Bennett C. L. et al., 2003a, ApJS, 148, 1 Bennett C. L. et al., 2003b, ApJS, 148, 97 Birkinshaw M., 1999, Phys. Rep., 310, 97 Bond J. R., Jaffe A. H., Knox L., 1998, Phys. Rev. D, 57, 2117 Borys C., Chapman S., Halpern M., Scott D., 2003, MNRAS, 344, 385 Bullock J. S., Kolatt T. S., Sigad Y., Somerville R. S., Kravtsov A. V., Klypin A. A., Primack J. R., Dekel A., 2001, MNRAS, 321, 559 109 110 Calabretta M. R., Greisen E. W., 2002, A&A, 395, 1077 Carlstrom J. E., Holder G. P., Reese E. D., 2002, ARA&A, 40, 643 Cen R., 1992, ApJS, 78, 341 Cen R., 2002, ApJS, 141, 211 Cen R., 2003, ApJ, 591, 12 Chiu W. A., Fan X., Ostriker J. P., 2003, ApJ, 599, 759 Chiu W. A., Ostriker J. P., 2000, ApJ, 534, 507 Cooray A., 2000, Phys. Rev. D, 62, 103506 Cooray A., Hu W., Tegmark M., 2000, ApJ, 540, 1 Copi C. J., Huterer D., Starkman G. D., 2004, Physical Review D, 70, 043515 da Silva A. C., Barbosa D., Liddle A. R., Thomas P. A., 2000, MNRAS, 317, 37 da Silva A. C., Barbosa D., Liddle A. R., Thomas P. A., 2001, MNRAS, 326, 155 de Oliveira-Costa A., Tegmark M., Zaldarriaga M., Hamilton A., 2004, Phys. Rev. D, 69, 063516 De Zotti G., Perrotta F., Granato G. L., Silva L., Ricci R., Baccigalupi C., Danese L., Toffolatti L., 2002, ArXiv Astrophysics e-prints Diego J. M., Silk J., Sliwa W., 2003, MNRAS, 346, 940 Eales S., Lilly S., Gear W., Dunne L., Bond J. R., Hammer F., Le F èvre O., Crampton D., 1999, ApJ, 515, 518 Ebeling H., Voges W., Bohringer H., Edge A. C., Huchra J. P., Briel U. G., 1996, MNRAS, 281, 799 Efstathiou G., 2003, MNRAS, 346, L26 111 Eke V. R., Navarro J. F., Steinmetz M., 2001, ApJ, 554, 114 Fan X. et al., 2004, AJ, 128, 515 Finkbeiner D. P., Davis M., Schlegel D. J., 1999, ApJ, 524, 867 Fosalba P., Gaztañaga E., 2004, MNRAS, 350, L37 Fosalba P., Gaztañaga E., Castander F. J., 2003, ApJ, 597, L89 Furlanetto S. R., Zaldarriaga M., Hernquist L., 2004, ApJ, 613, 1 Gao L., White S. D. M., Jenkins A., Frenk C., Springel V., 2005, ArXiv Astrophysics e-prints Gunn J. E., Peterson B. A., 1965, ApJ, 142, 1633 Hamilton A. J. S., 1997, MNRAS, 289, 285 Hernández-Monteagudo C., Rubiño-Martı́n J. A., 2004, MNRAS, 347, 403 Hinshaw G. et al., 2003, ApJS, 148, 135 Hirata C. M., Seljak U., 2003, Phys. Rev. D, 67, 43001 Holder G. P., 2002, ApJ, 580, 36 Holland W. S. et al., 1999, MNRAS, 303, 659 Huffenberger K. M., Seljak U., 2005, New Astronomy, 10, 491 Huffenberger K. M., Seljak U., Makarov A., 2004, Phys. Rev. D, 70, 063002 Hughes D. H. et al., 1998, Nature, 394, 241 Jain B., Seljak U., White S., 2000, ApJ, 530, 547 Jarrett T. H., Chester T., Cutri R., Schneider S., Skrutskie M., Huchra J. P., 2000, AJ, 119, 2498 112 Kaplinghat M., Knox L., Song Y.-S., 2003, Physical Review Letters, 91, 241301 Knox L., Holder G. P., Church S. E., 2004, ApJ, 612, 96 Kogut A. et al., 2003, ApJS, 148, 161 Komatsu E., Seljak U., 2001, MNRAS, 327, 1353 Komatsu E., Seljak U., 2002, MNRAS, 336, 1256 Kosowsky A., 2003, New Astronomy Review, 47, 939 Lacey C., Cole S., 1993, MNRAS, 262, 627 Loeb A., Barkana R., 2001, ARA&A, 39, 19 Maddox S. J., Efstathiou G., Sutherland W. J., Loveday J., 1990, MNRAS, 243, 692 Maselli A., Ferrara A., Ciardi B., 2003, MNRAS, 345, 379 Mason B. S. et al., 2003, ApJ, 591, 540 McQuinn M., Furlanetto S. R., Hernquist L., Zahn O., Zaldarriaga M., 2005, ApJ, 630, 643 Mellema G., Iliev I. T., Alvarez M. A., Shapiro P. R., 2005, ArXiv Astrophysics e-prints Miralda-Escudé J., Haehnelt M., Rees M. J., 2000, ApJ, 530, 1 Mo H. J., White S. D. M., 1996, MNRAS, 282, 347 Moore B., Governato F., Quinn T., Stadel J., Lake G., 1998, ApJ, 499, L5 Myers A. D., Shanks T., Outram P. J., Frith W. J., Wolfendale A. W., 2004, MNRAS, 347, L67 Navarro J. F., Frenk C. S., White S. D. M., 1997, ApJ, 490, 493 Olive K. A., Steigman G., Walker T. P., 2000, Phys. Rep., 333, 389 113 Page L. et al., 2003, ApJS, 148, 233 Pen U., 1998, ApJS, 115, 19 Perlmutter S. et al., 1999, ApJ, 517, 565 Peterson J. B., Pen U. L., Wu X. P., 2004, American Astronomical Society Meeting Abstracts, 205, Press W. H., Flannery B. P., Teukolsky S. A., 1986, Numerical recipes. The art of scientific computing. Cambridge: University Press, 1986 Press W. H., Teukolsky S. A., Vetterling W. T., Flannery B. P., 1992, Numerical recipes in C. The art of scientific computing. Cambridge: University Press, —c1992, 2nd ed. Readhead A. C. S. et al., 2004, ArXiv e-print astro-ph/0402359 Refregier A., Teyssier R., 2002, Phys. Rev. D, 66, 043002 Riess A. G. et al., 1998, AJ, 116, 1009 Rijkhorst E., Plewa T., Dubey A., Mellema G., 2005, ArXiv Astrophysics e-prints Ruhl J. et al., 2004, in Astronomical Structures and Mechanisms Technology. Edited by Antebi, Joseph; Lemke, Dietrich. Proceedings of the SPIE, Volume 5498, pp. 11-29 (2004). Rybicki G. B., Press W. H., 1992, ApJ, 398, 169 Santos M. G., Cooray A., Haiman Z., Knox L., Ma C., 2003, ApJ, 598, 756 Sasaki S., 1994, PASJ, 46, 427 Schlegel D. J., Finkbeiner D. P., Davis M., 1998, ApJ, 500, 525 Seljak U., 1998, ApJ, 503, 492 Seljak U., Burwell J., Pen U., 2001, Phys. Rev. D, 63, 63001 114 Seljak U. et al., 2005, Phys. Rev. D, 71, 103515 Seljak U., Zaldarriaga M., 1996, ApJ, 469, 437 Sheth R. K., Tormen G., 2002, MNRAS, 329, 61 Shimon M., Rephaeli Y., 2004, New Astronomy, 9, 69 Smail I., Ivison R. J., Blain A. W., 1997, ApJ, 490, L5 Snowden S. L. et al., 1997, ApJ, 485, 125 Sofue Y., Rubin V., 2001, ARA&A, 39, 137 Song Y., Cooray A., Knox L., Zaldarriaga M., 2003, ApJ, 590, 664 Spergel D. N. et al., 2003, ApJS, 148, 175 Springel V., White M., Hernquist L., 2001, ApJ, 549, 681 Tegmark M., 1997, Phys. Rev. D, 55, 5895 Tegmark M., de Oliveira-Costa A., Hamilton A. J., 2003, Phys. Rev. D, 68, 123523 Tegmark M., Eisenstein D. J., Hu W., de Oliveira-Costa A., 2000, ApJ, 530, 133 Toffolatti L., Argueso Gomez F., de Zotti G., Mazzei P., Franceschini A., Danese L., Burigana C., 1998, MNRAS, 297, 117 Totani T., Takeuchi T. T., 2002, ApJ, 570, 470 Trushkin S. A., 2003, Bull. Special Astrophys. Obs., 55, 90 White M., Hernquist L., Springel V., 2002, ApJ, 579, 16 White M., Majumdar S., 2004, ApJ, 602, 565 Wyithe J. S. B., Loeb A., 2003, ApJ, 588, L69 Wyithe J. S. B., Loeb A., 2004, Nature, 432, 194 115 Zahn O., Zaldarriaga M., Hernquist L., McQuinn M., 2005, ArXiv Astrophysics e-prints Zaldarriaga M., Seljak U., 1999, Phys. Rev. D, 59, 123507 Zhang P., Pen U., Trac H., 2004, MNRAS, 347, 1224 Zhang P., Pen U., Wang B., 2002, ApJ, 577, 555 Zhou W., Wu X., 2004, ApJ, 600, 501

1/--страниц