text:
# The evolution of the galaxy distribution

## 1 Introduction

Studies of the large-scale distribution of galaxies have traditionally focussed on problems such as testing hypotheses for the identity of the dark matter, the nature of the initial density perturbations and the mechanism of structure growth. Properties of the observed large-scale structure are also often used to estimate the values of fundamental cosmological parameters. Although none of these issues can be regarded as settled, there is now a growing consensus that cold dark matter (CDM) is the most likely candidate for the dark matter, that cosmic structure grew by the gravitational amplification of random-phase initial density fluctuations of inflationary origin, and that the fundamental cosmological parameters have the following values: density parameter, $`\mathrm{\Omega }_00.3`$, cosmological constant term, $`\mathrm{\Lambda }_00.7`$ and Hubble constant (in units of 100 km s-1 Mpc-1), $`h0.7`$.

Cosmological constraints reflect only one aspect of the information encoded in the pattern of galaxy clustering. Another, equally interesting aspect, concerns the processes responsible for the formation and evolution of galaxies. To extract this kind of information requires very extensive datasets and these are only now becoming available in the form of a new generation of galaxy surveys, like 2dF (?), Sloan (?) and 2MASS (?). The expectation is that these new datasets will provide, in addition to further cosmological constraints, some understanding of how the physics of galaxy formation manifests itself in the clustering of galaxies as a function of internal properties such as morphology, luminosity, colour, or star formation rate. Not only is this important for testing models of galaxy formation, but it is also required for extracting accurate cosmological information from the new surveys. Although it seems plausible that on very large scales the galaxy distribution traces the underlying mass in a simple way (??), complex biases are predicted to be present on small and intermediate scales (??).

To extract useful information from observational data of the quality and size of the new surveys, it is necessary to have detailed theoretical predictions. In this area too there have been significant advances in recent years, largely through the development of increasingly realistic ab initio calculations of galaxy formation and evolution. Two strategies have been developed for this purpose. In the first one, cosmological N-body/gasdynamics simulations are used to follow the coupled evolution of dark matter and gas, in particular, the cooling of gas in galactic dark matter halos (e.g. ??????). A phenomenological model is employed to decide when and where stars and galaxies form from this cooled gas and to include the associated feedback effects. In the second strategy, only the evolution of the dark matter component is simulated directly, or the assembly history of halos is obtained with a Monte-Carlo method, and the behaviour of the gas is calculated by solving a simple, analytical, spherically symmetric cooling-flow model. As in the direct simulation approach, star formation and feedback are included in a phenomenological way.

The two strategies offer different advantages. Direct simulations solve the evolution equations for gravitationally coupled dark matter and dissipative gas without imposing any restrictions on geometry. However, limited resolution restricts the range of length and mass scales that can be studied, and the expense of large simulations makes it impractical to carry out extensive parameter space explorations. Because of its simplified treatment of gas dynamics, semi-analytic modelling can follow an essentially unlimited range of length and mass scales, and is sufficiently flexible that the effects of varying assumptions and parameter values can be readily explored. Additional processes that cannot currently be easily investigated at the resolution available in direct simulations, such as those determining galaxy morphology, or the effects of dust obscuration, can be readily incorporated into the semi-analytic models by straightforward extensions to the phenomenological model of star formation and feedback. The numerical resolution and physical content of a typical N-body/gasdynamic simulation can be mimicked in a semi-analytic model and ? have shown that, at least in the case where only the simplest gas physics are modelled, the two techniques give reassuringly similar statistical results.

In this paper we combine large N-body simulations with the semi-analytic model of ? to investigate certain properties of the galaxy distribution that are relevant to the new generation of galaxy redshift surveys. We begin by displaying images that illustrate the evolution of the galaxy population in a representative volume of a simulated CDM universe. These images furnish some qualitative understanding of the mechanisms responsible for establishing the spatial distribution of galaxies of different kinds. We then focus specifically on the dependence of the two-point correlation function on galaxy luminosity, the morphology-density relation, and the topology of the galaxy distribution as measured by the genus, and on the evolution of these properties with redshift. The dependence of clustering on luminosity and colour have previously been considered, using similar techniques, by ?? and ?. The first and last of these studies found a weak increase in the correlation function with luminosity on large scales, but the second failed to detect any effect. These papers used simulations of relatively small volumes and so were unable to investigate the clustering of the brightest galaxies, for which luminosity-dependent effects are expected to be strongest. The new generation of redshift surveys will include large samples of very bright galaxies and may well be able to measure this kind of effect. In this paper we extend earlier work and investigate clustering at the bright end of the galaxy luminosity function. Closely related to the dependence of clustering strength on luminosity is the morphology-density relation which we also quantify in our simulations, both at the present day and at $`z=1`$. Finally, we provide the first determination of the genus curves predicted for galaxies in a CDM model; previous simulations had only been able to address the genus curves of the dark matter distribution. Although our model predictions are directed at the new surveys, we carry out limited comparisons with available observational data.

The remainder of this paper is organised as follows. In §2, we describe our simulation and modelling techniques. In §3 we present colour images of the evolution of a slice of our simulated volume (these images are publically available at `http://www.astro.caltech.edu/~abenson/Mocks/mocks.html`). In §4 we present quantitative estimates of clustering, namely the correlation length-luminosity and morphology-density relations and the genus statistic, and compare our results to observations. Finally, we present our conclusions in §5.

## 2 Method

The need for realistic modelling of galaxy formation as a prerequisite for deriving reliable clustering predictions has been emphasised by ? who showed, for example, that the form of the two-point correlation function on small scales is strongly influenced by the physical processes governing galaxy formation. Such processes are readily taken into account when the techniques of semi-analytic modelling are grafted into N-body simulations of the dark matter (?????). Monte-Carlo implementations of the semi-analytic model can also be used directly for clustering studies, without N-body simulations, but they only work well on scales larger than the Lagrangian radii of the dark matter halos which host galaxies, in practice on scales in excess of a few Mpc (although see ?), for which the bias can be calculated using the analytic formula of ?; see for example ?. Modelling the visible properties of galaxies explicitly allows simulated samples to be selected according to criteria closely patterned on observational selection procedures (e.g. by magnitude, colour, morphology, etc), thus allowing rigorous comparisons with observations to be made.

In this paper, we use the techniques introduced by ? and extended by ? to graft our semi-analytic model of galaxy formation onto N-body simulations. Full details of our semi-analytic model and the extensions required to study galaxy clustering are given in ? and ? respectively. Briefly, dark matter halos are identified in the simulation at the redshift of interest using the friends-of-friends algorithm with the standard linking length of 0.2 (?). After cleaning the halo catalogue in the manner described below, the mass of each halo is input into the semi-analytic model. By means of simple, physically-motivated prescriptions, described in detail in ?, the model calculates the amount of gas that cools in a virialized halo of that mass to make a galaxy, as well as the star formation rate, the reheating of left-over gas by stellar winds and supernovae, and the chemical evolution of the gas and stars. Galaxies are allowed to merge within common dark matter halos, producing elliptical galaxies and bulges from stellar disks. The spectrophotometric evolution of the galaxies is calculated using a standard stellar population synthesis model (??). Extinction by dust is included using the models of ?. The most massive galaxy in each halo is identified with the central galaxy and placed at the centre of mass (and given the peculiar velocity of the centre of mass). Other galaxies (satellites) are assigned the position and peculiar velocity of a randomly chosen dark matter particle within the halo. In this way, satellite galaxies always trace the dark matter within a given halo. Modelling of this kind has been successfully applied to study a large variety of properties of the galaxy population (e.g ????????.)

We consider only dark matter halos in the simulation containing ten or more particles. Since we are interested in galaxies of all luminosities, including the faint ones that occupy halos with masses close to the ten-particle limit, it is important to check that small halos are actually bound objects. For this, we calculate the total energy of each halo by summing the kinetic energy (measured relative to the centre of mass of the halo) and the gravitational energy (due to the interaction between all of the particles in the halo). If a halo is found to have positive energy (and so to be unbound), we remove the least bound particle and recompute the total energy. This process is repeated until either the energy becomes negative (in which case we now have a bound halo with a lower mass than the original) or there are fewer than ten particles left in the halo (in which case we discard it). In this way, we construct a new halo catalogue containing only bound objects.

Typically, approximately 10% of the halos from the original catalogue fail the binding energy test and are excluded. A slightly smaller fraction have particles removed but remain in the catalogue. Most of the excluded halos come from the low mass end of the distribution, with the excluded fraction dropping rapidly as the halo mass increases. We find many examples of halos contaminated by interlopers for which removal of a small number of the least bound particles results in a bound object. (Note that since the binding energy test only affects halos near the resolution limit of the N-body simulation, it does not alter any of the results of ? who considered only bright galaxies that form in halos well above the resolution limit.)

We adopt similar values for the parameters of the semi-analytic model as did ?, except that the parameters describing the normalisation and shape of the power spectrum ($`\sigma _8`$ and $`\mathrm{\Gamma }`$, see below) are set to the actual values used in the simulation. The parameters of the semi-analytic model are slightly different from those used in ?, but the differences in the predictions are negligible. Furthermore, since the dark matter halo mass function in the simulations differs somewhat from the Press-Schechter form assumed by ?, we find that the model works better if the value of $`\mathrm{{\rm Y}}`$ (the ratio of total to visible stellar mass, which depends on the fraction of brown dwarfs) is set equal to 1 rather than 1.4 as in ?. The recycled fraction in the calculation of the chemical enrichment is modified accordingly (?). ? showed that in a $`\mathrm{\Lambda }`$CDM cosmology, this model produces a real-space two-point galaxy correlation function which is remarkably similar to that measured in the APM survey by ?, in contrast to the $`\tau `$CDM cosmology which fails to match the observed two-point correlation function on all scales. Furthermore, ? showed that clustering predictions are robust to changes in the semi-analytic parameters provided that the model matches the bright end of the local galaxy luminosity function. The evolution of the galaxy correlation function with redshift is in good agreement with SPH simulations of galaxy formation (??). In this paper, we consider only the $`\mathrm{\Lambda }`$CDM model. We have checked that the parameters we have adopted do produce a correlation function of $`L_{}`$ galaxies identical to that of ?.

We use two different N-body simulations. The first one is the “GIF” $`\mathrm{\Lambda }`$CDM simulation, a full description of which may be found in ? and ?. This is a 17 million dark matter particle simulation in a cubic volume of side $`141.3h^1`$Mpc, with cosmological parameters $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Lambda }_0=0.7`$, $`h=0.7`$, $`\mathrm{\Gamma }=0.21`$ and $`\sigma _8=0.9`$ (where $`\mathrm{\Gamma }`$ is the power spectrum shape parameter and $`\sigma _8`$ is the linearly extrapolated rms mass fluctuation in a sphere of radius $`8h^1`$ Mpc). The mass of the smallest resolved halo in this simulation is $`1.4\times 10^{11}h^1M_{}`$. We also analyse the “$`512^3`$” simulation described by ? and ?, which has identical cosmological parameters to the GIF simulation (although a slightly different transfer function for the input power spectrum.) The particle mass in this simulation is larger, by a factor of roughly 5, than in the GIF simulation, so that only dark matter halos more massive than $`7\times 10^{11}h^1M_{}`$ are resolved, but the volume is approximately 40 times larger than in the GIF simulation. The $`512^3`$ simulation is ideal for studying the brightest galaxies which are only found in very massive halos and have low abundance. ? noted that the dark matter correlation function in the $`512^3`$ simulation differed slightly from that in the GIF simulation due to the large-scale power which is included in the larger volume but is missing in the smaller one. Throughout this paper, we apply a small correction to all the correlation functions determined from the GIF simulation by adding the quantity $`\mathrm{\Delta }\xi (r)=\xi _{512^3}(r)\xi _{\mathrm{GIF}}(r)`$ (where $`\xi _{512^3}`$ and $`\xi _{\mathrm{GIF}}`$ are the correlation functions of the dark matter, in real or redshift space as appropriate, in the $`512^3`$ and GIF simulations respectively). The required correction is not necessarily the same for galaxies and dark matter of course. However, the correction is at most 20% in $`\xi (r)`$ over the range of scales considered in this work and, furthermore, ignoring it does not alter any of our conclusions.

## 3 Images of the galaxy distribution

Fig. 1 shows the distributions of dark matter and galaxies in slices through the GIF N-body simulation volume at six different redshifts. Each slice has dimensions $`141.3\times 141.3\times 8h^3\text{Mpc}^3`$. The dark matter is represented by the greyscale, obtained by adaptively smoothing the N-body mass distribution. The shade intensity is proportional to the logarithm of the dark matter density in each pixel (the darker the pixel, the greater the projected density of the slice). Each galaxy brighter than $`M_\mathrm{B}5\mathrm{log}h=19`$ that formed in this region is depicted as a circle whose size is proportional to the rest-frame B-band absolute magnitude and whose colour indicates the rest-frame B-V colour, as given by the scales at the top of the figure. Fig. 2 displays zoomed images, at three redshifts, of the areas delineated by the coloured boxes in Fig. 1: the region around a supercluster (left), and a more typical region (right).

### 3.1 The biased galaxy distribution

The images in Fig. 1 illustrate what a realistic distribution of galaxies might look like. At $`z=0`$, this patch of universe is populated by galaxies with a wide range of colours. The dark matter has acquired the filamentary appearance characteristic of gravitational growth from cold dark matter initial conditions. The galaxies “light up” the filaments and superclusters of the dark matter distribution, but are conspicuously absent from regions in which the density of dark matter is low. The emptiness of the voids is quite striking (c.f. ?). The voids would not be as empty of galaxies as they are if, instead of using the semi-analytic model, “galaxies” were simply identified with randomly selected dark matter particles. Such a Poisson process inevitably places a small fraction of galaxies in voids, but the physics of galaxy formation do not allow this: at the magnitude limit appropriate to Figs 1 and 2, voids are truly voids: no galaxies form in them.

While to the eye the distribution of galaxies in the figures may appear to follow the broad features of the dark matter distribution, it is impossible to judge in this way how faithful this tracing really is. That galaxies in CDM models might not be perfect tracers of the mass has been suspected for some time (??) and was demonstrated explicitly for the GIF simulations by ? and ?. These studies showed, for example, that on scales below a few megaparsecs, galaxies in the simulations are less strongly clustered than the dark matter (or “anti-biased” relative to the mass). The upper panel in Fig. 3 shows the correlation functions of the galaxies for the different redshifts shown in Fig. 1 and demonstrates one aspect of the bias. The correlation functions of these galaxies have remarkably similar amplitudes, in contrast to the correlation function of dark matter (whose value at $`r=5h^1`$Mpc is shown as solid points in Fig. 3), the amplitude of which evolves rapidly with redshift. In Fig. 3, we also compare our determinations of the correlation function to those derived from an analytical approach to the bias (c.f. ?). For this, we compute the effective bias of the galaxy population as:

$$b(z)=\frac{1}{N}\underset{i=1}{\overset{N}{}}b(M_i,z),$$
(1)
where $`N`$ is the number of galaxies in the simulation, $`M_i`$ is the mass of the halo in which the $`i^{\mathrm{th}}`$ galaxy is found, and $`b(M,z)`$ is the bias of dark matter halos of mass $`M`$ and redshift $`z`$, which we calculate using the fitting formula of ? (which is based on the model of ?). An approximation to the galaxy correlation function is then $`b(z)^2\xi _{\mathrm{DM}}(r,z)`$, where $`\xi _{\mathrm{DM}}`$ is the correlation function of dark matter. (For the galaxy samples in our simulation, we find $`b(z)=1.07`$, 1.24, 1.47, 2.05, 2.91 and 5.21 for $`z=0.0`$, 0.5, 1.0, 2.0, 3.0 and 5.0 respectively.) In the lower panel of Fig. 3 we plot the ratio of the galaxy correlation functions in the simulation to the analytical approximation. On scales above a few megaparsecs the analytical bias approach works well (it does slightly underestimate the correlation functions, but given the large biases present in our high-redshift samples, the approximation is actually rather good), but, as expected, it fails on smaller scales where our model predicts a scale dependent bias.

Certain kinds of bias are readily apparent to the eye in Fig. 1. For example, the largest dark matter clumps at $`z=0`$ are preferentially populated by red galaxies, while the field contains a mixture of galaxy colours. Similarly, the brightest galaxies in the region are also preferentially found at the centres of rich clusters. The first of these biases, the “colour-density” relation is intimately related to a “morphology-density” relation and is a natural outcome of hierarchical clustering, as we will discuss in more detail below. It has been investigated before in these simulations, in a somewhat different form, by ? and ?, who found that the two-point correlation function of elliptical galaxies is higher than that of spirals on small scales. This kind of bias is, of course, known to occur in the real universe (e.g. ??).

Biases in the distribution of galaxies are the inevitable by-product of the complex physics of galaxy formation. They affect galaxies of different types in different ways and need to be understood before attempting to interpret cosmological observations.

### 3.2 The evolution of the galaxy population

The evolution of galaxies is driven by a number of processes. The most obvious one is the aging of the stellar populations. Even if galaxies lived in isolation, stellar evolution would cause their luminosities and colours to evolve. However, no galaxy is an island: accretion, mergers and interactions are common. Many galaxies are observed to be forming stars today, possibly by converting new gas supplied externally (although see ?), thus increasing their mass and size. Others are observed to be expelling gas through galactic winds and tidal encounters. Nucleosynthesis in stars causes the metallicity of galactic gas and stars to evolve and this, in turn, affects the integrated stellar spectra.

The nett effects of the various processes driving galaxy evolution are readily apparent in Fig. 1. Since young stars produce copious amounts of blue light, galaxies with high star formation rates relative to their total stellar mass appear bluer than galaxies with low relative star formation rates. At $`z=0`$, the majority of the galaxies in the simulation are quite red because they have low relative star formation rates and are made primarily of old stars. As we look into the past, the appearance of the population rapidly changes. By $`z=1`$, galaxies are much bluer because the typical relative star formation rate is higher than at $`z=0`$ and the average stellar age is younger. Beyond $`z=1`$, the galaxies remain blue, reflecting the youth of their stellar populations. The apparent star formation rate per unit volume declines at high redshift because fewer galaxies are seen above our magnitude selection, even though those that are seen still have high star formation rates. The detailed star formation history in our model is discussed in ?. It is broadly in agreement with the star formation history of the Universe as inferred from observations by, amongst others, ?, ? and ?. Although many quantitative details of the observations remain uncertain (due to complications arising, for example, from dust obscuration), the general behaviour seems to consist of a rapid rise in star formation rate between $`z=0`$ and $`z=1`$, followed by a slowly declining (or perhaps constant) star formation rate at higher redshifts. This is the kind of behaviour exhibited by our simulations.

As the images in Fig. 1 illustrate, in hierarchical models of galaxy formation the number of galaxies is constantly changing. Galaxies are born as new dark matter halos form and gas is able to cool in them and turn into stars. The population is depleted when galaxies merge together. Of course, the number of galaxies detected in a particular survey will depend crucially on the selection criteria. All these effects can be clearly seen in the images of our simulations. At $`z=5`$, there are very few galaxies present because only a handful of massive dark halos have had time to collapse. In those that have, galaxies have had little time to form stars, while feedback from supernovae has strongly suppressed star formation in small halos. The majority of the galaxies seen in the images at this epoch occur in halos of mass $`10^{1112}h^1M_{}`$ and have stellar masses of a few times $`10^9h^1M_{}`$; the very brightest galaxies are found in the tail of halos extending to masses close to $`10^{13}h^1M_{}`$. By $`z=3`$ the number of galaxies has increased significantly, as more halos have collapsed and more galaxies have been able to form. At this epoch several extremely bright (in the B-band) galaxies are visible. The increased number of galaxies in the image is due, in part, to our selection in the B-band which is sensitive to the star formation rate. At $`z=2`$ there is a noticeable increase in the abundance of galaxies as structures continue to form. The most obvious change from $`z=2`$ to $`z=1`$ is a substantial reddening of the galaxies, a trend which continues to $`z=0`$ as star formation rates decline and stellar populations age.

Many faint blue galaxies are formed in the filamentary network of the dark matter. In the comoving coordinates of Fig. 1, galaxies move rather little between $`z=5`$ and the present. For example, the progenitor of the large supercluster marked by a red box at $`z=0`$ is already clearly visible at $`z=3`$ as a concentration of young galaxies. In other words, by virtue of forming in the highest density regions, galaxies are strongly biased at birth. This is a fundamental outcome of hierarchical clustering (??). It underlies the results of ? and ?, who argued that Lyman-break galaxies at $`z=3`$ would be expected to be strongly clustered, as was subsequently found to be the case observationally (?).

## 4 Quantitative analysis and comparisons with observations

In this section, we compare the properties of our model with the limited observational data currently available. More stringent comparisons will be possible with the forthcoming 2dF and Sloan galaxy redshift surveys. We consider, in turn, the variation of the clustering length with the luminosity of a sample, the morphology-density relation and the topology of the galaxy distribution.

### 4.1 Correlation length versus separation

The evolution of the two-point correlation function of dark matter in CDM models is now well established. In the linear regime, it evolves according to the linear growth factor; in the non-linear regime, its evolution can be calculated accurately using N-body simulations (see, for example, ?). By contrast, the evolution of the two-point correlation function of galaxies has only begun to be investigated in detail recently. In Fig. 4 we plot the redshift-space correlation lengths, $`s_0`$ (top panel), and the real-space correlation lengths, $`r_0`$ (middle panel), of galaxies brighter than a particular rest-frame (dust-extinguished) B-band magnitude, as a function of their mean separation, $`d`$, at $`z=0`$ (solid line), $`z=1.0`$ (dashed line) and $`z=3`$ (dot-dashed line). For reference, the lower panel of Fig. 4 shows the relation between mean separation and absolute B-band magnitude. Statistical errors (obtained by assuming Poisson counting statistics to estimate the error in $`\xi (r)`$ and propagating this error through to the determination of $`r_0`$) are shown by the error bars. At small pair separations, we show results from the GIF simulation (applying the small correction for finite volume effects discussed above), and at large pair separations, we use the $`512^3`$ simulation which is more accurate on large scales. (The sudden decrease in the size of the errorbars at $`15h^1`$Mpc is due to this change of simulation volume.) For reference, we show the correlation lengths of the dark matter at $`z=0`$ and $`z=1`$ as horizontal arrows (at $`z=3`$ the dark matter correlation length is $`1.1`$ and $`1.5h^1`$Mpc in real and redshift-space respectively).

We compare our model predictions with results from the Southern Sky Redshift Survey 2 (SSRS2) (?), which are shown by filled squares in Fig. 4. These data exhibit a nearly constant redshift-space clustering length for galaxies fainter than $`M_\mathrm{B}5\mathrm{log}h=20`$, followed by a rapid increase at brighter magnitudes. ? compared their data to two simple models of galaxy bias. In the first, they assigned a dark matter halo mass to galaxies in their sample using the Tully-Fisher and Faber-Jackson relations, and then applied the techniques of ? to compute the bias of these halos (and hence of the galaxies which occupy them). This model provides a reasonable match to the observed behaviour of the faint galaxies in the survey, but it is unable to reproduce the strong luminosity dependent bias observed for galaxies brighter than $`L_{}`$. Their second model is based upon the work of ?, who developed a description of bias from the non-linear evolution of the density field. This model is able to match the luminosity dependent bias for bright galaxies, but predicts too strong a relation for faint galaxies and so is also ruled out by the data. Our model of galaxy clustering, on the other hand, does produce a trend similar to that observed, namely a relatively constant clustering length for small $`d`$, followed by a rise in correlation length for the rarest objects. From Fig. 4 it appears that our model is not consistent with the data over the whole range of separations. However, this discrepancy may result from sampling variance in the observations, as we discuss below. First, we explain how this trend arises in our models.

It is well known that dark matter halos are biased relative to the underlying mass, with the most massive halos being the most strongly clustered (e.g. ??). Thus, for the brightest galaxies to be the most strongly clustered, it is necessary that they should preferentially inhabit more massive dark matter halos than those occupied by their lower luminosity counterparts. In Fig. 5 we show mass functions of dark matter halos, weighted by the number of occupant galaxies (solid lines) and unweighted (dashed lines), for three values of $`d`$. In each panel, the solid histogram gives the galaxy-weighted halo mass function when galaxies are selected by their dust-extinguished B-band magnitude.The vertical arrow in each panel shows the location of $`M_{}`$, defined by $`\sigma (M_{})=\delta _\mathrm{c}`$, where $`\sigma (M)`$ is the mass variance in spheres containing a mass $`M`$ on average, and $`\delta _\mathrm{c}`$ is the critical overdensity for collapse in the spherical top-hat model. A simple understanding of the galaxy bias for each sample may be gained from this figure. For the two smaller values of $`d`$, the relative numbers of galaxies in highly-biased, cluster-sized halos ($`M_{\mathrm{halo}}>10^{14}h^1M_{}`$) and in weakly-biased, galaxy-sized halos ($`M_{\mathrm{halo}}10^{12}h^1M_{}`$) are comparable in the two samples. As a result, these two samples have quite similar correlation lengths. However, for the sparser sample with $`d=25h^1`$Mpc, the relative number in cluster-sized halos is much higher than in the two other cases. As a result, this sample has a larger correlation length than the other two. (For this very bright sample, we find that many galaxies occupying halos near the peak mass of $`M_{\mathrm{halo}}10^{14}h^1M_{}`$ have undergone a recent burst of star formation.)

As noted above, our model is not in perfect agreement with the data of ?. The SSRS2 survey, however, covers a relatively small volume and so sample variance may not be negligible. ? estimated the effects of sample variance on their results and concluded that the luminosity dependence of $`s_0`$ for faint (sub-$`L_{}`$) galaxies could well be due to sample variance, while that for brighter galaxies seemed to be a real effect. We can estimate the size of this uncertainty directly from our simulations, which cover a much larger volume than the real survey. For a given value of $`d`$, we extract from the simulation fifty randomly placed cubic regions of volume equal to that of a volume-limited SSRS2 catalogue cut at the same absolute magnitude. We then measure the correlation length in each of the fifty cubes. In Fig. 4 we plot the median correlation lengths from the cubes as open circles, with error-bars indicating the 10% and 90% intervals of the distribution. (We note that even for the largest value of $`d`$ shown, the number of independent SSRS2 volumes that fit within the $`512^3`$ simulation cube is still reasonably large, $`20`$, and so our estimates of sample variance at these separations should still be accurate.) Evidently, sample variance in a survey the size of the SSRS2 is large, and this can account for the differences with our model. At the smallest values of $`d`$, the median value of $`s_0`$ from the subsamples is biased low relative to that measured in the full simulation volume, because a large fraction of the clustering signal comes from galaxies in and around a few large clusters and these are often missing from small volumes cut out of the simulation box. At larger separations, the main source of sample variance is the low abundance of the brightest galaxies. Our estimate of sample variance confirm the conclusion reached by ?, namely that the luminosity dependence of $`s_0`$ observed for sub-$`L_{}`$ galaxies in the SSRS2 is due to sample variance, while that for brighter galaxies is real.

The evolution of the redshift-space correlation lengths of galaxies of different abundance is also illustrated in Fig. 4. For samples selected according to our specific criteria (i.e. according to rest-frame B-band absolute magnitude), the correlation lengths at $`z=1`$ and $`z=3`$ are slightly smaller than at $`z=0`$. This variation, however, is much smaller than the variation in the correlation length of the dark matter, as indicated by the horizontal arrows in Fig. 4. Thus, galaxies selected according to our criteria are more strongly biased at $`z=1`$ and $`z=3`$ than at $`z=0`$, even for quite small values of $`d`$. The images of the galaxy distribution displayed in Fig. 1 show exactly how this effect arises. The galaxies present at the highest redshift, $`z=3`$, have formed in regions that are destined to become incorporated into clusters or superclusters of galaxies by the present day. Such regions are amongst the most overdense at high redshift, and so small-scale density fluctuations tend to collapse earlier there than in less overdense regions (???). Forming as they do in the most highly biased regions of the universe, galaxies at high redshift naturally end up being strongly biased themselves. Our model predicts a strong luminosity dependence of clustering length for the brightest galaxies from $`z=0`$ to 3, an effect which is seen in both real and redshift-space

### 4.2 Morphology-density relation

The existence of a correlation between galaxy morphology and local environment is a remarkable feature of the galaxy population. ? showed that the fraction of galaxies of different morphological types is strongly correlated with the local galaxy density: elliptical and S0 galaxies are found preferentially in high density regions while spiral galaxies are found preferentially in low density regions.

Early N-body simulations suggested that a morphology-density relation is a natural outcome of hierarchical clustering from CDM initial conditions (?). This was explicitly shown to be the case by Monte-Carlo based semi-analytic modelling (??). These calculations had no information on the spatial distribution of galaxies and so the relation they established is one between morphology and cluster mass, rather than between morphology and galaxy density. The implementation of semi-analytic techniques in high resolution N-body simulations allowed the radial distributions of different kinds of galaxies inside large clusters to be calculated for the first time (??). This work has shown, for example, that galaxies of different colours are spatially segregated within the cluster.

The emergence of a colour-density relation is clearly illustrated in the images of Fig. 1. Redder galaxies (with $`BV>0.7`$) which, in our model, are primarily ellipticals or S0’s (see Fig. 12 of ?), are over-abundant in the most overdense regions relative to the field (e.g. the large supercluster in the middle of the left-hand edge, or the large cluster near the centre of the bottom edge). The correlation between galaxy colour (or morphology) and environment is a byproduct of the biases discussed in the preceding subsection: the oldest, reddest galaxies form in the highest density regions where the production of elliptical galaxies by mergers is also favoured.

To quantify the morphology-density relation apparent in Fig. 1, we proceed in a manner analogous to the analysis of an observational sample. We use a technique patterned on that employed by ?. First, we apply a friends-of-friends group-finding algorithm to the real-space distribution of model galaxies brighter than a particular B-band absolute magnitude, using different values of the linking length. (? applied their group finder in redshift-space, but used an anisotropic linking length to account for distortions in the redshift direction. Since we have a galaxy catalogue in real-space we perform group finding there, avoiding the complications of redshift-space distortions.) In this manner, we build up nested sets of groups as a function of the enclosed density (we consider three or more galaxies linked together to be a “group”). For very large linking lengths all galaxies will belong to a group, but as the linking length is decreased each galaxy will at some point no longer be a member of a group. We assign each galaxy a local density corresponding to the surface density of the group of which it was last a member. The local density at the surface of a group formed with a linking length $`r`$ is approximately $`n=3/2\pi r^3`$ (e.g. ?). The morphological type of the galaxy is assigned according to our standard definitions based on dust-extinguished B-band bulge-to-total luminosity ratios: galaxies with B/T$`{}_{\mathrm{B}}{}^{}<0.4`$ are labelled as spirals, while those with B/T$`{}_{\mathrm{B}}{}^{}0.4`$ are labelled as E/S0 (?). We constructed density-morphology relations for samples with different limiting absolute magnitudes. In each case, the measured densities were corrected to the density of galaxies at $`M_\mathrm{B}5\mathrm{log}h=17.5`$ using the CfA survey luminosity function, as was done by ?, namely we multiply the densities by a factor

$$f_\mathrm{n}=\frac{_{L_{\mathrm{ref}}}^{\mathrm{}}\mathrm{\Phi }(L)dL}{_{L_0}^{\mathrm{}}\mathrm{\Phi }(L)dL},$$
(2)
where $`L_0`$ and $`L_{\mathrm{ref}}`$ are the luminosities corresponding to $`M_\mathrm{B}5\mathrm{log}h19.5`$ and $`17.5`$ respectively and $`\mathrm{\Phi }(L)`$ is the luminosity function of the CfA survey.

The left-hand panel of Fig. 6 shows the model morphology-density relation at $`z=0`$ for three absolute magnitude cuts: $`M_\mathrm{B}5\mathrm{log}h=18.1`$ (dotted line; the completeness limit of the GIF simulation), $`19.5`$ (thin solid line; close to $`L_{}`$) and $`20.1`$ (dashed line; the completeness limit of the $`512^3`$ simulation). It is immediately apparent that our model does display a morphology-density relation with the correct trend: ellipticals/S0s are more common in high density environments. For high densities ($`>100`$ galaxies $`h^3`$Mpc-3) our model shows no relation since, by construction, no morphology-density relation can exist within individual halos because of the way in which we assign galaxies to dark matter particles. (We do not plot model results for densities greater than $`1000h^3`$Mpc-3 since these begin to probe single dark matter halos in our simulation, resulting in a poor determination of the morphology-density relation.) The right-hand panel of Fig. 6 shows that our model predicts a very similar morphology-density relation at $`z=1`$ as at $`z=0`$.

We can now elaborate a little on the cause of the morphology density relation in our model, focusing on the morphology-density relation at low densities (where the spatial distribution of morphological types within individual halos is unimportant and the local density is typically determined by averaging over regions containing many dark matter halos). The morphological mix of galaxies in a single dark matter halo can depend only upon the mass of that halo since this statistically determines the merger history and galaxy formation history in the halo. Therefore, for a morphology-density relation to exist: (i) there must be a relation between dark matter halo mass and local galaxy density and (ii) there must be a dependence of morphological fraction on halo mass. The first of these requirements is naturally met in hierarchical cosmologies since the most massive halos are preferentially found in the densest environments. The left-hand panel of Fig. 7 shows the local galaxy number density (for galaxies brighter than $`M_\mathrm{B}5\mathrm{log}h=19.5`$ but corrected to $`17.5`$ as before) at the centres of dark matter halos as a function of the halo mass (i.e. we plot the density of the three-particle group, as defined above, which contains the central galaxy). Clearly, the densest regions of the galaxy population are associated with cluster-sized halos. The right-hand panel of Fig. 7 shows that our model also meets the second requirement. This figure shows the dark matter halo mass function weighted by the number of galaxies of a particular morphological class per halo. The halo mass function for E/S0 galaxies is shifted to higher mass halos relative to that for S galaxies. The resulting dependence of morphological mix on halo mass merely reflects the fact that galaxy formation in clusters is accelerated relative to the field, allowing enough time for galaxy-galaxy mergers to produce a large population of elliptical galaxies.

Fig. 6 also shows that the model morphology-density relation depends upon the absolute magnitude at the which the galaxies are selected. This property is simply a reflection of a morphology-luminosity relation that is present in our model. This property complicates the comparison with observational data since existing analyses are usually based on apparent magnitude limited samples. In the figure we compare our model predictions with the density-morphology relation measured in the CfA survey by ?. For densities less than $`100h^3`$Mpc-3 (indicated by the vertical arrow), we plot the E/S0 fraction in “groups” (defined as associations of 3 to 9 galaxies by ?), and for larger densities, we plot the fraction in clusters (associations of 10 or more galaxies). This shows that our model displays qualitatively similar behaviour to the observational data. A more detailed comparison may be possible with the new generation of large redshift surveys.

Finally, another important prediction of our model is that there should be a strong morphology-density relation well-established already at $`z=1`$, as shown in the right-hand panel of Fig. 6. At this high redshift, the model relation is qualitatively similar to that at $`z=0`$.

### 4.3 The genus curve for the topology of the galaxy distribution

The two-point correlation function contains only low order information about the spatial distribution of galaxies. To fully specify this distribution requires determining its higher order clustering properties. Alternatively, the genus, a measure of the topology of a smooth density field, provides a statistic that is sensitive to all of the higher order moments of the distribution (??).

The genus is defined as the number of topological holes minus the number of isolated regions of an isodensity surface. By varying the density at which this surface is placed, a genus curve can be constructed. The genus curve has the interesting property that an exact, analytic expression exists for the special case of a Gaussian random density field, namely

$$g(\nu )=A(1\nu ^2)\mathrm{exp}\left(\frac{\nu ^2}{2}\right),$$
(3)
where $`g`$ is the genus per unit volume and $`\nu `$ is defined by

$$\nu =\sqrt{2}\mathrm{erf}^1(12f),$$
(4)
where $`f`$ is the fraction of the volume above the density threshold and $`\mathrm{erf}^1`$ is the inverse of the error function. (For a Gaussian random field, but not for any other field, this definition implies that $`\nu ^2`$ is the variance of the field.) The amplitude, $`A`$, depends only on the second moment of the power spectrum of the smoothed density field. In general, if the field is not Gaussian, the shape of the genus curve may differ from eqn. (3). With the above definition for $`\nu `$, the genus curve remains the same under any dynamical evolution or biasing in which the initial and final densities at each Eulerian point are related by a monotonic, one-to-one mapping. Thus, this topological measure provides, in principle, a method for testing whether or not density fluctuations in the early universe were originally random and Gaussian, as predicted in generic inflationary models.

The genus curve has been determined for both N-body simulations of dark matter and for surveys of galaxies in the real Universe by many authors (e.g. ????????.) The genus curve for galaxies would be identical to that of the underlying dark matter if the galaxy and dark matter density fields were related by a monotonically varying factor, the bias. The bias relation (i.e. the galaxy overdensity as a function of dark matter overdensity) in our model is shown in Fig. 8 for galaxies brighter than $`M_\mathrm{B}5\mathrm{log}h=18.5`$. It is qualitatively similar to the relations found by ? using similar techniques. The symbols show the median relation which is clearly monotonic. However, there is substantial scatter around this relation as a result of which it is no longer guaranteed that the genus curve for galaxies will be identical to that of the dark matter.

We have used the technique described by ? to measure the genus of both dark matter and galaxies in our simulation. Fig. 9 shows the results for galaxies (dashed lines) brighter than $`M_\mathrm{B}5\mathrm{log}h=18.5`$ at $`z=0`$ (upper panels) and $`z=1.0`$ (lower panels) and also for dark matter (solid lines). These curves were calculated by smoothing the dark matter and galaxies in redshift-space onto a $`128^3`$ grid using a Gaussian filter of the form

$$W(r)=\frac{1}{\pi ^{3/2}\lambda _\mathrm{e}^3}\mathrm{exp}\left(\frac{r^2}{\lambda _\mathrm{e}^2}\right),$$
(5)
as is conventional in the literature on this subject. Smoothing lengths of $`\lambda _\mathrm{e}=6.0`$ and $`\lambda _\mathrm{e}=8.0h^1`$ Mpc were chosen to match those used by ? in their analysis of the CfA surveys. These smoothing lengths are over five times larger than the size of the grid cells on which the fields are tabulated and so the finite resolution of the grid has no effect on the calculation of the genus curve (?). Note that the mean comoving separation of galaxies in our simulation is $`4.3`$ and $`5.6h^1`$ Mpc at $`z=0`$ and $`z=1.0`$ respectively (corresponding to 34000 and 16000 galaxies in the simulation volume). At $`z=1`$ this is just smaller than the minimum smoothing scale, thereby providing the greatest number of independent resolution elements without allowing discreteness effects to become too large (?). We estimate the errors on each genus curve by bootstrap resampling of the galaxy catalogues. ? find that this procedure produces slightly larger errors than those estimated by considering several realizations of a mock catalogue.

In a $`\mathrm{\Lambda }`$CDM cosmology, N-body simulations have shown that the genus curve for the dark matter displays both a “*bubble shift*” (i.e. a shift to the right) with respect to the Gaussian random phase genus curve and an amplitude reduction relative to a Gaussian random density field having the same power spectrum (?). These two effects can be seen by comparing the dark matter genus curve to the two random phase genus curves (i.e. with the shape given by eqn. 3) shown in Fig. 9 by dotted and dot-dashed lines. The amplitude of the dotted curve is chosen to best fit (in a least-squares sense) the measured dark matter genus curve, while the amplitude of the dot-dashed curve is that expected for a Gaussian random field with the same power spectrum as the dark matter. The dark matter genus curve is displaced to the right of the dotted curve, showing the bubble shift, and has smaller amplitude than the dot-dashed curve, showing the amplitude drop. The genus curve for galaxies brighter than $`M_\mathrm{B}5\mathrm{log}h=18.5`$ shows definite differences from that for the dark matter. Firstly, no bubble shift exists for the galaxy genus curve which, instead, exhibits a small “meatball shift” (i.e. shift to the left). The galaxy curve also has a systematically larger amplitude than the dark matter curve. ? advocate the amplitude drop (defined as the ratio of the amplitudes of the best-fit random phase genus curved for the actual density field and that of a Gaussian random field with the same power spectrum) as a useful measure of the degree of phase correlation in the galaxy density field. To measure the amplitude drop, we “Gaussianise” the galaxy density field (i.e. we take it to Fourier space, randomise the phases subject to the reality condition, $`\delta _k=\delta _k^{}`$, and then restore it to real space). The amplitude drop is then simply the ratio of amplitudes of the best-fit random phase genus curves for the original and Gaussianized density fields. At $`z=0`$ we find amplitude drops of $`R=0.70`$ and $`R=0.84`$ for the dark matter in redshift-space, for $`\lambda _\mathrm{e}=6.0`$ and $`8.0h^1`$Mpc respectively. (In real-space we measure amplitude drops of 0.60 and 0.67 for the same two smoothing lengths, in good agreement with the determinations of ?.) For galaxies at $`z=0`$ we find $`R=0.84\pm 0.02`$ and $`R=0.90\pm 0.02`$ also for these same two smoothing scales. At $`z=1`$ (lower panels in Fig. 9) the amplitude drops are somewhat smaller, $`R=0.90\pm 0.03`$ and $`R=0.93\pm 0.07`$ for $`\lambda _\mathrm{e}=6.0`$ and $`8.0h^1`$Mpc respectively, as phase correlations due to non-linear growth of structure have not had as long to develop as at $`z=0`$.

As noted above, galaxies would have exactly the same genus curve as the dark matter if there were a one-to-one mapping between dark matter and galaxy density fields which preserved the density ranking. However, we do see significant differences between the galaxy and dark matter genus curves. This must be due either to the scatter in the biasing relation between galaxies and dark matter, or to systematic biases arising from the relatively small number of galaxies $`(10^4)`$ compared to dark matter particles $`(10^7)`$ in our samples. To test this latter possibility, we extracted twenty random samples of dark matter particles with the same abundance as the galaxies in our catalogue and computed their genus curves. We find that this sparse sampling is the primary cause of the differences between the galaxy and dark matter genus curves. Just as for the galaxy sample, the sparsely sampled dark matter shows no evidence for a bubble shift (and agrees closely with the galaxy genus curve for $`\nu <0`$) and also shows a higher genus curve amplitude compared to the fully sampled dark matter distribution. We conclude therefore that any differences between the genus curves for dark matter and galaxies due to the stochastic bias of Fig. 8 are negligible compared to the effects of sparse sampling of the galaxy density field. This sparse sampling at present severely limits the usefulness of the genus statistic for quantifying the Gaussianity of the initial dark matter distribution.

We also show in Fig. 9 the genus curve measured for the CfA surveys by ?. While our model results are in reasonable agreement with these data, the pronounced features in the data suggest that the CfA surveys do not have a sufficiently large volume to avoid significant sample variance effects. Indeed, when we extract “CfA survey” volumes from random locations in our simulations and measure their genus curves, we find that excursions such as those seen in Fig. 9 are very common.

## 5 Discussion and conclusions

We have implemented a semi-analytic model of galaxy formation in high resolution N-body simulations of the $`\mathrm{\Lambda }`$CDM cosmology in order to study the spatial distribution of galaxies and its evolution. The semi-analytic model requires setting values for a number of parameters which describe the physical processes that are modelled, such as gas cooling, star formation and the associated feedback mechanisms and chemical evolution, galaxy merging, the evolution of stellar populations, etc. In keeping with the general philosophy of our work in the subject, we have fixed all the model parameters by requiring a match to a handful of global properties of the local galaxy population, with the largest number of constraints coming from the local B-band and K-band luminosity functions. No further adjustment to these parameter values was allowed in the clustering study carried out in this paper. Thus, our clustering results are genuine predictions of the model and offer an opportunity to test the validity of the physical assumptions it requires as well as the realism of the $`\mathrm{\Lambda }`$CDM model as a whole. In this paper, we have considered three specific statistical measures of clustering: the correlation length (in real- and redshift-space) of samples of galaxies of different luminosity, the morphology-density relation and the genus curve. At $`z=0`$, our model may be tested by forthcoming data from the 2dF and Sloan surveys; at $`z=1`$, it may tested by planned surveys such as DEEP (?) and VIRMOS (?).

The results presented here extend and complement those presented in earlier papers in this series (??), as well as in the series of papers by Kauffmann and collaborators (???) who analysed one of the N-body simulations that we have analysed here, but using their own semi-analytic model. In ?, we examined the physical and statistical processes that segregate the galaxies from the dark matter and we showed that the (real-space) two-point galaxy correlation function in a $`\mathrm{\Lambda }`$CDM model that produces an acceptable galaxy luminosity function is in remarkably good agreement with observations. In ?, we considered clustering in redshift space and analysed the resulting distortions of the two-point correlation function, as well as its dependence on galaxy luminosity, morphology and colour. Kauffmann and collaborators studied many of these properties too, at the present day (?), and at high redshift (?), as well as the clustering of groups (?). On the whole, these two independent analyses agree quite well and the differences that do exist can be readily understood in terms of differences in the detailed assumptions for the physics of galaxy formation (see ? for a detailed discussion of differences between the two models). The evolution of clustering has also been studied using similar semi-analytic/N-body techniques by ? and, more recently, by ?.

Images of our simulation clearly illustrate many of the salient features of galaxy growth by hierarchical clustering. They show that galaxies approximately trace the filamentary structure and avoid the lowest density regions of the dark matter distribution, that the redder galaxies tend to predominate in the most massive dark halos and that the brightest galaxies occur almost exclusively in regions of high dark matter density. A series of time slices shows how the galaxy population changes in abundance and colour with the passage of time and demonstrates the primary effect behind biased galaxy formation: the formation of the first bright galaxies in regions of exceptionally high dark matter density.

From quantitative studies of the galaxy distribution we reach the following three conclusions:

1. The correlation length of galaxies in real and redshift-space increases rapidly with galaxy luminosity for galaxies brighter than $`L_{}`$, both at $`z=0`$ and at high redshifts, $`z<3`$.
2. A strong morphology-density relation, in the same sense as observed, is a natural outcome of hierarchical clustering from CDM initial conditions, and radpidly develops in our simulation. A clear morphology-density relation is predicted to be already in place at least since $`z=1`$.
3. The topology of the galaxy distribution, as measured by the genus statistic, differs significantly from that of the dark matter. However, the differences are due almost entirely to sparse sampling effects; the stochastic biasing between galaxies and dark matter is, at most, a minor effect.

We now discuss these points more detail.

The variation of the correlation length with luminosity or, equivalently, with mean inter-galaxy separation is one of the most striking results of our analysis. The correlation length is virtually insensitive to the mean separation within the sample out to separations of around $`10h^1`$Mpc (corresponding to luminosities of $`M_\mathrm{B}5\mathrm{log}h20.5`$), but for brighter, sparser samples it increases very rapidly. Thus, the redshift-space clustering length of galaxies of luminosity $`7L_{}`$ is predicted to be over twice as large as that of $`L_{}`$ galaxies. The pattern is similar in real- and redshift-space. Owing to the large volume of our simulations, this is the first time that this rapid increase in correlation length at the brightest luminosities has been unambiguously demonstrated. The main cause of this behaviour is the preponderance of such galaxies at the centres of massive clusters. These galaxies experience enhanced merger rates at early times and are the beneficiaries of late accretion of cool gas which, in our model, is always funnelled onto the central galaxy in each halo. Our predictions for the dependence of clustering strength on galaxy luminosity are in broad agreement with existing datasets, but these are rather small and thus subject to considerable sampling uncertainties. A better test of these predictions should be forthcoming shortly from the 2dF and Sloan surveys.

Our simulations develop a strong morphology-density relation similar to that observed in the local universe: elliptical/S0 galaxies predominate in rich clusters while spirals predominate in the field. We have used a technique patterned after observational procedures to characterise the morphology-density relation in our simulations and find that it quantitatively agrees rather well with observations. The cause of the morphology-density relation is closely related to the reasons behind the luminosity dependence of clustering. In rich clusters, galaxy evolution is “accelerated” relative to more ordinary regions of space, thus allowing sufficient time for mergers and interactions to build up a large population of bright elliptical/S0 galaxies. Remarkably, a strong morphology-density relation is already well-established by $`z=1`$.

Finally, we have investigated the topology of the galaxy distribution, providing the first theoretical prediction for the genus curve of galaxies, rather than merely of dark matter. Of course, if galaxy density were monotonically related to dark matter density, the two curves would be the same. In the simulations there is on average a monotonic relation between the two, but it has such large scatter that it does not preclude differences in the respective genus curves. We do actually find a difference in our simulations: the dark matter genus curve has a “bubble” shift whereas the galaxy genus curve has a “meatball” shift and also a higher amplitude. It turns out, however, that the differences are not due to genuine topological differences but rather to the sparse sampling of the density field provided by galaxies. The confusing effects of sampling have been pointed out by previous authors (e.g ??). Our simulations show that, in order to measure an unbiased genus curve for a clustered galaxy distribution (particularly at low overdensity), several hundred galaxies per smoothing volume are required. This is unlikely to be practical even with the new generation of large surveys and so we will have to live with these biases. The best approach for comparing models and data is therefore to analyse each in identical ways, so as to cancel out any systematic effects. Mock galaxy catalogues such as those presented here will be crucial for this approach.

To summarise, the combination of high-resolution N-body simulations of dark matter and semi-analytic modelling of galaxies provides a powerful technique for turning cosmological and galaxy formation theory into realistic realizations of the galaxy population that can be compared in detail with observations. Tests of this sort will become increasingly common with the new generation of galaxy surveys. In this way, it will be possible to extract not only the cosmological information encoded in the clustering pattern, but also valuable information regarding the physics of galaxy formation.

## Acknowledgements

AJB, SMC and CSF acknowledge receipt of a PPARC Studentship, Advanced Fellowship and Senior Fellowship respectively. CSF also acknowledges a Leverhulme Research Fellowship. CGL acknowledges support at SISSA from COFIN funds from MURST and funds from ASI. This work was supported in part by a PPARC rolling grant, by a computer equipment grant from Durham University and by the European Community’s TMR Network for Galaxy Formation and Evolution. We acknowledge the Virgo Consortium and GIF for making available the simulations used in this study.


format: text

text:
# Open charm enhancement by secondary interactions in relativistic nucleus-nucleus collisions?supported by DFG/RFFI and the Landau-Heisenberg program

## 1 Introduction

In the last decade the interest in hadronic states with charm flavors ($`c,\overline{c}`$) has been rising continuously in line with the development of new experimental facilities . This relates to the charm production cross section in $`pN`$ and $`\pi N`$ reactions as well as to their interactions with baryons and mesons which determine their properties (spectral functions) in the hadronic medium. The charm quark degrees of freedom play an important role especially in the context of a phase transition to the quark-gluon plasma (QGP)  where charmonium ($`c\overline{c}`$) states should no longer be formed due to color screening . However, the suppression of $`J/\mathrm{\Psi }`$ and $`\mathrm{\Psi }^{}`$ mesons in the high density phase of nucleus-nucleus collisions  might also be attributed to inelastic comover scattering (cf.  and Refs. therein) provided that the corresponding $`J/\mathrm{\Psi }`$-hadron cross sections are in the order of a few mb . Present theoretical estimates here differ by more than an order of magnitude  especially with respect to $`J/\mathrm{\Psi }`$-meson scattering such that the question of charmonium suppression is not yet settled. On the other hand, the enhancement of ’intermediate-mass dileptons’ in $`Pb+Pb`$ collisions at the SPS has been tentatively attributed to an enhancement of ’open charm’ in nucleus-nucleus collisions relative to $`pA`$ reactions at the same invariant energy $`\sqrt{s}`$ . Thus ’charmonium suppression’ and ’open charm enhancement’ are present facets of relativistic heavy-ion collisions, which provide a theoretical  and experimental challenge for the future .

In this letter we will argue that the open charm excess reported in Ref.  might be due to secondary meson-baryon (quark-diquark) interactions. We will employ the HSD transport approach  for the overall reaction dynamics of nucleus-nucleus collisions using parametrizations for the elementary production channels for the charmed hadrons $`D,\overline{D},D^{},\overline{D}^{},D_s,\overline{D}_s,D_s^{},\overline{D}_s^{},`$ $`J/\mathrm{\Psi },\mathrm{\Psi }(2S),\chi _{2c}`$ from $`NN`$ and $`\pi N`$ collisions. We recall that in the HSD transport approach the initial stages of a $`\pi N`$, $`pp`$, $`pA`$ or $`AA`$ reaction (at high energy) are described by the excitation of color neutral strings, where the leading quarks and ’diquarks’ in a baryonic string (or quarks and antiquarks in a mesonic string etc.) are allowed to rescatter again (in case of nuclear targets) with hadronic cross sections divided by the number of constituent quarks and antiquarks in the hadrons, respectively .

In order to explore the kinematical situation for secondary ’meson’-’baryon’ interactions we show in Fig. 1 the distribution in the invariant collision energy $`\sqrt{s}`$ from such secondary interactions (solid histogram) for a central collision of $`Pb+Pb`$ at 160 A$``$GeV from the HSD approach. This distribution extends to about 15 GeV while the corresponding primary baryon-baryon collisions are centered around 17.3 GeV (dashed histogram) and become less frequent than the ’meson’-’baryon’ collisions for $`\sqrt{s}`$ 13 GeV. Since the threshold for open charm production in meson-baryon collisions is $`\sqrt{s_0}=m_D+m_{\mathrm{\Lambda }_c}`$ 4.15 GeV there are a lot of secondary collisions above threshold. The crucial question now is the magnitude of the open charm cross section in ’meson’-’baryon interactions for 4.15 GeV $`\sqrt{s}`$ 15 GeV, i.e from threshold up to roughly 10 GeV above threshold.

As argued in Refs.  the creation of a $`c\overline{c}`$ pair at high $`\sqrt{s}`$ is due to a hard process and dominated by gluon-gluon fusion (cf. Fig. 2a) for an illustration of a $`pp`$ reaction). On the other hand, the quark annihilation mechanism is found to contribute significantly or even to become dominant for $`\pi N`$ reactions at lower $`\sqrt{s}`$ . This mechanism is depicted in Fig. 2b) for the particular final state $`\overline{D}\mathrm{\Lambda }_c`$, where – by convention – $`\overline{D}`$ stands for a $`D`$-meson with a $`\overline{c}`$-quark. However, the perturbative $`q\overline{q}1\mathrm{gluon}\mathrm{c}\overline{\mathrm{c}}`$ process is included in standard PYTHIA calculations  and drops very fast with decreasing $`\sqrt{s}`$ close to the threshold of $`D`$meson production . We recall that in Ref.  the open charm cross sections for $`pN`$ and $`\pi N`$ reactions have been calculated within PYTHIA using MRS G (next to leading order) structure functions from the PDFLIB package  for the gluon distribution of the proton, a bare charm quark mass $`m_c=`$ 1.5 GeV and $`k_T`$ = 1 GeV. The results from PYTHIA (scaled to the available data ) for $`D,\overline{D},D^{},\overline{D}^{},D_s,\overline{D}_s,D_s^{},\overline{D}_s^{},`$ as a function of $`\sqrt{s}`$ 10 GeV have been parametrized by an expression of the form,

$$\sigma _X(s)=a_X(1Z)^\alpha Z^\beta ,$$
(1)
with $`Z=\sqrt{s_X^0}/\sqrt{s}`$ where $`\sqrt{s_X^0}`$ denotes the threshold for the channel $`X`$ in $`pN`$ or $`\pi N`$ reactions. The individual parameters $`a_X,\alpha `$ and $`\beta `$ from this fit are given in Table 1 for $`\pi N`$ reactions. As found in Ref.  the contribution to open charm production from ’secondary’ interactions – employing the parametrizations (1) – are about $``$ 9 % in central $`Au+Au`$ collisions at 160 A$``$GeV. This order of magnitude for the open charm enchancement is far below the observation of Ref. .

On the other hand, close to threshold, one would expect the exclusive reaction $`\pi N\overline{D}\mathrm{\Lambda }_c`$ to dominate. In case of $`S`$ -wave production this cross section then should rise as $`(1Z)^{0.5}`$ for a two-body final channel. Such a behaviour is well known experimentally  for the reactions $`\pi ^+p\rho ^+p`$, $`\pi ^{}p\omega n`$, $`\pi ^{}pK^0\mathrm{\Lambda }`$ or $`\pi ^{}p\varphi n`$ as displayed in Fig. 3. These experimental cross sections can be described by the expression

$$\sigma _{\pi NXB}(\sqrt{s})=a_X(1\frac{\sqrt{s_0}}{\sqrt{s}})^{0.5}(\frac{\sqrt{s_0}}{\sqrt{s}})^{\gamma _X}$$
(2)
with $`a_\rho =15`$ mb, $`a_\omega =12`$ mb, $`a_{K^0}=2.5`$ mb, $`a_\varphi =`$ 0.2 mb and $`\gamma _\omega =\gamma _K`$ = 6, $`\gamma _\rho `$ = 5 and $`\gamma _\varphi `$ = 10, respectively. In Eq. (2) $`\sqrt{s_0}`$ denotes the threshold for each reaction separately. For $`\sqrt{s}\sqrt{s_0}m_\pi `$ three-body final states (with an additional pion) become possible while multi-particle production dominates at high invariant energy. To demonstrate the relative importance of multi-particle final states the inclusive cross sections for $`\rho ^+,\omega `$ and $`\varphi `$ are shown in terms of the thick lines in Fig. 3 using the parametrizations from the review . The inclusive yield for all open charm mesons with a $`\overline{c}`$ quark is shown in terms of the lower thick solid line in Fig. 3 where the parametrisations (1) have been employed with the parameters from Table 1 (cf. ). The related data (full squares) have been taken from the review of Tavernier . We recall that for $`\sqrt{s}\sqrt{s_0}`$ 6 GeV the parametrization reflects the PYTHIA results (scaled to the available data) and at lower energies lacks a more fundamental justification as pointed out in Ref. .

Since the binary exclusive reactions dominate the threshold behaviour for $`\rho ,\omega ,K^0`$ and $`\varphi `$ production, one can expect a similar relation to hold also for the channels $`\overline{D}\mathrm{\Lambda }_c`$, $`\overline{D}^{}\mathrm{\Lambda }_c`$, $`\overline{D}\mathrm{\Sigma }_c`$, $`\overline{D}^{}\mathrm{\Sigma }_c,etc`$. Here the reaction mechanism is expected to be dominated by $`D^{}`$-meson or $`D^{}`$-Reggeon exchange as illustrated in Fig. 2c). The question that remains is the size of a related parameter $`a_D`$ in (2). Here we refer to the Regge model for an estimate.

In Ref.  the cross section for the reaction $`\pi N\overline{D}(\overline{D}^{})\mathrm{\Lambda }_c`$ has been estimated within the framework of the Quark-Gluon String Model (QGSM) developed in Ref. . The GGSM is a nonperturbative approach based on a topological $`1/N_f`$ expansion in QCD, where $`N_f`$ denotes the number of flavors, and on the color-tube model. This approach can be considered as a microscopic model describing Regge phenomenology in terms of quark degrees of freedom. It provides the possibility to establish relations between many soft hadronic reactions as well as masses and partial widths of resonances with different quark contents (see e.g. the review ). Recently, this model also has been successfully applied to the description of the nucleon and pion electromagnetic form factors  and deuteron photodisintegration .

The amplitude of the reaction $`\pi N\overline{D}(\overline{D}^{})\mathrm{\Lambda }_c`$ corresponding to the planar graph of Fig. 2c) with $`u`$ and $`\overline{c}`$ quark exchange in the $`t`$ channel – using the Mandelstam variables $`s`$ and $`t`$ – can be written as (see Ref. )

$$A_{\pi ^{}p\overline{D}\mathrm{\Lambda }_c}(s,t)=\sqrt{2}g_0^2F(t)(s/s_0^{u\overline{c}})^{\alpha _{u\overline{c}}(t)1}(s/\overline{s}),$$
(3)
where $`g_05.8`$ is a universal flavor independent coupling constant, $`\alpha _{u\overline{c}}(t)=\alpha _D^{}(t)`$ is the $`D^{}`$ Regge trajectory, $`\overline{s}=1`$ GeV2, $`s_0^{u\overline{c}}`$= 4.9 GeV2 is the flavor dependent scale factor and $`F(t)`$ is the form factor describing the $`t`$ dependence of the residue. We have considered two forms of the $`D^{}`$ Regge trajectory, i.e.
i) a linear trajectory (set $`a`$ from Ref. )

$$\alpha _D^{}(t)=0.86+0.5t$$
(4)
ii)a nonlinear trajectory, i.e. the square root trajectory from Ref.

$$\alpha _D^{}(t)=\alpha _D^{}(0)\gamma [\sqrt{T}\sqrt{Tt}]$$
(5)
with $`\alpha _D^{}(0)=1.02`$, $`\gamma =3.65`$ GeV-1, $`\sqrt{T}=`$3.91 GeV. Both choices lead to very similar total cross sections such that we do not discuss the two choices separately. The form factor $`F(t)`$ determining the $`t`$-dependence of the residue has been parametrized in Ref.  as

$$F(t)=\mathrm{\Gamma }(1\alpha _D{}_{}{}^{}(t))$$
(6)
and is motivated by dual models. The form factor (6) is convenient for an analytical continuation of the amplitude to the region of positive $`t`$, where the $`\mathrm{\Gamma }`$ function in (6) is an exponentially decreasing function of $`t`$. In particular, the form (6) was successfully used to relate the normalization of different planar graphs at $`t`$= 0 with the widths of the mesons lying on the corresponding Regge trajectories. For example, in Ref.  a width $`\mathrm{\Gamma }(D^+D^0\pi ^++D^+\pi ^0)=(3/2)\mathrm{\Gamma }(D^+D^0\pi ^+)`$= 20 keV was found. This value is not very different from the result based on the QCD sum rules  giving $`48`$ keV, while the experimental upper limit is 131 keV.

However, in the region of negative $`t`$ the parametrization (6) exhibits a factorial growth (which is faster than exponential) and therefore is not acceptable. For $`t0`$ we will use the convential parametrization

$$F(t)=\mathrm{\Gamma }(1\alpha _D{}_{}{}^{}(0))\mathrm{e}xp(R^2t)$$
(7)
with $`R^2=0÷0.2`$ GeV-2.

The differential cross section for the reaction $`\pi ^{}p\overline{D}\mathrm{\Lambda }_c`$ then is

$`{\displaystyle \frac{d\sigma _{\pi ^{}p\overline{D}\mathrm{\Lambda }_c}}{dt}}`$ $`={\displaystyle \frac{1}{64\pi s}}{\displaystyle \frac{1}{(p_\pi ^{\mathrm{cm}})^2}}|A_{\pi ^{}p\overline{D}\mathrm{\Lambda }_c}(s,t)|^2,`$ (8)
where $`p_\pi ^{\mathrm{cm}}`$ denotes the pion momentum in the cms, and the total cross section is obtained by integrating (8) over the kinematical allowed regime.

The result for the channel $`\pi ^{}p\overline{D}\mathrm{\Lambda }_c`$ (using $`R^2=0`$) is shown in Fig. 3 by the lower solid line as a function of $`\sqrt{s}\sqrt{s_0}`$, which can be well described by the function (2) using $`a_D`$ = 0.027 mb and $`\gamma _D=\gamma _\omega =\gamma _K`$ = 6 (lower dashed line). The ratio of the $`\omega `$ cross section to the ’estimated’ $`\overline{D}`$ cross section close to threshold ( $`\sqrt{s}\sqrt{s_0}0.4`$ GeV) is $``$ 500; a similar ratio holds for the inclusive cross sections at $`\sqrt{s}\sqrt{s_0}30`$ GeV where explicit data are available for the $`\pi ND\overline{D}+X`$ reaction. Thus our estimate in the Regge model sounds reasonable, however, has to be controlled by experiment in order to allow for definite conclusions.

We mention that similar exclusive cross sections can also be obtained from boson-exchange models with $`D^{}`$ or $`D_1`$ exchange (cf. Fig. 2c) employing monopole form factors at the vertices with cutoffs $`\mathrm{\Lambda }(1.61.9`$ GeV) while fixing the couplings $`g_{\overline{D}\pi \overline{D}^{}}`$ and $`g_{\overline{D}^{}\pi =D_1}`$ to the $`\pi \overline{D}`$ and $`\pi \overline{D}^{}`$ decay widths, respectively. Furthermore, the cross section for the reaction $`NN\overline{D}(D^{})\mathrm{\Lambda }_cN`$ has been calculated within the Reggezied one-pion exchange model using the pion form factor $`F_\pi =\mathrm{exp}(R_\pi ^2t)`$ with $`R_\pi ^2=0.1(GeV/c)^2`$ as in Ref. . Since the explicit cross sections for the $`NN`$ reactions at low $`\sqrt{s}`$ are of no major importance for our present study, we refer the reader to a forthcoming publication on these issues .

For the transport calculations (to be discussed below) we estimate the cross section for charmed hadron production in $`\pi N`$ reactions by employing the ansatz (2) for all charmed hadrons such as $`\overline{D},\overline{D}^{},\overline{D}_s,\overline{D}_s^{}`$ with associated charmed hyperons $`\mathrm{\Lambda }_c,\mathrm{\Sigma }_c,\mathrm{\Lambda }_c^s...`$. We explicitly adopt $`a_D=1/3a_D^{}=a_{D_s^{}}=3a_{D_s}`$ (as in Ref. ) with $`a_D0.027`$ mb. These estimates, of cause, only have exploratory character as stated above.

We now turn to the results of transport calculations for the reaction $`Pb+Pb`$ at 160 A$``$GeV. A description of the transport approach is given in Refs.  and the production and propagation of charmed mesons is described in Refs. . The novel phenomenon addressed here with respect to Ref.  is the secondary production of open charm mesons by ’meson’-’baryon’ collisions with the new parametrizations (2) for the low-energy $`\pi N\overline{D}(\overline{D}^{})\mathrm{\Lambda }_c(\mathrm{\Sigma }_c)`$ processes. We mention in passing that the role of secondary reactions on intermediate mass dilepton pairs via the $`q\overline{q}\gamma l^+l^{}`$ (Drell-Yan) mechanism has been investigated in the UrQMD transport model in Ref.  for nucleus-nucleus collisions at SPS energies before.

In Fig. 4 we show the open charm multiplicity (all $`D`$, $`\overline{D},D^{}`$, $`\overline{D}^{},D_s^{}`$, $`\overline{D}_s^{},D_s,\overline{D}_s`$) as a function of the impact parameter $`b`$ from the HSD approach for $`Pb+Pb`$ at 160 A$``$GeV. Here the dashed line gives the yield from primary baryon-baryon collisions, which is essentially the yield from $`pn`$ reactions at $`\sqrt{s}`$ 17 GeV times the number of ’hard’ $`pN`$ collisions (as given by Eq. (6) of Ref. ), while the solid line corresponds to the yield of open charm from secondary reactions as described above. Whereas for peripheral collisions the secondary production is slightly smaller than the direct channel, the yield from secondary meson-baryon channels becomes larger with decreasing impact parameter $`b`$. The tiny kinks in the curves in Fig. 4 are due to the limited statistics of the transport calculations. Without explicit representation we note that the final differential spectra of the $`D`$-mesons in transverse momentum and rapidity are very similar to those from primary and secondary channels, respectively, due to rescattering with the surrounding hadrons (cf. Ref. ).

In order to compare to the data from the NA50 Collaboration  we have to adopt a model to convert the number of participating nucleons $`A_{part}`$ to the impact parameter $`b`$ from the transport calculation. We use

$$A_{part}=2AN_0(b),$$
(9)
where $`N_0(b)`$ stands for the number of noninteracting nucleons (with no ’hard’ collisions) from the transport code at impact parameter $`b`$ while $`A`$ is the target (projectile) mass number. Denoting the expected number of $`D,\overline{D}`$ mesons from primary $`BB`$ collisions by $`N_{prim.}`$ and the number from secondary ’meson’-’baryon’ interactions by $`N_{sec.}`$ this leads to the ratio

$$R(b)=\frac{N_{prim.}+N_{sec.}}{N_{prim.}}(b)=R(A_{part}(b))=R(A_{part}),$$
(10)
which is displayed in Fig. 5 (solid histogram) in comparison to the data on the open charm enhancement from Ref. . The calculated ratio first increases fast with $`A_{part}`$ and becomes almost constant for $`A_{part}100`$. This is due to the fact that secondary semihard ’meson’-’baryon’ reactions with ’wounded’ nucleons, that have scattered at least once, already set in for large impact parameter $`b`$; their relative number then increases only slightly faster than $`A_{part}`$. Obviously, the general trend of the data can be roughly described with increasing centrality for the secondary cross sections specified above. The present statistical uncertainty of the data  does not allow for a final conclusion.

In summary: in this letter we have argued that for ’low’ energy ’meson’-’baryon’ reactions the dominant $`c\overline{c}`$ production is related to the two-body (or quasi two-body) reactions $`\pi N\overline{D}(\overline{D^{}})`$ $`\mathrm{\Lambda }_c(\mathrm{\Sigma }_c)`$. Estimates within the framework of the Quark-Gluon String model suggest cross sections of a few $`\mu b`$ for $`\pi N\overline{D}\mathrm{\Lambda }_c`$ in the region of about 1 GeV above threshold. The estimated order of magnitude for the open charm cross section (cf. lower part of Fig. 3) is found to be compatible with the ’open charm enhancement’ claimed by the NA50 Collaboration at the SPS  without employing the assumption of thermal and chemical or statistical equilibrium as advocated in Refs. . It should be stressed, however, that experimental investigations on open charm production in $`\pi N`$ reactions at invariant energies of 4.2 $`\sqrt{s}`$ 15 GeV are mandatory to confirm or disprove our suggestion.

We, furthermore, note that the cross section for charmonia such as $`J/\mathrm{\Psi }`$, $`\chi _c`$ or $`\mathrm{\Psi }^{}`$ are not substantially enhanced by such secondary reaction channels since their cross sections are small in $`\pi N`$ collisions  such that no substantial ’enhancement’ of charmonia relative to the primary $`NN`$ reaction channels is expected.

The authors are grateful to E.L. Bratkovskaya, A.B. Kaidalov and R. Vogt for clarifying discussions.