{'fraction_non_alphanumeric': 0.03522322024524668, 'fraction_numerical': 0.00972595686252074, 'mean_word_length': 4.5536962507943235, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Cosmological Perturbation Theory and Structure Formation

These lecture notes discuss several topics in the physics of cosmic structure formation starting from the evolution of small-amplitude fluctuations in the radiation-dominated era. The topics include relativistic cosmological perturbation theory with the scalar-vector-tensor decomposition, the evolution of adiabatic and isentropic initial fluctuations, microwave background anisotropy, spatial and angular power spectra, the cold dark matter linear transfer function, Press-Schechter theory, and a brief introduction to numerical simulation methods.


Introduction
The advent of precision cosmology with measurements of the Cosmic Microwave Background (CMB) radiation has brought a new focus on the physics of the z > 1100 universe. The connection between two fossil relics of the early universe-CMB fluctuations and galaxies-offers the possibility for improved understanding of both as well as tighter constraints on cosmological parameters. These lectures present an overview of the physics of structure formation starting after the inflationary epoch, going through recombination, and extending to the present day.
Cosmological perturbation theory is often regarded as a highly technical subject beyond the scope of graduate courses in cosmology. However, these lectures may suggest how, with minor simplification, it can be integrated into a treatment of large scale structure. The potential benefits are many-a detailed understanding of CMB anisotropy, a quantitative explanation for the key features of the cold dark matter transfer function, and an extension of the intuition developed from the Newtonian theory of large scale structure to the physics of the post-inflationary universe. The aim of the current notes is more modest, but it is hoped that they provide some direction for those wishing to delve more deeply into cosmological perturbation theory and structure formation.

Overview of Relativistic Cosmological Perturbation Theory
The goal of cosmological perturbation theory is to relate the physics of the early universe (e.g. inflation) to CMB anisotropy and large-scale structure and to provide the initial conditions for numerical simulations of structure formation. The physics during the period from the end of inflation to the beginning of 1 2 Edmund Bertschinger nonlinear gravitational collapse is complicated by relativistic effects but greatly simplified by the small amplitude of perturbations. Thus, an essentially complete and accurate treatment of relativistic perturbation evolution is possible, at least in the context of simple fluctuation models like inflation. Numerous articles have been written on this topic. See [1,2] for references and for more detailed treatments similar in style to the present lecture notes.
The starting point for cosmological perturbation theory is the metric of a perturbed Robertson-Walker spacetime, Spatial coordinates take the range 1 ≤ i, j ≤ 3; x i (or x for all three) is a comoving spatial coordinate; τ is conformal time; a(τ ) is the cosmic expansion scale factor; units are chosen so that the speed of light is unity; and γ ij ( x ) is the 3-metric of a maximally symmetric constant curvature space. Conformal time is related to the proper time measured by a comoving observer (i.e. one at fixed x) by dt = a(τ )dτ . For cosmological perturbation theory it is more convenient than proper time. The metric perturbations are given by h µν = g (1) µν /a 2 . The Robertson-Walker model is characterized by one function of time, a(τ ), and one constant, the curvature constant K. The expansion factor obeys the Friedmann equation, ȧ a 2 = 8π 3 Ga 2ρ (a) − K where a dot denotes d/dτ andρ(a) is the total mean energy density. Note that the usual Hubble parameter is H =ȧ/a 2 . The spatial part of the Robertson-Walker metric takes a simple form in spherical coordinates (χ, θ, φ): γ ij dx i dx j = dχ 2 + r 2 (χ) dθ 2 + sin 2 θdφ 2 . ( The angular radius r(χ) takes one of three forms depending on the sign of K:

Scalar-Vector-Tensor Decomposition
In linear perturbation theory, the metric perturbations h µν are regarded as a tensor field residing on the background Robertson-Walker spacetime. As a symmetric 4 × 4 matrix, h µν has 10 degrees of freedom. Because of the ability to make continuous deformations of the coordinates, 4 of these degrees of freedom are gauge-(coordinate-) dependent, leaving 6 physical degrees of freedom. A proper treatment of cosmological perturbation theory requires clear separation between physical and gauge degrees of freedom.
The metric degrees of freedom in linear pertubation theory were classified by Lifshitz in 1946 [3]. Lifshitz presented the scalar-vector-tensor decomposition of the metric. It is based on a 3 + 1 split of the components, which we rewrite as follows: Here, γ ij is the matrix inverse of γ ij . The trace part of h ij has been absorbed into φ so that S ij has only 5 independent components. With the space-time (3 + 1) split, we use γ ij (or γ ij ) to lower (or raise) indices of spatial 3-vectors and tensors. For convenience, spatial derivatives will be written using the 3-dimensional covariant derivative ∇ i defined with respect to the 3-metric γ ij ; it is a three-dimensional version of the full 4-dimensional covariant derivative presented in general relativity textbooks. For example, if K = 0, we can choose Cartesian coordinates such that γ ij = δ ij and ∇ i = ∂/∂x i .
The scalar-tensor-vector split is based on the decomposition of a vector into longitudinal and transverse parts. For any three-vector field w i ( x ), we may write The curl and divergence are defined using the spatial covariant derivative, e.g.
The longitudinal/transverse decomposition is not unique (e.g. one may always add a constant to w i ) but it always exists. The terminology arises because in the Fourier domain w i is parallel to the wavevector while w ⊥ i is transverse (perpendicular to the wavevector). Note that w i = ∇ i φ w for some scalar field φ w . Thus, the longitudinal/transverse decomposition allows us to write a vector field in terms of a scalar (the longitudinal or irrotational part) and a part that cannot be obtained from a scalar (the transverse or rotational part).
A similar decomposition holds for a two-index tensor, but now each index can be either longitudinal or transverse. For a symmetric tensor, there are three possibilities: both indices are longitudinal, one is transverse, or two are transverse. These are written as follows: where The first term in equation (8) is a longitudinal vector while the second term is a transverse vector. The divergence of the doubly-transverse part, S T ij , is zero.

Edmund Bertschinger
For a traceless symmetric tensor, the doubly and singly longitudinal parts can be obtained from the gradients of a scalar and a transverse vector, respectively: Although we will not prove it here (see [1] for the details), under an infinitesimal coordinate transformation S ij and S ⊥ ij can change while S T ij is invariant. Similarly, the longitudinal and transverse parts of w i = h 0i are both gauge-dependent. Now we have the mathematical background needed to perform the scalartensor-vector split of the physical degrees of freedom of the metric. The "tensor mode" represents the part of h ij that cannot be obtained from the gradients of a scalar or vector, namely S T ij . The tensor mode is gauge-invariant and has two degrees of freedom (five for a symmetric traceless 3 × 3 matrix, less three from the condition γ jk ∇ k S T ij = 0). Physically it represents gravitational radiation; the two degrees of freedom correspond to the two polarizations of gravitational radiation. Gravitational radiation is transverse: a wave propagating in the zdirection can have nonzero components h xx − h yy and h xy = h yx but no others. The tensor mode S T ij behaves like a spin-2 field under spatial rotations. Note that my definition of the tensor mode strain S ij here differs by a factor of 2 from the conventional h T T ij in [4]. The "vector mode" behaves like a spin-1 field under spatial rotations. It corresponds to the transverse vector parts of the metric, which are found in w ⊥ i and S ⊥ ij . Each part has two degrees of freedom. Although we will not prove it here (see [1] for the details), it is possible to eliminate two of these degrees of freedom by imposing gauge conditions (coordinate conditions). The synchronous gauge of Lifshitz [3] is one popular choice, with w ⊥ i = w i = 0. Another choice is the "Poisson" or transverse gauge [1] which sets S ⊥ ij = 0. In either case, there are two physical degrees of freedom and they correspond physically to gravitomagnetism. Although this effect is less well known than gravitational radiation or Newtonian gravity, it produces magnetic-like effects on moving and spinning masses, such as the precession of a gyroscope in the gravitational field of a spinning mass (Lense-Thirring precession). This phenomenon has not yet been discovered experimentally but should be measured by the Gravity Probe B satellite to be launched in 2002 [5].
The "scalar mode" is spin-0 under spatial rotations and corresponds physically to Newtonian gravitation with relativistic modifications. The scalar parts of the metric are given by φ, ψ, w i , and S ij . Any two of these may be set to zero by means of a gauge transformation. One popular choice is w i = S ij = 0, also known as the conformal Newtonian gauge. It corresponds to the scalar mode in the transverse gauge, defined by the gauge conditions This gauge in linearized general relativity is the gravitational analogue of Coulomb gauge in electromagnetism, where the magnetic vector potential is transverse, In the gravitational transverse gauge, w i is transverse and S ij is doubly transverse. This gauge is convenient for developing intuition although not necessarily the best for computation. The variables (φ, ψ, w ⊥ i , S T ij ) correspond to gauge-invariant variables introduced by Bardeen [6] as linear combinations of metric variables in other gauges.
In linear perturbation theory, the scalar, vector, and tensor modes evolve independently. The vector and tensor modes produce no density perturbations and therefore are unimportant for structure formation, although they do perturb the microwave background. In the remainder of these lectures, only the scalar mode will be considered.
The Einstein equations give the equations of motion for the metric perturbations in terms of the energy-momentum tensor, the source of relativistic gravity. Here we consider the case of a perfect fluid (or several perfect fluid components combined), for which where ρ and p are the proper energy density and pressure in the fluid rest frame and V µ is the fluid 4-velocity. We split T µν into time and space components as we did for the metric. It is convenient to use mixed components because the conformal factor a 2 then cancels, yielding T i j = [p(τ ) + δp( x, τ )] δ i j , whereρ(τ ) andp(τ ) are respectively the energy (or mass) density and pressure of the Robertson-Walker background spacetime, v i = γ ij dx j /dτ is the fluid 3-velocity (assumed nonrelativistic), and W is a velocity potential. We are assuming that the velocity field is irrotational or, if it is not, that the matter fluctuations are linear so that the longitudinal and transverse velocity components evolve independently. This is a good approximation prior to the epoch of galaxy formation.
We are also making another approximation, namely that the shear stress (the non-diagonal part of T i j ) is negligible compared with the pressure. This is a bad approximation for relativistic neutrinos, whose free-streaming leads to a non-isotropic momentum distribution after neutrino decoupling at a temperature of about 1 MeV. (A similar but smaller effect occurs for photons after they decouple at a temperature of about 0.3 eV.) However, the effect of this error is modest during the radiation-dominated era (photons contribute more pressure than neutrinos, and are isotropic) and is unimportant during the matterdominated era, when the nonrelativistic matter density is the major source of gravity. We introduce this approximation here in order to simplify the treatment; full calculations include the effects of neutrino shear stress [2].
When the stress tensor is isotropic, the Einstein equations in transverse gauge imply that the two scalar potentials φ and ψ are equal [1] and the metric takes the simple form This is a simple cosmological generalization of the standard weak-field limit of general relativity, which is recovered by setting a = 1 and γ ij = δ ij . Given the stress-energy tensor, we can now obtain field equations for φ( x, τ ) using the Einstein equations. The only independent equation for the Robertson-Walker background is the Friedmann equation (2). For the perturbations, we obtain and There are more equations than in Newtonian gravity because the Einstein equations have local energy-momentum conservation built into them. Equation (14) is remarkably like the Poisson equation of Newtonian gravity, especially when one realizes that the factor of a 2 is needed to convert the Laplacian from comoving to proper spatial coordinates. Otherwise there are three differences. First, the spatial curvature adds a scale to the Laplacian in a curved space. Second, the density source is δρ = ρ −ρ rather than ρ because φ represents a metric perturbation on top of the background Robertson-Walker spacetime. Third, momentum density (through W ) is a source for gravity in general relativity. However, this term is comparable with the density term only on length scales comparable to or larger than the Hubble length H −1 .
Two of our field equations are unfamiliar from Newtonian dynamics. Equation (15) is reminiscent of the Zel'dovich approximation [7] in that the velocity potential is proportional to the gravitational potential in the linear regime when φ( x, τ ) = φ 0 ( x )D φ (τ ). However, it is a more general result whose Newtonian limit may be derived by combining the continuity equation with the time derivative of the Poisson equation.
Equation (16) is more surprising, as it gives the gravitational potential through a second-order time evolution equation as opposed to the action-at-adistance Poisson equation (14). The reader may well wonder about the compatibility of the two equations as well as the physical validity of the Poisson equation in general relativity, but equations (14)-(16) are all correct. See [1] for a full discussion. In brief, causality is restored by the tensor mode. The time evolution equation for the potential can be regarded as arising from local energy-momentum conservation (which is built into the Einstein equations) combined with the Poisson equation. Indeed, an equation similar to equation (16) can be derived in Newtonian gravity from the fluid and Poisson equations. Its utility arises from the fact that when pressure perturbations are gravitationally unimportant (i.e. when matter and/or dark energy dominate over radiation), the time evolution of perturbations can be obtained from a single differential equation depending only on the background curvature, pressure, and density.

Perfect Fluid Model
To conclude this section, we present a simple model for understanding CMB anisotropy and large scale structure. We reduce the matter and energy content Cosmological Perturbation Theory and Structure Formation 7 of the universe to a vacuum energy (with ρ v = Λ/8πG contributing to the mean expansion rate but not to the fluctuations) and two fluctuating components, CDM and radiation (photons and neutrinos, the former coupled to baryons by electron scattering until recombination). The CDM and radiation are treated as perfect fluids, ignoring the free-streaming of neutrinos and the diffusion and streaming of photons during and after recombination.
The fluid equations follow from local energy-momentum conservation, ∇ µ T µν = 0. For a perfect fluid, with the scalar mode only (no vorticity and no gravitational radiation), the fluid equations in the metric of equation (13) are and These are familiar from Newtonian fluid dynamics, with some extra terms. The damping (ȧ/a) terms arise from Hubble expansion and the use of comoving coordinates. The ∂ τ φ term in the continuity equation (17) arises from the deformation of the spatial coordinates, i.e. the φδ ij contribution to the metric.
(The flux term in eq. 17 describes transport relative to the comoving coordinate grid, which is deforming when ∂ τ φ = 0.) The energy flux (or momentum density, they are equal) includes pressure in relativity because of the pdV work done by compression. Equations (17) and (18) apply separately to perfect fluid component.

Adiabatic and Isocurvature Modes
This section discusses the evolution of perturbations from the end of inflation (taken to be at τ = 0 for all practical purposes) through recombination using a two-fluid approximation. We approximate the matter and energy content in the universe as being two perfect fluids: Cold Dark Matter (a zero-temperature gas) and radiation (photons coupled to electrons and baryons before recombination, plus neutrinos which are approximated as behaving like photons). The net energy density perturbation is thus whereρ c ∝ a −3 and, if we neglect the contribution of baryons,ρ r ∝ a −4 . We could add a third component for the baryons, but this complicates the presentation without adding essential new behavior. In the interest of economy, we leave baryons out. We define the relative density contrasts δ c ≡ δρ c /ρ c and δ r ≡ δρ r /ρ r . The total pressure is simply p( x, τ ) = 1 3ρ r (a)(1 + δ r ) + p v where p v = −ρ v is the spatially constant negative pressure of vacuum energy (cosmological constant), should any be present.
As we noted in the previous section, density perturbations couple gravitationally only to the scalar mode of metric perturbations, and the scalar mode cannot generate transverse vector fields. As a result, the peculiar velocity fields of CDM and radiation are longitudinal and are fully characterized by their potentials W c and W r for CDM and radiation, respectively.
As a final simplifying assumption, we will assume that we are studying effects on distance scales less than the curvature length |K| −1/2 = |Ω 0 −1| −1/2 H −1 0 , which is at least 5 Gpc and possibly infinitely large.
With these assumptions, we can linearize the fluid equations (17) and (18) for δ 2 c ≪ 1 and δ 2 r ≪ 1 to obtain These equations are easy to understand. The ∂ τ φ terms have the same origin as in equation (17), namely the deformation of the comoving coordinate system. The velocity potential terms are related to the velocity divergence, The radiation density fluctuations grow at 4/3 the rate of matter fluctuations to compensate for the fact that the photon number density is proportional to ρ 3/4 r . As expected, the velocity equations show that matter velocities redshift (for φ = 0) with v c ∝ a −1 while radiation velocities (e.g. the speed of light) do not redshift. (Momenta redshift ∝ a −1 but the speed of light is unaffected by expansion.) Finally, radiation pressure gives rise to an acoustic restoring force − ∇p r = −(1/3)ρ r ∇δ r which shows up as the (1/4)δ r term in the radiation longitudinal velocity equation.
Equations (20) are coupled by gravity, which can be determined from the Poisson equation (14), which becomes We can solve equations (20) and (21) most easily by expanding the spatial dependence in eigenfunctions of the Laplacian ∇ 2 . In flat space these are plane waves exp(i k · x ), and this is a good approximation even if the background is curved, provided k ≫ |K| 1/2 . This approximation is valid for almost all applications except large angular scale CMB anisotropy. Note that k is a comoving wavevector; the physical wavelength is 2πa/k.
By expanding the spatial dependence in plane waves, we now have a fourthorder linear system of ordinary differential equations in time. We are free to define linear combinations of the variables (δ c , δ r , W c , W r ) as our fundamental variables. In terms of the mechanisms for generating primeval perturbations, the most natural variables are the metric perturbation φ and the specific entropy where is the effective one-fluid sound speed of the matter and radiation. (Acoustic signals actually do not propagate at this speed; they propagate with speed 3 −1/2 through the radiation and they do not propagate at all through CDM.) For a two-component radiation plus CDM universe, the solution of the Friedmann equation (2) gives The radiation-dominated era ends and matter-dominated era begins at redshift 1 + z eq = a −1 eq = 2.5 × 10 4 Ω m h 2 . A cosmological constant has no significant effect on y(τ ) provided that, during the times of interest, the vacuum energy density is much less than the radiation or cold dark matter density.
The fluid and Poisson equations can now be combined to give a pair of second-order ordinary differential equations, It is interesting to note that φ and η evolve independently aside from the source term each provides to the other. The coupling implies that entropy perturbations are a source for the growth of gravitational perturbations and vice versa. For a given wavenumber k, there are two key times in the evolution of φ and η: the sound-crossing time τ ≈ π/(kc s ) and the time of matter-radiation equality at τ ≈ τ e . [To be precise, a = a eq at τ /τ e = 2( √ 2 − 1).] For τ ≪ τ e , φ and η decouple and they both have solutions that decay rapidly with τ as well as solutions that are finite as τ → 0. The latter are conventionally called "growing modes" while the former are called decaying. The growing mode solution in the radiation era τ ≪ τ e is where ω ≡ k/ √ 3 is the phase of acoustic waves in the radiation fluid, C = 0.5772 . . . is the Euler constant, and Ci(x) is the cosine integral defined by In equation (26) we have neglected all terms O(τ /τ e ) aside from the lowest-order contribution of I( k ) to φ.
The general growing-mode solution contains two k-dependent integration constants, A( k ) and I( k ), which represent the initial (τ → 0) amplitude of φ and η, respectively. Conventionally they are called the "adiabatic" (so-called because the entropy perturbation η = 0 at τ = 0, although "isentropic" would be more appropriate) and "isocurvature" (so-called because φ = 0 at τ = 0) mode amplitudes. In addition, there are decaying modes whose amplitude at τ = τ e is negligible unless the initial conditions were fine-tuned or there are external sources such as topological defects to continually regenerate the decaying mode. The solution of equations (26) is valid only for inflation or other mechanisms that generate fluctuations exclusively in the early universe, but not for topological defects.
As a result of reheating, inflation produces fluctuations with spatially constant ratios of the number densities of all particle species. For CDM and radiation, for example, ρ 3/4 r /ρ c = constant hence I = 0 and thus η = 0 until radiation pressure forces begin to separate the radiation from the CDM, which takes a sound-crossing time or τ ∼ ω −1 . Thus, although inflation produces the adiabatic mode of initial fluctuations, the differing equations of state p(ρ) for the two-fluid system lead to entropy perturbations proportional to the initial adiabatic mode amplitude A( k ) (see the second of eqs. 26).
Physically, the number density of particles of all species are proportional to each other in the adiabatic mode as long as there has been too little time for particles to travel a significant difference via their thermal speeds. However, the differing thermal speeds of different particle species (e.g. photons and cold dark matter) eventually cause the particle densities to evolve differently, so that the number density ratios are no longer constant and the entropy perturbation is no longer zero.
First-order phase transitions, on the other hand, change the equation of state without moving matter or energy on macroscopic scales, hence δρ r + δρ c = θ c = θ r = 0 initially implying φ = 0 initially but η = 0. This isocurvature mode is characterized by the initial entropy perturbation amplitude I( k ). As the first of equations (26) shows, the isocurvature mode generates curvature fluctuations when the universe becomes matter-dominated.
The isocurvature CDM inflationary model has a fine balance such that δ r + yδ c = 0 initially. This may involve nonlinear fluctuations of δ c at early times, but δ r is negligible early in the radiation-dominated era y ≪ 1. For this reason, the isocurvature mode was originally called isothermal [8].
On small scales, kτ ≫ 1, equations (25) may be solved using a WKB approximation. For our purposes it suffices to take the ωτ ≫ 1 limit of equations (26): The radiation density fluctuation scales as which oscillates in proportion to cosωτ for the adiabatic mode and grows in proportion to y ∼ τ /τ e for the isocurvature mode.
Note from equation (28) that the isocurvature mode induces oscillations in the potential with a π/2 phase shift relative to the adiabatic mode. This shift translates into a shift in the positions of the acoustic peaks in the CMB angular power spectrum.

The Adbiabatic Growing Mode as a Standing Wave
It is instructive to consider the general adiabatic perturbation in the radiationdominated era including the decaying mode. This problem is easy to analyze and shows the role of acoustic waves in the radiation-dominated era.
The two linearly independent solutions of the first of equations (25) for y ≪ 1 and adiabatic initial conditions (η = 0) are where the subscripts label the growing and decaying modes, respectively. They are the real and imaginary parts of For a given spatial harmonic, exp(i k· x )φ ω represents a wave traveling in the k-direction. Thus, the growing mode φ g = 1 2 (φ ω +φ * ω ) represents a superposition of waves traveling along k and − k, i.e. a standing wave The waves travel at the sound speed ω/k = 1/ √ 3 in the radiation fluid. The entropy perturbation induced by a traveling wave potential perturbation is easily calculated in the radiation-dominated era using the Green's function method applied to the second of equations (25) with φ = φ ω for the source term. The result is The real part of this equation gives the A( k ) contribution to η in equations (26).

CMB Anisotropy
This section presents a simplified treatment of CMB anisotropy, with the aim of highlighting the essential physics without getting lost in all of the details. More complete treatments are found in [9,10]. The microwave background radiation is fully described by a set of photon phase space distribution functions. Ignoring polarization (a few percent effect), all the information is included in the intensity or in f ( x, p, τ ) where f d 3 p is the number density of photons of conjugate momentum p at position and time ( x, τ ). The conjugate momentum is related to the proper momentum measured by a comoving observer, p/a, so that p = constant along a photon trajectory in the absence of metric perturbations.
Remarkably, despite metric perturbations and scattering with free electrons, the CMB photon phase space distribution remains blackbody (Planckian) to exquisite precision: where T 0 = 2.728 K is the current CMB temperature and ∆( x, n, τ ) is the temperature variation at position ( x, τ ) for photons traveling in direction n.
The phase space density is blackbody but the temperature depends on photon direction as a result of scattering and gravitational processes occurring along the line of sight.
The phase space density may be calculated from initial conditions in the early universe through the Boltzmann equation where the right-hand side is a collision integral coming from nonrelativistic electron-photon elastic scattering. In the absence of scattering, the phase space density is conserved along the trajectories Note that the photon energy changes via gravitational redshift (the first term in dp/dτ ) or a time-changing potential, but in either case the energy of all photons is changed by the same factor. The photon direction changes by transverse deflections due to gravity; d n/dτ is twice the Newtonian value, a result wellknown in gravitational lens theory [11]. Because ∂f /∂ n = 0 for the unperturbed background (the anisotropy arises due to the metric perturbations), gravitational lensing affects the CMB anisotropy only through nonlinear effects on small scales [12]. The procedure for computing CMB anisotropy is to linearize the Boltzmann equation assuming ∆ 2 ≪ 1 [13]. Until recently, the temperature anisotropy ∆ was expanded in spherical harmonics and the Boltzmann equation was solved as a hierarchy of coupled equations for the various angular moments [2]. In 1996, Seljak and Zaldarriaga [9] developed a much faster integration method call CMBFAST based on integrating the linearized Boltzmann equation along the observer's line of sight before the angular expansion is made: where τ 0 is the present conformal time, χ is the radial comoving coordinate of equation (3), subscript "ret" denotes evaluation using retarded time τ = τ 0 − χ, δ γ is the relative density fluctuation in the photon gas (δ γ = δ r in our two-fluid approximation), v e is the mean electron velocity (i.e. the baryon velocity, which equals v r in our two-fluid approximation), σ T is the Thomson cross section, and the Thomson optical depth is We have left out small terms ("pol. term") due to polarization and the anisotropy of Thomson scattering in equation (36). The terms in equation (36) are easy to understand. The exp(−τ T ) factor accounts for the opacity of electron scattering, which prevents us from seeing much beyond a redshift of 1100, where τ T ≈ 1. The CMB anisotropy appears to come from a thin layer called the photosphere, just like the radiation from the surface of a star.
The two gravity terms give the effective emissivity due to the photon energy change caused by a varying gravitational potential. For a blackbody distribution, if all photon energies are increased by a factor 1 + ǫ, the distribution remains blackbody with a temperature increased by the same factor. Thus, the lineof-sight integration of the fractional energy change of equation (35) translates directly into a temperature variation.
The terms proportional to the Thomson opacity an e σ T are the effective emissivity due to Thomson scattering. Bearing in mind that the photon-baryon plasma is in nearly perfect thermal equilibrium (the temperature variations are only about 1 part in 10 5 ), the photons we see scattered into the line of sight from a given fluid element have the blackbody distribution corresponding to the temperature of that element. Recalling that the energy density of blackbody radiation is ρ γ ∝ T 4 , we see that 1 4 δ γ is simply the fractional temperature variation of the fluid element. In other words, if we carried a thermometer to the photosphere and measured its reading in direction − n, it would read a fraction 1 4 δ γ different from the mean. Now the energy density is defined in the fluid rest frame, while the fluid is moving with 3-velocity v e , so the temperature measured by a stationary thermometer (with fixed x i ) is changed by a Doppler correction v e · n. (If the photosphere is approaching us, v e n, the radiation is blue-shifted.) Comparing the CMB with a star, when we look at the surface of a star we see the temperature of the emitting gas (to the extent that local thermodynamic equilibrium applies), corresponding to the 1 4 δ γ term. For ordinary stars the Doppler boosting and gravitational redshift effects are negligible, although they are appreciable for supernovae and white dwarfs, respectively. For the CMB, on the other hand, all four emission terms in equation (36) are comparable in importance.
It is instructive to approximate the photosphere as an infinitely thin layer by adopting the instantaneous recombination approximation, according to which the free electron fraction and hence opacity drops suddenly at τ rec , the conformal time of recombination at z ≈ 1100: Substituting into equation (36) and integrating by parts the gravitational redshift term gives the result first obtained in another form by Sachs and Wolfe in where we have dropped the argument τ 0 from ∆, χ rec = τ 0 − τ rec , and W γ is the velocity potential for the photons (W γ = W r in the two-fluid approximation discussed above). We could have used the last of equations (20) to write the sum of the intrinsic ( 1 4 δ γ ) and gravitational (φ) contributions as ∂ τ W r . The combination of derivatives (∂ τ + ∂ χ )W r is the time derivative along the past light cone; the results are evaluated at recombination.
While it is interesting that the primary anisotropy depends on the velocity potential derivative along the line of sight, there appears to be no fundamental significance to the simple dependence in equation (39), aside from the fact that CMB anisotropy is produced by departures from hydrostatic equilibrium. (In hydrostatic equilibrium, v r = 0.) If the radiation gas were in hydrostatic equilibrium in a static gravitational potential, there would be no primary CMB anisotropy. Indeed, it can be shown from purely thermodynamic arguments that for a photon gas in equilibrium. Sachs and Wolfe [14] showed that for adiabatic perturbations in the matterdominated era, on scales much larger than the acoustic horizon (ωτ ≪ 1 in eq. 26), the sum of the photospheric terms in equation (39) (the terms evaluated at recombination) is 1 3 φ. Thus, on angular scales much larger than 1 • (roughly the size of the acoustic horizon), the CMB anisotropy is a direct measure of the gravitational potential on the photosphere at recombination.

Angular Power Spectrum
The angular power spectrum gives the mean squared amplitude of the CMB anisotropy per spherical harmonic component. Thus, we expand the anisotropy in spherical harmonic functions of the photon direction: Observations of our universe give definite numbers for the a lm (with experimental error bars, of course). Theoretically, however, we can only predict the probability distribution of the a lm . For statistically isotropic fluctuations (i.e. having no preferred direction a priori), the a lm are random variables with covariance where the Kronecker delta's make the a lm 's uncorrelated (δ ll ′ = 1 if l = l ′ and 0 otherwise). The variance of each harmonic is given by the angular power spectrum C l ; rotational symmetry ensures that, theoretically, it is independent of m.
To calculate the angular power spectrum, we expand φ( x, τ ) (and the other variables) in plane waves (assuming that the background is flat, K = 0), where µ = − k · n/k with a minus sign because − n is the radial unit vector at the origin. Note that many cosmologists insert a factor (2π) −3 in the Fourier integral going from k-space to x-space; take heed of the conventions when they matter.
In equation (42) we have used the spherical wave expansion of a plane wave in terms of spherical Bessel functions j l (x) and Legendre polynomials P l (x). This gives where, in the instantaneous recombination approximation, The temperature expansion coefficient takes a simple form in terms of ∆ l : wherek is a unit vector in the direction of k.
To get the angular power spectrum now we must relate ∆ l ( k ) to the initial random field of potential or entropy fluctuations that induced the CMB anisotropy. Let us assume we have adiabatic fluctuations, for which φ( k, τ → 0) = A( k ). We define the CMB transfer function which depends only on the magnitude of k because the equations of motion have no preferred direction. To compute ensemble averages of products, we will need the two-point function for a statistically homogeneous and isotropic random field, whose variance defines the power spectrum: Note that many cosmologists use an alternative definition of the power spectrum with a factor (2π) 3 inserted on the right-hand side. The power spectrum so defined is greater by a factor (2π) 3 than the one given here. [The extra factors of (2π) 3 in eq. 42 are then important.] The definition presented here is conventional in field theory and it ensures that P A (k)d 3 k is the contribution per unit volume in k-space to the variance of A( x ). Because of this interpretation, the power spectrum is also known as the spectral density. See [15] for further discussion.
For scale-invariant fluctuations (equal power on all scales, as predicted by the simplest inflationary models), P A ∝ k −3 .

Edmund Bertschinger
Combining equations (41), (45), and (46) gives the formal expression for the angular power spectrum as an integral over the three-dimensional power spectrum: It is difficult to get a simple approximation to D l (k) that makes this analytically tractable, except in the limit of large angular scales where the intrinsic and gravitational contributions to ∆( n ) dominate, with ∆ ≈ 1 3 φ. Then D l = 1 3 j l (kχ rec ) and the integral may be evaluated analytically for P A (k) ∝ k n−4 with fixed n. When n = 1 (the scale-invariant spectrum), the result gives l(l+1)C l = constant or equal power on all angular scales. The phenomenon of acoustic peaks in l(l + 1)C l is due to the acoustic oscillations in φ and W γ which modify the potentials at recombination from their scale-invariant primeval forms.
This presentation has been based on the traditional Fourier space representation of the potential and other fields. The physical interpretation of the CMB two-point function is somewhat clearer in real (position) space, although a more detailed analysis based on Green's functions is needed before one can reap the rewards. For an introduction to the position-space perspective see [16].

Structure Formation in ΛCDM Models
This section summarizes briefly the steps needed to go from fluctuation evolution in the early universe to predictions of galaxy formation and large scale structure. The presentation will focus on techniques rather than model results. However, in order to simplify the discussion, we will focus on the favored cold dark matter family of models (with curvature and vacuum energy as options) with adiabatic density fluctuations from inflation.

CDM Linear Transfer Function
During the radiation-dominated era, CDM density fluctuations grow only logarithmically. Once the universe becomes matter-dominated at τ ∼ τ e (eq. 24), equation (16) shows that the potential becomes frozen in as long as vacuum energy (or other sources of pressure) and curvature are dynamically unimportant. Thus, from equation (14), the CDM density fluctuation amplitude grows linearly with a since a eq : We have used subscript "m" to include all cold matter including baryons, which cluster with the CDM after recombination. We have also included the velocity potential term 3(ȧ/a)W in equation (14) as part of δ m as it is in Bardeen's gauge-invariant density perturbation variable ǫ [6]. (The velocity potential term is anyway negligible on scales much smaller than the Hubble distance H −1 .) To fully specify the CDM density fluctuation spectrum we now need to know the gravitational potential, or equivalently the density fluctuations, at the end of the radiation era. There is a crucial physical scale, the Jeans length for the radiation gas at τ = τ e , (Including the effects of baryons decreases this by up to 15% from what is given here.) Perturbations of wavelength longer than 2π/k e have always been Jeans unstable, so that φ( k, τ e ) ≈ A( k ) for k ≪ k e . Note that only the acoustic horizon-not the causal horizon or Hubble scale-enters as a length scale in the CDM transfer function.
To get the amplitude for shorter-wavelength perturbations which have suffered acoustic damping (i.e. they were Jeans stable for a period of time), it is not enough to use the radiation-era solution for the gravitational potential in equation (49), because the CDM density fluctuations dominate over those of radiation before the universe becomes matter-dominated. Instead, we turn to equations (14) and (26) for ωτ ≫ 1 in the radiation era. For the adiabatic mode this gives The gravitational impulse created by the acoustic oscillations of the radiation fluid causes the CDM to oscillate slightly but as ωτ grows the cosine integral quickly dies and the CDM grows logarithmically. This growth is important for structure formation as it allows the CDM perturbations to dominate over the radiation perturbations enough for galaxies to form in less than a factor 1000 expansion since recombination when the large-scale radiation fluctuations were only a few parts in 10 5 . Combining results, we obtain an approximate solution for the CDM transfer function T m (k) giving the amplitude of CDM density fluctuations relative to the primeval potential A( k ): where Bardeen et al. [17] give a fitting formula for T m (k)/k 2 , which, when normalized to 1 as k → 0, is conventionally called the CDM transfer function. However, given the measurements of the gravitational potential made possible from the CMB, it is preferable to retain the normalization of equation (52) with the potential A( k ) factored out in equation (53). The CDM linear power spectrum follows by analogy with equation (47) for the potential. The result, for a matter-dominated universe, is Theoretical CDM power spectra compared with the deconvolved power spectrum of galaxies from the Las Campanas Redshift Survey [19]. The error bars are strongly correlated at large scales.
The constant-curvature scale-invariant spectrum of potential fluctuations P A ∝ k −3 was first proposed by Harrison [18] for its naturalness properties long before inflationary cosmology gave it a concrete basis. With this spectrum, the CDM density spectrum has the well-known limiting behavior The power spectrum is useful for calculating the variance of the smoothed density field. If we convolve δ( x, τ ) with a spherical window function whose Fourier transform is W (k) normalized so that W (0) = 1, then From this, we see that the CDM density fluctuations have logarithmically divergent power at small scales. Thus, structure formation in a CDM model proceeds by hierarchical or bottom-up clustering with small scales collapsing by gravitational instability before being incorporated into larger objects. However, the slowness of the logarithm means that the structure formation proceeds rapidly from small to large scales. Figure 1 compares the linear power spectra extrapolated to τ = τ 0 (i.e. redshift 0) for several flat CDM models. The curves are all normalized so that the density smoothed with a 8 h −1 Mpc spherical "tophat" has standard deviation σ 8 = 0.9 in agreement with observed galaxy clustering for a "linear bias factor" b = σ −1 8 = 1.1. The ΛCDM model with (Ω b = 0.05, Ω c = 0.30, Ω Λ = 0.65, h = 0.65) has the most power on large scales, possibly too much (although an acceptable fit can be found with slightly higher Ω c ). The SCDM model with (Ω b = 0.05, Ω c = 0.95, h = 0.50) has too little power on large scales. When combined with the CMB anisotropy, this provides a strong argument for dark energy (e.g. a cosmological constant) [20].
The curves labeled "BBKS" show the linear power spectra for the ΛCDM and SCDM models in the limit Ω b = 0 (with Ω c = 0.35 and 1.0 for the two models, respectively). They are obtained from the fitting formula given by equation (G3) of [17] and agree well with full calculations using the COSMICS code [21], which was used for the other curves in Figure 1. Baryons reduce the amount of small-scale power relative to scales k < k e because they cannot experience gravitational instability on scales smaller than the photon-baryon Jeans length until after recombination. The oscillations in the ΛCDM power spectrum are also due to baryons and are similar to the cosine integral contributions to δ c (k) in equation (51). They arise from the acoustic oscillations of the baryons between τ e and τ rec . When the power spectra are normalized at small scales (as they effectively are with σ 8 ), increasing the baryon fraction increases the amount of large-scale power.
Aside from the baryon effect, the main parameter affecting the CDM spectrum (and the only one in the BBKS formula) is k e . As may be checked using equation (50), the peak of the CDM power spectrum is closely equal to k e . Measurements of galaxy clustering actually constrain k e /h = 0.1Γ Mpc −1 where Γ ≡ Ω m h is a commonly used parameter. However, as Figure 1 shows, a second parameter, the baryon fraction Ω b /Ω m , also significantly affects the CDM spectrum.

Evolution of Large Scale Structure
The results given above concerned the shape of the CDM power spectrum, which is frozen in after the universe becomes matter-dominated. However, the amplitude grows as a(τ ) as assumed in equation (53) only if Ω m = 1. While this is a good approximation soon after recombination, it breaks down at small redshift when Ω m < 1 today. In this case we must include dark energy, curvature, and any other important contributor to the expansion rate of the universe.
The evolution of CDM density fluctuations for τ ≫ τ e is given by equations (14) and (20) ignoring the radiation terms. On subhorizon scales for φ 2 ≪ δ 2 m ≪ 1 we have Because the coefficients are spatially homogeneous, the time and space dependences of the solution factor: For Ω m > 0, D + grows with time and D − decays. For Ω m = 1 the growing solution is D + = a. For Ω m < 1 the solutions depend on the curvature K (sometimes expressed as Ω K = −K/H 2 0 ) and on the equation of state of the dark energy. See Peebles [22] and Peacock [15] for examples.
Peculiar velocities offer another probe of large-scale structure that is sensitive directly to all matter, dark or luminous. The divergence of peculiar velocity v = d x/dτ depends on the growth rate of density fluctuations: Measurements of galaxy velocities actually give v/aH, so that the ratio of velocity divergence to density fluctuation gives the logarithmic growth rate of density fluctuations. An empirical fit to a range of models gives d ln D + /d ln a = f (Ω m ) ≈ Ω 0.6 . In practice, the galaxy number density may be a biased tracer of the total matter density, conventionally parameterized by the bias factor b: δ galaxies = bδ m . Comparing galaxy clustering and velocities then allows determination of Ω 0.6 m /b.

Numerical Simulation Methods
Linear theory breaks down once dark matter clumps become significantly denser than the cosmic mean. Once this happens, full numerical simulations are needed to accurately follow the nonlinear dynamics of structure formation. Moreover, radiative processes, star formation, and supernova feedback have a strong effect on the baryonic component once galaxies begin forming. Although these processes can only be included in simulations in a phenomenological way, numerical simulations are still the best way to advance our theoretical understanding of structure formation deep into the nonlinear regime.
Here we present only a brief synopsis of numerical simulation methods. For a comprehensive review of both the methods and their applications see [23].
Particle methods are used for representing the dark matter. The phase space of dark matter is represented by a sample of particles evolving under their mutual gravitational field: (59) where δ D is the Dirac delta function. In practice, the particles are spread over a softening distance to avoid unphysical two-body relaxation due to large-angle scattering. (A simulation may use particles of mass 10 9 M ⊙ , leading to unphysical tight binaries if the forces are not softened.) The art of dark matter simulation is mainly in calculating forces much faster than O(N 2 ) for a softening distance several orders of magnitude smaller than the simulation volume. A variety of algorithms and codes are available; see [23] for details.
The initial conditions for dark matter positions and velocities are given from linear theory by the Zel'dovich approximation [7], which relates the comoving position to its initial value q, a Lagrangian variable (i.e. one that labels each particle): x( q, τ ) = q + D + (τ ) ψ( q ) , ∇ q · ψ = −δ m .
The initial density fluctuation field δ m ( x, τ ) is obtained as a sample of a gaussian random field with appropriate power spectrum. Baryon physics is added by numerical integration of the fluid equations (17) and (18) with the addition of heating and cooling, ionization, etc. There are two classes of methods in widespread use in cosmology: smoothed-particle hydrodynamics (SPH) and grid-based methods. SPH is an extension of the N-body dynamics of equations (59) to include pressure forces and thermodynamics. It has the advantage of concentrating the spatial resolution and computational effort where the baryons are densest, with the drawback of relatively large numerical viscosity and diffusivity owing to the particle discreteness. Grid-based methods use finite-difference techniques to approximate the fluid equations on a (usually) regular lattice. High accuracy methods developed for aerodynamic and other applications provide excellent resolution of shock waves, but with a resolution limited by the grid. Adaptive mesh refinement techniques are now being used to refine the resolution where needed, resulting in the ability to follow structure formation to much higher densities than previously possible [24].
We conclude this introduction to numerical simulation methods with a mosaic of gravitational potential and dark matter density from a high-resolution simulation of the mixed cold+hot dark matter model. Figure 2 shows the linear theory and nonlinear simulation results. On the scales shown, ΛCDM would be almost indistinguishable.
The gravitational potential shown in Figure 2 is a statistical fractal at the end of inflation, viz. a gaussian random field with power spectrum P A ∝ k −3 . During the radiation era, acoustic damping smooths the potential on scales smaller than k −1 e . Because the density fluctuations are related to the potential by two spatial derivatives, the corresponding linear density field is dominated by high spatial frequency components. Indeed, if not for the smoothing imposed by finite spatial resolution, the variance of density fluctuations would be logarithmically divergent. It is instructive to compare the post-recombination density with the inflationary potential. The density is almost proportional to the potential but with opposite sign. This is because the density/potential transfer function (52) varies only logarithmically on small scales.
The nonlinear evolution of the density field causes the mass to cluster strongly into wispy filaments and dense clumps; most of the volume has low density. The intermittency of structure is a sign of strong deviations from gaussian statistics.

Press-Schechter Theory
Numerical simulations show that most of the mass falls into dense self-gravitating clumps. As a crude approximation, let us suppose that the mass profile of every clump is spherically symmetric. In this case every gravitationally bound shell will eventually fall into the center of the clump, at a time that is dependent on the enclosed mass and energy of the shell. For linear growing mode initial conditions, collapse occurs when linear theory would predict that the mean density contrast interior to the shell is δ c = 1.686 [15]. Thus, in the spherical model, the distribution of clump masses is simply dependent on the statistics of the smoothed linear density field.
Press and Schechter [26] devised an elegant method for computing the mass distribution of gravitationally bound clumps using the simple spherical model. Their procedure is first to smooth the linear density field by convolution with a spherical "tophat" of radius R enclosing mean mass M (R) = (4π/3)ρ m R 3 : According to the Press-Schechter ansatz, the regions of the initial density field withδ > δ c = 1. Potential and density for a mixed cold+hot dark matter model with Ω ν = 0.2 [25]. A slice of size 50h −1 Mpc is shown. Upper left: potential at the end of inflation. Upper right: potential at the end of recombination. Lower left: density at the end of recombination. Lower right: density at redshift 0.