Open Collections

UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Finding constraints on dark energy using hydrogen intensity mapping Fink, Allison 2017

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


24-ubc_2018_february_fink_allison.pdf [ 2.11MB ]
JSON: 24-1.0362240.json
JSON-LD: 24-1.0362240-ld.json
RDF/XML (Pretty): 24-1.0362240-rdf.xml
RDF/JSON: 24-1.0362240-rdf.json
Turtle: 24-1.0362240-turtle.txt
N-Triples: 24-1.0362240-rdf-ntriples.txt
Original Record: 24-1.0362240-source.json
Full Text

Full Text

Finding Constraints on Dark Energy Using HydrogenIntensity MappingbyAllison FinkB.A., Bryn Mawr College, 2011A THESIS SUBMITTED IN PARTIAL FULFILLMENTOF THE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinTHE FACULTY OF GRADUATE AND POSTDOCTORALSTUDIES(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)December 2017c© Allison Fink, 2017AbstractIn this MS.c. thesis, we demonstrate a method for estimating the expansionhistory of the universe using the hydrogen intensity map from the Canadian Hy-drogen Intensity Mapping Experiment (CHIME), which will be generated in thenear future. This map will be in angular and redshift space, where redshift of thehydrogen due to the expansion serves as a time variable. The expansion history,dependent on cosmological parameters via the Einstein equations, determines thedistance away from us in units of a grid comoving with the expansion at which lightof each redshift was emitted. We use knowledge of the fixed comoving distance,approximately 150 megaparsecs, that baryon acoustic oscillations, or primordialsound waves, traveled away from the centers of matter perturbations, where thereis a corresponding peak in the matter correlation function subject to uncertainty ofthe initial quantum mechanical fluctuations. We explain the method by which wefit the correlation function to a model for expansion determined by the equationof state of dark energy, to constrain this parameter. We test the method using athree-dimensional realization of the theoretical matter power spectrum calculatedfrom CAMB (Code for Anisotropies in the Microwave Background), providing anestimate of constraints obtained from a small redshift region spanning one soundhorizon diameter in redshift space assuming a constant equation of state of darkenergy and fixed values of the other cosmological parameters. We explain how togeneralize this method to a more complete analysis.iiLay SummaryIn this MS.c. thesis, we demonstrate a method to estimate the universe’s expan-sion history using a hydrogen intensity map to be made by the Canadian HydrogenIntensity Mapping Experiment (CHIME). This map will be in angular and redshiftspace, where redshift of the hydrogen due to the expansion serves as a time vari-able. The expansion history, dependent on cosmological parameters via the Ein-stein equations, determines the distance away from us in units of a grid comovingwith the expansion at which light of each redshift was emitted. Using knowledgeof the fixed comoving distance that primordial sound waves traveled away from thecenters of matter perturbations, where there is a corresponding peak in the mattercorrelation function subject to uncertainty of the initial quantum mechanical fluc-tuations, we fit the correlation function to a model for expansion determined by theequation of state of dark energy, to constrain this parameter.iiiPrefacePart 2 of this work is based on the work conducted in the Department of Physicsand Astronomy at the University of British Columbia by Dr. Gary Hinshaw, Dr.Kris Sigurdson, and Allison Fink. I was responsible for designing a method formeasuring the expansion history of the universe using simulated data. This in-volved writing computer codes in PYTHON and testing the methods.ivTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viiiAcknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Motivation and Brief Overview of Research Goals . . . . . . . . . 11.2 The Question of Dark Energy . . . . . . . . . . . . . . . . . . . . 21.3 Evidence for the Expansion of the Universe . . . . . . . . . . . . 41.3.1 The Cosmic Microwave Background . . . . . . . . . . . . 41.3.2 Big Bang Nucleosynthesis . . . . . . . . . . . . . . . . . 51.3.3 The Hubble Curve . . . . . . . . . . . . . . . . . . . . . 61.4 Probes of Dark Energy . . . . . . . . . . . . . . . . . . . . . . . 71.4.1 Type IA Supernovae . . . . . . . . . . . . . . . . . . . . 71.4.2 Weak Lensing . . . . . . . . . . . . . . . . . . . . . . . . 81.4.3 Galaxy Cluster Abundances . . . . . . . . . . . . . . . . 91.4.4 Baryon Acoustic Oscillations . . . . . . . . . . . . . . . 111.5 The Friedmann-LeMaıˆtre-Robertson-Walker Metric . . . . . . . . 15v1.6 Thermal History of the Universe . . . . . . . . . . . . . . . . . . 171.6.1 Statistics and Entropy in the Early Universe . . . . . . . . 181.6.2 The Boltzmann Equation . . . . . . . . . . . . . . . . . . 211.6.3 The Neutrino and Photon Energy Densities . . . . . . . . 221.6.4 Big Bang Nucleosynthesis . . . . . . . . . . . . . . . . . 241.6.5 The Free Electron Fraction From Recombination . . . . . 281.7 Cosmological Perturbation Theory . . . . . . . . . . . . . . . . . 351.7.1 Degrees of Freedom for Perturbations . . . . . . . . . . . 351.7.2 Choosing a Gauge for Scalar Perturbations . . . . . . . . 391.7.3 The Boltzmann Equation for Photons . . . . . . . . . . . 421.7.4 Perturbed Einstein Equations and Conservation Equations 471.8 Solutions for the Matter Power Spectrum . . . . . . . . . . . . . . 501.8.1 Initial Conditions for the Boltzmann Equations . . . . . . 511.8.2 Large Scale Solutions . . . . . . . . . . . . . . . . . . . 541.8.3 Cold Dark Matter Solution . . . . . . . . . . . . . . . . . 581.8.4 Small Scale Solutions . . . . . . . . . . . . . . . . . . . 591.8.5 The Power Spectrum and Transfer Function . . . . . . . . 661.8.6 Initial Conditions for the Potential . . . . . . . . . . . . . 67Use of Inflation to Solve the Horizon Problem . . . . . . . 67Variance of the Field Perturbations . . . . . . . . . . . . . 70The Power Spectrum Following Inflation . . . . . . . . . 721.9 The Hydrogen Brightness Temperature Power Spectrum . . . . . 741.10 Mapmaking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 762 Implementation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 812.1 Determining Cosmological Parameters From Maps . . . . . . . . 812.2 The Flat Sky Approximation . . . . . . . . . . . . . . . . . . . . 832.3 Calculation of the Matter Power Spectrum . . . . . . . . . . . . . 862.4 Calculation of the Clνν ′ . . . . . . . . . . . . . . . . . . . . . . . 872.5 Window Function . . . . . . . . . . . . . . . . . . . . . . . . . . 882.6 Generation of Sample Realizations . . . . . . . . . . . . . . . . . 902.7 Likelihood Analysis . . . . . . . . . . . . . . . . . . . . . . . . . 932.8 Calculating the Likelihood . . . . . . . . . . . . . . . . . . . . . 97vi2.9 The Covariance Matrix . . . . . . . . . . . . . . . . . . . . . . . 972.9.1 Model for the Covariance Matrix . . . . . . . . . . . . . . 982.9.2 Inverting the Covariance Matrix . . . . . . . . . . . . . . 992.10 Testing the Method for Bias . . . . . . . . . . . . . . . . . . . . 1002.11 Finding Best Estimate Constraints on w . . . . . . . . . . . . . . 1012.12 Conclusion and Direction for Future Work . . . . . . . . . . . . . 101Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104A Gauge Freedom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111B Spherical Harmonics . . . . . . . . . . . . . . . . . . . . . . . . . . 119viiList of FiguresFigure 2.1 Theoretical Correlation Function, 767 MHz (z=.851) . . . . . 90Figure 2.2 Theoretical Correlation Function, 767 and 780 MHz (z=.851and .821) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91Figure 2.3 Theoretical Radial Correlation Function, 767 MHz (z=.851) . 91Figure 2.4 Simulated 21 cm Brightness Temperature Anisotropies, 767MHz (z=0.851) . . . . . . . . . . . . . . . . . . . . . . . . . 92Figure 2.5 Simulated 21 cm Brightness Temperature Anisotropies, 780MHz (z=0.821) . . . . . . . . . . . . . . . . . . . . . . . . . 92Figure 2.6 Simulated 21 cm Brightness Temperature Anisotropies, 790MHz (z=0.797) . . . . . . . . . . . . . . . . . . . . . . . . . 93Figure 2.7 Simulated Angular Correlation Function, 767 MHz (z=.851) . 94Figure 2.8 Simulated Radial Correlation Function, 767 MHz (z=.851) . . 94Figure 2.9 Simulated Correlation Function, 767 and 780 MHz (z=.851and .821) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95Figure 2.10 Likelihood of Model Given Theoretical Correlation Function . 100Figure 2.11 Likelihoods of Model Given Respective Physical Realizationsof the Correlation Function . . . . . . . . . . . . . . . . . . . 102Figure 2.12 Likelihood of Standard Deviation in w for 10 Sample Realiza-tions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102viiiAcknowledgmentsI would like to thank my MS.c. supervisors, Gary Hinshaw and Kris Sigurdson,for their mentorship and guidance. I am extremely grateful for all I have learnedthrough working with them. Also, I would like to thank my parents, friends, andpeers for their support, understanding, and help in many ways.ixChapter 1Introduction1.1 Motivation and Brief Overview of Research GoalsAn understanding of dark energy, which makes up approximately 75 percentof the energy in the universe, is an important unresolved question in physics.The Canadian Hydrogen Intensity Mapping Experiment (CHIME) telescope, underconstruction near Penticton, B.C. at the time of writing, will take the measurementsneeded to estimate values of cosmological parameters, such as the equation of stateof dark energy, accurately enough to settle this question. As we will explain later,the values of cosmological parameters can be estimated and constrained from thevariance of the observed fluctuations in the matter density field. As originallyproposed in ([61]), at the time that dark energy is noticeable, the intensity of thecharacteristic 21 cm wavelength photon emitted by neutral hydrogen is an effectivetracer of this matter distribution. Due to the universe’s continuous expansion, thiswavelength gets increasing redshifted with the time since emission. The CHIMEtelescope will map the intensity of the neutral hydrogen signal over the sky as afunction of redshift, creating a three-dimensional map with redshift serving as aparametrization of the time variable. We can then fit the data from this map to acosmological model, thereby obtaining constraints on the cosmological parametersdetermining the expansion history.As explained in ([67]), ([41]), CHIME will measure the intensity of the neutralhydrogen signal over half the sky and at a range of frequencies f = 400-800 MHz,1or redshifts z=0.7 to z=2.5, which covers a significant range of time during whichdark energy is noticeable in the matter power spectrum. The map will be in 0.39MHz intervals and have an angular resolution of 0.52− 0.26◦. In this project, weoutline a specific method for finding constraints on the expansion history usingthe visible full three-dimensional shape of the matter correlation function obtainedfrom a map taken at a similar resolution as that of CHIME, and show examples ofconstraints obtained from simulated data. In our analysis of the matter correlationfunction, we use the method of baryon acoustic oscillations as a standard ruler, arobust probe of the expansion history ([36]), ([22]), ([47]), ([78]), see also ([37]),which we will explain shortly.In the next sections of the introduction, we first explain the evidence for theexpansion of the universe, followed by an overview of the question of dark energyand the methods that have been used so far to probe it. This includes more de-tails about the method of using baryon acoustic oscillations as a standard ruler andhydrogen intensity mapping. We next explain in detail the background physics rel-evant for the calculation of the matter power spectrum which is used to determinethe expansion history through the signal of baryon acoustic oscillations. Muchof the discussion and topics covered in the introduction are based on those givenin ([34]), as well as ([45]), ([50]), and related work by the same author cited inthe text, and supplemented by more explanations from other references as needed,such as ([86]) and ([9]). Note that in all references to ([34]) and ([86]) given in thetext, the reader should keep in mind the list of corrections to these sources givenin respectively ([6]) and ([85]). In the main body of the thesis, we explain ourmethod for measuring the expansion history, how we have tested the method, andour conclusions about the results. We end by summarizing steps for the future forobtaining a more complete analysis of the expansion history.1.2 The Question of Dark EnergyThere are many experimental sources of evidence that the universe’s expansionis accelerating, which we will discuss later. The source of this is presently not fullydetermined and so it is simply called dark energy. There are a few possibilities open(see ([29])). It is possible that dark energy can be accounted for by introducing an2unknown substance into the Einstein equations. The dark energy contribution tothe stress energy tensor in general relativity, assuming that the dark energy is onaverage homogeneous and isotropic, is characterized by a pressure and a density(see ([82])). The behavior of the density of dark energy can then be parametrizedin terms of an equation of state, w, which will be useful as follows ([34]). Asexplained in later sections, assuming that dark energy does not interact with othercomponents, it obeys a separate continuity equation, found as a consequence of theEinstein equations:∂ρDE∂ t+a˙a[3ρDE +3PDE ] = 0. (1.1)Integrating this from redshift 0 to a, we haveρDE α e−3∫ a0daa (1+w(a)), (1.2)wherew≡ PDEρDE. (1.3)If w(a) is constant, this isρDE α a−3(1+w). (1.4)Currently, there is room for freedom in this equation of state parameter. Using acombination of methods, the WMAP analysis ([56]) found, assuming a constantequation of state,w =−1.10±0.14 (1.5)to 68% confidence level. If w = −1, dark energy corresponds to a cosmologicalconstant Λ in the Einstein equations, Gµν+Λgµν = 8piG3 Tµν . When we put the cos-mological constant term on the right-hand side, we can interpret it as an additionalcontribution to the energy density which we call ρΛ = Λc48piG . This can be conve-niently written in terms of the critical density ρcr =3H208piGc2 , where H0 is the presentvalue of the expansion rate, as ρΛ =ΩΛρcr. The WMAP analysis is consistent withρΛ found from ΩΛ ≈ 0.725±0.016 and H0 = 70.2±1.4 kms Mpc . Recently, ([83]) hassuggested that the value of this cosmological constant contribution to the energydensity can be satisfactorily explained as being due to that of the vacuum energy ofquantum field theory, showing that this vacuum contains large fluctuations which3largely cancel each other, so that a small net energy contribution to the expansionis achieved. However, there is not yet enough experimental information to rule outother possible contributions, such as another substance with a different equation ofstate parameter or a breakdown of Einstein’s equations.1.3 Evidence for the Expansion of the UniverseThe universe is believed by most physicists today to have begun with a space-time singularity called the Big Bang and to have been expanding ever since. Herewe explain the three pieces of evidence for this, as given in ([34]).1.3.1 The Cosmic Microwave BackgroundThe first piece of evidence put in favor of the big bang theory, as explainedin ([34]), ([86]), see also ([9]), ([82]), is the observation of the cosmic microwavebackground, which is the radiation predicted to be left over from the big bang.One possible assumption for the big bang initial conditions is that at early times,matter and radiation were in thermal equilibrium, which means that they were at acommon temperature. This implies that photons followed Bose-Einstein distribu-tions. As will be explained later, given a knowledge of the number of relativisticdegrees of freedom at a given temperature of the universe, from entropy conserva-tion, which holds in equilibrium, one can then find the scale factor as a functionof the temperature. From this, we conclude that there was a time in the early uni-verse when temperatures were high enough such that there were no neutral atoms,only free electrons, and the scattering rate for electrons interacting with photonsby Compton scattering was in much higher than the expansion rate, so in equilib-rium. This persisted until a temperature when the universe’s expansion rate be-came greater than the scattering rate and the photons decoupled from the electrons,allowing the observation of a cosmic microwave background. The thermal equilib-rium condition for the photons and matter up to the time of last scattering impliesthat the photons would have a blackbody spectrum, with intensity as a function offrequency given byI(ν) =4pi h¯ν3/c2e2pi h¯ν/kBT −1 . (1.6)4This prediction of the big bang under equilibrium initial conditions has been ex-perimentally verified. The cosmic microwave background was observed in the1960s by Penzias and Wilson, who found that measurements of one frequency andtemperature agreed with a blackbody spectrum ([71]). The full intensity versusfrequency graph was measured by the COBE satellite with results published in1994 ([64]). It was found to be an almost perfect blackbody with a temperature ofapproximately T = 3K.1.3.2 Big Bang NucleosynthesisAs explained in ([34]), after the universe cooled below the binding energy ofnuclei, it was possible for nucleons to fuse together. By solving the Boltzmannequations in an expanding universe, one can find that the main elements that formedwere 2H, also known as deuterium D, and 4He, with small amounts of 3He and 7Li.One can make definite predictions for the relative abundances of light elementswith respect to hydrogen, as a function of the total baryon density. We explain themethods of measuring these light elements as in ([27]), ([34]) below.The relative abundances of deuterium relative to hydrogen have been con-strained in ([26]) by observing the spectrum of light from distant quasars, whichwill be damped at the corresponding two wavelengths corresponding to the transi-tion from the n = 1 to n = 2 level for hydrogen and deuterium. Due to the smallerabundance of deuterium, the wavelength corresponding to absorption by deuteriumhas a smaller amount of damping, and this can be used to determine their fractionalabundances.Measurements of the other primordial element abundances can be inferred frommeasurements in stars. The fractional abundance of 4He with respect to hydrogen,known as Y, at the time of big bang nucleosynthesis, is not the same as it is todaybecause this abundance increases with time due to processing. However, the de-gree of this processing is indicated by the fractional abundance of oxygen relativeto hydrogen, O/H. Y can be measured as a function of O/H and extrapolated backto the values at O/H = 0, which indicates the abundances at the time of big bangnucleosynthesis, as done in ([52]). 7Li has been measured in the oldest stars inthe halo of our galaxy, where it is believed with support from observations that the5stars with the higher surface temperature preserve their lithium. Lastly, althoughin stars, deuterium is converted to 3H, it is believed, also with support from ob-servations, that the total abundance of D plus 3He in stars is constant, so that the3He abundances can be inferred from the total measurements combined with the Dmeasurements.Based on the measurements of the relative abundances above, and the predic-tions of the relative abundances as a function of baryon density, one can find acommon range of baryon densities in which the theoretical predictions are withinthe range of the error bars of the measurements. This result is considered to beone of the strongest pieces of evidence for the big bang theory since given that theabundances span about nine orders of magnitude, it would otherwise be hard toexplain the satisfaction of the precision required for this agreement.1.3.3 The Hubble CurveThe last piece of evidence in favor of the big bang theory, explained in ([34]),is given by the Hubble curve, which shows that galaxies appear to be moving awayfrom us. Hubble inferred the velocities of galaxies from their redshifts, z, definedbyz =λobs−λemitλemit=1a−1, (1.7)up to approximately z ≈ 0.003. For small redshifts (up to z ≈ .0003, as measuredby Hubble), the Doppler formulaz =√1+(v/c)1− (v/c) −1, (1.8)can be treated non-relativistically and becomesz≈ vc(1.9)and is then used to obtain v. Hubble then measured the redshifts and distances tothe galaxies D and plotted v as a function of D, and found that the graph could befitted to a straight line with a positive slope. This curve could be explained as beingdue to the expansion of the universe. If the universe were expanding, the distances6to galaxies could be written asD = a(t)χ (1.10)where a(t) is the scale factor for the distance (equal to one in the present day), andχ is the comoving length out to the galaxy. Assuming that there is no significantmotion independent of expansion, one can treat χ as constant in time, and thevelocity is given byv =ddt(aχ) =dadtχ = HD, (1.11)where H ≡ da/dta is the Hubble rate of expansion. Assuming the data was over a suf-ficiently small range of redshifts for H to be approximately constant, the straight-line trend could then be explained, and this slope could be used to approximate theexpansion rate in the recent past.1.4 Probes of Dark EnergyWe now explain methods that have been applied to constrain dark energy in thepast, ending with baryon acoustic oscillations, the method used in this work, using([34]) and other sources referenced. More details about each of these methods canbe found in ([84]).1.4.1 Type IA SupernovaeAs explained in ([34]), a type IA supernova serves as a standard candle, orobject with an intrinsic brightness, or luminosity L. The measured flux of the su-pernova is given byF =L4pid2L(1.12)where L is the known luminosity at the source, anddL ≡ χa (1.13)is the luminosity distance. Here,χ(z) =∫ t0tcdt′a(t′) =∫ z0cdz′H(z′) , (1.14)7dependent on the expansion history, is the comoving distance out to redshift z, orin other words the comoving distance away from us at which that the light, nowobserved by us at time t0 to have redshift z, was emitted. The high luminosities ofthe supernovae make it possible to measure χ(z) to relatively high redshifts of z =1.7, from which H(z) can be estimated, extending our knowledge of the expansionhistory past the points found for the original Hubble diagram. Type Ia supernovaeserve as the first experimental confirmation of the existence of the acceleratingexpansion of the universe due to a cosmological constant or dark energy, publishedin ([74]).1.4.2 Weak LensingAs explained in ([34]), according to general relativity, perturbations to the met-ric influence the path of a light ray according to the geodesic equation. We canwrite the geodesic equation in terms of the metric perturbationΦ, the source angle,θ S, and the observed angle θ , for a galaxy observed a distance χ(z) away from us,assuming a flat sky approximation in the limit of small angles, to obtain∂θ Si∂θ j−δi j ≡(−κ− γ1 −γ2−γ2 −κ+ γ)= 2∫ χ0dχ ′Φ,i j(~x(χ ′))(1− χ ′χ). (1.15)In the matrix elements on the left-hand side, κ is the convergence, and γ1 and γ2are the two components of the shear.The ellipticities of galaxies are related to the source angle θ S according toε1 ≡ qxx−qyyqxx+qyy (1.16)ε2 ≡ 2qxyqxx+qyy , (1.17)where qi j are the quadrupole moments of the of the observed intensity, given byqi j =∫d2θ Iobs(~θ)θiθ j =∫d2θ I(~θS)θiθ j. (1.18)8Assuming small angles, we can therefore substituteθi = (A−1)i jθ Sj (1.19)where Ai j ≡ ∂θSi∂θ j with the matrix elements given in the left-hand side of (1.15), in(1.18), and integrate over ~θ S. Substituting the results into (1.16) and (1.17), andassuming small values for the distortions, we obtainε1 ≈ 2γ1 (1.20)ε2 ≈ 2γ2. (1.21)Hence the observations of the ellipticities of galaxies tell us the shear field. We canthen calculate the two-dimensional shear power spectrum, 〈ψi j(~l)ψlm(~l′)〉 whereψi j ≡ Ai j−δi j = 2∫ χ∞0dχW (χ)∫ χ0dχ ′Φ,i j(~x(χ ′))χ ′(1− χ ′χ), (1.22)where the right-hand side follows from (1.15) this time averaging over the distribu-tion of redshifts of the survey using the normalized weighting factor W (χ ′). Takingthe derivatives in l space, where l is the conjugate variable to θ on the sky, we thenobtain for the variance〈ψi j(~l)ψlm(~l′)〉= (2pi)2δ 2(~l−~l′)∫ χ∞0dχg2(χ)χ2lil jlllmχ4PΦ(lχ), (1.23)whereg(χ)≡ 2χ∫ χ∞χdχ ′(1− χχ ′)W (χ ′). (1.24)This will allow us to infer the matter power spectrum as a method to constrain darkenergy.1.4.3 Galaxy Cluster AbundancesAs we will explain below (see ([34]), ([65]), and references therein), one ofthe ways to infer the matter power spectrum and therefore determine constraintson dark energy is by measuring n(M,z), the number density of galaxy clusters as9a function of their mass and redshift. After gravitational collapse has occurred,the average mass enclosed in a sphere of comoving radius R is given by a constantvalueM(R) =Ωmρcr43piR3 = 1.16×1015Ωmh−1M(R10h−1MpC)3, (1.25)where ρcr is explained in (1.52), and ρ¯m/ρcr = Ωm/a3, and h ≡ H0sMpC100km ≈ 0.7.One can show that the matter power spectrum becomes nonlinear at a scale ofR . 10h−1 MpC, and so galaxy clusters, which can have masses up to 1015M,assuming their average densities are on the order of the average density of theuniverse, are enclosed in large enough radii to have their abundances predicted bylinear theory. This can be done using Press-Schechter theory.Press-Schechter theory, introduced in ([72]), assumes that the fraction f (R,z)of the total matter existing in collapsed regions with mass M(R) > 43piR3Ωmρcr isgiven byf (R,z) =2√2piσ(R,z)∫ ∞δcdδRe−δ2R/2σ2(R,z), (1.26)whereδR(~x)≡∫d3x′δ (~x′)WR(~x−~x′) (1.27)is the overdensity field smoothed over a scale of comoving radius R using the tophatwindow functionWR(r) =(43piR3)−1Θ(1− rR), (1.28)and where σ2(R,z) is the variance of δR, which can be written in terms of the linearmatter power spectrum, P(k,z), asσ2(R,z) =∫ ∞0dkk∆2(k,z)W˜ 2R (k), (1.29)where∆2(k,z)≡ (2pi)−34pik3P(k,z), (1.30)andδc ≈ 1.686 (1.31)10is the calculated critical overdensity for collapse.This formula for f (R,z) does not have an initial clear justification, and wasnot fully justified in the original paper. It also uses the overdensities found fromlinear perturbation theory, whereas the mass inside these regions has reached thescale of galaxies which are at approximately 1 MpC, where nonlinearities becomeimportant. However, after comparison with numerical simulations, the fraction ofcollapsed galaxy clusters of a given mass range can be accurately found from thisformula.Then, since− ∂ f (R,z)∂R dRdM dM is equal to M times the total number of objects withmasses between M and M+ dM divided by the total mass, the number density ofcollapsed galaxy clusters with masses between M and M+dM is given bydn(M,z) =−ρ¯mM∂ f (R,z)∂RdRdMdM=√2piρ¯mδc3M2σe−δ2c2σ2(−RσdσdR)dM.(1.32)The matter power spectrum can then be inferred by measuring n(M,z) of collapsedgalaxy clusters as a function of mass and redshift, and using the equations (1.29)and (1.32) above. For a discussion of the ways of determining the masses of ob-served clusters, see ([34]).1.4.4 Baryon Acoustic OscillationsLastly, we give an overview of baryon acoustic oscillations as a probe of cos-mic expansion, our focus in this work (see ([37]), ([34]), ([50]), ([3]), ([45])). Inthe early universe, cold dark matter and baryons had a common origin in perturba-tions following inflation. Before the universe became transparent to light, baryonsand photons were strongly coupled due to Compton scattering and participatedin baryon acoustic oscillations, or longitudinal waves in the baryon-photon fluidcaused by gravity and pressure, and traveled a comoving distance of approximately150 MpC, called the sound horizon, away from the centers of the perturbations be-fore being released from the drag of the photons during recombination. This madea positive contribution to the correlation function of the matter fractional overden-sities at the comoving separation given by the sound horizon. This contribution11remained following the drag epoch, when the baryons were no longer coupled tothe photons and like the dark matter were only subject to gravity. Because of noisein the correlation function from cosmic variance set up by inflation, intermediatecorrelations near the sound horizon do not entirely cancel, so a peak will not al-ways be observed in this precise location; however, as stated earlier, the BAO peakhas proven to be a robust feature ([36]), ([22]), ([47]), ([78]), ([37]). Below weexplain more about how the expansion history has been estimated in the past usingbaryon acoustic oscillations and how it will be measured in experiments such asCHIME.It is possible to estimate the matter correlation function using galaxy redshiftsurveys, or, as a more recent development, hydrogen intensity mapping. As ex-plained in ([34]), galaxy redshift surveys measure the angular positions and red-shifts of galaxies. This can be used to estimate the matter correlation function. Aswe will explain later, when temperatures in the early universe became low enough,protons and free electrons were able to combine to form neutral hydrogen duringthe period of recombination. Later, as explained in ([90]), galaxies formed in thematter overdensities, which emitted light which reionized the atoms, with someremaining neutral hydrogen existing in dense clouds in galaxies. These are knownas damped Lyman alpha systems ([25]). Therefore, the neutral hydrogen fractionaloverdensities can be used to determine the total matter fractional overdensities, upto a possibly scale and time-dependent biasing factor ([25]). Due to the proton-electron spin interaction, there are two 1s hyperfine levels in neutral hydrogen, onetriplet state where the nuclear and electron spins are parallel, and the other singletstate where they are antiparallel, where the antiparallel spin has the lower energy([90]). As the neutral hydrogen transitions from the higher to lower of these energystates it emits a 21 cm wavelength, or 1420.405 MHz, photon, whose redshift canthen be determined from the relationν =1420.4051+ z. (1.33)Since baryon acoustic oscillations occur on significantly larger scales than galax-ies, it is useful to instead measure over larger scales the intensity of this neutralhydrogen signal as a function of frequency to avoid having to image the galaxy12distribution ([41]). We will explain later how the hydrogen intensity is related tothe matter power spectrum and how the map is constructed.Given a three-dimensional map in angular and redshift space, the correlationfunction, centered on a given redshift z, is parametrized by a separation ∆z in theline of sight direction and a separation θ in the transverse direction, where we onlytake the magnitude of the angular separation, and not its direction, to be important.The comoving separation, in terms of these parameters, is given by ([77])r(z,z′,θ) =√χ(z)2+χ(z′)2−2χ(z)χ(z′)cos(θ), (1.34)where χ(z) is given by (1.14). Assuming a value for the standard ruler comovingdistance rs, the correlation function near the peak can be fit to a model, to constrainthe parameters determining the expansion history.This can be done in a few different ways. In earlier galaxy surveys, such asthe Sloan Digital Sky Survey, it was not possible to resolve the three-dimensionalshape of the correlation function accurately. Therefore, a volume averaged samplewas taken assuming a fiducial model for translating angles and redshifts into dis-tances in order to calculate the volume. The SDSS team then fit C f it(r f id), the 3Dcorrelation function averaged over all directions according to the fiducial model, toa model as ([17]), ([89])C f it(r f id) = B2C f id(αr f id)+A(r f id). (1.35)Here, C f id is the correlation function of the fiducial model,A(r) =a1r2+a2r+a3 (1.36)with free parameters a1, a2, and a3, and the information about the true cosmologywas determined by the parameter α, known as the dilation, or change in scale,parameter, given byα =DV/rs(DV/rs) f id. (1.37)In this expression, DV is defined as the volume averaged distance to redshift z, andwas first used in ([37]). The best estimate of the distance to redshift z assuming only13knowledge of the correlation function in the radial direction at that redshift is czH(z) ,and the best knowledge of the distance to this redshift assuming only knowledgeof the correlation function in the transverse direction is χ(z) = (1+ z)DA(z). Sinceone averages over a volume, one mixes radial and angular information about thedistance to redshift z. One must then weight the dependence of this distance onthese parameters according the contribution that each makes to the volume overwhich the correlation function is measured. As the angular direction contains twodimensions whereas the radial direction contains only one, the volume averaginggivesDV (z) =((1+ z)2D2A(z)czH(z))1/3. (1.38)The dilation parameter α in (1.35) ensures that when r f id = rtrues , the correla-tion function can be expressed in terms of a change in scale of the old correlationfunction determined by the ratio of the volume averaged distances to redshift z:Ctrue(rtrues ) = B2C f id(DtrueVD f idVr f ids)+A(rtrues ). (1.39)With this value of α, they thereby fit to a model for cosmological parameters de-termining DV/rs, or, assuming a value for rs, such as the WMAP values, DV alone.There are some weaknesses to this method, as pointed out by ([76]). The firstis that it mixes angular and radial scales which depend on the cosmological param-eters in different ways, allowing degeneracy among cosmological parameters. Thesecond is that it relies on a fiducial model for the correlation function, which couldinduce cosmology dependence. As an alternative method, ([76]), ([77]) tested amethod which combined constraints from the purely radial and transverse mea-sured correlation functions. The transverse correlation function was fit to a modelwith parameters w and Ωm determining the angle θBAO(z) at which the BAO peakwas located, and similarly the radial correlation function was fit to a model withthese parameters determining a peak value ∆zBAO(z). The resulting constraints werecomparable to the previous method applied to the same data set. This method didnot consider the full shape of the correlation function for all intermediate directionshowever, which would provide more information. As we will explain, this is what14we hope to achieve in the method that we will use to analyze the CHIME data.1.5 The Friedmann-LeMaıˆtre-Robertson-Walker MetricA Friedmann- LeMaıˆtre-Robertson-Walker (FLRW) cosmological model, de-scribed as in ([82]), can be derived assuming homogeneity and isotropy. This as-sumption implies that the universe has a constant curvature everywhere. If thecurvature is positive, the spatial geometry of the universe is that of a 3-sphere, ifthe curvature is zero, the spatial geometry is flat in 3 dimensions, and if the curva-ture is negative, the spatial geometry is that of a 3-hyperboloid. These possibilitieshave metrics given by respectivelyds2 =−dt2+a2(t)(dΨ2+ sin2Ψ(dθ 2+ sin2θdφ 2))ds2 =−dt2+a2(t)(dΨ2+Ψ2(dθ 2+ sin2θdφ 2))ds2 =−dt2+a2(t)(dΨ2+ sinh2Ψ(dθ 2+ sin2θdφ 2)) (1.40)where t is the proper time of the observer in a frame where the spacelike hypersur-faces are homogeneous and isotropic.Our goal is to solve the Einstein equation,Gµν = 8piGTµν . (1.41)The left-hand side is a function of the metric, given by Gµν = Rµν − 12 gµνR whereRµν is the Ricci tensor, which can be written in terms of Christoffel symbols,Γαµν =12gασ (gµσ ,ν +gνσ ,µ −gµν ,σ ), (1.42)asRµν = Γαµν ,α −Γαµν ,α +Γαµν ,αΓαµν ,α −Γαµν ,αΓαµν ,α . (1.43)We now explain the right-hand side, Tµν , the stress energy tensor, as given in ([34]).We will define it for a particle species i in terms of the four-momentum componentspµ and the metric. The norm of the physical momentum is given bygµν pµ pν =−m2, (1.44)15where m is the rest mass.In terms of these, the stress energy tensor contribution from a species i is givenby ([44]), ([34]), ([86])T µν (~x, t) =−gi√| g |∫d p1d p2d p3(2pi)3fi(~x,~p, t)pµ pνp0, (1.45)where gi is the spin degeneracy, | g | is the absolute value of the determinant ofthe metric, and fi is the phase space distribution, equal to a Fermi-Dirac or Bose-Einstein function.The most general form of the stress energy tensor consistent with homogeneityand isotropy is given by that of a perfect fluid ([82]),T µν =−ρPPP , (1.46)where ρ and P are energy density and pressure. In a coordinate system where theoff-diagonal components of the metric are zero, the physical magnitudes of thespatial part of the four-momentum is given byp =√gi j pi p j. (1.47)For each species i, the contributions to ρ and P at a given temperature can then bewritten as ([34])ρi = gi∫ d3 p(2pi)3fi(~x,~p, t)E(p) (1.48)Pi = gi∫ d3 p(2pi)3fi(~x,~p, t)p23E(p). (1.49)We now explain the final form of the Einstein equations in an FLRW metric,as in ([34]), ([82]). To lowest order assuming spatial isotropy and homogeneity,evaluating the left and right-hand sides of the Einstein equations, we are left with16two equations:3a˙2a2= 8piGρ− 3Ka2(1.50)and3a¨a=−4piG(ρ+3P), (1.51)where K = 1,0,−1 depending on the sign of the curvature.(1.50) is called the Friedmann equation. Defining the Hubble parameter, thecritical density, and the density due to curvature, respectively byH ≡ da/dta, ρcr ≡ 3H208piG, ρk ≡ −3K8piGa2 , (1.52)the Friedmann equation can be written asH2H20=ρ+ρkρcr. (1.53)Here, ρcr has a measured valueρcr = 8.098×10−11h2 eV 4. (1.54)From here on, we will assume zero curvature, which is consistent with recent ob-servations ([56]). Therefore, ρk = 0 in (1.52), and the metric is given by the secondexpression of (1.40), which in Cartesian coordinates isds2 =−dt2+a2(t)(dx2+dy2+dz2) . (1.55)Later we will discuss the contributions to ρ, including radiation, matter, and darkenergy, and how they behave with the scale factor. These, together with the signof the curvature, determine the evolution of the scale factor through the Friedmannequation (see ([82]).1.6 Thermal History of the UniverseIn the next few sections, we will explain the thermal history of the early uni-verse, using ([34]), ([9]), ([10]), ([86]), including important events such as neu-17trino decoupling, big bang nucleosynthesis, and recombination. From the physicsof neutrino decoupling, we can determine the fractional contribution of the criticaldensity of radiation, which will be useful for the sound horizon calculation. Theinitial results from big bang nucleosynthesis will tell us the abundance of hydrogen,from which we can calculate the free electron fraction during recombination. Fromthis we can determine the optical depth which is used in the matter power spectrumcalculation. From now on, with only a few exceptions, we will use natural unitswhereh¯ = c = kB = 1. (1.56)Also, we note that from now on, overdots will refer to derivatives with respect tothe conformal time,η =∫ t0dta(t). (1.57)Finally, it is useful to keep in mind that our notation will be generally compatiblewith the conventions of ([45]) and related work by the same author.1.6.1 Statistics and Entropy in the Early UniverseIn this section, we will follow the discussion and derivations given in ([9]).As explained above, we have assumed that under the conditions of the big bang,all particles in the universe were in thermal equilibrium. A system is in ther-mal equilibrium when all the particles in the system are at a common tempera-ture, and is achieved when there is both kinetic equilibrium and chemical equi-librium. As explained in ([34]), ([9]), kinetic equilibrium implies that particles iscatter rapidly enough to take either Fermi-Dirac or Bose-Einstein distributions,fi = 1e(Ei−µi)/Ti±1 , where Ti can depend on the species. Chemical equilibrium im-plies that the ΣiµidNi = 0, and that the rate of a scattering reaction and its reversereaction are the same. In all calculations, we will assume conditions of kineticequilibrium for scattering reactions, whereas at later times we cannot always as-sume chemical equilibrium.A particle species i will behave relativistically when mi << Ti, so we can useE = p in the distribution functions. Photons always have a chemical potential ofzero. Although we will not discuss this here, it is considered a reasonable assump-18tion that µ in the early universe for all other particle species has always been small([86]), ([9]), ([10]). If the relativistic particles have negligible chemical potentials,we then haveni = gi∫ d3 p(2pi)31ep/Ti±1 =ζ (3)pi2giT 3i , bosons (1.58)=ζ (3)pi2giT 3i ×34, f ermions, (1.59)where gi is the spin degeneracy of the species.If T << E − µ, quantum statistics can be ignored, and both the Fermi-Diracand the Bose-Einstein distributions reduce to the Maxwell-Boltzmann distribution,given by f (E,µ,T ) = eµ/T e−E/T . We then haveni = gieµi/Ti∫ d3 p(2pi)3e−(mi+p22mi)/Ti= gie(µi−mi)/Ti(miTi2pi)3/2,(1.60)which shows that, for small µ, the particles in equilibrium become suppressed withdecreasing temperature after they become nonrelativistic.The first law of thermodynamics states that the change in entropy of a systemat temperature T is given bydS =dU +PdV −ΣiµidNiT(1.61)where, in chemical equilibrium, the last term vanishes. For relativistic particles,which have mi << Ti− µi, and assuming chemical potentials are negligible, wehaveρi = gi∫ d3 p(2pi)31ep/Ti±1E(p) = gipi230T 4i , bosons=78gipi230T 4i , f ermionsPi = gi∫ d3 p(2pi)31ep/Ti±1p23E(p)=ρi3.(1.62)19From (1.62), it follows that for relativistic species i,∂Pi∂Ti=ρi+PiTi. (1.63)For a system i in thermal equilibrium at temperature Ti, (1.61) can then bewritten asdSi =d((ρi+Pi)V )−V dPiTi=1Ti(d((ρi+Pi)V )−V (ρi+Pi)Ti dTi)= d(ρi+PiTiV).(1.64)The entropy density contribution, si, is then given bysi =ρi+PiTi. (1.65)We suppose that the universe consists of noninteracting systems i each in equi-librium at temperatures Ti. Each system obeys the continuity equation (1.192), and(1.63). We then havedSidt=ddt(ρi+PiTiV)=VTi(dρidt+1VdVdt(ρi+Pi))+VTi(dPidt− ρi+PiTidTidt)= 0.(1.66)If it is the case that the universe consists of noninteracting systems i each inequilibrium at temperatures Ti, so obeying (1.66), we can add all the contributionssi (1.65) to obtaina3s =2pi245g∗S(T )T 3a3 = constant, (1.67)whereg∗S(T ) = ∑i=bosonsgi(TiT)3+78 ∑i= f ermionsgi(TiT)3. (1.68)20This equation will be important for relating the photon and neutrino temperatures.It will also be useful to express the total energy density from relativistic contri-butions in terms of the temperature of the photons T, using (1.62), asρr =pi230g∗(T )T 4, (1.69)where g∗(T ) is the effective number of relativistic degrees of freedom, given byg∗(T ) = ∑i=bosonsgi(TiT)4+78 ∑i= f ermionsgi(TiT)4. (1.70)1.6.2 The Boltzmann EquationWe now explain the Boltzmann equation, as given in ([9]), see also ([34]),([86]), which gives the evolution of the number density n1 of particles undergoingthe scattering reaction 1+2↔ 3+4 between species 1, 2, 3, and 4. The Boltzmannequation is given byddt(n1a3) =−αn1n2a3+βn3n4a3. (1.71)It says that the rate of change in the comoving number density of particles, n1a3,is equal to the rate of change due to destruction, −αn1n2a3, where α, also writtenas 〈σv〉, is the thermally averaged cross section, plus the rate of change due tothe reverse process from creation, βn3n4a3. This has used the fact that the rateof a reaction per volume is proportional to the physical number densities of thereactants. As mentioned above, the condition of chemical equilibrium implies thatthe right-hand side vanishes. Therefore,β =(n1n2n3n4)equilibriumα, (1.72)and we have1a3ddt(n1a3) =−α(n1n2−(n1n2n3n4)equilibriumn3n4). (1.73)21It is easier to see how this equation behaves when we rewrite it asdln(n1a3)dlna=−αn2H(1−(n1n2n3n4)equilibriumn3n4n1n2). (1.74)From here we see that the tendency is for n1a3 at a given temperature to movein the direction of the chemical equilibrium ratio. If the reaction rates niα arelarge compared to H, then any deviations of nia3 from equilibrium will quicklygo towards equilibrium. However, when reaction rates are small compared to H,then the right-hand side will approach zero before the equilibrium values can bereached. The reaction rate is then said to freeze out, leaving a constant abundanceof each species per comoving volume corresponding to the time of freezeout. Inthe early universe during radiation domination, H may be computed in terms of thephoton temperature T using the Friedmann equation and (1.69) asH =pi3√8piG10g∗(T )T 2. (1.75)By calculating α(T )n2(T )/H(T ) for a particular reaction, we can determine therange of temperatures at which it was in equilibrium.1.6.3 The Neutrino and Photon Energy DensitiesIn this section, we will explain how to find the quantities Ωγ and Ων , fromconsiderations of early scattering events in the universe, and the subsequent contri-butions to the density, as explained in ([9]), ([34]). In the early universe, neutrinoswere in equilibrium with the baryon-photon plasma at a common temperature andinteracted with it through weak force mediated reactions such as ν+ ν¯ ↔ e++ e−and e−+ ν¯ ↔ e−+ ν¯ . From a Boltzmann equation calculation, one can find thatthe number density of the neutrinos froze out at about T = 0.8 MeV. Immediatelyfollowing this, when the reaction rate for electron-positron creation and annihi-lation was still greater than the Hubble rate of expansion, electrons, positrons,and photons followed their equilibrium distributions at the common temperatureof the plasma. As the temperature decreased and the electrons and positrons be-came nonrelativistic, the ratio of their equilibrium number densities, (1.60), to that22of photons, (1.58), got increasingly smaller due to the exponential factors, so theelectrons and positrons annihilated at a faster rate than they were created, and theirnumber densities were suppressed. After neutrino decoupling, we can approxi-mate the neutrinos and the rest of the plasma as noninteracting systems each inequilibrium, so their entropies were separately conserved.At the time of neutrino decoupling, before electron-positron creation and anni-hilation had frozen out, the neutrinos and plasma were at a common temperature,T1. Evaluating the total entropy (1.67) givess(a1)a31 =43pi290T 31 a31, (1.76)accounting for the fact that there were 2 spin states for photons, 2 spin states foreach of electrons and positrons, and 1 spin state for each of 3 species of relativisticneutrinos and antineutrinos.Equating the contributions to the neutrino entropy before and after annihilationusing (1.76), we can determine Tν after annihilation from the condition a1T1 =a2Tν . Doing the same for the plasma, we see that due to the change in g∗S ofthe plasma after annihilation, we cannot have T 3γ a3 = constant as is the case forneutrinos, implying that the plasma and neutrinos had different temperatures Tγ andTν . Evaluating the total entropy of the plasma plus neutrinos following electron-positron annihilation, we then haves(a2)a32 =2pi245(2T 3γ +6×78T 3ν )a32. (1.77)Equating s(a1)a31 with s(a2)a32 and eliminating T1 by substituting a1T1 = a2Tν , wefindTνTγ=(411)1/3. (1.78)Using this and (1.62), we haveρνργ=6×7/82(411)4/3. (1.79)The present-day energy density of photons is given by (1.62) using the present day23temperature of Tγ ,which is Tγ0 = 2.725K, and g= 2.We can defineΩγ byΩγ =ργρcrwhere ρcr is given by the value (1.54), and similarly Ων = ρνρcr . We therefore findΩγ =2.47×10−5h2, (1.80)so (1.79) givesΩν =1.68×10−5h2. (1.81)Following neutrino decoupling and electron-positron annihilation, g∗(T ) stayedconstant, with photons and neutrinos being the only remaining components to ra-diation. One can show, using the geodesic equation ([34]) that the energy of arelativistic particle changes proportional to a−1, whereas volume is proportional toa−3. Therefore, after, after neutrino decoupling and electron-positron annihilation,the contributions to radiation are given by ργ(a)/ρcr = Ωγ/a4, and ρν(a)/ρcr =Ων/a4, and for matter, ρm(a)/ρcr = Ωm/a3. We similarly define ΩDE , and Ωk asthe fractions of ρcr taken up today by the density due to dark energy and curvature.From (1.52), curvature scales as a−2. We can then write the Friedmann equation(1.53) asH2H20=Ωma3+Ωγa4+Ωνa4+ΩDEe−3∫ a0daa (1+w(a))+Ωka2. (1.82)For the Friedmann equation to hold at a = 1 this implies, for a constant w(a), that1 =Ωm+Ωr +Ων +ΩDE +Ωk. (1.83)1.6.4 Big Bang NucleosynthesisIt will be important for solving for the power spectrum to have a knowledge ofthe fraction of the total number density of baryons occupied by protons and hydro-gen atoms at the time of recombination. The hydrogen fraction will be importantfor our calculation of the free electron fraction for the optical depth during recom-bination. For this section, we will rely primarily on the account given in [86], withsome content from [34].We can find an expression for the fractional abundance Xi of the total nucleons24occupied by an element in equilibrium. Using the equilibrium Maxwell-Boltzmanndistributions, ignoring the binding energy, Bi = Zimp +(Ai− Zi)mn−mi and themass difference between protons and neutrons except in the exponential and there-fore using the mass mN of a nucleon which could be taken to be either a proton orneutron, we haveninZip nAi−Zin=gi2AiA3/2i(mNT2pi)3(1−Ai)/2eBi/T . (1.84)Writing this in terms of the fractional abundancesXi ≡ ninN , Xp ≡npnN, Xn =nnnN, (1.85)and definingε ≡ 12nN(2pi)3(2pimNT )−3/2, (1.86)we haveXi =gi2XZip XAi−Zin A3/2i εAi−1eBi/T . (1.87)Since ε is then very small in the early universe, we find that, assuming equilibriumin the early universe, considering the dominant factors in the above expression, aspecies will only occupy a significant fraction of the total nucleons ifTi ∼ Bi(Ai−1)|lnε| , (1.88)which ranges from T ∼ 0.75× 109 K for D to 3.1× 109 K for 4He, with similardependencies as that of 4He for more massive elements. However, by the timethe universe reached 3.1× 109 K, only two-body processes had appreciable reac-tion rates in the Boltzmann equation. The simplest of these was n+ p→ D+ γ.Therefore, D started its formation first, rather than the heavier elements. Otherreactions with appreciable rates included D+D→ 3H + p and D+D→ 3He+n,which froze out at T ≈ 109 K leaving only a small amount of deuterium. Also,D+ 3H→ 4He+n and D+ 3He→ 4He+ p took place rapidly enough that we canassume that by 109 K, almost all of the H3 and He3 had converted to He4. Thereare no stable nuclei with A = 5 or A = 8, so two-body processes could not produce25any heavier elements.Since almost all of the composite nuclei ended up as 4He,we can determine thefraction of baryons at the end of big bang nucleosynthesis occupied by hydrogenfrom the fraction occupied by helium. To a good approximation, using atomic massunits to relate mass to number densities,Y ≡ mHemtotal=4nHe4nHe+nH, (1.89)whereasXn ≡ nnnn+np , (1.90)wherenn = 2nHe, np = 2nHe+nH , (1.91)soY = 2Xn. (1.92)In the early universe, neutrons and protons were coupled through the weakforce mediated reactionsn+ν ↔ p+ e−, n+ e+↔ p+ ν¯ , (1.93)andn↔ p+ e−+ ν¯ , (1.94)or neutron decay.At these temperatures, which are before electron-positron annihilation had frozenout, leptons e+ and e− were in equilibrium and had the same abundance n(0)l . Asa first approximation we neglect neutron decay and therefore solve for the neutronfraction using the Boltzmann equation (1.73) for either of the first two reactions,which becomesa−3d(nna3)dt= n(0)l 〈σv〉{npn(0)nn(0)p−nn}, (1.95)wheren(0)i ≡ gi∫ d3 p(2pi)3e−Ei/T (1.96)26is the number density divided by e−µi/T , assuming Maxwell-Boltzmann statistics,and we have used the fact that in thermal equilibrium, µ1+µ2 = µ3+µ4. This maybe rewritten in terms of Xn ≡ nnnn+np and the neutron-proton conversion rateλnp ≡ n(0)l 〈σv〉=255τnx5(12+6x+ x2) (1.97)asdXndt= λnp{(1−Xn)e−Q/T −Xn}, (1.98)where e−Q/T = n(0)nn(0)pand Q≡mn−mp = 1.30 MeV, where x≡ QT , and τn = 886.7sec−1is the neutron lifetime. By comparing the reaction rate to the expansion rate, onefinds that the reaction froze out at a value of T ≈ 0.5 MeV, whereXn(t f reeze)≈ 0.15. (1.99)Next, we wish to account for the effects of neutron decay, (1.94). Since all neutronshave been decaying with times since production, this would require that Xn(t ≥t f reeze) ∝ e−t/τn . However, from the Friedmann equation and (1.75), we can findthat the time-temperature relation at the time after electron-positron annihilation isgiven by t = 132sec(0.1MeVT)2. Therefore, at freezeout, when T = 0.5MeV, therehad been negligibly small times available for decay compared to τn over the courseof production, so decays could not have significantly affected the result. Therefore,Xn(t f reeze)e−t/τn is a good approximation. After fitting the numerical results to adecaying exponential curve, we haveXn(t ≥ t f reeze)≈ 0.1609e−t/885. (1.100)At T ≈ 1 MeV, or 109K, the reaction for deuterium conversion fell out of equilib-rium as the reaction rate became smaller than H. At this time neutrons and protonswere nonrelativistic, decreasing g∗. Using the corresponding relation, this occursat a time of 168 seconds, during which, assuming that the 3H and 3He were imme-diately converted to 4He, we have, from (1.92),Yp ≈ 2×0.1609e−168/885 ≈ 0.27, (1.101)27with small adjustments given more detailed calculations. To be consistent with thecalculations in [86] of the free electron fraction, we will assume Yp = 0.76, so that,at recombination, using the expression for the baryon fraction,n≡ np+nH = 0.76nB = 0.76(3H20Ωb8piGmp)(TγTγ0)3. (1.102)1.6.5 The Free Electron Fraction From RecombinationIt will also be important to calculate the free electron fraction in order to cal-culate the optical depth τ, which, given a present time, parametrizes the amountof time looking back to a particular event in terms of an integral over the rate ofCompton scattering. The Compton optical depth is given by (see eg. ([50]))τ(η) =∫ η0dη ′τ˙, (1.103)whereτ˙ = neσT a′, (1.104)with ne the free electron number density. One can also specify the optical depthfrom the present, η0, to an earlier time η , given by, as in ([50]),τ(η) =∫ η0ηdη ′τ˙(η ′), (1.105)where we have τ˙ =−neσT a.For the calculation of the free electron fraction, we will follow the discussiongiven in ([86]). As the universe expanded, the electrons became bound into hydro-gen atoms. This period is called recombination, although it really refers to the firsttime that the electrons combined with protons to form hydrogen atoms. We willnow find the free electron fraction as a function of temperature. The initial con-ditions for this calculation depend on the initial fractional abundance of hydrogenfrom the results of Big Bang nucleosynthesis, (1.102). We will consider scatteringat T < 4400 K. At the beginning of this period, the reaction p+e↔H1s took placein thermal equilibrium. In general, the energy of a photon emitted in the transition28from a free electron to the ground state is such that it can go on to ionize anotherhydrogen atom, and similarly that of a photon emitted from a transition n ≥ 3 tothe ground state is such that it can go on to excite another atom in the ground stateto an n = 2 state, so that in both cases, these transitions produce no net effect. Theonly way the ground state can effectively be reached is through a series of radiativedecays, and the states of these atoms are therefore by our assumption in thermalequilibrium.If equilibrium were maintained as the ionization fraction decreased, so that theSaha equation gave the correct free electron fraction, then there would be no freeelectrons remaining in the universe, so to explain the existing free electron fraction,one must consider nonequilibrium conditions ([9]). A consideration of nonequilib-rium conditions can be done using some approximations explained below as in([86]), based on the original papers ([68]) and ([91]).During these times, electrons and hydrogen atoms were already nonrelativisticand we can use the Maxwell-Boltzmann distribution (1.60). Ignoring the massdifference between the proton and hydrogen except in the exponential, using thefact that due to equilibrium, the chemical potentials cancel in the numerator anddenominator, and assuming charge neutrality so that the number of protons andfree electrons are equal, we haven1snpne=n1sn2p=(meT2pi)−3/2eB1/T , (1.106)where B1 ≡mp+me−mH = 13.6 eV is the binding energy of hydrogen. Similarly,during thermal equilibrium, since µ1s = µnl, we have for the photon transition pro-cessesnnln1s∝ e−(Enl−E1s)/T , (1.107)where the energy difference in the exponent is at least 10.6 eV as for n = 2, so wecan assume to a good approximation that for temperatures below the mentioned4400 K, nH = n1s. We can then solve for the free electron fraction through thequantity Xe≡ npnp+n1s ≈npnp+nH, as a function of temperature using the Saha equation,29which holds in equilibrium, given byXe(1+SXe) = 1, (1.108)whereS≡((np+n1s)n1sn2p)∣∣equilibrium≈ 0.76nb(meT2pi)−3/2eB1/T ,(1.109)where we have used (1.102) with np + n1s ≈ n and (1.106), and n is known as afunction of temperature using (1.102). The Saha equation holds at early times, sothat one can get an estimate of the time of recombination, defined as the time whenXe = 0.1, as 0.3 eV ≈ 3600 K.Of the radiative decays to the ground state, the final transition from n = 2to n = 1 happens especially slowly, so it fell out of equilibrium faster than theother reactions. Therefore, at the time period we are considering, we assume thatonly the ground state, 1s, was not in equilibrium, whereas the other energy levelsof hydrogen were in equilibrium with the photons and each other and could beapproximated as taking their equilibrium number densities at the temperature ofthe photons.With the approximation of only 1s not in equilibrium, we can use the equi-librium number densities to find the number densities of nnl, n ≥ 2, in terms ofn2s. Ignoring the fine, hyperfine, and Lamb energy shifts, we can assume that thebinding energy is independent of l, resulting innnl = (2l+1)n2se(Bn−B2)/T . (1.110)We now consider the Boltzmann equations for the recombination and reionizationreactions e+ p↔ nnl for all n ≥ 2 and l. This is the same as the discussion onthe Boltzmann equation given earlier, except in this case there is only one par-ticle on the right-hand side. Again using the fact that the rate of a reaction pervolume is proportional to the physical number densities of the reactants, the net30rate of increase in the number of free electrons in a comoving volume a3 is givenby −αn2ea3, the rate of recombination, which is obtained after collisions betweenelectrons and protons, plus βn2sa3, the rate of reionization resulting from excitedstates nnl, proportional to n2s from (1.110), where α and β are coefficients whichdepend on temperature. The Boltzmann equation is thenddt(nea3) =−αn2ea3+βn2sa3. (1.111)Dividing by the constant na3, where n≡ np+nH is given by (1.102), then gives anequation for the evolution of the free electron fraction,ddt(nen)=−αn2en+βn2sn, (1.112)which we will use to replace the equilibrium version of the Saha equation. Thecoefficient α can be calculated to beα = 2.84×10−11T−1/2cm3s−1. (1.113)As before, the ratio of α, the thermally averaged cross section, to β can be fixedusing the fact that if the electrons, protons and excited states are in chemical equi-librium then the right-hand side of (1.111) must vanish, so we haveβα=(n2en2s)equilibrium=(meT2pi)3/2e−B2/T . (1.114)In order to find a complete equation for the free electron fraction, we must thenfind an expression for n2s in (1.112) in terms of the free electron fraction, or equiv-alently ne and nH . For this, we make the additional assumption that the number ofexcited states of hydrogen following recombination changes very slowly comparedto the number in the 1s state. In other words, we assume that the change in the num-ber density of atoms in an excited state nnl due to recombination and reionizationexchanges with free electrons cancels with its change due to exchanges with the 1s31state. Therefore, we use the equalityαn2e−βn2s = κn2s− εn1s, (1.115)where κn2s is the rate of decays from n = 2 to the ground state, and εn1s is the rateof the reverse process. In a similar way as for the ratio β/α, since if the 1s andexcited states are in chemical equilibrium, there is no net rate of change in the n1sstate, we haveε = κ(n2sn1s)equilibrium= κe−(B1−B2)/T . (1.116)As explained above, since the temperatures at which we consider the Boltz-mann equation are much less than the binding energy difference between the sec-ond and higher states, and these states are assumed to be in equilibrium, (1.110)shows that the number density of higher excited states is small, so in (1.115) wecan use the approximationn1s ≈ nH −n2s−n2p = nH −4n2s, (1.117)which will allow us to find an expression for n2s in terms of nH and ne.We can find an expression for κ as follows. One of the ways the final transitionfrom n= 2 to the ground state 1s can occur is the transition 2p→ 1s+γ, which hasa rate Γ2p, which can be calculated from the solution to the quantum mechanicalscattering problem. The emitted photon in this transition, known as a Lyman αphoton, can bring another hydrogen atom in the 1s state to a 2p state, but hasbarely enough energy to do so, and therefore, due to the cosmological redshift ofits frequency, it will not succeed unless it quickly reaches another 1s hydrogenatom. Also considered is the transition 2s→ 1s+ γ + γ, occurring at a rate Γ2s,where the two emitted photons do not have enough energy to produce an excitedstate from another 1s state. The effective decay rate from excited states to theground state is then given byκn2s = (Γ2s+3PΓ2p)n2s, (1.118)where we have used (1.110) to express n2p in terms of n2s, and where P(t) is the32probability that the Lyman alpha photon emitted at time t does not re-excite another1s state to a 2p state. By the Breit-Wigner formula, the cross section for the Lymanalpha photon exciting an n1s state to an n2p state is given, in terms of the meanfrequency ωα = B1−B2 of the Lyman α photon, byσ(ω) =32(2pi2Γ2p(ωα)2)P(ω), (1.119)whereP(ω)dω is the probability that the Lyman alpha photon has frequency be-tween ω and ω+dω, which is given byP(ω)dω =Γ2p2pi1(ω−ωα)2+Γ22p/4dω. (1.120)In terms of these, the probability that a Lyman alpha photon emitted at time t willnot excite a 1s goes as the negative exponential of the number of scattering eventsafter a time t, and is given byP(t) =∫ ∞−∞dωP(ω)exp(−∫ ∞tdt ′n1s(t ′)σ(ωa(t)/a(t ′))), (1.121)where the factor of a(t)/a(t′) in (1.121) accounts for cosmological redshift of thefrequency ω at a time t′. Performing the integral with some approximations, andconsidering the relation between the scale factor, time, and temperature, one canfind that this expression can be simplified at temperatures less than the temperatureswe have considered to approximatelyP(t) =8piH3λ 3αΓ2pn1s, (1.122)where λα = 2piω−1α . Therefore, (1.115) becomesαn2e−βn2s = (Γ2s+3PΓ2p)n2s− εn1s, (1.123)where from (1.116),ε = (Γ2s+3PΓ2p)e−(B1−B2)/T , (1.124)and β can be related to the known value of α by (1.114). Substituting (1.117) into33(1.123) and solving for n2s, we obtainn2s =αn2e + εnHΓ2s+3PΓ2p+β +4ε. (1.125)Substituting this into the Boltzmann equation (1.112) along with the results (1.114),(1.121), and (1.116), and using the approximation that T << B1−B2 in the expo-nent of the last expression, we obtainddt(nen)=(Γ2s+3PΓ2pΓ2s+3PΓ2p+β)(−αn2en+α(meT2pi)3/2e−B1/TnHn). (1.126)Finally, substituting Xe, this time defined as Xe ≡ nen = 1− nHn , andS = n(meT2pi)−3/2eB1/T , (1.127)as in (1.109), the Boltzmann equation can be entirely expressed in terms of oneunknown variable Xe asdXedt=(Γ2s+3PΓ2pΓ2s+3PΓ2p+β)αn(−X2e +S−1(1−Xe)). (1.128)The equation can be written with the temperature as the time variable using the factthat a ∝ T−1 during radiation domination, sodtdT=− 1HT, (1.129)to obtaindXedT=αnHT(1+βΓ2s+8piH/λ 3αn(1−Xe))−1(X2e − (1−Xe)/S). (1.130)From (1.102), we can then find ne(T ). ne(z) can the be found given ne(T ),using (1.75) and H(z). From (1.104), we can then determine τ(z). A fitting formulafor the optical depth from the present to a time η during recombination is given by34([45])τ(z) =Ωc1b( z1000)c2, (1.131)wherec1 = 0.43, c+2 = 16+1.81lnΩb. (1.132)1.7 Cosmological Perturbation TheoryBelow, we explain the theoretical predictions for the evolution of the matterdensity perturbations to linear order, from which the power spectrum and the mat-ter correlation function are obtained. The perturbations to the metric and the contri-butions to the stress energy tensor obey the Boltzmann equations and the Einsteinequations. From these equations, the variances of the solutions to the fluctuationsare obtained in terms of their initial values, which are believed to have arisen fromquantum mechanical fluctuations in the inflaton field. To give the reader an under-standing for the physics processes involved, we will discuss an analytical solutionto the power spectrum, based on those given in ([45]) and related work by thesame author referenced in the text, and references therein. For an alternative ac-count of these equations, see also ([62]). However, to solve the equations for thepower spectrum to the required accuracy for our project, numerical calculationsare required, and we will use the CAMB (Code for Anisotropies in the MicrowaveBackground) online resource, available at the link given in ([4]).1.7.1 Degrees of Freedom for PerturbationsWe now explain a general decomposition for the perturbations to the FLRWmetric, originally shown by Lifshitz ([60]), using ([21]) and ([30]). A general per-turbed FLRW metric in the form below can be written as a sum of an unperturbedpart plus a perturbation hµν ,ds2 = a2(η)(−dη2+ γi j(~x)dxidx j +hµν(η ,~x)dxµdxν), (1.133)where γi j is the unperturbed 3-metric without factors of a2, here equal to δi j sincewe are assuming zero curvature. The perturbations hµν can be written as a sum ofparts with nonzero components which transform only among themselves under a35spatial rotation as a scalar for the time-time component, vector for the time-spacecomponents, and a spatial tensor for the spatial components. The spatial tensor canbe then be decomposed again into self-transforming parts as the sum of a trace orfunction of a scalar, and a traceless part with respect to the three-metric, γ i jSi j = 0,where γ i j is the inverse of γi j. We write this as follows:h00 = 2Ψ, (1.134)h0i = wi, (1.135)hi j = 2(Φγi j +Si j). (1.136)We can further decompose the vector wi and tensor modes Si j into functions ofscalars, vectors, and tensors each of which transform only among themselves. Bythe Helmholtz theorem, a smooth vector field wi can be decomposed aswi = w⊥i +w‖i , (1.137)where, using the covariant derivative ∇i in terms of γi j,w‖i = ∇iφ (1.138)for a scalar field φ , and∇ jw⊥j = 0. (1.139)As an extension of this theorem, a symmetric tensor Si j can be decomposed intoscalar, vector and tensor parts, asSi j = S‖i j +S⊥i j +STi j, (1.140)where S‖i j is obtained from the divergence of a scalar field φs according toS‖i j = (∇i∇ j−13γi j∇2)φs, (1.141)36S⊥i j is obtained from a transverse, or divergence free, vector field S⊥i according toS⊥i j = ∇iS⊥j +∇ jS⊥i , (1.142)and finally, STi j is a tensor mode which obeys∇iSTi j = 0. (1.143)The ten degrees of freedom of a general symmetric four -by-four matrix arethen accounted for by 4 scalars, two vectors each with two independent compo-nents, and one tensor with two independent components. The decomposition wehave written down here will be important because in linear perturbation theory,we can decompose Einstein’s equations in such a way that the scalar, vector, andtensor perturbations can be solved for independently of each other. Physically,the equations for the vector modes of the metric perturbation correspond to thoseof gravitomagnetism, and those of the tensor modes correspond to gravitationalwaves. In this thesis, we will only be interested in the scalar perturbations, whichdetermine the matter power spectrum.The remaining part of this discussion will be based on ([45]) and referencestherein, which uses notation consistent with the original paper ([53]). At any giventime, a perturbation as a function of comoving coordinates ~x can be decomposedinto Fourier modes~k. For linear perturbation theory, we find that all modes~k ofthe perturbations evolve independently ([34]), which can be verified by Fouriertransforming the equations of motion. Given this independent evolution, we willwrite the metric used to obtain the equation of motion for the ~kth mode as theunperturbed metric plus the perturbation due to the contribution from this mode.For an FLRW metric, this equals the Fourier coefficient times a function Q whichis an eigenfunction of the Laplacian,∇2Q = γ i j∇i∇ jQ =−k2Q, (1.144)where the covariant derivatives are taken using the 3-metric γi j. Therefore, if thereis zero curvature,Q = ei~k·~x. (1.145)37k is known as the comoving wave number. It is useful because kη equals the ratioof the comoving horizon to the comoving wavelength of the perturbation ([34]).The~kth mode contributions to the scalar perturbations above, w‖i and S‖i j, arethen proportional to the operations acted on the scalar Q given byQi =−1k∇iQ (1.146)andQi j =1k2∇i∇ jQ− 13γi jQ. (1.147)Indices on Qi and Qi j can be raised and lowered by the 3-metric γi j.Therefore, to obtain the equation of motion for the mode ~k, we can use themetric with scalar perturbations written asg00 =−a2(1+2AGQ)g0 j =−a2BGQ jgi j = a2(γi j +2HGL Qγi j +2HGT Qi j).(1.148)The perturbations to the stress energy tensor follow the same decompositioninto self-transforming parts under rotations as the perturbations to the metric ([86]).We can write the stress energy tensor contribution from a fluid or fluid combinationx with scalar perturbations as ([53]), ([45])(T 00 )x =−(1+δGx Q)ρx(T 0i )x = (ρGx + pGx )VGx Qi(T j0 )x =−(ρGx + pGx )V Gx Q j(T ij )x = pGx (δij +δ pGxpGxδ ijQ+ΠGx Qij).(1.149)Here, δx = δρρ , Πx is the anisotropic stress, and, for consistency with a later expla-nation,V Qi =V Qi = adxidτ, (1.150)where τ is the proper time. We can then relate the components of the stress energy38tensor of all the matter, (T µν )T , to the individual contributions (Tµν )x, by writing(T µν )T = Σx(Tµν )x summing over all the fluids x.In the above expressions for the perturbed metric and stress energy tensor, thesuperscripts G indicate that the choice of perturbations is not unique to the physics,but depends on the choice of gauge. This is a consequence of the Bianchi identity∇µGµν = 0, which places constraints on Rµν , and therefore introduces redundancyin the Einstein equations ([30]). Any gauge transformation capable of influencingscalar perturbations can be given by ([45])η˜ = η+T Qx˜i = xi+LQi,(1.151)recalling the definitions (1.145) and (1.146) for Q and Qi. Since there are two de-grees of freedom that cause the coordinate changes in the scalars, we can eliminatetwo of the four scalar coordinate degrees of freedom by a change of coordinates. Ina similar way when dealing with vector perturbations, one of the two-componentvector fields can be eliminated, leaving only 6 degrees of freedom in total, consis-tent with the four additional constraints of the Bianchi identity. The gauge freedomand the choices for it used in this thesis will be explained in the next subsection.1.7.2 Choosing a Gauge for Scalar PerturbationsWe now wish to relate the changes in the components of the scalar metric per-turbations to the coordinate transformation. In what follows, we will continue touse the explanations and notation of ([45]), originally given in ([53]); see also([86]). As explained in Appendix A, under a gauge transformation, xλ ′ = xλ +ξ λ ,the metric and energy-momentum tensors transform by the Lie derivative with re-spect to a four-vector ξ λ , which to first order in hµν and ξ λ gives∆hµν(x) =−g¯λµ(x)∂νξ λ (x)− g¯λν(x)∂µξ λ (x)−∂λ g¯µν(x)ξ λ (x), (1.152)and similarly to first order in the perturbations to the stress energy tensor and39ξ µ , we have∆δTµν(x) =−LξTab=−T¯λµ(x)∂νξ λ (x)− T¯λν(x)∂µξ λ (x)−∂λ T¯µν(x)ξ λ (x)(1.153)([86]).Evaluating (1.152), we can therefore find that the scalar metric components(1.148) change under the coordinate transformation (1.151) byAG˜ = AG− T˙ − a˙aTBG˜ = BG+ L˙+ kTHG˜L = HGL −k3L− a˙aTHG˜T = HGT + kL,(1.154)and also, from (1.153),V G˜x =VGx + L˙δ G˜x = δGx +3(1+wx)a˙aTδ pG˜x = 3c2xρx(1+wx)a˙aTΠG˜x =ΠGx ,(1.155)wherec2x ≡p˙xρ˙x. (1.156)One very commonly used gauge is the conformal Newtonian, also known aslongitudinal, gauge. The coordinates used for this gauge are determined from theconditionBN = HNT = 0. (1.157)From now on, we will also relabel the remaining conformal Newtonian gauge vari-ables, now given by functions of k, asA =Ψ, HL =Φ. (1.158)40The total matter gauge is defined by the condition that HTT = 0, as in the New-tonian gauge, and also the condition that the total four-velocity of all the mattercomponents of the universe is orthogonal to the constant time hypersurfaces. Asexplained in ([1]), the four-velocity of a matter component is given byuµ =(dηdτ,dxidτ), (1.159)where τ is the proper time in the rest frame of the matter, and satisfies uµuµ =−1.The second condition of the total matter gauge then implies that (see ([44]))uµdxµ = 0 (1.160)if dη = 0. By the four-velocity normalization, we then haveuµ = (−a(1+AQ),0,0,0), (1.161)so that raising an index on uµ gives the general expressionuµ =(1a(1−AQ), 1aBQi). (1.162)Therefore, comparing with (1.159), as stated in ([53]), ([44]),V = B, (1.163)where V is the magnitude of the spatial component of the velocity perturbation,V Qi = adxidτ. (1.164)We can convert the variables in the total matter gauge, denoted with a super-script T, from the Newtonian gauge, denoted with a superscript N, using (1.155)and the conditions of the Newtonian gauge (1.157), as∆x ≡ δ Tx = δNx +3(1+wx)a˙aV NTk, (1.165)41V Tx =VNx . (1.166)One of the advantages of using the variable ∆T is that, in addition to having theinterpretation as an overdensity in the total matter gauge, it is also a gauge invari-ant quantity ([34]), ([53]). One can verify, using the relation δT = 1ρT (Σiρiδi), thatat late times when ρr → 0, we can neglect the radiation contributions to the over-density so that ∆T → ∆m and we can therefore use ∆T to define the matter powerspectrum in the total matter gauge, as done in ([45]). Also, for modes k far belowthe horizon, ∆T reduces to δNT , and so we can use this when solving for small scalesolutions ([53]), ([34]).When solving for the matter power spectrum, most of the time we will followthe method given in ([45]), which is to work in the Newtonian gauge, using thecorresponding metric perturbation variables Φ and Ψ, as well as temperature per-turbations found from Θ= 4δNγ , but to write overdensities of matter components xin terms of a gauge invariant variable ∆x, which represents overdensites in the totalmatter gauge. Anisotropic stress and velocities are the same in both gauges so donot need to be distinguished.However, one additional gauge, which we will use when solving for initialconditions, will be the gauge with spatially flat slicing, defined by Φ= 0, HT = 0([34]). Finally, also when solving for initial conditions, we will make use of twogauge-invariant quantities, given in ([34]), originally demonstrated in ([18]),v≡ ikB+ kˆiδT 0i(ρ+P)(1.167)andΦH ≡−Ψ+aH(B+ H˙T ). (1.168)1.7.3 The Boltzmann Equation for PhotonsPhoton temperature and polarization perturbations obey a Boltzmann equationwith collision terms due to Compton scattering with free electrons. The results ofthese calculations were first given in ([24]). We will not go into a full derivationhere; however, a good account is given in ([86]) which explains the basic form andmotivation of the equations. Another thorough derivation is given in the Newtonian42gauge in ([57]); see also ([34]), ([62]), and ([66]) for more discussion. The equa-tions quoted in this section are those given in ([45]) which writes them in terms ofmultipole moments.As explained in ([7]), if the axes are chosen such that the electromagnetic wavetravels in the direction zˆ, the electric field part of the wave can be expressed asa linear combination of wave components polarized in the x direction and the ydirection with an arbitrary choice of phase for each:~E =(ExeiδxEyeiδy)eikz−ωt = E(cos(φ)eiδxsin(φ)eiδy)eikz−ωt , (1.169)for some angle φ , where E =√E2x +E2y and δx and δy are the phases for the wavecomponent polarized in the x and y directions. The normalized two-dimensionalvector in parentheses in the second equality is known as a Jones vector and de-scribes an arbitrary state of electromagnetic wave polarization. The followingdiscussion of Stokes parameters is based on ([57]). The polarization state of anelectromagnetic wave can be specified in terms of time averages depending on theindividual components and phases in terms of Stokes parameters,I = 〈E2x 〉+ 〈E2y 〉Q = 〈E2x 〉−〈E2y 〉U = 〈2ExEycos(δx−δy)〉V = 〈2ExEysin(δx−δy)〉.(1.170)Left and right hand circularly polarized light corresponds to a Jones vector 1√2(1, i)and 1√2(1,−i), respectively ([7]). As explained in ([86]), linear polarization cor-responds to a linear combination of left and right hand circular polarization withequal absolute values for each component, and hence a real Jones vector up to anoverall phase; it is expected that all photons are linearly polarized from Comptonscattering, and V = 0. I determines the intensity. The parameters Q and U are co-ordinate system dependent. Under a counter clockwise rotation by φ of the x and y43axes, they transform asQ′ = Qcos(2φ)+Usin(2φ)U ′ =−Qsin(2φ)+Ucos(2φ),(1.171)and also under a rotation by φ , the angleα =12tan−1(UQ)(1.172)transforms asφ ′= α−φ , (1.173)and is the same for α + pi. This can be generalized for an arbitrary direction ofpropagation nˆ for axes perpendicular to this direction. Therefore, the polarizationfield in a direction nˆ can be represented by a line segment symmetric about theorigin in the axes formed from this plane, of lengthP =√Q2+U2 (1.174)in the direction α with respect to the current axes.Quantum mechanically, photons possess polarization degrees of freedom intheir state space. If the axes are chosen such that the photon travels in the directionzˆ, we can express the polarization state of a photon as a linear combination oforthogonal basis states |εx〉 and |εy〉 polarized in the x and y directions respectivelyas| ε〉= axeiδx |εx〉+ayeiδy |εy〉. (1.175)In the quantum mechanical description, the polarization state in terms of an orthog-onal basis takes the same form as the Jones vector, and the values for the observableStokes parameters are replaced by expectation values 〈Ψ | Oˆ |Ψ〉 of the following44corresponding operatorsIˆ = |εx〉〈εx|+|εy〉〈εy|Qˆ = |εx〉〈εx|−|εy〉〈εy|Uˆ = |εx〉〈εy|+|εy〉〈εx|Vˆ = i(|εy〉〈εx|−|εy〉〈εx|).(1.176)Comparing the coefficients in the expansion of ρ = |Ψ〉〈Ψ| in terms of | εi〉〈ε j |with the expectation values 〈Ψ | Iˆ + Qˆ | Ψ〉, and likewise for other operators, wefind that the density matrix for a system of photons with polarization state | Ψ〉 isgiven byρ = IIˆ+QQˆ+UUˆ +VVˆ =12(I+Q U− iVU + iV I−Q), (1.177)where I is the expectation value of the operator Iˆ, and likewise for the other opera-tors. Here, the trace of ρ corresponds to the intensity. We note that in the quantummechanical case, as explained in ([86]), the polarization state of a photon now tellsthe information about its spin state, or helicity. The polarization states of photonswith the Jones vectors corresponding to left and right hand circularized light givenabove are states of definite corresponding helicity, or spin about their direction ofmotion.To zeroth order, there is no polarization for photons in the early universe, sowe must have Q =U = 0 as well as V = 0. To first order, introducing additionalfactors, we can write the perturbation to the components of the density matrix for asystem of photons as ΘT = ∆TT , ΘQ =QT , and ΘU =UT . In terms of the same units,we also have ΘP =√Θ2Q+Θ2U .For a given direction of propagation nˆ, the Boltzmann equation for photons inthe presence of Compton scattering, taking into account polarization and the pres-ence of gravitational field perturbations, is given in Fourier space by two coupledequations for the fractional temperature perturbations ΘT and fractional polariza-45tion perturbation strength ΘP ([24]), see ([34]), ([86]) ([45]), ([62]):Θ˙T + ikµΘT =−Φ˙− ikµΨ+ τ˙(−ΘT +ΘT 0− iµVb− 12P2(µ)Π) (1.178)Θ˙P+ ikµΘP = τ˙(−ΘP+ 12(1−P2(µ))Π), (1.179)where µ ≡ kˆ · nˆ, τ˙ = neσT a here taking τ to be the Compton optical depth, and neis the free electron number density, Vb is the electron velocity perturbation, andΠ=ΘT 2+ΘP2+ΘP0. (1.180)For an arbitrary wave vector~k and nˆ, we can orient the axes such that ΘQ =ΘPand ΘU = 0. Since ΘP depends on nˆ only through µ, we can then expand thevariables Θ and ΘQ into Legendre polynomials Θl(k) usingΘl(k) =il(2l+1)∫ 1−1dµ ′2Pl(µ ′)Θ(k,µ ′), (1.181)take moments, and substitute the recursion relation(l+1)Pl+1(µ) =−i(2l+1)µPl(µ)− lPl−1(µ) (1.182)to find the following equations for Θ and ΘQ respectively ([24]), ([50]):Θ˙0 =−k3Θ1− Φ˙ (1.183)Θ˙1 = k(Θ0+Ψ− 25Θ2)− τ˙(Θ1−Vb) (1.184)Θ˙2 = k(23Θ1− 37Θ3)− τ˙(910Θ2− 110ΘQ2 −12ΘQ0)(1.185)Θ˙l = k(l2l−1Θl−1−l+12l+3Θl+1)− τ˙Θl, l > 2 (1.186)andΘ˙Q0 =−k3ΘQ0 − τ˙(12ΘQ0 −110(Θ2+ΘQ2 ))(1.187)46Θ˙Q1 = k(ΘQ0 −25ΘQ2)− τ˙ΘQ1 (1.188)Θ˙Q2 = k(ΘQ1 −37ΘQ3)− τ˙(910ΘQ2 −110Θ2− 12ΘQ0)(1.189)Θ˙Ql = k(l2l−1ΘQl−1−l+12l+3ΘQl+1)− τ˙ΘQl , l > 2. (1.190)Using the notation Nl in place of Θl, we can also write analogous equations formassless neutrino perturbations, which then behave the same as photons except forthe absence of polarization, and with τ˙ = Perturbed Einstein Equations and Conservation EquationsThe solution for the power spectrum depends on the following equations asgiven in ([45]). As a consequence of the Einstein equation and the Bianchi identity∇µGµν = 0, we have the conservation law∇µT µν = 0 (1.191)([30]). The ν = 0 component of this equation is the continuity equation. To zerothorder in perturbations, this is given byρ˙xρx=−3(1+wx) a˙a . (1.192)As in the expression for the stress energy tensor (1.149), here, and from now on,x can represent either the sum total of the matter, or a single fluid if there are noeffects from interactions. This equation to first order in perturbations isδ˙x =−(1+wx)(kVx+3Φ˙)−3 a˙aδwx, (1.193)whereδwx =px+δ pxρx+δρx−wx=(δ pxδρx−wx)δx.(1.194)47In the tight coupling limit, when Compton scattering was taking place, and be-fore recombination, photon and baryon number was conserved separately, and wetherefore have a separate continuity equation for photon and baryon perturbations:δ˙b =−kVb−3Φ˙, (1.195)Θ˙0 =−k3Vγ − Φ˙ (1.196)([46]), see also ([43]), ([45]). Comparing (1.196) with the monopole equation(1.183) then gives, as in ([45]),Θ1 =Vγ . (1.197)Also, cold dark matter obeys the same continuity equation as for baryons, andneutrinos obey the same continuity equation as photons.The Euler equation, or ν = i component of (1.191), is given byV˙x =(− a˙a(1−3wx)− w˙x1+wx)Vx+δ px/δρx1+wxkδx− 23wx1+wxkΠx+ kΨ. (1.198)From the continuity equation (1.192), we havew˙x =ρ˙xρx(pxρ˙x−wx)(1.199)=−3(1+wx)(c2x−wx)a˙a, (1.200)wherec2x ≡p˙xρ˙x. (1.201)We can definewxΓx =(δ pxδρx− c2x)δx. (1.202)Using (1.199) and the definition (1.202), we can rewrite (1.198) asV˙x =− a˙a(1−3c2x)Vx+c2x1+wxkδx+wx1+wxkΓx− 23wx1+wxkΠx+ kΨ. (1.203)To allow for interactions for baryons, we can use the Boltzmann equation, obtained48in ([34]) following a similar method as for photons. However, since we showed thatthe velocity evolution for photons is given by (1.184) in the Boltzmann equationsfor photons, where the interaction terms have been previously calculated, we canalso find the baryon evolution using this and the fact that Compton scattering en-forces momentum conservation ([50]), ([45]). The momentum densities for thephotons and baryons are given by(ργ + pγ)Vγ =43ργVγ (1.204)and(ρb+ pb)Vb ≈ ρbVb, (1.205)so assuming the densities vary slowly compared to the velocity,1R(V˙γ)interaction+(V˙b)interaction = 0. (1.206)For the l = 1 component of the photon temperature Boltzmann equations, the in-teraction term was(V˙γ)interaction =−τ˙(Vγ −Vb). (1.207)Therefore, adding the compensating term to the non-interacting baryon terms in(1.203) gives the baryon Euler equationV˙b =− a˙aVb+ kΨ+τ˙R(Θ1−Vb). (1.208)Also, we note by comparing the non-interacting terms from the Euler equation(1.198) for photons and the photon dipole equation (1.184), and analogous equa-tions for the neutrinos, that Πγ = 125 Θ2, Πν =125 N2, so, adding the contributionsto the stress energy tensor (1.149), we havepTΠT =125(pγΘ2+ pνN2). (1.209)The cold dark matter does not interact with the other fluids, so the continuityand Euler equations are the same as for the baryons except without the coupling49term. Combining them gives the equation for the evolution given by ([50])δ¨c+a˙aδ˙c =−k2Ψ−3Φ¨. (1.210)As explained in ([45]), we can find equations for Ψ and Φ using the Einsteinequations. The time-time component, to first order in the perturbations, is3(a˙a)2Ψ−3 a˙aΦ− k2Φ=−4piGa2ρTδT . (1.211)The time-space component isa˙aΨ− Φ˙= 4piGa2(1+wT )ρT VTk (1.212)where the T subscript indicates the net contribution of all the components. Com-bining the time-time and time-space equations gives the Poisson equationk2Φ= 4piGa2ρT(δT +3a˙a(1+wT )VTk)= 4piGa2ρT∆T(1.213)when written in terms of the total matter gauge overdensity. The traceless space-space component of the Einstein equations givesk2(Φ+Ψ) =−8piGa2 pTΠT (1.214)where ΠT is given by (1.209).1.8 Solutions for the Matter Power SpectrumIn the solutions for the matter power spectrum, we will distinguish betweensuperhorizon and subhorizon scales. The comoving horizon is defined asη =∫ ttedt′a(t′) =∫ aaeda′a′2H(a′) (1.215)50as the total comoving distance that light could have traveled since the end of in-flation. We see that it is the integral over the logarithm of the comoving Hubbleradius, 1aH ([34]). If perturbations are on a scale greater than the Hubble radius,no causal physics will have influenced them in the time that it takes for the scalefactor to increase by a factor of e. During radiation and matter domination, 1aH ≈ ηand we will use η in place of 1aH ([45]).1.8.1 Initial Conditions for the Boltzmann EquationsWe now trace back the equations to initial conditions following inflation, using([34]), ([58]), ([45]), ([50]). We will generally assume adiabatic initial conditionsfor all the components. Adiabatic initial conditions require that the entropy pertur-bations, given bySxy ≡δ(nxny)(nx/ny)(1.216)are zero for all components x, y in the universe, and also that the bulk velocities ofall components are essentially zero initially at the end of inflation. We then have,in terms of the total matter gauge variables,Sxy =δnxnx− δnyny=∆x1+wx− ∆y1+wy= 0. (1.217)The second equality can be verified for the case of nonrelativistic particles wherew = 0, and for radiation using the expression ρi = gipi230 T4 given by (1.62) and thefact that ni =giT 3pi2 after evaluating the integral (1.96). Therefore, we have∆c = ∆b = ∆m, (1.218)and∆γ = ∆ν = ∆r =43∆m. (1.219)We will now relate Θ0(k,0) to the potential, Φ(k,0). Using the fact that duringthis time, radiation dominates, and so a∝η , the Einstein equation (1.211) becomes51during radiation dominationΦ˙η−Ψ= 2Θr,0. (1.220)Here,Θr,0 = (1− fν)Θ0+ fνN0, (1.221)where fν ≡ ρνρν+ργ ≈ .405. As explained in ([34]), at times early enough followinginflation, wavelengths of all perturbations could be taken to be arbitrarily greaterthan the horizon, allowing terms containing k in the Boltzmann equations, as wellas higher moments of the Boltzmann equations, to be neglected. Also, couplingterms can be neglected since no causal physics can occur for scales above the hori-zon. The Boltzmann continuity equations then become, in terms of Newtonianoverdensities,Θ˙0 = ˙N0 =−Φ˙, (1.222)δ˙c = δ˙b =−3Φ˙. (1.223)Differentiating (1.220) and using the Boltzmann equations (1.222), we haveΦ¨η+ Φ˙− Ψ˙=−2Φ˙. (1.224)We note from the Einstein equation (1.214) that neglecting anisotropic stressleads to the condition Ψ=−Φ, so this becomesΦ¨η+4Φ˙= 0, (1.225)which has a power law solution p = 0 or p = −3. Ignoring the p = −3 solutionwhich is a decaying mode, (1.220) implies that we then haveΦ=−Ψ= 2Θr,0 (1.226)which is constant. By the earlier assumption, we have Θ0(k,0) =N0(k,0) whichis then a constant, and so the above equation impliesΘ0(k,0) =12Φ(k,0). (1.227)52Also, using the Boltzmann equation (1.222) on superhorizon scales, Θ˙0 =−Φ˙,and the initial conditions (1.227), and using the fact that the potential decays afterhorizon crossing, we obtain the approximationΘ0(k,η)≈ 12Φ(k,0)−Φ(k,η)+Φ(k,0) =32Φ(k,0), (1.228)which has been shown to hold in the limit that Π= 0 ([50]), ([49]), ([54]).In the next section on large scale solutions, we will show that before the epochof matter radiation equality, when anisotropic stress is accounted for, that we ac-tually have (1.245)−Ψ = (1+ 25 fν)−1Φ. Therefore, following the same steps asabove, we have, correcting for anisotropic stress, ([50])Θ0(k,0) =−12Ψ(k,0) =12(1+25fν)−1Φ(k,0). (1.229)It then turns out that the approximation (1.228) can be replaced to good accuracyby ([50])Θ0(k,0)≈ 32(1+25fν)−1Φ(k,0). (1.230)This will be important for determining the small-scale solutions which enter thehorizon before matter radiation equality.Velocities for matter in the early universe can be entirely specified in terms ofthe initial conditions for the potential using the Euler equation (1.208), neglect-ing coupling terms above the horizon. The Euler equations are then identical forbaryons and cold dark matter, and we have that velocities started off behaving as([34])Vb =Vc =kΨ2aH. (1.231)Neglecting anisotropic stress, and assuming the velocities started off zero afterinflation for our adiabatic initial conditions, since (1.231) also satisfies the Eulerequation for photons and neutrinos (1.184) we also haveVν =Vγ =Vm. (1.232)This initial velocity condition will be used when solving for the variance of the53potential at horizon crossing in terms of variables describing the scalar field frominflation. In the next section, we will solve for the matter power spectrum in termsof the initial variance of the potential, and show that correction factors to the powerspectrum are introduced when anisotropic stress is perturbatively accounted for.The initial condition in terms of the inflaton field variables that we will later use,however, as shown in ([34]), has not accounted for anisotropic stress.1.8.2 Large Scale SolutionsWe now explain the large scale solutions for the Einstein equations followingthe initial conditions, using ([45]), ([49]), ([48]) (see also ([34])); the results inthe absence of anisotropic stress are based on the original solution in ([53]). Asexplained above (see (1.165)), for modes above the horizon at late times, in the totalmatter gauge, we have ∆m → ∆T . Since this second quantity is a gauge invariantvariable, we will use it to define the large-scale power spectrum.It follows from the definition (1.202) thatpTΓT = δ pT − p˙Tρ˙T δρT=∑ipiΓi+(c2i − c2T )δρi.(1.233)We then find that, for a universe containing matter and radiation,ΓT =−431−3wT1+wTS, (1.234)where (see (1.216), (1.217))S≡ δ (nm/nr)nm/nr= ∆m− 34∆r (1.235)is the fractional perturbation in the ratio of the number densities of matter andradiation. Using this, the relations (1.165) and (1.166) and the Einstein equations,we can then write (1.193) and (1.203) as ([45]), ([53])∆˙T −3wT a˙a∆T =−(1+wT )kVT −2a˙awTΠT , (1.236)54V˙T +a˙aVT =43wT(1+wT )2k(∆T − (1−3wT )S)+ kΨ− 23kwT1+wTΠT . (1.237)As explained previously, if we’ve assumed adiabatic initial conditions then S =0. Also, from combining the continuity equations (1.195)-(1.196) for matter andradiation, and using the definition (1.235) and the previous discussion, we have asstated in ([49]), ([45]), S˙= k(Vr−Vm), so S˙= 0 initially. S˙ also has essentially zerotime derivative by the baryon and photon Euler equations, so we will set S = 0.We will neglect anisotropic stress as a first approximation. Then from (1.214),Ψ≈−Φ, and Φ is related to ∆T by the Poisson equation (1.213), where a2ρT canbe written in terms of a˙a by the Friedmann equation. We also use the fact that duringmatter or radiation domination, we havea˙a=(1+a)1/2akeq (1.238)where keq = aeqHeq, is given byk2eq =2ΩmH20aeq. (1.239)is the scale that enters the horizon at matter radiation equality. The equations(1.236) and (1.237) can then be combined to form a second order differential equa-tion for ∆T . The solution for the growing mode can be written in terms of y ≡ aaeqas ([53]), ([45])∆T = A(k)UG(y) = A(k)(y3+29y2− 89y− 169+169√y+1)1y(1+ y)(1.240)which reduces to ([49])∆T = A(k)109y2, y << 1= A(k)y, y >> 1,(1.241)where we will later find A(k) by initial conditions.Next it is possible to account for anisotropic stress, as in ([48]), ([45]). Beforerecombination, due to Compton scattering, the photon distribution became more55isotropic than neutrinos, so the main contribution to the anisotropic stress camefrom neutrinos. Recalling the equation for anisotropic stress (1.209), we then haveΠT ≈ 125 fνN2, (1.242)where as mentioned previously fν =ρνρν+ργ ≈ .405. From the neutrino Boltzmannequation for the quadrupole, neglecting higher moments on superhorizon scales,and recalling the interpretation of the dipole as a velocity and the initial conditions(1.231), we have˙N2 =23kN1 =23kVν =23kVT . (1.243)To correct for anisotropic stress, one can use an iterative approach. The first stepis to substitute the growing mode (1.240) into (1.236) to find the zeroth order so-lution for VT in the absence of stress. With this, one solves for N2 in (1.243) andhence ΠT through (1.242), then substitutes ΠT back into the equations (1.236) and(1.237), to re-solve for ∆T . The corrected solution picks up extra factors dependingon fν as∆T (a) = (1+25fν)A(k)109y2, y << 1,= (1+415fν)A(k)y, y >> 1.(1.244)Then substituting this result into (1.236) and (1.237) and using the expressionfound for ΠT , one can find an expression for VT and Ψ in terms of ∆T , whereas ∆Tcan be related to Φ using the Poisson equation. From this we findΨ=−(1+ 25fν)−1Φ, y << 1=−Φ, y >> 1.(1.245)Next we must find the normalization factor A(k). From the Poisson equation(1.213), we have in the radiation dominated limit910(1+25fν)−1Φ(k,η) =4piGρT a4k2a2eqA(k), (1.246)56and in the matter dominated limit,(1+415fν)−1Φ(k,η) =4piGρT a3k2aeqA(k), (1.247)from which we have, as mentioned in ([49]), that the potentials are then approx-imately constant values at the respective times. Φ(k,η) in (1.246) can then beapproximated by Φ(k,0), whereas at aeq, we must have the right-hand sides of theabove two equations equal, soΦ(k,ηeq) =910(1+25fν)−1(1+415fν)Φ(k,0). (1.248)Since it stays constant after radiation is negligible it may continue to be approx-imated during matter domination by this value. Therefore, substituting this into(1.247), and using the fact that during matter domination, 4piGρa3 = 32 H20Ωm, wehave, as in ([45])A(k) =35ΩmH20(1+25fν)−1k2aeq =65(kkeq)2(1+25fν)−1Φ(k,0). (1.249)As explained in ([45]), at late times a>> aeq, and for superhorizon modes, thecontinuity and Euler equations can be combined to give∆¨T +a˙a∆˙T = 4piGρT a2∆T . (1.250)By the Friedmann equations, H obeys this same equation. Therefore, we can sub-stitute a solution ∆T ∝ D(a) ∝ H to obtainD(a) ∝ H∫ a0da′(a′H(a′))3. (1.251)The constant of proportionality will be fixed to initial conditions, where from(1.244) we must require thatD(a) = a (1.252)57in the matter dominated era, consistent with the notation of ([34]). This givesD(a) =5ΩmH(a)2∫ a0da′(a′H(a′))3. (1.253)We then make the replacement a→ D in the general solution. We therefore have,using (1.244) and (1.249),∆T =35H20Ωm(1+415fν)(1+25fν)−1k2Φ(k,0)D(a)=65(kkeq)2(1+415fν)(1+25fν)−1Φ(k,0)D(a)(1.254)at late times after matter domination.1.8.3 Cold Dark Matter SolutionWe now explain the cold dark matter solution, following ([45]), ([50]). In theradiation dominated era we obtain as a solution of the equation (1.210)δc = 3Θ0(k,0)+∫ η0(lna′− lna)a′a˙′(k2Ψ+3Φ¨)dη ′. (1.255)After the potentials decay following entry into the horizon we haveδc ≈ I1Φ(k,0)ln(I2aaH), aH << a << aeq, (1.256)where aH is the scale factor at horizon crossing, which can be found to beaH(k)aeq=1+√1+8(k/keq)24(k/keq)2≈√22keqk, k >> keq, (1.257)andI1 = 9.11(1+0.128 fν +0.029 f 2ν ), I2 = 0.594(1−0.631 fν +0.284 f 2ν ). (1.258)This radiation dominated solution must then be joined onto the growing mode inthe matter dominated era. After radiation domination but before the drag epoch,58again using y ≡ aaeq , the evolution of the cold dark matter can be rewritten as theequationd2dy2δc+(2+3y)2y(1+ y)ddyδc =32y(1+ y)ΩcΩmδc. (1.259)The solution isδc(k,η) = I1Φ(k,0)[A1U1(η)+A2U2(η)], (1.260)where U1 and U2 are written in terms of hypergeometric functions F byUi = (1+ y)−αiF(αi,αi+12,2αi+12;11+ y), (1.261)whereαi =1±√1+24Ωc/Ωm4, (1.262)andA1 =−Γ(α1)Γ(α1+ 12)Γ(2α1+ 12)(ψ(α1)+ψ(α1+12)−ψ(α2)−ψ(α2+ 12))×(ln(I2aeqaH)+2ψ(1)−ψ(α2)−ψ(α2+ 12)),(1.263)whereψ(x) = Γ′(x)/Γ(x), (1.264)and with a corresponding expression for A2 with 1→ Small Scale SolutionsWe now explain the solutions to the Boltzmann equations and the baryon frac-tional overdensities and velocities on small scales, as described in ([50]), ([34]),([70]) and references therein, where the last reference is the original paper withthese calculations. See also ([65]). The tight coupling approximation assumes thatcollisions between photons and electrons are rapid enough such that τ˙ = neσT a>>1, as was true at early times when the free electron fraction was higher. If we as-sume that there exists some l > 2 beyond which the Boltzmann hierarchy may betruncated, we can rewrite the photon Boltzmann equations with theΘl terms on the59right-hand side, for this l, ignoring Θl+1, as (see [34])Θl = τ˙−1(−Θ˙l + k(l2l−1Θl−1))≈ 12kτ˙Θl−1,(1.265)where in the second equality we have neglected the first term on the right-hand sidebecause by iteration it is to second order in τ˙−1. Therefore, the main contributionto moments l > 2 is proportional to successively higher powers of kτ˙−1. Since τ˙is large, it is therefore a good approximation to only consider the first 3 moments,Θ0, Θ1, and Θ2.In the following discussion, based on that given in ([50]) and references therein,we will explain the dispersion relation for the Boltzmann equations. One can showfrom the polarization monopole and quadrupole equations (1.187) and (1.189) thatto zeroth order in τ˙−1,ΘQ2 =ΘQ0 =14Θ2. (1.266)Substituting this into the temperature quadrupole equation (1.185) neglecting highermoments we haveΘ2 = τ˙−1 f−12(23kΘ1− Θ˙2)≈ τ˙−1 f−1223kΘ1 (1.267)wheref2 =34, (1.268)and in the second equality we have neglected the Θ˙2 term because by iterating theright-hand side it is to second order in τ˙−1. This shows that the quadrupole scaleswith the dipole by an additional power of τ˙−1 and may also be neglected, leavingonly the monopole and dipole, as a first approximation.We can now find the dispersion relation assuming that ω is large enough suchthat Φ, Ψ and R can be treated as constants, and effects of expansion can be ne-glected. We will assume a solution of the form Θ1 ∝ ei∫dηω . From the baryonEuler equation (1.208), recalling as established earlier the interpretation Vγ = Θ1,60we haveVb =Θ1− τ˙−1R(V˙b+ a˙aΘ1− kΨ). (1.269)Therefore iterating (1.269) twice with these approximations, we obtainVb ≈Θ1− τ˙−1R(iωΘ1− kΨ)− τ˙−2R2ω2Θ1. (1.270)Substituting this and the quadrupole expression (1.267) into the photon dipoleequation (1.184) and keeping terms to first order in τ˙−1 givesiω(1+R)Θ1 = k(Θ0+(1+R)Ψ)− τ˙−1R2ω2Θ1− 415 τ˙−1 f−12 k2Θ1. (1.271)Differentiating (1.271) and using the photon monopole equation (1.183) to re-late Θ0 and Θ1, ignoring time derivatives of potentials, we obtain(1+R)ω2 =k23+ iτ˙−1ω(R2ω2+415k2 f−12 ). (1.272)Iterating ω2 on the right-hand side in parentheses, and solving the quadratic equa-tion, the dispersion relation to first order in τ˙−1 is approximatelyω =± k√3(1+R)+i6k2τ˙−1(R2(1+R)2+45f−1211+R). (1.273)This assumes ([70]) that the damping term is small compared to the first term.To go farther, it will now be useful to understand the equation for the monopoleto zeroth order in τ˙−1. As explained in ([34]), substituting (1.269) into the photondipole equation (1.184) givesΘ˙1 =− R˙1+RΘ1+11+RkΘ0+ kΨ; (1.274)if we differentiate the photon monopole equation (1.183) and substitute into it(1.274) and (1.183) again, we obtain the forced oscillator equationc2sddη(c−2s Θ˙0)+ c2s k2Θ0 =−k23Ψ− c2sddη(c−2s Φ˙), (1.275)61([34]), ([46]), ([50]), wherec2s =13(1+R)=13(1+ 3Ωb4Ωγ (1+z)). (1.276)As explained in ([50]), this is then a harmonic oscillator with a time vary-ing mass, me f f = (1+R). The dispersion relation for an oscillator with a massvarying slowly over a period of the oscillation is approximately ω = kcs. The en-ergy, E = 12 me f fω2A2, of this harmonic oscillator changes slowly over a period ofthe oscillation, so the oscillator will then have an approximate adiabatic invariant(see ([14])) given by its energy E = 12 me f fω2A2, divided by its frequency ω, fromwhich we conclude A ∝ (1+R)−1/4. We then find the zeroth order solutions arewell approximated byΘˆ0+Ψ≈ Θˆ0 = (1+R)−1/4CAcos(krs) (1.277)δˆb = 3(1+R)−1/4CAcos(krs). (1.278)andΘˆ1 =−√3(1+R)−3/4CAsin(krs), (1.279)wherers(z) =∫ ∞zcsdz′H(z′) , (1.280)and in the first equation we have used the fact that Ψ decays after it crosses thehorizon. The coefficient CA, referring to the assumption of adiabatic initial con-ditions, is determined from the initial solution for the mode before it entered thehorizon. From (1.230) and the fact that at early times, R→ 0, we have that this isgiven byCA(k)≈Θ0(k,0)|adiabatic= 32(1+25fν)−1Φ(k,0). (1.281)These oscillating solutions then represent the effect of baryon acoustic oscillations,with sound horizon given by (1.280).Next, we notice that that the first order dispersion relation (1.273) adds a damp-ing factor from the dispersion relation of the zeroth order solutions, ω = csk, wherein finding these dispersion relations in both cases we have assumed frequencies are62high enough such that the effects of expansion are small in comparison. It is there-fore possible to have the first order solutions reduce to the zeroth order solutions inthe absence of the new corrections by writing ([50])Θ0 = e−(k/kD)2(1+R)−1/4CAcos(krs), (1.282)δb = e−(k/kD)23(1+R)−1/4CAcos(krs), (1.283)Θ1 =−e−(k/kD)2√3(1+R)−3/4CAsin(krs), (1.284)where (see also ([34]),k−2D ≡16∫ η0dη ′τ˙−1(R2+ 45 f−12 (1+R)(1+R)2). (1.285)We can confirm that these are approximate solutions to the small scale equations,in which large k makes ω large, when terms to first order in τ˙−1 are included (see([70])).The general solution to the baryon Euler equation is given byaVb =∫ η0dη ′a′(τ˙dΘ1+ kΨ)e−τd , (1.286)where as defined in ([50]),τ˙d ≡ τ˙R =neσT aR, (1.287)in this case with τ the Compton optical depth, and where τd(η ′) =∫ η0η ′neσT a′′dη ′′Rin the exponent is the optical depth from the present to an earlier time η ′ duringrecombination. The first term of (1.286) is the effect of the baryons following thepath of the photons due to Compton drag, an effect which diminishes with time asthe free electron fraction decreases. The second is the effect of infall into potentialwells, which will start to take effect around the time when τd has decreased tounity. Therefore, we can define this transition redshift, zd byτd(zd) = 1 (1.288)as the time when the baryons are released from the drag of the photons and begin to63behave identically to the cold dark matter. We use the fitting formula for zd givenby ([35])zd = 1291(Ωmh2)0.2511+0.659(Ωmh2)0.828(1+b1(Ωbh2)b2), (1.289)whereb1 = 0.313(Ωmh2)−0.419(1+0.607(Ωmh2)0.674), b2 = 0.238(Ωmh2)0.223, (1.290)with Ωm =Ωb+Ωc.Immediately following the drag epoch, we then have the approximate initialsmall scale velocity solution, given in ([50]), taking into account damping,Vb(k,ηd) = Θˆ1(k,ηd)Db(k) =−Db(k)√3(1+R)−3/4CAsin(krs(ηd)), (1.291)where Db(k) is given byDb(k) =∫ η00dηVb(η)e−(k/kD)2, (1.292)where Vb is the normalized drag visibility function,Vb =aτ˙de−τd∫ η00 dηaτ˙de−τd. (1.293)In making this approximation, one has assumed that Θˆ1 varies slowly comparedto the rest of the integral, and used the fact that τ˙de−τd is approximately a deltafunction, peaked at ηd , in comparison with variations of η .As explained previously in the section on baryon acoustic oscillations (see([37]), ([34]), ([50]), ([3]), ([45])), following the drag epoch, recalling that theoverdensities of baryons and cold dark matter were initially the same, the overden-sities of baryons have traveled a comoving distance ofrs(zd) =∫ ∞zdcsdzH(z)(1.294)away from the center of the matter perturbations. The cold dark matter did notparticipate in the acoustic oscillations, but stayed near the centers of the pertur-64bations. This then made a positive contribution to the correlation function of thematter fractional overdensities at a comoving separation of rs(zd)≈ 150 MpC. Thiscontribution remained following the drag epoch, when the baryons, like the colddark matter, were only affected by gravity. This allows BAO’s to serve as a stan-dard ruler, subject to noise in the correlation function from cosmic variance set upby inflation.Following the drag epoch, coupling to photons ends and the baryons and thecold dark matter obey the same equations. The solution for each is then found bymatching to the initial solutions for δ and V at the drag epoch, where we can usethe fact that ([8]) δ˙ =−kV from (1.195) neglecting Φ˙. The result is ([50])δb(k,η) = (G1(ηd)δb(k,ηd)+G2(ηd)kVb(k,ηd))D1(η), (1.295)with an identical equation for cold dark matter, whereG1(η) =D˙2D1D˙2− D˙1D2, G2(η) =D2D1D˙2− D˙1D2, (1.296)where D1 and D2 are respectively the growing and decaying dark matter solutionsfound in ([69])D1 =23+ y, D2 =158(2+3y)ln((1+ y)1/2+1(1+ y)1/2−1)− 454(1+ y)1/2. (1.297)As explained in ([34]), at late times, for subhorizon modes, the baryons andcold dark matter obey the same equations, given by (1.195), (1.208) without thecoupling term, and (1.213) in the limit where kη >> 1. We write these asδ˙m =−kVm−3Φ˙, (1.298)V˙m =− a˙aVm+ kΦ, (1.299)k2Φ= 4piGa2ρmδm. (1.300)These can be combined to obtain an identical equation for δm as for ∆T in (1.250).Consistent with the notation of ([45]), from the form of D1 in (1.297), we see that65the initial conditions for the solution, taken for a >> aeq, is that it is proportionalto a, as for the large scale solution. Therefore, following the same argument as forthe large scale solutions, we replace a→ D, given by (1.253), in (1.295) and thecorresponding equation for cold dark matter, as done in ([45]).1.8.5 The Power Spectrum and Transfer FunctionIt is convenient to define the matter power spectrum in terms of the variance offractional overdensities in Fourier space as ([34])〈δ ∗m(k,z)δm(k′,z′)〉 ≡ (2pi)3δ 3(~k−~k′)Pm(k,z,z′). (1.301)For the following discussion, see ([45]), ([50]), ([35]). Given the expressions forthe total matter overdensity, which can be found from the results of the previoussection, and are proportional to Φ(k,0), we can now define the power spectrum ofthe overdensities of matter at late times for an arbitrary scale in terms of the largescale solution (1.254) asPm(k,a) =925k4Ω2mH40(1+415fν)2(1+25fν)−2PΦ(k,0)T 2(k)D(a)2. (1.302)Here, 〈Φ∗(k,0)Φ(k′,0)〉 = (2pi)3δ 3(~k−~k′)PΦ(k,0) (see ([34])) and will be dis-cussed in the section on initial conditions from inflation. T (k) is the transfer func-tion, by construction unity on large scales, defined as ([35])T (k)≡ δ (k,z = 0)δ (0,z = 0)δ (0,z = ∞)δ (k,z = ∞). (1.303)The transfer function is defined as a weighted sum of contributions for baryons andcold dark matter, asT (k) =ΩbΩmTb(k)+ΩcΩmTc(k). (1.304)We can substitute the adiabatic solution for baryon fractional overdensity (1.283)and velocity (1.291) into (1.295) to obtain the expression for δb. Then using thisand the limit of the large scale solution (1.254) in the definition (1.303), the transfer66function for baryons isTb(k)=154(1+415fν)−1(keqk)2Db(k)(1+R)−1/4(cos(krs)− D2D˙2kcssin(krs))G1 |η=ηd ,(1.305)where the first term is the density contribution and the second the velocity contri-bution.Similarly, the dark matter contribution to the transfer function can be foundfrom the solution to δm to evaluate (1.295) applied to dark matter, to findTc(k)= I1(1+25fν)(1+415fν)−1 56(keqk)2(G1[A1U1+A2U2]−G2[A1U˙1+A2U˙2]) |η=ηd .(1.306)1.8.6 Initial Conditions for the PotentialFor the following explanation of inflation, first proposed in ([39]), we will relymostly on ([34]), with a few other citations given as necessary. See also the refer-ences therein for the other early papers on inflation.Use of Inflation to Solve the Horizon ProblemIf we assume the same expression for H before z∗, the redshift of the CMBwhen the universe first became transparent to light, as after it, we see that thecomoving Hubble radius has been increasing since the big bang, and one can cal-culate that the maximum comoving distance that light could have traveled from thebeginning of the universe to the time of transparency is much less than the distancethat photons from the CMB reaching us from widely separated patches of the skywere when the light was emitted. However, the photons of the CMB have approxi-mately the same temperature. This is called the horizon problem. To explain this,physicists appeal to new physics that operated before the universe became transpar-ent to light. According to this suggestion, the comoving Hubble radius was oncelarge, then decreased in a period called inflation, allowing the uniform temperaturephotons we observe to have once been in contact with each other. If we assumethat the particles that we see today were in causal contact with each other beforeinflation, then the comoving Hubble radius at that time must have been greater than67it is today (see earlier discussion):1aiHi>1H0. (1.307)For a model of inflation, we would then choose values of ae, ai, and H thatwould satisfy this condition. It is usually assumed that inflation ended at aroundTe = 1015 GeV before other known physical processes. This can be related to ascale factor using the relation T = Tγ0a where Tγ0 is the observed temperature of theCMB, to obtainae =Tγ0Te=Tγ01015GeV≈ 10−28. (1.308)This can be related to HeH0 by substituting this value for ae the Friedmann equation,and we can then find a ratio of the comoving Hubble radii before and after inflationusing (1.307). To satisfy the resulting extremely large ratio, one generally assumesthat H is approximately constant during inflation, so that the ratio is simply theratio of the scale factors, and thata = aeeH(t−te), t < te. (1.309)From this, one concludes thatH(ti− te)≈ 60, (1.310)so the scale factor increases by about 60 orders of magnitude during inflation.A scalar field is capable of making H approximately constant and thereforesolving the horizon problem in this way. A classical scalar field φ in curved space-time has an action ([12])S =∫dη d3x√−|g|(12gµν∂φ∂xµ∂φ∂xν−V (φ))(1.311)and a stress energy tensor ([34])Tαβ = gαν ∂φ∂xν∂φ∂xβ−gαβ(12gµν∂φ∂xµ∂φ∂xν+V (φ))(1.312)68which is a divergenceless quantity as a consequence of the field’s equations ofmotion ([86]). When quantum mechanics and general relativity are attempted tobe combined (see ([82])), classical field values are replaced by field operators φˆ ,acting on a Hilbert space, and the stress energy tensor will be a function of theseoperators. However, it is assumed that on large scales, the field can be treatedclassically, and only fluctuations are treated quantum mechanically. Therefore, wewriteφˆ = φ (0)+ ˆδφ , (1.313)andTˆ µν = Tµ(0)ν (φ (0))+ ˆδTµν (φ(0), ˆδφ). (1.314)The field is assumed to zeroth order in perturbations to be homogeneous. There-fore, the homogeneous part of the stress tensor isρ =−T 0(0)0 =12(dφ (0)dt)2+V (φ (0)), (1.315)P = T i(0)i =12(dφ (0)dt)2−V (φ (0)). (1.316)The Einstein equations for a perfect fluid (1.50) and (1.51) gived2a/dt2a=−4piG3(ρ+3P). (1.317)Therefore, to produce a positive acceleration of the scale factor a, P must be neg-ative, and therefore the kinetic energy term must be less than the potential energyterm. This implies that φ (0) is trapped at a value which is a false minimum ofthe potential V (φ (0)). Inflation ends when φ (0) reaches its true ground state. Weassume that it is slowly rolling towards its ground state so that the density is ap-proximately constant, allowing H to be constant. Two quantities that characterizethe slow roll areε =116piG(V ′V)2(1.318)69andδ = ε− 18piGV ′′V, (1.319)where the primes are derivatives with respect to φ (0). The zeroth order equation ofthe field is¨φ (0)+2aHφ˙ (0)+a2V ′ = 0, (1.320)which follows from the time-time Einstein equation.Variance of the Field PerturbationsNext we solve for the power spectrum of the perturbations to the inflaton fieldδφ . When the operator ˆδφ is found as a function of time, given an initial state,the expectation value of the fluctuation can be predicted for an ensemble of uni-verses. Using the gauge with spatially flat slicing, the zeroth component of thedivergenceless of the perturbed stress energy tensor equation becomes, in Fourierspace,¨ˆδφk+2aH˙ˆδφ + k2 ˆδφk = 0. (1.321)Making the change of variablesˆ˜δφk = a ˆδφk (1.322)we have¨˜ˆδφk+(k2− 2η2)ˆ˜δφk = 0, (1.323)which is the equation of a harmonic oscillator.Writing out the perturbed part of the action (1.311) in terms of ˜δφ in realspace, our conjugate variable to ˜δφ is ∂L∂ ( ˙˜δφ)= ˙˜δφ(η ,~x) ([12]). For the followingdiscussion, see ([12]), ([20]), ([15]), ([34]), but note some differing conventions;our convention will be consistent with ([20]). The conjugate variables must obeythe canonical commutation relations[ ˆ˜δφ(η ,~x),˙˜ˆδφ(η ,~x′)] = iδ 3(~x−~x′). (1.324)70In Fourier space, we haveˆ˜δφ(η ,~x) =∫ d3k(2pi)3ei~k·~x ˆ˜δφk. (1.325)By expanding ˆ˜δφ(η ,~x) and˙˜ˆδφ(η ,~x′) on the left-hand side of (1.324) in terms oftheir Fourier coefficients as above, we then find that (1.324) will be obeyed if[ ˆ˜δφk,˙˜ˆδφ †k′ ] = i(2pi)3δ 3(~k−~k′). (1.326)The solution to (1.323) is thenˆ˜δφk = aˆke−ikη√2k(1− ikη)+ aˆ†−keikη√2k(1+ikη)(1.327)where aˆk and aˆ†−k are time independent operators. If aˆk and aˆ†k obey commutationrelations[aˆk, aˆ†k′ ] = (2pi)3δ 3(~k−~k′), (1.328)with all others zero, we can check that the canonical commutation relations (1.326)will be obeyed. Also, by writing the Hamiltonian, given in terms of the field op-erators, now in terms of the ak and a†k operators, and requiring that the vacuum bethe lowest energy eigenstate, we require that ak|0〉= 0 for all k. ([15])We can use the solution for the field operator above to solve for the energyeigenstates in the | ˜δφk〉 basis to determine the probability that the field takes thevalues ˜δφk. We will assume that the perturbation to the inflaton field has decayedto the ground state by the time it leaves the horizon (for a discussion of this as-sumption, see ([86])). Given an ensemble of universes, the statistics due to ˜δφ canthen be predicted. Using (1.322) and (1.327), the fluctuations have a mean of zero.Computing the variance, diagonal in k space, in terms of the definition of the power71spectrum (1.301), we havePδφ (k) =1(2pi)31a2〈0| ˆ˜δφ †k ˆ˜δφk|0〉=1a212k3η2(1+ k2η2)=H22k3(1+k2a2H2).(1.329)During inflation, we redefine the conformal time to be a negative quantityη =∫ ttedta(t), t < te (1.330)which approaches zero closer to the end of inflation. In the limit that k|η |→ 0,or aH >> k, when the mode is far outside the horizon, this approaches a constantvalue H22k3 . From the equation above, the power spectrum at the point the modeleaves the horizon can be approximated by, up to a factor of two,Pδφ |aH=k=H22k3. (1.331)The Power Spectrum Following InflationIn this section, continuing the discussion in ([34]), we obtain an expression forthe power spectrum due to the metric perturbation Φ initially after inflation. Inthe spatially flat slicing gauge, at horizon crossing, evaluating the gauge invariantquantity v (1.167) using the expression for the stress energy tensor, for δT 0i , wehavev = ikB− ikφ˙(0)δφ(ρ+P)a2= ikB− ikδφφ˙ (0).(1.332)The gauge invariant quantity (1.168), also in the spatially flat slicing gauge, is givenbyΦH = aHB. (1.333)72Then using the gauge invariant linear combinationζ ≡−ΦH − iaHk v, (1.334)in the spatially flat slicing gauge, at horizon crossing, we haveζ =−aHδφ˙φ (0)|aH=k. (1.335)ζ can be shown to be conserved outside the horizon ([34]) from the vanishingdivergence of the stress energy tensor. Therefore, we havePζ =(aHφ˙ (0))2Pδφ=4piGεPδφ .(1.336)In the conformal Newtonian gauge,ζ =− ikiδT0i aHk2(ρ+P)−Ψ. (1.337)We can find an expression for this conserved quantity after inflation. On the otherhand, we know the stress energy tensor after inflation during the period of radiationdomination,ikiδT 0i = k(ρ+P)Vγ , (1.338)so evaluating (1.337) after inflation, we haveζ =−aHΘ1k−Ψ=−32Ψ(1.339)after using the initial condition for the velocity (1.231), implyingPΦ =49Pζ . (1.340)Then substituting (1.331) and (1.336) into (1.340), we have an expression for PΦ73at horizon crossing:PΦ =8piGH29εk3|aH=k . (1.341)We can rewrite this as ([34])PΦ(k) =50pi29k3(kH0)n−1δ 2H(ΩmD(a = 1))2, (1.342)where n is the scalar spectral index, given byn−1 =−4ε−2δ , (1.343)and δH is referred to as the scalar amplitude. In terms of this, the matter powerspectrum (1.302) becomesP(k,a) = 2pi2δ 2HknHn+30(1+415fν)2(1+25fν)−2T 2(k)(D(a)D(a = 1))2. (1.344)1.9 The Hydrogen Brightness Temperature PowerSpectrumThe spectral radiance, or energy per time per solid angle per unit emitting areaper unit frequency, of the neutral hydrogen signal will be proportional to its densityin galaxies. As in ([80]), we will denote this quantity by I(ν ,T ). As explained in([11]), at low temperatures, where hν << kBT, now no longer using natural units,the spectral radiance for a blackbody, given by Planck’s law,I(ν ,T ) =2hν3c21ehν/kBT −1 , (1.345)reduces to the Raleigh-Jeans limitI(ν ,T ) = 2kTν2/c2. (1.346)The spectral radiance of the neutral hydrogen, a non-perfect emitter with emissivityε, is given by ε times this, which may then be thought of as a blackbody spectrum74with an effective, or brightness, temperatureTb = εT. (1.347)As explained in ([25]), the hydrogen brightness temperature can be obtained fromthe measured intensity of the emitted neutral hydrogen signal byTb(ν , nˆ) =c2I(ν , nˆ)2ν2kB. (1.348)Since I(ν , nˆ) is proportional to the hydrogen densities in galaxies, we can writethe total observed brightness temperature asTb(z, nˆ) = T¯b(z)(1+δHI(nˆ,z)), (1.349)where δHI is the fractional overdensity of neutral hydrogen. The signal, or fluctua-tion, in the brightness temperature field is thenδTb(z, nˆ) = T¯b(z)δHI(nˆ,z). (1.350)Here,δHI = bδm, (1.351)where b is the bias giving the ratio of neutral hydrogen to matter fractional over-densities.We will assume that b= 1. Then, using (1.350) and the definition of the matterpower spectrum (1.301), it follows that the brightness temperature power spectrumis given by〈 ˜δT ∗b(~k,z) ˜δT b(~k′,z′)〉= (2pi)3δ 3(~k−~k′)PTb(k,z,z′) (1.352)where ([80])PTb(k,z,z′) = T¯b(z)T¯b(z′)Pm(k,z,z′). (1.353)75The mean brightness temperatures are given by ([31]), ([80]), ([19])T¯b(z) = 0.3(ΩHI√0.29√2.510−3)((1+ z)2H(z)/H0)(1.354)in mK, where ΩHI is the comoving mass density of neutral hydrogen in units ofρcr. For this, we use the value 5×10−4 given in ([63]). See also ([40]), ([90]) fordiscussions of the brightness temperature.Also, in (1.351), we have neglected redshift space distortions (see ([34]) for aderivation of the following expressions), where b would be replaced by (b+ fµ2),wheref (a) =aD1dDda, (1.355)and µ = zˆ · kˆ, where zˆ is the direction along the line of sight. To measure red-shift space distortions, we would measure the power spectrum P(k) and the powerspectrum in redshift space Ps(k) and use the formulaPs(k) = P(k)(1+βµ2k )2, (1.356)where β = fb . From f (a), we can solve for the growth factor, D. We can then com-pare this with the growth factor determined from the parameters in the expansionhistory, (1.253). This would then enable us to test the equations of general relativityand should be taken into account in the future.1.10 MapmakingWe will now describe the mapmaking process that CHIME will use, as ex-plained in ([80]) and ([79]). The CHIME telescope will focus radiation from thesky through five separate 20 m by 100 m cylinders, each containing an array ofdual polarization feeds where the focused data is collected ([67]). The electricfield density received from the sky within a frequency interval dν and solid angleinterval dnˆ is given by ([79])d~E = (µ0c)1/2~ε(nˆ,ν)d2nˆdν . (1.357)76We will define our coordinate system with an angle φ relative to the fixed sky asthe earth turns, with basis vectors θˆ and φˆ . For a feed located at position ~ri, theoutput in response to the electric field from the sky is given byFi(φ ,ν) =∫d2nˆAai (nˆ,φ)√Ωiεa(nˆ,φ)e2pi nˆ·~ri/λ , (1.358)where Ωi =∫d2nˆ|~Ai(nˆ,φ)|2, so that Aai (nˆ,φ)√Ωiis a normalized antenna reception pat-tern in the direction of coordinate a orthogonal to the line of sight.At each point φ of the rotation, the telescope measures the correlation, aver-aged over all incoming directions nˆ, between the signals of the two feeds i andj, described by the visibility function Vi j. Allowing for noise in the signal, this isgiven as a sum of a signal and a noise term ni j. The correlation between the electricfield densities in the polarization directions of the coordinate system is determinedfrom〈εa(nˆ,ν)ε∗b (nˆ′,ν ′)〉=2kBλ 2δ (nˆ− nˆ′)δ (ν−ν ′)(PT T (nˆ)+PQQ(nˆ)+PUU(nˆ)+PVV (nˆ)),(1.359)wherePT =12(1 00 1), PQ =12(1 00 −1), PU =12(0 11 0), PV =12(0 −ii 0),(1.360)and T, Q,U, and V are polarization brightness temperatures, with T correspondingto the expression Tb, (1.348). In terms of Fi, the visibilities can then be defined asVi j =λ 2kB〈FiF∗j 〉=∫d2nˆ(BTi jT +BQi jQ+BUi jU +BVi jV )+ni j, (1.361)whereBXi j = 2Aai A∗bj√ΩiΩ jPXabe2piinˆ·(~ri−~r j)/λ . (1.362)This ensures that for a feed i, we have Vii = Tb.However, we are only interested in the T contribution as this determines theneutral hydrogen signal. From here on, we will assume that this contribution is77already separated out in the calculation of the visibility, and will follow the discus-sion in ([80]); the full treatment is covered in ([79]). In this case we then have, interms of the brightness temperature (1.348), for each wavelength λ ,Vi j(φ ,λ ) =∫d2nˆBi j(nˆ,φ ,λ )Tb(nˆ,λ )+ni j(φ ,λ ) (1.363)whereBi j(nˆ,φ ,λ ) =1√ΩiΩ jAi(nˆ,φ)A∗j(nˆ,φ)e2piinˆ·(~ri−~r j)(φ)/λ , (1.364)where Ai = |~Ai|.For a given ν = c/λ , the total visibility function Vi j is periodic in φ so must bewritten in terms of nonzero Fourier modes eimφ for integer m. Therefore, its Fouriercomponents are given byV i jm =∫ dφ2piVi j(φ)e−imφ . (1.365)ExpandingBi j(nˆ,φ) =∑lmBi jlm(φ)Y∗lm(nˆ) (1.366)andTb(nˆ) =∑lm′almYlm(nˆ), (1.367)we then haveV i jm =∑lm′∫ dφ2piBi jlm′(φ)alm′e−imφ +ni jm. (1.368)To obtain the transformed coefficients Bi jlm(φ) with respect to the earth on its rota-tion by φ from φ = 0, we replace Y ∗lm(nˆ) with Y∗lm(θ(nˆ),φ(nˆ)− φ) and rewrite interms of the original basis Y ∗lm(nˆ). Using the fact that Y∗lm has φ dependence e−imφ ,we then haveBi jlm′(φ) = Bi jlm′(φ = 0)eimφ . (1.369)78Therefore, defining Bi jlm(φ = 0)≡ Bi jlm, (1.368) becomesV i jm =∑lBi jlmalm+ni jm. (1.370)Using the fact that for a real field, alm = a∗l−m,V i j∗−m =∑l(−1)mBi j∗l,−malm+ni j−m. (1.371)We then define for m≥ 0Bi j+lm = Bi jlm, ni j+m = ni jmBi j−lm = (−1)mBi j∗l−m, ni j−m = ni j∗−m(1.372)with Bi j−l0 = ni j−0 = 0, so we can write Vi jm in terms of only the alm with m ≥ 0,obtainingV i j±m =∑lBi j±lm alm+ni j±m . (1.373)Combining the indices ± and feed pairs i, j into a single index α, and accountingfor the different frequencies ν , this can be written asv = Ba+n (1.374)where(vm)αν =Vανm , (1.375)(Bm)αν ,lν ′ = Bανlm δνν ′, (1.376)and(am)lν = almν (1.377)for each individual m, and similar for the noise term.The data consists of a combination of foreground, noise, and signal, with thenoise, foreground, and signal covariance matrices assumed known, with modelsfor the foregrounds given in ([80]). Using Karhunen-Loe`ve (KL) transforms (see([34]), ([80]) for a discussion), one identifies and removes the modes of the data79with a low signal to foreground ratio, followed by the modes that have a low signalto foreground plus instrumental noise ratio. Finally, to make the map, the almv ofthe signal is given by solving for the maximum likelihood a of the brightness tem-perature given the data for the cleaned visibilities vclean, according to the equationL(a | vclean) = 1| piN |e− 12 (vclean−Ba)†N−1(vclean−Ba) (1.378)which has the solutiona = (N−1/2B)†N−1/2vclean. (1.379)80Chapter 2Implementation2.1 Determining Cosmological Parameters From MapsFor this discussion, we will use ([33]). Since it is a function on a sphere, thebrightness temperature anisotropy map at a given frequency can be expanded inspherical harmonics asδTb(ν ,θ ,φ) = Σ∞l=0Σlm=−lalmνYlm(θ ,φ). (2.1)The correlation function corresponding to the redshift slices that we observe canbe written asC(nˆ1,ν , nˆ2,ν ′) = Σlm〈a∗lmνY ∗lm(nˆ1,ν)Σl′m′al′m′ν ′Ylm(nˆ2,ν ′)〉= ΣlmΣl′m′〈a∗lmνalmν ′〉Y ∗lm(nˆ1)Yl′m′(nˆ2).(2.2)Using the expressionalmν =∫dnˆYlm(nˆ)δTb(nˆ,ν) =∫dnˆYlm(nˆ)∫ d3k(2pi)3˜δTb(k,ν)eikχ(ν)(kˆ·nˆ), (2.3)we have〈a∗lmνal′m′ν ′〉= 〈∫dnˆ1Y ∗lm(nˆ1)∫ d3k(2pi)3˜δTb∗(k,ν)e−ikχ(ν)(kˆ·nˆ1)∫dnˆ2Ylm(nˆ2)∫ d3k′(2pi)3˜δTb(k′,ν)eik′χ(ν ′)(kˆ′·nˆ2)〉.(2.4)81Using (B.14) and the definition of the Legendre polynomial (B.13), and theorthogonality of the spherical harmonics, we have∫dΩY ∗lm(nˆ)e−ikχ(ν)kˆ·nˆ = 4pi(−i)l jl(kχ(ν))Y ∗lm(kˆ). (2.5)Then using the relations (1.352) and (1.353), and orthogonality of the sphericalharmonics again, we obtain the exact expression〈a∗lmνal′m′ν ′〉= δll′δmm′Clνν ′ (2.6)whereClνν ′ = δll′δmm′2piT¯b(ν)T¯b(ν ′)∫dkk2 jl(kχ(ν)) jl(kχ(ν ′))Pm(k,ν ,ν ′). (2.7)Then, using (2.6) and the definition of Legendre polynomials (B.13) in (2.2), wehaveC(nˆ1,ν , nˆ2,ν ′) = Σl 2l+14pi Clνν ′Pl(nˆ1 · nˆ2). (2.8)We notice that this theoretical correlation function is dependent only on angularseparation θ = cos−1(nˆ1 · nˆ2), so we can express it as C(ν ,ν ′,θ).Now given our data δTb(nˆ,ν) or equivalently the almν , we must find the bestestimate of the cosmology. We can estimate the true covariance matrix Clνν ′ ofthe almν by using the fact that for all m, the almν are independently drawn, due toquantum mechanical fluctuations from inflation, from a distribution with the samecovariance Clνν ′ ([34]) and hence, the best estimate of Clνν ′ is the average productCdatalνν ′ =12l+1Σma∗lmνalmν ′. (2.9)The best estimate of the correlation function is thenCdata(ν ,ν ′,θ) = Σl 2l+14pi Cdatalνν ′Pl(cos(θ)). (2.10)The theoretical power spectrum can be straightforwardly related to the theoretical82correlation function by a Fourier transform, as follows: ([86])C(z, nˆ,z′, nˆ) = 〈∫ d3k(2pi)3˜δT ∗b(k,z)e−i~k·r(z)nˆ∫ d3k′(2pi)3˜δT b(k′,z′)ei~k′·r(z′)nˆ′〉=∫ d3k(2pi)3∫ d3k′(2pi)3〈 ˜δT ∗b(k,z) ˜δT b(k′,z′)〉ei(~k′·r(z′)nˆ′−~k·r(z)nˆ)=∫ d3k(2pi)3PTb(k,z,z′)ei~k·(r(z)nˆ−r(z′)nˆ′)(2.11)using the definition of the brightness temperature power spectrum (1.352), andtherefore the Clνν ′ contain the same information as the power spectrum. We arenot able to take the Fourier transform of the map or correlation function becausewe do not know r(z), so we cannot estimate the power spectrum directly from thedata. Therefore, we must constrain the cosmology using the Cdatalνν ′ or Cdata(ν ,ν ′,θ)as our information.We first calculate the theoretical Clνν ′ from the power spectrum and then useit to generate simulated almν . This allows us to generate the data correspondingprecisely to the position-time coordinates that we observe, namely the points onthe sky with spherical coordinates (r(z)nˆ, t(z)). Generating the data in l space hasthe advantage over angular space in that the covariance matrix is diagonal withrespect to different l.2.2 The Flat Sky ApproximationThe calculation of the Clνν ′ is very computationally expensive to perform forlarge l. Fortunately, for large values of l which span small values of angle acrossthe sky, the flat sky approximation is a good approximation, and agrees with thefull-sky approximation to the one percent level for l > 10 ([33]), ([80]). The flatsky approximation can be derived as follows (([33]), see also ([80]) for the finalform of the expression). We write position vectors on the sky asnˆ = ~m+~θ (2.12)where ~m is a position vector pointing in the center of the field of view, and θ << 1is a two-dimensional vector in the plane tangent to the position vector ~m. Under83the flat sky approximation, we can relate the almν to the Fourier transform of theangular part of δTb(ν , nˆ), given by˜δT b(ν ,~U)≡∫d~θe2pii~U ·θδTb(ν , nˆ) (2.13)where ~U is a two-dimensional vector in the plane of the sky. Given ~U =(UcosφU ,UsinφU)and ~θ = (θcosφ ,θsinφ), using the expansione−2pii~U ·θ = Σm(−i)mJm(2piUθ)eim(φU−φ), (2.14)where Jm is an ordinary Bessel function, and using the small angle approximationYlm(θ ,φ)≈ Jm(lθ)√l2pieimφ , θ → 0, (2.15)(2.13) becomes˜δT b(ν ,~U) =√1UΣm(−i)meimφU∫dθY ∗2piU,m(θ ,φ)δTb(ν , nˆ). (2.16)Then using the approximationalmν ≈∫d~θY ∗lm(θ ,φ)δTb(ν , nˆ) (2.17)under the assumption that l is large, we have˜δT b(ν ,U)≈√1UΣm(−i)meimφU a2piU,m(ν). (2.18)Then taking the variance, usingΣmeim(φU−φU ′) = 2piδ (φU −φU ′), (2.19)and the definition〈a∗lmνalmν ′〉=Clνν ′δl,l′δm,m′, (2.20)84we have〈 ˜δT b(ν ,~U) ˜δT ∗b(ν ′, ~U ′)〉= 2piC2piU(ν ,ν ′)δU,U ′Uδ (φU −φU ′). (2.21)We can simplify the right-hand side by determining that, again using the smallangle expansion (2.15) and the relation (2.19),δ 2(~U− ~U ′) =∫d~θe−2pii(~U−~U ′)·~θ≈∫d~θ√1UU ′Σmm′(−i)m−m′Y ∗2piU,m(θ ,φ)Y2piU ′,m′(θ ,φ)eimφU e−im′φU ′= 2piδU,U ′Uδ (φU −φU ′),(2.22)obtaining〈 ˜δT ∗b(ν ,~U) ˜δT b(ν ′,~U ′)〉 ≈Cl=2piU,ν ,ν ′δ 2(~U− ~U ′). (2.23)With this correspondence established, we can estimate Clνν ′ from the varianceof (2.13), which we now relate to the matter power spectrum. SubstitutingδTb(ν , nˆ) = T¯b(ν)∫ d3k(2pi)3eiχ(ν)(k‖+k⊥˙~θ)δ˜m(ν ,~k) (2.24)into (2.13) and integrating over~k⊥ = 2piχ(ν)(U′i ,U′j) gives˜δT b(ν ,~U) =T¯b(ν)2piχ(ν)2∫ ∞−∞dk‖eik‖χ(ν)δ˜m(k‖mˆ+2pi~Uχ(ν)). (2.25)Computing the variance, we have for l = 2piU, ([33]), ([80]), neglecting red-shift space distortions,Clνν ′ ≈ 1piχ(ν)χ(ν ′)∫ ∞0dk‖cos(k‖∆χ)Pm(k)T¯b(ν)T¯b(ν ′), (2.26)where χ and χ ′ are the comoving distances to redshift z(ν) and z(ν ′), ∆χ = χ−χ ′,85andk =√k2‖+(lχ¯)2(2.27)where χ¯ = χ+χ ′2 .2.3 Calculation of the Matter Power SpectrumIf we assume that Ωνh2 = 0.00064 and Ωk = 0, our present-day model of thepower spectrum depends the equation of state of dark energy, w, and six other in-dependent input cosmological parameters. The first three parameters are the scalaramplitude As, scalar spectral index ns, and reionization optical depth, τ, defined asin ([81]) by τ =∫ η0ηr dη ′neσT a′. In terms of h ≡ H0sMpC100km ≈ 0.7, the other three areeither Ωbh2, Ωch2, ΩΛ, or equivalently, Ωbh2, Ωch2, and h, where ΩΛ is fixed to1−Ωb−Ωc−Ων−Ωk. The fiducial values are given by w=−1,Ωbh2 = 0.02260,Ωch2 = 0.112000, h = 0.7, As = 2.46e−09, ns = 0.96, and τ = 0.09.For our power spectrum, we choose kmax = 100h/MpC. This high value of kwas needed to achieve an approximate convergence from the correlation functionalong the line of sight direction when performing the integrals over the oscilla-tory functions. To test convergence, we took the Fourier transform of this powerspectrum (the correlation function), with kmax ≈ 100h/MpC, as done in ([75]), andfound that the shape stayed constant with increasing k.The power spectrum from CAMB sets the values kn for the power spectrum askn = k0en/m (2.28)where m is the value chosen for k per log int. Therefore, we chose to perform logintegration over the power spectrum when evaluating the Clνν ′. We chose a k perlog int of 2000 to achieve good accuracy.CAMB also has the option of taking into account nonlinearities, which in arealistic case should be accounted for in the correlation function. However, weused the linear power spectrum settings on CAMB as a first approximation to avoidcomplications and to shorten computation time.Because the transfer functions are identical at late times, to get the matter power86at arbitrary z and z′ we can use the matter power spectrum generated by CAMB atz = z′= 0, divide out the growth functions at z=0, and multiply by the new growthfunctions, to obtainPm(k,z,z′) = Pm(k,0,0)D(0)2 D(z)D(z′). (2.29)2.4 Calculation of the Clνν ′Given that the observables in the experiment are frequencies from 400 to 800MHz, we chose an evenly spaced frequency grid in 1 MHz intervals. Using thefrequency redshift relation (1.33), the redshift values z[i] along this frequency gridare, going to lower redshift, in terms of ν = 400+ i MHz,z[i] =1420.405400+ i−1. (2.30)Then for our chosen theory, corresponding to w = −1, we evaluated the Clνν ′for this grid of z values. The only significant contributions to the covariance matrixwere those that differed by a given separation of frequencies up to slightly morethan a BAO range on each side of a central value. Using the relation∆z =1cH(z1)rs (2.31)the range of frequencies that a BAO centered at a given redshift z would cover isgiven byNz = ∆ν(z)≈ dνdz∆z =− 1420rsc(1+ z)2H(z)≈ 50(1+ z)2√Ωm(1+ z)3+ΩDE(1+ z)3(1+w).(2.32)When calculating the Clνν ′, we use the full sky approximation for approximatelyl ≤ 70 and the flat-sky approximation for l > 70. We chose a maximum l of ap-proximately 2000 to evaluate the correlation function accurately.872.5 Window FunctionSince the CHIME data will only probe structure on scales in the range of theBAO radius, we must use window functions on our data. Window functions elim-inate ringing effects from a higher and lower bound to the number of l that weinclude in our correlation function. We will use a Gaussian beam pattern, whichwe explain below as given in ([87]), ([88]). The shape of the window function isdetermined by the beam function, Bi(xˆ · nˆ), where nˆ and xˆ are radial vectors spec-ifying the position on the sky at a redshift slice. In terms of the beam pattern, wehaveδTb(xˆ) =∫dnˆB(xˆ · nˆ)ΣlmalmYlm(nˆ). (2.33)We can write our beam pattern in terms of relative angle θ = cos−1(xˆ · nˆ) on thesky as ([87]), ([88])B(xˆ · nˆ) = 12piσ2e−θ22σ2 , (2.34)so thatB(xˆ · nˆ)dnˆ = θdθdφ2piσ2e−θ22σ2 . (2.35)Here, σ << 1, meaning that it falls off very fast with θ , allowing for the smallangle approximation over the angular integral. Since the beam pattern is a functionof xˆ · nˆ, we can expand it asB(xˆ · nˆ) = 14piΣl(2l+1)blPl(xˆ · nˆ) = ΣlmBlYlm(xˆ)Y ∗lm(nˆ). (2.36)Substituting in this expansion and using the orthogonality relation∫dnˆYlm(nˆ)Y ∗l′m′(nˆ) = δl,l′δm,m′, (2.37)the integral over dnˆ is eliminated, givingδTb(xˆ) = ΣlmblalmYlm(xˆ). (2.38)88The coefficients bl are given by inverting the expansion using the orthogonality ofLegendre polynomials∫ 1−1Pl(cosθ)Pl′(cosθ)d(cosθ)=∫ ∞−∞Pn(1− θ22)Pm(1− θ22)θdθ =22l+1δl,l′(2.39)using the small angle approximation cosθ ≈ 1− θ 22 , obtainingbl =∫ ∞−∞θdθσ2e−θ22σ2 Pl(1− θ22)=∫ ∞0dxPl(1−σ2x)e−x.(2.40)We can approximate this as∫ ∞0dxPl(1−σ2x)e−x = Σ∞n=0(−σ2)nPnl (1)= Σln=0(−σ2)n12nn!(l+n)!(ln)!≈ Σln=0(−σ2)n12nn!(l+12)2n≈ e l(l+1)σ22 .(2.41)We therefore haveδTb(xˆ) = e−l(l+1)σ2almYlm(xˆ) (2.42)andCl,beam = e−l(l+1)σ2Cl. (2.43)We window our data approximating the effect of (2.43) as ([34])Cl,beam = e−l2σ2Cl. (2.44)The full width at half maximum of a Gaussian beam is related to the σ by ([34])σ =FWHM√8ln2. (2.45)89Figure 2.1: Theoretical Correlation Function, 767 MHz (z=.851)We choose a Gaussian beam window function with a full width at half maximumof 0.2 degrees. As the frequency bins are small, we neglect the use of the windowfunction along the line of sight.Sample views of the theoretical correlation function plotted using these meth-ods are shown in Figs. 2.1 , 2.2 , and 2.3 .2.6 Generation of Sample RealizationsNext we enter this covariance matrix into the computer and generate a set ofcomplex Gaussian almv that has the required covariance matrix, subject to the con-straint thata∗lmν = al−mν . (2.46)Sample realizations, plotted using Healpix ([38]), are shown in Figs. (2.4), (2.5),and (2.6).From our simulated almν ,we computed the cosmic variance limited Clνν ′, givenbyCdatalνν ′ =12l+1Σlm=−la∗lmνalmν ′ (2.47)90Figure 2.2: Theoretical Correlation Function, 767 and 780 MHz (z=.851 and.821)Figure 2.3: Theoretical Radial Correlation Function, 767 MHz (z=.851)91Figure 2.4: Simulated 21 cm Brightness Temperature Anisotropies, 767 MHz(z=0.851)Figure 2.5: Simulated 21 cm Brightness Temperature Anisotropies, 780 MHz(z=0.821)92Figure 2.6: Simulated 21 cm Brightness Temperature Anisotropies, 790 MHz(z=0.797)and calculated the correlation function. The angular, radial, and one side viewresult are plotted in Figs. (2.7), (2.8), and (2.9).2.7 Likelihood AnalysisFor the following discussion introducing likelihood analysis, see ([34]). Theformula for the likelihood function for a theory comes from Bayes’ law, which saysP[B∩A] = P[B | A]P[A] = P[A | B]P[B]. (2.48)In other words, the probability of both B and A equals the probability of B givenA times the probability of A, or equivalently the probability of A given B times theprobability of B. The second equality implies that P[{pi},C | {di}], the probabilityof the theory, defined by a set of parameters pi, with a covariance matrix C, whereCIJ = 〈(dI− d¯I)∗(dJ− d¯J)〉, (2.49)93Figure 2.7: Simulated Angular Correlation Function, 767 MHz (z=.851)Figure 2.8: Simulated Radial Correlation Function, 767 MHz (z=.851)94Figure 2.9: Simulated Correlation Function, 767 and 780 MHz (z=.851 and.821)given a set of data {di}, is given byP[{pi},C | {di}] = P[{di} | {pi},C]P[{pi},C]P[{di}] . (2.50)The denominator of the right-hand side will be the same for all theories, being equalto the integral of P[{pi},C | {di}] over all possible values of {pi} and C, and sodoes not depend on the parameters in the theory. If we assume no prior informationabout the correct theory, we can assume that the second term in the numerator ofthe right-hand side, P[{pi},C], is the same for all theories. Therefore, factoringout these two terms, we haveP[pi,C | {di}] ∝ P[{di} | {pi},C]. (2.51)In other words, we have the statement that the probability of a theory given a set ofdata is proportional to the likelihood of the data given the theory. We will use thisstatement to determine the constraints on our parameters.95Given a set of complex valued data drawn from a multivariate Gaussian distri-bution, the likelihood function, defined asL≡ P[{di} | pi,C], (2.52)is given byL =1(2pi)N/2√| C | exp(−12(d− d¯)∗I C−1IJ (d− d¯)J)(2.53)where N is the number of data points. The dI are the values of the data in the array,and d¯I is the mean value of the data according to the theory. C is the covariancematrix of the data predicted by the theory, given by (2.49).In our case, we wish to maximize the likelihood, or minimize the negative ofthe exponential,lnL =−12(ln(| C |)+χ2), (2.54)whereχ2 = (d− d¯)∗I C−1IJ (d− d¯)J. (2.55)Since the determinant generally changes very slowly over the range of parameterspace compared to χ2, we will neglect the determinant in our likelihood function.The points on the noisy correlation function are quadratic in Gaussian randomvariables, and so follow a χ2 distribution. Our correlation function is calculatedfrom the expression (2.10) which makes use of the almν as in (2.9). Since the almνare complex, there are 2(2l+ 1) numbers in total for each l,ν ,ν ′; however, sincethe Clνν ′ are real, this implies almν = a∗l−mν , so the χ2 distribution has 2l + 1 in-dependent degrees of freedom for each l and ν , needed to specify the correlationfunction. By the central limit theorem, a χ2 distribution with N degrees of free-dom approaches a Gaussian distribution for large N ([5]) with the same mean andvariance. Therefore, in our project, we will use a Gaussian form of the likelihoodfunction for the correlation function points C. This is given byL(w) ∝ exp(−12(Cmodel−Cdata)I(C−1)IJ(Cmodel−Cdata)J), (2.56)96normalized to have integral value one. In practice, the covariance matrix we willuse will approximate the covariance of the theoretical points on the correlationfunction, and Cmodel will be a model for the correlation function in the vicinity ofthe peak. We will explain this model below.2.8 Calculating the LikelihoodFollowing a method similar to ([77]) and ([76]), we model the theory as aparabola centered on the sound horizon rs. Our model is:Cmodel(ν ,ν ′,θ ,w) = A−B(r(ν ,ν ′,θ ,w)− rs)2 (2.57)where r(ν ,ν ′,θ ,w) is determined from (1.34) and (1.14), with ν determining theredshift as in (1.33). The sound horizon rs, largely independent of w, is calculatedusing the formula as given in (1.294) using the fitting form for zd (1.289), dividedby a factor of 154.66/150.82 to account for the differences between the CAMBnumerical values and the fitting form above, as has been shown to account for thedifferences in past projects ([55]). We tested that the linear correlation function hada peak at this sound horizon in all directions to a very good approximation whenconverting redshifts and angles to comoving distances. Due to a copying errordiscovered afterwards we had used Ωγh2 = 2.57× 10−5 rather than the correctvalue of 2.47× 10−5 in the code when calculating the sound horizon, but this didnot significantly affect the results. We choose only frequencies in a small range ofthe peak (covering only those frequencies and angles that correspond to 145 to 155MpC if w =−1). After this, we calculate a grid of χ2[A,B,w] and exponentiate toform L[A,B,w] as in (2.56). We marginalize over, or integrate the likelihood ove,rall possible values of, A and B, and finally normalize by dividing by∫L(w)dw. Inthe next section, we explain the steps of the calculation of the constraints from thismodel.2.9 The Covariance MatrixBelow, we explain a formula for the covariance matrix of the points on thecorrelation function based on the values of the theoretical correlation function,97using the input theory of w = −1. In an actual experiment we would only haveaccess to the noisy Clνν ′, in which case it might be useful to use other methods,which have not been explored here. See ([28]), ([32]) for discussion material.2.9.1 Model for the Covariance MatrixThe points of the correlation function, like the Clνν ′, are drawn from a χ2 dis-tribution which approximates a Gaussian distribution for a large number of degreesof freedom. We neglect the determinant. We assume when applying Baysess the-orem that the covariance matrix does not change much with respect to variationsabout the correct theory, and therefore we use the covariance matrix calculated forthe correct theory. We start with the calculation of the covariance matrix of theClνν ′ :〈∆Clνν ′∆Clν ′′ν ′′′〉= 〈(Clνν ′−C¯lνν ′)(Clν ′′ν ′′′−C¯lν ′′ν ′′′)〉= 〈Clνν ′Clν ′′ν ′′′〉−C¯lνν ′C¯lν ′′ν ′′′.(2.58)With the assumption that the Clνν ′ are Gaussian and not χ2 variables, we can applyWick’s theorem, which applies to expectation values of n-point correlation func-tions between Gaussian variables, in the case of n=2. The first term of (2.58) isthen〈Clνν ′Clν ′′ν ′′′〉= 1(2l+1)(2l′+1)Σl−lΣl′−l′〈a∗lmνalmν ′a∗l′m′ν ′′al′m′ν ′′′〉=1(2l+1)(2l′+1)Σl−lΣl′−l′(〈a∗lmνalmν ′〉〈a∗l′m′ν ′′almν ′′′〉+ 〈a∗lmνa∗l′m′ν ′′〉〈almν ′al′m′ν ′′′〉+ 〈a∗lmνal′m′ν ′′′〉〈a∗l′m′ν ′′almν ′〉).(2.59)Using the first term of this expression to cancel the second term in (2.58), thensubstitutinga∗lmν = (−1)lal−mν (2.60)in the second term, and the definition〈a∗lmνal′m′ν ′〉= δmm′Clνν ′, (2.61)98our expression (2.58) becomes〈∆Clνν ′∆Cl′ν ′′ν ′′′〉= δll′ 12l+1(C¯lνν ′′C¯lν ′ν ′′′+C¯lνν ′′′C¯lν ′ν ′′). (2.62)Therefore, evaluating the covariance of the points of the correlation function andsubstituting (2.62), we obtain (see [32]):〈∆C(ν ,ν ′,θ)∆C(ν ′′,ν ′′′,θ ′)〉= 〈∆∑l2l+14piClνν ′Pl(cos(θ))∆∑l2l+14piClν ′′ν ′′′Pl(cos(θ ′))〉=∑l(2l+1)2(4pi)2〈∆Clνν ′∆Clν ′′ν ′′′〉Pl(cos(θ))Pl(cos(θ ′))=∑l2l+1(4pi)2Clνν ′′Clν ′ν ′′′Pl(cos(θ))Pl(cos(θ ′)).(2.63)We calculate this from the theory for the case of w =−1. in a realistic scenario wewill not have access to this, so we would either estimate it using the noisy Clνν ′ orcalculate the deviation from the mean of many realizations of the Clνν ′ from thedata as input. This is frequently done, under the assumption is that the distributionof Clνν ′ about the data is similar to the distribution about the correct theory ([73]).2.9.2 Inverting the Covariance MatrixThe covariance matrix was very ill-conditioned, meaning that the inverse couldnot be calculated accurately to the existing precision due to pile-up error. To do getaround this, we calculated a pseudo inverse using the singular value decomposition([73]). The singular value decomposition is given by decomposing the covariancematrix asC =UWV T (2.64)where U and V are orthogonal matrices and W was the eigenvalue decompositionof the matrix with the very small eigenvalues removed, and inverting to getC−1 =VW−1UT . (2.65)99Figure 2.10: Likelihood of Model Given Theoretical Correlation FunctionWe used a range of usable eigenvalues that ranged about 18 orders of magnitude.This was close to what one should expect since it matched the floating point pre-cision of the python settings we were using to calculate the existing covariancematrix. However, we found that to obtain sensible results for the likelihood, takingthe pseudo inverse itself was not enough; it was necessary to use some increasedprecision, using the python package mpmath, to calculate the pseudo inverse.2.10 Testing the Method for BiasTo determine whether this method is biased, we put the theoretical correlationfunction, calculated for w =−1, in the place reserved for the data in the likelihoodfunction. In Fig. (2.10) we show the likelihood obtained from one BAO region.The maximum of the likelihood function occurs very near w = −1. We concludethat there is a very small bias in the constraints on w from this method.1002.11 Finding Best Estimate Constraints on wWe also had to test whether the likelihood function accurately predicts w towithin the uncertainty posed by cosmic variance. In Fig. (2.11) we show 10 ex-amples of the constraints on w from our likelihood function, each obtained from aseparate realization of data generated from the theoretical Clνν ′. We see that eachof these realizations contains a reasonable likelihood of w =−1. We found an av-erage sample width of a likelihood to be σ¯w,i ≈ .03. We take the best fit parametervalues of w obtained from each of the ten likelihood curves, wmax,i, and calculatethe distribution around the mean,σwmax =√Σi(w¯max−wmax,i)2N−1 . (2.66)We find σwmax ≈ 0.02. To investigate this possibility of there being a systematicdisagreement, we repeatedly generated 10 sets of random numbers from the dis-tribution of the sample likelihood plot approximated as a Gaussian of width 0.03.We then calculated the standard deviation of the 10 values each time, which wecalled σ10,i. We made a histogram of these values to determine the spread of theσ10,i, shown in Fig. (2.12) . We find that our obtained value σwmax = 0.02 was stillwithin the spread of σ10,i to a reasonable confidence level.2.12 Conclusion and Direction for Future WorkIn this thesis, we have shown a method for constraining the equation of state ofdark energy using the full three-dimensional shape of the correlation in the vicinityof the BAO peak and tested the method for one full span of the BAO centered ata sample redshift. This method does not make use of the full shape of the correla-tion function; however, given experimental noise, the full shape of the correlationfunction is difficult to determine precisely, whereas the BAO scale, on which themethod exclusively relies, is considered a robust probe of cosmological parameters([36]), ([22]), ([47]), ([78]), see also ([37]). This method also does not rely on afiducial model for the correlation function, and so is guaranteed to give constraintsin a cosmology-independent way. The computation time is also short enough that101Figure 2.11: Likelihoods of Model Given Respective Physical Realizationsof the Correlation FunctionFigure 2.12: Likelihood of Standard Deviation in w for 10 Sample Realiza-tions102one can marginalize over multiple values of the cosmological parameters in a morecomplete analysis using a program which uses the Levenburg-Marquardt method,such as COSMO MC ([59]).In a complete analysis, to be handled in the future, we will need to use a greaterregion of redshift space to evaluate our likelihood function. Given more time andcomputational resources, we would combine the constraints from BAO’s centeredon widely separated redshift regions. If computationally feasible, we could alsoinclude correlations between BAO’s in neighboring redshift regions; thus the ex-ponential of the likelihood function would be− 12∆C(ν ,ν ′,θ)C−1∆C(ν ′′,ν ′′′,θ ′) (2.67)for all possible values of the arguments in parentheses, to make full use of all theinformation obtainable from CHIME’s 3D map using the BAO as a standard ruler.We would also need to marginalize over the values of the parameters which appearin the argument of the expression for r and rs in (1.294) and (1.34), includingΩm and ΩDE , allow for nonzero curvature, and to allow for the possibility that wdepends on redshift. The most straightforward way to generalize to a non-constantw assuming nothing about the nature of dark energy is given in ([51]). Here, wewould divide the region from z = 0 to z into small redshift intervals(0,z1+∆z1),(z2−∆z2 ,z2+∆z2), ...,(z j−∆z j ,z j +∆z j), (2.68)and assume that w = wi is constant in the region (zi−∆zi ,zi +∆zi), then treat thewi as free parameters to be constrained. The equation of state of dark energy in aregion (z j−∆z j ,z j +∆z j) is given byρDE(z) =ΩDEe−3∫ a0daa (1+w(a))= ρDE(z = 0)(1+ z1+ z j−∆z j/2)3(1+w j)Π j−1i=1(1+ zi+∆zi/21+ zi−∆zi/2)3(1+wi).(2.69)We would also need to take into account redshift space distortions, a nonlinearcorrelation function, and experimental noise.103Bibliography[1] Four-velocity. [Online; accessed05-December-2017]. → pages 41[2] Associated Legendre Polynomials. Legendre polynomials. [Online;accessed 05-December-2017]. → pages 119[3] Baryon Acoustic Oscillations. acoustic oscillations. [Online; accessed05-December-2017]. → pages 11, 64[4] CAMB Web Interface. camb form.cfm. [Online; accessed07-December-2017]. → pages 35[5] Chi-Square: Testing for Goodness of Fit.∼drip/133/ch4.pdf. [Online; accessed05-December-2017]. → pages 96[6] Modern Cosmology: Errata.∼dodelson/errata.html.[Online; accessed 05-December-2017]. → pages 2[7] Jones Calculus. calculus. [Online;accessed 05-December-2017]. → pages 43[8] 4: Cosmological Perturbation Theory. [Online;accessed 05-December-2017]. → pages 65[9] 3: Thermal History., .[Online; accessed 05-December-2017]. → pages 2, 4, 17, 18, 19, 21, 22, 29104[10] Cosmology Part III Mathematical Tripos., .[Online; accessed 05-December-2017]. → pages 17, 19[11] Brightness Temperature. temperature. [Online; accessed05-December-2017]. → pages 74[12] 6: Initial Conditions from Inflation. [Online;accessed 05-December-2017]. → pages 68, 70[13] Plane Wave Expansion. wave expansion.[Online; accessed 05-December-2017]. → pages 121[14] Adiabatic Invariant. invariant, .[Online; accessed 05-December-2017]. → pages 62[15] Canonical Quantization. quantization, . [Online; accessed05-December-2017]. → pages 70, 71[16] Spherical Harmonics. harmonics, .[Online; accessed 05-December-2017]. → pages 119[17] L. Anderson et al. MNRAS, 427:3435–3467, 2012. → pages 13[18] J. M. Bardeen. Physical Review D, 22:1882, 1980. → pages 42[19] R. Barkana and A. Loeb. Rep. Prog. Phys., 70:627, 2007. → pages 76[20] D. Baumann. TASI Lectures on Inflation., 2012. → pages 70[21] E. Bertschinger. Cosmological Perturbation Theory, Nonlinear Models andNumerical Simulations of Cosmic Structure Formation. In Cosmology 2000,page 1.1, 2000. → pages 35[22] C. Blake and K. Glazebrook. Astrophys. J., 594:665, 2003. → pages 2, 12,101[23] M. L. Boas. Mathematical Methods in the Physical Sciences. John Wileyand Sons, 3 edition, 2006. → pages 119105[24] J. Bond and G. Efstathiou. Astrophys. J. Letters, 285:L45, 1984. → pages42, 46[25] P. Bull, P. Ferreira, P. Patel, and M. Santos. Astrophys. J, 803:21, 2015. →pages 12, 75[26] S. Burles and D. Tytler. Astrophys. J., 507:732–744, 1998. → pages 5[27] S. Burles, K. Nollett, and M. Turner. Big-Bang Nucleosynthesis: LinkingInner Space and Outer Space.,1999. → pages 5[28] A. Cabre´, P. Fosalba, E. Gaztan˜aga, and M. Manera. MNRAS, 381:1347–1368, 2007. → pages 98[29] R. Caldwell and M. Kamionkowski. Ann. Rev. Nucl. Part. Sci., 59:397–429,2009. → pages 2[30] S. Carroll. Spacetime and Geometry: An Introduction to General Relativity.Addison Wesley, 2004. → pages 35, 39, 47[31] T.-C. Chang, U.-L. Pen, J. B. Peterson, and P. McDonald. Phys. Rev. Lett.,100:091303, 2008. → pages 76[32] M. Crocce, A. Cabre´, and E. Gaztan˜aga. MNRAS, 414:329–349, 2011. →pages 98, 99[33] K. K. Datta, T. Choudhury, and S. Bharadwaj. MNRAS, 378:119, 2007. →pages 81, 83, 85[34] S. Dodelson. Modern Cosmology. Academic Press, 2003. → pages 2, 3, 4,5, 6, 7, 8, 9, 11, 12, 15, 16, 17, 18, 21, 22, 24, 37, 38, 42, 43, 46, 49, 51, 52,53, 54, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 70, 72, 73, 74, 76, 79, 82,89, 93[35] D. Eisenstein and W. Hu. Astrophys. J., 496:605, 1998. → pages 64, 66[36] D. Eisenstein, W. Hu, and M. Tegmark. Astrophys. J., 504:L57–L61, 1998.→ pages 2, 12, 101[37] D. Eisenstein et al. Astrophys. J., 633:560–574, 2005. → pages 2, 11, 12,13, 64, 101[38] K. Go´rski, B. Wandelt, E. Hivon, F. Hansen, and A. Banday. The HEALPixPrimer., 2010. Version 2.15a.[Online; accessed 05-December-2017]. → pages 90106[39] A. Guth. Phys. Rev. D., 23:347, 1981. → pages[40] A. Hall, C. Bovin, and A. Challinor. Phys. Rev. D, 87:064026, 2013. →pages 76[41] M. Halpern. CHIME: a Progress Report. Canadian HydrogenIntensity-Mapping Experiment., 2017.[Online; accessed 05-December-2017]. → pages 1, 13[42] T. Harmark. Notes on Lie Derivatives and Killing Vector Fields .∼harmark/Killingvectors.pdf, 2008. [Online; accessed05-December-2017]. → pages 114[43] W. Hu. Astro 448: The Cosmic Microwave Background.∼whu/Presentations/ast448n.pdf, . [Online;accessed 05-December-2017]. → pages 48[44] W. Hu. Ast 448 Set 2: Perturbation Theory.∼whu/Courses/Ast448 15/ast448 2.pdf, .[Online; accessed 05-December-2017]. → pages 16, 41[45] W. Hu. Wandering in the Background: A Cosmic Microwave BackgroundExplorer. Ph.D. Thesis, University of California at Berkeley., 1995. → pages 2, 11, 18, 35, 37,38, 39, 42, 43, 46, 47, 48, 49, 50, 51, 54, 55, 57, 58, 64, 65, 66[46] W. Hu. Damping.∼whu/araa/node12.html,2001. [Online; accessed 05-December-2017]. → pages 48, 62[47] W. Hu and Z. Haiman. Physical Review D., 68:3004, 2003. → pages 2, 12,101[48] W. Hu and N. Sugiyama. Astrophys. J., 444:489–506, 1995. → pages 54, 55[49] W. Hu and N. Sugiyama. Phys. Rev. D, 51:2599, 1995. → pages 53, 54, 55,57[50] W. Hu and N. Sugiyama. Astrophys. J., 471:542–570, 1996. → pages 2, 11,28, 46, 49, 50, 51, 53, 58, 59, 60, 62, 63, 64, 65, 66[51] D. Huterer, E. Linder, and J. Weller. Constraints on Dark Energy and itsModels.∼huterer/Papers/hutwel.pdf.[Online; accessed 05-December-2017]. Published in: Yellow Book on DarkEnergy, SNOWMASS 2001, eds. E. V. Linder. → pages 103107[52] Y. I. Izotov and T. Thuan. Astrophys. J., 500:188, 1998. → pages 5[53] H. Kodama and M. Sasaki. Prg. Theor. Phys. Suppl., 88:1, 1984. → pages37, 38, 39, 41, 42, 54, 55[54] H. Kodama and M. Sasaki. Int. J. Mod. Phys., A1:265, 1986. → pages 53[55] E. Komatsu. Email correspondence with G. Hinshaw. → pages 97[56] E. Komatsu et al. Astrophys. J. Suppl., 192:18, 2011. → pages 3, 17[57] A. Kosowsky. Ann. Phys., 246:49, 1996. → pages 43[58] D. Langlois. C.R. Physique, 4:953–959, 2003. → pages 51[59] A. Lewis and S. Bridle. Phys. Rev. D, 66:103511, 2002. → pages 103[60] E. Lifshitz. J. Phys. USSR, 10:116, 1946. → pages 35[61] A. Loeb and J. Wyithe. Physical Review Letters, 100:161301, 2008. →pages 1[62] C. Ma and E. Bertschinger. Astrophys. J., 455:7–25, 1995. → pages 35, 43,46[63] K. Masui et al. Astrophys. J., 763:L20, 2013. → pages 76[64] J. Mather et al. Astrophys. J., 420:439–444, 1994. → pages 5[65] H. Mo, F. van den Bosch, and S. White. Galaxy Formation and Evolution.Cambridge University Press, 2010. → pages 9, 59[66] V. Mukhanov. Physical Foundations of Cosmology. Cambridge UniversityPress, 2005. → pages 43[67] L. Newburgh et al. Calibrating chime: a new radio interferometer to probedark energy. volume 9145, page 9145, 2014. doi:10.1117/12.2056962. →pages 1, 76[68] P. Peebles. Astrophys. J., 153:1, 1968. → pages 29[69] P. Peebles. Large Scale Structure of the Universe. Princeton UniversityPress, 1980. → pages 65[70] P. Peebles and J. Yu. Astrophys. J., 162:815, 1970. → pages 59, 61, 63[71] A. Penzias and R. Wilson. Astrophys. J, 142:419, 1965. → pages 5108[72] W. Press and P. Schechter. Astrophys. J, 187:425–438, 1974. → pages 10[73] W. H. Press, B. P. Flannery, S. A. Teukolsky, and W. Vetterling. NumericalRecipes: The Art of Scientific Computing. Cambridge University Press,1986. → pages 99[74] A. G. Riess et al. Astron. J., 116:1009–1038, 1998. → pages 8[75] A. Sa´nchez, C. Baugh, and R. Angulo. MNRAS, 390:1470–1490, 2008. →pages 86[76] E. Sa´nchez, D. Alonso, F. Sa´nchez, J. Garcı´a-Bellido, and I. Sevilla.MNRAS, 434, Issue 3:2008, 2013. → pages 14, 97[77] E. Sa´nchez et al. MNRAS, 411:277–288, 2011. → pages 13, 14, 97[78] H. Seo and D. Eisenstein. Astrophys. J., 598:720, 2003. → pages 2, 12, 101[79] J. Shaw, K. Sigurdson, U. Pen, A. Stebbins, and M. Sitwell. Astrophys. J.,781:57, 2014. → pages 76, 78[80] J. Shaw, K. Sigurdson, M. Sitwell, A. Stebbins, and U.-L. Pen. Astrophys. J.,781:57, 2014. → pages 74, 75, 76, 78, 79, 83, 85[81] M. Shull and A. Venkatesan. Optical Depth of the Cosmic MicrowaveBackground and Reionization of the Intergalactic Medium., 2007. → pages 86[82] R. M. Wald. General Relativity. University of Chicago Press, 1984. →pages 3, 4, 15, 16, 17, 69, 111, 114[83] Q. Wang, Z. Zhu, and W. Unruh. Phys. Rev. D, 95:103504, 2017. → pages 3[84] D. Weinberg et al. Physics Reports, 530, Issue 10:87–255, 2013. → pages 7[85] S. Weinberg. Corrections to the First Edition of Cosmology.∼weinberg/swcorrections.pdf. [Online; accessed05-December-2017]. → pages 2[86] S. Weinberg. Cosmology. Oxford University Press, 2008. → pages 2, 4, 16,17, 19, 21, 24, 28, 29, 38, 39, 40, 42, 43, 45, 46, 69, 71, 83, 118[87] M. White. Phys. Rev. D, 46:4198–4205, 1992. → pages 88[88] M. White and M. Srednicki. Astrophys. J., 443:6, 1995. → pages 88109[89] X. Xu et al. MNRAS, 427:2146, 2012. → pages 13[90] S. Zaroubi. The Epoch of Reionization. in: The First Galaxies - Theoretical Predictions and ObservationalClues”, Springer, 2012, eds. T. Wiklind, B. Mobasher and V. Bromm. →pages 12, 76[91] Y. B. Zel’dovich, V. Kurt, and R. Sunyaev. JETP, 28:146, 1969. → pages 29110Appendix AGauge FreedomWe will now explain, as given in ([82]), the motivation for the statement thatthe gauge freedom of the metric and stress energy tensor in general relativity isdetermined by a Lie derivative. In general relativity, spacetime is described as ann dimensional manifold. This is a set of points and subsets Oα such that eachpoint is in at least one subset, for each Oα , there is a one to one and onto mapψα from this subset to an open subset of Rn, and for any overlapping subsets Oαand Oβ on the manifold, for the corresponding maps ψα and ψβ , ψα [Oα ∩Oβ ] andψβ [Oα ∩Oβ ] are open sets, and the map between them Ψβ ◦Ψ−1α is C∞ (havinginfinitely many derivatives) . A manifold M consists of a tangent space Vp andcotangent space V ∗p at each point p, on which we can define, respectively, vectorsand dual vectors. Tangent space vectors v are defined as maps from the space ofC∞ functions f (p) : M → R, where R denotes real numbers, to R, which satisfylinearity: for a and b in R,v(a f +bg) = av( f )+bv(g), (A.1)and the Leibnitz rule: for two C∞ functions f and g : M→ R,v( f g) = f (p)v(g)+g(p)v( f ). (A.2)111For each coordinate mapping ψ, we can construct a basis for the tangent space,given byXµ( f ) =∂∂xµ( f ◦ψ−1)|ψ(p), (A.3)for µ = 1 to n. In terms of this, we can specify all vectors v as a linear combinationof the basis vectors with components vµ asv( f ) =n∑µ=1vµXµ( f ). (A.4)It will also be useful to define a tangent vector T of a curve C : R→M, by, for allfunctions f in the space of C∞ functions from M to R,T ( f ) =ddt( f ◦C). (A.5)Given a basis vµ for the tangent space, we can also define a basis of the cotan-gent space v1∗, ...,vn∗ of dual vectors, which are linear maps f from vectors to realnumbers, by the requirement thatvµ∗(vν) = δµν . (A.6)With respect to the tangent and cotangent space, we can define (m,n) type tensors,or multilinear maps V ∗p × ...×V ∗p ×Vp× ...×Vp→ R, which act on m copies of thecotangent space and n copies of the tangent space. These are denoted by T µ1...µmν1....νn ,using abstract index notation where the superscripts and subscripts are used tospecify that the tensor is a type (m,n) tensor, rather than to refer to the value of itscomponents in a coordinate system.Two spacetimes, considered as manifolds M and N, are equivalent if and only ifthey are related by a diffeomorphism φ , or injective (one to one) C∞ map betweenpoints on the manifold with an inverse that is C∞. Assuming that this holds, if werequire that a vector on N is equivalent at each point translated by the diffeomor-phism to the vector on the original point at M, then for all C∞ functions f1, f2 suchthat f2(φ(p)) = f1(p) for all points p on M,(φ ∗v)( f2(φ(p)) = v( f1(p)), (A.7)112which requires that f2 = f1 ◦ φ−1. Therefore, the corresponding map φ ∗ : Vp →Vφ(p) between tangent vectors v on the two manifolds is defined, for f on N, by(φ ∗v)( f ) = v( f ◦ φ). (A.8)We can at the same time define a mapping of dual vectors by φ∗ : V ∗φ(p)→ V ∗p byrequiring that for all vectors v in Vp,(φ∗µ)ava = µa(φ ∗v)a. (A.9)Similarly, we can extend φ ∗ : Vp→Vφ(p) to general tensors, as(φ ∗T ) (µ1)b1 ...(µk)bk(t1)a1 ...(tl)al =T b1...bka1...ak (φ∗µ1)b1 ...(φ∗µk)bk([φ−1]∗t1)a1 ...([φ−1]∗tl)al ,(A.10)and to extend φ∗ to general tensors we would use the fact that φ∗ = (φ−1)∗. Thesemaps again are necessary to ensure that a tensor on N is equivalent at each pointtranslated by the diffeomorphism to the tensor on the original point at M.Two physical theories will then be equivalent if the spacetimes are related bya diffeomorphism φ and the tensors found by the equations of the theory are re-lated by the transformations T ′ = φ ∗T given above. Conversely, as mentionedearlier, the only way that two spacetimes are equivalent is if they are related by adiffeomorphism φ . Also, if physical tensors are not related, for one of these dif-feomorphisms, by the corresponding transformation T ′ = φ ∗T, for some point thetensor will differ and therefore so will the theories. Therefore, the diffeomorphismrelations above account for all gauge, or unphysical, degrees of freedom.A one-parameter family of diffeomorphisms φt is a C∞ map from R×M→M,such that for each parameter value t in R, φt : M→M is a diffeomorphism and forall t,s in R, φt ◦ φs = φt+s. For an arbitrary point p ε M, φt(p) : R→M defines acurve, called the orbit of φt(p), which passes through p at t = 0 and has tangentvector v |p at t = 0. Since φt(p) is a well defined function, each point is on one andonly one such curve with this choice of parametrization. The requirement that, foreach point p, the orbit of φt(p) have tangent vector v |p at t = 0 then allows φt tospecify a unique vector field va, which is called the generator of φt(p). It can alsobe proven that for each possible smooth vector field va defined on M, there is only113one possible set of curves, considered as maps R→M, known as integral curves,given their starting points on the manifold at parameter value t = 0, such that foreach p, v |p is the tangent to a curve at p, and one and only one curve in the setpasses through each point p ε M. Given the definition of a tangent vector (A.5),and the expressions (A.3) and (A.4), this curve can be specified in a coordinatesystem by solvingvµ =dxµdt(A.11)for xµ(t) at each point p subject to specifying which point the curve runs throughat t = 0. Assuming the integral curves for this vector field can take all possiblecurve parametrization values, we can require that for a choice of parametrizationof the integral curve of v |p containing p such that t = 0 at p, φt(p) is the pointat parameter value t along this curve. We can then verify that φt(p) is a one-parameter family of diffeomorphisms. Given an arbitrary one-parameter familyof diffeomorphisms φt(p) generated by va, since for all s, φt ◦φs(p) = φt+s(p), wehave that φt ◦φs(p) covers the same points on the curve as φt(p) except with startingpoint at φs(p) instead of p, and with a parameter value differing by s at all points,and with tangent vector v |φs(p) at t = 0. This can then be repeated for all otherpoints belonging to different orbits under the diffeomorphism, so we have a set ofcurves such that one and only one curve in the set passes through each point, up toa choice of starting point. However, from above, there is only one set of curves, ormaps R→M, given arbitrary starting points at t = 0,with tangent vectors v |p, suchthat one and only one curve in the set passes through each point, so φt(p) is thencompletely specified. Therefore, given a generator, the diffeomorphism generatedis unique.As explained in ([42]), ([82]), we will use a Lie derivative to quantify the rateof change under a diffeomorphism φt with respect to a change in parameter valuet. Since φt determines a generator va as described above, and since the generatorva of φt specifies the diffeomorphism uniquely, we have a well defined expressionif we write the Lie derivative as Lv, so that we = limt→0φ ∗−t(−Ta1...anb1....bmt. (A.12)114Equivalently, the Lie derivative of a tensor gives the rate of change of the tensoralong the integral curve of the generator va of the diffeomorphism, or in the direc-tion of the tangent vector of the integral curve. Similarly, we can define the Liederivative of a function by the rate of change along the integral curve of va withrespect to t, or with respect to the diffeomorphism φt , asLv( f ) = limt→0f (φt(p))− f (p)t. (A.13)Therefore, from (A.5) Lv( f ) = T ( f ) for the curve φt(p), so we haveLv( f ) = v( f ). (A.14)We now wish to find an expression for the Lie derivative of a tensor for a givenone-parameter family of diffeomorphisms φt in terms of the generating vector fieldva for φt . As a starting point for this, we can use a coordinate system on M where talong the integral curves of va is chosen as one of the coordinates x1. Then given theproperties of diffeomorphisms on tensors described above, φ−t enacts a coordinatetransformation x1→ x1+ t, so from (A.12) we haveLvTµ1...µiν1...ν j =∂T µ1...µiν1...ν j∂x1. (A.15)As a partial derivative operator, the Lie derivative satisfies the Leibnitz rule forproducts of components of tensors. Also, from the above expression (A.11) for thetangent vector, and using (A.3) and (A.4), we have for this coordinate systemva =(∂∂x1)a, (A.16)where we have used abstract index notation. Next, we wish to define the Lie deriva-tive of a tensor, given by (A.15) in a specific coordinate basis, in a coordinate in-dependent manner in terms of va. An arbitrary derivative operator ∇b on tensors ingeneral relativity is defined as an operator that satisfies the conditions of linearity,the Leibnitz rule, commutativity with contraction, for all functions f and vectorsta, t( f ) = ta∇a f , and the torsion free condition (∇a∇b−∇b∇a) f = 0. Using these115properties we can find an expression for the commutator vector, for an arbitraryderivative operator ∇b,[v,w]( f )≡ v{w( f )}−w{v( f )}= (va∇awb−wa∇avb)∇b f ,(A.17)so that we have[v,w]a = vb∇bwa−wb∇bva. (A.18)We are given that in our choice of coordinate system,Lvwa =∂wa∂x1. (A.19)With respect to the derivative operator ∂∂xµ we have[v,w]µ = Σv(vv∂wµ∂xν−wν ∂vµ∂xν). (A.20)From (A.16), this evaluates to the right-hand side of (A.19), so we then can definethe Lie derivative on a vector field asLvwa = [v,w]a. (A.21)We can then extend this to dual vectors and tensors. By (A.13), we haveLv(uawa) = v(uawa). (A.22)However, by the properties of an arbitrary derivative operator ∇a, as stated above,we havev(uawa) = vb∇b(uawa)= vbwa∇bua+ vbua∇bwa.(A.23)116Using the Leibnitz rule (A.2) and (A.21) , we havev(uawa) = waLvua+ua[v,w]a. (A.24)Using (A.23) and the expression (A.18) for the commutator we findLvua = vb∇bua+ub∇avb. (A.25)In a similar way, we can find by applying induction vc∇cT −Σ∇cvai +Σ∇b j vc. (A.26)Given a mean unperturbed metric g¯ab, it is possible to determine a one-parameterfamily of spacetimes gab(λ ), if these metrics can be expanded asgab = g¯ab+ γabλ , (A.27)where γab = dgabdλ |λ=0, so that λ quantifies the strength of the perturbation aboutg¯ab. As explained above, two metric perturbations about the same manifold M areequivalent if and only if they are related by a diffeomorphism. We also wish toexpandφ ∗gµν = φ ∗gµν +ddλ(φ ∗gµν) |λ=0 λ , (A.28)so that, assuming the first term is equivalent after the diffeomorphism φ ∗ has beenapplied, ddλ (φ∗gµν) |λ=0 represents an equivalent physical perturbation. This thenrequires that we are able to express φ ∗ in terms of λ . The most general diffeomor-phism from M to M depending on λ for all possible values of λ is provided by aone-parameter family of diffeomorphisms φλ , which was defined above. We thenhaveγab =dgabdλ∣∣∣∣λ=0, γ ′ab =d(φ ∗λgab)dλ∣∣∣∣λ=0. (A.29)We can then verify, for the metric perturbations hµν = γµνλ and h′µν = γ ′µν , using117(A.12) and (A.26), thath′µν = hµν −Lξgµν= hµν −ξ λ∇λgµν +gλν∇µξ λ +gµλ∇νξ λ ,(A.30)where ξ a is the vector field generating φλ . Evaluating the Lie derivative to firstorder in hµν and ξ λ , using ordinary differentiation for the derivative operator, gives([86])∆hµν(x) =−g¯λµ(x)∂νξ λ (x)− g¯λν(x)∂µξ λ (x)−∂λ g¯µν(x)ξ λ (x). (A.31)As explained above, all other tensors, including the energy-momentum tensorsmust also undergo the same transformation if the physics is to be equivalent. Tofirst order in the perturbations to the stress energy tensor and ξ λ , we have ([86])∆δTµν(x) =−LξTab=−T¯λµ(x)∂νξ λ (x)− T¯λν(x)∂µξ λ (x)−∂λ T¯µν(x)ξ λ (x).(A.32)118Appendix BSpherical HarmonicsHere we give an overview of spherical harmonics, relying on content from([23]), ([16]), ([2]). Any function on a sphere can be written in terms of the solu-tions to Laplace’s equation, which is given by∇2 f (r,θ ,φ) = 0. (B.1)This is solved using the technique of separation of variables by writing the solutionin the most general form asf (r,θ ,φ) = R(r)Θ(θ)Φ(φ). (B.2)The equation is the sum of a φ dependent function and a jointly θ and r dependentequation. We can then separate out the φ dependence as1Φd2Φdφ 2=−m2. (B.3)The solution to the φ∗ dependent equation isΦ(φ) = eimφ , (B.4)119where to ensure spherical symmetry, m is an integer. For this m, the equation canthen be further decomposed into a radial equation1Rddr(r2dRdr)= k, (B.5)and the angular equation, for this solution k and m, which is given by1sin(θ)ddθ(sin(θ)dΘdθ)− m2sin2θΘ+ kΘ= 0. (B.6)The solutions to the second equation, if we assume that k = l(l + 1) for a giveninteger value of l, are the associated Legendre polynomialsΘ(θ) = Plm(cos(θ), (B.7)wherePlm(x) =(−1)m2ll!(1− x2)m/2 dl+mdxl+m(x2−1)l, (B.8)which are zero for |m|> l and obeyPlm(x) = (−1)m (l−m)!(l+m)!Pl−m(x). (B.9)The complete angular solution Θ(θ)Φ(φ) is therefore given by the spherical har-monics,Ylm(θ ,φ) =√2l+14pi(l−m)!(l+m)!Plm(cos(θ))eimφ (B.10)which obeyY ∗lm = (−1)mYl−m. (B.11)The Ylm satisfy an orthogonality condition∫ piθ=0∫ 2piφ=0dΩYlm(θ ,φ)Yl′m′(θ ,φ) = δl,l′δm,m′, (B.12)120and a completeness relation, and so serve as a basis for functions on a sphere. Fromthe spherical harmonics, we can form Legendre polynomials:Pl(xˆ · yˆ) = 4pi2l+1Σlm=−lY∗lm(xˆ)Ylm(yˆ), (B.13)where xˆ and yˆ are radial unit vectors pointing in the directions xˆ and yˆ. Finally,it will also be useful to know that a plane wave can be expanded in Legendrepolynomials using the Raleigh expansion ([13]),ei~k·~r = Σl(2l+1)il jl(kr)Pl(kˆ · rˆ). (B.14)121


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items