Cosmology. The Origin and Evolution of Cosmic Structure - Coles P., Lucchin F

288 Nonlinear Evolution

(as usual, the symbol ‘+’ indicates the growing mode, while ‘−’ denotes the decaying mode). The combination of growing and decreasing modes in Equations (14.1.1) is necessary to satisfy the correct boundary condition on the velocity: Vi = 0 requires that δ+(ti) = 35 δi. One can assume that, after a short time, the decaying mode will become negligible and the perturbation remaining will just be δ δ+(ti). Let us take the initial value of the Hubble expansion parameter to be Hi. Assuming that pressure gradients are negligible, the sphere representing the perturbation evolves like a Friedmann model whose initial density parameter is given by

where the su x ‘p’ denotes the quantity relevant for the perturbation, while ρ(ti) and Ω(ti) refer to the unperturbed background universe within which the perturbation resides. Structure will be formed if, at some time tm, the spherical region ceases to expand with the background universe and instead begins to collapse. This will happen to any perturbation with Ωp(ti) > 1. From Equations (14.1.2) and (2.6.4) this condition can easily be seen to be equivalent to

where Ω is the present value of the density parameter. In universes with Ω < 1, however, the fluctuation must exceed the critical value (1 − Ω)/Ω(1 + zi); it is interesting to note that in this case the condition (14.1.3) implies that the growing perturbation reaches the nonlinear regime before the time t at which the universe becomes curvature dominated and therefore enters a phase of undecelerated free expansion. For Ω 1, on the other hand, there is no problem.

from which we easily obtain that the density of the perturbation at time tm is

the value of tm, from Equation (2.4.9) (where t0 is replaced by ti) and Equation (14.1.5), is just

In an Einstein–de Sitter universe the ratio χ between the background density, ρ(tm), and the density inside the perturbation, ρp(tm), is obtained from the previous equation and from

which corresponds to a perturbation δ+(tm) 4.6; the extrapolation of the linear growth law, δ+ t2/3, would have yielded, from (14.1.6),

corresponding to the approximate value ρp(tm)/ρ(tm) 1 + δ+(tm) 2.07. The perturbation will subsequently collapse and, if one can still ignore pressure e ects and the configuration remains spherically symmetric, in a time tc of order 2tm, one will find an infinite density at the centre. In fact, when the density is high, slight departures from this symmetry will result in the formation of shocks and considerable pressure gradients. Heating of the material will occur due to the dissipation of shocks which converts some of the kinetic energy of the collapse into heat, i.e. random thermal motions. The end result will therefore be a final equilibrium state which is not a singular point but some extended configuration with radius Rvir and mass M. From the virial theorem the total energy of the fluctuation is

If in the collapsing phase we can ignore the possible loss of mass from the system due to e ects connected with shocks, and possible loss of energy by thermal radiation, the energy and mass in (14.1.10) are the same as the fluctuation had at time tm,

where Rm is the radius of the sphere at the moment of maximum expansion. Having assumed that the pressure is zero, in Equation (14.1.11) no account is taken of the contribution of thermal energy; the kinetic energy due to the expansion is zero by definition at this point. From Equations (14.1.10) and (14.1.11) we therefore have Rm = 2Rvir, so that the density in the equilibrium state is ρp(tvir) = 8ρp(tm). One usually assumes that at tc, the time of maximum compression, the density is of order ρp(tvir). Numerical simulations of the collapse allow an estimate to be made of the time taken to reach equilibrium: one finds that tvir 3tm. If at times tc and tvir the universe is still described by an Einstein–de Sitter model, the ratios

290 Nonlinear Evolution

between the density in the perturbation and the mean density of the universe at these times are

respectively. An extrapolation of linear perturbation theory would give

which correspond to values of 2.68 and 3.20 for the ratio of the densities, in place of the exact values given by Equations (14.1.12 a) and (14.1.12 b).

14.2 The Zel’dovich Approximation

The model discussed in the previous section, though very instructive in its conclusions, su ers from some notable defects. Above all, reasonable models of structure formation do not contain primordial fluctuations at ti trec, which are organised into neat homogeneous spherical regions with zero peculiar velocity at their edge. Moreover, even if this were the case at the beginning, such a symmetrical configuration is strongly unstable with respect to the growth of non-radial motions during the expansion and collapse phases of the inhomogeneity. In fact, the classic work of Lin et al. (1965) showed that, for a generic triaxial perturbation, the collapse is expected to occur not to a point, but to a flattened structure of quasi-two-dimensional nature. The usual descriptive term for such features is pancakes.

The spherical top-hat model is only reasonably realistic for perturbations on scales just a little larger than MJ(i)(zrec). In this case, however, pressure is not negligible and dissipation can be significant during the collapse. Presumably what form in such a situation are more or less spherical protoobjects in which gravity is balanced by pressure forces.

It is more complicated to study the development of perturbations on scales M MD(a)(zrec). Of course, one could simply resort to numerical methods like those we shall discuss in Section 15.5. However, some simplifying assumptions are possible. For example, in this situation, pressure would be e ectively zero and the fluid can be treated like dust. Under this assumption it is in fact possible to understand the growth of structure analytically using a clever approximation devised by Zel’dovich (1970). This approximation actually predicts that the density in certain regions – called caustics – should become infinite, but the gravitational acceleration caused by these regions remains finite. Of course, in any case one cannot justify ignoring pressure when the density becomes very high, for much the same reason as we discussed in Section 15.1 in the context of spherical

collapse: one forms shock waves which compress infalling material. At a certain point the process of accretion onto the caustic will stop: the condensed matter is contained by gravity within the final structure, while the matter which has not passed through the shock wave is held up by pressure. It has been calculated that about half the material inside the original fluctuation is reheated and compressed by the shock wave. An important property of the structures which thus form is that they are strongly unstable to fragmentation. In principle, therefore, one can generate structure on smaller scales than the pancake.

Let us now describe the Zel’dovich approximation in more detail, and show how it can follow the evolution of perturbations until the formation of pancakes. Imagine that we begin with a set of particles which are uniformly distributed in space. Let the initial (i.e. Lagrangian) coordinate of a particle in this unperturbed distribution be q. Now each particle is subjected to a displacement corresponding to a density perturbation. In the Zel’dovich approximation the Eulerian coordinate of the particle at time t is

where r = a(t)x, with x a comoving coordinate, and we have made a(t) dimensionless by dividing throughout by a(ti), where ti is some reference time which we take to be the initial time. The derivative on the right-hand side is taken with respect to the Lagrangian coordinates. The dimensionless function b(t) describes the evolution of a perturbation in the linear regime, with the condition b(ti) = 0, and therefore solves the equation

This equation corresponds to (10.6.14), with vanishing pressure term, which describes the gravitational instability of a matter-dominated universe. For a flat matter-dominated universe we have b t2/3 as usual. The quantity Φ0(q) is proportional to a velocity potential, i.e. a quantity of which the velocity field is the gradient, because, from Equation (14.2.1),

this means that the velocity field is irrotational. The quantity Φ0(q) is related to the density perturbation in the linear regime by the relation

which is a simple consequence of Poisson’s equation.

The Zel’dovich approximation is therefore simply a linear approximation with respect to the particle displacements rather than the density, as was the linear solution we derived above. It is conventional to describe the Zel’dovich approximation as a first-order Lagrangian perturbation theory, while what we have dealt

292 Nonlinear Evolution

with so far for δ(t) is a first-order Eulerian theory. It is also clear that Equation (14.2.1) involves the assumption that the position and time dependence of the displacement between initial and final positions can be separated. Notice that particles in the Zel’dovich approximation execute a kind of inertial motion on straight line trajectories.

The Zel’dovich approximation, though simple, has a number of interesting properties. First, it is exact for the case of one-dimensional perturbations up to the moment of shell crossing. As we have mentioned above, it also incorporates irrotational motion, which is required to be the case if it is generated only by the action of gravity (due to the Kelvin circulation theorem). For small displacements between r and a(t)q, one recovers the usual (Eulerian) linear regime: in fact, Equation (14.2.1) defines a unique mapping between the coordinates q and r (as long as trajectories do not cross); this means that ρ(r, t) d3r = ρ(ti) d3q or

where |J(r, t)| is the determinant of the Jacobian of the mapping between q and r: ∂r/∂q. Since the flow is irrotational, the matrix J is symmetric and can therefore be locally diagonalised. Hence

the quantities 1 +b(t)αi are the eigenvalues of the matrix J (the αi are the eigenvalues of the deformation tensor). For times close to ti, when |b(t)αi| 1, Equation (14.2.6) yields

which is the law of perturbation growth in the linear regime.

Equation (14.2.6) indicates that at some time tsc, when b(tsc) = −1/αj, an event called shell-crossing occurs such that a singularity appears and the density becomes formally infinite in a region where at least one of the eigenvalues (in this case αj) is negative. This condition corresponds to the situation where two points with di erent Lagrangian coordinates end up at the same Eulerian coordinate. In other words, particle trajectories have crossed and the mapping (14.2.1) is no longer unique. A region where the shell-crossing occurs is called a caustic. For a fluid element to be collapsing, at least one of the αj must be negative. If more than one is negative, then collapse will occur first along the axis corresponding to the most negative eigenvalue. If there is no special symmetry, one therefore expects collapse to be generically one dimensional, i.e. to a sheet or ‘pancake’. Only if two (or three) negative eigenvalues, very improbably, are equal in magnitude can the collapse occur to a filament (or point). One therefore expects ‘pancake’ formation to be the generic result of structure collapse.

The Zel’dovich approximation matches very well the evolution of density perturbations in full N-body calculations until the point where shell crossing occurs

Figure 14.1 Comparison of the Zel’dovich approximation (b) and an N-body experiment

(a) for the same initial conditions. Agreement is good, except for the ‘fuzzy’ appearance of the pancake regions which is due to the motion of particles after shell-crossing.

(Coles et al. 1993a); we shall discuss N-body methods later on. After this, the approximation breaks down completely. According to Equation (14.2.1) particles continue to move through the caustic in the same direction as they did before. Particles entering a pancake from either side merely sail through it and pass out the opposite side. The pancake therefore appears only instantaneously and is rapidly smeared out. In reality, the matter in the caustic would feel the strong gravity there and be pulled back towards it before it could escape through the other side. Since the Zel’dovich approximation is only kinematic it does not account for these close-range forces and the behaviour in the strongly nonlinear regime is therefore described very poorly. Furthermore, this approximation cannot describe the formation of shocks and phenomena associated with pressure. The problem of shell-crossing is inevitable in the Zel’dovich approximation. In order to prevent this from interfering too much in calculations, one can filter out the small-scale fluctuations from the initial conditions which give rise to shell-crossing. If the power spectrum is a decreasing function of mass, then the large scales can be evolving in the quasilinear regime (i.e. before shell-crossing) even when a higher resolution would reveal considerable small-scale caustics. By smoothing the density field one removes these small-scale events but does not alter the kinematical evolution of the large-scale field. The best way to implement this idea appears to be to filter the initial power spectrum according to

where knl < kG < 1.5knl and knl is the characteristic nonlinear wavenumber given approximately by

294 Nonlinear Evolution

so that the RMS density fluctuation σM on a scale R 2π/knl is of order unity. The performance of the Zel’dovich approximation, the ‘smoothed’ Zel’dovich approximation and a full N-body simulation from a realisation of Gaussian initial conditions is shown in Figure 14.1.

14.3 The Adhesion Model

The smoothed Zel’dovich approximation merely ignores the problem of shellcrossing. If one is forced to deal with it, in other words if one wants to study the mass distribution on scales where σM > 1, then one must come up with some other approach. One relatively straightforward way to extend the Zel’dovich approximation is through the so-called adhesion model.

In the adhesion model one assumes that the particles stick to each other when they enter a caustic region because of an artificial viscosity which is intended to simulate the action of strong gravitational e ects inside the overdensity forming there. This ‘sticking’ results in a cancellation of the component of the velocity of the particle perpendicular to the caustic. If the caustic is two dimensional, the particles will move in its plane until they reach a one-dimensional interface between two such planes. This would then form a filament. Motion perpendicular to the filament would be cancelled, and the particles will flow along it until a point where two or more filaments intersect, thus forming a node. The smaller the viscosity term is, the thinner the sheets and filaments will be, and the more point-like the nodes will be. Outside these structures, the Zel’dovich approximation is still valid to high accuracy. Comparing simulations made within this approximation with full N-body calculations shows that it is quite accurate for overdensities up to δ 10.

Let us begin by rewriting the Euler and continuity equations, together with the Poisson equation (all ignoring the e ects of pressure), in a slightly altered form

and r = a(t)x; x is a comoving coordinate. The Equation (14.3.1 c) is not needed in this section, but we have included it here for the sake of completeness. The Zel’dovich approximation is equivalent to putting the right-hand side of (14.3.1 a)

∂U

∂b

The adhesion model involves modifying the Equation (14.3.2 b) by introducing a viscosity term ν, which allows the particles to stick together:

The e ect of this term is to make the particles ‘feel’ the inside of collapsed structures. It remains negligible outside these regions. The viscosity ν has the dimensions of ‘length squared’ in this representation because our ‘time’ coordinate is actually dimensionless, so the model basically requires that d √ν should be much less than the typical dimension of the structures forming. Equation (14.3.3) is well known in the mathematical literature as the Burgers equation. In many cases, and this is true in our case, this equation has an exact solution. With the so-called Hopf–Cole substitution,

For small values of ν the main contribution to the integral in Equation (14.3.6) comes from regions where the function G has a maximum. This property allows a simplified treatment of the problem. The Eulerian position of the particle can be found by solving the integral equation

The adhesion model furnishes results in accord with the Zel’dovich approximation at distances l d from the structure, but allows one to follow the formation of structure insofar as it prevents structure from being erased by shell-crossing. It also allows one to avoid the singularities which occur in the usual Zel’dovich approximation. In many simple cases the solution (14.3.6) does indeed allow one to study the formation of structure to high accuracy even in a highly advanced phase of nonlinearity.

The spatial distribution of particles obtained by letting the parameter ν tend to zero represents a sort of ‘skeleton’ of the real structure: nonlinear evolution generically leads to the formation of a quasicellular structure, which is similar

296 Nonlinear Evolution

to a ‘tessellation’ of irregular polyhedra having pancakes for faces, filaments for edges and nodes at the vertices. This skeleton, however, evolves continuously as structures merge and disrupt each other through tidal forces; gradually, as evolution proceeds, the characteristic scale of the structures increases. In order to interpret the observations we have already described in Chapter 4, one can think of the giant ‘voids’ as being the regions internal to the cells, while the cell nodes correspond to giant clusters of galaxies. While analytical methods, such as the adhesion model, are useful for mapping out the skeleton of structure formed during the nonlinear phase, they are not adequate for describing the highly nonlinear evolution within the densest clusters and superclusters. In particular, the adhesion model cannot be used to treat the process of merging and fragmentation of pancakes and filaments due to their own (local) gravitational instabilities.

14.4 Self-similar Evolution

A possible way to treat highly nonlinear evolution in the framework of ‘bottom-up’ scenarios is to introduce the concept of self-similarity or hierarchical clustering. As we have already explained, in the isothermal baryon model or in the more modern CDM model, the first structures to enter the nonlinear regime are expected to be on a mass scale of order MJ(i)(zrec). Galaxies and larger structures then form by merging of such objects into objects of higher mass. This process is qualitatively di erent from that described by the Zel’dovich and adhesion approximations, which are more likely to be accurate on scales relevant to clusters and superclusters, while we need something else to describe the formation of structure on scales up to this.

14.4.1 A simple model

To illustrate some of these ideas, let us assume that the Universe is well-described by an Einstein–de Sitter model. A perturbation with mass M > MJ, which we use from now on to mean MJ(i)(zrec), arrives in the nonlinear regime, approximately, at a time tM such that

where σM (trec) is the RMS mass fluctuation on the scale M at t = trec. One therefore has the relationship

where the quantity αrec is defined in Section 13.4. From Equation (14.4.2) it follows that

where tJ = tM for M MJ. As we explained in Section 14.1, if we think of the perturbation as a spherical ‘blob’, then the time tM practically coincides with the moment at which the perturbation ceases to expand with the background Universe and begins to collapse. In the general case expressed by (14.4.2), one can apply the simple scheme described in Section 14.1: one can easily obtain from Equations (14.1.6) and (14.4.2) that, at virial equilibrium, the perturbation has a density

where we have put ρM(MJ) = ρJ. If rM is the radius of a (collapsed) perturbation of mass M, from (14.4.4) and from the fact that M ρMrM3 , one finds

where the meaning of rJ is clear; the exponent γvir is given by the relation

We can also relate the mass M to the virial velocities generated by it, VM , in this model. The result is

If nrec = −2, then this can explain the M VM4 relationship implied by the observed correlation between L and V for galaxies, known as the Tully–Fisher relationship, Equation (4.3.2).

A simple interpretation of the model just described, which is called the hierarchical clustering model, is the following. The Universe at time tM on a scale r < rM contains condensed objects of various masses M and corresponding sizes rM according to a hierarchical arrangement, in which the objects of one scale are the building blocks from which objects on higher scales are made. This arrangement holds up to the scale M which is the largest mass scale to have reached virial equilibrium. For masses greater than M , fluctuations are small and still evolving in the linear regime so that, for r > rM , we have δρm(r) σM M−αrec r−3αrec = r−(3+nrec)/2. These small fluctuations will grow and, when t > tM , objects on a higher mass scale than M will collapse and form a higher level of the hierarchy. Simple though it is, this description seems to provide a fairly accurate representation of the behaviour of N-body simulations of hierarchical clustering in the highly nonlinear phase.