Including star formation and supernova feedback within cosmological simulations of galaxy formation

We investigate phenomenological models of star formation and supernova feedback in N-body/SPH simulations of galaxy formation. First, we compare different prescriptions in the literature for turning cold gas into stars neglecting feedback effects. We find that most prescriptions give broadly similar results: the ratio of cold gas to stars in the final galaxies is primarily controlled by the range of gas densities where star formation is allowed to proceed efficiently. In the absence of feedback, the fraction of gas that cools is much too high resulting, for example, in a K-band luminosity function that is much brighter than observed. This problem is ameliorated by including a feedback model which either imparts radial kinetic perturbations to galactic gas or directly reheats such material and prevents it from cooling for a certain period of time. In both these models, a significant fraction of cold gas is heated and expelled from haloes with an efficiency that varies inversely with with halo circular velocity. Increasing the resolution of a simulation allows a wider dynamic range in mass to be followed, but the average properties of the resolved galaxy population remain largely unaffected. However, as the resolution is increased, more and more gas is reheated by small galaxies; our results suggest that convergence requires the full mass range of galaxies to be resolved.


INTRODUCTION
The first detailed modelling of galaxy formation more than twenty years ago (White & Rees 1978) revealed a significant difficulty in hierarchical clustering cosmologies. These assume that galaxies form when gas cools and condenses in dark matter haloes. However, at early times, pre-galactic gas clouds are so dense that the cooling time is extremely short. Consequently, in the absence of other processes, baryonic matter cools catastrophically in sub-galactic haloes and forms stars, leaving no gas available to make up the intergalactic medium observed at high redshift, or the intracluster medium observed at low redshift. White & Rees proposed a solution to this problem: energy released from stars in the course of their evolution would act as negative feedback on the gas, limiting its cooling rate and associated star formation.
Recent modelling of galaxy formation, including the effects of feedback, has made use of two related techniques, semi-analytic modelling and direct numerical simulation. The first is a direct descendant of the methodology developed by White & Rees (1978) and provides a computationally efficient way to calculate a detailed model of galaxy formation. The merger history of dark matter haloes is followed using either statistical techniques (e.g. Kauffmann et al. 1994;Cole et al. 1994Cole et al. , 2000Somerville & Primack 1999) or N -body simulations (e.g. Kauffmann, Nusser & Steinmetz 1997;Benson et al. 2001;Somerville et al. 2000), the evolution of gas is treated by means of a spherically symmetric cooling flow model, star formation and the effects of feedback are parameterised according to simple, observationally motivated rules and stellar evolutionary effects are incorporated by means of a stellar population synthesis model. Chemical evolution and dust extinction or emission are readily incorporated in this scheme (e.g. Granato et al. 2000). Semianalytic modelling has proved to be extremely successful for the interpretation of observational data at both low and high redshift, as reviewed recently by Baugh et al. (2001).
In this paper, we are concerned with the treatment of star formation and feedback in gasdynamical simulations of galaxy formation. Numerical simulations have the advantage that some of the assumptions made in the semi-analytic treatment of gas dynamics can be relaxed. In particular, no artificial symmetries are imposed on a situation which is inherently asymmetric. However, the application of simulation techniques to galaxy formation is still in its infancy because it places strenuous demands on computational resources, consequently limiting the dynamic range of resolved structures. Early work, while making brave first attempts at simulating cosmological galaxy formation, was hampered to varying degrees by relatively poor resolution (e.g. Cen & Ostriker 1992;Katz, Hernquist & Weinberg 1992;Evrard, Summers & Davis 1994;Katz, Weinberg, & Hernquist 1996;. It is only over the last few years that parallel supercomputers have enabled simulations large enough to resolve present-day bright galaxies adequately within statisticallyrelevant volumes (e.g. Pearce et al. 1999;Fardal et al. 2000;Pearce et al. 2000). Detailed tests of varying assumptions and parameter values in simulations based on the Smooth Particle Hydrodynamics technique are presented in Kay et al. (2000 It is unfeasible to consider an explicit treatment of star formation and feedback within a cosmological simulation because these processes occur on much smaller scales than can be resolved with present techniques. The only approach (which is currently also used in semi-analytic models) is to develop phenomenological prescriptions and test their validity by comparing their effects to as many observables as possible. The main aim of this paper is to compare the outcome of adopting various star formation and feedback prescriptions in a series of small test simulations, similar to those analysed by K2000. We assume that feedback is produced only as a result of energy being injected into gas through supernovae occurring in the star forming regions that develop in the simulations. We also investigate the effects of varying the numerical resolution, albeit over a small dynamic range. Our primary concern is the effect of feedback on the properties of the resolved galaxy population, mainly at z = 0, but we also examine the properties of reheated gas, since it is these predictions that may provide some of the most stringent observational tests.
We begin with a short description of our numerical methods in Section 2. We compare various star formation models in Section 3 before investigating feedback in Section 4, focusing on both the galaxies and the distribution of reheated gas. The sensitivity of our results to numerical resolution is explored in Section 5. Finally, we draw conclusions in Section 6.

The code
The simulations were carried out using the hydra code † (Couchman, Thomas & Pearce 1995). This combines an adaptive particle-particle/particle-mesh (AP 3 M ) algorithm to calculate gravitational forces (Couchman 1991) with Smoothed Particle Hydrodynamics (SPH;e.g. Benz 1990;Monaghan 1992 and references therein) to estimate hydrodynamical quantities. The SPH implementation differs from that described in Couchman, Thomas & Pearce (1995), as discussed by Thacker & Couchman (2000). In addition to these changes, the SPH neighbour search algorithm in hydra has been modified in order to avoid setting an artificial lower limit on the gas density, present in older versions of the code (which used a single chaining mesh to find neighbours for both the PP gravity and SPH calculations). In the new version, a separate neighbour list for the gas particles is constructed using the algorithm described in Appendix A2 of Theuns et al. (1998).
Radiative cooling is modelled in the manner described by Thomas & Couchman (1992), except that tabulated cooling rates (as a function of temperature and metallicity) from Sutherland & Dopita (1993) are implemented. This assumes that the gas is in collisional ionization equilibrium and is applicable only to radiative processes at temperatures above 10 4 K . Contributions to the cooling rate at temperatures below 10 4 K are ignored, as are the radiative heating contributions from a photo-ionizing background. The former are only important in structures on scales significantly below the resolution of our simulations. Although photo-heating is important in determining the structure of the low density gas responsible for the Ly-α forest, previous studies have shown that the process is unimportant for galaxies forming in haloes with circular velocities above ∼ 50 km s −1 (e.g. Quinn, Katz & Efstathiou 1996;Navarro & Steinmetz 1997;Weinberg, Hernquist & Katz 1997), as is the case in our simulations.
The gas particle mass was held fixed for all simulations within each cosmological model, at ∼ 5 × 10 8 h −1 M⊙ in the case of SCDM and to ∼ 3×10 8 h −1 M⊙ in the case of ΛCDM. Both models assumed an initial (equivalent Plummer) softening length of ∼ 20 h −1 kpc which was fixed in comoving coordinates until z ∼ 1, after which it was fixed in physical coordinates to 10 h −1 kpc until z = 0. The SCDM simulations were started at a redshift of 24 and the ΛCDM simulations at a redshift of 49 and both sets were run to z = 0.

STAR FORMATION
In this section we concentrate on comparing several star formation models drawn from the literature. We present results for SCDM simulations with 32768 gas and 32768 dark matter particles, initially perturbed from a regular grid, within a 10 h −1 Mpc comoving box. The initial conditions and parameters are identical to those of the fiducial simulation of K2000, where further details may be found. Although the box-size is too small to allow statistically robust results for the galaxy population, the simulations are only used for comparative purposes.

Star formation models
We now summarise each star formation model studied. In all cases, the basic procedure is the same: at the end of each timestep, a subset of gas particles is identified as being eligible to form stars. A fraction of them are then converted into stars, which are subsequently subject only to gravitational interactions. Details of each method are given in Table 2.

The Pearce method
This is based on the method used by Pearce (1998). Potential star particles are identified with gas particles having overdensity, δ > 550 and temperature, T < 10 5 K . We carried out two simulations, both using the same temperature threshold as Pearce but one with an overdensity threshold of 5000 (labelled pearce1) and the other with a threshold of 50,000 (labelled pearce2). Our chosen overdensity thresholds are higher than that used by Pearce, who selected a value in order to ensure that a certain fraction of gas was converted into stars by z = 0. However, Pearce's value was also close to the maximum resolved SPH gas overdensity in his simulations. This limit is estimated by considering a sphere of radius 2hmin containing NSPH neighbours. (The minimum SPH smoothing length, hmin is set equal to 1.17 times the equivalent Plummer softening length.) For our choice of parameters (2hmin = 0.0234 h −1 Mpc and NSPH = 32), this limit is δ ∼ 20, 000 at z = 0. We deliberately set our choice of overdensity thresholds to straddle this resolution limit.

The Summers method
In the method used by Summers (1993) gas particles are converted into stars if they satisfy a temperature criterion, T < 3 × 10 4 K and a physical density criterion, nH > 0.1 cm −3 . Additionally, eligible gas particles must be within overdense regions, δ > 10, so that star formation does not occur within shallow potential wells at high redshift. We performed two simulations using the Summers method, labelled summers1 and summers2 respectively. For summers1, we used a density criterion, ρ > 10 −25 g cm −3 (ρ = mH nH/X; we assume X = 0.76), approximately a factor of 2 lower than used by Summers. This threshold corresponds to a baryon overdensity, δ ∼ 3.5 × 10 5 at z = 0 which is over an order of magnitude higher than our resolution limit. Because of this, we carried out a second simulation (summers2) for which we chose a density threshold of ρ > 10 −26 g cm −3 .

The Navarro & White method
We have also studied the method of Navarro & White (1993) who applied it to follow star formation in disk galaxies. This is a popular method that has since been used by several other authors, with minor modifications (e.g. Steinmetz & Müller 1995;Tissera, Lambas, & Abadi 1997;Carraro et al. 1998;Buonomo et al. 2000). To summarise, regions of star formation are identified in which the gas is in a convergent flow (∇.v < 0) and has a density, ρ > 7 × 10 −26 g cm −3 . Gas particles satisfying these criteria are split into equal mass gas and star particles, after a time interval, ∆t * = ln2 τ dyn , where τ dyn ∝ 1/ √ ρ is the local dynamical time of the gas. When the gas mass is less than 5 per cent of its original value, the remaining fraction is converted into stars. We label the simulation using this method navarro.

The Katz, Weinberg & Hernquist method
The fourth method investigated is that used in the cosmological simulations of Katz, Weinberg & Hernquist (1996; see also Katz 1992). Four criteria are considered to decide if a gas particle is eligible to form stars: the divergence criterion, as used in the Navarro & White method; the physical density threshold used as in the Summers method; an overdensity threshold, δ > 55.7, and a Jeans instability condition, τJ > τ dyn , where τJ ≡ h/cs is a "Jeans" timescale, h is the particle's SPH smoothing length and cs = (γ(γ − 1)u) 1/2 ‡ is the adiabatic sound-speed of the gas. (Note that this condition depends on the resolution of the simulation since h is limited by the gravitational softening.) If a gas particle is eligible to form stars, a third of its mass is converted if r < 1 − exp(−∆t/min(τ dyn , τ cool )), where r is a random number between 0 and 1, ∆t is the timestep and τ cool ∝ 1/ρ is the cooling time. This method gives each particle a dual identity (we refer to it as a schizophrenic particle), allowing it to carry both a gas mass and a stellar mass. Gravitational forces are calculated using the total mass and hydrodynamical forces using only the gas mass. Once the gas mass has dropped below 5 per cent of its original value, the particle becomes a purely collisionless star particle. We refer to the simulation using this method as katz.

The Mihos & Hernquist method
The final model studied is based on the method used by Mihos & Hernquist (1994). This uses the Schmidt law (Schmidt 1959; see also Kennicutt 1998 and references therein), an empirical (power-law) relation between the star formation rate per unit area and the surface density of Hi gas in late-type galaxies. Mihos & Hernquist convert this to a 3-dimensional relation, such that the fractional change in mass converted into stars during a timestep is ∆M * /Mgas = C sfr ρ N−1 ∆t. This equation is applied to gas particles with δ > 10 and ∇.v < 0, although rather than converting part of a particle's mass into stars, we convert all of it if ∆M * /Mgas > r, where r is a random number between 0 and 1. We study a simulation using this method, labelled mihos, in which we have set N = 3/2, consistent with the observed value, and C sfr = 5 × 10 −4 cm 3/2 g −1/2 s −1 .

Results
3.2.1 General properties at z = 0 Table 3 lists general properties of the 7 simulations at z = 0, and of the equivalent simulation without star formation (labelled gas). Column 2 gives the CPU time required to complete each run, relative to gas. (To provide a fair compar- ‡ u = kT /(γ − 1) µm H ; we take the ratio of specific heats at constant volume and constant pressure for a monatomic ideal gas to be γ = 5/3, and the mean molecular weight of a gas of primordial composition to be µ = 0.59 ison, differences in the machine-dependent clock rate were removed by dividing the total time taken for each simulation by the average time per step over the first 10 steps.) The least efficient method is navarro, which takes around a third longer to complete than gas because of the increased number of particles produced by star formation. The rest of the simulations all take around half of the time to complete than gas, with the shortest being pearce1. Including a prescription for star formation reduces the CPU requirement because of the reduced SPH workload. Column 3 gives the fraction of baryons at z = 0 in stars or gas with ρ > 500 ρ and T < 30, 000K , labelled f cool . Gas that satisfies these criteria (i.e. it is cold and dense) is readily associated with galaxies, since the two thresholds enclose a region of the ρ−T plane containing material which has radiatively cooled (see, for example, §3.2 of K2000). For all simulations, including gas, fg = 0.34 ± 0.01. As found by Pearce (1998), the conversion of cold dense gas into stars does not significantly affect the amount of gas that cools.
Column 4 gives the total fraction of baryons in stars, f * , for each run. The range of f * is much wider than for f cool , varying from 15 per cent (katz) to 28 per cent (pearce1). Increasing the overdensity threshold from 5000 ( pearce1) to 50,000 (pearce2) decreased f * by only 0.03. This small difference is due to the extra time it takes for cooled material to collapse to a higher overdensity before forming stars. Summers2 produced a similar fraction of stars as pearce1 and pearce2, but 40% more than summers1. The difference in the summers simulations is also due to the choice of density threshold which, at z = 0, is about an order of magnitude higher for summers1 than for pearce2. katz and navarro produce similar values of f * (to within 20%) to summers1 because they use the same density threshold. The number of stars forming at any given time in mihos is controlled by the choice of normalisation constant, C sfr . The value, C sfr = 5 × 10 −4 cm 3/2 g −1/2 s −1 , was chosen (by trial and error) to give a similar fraction of stars at z = 0 as pearce2.
The final column in Table 3 lists the number of galaxies in each simulation. These were identified using a friends-offriends group-finder (Davis et al. 1985) on the cold gas and star particles (i.e. on those used to calculate f cool above), with a linking length, b = 0.1 (comparable to the gravitational softening length.) Any objects with fewer star particles than NSPH (32 for the simulations studied in this paper) were discarded. All runs contain ∼ 70 galaxies at z = 0 to within 10%.  . The cumulative fraction of cooled gas (ρ > 500 ρ and T < 30, 000K ) for the simulation without star formation, is also plotted. The instantaneous star formation rate is proportional to the gradient of each star fraction curve.

Star formation histories
All simulations started cooling material and forming stars within the first 2 Gyr (z > 2). The exact time when the first stars form depends upon two factors: the resolution of the simulation, which determines the formation time of the first resolved dark matter haloes, and the range of gas densities where star formation occurs. Since the resolution is fixed in all our simulations, only the second factor is relevant.
The star formation history in pearce1 is similar to that of cooled material, i.e. it rises sharply from z = 5 until z ∼ 2, before rolling over, implying a decreasing star formation rate at lower redshift. However, in pearce2, the delay caused by a higher overdensity threshold causes star formation to begin later and to proceed at a slower rate until around 5 Gyr (z ∼ 1), after which the z = 0 normalised star formation rate is higher than in pearce1. summers1, navarro and katz have the same histories at all times, demonstrating that the additional features of the latter two models are irrelevant in determining the average star formation rate within the simulations. (We exclude the latter two simulations from further discussion since they produce nearly identical results to summers for all quantities studied in this paper.) These curves are very similar to those from the pearce simulations until z = 1 − 2, where a "kink" develops, due to the behaviour of the gravitational softening. At early times, the softening is constant in comoving coordinates, and so the minimum resolved length scale in physical coordinates grows in proportion to the expansion factor, a = 1/(1 + z), leading to a decrease (as (1 + z) 3 ) in the maximum resolved physical density. By z ∼ 1.7, this density has dropped to ρ ∼ 10 −26 g cm −3 , the summers1 density threshold. This suppresses further star formation until z = 1, when the softening becomes constant in physical coordinates, allowing the star formation rate to increase again. summers2 exhibits a similar star formation history to pearce1. In these cases no kink is evident since the density thresholds are always lower than the maximum resolved density.
Finally, in mihos, the relative fraction of stars rises more sharply than other curves at high redshift but flattens off by z = 0. For this run, an approximate minimum density threshold can be calculated by first equating the relevant star formation timescale, τ sfr = Mgas∆t/∆M * = 1/C sfr ρ 1/2 gas (c.f. §3.1.5) to the lookback time for an Ω = 1 universe, τ look = (2/3H0)(1 − (1 + z) −3/2 ); at higher densities, τ sfr < τ look . Assuming that δ ∼ ρ/ ρ (i.e. for δ >> 1), the minimum overdensity scales as δmin ∝ ((1 + z) 3/2 − 1) −2 , a decreasing function of redshift. At high redshift, this overdensity threshold was sufficiently low to allow efficient star formation to take place. However, at lower redshift, the decrease in the density of the Universe and in the available time left to form stars before z = 0, causes the number of gas particles above this threshold to drop rapidly.

Cold gas fractions
In Fig. 2, we plot the median mass ratio of cold gas to stars against the stellar mass of each galaxy with 32 or more star particles and a non-zero gas mass, in bins of width, ∆logM * = 0.35. Error bars represent 10 and 90 percentile points of the distribution in each bin. There is a clear trend for more massive galaxies to contain proportionally less cold gas. This trend reflects the distribution of mean stellar ages with galaxy mass: the more massive galaxies in all star formation models harbour, on average, older stellar populations than less massive galaxies. Therefore, cooled baryons have had more time to collapse and cross the imposed threshold for star formation.
The cold gas to star mass fraction in the brightest galaxies (above M * = 10 11 h −1 M⊙, approximately the mass of an L * galaxy) is less sensitive to the choice of star formation prescription than for less massive galaxies. However, below M * = 10 11 h −1 M⊙, simulations with a higher density threshold (e.g. summers2) contain significantly more cold gas than simulations with a lower threshold. For example, the median ratio of cold gas to stars at M * = 2.4 × 10 10 h −1 M⊙ is ∼ 0.4 for summers2 but is 1.6 for summers1. Again, this trend reflects the reduction in star formation activity when the density threshold is high.
Observations of atomic and molecular gas in disk galaxies can, in principle, be used to place constraints on plausible star formation models. However, this can only be done in a rather crude manner because the simulations do not have enough resolution to distinguish between morphological types. By selecting those galaxies with cold gas we are at least able to identify candidate disk galaxies in the simulations. From the data plotted in Fig. 9 of Cole et al. (2000), L * disk galaxies (M/LB ∼ 2 and M * B ∼ −19.5) have Mcg/M * ∼ 0.1 − 0.3, values which do not increase significantly for galaxies an order of magnitude less massive. This is broadly consistent with our results for simulations with low density thresholds, but inconsistent with those with high thresholds (extrapolating the result for summers1 down to M * = 10 10 h −1 M⊙ implies a ratio of at least 2).

Galaxy mass and luminosity functions
Cumulative mass functions for the galaxy population at z = 0 are plotted in Fig. 3. For comparison, the mass function for simulation gas is also displayed. All the mass functions are very similar, demonstrating that the abundance of galaxies over the mass range resolved in these simulations is insensitive to the inclusion of star formation.
We can assign luminosities to the galaxies in the simulations using a stellar population synthesis model. We calculate K-band luminosities since the light in this band is less affected by dust than at other bands and, since it is dominated by late-type stars, it provides a more robust estimate of the underlying stellar mass of the galaxy. To calculate galaxy luminosities, the formation time, t f , of each star particle, m * , was used to calculate the luminosity per unit mass at the present day, l(t0 − t f ). The total magnitude is then  where MK(0) is the zero-point and N * is the number of star particles in the galaxy. The function, l(t0 −t f ) was generated using the stellar population synthesis models of Bruzual & Charlot (1993), assuming solar metallicity and the stellar IMF proposed by Kennicutt (1983; see also Cole et al. 2000). The total magnitude of each galaxy can be rescaled to take into account non-visible stellar material (i.e. a population of brown dwarfs) by adding the term 2.5logΥ, where Υ ≡ 1 + M (brown dwarfs)/M (visible stars) (Cole et al. 2000).
Cumulative luminosity functions for galaxies in our simulations at z = 0 are shown in Fig. 4. Only galaxies with 32 or more star particles (or equivalent) are included. The best-fit Schechter function to the K-band galaxy luminosity function from the 2dF and 2-Mass surveys is also plotted as a solid curve , with the value of M * K marked with a vertical arrow. The model luminosity functions are plotted assuming Υ = 1; the shift required for simulation pearce1 to agree with the data at M * K is indicated by a horizontal arrow. This implies a value of Υ = 3, i.e. twice as much mass is required in brown dwarfs than in visible stars, which is in clear contradiction with observations (see Balogh et al. 2001).

Summary
We have examined various prescriptions in the literature for turning gas into stars in cosmological simulations. The determinant factor in each model is the range of density (usually defined by a lower limit, or threshold) at which star formation is allowed to proceed efficiently: simulations with lower thresholds naturally form more stars than simulations with higher thresholds. The cooling rate of the gas is largely unaffected by the inclusion of a star formation mechanism. In these models, the star formation efficiency determines the fraction of galactic baryons that end up in the form of cold gas or stars. Simulations with a high density threshold contain more cold gas in (particularly sub-L * ) galaxies than simulations with low density thresholds; the former are discrepant with observational data. All of the models overpredict the abundance of galaxies at a given magnitude, or equivalently produce galaxies that are (about 0.5-1 magnitudes) too bright at a given abundance. In the next section, we consider feedback mechanisms that act to reduce galaxy luminosities.

FEEDBACK
In this section we study the effects of energy injection associated with Type II supernovae explosions in a cosmological simulation. On energetic grounds alone, these events must be important during galaxy formation, and their role as a negative feedback mechanism allows ab-initio hierarchical models to match the observed galaxy luminosity function (White & Rees 1978;White & Frenk 1991;Cole 1991). Ideally, we would like to be able to resolve individual supernova events and model the multi-phase intergalactic medium in our simulations. However, this is impractical, both because of computational limitations and because of the lack of a detailed understanding of the processes involved. At present, the only practical strategy is to include the effects of supernova feedback phenomenologically.
Phenomenological feedback models within numerical simulations have already been considered by a number of authors studying galaxy formation/evolution. For example, supernova feedback was first included within SPH simulations by Katz(1992), who studied the collapse of a galaxysized density perturbation, and used a scheme where every star particle supplied neighbouring gas particles with thermal energy (heat). However, this had very little effect on the galaxy's properties, because the high gas density at the locations where energy is injected implies a cooling time which is much shorter than the dynamical time. As a result, the energy injected is simply radiated away. For future reference we will refer to this scheme for energy injection as "thermal feedback".
An alternative feedback scheme was proposed by Navarro & White (1993). In this model, a fraction, fv, of supernovae energy is supplied to neighbouring gas particles in the form of kinetic energy, with the rest being added as heat. Like Katz, Navarro & White found that for low values of fv (i.e. almost all thermal energy), very little feedback occurs and consequently the star formation rate is very high. However, for higher values of fv the star formation rate is reduced because kinetic perturbations are able partially to halt the local collapse of a density perturbation. We refer to this method as "kinetic feedback". Kinetic feedback has since been used by other authors (e.g. Mihos & Hernquist 1994).
Subsequent attempts have been made to boost the effects of thermal feedback by preventing the gas from immediately re-radiating the supernova energy. This was first considered by Gerritsen (1997) who argued that previous feedback schemes did not produce an overall disk morphology resembling those observed, with spherical cavities, or bubbles, filled with hot gas. They found that injecting all the energy into one particle as heat and preventing it from cooling for a finite length of time, had the desired effect. A more elegant approach was taken by Springel (2000), who created a second pressure term in the equations of motion, representing turbulent pressure caused by supernovae. Yet another example is the study by Thacker & Couchman (2000), who considered a scheme in which a "cooling" density for each gas particle injected with energy is obtained from a simple pressure equilibrium argument. This density is lower than the standard SPH density, and so lengthens the particle's cooling time, before resorting back to the SPH density after a finite period.
Attempts to simulate galaxy formation in a large region, starting from realistic cosmological initial conditions, and including a model for feedback have been made by several groups using Eulerian grid codes (e.g. Cen & Ostriker 1992Yepes et al. 1997) and SPH codes (Summers 1993;Katz, Weinberg, & Hernquist 1996). The former technique is hindered by poor spatial resolution, with galaxies being several factors smaller than the minimum resolved scales (although this situation is improving, particularly with adaptive mesh codes). Studies using SPH simulations have focussed on thermal feedback, and like Katz (1992), have also found that adding thermal energy to the gas makes little difference to the resulting galaxy population. In the next subsection we will describe our own implementations of thermal and kinetic feedback in our simulations, which we regard as a natural follow-up to the work of Summers (1993) and Katz et al. (1996). As we shall see, our models are capable of producing strong feedback, reducing the amount of material which cools and forms stars, and resulting in a more realistic K-band luminosity function.

Cooling and Star formation
We have studied the effects of feedback in the ΛCDM and SCDM cosmological models described in § 2.2. Similar trends with respect to variation of the parameters are present in both models, so we will focus mainly on ΛCDM. We implement the "decoupling" technique described by Pearce et al. (2000) whereby gas particles with temperature below T = 12, 000K are ignored when estimating the gas density of hot particles with temperature above T = 10 5 K . This technique effectively prevents the unphysical cooling of hot gas induced purely by the presence of nearby cold, dense galactic gas. (Note, however, that this problem is largely circumvented due to the inclusion of a star formation prescription which removes the majority of galactic material from the SPH calculations.) This modification is a crude attempt to model a multiphase gas with SPH and has been discussed further elsewhere (e.g. K2000; Croft et al. 2000; see also Ritchie & Thomas 2001).
Star formation was incorporated using the Pearce method (c.f. § 3.1.1), assuming an overdensity threshold of 5000 and a temperature threshold of 30, 000K (essentially the pearce1 criteria). When a gas particle satisfies both criteria it is instantaneously converted into a star particle. In the ΛCDM simulations we assumed an unevolving gas metallicity of Z = 0.3Z⊙ while for the SCDM simulations, we retained the metallicity of Z = 0.5Z⊙ assumed in the previous section. Note that our choice of metallicity results in an overestimate of the cooling rate of gas at high redshift, where Z is considerably smaller. Self-consistent models of feedback and metal enrichment within simulations will be presented in a future paper.

Supernova energetics
We adopt the same parameter values for the supernovae energetics as Katz (1992) and Katz et al. (1996). To summarise, the amount of energy available from supernovae is calculated assuming that stars with mass above 8M⊙ release 10 51 erg of energy into the surrounding gas. A Miller & Scalo (1979) IMF is assumed, with lower and upper mass limits of 0.1M⊙ and 100M⊙ respectively, implying a specific energy released in supernovae of ∼ 3.7 × 10 15 erg g −1 , equivalent to a temperature of 1.8 × 10 7 K . The size of the timestep in the simulations (typically 3 Myr) is about one order of magnitude smaller than the lifetime of an 8M⊙ star (τ8 ∼ 20Myr). For this reason, the energy was released gradually, using the following functioṅ uSN(t) = 3.7 × 10 15 exp (−(t − t0)/τ8) τ −1 where t0 is the time when the star formation event occurred. Energy was released until t − t0 = 200Myr, whenuSN becomes negligible.

Feedback models
Our prescription for thermal feedback works as follows. When a gas particle is to be converted into a star particle, it first passes through an intermediate phase, in which it is labelled a supernova particle. Supernova particles feel the gravitational force only, but the SPH algorithm maintains both their densities and neighbour list. Smoothing lengths are fixed at the minimum value, hmin, and each gas neighbour within the SPH smoothing radius (2hmin) receives a proportion of the thermal energy available at that time. For supernova particle i, the energy given to gas particle j, over the timestep ∆t, is where 0 ≤ ǫ ≤ 1 is the feedback efficiency parameter, and H is the Heaviside step function. We choose not to weight the contributions using the SPH kernel since this as an unnecessary complication in cosmological simulations, where the forces within galaxies are softened. Furthermore, the efficiency of energy transfer will vary depending on the number of neighbours of each supernova particle. Our method ensures that the same amount of energy per supernova is always transferred to the gas and so the efficiency is controlled solely by the choice of ǫ.
As discussed above, this scheme is expected to have little effect on the galaxy properties, since the gas will reradiate the energy in a short time. To circumvent this problem, we follow the method of Gerritsen (1997) and prevent the gas from cooling over a period of time. This timescale is controlled using the parameter τ , which is the length of the non-radiative period in units of the age of the universe at z = 0.
We also study kinetic feedback, where instead of heating the gas, it receives a velocity perturbation radially away from the supernova particles. The total magnitude of the velocity added to gas particle j, from supernova particle i, is This scheme is the "momentum" implementation of Navarro & White 1993.

Results
We present results from 6 different parameter choices for the thermal feedback model and 3 for the kinetic model. Specifically, for simulations with thermal feedback, we consider τ = 0.01, 0.1, 1.0 for ǫ = 0.1, 1.0. For simulations with kinetic feedback, we consider ǫ = 0.001, 0.01, 0.1. Details of the feedback models for ΛCDM and SCDM simulations are summarised in Table 4.

General properties at z = 0
We identify galaxies by applying a friends-of-friends groupfinder to the gas and star particles as in Section 3.2.1, except that we use a linking length b = 0.2. This gives a slightly larger number of galaxies than b = 0.1. (In simulations without feedback, both values of b pick out the same galaxies, but when feedback is included a larger linking length is required because the galaxy material becomes slightly distended by the disturbance produced by the injection of supernova energy. Since the internal structure of galaxies is unresolved, the exact choice of linking length is unimportant provided it captures the majority of galaxy material.) We also demand that each galaxy should contain a minimum of 32 star particles, regardless of the number of gas particles, in order to provide a reasonable estimate of its luminosity. The minimum baryonic mass of a galaxy in these simulations is 1.0×10 10 h −1 M⊙ for the ΛCDM model and 1.6×10 10 h −1 M⊙ for the SCDM model. Column 5 of Table 4 gives the fraction of baryons (cold gas and stars) in galaxies, labelled f gal . For simulations with the same feedback parameters, the fraction of galaxy material is always higher in the ΛCDM case than in the SCDM case though the trends with varying ǫ and τ are the same. A larger feedback efficiency, ǫ, causes a decrease in the fraction of galactic material at z = 0 (for both thermal and kinetic feedback). In the thermal case, increasing the efficiency increases the pressure of reheated gas particles, with the result that more material can rise out of galaxies, before it can recool and form stars. In the kinetic case, a larger ǫ causes an increase in the fraction of gas particles with velocities greater than the escape velocity of each galaxy. For a given value of ǫ in the simulations with thermal feedback, a larger τ results in a smaller f gal . This is because the longer the reheated gas is prevented from cooling, the easier it is for this gas to escape from the galactic potential. (Note that for the largest value of τ considered, τ = 1, reheated gas can never re-radiate its energy before the end of the simulation.) The number of resolved galaxies in each simulation at z = 0 is listed in Column 6. This number is sensitive to the feedback mechanism. In general, SCDM simulations contain more galaxies than the equivalent ΛCDM simulations, although the exact values depend on both the mass resolution of the simulation and the shape of the galaxy mass function near this threshold.

Star formation rates
In Fig. 5, we illustrate the effect of feedback on the global star formation rate, by plotting the cumulative fraction of baryons in stars as a function of time for ΛCDM models.
In each panel we compare results between simulations with the same feedback prescription, but with different values of one parameter. For thermal feedback, dotted lines represent results from simulations with τ = 0.01, short-dashed with τ = 0.1 and long-dashed with τ = 1; for kinetic feedback, the same sequence represents results for ǫ = 0.001, 0.01 & 0.1. The equivalent simulation without feedback is also plotted as a solid line. Consider first the case of thermal feedback. For ǫ = 0.1, the efficiency of star formation is sensitive to the value of τ (when varied over 2 orders of magnitude). However, increasing ǫ by a factor of 10 produces only small differences in the star formation rate when τ is varied. In the latter case, enough energy is deposited into the gas that even on the shortest timescale considered, τ = 0.01 (of order 100 Myr), much of the reheated material is unable to cool and form stars before the end of the simulation. Kinetic feedback can be very effective at reducing the amount of star formation, even though gas above 10 4 K is able to cool radiatively at all times. Although supplying 0.1 per cent of the available Figure 6. Thermal history (density, temperature and entropy as a function of time) of a gas particle in the ΛCDM simulation with thermal feedback (ǫ = 1, τ = 1) that would otherwise have cooled and formed stars. The density threshold for star formation is plotted as a dotted line. Redshift intervals are labelled along the top of the figure. energy makes no significant difference to the star formation rate, increasing the efficiency to 1 and 10 per cent lowers the fraction of stars formed by z = 0 to around 2/3 and 1/4 of the no-feedback values respectively.

The fate of reheated gas
Virtually all the gas that is reheated, either thermally or kinetically, has a temperature in the range 10 4 − 10 5 K . The temperature to which it is reheated and its subsequent thermal evolution depends on the specific feedback mechanism. For thermal feedback, at least 80% of the remaining gas that was once reheated is hotter than 10 5 K since any gas that can cool again is rapidly turned into stars. In the case of kinetic feedback, the final temperature of the reheated gas depends on the value of ǫ: a greater efficiency leads to hotter gas. Unlike in the thermal models, gas in the kinetic models is not reheated instantaneously but has to thermalise its energy, a process which will be particularly susceptible to numerical resolution .
To investigate the fate of reheated gas, we first illustrate the typical history of a reheated gas particle that would otherwise have cooled within a galaxy and been converted into a star particle. We have selected a particle in the ΛCDM simulation with thermal feedback, for ǫ = 1 and τ = 1. The density (n/ cm −3 ), specific "entropy" (defined as s = kT /n 2/3 ) and temperature of this particle are plotted in Fig. 6. Also shown (dotted line) is the density at each redshift equivalent to the overdensity threshold for star formation, i.e. n > 5000 ρ /µmH. Initially (z > 5), the particle cools adiabatically (i.e. isentropically) from its initial temperature of 100K due to the expansion of the Universe. (Note that our choice of initial temperature is arbitrary but is too low to have any effect on the subsequent evolution of the particle. The true temperature of the gas at high redshift is determined by the intensity of the UV background, which we ignore for this study.) At z = 5 the particle has fallen into a halo where it is shock-heated to a temperature of 10 4 K : the density is sufficiently high that the particle is able to radiate any excess thermal energy within a timestep (K2000). At this point the density does not increase (as expected during shockheating) since the region is not fully resolved by the SPH algorithm. A few of the particle's neighbours are not part of the overall collapse (leading to an underestimate of the gas density) even though the net flow is converging, which results in viscous heating. Shortly after this time, however, the density of the particle starts increasing as it condenses to form part of a "galaxy" within a forming dark matter halo. (The halo, with Vc ∼ 60 km s −1 , contains approximately 25 particles and so is also poorly resolved at this stage.) The density reaches a peak value of 0.1 cm −3 at z ∼ 3. However, before the particle crosses the threshold for star formation, it absorbs supernovae energy which raises its temperature by about 2 orders of magnitude, to 10 6 K . Correspondingly, the density decreases by 3 orders of magnitude (aided by the fact that since the particle is at T > 10 5 K , its SPH density is determined only by other hot particles) and its entropy jumps to ∼ 100keVcm 2 . Apart from a few minor heating events, the particle evolves approximately isentropically for the remaining time (z < 3) as it moves into less dense environments.
In order to quantify the fate of gas reheated by supernovae energy, we plot in Fig. 7 the fraction of reheated gas particles in the ΛCDM simulations as a function of time, which lie in regions of different overdensity. The solid lines show the entire reheated fraction, the long-dashed lines the fraction with δ < 500, the short-dashed lines the fraction with δ < 50 and the dotted lines the fraction with δ < 10. The first cut approximately selects gas that has been expelled from galaxies (blow-out), the intermediate overdensity preferentially picks out gas expelled from haloes altogether (blow-away), while the final criterion indicates particles that have been ejected into weakly overdense regions such as those responsible for the Ly-α forest observed in quasar spectra.
The first reheated particles appear at around z = 5, when the first stars form and the associated supernovae inject energy into the gas. The reheated fraction increases steadily to z = 0 in all simulations except in the case of thermal feedback with ǫ = 0.1 and τ = 0.01. In this simulation, the amount of reheated gas that can cool again and form stars at z < 1 outweighs the total amount of gas that is reheated. As expected, at low redshift, the feedback models with the highest star formation rates have the lowest fraction of reheated gas. However, the amount of reheated material is lower in the kinetic feedback models at any given time, than in a thermal model which produces approximately the same fraction of stars. For example, both the thermal feedback model with ǫ = 1, τ = 1 and the kinetic feedback model with ǫ = 0.1 predict that 5-6% of the baryons are locked into stars in resolved galaxies (M > 10 10 M⊙) at z = 0. However, ∼ 34% of the baryons are reheated in the thermal case, whereas only ∼ 20% are reheated in the kinetic case. In all our simulations, a significant fraction of the gas that absorbs supernova energy escapes from galaxies and haloes. Kinetic feedback is more efficient at transporting the gas into the intergalactic medium than thermal feedback. In the latter models, as much as 10% of the baryons are reheated and expelled into weakly overdense regions (δ < 10) by z = 0, so that the ratio of diffuse reprocessed gas to stars is as high as 2. Kinetic feedback models predict a similar ratio, even though the total fraction of reheated gas is much lower.

Galactic winds
One of the effects of energy injection is to expel gas from haloes in the form of a galactic wind. The time at which the gas is expelled is only approximately known in the simulations because of the relatively sparse output times. To estimate it, we first identify dark matter haloes at the simulation output times (z = 5, 3, 2, 1, 0.5, 0). Haloes are defined as overdense spheres enclosing an average density ∆c ρcr(z), where ρcr(z) = 3H 2 (z)/8πG; we set ∆c = 178 Ω 0.45 (z) for the ΛCDM model (Eke, Cole & Frenk 1996). We then tag each particle with the latest output time when it was inside a halo, and note the corresponding halo mass. The circular velocity of a halo is related to its mass and redshift by where M14 is the total mass of the halo in 10 14 h −1 M⊙. A reheated gas particle is deemed to have been expelled from the halo when it first moves beyond the virial radius, providing it was identified with a resolved halo (32 or more particles) before the reheating event and also prior to expulsion. All expelled particles are then grouped together to calculate the expulsion mass for each halo. In Fig. 8, we plot the median mass fraction of expelled gas as a function of the circular velocity of the halo from which the gas was blown away. (We define the halo gas mass as M halo Ω b /Ω0, where M halo is the total halo mass.) Error bars illustrate 10 and 90 percentile points of the distribution within each bin. Stars, triangles and squares represent results from ΛCDM thermal feedback simulations with τ = 0.01, 0.1 & 1.0, and kinetic feedback simulations with ǫ = 0.001, 0.01 & 0.1 respectively. Although the range of resolved circular velocities is limited in these simulations, Fig. 8 shows a clear trend of decreasing expulsion fraction with increasing circular velocity for nearly all cases, with only a few per cent of gas escaping from Milky Way sized haloes (Vc ∼ 220 km s −1 ) and larger. Increasing the available energy (large ǫ) causes a larger fraction of the gas to be expelled, particularly in the haloes with the lowest values of Vc. A larger value of ǫ raises the mean temperature of the reheated gas, allowing a greater fraction to escape from the smaller systems. However, the average cooling times for these systems are considerably shorter than their dynamical times, and so the amount of gas escaping in the thermal feedback simulations is lower for smaller values of the delay timescale, τ .

Galaxy properties
We now investigate the effects of feedback on the present day galaxy population. Firstly, we examine its effect on the amount of cold gas in galaxies, by plotting the mean ratio of cold gas mass to stellar mass, Mcg/M * , against M * for galaxies (with non-zero gas content) in the ΛCDM simulations with 32 or more particles. The points in Fig. 9 represent the median values within bins of width, ∆logM * = 0.5, and error bars represent 10 and 90 percentile points within each bin. The trend of decreasing cold gas fraction with galaxy mass is similar to that seen in Fig. 2 for simulations without feedback, but all our feedback models predict median cold gas fractions below 0.5 for galaxies with M * > 10 10 h −1 M⊙. This agrees well with the observational data summarised by Cole et al. (2000;see Section 3.2.3) The overall reduction in star formation activity caused by feedback has a substantial impact on the galaxy luminosity function. Cumulative K-band luminosity functions are shown in Fig. 10 for the ΛCDM simulations. Luminosities for each galaxy (consisting of at least 32 star particles) were calculated using the method described in Section 3.2.4. The best cumulative Schechter function fit to the data from Cole et al. (2001) from the 2dF and 2-Mass surveys is plotted as a solid curve, with the value of M * K indicated with a vertical arrow.
Results from the thermal feedback simulation with ǫ = 0.1, τ = 0.01 and the kinetic feedback simulation with ǫ = 0.001 are similar to that from a simulation with no feedback. These all produce a luminosity function which is ∼ 2 magnitudes brighter than observed. As the strength of the feedback is increased (by increasing ǫ and/or τ ), the stellar mass and hence the luminosity of the galaxies are reduced. Both the thermal model with ǫ = 1 & τ = 1 and the kinetic model with ǫ = 0.1 predict about the correct number of L * K galaxies, although the latter model underestimates the abundance of fainter galaxies.

Summary
We have investigated two ways of including feedback from supernovae energy in cosmological simulations, in order to reduce the amount of gas that cools and forms stars. The first method, thermal feedback (Katz 1992), adds thermal energy to surrounding gas (with efficiency controlled by the parameter ǫ) and prevents such material from recooling for a period (controlled by the parameter τ ). The second method, kinetic feedback (Navarro & White 1993), adds velocity perturbations to the gas radially away from the stars, and allows the gas to cool at all times. Increasing ǫ and/or τ increases the strength of feedback, reducing the global star formation rate, but without significantly affecting the average ratio of cold gas to stars in galaxies. The feedback works in both models by reheating gas which has already cooled, with stronger feedback resulting in more gas escaping into regions of low overdensity. The efficiency with which gas is expelled from dark matter haloes is higher in haloes with Vc ∼ 100 km s −1 than in haloes with higher circular velocity. By limiting the star formation rate, feedback lowers the K-band luminosity of galaxies and, with a suitable choice of parameter values, it is possible to reproduce the bright end of the galaxy luminosity function.

RESOLUTION EFFECTS
The simulations discussed so far have a fixed mass resolution, with the smallest resolved galaxies having a baryonic mass ∼ 10 10 h −1 M⊙, roughly one tenth that of an L * galaxy. In the absence of feedback, resolution effects completely determine the outcome of a simulation: as smaller objects are resolved, increasingly large amounts of gas cool, raising both the global fraction of cooled material as well as the mass of the majority of the galaxies. Feedback limits the amount of gas that cools and so it seems plausible that an appropriate feedback prescription might lead to results that do not depend strongly on resolution. In this case, an increase in resolution would lead to smaller galaxies being resolved without affecting the properties of the larger objects.
To address the issue of resolution, we performed two additional simulations, using both thermal and kinetic feedback models. For thermal feedback, we set ǫ = 1 and τ = 0.1; for kinetic feedback, we set ǫ = 0.1. The cosmological parameters and initial conditions were identical to those used for the SCDM model in §3, except that the number of gas and dark matter particles were increased from 32768 (32 3 ) of each to 110592 (48 3 ). This results in a decrease in the minimum baryonic mass of a resolved galaxy by about a factor of 3, from ∼ 1.6 × 10 10 h −1 M⊙ to ∼ 4.8 × 10 9 h −1 M⊙. When constructing the initial density field, the phases of the large-scale waves were preserved, so that the same overall large-scale structure is present in all the simulations. We also increased the spatial resolution by decreasing the softening length in proportion to m 1/3 , where m is the dark matter particle mass. This corresponds to a reduction by a factor of 1.5, for the duration of the simulation, to a final (Plummer) of ∼ 6.7 h −1 kpc. All other parameters were fixed at the values used in the low resolution simulations. Table 5 compares the properties of the high resolution simulations with their low resolution counterparts. The last two columns give the fraction of baryons in resolved galaxies, f gal , and the total number of galaxies, N gal . Galaxies are identified as before, using a friends-of-friends group-finder on the distributions of gas and star particles with a linking length of b = 0.2, and only those with 32 or more star particles are retained. As expected, an increase in resolution leads to an increase in both the number of resolved galaxies and the galaxy baryon fraction.
The main differences in the simulations with low and high resolution occur at early times and in small haloes. For example, the star formation histories (Fig. 11) are quite similar except within the first gigayear, when smaller objects are resolved in the high resolution simulation. Increasing the resolution serves to extend the range of resolved galaxies to smaller masses, but has relatively small effects on more massive galaxies. For example, in Fig. 12, we plot the median ratio of cold gas to stars against stellar mass for galaxies with non-zero gas content and 32 or more particles (error bars represent 10 and 90 percentile points within each bin, of width ∆logM * = 0.5 for low resolution simulations and  Median ratio of cold gas to stars against stellar mass for galaxies with non-zero gas content and 32 or more particles, in feedback simulations with varying mass resolution. Error bars represent 10 and 90 percentile points in each bin. 0.35 for high resolution simulations). The cold gas fractions in the high resolution simulations are consistent with the low resolution simulations for galaxies resolved in both.
The extension of the resolved galaxy population to smaller masses as the resolution is increased is clearly illustrated by the K-band galaxy luminosity functions, plotted in Fig. 13. In the thermal feedback case, the luminosity Figure 13. K-band galaxy luminosity function for the simulations with thermal and kinetic feedback and varying mass resolution. Also plotted is the best-fit Schechter function to the result from the 2dF and 2-Mass surveys . The value of M * K is illustrated with a vertical arrow.
function in the higher resolution simulation matches well on to the low resolution luminosity function. In the kinetic feedback case, there is some mismatch around M * K , but the number of galaxies is so small here that these differences are not significant. We conclude that, overall, the increase in the resolution of the simulation has little effect on the population of massive galaxies.
Resolution affects are more noticeable in the amount and distribution of reheated gas. Fig. 14 shows the fraction of gas at different epochs lying in regions of different overdensity (c.f. Fig 7). Increasing the resolution results in an increase in these fractions at all redshifts. The hydrogen column density distribution of reheated gas is illustrated in Fig.15 for the four simulations, at z = 3, 2, 1 & 0. Column densities were calculated by smoothing the density of local particles onto a 64 × 64 grid, using the SPH kernel with width set by the particles' smoothing lengths. Overlaid as circles (with radii proportional to M 1/3 gal ) are the positions of galaxies with 32 or more star particles. At z = 3 the gas in the low resolution simulations is confined to the few regions where the first galaxies have formed (most of these are unresolved.) In the simulations with higher resolution the gas has a much larger covering factor because many more galaxies are resolved and so more of the volume is occupied by galaxies. The range of column densities is high in these regions (reaching ∼ 10 20 cm −2 ). At lower redshifts, the distribution of gas becomes much more skewed towards regions of higher overdensity, as the galaxies and haloes merge into larger units. By z = 0, the variation of density contrast in the gas is prominent, although the distribution is smoother in kinetic than in thermal models. As was inferred from Fig. 7, kinetic feedback is able to transport gas into regions of lower den- Figure 15. Column density of reheated hydrogen gas in SCDM simulations, with the positions of the galaxies superimposed. From top to bottom, the rows are for redshifts 3,2,1 & 0 respectively. From left to right, the columns are for thermal feedback in the 32 3 & 48 3 simulations, followed by kinetic feedback for 32 3 and 48 3 simulations. Radii of the circles are scaled proportional to M 1/3 gal . Figure 14. Fraction of baryons supplied with supernova energy, at a given time, for thermal and kinetic feedback simulations with varying resolution (solid lines). Also plotted are subsets with δ < 500 (long-dashed lines), δ < 50 (short-dashed lines) and δ < 10 (dotted lines). Redshift intervals are marked along the top of the figure.
sity faster than thermal feedback, even though both cause a similar reduction of star formation by z = 0.

CONCLUSIONS
This paper has surveyed various phenomenological models for incorporating the effects of star formation and associated energy feedback (from supernovae) within cosmological simulations of galaxy formation, using the SPH code hydra). Our main conclusions are as follows.
(i) Most prescriptions for turning cold gas into stars used in previous SPH simulations give similar results when adjusted to give the same present day stellar density, Ω * . The mass in stars formed as a function of time is primarily determined by the gas density at which star formation is assumed to proceed efficiently. Thus, adopting a global density threshold (either in comoving or in physical coordinates), a Schmidt law, or a more complicated prescription involving, for example, the Jeans mass, makes little difference to the global properties of the galaxy population that forms. Given the limited resolution of current simulations (in which the internal properties of galaxies are unresolved), it is simplest to adopt a global density threshold.
(ii) Incorporating star formation within a cosmological simulation makes no significant difference to the rate at which gas cools. The amount of gas cooled by z = 0 only varies up to 1 per cent in simulations with and without star formation and the present day baryonic (gas+stars) mass function of galaxies also remains unchanged. However, when cold, dense gas particles are replaced by collisionless particles, a saving of around 50 per cent of the CPU time is achieved from reducing the (computationally-expensive) work of calculating hydrodynamical forces.
(iii) Simulations without feedback predict cooled gas fractions of around 35 per cent. Not only are these values resolution-dependent, but they exceed observational measurements by around a factor of 5. The problem of "overcooling" is further highlighted in the K-band galaxy luminosity function, which requires the majority of cooled gas to be "hidden" in low-mass non-visible stars (such as browndwarfs), in order to provide a reasonable match to the observed luminosity function.
(iv) Injecting supernova energy to surrounding gas particles, as thermal or kinetic energy, can reduce the star formation rate appreciably. However, thermal models can only achieve this if the reheated gas is prevented from cooling for a considerable fraction of a Hubble time. For maximal thermal feedback (i.e. all available energy is injected into the gas), the results are less sensitive to the choice of delay timescale than for models with a lower efficiency. Average fractions of cold gas in galaxies are unaffected by the inclusion of feedback.
(v) A significant fraction of the gas reheated by supernovae is able to escape galaxy haloes and move into the low density regions (δ < 10) responsible for the Ly-α forest. If this gas is viewed as a tracer of reprocessed material, this mechanism efficiently transports metals into low-density regions.
(vi) The efficiency with which supernovae expel gas from haloes depends on the circular velocity of the halo ( 100 − 300 km s −1 in these simulations) -more gas is expelled from smaller haloes, particularly for models capable of strong feedback. Larger simulations are required to establish robustly the form of this relation.
(vii) Increasing the mass resolution by about a factor of ∼ 3 increases the total amount of resolved galactic material primarily because smaller galaxies are then resolved. Consequently, the K-band galaxy luminosity function is very similar over the range of magnitudes resolved in our high and low resolution simulations. The ratio of cold gas to stars is also essentially unaffected, as is distribution of gas perturbed by supernovae, although increasing the resolution increases the amount of this material, due to the higher number of resolved galaxies. Hence, until numerical convergence is established, the simulations underpredict the amount of this material.
To fully test the effectiveness of the feedback schemes that we have explored in this paper will require larger simulations to be performed in statistically significant volumes. This will allow definite predictions to be made for e.g. the galaxy luminosity function and correlation function, which can readily be compared with observations.

ACKNOWLEDGMENTS
We would like to thank Andrew Benson, George Efstathiou, Cesario Lia, Rob Thacker, Peter Thomas, Tom Theuns, David Weinberg and Simon White for helpful discussions. We also thank the referee for constructive comments that have improved the clarity of this paper. STK acknowledges PPARC for financial support. 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/). This work was supported by the EC TMR network for "Galaxy formation and evolution".