text:
# Simulating the Formation of the Local Galaxy Population

## 1 Introduction

Over the last decade phenomenological modelling has made it possible to follow many aspects of the formation and evolution of galaxies within the currently favored hierarchical paradigm for the growth of cosmic structure. Recently, the grafting of techniques originally developed by White & Frenk (1991); Kauffmann et al. (1993); Cole et al. (1994) onto high resolution N-body simulations has allowed the spatial and kinematic distributions of galaxies to be predicted in detail as a function of their intrinsic properties (Kauffmann et al. 1997, 1999a, 1999b; Diaferio et al. 1999, 2001; Benson et al. 2000a, b, 2001a, 2001b; Springel et al. 2000). This work clarifies many aspects of the problem of ‘galaxy biasing’ and supersedes the heuristic models previously used to relate the galaxy and mass distributions in CDM cosmogonies. Its goal is twofold: to better understand the physical processes driving the formation, evolution and clustering of galaxies, and to test the standard structure formation paradigm.

With the steady improvement of computer performance and of simulation codes, dissipationless simulations are able to achieve ever higher mass resolution. In a recent example, Springel et al. (2000, hereafter S00) followed the formation of all galaxies brighter than the Fornax dwarf spheroidal within a cluster similar in mass to Coma. In the simulations presented below, the use of $`7\times 10^7`$ particles allows us to follow the formation of all galaxies more luminous than the LMC out to 8000 km s-1 from the Milky Way. Improved computers have also greatly enhanced the ability of cosmic gas-dynamics codes to simulate galaxy formation (e.g. Pearce et al. 2001; Nagamine et al. 2000; White et al. 2000). Although such simulations remove the uncertainties due to phenomenological modelling of the dynamics and cooling of diffuse gas, they retain phenomenological models for the much more uncertain processes of star formation and supernova feedback. Moreover, their greatly increased computational cost makes it impossible for them to resolve galaxy formation over a volume as large as that modelled in this paper. A further major advantage of the techniques used here is that the efficiencies for uncertain processes like star formation and feedback can be varied to study their influence and then adjusted to fit observation without the need to run a new simulation for each new parameter set.

The simulations we present below use initial conditions generated using a technique developed by Kolatt et al. (1996, see () for a related technique). An all-sky redshift survey is smoothed heavily to produce an estimate of the galaxy overdensity field in a spherical volume centred on the Milky Way. This is assumed to be a known constant times the similarly smoothed local mass overdensity field. One then solves for the linear overdensity field at high redshift which would evolve into this quasi-linear local field. The initial conditions for the simulation are taken to be a random realisation of a Gaussian random field with a suitable CDM (or other) power spectrum, but constrained so that when suitably smoothed the overdensity field is equal to that inferred from the local galaxy distribution. Simulations run from such initial conditions reproduce the large-scale structure of the local universe but have the characteristics of the assumed CDM model on small scales where the smoothed local galaxy distribution imposes no constraints.

Although galaxy formation should occur in such simulations exactly as in random realisations of the underlying CDM model, there are a number of advantages to having the large-scale structure of the simulation correspond in detail to that of the local universe. For certain types of galaxies, in particular dwarfs, surveys are restricted to our local neighborhood. For many types of galaxies, surveys are most complete and the properties of the galaxies best characterised in this region. When interpreting such surveys one must be wary of biases introduced by the particular structure of our neighborhood. Such biases are clearly minimised in models which reproduce the local structure. Distances to galaxies can only be measured with sufficient accuracy to estimate their peculiar motions out to redshifts of about 10,000 km s-1. As a result, detailed studies of large-scale flows are only possible in our local neighborhood. Such studies aim to verify that flows are gravitationally induced and to use them to measure the cosmic density parameter. Simulations of the nearby universe are ideal for calibrating such studies and for checking that they produce unbiased estimates of $`\mathrm{\Omega }_m`$.

The present paper is organized as follows. In Section 2 we describe both how we construct constrained initial conditions from the density field of the *IRAS* 1.2 Jy survey and how we carry out dark matter simulations from these initial conditions. Section 3 explains how halo catalogues and halo merging trees are built from the simulation outputs and summarises the phenomenological treatment of galaxy formation which we graft onto these trees, emphasising points where it differs from the treatment in Kauffmann et al. (1999a, hereafter K99). This section also compares our dark matter distributions to those found in unconstrained simulations and our simulated distributions of galaxy luminosity, colour and morphology to those observed in large surveys. Section 4 begins our detailed comparison with the local universe by matching simulated rich clusters object by object with real rich clusters. In Section 5 we explain how we generate ‘mock’ catalogues for direct comparison with the *IRAS* *PSCz* and *UZC* surveys. In Section 6 we use these catalogues to carry out a point-by-point comparison of the smoothed galaxy and mass density fields, while Section 7 extends the comparison to smaller scales using cross-correlation statistics. Section 8 displays mock versions of the *Mark III* catalogue of galaxy peculiar velocities to illustrate the application of our simulations to cosmic flow problems. Finally, we give a brief summary and discussion of our results in Section 9.

## 2 Dark Matter

Our modelling is based on large high resolution simulations of the evolution of the dark matter distribution in a spherical region surrounding the Milky Way. The first half of this section explains how we obtain initial conditions for these simulations; the second half describes the code used to follow their evolution.

### 2.1 Constructing initial conditions

We use the techniques developed by Kolatt et al. (1996, hereafter K96) to set up initial conditions such that at $`z=0`$ the evolved mass overdensity in the simulation is a suitably scaled version of the galaxy overdensity in the *IRAS* 1.2 Jy survey of Fisher et al. (1994, 1995) once both are smoothed on the same large scale. We take all the observed galaxies out to a redshift of 12,000 km s-1, weight each by the inverse of the survey selection function at its redshift, and smooth with a Gaussian of 1-D dispersion $`5h^1\text{ Mpc}`$ to obtain the galaxy density field at all points within a sphere of radius 10,000 km s-1. When smoothing we take care to account properly for the regions of space not sampled by the observational survey. These are primarily behind the Galactic Plane and in the strip which was not scanned by the *IRAS* satellite.

For each cosmology ($`\mathrm{\Lambda }`$CDM and $`\tau `$CDM as in K99) we obtain a ‘target’ mass overdensity field by scaling the galaxy overdensities so that their rms within the region out to 10,000 km s-1 is equal to the value obtained by smoothing the $`z=0`$ linear power spectrum of the particular cosmological model with a Gaussian of dispersion $`5h^1\text{ Mpc.}`$ As described by K96 the Eulerian Zel’dovich-Bernoulli equation can then be integrated back in time to give the linear density field at $`z=50`$ which gives rise to the target smooth density field at the present day. This linear field is then “Gaussianised”, i.e. it is mapped onto a new field such that the new density at each point is a monotonic function of the original overdensity and is Gaussian on the observed region with rms equal to that expected at $`z=50`$ from the smoothed theoretical power spectrum.

We now apply the Hoffman-Riback algorithm (Hoffman & Ribak 1991, 1992; Ganon & Hoffman 1993) to generate an initial displacement field on a 2563 grid which (i) is periodic on a simulation cube of side $`L=240`$ $`h^1`$ Mpc, (ii) is a Gaussian random field with power per mode equal to that expected at $`z=50`$ for the chosen cosmology for all wavenumbers between the fundamental, $`k_0=2\pi /L`$ and $`64k_0`$, and (iii) is constrained so that when smoothed with a Gaussian of comoving dispersion $`5h^1\text{ Mpc}`$ it reproduces the target initial density field everywhere within a sphere of radius 80 $`h^1`$ Mpc centred on the cube.

We supplement this constrained low-frequency displacement field with an unconstrained high-frequency displacement field constructed as follows. On a periodic cube of side 80 $`h^1`$ Mpc we set up an unconstrained Gaussian random field with power per mode equal to that expected in the chosen cosmology for all wavenumbers between $`64k_0`$ and $`2\pi /(0.7h^1\text{ Mpc})=343k_0`$. The latter corresponds to a wavelength equal to twice the mean interparticle separation of our high resolution region (see below). We replicate this smaller cube 27 times to obtain a high-frequency displacement field everywhere within the 240 $`h^1`$ Mpc cube.

Because the version of the integration code we use has vacuum boundary conditions, it works best on near-spherical regions. We therefore create an unperturbed but variable resolution particle load within a sphere of radius $`240h^1\text{ Mpc}\times \sqrt{3}/2=207.85h^1\text{ Mpc.}`$ We first fill the 240 $`h^1`$ Mpc cube with a uniform but irregular “glass” of equal mass particles. (See White (1995) for a discussion of such glass distributions.) We then border this cube with periodic replications and select the sphere which just encloses the central cube as the final region to be simulated. Within this large sphere, we excise all particles which will end up less than 80 $`h^1`$ Mpc from the center by $`z=0`$, and we replace them with a region of the same size and shape cut from a glass distribution with a mean interparticle separation of 0.35 $`h^1`$ Mpc. The particle masses are adjusted to ensure the correct mean mass density in each of the two regions. Table 1 lists the cosmological parameters and the particle numbers and masses for each of our two simulations.

The final initial condition (at $`z=50`$) is created by interpolating the high- and low-frequency displacement fields to the unperturbed position of each particle, displacing it, and assigning it a peculiar velocity proportional to its displacement in accordance with usual Zel’dovich approximation. We move the low-mass, high-resolution particles in the inner region using the sum of the two displacement fields. The high-mass, low-resolution particles in the outer region are perturbed using the low-frequency constrained displacement field only. The latter field must be periodically replicated in order to displace the low resolution particles outside the 240 $`h^1`$ Mpc central cube.

### 2.2 Simulating the Dark Matter

We carried out simulations from the above initial conditions using the parallel tree-code gadget (Springel et al. 2001). Both simulations were started with fixed comoving softenings of 70 $`h^1`$ kpc and 420 $`h^1`$ kpc for high- and low-resolution particles respectively. Once the corresponding physical comoving softenings reached 20 $`h^1`$ kpc and 120 $`h^1`$ kpc they were kept constant at these values.

For the $`\mathrm{\Lambda }`$CDM simulation, this resulted in particles in the densest regions of the high resolution region requiring roughly 75000 adaptive timesteps to reach redshift zero. The number of equivalent full timesteps (the total number of forces computed divided by the number of particles) was much smaller however, of order 400. Both runs were performed on 512 processors of the CRAY T3E at the Computer Center (RZG) of the Max-Planck Society at Garching. The $`\mathrm{\Lambda }`$CDM simulation required about 104 processor-hours to complete (the total CPU time consumption was about 50000 hours), while the $`\tau `$CDM simulation was somewhat faster. We saved the particle positions and velocities at 41 times logarithmically spaced in expansion parameter from $`z=10.3`$ to $`z=0`$.

In Fig. 1 we show the present-day dark matter distribution in our $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models. The projected region is a 30 $`h^1`$ Mpc thick slice parallel to the supergalactic plane and extending from -15 $`h^1`$ Mpc to 15 $`h^1`$ Mpc in SGZ. The side of the slice is 180 $`h^1`$ Mpc long, and the Milky-Way is at its centre. Well-known local structures are easily identified, in particular, the Great Attractor region, the Coma-Great Wall structure, the Pisces-Perseus supercluster filament, and the Local Void.

Many of the known nearby rich clusters, including Virgo, can also be identified, even though the smoothing kernel used to generate the initial conditions is substantially larger than a rich cluster and such clusters are actually quite poorly represented in the *IRAS* 1.2 Jy survey.

One can start to see the low resolution region at the corners of these images which extend more than 80 $`h^1`$ Mpc from the centre. In the remainder of this paper we shall concentrate on analysis of the region within 80 $`h^1`$ Mpc and we will exclude from consideration the few dark haloes near the boundary which are contaminated with high-mass particles.

As a test that our simulations are correctly reproducing the small scale structure expected in each cosmology, Fig. 2 compares their halo mass functions at $`z=0`$ with those expected according to the large ensemble of unconstrained simulations analysed by Jenkins et al. (2001). As discussed in the next section, we identify haloes in our simulations using a friends-of-friends group-finder with $`b=0.2`$ (Davis et al. 1985). The analytic fit which Jenkins et al. (2001) use to summarise their data for such FOF mass functions is also an excellent description of our own halo data. Note that no parameters were adjusted in making this comparison. Significant discrepancies occur only for the largest mass haloes. These are plausibly a consequence of the constraints we impose on our initial conditions.

## 3 Putting in the Galaxies

We follow the assembly, the evolution and the merging of galaxies within our dark matter simulations using phenomenological prescriptions based closely on those of K99. In this section we give an outline description of the scheme as a whole. We then discuss particular aspects of the modelling where we have found it advantageous to modify the K99 prescriptions. Finally we describe how we set the parameters which describe uncertain physical processes (for example, the efficiencies of star formation and feedback) using observed properties of galaxies such as the Tully–Fisher relation and the luminosity function.

Our galaxy formation scheme is based on the idea, originally introduced by White & Rees (1978), that galaxies result from the condensation of diffuse gas at the centres of dark haloes. As a result all galaxies are found within dark haloes, either as the principal central object or as satellites in orbit about it. The first step is thus to identify all the dark haloes in each snapshot of the simulation and to construct a tree which describes how they merge together from snapshot to snapshot. The galaxy population can be reconstructed by following evolution down the tree from the earliest times to the present. Within each halo at each time diffuse gas cools onto the disk of the central galaxy, stars form from cold gas both in this galaxy and in its satellites, feedback from this star formation reheats some of the cold gas, and satellites occasionally merge with the central galaxy, modifying its morphology if they are sufficiently massive. Galaxy motions are followed in a dynamically consistent way by attaching each galaxy to a dark matter particle. For a central galaxy this identification is with the most bound particle of the halo. The particle identification remains fixed if this halo merges with a more massive one so that the galaxy becomes a satellite.

We now briefly discuss aspects of this modelling, emphasising aspects where our implementation differs in detail from that of K99.

### 3.1 Tree structure

As in K99 we build a halo catalogue for each of the 41 snapshots of each simulation using a friends-of-friends algorithm with linking length 0.2 times the mean interparticle separation (Davis et al. 1985). For further analysis we consider only haloes with 10 or more particles; smaller haloes frequently disappear in later snapshots. The more massive haloes are linked from each snapshot to the next to form a tree exactly as in K99. Note that although we can pinpoint the location and mass of haloes quite well down to our limit of 10 particles, we can only reliably track the merging history of haloes more massive than about 100 particles. Lower mass haloes can typically be traced back through a few snapshots only before they vanish below the resolution limit. As a result, we can get reasonable estimates for the mass of the central galaxy in a 10 to 20 particle halo, but we cannot infer anything about its morphology for halo masses below 100.

This leads us to define the *luminosity* resolution limit of our simulated galaxy catalogues to be the typical luminosity of the central galaxy in a 10 particle halo, and the *morphology* resolution limit to be the typical luminosity of the central galaxy in a 100 particle halo. At $`z=0`$ the circular velocities of 10 and 100 particle haloes are 45 km s-1 and 103 km s-1 respectively in our $`\mathrm{\Lambda }`$CDM simulation. For $`\tau `$CDM the corresponding numbers are 76 km s-1 and 160 km s-1. (We define the circular velocity $`V_{200}`$ of a halo to be the value measured at $`R_{200}`$, the radius of the sphere centred on the most bound particle within which the mean density is 200 times the critical density of the Universe.) In the galaxy catalogues we construct below, the corresponding *luminosity* and *morphology* limits are $`M_B=16.27`$ and $`18.46`$ for $`\mathrm{\Lambda }`$CDM, $`M_B=18.50`$ and $`20.60`$ for $`\tau `$CDM. In this paper, “resolution limit” without a qualifier will always denote the *luminosity* resolution limit.

Galaxies are of three types in our evolution scheme: central galaxies which are always identified with the most bound particle (MBP) of their halo, satellite galaxies which are identified with other halo particles, and field galaxies whose particle is not a member of any catalogued halo. When setting up the galaxy population at redshift $`z_0`$, we adopt the following procedures which differ slightly from those of K99.

1. If no particle in a halo was a galaxy in the previous snapshot at $`z_1>z_0`$, then we create a new galaxy with default properties based on a simple infall model for the earlier evolution of the halo up to $`z_0`$.
2. If no particle in the halo was the central galaxy of a halo at $`z_1`$, we identify the most massive (in terms of total stellar and cold gas mass) of the galaxies which do populate the halo as its central galaxy, we transfer the properties of this galaxy to the halo MBP, and we delete the particle previously assigned to it from the galaxy list.
3. If at least one particle was a central galaxy at $`z_1`$, we identify the new central galaxy as the one whose halo donated the most mass to the $`z_0`$ halo. Its properties are transferred to the MBP and its previous particle identification is deleted from the galaxy list. The other ‘old’ central galaxies become satellite galaxies in the new halo.
4. Galaxies from $`z_1`$ which are part of no halo at $`z_0`$ are designated as field galaxies. These objects are always few in number and are near the resolution limit of the simulation. They play no significant role in the results we present below.

This scheme is designed to ensure that a halo central galaxy always inherits the properties of the most plausible (usually the most massive) candidate for its progenitor in the previous snapshot. We have tested a variant in which inheritance is based on luminosity rather than mass and find it to make no discernable difference. The scheme corrects a loophole in the K99 procedures which allowed relatively massive haloes occasionally to find no progenitor for their central galaxies. As a result they were assigned relatively low-mass central galaxies and they appeared as outliers in the Tully–Fisher relation.

### 3.2 Cooling, star formation and feedback

The 41 snapshots stored from $`z=10.42`$ to $`z=0`$ sample the evolution of the system more coarsely than is desirable when simulating the evolution of the galaxy population. As in K99 we circumvent this by dividing the time between each pair of snapshots into 50 small steps, which we use for integrating the simple phenomenological equations which describe gas infall and cooling, star formation, feedback and satellite merging. We have checked that doubling the number of steps does not change our results.

We treat cooling using exactly the same prescriptions as K99. These assume that gas cools out of the hot atmosphere of a halo onto the disk at its centre at a rate which depends only on redshift, on halo mass and on halo hot gas content. For haloes in which cooling is rapid (at high redshift, relatively small mass, and large hot gas fraction) the atmosphere is assumed to condense out on the dynamical time of the halo. When cooling times are longer, the condensation rate is taken from a simple cooling flow model. We calculate cooling times using the collisional ionisation cooling curve of Sutherland & Dopita (1993) assuming primordial abundances. We do not attempt to follow chemical enrichment. We treat the global baryon fraction of the models $`f_{\mathrm{bar}}\mathrm{\Omega }_{\mathrm{bar}}/\mathrm{\Omega }_0`$ as an adjustable parameter.

These prescriptions give rise to a problem that has been noticed at least since the work of White & Frenk (1991). In galaxy cluster haloes they predict strong cooling flows (similar to those apparently observed in clusters like Perseus) which deposit large amounts of cold gas at cluster centres. The same star formation prescriptions which produce reasonable galaxies at the centre of $`10^{12}\mathrm{M}_{}`$ haloes then predict massive, luminous, star-forming galaxies which are quite inconsistent with the observed properties of the central galaxies in clusters. This is just another facet of the well-known problem of the fate of the gas apparently deposited by cooling flows (Fabian et al. 1991; Allen & Fabian 1997).

A number of authors have tried to deal with this problem in the present context by building more elaborate models for the structure and thermodynamic history of the gas in dark matter haloes (Somerville & Primack 1999; Cole et al. 2000, hereafter SP99, C00). We have tried implementing the prescriptions advocated by SP99 and find that in our own models they produce only a rather limited improvement. As a result, we have preferred to retain the simpler but cruder ad hoc solution of K99 – we switch off all cooling of gas in haloes with circular velocities exceeding 350 km s-1. This is equivalent to assuming that gas processed through cluster-like cooling flows does not end up being available for normal star formation. It also ensures that no disks are formed with circular velocities exceeding the largest observed for real disk galaxies.

Our models for star formation and feedback are also identical to those of K99. The star formation rate in a galaxy is assumed to be proportional to its cold gas mass divided by its dynamical time:

$$\dot{M}_{}=\alpha \frac{M_{\mathrm{cold}}}{t_{\mathrm{dyn}}}=100\alpha H(z)M_{\mathrm{cold}}$$
(1)
where $`\alpha `$ is a free parameter to be fitted against observation, $`H(z)`$ is the Hubble parameter, and the second equality follows from a very simple model for the typical dynamical time of galaxy disks, taken to be $`0.1R_{200}/V_{200}`$. We use the population synthesis models of Bruzual & Charlot (1993) to compute the luminosities of our galaxies from their star formation histories assuming a Scalo IMF.

Feedback is assumed to act by reheating some of the cold gas in the disk. We write the reheating rate as:

$$\dot{M}_{\mathrm{reheat}}=ϵ\frac{4}{3}\frac{\dot{M}_{}\eta _{\mathrm{SN}}\mathrm{E}_{\mathrm{SN}}}{\mathrm{V}_{200}^2}$$
(2)
where $`\eta _{\mathrm{SN}}=5\times 10^3\mathrm{M}_{}^1`$ is the number of supernovae produced per solar mass of stars in a Scalo IMF, $`\mathrm{E}_{\mathrm{SN}}10^{51}\mathrm{erg}`$ is the mean energy released per supernova, $`\mathrm{V}_{200}`$ is the circular velocity of the halo at $`R_{200}`$, and $`ϵ`$ is a free parameter to be tuned to match observations. Since it is unclear whether the reheated gas will remain in the halo or will be expelled, a model of each type was discussed by K99. We have tried both schemes in our own simulations and, as in K99, we find that $`\mathrm{\Lambda }`$CDM fits the observations better with retention feedback. For $`\tau `$CDM , however, we are able to fit the Tully–Fisher relation with less scatter using the ejection scheme. This may appear surprising in view of fig. 6 of K99, but it is accounted for by our different prescription for computing $`V_{\mathrm{disk}}`$ (see below).

### 3.3 Dust

We use the simple model of dust extinction proposed by SP99. We compute the face–on $`B`$-band optical depth as:

$$\tau _B=\tau _{B,}\left(\frac{L_B}{L_{B,}}\right)^\beta $$
(3)
with $`\tau _{B,}=0.8`$, $`L_{B,}=6\times 10^9\mathrm{L}_{}`$ and $`\beta =0.5`$, as given by Wang & Heckman (1996). We then use the galactic extinction curve of Cardelli et al. (1989) to calculate the optical depth in other bands. We pick a random inclination for each galaxy, and apply a dust correction to its disk component only. This is calculated using a “slab” model as:

$$A_\lambda =2.5\mathrm{log}_{10}\left(\frac{1\mathrm{e}^{\tau _\lambda \mathrm{sec}\theta }}{\tau _\lambda \mathrm{sec}\theta }\right)$$
(4)
K99 apply a dust correction only to the disks of galaxies with star formation rates $`\mathrm{SFR}>0.5M_{}/\mathrm{yr}`$. Here, we correct all disks, regardless of their SFR. In practice, the differences induced by the change in prescription are not large, and although the new prescription gives somewhat better results, it is still unable to eliminate the apparent excess of galaxies at the bright end of the luminosity function.

### 3.4 Merging

As in K99, we allow a satellite galaxy to merge only with the central galaxies of its halo. Each satellite galaxy carries a merging time counter which is initialized with an estimate of the timescale for orbital decay through dynamical friction. This counter is reinitialized each time the satellite’s halo undergoes a major merger or is accreted by a more massive halo. Otherwise it decreases steadily and the satellite is assumed to merge with the central galaxy when it reaches zero.

The initial timescale is based both on analytic estimates (Binney & Tremaine 1987) and on fits to numerical simulations (Navarro et al. 1995). Some problems with these formulae have been pointed out by S00, who note that when haloes of comparable mass merge the inferred merging times for their central galaxies are typically much shorter than the time actually taken in high resolution simulations of the process. This could, in principle, lead to an overmerging problem. There is also disagreement in the literature on the appropriate formula for the Coulomb logarithm when estimating dynamical friction timescales. Somerville & Primack (1999) take $`\mathrm{ln}\mathrm{\Lambda }=\mathrm{ln}(1+(M_{\mathrm{halo}}/M_{\mathrm{sat}})^2)`$ while K99 and S00 use $`\mathrm{ln}\mathrm{\Lambda }=\mathrm{ln}(M_{\mathrm{halo}}/M_{\mathrm{sat}})`$ and $`\mathrm{ln}\mathrm{\Lambda }=\mathrm{ln}(1+M_{\mathrm{halo}}/M_{\mathrm{sat}})`$ respectively. Here, we adopt the formula of SP99 since it improves the bright–end behaviour of our luminosity functions.

We attempt to avoid any overmerging problem by ensuring that galaxies never merge on timescales shorter than their halo crossing times. If $`t_{\mathrm{fric}}`$ is the merging timescale derived from dynamical friction considerations, we initialise the merging time counter using:

$$t_{\mathrm{merging}}=t_{\mathrm{dyn}}(2+\frac{t_{\mathrm{fric}}}{t_{\mathrm{dyn}}})$$
(5)
This encapsulates the physical assumption that the satellite cannot merge with the central galaxy fewer than two crossing times after it first enters the halo.

When the merging time counter vanishes, we add the properties of the satellite to those of the central galaxy, and we remove the satellite from the galaxy list. We introduce a bulge formation threshold $`f_{\mathrm{bulge}}`$ as a free parameter, used to match the observed abundance of galaxies by morphology. Mergers in which $`M_{\mathrm{sat}}/M_{\mathrm{central}}>f_{\mathrm{bulge}}`$ are considered major mergers and create an elliptical galaxy containing the stars of both progenitors, as well as stars formed from their cold gas component in a burst assumed to last $`10^8`$ years. A merger with $`M_{\mathrm{sat}}/M_{\mathrm{central}}<f_{\mathrm{bulge}}`$ is considered minor merger. We add the cold gas of the satellite to that of the central galaxy. If no bulge already exists in the central galaxy, we create a new bulge with the stars of the merging satellite, otherwise we simply add the stars to the existing bulge. Such a scheme is sufficient for our current purposes, more realistic treatments include that of Somerville et al. (2001).

### 3.5 Setting parameters

The models set out above have four free parameters: the baryon fraction $`f_{\mathrm{bar}}`$, the star formation efficiency $`\alpha `$, the supernova feedback efficiency $`ϵ`$, and the bulge formation threshold $`f_{\mathrm{bulge}}`$. We choose these parameters within the range which seems physically plausible in such a way as to maximise the agreement of the resulting galaxy populations with observation. Notice that this procedure is only feasible in a scheme of the type we have implemented where changes in such physical parameters do not require the supercomputer simulations to be repeated. We attempt to match the following observations:

1. The velocity–luminosity relation for simulated spirals *before dust correction* should match the published $`I`$-band relation of Giovanelli et al. (1997), which is already corrected for internal extinction.
2. The $`b_j`$-band luminosity function should be in rough agreement with those of the 2dF and SDSS galaxy surveys (Folkes et al. 1999; Blanton et al. 2001). Note that our poor resolution prevents meaningful comparison of the faint end slopes so we concentrate on getting reasonable fits for the brighter galaxies.
3. Galaxies with $`V_{\mathrm{disk}}`$ $`220\text{km s}\text{-1}`$ should host $`10^{11}\mathrm{M}_{}`$ of stars and a few times $`10^9\mathrm{M}_{}`$ of cold gas (c.f. K99, SP99). The latter includes both molecular and atomic gas.
4. The morphology distribution should agree with the observation: this is discussed further in section 3.5.4.

In Table 2 we list our estimated best parameters given these constraints, which are described in more detail in later subsections. Note that the baryon fraction in the $`\tau `$CDM model is high and inconsistent with standard big–bang nucleosynthesis (Burles et al. 1999), although it matches the baryon fraction observed in massive clusters (Ettori & Fabian 1999). Interestingly, if we adopt the preferred cooling model of SP99, we can fit the Tully–Fisher relation with baryon fractions of order $`f_{\mathrm{bar}}0.08`$ and $`f_{\mathrm{bar}}0.15`$ for $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM respectively, albeit with overbright ”cD” galaxies. This same trend in the required $`f_{\mathrm{bar}}`$ was also found by SP99. Note that taking into account the chemical enrichment of the hot halo gas would boost the cooling rate at late times, and allow the observed Tully–Fisher relation to be matched for a lower baryon fraction. (See, for example, models n and c of van Kampen et al. 1999.) We also find that lower baryon fractions are allowed if we assume a substantial fraction ($`20\%`$) of the mass of newly formed stars is returned to the cold gas of the galaxy, as expected from SNe and stellar winds (see C00). We did not include this possibility in the models described here, as it makes more sense to use it in combination with a full chemical enrichment scheme.

In Fig. 3 we show the present-day distribution of galaxies in our $`\mathrm{\Lambda }`$CDM model for this choice of parameters. Galaxies are overplotted on the dark matter distribution using symbols whose size increases with luminosity and whose colour represents $`BV`$ index. The geometry of the slice is the same as in Fig. 1. Note the predominance of red galaxies in clusters and that of blue galaxies in the field. Notice also that the local voids are indeed nearly empty of galaxies.

#### 3.5.1 Tully–Fisher relations

In K99 the Tully–Fisher relation was plotted for central galaxies with ‘spiral morphology’ assuming the disk rotation velocity $`V_{\mathrm{disk}}`$ to be equal to $`V_{200}`$ the circular velocity of the host halo. When we used the same prescription, we found a number of low luminosity galaxies, off by more than half a magnitude both from the observed relation and from the mean relation for the simulated galaxies. Similar outliers can be seen in fig. 6 of K99 – at $`V_{200}`$ $`300\text{km s}\text{-1}`$ some of their ‘Sb/Sc’ galaxies are too faint by more than one magnitude. This problem arises because $`V_{200}`$ $``$$`V_{\mathrm{disk}}`$ is a poor approximation for a subset of galaxies with unusual halo assembly histories. To get a better estimate for $`V_{\mathrm{disk}}`$, we proceed as follows. Throughout the evolution of our galaxy population we accumulate $`\mathrm{\Delta }M_{\mathrm{stars}}V_{200}^2`$ for each star-forming central disk where $`V_{200}`$ is the circular velocity of its current halo. The circular velocity of the disk at any redshift is then defined by

$$V_{\mathrm{disk}}^2=\frac{\underset{z_i}{}\mathrm{\Delta }M_{\mathrm{stars},\mathrm{z}=\mathrm{z}_\mathrm{i}}V_{200}^{2}{}_{\mathrm{z}=\mathrm{z}_\mathrm{i}}{}^{}}{M_{\mathrm{stars},\mathrm{disk}}},$$
(6)
where the sum is over all redshift intervals since the last major merger and $`M_{\mathrm{stars},\mathrm{disk}}`$ is the stellar mass of the disk. Satellite galaxies keep the disk rotation velocity they had when they last were central galaxies.

In Fig. 4 we plot Tully–Fisher relations for simulated galaxies with $`100<`$$`V_{\mathrm{disk}}`$ $`<300\text{km s}\text{-1}`$ and $`1.5<M_{B,\mathrm{bulge}}M_{B,\mathrm{total}}<2.2`$ (appropriate for Sb/Sc galaxies according to Simien & de Vaucouleurs 1986). This is the same criterion used by K99 except that here we do not restrict ourselves to central galaxies. We compare our results with the fit to the observational data given by Giovanelli et al. (1997):

$$M_I5\mathrm{log}(h)=21.00\pm 0.027.68\pm 0.24(\mathrm{log}(W)2.5),$$
(7)
where $`W=2\times V_{\mathrm{disk}}\text{.}`$ The dashed lines show $`\pm 1\sigma `$ of what Giovanelli et al. term the “intrinsic scatter” of their relation. This relation has been corrected for internal extinction so we have not included any effects of dust in our model predictions.

Even with our revised definition of $`V_{\mathrm{disk}}`$, there are still a few galaxies with low $`I`$-band luminosity which show up below the observed relation. We could remove the five low-luminosity galaxies at $`\mathrm{log}(W)2.4`$ in $`\mathrm{\Lambda }`$CDM by adopting a more stringent morphological selection criterion, for example $`1.6<M_{B,\mathrm{bulge}}M_{B,\mathrm{total}}<2`$.

We have checked that the galaxies lying below (more than 0.7 magnitude) the expected relation are redder than ”typical” spirals: their mean $`BV`$ index is 0.9 and 1.0 in $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM respectively, instead of the 0.7 found for “normal” spirals. These galaxies are older systems, for example, halo satellites which still satisfy the morphological selection criterion in the $`B`$-band, have kept the disk rotation velocity when they last were a central galaxy, but where star formation has stopped because of gas exhaustion, resulting in a relatively low stellar mass and $`I`$-band luminosity. They might thus correspond to a subset of observed S0’s.

Giovanelli et al. (1997) quote an intrinsic scatter of $`ϵ_{int}=0.28x+0.26`$ for their Tully–Fisher relation, with $`x=\mathrm{log}W2.6`$. This is consistent with our $`\mathrm{\Lambda }`$CDM model, but is exceeded by the $`\tau `$CDM model at the bright end. We find the dependence of the scatter on our rather uncertain morphological selection criterion to be quite strong. A larger scatter is expected in high-density cosmologies because of their late structure formation, as pointed out by Buchalter et al. (2001).

#### 3.5.2 Luminosity functions

In Fig. 5 we compare the $`b_j`$-band luminosity functions (LF) of our simulations with those derived from the 2dF and SDSS redshift surveys. The vertical lines correspond to our *luminosity* and *morphology* resolution limits. As already stressed, a major difficulty in reproducing observed luminosity functions with the kind of modelling employed here is in establishing a sufficiently strong cut-off at high luminosities by avoiding the formation of overly bright central galaxies in groups and clusters. Despite our best efforts and our ad hoc cooling switch at $`V_{200}`$ $`=350`$ km s-1, this has not been completely successfully achieved in the models, where we still have a bright galaxy at $`M_{b_j}24.5`$, possibly due to the overmerging issue discussed in Section 3.4. A second major difficulty is in suppressing the luminosity of galaxies in low mass haloes sufficiently to reproduce the flat faint-end slope of observed LFs. We have not stressed this issue here because our resolution does not allow us to address it adequately.

The $`\mathrm{\Lambda }`$CDM LF shows a deficit of $`L_{}`$ galaxies, with respect to the observations, whereas the $`\tau `$CDM LF produces somewhat too many such galaxies. As pointed out by Kauffmann et al. (1993), this difference is a consequence of the differing number of $`V_{200}200`$ km s-1 haloes in the two cosmologies, given our requirement that both should reproduce the observed Tully–Fisher relation. The global shape of the $`\tau `$CDM LF agrees somewhat better with the observations. The shapes of both our luminosity functions are closer to the observations than those given in K99. This is a result of our tuning our cooling and merging prescriptions as noted in previous sections (see also SP99 and C00).

The results of S00 suggest that the luminosity function shape would be improved if we used simulations of sufficiently high resolution to follow the merging of galaxies directly, rather than having to use a phenomenological model. The “knee” in our simulated $`b_j`$-band LFs occurs at $`M_{b_j}22`$ in both our models. 93% and 96% of the galaxies brighter than these limits are central galaxies of haloes in $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM respectively. Of these central galaxies, 73% and 16% belong to haloes with $`V_{200}>350`$ km s-1, showing that our LF results in this range are being affected strongly by the cooling cut-off and by the merging which can brighten galaxies above it.

#### 3.5.3 Masses for Milky Way look-alikes

We select Milky Way look-alikes as the central galaxies of haloes with $`200<`$$`V_{\mathrm{disk}}`$ $`<240\text{km s}\text{-1}`$ which also satisfy the “Sb/Sc” criterion we used when plotting the Tully–Fisher relation. Table 3 presents their mean stellar and cold gas masses together with their $`B`$-band and $`I`$-band absolute magnitudes, their colours, their disk SFR’s and their abundance in the simulations.

Our two models show similar cold gas and stellar masses, quite comparable with those estimated for the Milky Way, although their present-day star formation rates differ by a factor of two. We find an average SFR of 1.8 $`M_{}/\mathrm{yr}`$ for disks of Milky Way–like galaxies in $`\mathrm{\Lambda }`$CDM and 3.4 $`M_{}/\mathrm{yr}`$ for such disks in $`\tau `$CDM. This reflects the later formation of structure in $`\tau `$CDM which requires more star formation at late times to get enough stars to reproduce the Tully–Fisher relation. Both values are in reasonable agreement with the 2 to 3 $`M_{}/\mathrm{yr}`$ estimated for the mean star formation rate in the Milky Way’s disk over the last few Gyr (Rocha-Pinto et al. 2000b, a).

The difference in the number of galaxies with $`200<`$$`V_{\mathrm{disk}}`$ $`<240`$ km s-1 in $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM is substantial but is fully accounted for by the difference in the halo mass functions.

#### 3.5.4 Morphologies

As in Kauffmann et al. (1993), we assign a morphological type based on the $`B`$-band bulge–to–disk ratio: galaxies with $`L_{\mathrm{B},\mathrm{bulge}}/L_{\mathrm{B},\mathrm{disk}}>1.52`$ are ellipticals, S0’s have $`0.68<L_{\mathrm{B},\mathrm{bulge}}/L_{\mathrm{B},\mathrm{disk}}<1.52`$, spirals have $`L_{\mathrm{B},\mathrm{bulge}}/L_{\mathrm{B},\mathrm{disk}}<0.68`$ and irregulars have no bulge.

We adjust our major merger threshold parameter so that the overall morphology fractions in our simulations match those derived by Baugh et al. (1996) from the observations of Loveday (1996): 13%, 20% and 67% for Es, S0s, and Sps + Irrs. We estimate these ratios in the simulations using all galaxies brighter than our morphology resolution limit. Taking $`f_{\mathrm{bulge}}=0.1`$ we obtain 12%, 19% and 69% as the relative fraction of the three morphological classes in $`\mathrm{\Lambda }`$CDM, and 16%, 14% and 70% in $`\tau `$CDM. For $`\tau `$CDM our S0 fraction is somewhat smaller than the observed 20%. This may be another hint for overmerging: central galaxies generally spend little time within the S0 bulge–to–disk ratio limits, either growing bigger disks or evolving into Es via mergers. Both the cooling and the merging processes are more rapid in $`\tau `$CDM, leading to a reduced number of S0s. Given the different merging prescriptions used by different authors, one must be careful when comparing the $`f_{\mathrm{bulge}}`$ values they derive. For example, SP99 and K99 both fit the observed morphological ratios with $`f_{\mathrm{bulge}}=0.25`$, but with a different formula for the dynamical friction timescale than is used here. As C00 emphasize, the merger simulations of Walker et al. (1996) and Barnes (1998) suggest $`0.1<f_{\mathrm{bulge}}<0.3`$.

In Fig. 6, we plot mean morphology fractions as a function of the maximum $`B`$ absolute magnitude of the galaxies considered. The fraction of bulge–dominated galaxies increases with luminosity. The most luminous galaxies, which are mainly the central galaxies of clusters, are all ellipticals. The variation of the morphological mix with luminosity is smooth, except for a steep rise in the elliptical fraction at $`M_\mathrm{B}22`$. This is a consequence of our ad hoc suppression of cooling onto central galaxies in haloes with $`V_{200}>350\text{km s}\text{-1}\text{.}`$ This suggests that the observed predominance of bright ellipticals in clusters is related to the observed inability of cooling flows to generate new disks in such systems.

#### 3.5.5 Colour distributions

Although we do not use colours when setting the free parameters of our models, we briefly discuss here the colour distributions of our galaxies. The solid histogram in Fig. 7 shows the distribution in $`BV`$ colour of all the galaxies brighter than our morphology limit. Two peaks are evident in these distributions which are separated when we plot individual distributions for ellipticals and spirals, as defined in 3.5.4. There is good agreement with K99 and with the colours derived from the RC3 catalogue of de Vaucouleurs et al. (1991) (shown as the dashed histogram). If we restrict ourselves to galaxies residing in massive clusters (those with $`V_{200}>500\text{km s}\text{-1}`$), we obtain the dotted plots, showing that cluster galaxies are redder than their field counterparts (see also fig. 16 of S00).

## 4 Nearby Clusters

This section begins the detailed comparison of our simulations with the individual structures seen in the local universe. We start with a one–to–one comparison of the largest simulated clusters with their nearby counterparts. It is surprising that this comparison is possible at all because rich clusters are systematically under-represented in the *IRAS* 1.2 Jy survey from which we constructed our initial conditions. The next subsection shows that a convincing identification of simulated and real clusters is, nevertheless, possible, and we go on in subsequent subsections to explore this correspondance in terms of cluster dark matter properties, galaxy luminosity functions, and cluster mass–to–light ratios.

### 4.1 Observed versus simulated clusters

Among nearby rich clusters we select Coma, Virgo, Centaurus, Hydra and Perseus as particularly interesting for comparison with our simulations. All are relatively massive and lie within $`8000\text{km s}\text{-1}`$. We identify them with their simulated counterparts as follows. We take the observed coordinates and redshift from the nasa/ipac extragalactic database 111http://nedwww.ipac.caltech.edu/forms/byname.html, and then look in the simulations for the most massive group less than 7 $`h^1`$ Mpc away from the expected redshift–space position. This search radius has been chosen arbitrarily to ensure that we recover all five selected clusters.

The positional agreement between real and simulated clusters, given in Table 4, is quite good except for our $`\mathrm{\Lambda }`$CDM Centaurus cluster, which is 6.7 $`h^1`$ Mpc away from the expected position. The very good match for Coma is particularly surprising: the error in position is less than 2 $`h^1`$ Mpc in both models, despite the fact that the cluster is close to the edge of the high-resolution region, where the constraints are noisy due to the relatively poor sampling of the 1.2 Jy survey.

### 4.2 Cluster dark matter

In Table 4 we list for each cluster the total dark matter mass (as given by our group-finder), the virial mass (defined as the mass within $`R_{200}`$), both in units $`10^{14}`$ $`h^1`$ $`M_{}`$, the virial radius in $`h^1`$ Mpc, and the line–of–sight galaxy velocity dispersion in units of km s-1. We take the real data from Girardi et al. (1998). In the observed virial mass column, the first figure is $`M_{200}`$, inferred from the virial theorem applied to galaxies within $`R_{200}`$ but without including a surface term (the traditional virial estimate). The second figure is $`M_{200,\mathrm{corr}}`$, which includes a term to correct for the surface “pressure” and should be more accurate (see Girardi et al. 1998 for details). We compute line–of–sight velocity dispersions in our simulations using all galaxies above our resolution threshold. The first figure is obtained by considering all galaxies in the cluster halo, the second only those galaxies within the virial radius. There is no significant difference.

It is striking that the ranking of clusters by mass in the simulations agrees with their ranking in the real universe. For individual clusters the simulated $`M_{200}`$ is usually within a factor of 2 of the observed value and is often better. If we assume that $`M_{200,\mathrm{corr}}`$ is the better observational estimator, then for $`\mathrm{\Lambda }`$CDM Coma, Hydra and Perseus agree well, but Virgo is too massive by a factor of two, while for $`\tau `$CDM Virgo agrees well but Perseus and Coma are too massive by a factor of two whereas Hydra is undermassive by a factor exceeding 2. This casts some doubt on the identification of Hydra in $`\tau `$CDM – however, it lies only $`4h^1\text{ Mpc}`$ away from the expected location and the next major (and more massive) cluster is $`7h^1\text{ Mpc}`$ away.

For the line–of–sight velocity dispersion the ranking of all five clusters agrees with observation in both cosmologies. The individual values also agree quite well in most cases. The larger deviations correspond to the discrepant masses just discussed. Recall that Centaurus has a well–known bimodal velocity structure (Lucey et al. 1986), reflecting the fact that the groups Cen30 and Cen45 are currently merging (Churazov et al. 1999). The quoted observed value is inferred if the two groups are considered together.

### 4.3 Luminosities and M/L ratios

For each cluster, we identified the most luminous galaxy (in $`B`$) and calculated its $`BV`$ colour. In Table 5, we compare the results with the data of Sandage (1972, hereafter S72). For each cluster the first line refers to the brightest cluster galaxy in our $`\mathrm{\Lambda }`$CDM model, the second to the BCG in $`\tau `$CDM, and the last to the observed BCG. In (almost) all cases our simulated BCGs are brighter than observed. (The BCG in the low mass $`\tau `$CDM ‘Hydra’ is an exception.) The discrepancy is less than about a magnitude for all objects except the well-known peculiar galaxy NGC 1275 in Perseus, which is more than one and a half magnitude fainter than our prediction. The observed values quoted have been corrected for galactic extinction as follows. Most corrections computed by S72 used a simple, latitude dependent model, except for Perseus where S72 estimated the correction based on the $`BV`$ index of NGC 1275, and found $`A_V=0.3`$ and $`A_B=0.4`$. Given the peculiar spectrum of the galaxy, these values are quite uncertain. The $`V`$-band magnitude of the BCG of Centaurus found in Abell et al. (1989) has not been corrected for extinction. Recall that, owing to their low galactic latitude, extinction is a major issue for the $`B`$-band magnitude of Centaurus and Perseus. For these two clusters, we have therefore taken the estimation of galactic extinction given by Schlegel et al. (1998) from their $`100\mu \mathrm{m}`$ all-sky map of dust emission: they find $`A_V=0.38`$ and $`A_B=0.49`$ for Centaurus and $`A_V=0.57`$, $`A_B=0.74`$ for Perseus.

The discrepancy in brightness between the simulated BCGs and the data may be due in part to the luminosities of the observed BCGs being underestimated (c.f. Uson et al. 1991; Gonzalez et al. 2000), and in part to the overmerging problem discussed in 3.4. Note that Coma has two very bright galaxies in a binary in its core. The luminosities of our simulated Coma BCGs are similar to the sum of those of NGC 4889 and NGC 4874.

The colour indices of our simulated BCGs are all very red, reflecting their old stellar populations. (Note that we use solar metallicity population synthesis models when estimating the observed properties of our stellar populations.) They agree well with the observed BCG colours in all cases except that of NGC 1275. The anomalous A-type spectrum of this galaxy is thought to result from star formation associated either with the strong cooling flow in Perseus or with an apparent ongoing merger with a gas-rich galaxy.

In Fig. 8, we give $`m_V`$ for the ten most luminous galaxies in our simulated Coma and Virgo clusters. They are compared to the data of Gudehus (1995). The BCG has been discussed in the previous section, and is typically brighter in the models than in the simulations. The remaining nine galaxies are fainter than observed by about one magnitude. It appears that choosing parameters to fit the Tully–Fisher relation has resulted in too few stars forming in cluster galaxies. Possible reasons might be a systematic underestimate of $`V_{\mathrm{disk}}`$ (and thus the assignment of low $`L`$) for the central galaxies of our haloes, or an underestimate of the effects of dust (and thus of the stellar masses) in these same ‘Tully–Fisher’ galaxies. The first explanation would affect the luminosities of all galaxies and would so shift the overall LFs of the models, the second would affect primarily cluster galaxies and would have little effect on the global LFs. The fact that Diaferio et al. (2001) get a good fit (actually slightly overluminous) to the bright end of the LF’s of the CNOC1 clusters, while S00 slightly underpredict the total luminosity of Coma, suggests that the problem is not fundamental but lies in the details of the phenomenological modelling. The discrepancy in our own models is also evident in the cluster mass-to-light ratios as we now show.

We compare simulated mass–to–light ratios in the $`b_j`$-band with those derived by Girardi et al. (2000). We approximate the $`b_j`$ magnitude by taking the “galaxy colour relation”: $`b_j=B0.35(BV)`$ (Blanton et al. 2001). We estimate two ratios for our simulated clusters, one using the total mass and light in our simulated haloes, the other using values within their virial radii:

$$\mathrm{{\rm Y}}_1=\frac{M_{\mathrm{tot}}}{L_{b_j,\mathrm{tot}}},$$
(8)
$$\mathrm{{\rm Y}}_2=\frac{M_{\mathrm{vir}}}{L_{b_j,<\mathrm{Rvir}}}.$$
(9)
We give all ratios in solar units for $`h=1`$.

If $`\varphi _{b_j}`$ is our luminosity function for the $`b_j`$-band, we include an estimated contribution to the total luminosity from galaxies fainter than our luminosity resolution limit as follows:

$$L_{b_j,\mathrm{tot}}=L_{b_j,\mathrm{faint}}+L_{b_j,\mathrm{bright}}$$
(10)
$$L_{b_j,\mathrm{bright}}=\underset{L_{b_j}>\mathrm{resolution}}{}L_{b_j,i}$$
(11)
$$L_{b_j,\mathrm{faint}}\frac{N_{bright}}{_{L_{b_j,\mathrm{res}}}^{\mathrm{}}\varphi _{b_j}(L)𝑑L}\times _0^{L_{b_j,\mathrm{res}}}L\varphi _{b_j}(L)𝑑L$$
(12)
Table 6 lists the results: the two numbers cited for the data are the minimum and maximum values found by Girardi et al. (2000) when varying their sample of galaxies, their maximum clustercentric radius, and their method to estimate luminosities. For Centaurus, Girardi et al. (2000) did not provide any data since the cluster is multipeaked. The observed ratios of Girardi et al. (2000) in the $`b_j`$-band scatter around a value of 250.

For our $`\mathrm{\Lambda }`$CDM model, we find values scattering around $`\mathrm{{\rm Y}}_1500`$ and $`\mathrm{{\rm Y}}_2600`$, and around $`\mathrm{{\rm Y}}_1600`$ and $`\mathrm{{\rm Y}}_2700`$ for the $`\tau `$CDM model, although with a noticeably larger scatter.

The higher value within the virial radius is a consequence of the radial variation of the morphology and colour of the galaxies. The typical discrepancy with observation appears to be factor of 2 and 2-3 for the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models respectively. This is somewhat larger than expected from the results presented above since our simulated cluster masses are not systematically in error, the brighter galaxies in our clusters are typically too faint by a magnitude or less, and our cluster luminosity functions $`\varphi _{b_j}`$ have faint end slopes (which are important for the luminosity correction discussed above) which are steeper than those of most observational determinations (Folkes et al. 1999; Blanton et al. 2001).

Further work, including careful analysis of the observational mass determination techniques is needed to understand how much of this discrepancy comes from incorrect galaxy modelling and how much from systematics introduced by the mass and luminosity measurement techniques. For example, the higher M/L ratios generically obtained in the $`\tau `$CDM model compared to $`\mathrm{\Lambda }`$CDM are not reflected in the Johnson $`B`$-band, where the figures are similar. At these wavelengths, our $`\mathrm{\Lambda }`$CDM simulation gives total M/L ratios of 536, 501, 512, 490 and 402 for Coma, Virgo, Hydra, Centaurus and Perseus respectively. Using modelling techniques very similar to ours, K99 found the $`B`$-band M/L ratios of clusters in $`\mathrm{\Lambda }`$CDM to be biased with respect to that of the Universe as a whole by an amount which is quite consistent with observation. In an independent analysis Benson et al. (2000b) obtained $`\mathrm{{\rm Y}}_\mathrm{B}400`$ for massive haloes in $`\mathrm{\Lambda }`$CDM. For comparison, Kent & Gunn (1982) measured $`\mathrm{{\rm Y}}_\mathrm{B}=360`$ for the Coma cluster, perhaps the best observed of local clusters.

It is also possible to use the galaxy populations in our simulations to study the radial distribution of galaxies within clusters. Detailed analyses of this kind were carried out by Diaferio et al. (1999, 2001) comparing their GIF simulations to CfA and CNOC clusters. At the resolution of our (and their) simulations, it is only possible to get useful results by stacking many clusters. We have done this with our data and find results which are very similar to that of the earlier work. In the mean the radial distribution of galaxies within clusters is similar to but slightly less concentrated than that of the dark matter. Blue and star-forming galaxies are much less concentrated towards the cluster centres than redder galaxies. The proportion of elliptical galaxies rises (and that of spirals drops) near the cluster centre. We do not repeat the plots of the earlier paper here since we have little to add. Further progress on these issues requires simulations of significantly higher resolution, such as those of S00.

## 5 Comparison catalogues

In this section, we recall briefly the features of the two local galaxy catalogues with which we compare our simulations. We use the far-infrared selected *PSCz* catalogue since it includes the IRAS 1.2 Jy survey which was used as the density constraint on our initial conditions. In addition it has near full-sky coverage and the highest available IR source density in the region we have modelled. A complication is that comparison with this survey requires us to model the FIR luminosity of simulated galaxies. At optical wavelengths, we have chosen the recently completed *UZC* as our reference catalogue. Its sky coverage is substantially smaller than that of the *PSCz* but it has a higher source density and we can compare the observed $`B`$-band luminosities directly to our simulations.

### 5.1 *PSCz* catalogue

The *PSCz* catalogue has been described by Saunders et al. (2000). It contains 15411 *IRAS* galaxies and covers some 84% of the sky, excluding regions of low galactic latitude and cirrus, and regions unobserved by the *IRAS* satellite. 1.2% of *IRAS* galaxies within the *PSCz* region and with $`b_j<19.5`$ have unknown redshift, and we discard them from the catalogue. We then consider the *PSCz* catalogue is limited both at $`b_j<19.5`$ and at $`f_{60\mu \mathrm{m}}>0.6`$ Jy. Also, we only consider galaxies with $`cz<8000\text{km s}\text{-1}\text{.}`$

To compare the results of our simulations with the *PSCz,* we need to estimate the IR luminosity of the galaxies at $`60\mu \mathrm{m}`$. Our model as described so far does not provide us with thisg information. We assume that the $`60\mu \mathrm{m}`$ luminosity of a star-forming galaxy has two components, one coming from star-forming regions and proportional to its star formation rate, the other representing the re-emission of obscured light from older stars and proportional to its $`I`$-band luminosity. For very low star formation rates we assume a negligible FIR luminosity. This threshold in star formation rate is set empirically. Formally, we write:

$$L_{\mathrm{FIR}}=a_{\mathrm{SFR}}\dot{M}_{}+\alpha _\mathrm{I}L_\mathrm{I}\mathrm{if}\dot{M}_{}>0.01\frac{M_{}}{t_{\mathrm{Hubble}}}$$
(13)
$$L_{\mathrm{FIR}}=0\mathrm{otherwise}$$
(14)
This is known to be at best a rough approximation for *IRAS* galaxies (Helou 1986), but it gives sufficiently accurate results for the purposes of this paper. A more detailed analysis would require consistent inclusion of the IR band in our SA model (see Devriendt & Guiderdoni 2000). We determine the proportionality coefficients $`a_{\mathrm{SFR}}`$ and $`\alpha _\mathrm{I}`$ for each cosmology by matching the IR luminosity function of simulated galaxies brighter than the resolution limit to the *IRAS* luminosity function given in Saunders et al. (1990). Because of resolution effects we are able to fit the observed *IRAS* luminosity function only above a limiting $`60\mu \mathrm{m}`$ luminosity of $`L_{60}=\mathrm{2.7\hspace{0.25em}10}^8\mathrm{L}_{}`$ and $`L_{60}=\mathrm{1.50\hspace{0.25em}10}^8\mathrm{L}_{}`$ in $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM respectively, which we define as our FIR resolution limits. Here, $`L_{60}`$ is expressed in units of the bolometric luminosity of the sun, and is consistent with the assumed galaxy FIR spectrum of Saunders et al. (1990). To avoid confusion we will in the following refer to our previous luminosity resolution limits as the “optical resolution limits” to distinguish them from the “FIR resolution limits”. Interestingly, the fitting values of $`a_{\mathrm{SFR}}`$ that we obtain: $`a_{\mathrm{SFR}}=\mathrm{1.68\hspace{0.25em}10}^9`$ and $`a_{\mathrm{SFR}}=\mathrm{6.4\hspace{0.25em}10}^8`$ in $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM respectively are close to the one estimated by Rowan-Robinson (2000) for bright IR galaxies. As a check, we compare the number of simulated and observed galaxies with $`f_{60\mu \mathrm{m}}>0.6`$ Jy, $`cz<8000\text{km s}\text{-1}`$, $`b_j<19.5`$ and brighter than both the optical and FIR resolution limits. The figures are reported in Table 7. Note that they are different for our two cosmologies because of the differing simulation resolution limits.

### 5.2 *UZC* catalogue

According to Falco et al. (1999), the *UZC* is $`98\%`$ complete for $`m_{\mathrm{Zwicky}}15.5`$ in the ranges $`20^\mathrm{h}\alpha _{1950}4^\mathrm{h}`$ and $`8^\mathrm{h}\alpha _{1950}17^\mathrm{h}`$, for declinations $`2^{}.5\delta _{1950}50^{}`$. This covers roughly one third of the sky. We again consider only galaxies with $`cz<8000\text{km s}\text{-1}\text{,}`$ and with absolute magnitudes above the optical resolution limit for our simulated galaxies.

To compare with the *UZC,* we need to cut the simulated galaxy catalogues in the appropriate $`B`$-band magnitude. However, the assumption of $`M_{Zwicky}M_B`$ is known to be rather approximate (Huchra 1976; Bothun & Cornell 1990), with a $`1\sigma `$ error of 0.3 mag. We adjust our threshold in simulated $`m_B`$ to obtain good agreement between the predicted and observed number of galaxies in the *UZC* region. This process is independently applied to both our models. It leads us to adopt $`M_{Zwicky}M_B=0.4`$ and $`+0.9`$ for the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models respectively. These shifts are larger than can plausibly be attributed to the difference in photometric systems, but allow us to compensate in part for the systematic offsets in luminosity function that are visible in Fig. 5. Again, the resulting number of galaxies in our real and mock *UZC* catalogues are reported in Table 7.

## 6 Smoothed density fields

The constrained density field used to set up the initial conditions for our two simulations was derived from the *IRAS* 1.2 Jy catalogue, smoothed with a Gaussian of dispersion 5 $`h^1`$ Mpc. As a result we expect the galaxy density fields in our simulation to resemble each other and to resemble the real galaxy field when smoothed on similar scales. In this section we will check how well this works out. In the following, we will denote by “optical galaxies” the sample of all simulated galaxies brighter than our $`B`$-band resolution limit. Similarly, we will call “FIR galaxies” all simulated galaxies brighter than both the optical and the FIR resolution limits. We will apply apparent magnitude limits to these samples only in section 6.4 where we compare with the observed *PSCz* sample.

### 6.1 Smoothed maps

As a qualitative example, Figs. 9 and 10 show the overdensities of optical galaxies in the supergalactic plane, using Gaussian kernels with 5 and 10 $`h^1`$ Mpc smoothing length respectively (the symbol $`R`$ in Equation 15 below). As expected, there is already good agreement between the two cosmologies on the 5 $`h^1`$ Mpc smoothing scale, and the match is even better at 10 $`h^1`$ Mpc. The three major overdensity peaks are Coma at the top of the plot, the Great Attractor at the center–left, and the Pisces-Perseus complex at the center–right. The Virgo cluster/Local Supercluster complex is also evident just above the centre of the plots. These maps can be compared with galaxy density maps reconstructed directly from various *IRAS* surveys (e.g. Branchini et al. 1999). The resemblance is quite close and is encouraging because of the many steps between using the observed density fields to constrain the initial dark matter distributions in our models and constructing simulated galaxy density fields from the galaxies which are formed during their evolution.

To obtain more quantitative results, we now compare various smoothed density fields point-by-point. To do this we assign the dark matter particles and the galaxies to a regular grid with a spacing of 10 $`h^1`$ Mpc using a CIC scheme. We keep only those grid points within $`cz<8000\text{km s}\text{-1}\text{,}`$ resulting in some $`2100`$ nodes. We further smooth this density field using a Gaussian kernel of the form

$$W(r)=\frac{1}{(2\pi R^2)^{3/2}}\mathrm{exp}\frac{r^2}{2R^2}$$
(15)
We take $`R=10`$$`h^1`$ Mpc, corresponding to twice the smoothing used when generating the initial conditions. After this smoothing, we correct for artificial effects introduced by the spherical boundary of our galaxy distribution. This is done by multiplying the smoothed field by the ratio of fields obtained by applying similar smoothing procedures to two Poisson distributed catalogues, one uniform within a box of side $`220h^1\text{ Mpc}`$ centred on the simulation sphere, the second uniform within the region where our data lie ($`cz<8000\text{km s}\text{-1}\text{,}`$ together with the *PSCz* mask when appropriate). We also evaluate the volume of the intersection between the smoothing kernel and our selection region using this simple Monte-Carlo scheme and we discard nodes where this volume is less than half of the volume of the smoothing kernel.

Note that we have chosen the spacing of our grid so that the smoothed density fields are not heavily oversampled, while retaining a reasonable number of data points. For consistency, we define the mean density of each smoothed field as the straight average of its values at our final set of nodes. We can then define the smoothed overdensity at each node through

$$\delta _{\mathrm{s},10}=\frac{\rho _{\mathrm{s},10}1}{\overline{\rho }_{\mathrm{s},10}}$$
(16)
We now discuss scatter plots which show point-by-point comparisons of various of these density fields. We begin by focussing on different components within a given simulation. We then compare each component between our two simulations. Finally we compare our simulated FIR galaxy density fields with the real *PSCz* data.

### 6.2 Optical galaxies vs. FIR galaxies vs. mass

In Fig. 11, we compare the overdensity distributions of optical galaxies, of FIR galaxies and of mass in each of our simulations.

The overdensities of optical galaxies and of dark matter are remarkably tightly correlated in both simulations. In each case the galaxies are almost unbiased with respect to the underlying DM except in the highest density regions ($`\delta _{\mathrm{s},10}1\mathrm{to}\mathrm{\hspace{0.33em}3}`$), where the optical $`\mathrm{\Lambda }`$CDM galaxies are significantly antibiased.

Biases are more evident and the scatter is larger when the FIR overdensity fields are compared with those for the dark matter. For $`\mathrm{\Lambda }`$CDM the trends are similar to but stronger than those for the optical galaxies. At DM overdensities of order $`\delta _{\mathrm{s},10}3`$, galaxies are antibiased by a factor of almost 3. This is expected given our model for the FIR luminosity: cluster galaxies have low star formation rates and so are assigned small or zero $`L_{\mathrm{FIR}}`$ (consistent with observations of most E/S0 galaxies) and so are almost all excluded from our simulated FIR catalogues. This substantially reduces the number of galaxies counted in the densest regions. On the other hand, in mean density and underdense regions star–forming spirals predominate, and these are included in both the optical and the FIR samples. As a result the scatter plots for the two populations are quite similar in this density range.

A similar but weaker pattern is observed for the $`\tau `$CDM FIR galaxies and for the same reasons. The weakening in this case reflects the fact that the present day star formation rate in the $`\tau `$CDM cosmology is five times that in the $`\mathrm{\Lambda }`$CDM cosmology, and a significant number of galaxies in dense regions have fallen in recently enough for their ongoing star formation to take them above our $`60\mu \mathrm{m}`$ resolution limit.

The larger scatter in the FIR plots may well be a real effect, reflecting, perhaps, the stochasticity introduced by the fact that FIR luminosity is linked to star formation rate rather than stellar mass. The resulting underweighting of all group and cluster environments introduces scatter as well as bias into the relation between smoothed mass and smoothed light. We checked that the lower scatter in the optical case is not a result of the larger number of galaxies in our simulated optical samples. Randomly picking one third of the optical galaxies and using this subset to construct the overdensity field does not lead to a noticeable increase in the scatter in the optical plot.

As a consistency check, we plot in the two bottom panels of Fig. 11 the overdensities of FIR galaxies against those of optical galaxies. The effect of our FIR modelling is again very clear: in both models, but particularly in the $`\mathrm{\Lambda }`$CDM model, significantly overdense regions in terms of optical galaxies become moderately overdense in terms of FIR galaxies.

### 6.3 $`\mathrm{\Lambda }`$CDM vs. $`\tau `$CDM

In this section, we compare the distributions of dark matter, of optical and of FIR galaxies between our two cosmologies. Recall that although the initial conditions of the two simulations are nearly identical when smoothed on scales of 5 $`h^1`$ Mpc (apart from a difference in fluctuation amplitude) they differ on smaller scales. As a result, the build-up of galaxies is almost uncorrelated between them.

The upper left plot of Fig. 12 compares dark matter overdensities in the two models. The shape is as expected from the initial set-up of the simulations. Initial amplitudes were chosen so that the abundance of objects of rich cluster mass would agree in the two cases. This requires the amplitude of linear overdensity fluctuations to be larger in $`\mathrm{\Lambda }`$CDM than in $`\tau `$CDM. The difference in amplitude carries over to the final time and shows up as a clear “antibias” of the $`\tau `$CDM density relative to the $`\mathrm{\Lambda }`$CDM field. On the other hand, the correlation between the two overdensity fields is strong, showing that the differences in the initial conditions (and so in the evolved mass fields) on small scales have little influence on larger scale fluctuations for the 10 $`h^1`$ Mpc smoothing used to make these plots.

The upper right plot shows that the galaxy formation recipes we have implemented, which are tuned to reproduce the observed Tully–Fisher relation and, to a lesser extent the observed luminosity function, are only partially successful in compensating for this difference in dark matter clustering. The distribution of optical galaxies in $`\tau `$CDM is still significantly antibiased relative to that in $`\mathrm{\Lambda }`$CDM, although the effect is weaker than for the dark matter.

Somewhat surprisingly, this difference in clustering amplitude is almost absent in the simulated FIR galaxy distributions. The lower left plot in Fig. 12 shows that the overdensity fields of the two simulations are very similar for this population. The higher star-formation rates of non-cluster galaxies in the $`\tau `$CDM model clearly boost the relative amplitude of the density fluctuations in FIR population sufficiently to compensate for the weaker clustering of galaxies by optical luminosity (or stellar mass).

In the lower right plot, we compare the *redshift–space* distributions for this same FIR galaxy population. As expected, the scatter is greater in redshift space and the overall amplitude of the density fluctuations is increased. In addition, an extended plume appears at high densities. These points correspond to the environments of massive clusters, in particular the Great Attractor and Pisces-Perseus regions. It is interesting that peculiar velocities appear to enhance the density contrast of these structures more strongly in $`\mathrm{\Lambda }`$CDM than in $`\tau `$CDM.

From this series of four plots we conclude that the point-by-point agreement between our two simulations is quite good for all three components, at least for the rather large smoothing employed here. As expected, the mass is more weakly clustered in $`\tau `$CDM, and this effect carries over to the optical galaxy population. The smoothed FIR galaxy populations in the two simulations are very similar.

### 6.4 Simulated *PSCz* galaxies vs. observations

We now compare the distribution of the simulated FIR galaxies directly to the observed distribution in the *PSCz* catalogue. This requires a slightly different FIR sample than that used in the last section. In addition to the redshift and simulation resolution limits already enforced, we need to apply the sky mask of the observed catalogue and its flux limits at 0.6 Jy and $`b_j<19.5`$. The effects of these additional selection criteria are consistently taken into account when constructing overdensity fields both for the observations and for the simulated data.

The left and right plots of Fig. 13 compare simulated and observed overdensity fields for the $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM cosmologies respectively. These fields are, of necessity, constructed in redshift–space. The agreement is reasonably good, with greater scatter in the $`\mathrm{\Lambda }`$CDM case. In neither cosmology is there any obvious bias between simulation and observation. This is a nice confirmation that our SA and FIR luminosity modelling schemes produce an FIR galaxy population with a large-scale distribution which quite closely resembles that used to define our initial conditions.

## 7 Correlation functions

In this section we use correlation techniques to extend both our clustering measurements and our comparison with observations to smaller scales. We begin by using auto- and cross-correlations of the various simulated populations to quantify small-scale bias. We then cross-correlate observed and mock catalogues to further explore how well we have reproduced our own cosmic neighborhood.

### 7.1 Autocorrelations

In Fig. 14 we plot real–space autocorrelation functions for the dark matter and for galaxies brighter than our $`B`$-band resolution limit (“optical galaxies”) for each of our simulations. For reference, we also plot the dark matter autocorrelations for the *GIF* simulations analysed by K99 and the autocorrelation function of real optical galaxies as inferred from inversion of the angular correlation data for the *APM* survey (Baugh 1996). We compute our correlation functions up to a scale of 10 $`h^1`$ Mpc. Note that the analysis in Jenkins et al. (1998) shows that the *GIF* results are close to the ensemble-averaged expectations for the two cosmologies.

A comparison of the dark matter correlations between our simulations and the *GIF* simulations shows similar behaviour in the two cosmologies. There is reasonably good agreement on small scales, but our simulations have more power than the *GIF* models on scales of a few Mpc. This is clearly a reflection of the particular way in which we have constrained the large-scale density field in our initial conditions using the observed distribution of 1.2 Jy galaxies.

On small scales the optical galaxy distribution is quite strongly antibiased relative to the mass in both our cosmologies, while on scales of a few Mpc the difference is much smaller. Indeed, in $`\tau `$CDM the correlations of dark matter and optical galaxies are almost equal beyond 1.5 $`h^1`$ Mpc. For neither case are the autocorrelations of the optical galaxies close either to a power law or to the *APM* data. It is unclear whether this is a problem, since correlations calculated for the relatively local region we are modelling should not necessarily reproduce those found for much larger “representative” regions. We will see below that our $`\mathrm{\Lambda }`$CDM model does seem to reproduce the correlation statistics of local galaxies quite well, while our $`\tau `$CDM model does not. On scales of 5 to 10 $`h^1`$ Mpc the correlation amplitude for optical galaxies is quite close to the *APM* values in both our simulations.

### 7.2 Cross-correlations

To make a more quantitative comparison between our simulations and the observed *PSCz* and *UZC* catalogues we have made mock *PSCz* and *UZC* catalogues which reproduce in detail the sky coverage and the apparent luminosity limits of the observational data. We also limit both the mock catalogues and the real catalogues to galaxies with absolute luminosities brighter than the relevant resolution limits of our simulated catalogues (see above). We can then compare real and simulated distributions in detail using auto- and cross-correlations. We have applied the same sky masks and depth selection functions to our simulated mass distributions to produce mass catalogues which can be compared directly with the observed galaxy distributions using the same techniques.

To compute auto- and cross-correlations we use the estimator suggested by Hamilton (1993):

$$\xi _{12}(r)=\frac{D_1D_2R_1R_2}{D_1R_2D_2R_1}1,$$
(17)
where $`DD`$, $`RR`$ and $`DR`$ refer to the number of data-data, random-random, and data-random pairs, and the subscripts refer to the two catalogues. We compute these correlations from 0.7 $`h^1`$ Mpc up to 15 $`h^1`$ Mpc. The random catalogues used here are generated using the same angular mask as the corresponding galaxy (or mass) catalogue and a selection function in depth derived from the relevant luminosity function (that of Saunders et al. 1990 for the *PSCz,* mock *PSCz* and associated mass catalogues; that of Marzke et al. 1994 for the *UZC* and its associated mass catalogue; those of 3.5.2 for the mock *UZC* catalogues). We ensure that each random catalogue contains at least ten times as many points as its corresponding “data” catalogue. This ensures that uncertainties in pair counts are dominated by the number of available galaxies (or mass particles).

#### 7.2.1 Cross–correlations with FIR galaxies

In Fig. 15, we give *redshift–space* auto– and cross–correlation functions for our observed *PSCz,* mock *PSCz* and *PSCz* mass catalogues.

A surprising result in this figure is that the autocorrelation of the mock *PSCz* catalogue agrees almost perfectly with that of the real data for $`\mathrm{\Lambda }`$CDM and is quite close to it for separations above 1.5 $`h^1`$ Mpc in the $`\tau `$CDM case also. In both cases the model autocorrelation is close to a power-law of index $`1.28`$, the value found by Fisher et al. (1994) for the 1.2 Jy sample as a whole. Over the range we plot, the bias of the simulated *PSCz* galaxies relative to the mass is almost constant; the $`\mathrm{\Lambda }`$CDM galaxies are substantially antibiased while the $`\tau `$CDM galaxies are almost unbiased.

As expected, cross-correlations of the observed *PSCz* galaxies with either our mock *PSCz* galaxies or the simulated mass distribution show a different behaviour. On scales larger than about 5 $`h^1`$ Mpc the cross-correlation between real and mock galaxies is as strong as (and effectively equal to) the autocorrelation of either population. This is a striking confirmation of the effectiveness of our constraint and galaxy formation procedures. On scales below 5 $`h^1`$ Mpc this cross-correlation flattens out to a constant value, reflecting the fact that the small-scale structure in our simulations is unrelated to that in the real Universe.

The cross-correlation between the observed *PSCz* galaxies and the simulated mass distribution has a similar shape, and indeed, for $`\tau `$CDM (where the mock *PSCz* galaxies are almost unbiased) it is almost identical to the cross-correlation with the mock galaxy catalogue. For $`\mathrm{\Lambda }`$CDM the cross-correlation with mass is stronger than that with the mock galaxies. On large scale this cross-correlation is approximately the geometric mean of the autocorrelations of the two populations. This would be expected for a pure linear bias model but may be surprising in view of the scatter and t¡he nonlinearity evident in Figs. 11 and 13.

#### 7.2.2 Cross–correlations with optical galaxies

Fig. 16 shows auto- and cross-correlation functions in the same format as Fig. 15 but for observed and mock galaxy catalogues and for mass catalogues with the sky mask and depth distribution of the *UZC* catalogue.

It is interesting to compare the autocorrelation functions of this plot with the real-space autocorrelations in Fig. 14. As was also the case for the mock FIR galaxies, the redshift-space autocorrelation functions of the mock optical galaxies are parallel to those of the dark matter; for $`\mathrm{\Lambda }`$CDM the galaxies are antibiased, although less strongly so than the FIR galaxies, while for $`\tau `$CDM the galaxies show a slight positive bias. In both cases the dark matter and galaxy curves are much more nearly parallel than are the corresponding real-space functions in Fig. 14. Note also the substantial differences between the mass autocorrelations plotted in Fig. 16 and those shown for the same simulations in Fig. 15. These reflect the different sky coverage and depth distribution of our *PSCz* and *UZC* catalogues and emphasise that neither should be considered a fair sample of the Universe as a whole.

When we compare our mock catalogue autocorrelations with those of the real *UZC* galaxies we find excellent agreement below 5 $`h^1`$ Mpc for $`\mathrm{\Lambda }`$CDM with a slight underprediction of the observed amplitude on larger scales. For $`\tau `$CDM the predicted autocorrelations are low on all scales, with the difference in $`\xi (s)`$ ranging from 30% on large scales to a factor of 2 around 1 $`h^1`$ Mpc. Comparing with Fig. 14 and Fig. 15, we see that while the autocorrelations of observed optical galaxies in our *UZC* region are significantly weaker than those of observed FIR galaxies in our *PSCz* region, they nevertheless substantially exceed those measured for galaxies in the *APM* survey (even allowing for the difference between real and redshift space). This again emphasises that the regions of space we are analysing are quite small and should not be thought representative.

The cross–correlations of the observed *UZC* galaxies with the mass and mock *UZC* galaxy distributions in our simulations are almost flat below $``$ 5 $`h^1`$ Mpc, especially in the $`\tau `$CDM model. On larger scales they do approach the autocorrelation amplitude of the observed galaxies but the convergence is less compelling than was the case when we compared our mock FIR galaxies with the *PSCz* data. This difference may in part reflect the smaller volume covered by our *UZC* catalogues. It may also be due to the fact that we used an FIR galaxy catalogue when setting up our constrained initial conditions.

To summarise these sections on correlation functions, our initial condition generation and galaxy evolution procedures have allowed us to build physically consistent $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM models which reproduce well both the individual structures and the statistical clustering properties of star-forming galaxies in our local neighborhood. For $`\mathrm{\Lambda }`$CDM the statistical properties of optically selected galaxies from the *UZC* catalogue are also well reproduced, but this is not the case in our $`\tau `$CDM simulation where the optical galaxies are significantly more weakly clustered.

## 8 Mock Catalogues of Peculiar Velocities

As an example of how our simulations can be used to test and calibrate the methods of cosmic-flow analysis, we have generated Monte-Carlo mock samples mimicking real catalogues of peculiar velocities. In particular, we show here mock catalogues corresponding to the *Mark III* Catalogue of Peculiar Velocities (Willick et al. 1995, 1996, 1997). These were made with an improved version of the procedure described in K96, which provided the basis for earlier testing of methods. The major advance allowed by our new simulations comes in part from their ability to follow the nonlinear evolution of the mass distribution on sub-Mpc scales, but more importantly from the fact that they simulate directly the formation of the galaxies themselves. In K96, “galaxies” were inserted by applying a statistical biasing scheme to the smoothed mass density field, and their magnitudes were drawn at random from an assumed luminosity function. Here, galaxies form in a physically consistent way through the condensation of gas at the centres of dark haloes and their morphologies reflect their actual merging histories. This more realistic treatment should better capture possible correlations between galaxy properties and the underlying density and velocity fields. Our simulations, based on constrained realizations, then allow a full reproduction of the correlation between systematic errors and the actual signal of density and velocity fields.

For the purpose of this section we distinguished between “spiral” and “elliptical” galaxies by defining ellipticals as galaxies with total stellar mass $`>4.5\times 10^{10}\mathrm{M}_{}`$, and bulge-to-total $`V`$-band luminosity ratio $`>0.4`$. Small adjustments were made to the semi-analytic absolute magnitudes in order to fine-tune their match with a Schechter luminosity function. We then assigned to each galaxy an “observed” linewidth based on an assumed Tully–Fisher relation and scatter as in K96 – in practice these linewidths were always close to those assigned by the algorithms described earlier which are based on the actual circular velocities of the haloes. To simulate the selection based on magnitude in each dataset of the *Mark III* catalogue, we assumed Galactic extinction as a function of Galactic latitude as in K96, and we slightly smeared the magnitude limits to take into account scatter in the relation between the magnitudes used for selection and those appearing in the Tully–Fisher relation. Fig. 17 demonstrates our success in matching the final distribution of redshift and apparent magnitude in our mock catalogues with that of the real data from the largest single dataset within the *Mark III* catalogue.

Each mock dataset was then diluted at random (trying to mimic selection by other independent properties such as inclination or lack of a strong bar) to match the number of galaxies in the real catalogues. This random sampling, along with the random distance errors (introduced by the TF scatter) was repeated 10 times to generate 10 mock catalogues. The mock data were then grouped, and statistically corrected for Malmquist bias in exactly the same way as the real data (details are in the *Mark III* papers and are summarised in the POTENT analysis, Dekel et al. 1999).

Mock and real maps of radial peculiar velocities in a slice about the Supergalactic plane arEe shown in Figs. 18 to 20, illustrating the degree of agreement between simulated and observed peculiar velocity fields and between the real-space density fields reconstructed from them. For the simulations, the reconstructed density fields can also be compared with the real density fields which are plotted in Fig. 10. Overall the level of agreement is impressive. Systematic differences are visible in certain areas and presumably reflect noise effects together with systematics arising because the smoothed *IRAS* 1.2 Jy density field which constrained our initial conditions does not perfectly trace the true mass density field. Differences may also arise because the observed peculiar velocities are generated in part by the mass distribution outside the region where we constrained our initial conditions.

Similar mock catalogues provided the basis for a revised likelihood analysis of peculiar velocities by Silberman et al. (2001). The proper incorporation of nonlinear effects revealed a systematic overestimate of $`\mathrm{\Omega }_\mathrm{m}`$ in earlier linear analyses. It then allowed the development of unbiased nonlinear methods, which brought the best estimates from the *Mark III* and *SFI* data to $`\mathrm{\Omega }_\mathrm{m}=0.35\pm 0.1`$.

We publically release the mock catalogues for *Mark III.* They are available upon request from AE or AD (eldar@phys.huji.ac.il). We can also provide similar mock catalogues tailor-made for other peculiar-velocity catalogues upon request.

## 9 Summary and Conclusions

The goal of the current study has been to carry out physically consistent simulations of the formation and evolution of the local galaxy population within the currently dominant CDM structure formation paradigm. We have attempted to reproduce not only the statistics of galaxy properties (luminosity functions, colour and morphology distributions, Tully–Fisher relations…) and of galaxy clustering (correlation functions, cluster luminosities…) but also the actual spatial distribution of clusters and superclusters within 80 $`h^1`$ Mpc of the Milky Way. This is an ambitious undertaking and clearly we have been more successful in some aspects of it than in others.

On small scales the statistics of the dark matter distributions in our $`\mathrm{\Lambda }`$CDM and $`\tau `$CDM simulations are a very good match to those expected. On scales beyond a Mpc or two, discrepancies are apparent which reflect both the relatively small volume simulated and the constraints we imposed on the initial large-scale fluctuations. For both simulations the morphology of the dark matter structure on large scales is a good match to that of the observed 1.2 Jy *IRAS* sample which we used as a constraint. In addition, the most massive clusters agree quite well with observed clusters both in position and in mass.

The distributions of properties of individual galaxies are in fair agreement with those observed. The Tully–Fisher relations we find for our simulated spiral galaxies are, by construction, in excellent agreement with observation. Our distributions of morphology and colour also resemble those observed and show the correct dependence on galaxy environment. Our luminosity functions, on the other hand, show significant discrepancies with those emerging from the 2dF and SDSS surveys. This problem is substantially reduced from K99, the paper on which our modelling is based, because of a number of minor improvements we have introduced, but it is still worse than found by several other groups carrying out similar modelling (e.g. Cole et al. 2000; Somerville et al. 2001, S00). We conclude that this disagreement probably reflects an inadequacy in our galaxy modelling, rather than any fundamental problem with the paradigm.

The statistics of galaxy clustering in our $`\mathrm{\Lambda }`$CDM simulation is in good agreement with that seen both in the far-infrared selected *PSCz* survey and in the optically selected *UZC* survey, at least as inferred from the autocorrelation functions of the galaxy populations. For our $`\tau `$CDM simulation the agreement with the *PSCz* survey is still quite good, but the amplitude of optical galaxy correlations is about 40% lower than in the *UZC* survey.

This latter problem may be related to the fact that the galaxy populations of our simulated galaxy clusters are in poor agreement with observation. Our brightest cluster galaxies are mostly about a magnitude too bright, while other bright cluster galaxies are about a magnitude too faint. Overall our clusters have mass-to-light ratios at least a factor of 2 larger than observed. This difficulty is apparently a consequence of the details of our galaxy formation modelling. Diaferio et al. (2001), using the simulations of K99, find a similar problem with the brightest cluster member but find their other cluster galaxies to be slightly too bright, while S00 find cluster luminosity functions which agree well in shape with observation, but mass-to-light ratios which are somewhat high. S00 show the improvement in luminosity function shape and the reduction in central galaxy luminosity to be a consequence of the greatly improved resolution of their simulation. This allows galaxy merging to be followed explicitly rather than inserted “by hand” using a dynamical friction model.

Point-by-point comparison of the smoothed density fields of mass and of optical- and FIR-selected galaxies in our two simulations, as well as comparison with the observed density field of *PSCz* galaxies, illustrates both the nature of the biases which relate these various components and the degree to which our techniques have been successful in reproducing the observed spatial distribution of star-forming galaxies. Our simulations produce a remarkably tight relation between the densities of optically selected galaxies and of dark matter. Biases are present but are relatively weak, too weak in the $`\tau `$CDM case to be consistent with observation, as already noted above. The fact that star-forming galaxies avoid rich groups and clusters results both in more scatter and in a strong bias at high density in the relation between FIR galaxy density and either mass or optical galaxy density. Comparing density fields between simulations or with the observations shows the strong relations introduced by our initial condition constraints as well as biases caused by differing amounts of dynamical evolution. The excellent agreement between our simulated FIR galaxy distributions and the *PSCz* data is best illustrated by the cross-correlation between the two which is within 20% of the autocorrelation of the *PSCz* data on scales above 5 $`h^1`$ Mpc.

Finally our mock catalogues of peculiar velocity data agree well with the real *Mark III* data. The reconstructed density fields from the mock catalogues show a close correspondance to each other, to the true density fields in the simulations, and to the reconstructed density field from the real data. This is both a reassuring demonstration that our techniques have achieved their primary goal, and a signpost to the way our simulations may be used to calibrate quantitative scientific analysis of the local galaxy distribution. We will present more detailed applications in future work, and we release a variety of galaxy, halo and dark matter catalogues in order that others can use them also. The data are available at the URL: http://www.mpa-garching.mpg.de/NumCos/CR/index.html

## Acknowledgements

The simulations presented in this paper were carried out on the T3E supercomputer at the Computing Center of the Max-Planck-Society in Garching, Germany. This work was supported in part by a grant from the German-Israel Science Foundation.