text:
# The temperatures of dust-enshrouded active galactic nuclei

## 1 Introduction

There is increasing evidence (Magorrian et al. 1998; van der Marel 1999, Ferrarese & Merritt 2000, Gebhardt et al. 2000) for massive dark objects (MDOs), probably supermassive black holes (Kormendy & Bender 1999), in the centers of most large nearby galaxies. The comoving mass density of these MDOs is higher by at least a factor of two (Salucci et al. 1999), and possibly more (Phinney 1997; Haehnelt, Natarajan & Rees 1998), than can be explained by accretion onto optically-luminous AGNs.

Haehnelt et al. (1998) suggested one possible solution to this problem: that much of the accretion is obscured by dust, and that this population of dust-obscured AGNs is missing from optically selected samples (see also Fabian 1999 for a specific model). Dust-enshrouded AGNs are certainly known at high redshift: see Table 1 for some examples, and McMahon et al. (1999) and Benford et al. (1999). However, note that the fraction of their bolometric luminosity that originates from an AGN as opposed to a starburst has not been well established in most cases. In an important recent development, analysis of the hard X-ray background (Fabian & Iwasawa 1999; Salucci et al. 1999) has revealed that if the sources responsible for producing the X-ray background are obscured (type-2) quasars, then they may also be the class of sources that are responsible for producing the high local MDO density, as dust and the associated density of molecular and neutral gas is transparent to hard X-rays.

Fabian & Iwasawa (1999) suggest that the bulk of the energy absorbed at optical, ultraviolet and soft X-ray wavelengths is reradiated at far-infrared and submillimetre wavelengths by obscured AGNs. We investigate this scenario using constraints from far-infrared and submillimetre source counts and backgrounds. If we assume that interstellar dust is responsible for the absorption and reradiation, then the question becomes: what is the dust temperature of these sources?

Almaini et al. (1999) and Gunn & Shanks (2000) recently proposed that the obscured quasars contributing to the X-ray background have cold ($`40`$ K) dust shrouds, and that they account for about ten to twenty per cent of the mJy sources detected at 850 $`\mu `$m using the SCUBA camera at the JCMT (Blain et al. 2000 and references therein). Their motivation for suggesting this temperature comes from observations of high-redshift radio-quiet AGNs (Benford et al. 1999). Interestingly this is the most favoured temperature for the SCUBA sources derived from current measurements of the far-infrared and submillimetre background, source counts and redshift distribution (Blain et al. 1999a,b; Trentham, Blain & Goldader 1999).

This idea is very attractive. The dust in type-2 AGNs, which are required to produce the hard X-ray background, is observed directly by SCUBA reradiating away the absorbed ultraviolet and soft X-ray energy in the submillimetre waveband. Two of these SCUBA sources, in the field of the lensing cluster Abell 370 (Smail, Ivison & Blain 1997), which are known from optical spectroscopy to contain AGNs, have in fact been seen at X-ray wavelengths with Chandra by Bautz et al. (2000). On the other hand, much of the accretion may take place surrounded by hotter ($`60`$ K) dust shrouds, as would be required by radiative transfer models of the spectral energy distributions (SEDs) of luminous infrared galaxies (Rowan-Robinson 2000 and references therein). These objects would then on the whole have SEDs skewed to wavelengths too short to be seen easily with SCUBA, yet too long to have been seen in ground-based near-infrared surveys. They may be present in substantial numbers but until recently have been very difficult to detect at all. The first such objects to be seen in a systematic survey were recently described by Wilman, Fabian & Gandhi, who found $`ISO`$ 6.7-$`\mu `$m and 15-$`\mu `$m counterparts of two hard X-ray sources detected by Chandra. Many local infrared-luminous AGNs, for example the ‘Warm ULIGs’ of Sanders & Mirabel (1996), have temperatures close to 60 K. Higher temperatures in excess of 100 K are observed for a few, rare high-redshift ULIGs, such as P09104+4109 (Kleinmann et al. 1988) and APM 08279+5255 (Lewis et al. 1998), with temperatures in excess of 100 K. The fortunate circumstances under which APM 08279+5255 (Irwin et al. 1998) was discovered (the infrared source happens to lie close to the caustic of a lensing galaxy so is highly magnified) are perhaps suggestive that many such objects (without such a large magnification) exist but are as yet unidentified. Even higher dust temperatures (150–200 K) are implied by the observed bumps in the SEDs of the most luminous Seyfert galaxies at about 25 $`\mu `$m (Sanders & Mirabel 1996), although the bolometric luminosities of these galaxies are not dominated by these hot components.

Therefore it is worth investigating the constraints on the temperatures of dust-enshrouded AGNs, the observations of which are described in Section 2. In Section 3 we demonstrate the effect of the dust temperature on both the observable source counts and the efficiency of accretion. We conduct more detailed calculations in Section 4. In Section 5 we discuss our results and indicate the most promising lines of enquiry for the future. Throughout this paper we assume that $`H_0=50`$ km s-1 Mpc-1 and $`\mathrm{\Omega }_0=1`$.

## 2 Dust-enshrouded AGNs

It is plausible that dust-enshrouded AGNs are responsible for the production of the hard X-ray background radiation and constitute a subset of the SCUBA galaxies. There is an additional argument to support the presence of obscured AGNs. If the SCUBA sources have similar gas density structures as local ultraluminous infrared galaxies (ULIGs), specifically that the bolometric luminosity originates from a very dense compact central region (Downes & Solomon 1998; Sakamoto et al. 1999), as suggested by their similar SEDs (Iviso et al. 1998, 2000), then a very efficient power source is required. The only $`z=0`$ stellar systems that exist at the same densities, 100 M pc-3 or greater, are the cores of elliptical galaxies (Kormendy & Sanders 1992); however, because only about one per cent of the cosmological density in stars (Fukugita, Hogan & Peebles 1998) is in these systems, it is unlikely that their formation alone could have generated the large fraction of the bolometric energy density of the Universe ($`>10`$ per cent) that is accounted for by the SCUBA sources (Trentham 2000). High-mass stars would be another alternative, but this would then require the MDOs to be clusters of stellar remnants, not black holes, an idea which seems to be ruled out given the high angular resolution of recent stellar-kinematical measurements of the centers of nearby galaxies (e.g. Kormendy & Bender 1999) unless the clusters of stellar remnants coalesce to form supermassive black holes.

An alternative efficient power source would be accreting AGNs. However, it is unlikely that sources with bolometric luminosities that originate from these compact central regions, with cold rotating nuclear gas disks on scales $`<300`$ pc, can be powered by AGNs alone, if their only source of fuel is these gas disks. This is because it is difficult to remove sufficient angular momentum from the gas to allow it to reach the event horizon of a black hole without converting the gas into stars (Begelman 1994). This is a long-standing astrophysical problem, and it is now confirmed by observations on scales where the central black hole contributes siginificantly to the local gravitational potential (e.g. Sakamoto et al. 1999).

## 3 The observational effects of different dust temperatures in accreting AGNs

A luminosity function can be adopted to describe the population of submillimetre-luminous galaxies, for example the Gaussian function chosen by Trentham et al. (1999). This can be integrated over redshift and luminosity to make predictions of the source counts and background radiation intensities in the submillimetre and far-infrared wavebands, if an (SED) is assumed, (which is normalised to the spectral luminosity at an arbitrary wavelength of 60 $`\mu `$m, $`L_{60}`$, chosen for convenient comparison with IRAS studies of the local Universe. This luminosity function is consistent with the observed source counts of galaxies at 850 and 175 $`\mu `$m and the intensity of the far-infrared background radiation (Puget et al. 1996, Hauser et al. 1998; Schlegel et al. 1998, Finkbeiner et al. 2000) if the dust temperature $`T40`$ K, assuming that the dusty soures all radiate as modified blackbodies with a single temperature, and that the dust emissivity varies with frequency as a power law $`\nu ^p`$ where $`p=1.5`$.

In Fig. 1, we show the effects of changing the dust temperature, and thus modifying the SED, on the counts of galaxies at 850 and 175 $`\mu `$m. The results are presented in the form of ratios as compared with the values at $`T=40`$ K, and so details of the normalization of the counts and backgrounds need not be considered here. The Gaussian luminosity function used has a characteristic 60-$`\mu `$m luminosity $`L_{60}^{}=3\times 10^{13}`$ L and a width $`\sigma =0.5`$. No luminosity evolution is assumed and the comoving density of sources is assumed to be constant in the redshift range $`1<z<3`$ and zero elsewhere (Trentham et al. 1999).

At 850 $`\mu `$m the counts are less at greater dust temperatures, because the hotter sources are being probed further down the Rayleigh–Jeans tail of their SEDs, and so only sources with greater bolometric luminosity can exceed the flux density limit required to appear in the count. The temperature dependence of the 175-$`\mu `$m source counts, and of the background radiation intensity at 140 and 240 $`\mu `$m, is much weaker as these quantities are derived from the properties of the SED much closer to its peak.

If we assume that the bolometric far-infrared luminosity of the AGN, $`L_{\mathrm{bol}}`$, is powered entirely by accretion onto a black hole at a fraction $`f_{\mathrm{edd}}`$ of the Eddington luminosity and that the accretion efficiency is independent of luminosity, then the mass of the central black hole (the MDO) is

$$\frac{M}{\mathrm{M}_{}}=3.5\times 10^5f_{\mathrm{edd}}^{}{}_{}{}^{1}\left(\frac{L_{\mathrm{bol}}}{\mathrm{L}_{}}\right).$$
(1)
In this accretion model, a greater bolometric luminosity must derive from a greater accretion rate, and so correspond to a greater black hole mass. As shown in Fig. 1, the comoving MDO density generated in order to power the 850- and 175-$`\mu `$m sources depends on the dust temperature, in the sense that very cold and very hot sources have greater bolometric luminosities per unit spectral luminosity at the wavelength of 60 $`\mu `$m, at which the luminosity function is defined, or indeed per unit bolometric luminosity.

It is clear from Fig. 1 that a greater comoving MDO density is produced by hotter sources per unit contribution to the 850-$`\mu `$m source counts. If we assume that the bolometric far-infrared luminosities of the sources are proportional to their absorbed soft X-ray flux densities, that is the fraction of the energy absorbed by dust and gas at wavelengths outside the soft X-ray window is either negligible or independent of luminosity, then the absolute value of the hard X-ray background intensity is almost proportional to the comoving MDO density. The difficulties in normalizing this relation are discussed in Section 4. The dust temperature could thus have a significant effect on the submillimetre/far-infrared observability of the sources responsible for the generation of the hard X-ray background. Because the numbers plotted in Fig. 1 are ratios, when the models are matched to observations in detail, the absolute values of these (and other) quantities must also be considered as a function of temperature.

## 4 Specific normalized models of accretion

We now investigate the normalization of specific models of dust-enshrouded accretion onto black holes at high redshift. As discussed above, we can integrate over redshift and luminosity to derive the comoving mass density of the black holes that are accreting at any redshift $`z`$ using equation (1). By including the density of black holes that have accreted in the past, but are quiescent at $`z`$, we can compute the total mass density of MDOs at $`z`$.

We can constrain the redshift-dependant luminosity function of $`\varphi (L_{60},z)`$ of accreting sources, their temperatures $`T`$, and $`f_{\mathrm{edd}}`$ in combination if we know the mass density of MDOs at $`z=0`$. Further constraints are imposed by measurements of the counts and backgrounds if we assume an SED for the sources, and by the X-ray background if we include some infrared-to-X-ray conversion factor. A number of assumptions must first be made: 1) All sources are assumed to have the same temperature. In reality sources probably have a range of temperatures; consequently $`T`$ should be regarded as a luminosity-weighted average temperature. 2) A specific form of $`\varphi (L_{60},z)`$ must be assumed. The dependence on $`L_{60}`$ is assumed to be a Gaussian function with width $`\sigma =0.5`$ and characteristic luminosity $`L_{60}^{}`$. No luminosity evolution is assumed, that is $`L_{60}^{}`$ does not vary with redshift, and the evolution of the comoving space density of sources is chosen to ensure that the integral $`_0^{\mathrm{}}L_{60}\varphi (L_{60},z)dL`$ varies with redshift in the same way as the ultraviolet AGN luminosity density (Boyle & Terlevich 1998). 3) Any mode of accretion that is either not dust-enshrouded or not generating powerful emission, for example advection-dominated accretion, is assumed to be negligible. 4) When constraining $`\varphi (L,z)`$, $`T`$, and $`f_{\mathrm{edd}}`$ using measurements of the backgrounds and source counts at various wavelengths, we require that the backgrounds or source counts are never overproduced. It is not a concern if the model underpredicts the observed quantities, as there is certainly a significant contribution from infrared-luminous galaxies that are powered by star-formation activity as opposed to AGN accretion. The models are constrained by the upper limit to the 2.8-mm source counts (Wilner & Wright 1997), the 850-$`\mu `$m and 450-$`\mu `$m source counts (Blain et al. 2000), the 175-$`\mu `$m source counts (Kawara et al. 1998; Puget et al. 1999), the submillimetre and far-infrared background radiation intensity (Puget et al. 1996, Fixsen et al. 1998; Hauser et al. 1998; Schlegel et al. 1998, Finkbeiner et al. 2000) and the MDO density at $`z=0`$, $`\rho _{\mathrm{MDO}}(0)`$. The methodology used to make these comparisons is described in detail by Trentham et al. (1999).

## 5 Results and discussion

In all the models, very low dust temperatures are ruled out to avoid breaching both the limit to the observed 2.8-mm counts (Wilner & Wright 1997) and the submillimetre background intensity (Fixsen et al. 1998). For most values of $`L_{60}^{}`$, the most important constraint becomes the requirement that the 175-$`\mu `$m source counts are not overproduced, a constraint which also depends on $`T`$. The 850- and 450-$`\mu `$m counts are exceeded for some values of $`L^{}`$ and $`T`$, but these values are always also ruled out by the observed intensity of either the 850- or 240-$`\mu `$m backgrounds. As the value of $`L^{}`$ increases, the count constraints become steadily more important relative to the background constraints, but never become more significant for the values of $`L_{60}^{}<10^{14}`$ L.

The results can be summarized by a simple equation:

$$\frac{t}{\mathrm{K}}\left[\frac{\rho _{\mathrm{MDO}}(0)}{\mathrm{M}_{}\mathrm{Mpc}^3}f_{\mathrm{edd}}\right]^{0.36}.$$
(2)
This equation is an equality if the submillimetre and far-infrared backgrounds are produced entirely by dust-enshrouded AGNs; it is an inequality if a significant fraction of the background comes from another population of objects.

From equation (2), $`T70`$ K is required for $`f_{\mathrm{edd}}>0.1`$, if we assume the higher value of the local MDO density implied by the MDO mass–galaxy luminosity correlation of Magorrian et al. (1998) i.e. $`\rho _{\mathrm{MDO}}(0)=1.5\times 10^6\mathrm{M}_{}\mathrm{Mpc}^3`$. Such objects are probably not the SCUBA-selected galaxies, which probably have temperatures of about 40 K (Ivison et al. 1998, 2000; Blain et al. 1999b; Trentham et al. 1999). An additional warm population of dusty galaxies must be hypothesized to explain the X-ray background and the local MDO density in this case.

This scenario is interesting, because the observed 140- and 240-$`\mu `$m background radiation intensity can be explained entirely by this population, as shown in Fig. 2. The cooler SCUBA sources in the models of Blain et al. (1999a,b) and Trentham et al. (1999) tend to underpredict the 140-$`\mu `$m background and the 100-$`\mu `$m and 60-$`\mu `$m backgrounds of Finkbeiner et al. (2000). Hence, a hotter additional population could be argued for a priori, and it is attractive that such a population can simultaneously provide the accretion necessary to generate the MDO density at $`z=0`$.

If, on the other hand, we assume the lower value of the local MDO density implied by van der Marel (1999) i.e. $`\rho _{\mathrm{MDO}}(0)=2\times 10^5\mathrm{M}_{}\mathrm{Mpc}^3`$, and that MDOs formed in only the dust-enshrouded high-redshift AGNs that contribute to the hard X-ray background, then sources with temperatures as low as 30 K are consistent with the observations. Such objects could easily make up a subset of the SCUBA sources. This is essentially the scenario proposed by Almaini et al. (1999) and studied by Gunn & Shanks (2000).

### 5.1 Constraints from X-ray observations

It is difficult to use measurements of the hard X-ray background as a strong constraint on the history of accretion, because we do not know the SED of a typical obscured AGN over the entire spectrum from submillimetre to hard X-ray wavelengths. For example, the dust-enshrouded AGNs listed in Table 1 have very different SEDs.

Fabian & Iwasawa (1999) estimated that the total X-ray flux density that is absorbed and reradiated by the dust-enshrouded AGN population is about 3 nW m-2 sr-1. This is the integrated value over all wavelengths and is thus independent of the dust temperature. This is 58, 7, and 7 per cent of the background radiation integrated over all wavelengths for the three solid lines on Fig. 2 in order of increasing thickness, and so in these three models, 42, 93, and 93 per cent repsectively of the energy absorbed by the dust shrouds of the AGNs must come from either a population of Compton-thick sources that are not present in X-ray surveys, or from outside the X-ray window defined by Fabian & Iwasawa (1999).

One possibility is that very luminous heavily obscured AGNs, perhaps like Markarian 273 (that looks like a starburst at near-infrared and optical wavelengths but like a dusty AGN at mid-infrared wavelengths; Lutz et al. 1999), or maybe even more obscured, are common at high redshift. These could well be Compton-thick. This possibility is suggested by the non-detection of six SCUBA sources with the $`Chandra`$ X-ray satellite (Fabian et al. 2000; see Table 2 and the accompanying list of energy budgets). Such a scenario also would be suggested by the spectroscopy of the SCUBA sources by Barger et al. (1999), who find no Seyfert 1’s like Markarian 231, so that if their sources are AGN-powered they must be heavily obscured.

Table 2 should be considered in conjunction with the following energy contributions to the extragalactic background light. $`ϵ_{\mathrm{cc}}+ϵ_{\mathrm{hc}}=2\mathrm{nW}\mathrm{m}^2\mathrm{sr}^1`$ from Fabian & Iwasawa (1999), after removing their correction for Compton-thick sources. $`ϵ_{\mathrm{cc}}+ϵ_{\mathrm{cC}}7\mathrm{nW}\mathrm{m}^2\mathrm{sr}^1`$ from Trentham et al. (1999), summing all the energy in the SCUBA sources, the temperatures of which are constrained to not overproduce the far-infrared background at the high-$`T`$ end, and to produce about one-third of the 450-$`\mu `$m background at the low-$`T`$ end (Blain et al. 1999a). This is a lower limit due to the possibility of a contribution from sources powered by star formation. $`ϵ_{\mathrm{cc}}+ϵ_{\mathrm{hc}}+ϵ_{\mathrm{cC}}+ϵ_{\mathrm{hC}}=5\mathrm{nW}\mathrm{m}^2\mathrm{sr}^1`$ for the thin solid line in Fig. 2.

The possibility of there existing many hot Compton-thick sources is also realistic as both more cold sources appear to be Compton-thick than Compton-thin (see the first two lines of Table 2), and more Compton-thin sources appear to be hot than cold (see the second and fourth lines of Table 2). Note that caution is required when directly interpreting the numbers in the fourth, fifth and sixth columns of Table 2 since the sensitivities of Chandra and SCUBA to the same template object in each survey are different). Such sources may appear at hard X-ray wavelengths, like those probed more sensitively by the XMM-Newton satellite.

If the entire observed intensity of the 140-$`\mu `$m background radiation, about 30 nW m-2 sr-1, ten times greater than Fabian & Iwasawa (1999) estimated, is to be explained by dust-enshrouded AGNs, then higher-temperature sources radiating at a high efficiency are required. However, other mechanisms, for example warm dust reradiating absorbed starlight from normal star-forming galaxies at $`z2`$ (see Appendix A), or an isotropic warm Galactic dust component (Lagache et al. 1999), could be responsible for a significant fraction of the 140-$`\mu `$m background, thus reducing the background radiation intensity that is required from this additional AGN population.

In summary, we cannot distinguish between these scenarios at present. The properties of high-redshift dust-enshrouded AGNs remain only poorly constrained. Nevertheless using the results presented here we can suggest productive lines of enquiry for the future.

### 5.2 Making observational progress in understanding the growth of MDOs and dust-enshrouded AGNs

More information about both the properties of distant dusty galaxies and the density of MDOs at low redshifts is required to better understand the dominant mode of AGN accretion at high redshifts. Based on the models preesented here, we now list the most promising directions for future progress. (1) Multiwaveband follow-up observations of high-redshift far-infrared and submillimetre selected galaxies will be required to identify candidate dusty AGNs that are undergoing significant amounts of accretion. Optical and near-infrared spectroscopy is currently available for some SCUBA sources (Ivison et al. 1998, 2000; Barger et al. 1999). This can be used to provide some evidence for the presence or absence of an AGN, but not to unambiguously determine the nature of the power source. The likely extreme optical faintness of some SCUBA galaxies (Barger, Cowie & Richards 2000; Smail et al. 2000) will make spectroscopic observations of complete samples very difficult. Mid-infrared measurements of PAH features (for example Genzel et al. 1998) and the shapes and ratios of near-infrared emission lines (for example Veilleux et al. 1999) may allow the unambiguous diagnosis of AGNs using future instruments like $`SIRTF`$, $`FIRST`$ and $`SPECS`$. (2) The best wavelength at which to make an efficient search for more examples of dust-obscured AGNs depends on both the dust temperatures of the objects, and the sensitivity of detectors. The SED of an object with a dust temperature $`T`$ at redshift $`z`$ peaks at a wavelength of about $`65(1+z)(T/40\mathrm{K})^1`$$`\mu `$m in the observers frame. Hence, cooler objects are probably best searched for at submillimetre wavelengths, while warmer objects will probably be detected more efficiently in far-infrared surveys. The properties of galaxies detected using SCUBA at 850-$`\mu `$m, ISO at 175-$`\mu `$m, and IRAM at 1.25 mm (e.g. Bertoldi et al. 2000) are already revealing the best strategies for future surveys (Blain 1999). (3) Improved modeling of the velocity dispersion anisotropy in observations of the stellar kinematics in nearby galaxies (see Section 4.6 of van der Marel 1999) will hopefully allow a more precise determination of MDO masses, and thus of the local black hole density $`\rho _{\mathrm{MDO}}(0)`$. The improved form of the relation between galaxy luminosity and MDO mass that will result may also provide another constraint on the accretion history, since the local early-type galaxy luminosity function is well known (e.g. Binggeli, Sandage & Tammann 1988). Correlating MDO masses with other properties, like velocity dispersion (Ferrarese & Merritt 2000; Gebhardt et al. 2000) will help too, particularly if the correlation is tight. (4) More sensitive measurements of the hard ($``$ 10 keV) X-ray source counts will soon be available from XMM-Newton. By cross correlating sensitive hard-X-ray maps with deep submillimetre and far-infrared observations, the nature of dust-enshrouded AGNs can probably be pieced together slowly on an object by object basis.

Ultimately, these pieces of information will allow a more rigorous calculation similar to that described in Section 4 to carried out. This will allow the nature of the dust-enshrouded AGN population to be clarified, providing a strong constraint on $`f_{\mathrm{edd}}`$ and thus telling us about the fundamental physics of accretion.

## 6 Conclusions

It is not yet possible to determine a unique description of the SEDs of the dust-enshrouded AGNs that are hypothesized to be responsible for both the hard X-ray background radiation intensity and the density of massive dark objects (MDOs) in the cores of low-redshift galaxies. If the density of MDOs is at the bottom of the current observational range, then the low-temperature dust-enshrouded AGNs described by Almaini et al. (1999), a subset of the high-redshift SCUBA galaxies, can account for all the accretion. However, a hotter population of dust-enshrouded AGNs is required if the MDO density is towards the top of the range permitted by current observations. This second scenario has the advantage of accounting for most of the far-infrared background radiation. New clues have been provided by the lack of strong cross-correlation between SCUBA and $`Chandra`$ sources, implying that most cold sources, if they are AGNs, are Compton-thick. There are more hot than cold Compton-thin sources, and so one might expect there to be significant numbers of hot Compton-thick sources seen by neither $`Chandra`$ nor SCUBA. A hotter population of dust-enshrouded AGNs does indeed seem to be a realistic proposition. Deeper observations at higher X-ray energies with XMM - Newton will help to uncover such a population, if it exists.

## Acknowledgments

We thank Priya Natarajan, Max Pettini, and Dave Sanders for helpful discussions, and acknowledge PPARC and the Raymond & Beverly Sackler Foundation for financial support.

APPENDIX A An alternative source of the submillimetre and far-infrared background radiation, and one which must contribute to some degree, is the reprocessed radiation from young stars that is absorbed and reradiated by dust in normal star-forming galaxies at redshifts $`0z4`$ (Steidel et al. 1999). The contribution to the background light at four different wavelengths is shown in Fig. 3 as a function of the temperature of the dust that is radiating the absorbed energy. If this temperature is between 50 and 100 K, then much of the 140-$`\mu `$m background radiation intensity can be explained, while there is a negligible contribution to the 850-$`\mu `$m background, which can be produced almost entirely by the known population of more luminous SCUBA sources. Dust temperatures of 50-100 K are similar to those of the dust clouds that dominate the bolometric luminosity of Galactic star-forming regions (Wynn-Williams 1982). For $`T=75`$ K, the 850-$`\mu `$m flux of even a very luminous $`M_{\mathrm{AB}}`$(1700 Å) $`=23`$ Lyman break galaxy at $`z=3`$ would be only about 0.7 mJy, lower than the $`>2`$ mJy fluxes of the SCUBA sources. See also the analyses of Ouchi et al. (1999) and Peacock et al. (2000).


format: text

text:
# Simulations of galaxy formation in a cosmological volume

## 1 Introduction

A detailed understanding of galaxy formation is one of the central goals of contemporary astrophysics. Over the past decade, this goal seems, finally, to have come within reach. On the observational side, data from the Keck and Hubble Space telescopes have revolutionised our view of the high redshift Universe. From the theoretical point of view, it is clear that studying formation involves a synthesis of ideas from a wide range of specialities. A full theoretical treatment requires consideration of the early Universe processes that created primordial density fluctuations, the non-linear dynamics that result in the formation of dark matter halos, the dissipational processes of cooling gas, the microphysics and chemistry that precipitate star formation, the feedback due to energy exchange between supernovae and interstellar gas, the merging of galaxy fragments and the interaction of galaxies with their large-scale environment.

Because of its non-linear character, lack of symmetry and general complexity, galaxy formation is best approached theoretically using numerical simulations. As a minimum, realistic simulations must be able to follow the evolution of the dark matter from appropriate initial conditions, together with the coupled evolution of a dissipative gas component that will eventually produce the visible regions of galaxies. The advent of large computers, particularly parallel supercomputers, together with the development of efficient algorithms, has given considerable impetus to this approach in recent years.

Both Eulerian and Lagrangian gas dynamics techniques have been used to simulate galaxy formation. Currently, only the latter have achieved sufficient dynamic range to resolve galaxy formation within cosmological volumes. For example, the recent large Eulerian simulations of Blanton et al. (1999) have spatial resolution of order $`300500\mathrm{kpc}`$, whereas typical Lagrangian techniques can readily achieve spatial resolutions of order $`20\mathrm{kpc}`$ or better (e.g. Carlberg, Couchman & Thomas 1990, Katz, Hernquist & Weinberg 1992, 1999, Evrard, Summers & Davis 1994, Frenk et al. 1996, Navarro & Steinmetz 1997). Recently, the power of the traditional, fixed-grid Eulerian approach has been greatly enhanced by the inclusion of adaptive mesh refinement techniques which have been applied to the study of the early phases of galaxy formation (Abel et al. 1998).

In this paper, we carry out simulations of galaxy formation using the Lagrangian technique of “Smooth particle hydrodynamics” or SPH (Gingold & Monaghan 1977; Lucy 1977). We focus on the thermodynamic history of the gas that cools to make galaxies, on the rate at which galaxies build up with time, on their abundance at the present day and on their spatial and velocity distributions. Early results from one of the simulations that we analyse here were presented in a previous paper (Pearce et al. 1999), where we focused on the spatial distribution of galaxies. Our work follows on from earlier attempts to grow virtual galaxies in a cosmological volume using the SPH technique, most notably by Carlberg, Couchman & Thomas (1990), Katz et al. (1992, 1996, 1999), Evrard et al. (1994) and Frenk et al. (1996). These simulations succeeded in resolving individual galaxies and allowed the first investigations of their distribution and of the dynamics of galaxies in clusters. However, the volumes modelled in these early simulations were rather small, permitting only limited statistical studies. Our simulations contain about one order of magnitude more particles and produce about 40 times more galaxies than the largest previous simulation of this kind by Katz, Weinberg & Hernquist (1996). SPH simulations with even higher resolution have been used to study the formation of individual galaxies (e.g. Navarro & Steinmetz 1997, Weil, Eke & Efstathiou 1998).

Simulations of the kind studied here are complementary to semi-analytic studies of galaxy formation (e.g. Kauffmann et al. 1997, 1999a,b; Diaferio et al. 1999; Guiderdoni et al. 1998; Benson et al. 1999a,b; Somerville & Primack 1999, Cole et al. 2000.) The main advantage of the hydrodynamic simulations is that they follow the dynamics of diffuse cooling gas in full generality whereas the semi-analytic treatment involves a variety of approximations such as spherical symmetry and quasistatic evolution. A comparison of results from the two methods has recently been performed by Benson et al. (2000c). Both approaches require a phenomenological treatment of critical processes such as star formation and feedback since these involve scales far below the resolution limit of all current simulations.

The remainder of this paper is laid out as follows. In Section 2 we describe our simulations and the physical approximations we have made. In section 3 we discuss the evolution of the global properties of the gas, the extraction of a dark matter halo catalogue, and the properties of the gas within these halos. In Section 4 we focus on the properties of galaxies, their formation histories, luminosity function, age distribution, etc. In Section 5 we consider the correlation function of the galaxies as a function of redshift and galaxy mass, the relative pairwise velocity distributions of galaxies and dark matter, and we determine how well the galaxies trace the mass distribution within the largest dark matter halos. We conclude in Section 6. This work is part of the programme of the Virgo consortium for cosmological simulations.

## 2 Simulations

Our simulations were performed using a parallel, adaptive, particle-particle, particle-mesh SPH code (Pearce & Couchman 1997) which can be run in parallel on the Crays T3D and T3E, or serially on a single processor workstation. This code is essentially identical in operation to the publicly released version of hydra, described in detail by Couchman, Thomas & Pearce (1995).

We have simulated galaxy formation in two cold dark matter (CDM) models (e.g. Davis et al. 1985) with the same values of the cosmological parameters as the SCDM and $`\mathrm{\Lambda }`$CDM models of Jenkins et al. (1998). For SCDM these are as follows: mean mass density parameter, $`\mathrm{\Omega }_0=1`$; cosmological constant, $`\mathrm{\Lambda }/(3H_0^2)=0`$; Hubble constant (in units of 100 km s-1 Mpc-1), $`h=0.5`$; and rms linear fluctuation amplitude in 8 $`h^1\mathrm{Mpc}`$ spheres, $`\sigma _8=0.6`$. For $`\mathrm{\Lambda }`$CDM , the corresponding parameter values are: $`\mathrm{\Omega }_0=0.3`$; $`\mathrm{\Lambda }/(3H_0^2)=0.7`$; $`h=0.7`$; and $`\sigma _8=0.9`$. Each simulation followed 2097152 dark matter particles and 2097152 gas particles in a cube of side $`100\mathrm{Mpc}`$ and required around 10000 timesteps (and $`10^5`$ processor hours on a Cray-T3D) to evolve from $`z=50`$ to $`z=0`$. In both cases, we employed a $`\beta `$-spline gravitational softening which remained fixed at the Plummer-equivalent comoving value of $`50\mathrm{kpc}`$ until $`z=1.5`$ ($`z=2.5`$ for $`\mathrm{\Lambda }`$CDM ). Thereafter, the softening remained fixed at $`50/(1+z)=10h^1\mathrm{kpc}`$ in physical coordinates, and the SPH smoothing length was set to match this value.

We deliberately chose to set the gas mass per particle to be $`2\times 10^9\mathrm{M}_{}`$ in both cosmologies so as to have identical mass resolution. As we typically smooth over 32 SPH neighbours, the smallest resolved objects have a gas mass of $`6.4\times 10^{10}\mathrm{M}_{}`$. The baryon fraction, $`\mathrm{\Omega }_b`$, was set from the nucleosynthesis constraint, $`\mathrm{\Omega }_bh^2=0.015`$ (Copi, Schramm & Turner 1995) and although this is somewhat lower than the more recent limits given by Tytler et al. (2000), it was the current value at the time these simulations were carried out. We also assumed a constant gas metallicity of 0.3 times the solar value in both cases. We expect more gas to cool globally in the $`\mathrm{\Lambda }`$CDM model because structure forms earlier in this model. With our chosen parameters, our simulations were able to follow the cooling of gas into galactic-sized dark matter halos. The resulting “galaxies” are typically made up of $`501000`$ particles and so have cold gas masses in the range $`10^{11}\mathrm{M}_{}2\times 10^{13}\mathrm{M}_{}`$. With a spatial resolution of $`10h^1\mathrm{kpc}`$, we cannot resolve the internal structure of these galaxies and we must be cautious about the possibility of enhanced tidal disruption, drag, and merging within the largest clusters. However, as we argue below, there is no evidence that this is a major problem.

As in all numerical simulation work, approximations and compromises are necessary if the calculation is to completed within a reasonable amount of time. An important choice is the effective resolution of the simulation. Increasing the spatial resolution requires increasing both the number of particles (to prevent 2-body effects) and the number of timesteps (to follow smaller structures). Our chosen value of $`10h^1\mathrm{kpc}`$ for the gravitational softening is larger than the scalelength of typical galaxies, but is the best that we can do with our current computer resources. For a given simulation volume, the number of particles determines the mass resolution. In cosmological simulations that involve only collisionless dark matter, the choice of particle number is often driven by computer memory limitations. Cosmological gas dynamics simulations are more computationally demanding than collisionless simulations and memory considerations do not, at present, limit their size. Since approximately 10% of all the baryons in the universe have cooled by the present to form stars and dense galactic gas clouds (Fukugita, Hogan & Peebles 1998), and around 50 particles are required to resolve a galaxy in a simulation, a simple estimate of the number of gas particles required to obtain $`N_{gal}`$ galaxies is: $`10\times N_{gal}\times 50`$, or roughly half a million particles for $`1000`$ galaxies. In practice, many galaxies contain more than the minimum number of particles and so this is an underestimate. Since galaxies of the characteristic luminosity, $`L_{}`$ or greater, have a space density of $`0.01/(h^1\mathrm{Mpc})^3`$, a simulation of a volume like those studied by Jenkins et al. (1998), of side $`250h^1\mathrm{Mpc}`$, would produce over $`100,000L_{}`$ galaxies and require around $`10^8`$ gas particles. For this work, we have chosen to use just over $`2\times 10^6`$ million gas particles in a volume of side only 100 Mpc.

In addition to the choice of simulation parameters, other choices need to be made to model physical processes that operate on sub-resolution scales. Star formation and the associated feedback from energy released in supernovae and stellar winds are the most obvious sub-resolution processes relevant to simulations of galaxy formation. As the resolution is increased, progressively smaller, denser halos are resolved whose gas has a short cooling time. In hierarchical clustering theories, a simple model predicts that most of the baryonic material should have cooled (and presumably turned into stars), in small, subgalactic objects at high redshift (White & Rees 1978, Cole 1991, White & Frenk 1991). Since this does not appear to be the case in the universe, some mechanism must have prevented much of the gas from cooling. Reheating of the gas by the energy released in the course of stellar evolution is a commonly invoked mechanism to quench star formation in small galaxies. This, however, is still a very poorly understood process. Various simplified models have been implemented in simulations of galaxy formation (e.g. Katz 1992; Navarro & White 1993; Metzler & Evrard 1994; Steinmetz & Muller 1995; Evrard, Metzler & Navarro 1996; Katz et al. 1996; Mihos & Hernquist 1994, 1996; Gerritsen & Icke 1997; Navarro & Steinmetz 1997).

In this paper, we have taken the simpler route of neglecting feedback effects altogether. Instead, we appeal to the resolution threshold of the simulation to prevent all the gas from cooling. This is analogous to assuming a model in which gas cooling is suppressed within dark matter halos of mass below $`6.4\times 10^{10}\mathrm{M}_{}`$), whilst leaving larger halos unaffected. Although relatively simple, this model gives results that are not too dissimilar from those produced by detailed semi-analytic models in which feedback is treated with greater care. In semi-analytic modelling, the star formation and feedback prescriptions are tuned by the requirement that the model should reproduce observable properties of the local galaxy population such as the luminosity function. The difference in the amounts of gas that cool in our simulations and in the semi-analytic model of Cole et al. (2000) are illustrated in Fig. 1 of Pearce et al. (1999). Near the resolution threshold, approximately 50% more gas cools in the simulation than in the semi-analytic model and the difference falls rapidly with increasing galaxy mass. A detailed comparison of results from these simulations and semi-analytic modelling may be found in Benson et al. (2000c).

The second area in which sub-resolution physics is important is the interaction of cool gas clumps with the surrounding hot gas. In reality, the gas within such clumps would presumably become incorporated into molecular clouds or turned into stars. These processes occur on scales that are several orders of magnitude below the resolution limit of current simulations and so must be modelled phenomenologically. Different authors have chosen to do this in different ways. For example, Evrard et al. (1994) and Frenk et al. (1996) opted to leave the cold gas alone. Although this is the simplest possible approach, it has the disadvantage that the smoothing inherent in the SPH method can artificially boost the cooling rate of hot gas that comes into contact with the gas that has cooled to very high density inside a galaxy. This process is illustrated in the lower panels of Fig. 1 which show the temperature and density profile of a galaxy in the simulation. Model galaxies are typically surrounded by hot gaseous coronae. At the interface with the cold galactic gas, the density of the hot coronal gas is overestimated and its cooling rate is artificially enhanced. This gives rise to the rapidly cooling material which is clearly visible in the region between the effective search lengths of the cold and hot material, indicated by the dashed and dotted lines in the figure respectively. The artificial cooling at this contact interface arises because SPH is not designed to deal with the steep density gradients that develop in a multiphase medium. As Thacker et al. (2000) have shown, this limitation is present to some extent in all standard implementations of SPH, although recently Ritchie & Thomas (2000) have proposed an alternative formulation of SPH utilising the gas pressure rather than the gas density which overcomes this limitation.

An alternative approach to leaving the galactic cold dense gas alone is to replace it with collisionless “star” particles, thus decoupling it entirely from the rest of the baryonic material (Navarro & White 1993, Steinmetz & Muller 1995, Katz et al. 1996). While this might appear as a more realistic solution than leaving cold material as gas, in practice, the relatively large softening lengths normally used in SPH simulations imply unrealistically low binding energies for the star particles. As a result, stellar objects are prone to disruption by collisions and tidal encounters, particularly in strongly clustered regions.

In this study, we have chosen an approach which is intermediate between the two extremes just discussed. In our phenomenological model, hot gas, defined to be hotter than $`10^5\mathrm{K}`$, is prevented from interacting with (galactic) material at a temperature below $`12000\mathrm{K}`$ (but not vice-versa). All other SPH forces remain unchanged. The threshold temperature of $`10^5\mathrm{K}`$ is less than the virial temperature of the smallest resolved halos in our simulations. In our implementation of SPH, this assumption preserves force symmetry and momentum conservation, and effectively prevents overcooling of hot gas induced purely by the presence of a nearby clump of cold, dense gas. This is clearly seen in the top panels of Fig. 1. In contrast to the standard model in which the hot and cold phases remain strongly coupled, there is now no interface discontinuity. As a result, the density profile in the vicinity of a cold clump remains smooth and a “cooling flow” develops around it. Recently this effect has been confirmed by Croft et al. (2000) who suggest several alternative routes to circumvent it.

In our scheme, the cold gas “feels” all the usual SPH forces, including viscous drag when it moves in a hot intergalactic medium, and interacts with other cold gas clouds in the usual way. The outer accretion shock of a galaxy cluster is still properly modelled as the infalling material can “see” the hot gas in the cluster. Our approximation has only a small effect on the masses of small galactic objects as may be seen in Fig. 2, where we compare the masses of galaxies formed in two small test simulations. These had similar parameters and mass resolution to those employed in our main simulations, but in one case the hot gas was decoupled from the cold gas and in the other case it was not. For most galactic-sized objects, this change produced a slight reduction in mass, but for the largest object, which has over 1000 particles without decoupling, the mass is considerably reduced to 370 particles. This effect, as well as the results of varying several simulation parameters are explored further by Kay et al. (2000).

## 3 The properties of the gas

In this section, we examine both the evolution of the gas in the simulation as a whole and the properties of gas within virialised dark matter halos.

### 3.1 Global properties of the gas

Following Kay et al. (2000), we define three different gas phases: (i) cold galactic gas (temperature below $`12000\mathrm{K}`$ and overdensity greater than 10); (ii) hot halo gas (temperature above $`10^5\mathrm{K}`$); and (iii) cold, uncollapsed gas (overdensity less than 10; this is essentially diffuse gas in “voids.”) We ignore the effects of photoionization because at this mass resolution a photoionising field has a negligible effect upon galactic or halo properties. In practice, little gas is both cold and at overdensities between 5 and 1000, so that effectively the uncollapsed and galaxy phases are completely disjoint. The evolution of the fraction of gas in each of these phases is illustrated in Fig. 3 for the two cosmological models that we have calculated, SCDM (bold lines) and $`\mathrm{\Lambda }`$CDM (thin lines). These gas fractions are very similar to those derived using a semi-analytic model of cooling gas implemented in these simulations by Benson et al. (2000c).

The properties of cold, uncollapsed material in voids are only followed in an approximate manner in our simulations. The SPH method generally has poor spatial resolution in regions in which the gas is diffuse because it must smooth over a sparse particle distribution. This limitation is exacerbated in our specific implementation which imposes a maximum search length, resulting in a minimum resolvable density of a few times the mean density in the computational volume. Test simulations without such a minimum resolvable density demonstrate that these inaccuracies do not affect the properties of the galaxies because the gas forces are small in void regions and have no effect on the dense regions that ultimately form galaxies (Kay et al. 2000). However, outwardly propagating shocks generated by the formation of large clusters are poorly followed in underdense regions and so we can say little about the detailed properties of voids (beyond their mere presence). This is a regime in which Eulerian techniques are more accurate than SPH (e.g. Cen & Ostriker 1996).

As the evolution proceeds the cold, diffuse phase is depleted as gas falls into dark halos where it is shocked into the hot phase, thereafter cooling into galaxies. The time derivative of the dashed curve that tracks the galaxy phase in Fig. 3 effectively defines the global star formation rate in the computational volume. This is plotted in Fig. 4. The star formation rate rises rapidly at early times, has a broad maximum around $`z12`$ and then gently declines to the present. This behaviour is broadly similar to that seen in real data (Steidel et al. 1999), indicating that the limited resolution of the simulations and other approximations do not result in too unrealistic a model. At late times, Fig. 3 shows that there is a substantial transfer of gas from the uncollapsed phase into the hot phase, with around $`50\%`$ of the gas ending up at temperatures above $`10^5\mathrm{K}`$ in both cosmological models. At $`z1`$, when only around $`30\%`$ of the gas is above $`10^5\mathrm{K}`$, the gas fractions in the different phases in our simulations agree well with those in the simulations of Evrard et al. (1994) (which were stopped at this point.)

As expected, more gas cools into the galaxy phase in the $`\mathrm{\Lambda }`$CDM than in the SCDM model. This is due to our adoption of a single gas particle mass in the two simulations which results in the same effective resolution in both cases. More gas cools at early times in the $`\mathrm{\Lambda }`$CDM simulation because structure forms earlier in this case. If all the gas that cools into galaxies were assumed to turn into stars, then the SCDM model would have roughly the observed mean mass density in stars. In the $`\mathrm{\Lambda }`$CDM model, on the other hand, about twice as much gas cools and, if all this gas were assumed to turn into stars, an uncomfortably large value of the mass-to-light ratio would then be required to match the observed mean stellar mass density. The amount of gas that cools in our simulations is determined by the mass resolution and a small decrease in the resolution of the $`\mathrm{\Lambda }`$CDM model would reduce the final amount of cold gas to the same level as in the SCDM model.

### 3.2 Properties of gas in halos

To construct a catalogue of dark matter halos from the simulations, we first identified cluster centres using the friends-of-friends grouping algorithm (Davis et al. 1985) with a small value of the linking length, $`b=0.05`$. We then grew a sphere around the centre of mass of each of these halos until the mean overdensity within it reached a value of 178 in the case of SCDM or 324 in the case of $`\mathrm{\Lambda }`$CDM . These are the overdensities corresponding to the virial radius according to the spherical top-hat model for clusters (Eke, Cole & Frenk 1996.) The resulting set of halos was cleaned by removing the smaller of any two overlapping halos. This procedure produced catalogues of 1320 halos in SCDM and 1795 halos in $`\mathrm{\Lambda }`$CDM , with more than 28 particles per halo in each case.

To check whether the presence of baryons affects the halo catalogue, we ran a pure N-body simulation of the $`\mathrm{\Lambda }`$CDM model and extracted a halo catalogue from both it and the baryonic simulation with gas cooling. The mass function of halos in the test simulation is almost indistinguishable from that in the original simulation, as illustrated in Fig. 5. This figure also compares the mass functions in the simulations with the mass function given in Jenkins et al. (2000). This gives a better fit than the predictions of the Press-Schechter model (Press & Schechter 1974). Although for the values of $`\mathrm{\Omega }_b`$ that we have assumed here, the presence of gas has little effect on the halo mass function, it does affect the central structure of the halo into which it cools (Pearce et al. 2000).

In Fig. 6 we plot the mass fraction of gas within the virial radius of the halos normalized to the mass fraction in the simulation as a whole, together with the proportions in the hot and cold phases. The total gas mass fraction varies little with halo mass over the range $`10^{12}\mathrm{M}_{}<M<10^{15}\mathrm{M}_{}`$, although there is some indication of an upturn at low masses in the $`\mathrm{\Lambda }`$CDM case and at high masses in the SCDM case. The mean value of the gas fraction is $`0.85`$ for $`\mathrm{\Lambda }`$CDM and $`0.8`$ for SCDM with a scatter of $`0.05`$, indicating that the gas is slightly more extended than the dark matter in both models. A few small mass halos have a gas mass fraction greater than unity. These halos, however, are near the resolution limit of the simulations. The crosses show the mass fraction of hot ($`T>12000\mathrm{K}`$) gas in the halos. This increases rapidly with halo mass, from $`1020\%`$ for the smallest halos in the simulations to $`6080\%`$ for the largest. The fraction of cold ($`T<12000\mathrm{K}`$) gas shows the opposite trend and varies from $`7080\%`$ at the small mass end to $`1020\%`$ at the high mass end.

The fraction of baryonic material that cools within dark matter halos can also be calculated using semi-analytical techniques (White & Frenk 1991). As we discussed in the preceeding section, in semi-analytic models the cooling of gas is regulated by feedback processes due to energy liberated by stars. In these simulations we ignore feedback altogether and cooling in small mass objects is then limited by resolution effects. This leads to some differences in the gas fractions derived with these two techniques. For example, in the semi-analytic model of Cole et al. (2000), the gas fraction of cold gas has a maximum of $`40`$% at a halo mass of $`10^{12}\mathrm{M}_{}`$. This is close to the resolution limit of our simulations. At small masses, the cold fraction is suppressed by feedback and at larger masses, it is suppressed by the long cooling time of the diffuse gas in halos. The first of these effects is not present in our simulations (which do not resolve the regime where the feedback is important in the semi-analytic models), but the second one is. On these large mass scales, the simulations agree well with the semi-analytic results and even over the entire range of mass, the difference is only of the order of 50% (Benson et al. 2000c). Overall, somewhat more mass cools in the simulations than in the semi-analytic model, suggesting that feedback effects are not negligible even above the resolution limit of our simulations.

The hot gas that falls into a halo is quickly thermalized by shocks and heated to the virial temperature. This is illustrated in Fig. 7 which shows the relation between the mean, mass-weighted gas temperature (for gas with $`T>10^5\mathrm{K}`$) and the virial mass of the host halo. For masses above a few times $`10^{12}\mathrm{M}_{}`$, there is a very tight correlation between gas temperature and halo mass, indicating that the gas is relaxed and close to equilibrium. At the small mass end, there are some halos with anomalously large temperatures. These tend to be in the vicinity of larger halos whose own hot halos contaminate the smaller ones, boosting their temperature. The dashed line is the relation $`MT^{\frac{3}{2}}`$ predicted by a simple, spherically symmetric, equilibrium model for the halo and its gas. The simulations follow this relation well over most of the mass range plotted in the figure. The X-ray properties of the largest halos in our simulations have been studied by Pearce et al. (2000).

## 4 The properties of galaxies

We now consider the properties of the population of galaxies in the simulations, focussing on their abundance, ages, merger histories, and individual star formation rates.

### 4.1 Demographics of the galaxy population

We identify “galaxies” in the simulations with dense clumps of cold gas. Locating them is straightforward except in the small number of cases when they are undergoing a merger or are being tidally disrupted within a large cluster. To find galaxies, we use the same friends-of-friends algorithm that we employed for finding dark halos, but with a linking length of only $`b=0.02`$ for SCDM and $`b=0.0134`$ for $`\mathrm{\Lambda }`$CDM, values that are 10% of those required to obtain virialised halos in these cosmologies. When calculating the mass of each galaxy, we consider only particles with temperature below $`12000\mathrm{K}`$, a condition that rejects only very few particles. In practice, as the cold clumps of gas that make up the galaxies have contracted by a large factor the galaxy catalogues and the masses and positions of each are insensitive to these choices. We find a total of about 2000 resolved galaxies in each of the two simulations at the present day.

The evolution in the abundance of galaxies is shown in Fig 8. Here, we characterize each galaxy by its cold gas mass and plot the cumulative abundance at various epochs in our two simulations. In the $`\mathrm{\Lambda }`$CDM model, the first substantial galaxies begin to form at a redshift $`z5`$, although some fairly massive objects are present even before this. The massive galaxy population builds up rapidly between $`z5`$ and $`z3`$ when it increases approximately ten-fold. A similar increase occurs between $`z3`$ and $`z1`$. As time progresses, the mass of the largest galaxies also increases. After $`z=0.5`$, the growth in the population slows down somewhat because the rate at which new galaxies are forming is more or less balanced by the rate at which existing galaxies are destroyed by mergers and tidal disruption. At $`z=0`$, the abundance of galaxies with gas mass $`M>10^{11}\mathrm{M}_{}`$ approximately matches the observed abundance of $`L_{}`$ or brighter galaxies of $`0.003`$ per $`h^3\mathrm{Mpc}^3`$. The evolution in the SCDM model is qualitatively similar, but shifted towards lower redshifts.

The creation sites of galaxies at recent times are illustrated in Fig. 9 which shows the locations where new gas cools between $`z=0.3`$ and the present in the SCDM model. (Broadly similar behaviour is seen in the $`\mathrm{\Lambda }`$CDM model.) The radius of the circle around each of these locations is scaled by the mass of new cold material. Some gas cools onto preexisting objects (note the small circles at the centre of some clusters), but 383 new resolved galaxies, about 15% of the total, formed during this period. The new galaxies are small and tend to form along filaments or in low density regions, away from the larger clusters.

Galaxies are destroyed preferentially in large clusters, where tidal forces are largest. For example, about $`20\%`$ of the cold material within the virial radius of the largest halo at $`z=0.3`$ is reheated before $`z=0`$, and is associated with the destruction of about half the small (around 1/3$`L_{}`$) galaxies in this halo. As discussed earlier, tidal disruption of galaxies in clusters is likely to be enhanced by the limited spatial resolution of the simulations which results in artificially low binding energies for small objects. Larger simulations are required to investigate this important problem in detail.

The major consequence of “decoupling” the hot gas from the cold phase, as discussed in the preceeding section, is that no supermassive objects form. The biggest galaxy in the $`\mathrm{\Lambda }`$CDM model contains $`5\times 10^{12}\mathrm{M}_{}`$ of cold gas. The turnover in the abundance of galaxies at small masses in the two simulations is a numerical artifact. At least 32 particles (corresponding to a mass of $`6.4\times 10^{10}\mathrm{M}_{}`$) of gas must be present for cooling to be efficient; objects below this are in the process of forming or are being ablated as they move in a hot halo. The effective mass threshold in both simulations is actually slightly higher than the nominal value and corresponds to a mass threshold around $`10^{11}\mathrm{M}_{}`$, or 50 gas particles.

The number of galaxies within the virial radius of each halo at the present day is plotted against the halo mass in Fig. 10. The largest cluster contains nearly 30 large galaxies and has a virial radius of over $`1.5h^1\mathrm{Mpc}`$. Its galaxy content is consistent with those of $`10^{15}h^1\mathrm{M}_{}`$ clusters which typically contain 30 - 100 $`L_{}`$ galaxies. This is, in fact, the only large Abell-type cluster that formed in our simulations. All dark halos in the simulations that contain more than a few hundred particles have at least one galaxy within them. There are no empty halos of mass greater than $`10^{13}\mathrm{M}_{}`$ in either cosmology.

### 4.2 Star formation and the age of galaxies

Although, for simplicity, we have not included any prescription for star formation in our simulations, we can still derive some general conclusions regarding the properties of the stars expected to form. Stars form from gas that has cooled into galactic dark halos. The approximate time at which an individual gas particle is available for star formation is easy to estimate in the simulation. We take this to be halfway between the output time at which the particle first has an overdensity greater that 10 and a temperature smaller than $`12000\mathrm{K}`$, and the previous output time. These criteria have previously been shown to pick out solely baryonic material that has cooled into galactic objects (e.g. Kay et al. 2000). We assume that as soon as a particle becomes available for star formation, a mass of stars equal to the mass of the particle forms.

By calculating a star formation time for each particle that ends up in galaxy, we can derive a mean, mass-weighted stellar age or formation time, $`t_{\mathrm{form}}`$, for each galaxy as

$$t_{\mathrm{form}}=\underset{i=1}{\overset{N}{}}\frac{t_i}{N},$$
(1)
where the sum is over all particles found within the final galaxy. This definition of formation time is based exclusively on the age of the stellar population of a galaxy and takes no account of when the galaxy was assembled. An alternative definition of formation epoch is the time at which half the final galaxy mass was assembled into a single progenitor. As illustrated in Fig. 11, these two times can be very different. For the largest galaxy in the SCDM simulation volume (which has a gas mass of $`2.4\times 10^{12}\mathrm{M}_{}`$ and for which we have enough resolution to follow many progenitors), the average formation time of its stars corresponds to redshift $`z2.5`$, but half of the stellar material is only assembled into a single distinct object at $`z0.7`$. The formation history of this galaxy bears little resemblance to the classic spherical top-hot collapse model. Rather, the galaxy is assembled through repeated mergers of sub-units, many of which have a mass close to the resolution threshold of the simulation and which are originally spread over quite a large comoving volume. These fragments tend to line up with filaments of the dark matter distribution, giving rise to an anisotropic fragment distribution which is particularly visible at $`z=3`$ in Fig. 11. The mass distribution of the progenitors of this massive galaxy at various epochs is shown in Fig. 12.

Fig. 13 shows the distribution of the mean stellar ages (as a fraction of the lookback time) for galaxies in the SCDM simulation. Contrary to the naive expectation for a hierarchical scenario, there is a weak trend for the largest galaxies to be the oldest. This may be an artificial consequence of the resolution limit of our model which is at a fixed mass. Models with more realistic feedback and no lower mass limit produce a weak trend in the opposite sense for galaxies in the field and no trend at all for galaxies in clusters (Kauffmann 1996).

An estimate of the star formation rate (SFR) in our simulated galaxies can be derived from the gas mass that cools between successive output times. This crude measure (we only have 10 output times) is shown for the 20 largest galaxies selected at the end of the SCDM simulation in Fig. 14. Clearly, the inferred aggregate SFR within each of these objects is very high but, as Fig. 11 shows, these galaxies are broken up into many small precursor objects at early times each of which has a modest SFR. The aggregate SFR for the most massive objects today is sharply peaked at $`z3`$. This contrasts with the star formation history of more typical galaxies, illustrated in Fig. 15. This shows SFRs for every $`50^{\mathrm{th}}`$ galaxy ordered by mass. These curves are more extended in time and the mean SFR for the average galaxy is much lower than the aggregate star formation rate of the largest objects. The SFRs generally decline towards low redshift as galaxies consume the gas that is able to cool by the present day. The SFR integrated over all the galaxies in the simulation is shown in Fig. 4.

By convolving the SFR in each galaxy with the Bruzual & Charlot (1993) stellar population synthesis model, we have calculated the $`K`$band luminosity of each galaxy. We have assumed a Salpeter Initial Mass Function with a low mass cutoff of $`0.1\mathrm{M}_{}`$ and solar metallicity. (This is higher than our assumed metallicity of 0.3 solar for the ambient gas, since the galactic material will have been enriched by reprocessed metals.) As shown in Fig. 16, the model luminosity function in both models has a similar shape to the $`K`$band luminosity function measured by Gardner et al. (1997). Furthermore, the model luminosity function can be made to agree quite well with the data by assuming that only a fraction $`1/\mathrm{{\rm Y}}`$ of the stellar mass is luminous (and the rest is in brown dwarfs), with $`\mathrm{{\rm Y}}=2.8`$ for $`\mathrm{\Lambda }`$CDM  and $`\mathrm{{\rm Y}}=0.83`$ for SCDM. The former is somewhat larger and the latter somewhat smaller than the values inferred from observations (Cole et al. 2000 and references therein) because, as discussed earlier, we have used the same mass per gas particle in the two simulations and, as a result, more gas cools globally in the $`\mathrm{\Lambda }`$CDM model than in the SCDM case. Due to limited resolution, only the bright end of the luminosity function ($`L>L_{}`$) is accessible in our simulations.

## 5 Clustering properties of the Galaxies

In this section we examine both the two-point spatial clustering and the pairwise velocity dispersions of the galaxies as a function of separation and we compare them to observational determinations. We also investigate the dependence of spatial clustering on galaxy mass. We test a virial mass estimator, that of Heisler, Tremaine, & Bahcall (1985, hereafter HTB), on the galaxy samples contained within a few of the largest dark matter halos to see how well the true mass within the virial radius is recovered.

Before proceeding with our analysis, we need to define a suitable galaxy sample. As discussed earlier, the galaxies are easily identified because of their very high gas density compared to the mean. The exact membership of the galaxy catalogues obtained using a group finder is relatively insensitive to the input parameters used to select them. In this section, galaxy catalogues for both cosmologies were generated by applying a “friends-of-friends” group finder with a linking length of $`0.02(1+z)`$ on a subset of gas particles with temperature less than $`12000\mathrm{K}`$. All gas particles were counted for the purposes of defining the linking length itself. Catalogues were constructed at $`z`$=0, 1 and 3. They contain 2263, 2360 and 1502 galaxies respectively for $`\mathrm{\Lambda }`$CDM and 2089, 1647 and 550 galaxies for SCDM, counting objects of 8 particles or more.

### 5.1 The galaxy autocorrelation function

Figure 17 shows the redshift evolution of the two-point correlation function of galaxies in the $`\mathrm{\Lambda }`$CDM and SCDM models. In both cases, the galaxy correlation function evolves little over the range $`z=`$0–3, in contrast to the mass autocorrelation functions also shown in the figure. The galaxy correlation functions are essentially unbiased on the largest scales relative to the $`z=0`$ mass correlation function but show a distinct anti-bias on scales below 3 $`h^1\mathrm{Mpc}`$ for $`\mathrm{\Lambda }`$CDM and 2 $`h^1\mathrm{Mpc}`$ for SCDM, and a positive bias on scales below $`100`$$`h^1\mathrm{kpc}`$ in both cases. Also plotted is Baugh’s (1996) measured galaxy autocorrelation function determined by deprojecting the angular correlation function of galaxies in the APM catalogue. (The curve plotted assumes that clustering is fixed in comoving coordinates.) On small scales, the galaxy correlation function in the simulations is encouragingly close to the real galaxy correlation function. On large scales, the SCDM model fails because it lacks sufficient large-scale power, but the $`\mathrm{\Lambda }`$CDM model continues to do well. The departure of the $`\mathrm{\Lambda }`$CDM correlation function of both galaxies and dark matter below the APM result at the largest scales plotted is due to the small volume of the simulation. A better comparison between the APM data and the dark matter correlation function at large separations is shown in Fig. 5 of Jenkins et al. (1998).

To examine the dependence of spatial clustering on galaxy mass, we subdivided the $`\mathrm{\Lambda }`$CDM $`z=0`$ sample into low-, middle- and high-mass subsamples, each with approximately 750 galaxies. In terms of the gas particle number, $`N_{part}`$, the three subsamples have $`8N_{part}<50`$, $`50N_{part}<95`$ and $`N_{part}95`$ respectively, corresponding to $`K`$-band luminosities in the ranges ($`20`$, $`22.3`$), ($`22.3`$, $`22.8`$) and brighter than $`22.8`$ . Figure 18 shows the 2-point correlation function for the three galaxy subsamples. The clustering strength shows a clear trend, with the most massive galaxies having the strongest clustering. The difference between the medium- and low-mass samples is much less pronounced than the difference between the medium- and high-mass samples.

To conclude, we find that the clustering properties of the galaxies in the simulations are markedly different from those of the dark matter. In particular, the galaxy correlation function evolves weakly with redshift and, coincidentally, appears unbiased on large-scales at the present epoch. At the present day, the correlation amplitude of the galaxies in both models is weaker than that of the dark matter at intermediate scales around 1 $`h^1\mathrm{Mpc}`$ and this is particularly pronounced for $`\mathrm{\Lambda }`$CDM . The galaxy correlation function in the $`\mathrm{\Lambda }`$CDM model is close to that determined observationally. Its shape is rather noisy because of the small sample, but it is much better fit by a power-law than the mass distribution. The clustering strength of galaxies increases with galaxy luminosity, but this effect is only strong at very bright luminosities.

### 5.2 The galaxy projected pairwise velocity dispersion.

In this subsection, we compare the projected pairwise velocity dispersions of galaxies and dark matter. As is well-known, a reliable estimate of this statistic requires a rather large sample volume (Marzke et al. 1995, Mo, Jing & Boerner 1997). The sample volumes that we analyse in this paper are probably too small to provide an accurate estimate of the global pairwise velocity dispersions, and we will therefore concentrate on the relative dispersions of galaxies and dark matter. Nonetheless, with this caveat in mind, we plot, for interest, the pairwise velocity dispersions determined from the LCRS catalogue by Jing, Mo & Boerner (1998).

Fig. 19 shows the measured projected pairwise velocity dispersions of both dark matter and galaxies in the $`\mathrm{\Lambda }`$CDM and SCDM models. The detailed definition of the projected pairwise dispersion is given in Jenkins et al. (1998). Because of our smaller simulation volumes, we have decreased the length-scale of the projection from $`\pm `$25 $`h^1\mathrm{Mpc}`$ to $`\pm `$10 $`h^1\mathrm{Mpc}`$ and $`\pm `$8 $`h^1\mathrm{Mpc}`$ for the $`\mathrm{\Lambda }`$CDM and SCDM models respectively. We plot the same LCRS data points (taken from Jing et al. 1998), as plotted in Fig. 11 of Jenkins et al. (1998).

There is a remarkable difference between the pairwise velocity correlations of galaxies and those of the dark matter on scales below $`2`$$`h^1\mathrm{Mpc}`$ in both cosmological models. For $`\mathrm{\Lambda }`$CDM , the galaxy pairwise velocities turn out to be not much higher than the observational data points. In fact, the dark matter pairwise dispersion in our 70 $`h^1\mathrm{Mpc}`$ $`\mathrm{\Lambda }`$CDM simulation cube is around 100km s$`^1`$higher than the value obtained by averaging over larger 239.5 $`h^1\mathrm{Mpc}`$ boxes in Jenkins et al. (1998; Fig. 11). One might infer from this that the global average for the projected galaxy pairwise dispersion is also overestimated. The difference between the behaviour of the galaxies and dark matter becomes smaller with increasing pair separation and is similar on the largest scales probed by these simulations.

Why do the galaxies have a lower pairwise dispersion than the dark matter? To help answer this question we have constructed a “shadow galaxy” catalogue for $`\mathrm{\Lambda }`$CDM by selecting the nearest dark matter particle to each galaxy. The shadow catalogue has, by construction, virtually an identical 2-point correlation function to the galaxies. In Fig. 19 we also plot, as a dotted line, the projected pairwise dispersion for the shadow catalogue. This is almost identical to the projected pairwise dispersion of the galaxies. This shows that the strong difference between the galaxies and the dark matter arises mostly from the way in which galaxies populate dark matter halos rather than from a strong intrinsic bias in the velocities of galaxies relative to the dark matter particles at the same location. The fact that the dispersion of the shadow catalogue is everywhere higher than that of the galaxies indicates that some residual is present, but the closeness of the curves shows that this is not a large effect. Thus, the low galaxy velocity dispersion is a statistical effect which may, in principle, be driven by differences in the ratio of the number of galaxies to total mass in high- and low-mass halos and also by differences in the way in which galaxies are distributed within halos compared to the dark matter.

Our results for $`\mathrm{\Lambda }`$CDM are very similar to those in Benson et al. (2000a), but differ in some respects from those of Kauffmann et al. (1999b). In both papers, semi-analytical modelling of galaxy formation was used to create synthetic galaxy catalogues from the same $`\mathrm{\Lambda }`$CDM N-body simulation, but the detailed placement of galaxies in halos was different. Generally, in the semi-analytic models the efficiency of galaxy formation per unit mass of dark matter is a strong function of dark halo mass and peaks for halos with masses around $`10^{12}\mathrm{M}_{}`$. Benson et al. (2000a) found that by populating the dark matter halos unevenly as their semi-analytical model predicts, not only does one obtain a galaxy correlation function which matches the APM result of Baugh (1996) well, but also a comparable difference as found here in the projected pairwise dispersion of galaxies relative to the dark matter. Kauffmann et al. (1999b) find a very similar two-point correlation function, but rather little difference between the pairwise dispersions of galaxies and dark matter. As the analysis of Benson et al. (2000a) shows, this discrepancy can be traced back to relatively small variations in the precise form of the halo occupation number predicted in the two models, particularly for halos with mass $`>10^{13}\mathrm{M}_{}`$. These generate a large difference in the pairwise dispersions but have a much weaker effect on other statistics such as the two-point correlation function. A detailed comparison between our simulations and semi-analytical modelling techniques applied to dark matter only versions of our simulations is given in Benson et al. (2000c).

In summary, we find significant differences in the pairwise velocity dispersion of galaxies relative to the dark matter on scales below 2 $`h^1\mathrm{Mpc}`$ . The amplitude of this statistical effect appears sufficient to reconcile $`\mathrm{\Lambda }`$CDM with current observations. The differences in the velocity dispersions of galaxies and dark matter arise because the galaxies populate dark halos with a variable efficiency which decreases for the highest mass halos.

### 5.3 Virial theorem mass estimates of Galaxy clusters

In this section we employ a commonly used virial estimator to see how well the virial masses can be determined for a few of the most massive dark halos using the projected positions and velocities of the galaxies located within the halo. The estimator we have selected is one of those discussed by Heisler et al. (1985), specifically,

$$M_{VT}=\frac{3\pi N}{2G}\frac{\underset{i}{}V_{r,i}^2}{_{i<j}1/R_{ij}},$$
(2)
where $`V_{r,i}`$ is the line-of-sight velocity of the $`i`$th galaxy (in the observer frame where the cluster is at rest) and $`R_{ij}`$ is the projected distance between galaxies $`i`$ and $`j`$. HTB show that this estimator is fairly robust. It does assume, however, that the galaxies trace the dark matter. If the galaxies are, for instance, more concentrated within a dark halo than the dark matter then the kinetic energy-like upper sum in equation (2) will be underestimated, whilst the potential energy-like sum will be overestimated. The result then is a virial mass estimate which will be lower than the true mass. Additional sources of uncertainty in applying the estimator to galaxy clusters are that the dark halos are not that well isolated from their surroundings and may not be in close virial equilibrium particularly if viewed at an epoch of major mass accretion.

As a sample, we take the 9 most massive clusters from both $`\mathrm{\Lambda }`$CDM and SCDM models. The average number of galaxies in our $`\mathrm{\Lambda }`$CDM halos is N=12.4 and N=10.7 for SCDM. We project each cluster from a large number of different directions to get the full distribution of the estimated virial mass. We then average the ratio of the estimated mass to the true mass for the cluster samples. Using $`\mathrm{log}_{10}(M)`$ as the mass variable we find that the distribution of the estimated mass about the true mass is very similar in both cosmologies but it is slightly more peaked in $`\mathrm{\Lambda }`$CDM than SCDM. With only nine halos for each model, the conclusions that we can draw are not very strong. The distributions have significant tails towards low masses, but the exact form is dependent on only one or two of the halos. The estimator shows a bias in $`\mathrm{log}_{10}(M)`$ of $`0.08\pm 0.08`$ for $`\mathrm{\Lambda }`$CDM and $`0.02\pm 0.2`$ for SCDM which is consistent with zero. A more reliable statistic which does not depend so much on the tail of the distributions is the range $`\pm \mathrm{\Delta }\mathrm{log}_{10}M`$ about the true mass which includes 3/4 of the distribution. For $`\mathrm{\Lambda }`$CDM this is about 0.23 and for SCDM 0.2. In Figure 20 we combine the results for $`\mathrm{\Lambda }`$CDM and SCDM. For the combined sample, 3/4 of the distribution is included within a range of $`\pm 0.21`$. This spread of values is a little larger than the corresponding spread found by HTB, who quote a range of 0.2 for N=5, and 0.15 for N=10 but these are for realisations of virialised Mitchie models. This is not perhaps surprising given that the galaxies do not trace the mass exactly nor are the clusters in complete virial equilibrium.

In Figure 21 we compare the average radial distribution of the galaxies and mass within the virial radius for the 9 halos from each cosmology analysed above. The galaxy radial distribution is shown with a solid line and the mass radial distribution with a dotted line. There are 112 galaxies contributing to the $`\mathrm{\Lambda }`$CDM histogram and 97 to SCDM. All the histograms are normalised so as to have unit area underneath. Broadly speaking the galaxies and mass have a similar distribution. A $`\chi ^2`$ analysis of the galaxies and mass distributions yields $`\chi ^214`$ for 8 d.o.f. for both $`\mathrm{\Lambda }`$CDM and SCDM. Taken individually this is unremarkable with a higher value occurring approximately 8% of the time. The joint $`\chi ^2`$ has a probability which is only marginally more significant at 3%. Without a larger sample it is difficult to draw any firmer conclusion about the relative distributions except to note that several hundred galaxies are needed in order to distinguish the distributions.

In summary, we find that applying the estimator, eqn. 2, to the galaxies in the largest dark matter halos taken from either of our two simulations, one can infer the masses of the dark halos in a relatively unbiased way, albeit with a significant uncertainty for a single determination of a particular halo. Consistent with this conclusion, the radial distributions of galaxies and mass are similar. These mass estimates, of course, neglect projection effects in the identification of cluster galaxies which can add significantly to the uncertainty (eg van Haarlem, Frenk & White 1997). It is interesting to note that for clusters more massive that $`1.3\times 10^{14}`$ $`\mathrm{M}_{}`$, the mean baryon mass fraction in the form of galaxies is 17% and 8.4% for $`\mathrm{\Lambda }`$CDM and SCDM respectively. These values differ from the global fractions which are 12% for $`\mathrm{\Lambda }`$CDM and 7% for SCDM. Thus, the mass-to-light ratios of clusters cannot, in general, be used to estimate $`\mathrm{\Omega }_0`$ reliably without accurately accounting for the difference in the efficiency of galaxy formation in clusters compared to the universe as a whole.

## 6 Discussion and conclusions

Our simulations have many limitations and should be regarded primarily as illustrative of general trends likely to characterize galaxy formation in hierarchical CDM models. Most of these limitations stem from the limited resolution that is attainable with current computing resources if one wishes to simulate a representative volume of the universe. Our models do not include a self-consistent treatment of star formation and feedback. Instead, the ability of gas to cool in small objects is limited by the resolution of the simulation. This is clearly a crude approximation. The one virtue of the strategy we have adopted is that the resolution limit can be adjusted to ensure that approximately the right amount of gas cools by the end of the simulation. We have done this for the SCDM simulation and kept the same resolution for the $`\mathrm{\Lambda }`$CDM simulation with the result that more gas cooled in the latter case than is observed in the form of stars in the real universe. Clearly, further progress will require more realistic, resolution-independent and physically-based treatments of star formation and feedback.

Another limitation of our simulations, stemming also from limited resolution, in this case spatial resolution, is the inadequate modelling of the internal structure of galaxies. It seems clear that a gravitational softening of $`10h^1\mathrm{kpc}`$ is much too large to prevent the artificial disruption of galaxies in rich clusters. Our strategy of retaining gas that has cooled into galaxies as a gaseous component rather than turning it into stars as is sometimes done, partially compensates for this, but it is clear that simulations with much higher spatial resolution are required. Unlike the problem of star formation and feedback, this is a limitation that can, in principle, be overcome by increased computing power. For example, a simulation similar to ours but with a gravitational softening of only $`2h^1\mathrm{kpc}`$ would require $`100,000`$ timesteps at least for the particles in the densest regions.

In spite of these limitations, our simulations demonstrate the potential of direct simulation for realistic modelling of galaxy formation. This study is the first to follow the formation and evolution of a large number of galaxies in a volume, a cube of side $`100\mathrm{Mpc}`$, which may be regarded as representative. Over 2000 galaxies brighter than $`L_{}`$ formed in each of our simulations. Both in this paper and elsewhere (Pearce et al. 1999) we have shown that the galaxies that form in the simulations have a spatial distribution that is consistent with observations not only locally, but also at high redshift.

Here, we have used our simulations to illustrate the hierarchical build up of galaxies and to give an indication of the rate at which gas is expected to cool into fragments which subsequently merge to form larger galaxies. The aggregate star formation rate in the fragments that end up in a massive galaxy today adds up to a large value and is sharply peaked at high redshift. By contrast, the star formation rate of a typical galaxy today is more modest and spread out in time. The mean star formation rate in the simulation as a whole rises at early times, reaches a broad maximum between $`z=2`$ and $`z=1`$ and declines towards the present, in a manner reminiscent of the data of Steidel et al. (1999). Although this behaviour is, to a large extent, determined by resolution effects, the fact that it resembles the data suggests that the bulk properties of the model galaxy population today are not too unrealistic. In the $`\mathrm{\Lambda }`$CDM simulation, the first substantial galaxies form at $`z5`$. The galaxy population builds up rapidly until $`z=1`$ and there is a marked decline in the rate of change of the galaxy mass function after $`z=0.5`$. In the SCDM simulation there is more evolution at recent times. All these trends have been previously emphasized in semi-analytic models of galaxy formation (e.g. Cole et al. 1994, Kauffmann 1996). At the present day both simulations match the bright end of the observed K-band luminosity function for suitable values of the mass-to-light ratio.

The galaxy autocorrelation functions evolve little with redshift over the range $`z=03`$. The galaxy correlation function agrees closely with that of the dark matter at separations of $`10h^1\mathrm{Mpc}`$ at $`z=0`$ but differs significantly on smaller scales. For both, the $`\mathrm{\Lambda }`$CDM and SCDM models, the galaxy correlation function appears closer to a power law than the mass correlation function. The amplitude of the galaxy correlation function increases with the mass of the galaxies. The projected pairwise velocity distribution of the galaxies is significantly lower, particularly in $`\mathrm{\Lambda }`$CDM , than that of the dark matter.

Within the most massive dark matter halos the galaxies trace the dark matter faithfully. A virial mass estimator applied to the clusters correctly infers the appropriate dark matter mass albeit with a large dispersion for a single determination. However, galaxy formation is more efficient in such clusters than in the simulation as a whole, suggesting that the cosmic density parameter cannot safely be inferred from the mass-to-light ratios of clusters and the mean luminosity density of the Universe.

## Acknowledgements

We would like to thank Shaun Cole, Eric Tittley, Scott Kay and Andrew Benson for providing analysis software and data. The work presented in this paper was carried out as part of the programme of the Virgo Supercomputing Consortium (http://star-www.dur.ac.uk/ frazerp/virgo/) using computers based at the Computing Centre of the Max-Planck Society in Garching and at the Edinburgh Parallel Computing Centre. This work was supported in part by grants from PPARC, EPSRC, and the EC TMR network for “Galaxy Formation and Evolution”. CSF acknowldeges a Leverhulme Research Fellowship and PAT holds a PPARC Lecturer Fellowship. HMPC is supported by NSERC of Canada.