вход по аккаунту


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
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.
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.
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.
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
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
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.
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 Introduction
2 Sunyaev-Zel’dovich effect in WMAP and its effect on cosmological parameters
2.1 Original abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Cross-spectra combination for SZ . . . . . . . . . . . . . . . . . . . . . . .
2.3.1 Estimator and data . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.4 Conclusions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3 ACT simulations and power spectrum analysis
3.1 Original Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 Simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.1 Hydrodynamical simulations . . . . . . . . . . . . . . . . . . . . .
3.3.2 Lensed CMB . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.3 Sunyaev-Zel’dovich effect . . . . . . . . . . . . . . . . . . . . . . . .
3.3.4 Kinetic Sunyaev-Zel’dovich effect . . . . . . . . . . . . . . . . . . .
3.3.5 Point sources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.6 Detector noise . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3.7 Caveats and other signals not included . . . . . . . . . . . . . . .
3.3.8 Summary of simulations . . . . . . . . . . . . . . . . . . . . . . . .
3.4 Power spectrum analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.4.1 Multifrequency estimator . . . . . . . . . . . . . . . . . . . . . . .
3.4.2 Relation to Wiener filter . . . . . . . . . . . . . . . . . . . . . . . .
3.4.3 Frequency dependence and scatter . . . . . . . . . . . . . . . . . .
3.4.4 Measured signal power spectrum . . . . . . . . . . . . . . . . . . .
3.4.5 Measuring primordial amplitude and slope . . . . . . . . . . . . .
3.5 One-point distribution function analysis . . . . . . . . . . . . . . . . . . .
3.5.1 Signal one-point distribution functions . . . . . . . . . . . . . . . .
3.5.2 Likelihood . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.5.3 Separation by distribution function . . . . . . . . . . . . . . . . . .
3.6 Conclusions
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4 Progress towards simulation of reionization with radiation transport 71
4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2 Source model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.1 Halo population model . . . . . . . . . . . . . . . . . . . . . . . . .
4.2.2 Source luminosity . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.3 Intergalactic Medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4 Radiation transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.1 Cen’s Flux Method . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.4.2 Photon Absorption . . . . . . . . . . . . . . . . . . . . . . . . . . .
4.5 Tests and future directions . . . . . . . . . . . . . . . . . . . . . . . . . . .
A Semi-analytic simulations of the Sunyaev-Zel’dovich effect
A.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.2 Analytic halo profile . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.3 Halo catalog . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
A.4 Semi-analytic Sunyaev-Zel’dovich maps . . . . . . . . . . . . . . . . . . .
B Importing flat-sky simulations into the ACT geometry
C Signal filtering in ACT maps
C.1 Linear (Wiener) filter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
C.2 An example: point source identification . . . . . . . . . . . . . . . . . . . 101
D Cross-spectrum contaminant estimator
Chapter 1
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
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
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?”
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
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,
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
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
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–
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,
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
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
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
ν (GHz)
∆T SZ /(T CMB y)
Cl /ClRJ
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
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.
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
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
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.
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 ,
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
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 =
wlli 0 ClCMB
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
Sl =
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.
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
∝ −2f (x),
where the frequency dependence relative to the RJ regime is given by
f (x) = 2 −
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
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
L ∝ exp − [D − hDi] Σ [D − hDi]) ,
where the covariance Σ = h(D − hDi)(D − hDi) † i can be written as Σ = {Σii
ll0 }. We
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 ,
where we have defined the auxiliary matrices
F ≡ E − ES S ES
S† E
w † Σ−1
E ≡ Σ−1 − Σ−1 w w † Σ−1 w
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 power spec.cfm
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.
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
l(l+1) Cl / 2π (µK )
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.
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.
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
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.
l(l+1) Cl / 2π (µK )
in 9
ov C
a 95
% li
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 × 102
no SZ
−0.13 −0.26
+0.07 +0.15
0.245−0.06 −0.10
−0.015 −0.029
+0.07 +0.10
0.19−0.08 −0.14
−0.11 −0.20
+0.06 +0.12
0.74−0.05 −0.09
−0.04 −0.07
with SZ, no SZ prior
< 7.1 (95%)
−0.18 −0.33
+0.07 +0.15
0.234−0.06 −0.10
−0.016 −0.030
+0.07 +0.10
0.20−0.08 −0.14
−0.11 −0.22
+0.07 +0.14
0.76−0.06 −0.11
−0.04 −0.07
with SZ prior
< 2.9 (95%)
−0.15 −0.29
+0.07 +0.15
0.243−0.06 −0.10
−0.016 −0.031
+0.07 +0.10
0.19−0.08 −0.14
−0.11 −0.22
+0.06 +0.13
0.75−0.05 −0.10
−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.
Point source amplitude A (µ K )
Best fit
WMAP est.
unresolved src.
SZ amplitude B
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.
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
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
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
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
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
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,
Band (GHz)
Beam FWHM (arcmin.)
Sensitivity per pixel (µK)
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
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,
ΩΛ = 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̂) =
d2 l Θl exp(in̂ · l).
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.
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
kSZ x 10
Rad. pt. src. x 10
IR pt. src.
(100 µΚ)
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).
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.
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̂) =
r(χ)r(χLS − χ) δ(n)
r(χLS )
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
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
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)
Θ (n̂) = −2y(n̂) 2 −
l(l+1) Cl / 2π
l(l+1) Cl / 2π
IR Sources
Radio Sources
10000 1000
10000 1000
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.
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̂) =
σT ne (n)
kB Te (n)
me c2
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.
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.
(n̂) =
σT Xe (n)ne (n)
v(n) · n̂
exp[τ (χ)]adχ.
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).
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
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
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
S dS
Cl = 2
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.
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,
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
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
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
The differential counts for radio sources at the level of our flux cut are shallower,
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.
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
scales. Assuming the noise and Gaussian beams in Table 3.1, the power spectrum of
noise is given by:
= (σα bα ) exp
(lb )2
√ α
8 log 2
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.
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
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.
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.
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
145 GHz
217 GHz
265 GHz
(100 µΚ)
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.
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.
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.
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
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
C = hee† i = Cαβ
(Z) + Nαβ
= Sαβ
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
L(e) = (2π)−N/2 det(C)−1/2 exp(− e† C−1 e).
A second order expansion in Zp gives
ln L(Z + δZ) = ln L(Z) +
X ∂ ln L(Z)
1 X ∂ 2 ln L(Z)
δZp δZp0 ,
2 0 ∂Zp ∂Zp0
∂ ln L(Z)
∂ ln L(Z)
∂Zp ∂Zp0
∂C −1
C e − C−1
∂C −1 ∂C −1
= e† C−1
C e
∂C −1 ∂C
− tr (C−1
= tr (e† C−1
The ensemble average of the second expression is the Fisher matrix
∂C −1 ∂C
Fpp0 = tr (C−1
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
bp = tr [C−1
∂C −1
C N].
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:
Cpp0 = hẐp Ẑp0 i = Fpp
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.
We will also use the (cross) power spectrum of signals, denoted S mn
= hslm sln iδll0 . We
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
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:
Rαβmn Smn
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
m0 n0
Rαβm0 n0 δll0
× [δmm0 δnn0 + δmn0 δnm0 − δmm0 δnn0 δmn0 δnm0 ]
The final expression in the brackets is shorthand for:
 1 if (m = m and n = n ) or (m = n and n =
[. . .] =
m0 ),
 0 otherwise.
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
Rαβmn Amn Smn
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
m0 n0
Rαβm0 n0 Smn
× [δ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.
Relation to Wiener filter
If all signals have a spatially uniform frequency dependence, then s lα =
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
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 =
Smm0 ym0
m C e.
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
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.)
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
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.
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
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
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
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 .
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
realizations χ2 =
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
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
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
For case 2, assuming the shape of the point source power spectrum improved the
l(l+1) ∆Cl / 2π
l(l+1) ∆Cl / 2π
Monte Carlo
Case 1
Case 2
Case 3
Case 4
Case 1
Case 2
Case 3
Case 4
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.
l(l+1) Cl / 2π
l(l+1) Cl / 2π
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.
l(l+1) Cl / 2π
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
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
l(l+1) Cl / 2π
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-
l(l+1) Cl / 2π
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.
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
qq 0 =
X dCl
dq 0
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
C(l < 2000) =
A 2 × 10−2 −4 × 10−3
 n
1 × 10−3
However, if we can extract the CMB, the covariance improves.
A 4 × 10
−9 × 10
2 × 10−4 
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
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
Total @ 217 GHz
IR Sources
Radio Sources
l(l+1) Cl / 2π
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.
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
IR sources
-3 -2 -1 0
∆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.
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) =
Pa (∆Ti ),
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
over bins:
log L(a) =
Nj log Pa (∆Tj ),
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.
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
. 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
+ σ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
σkSZx10 (K)
σCMBx10 (K)
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 .
+ σ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
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
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
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
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
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
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.
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,
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
(δ1 − δ0 )2
exp −
2(s − S) (s − S)
2π(s − S)
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 =
f (s, δ1 |S, δ0 )
casts the conditional mass function in more familiar form, as a number per unit
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.68647 −
(1 + δ)2/3
(1 + δ)0.58661
(1 + δ)1/2
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
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 =
[nt+∆t (m) − nt (m)] dm
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 )
Thus the mean number of halos destroyed (that is, incorporated) during a time step is
h∆Ndest i = hN i × 1 −
D(t + ∆t)
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.
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
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.
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)
`(m, t, tf ) = ε
t − tf
t − tf
exp −
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) =
∆Nform (mi , tf ) `(mi , t, tf ) .
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
10 × L/εVc (solar masses/second/Mpc )
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
100 h-1 Mpc box, 323 resolution
z = 15
few sources
Fraction of light
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
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
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
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:
= P (x) − R(x).
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
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
−0.5 T
× 1+
106 K
103 K
0.7 !−1
cm3 s−1
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
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.
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) =
I(~x, n̂0 , ν, t) dn̂0 ,
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
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
One objection to this method is that at large distance shadows are cast imperfectly,
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
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
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
= −κI + j,
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/κ,
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) =
= (I0 − j/κ)[1 − exp(−κs)] + js.
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 κ.
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:
A∆Ω =
dΩ a(I0 (~x), s(~x)),
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,
3 Z
A∆Ω =
dΩ a (I0 (~x), s(~x)) |êi · n̂| ,
face i
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
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∆Ωk .
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
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
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
Θ (n̂) = −2y(n̂) 2 −
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̂) =
σT ne (n)
kB Te (n)
me c2
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
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)
Tgas (r) = Tgas (0) ygas
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-
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 )3+α
)−α (1
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 ≡
ln(1 + x)
1−B× 1−
1 −1
3 γ − 1 ln(1 + c)
η(0) γ
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
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
The temperature normalization is given by
Tgas (0) = 8.80 keV × η(0)
M/1015 h−1 M
rvir /1 h−1 Mpc
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.
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:
c(M, z) =
M∗ (z = 0)
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
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
Celestial Sphere
δ = −55 deg
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
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
Each redshift box is then tiled onto the cone surface with a random offset and projection direction.
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
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
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 ,
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].
The covariance of the estimate is also easy to write down:
hs̃s̃† i = WCW† .
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 ,
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)
s (k) + Pn (k)
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
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
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
for y, by minimizing
f (y) =
1 †
y Cy − d† y,
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.
145 GHz
217 GHz
265 GHz
145 GHz
217 GHz
265 GHz
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.
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
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:
−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
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
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.
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)] .
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,
−2 log L + const. = C† aC + (b† C + C† b) + c† Σ−1 c,
where we have introduced the shorthand
a ≡ w † Σ−1 w
b ≡ −2w † Σ−1 (D − SA − ZB)
c ≡ D − SA − ZB.
We complete the square for variable C:
1 −1 †
1 −1
−2 log L + const. = C + a b a C + a b
+c† Σ−1 c − b† a−1 b
Note that the final term may be rewritten in terms of our auxiliary variable c,
1 † −1
† −1
† −1
b a b=c Σ w w Σ w
w Σ
Thus if we define the symmetric matrix
E ≡ Σ−1 − Σ−1 w w † Σ−1 w
w † Σ−1 ,
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 †
−2 log L + const. = C + a b a C + a b + c† Ec.
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,
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],
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
We wish to repeat this sequence to marginalize the variable A. Thus we write
−2 log LC + const. = A† fA + (g† A + A† g) + h† Eh
where we have introduced
≡ S† ES
g ≡ −2S† E(D − ZB)
h ≡ D − ZB.
Again we complete the square, now for variable A:
−2 log LC + const. =
† A + f −1 g f A + f −1 g
+h† Eh − g† f −1 g.
Noting the last term may be re-written in terms of h,
1 † −1
g f g = h ES S ES
S E h,
F ≡ E − ES S† ES
S† E.
we write with analogy to equation (D.8),
Define LC,A ≡
dCdA L =
dA LC , and as before we integrate out our contaminants:
−2 log LC,A + const. = h† Fh
= [D − ZB]† F [D − ZB].
To estimate the amplitudes B, we express the marginalized likelihood as
−2 log LC,A + const. = B† pB + (q† B + B† q) + r
where we have introduced
p ≡ Z† FZ
q ≡ −2Z† FD
r ≡ D† FD.
We complete the square one final time for variable B:
−2 log LC,A + const. = [B − p−1 q]† p[B − p−1 q]
1 † −1
+r − q p q.
Now LC,A (B) ∝ exp[− 21 (B − hBi)† (ΣB )−1 (B − hBi)]. If we estimate
B̄ ≡ (Z† FZ)−1 Z† FD,
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 .
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
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 ,
then substitute into E,
E ≡
wΣCw † + Σ0
− wΣCw † + Σ0
−1 −1
C †
w wΣ w + Σ
w † wΣC w † + Σ0
−1 −1 † 0 −1
w Σ
w w † Σ0
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
− Σ0
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 .
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
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
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
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
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,
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
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
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
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
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 051 Кб
Пожаловаться на содержимое документа