text:
# 1 Introduction

## 1 Introduction

Proteins consist of 20 different amino acids with a great diversity with regard to size, polarity and charge. The understanding of water and interactions with water seems to be important in order to understand protein folding in general, and the special feature of cold unfolding of several small globular proteins in particular.

We in the present work will represent the energy difference between the unfolded and folded interior, with regard to the water, by mimicking additional hydrogen bonds from which we calculate the hydration heat capacity change upon protein unfolding. A justification of the model is the ability for water molecules to form an “ice-like” shell (“iceberg” in the terminology of Frank and Evans) around apolar surfaces and thus create more hydrogen bonds. Reduction of both enthalpy and entropy upon apolar hydration seems to be well established.

However, the protein interior that becomes hydrated upon unfolding, also consists of surfaces that has polarity, which means that the surface has permanent dipoles and charges. The heat capacity upon purely polar hydration becomes surprisingly negative. For apolar surfaces experiments show that the hydration contribution to the heat capacity upon solvation is positive. Also for proteins where part of the surface is polar this heat capacity is positive. Thus for simplicity we in this work will use the apolar “ice-like” shell picture to make a effective model for the hydration effect upon protein unfolding. In this way we may neglect some crucial features of polar solvation.

Finally, we apply equilibrium statistical mechanics to the model and calculate the hydration heat capacity increment, which we compare with experimental data from Privalov and Makhatadze on four different proteins.

## 2 Hydration upon protein unfolding

We will use a refined version of a model first proposed by Hansen et al.. The model studied here was applied by Bakk et al. on a complete protein folding model, but they did not study the hydration effect separately. In this work we will study specifically this hydration upon protein unfolding.

Protein unfolding involves a cavity formation in water with a rearrangement of the water molecules surrounding the unfolded protein. When estimating the solvation energy of exposing the interior of a protein to water, one has to calculate the energy difference between hydrated water, associated to the protein, and bulk water. More precisely, the hydration is defined as the transfer of a solute from a fixed position in the ideal gas phase to a fixed position in the solvent, i.e., water in the present case.

In order to model the effect upon inserting a surface into water, i.e., unfolding of a protein, we use the simplified analogy of a classical electrical dipole in an external electrical field whose energy is

$$E_1=ϵ\mathrm{cos}\vartheta ,$$
(1)
where $`ϵ`$ is a bending distortion constant. The angle $`\vartheta `$ is the polar angle. Eq. (1) is the hydration model used in the works by Bakk et al., and it extends the interpretation of the hydration model applied by Hansen et al. and Bakk et al..

The idea of representing the solvent by dipoles in protein folding was introduced by Warshel and Levitt, and later applications by Russell and Warshel, Fan et al., and Avbelj.

In addition to the energy due to the external field \[Eq. (1)\] we will add a coupling term

$$E_2=\frac{1}{2}\underset{i,j}{}J_{ij}𝐬_i𝐬_j,$$
(2)
thus modeling pair interactions between the water molecules, where $`J_{ij}`$ is the coupling constant between water molecules $`i`$ and $`j`$, and $`𝐬_i`$ is the dipole moment of water molecule $`i`$. For simplicity we put $`|𝐬_i|=1`$. It can be shown that the energy $`E_1+E_2`$ \[Eqs. (1) and (2)\] in a mean field solution can be represented as

$$E(\vartheta )=E_1+E_2=(ϵ+bm)\mathrm{cos}\vartheta +\frac{1}{2}bm^2,$$
(3)
where $`bm`$, with $`b=_jJ_{ij}`$, is the mean field coupling between a water molecule and its surrounding water molecules, has to be determined self consistently. The average dipole moment is defined by $`m=cos\vartheta `$. The standard mean field solution in an effective field $`ϵ_e=ϵ+bm`$ is

$$m=\mathrm{coth}\left(\frac{ϵ_e}{RT}\right)\frac{RT}{ϵ_e}.$$
(4)
Note that here and below $`ϵ`$ and $`b`$ are energies per mole, as the gas constant $`R`$ is used.

For $`N`$ dipoles per protein the partition function for the total hydration contribution upon protein unfolding is

$$\begin{array}{cc}\hfill Z& =\left[_0^{2\pi }𝑑\phi _0^\pi 𝑑\vartheta sin\vartheta \mathrm{exp}\left(\frac{E(\vartheta )}{RT}\right)\right]^N\hfill \\ & =\left[\frac{4\pi RT\mathrm{sinh}\mathbf{(}ϵ_e/(RT)\mathbf{)}}{ϵ_e}\mathrm{exp}\left(\frac{bm^2}{2RT}\right)\right]^N.\hfill \end{array}$$
(5)
The total hydration heat capacity change per mole of proteins is

$$\mathrm{\Delta }C=RT^2\frac{^2}{T^2}\mathrm{ln}Z.$$
(6)
To obtain $`\mathrm{\Delta }C`$, the self consistent Eq. (4) has to be solved numerically with respect to $`m`$. Eq. (4) can be obtained from Eq. (5) via

$$m=\frac{1}{N}\left(\frac{\mathrm{ln}Z}{ϵ_e}\frac{ϵ_e}{ϵ}+\frac{\mathrm{ln}Z}{m}\frac{m}{ϵ}\right).$$
(7)
## 3 Discussion

We want to compare the heat capacity change upon unfolding solvation of the protein interior with experiments. The proteins considered are myoglobin (Mb), lysozyme (Lys), cytochrome c (Cyt), and ribonuclease (Rns) which we compare with experimental data from Privalov and Makhatadze on the hydration contribution to the heat capacity change upon protein unfolding.

The total hydration heat capacity change is shown in Figure 1. The parameter fit to the experimental data agrees quite well with these data.

Figure 1 shows that the heat capacity has a maximum around 25C for Mb, Lys, and Cyt, while this maximum is shifted to around 50C for Rns. Also in Table I a similar relation is reflected.

With some interest we note from Table I that both the “electric field” constant $`ϵ`$ and the coupling constant $`b`$ are essentially the same for all of the four proteins, but there is a small deviation for Rns.

This small deviation for Rns reflects itself in the ratios between accessible polar and total surface area (reported from Makhatadze and Privalov), $`\mathrm{\Delta }A_p/\mathrm{\Delta }A_t`$, which are almost equivalent for Mb, Lys, and Cyt, while this ratio is significantly larger for Rns. Hence, the parameters $`ϵ`$ and $`b`$ may be regarded as effective ones for the combined effect of apolar and polar surfaces as discussed in Section 1.

One notes that Rns differ a bit from the other three proteins considered. This may reflect its larger fraction of polar surface which then also can affect qualitative properties. Our model is more like an effective one for a mixed polar and apolar surface. Thus features specific for polar surfaces are not properly taken into account, but are more or less taken into account by adjusting available parameters. E.g., our present model may seem to give too small curvature on figure 1 for Rns. A reason for this may be the negligence of quantization, which will lower the specific heat and thus increase its curvature for decreasing temperatures. The polar (ionic) forces are relatively strong and the hydrogen atom is light, which both favor quantum effects. Use of a two-level system for this kind of problem can thus reflect such quantization.

## 4 Conclusion

We have proposed a thermodynamical model for the hydration of the protein interior that becomes exposed to water upon unfolding. To our knowledge this is the first model studied and compared to experimental data on the pure hydration heat capacity increment over such broad temperature range (5C to 125 C).

Hydration is modeled in an “ice-like” shell analogy, where the water molecules are represented by interacting dipoles in an external field. Compared with experimental data from Privalov and Makhatadze for the four proteins myoglobin (Mb), lysozyme (Lys), cytochrome c (Cyt), and ribonuclease (Rns) the model fits quite well. The specific values of the the field coupling constant ($`ϵ`$) and dipole coupling constant ($`b`$) \[see Eq. (5)\] are the same for Mb, Lys, and Cyt, while it is slightly different for Rns.

In a future model one should take into account qualitative opposite effects, with respect to the heat capacity, by apolar and polar hydration. Experimentally one finds that the heat capacity change is negative for hydration of purely polar surfaces, in contrast to the apolar surfaces where this heat capacity change is positive.

We have reason to expect that the present model, due to its simplicity, also may be useful in more complete protein folding models.

## Acknowledgments

A. B. thanks the Research Council of Norway for financial support (Contract No. 129619/410).

## Figure captions

Fig. 1. The hydration heat capacity change upon unfolding of four different proteins. Experimental data from Privalov and Makhatadze. Parameters, corresponding to Eq. (5), for the fit to the experimental data are listed in Table I.

A. Bakk, A. Hansen, and J. S. Høye
Specific heat upon aqueous unfolding of …


format: text

text:
# Including star formation and supernova feedback within cosmological simulations of galaxy formation

## 1 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.

These ideas carry through immediately to modern cosmological theory, which is based upon the cold dark matter paradigm (\[White & Frenk 1991\]). This provides a natural motivation for hierarchical clustering and so it inherits the “cooling catastrophe”, generic to this kind of models. There is now a growing body of work attempting to understand likely sources of feedback (e.g. \[Dekel & Silk 1986\]; \[Efstathiou 1992\]; \[Thoul & Weinberg 1996\]; \[Silk & Rees 1998\]; \[Mac Low & Ferrara 1999\]; \[Efstathiou 2000\]; \[Madau, Ferrara & Rees 2000\]), their effects on models of galaxy formation (e.g. \[White & Frenk 1991\]; \[Cole 1991\]; \[Katz 1992\]; \[Navarro & White 1993\]; \[Kauffmann et al. 1994\]; \[Cole et al.1994\]; \[Katz, Weinberg, & Hernquist 1996\]; \[Yepes et al. 1997\]; \[Somerville & Primack 1999\]; \[Thacker & Couchman 2000\]; \[Cole et al. 2000\]), and on the chemical enrichment of the intergalactic medium at high redshift (e.g. \[Ellison et al. 1999\]; \[Schaye et al. 2000\]). Efforts to detect observational signatures of feedback are also being made (e.g. \[Martin 1999\]; \[Theuns, Mo & Schaye 2001\]; \[Pettini et al. 2001\]).

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 (\[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. \[Cole et al.1994\], \[Cole et al. 2000\]; \[Somerville & 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\]). Semi-analytic 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. (\[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\]; \[Frenk et al. 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 statistically-relevant 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. (\[Kay et al 2000\], hereafter K2000)

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.

## 2 Numerical method

### 2.1 The code

The simulations were carried out using the hydra code 111hydra is publicly available and can be downloaded from http://hydra.mcmaster.ca(\[Couchman, Thomas & Pearce 1995\]). This combines an adaptive particle-particle/particle-mesh (AP3M ) 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 (\[Couchman, Thomas & Pearce 1995\]), as discussed by Thacker & Couchman (\[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. (\[Theuns et al. 1998\]).

Radiative cooling is modelled in the manner described by Thomas & Couchman (\[Thomas & Couchman 1992\]), except that tabulated cooling rates (as a function of temperature and metallicity) from Sutherland & Dopita (\[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\text{K}`$. Contributions to the cooling rate at temperatures below $`10^4\text{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-$`\alpha `$ forest, previous studies have shown that the process is unimportant for galaxies forming in haloes with circular velocities above $`50\text{km}\text{s}^1`$ (e.g. \[Quinn, Katz & Efstathiou 1996\]; \[Navarro & Steinmetz 1997\]; \[Weinberg, Hernquist & Katz 1997\]), as is the case in our simulations.

### 2.2 Cosmological models

We carried out simulations using two cold dark matter models: SCDM and $`\mathrm{\Lambda }`$CDM. The former assumes a critical density universe, $`\mathrm{\Omega }_0=1`$ and $`\mathrm{\Lambda }_0=0`$, Hubble constant, in units of $`100\text{km}\text{s}^1\text{Mpc}^1`$, $`h=0.5`$ and baryon density, $`\mathrm{\Omega }_b=0.06`$. The $`\mathrm{\Lambda }`$CDM model assumes a flat, low density universe, with $`\mathrm{\Omega }_0=0.35`$ and $`\mathrm{\Lambda }_0=0.65`$, Hubble constant, $`h=0.71`$, and baryon density, $`\mathrm{\Omega }_b=0.038`$. Current determinations of cosmological parameters (from, for example, microwave background anisotropies, high redshift supernovae, and the baryon fraction in galaxy clusters) favour the $`\mathrm{\Lambda }`$CDM model. However, we do not expect our conclusions to differ significantly in other hierarchical models.

Initial density fluctuations were calculated for the SCDM model using the approximate transfer function given by Bond & Efstathiou (\[Bond & Efstathiou 1984\]) with shape parameter, $`\mathrm{\Gamma }=\mathrm{\Omega }_0h=0.5`$ and amplitude set by $`\sigma _8=0.6`$ (as required for agreement with the present day abundance of rich clusters; e.g. \[Vianna & Liddle 1996\]; \[Eke, Cole & Frenk 1996\]). For the $`\mathrm{\Lambda }`$CDM model, the transfer function of Bardeen et al. (\[Bardeen et al. 1986\]) was adopted and $`\mathrm{\Gamma }`$ was set to $`\mathrm{\Gamma }=\mathrm{\Omega }_0h\mathrm{exp}(\mathrm{\Omega }_b/\mathrm{\Omega }_0(\mathrm{\Omega }_0+\sqrt{2h}))=0.21`$ (\[Sugiyama 1995\]). The amplitude of density fluctuations was normalised by setting $`\sigma _8=0.9`$, as required by the cluster abundance argument in this case.

The gas particle mass was held fixed for all simulations within each cosmological model, at $`5\times 10^8h^1\text{M}_{}`$ in the case of SCDM and to $`3\times 10^8h^1\text{M}_{}`$ in the case of $`\mathrm{\Lambda }`$CDM. Both models assumed an initial (equivalent Plummer) softening length of $`20h^1\text{kpc}`$ which was fixed in comoving coordinates until $`z1`$, after which it was fixed in physical coordinates to $`10h^1\text{kpc}`$ until $`z=0`$. The SCDM simulations were started at a redshift of 24 and the $`\mathrm{\Lambda }`$CDM simulations at a redshift of 49 and both sets were run to $`z=0`$.

## 3 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 $`10h^1\text{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.

### 3.1 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.

#### 3.1.1 The Pearce method

This is based on the method used by Pearce (\[Pearce 1998\]). Potential star particles are identified with gas particles having overdensity, $`\delta >550`$ and temperature, $`T<10^5\text{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 $`2h_{\mathrm{min}}`$ containing $`N_{\mathrm{SPH}}`$ neighbours. (The minimum SPH smoothing length, $`h_{\mathrm{min}}`$ is set equal to 1.17 times the equivalent Plummer softening length.) For our choice of parameters ($`2h_{\mathrm{min}}=0.0234h^1\text{Mpc}`$ and $`N_{\mathrm{SPH}}=32`$), this limit is $`\delta 20,000`$ at $`z=0`$. We deliberately set our choice of overdensity thresholds to straddle this resolution limit.

#### 3.1.2 The Summers method

In the method used by Summers (\[Summers 1993\]) gas particles are converted into stars if they satisfy a temperature criterion, $`T<3\times 10^4\text{K}`$ and a physical density criterion, $`n_\mathrm{H}>0.1\text{cm}^3`$. Additionally, eligible gas particles must be within overdense regions, $`\delta >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, $`\rho >10^{25}\text{g}\text{cm}^3`$ ($`\rho =m_\mathrm{H}n_\mathrm{H}/X`$; we assume $`X=0.76`$), approximately a factor of 2 lower than used by Summers. This threshold corresponds to a baryon overdensity, $`\delta 3.5\times 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 $`\rho >10^{26}\text{g}\text{cm}^3`$.

#### 3.1.3 The Navarro & White method

We have also studied the method of Navarro & White (\[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 ($`.𝐯<0`$) and has a density, $`\rho >7\times 10^{26}\text{g}\text{cm}^3`$. Gas particles satisfying these criteria are split into equal mass gas and star particles, after a time interval, $`\mathrm{\Delta }t_{}=\mathrm{ln2}\tau _{\mathrm{dyn}}`$, where $`\tau _{\mathrm{dyn}}1/\sqrt{\rho }`$ 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.

#### 3.1.4 The Katz, Weinberg & Hernquist method

The fourth method investigated is that used in the cosmological simulations of Katz, Weinberg & Hernquist (\[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, $`\delta >55.7`$, and a Jeans instability condition, $`\tau _\mathrm{J}>\tau _{\mathrm{dyn}}`$, where $`\tau _\mathrm{J}h/c_s`$ is a “Jeans” timescale, $`h`$ is the particle’s SPH smoothing length and $`c_s=(\gamma (\gamma 1)u)^{1/2}`$ 222$`u=kT/(\gamma 1)\mu m_\mathrm{H}`$; we take the ratio of specific heats at constant volume and constant pressure for a monatomic ideal gas to be $`\gamma =5/3`$, and the mean molecular weight of a gas of primordial composition to be $`\mu =0.59`$ 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\mathrm{exp}(\mathrm{\Delta }t/\mathrm{min}(\tau _{\mathrm{dyn}},\tau _{\mathrm{cool}}))`$, where $`r`$ is a random number between 0 and 1, $`\mathrm{\Delta }t`$ is the timestep and $`\tau _{\mathrm{cool}}1/\rho `$ 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.

#### 3.1.5 The Mihos & Hernquist method

The final model studied is based on the method used by Mihos & Hernquist (\[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 $`\mathrm{\Delta }M_{}/M_{\mathrm{gas}}=C_{\mathrm{sfr}}\rho ^{N1}\mathrm{\Delta }t`$. This equation is applied to gas particles with $`\delta >10`$ and $`.𝐯<0`$, although rather than converting part of a particle’s mass into stars, we convert all of it if $`\mathrm{\Delta }M_{}/M_{\mathrm{gas}}>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_{\mathrm{sfr}}=5\times 10^4\mathrm{cm}^{3/2}\mathrm{g}^{1/2}\mathrm{s}^1`$.

### 3.2 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 comparison, 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 $`\rho >500\rho `$ and $`T<30,000\text{K}`$, labelled $`f_{\mathrm{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 $`\rho T`$ plane containing material which has radiatively cooled (see, for example, §3.2 of K2000). For all simulations, including gas, $`f_\mathrm{g}=0.34\pm 0.01`$. As found by Pearce (\[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_{\mathrm{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_{\mathrm{sfr}}`$. The value, $`C_{\mathrm{sfr}}=5\times 10^4\mathrm{cm}^{3/2}\mathrm{g}^{1/2}\mathrm{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-of-friends group-finder (\[Davis et al. 1985\]) on the cold gas and star particles (i.e. on those used to calculate $`f_{\mathrm{cool}}`$ above), with a linking length, $`b=0.1`$ (comparable to the gravitational softening length.) Any objects with fewer star particles than $`N_{\mathrm{SPH}}`$ (32 for the simulations studied in this paper) were discarded. All runs contain $`70`$ galaxies at $`z=0`$ to within 10%.

#### 3.2.2 Star formation histories

Fig. 1 illustrates the star formation histories in the simulations by means of the cumulative fraction of baryons turned into stars as function of time, relative to the $`z=0`$ fraction (the corresponding redshift intervals are marked along the top of the figure). The cumulative fraction of cooled gas ($`\rho >500\rho `$ and $`T<30,000\text{K}`$) for the simulation without star formation, is also plotted. The instantaneous star formation rate is proportional to the gradient of each star fraction curve.

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 $`z2`$, 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 ($`z1`$), 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=12`$, 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 $`z1.7`$, this density has dropped to $`\rho 10^{26}\text{g}\text{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, $`\tau _{\mathrm{sfr}}=M_{\mathrm{gas}}\mathrm{\Delta }t/\mathrm{\Delta }M_{}=1/C_{\mathrm{sfr}}\rho _{\mathrm{gas}}^{1/2}`$ (c.f. §3.1.5) to the lookback time for an $`\mathrm{\Omega }=1`$ universe, $`\tau _{\mathrm{look}}=(2/3H_0)(1(1+z)^{3/2})`$; at higher densities, $`\tau _{\mathrm{sfr}}<\tau _{\mathrm{look}}`$. Assuming that $`\delta \rho /\rho `$ (i.e. for $`\delta >>1`$), the minimum overdensity scales as $`\delta _{\mathrm{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.

#### 3.2.3 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, $`\mathrm{\Delta }\mathrm{log}M_{}=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\text{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\text{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\times 10^{10}h^1\text{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. (\[Cole et al. 2000\]), $`L^{}`$ disk galaxies ($`M/L_B2`$ and $`M_B^{}19.5`$) have $`M_{\mathrm{cg}}/M_{}0.10.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\text{M}_{}`$ implies a ratio of at least 2).

#### 3.2.4 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(t_0t_f)`$. The total magnitude is then

$$M_K=5\mathrm{l}\mathrm{o}\mathrm{g}h+M_K(0)2.5\mathrm{log}\underset{i=1}{\overset{N_{}}{}}l(t_0t_f)m_{},$$
(1)
where $`M_K(0)`$ is the zero–point and $`N_{}`$ is the number of star particles in the galaxy. The function, $`l(t_0t_f)`$ was generated using the stellar population synthesis models of Bruzual & Charlot (\[Bruzual & Charlot 1993\]), assuming solar metallicity and the stellar IMF proposed by Kennicutt (\[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.5\mathrm{log}\mathrm{{\rm Y}}`$, where $`\mathrm{{\rm Y}}1+M(\mathrm{brown}\mathrm{dwarfs})/M(\mathrm{visible}\mathrm{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 (\[Cole et al. 2001\]), with the value of $`M_K^{}`$ marked with a vertical arrow. The model luminosity functions are plotted assuming $`\mathrm{{\rm Y}}=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 $`\mathrm{{\rm Y}}=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\]).

### 3.3 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.

## 4 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(\[Katz 1992\]), who studied the collapse of a galaxy-sized 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 (\[Navarro & White 1993\]). In this model, a fraction, $`f_v`$, 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 $`f_v`$ (i.e. almost all thermal energy), very little feedback occurs and consequently the star formation rate is very high. However, for higher values of $`f_v`$ 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 (\[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 (\[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 (\[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 1992\], \[Cen & Ostriker 1999\]; \[Yepes 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 (\[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 (\[Summers 1993\]) and Katz et al. (\[Katz, Weinberg, & Hernquist 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.

### 4.1 Simulation details

#### 4.1.1 Cooling and Star formation

We have studied the effects of feedback in the $`\mathrm{\Lambda }`$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 $`\mathrm{\Lambda }`$CDM. We implement the “decoupling” technique described by Pearce et al. (\[Pearce et al. 2000\]) whereby gas particles with temperature below $`T=12,000\text{K}`$ are ignored when estimating the gas density of hot particles with temperature above $`T=10^5\text{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,000\text{K}`$ (essentially the pearce1 criteria). When a gas particle satisfies both criteria it is instantaneously converted into a star particle. In the $`\mathrm{\Lambda }`$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.

#### 4.1.2 Supernova energetics

We adopt the same parameter values for the supernovae energetics as Katz (1992) and Katz et al. (\[Katz, Weinberg, & Hernquist 1996\]). To summarise, the amount of energy available from supernovae is calculated assuming that stars with mass above $`8\text{M}_{}`$ release $`10^{51}\text{erg}`$ of energy into the surrounding gas. A Miller & Scalo (\[Miller & Scalo 1979\]) IMF is assumed, with lower and upper mass limits of $`0.1\text{M}_{}`$ and $`100\text{M}_{}`$ respectively, implying a specific energy released in supernovae of $`3.7\times 10^{15}\text{erg}\mathrm{g}^1`$, equivalent to a temperature of $`1.8\times 10^7\text{K}`$. The size of the timestep in the simulations (typically $`3\mathrm{Myr}`$) is about one order of magnitude smaller than the lifetime of an $`8\text{M}_{}`$ star ($`\tau _820\mathrm{M}\mathrm{y}\mathrm{r}`$). For this reason, the energy was released gradually, using the following function

$$\dot{u}_{\mathrm{SN}}(t)=3.7\times 10^{15}\mathrm{exp}((tt_0)/\tau _8)\tau _8^1\text{erg}\mathrm{g}^1,$$
(2)
where $`t_0`$ is the time when the star formation event occurred. Energy was released until $`tt_0=200\text{Myr}`$, when $`\dot{u}_{\mathrm{SN}}`$ becomes negligible.

#### 4.1.3 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, $`h_{\mathrm{min}}`$, and each gas neighbour within the SPH smoothing radius (2$`h_{\mathrm{min}}`$) 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 $`\mathrm{\Delta }t`$, is

$$\mathrm{\Delta }u_{\mathrm{SN},j}=ϵ\dot{u}_{\mathrm{SN},i}\mathrm{\Delta }t\frac{H(2h_{\mathrm{min}}|𝐫_i𝐫_j|)}{_{j=1}^{N_{\mathrm{gas}}}H(2h_{\mathrm{min}}|𝐫_i𝐫_j|)},$$
(3)
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 re-radiate the energy in a short time. To circumvent this problem, we follow the method of Gerritsen (\[Gerritsen 1997\]) and prevent the gas from cooling over a period of time. This timescale is controlled using the parameter $`\tau `$, 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

$$\mathrm{\Delta }v_{\mathrm{SN},\mathrm{j}}=\sqrt{2ϵ\dot{u}_{\mathrm{SN},i}\mathrm{\Delta }t\frac{H(2h_{\mathrm{min}}|𝐫_i𝐫_j|)}{_{j=1}^{N_{\mathrm{gas}}}H(2h_{\mathrm{min}}|𝐫_i𝐫_j|)}}.$$
(4)
This scheme is the “momentum” implementation of \[Navarro & White 1993\].

### 4.2 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 $`\tau =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 $`\mathrm{\Lambda }`$CDM and SCDM simulations are summarised in Table 4.

#### 4.2.1 General properties at $`z=0`$

We identify galaxies by applying a friends-of-friends group-finder 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\times 10^{10}h^1\text{M}_{}`$ for the $`\mathrm{\Lambda }`$CDM model and $`1.6\times 10^{10}h^1\text{M}_{}`$ for the SCDM model.

Column 5 of Table 4 gives the fraction of baryons (cold gas and stars) in galaxies, labelled $`f_{\mathrm{gal}}`$. For simulations with the same feedback parameters, the fraction of galaxy material is always higher in the $`\mathrm{\Lambda }`$CDM case than in the SCDM case though the trends with varying $`ϵ`$ and $`\tau `$ 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 $`\tau `$ results in a smaller $`f_{\mathrm{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 $`\tau `$ considered, $`\tau =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 $`\mathrm{\Lambda }`$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.

#### 4.2.2 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 $`\mathrm{\Lambda }`$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 $`\tau =0.01`$, short-dashed with $`\tau =0.1`$ and long-dashed with $`\tau =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 $`\tau `$ (when varied over 2 orders of magnitude). However, increasing $`ϵ`$ by a factor of 10 produces only small differences in the star formation rate when $`\tau `$ is varied. In the latter case, enough energy is deposited into the gas that even on the shortest timescale considered, $`\tau =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\text{K}`$ is able to cool radiatively at all times. Although supplying 0.1 per cent of the available 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.

#### 4.2.3 The fate of reheated gas

Virtually all the gas that is reheated, either thermally or kinetically, has a temperature in the range $`10^410^5\text{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\text{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 (Katz et al. \[Katz, Weinberg, & Hernquist 1996\]).

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 $`\mathrm{\Lambda }`$CDM simulation with thermal feedback, for $`ϵ=1`$ and $`\tau =1`$. The density ($`n/\text{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\rho /\mu m_\mathrm{H}`$.

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\text{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 shock-heating) 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 $`V_c60\text{km}\text{s}^1`$, contains approximately 25 particles and so is also poorly resolved at this stage.) The density reaches a peak value of $`0.1\text{cm}^3`$ at $`z3`$. 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\text{K}`$. Correspondingly, the density decreases by 3 orders of magnitude (aided by the fact that since the particle is at $`T>10^5\text{K}`$, its SPH density is determined only by other hot particles) and its entropy jumps to $`100\mathrm{k}\mathrm{e}\mathrm{V}\mathrm{c}\mathrm{m}^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 $`\mathrm{\Lambda }`$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 $`\delta <500`$, the short-dashed lines the fraction with $`\delta <50`$ and the dotted lines the fraction with $`\delta <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-$`\alpha `$ 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 $`\tau =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,\tau =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}\text{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 ($`\delta <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.

#### 4.2.4 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 $`\mathrm{\Delta }_c\rho _{\mathrm{cr}}(z)`$, where $`\rho _{\mathrm{cr}}(z)=3H^2(z)/8\pi G`$; we set $`\mathrm{\Delta }_c=178\mathrm{\Omega }^{0.45}(z)`$ for the $`\mathrm{\Lambda }`$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

$$V_c=740M_{14}^{1/3}\left(\frac{H(z)}{H_0}\right)^{1/3}\left(\frac{\mathrm{\Delta }_c}{178}\right)^{1/6}\text{km}\text{s}^1,$$
(5)
where $`M_{14}`$ is the total mass of the halo in $`10^{14}h^1\text{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_{\mathrm{halo}}\mathrm{\Omega }_\mathrm{b}/\mathrm{\Omega }_0`$, where $`M_{\mathrm{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 $`\mathrm{\Lambda }`$CDM thermal feedback simulations with $`\tau =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 ($`V_c220\text{km}\text{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 $`V_c`$. 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, $`\tau `$.

#### 4.2.5 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, $`\mathrm{M}_{\mathrm{cg}}/\mathrm{M}_{}`$, against $`M_{}`$ for galaxies (with non-zero gas content) in the $`\mathrm{\Lambda }`$CDM simulations with 32 or more particles. The points in Fig. 9 represent the median values within bins of width, $`\mathrm{\Delta }\mathrm{log}M_{}=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\text{M}_{}`$. This agrees well with the observational data summarised by Cole et al. (\[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 $`\mathrm{\Lambda }`$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,\tau =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 $`\tau `$), the stellar mass and hence the luminosity of the galaxies are reduced. Both the thermal model with $`ϵ=1`$ & $`\tau =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.

### 4.3 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 $`\tau `$). 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 $`\tau `$ 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 $`V_c100\text{km}\text{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.

## 5 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\text{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 $`\tau =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\times 10^{10}h^1\text{M}_{}`$ to $`4.8\times 10^9h^1\text{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.7h^1\text{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_{\mathrm{gal}}`$, and the total number of galaxies, $`N_{\mathrm{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 $`\mathrm{\Delta }\mathrm{log}M_{}=0.5`$ for low resolution simulations and 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 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\times 64`$ grid, using the SPH kernel with width set by the particles’ smoothing lengths. Overlaid as circles (with radii proportional to $`M_{\mathrm{gal}}^{1/3}`$) 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}\mathrm{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 density faster than thermal feedback, even though both cause a similar reduction of star formation by $`z=0`$.

## 6 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, $`\mathrm{\Omega }_{}`$. 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 brown-dwarfs), 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 ($`\delta <10`$) responsible for the Ly-$`\alpha `$ 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 ( $`100300\text{km}\text{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”.


text:
# 1 Introduction

## 1 Introduction

An essential property of gravitational dynamics is that the canonical Hamiltonian vanishes on the physical trajectories of the system, which then includes a constraint $`0`$. This reflects the fact that the evolution of the gravitational field is given in terms of a parameter $`\tau `$ which does not have physical significance. This feature leads to a fundamental difference between ordinary quantum mechanics and the quantization of gravitation, because the existence of a unitary quantum theory is related to the possibility of defining the time as an absolute parameter. The identification of a global phase time  can therefore be considered as the previous step before quantization .

For minisuperspace models we have an action functional of the form

$$S[q^i,p_i,N]=_{\tau _1}^{\tau _2}\left(p_i\frac{dq^i}{d\tau }N\right)𝑑\tau $$
(1)
where $`N`$ is a Lagrange multiplier enforcing the quadratic Hamiltonian constraint

$$=G^{ij}p_ip_j+V(q)0,$$
(2)
with $`G^{ij}`$ the reduced version of the DeWitt supermetric. The extremal condition $`\delta S=0`$ gives the canonical equations

$$\frac{dq^i}{d\tau }=N[q^i,],\frac{dp_i}{d\tau }=N[p_i,].$$
(3)
The solution of these equations describes the evolution of a spacelike hypersurface along the timelike direction; the presence the of the multiplier $`N`$ introduces an arbitrariness in the evolution which is associated to a multiplicity of times. From a different point of view, the constraint $`0`$ acts as a generator of gauge transformations which can be written

$$\delta _ϵq^i=ϵ(\tau )[q^i,],\delta _ϵp_i=ϵ(\tau )[p_i,],\delta _ϵN=\frac{ϵ(\tau )}{\tau }.$$
(4)
From (3) and (4) we see that the dynamical evolution can be reproduced by a gauge transformation progressing with time, that is, any two succesive points on each classical trajectory are connected by a gauge transformation . Hence, the gauge fixation can be thought not only as a way to select one path from each class of equivalent paths in phase space, but also as a reduction procedure identifying a time for the system.

Admissible gauge conditions are those which can be reached from any path in phase space by means of gauge transformations leaving the action unchanged. Under a gauge transformation defined by the parameters $`ϵ^m`$ the action of a system with constaints $`C_m`$ changes by

$$\delta _ϵS=\left[ϵ^m(\tau )\left(p_i\frac{C_m}{p_i}C_m\right)\right]_{\tau _1}^{\tau _2}.$$
(5)
For ordinary gauge systems, which include constraints that are linear and homogeneous in the momenta, it is $`\delta _ϵS=0`$, and gauge conditions of the form $`\chi (q,p,\tau )=0`$ (canonical gauges) are admissible. In the case of gravitation, instead, the Hamiltonian constraint is quadratic in the momenta, and we would have $`\delta _ϵS0`$ unless $`ϵ(\tau _1)=ϵ(\tau _2)=0`$ ; then gauge conditions involving derivatives of Lagrange multipliers as, for example, $`\chi dN/d\tau =0`$ (derivative gauges) should be used  . These gauges cannot define a time in terms of the canonical variables. At the quantum level this has the consequence that the usual Fadeev–Popov path integral procedure for quantizing gauge systems could not be applied .

However, in recent papers - we have shown that if the system results separable the action can be provided with gauge invariance at the end points by means of a canonical transformation $`(q^i,p_i)(Q^i,P_i)`$ which matches the constraint $``$ with one of the new momenta. Then canonical gauges are admissible and a global phase time can be identified by imposing $`\tau `$dependent gauge conditions, and, simultaneously, the transition amplitude can be obtained by means of the usual path integral procedure for gauge systems. Moreover, we have shown that a transition amplitude between two given quantum configurations in terms of only the original coordinates can be found only in the case that an intrinsic time exists.

A function $`t(q^i,p_i)`$ is a global phase time if $`[t,]>0.`$ Because the supermetric $`G^{ik}`$ does not depend on the momenta, a function $`t(q)`$ is a global phase time if the bracket

$`[t(q),]`$ $`=`$ $`[t(q),G^{ik}p_ip_k]`$
$`=`$ $`2{\displaystyle \frac{t}{q^i}}G^{ik}p_k`$
is positive definite. If the supermetric has a diagonal form and one of the momenta vanishes at a given point of phase space, then no function of its conjugate coordinate only can be a global phase time. This is something to be remarked, as in some early works  the scale factor was used as time parameter in the obtention of a Schrödinger equation for the minisuperspaces, and this is not licit in general; examples of this are the cases of Kantowski–Sachs anisotropic universe  and isotropic models yielding from the low energy limit of closed bosonic string theory . For a constraint whose potential can be zero for finite values of the coordinates, the momenta $`p_k`$ can be all equal to zero at a given point, and $`[t(q),]`$ can vanish. Hence an intrinsic global phase time $`t(q)`$ can be identified only if the potential in the constraint has a definite sign. In the most general case a global phase time should be a function including the canonical momenta; in this case it is said that the system has an extrinsic time $`t(q,p)`$ , because the momenta are related to the extrinsic curvature. When only an extrinsic time exists we shall have to revise some points of the path integral quantization to which we are used: as we shall see below, a quantum description in terms of only the original coordinates may be impossible if we want to work in a theory with a clear notion of time.

## 2 A simple model: The closed de Sitter universe.

In order to discuss this topic we shall use a simple model which presents all the peculiarities of the models of physical interest; this model has the special feature that the transition probability is trivially known since it has one constraint and only one degree of freedom, so that there is only one physical state for the system. Despite its simplicity, it is a good example to understand the quantization with extrinsic time.

Consider the Hamiltonian constraint of the most general empty homogeneous and isotropic cosmological model:

$$=\frac{1}{4}e^{3\mathrm{\Omega }}\pi _\mathrm{\Omega }^2ke^\mathrm{\Omega }+\mathrm{\Lambda }e^{3\mathrm{\Omega }}0.$$
(6)
This Hamiltonian corresponds to a universe with arbitrary curvature $`k=1,0,1`$ and non zero cosmological constant; we shall suppose $`\mathrm{\Lambda }>0`$. In the case $`k=0`$ we obtain the de Sitter universe; although the absence of matter makes this universe basically a toy model, it has received considerable attention because it reproduces the behaviour of models with matter or with non zero curvature when the scale factor $`ae^\mathrm{\Omega }`$ is great enough. The classical evolution is easy to obtain, and corresponds to an exponential expansion. In fact, for both $`k=0`$ and $`k=1`$ the potential is never zero, and then $`\pi _\mathrm{\Omega }`$ cannot change its sign. Instead, for the closed model $`\pi _\mathrm{\Omega }=0`$ is possible.

It is convenient to work with the rescaled Hamiltonian $`H=e^\mathrm{\Omega }:`$

$$H=\frac{1}{4}e^{4\mathrm{\Omega }}\pi _\mathrm{\Omega }^2k+\mathrm{\Lambda }e^{2\mathrm{\Omega }}0.$$
(7)
The constraints $`H`$ and $``$ are equivalent because they differ only in a positive definite factor.

We shall turn the parametrized system asociated to the Hamiltonian (7) into an ordinary gauge system by means of two succesive canonical transformations. The $`\tau `$independent Hamilton–Jacobi equation for the Hamiltonian $`H`$ is

$$\left(\frac{W}{\mathrm{\Omega }}\right)^24ke^{4\mathrm{\Omega }}+4\mathrm{\Lambda }e^{6\mathrm{\Omega }}=4e^{4\mathrm{\Omega }}E.$$
(8)
Matching $`E=\overline{P}_0`$ we obtain the solution

$$W(\mathrm{\Omega },\overline{P}_0)=2sign(\pi _\mathrm{\Omega })𝑑\mathrm{\Omega }e^{2\mathrm{\Omega }}\sqrt{\mathrm{\Lambda }e^{2\mathrm{\Omega }}k\overline{P}_0},$$
(9)
which is the generating function of the canonical transformation $`(\mathrm{\Omega },\pi _\mathrm{\Omega })(\overline{Q}^0,\overline{P}_0)`$ defined by

$$\overline{Q}^0=\frac{W}{\overline{P}_0}=sign(\pi _\mathrm{\Omega })\mathrm{\Lambda }^1\sqrt{\mathrm{\Lambda }e^{2\mathrm{\Omega }}k\overline{P}_0},\overline{P}_0=H.$$
(10)
Then we define the function $`F=\overline{Q}^0P_0+f(\tau )`$ which generates the second canonical transformation yielding a non vanishing true Hamiltonian $`h=f/\tau `$ and $`Q^0=\overline{Q}^0`$, $`\overline{P}_0=P_0.`$ This second transformation may seem to be an unnecesary sofistication, but it plays a central role in the case of cosmological models with true degrees of freedom, as it allows the new coordinates $`Q`$ associated to observables to be fixed at arbitrary values at the boundaries. If possible, $`f`$ should be chosen in such a way that the reduced Hamiltonian is conservative .

The variables $`Q^0`$ and $`P_0`$ describe the gauge system into which the model has been turned. Therefore the gauge can be fixed by means of a $`\tau `$dependent canonical condition like $`\chi Q^0T(\tau )=0`$ with $`T`$ a monotonic function of $`\tau `$. Then as $`[Q^0,P_0]=1`$ we can define the global phase time as

$$t=Q^0|_{P_0=0}=sign(\pi _\mathrm{\Omega })\mathrm{\Lambda }^1\sqrt{\mathrm{\Lambda }e^{2\mathrm{\Omega }}k}$$
(11)
As on the constraint surface $`P_0=0`$ we have

$$\pi _\mathrm{\Omega }=2sign(\pi _\mathrm{\Omega })e^{2\mathrm{\Omega }}\sqrt{\mathrm{\Lambda }e^{2\mathrm{\Omega }}k},$$
(12)
(so that in the case $`k=1`$ the natural size of the configuration space is given by $`\mathrm{\Omega }\mathrm{ln}(\sqrt{\mathrm{\Lambda }})`$ ) we can write

$$t(\mathrm{\Omega },\pi _\mathrm{\Omega })=\frac{1}{2}\mathrm{\Lambda }^1e^{2\mathrm{\Omega }}\pi _\mathrm{\Omega },$$
(13)
which is in agreement with the time obtained by matching the model with the ideal clock . It is interesting to notice that this expresion for the time can be also obtained by demanding that its functional form presents no ambiguities under different ways of factorizing the Hamiltonian constraint.

Now an important difference between the cases $`k=1`$ and $`k=1`$ arises: for $`k=1`$ the potential has a definite sign, and the constraint surface splits into two disjoint sheets given by (12). In this case the evolution can be parametrized by a function of the coordinate $`\mathrm{\Omega }`$ only, the choice given by the sheet on which the system remains: if the system is on the sheet $`\pi _\mathrm{\Omega }>0`$ the time is $`t=\mathrm{\Lambda }^1\sqrt{\mathrm{\Lambda }e^{2\mathrm{\Omega }}+1}`$, and if it is on the sheet $`\pi _\mathrm{\Omega }<0`$ we have $`t=\mathrm{\Lambda }^1\sqrt{\mathrm{\Lambda }e^{2\mathrm{\Omega }}+1}`$. The deparametrization of the flat model is completely analogous. For the closed model, instead, the potential can be zero and the topology of the constraint surface is no more equivalent to that of two disjoint planes. Although for $`\mathrm{\Omega }=\mathrm{ln}(\sqrt{\mathrm{\Lambda }})`$ we have $`V(\mathrm{\Omega })=0`$ and $`\pi _\mathrm{\Omega }=0,`$ it is easy to verify that $`d\pi _\mathrm{\Omega }/d\tau 0`$ at this point. Hence, in this case the coordinate $`\mathrm{\Omega }`$ does not suffice to parametrize the evolution, because the system can go from $`(\mathrm{\Omega },\pi _\mathrm{\Omega })`$ to $`(\mathrm{\Omega },\pi _\mathrm{\Omega })`$; therefore we must necessarily define a global phase time as a function of both the coordinate and the momentum.

The system has one degree of freedom and one constraint, so that it is pure gauge. In other words, there is only one physical state, in the sense that from a given point in the phase space we can reach any other point on the constraint surface by means of a finite gauge transformation. For this reason, if the deparametrization procedure is consistent it should be possible to verify that the transition probability written in terms of the variables which include a globally well defined time is equal to unity

The quantization proceeds as in references -, and the observation above is reflected in the fact that we obtain

$`<Q_2^0,\tau _2|Q_1^0,\tau _1>`$ $`=`$ $`{\displaystyle DQ^0DP_0DN\delta (Q^0\tau )\mathrm{exp}\left(i_{\tau _1}^{\tau _2}\left[P_0\frac{dQ^0}{d\tau }NP_0\frac{f}{\tau }\right]𝑑\tau \right)}`$ (14)
$`=`$ $`{\displaystyle DQ^0DP_0\delta (P_0)\delta (Q^0\tau )\mathrm{exp}\left(i_{\tau _1}^{\tau _2}\left[P_0\frac{dQ^0}{d\tau }\frac{f}{\tau }\right]𝑑\tau \right)}`$
$`=`$ $`\mathrm{exp}\left(i{\displaystyle _{\tau _1}^{\tau _2}}{\displaystyle \frac{f}{\tau }}𝑑\tau \right),`$
and then the probability for the transition from $`Q_1^0`$ at $`\tau _1`$ to $`Q_2^0`$ at $`\tau _2`$ is

$$\left|<Q_2^0,\tau _2|Q_1^0,\tau _1>\right|^2=1.$$
(15)
When the model is open or flat the coordinates $`\mathrm{\Omega }`$ and $`Q^0`$ are uniquely related, and the result can be easily understood in the sense that once a gauge is fixed there is only one possible value of the scale factor $`ae^\mathrm{\Omega }`$ at each $`\tau `$. But in the case of the closed model we have seen that this is not true: at each $`\tau `$ there are two possible values of $`\mathrm{\Omega }`$; instead, there is only one possible value of $`\pi _\mathrm{\Omega }`$ at each $`\tau `$. Hence the transition probability in terms of $`Q^0`$ does not correspond to the evolution of the coordinate $`\mathrm{\Omega }`$, but rather of its derivative; then we conclude that the amplitude $`<Q_2^0,\tau _2|Q_1^0,\tau _1>`$ corresponds to an amplitude $`<\pi _{\mathrm{\Omega },2}|\pi _{\mathrm{\Omega },1}>`$, and we have

$$|<\pi _{\mathrm{\Omega },2}|\pi _{\mathrm{\Omega },1}>|^2=1.$$
(16)
The fact that $`<Q_2^0,\tau _2|Q_1^0,\tau _1>`$ is not equivalent to $`<\mathrm{\Omega }_2|\mathrm{\Omega }_1>`$ is natural, as the nonexistence of an intrinsic time makes impossible to find a globally good gauge condition giving $`\tau `$ as a function of $`\mathrm{\Omega }`$ only . But precisely for this reason, this should not be taken as a failure of the quantization procedure, because a characterization of the states in terms of only the original coordinates is not correct if we want to retain a clear notion of time on the whole evolution.

## 3 Discussion

It is common to regard an intrinsic time as more natural, and the necessity of defining an extrinsic time as a problematic peculiarity. However, this is perhaps only a consequence of working with simple parametrized systems like, for example, the relativistic particle; the formalism for these systems, when put in a manifestly covariant form, has the time included among the coordinates, and the evolution is given in terms of a physically meaningless time parameter. But while for these systems the time coordinate always refers to an external clock, this is clearly not the case in cosmology; for example, in the case of pure gravitational dynamics the coordinates are the elements of the metric $`g_{ab}`$ over spatial slices, and in principle there is not necessarily a connection between $`g_{ab}`$ and anything external. Rather, such a relation can be thought to exist for the derivatives $`dg_{ab}/d\tau `$ of the metric, as they appear in the expression for the extrinsic curvature $`K_{ab}`$ which describes the evolution of spacelike 3-dimensional hypersurfaces in 4-dimensional spacetime. If no matter fields are present the canonical momenta are given by

$$\pi _ip^{ab}=2G^{abcd}K_{cd},$$
(17)
and then one must expect the momenta to appear in the definition of a global phase time. The existence of a time in terms of only the coordinates should therefore be understood as a sort of an accident related to the fact that, in some special cases which do not represent the general features of gravitation, there exists a relation that enables to obtain the coordinates in terms of the momenta with no ambiguities.

This means that we must revise some points of the path integral quantization to which we are used. Indeed, the system considered in this brief note represents an example of a cosmological model for which the characterization of the quantum states must include the momenta.

It is worth noting that at the quantum level the definition of an extrinsic time in terms of the functional form (13) is not sufficient; in fact, it is also necessary to propose a prescription for the operatorial order between coordinates and momenta to give a precise definition of time. Indeed, it is possible to verify that an ordering which leads to define an extrinsic time for the closed de Sitter universe is given by

$$\widehat{t}\widehat{\pi }_\mathrm{\Omega }e^{2\widehat{\mathrm{\Omega }}}.$$
(18)
We must remark that different orderings generate in the conmutator $`[\widehat{t},\widehat{H}]`$ linear terms in $`\pi _\mathrm{\Omega }`$, which lead to a non definite sign.

At this point, it could be interesting to comment that the description of the evolution of cosmological models in terms of the momenta is not a particular property of quantization; in fact, this is a common feature of classical cosmology where it is usual to deal with cosmological dynamics where the time is given in terms of the Hubble constant H, namely

$$t_\mathrm{H}=\mathrm{H}^1\frac{e^{3\mathrm{\Omega }}}{\pi _\mathrm{\Omega }}.$$
(19)
Note that on the constraint surface the time (13) can be rewritten as

$$t\frac{ke^{2\mathrm{\Omega }}}{\pi _\mathrm{\Omega }}\frac{\mathrm{\Lambda }e^{4\mathrm{\Omega }}}{\pi _\mathrm{\Omega }},$$
(20)
which has an analogous form.

The problem of time in quantum gravity is an important open question. In this paper we have identified time variables for a particular minisuperspace model which is of interest in quantum cosmology by means of a systematic procedure to deparametrize the constrained system, and we have given a consistent quantization with a clear notion of time. We have pointed out the relevance of the study of the structure of the phase space in terms of the extrinsic curvature within the context of the procedure to find a global phase time for deparametrizing cosmological models which do not admit an intrinsic time.

Acknowledgement

This work was supported by CNEA and CONICET.