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

308 Nonlinear Evolution

Figure 14.3 Numerical simulations from scale-free initial conditions with spectral index n = 0. The time sequence runs from left to right and top to bottom. The development of a filament–cluster–void network with an increasing characteristic size is clearly seen.

or P3M in general depends upon the degree of clustering one wishes to probe. Strongly nonlinear clustering in dense environments probably requires the force resolution of P3M. For larger-scale structure analyses, where one does not attempt to probe the inner structure of highly condensed objects, PM is probably good enough. One should, however, recognise that the short-range forces are not computed exactly, even in P3M, so the apparent extra resolution may not necessarily be saying anything physical.

Some simulations of structure formation in models with scale-free (i.e. n = const.) initial conditions are shown in Figure 14.3. One can see that not only does one form isolated ‘blobs’ which resemble those handled by the hierarchical model, the appearance of pancakes and filaments is also generic. In the CDM

and HDM models, which are not scale free, the behaviour is rather simpler than the scale-free simulations which can be analysed with the techniques of Section 14.4 and 14.5. In the HDM model, where the initial spectrum is cut o on small scales, Zel’dovich pancakes form readily on supercluster scales, but that nonlinear processes do not create galaxy-size fluctuations rapidly enough to agree with the observations. The structure in a CDM model is much more clumpy on small scales but smoother on large scales.

14.6.3 Tree codes

An alternative procedure for enhancing the force resolution of a particle code whilst keeping the necessary demand on computational time within reasonable limits is to adopt a hierarchical subdivision procedure. The generic name given to this kind of technique is ‘tree code’. The basic idea is to treat distant clumps of particles as single massive pseudo-particles. The usual algorithm involves a mesh which is divided into cells hierarchically in such a way that every cell which contains more than one particle is divided into 23 sub-cells. If any of the resulting sub-cells contains more than one particle, that cell is subdivided again. There are some subtleties involved with communicating particle positions up and down the resulting ‘tree’, but it is basically quite straightforward to treat the distant forces using the coarsely grained distribution contained in the high level of the tree, while short-range forces use the finer grid. The greatest problem with such codes is that, although they run quite quickly in comparison with particle–mesh methods with the same resolution, they do require considerable memory resources. Their use in cosmological contexts has so far therefore been quite limited, one of the problems being the di culty of implementing periodic boundary conditions in such algorithms.

14.6.4 Initial conditions and boundary e ects

To complete this section, we make a few brief remarks about starting conditions for N-body simulations, and the e ect of boundaries and resolution on the final results.

Firstly, one needs to be able to set up the initial conditions for a numerical simulation in a manner appropriate to the cosmological scenario under consideration. For most models this means making a random-phase realisation of the power spectrum – see Section 14.8. This is usually achieved by setting up particles initially exactly on the grid positions, then using the Zel’dovich approximation, Equation (14.2.1), to move them such as to create a density field with the required spectrum and statistics. The initial velocity field is likewise obtained from the primordial gravitational potential. One should beware, however, the e ects of the poor k-space resolution at long wavelengths. The assignment of k-space amplitudes requires a random amplitude for each wave vector contained in the reciprocal-space version of the initial grid. As the wave number decreases,

310 Nonlinear Evolution

the discrete nature of the grid becomes apparent. For example, there are only three (orthogonal) wave vectors associated with the fundamental mode of the box. When amplitudes are assigned via some random-number generator, one must take care that the statistically poor sampling of k-space does not lead to spurious features in the initial conditions. One should use a simulation box which is rather larger than the maximum scale at which there is significant power in the initial spectrum.

At the other extreme, there arises the question of the finite spacing of the grid. This puts an upper limit, known as the Nyquist frequency, on the wavenumbers k which can be resolved, which is defined by kN = 2π/d, where d is the mesh spacing. Clearly, one should not trust structure on scales smaller than k−N1.

One is therefore warned that, although numerical methods such as these are the standard way to follow the later nonlinear phases of gravitational evolution, they are not themselves ‘exact’ solutions of the equations of motion and results obtained from them can be misleading if one does not choose the resolution appropriately.

14.7 Gas Physics

So far we have dealt exclusively with the behaviour of matter under its self-gravity. We have ignored pressure gradient terms in the equation of motion of the matter at all times after recombination. While this is probably a good approximation in the linear and quasilinear regimes, when the Jeans mass is much smaller than scales of cosmological interest, it is probably a very poor representation of the late nonlinear phase of structure formation. As we shall see, hydrodynamical e ects are clearly important in determining the behaviour of the baryonic part of galaxies, even if the baryons are only a small fraction of the total mass. Nonlinear hydrodynamical e ects connected with the formation of shocks are also very important in determining how a collapsing structure reaches virial equilibrium.

14.7.1 Cooling

One of the important things to explain in hierarchical clustering scenarios is the existence of a characteristic scale of 1011M in the mass spectrum of galaxies. Because gravity itself does not pick out any scale, some other physical mechanism must be responsible. Since only the baryonic part of the galaxy can be seen, and it is only this part which is known to possess characteristic properties, it is natural to think that gas processes might be involved. A good candidate for such a process is the cooling of the gas forming the galaxy.

Following Rees and Ostriker (1977), let us consider a simple model of a galaxy as a spherical gas cloud (i.e. no non-baryonic material) in the manner of Section 14.1. After collapse and violent relaxation (the process which converts the radial collapse motion into random ‘thermal’ motions) this cloud will be supported in equilibrium at its virial radius R and will have a temperature T GMµ/R, where µ is the mean molecular weight. If this temperature is high, as it will be for interesting mass scales, the cloud will be radiating and therefore cooling. The balance

between pressure support and gravity which determines the size of the object depends on two characteristic timescales: the cooling time

and the dynamical time, defined to be the free-fall collapse time for a sphere of mass M and radius R,

where ρ is the mean baryon density and Λ(T) in Equation (14.7.1) is the cooling rate (energy loss rate per unit volume per unit time) for a gas at temperature T (Λ is tabulated in standard physics texts for di erent kinds of gas). There are three main contributions to cooling in a hydrogen–helium plasma which is what we expect to have in the case of galaxy formation: free–free (bremsstrahlung) radiation, recombination radiation from H and He, and Compton cooling via the cosmic microwave background. This last one is e cient only if z > 10 or so. Since it is not known whether galaxy formation might have taken place at such high redshifts, this may play a role but for simplicity we shall ignore it here.

The two timescales tdyn and tcool, together with the expansion timescale τH = H−1, determine how the protogalaxy cools as it collapses. If tcool > τH, then cooling cannot have been important and the cloud will have scarcely evolved since its formation. If τH > tcool > tdyn, then the gas can cool on a cosmological timescale, but the fact that it does so more slowly than the dynamical characteristic time means that the cloud can adjust its pressure distribution to maintain the support of the cooling matter. There is thus a relatively quiescent quasi-static collapse on a timescale tcool. The last possibility is that tcool < tdyn. Now the cloud cools so quickly that dynamical processes are unable to adjust the pressure distribution in time: pressure support will be lost and the gas undergoes a rapid collapse on the free-fall timescale, accompanied by fragmentation on smaller and smaller scales as instabilities develop in the cloud which is behaving isothermally.

It is thought that the condition tcool < tdyn is what determines the characteristic mass scale for galaxies. Only when this criterion is satisfied can the gas cloud collapse by a large factor and fragment into stars which allow the cloud to be identified as a galaxy. Furthermore, if structure formation proceeds hierarchically, the gas must cool on a timescale at least as small as tdyn, otherwise it will not be confined in a bound structure on some particular scale but will instead be disrupted as the next level of the hierarchy forms.

Let us now add non-baryonic matter into this discussion. What changes here is that the dynamical timescale for a collapsing cloud will be dominated by the dark matter while cooling is enjoyed only by the gas. Let us assume a spherical collapse model again. Notice that the dynamical timescale (14.7.2) is essentially the time taken for a perturbation to collapse from its maximum extent which can be identified as the turnaround radius Rm in Section 15.1. Putting in some

312 Nonlinear Evolution

numbers one finds that

One can estimate the cooling timescale by assuming that gas makes up a fraction Xb of the total mass M and that it is uniformly distributed within the virial radius which will be Rm/2. We then take the gas temperature to be the same as the virial temperature of the collapsed object: T 2GMµ/5kBRm. We also assume that the gas has not been contaminated by metals from an early phase of star formation (metals can increase the cooling rate and thus lower the cooling time considerably), and therefore adopt the appropriate value of Λ(T) for a pure hydrogen plasma at temperature T. Using Equation (14.7.1) we find that

which, for Xb 0.05, gives M 3 × 1011M . While this theory therefore gives a plausible account of the characteristic mass scale for galaxies, it is obviously extremely simplified. Hydrodynamical e ects may be important in many other contexts, such as cluster formation, the collapse of pancakes and also the feedback of energy from star formation into the intergalactic medium. A detailed theory of the origin of structure including gas dynamics, dissipation and star formation is, however, still a long way from being realised.

14.7.2 Numerical hydrodynamics

In the above we discussed an example where gas pressure forces are important in the formation of cosmic structure. Understanding of these e ects is highly qualitative and applicable only to simple models. In an ideal world, one would like to understand the influence of gas pressure and star formation in a general context. E ectively, this means solving the Euler equation, including the relevant pressure terms, self-consistently. The appropriate equation is

The field of cosmological hydrodynamics is very much in its infancy, and it is fair to say that there are no analytic approximations that can be implemented with any confidence in this kind of analysis. The only realistic hope for progress in the near future lies with numerical methods, so we describe some of the popular techniques here.

In smoothed-particle hydrodynamics (SPH) one typically represents the fluid as a set of particles in the same way as in the N-body gravitational simulations described in Section 14.6. Densities and gas forces at particle locations are thus calculated by summing pairwise forces between particles. Since pressure forces are expected to fall o rapidly with separation, above some smoothing scale h (see below), it is reasonable to insert the gas dynamics into the part of a particle code that details the short-range forces such as the particle–particle part of a P3M code. It is, however, possible to include SPH dynamics also in other types of simulation, including tree codes.

One technique used to insert SPH dynamics into a P3M code is to determine local densities and pressure gradients by a process known as kernel estimation. This is essentially equivalent to convolving a field f(x) with a filter function W to produce a smoothed version of the field:

fs(r) = f(x)W(x − r) d3x, (14.7.7)

where W contains some implicit smoothing scale; one possible choice of W is a Gaussian. If f(x) is just the density field arising from the discrete distribution of particles, then it can be represented simply as the sum of delta-function contributions at each particle location xi and one recovers Equation (14.6.2). We need to represent the pressure forces in the Euler equation: this is done by specifying the equation of state of the fluid p = (γ − 1)Hρ, where H is the thermal energy, ρ the local density and p the pressure. Now one can write the pressure force term in Equation (14.7.6) as

The gradient of the smoothed function fs can be written

fs(r) = f(x) W(x − r) d3x, (14.7.9)

so that the gas forces can be obtained in the form

The form of Equation (14.7.10) guarantees conservation of linear and angular momentum when a spherically symmetric kernel W is used. The adiabatic change in the internal energy of the gas can similarly be calculated:

where vij is the relative velocity between particles. For collisions at a high Mach number, defined as the ratio of any systematic velocity to the thermal random

314 Nonlinear Evolution

velocity, thermal pressure will not prevent the particles from streaming freely, but in real gases there is molecular viscosity which prevents interpenetration of gas clouds. This is modelled in the simulations by introducing a numerical viscosity, the optimal form of which depends upon the nature of the simulation being attempted.

The advantage of particle-based methods is that they are Lagrangian and consequently follow the motion of the fluid. In practical terms, this means that most of the computing e ort is directed towards places where most of the particles are and, therefore, where most resolution is required. As mentioned above, particle methods are the standard numerical tool for cosmological simulations. Classical fluid dynamics, on the other hand, has usually followed an Eulerian approach where one uses a fixed (or perhaps adaptive) mesh. Codes have been developed which conserve flux and which integrate the Eulerian equations of motion rapidly and accurately using various finite-di erence approximation schemes. It has even proved possible to introduce methods for tracking the behaviour of shocks accurately – something which particle codes struggle to achieve. Typically, these codes can treat many more cells than an SPH code can treat particles, but the resolution is usually not so good in some regions because the cells will usually be equally spaced rather than being concentrated in the interesting high-density regions.

An extensive comparison between Eulerian and Lagrangian hydrodynamical methods has recently been performed, which we recommend to anyone thinking of applying these techniques in a cosmological context. Each has its advantages and disadvantages. For example, density resolution is better in the state-of-the- art Lagrangian codes, and the thermal accuracy better in the Eulerian codes. Conversely, Lagrangian methods have poor accuracy in low-density regions, presumably due to statistical e ects, while the Eulerian codes usually fail to resolve the temperatures correctly in high-density regions due to the artificially high numerical viscosity in them.

14.8 Biased Galaxy Formation

It should be obvious by now that the complexities of nonlinear gravitational evolution, together with the possible influence of gas-dynamical processes on galaxy formation, mean that a full theory of the formation of these objects is by no means fully developed. Structure on larger scales is less strongly nonlinear, and therefore is less prone to hydrodynamical e ects, so may be treated fairly accurately using linear theory as long as σM 1 or, better still, using approximation methods such as the Zel’dovich and adhesion approximations. The problem is that, when one seeks observational data with which to compare theoretical predictions, these data invariably involve the identification of galaxies. Even if we give up on the task of understanding the details of the galaxy-formation process, we still need to know how to relate observations of the large-scale distribution of galaxies to that of the mass.

In Section 13.9 we discussed the Poisson clustering model, which is a statistical statement of the form ‘galaxies trace the mass’. In this model the two-point cor-

relation function of galaxies is equal to the covariance function of the underlying density field. In recent years, however, it has become clear that this is probably not a good representation of reality. In the spirit of the spherical collapse model one might imagine that galaxies should form not randomly sprinkled around according to the local density of matter, but at specific locations where collapse, cooling and star formation can occur. Obvious sites for protostructures would therefore be peaks of the density field, rather than randomly chosen sites. This simple idea, together with the assumption that the large-scale cosmological density field is Gaussian (see Section 14.8), led Kaiser (1984) (in a slightly di erent context; see Section 16.5) to suggest a biased galaxy formation, so that the galaxy correlation function and the matter autocovariance function are no longer equivalent. The way such a bias might come about is as follows. Suppose the density field δM, smoothed on some appropriate mass scale M to define a galaxy, is Gaussian and has variance σM2 . The covariance function ξ(r) of δM is

where the average is taken over all spatial positions x and x such that |x−x | = r. If galaxies trace the mass, then the two-point correlation function of galaxies ξgg(r) coincides with ξ(r). If galaxies do not trace the mass, this equality need not hold. In particular, imagine a scenario where galaxies only form from highdensity regions above some threshold δc = νσM , where ν is a dimensionless threshold. The existence of such a threshold is qualitatively motivated by the spherical model of collapse, described in Section 14.1, within which a linear value of δc 1.68 would seem to be required for structure formation. To proceed we need to recall that, for such a Gaussian field, all the statistical information required to specify its properties is contained in the autocovariance function ξ(r). It is straightforward to calculate the correlation function of points exceeding δc using the Gaussian prescription because the probability of finding two regions separated by a distance r both above the threshold will be just

Now, as explained in Section 13.7, the N-variate joint distribution of a set of δi can be written as a multivariate Gaussian distribution: for the case where N = 2, which is needed in Equation (14.8.2), using the substitution δi = νiσ and w(r) =

ξ(r)/σ2, we find

The two-point correlation function for points exceeding νc = δc/σ is then

316 Nonlinear Evolution

where Q1 = P>δc ; see Equation (14.5.4). The exact calculation of the integrals in this equation is di cult but various approximate relations have been obtained. For large νc and small w we have

while another expression, valid when w is not necessarily small, is

Kaiser initially introduced this model to explain the enhanced correlations of Abell clusters compared with those of galaxies; see Section 16.5. Here the field δ is initially smoothed with a filter of radius several Mpc to pick out structure on the appropriate scale. If galaxies trace the mass, and so have ξgg(r) ξ(r), then the simple relation (14.8.5) explains qualitatively why cluster correlations might have the same slope, but a higher amplitude than the galaxy correlations. This enhancement is natural because rich clusters are defined as structures within which the density of matter exceeds the average density by some fairly well-defined factor in very much the way assumed in this calculation.

This simple argument spawned more detailed analyses of the statistics of Gaussian random fields, culminating in the famous ‘BBKS’ paper of Bardeen et al. (1986), which have refined and extended, while qualitatively confirming, the above calculations. The interest in most of these studies was the idea that galaxies themselves might form only at peaks of the linear density field (this time smoothed with a smaller filtering radius). If galaxies only form from large upward fluctuations in the linear density field, then they too should display enhanced correlations with respect to the matter. This seemed to be the kind of bias required to reconcile the standard CDM model with observations of galaxy-peculiar motions and also the cause of the apparent discrepancy between dynamical estimates of the mass density of the Universe of around Ω0 0.2 when the theoretically favoured value is Ω0 1. We shall discuss the question of velocities in detail in Chapter 18 and we have referred to it also in Chapter 4. Nevertheless, some comments here are appropriate. The velocity argument can be stated simply in terms of a sort of cosmic virial theorem. If galaxies trace the mass, and have correlation function ξ(r) and mean pairwise velocity dispersion at a separation r equal to v2(r), then this theorem states that

with a calculable constant of proportionality; see Section 18.5 for details.

There are problems with this theorem in the context of standard CDM. First, if one runs a numerical simulation of CDM to the point when the correlation function of the mass has the right slope compared with that of the observations, then the accompanying velocities v are far too high. A low-density CDM seems to be a much better bet in this respect, but this may be because the slope of the correlation function is not a very good way to determine the present epoch in a simulation. The same thing, however, seems to happen in our Universe, where the

observed correlation function and the observed pairwise peculiar motions give Ω 0.2. One way out of this, indeed the obvious way out apart from the fact that it appears to contradict inflation, is to have Ω 0.2 and leave it at that. There is another way out, however, which involves bias of the sort discussed above. Taking (14.8.6) as a qualitative model, one might argue that in fact ξ(r) is wrong by a factor ν2/σM2 and, if this bias is large, one can reconcile a given v with Ω = 1. A bias factor b, defined by

of around b 1.5–3 seems to be required to match small-scale clustering and peculiar velocity data with the standard CDM model. Notice also that true density fluctuations are smaller than the apparent fluctuations in counts of galaxies, so that fluctuations in the microwave background are smaller by a factor 1/b in this picture than they would be if galaxies trace the mass.

The parameter b often arises in the cosmological literature to represent the possible di erence between mass statistics and the statistics of galaxy clustering. The usual definition is not (14.8.8) but rather

b2 = σ2(galaxies), (14.8.9)

8

2

σ8 (mass)

where σ82 represents the dimensionless variance in either galaxy counts or mass in spheres of radius 8h−1 Mpc. This choice is motivated by the observational result that the variance of counts of galaxies in spheres of this size is of order unity, so that b 1/σ8(mass). Unless stated otherwise, this is what we shall mean by b in the rest of this book. Many authors use di erent definitions, e.g.

which is called the linear bias model. While a relation of the form (14.8.10) clearly entails (14.8.9) and (14.8.8), it does not follow from them, so these definitions are not equivalent. While there is little motivation, other than simplicity, for supposing the bias parameter to be a simple constant multiplier on small scales, it can be shown that, as long as the bias acts as a local function of the density, the form (14.8.8) should hold on large scales, even if the biasing relationship is complicated (Coles 1993).

Alternatives to (14.8.10), which are not equivalent, include the high-peak model and the various local-bias models (Coles 1993). Non-local biases are possible, but it is rather harder to construct such models (Bower et al. 1993). If one is prepared to accept an ansatz of the form (14.8.10), then one can use linear theory on large scales to relate galaxy-clustering statistics to those of the density fluctuations, e.g.