UBC Theses and Dissertations

UBC Theses Logo

UBC Theses and Dissertations

Warm dark matter collapse : real space analysis methods Loranger, Jonathan 2013

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

Item Metadata


24-ubc_2014_spring_loranger_jonathan.pdf [ 337.23kB ]
JSON: 24-1.0085565.json
JSON-LD: 24-1.0085565-ld.json
RDF/XML (Pretty): 24-1.0085565-rdf.xml
RDF/JSON: 24-1.0085565-rdf.json
Turtle: 24-1.0085565-turtle.txt
N-Triples: 24-1.0085565-rdf-ntriples.txt
Original Record: 24-1.0085565-source.json
Full Text

Full Text

Warm Dark Matter Collapsereal space analysis methodsbyJonathan LorangerB.Sc., The University of Guelph, 2010A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF SCIENCEinThe Faculty of Graduate and Postdoctoral Studies(Physics)THE UNIVERSITY OF BRITISH COLUMBIA(Vancouver)October 2013c? Jonathan Loranger 2013AbstractWe study the evolution of a warm dark matter and perfect fluid systemto determine its behaviour in the linear regime. Comparative analysis isperformed between cold dark matter, hot dark matter and warm dark matterapproximating each case. Numerical issues causes differences between thewarm dark matter approximations and the respective case. Numerical issuesthat we have been unable to solve prevent the calculation of sufficient k-spacemodes to study interesting scales. Analytic methods to obtain the real spaceperturbations and distribution functions are derived.iiPrefaceFor my Master?s program, I was tasked with of solving the warm dark mattercollapse problem, assigned to me by my supervisor, Dr Kris Sigurdson. Inmy research, my supervisor, and Dr Annika Peter, mainly helped by pointingme towards possible fixes for programming bugs I encountered, and suggest-ing literature to review, but overall, I was in charge of writting the code,finding the numerical problems, and fixing them. Once we obtained data, Iwas responsible for performing the analysis and determining what methodswould be used, under the direction of my supervisor and Dr Peter.None of the details of this thesis have been published.iiiTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . vii1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11.1 Cosmological Definitions . . . . . . . . . . . . . . . . . . . . 21.2 Units . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 Unperturbed System . . . . . . . . . . . . . . . . . . . . . . . . 52.1 Metric Defnitions . . . . . . . . . . . . . . . . . . . . . . . . 52.2 Metric Calculations . . . . . . . . . . . . . . . . . . . . . . . 73 Derivation from the Metric . . . . . . . . . . . . . . . . . . . 104 Perfect Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164.1 Radiation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17ivTable of Contents4.2 Cold Dark Matter . . . . . . . . . . . . . . . . . . . . . . . . 185 Derivation from the Boltzmann Equation . . . . . . . . . . 195.1 Massless Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . 245.2 Massive Neutrinos . . . . . . . . . . . . . . . . . . . . . . . . 276 Initial Conditions and Normalization . . . . . . . . . . . . . 297 Numerical Details . . . . . . . . . . . . . . . . . . . . . . . . . 357.1 Energy Density Normalization . . . . . . . . . . . . . . . . . 367.2 Numerical Routines . . . . . . . . . . . . . . . . . . . . . . . 378 Simulation Results . . . . . . . . . . . . . . . . . . . . . . . . . 398.1 Perfect Fluid Comparison . . . . . . . . . . . . . . . . . . . . 398.2 Cold Dark Matter Comparison . . . . . . . . . . . . . . . . . 408.3 Massless Neutrino Comparison . . . . . . . . . . . . . . . . . 428.4 Warm Dark Matter Data . . . . . . . . . . . . . . . . . . . . 449 Data Analysis Routines . . . . . . . . . . . . . . . . . . . . . . 489.1 Transfer Function . . . . . . . . . . . . . . . . . . . . . . . . 489.2 Real Space Perturbations . . . . . . . . . . . . . . . . . . . . 499.3 Distribution Function in Real Space . . . . . . . . . . . . . . 5210 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57vList of Figures8.1 Perfect Fluid Comparison . . . . . . . . . . . . . . . . . . . . 408.2 High Mass Warm Dark Matter Comparison . . . . . . . . . . 418.3 Low Mass Warm Dark Matter, Absolute Difference . . . . . . 438.4 Warm Dark Matter: Evolution of ? . . . . . . . . . . . . . . . 458.5 Warm Dark Matter: Evolution of ? . . . . . . . . . . . . . . . 458.6 Warm Dark Matter: Evolution of ? . . . . . . . . . . . . . . . 46viAcknowledgementsI would like to thank my supervisor, Dr Kris Sigurdson, for the advice anddirection given during this thesis. I would also like to thank Dr AnnikaPeter for her help in problem solving and debugging the code. Thanks toDr Douglas Scott for ensuring my thesis was properly formatted. I thankmy parents for supporting me throughout my education. And lastly, I thankmy fiance, Sam Snobelen, for pushing me to continue when things seemedbleak.viiChapter 1IntroductionUnderstanding how planets, stars, galaxies formed has been a point of inter-est throughout history. With advances in general relativity, statistical me-chanics and electromagnetism from the last century, we are now in a positionto attempt to explain how these structures are formed. The current modelused for structure formation is the ?CDM model [3]. In this model, thereexists the standard model particles, along with cold dark matter (CDM) anda cosmological constant, ?, acting as dark energy. The presence of CDM anda cosmological constant in the model are due to observations that suggestthat visible matter does not account for the all gravitational force observablein our Universe. In effect, visible matter accounts for only roughly 5% of allmatter, with dark matter being approximately 27% and dark energy repre-senting the remaining 68% [1]. In cosmology,?dark? signifies that there isno interaction with electro-magnetic force. [9] This translates to no photonproduction/absorption during interactions with dark matter/energy. Thislack of interaction causes dark matter/energy to be elusive in nature, andtheir existence is inferred from gravitational interactions, which are muchweaker then electro-magnetic interactions. The ?cold? descriptor for colddark matter relates to the matter?s momentum. A cold particle exhibits a11.1. Cosmological Definitionssmall p/m ratio, where p is the momentum andm, the mass [9]. We will lookat the detailed implications of this in further chapters, but in short, it causesmany simplifications to the system. There has been much work done to un-derstand how cold dark matter collapses to form structures [5][15][14][2][10].This thesis will look at the effect of making the dark matter ?warm?. In thiscontext, warm signifies particles for which p/m begins relatively large andbecomes insignificant at late times [9]. Before we continue the discussion ofwarm dark matter (WDM) collapse, we will define some basic concepts incosmology, and look at the collapse of CDM in more detail.1.1 Cosmological DefinitionsIn the field of cosmology, time and distance are not uniquely defined quanti-ties. There exist numerous different ways to measure both of these quantities[3]. For distances, the curvature of space-time itself and the expansion ofthe Universe cause some issues of definition. When considering expansion,distances measured today were not the same as in the past, and this changein distance is significant. In order to deal with this change, the scale factor,a, is defined. If a distance x is measured today, it is at a distance of ax atanother time. This definition sets a0 = 1, where the subscript 0 indicatestoday. As our Universe is expanding, distances were smaller in the past, andhence past values of a are smaller then 1 [3]. The exact behaviour of a isbased on the composition of the Universe[3].In order to find a definition of time, we must decide if we desire propertime or coordinate time[3]. Proper time is measured based on a comoving21.1. Cosmological Definitionsobserver, while coordinate time is measured by an outside observer. Let usconsider the distance light can travel in an expanding universe since t = 0from the view of a comoving observer. The velocity of light can be takento be 1, as we are free to set c = 1. This choice of value for c gives usthat d/t = 1. If we consider only an infinitesimal time step dt, we obtaindt = a(t)dx. As we wish to calculate the distance light has travelled we areinterested in?(t) =? x(t)x(t=0)dx? =? t0dt?a(t?) . (1.1)As dt and a are always positive, ? is an increasing function. This propertymakes ? an ideal candidate for a time variable. As a time coordinate, it isrefered to as ?conformal time?[3]. This qunatity contains some importantinformation. As it represents the distance light has travelled, it also gives themaximum distance two objects can be apart to have been in causal contactin the past. We call this distance the ?horizon?[3]. As will be seen later, wewill be working in Fourier transformed k-space for much of our calculations.In k-space, modes for which k? < 1 have not yet crossed the horizon andthese are called ?super-horizon? modes, while those with k? > 1 are ?sub-horizon? modes. These two ranges of k will have distinct behaviour, as wewill see once we have derived the relevant equations for the evolution ofk-modes.31.2. Units1.2 UnitsIn this paper, we will be using only Mpc as units. In order to achieve this,we set the following constants to 1 : c, kB, ~. In this system, quantites areusually expressed in terms of eV. We use G to change from eV to Mpc. Thevalue of G in the appropriate units is G = 163.9964Mpc?2eV?4.4Chapter 2Unperturbed SystemWe begin our study with an unperturbed universe. It is isotropic and homo-geneous. All directions and all parts of space are the same. We will relax thehomogeneous condition as we progress, as we need inhomogeneities in thematter distribution to cause matter to gravitationally collapse and producethe structure we see in the Universe. We will need to determine the metricto use, derive the general relativistic equations, and define the quantities ofinterest. Once this framework is built, we will add perturbations to it toeventually obtain our desired warm dark matter system.2.1 Metric DefnitionsThe metric contains information about the curvature of space-time. This isone of two components necessary to determine the general relativity equa-tions, the other component being the stress-energy tensor. For an isotropicand homogeneous expanding Universe, the metric to use is[9]ds2 = ?a2d?2 + a2dxidxi. (2.1)We use conformal time as our time variable. ds2 is the proper distance52.1. Metric Defnitionsin space-time in this metric. Roman indices are used to sum over the spatialindices (1,2,3), while greek indices will also sum over the temporal indice(0). The metric may also be written asds2 = g??dx?dx? . (2.2)We must work our way from this metric to solving for the Einstein tensor,G?? :G?? ? R?? ?12g??R = 8?GT?? , (2.3)where, R?? is the Ricci tensor, R is the Ricci scalar and T?? is the stress-energy tensor[3]. We will go through the process of defining and deriving allthese quantities for this system. We begin withR?? = ????,? ? ????,? + ???????? ? ????????, (2.4)where ???? is the Christoffel symbol, while the comma represents a derivative,for example R,? is the derivative of R with respect to x?[3]. The Christoffelsymbol is defined as[3]???? =g??2 [g??,? + g??,? ? g??,? ]. (2.5)Continuing our definitions,R ? g??R?? (2.6)62.2. Metric CalculationsandT?? = Pg?? + (?+ P )U?U? , (2.7)where U? = dx?/??ds2 is the four-velocity[9].2.2 Metric CalculationsWe now have all the pieces to solve the unperturbed system. We obtain thefollowing Christoffel symbols :?000 =( a?a); (2.8a)?0ij = ?i0j = ?ij0 = ?ija?a . (2.8b)All other ?s are zero. Here the dot represents a derivative with respectto ? .We now use the Christoffel symbols to obtain the Ricci tensor :R00 = ?3dd?a?a ; (2.9a)Rij = ?ij{dd?( a?a)+ 2( a?a)2}. (2.9b)The other Ricci tensor coefficients are zero.The Ricci scalar can now be calculated. We obtainR = 6a2( a?a). (2.10)72.2. Metric CalculationsThe stress-energy tensor, T?? , is the last piece we need. It is defined asT 00 = ??, (2.11a)T 0i = T i0 = 0, (2.11b)T ij = ?ijP, (2.11c)where ? and P are the energy density and the pressure, respectively.It will be noted that we used a different configuration of indices for T??then we used in the Einstein equation (2.3). This is done to avoid scalefactors appearing into the stress-energy tensor. We will now need to changethe Einstein tensor to this form in order to solve the Einstein equation. Tochange indices, we use T?? = g??T??.After the change in indices, we obtain the final Einstein equations :( a?a)2= 8?G3 a2?; (2.12)dd?( a?a)= ?4?G3 a2(?+ 3P ). (2.13)We will redefine the energy density and pressure as follows ?? = 8piGa43 ?and P ? = 8piGa43 P , changing the above equations toa? =???, (2.14)dd? (a?a) = ??? + 3P ?2a2 . (2.15)We may use these two equations to solve for a?, and correspondingly a.Once we have a, we know how the unperturbed universe evolves, solving it82.2. Metric Calculationscompletely. These equations are valid for a flat space-time. For a space-time that is positively or negatively curved, an additional term appears,modifying the first equation by the addition of a ?? term. ? represents theoverall curvature of space-time[3]. Hence, for a positively curved universe,? > 0, and ? < 0 for a negatively curved universe. We will be mainlyconcerned with a universe with ? = 0 but the spherical collapse case willuse ? > 0.The next step in our derivation is to add perturbations to the energy-stress tensor and determine their evolution equations.9Chapter 3Derivation from the MetricThere exists two common gauges used in cosmology for the evolution of per-turbations, the synchronous gauge and the conformal Newtonian gauge.[9]In this thesis, the focus will be solely on the conformal gauge. This choice ismade due to the physical interpretation of the conformal Newtonian gauge.The gravitational perturbation ? is equivalent to the Newtonian gravita-tional potential, and as such, previous intuition can be used when analyzingthe gravitational perturbation. It should be noted that the majority of thederivations in Chapters 3,4 and 5 fill in the details of the derivations from[9].For the conformal Newtonian gauge, the perturbed metric equation canbe written as [9]ds2 = a2(?){?(1 + 2?)d?2 + (1? 2?)dxidxi}. (3.1)We may perform a similar analysis for this metric as was performedfor the unperturbed state. The equation for the Einstein tensor can bedecomposed as follows [3]:G?? + ?G?? = 8?G(T?? + ?T?? ). (3.2)10Chapter 3. Derivation from the MetricAs we have already solved the unperturbed parts, these will cancel, and weobtain?G?? = 8?G?T?? . (3.3)This will prevent us from repeating the derivations of the unperturbed sys-tem.From Ref. [7], the perturbed Christoffel symbols are:??000 = ??; (3.4a)??00i = ??0i0 = ?,i; (3.4b)??0ij = ??ij{2H(?+ ?) + ??}; (3.4c)??i00 = ?,i; (3.4d)??i0j = ????ij ; (3.4e)??ijk = (?,k?ij + ?,j?ij) + ?,i?jk. (3.4f)Having these perturbed Christoffel symbols, we can now calculate theperturbed Ricci tensor coefficients:?R00 = a?2{?3????2? ? 3 a?a(?? + ??)? 6dd?( a?a)?}; (3.5a)?R0i = ?2a?2(??+ a?a?)i= ?Ri0; (3.5b)?Rij = a?2[???+?2?? a?a(?? + 2??)?{2 dd?( a?a)+ 4( a?a)2}?]?ij+ a?2(?? ?),ij .(3.5c)11Chapter 3. Derivation from the MetricIt should be noted that due to the nature of index raising/lowering, theperturbed portion of the Ricci tensor is ?g??R?? + g???R?? .This now allows us to find the perturbed curvature scalar,R = a?2[?6??+ 2?2(2?? ?)? 6 a?a(?? + 3??)?12{dd?( a?a)+( a?a)2}?].(3.6)Using the metric for the conformal Newtonian gauge, (3.1), we can de-termine ?T?? using (2.7):?T 00 = ???; (3.7a)?T 0i = (?+ P )vi = ??T i0; (3.7b)?T ij = ?P?ij + ?ij ; ?ii = 0. (3.7c)We can now combine the perturbed equations in order to obtain ourEinstein equations:a?2{?2?2?+ 6 a?a ??+ 6( a?a)2?}= a?23???, (3.8a)?2a?2(??+ a?a?),i= ?a?23(?? + P ?)vi, (3.8b)a?2([2??+?2(? ? ?) + a?a(2??+ 4??)+{4 dd?( a?a)+ 2( a?a)2}?]?ji + (?? ?),ij)= a?2(12?P?)+ ?ji .(3.8c)The third equation in this set is rather cumbersome. We may decompose12Chapter 3. Derivation from the Metricit into the trace and traceless components. Before doing this, we will convertto momentum space via a Fourier transform. It will be simpler to decompose(3.8c) in momentum space, and we will also be working in this space once weadd interactions, as in linear theory it will decouple our equations [3], thusit will be good exercise to convert our equations now. The primary effect ofthis is the conversion of , i into ?iki. This gives us ?2 = ?k2. Also, sincethe ki vectors are Euclidean constructs, you can interchange their indicesfreely[3]. Focusing only on (3.8c), we obtaina?2([2??+ k2(?? ?) + a?a(2??+ 4??) +{4 dd?( a?a)+2( a?a)2}?]?ji ? kikj(?? ?))= a?2(12?P??ji)+ ?ji .(3.9)We will decompose this equation into two components, the trace andtraceless parts. For the trace, ?ii ? 3 and kiki ? k2. The trace is thena?2([2??+ k2(?? ?) + a?a(2??+ 4??) +{4 dd?( a?a)+ 2( a?a)2}?]3?k2(?? ?))= a?2(32?P?)+ ?ii.(3.10)Since ?ji is defined as traceless, ?ii = 0. Dividing both sides by 3a?2, weobtain2??+ k2(?? ?) + a?a(2??+ 4??) +{4 dd?( a?a)+ 2( a?a)2}??a2k23 (?? ?) = ?P?.(3.11)13Chapter 3. Derivation from the MetricNow, we deal with the traceless component of (3.8c). We can accomplishthis by applying the traceless operator k?ik?j ? 13?ij [3]:(k?ik?j ?13?ij)kikj(?? ?) =(k?ik?j ?13?ij)?ji . (3.12)Let us define(?+ P )? = ikj?T 0j , (3.13a)(?+ P )? =(k?ik?j ?13?ij)?ij . (3.13b)Simplifying the left hand side, we obtain our final equation23k2(?? ?) = 3(?? + P ?)?. (3.14)For ease, we will now restate the four basic equations ask2?+ 3 a?a(??+ a?a?)= 32???, (3.15a)k2(??+ a?a?)= 32(?? + P ?)?,(3.15b)??+ a?a(??+ 2??) +{2 dd?( a?a)+( a?a)2}? ? k23 (?? ?) =12?P?, (3.15c)k2(?? ?) = 92(?? + P ?)?,(3.15d)where, ? = kivi[3].This completes the main derivations for the evolution of the background14Chapter 3. Derivation from the Metricand gravitational perturbation quantities. In the subsequent chapter, wewill derive the equations for the perturbations for various types of particles,but beforehand, we will solve a simple system, the perfect fluid.15Chapter 4Perfect FluidBefore continuing on to consider the effects of particle interactions on theevolution of matter, we will look at the simplest type of matter, a perfectfluid. A perfect fluid is defined solely by its energy density, ? and pressure,P[9]. This simplifications allow us to solve the entire system using energy-momentum conservation.In general relativity, this principle is embodied by the following [3]T??;? = ??T?? + ????T?? + ????T ?? = 0, (4.1)where ; represents the covariant derivative, which means it is invariant undergauge transforms.[9] Looking at the perturbed part of these equations, weobtain the following 2 equations:?? = ?(1 + w)(? ? 3??)? 3 a?a(?P?? ? w)?; (4.2a)?? = ? a?a(1? 3w)? ?w?1 + w? +?P/??1 + w k2? ? k2? + k2?; (4.2b)Here w = P/?.We will not cover these derivations in detail here, but aside from a fewminor details, they are fairly straightforward calculations. The two main164.1. Radiationdetails that must be used to obtain them is that ?? = ???? ??????2 , and that?a2 =( a?a)2 using our units. Using these two substitutions, it should be anexercise in algebra to obtain the previous equations.There are two types of perfect fluids that we will be interested in, w = 1/3and w = 0 fluids. w = 1/3 fluids are radiation fluids while w = 0 fluidsrepresent cold dark matter (CDM)[9]. A radiation fluid will be included inour system to keep matter-radiation equality near to that of our Universe.The CDM is being considered as a test case for high mass WDM.For both of these perfect fluids, w? = 0 and ? = 0[9].4.1 RadiationSubstituting w = 1/3 into (4.2), we obtain?? = ?(4/3)(? ? 3??), (4.3a)?? = k2(? + ?). (4.3b)If only this perfect radiation fluid is assumed to exist, and (3.15a) or (3.15b)are used to evolve ?, this set of equations has an analytic solution of theform [4]? = 9?0(k?)3[?3 sin( k??3)? k? cos( k??3)], (4.4a)? = 6?0(k?)3[2?3((k?)2 ? 3) sin( k??3)? k?((k?)2 ? 6) cos( k??3)], (4.4b)? = ?3?3k?02(k?)2[2?3k? sin( k??3)? ((k?)2 ? 6) cos( k??3)], (4.4c)174.2. Cold Dark Matterwhere ?0 is the initial value of ?.This set of analytic solutions will be used to verify the numerical accuracyof the code, in the special case where only the perfect fluid is present. Asthe system of interest will have WDM mixed in with the perfect radiationfluid, the actual solution will differ from this analytic solution.4.2 Cold Dark MatterSubstituting w = 0 into (4.2), we obtain?? = ?? + 3??, (4.5a)?? = ? a?a? + k2?. (4.5b)Again, this set of equations will be used to verify the code for the high masscase, where we expect the WDM to behave like CDM.18Chapter 5Derivation from theBoltzmann EquationHaving now determined how the background and the gravitational potentialbehaves, we must now derive the evolution of the matter. Although weare primarily concerned with WDM, a bottom-up approach was used tobuild the code, where we started with the simplest systems and added tothem until we reached WDM. For this reason, we have covered the evolutionof cold dark matter, and will now consider massless neutrinos and massiveneutrinos, where massive neutrinos will be the basis for our WDM [9]. Beforewe begin to think about the evolution for a specific case of matter though,we need to examine the distribution function and the Boltzmann equation.It will be noted that this chapter is strongly based on Ref. [9].The Boltzmann equation utilizes phase space, which is a six dimensionalspace, having 3 positions, xi and 3 momenta, Pi = mUi. The positionsmay be represented by xi, but the momenta will need to be examined moreclosely.19Chapter 5. Derivation from the Boltzmann EquationIn the conformal gauge, we havePi = (1? ?)api, (5.1a)P i = ?(1 + ?)a3pi. (5.1b)Let us define the phase space distribution function asf(xi, Pj , ?)dx1dx2dx3dP1dP2dP3 = dN, (5.2)where dN is the number of particles in the volume dx1dx2dx3dP1dP2dP3.We will decompose the distribution function,f , into an unperturbed andfirst order component, as followsf(xi, Pj , ?) = f0(1 + ?(xi, Pj , ?)). (5.3)In our previous definition of P , the metric perturbation appeared. Thiswould complicate many calculations, but we may replace Pj by qj?apj [9]. Asqj is not the conjugate momentum, due to the lack of the (1??) term, we willneed to be careful when working with the energy-momentum tensor, and thevolume element. The volume element becomes dx1dx2dx3dq1dq2dq3(1? 3?)to first order. We can also define the comoving energy as ? =?q2 + a2m2.As this is the energy measured by a comoving observer, it is related to P0,20Chapter 5. Derivation from the Boltzmann Equationthe energy component of the momentum four-vector, throughP0 = ?(1 + ?)?, (5.4a)P 0 = (1? ?)a?2?. (5.4b)Another desired modification is the separation of the magnitude of q fromits direction. This can be easily accomplished with the following definitionqj = qnj where nini = ?ij . Our phase space distribution now depends on theposition, q, nj and time.The zeroth order distribution function, f0 will be the Fermi-Dirac func-tion for fermions(+) and Bose-Einstein for bosons (?). As the particlesdecouple while highly relativistic, q ? am in the distribution function, andas such, we may approximate ? as q in the distribution function only:f0(q) =gsh3P1eq/kBT0 ? 1 (5.5)where gs is the number of spin degrees of freedom, hP is the Planck constant,kB is the Boltzmann constant, and is set to 1, and T0 is the temperaturetoday. For the remainder of the discussion, we will redefine ?, q and m asfollows:q ? qT0; m? mT0; ?? ?T0. (5.6a)This is the notation used by CLASS [8].We may decompose the distribution function into a zeroth order and21Chapter 5. Derivation from the Boltzmann Equationhigher order components:f(xi, q, nj , ?) = f0(q)(1 + ?(xi, q, nj , ?) + ...), (5.7)where ? represents the first order perturbation. We will ignore the higherorder perturbations.We have previously written the energy-momentum tensor in terms of themetric (2.7). We can also express it in terms of the distribution functionT?? =?dP1dP2dP3(?g)?1/2P?P?P 0 f(xi, Pj , ?), (5.8)where g is the determinant of the g?? . For the conformal gauge, g = ?a8(1?6? + 2?). It should also be noted that, as discussed earlier, dP1dP2dP3 =T 30 dq1dq2dq3(1? 3?).Combining all of these in equation (5.8), we can obtain the energy-momentum tensor:T 00 = ?(T0a)4 ?q2dqd??f0(q)(1 + ?); (5.9a)T 0i =(T0a)?q3dqd?nif0(q)?; (5.9b)T ij =(T0a)?q3dqd?ninjq? f0(q)(1 + ?). (5.9c)It becomes easy to split this energy-momentum tensor into perturbedand unperturbed components as the only first order quantity present is ?.Because there exists simplifications to the integrand for various types ofmatter, we will wait to determine equivalence between this equation and22Chapter 5. Derivation from the Boltzmann Equationthe quantities from equation (3.7).As ? depends on ? , its evolution must be determined. We will use theBoltzmann equation to calculate this. The Boltzmann equation isDfd? =?f?? +dxid??f?xi +dqd??f?q +dnid??f?ni=(?f??)C, (5.10)Where right hand side is a collision term. As we will be dealing withdark matter, this term will be 0 for all cases we are interested in.There is one simplification that can be done immediately. Since only? depends on ni, ?f?ni is a first order quantity. Also, only the gravitationalpotentials, ? and ?, may cause changes to ni, and so dnid? is also a first orderquantity, making the combined term a second order perturbaion[3]. As weare only interested in first order quantities, we may ignore this term.We will now determine expression for dxid? . Multiplying by??ds2 weobtain [3]dxi??ds2??ds2d? =P iP 0 , (5.11a)= ?(1 + ? + ?)q? . (5.11b)There is now only one term left to calculate, dqd? . We will use the geodesicequation,[3]P 0 dp?d? + ????P?P ? = 0. (5.12)Setting ? = 0, we obtaindqd? = q??+ i?niki?. (5.13)235.1. Massless NeutrinosUtilizing all of these terms, we can obtain our general expression for theBoltzmann equation:???? + iq? (~k ? n?)? + df0dqqf0[?? ? i ?q (~k ? n?)?]= 0. (5.14)It should be noted that the only dependance on n? is through its dotproduct with ~k. As only the angle between both terms is significant, we mayset the direction of n? and only consider the angle between both vectors.There remains one more detail to iron out before utilizing equation(5.14). We must determine the k-space expression for ?. This will be donedifferently for the massless neutrino and massive neutrino cases that we willbe considering,i.e. the hot dark matter, and warm dark matter candidates,respectively.5.1 Massless NeutrinosFor massless neutrinos, m = 0, causing ? = q. This causes massless neutri-nos to be a type of hot dark matter (HDM) as they remain relativistic atall times. With this simplification, equations (5.9) and (3.7) may be usedtogether to obtain expressions for ? and P and their perturbations. We245.1. Massless Neutrinosobtain,? = 3P =(T0a)4 ?q3dqd?f0(q), (5.15a)?? = 3?P =(T0a)4 ?q3dqd?f0(q)?, (5.15b)?T 0i =(T0a)4 ?q3dqd?nif0(q)?, (5.15c)?ij =(T0a)4 ?q3dqd?(ninj ?13?ij)f0(q)?. (5.15d)We will integrate out the q dependance, and expand the angular depen-dence into a series of Legendre polynomials, Pl(k? ? n?)F (~k, n?, ?) ??q3dqf0(q)?q3dqf0(q)???l=0(?i)l(2l + 1)Fl(~k, ?)Pl(k? ? n?). (5.16)With this new definition, we have removed the dependance on q, andisolated the dependence on n?. Combining (5.15), (5.16) and (5.9), we mayobtain? = 14??d?F = F0, (5.17a)? = 3i16??d?k(k? ? n?)F = 34kF1, (5.17b)? = ? 316??d?[(k? ? n?)2 ? 13]F = 12F2. (5.17c)It is now time to simplify the Boltzmann equation for the massless neu-trinos. We can integrate (5.14) over?q3dqf0(q) and divide by?q3dqf0(q).The only term that is not evident is? d ln f0d ln q q3f0dq?q3f0dq . Using integration by parts,255.1. Massless Neutrinosit can be shown that this term is equal to ?4. With these, we obtain?F?? + ik?F = 4(??? ik??). (5.18)where ? = k? ? n?.Using the orthonormality of the Legendre polynomials, along with therecursion relation ?Pl(?) = 12l+1((l + 1)Pl+1 + lPl?1), we may decompose(5.18) into individual moments to obtain?? = ?43? + 4??, (5.19a)?? = k2(14? ? ?)+ k2?, (5.19b)F?l =k2l + 1[lFl?1 ? (l + 1)Fl+1]. (5.19c)It should be noted that each mode is only coupled to the neighbouringmodes. Also, we have now decoupled our single real space equation into aninfinite set of harmonic equations. We will choose some maximum mode,lmax, as a cutoff. We must be careful with this truncation though, as theincorrect lmax+1 mode will affect lmax and then trickle down the hierarchy. Inorder to minimize this error, we will use the following recurrence relationshipFlmax+1 =2lmax + 1k? Flmax ? Flmax?1, (5.20)which is inspired by the spherical Bessel behaviour exhibited by numericalsolutions.[9]We will choose to truncate at l = 2000, as is suggested in Ref. [9]. Much like265.2. Massive Neutrinosthe CDM case will be a test for a high mass WDM, the massless neutrinoswill be used as a test case for the low mass WDM.5.2 Massive NeutrinosUnlike with the massless neutrinos, the massive neutrinos do not allow forq = ?, as m 6= 0. This allows for no simplifications of equation (5.9). Also,since ? depends on ? and q, we may not integrate over q. As such, we willexpand ? itself as the following Legendre series?(~k, n?, q, ?) =??l=0(?i)l(2l + 1)?l(~v, q, ?)Pl(k? ? n?). (5.21)Applying this expansion to (5.9), we obtain?? = 4?(T0a)4 ?q2dq?f0(q)?0, (5.22a)?P = 4?3(T0a)4 ?q2dq q2? f0(q)?0, (5.22b)(?+ P )? = 4?k(T0a)4 ?q3dqf0(q)?1, (5.22c)(?+ P )? = 8?3(T0a)4 ?q2dq q2? f0(q)?2. (5.22d)Using the properties of the Legendre polynomials, as with the massless275.2. Massive Neutrinoscase, we may decouple the moments of ? in (5.14) as follows:??0 = ?qk? ?1 ? ??d ln f0d ln q ; (5.23a)??1 =qk3? (?0 ? 2?2)??k3q?d ln f0d ln q ; (5.23b)??l =qk(2l + 1)?(l?l?1 ? (l + 1)?l+1), l ? 2. (5.23c)As with the massless case, a truncation is necessary. We will use??(lmax+1) =(2lmax + 1)?qk? ??lmax ???(lmax?1). (5.24)We may trunctate at a much lower value of l than in the massless neutrinocase. This is due to the fact that the massive neutrinos will become non-relativistic, causing high multipole moments to decay rapidly. We choosel = 50, as suggested by [9].28Chapter 6Initial Conditions andNormalizationNow that we have derived all the equations to evolve, we must determine theinitial conditions. The first step in this process is to define a starting pointfor the simulation. At early times, during the radiation domination, ? = aa? .We can see this fact by looking at equation (2.14). ? ? a?4 for a radiationdominated Universe. As such, a? is constant. Hence, a = a?? , giving us thataa? = ? .[9]We set a such that our WDM begins radiation dominated, and useequation (2.14) to calculate a?. From equation (2.14), we also obtain that( a?a)2 = ??2 =?a2 . The value of ? obtained must also satisfy k??1. Ad-ditionally, since our WDM is radiation dominated, we may treat it as amassless neutrino at early times.There exists some freedom in choosing an initial condition. For a detaileddescription on the different options, see Ref. [7]. We choose the adiabaticinitial conditions for our case. This means that we set entropy perturbationsto 0. This is a similar case to the analytic solution in equation (4.4) thatwe studied earlier. By examining those solutions, we find that ?? = ?? = 0.29Chapter 6. Initial Conditions and NormalizationWe can use this, along with our early time details from earlier to transformequation (3.15), (4.3) and (5.19) intok2?+ 3?2? = ?32??, (6.1a)k2? ? = 2??, (6.1b)k2(?? ?) = 6??, (6.1c)??wdm = ??pf = ?43?, (6.1d)??wdm = k2(14?wdm ? ?wdm + ?), (6.1e)??pf = k2(14?pf + ?), (6.1f)??wdm =415?, (6.1g)where wdm and pf are used to denote the WDM and the perfect fluid cases,respectively. Also, substituting ? = ??2 into (6.1), and ignoring the (k?)2factors, we obtain?2? = ?, (6.2a)k2?2 ? = ?, (6.2b)k?26 (?? ?) = ?, (6.2c)??wdm = ??pf = ?43?, (6.2d)??wdm = k2(14?wdm + ?), (6.2e)??pf = k2(14?pf + ?), (6.2f)??wdm =415?wdm. (6.2g)30Chapter 6. Initial Conditions and NormalizationWe have removed the ?wdm term from equation (6.2d) as it is propor-tional to (k?)2 as seen from equation (6.2c). We can see from these equa-tions, that ??wdm and ??pf follow the same evolution. The same can be saidabout the ? equations. With this, we know that ?wdm = ?pf and ?wdm = ?pf.We must relate ? to ?wdm and ?pf, and similarly for ? and ?. This is donethrough? = (1?Rwdm)?pf +Rwdm?wdm, (6.3a)? = (1?Rwdm)?pf +Rwdm?wdm, (6.3b)? = Rwdm?wdm, (6.3c)whereRwdm =?wdm?wdm + ?pf. (6.4)We see from these, that if ?wdm = ?pf, ? = ?pf. The same holds for ?. Wenow have expressions for ? and ?. We wish to find an expression for ? now.Let us consider equation (6.2g). If we substitute our result from equation(6.2b) and integrate over ? , we obtain?wdm =(k?)215 . (6.5)Using this with equations (6.3c) and (6.2c), we obtain? = (1 +Rwdm)?. (6.6)We have now expressed all quantities of interest in terms of ?. As we31Chapter 6. Initial Conditions and Normalizationare working in the linear regime, we are free to choose an initial value of ?.We will choose ? = ?12 , such that ? = 1 initially.Although we have all physical quantities determined, we need to findinitial conditions for ?l in order to evolve our WDM past the radiationdominated domain. To do this, we will observe the evolution equations forthe moments. These reduce to??0 = ?qk? ?1, (6.7a)??1 =qk3? (?0 ? 2?2)??k3q?d ln f0d ln q , (6.7b)??2 =qk5 (2?1). (6.7c)Again, we will ignore the second moment in ??1, as it is proportional to(k?)2. We know that at early time, when am? q, ? = q. This allows us torewrite equation (5.22a) as?? =?q3f0(q)?0dq?q3f0(q)dq. (6.8)When deriving the massless neutrino case, we saw that if there was afactor of d ln f0d ln q , the integral fraction would simplify to ?4. We can see that?0 = ?14?wdmd ln f0d ln q (6.9)satisfies the integral equation. Putting this solution into our evolution equa-32Chapter 6. Initial Conditions and Normalizationtions, we obtain the following [9]:?1 = ??3qk?wdmd ln f0d ln q ; (6.10a)?2 = ?12?wdmd ln f0d ln q . (6.10b)With these initial conditions, we have all the necessary initial conditions forthe WDM+pf system we will be interested in, but we need to determine theinitial conditions for the CDM case.For the pure CDM case, we begin with a matter dominated universe,which gives us a?a = 2??1 = ?1/2. There are only three variables to solve forin the CDM system, ?, ?, ?, since ? = ?. We will begin our derivation of theinitial conditions by rewritting our equations, and substituting in our valuefor a?a . This gives us the following set of equations:k2?+ 12?2? = ?32??; (6.11a)2k2? ? =32??; (6.11b)? = ?; (6.11c)??cdm = ??cdm; (6.11d)??cdm = ?2? ?cdm + k2?; (6.11e)Here we set ?? = ?? = 0, as in the WDM case.By substituting ? = 4?2 , and eliminating terms proportional to (k?)2, we33Chapter 6. Initial Conditions and Normalizationobtain?cdm = ?2?, (6.12a)?cdm = k2??3 . (6.12b)We now have our CDM initial conditions. We are again free to set an initialvalue for ?.When we run a high mass WDM case, we will need initial conditionsfor ?l. We will follow a similar derivation of these as in the radiationdominated WDM. The main difference exists in the integral performed overq to obtain ?. Where we previously multiplied by q3f0dq before integrating,we now multiply by amq2f0dq. The difference is due to the estimate thatwas previously ? = q, but now is ? = am. By performing this new integraland dividing by ? as before, we obtain ?3 instead of ?4. This gives us?0 = ?13?cdmd ln f0d ln q , (6.13a)?1 = ??3qk?cdmd ln f0d ln q , (6.13b)?2 = 0. (6.13c)This concludes the derivation of the initial conditions. We may now begindiscussing details of the numerical routines to be used to solve the equationsderived.34Chapter 7Numerical DetailsThere are two major numerical elements that must be discussed here, thevalues to be used for the energy densities, and the numerical routines utilizedduring the simulation.It should be noted at this point that we have not defined what our perfectfluid or WDM will consist of, namely, we have not determined what f0 isfor either component. For the WDM, we will choosef0(q) =gs(2?)31eq + 1 , (7.1)while the perfect fluid will bef0(q) =3(2?)31eq + 1 +1(2?)31eq ? 1 . (7.2)The perfect fluid distribution function is made to represent a mix of photonsand neutrinos, while the WDM distribution is chosen to be comparable tothe massive neutrino case studied earlier.357.1. Energy Density Normalization7.1 Energy Density NormalizationNormally gs represents the number of degrees of freedom, but since we arenot specifying a specific WDM candidate, we shall leave gs as a free param-eter, giving us the freedom to set ?wdm to a desirable quantity. Looking atour equations, the only time the factor of gs will appear will be in conjecturewith T 40wdm. As this is also the only time T0wdm explicitly enters the equa-tions due to our scaling, this allows us to set the full product of gsT 40wdm todetermine the density. As we are not setting a specific value for T0wdm, wewill need to set mT0wdm instead of m.For the perfect fluid, we shall consider a mix of fermions and bosons.From Ref. [3], we can determine that?pf =?215{1 + 3(78)( 411)4/3}T 4pf. (7.3)We have assumed that the fermions are massless neutrinos and that thebosons are photons. This is done to attempt to simulate the real Universe.With this goal in mind, we shall also set ?crit to be that seen in the realniverse. We will also set ?pf/?wdm at a = 1 to be equal to the ?rad/?matterobserved.Using the radiation to matter ratio and the critical density, we maydetermine ?pf and ?wdm and use these to determine T0pf and gsT 40wdm. Forthe perfect fluid, we use equation (7.3), and for the WDM, we compare ?wdmwith 4?gs(T0wdma )4?q2dq?f0(q)/gs. We have now discussed how to set up thesystem and the equations that will be evolved within it.367.2. Numerical Routines7.2 Numerical RoutinesThere are two primary numerical calculations that are accomplished by thesimulation, solving the differential equations and performing the integrals.For the differential equation solver, a Runge-Kutta routine will be used,with adaptive step size. We will use the routine odeint to drive the solver.For more details on these routines consult Ref. [12].For the integration routine, we will use the Gauss-Laguerre quadratureroutine described in Ref. [8]. This is a routine that works well for exponen-tially decaying functions that must be integrated from 0 to ?. The basicquadrature expression isI =? ?0dqf0(q)g(q) ?n?i=1Wig(qi). (7.4)For the Gauss-Laguerre quadrature, we have this basic expression? ?0dqe?qh(q) ?n?i=1wih(qi), (7.5)where qi are the ith roots of Ln, the nth Laguerre polynomial. The weightswi are determined fromwi =qi(n+ 1)2[Ln+1(qi)]2. (7.6)If we set h(q) = eqf0(q)g(q), we may relate wi to Wi throughWi = wieqif0(qi). (7.7)377.2. Numerical RoutinesThis allows us to calculate the weights once for a set of q values, and thenmultiply them by the appropriate function g, depending on the integral.Another advantage of this routine is that one only needs a few q modes toobtain suitable precision. We choose to use 20 q modes for our case, as thisleads to convergence.Using both of these numerical routines, it is now possible to run thesimulation, and obtain results.38Chapter 8Simulation ResultsBefore looking at out results for the full simulation of the WDM+pf model,we will run the checks that have been discussed in previous chapters. Wewill verify first that the perfect fluid matches the analytic solution. We willthen compare the high mass WDM to the CDM cases and the low massWDM to the massless neutrino solution.8.1 Perfect Fluid ComparisonWe will be comparing the results from our numerical analysis and the ana-lytic solution, (4.4).We can see from Fig. 8.1 that there is no noticeable difference betweenboth analytic and numerical solutions in the pure perfect fluid case at latetimes. When k? ? 1, numerical issues cause the analytic solution to beinaccurate, but our computed solution matches with the expected value forthe analytic solution. Once k? becomes larger, there is a perfect matchbetween both cases. As the analytic solution is only valid for a radiationdominated universe, we will not be able to use the analytic solution for theWDM+pf system once there is a significant matter component.398.2. Cold Dark Matter Comparison-0.015-0.01-0.00500.0050.010.0150.020.0250.030.0350.040.001 0.01 0.1 1 10 100 1000 10000%DifferenceTime (Mpc)? comparaison? comparaison? comparaisonFigure 8.1: Perfect fluid comparison: Percentage difference between theanalytic solution and numerical solution for a perfect fluid. At early times,the solution does not match, but this is due to the numerical issues whenevaluating the analytic solution when k? ? 1.8.2 Cold Dark Matter ComparisonHaving now confirmed that the perfect fluid component is properly cal-culated, we will now focus on making sure that the WDM component isworking properly. We will begin with comparing a very high mass WDMscenario with the CDM case. We choose our mass such that the particles arenon-relativistic at the start of the simulation. In order for both these casesto be comparable, we must evolve the high mass WDM alone, without theperfect fluid component. This ensures that we are in a matter dominateduniverse.We can see from Fig. 8.2 that there are numerical differences betweenboth cases. The difference remains small during the time of integration. The408.2. Cold Dark Matter Comparison00.511.522.50.01 0.1 1 10 100 1000 10000%DifferenceTime (Mpc)? comparaison? comparaison? comparaisonFigure 8.2: High mass warm dark matter comparison: Percentage differencebetween CDM and WDM with m/T0 = 1015. At early times, the differenceis caused by a numerical inaccuracy in the integral to determine ?. At latetimes, the inaccuracy grows. The magnitude of the difference remains smallfor during the time of interest.error is caused by the Laguerre quadrature integration method being inac-curate. One possible cause of this inaccuracy are the roots of the Laguerrepolynomials used to calculate the weights. The values obtained throughusing the Gnu Scientific Library (GSL) Laguerre routine do not agree withthe ones from Mathematica.[16] The values from the GSL routine matchedthe roots found when using the recursion relation,L0(x) = 1,L1(x) = 1? x,L(k + 1)(x) =11 + k ((2k + 1? x)Lk(x)? kLk?1(x)).(8.1)418.3. Massless Neutrino ComparisonHence the GSL roots were used. It is possible that the GSL routine uses thesame recursion relation and that the same numerical inaccuracies appeared,causing both sets of roots to match, and causing the Mathematica roots tobe the proper solution.Another possible cause for the inaccuracy, is the number of momentspresent. In the CDM case, only the first two modes are considered as theothers should be 0. When evolving the WDM case, there is the possibilityof power leaking into higher moments as we do not automatically set themto 0. Even when setting the WDM case to only have 2 moments, thereare similar innacuracies. Without having precise integration methods, it ishence difficult to judge if there are significant errors introduced by powerleakage.There now remains only one more area where the code may not be func-tional. We must now look at the low mass extreme limit.8.3 Massless Neutrino ComparisonWe will now look at the final test case, the low mass WDM. This regimeshould behave like the massless neutrino case. In order to make sure thatour WDM has low enough mass, we set mT0 ? 1. This ensures that ? ? q.It should be noted that Fig. 8.3 shows the actual difference and notpercentile difference like the other plots. This is due to the oscillatory natureof the solution causing division by 0 if we tried to plot a relative difference.From Fig. 8.3, we see a good match between both cases at early times,given that ? varies between ?2 and 3, while ? is always between ?1 and428.3. Massless Neutrino Comparison-0.5-0.4-0.3-0.2- 1e-08 1e-06 0.0001 0.01 1 100 10000DifferenceTime (Mpc)? comparaison? comparaison? comparaisonFigure 8.3: Low mass dark matter, absolute difference: Difference betweenmassless neutrinos and WDM withm/T0 = 10?9. At early times, both casesare seen to agree, but as k? > 1, these cases start diverging.1, and ? oscillates between ?0.0015 and 0.0025. Both cases start divergingwhen k? > 1. The integration method does not always perform the bestaround origin crossing, due to the convergence check relying on the valueitself. As our solution crosses 0, the integrator may be causing numericalinaccuracies, but it is difficult to ascertain due to the large number of equa-tions, causing numerical issues difficult to track. Another possible source ofinaccuracy is the number of moments present. As the low mass WDM needsto be solved for an array of q modes, it was not computationally feasibleto solve it for the recommended 2000 modes or more with the equipmentI had access to. As such, there may be power reflecting back to the lowermoments instead of propagating.438.4. Warm Dark Matter Data8.4 Warm Dark Matter DataAlthough some results from our simulation may be questionable, we willnevertheless perform the WDM+pf simulation. We choose m such thatam ? q initially, but such that m ? q. The behaviour of the solution willdepend on the properties of the system when the modes cross the horizon,which occurs when k? = 1. Hence, small values of k will enter the horizonwhen ? ? am, and resemble CDM, while large k modes will cross when? ? q, behaving like massless neutrinos when crossing the horizon. Sinceam will eventually dominate ?, we expect it will begin behaving like CDMat late times, but with a diminished amplitude, due to the decay presentduring the transition between the HDM and CDM regimes. The differentialequation solver does not seem to converge when the modes cross while stillin the HDM regime, or nearing transition to the CDM regime. As such, wemay only look at purely CDM modes, and modes that transition to CDMshortly after horizon crossing. We will see the behaviour of the low k modesin Figs. 8.4, 8.5 and 8.6, where we look at the perturbations for a numberof different k values.Current versions of the code have difficulty with large k modes. Weare limited to k < 3Mpc?1 as this is the highest mode that will completeits evolution. Also, it should be noted that increasing values of k takesincreasing time to calculate. The inaccuracies mentioned earlier may causeerrors to accumulate, causing convergence of the differential equations to beslower and to also eventually cause the system to stop converging.Once the numerical issues are solved, we may transform the k space data448.4. Warm Dark Matter Data0204060801000.001 0.01 0.1 1 10 100 1000 10000?Time (Mpc)k = 0.0030k = 0.2060k = 1.2623Figure 8.4: Warm dark matter: Evolution of ? for k = 0.0030 Mpc?1,k = 0.2060 Mpc?1, and k = 1.2623 Mpc?1 with m/T0 = 107. We see thatthese modes lead to an increasing amlitude as expected.-1.2-1-0.8-0.6-0.4-0.200.001 0.01 0.1 1 10 100 1000 10000?Time (Mpc)k=0.0030k = 0.2060k = 1.2623CDMFigure 8.5: Warm dark matter: Evolution of ? for k = 0.0030 Mpc?1,k = 0.2060 Mpc?1, and k = 1.2623 Mpc?1 with m/T0 = 107. We seethat the higher k modes are not pure CDM. The increase in ? shortly afterk? = 1 represents the momentum of the particles inhibiting the collapse ofthe system, as HDM would.458.4. Warm Dark Matter Data-0.6-0.5-0.4-0.3-0.2- 0.01 0.1 1 10 100 1000 10000?Time (Mpc)The evolution of ? for different k modes for the WDM+pf fluidk = 0.0030k = 0.2060k = 1.2623massless neutrinoFigure 8.6: Warm dark matter: Evolution of ? for k = 0.0030 Mpc?1, andk = 0.2060 Mpc?1, k = 1.2623 Mpc?1 with m/T0 = 107. For the larger kmode, we see some oscillations after horizon crossing. These are due to thepresence of some mass as seen when comparing with the smooth masslessneutrino case.468.4. Warm Dark Matter Datainto real space to see how the WDM is collapsing. In the following chapter,we will discuss how we perform this real space transform.47Chapter 9Data Analysis RoutinesThere are two different real space transforms that we wish to look at. Wewill begin by examining ? and ? in real space. Afterwards, we will lookat the distribution functions in real space. We may use the informationobtained from the perturbations to verify the distribution function.9.1 Transfer FunctionThere is one important detail that has not been discussed yet that will al-low us to perform these transforms. Upon looking at the equations that weevolve, all terms are proportional, to first order, to a perturbation. As such,we may multiply all the perturbations by some function of k, without affect-ing the equations. What this allows us to do is rewrite our perturbationsas T (k)F (k), where T (k) is the quantity calculated in the previous sections(?, ?, ?, etc) and called the transfer function and F (k) is some function ofk[3]. We should note that despite T (k) being different for each perturba-tion, F (k) must be the same for all perturbations. What this decompositionallows us to do is set the initial distribution of matter. Since ? is constantat early times for all k, if we set F (k) to the Fourier transform of some f(x),?(x, ?0) will have the same shape as f(x). We may also use F (k) to set the489.2. Real Space Perturbationsproperties of our system. By modifying the magnitude of F (k), this willallow us to set either the magnitude of our initial conditions, or fix someend point quantity to a desired value. In our case, we chose to set ?8, theRMS overdensity in a sphere of 8 Mpc[9], such that it was comparable withthe real Universe, where?28 = ?[ d3k(2?)3 ?(~k)W8(~k)]2?, (9.1)where W8(~k) is the Fourier transformed 3 dimensional top hat function fora 8 Mpc sphere.9.2 Real Space PerturbationsWe can now discuss the real space perturbations. We wish to work in threedimensional real space. Usually a Fast Fourier Transform (FFT) algorithmis used to numerically calculate Fourier transforms. Although these canbe adapted to work in more then one dimension, we instead use the factthat our system is spherically symmetric to utilize a Hankel transform[11].The Hankel transform is similar to a Fourier transform but it utilizes theBessel functions as a basis instead of the complex exponential. The Hankeltransform is also only valid for functions which only depend on r, whetherit be in two dimensions or three. The transform is defined as [11]H?(k) =? ?0h(r)J?(kr)rdr, (9.2)499.2. Real Space Perturbationswhere J? is the Bessel function of first kind of order ?, with ? ? ?1/2. Theinverse is simplyh(r) =? ?0H?(k)J?(kr)kdk. (9.3)The integer values of ? represent the 2D transforms while the half integersare the 3D transforms. These are the ones we are interested in. We will belooking at the ? = 1/2 transform. The corresponding Bessel function isJ1/2(kr) =?2kr?sin(kr)kr . (9.4)We cannot simply Hankel transform our calculated values with equation(9.3) as our transfer functions where specifically calculated in Fourier space.As such, we wish to equate the Fourier transform to the Hankel transform.We will use f(r) and F (k) to denote the Fourier functions, and h(r) andH?(k) for the Hankel functions. If we perform the Fourier transform for aradially symmetric system, we obtainF (k) = 4?? ?0f(r)sin(kr)kr r2dr. (9.5)Comparing F (k) to H1/2(k), we can see that if we seth(r) = f(r)(2?)3/2?r?k, (9.6)then both F (k) and H1/2(k) are equal. The same analysis can be done forthe inverse transforms, where we find that we must take the inverse of r andk to obtain equivalence.509.2. Real Space PerturbationsThere is one subtlety that we now must consider before applying thisprocedure to ?. When calculating ?, we used the fact that in Fourier space,the divergence in real space can be replaced by ?ik in Fourier space. Thisis a property of the Fourier transform and does not hold for the Hankeltransform. As such, we must obtain a new method of transforming ?(k)using the Hankel transform. Beginning with the inverse Fourier transform,using the substitution procedure derived above, we will take the divergenceon both sides to obtain?r ? ~v(r) = 4?? ?0?r ? (~v(k)sin(kr)kr k2dk). (9.7)Performing the derivative on the RHS, we obtain?(r) = 4?? ?0~v(k)(sin(kr)kr2 +cos(kr)r)k2dk. (9.8)We can compare this expression with J3/2 and J1/2 to see that?(r) = (2?)3/2? ?0~v(k)k1/2r1/2(J1/2kr ? J3/2)kdk. (9.9)We can now transform our perturbations from k space into real space usinga Hankel transform.Unfortunately, due to the numerical restrictions discussed in the previouschapter, we do not cover a sufficient k range to apply the Hankel algorithmwith enough resolution to observe interesting perturbations. The maximumresolution is proportional to 2?/kmax and as we saw earlier, k < 3Mpc?1,and as such, x can only be sampled roughly every 2 Mpc, and we wish for519.3. Distribution Function in Real Spaceperturbations that are between 10-100 kpc.9.3 Distribution Function in Real SpaceThe other real space quantity of interest is the distribution function. Toobtain it, we must transform ?(~k, ?, qn?). Unlike the perturbations that wetransformed earlier, ? is not spherically symmetric due to the presence ofPl(k? ? n?). As such, we cannot perform the Hankel transform. We will needto perform a full 3D Fourier transform. Before we attempt to numericallyperform this, let us see if we may simplify the integrand.We begin by defining the three vectors of interest, ~k, ~q, and~x. We havesome freedom on how to define these vectors. Due to azimuthal symmetry,we know that our final solution must depend on only x? ? n?. As such, we canset the direction of x?, and place n? in the same plane. The simple way ofdoing this is by setting x? in the z? direction and place n? in the xz plane. Weleave k? as a free vector, as we will be integrating over all of k-space in ourFourier transform. The final form of our three vectors arex? = z?, (9.10a)n? = sin(?)x?+ cos(?)z?, (9.10b)k? = cos(?) sin(?)x?+ sin(?) sin(?)y? + cos(?)z?, (9.10c)where ? = x? ? n?.We may perform the appropriate dot products and rewrite our Fourier529.3. Distribution Function in Real Spacetransform as?(~x, ~q, ?) =?e?ikx cos???l=0(?i)l(2l + 1)?l(~k, ~q, ?)Pl(sin(?) cos(?) sin(?) + cos(?) cos(?))d?d? sin(?)F (k)k2dk,(9.11)where F (k) is again the k space representation of the initial conditions.We may now perform the angular integral analytically. We obtain?(~x, ~q, ?) = 4?? ??l=0(2l + 1)Pl(?)1(kx)l+1??l?m=0,even?m(kx) sin(kx)(?1)m/2+l?m=1,odd?mkx cos(kx)(?1)(m+1)/2??F (k)k2dk,(9.12)where ?m(x) are the reverse Bessel polynomials and are defined as[6]?n(x) =n?m=0(2n?m)!(n?m)!m!xm2n?m , (9.13a)?n(x) =n?m=0Cnmxm2n?m . (9.13b)The closed form of this expression was found by inputting the numer-ical results of integrating over Pl into the Online Encyclopedia of IntegerSequences [13].Having determined an analytic equation, we must now consider howto numerically compute this integral. The analytic solution works well for539.3. Distribution Function in Real Spacekx > 0.3 but begins to have numerical issues for small values of kx. ByTaylor expanding the solution, we obtain?(~x, ~q, ?) ?2?? ??l=0Pl(?)(2l + 1)??n=0[kxn(?1)nCnn+??m=1kxn+2m(?1)n+mCnn?ml=1 2l(2m+ 1 + 2l)]k2F (k)dk,(9.14)which converges rapidly for small values of kx. We will hence use the analyticsolution for kx > 0.3 and the Taylor expansion for kx ? 0.3. We can nownumerically perform the k space integral to obtain the distribution functionin real space. To perform this integral, we integrate starting at a small valueof k and end at a large value of k instead of integrating from 0 to ?. Thiscan be done since for a small k, F (k) approaches 0, and the Taylor expansionportion approaches one. This gives us a k2 dependance overall, causing smallvalues to contribute little to the integral. As for the high k regime, we seethat the highest power in k is proportional to kF (k)(sin(k) or cos(k)). For aspherically symmetric top hat of radius R, F (k) = 4? sin(kR)?Rk cos(kR)k3 . Assuch, for large k values, the integrand is proportional to k?1 or lower powers,and, we may stop at some maximum k without incurring large errors. Toobtain the distribution function, we simply substitute ? we calculate intoequation (5.3).In order to determine that this routine is valid, one may compare thedistribution functions with the perturbations. There should be a peak inthe distribution function related to ? and the peak should be displaced from549.3. Distribution Function in Real Space~q = 0 by ~v.55Chapter 10ConclusionWith the amount of matter that is dark matter, it is important to fullyunderstand the behaviour of dark matter in structure formation. Using theConformal Newtonian gauge, we have obtained a system of equations toevolve a variety of particle types through time. Using such equations, wehave been able to create a WDM+pf system and evolve it. Due to numericalcomplications we have a limited range of k values that may be calculated.Also, through the study of test cases, we have observed that some numericalissues exist in our solver. After obtaining accurate transfer functions, wewould be able to transform our data into real space.Once we have real space data, it would be possible to use it to generateinitial conditions for an N-body simulation. Our current set of equationsis linear, and as such, once the over density approaches 1, the equationsbecome inaccurate. Hence, the N-body simulation is necessary to obtainthe later evolution. We should note that due to the transitional propertiesof WDM, we must evolve it until it is non-relativistic, before inputting it intoan N-body simulation. This N-body simulation would allow us to observehow the WDM collapses and give us new insight in structure formation inour Universe.56Bibliography[1] Planck Collaboration. Planck 2013 results. XVI. cosmological parame-ters. Mar 2013, 1303.5076.[2] A. Del Popolo. Non-baryonic dark matter in cosmology. 2013,1305.0456.[3] S. Dodelson. Modern Cosmology. Academic Press, San Diego: London:Burlington, first edition, 2003.[4] Adrienne L. Erickcek and Kris Sigurdson. Reheating effects in thematter power spectrum and implications for substructure. Phys. Rev.D, 84:083503, Oct 2011.[5] J. R. Gott, III and J. E. Gunn. The coma cluster as an x-ray source:Some cosmological implications. Astrophys. J., 169:L13, Oct 1971.[6] HL Krall and Orrin Frink. A new class of orthogonal polynomials:The bessel polynomials. Transactions of the American MathematicalSociety, 65(1):100?115, 1949.[7] Hannu Kurki-Suonio. Introduction to cosmological perturbation theory.Lecture notes, Apr 2011.57Bibliography[8] Julien Lesgourgues and Thomas Tram. The cosmic linear anisotropysolving system (class) IV: efficient implementation of non-cold relics.JCAP, 1109:032, 2011, 1104.2935.[9] Chung-Pei Ma and Edmund Bertschinger. Cosmological perturbationtheory in the synchronous versus conformal newtonian gauge. Astro-phys. J., 1994, astro-ph/9401007.[10] Shmuel Nussinov. Some aspects of new cdm models and cdm detectionmethods. Mod.Phys.Lett., A24:2213?2223, 2009, 0907.3866.[11] A. Papoulis. Systems and Transforms With Applications in Op-tics. McGraw-Hill series in systems science. Malabar, Florida: RobertKrieger Publishing Company, 1968.[12] William H. Press, Saul A. Teukolsky, William T. Vetterling, andBrian P. Flannery. Numerical Recipes in C. Press Syndicate of theUniversity of Cambridge, Cambridge: New York: Port Melbourne, sec-ond edition, 1988.[13] N. J. A. Sloane. Online encyclopedia of integer sequences, Mar 2013.[14] Shruti Thakur and Anjan A Sen. Can structure formation distinguish? cdm from non-minimal f(r) gravity? 2013, 1305.6447.[15] David H. Weinberg, James S. Bullock, Fabio Governato, Rachel Kuziode Naray, and Annika H. G. Peter. Cold dark matter: controversies onsmall scales. 2013, 1306.0913.58Bibliography[16] Inc. Wolfram Research. Mathematica. wolfram Research, Inc, Cham-paign, Illinois, version 8.0 edition, 2010.59


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