text:
# What Shapes the Luminosity Function of Galaxies?

## 1 Introduction

The luminosity function of galaxies is one of the most basic properties of the galaxy population; yet it contains many valuable clues to the process of galaxy formation. The basic physical mechanisms which determine the form of the luminosity function were first described by Rees & Ostriker (1977) and White & Rees (1978). In this picture, galaxy formation is regulated by the rate at which gas is able to cool in the parent dark matter halos. These authors suggested that the sharp cut-off in the galaxy luminosity function arose from the long cooling times of gas in high mass halos (or high mass protogalaxies in the case of Rees & Ostriker). The model has been developed by many authors to follow in great detail the formation of galaxies in a hierarchical universe. Key improvements are the inclusion of galaxy merging and the evolution of stellar populations (White & Frenk, 1991; Cole, 1991; Kauffmann, White & Guiderdoni, 1993; Lacey et al., 1993; Cole et al., 1994; Kauffmann et al., 1999; Somerville & Primack, 1999; Cole et al., 2000; Benson et al., 2002). Such models are now being strongly tested by high-precision measurements of the galaxy luminosity function from large redshift surveys such as the 2dFGRS and 2MASS (Cole et al., 2001; Kochanek et al., 2001).

While the key physics of gas cooling and merging are now thought to be modeled with reasonable accuracy (Benson et al., 2001, 2002; Yoshida et al., 2002; Helly et al., 2003; Voit et al., 2002), other physics crucial to establishing the shape of the luminosity function remain poorly understood. The first uncertainty is the “feedback” needed to regulate the formation of dwarf galaxies, and hence reconcile the rather shallow slope of the faint end of the observed luminosity function with the relatively steep mass function of dark matter halos. While outflows of gas from galaxies have been observed at both low and high redshift (Martin, 1999; Pettini et al., 2002), the complex physics at work has not yet been understood in detail, and most models of galaxy formation simply adopt phenomenological rules to describe their effects. Previous work has typically assumed that the dominant feedback mechanism is the reheating of cold gas in the disk to the temperature of the diffuse gas halo (White & Frenk, 1991; Cole et al., 1994; Kauffmann et al., 1999; Efstathiou, 2000), although complete expulsion of disk gas from the halo was considered by Somerville & Primack (1999). However, although this form of feedback solves the problem with the faint-end slope of the luminosity function (as it was originally designed to do), it creates a second problem matching the very sharp (exponential with luminosity) cut-off seen at the bright end. In addition, the effectiveness with which the gas needs to be reheated seems excessive compared to observations of galactic outflows (Martin, 1999) and to what is seen in simulations of this process in individual galaxies (Strickland & Stevens, 2000). These papers suggest that the mass in the outflow is comparable to the gas mass that is turned into stars.

All current models of galaxy formation, calculated using either gas-dynamical simulations (Pearce et al., 2001; Kay et al., 2002; Katz et al., 2002; Weinberg et al., 2003) or semi-analytic techniques (e.g. Kauffmann et al. 1999; Somerville & Primack 1999; Cole et al. 2000), exhibit strong gas cooling in the central regions of groups and clusters. This leads to the formation of extremely bright galaxies (which are never seen in reality) unless some additional suppression of the gas cooling is assumed. The suppression mechanisms which have been considered in semi-analytical models are: simply switching off cooling and/or star formation in the most massive halos, e.g. Kauffmann et al. (1999); redistributing the gas within the halo so that it becomes so rarefied that its cooling time is longer than the age of the Universe, e.g. Cole et al. (2000); or suppressing cooling until the halo is completely virialized, e.g. van Kampen et al. (1999). The model presented in Cole et al. (2000) adopted a low value for the cosmic baryon fraction ($`\mathrm{\Omega }_\mathrm{b}=0.02`$), and a model for the halo gas distribution in which the core radius was a function of the amount of gas that had cooled in previous generations of halos. In combination, these assumptions were able to produce a good match to the galaxy luminosity function; however, both are now disfavored by current observational data (e.g. O’Meara et al. 2001; Allen et al. 2001; Ettori et al. 2002). If the more recent, higher value of the baryon abundance ($`\mathrm{\Omega }_\mathrm{b}0.04`$) is adopted, the growing core radius has little effect on the cooling rates (since the cooling radius is then significantly beyond the typical core radii).

In this paper, we investigate physical processes that may be responsible for shaping the bright end of the luminosity function. In addition to the conventional feedback process, in which cold gas in the disk is reheated to the temperature of the diffuse gas halo, the feedback energy may be used to regulate the formation of the galaxy in two further ways. While the wind flowing from the gas disk may contain a relatively small mass, the energy in the disk outflow may be transferred to the existing hot gas halo, causing it to expand within the halo potential. This makes the central gas more diffuse, lengthens its cooling time and so reduces the rate at which gas is supplied to the cold disk. This type of process is seen in simulations of the effect of injection of relativistic radio-emitting plasma in the centres of clusters (Quilis et al., 2001; Brüggen et al., 2002). An alternative is to assume that the gas expelled from the disk does not mix with the virialized diffuse gas halo, and that the energy per particle is sufficiently high that the material escapes the confining potential and is never recaptured.

A further means by which the supply of cold gas to the galaxy disk may be reduced is suggested by recent Chandra and XMM observations of galaxy clusters (Peterson et al., 2001; Tamura et al., 2001; Fabian et al., 2001; Johnstone et al., 2002; Böhringer et al., 2002; McNamara et al., 2001; Nulsen et al., 2002). These have led to a revival of the idea that thermal conduction may be an important source of heating in the central regions of clusters (Narayan & Medvedev, 2001; Gruzinov, 2002; Voigt et al., 2002). Indeed, the conductivity has been inferred to be close to the Spitzer rate expected for a fully ionized plasma (Fabian, Voigt & Morris, 2002). As first suggested by Fabian, Voigt & Morris (2002), heating due to conduction could plausibly counteract the excessive cooling predicted to occur in the most massive halos and thereby explain the dearth of highly luminous galaxies in the Universe.

In the remainder of this paper, we will demonstrate quantitatively the importance of each physical mechanism and its role in setting the form of the luminosity function. We will demonstrate that the milder forms of feedback are unable to account simultaneously for the sharp break in the luminosity function and for the flat faint-end slope. The two processes discussed above, thermal conduction and superwinds, may provide the answer to this problem. Thermal conduction is capable of suppressing the formation of the brightest objects if the conductive efficiency is high. Alternatively, a model that includes energetic gas expulsion is also able to produce reasonable fits to the observed luminosity function, but the energy required to power this “superwind” is larger than that available from supernova explosions and stellar winds.

## 2 The Model

We compute the luminosity function of galaxies in a cold dark matter Universe using a development of the semi-analytic model, galform, described by Cole et al. (2000). Our changes to the model concern the treatment of the gas distribution in halos, and the treatment of feedback.

### 2.1 Feedback

In Cole et al. (2000), the star formation rate was affected by feedback generated when cold gas is reheated by energy injected by supernovae explosions and stellar winds. This process transfers gas from the star-forming disk of the galaxy to a diffuse corona or halo. The rate at which cold gas is reheated is assumed to be related to the star formation rate, $`\dot{M}_{}`$, by:

$$\dot{M}_{\mathrm{reheat}}=\left(\frac{V_{\mathrm{disk}}}{V_{\mathrm{hot}}}\right)^\alpha \dot{M}_{},$$
(1)
where $`V_{\mathrm{disk}}`$ is the circular speed at the half-mass radius of the gas disk, and $`V_{\mathrm{hot}}`$ and $`\alpha `$ are adjustable parameters. Here we include this feedback mechanism, but we also consider two further schemes which correspond to different processes by which the energy from supernova explosions and stellar winds can regulate the gas cooling and star formation rates. A further fraction of the total energy will be radiated away and so is not available for feedback. We assume that the way the energy is apportioned amongst the three processes is fixed by the local supernova environment, and does not depend on the mass of the galaxy or on the instantaneous star formation rate.111In principle the supernova’s local environment—the interstellar medium—could depend on global properties of the galaxy, such as its mass, but we neglect any such dependence here. Thus, feedback in our model is specified by three parameters: $`ϵ_{\mathrm{reheat}},ϵ_{\mathrm{halo}},ϵ_{\mathrm{sw}}`$. It is useful to note that the total energy available from supernovae is of order $`10^{49}\text{erg}`$ for each Solar mass of stars that is formed, for a standard IMF.

1. Disk reheating. This is the feedback scheme explored by Cole et al. (2000). Gas is removed from the disk of the galaxy at a rate proportional to the star formation rate. This gas can be thought of as either being ejected at close to, but below the escape velocity of the halo, or as being heated to the halo virial temperature. In the first case, the gas will subsequently be heated back to the halo virial temperature, through mixing with the halo gas, or through further shocks when a new halo forms. In either case, the gas does not leave the halo and so is available for cooling in the future. We parameterize the energy invested in reheating the cold disk gas as $`ϵ_{\mathrm{reheat}}\mathrm{\hspace{0.17em}10}^{49}\text{erg}`$ per Solar mass of stars formed. If we approximate the gravitational potential wells of galaxies as being self-similar, then the energy required to eject the gas from the disk should scale as $`V_{\mathrm{disk}}^2`$, and one obtains formula (1) above with $`\alpha =2`$. Specifically, if we assume that gas is ejected with energy per unit mass equal to $`V_{\mathrm{disk}}^2`$, then the rate at which gas is reheated is given by:

$$\dot{M}_{\mathrm{reheat}}=\frac{5.6ϵ_{\mathrm{reheat}}}{\left(V_{\mathrm{disk}}/300\text{km}\text{s}^1\right)^2}\dot{M}_{}.$$
(2)
We therefore have the following relation between $`ϵ_{\mathrm{reheat}}`$ and $`V_{\mathrm{hot}}`$:

$$ϵ_{\mathrm{reheat}}=0.18\left(\frac{V_{\mathrm{hot}}}{300\text{km}\text{s}^1}\right)^2,$$
(3)
If the circular speed of the disk is equal to $`V_{\mathrm{hot}}`$, the mass of gas reheated is equal to the mass formed into stars. Feedback of this form is motivated by the study of Efstathiou (2000), although the efficiency assumed here is somewhat higher.
2. Energy injection into the gaseous halo. The feedback energy is assumed to redistribute the hot, diffuse gas so that the central density of the gas is reduced, and the central cooling time is increased. The energy injected per unit mass of stars formed is parameterized as $`ϵ_{\mathrm{halo}}10^{49}\text{erg}M_{}^1`$. The nature of this type of feedback process can be understood as follows: the total energy injected into the halo grows with time until it becomes comparable to the gravitational binding energy of the gaseous halo. At this point, the cooling rate drops dramatically. Since the binding energy of a halo is proportional to $`M_{\mathrm{halo}}^{5/3}`$, the fraction of the total gaseous mass which must be turned into stars to achieve this balance is strongly dependent on halo mass. Larger halos are therefore able to cool a larger fraction of their baryon reservoir to form galaxies. This mechanism is described in more detail in Bower et al. (2001) and, as noted above, is motivated by simulations of the injection of radio plasma in cluster centres (Quilis et al., 2001; Brüggen et al., 2002).
3. Gas expulsion from the halo. In this scheme, an energy of $`ϵ_{\mathrm{sw}}\mathrm{\hspace{0.17em}10}^{49}\text{erg}M_{}^1`$ is invested in heating a small mass of the cold gas disk to an energy much greater than the binding energy of the halo. If the energy is contained within this super-heated phase (and is not shared with the diffuse gas halo, as assumed in mechanism 2 above), the wind may be able to escape completely from the halo. This mechanism is motivated by the strong, high velocity winds inferred to exist around Lyman-break galaxies (Pettini et al., 2002; Adelberger et al., 2003) and has recently been studied theoretically by Shu, Mo & Mao (2003). We assume that the energy is injected with a thermal distribution of particle energies, i.e. so that the fraction of mass ejected with energy in the range $`E`$ to $`E+\mathrm{d}E`$ is proportional to $`\mathrm{exp}(E/E_{\mathrm{av}})\mathrm{d}E`$, where $`E_{\mathrm{av}}`$ is the mean ejection energy per particle. The fraction of the superwind gas with sufficient energy to escape the halo is computed using this distribution. Some of the expelled gas will be recaptured when a deeper potential well forms. To estimate this, we calculate the fraction of the expelled gas with energy less than the depth of the new potential well, and allow this fraction to be recaptured. This process continues as deeper potential wells form until either all of the gas is recaptured, or $`z=0`$ is reached.

We assume that the mass flux of material in the wind as it leaves the disk, $`\dot{M}_{\mathrm{superwind}}`$, is proportional to the star formation rate, with coefficient of proportionality $`\beta _{\mathrm{SW}}`$. The energy in the superwind as it leaves the disk is given by

$$\dot{E}_{\mathrm{superwind}}=10^{49}ϵ_{\mathrm{SW}}\dot{M}_{}\text{ ergs},$$
(4)
where $`10^{49}ϵ_{SW}`$ ergs is the energy per unit mass of star formation. The characteristic specific energy of the wind is simply $`E_{\mathrm{av}}=\dot{E}_{\mathrm{superwind}}/\dot{M}_{\mathrm{superwind}}\frac{1}{2}V_{\mathrm{SW}}^2`$, and so the characteristic velocity of the wind is

$$V_{\mathrm{SW}}=1002\sqrt{ϵ_{\mathrm{SW}}/\beta _{\mathrm{SW}}}\text{ km/s}.$$
(5)
We assume that material ejected from the disk requires an energy per unit mass $`\lambda _{\mathrm{SW}}V_{\mathrm{disk}}^2`$ in order to escape also from the halo. Therefore, in the case in which all of the wind material has the same energy, superwinds will be driven from halos with $`V_{\mathrm{disk}}<709\sqrt{ϵ_{\mathrm{SW}}/\beta _{\mathrm{SW}}\lambda _{\mathrm{SW}}}`$ km/s, and not from deeper potential wells. For typical halos in our model $`\lambda _{\mathrm{SW}}=2.9`$ gives a good estimate of the energy required to escape the combined gravitational pull of the galaxy and its dark matter halo, and so we use this value throughout this work. If we assume a distribution of energies in the wind, then there will be a smooth transition, as $`V_{\mathrm{disk}}`$ increases, between the regime in which the wind escapes and the regime in which it is retained. Thus, the actual mass flux in the superwind escaping from the halo is given by

$$\dot{M}_{\mathrm{superwind}}=f_{\mathrm{eject}}(E_{\mathrm{av}}/\lambda _{\mathrm{SW}}V_{\mathrm{disk}}^2)\beta _{\mathrm{SW}}\dot{M}_{},$$
(6)
where, for our chosen energy distribution, $`f_{\mathrm{eject}}(x)=\mathrm{exp}(x)`$. The material which does escape may be recaptured by larger halos forming at later times. Our estimate of the fraction of the wind which escapes and that which can be recaptured by larger halos is described in Appendix A.

We regard this mechanism as being more uncertain than mechanisms 1 and 2 because the expelled gas must leave the halo without sharing its energy with the diffuse hot component. It must therefore punch a well collimated hole through the halo, or remain contained within a buoyant bubble which is convected out of the halo (Quilis et al., 2001; Springel & Hernquist, 2003; Kay et al., 2002).

In general, it seems plausible that some fraction of the total feedback energy will be processed into each one of these forms, and we consider models which include a combination of these processes on an equal footing with models that include only one.

### 2.2 Halo Gas Distribution

The hot gaseous component in dark matter halos is assumed initially to have a density profile, $`\rho _\mathrm{g}(r)`$, given by the $`\beta `$-model, ie, at radius $`r`$ in the halo (whose dark matter density profile is assumed to have the NFW profile, Navarro, Frenk & White 1996, 1997),

$$\rho _\mathrm{g}(r)=\frac{\rho _0}{[1+(r/r_\mathrm{c})^2]^{3\beta /2}},$$
(7)
where $`\rho _0`$ is the density at the centre of the halo, $`r_\mathrm{c}`$ is the radius of the “core” and $`\beta `$ is a parameter which sets the slope of the profile on scales larger than $`r_\mathrm{c}`$.

Departing from the prescription of Cole et al. (2000), we adopt a gas density profile in the absence of energy injection with fixed $`r_\mathrm{c}=0.07r_{\mathrm{vir}}`$ (and $`\beta =2/3`$) in all halos. This provides a reasonable match to gas-dynamic simulations of non-radiative gas in clusters (e.g. Eke, Navarro & Frenk 1998) and to the observed X-ray profiles of relaxed clusters (e.g. Allen et al. 2001). The normalization of the profile is determined by the total diffuse gas mass remaining in the halo, and the temperature of the gas is set assuming hydrostatic equilibrium. As a boundary condition, we set the temperature at the virial radius equal to the virial temperature (eqn. 4.1 in Cole et al. 2000). This “default” profile is modified if the diffuse gas gains further energy (“excess energy”) as a result of energy injection (process 2 in Section 2.1). We use the prescription described in Bower et al. (2001), in which the excess energy first causes the slope of the gas profile to decrease down to a minimum value of $`\beta _{\mathrm{min}}=0.2`$, after which it increases the boundary temperature of the gas halo. The effect of the excess energy is to decrease the central density of the gas, lengthening its cooling time. Mass is conserved by pushing some of the diffuse gas outside the halo; however, this gas can be recaptured as the total halo mass (and thus the gravitational binding energy) increases. We assume that the excess energy is conserved during mergers between halos, although the results are not qualitatively affected by a small dilution or amplification of energy during halo mergers.

We define the effective cooling time at radius $`r`$, $`t_{\mathrm{cool}}^{}(r)`$, as the maximum true cooling time (i.e. that defined by Cole et al. 2000) occurring at smaller radii, plus the free-fall time from $`r`$ to the halo centre. This ensures that $`\mathrm{d}t_{\mathrm{cool}}^{}/\mathrm{d}r>0`$, so the cooling radius is always a smoothly increasing function of radius. Experiments with different approaches show that the results we present here are not sensitive to the details of this prescription.

### 2.3 Conduction

Conduction in the ionized gas can transport energy into the inner regions of the halo, effectively increasing the cooling time of the gas there. The rate at which energy is deposited into the shell between radii $`r`$ and $`r+\mathrm{d}r`$ is given by:

$$\mathrm{\Sigma }\mathrm{d}r=4\pi \frac{\mathrm{d}}{\mathrm{d}r}\left(\kappa r^2\frac{\mathrm{d}T}{\mathrm{d}r}\right)\mathrm{d}r,$$
(8)
where the conductivity, $`\kappa `$, may depend on radius through its temperature dependence. We approximate $`\mathrm{\Sigma }`$ as

$$\mathrm{\Sigma }=\alpha _{\mathrm{cond}}4\pi \kappa _\mathrm{s}T,$$
(9)
where $`\kappa _\mathrm{s}`$ is the Spitzer conductivity (Spitzer, 1962), and $`\alpha _{\mathrm{cond}}`$ is a parameter that absorbs the dependence on the shape of the temperature profile as well as any difference between the actual conductivity and the Spitzer rate. For a power-law temperature profile, $`Tr^a`$, and conductivity of the Spitzer form, $`\kappa _\mathrm{s}T^{5/2}`$, $`\alpha _{\mathrm{cond}}=f_{\mathrm{sp}}a(1+7a/2)`$ where $`f_{\mathrm{sp}}`$ is the ratio of the conductivity to the Spitzer value. Adopting a temperature profile with $`a=0.4`$ as suggested by recent X-ray observations (Voigt et al., 2002) gives $`\alpha _{\mathrm{cond}}=0.96f_{\mathrm{sp}}`$. Adopting a linear temperature gradient gives $`\alpha _{\mathrm{cond}}=4.5f_{\mathrm{sp}}`$.

The heating rate due to conduction is subtracted from the radiative cooling rate to give a net cooling rate for the gas. This net cooling rate is used to compute the cooling time. Conduction causes the cooling radius to become smaller than in the standard model. The result is a suppression of cooling in hot halos, and the mass at which this effect becomes important is determined by the parameter $`\alpha _{\mathrm{cond}}`$.

### 2.4 Photoionization and Merging

While we employ the detailed model of galaxy merging developed by Benson et al. (2002), we choose to use a much simpler model of the effects of reionization than used in that paper, in order to incorporate them easily into our calculation of the galaxy luminosity function. We simply assume that galaxy formation is completely suppressed by reionization in dark matter halos with circular velocities below $`V_{\mathrm{reion}}`$ after redshift $`z_{\mathrm{reion}}`$. Unless otherwise stated, we adopt $`V_{\mathrm{reion}}=50\text{km}\text{s}^1`$ and $`z_{\mathrm{reion}}=6`$. With this choice, this simple model matches the results of the full calculation of Benson et al. (2002) quite well.

## 3 Results

Throughout this paper, we use the same parameter values adopted by Benson et al. (2002), with the exception of a larger baryon fraction corresponding to $`\mathrm{\Omega }_\mathrm{b}=0.045`$, consistent with constraints from Big Bang nucleosynthesis (O’Meara et al., 2001). We assume a $`\mathrm{\Lambda }`$CDM universe with mean matter density $`\mathrm{\Omega }_0=0.3`$, cosmological constant term $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, Hubble constant222Here and below $`h`$ denotes the Hubble constant in units of 100km s-1 Mpc-1 $`H_0=70`$km s-1 Mpc-1 and linear fluctuation amplitude on spheres of radius 8h-1 Mpc, $`\sigma _8=0.9`$. Of these parameters, the uncertainty in the value of $`\sigma _8`$ has most impact on the model results. While several studies support a value of $`\sigma _80.9`$ (e.g., Bacon et al. (2002), Hoekstra et al. (2002) \[gravitational lensing\]; Eke, Cole & Frenk (1996), Vianna et al. (2002) \[cluster abundance\]; Spergel et al. (2003) \[cosmic microwave background (CMB)\]; Sievers et al. (2003) \[Sunyaev-Zeldovich effect\]), other recent analyses have suggested lower values, $`\sigma _80.7`$ (e.g., Peacock (2003) \[large scale structure\]; Melchiori et al. (2002) \[CMB\]; Jarvis et al. (2003) \[gravitational lensing\]; Allen et al. (2003), Smith et al. (2003) \[cluster abundance\]). A discussion of recent results may be found in Wang et al. (2003). Except where specified, we show models based on $`\sigma _8=0.9`$; however, (as we shall see) taking the lower value, $`\sigma _8=0.7`$, considerably eases the energy budget and/or reduces the conduction efficiency required to match the galaxy luminosity function. Throughout, we use the halo mass function derived from N-body simulations by Jenkins et al. (2001)333Note that we use the Jenkins et al. (2001) mass function to compute the abundances of halos at $`z=0`$, but retain the Press-Schechter approximation when computing halo merger trees. An improved calculation would use a self-consistent computation of halo merging histories., instead of the Press-Schechter mass function used by Cole et al. (2000).

We compare our model with recent determinations of the K-band luminosity function assuming a Kennicutt stellar initial mass function (Kennicutt, 1983). In order to facilitate comparison between models, we have kept the IMF fixed, and assumed a negligible fraction of brown dwarf stars444We set $`\mathrm{{\rm Y}}=1`$ in the notation of Cole et al. (2000). We choose the K-band in order to minimize the sensitivity of our results to recent star formation and to dust obscuration. The model of Cole et al. (2000) includes a detailed and fully self-consistent calculation of dust extinction which is used in this work. We find, however, that dust-obscuration has a negligible effect on our results for the K-band luminosity function (typically shifting the bright end faintwards by less than one tenth of a magnitude.) For the observational comparison, we use the local K-band luminosity functions measured by Cole et al. (2001) and Kochanek et al. (2001) (both based on the 2MASS survey) and the local ($`z<0.1`$) luminosity function derived from the much deeper K-band survey of Huang et al. (2002). The analysis by Huang et al. suggests a faint-end slope ($`\alpha _\mathrm{K}=1.37\pm 0.10`$), steeper than the values found by Cole et al ($`\alpha _\mathrm{K}=0.93`$) and by Kochaneck et al. ($`\alpha _\mathrm{K}=1.09`$). These latter two are also in good agreement with the faint end slope of the z-band luminosity function measured by Blanton et al. ($`\alpha _\mathrm{z}=1.08`$) from the SDSS survey. The z-band data should also be little affected by residual star formation and dust extinction, but have a deeper surface brightness limit. These discrepancies indicate that there remain significant systematic uncertainties in current measurements of the faint end of the K-band luminosity function, perhaps due to luminosity from low-surface brightness regions of galaxies being missed as recently suggested by Andreon (2002).

### 3.1 A Model with constant mass-to-light ratio

In Fig. 3 we show the simplest possible model of the luminosity function which we call Model 1 (shown as the dashed line). In this model, the mass function of dark matter halos (Jenkins et al., 2001) has been converted into a luminosity function simply by assuming a fixed mass-to-light ratio ($`M/L_\mathrm{K}=11M_{}/L_{\mathrm{K},}`$), chosen so as to match the knee of the observed luminosity function. As is well known, this produces a luminosity function which is much steeper at the faint end than is observed, and also fails to cut off at bright magnitudes (the halo mass function does possess a cut-off, but it occurs at much higher mass and lower abundance than shown in the plot).

### 3.2 Model including Cooling only

White & Rees (1978) argued that the difference between the halo mass and the galaxy luminosity functions is due to the dependence of the gas cooling time on halo mass and to feedback processes. We use the galform semi-analytic model to follow gas cooling and star formation in a merging hierarchy of dark matter halos in the $`\mathrm{\Lambda }`$CDM cosmology. In order to illustrate the simplest possible model first, we do not include photoionization suppression, feedback, galaxy merging or conduction. The result is Model 2 (shown as a dotted line in Fig. 3). It clearly displays the “overcooling problem”: gas has cooled into the smallest halos resolved in the calculation, producing an overabundance of faint galaxies. As a result, the faint-end slope of the luminosity function is much too steep. It should be noted that the results for this cooling-only model are sensitive to the mass resolution of the galform merger trees. The cooling time for halo gas decreases with decreasing halo mass down to halo virial temperatures around $`10^4`$K. Below this scale, cooling becomes ineffective unless molecular hydrogen is abundant. This means that for the luminosity function to be fully converged at all galaxy luminosities, we would need to resolve halos with virial temperatures as low as $`10^4`$K at all redshifts at which there is significant cooling. To obtain the results shown for Model 2, we ran galform at the highest mass resolution that was computationally feasible, which comes close to resolving $`10^4`$K halos at $`z=0`$. We believe that our results are substantially converged, but we cannot be certain that they would not change if the mass resolution were increased further.

### 3.3 Model with photoionization and merging

The problem of mass resolution encountered in Model 2 is effectively removed once we include feedback processes which are effective in low-mass halos. Model 3 (long-dashed line in Fig. 3) shows the effects of including photoionization suppression. As described by Benson et al. (2002), the formation of low-mass galaxies is suppressed after the Universe is reionized and, as we noted in §2, we adopt a simplified model of this process which nevertheless is a good approximation to a full calculation. The resolution of the galform merger trees is then sufficient to follow the lowest mass halos that are able to form galaxies, resulting in a converged solution. However, photoionization suppression alone is unable to produce a sufficiently flat faint end in the luminosity function. In addition, Model 3 also produces more bright galaxies, since less gas has been locked into stars in the smallest halos. When both galaxy merging (which was artificially switched off in Models 1 and 2) and photoionization are included, as in Model 4 (solid line in Fig. 3), the faint end remains too steep and even more bright galaxies are produced. In the models that follow, we investigate the effects of including feedback processes in addition to photoionization and merging.

### 3.4 Model with feedback—disk reheating

A solution to the faint-end problem is illustrated by Model 5.3 in Fig. 3.3. Here, feedback is included through the reheating of disk gas in star-forming galaxies. We use the standard prescription of Cole et al. (2000), but with a larger value of $`ϵ_{\mathrm{reheat}}=0.41`$ (equivalent to $`V_{\mathrm{hot}}=450\text{km}\text{s}^1`$) which is required in order to obtain a similar faint end slope to that in Cole et al. for the larger value of $`\mathrm{\Omega }_\mathrm{b}`$ assumed in this work. This form of feedback flattens the luminosity function considerably, resulting in reasonably good agreement with the observed faint end. While the slope is not as flat as that measured by Cole et al. (2001) or Blanton et al. (2003), it is in good agreement with the steeper slope reported by Huang et al. (2002). This achievement carries a price, however—the overabundance of bright galaxies (formed through excessive cooling in massive halos) is exacerbated, as there is now a much greater mass of diffuse hot gas remaining in the larger halos. This gas is sufficiently dense that the central regions are able to cool; consequently, Model 5.3 produces far too many bright galaxies. This result depends little on the choice of $`\sigma _8`$. Adopting $`\sigma _8=0.7`$ makes the brightest galaxies only 0.5 mag fainter. This clearly demonstrates a longstanding problem in semi-analytic models: previous calculations have either assumed low values of $`\mathrm{\Omega }_\mathrm{b}`$, or have invoked rather artificial ways to prevent the cooling which forms these overluminous objects.

The energy requirements of Model 5.3 are substantial, but not excessive. The reheating energy of $`0.41\times 10^{49}\text{erg}M_{}^1`$ should be compared to the total energy available from supernovae explosions which is approximately $`0.7\times 10^{49}\text{erg}M_{}^1`$ for a Salpeter IMF, or $`0.9\times 10^{49}\text{erg}M_{}^1`$ for the Kennicutt IMF adopted in galform. In a halo of circular velocity 250 km/s, the mass of gas reheated is more than 3 times the mass of gas formed into stars. If the level of reheating is reduced, as illustrated by models 5.1 ($`ϵ_{\mathrm{reheat}}=0.03`$) and 5.2 ($`ϵ_{\mathrm{reheat}}=0.13`$), then the formation of small galaxies is not suppressed sufficiently to match the observed luminosity function.

### 3.5 Models with energy injection

In Models 6.1–6.3 (Fig 3.4), we investigate the effect of assuming that a fraction of the supernova and stellar wind energy heats the diffuse gas halo, causing it to expand. This form of feedback suppresses the formation of both bright and faint galaxies, but it does not produce a sharp break in the luminosity function. Models 6.1 and 6.2 illustrate how as the energy spent in heating the diffuse halo increases, the break in the luminosity function becomes less pronounced. If the heating is made even stronger, the resulting luminosity function approaches a power-law (Model 6.3). This form of feedback clearly cannot solve the problem of overproduction of bright galaxies.

### 3.6 Model with conduction

We now consider two possible schemes that are capable of producing a good match to the observed luminosity function. The first involves balancing radiative cooling with thermal conduction. The second involves expelling gas from dark matter halos at such high energies that it is subsequently unable to cool.

Thermal conduction is expected to imprint a special scale on the galaxy population because of the strong temperature dependence of the Spitzer conductivity rate. Model 7.3 (solid line in Fig. 3.6) shows the result of including conduction with $`\alpha _{\mathrm{cond}}=25`$ in Model 5.3 ($`ϵ_{\mathrm{reheat}}=0.41`$, $`ϵ_{\mathrm{halo}}=0`$). Such a high value of the conductivity is indeed effective in suppressing the formation of the most massive galaxies since it prevents efficient gas cooling in group and cluster-sized halos. The brightest galaxies in the luminosity function are instead built through mergers. The result is a rather good match to the observed galaxy luminosity function. However, the conduction efficiency that we have assumed is extremely high. To operate at the required level, we must assume both that the conductivity is not suppressed below the Spitzer value, and that the effective temperature gradient is steeper than $`Tr^{2.5}`$ in the region of the cooling radius. Note that the Spitzer formula for thermal conductivity in an ionized plasma breaks down if the conductivity becomes too high, i.e. if the conduction “saturates” as described by Cowie & McKee (1977). The models shown in this paper do not take account of this saturation limit. However, using the estimate of the saturated heat flux from Cowie & McKee (1977), we have checked that our results are not significantly affected by saturation (for model 7.3, which has the most extreme conduction of all our models, there is only a small increase in the number of the very brightest galaxies).

In Models 7.1 and 7.2, we show the effect of assuming a more modest conduction efficiency ($`\alpha _{\mathrm{cond}}=1.0`$ and 0.1 respectively). In these models, conduction is not sufficient to suppress cooling in the larger halos adequately.

If we adopt a lower value for $`\sigma _8`$, however, a lower conduction efficiency gives a reasonable match to the observed luminosity function. Model 7.4 shows the luminosity function for the case $`\sigma _8=0.7`$ and $`\alpha _{\mathrm{cond}}=7`$. This conduction efficiency could be achieved if the temperature gradient was $`Tr^{1.3}`$ and the conduction was only slightly suppressed below the Spitzer value. Although this is still a high rate of conduction, it offers a promising route for explaining the bright end of the galaxy luminosity function.

### 3.7 Model with superwinds

The expulsion of gas from halos at high energy can, in principle, strongly suppress the formation of later generations of galaxies, hence affecting the shape of the luminosity function. Starting from Model 5.2, we add further feedback energy that expels cold gas completely, not only from the disk but also from the halo. The superwind must have high energy in order that the expelled material not be recaptured by more massive halos. The effect of a low power superwind is illustrated by Model 8.1 (dashed line in Fig. 3.6). This model, with $`ϵ_{\mathrm{sw}}=0.27`$ and $`\beta _{\mathrm{SW}}=3`$, has a relatively weak superwind. This gas expulsion is in addition to the reheating of cold disk gas ($`ϵ_{\mathrm{reheat}}=0.13`$); we have assumed that there is no heating of the diffuse halo ($`ϵ_{\mathrm{halo}}=0.0`$). Although the winds eject a large amount of gas, most of the material is recaptured as larger halos collapse and the luminosity function differs little from Model 5.2.

In Model 8.2 (dotted line), we have set $`ϵ_{\mathrm{sw}}=5.0`$ and $`\beta _{\mathrm{SW}}=1`$, corresponding to a mean energy per superwind particle of $`E_{\mathrm{av}}=15`$keV. Such an energetic wind is required to ensure that very little material is recaptured by group halos. In this model, the superwind dominates the feedback energy budget; indeed, the total energy required ($`5.13\times 10^{49}\text{erg}M_{}^1`$) significantly exceeds that available from supernovae alone. The model comes much closer to matching the luminosity function, but still overproduces bright galaxies. We can increase the superwind energy further, but a factor of 2 increase only results in a small improvement in the luminosity function. If we include conduction as well as superwinds, it is, of course, possible to improve the match but a high conduction efficiency ($`\alpha _{\mathrm{cond}}1`$) is still needed. Increasing the mass loading of the wind substantially (by increasing $`\beta _{\mathrm{SW}}`$) results in too few galaxies around the knee of the luminosity function.

As we found in the case of conduction, the luminosity function can be matched more easily if we adopt a lower value for $`\sigma _8`$. The case $`ϵ_{\mathrm{sw}}=5.0`$, $`\beta _{\mathrm{SW}}=1`$, $`\sigma _8=0.7`$ is illustrated by Model 8.3 (solid line). Given the uncertainties of our recapture prescription, this model gives a reasonable match to the luminosity function; it has a strong break at the correct luminosity, and overproduces bright objects only marginally. There are a variety of ways to further improve the match to observations: we could increase the energy injection (but $`ϵ_{\mathrm{sw}}>10.0`$ is required) or increase the mass loading of the wind so that the curve is shifted faintwards. An alternative strategy is to allow for a low level of conduction: $`\alpha _{\mathrm{cond}}1`$ is sufficient to produce a significant improvement in the match to the luminosity function when $`\sigma _8=0.7`$ and superwinds are present.

The parameters $`\alpha _{\mathrm{cond}}`$ and $`ϵ_{\mathrm{SW}}`$ are highly degenerate in their effects on the luminosity function (i.e. increasing either suppresses the bright end). While current computational limitations make it prohibitively expensive to perform an accurate $`\chi ^2`$ fit of the model parameters to the data, a crude estimate of the $`\chi ^2`$ surface for these two parameters shows that the data prefer models with strong superwinds ($`ϵ_{\mathrm{SW}}6`$) and high conductivity ($`\alpha _{\mathrm{cond}}30`$) for a model with $`\sigma _8=0.93`$. Lowering $`\sigma _8`$ to $`0.7`$ reduces the requirements to $`ϵ_{\mathrm{SW}}3`$ and $`\alpha _{\mathrm{cond}}25`$. Further investigation of the $`\chi ^2`$ surface would require a Bayesian prior to specify formally a physically allowed range for this parameter.

### 3.8 Other Considerations

While our primary goal in this paper is to examine the luminosity function of galaxies, it is prudent to check whether our models are in reasonable agreement with other basic properties of the galaxy population such as the Tully-Fisher relation and the galaxy autocorrelation function. We compare the models that best fit the galaxy luminosity function to these observables. We retain the model parameters found earlier and we do not attempt to adjust these or any other parameters in this comparison. We defer a more detailed comparison of our models with a wider range of observational constraints to a future paper.

#### 3.8.1 Tully-Fisher Relation

Simultaneously matching the galaxy luminosity function and the Tully-Fisher relation has been a long-standing problem for CDM-based semi-analytic models (see, for example, White & Frenk 1991, Kauffmann, White & Guiderdoni 1993 and Cole et al. 2000). In particular, Cole et al. (2000) found that while their best model reproduced the observed slope and scatter in the Tully-Fisher relation, the predicted circular velocities for galaxy disks were about 30% larger than the values measured for the data.

Figure 3.8.1 shows the Tully-Fisher relation predicted by four of our models (including those which most successfully match the observed luminosity function). It is readily apparent that our models all perform reasonably well in that they reproduce the slope and scatter of the observed relation, although they do not match the zero-point exactly. Model 8.3 does best, being offset by only approximately 10% in circular velocity over most of the range shown. This model, in fact, performs rather better than that of Cole et al. (2000), although an even better match would be desirable. This may involve revising the assumption of adiabatic invariance made in calculating the effects of the disk on the inner structure of the dark matter halo (see Cole et al. 2000) or considering the effects of the angular momentum of the superwind material or the details of our star formation prescription.

#### 3.8.2 Galaxy Correlation Function

One of the most remarkable successes of semi-analytic models of galaxy formation in the $`\mathrm{\Lambda }`$CDM cosmology is their ability to match the observed galaxy autocorrelation function (Kauffmann et al. 1999; Benson et al. 2000). As pointed out by Benson et al. (2000), the clustering of galaxies may be understood physically by reference to the halo occupation distribution (see, for example, Peacock & Smith 2000, Cooray & Sheth 2002), which specifies the number of galaxies that occupy dark matter halos of a given mass. We show, in Fig. 3.8.2, the mean number of galaxies brighter than $`M_\mathrm{B}5\mathrm{log}h=19.5`$ per halo as a function of halo mass for four of our models. For reference we also plot the fits obtained by Berlind et al. (2003) to the semi-anaytic model of Cole et al. (2000) and to a smooth-particle hydrodynamics calculation of galaxy formation as solid and dotted lines respectively. The two models considered by Berlind et al. (2003) both gave reasonable fits to the galaxy correlation function.

Our new models give rise to a mean halo occupancy which is qualitatively similar to those found by Berlind et al. (2003) – a power-law at high masses, which flattens before cutting off sharply at some critical mass. Quantitatively, the models with conduction (7.3 and 7.4) produce significantly fewer galaxies in low and high mass halos, while agreeing with the Berlind et al. (2003) fits at a halo mass of a few times $`10^{12}h^1M_{}`$. The models with superwinds (8.1 and 8.3) agree more closely with the Berlind et al. (2003) fits, although model 8.3 produces somewhat fewer galaxies in clusters. Since none of our models differ greatly from those of Berlind et al. (2003), we may expect them also to produce correlation functions in reasonable agreement with the observational data, at least on intermediate and large scales. On smaller scales, the mean number of galaxies per halo is not sufficient to specify the correlation function; this requires knowledge of the distribution of galaxy pairs per halo.

Figure 3.8.2 shows correlation functions (computed using the techniques of Benson et al. 2000) for the same four models, compared to both the predicted dark matter correlation function and the observed galaxy correlation function. As expected, all models do well at matching the data on scales larger than about $`1h^1`$Mpc. On smaller scales, all models are ‘antibiassed’ with respect to the dark matter, but they produce somewhat too much clustering. Benson et al. (2000) found similar results for their semi-analytic model. Model 8.3 gives, in fact, very similar clustering results to theirs (for the same $`\mathrm{\Lambda }`$CDM cosmology). The other models produce even stronger correlations on small scales, providing a worse match to the data. Model 8.3 does better than Model 8.2 because it places relatively fewer galaxies in clusters, thereby having many fewer pairs on small scales. We speculate that a model which performs even better at matching the observed luminosity function (e.g. producing a sharper cut-off at bright magnitudes) may also do even better at matching the observed correlation function.

## 4 Discussion

We have examined the key physics thought to be necessary to explain the shape of the galaxy luminosity function in a cold dark matter universe. While the cooling and condensation of gas in a merging hierarchy of dark matter halos remains the fundamental process through which galaxies form, we have demonstrated that at least two other processes must act to shape the luminosity function. A *minimal* model requires the inclusion of feedback mechanism(s) (beyond the heating resulting from the photoionization of the pregalactic gas) to flatten the faint end of the luminosity function and to suppress cooling at the centres of the massive halos of groups and clusters. If a fraction of the energy liberated by supernovae and winds goes into reheating disk gas and/or heating the diffuse gas halo, then the formation of faint galaxies is suppressed, resulting in a flattened faint-end slope that matches the available observational data adequately. These mechanisms on their own, however, are unable to produce a sharp cut-off at the bright end of the luminosity function. We have shown that there are two possible (but quite extreme) processes that can achieve this: (1) thermal conduction at about or above the Spitzer rate; (2) expulsion of gas from halos in superwinds at temperatures high enough to prevent its subsequent recapture.

The high value of $`\alpha _{\mathrm{cond}}`$ required to suppress the bright end of the luminosity function is discouraging since it implies both that the conductivity must be close to the Spitzer value (despite the presence of $`\mu G`$ magnetic fields in clusters; Narayan & Medvedev 2001; Taylor et al. 2002), and that the effective temperature gradient must be somewhat steeper than the gradients observed in galaxy clusters (Allen et al., 2001; Fabian, Voigt & Morris, 2002).

However, our method for calculating the effect of conduction is highly simplified. In particular, there are two issues that are not well addressed. Firstly, our calculation is based on the heat flux through a shell. We have not considered the total extent of the region in which conduction must be effective. We can define a radius $`r_{\mathrm{out}}`$, such that the initial thermal energy in the region $`r_{\mathrm{cool}}<r<r_{\mathrm{out}}`$ is equal to the total energy radiated from the region $`r<r_{\mathrm{cool}}`$, where $`r_{\mathrm{cool}}`$ is the cooling radius in the presence of conduction. Thermal energy needs to be conducted from a radius at least as large as $`r_{\mathrm{out}}`$ in order to balance the radiative losses from smaller radii. For an isothermal halo with an $`r^2`$ density profile, we find $`r_{\mathrm{cool}}r_{\mathrm{out}}=r_{\mathrm{cool}}^2`$, where $`r_{\mathrm{cool}}^{}`$ is the cooling radius calculated in the absence of conduction. The radius $`r_{\mathrm{out}}`$ can be large, but the temperature gradient that is required to transport heat effectively across this region is a strongly declining function of radius. Thus, so long as $`r_{\mathrm{out}}`$ lies within the virialized region of the cluster, this is unlikely to modify our solution significantly. Secondly, our method assumes that the logarithmic temperature gradient within $`r_{\mathrm{cool}}`$ is constant. If, instead, the gradient is peaked in this region, our method would underestimate the heat deposited in the shell at $`r_{\mathrm{cool}}`$. These two simplifications (which affect the required conduction coefficient in opposite senses) are difficult to model accurately and without ad-hoc assumptions. Clearly, this is a problem that needs to be addressed by numerical simulations of cooling in conductive plasmas. These would allow us to calibrate our simple model of conductive heating and to infer what conductivity is required, by removing the degeneracy between the shape of the halo temperature profile and the suppression factor $`f_{\mathrm{sp}}`$ for the conductivity.

The alternative model invokes highly energetic “superwind” outflows. The expulsion of gas from the halo reduces the reservoir of baryons available for cooling at later times, and reduces the number of stars formed by each parcel of cold gas. However, the fraction of baryons required to be in this “super-heated” phase at $`z=0`$ is quite high, and the energy budget is exceptional. We found the energy budget of a promising superwind model to be in region of $`5\times 10^{49}`$ergs per Solar mass of stars formed. It is quite possible our model may not treat the recapture process adequately and that we may therefore be overestimating the total energy requirement. It is difficult to model the process better, however, because the expelled gas may cool adiabatically and then be accreted by a completely different halo; alternatively, some material may escape into void regions and evade recapture altogether. A further problem is the need to assume a particular distribution for the wind particle energies. Our assumption that the distribution is thermal may lead us to overestimate the energy requirement.

On average, the mass of baryons ejected from halos amounts to 7% of mean cosmic value. This corresponds to about 1.3 times the fraction of baryons turned into stars. A compilation of observations by Martin (1999) suggests that the total mass loss rate from disk galaxies could certainly be as high as twice the star formation rate. If all of this mass loss were energetic enough to produce a superwind, these observations would provide some justification for the high value of $`ϵ_{\mathrm{sw}}`$ required in our model. The total energy budget, however, exceeds that available from supernovae ($`0.9\times 10^{49}\text{erg}M_{}^1`$ for the Kennicutt IMF) even though we have not allowed for radiative losses. Although the exact energy output from supernovae is uncertain, and the conditions we have assumed for superwind escape are also uncertain, this suggests that the supernova energy must need to be supplemented by energy released during black hole formation (Enßlin et al., 1998; Bower et al., 2001; Cavaliere et al., 2002).

Finally, we have briefly examined the ability of some of our most successful models to match the observed Tully-Fisher relation and galaxy correlation function. Without adjusting any parameters, we find, as in many previous semi-analytic calculations, that our models adequately match the slope and scatter of the Tully-Fisher relation, but tend to overpredict the circular velocity at a given luminosity (or, equivalently, underpredict the luminosity at given circular velocity.) In most cases, the agreement is at least as good as in previous studies, but for one of our superwind models the discrepancy in the circular velocity zero-point is reduced to only 10%. At least one of our models (Model 8.3) matches the galaxy autocorrelation function reasonably well (as well as the semi-analytic model of Benson et al. (2000)). Importantly, this success reflects the fact that clusters in this model contain relatively fewer galaxies, and thus many fewer pairs on small scales, than clusters in other models. This demonstrates how measurements of galaxy clustering place important constraints on halo mass-dependent feedback processes such as those considered in this study. Our results encourage a more detailed investigation of how our models fare when compared with a wider range of observables, a programme that we intend to pursue in future work.

## 5 Conclusions

One of the simplest properties of the galaxy population, the luminosity function, is now reasonably well determined observationally (although, as we have noted, some significant disagreements about the faint-end slope remain). We have shown that our theoretical understanding of this fundamental property has progressed enormously in recent years to the extent that it is now possible to explain the form of the luminosity function in quantitative detail. Nevertheless, two crucial ingredients—feedback and conduction—are understood only poorly at best. Further study of these physical processes is required before we can truly claim a satisfactory understanding of the luminosity function of galaxies.

## Acknowledgments

We thank J.-S. Huang and collaborators for providing their K-band luminosity function data in electronic form and Simon White and Peder Norberg for helpful conversations. AJB acknowledges the hospitality of the Institute of Computational Cosmology of the University of Durham where much of this work was completed. RGB acknowledges support from a Leverhulme Fellowship. This work was supported in part by the UK Particle Physics and Astronomy Research Council, the Royal Society and the EC Research Network for Research into the Physics of the Intergalactic Medium.

## Appendix A Ejection and Recapture of Superwinds

The wind is ejected at speed $`V_{\mathrm{SW}}`$, so has mean energy per unit mass of $`E_{\mathrm{av}}=\frac{1}{2}V_{\mathrm{SW}}^2`$. We assume that to escape the halo it needs energy per unit mass $`\lambda _{\mathrm{SW}}V_{\mathrm{disk}}^2`$. Assuming a thermal distribution of energies in the wind (i.e. $`f(E)\mathrm{exp}(E/E_{\mathrm{av}})`$), the fraction of the mass flowing out of the disk which also escapes from the halo is given by

$$f_{\mathrm{eject}}=_{x_{\mathrm{SW}}}^{\mathrm{}}\mathrm{exp}(x)dx=\mathrm{exp}(x_{\mathrm{SW}}),$$
(A1)
where $`x=E/E_{\mathrm{av}}`$ and $`x_{\mathrm{SW}}=\lambda _{\mathrm{SW}}V_{\mathrm{disk}}^2/E_{\mathrm{av}}`$. For the material which escapes the halo, the mean energy excess over that needed to eject it from the halo is

$$E_{\mathrm{esc}}=E_{\mathrm{av}}_{x_{\mathrm{SW}}}^{\mathrm{}}(xx_{\mathrm{SW}})\mathrm{exp}(x)dx/_{x_{\mathrm{SW}}}^{\mathrm{}}\mathrm{exp}(x)dx=E_{\mathrm{av}},$$
(A2)
i.e. the mean (kinetic+thermal) energy per unit mass for the gas escaping from the halo is the same as the mean (kinetic+thermal) energy per unit mass for the gas flowing out of the disk (since only the highest energy particles escape). Using therefore the same energy distribution, we assume that the fraction of particles recaptured in a halo of virial velocity $`V_{\mathrm{halo}}`$ is

$$f_{\mathrm{cap}}=_0^{x_{\mathrm{cap}}}\mathrm{exp}(x)dx=1\mathrm{exp}(x_{\mathrm{cap}})$$
(A3)
where $`x_{\mathrm{cap}}=C_{\mathrm{cap}}V_{\mathrm{halo}}^2/E_{\mathrm{av}}`$, and $`C_{\mathrm{cap}}`$ is a numerical coefficient of order unity. We have used $`C_{\mathrm{cap}}=1`$.


format: text

text:
# Non-classical nucleation in supercooled nickel

## I Introduction

Solidification of liquids usually begins by the formation of small clusters or nuclei. Understanding the structure of nuclei and the dynamics of nucleation and growth is of great scientific and technical interest. Nucleation is a highly non-linear process and the interplay between the non-linearity and a change in symmetry is not understood in, for example, the crystallization of super-cooled liquids. In addition to the interesting scientific questions such as the structure of the critical or nucleating droplet and the nature of the initial growth, the form of the nucleation process can have a profound effect on the morphology of the final state. Classical nucleation theoriesgunt  predict that the interior of the critical droplet will have the same structure as the stable phase. However Klein and Leyvraz(KL)kle86 ; kle01  have argued that the nucleating droplet in supercooled liquids with long range potentials is non-classical and can have a structure which is significantly different than that of the stable phase due to the presence of a spinodal. They suggest that metals with effective long range interactions due to elastic forces would be a candidate for exhibiting this process which takes place arbitrarily close to a spinodal. The non-classical droplet is predicted to have a density close to that of the supercooled liquid and a triangular structure in 2 dimensions. In 3 dimensions either randomly stacked triangular planes or a bcc structure is predicted depending on the details of the potential. In the initial growth phase of the droplet it is predicted that a core will develop that has the density, and symmetry, of the stable solid phase, or of a metastable solid phase.kle01 ; uk85

Simple models of elastic systems exhibiting an Austenite to Martensite transition have indicated a non-classical nucleation pattern.kle02  However, there has been no test of the prediction of non-classical nucleation in either experiment or simulations with realistic metal potentials. Here we present a computational study examining the nucleation and growth in metastable liquid nickel using state-of-the-art semi-empirical potentials. Our results indicate that the nucleation is non-classical and is consistent with nucleation near a spinodal. It is important to stress that by nucleating or critical droplet we refer to the saddle point structure. Physically this is the droplet that just reaches the critical size. This appears in the theory as a saddle point.kle01 ; ch ; lang

## II Method

To perform our simulations we use the modified embedded atom method (MEAM) potential, referred to as potential 4 in Ref. bas97 . We extended the radial cutoff, rcut, from 4.0 to 4.8 Åto conserve energy in the liquid phase. This alteration does not affect the solid properties. A thorough discussion of the MEAM formalism can be found in the literature bas97 ; bas92 ; bas94  and will not be repeated here. Recently, this particular potential was shown to describe the properties of liquid nickel quite well che02 . The calculated melting point was within 9 % of experiment. The calculated viscosity of 3.84 mPa-s at 1775 K was in agreement with the lower of the experimental range of values $``$4.0-6 mPa-s.

We perform a molecular dynamics quenching sequence using an isobaric-isothermal (NPT-ensemble). The MD calculations use a time step of 1 fs. Constant temperature is maintained by using a Nosé-Hoover nos91  thermostat. The Parinello-Rahman method for pressure control is used maintaining the pressure to within $`\pm `$100 bar. The sequence is begun by equilibrating a 3D periodic cell of (e.g.) 1500 atoms in the liquid state at a temperature slightly above the melting point of the system for 25 ps. The samples are cooled at 10 K increments and held at that temperature for the length of time corresponding to the desired cooling rate. For example, a 1012 K/s cooling rate would be obtained by holding the sample at each temperature step for 10 ps. The quench rates for this study ranged from 1011 K/s to 1013 K/s. A few calculations are performed using a larger cell of 2048 and 12000 atoms.

Averages of the potential energy and volume are determined at each temperature after 1 ps. The only exception to this is for the quench rate of 1013 K/s, where the average temperature and volume are determined for the last 0.5 ps. From this data, we determine whether the specimen is undergoing a liquid-to-crystal or a liquid-to-amorphous transition. The crystallization produces an inflection on the cooling curve and the inflection point determines the solidification temperature, T$`{}_{}{}^{^{}}{}_{m}{}^{}`$. The liquid-to-amorphous transition produces a change of slope in the cooling curve. The glass transition temperature, Tg, is determined from the intersection of straight lines fit to the cooling curve above and below Tg.

We use common neighbor analysis (CNA) to identify the local crystalline environment of each atom fak94 ; bar84 ; hon87 ; jon88 . This method provides an efficient way to identify the local environment of each atom in the system. The results are quite insensitive to small displacements of the atoms. Since our CNA is not configured to identify planar structures we also perform a visual analysis of possible planar structures both surrounding the 3D Bravais lattice identified by the CNA and free standing(i.e. without the Bravais core).

## III Results and Discussion

In Fig. 1, we present the average atomic energies and atomic volumes as a function of temperature for various quench rates. For quench rates below 5 x 1012 K/s a rapid decrease in energy and volume is seen in the region of temperature between 700-900 K. The equilibrium freezing temperature is 1540Kche02 . The decrease occurs over a smaller range of temperature as the quench rate is lowered. In contrast, for quench rates of 5 x 1012 K/s and higher (not shown) only a subtle change of slope in the energy or volume is seen. We attribute the rapid decrease in energy and volume to crystallization and the change of slope to solidification into an amorphous state. Below the transition, the curves for both energy and volume are parallel, but offset from each other. Thus the specific heat and thermal expansion of the crystallized material are independent of quench rate (which determines local structure) but the actual energy and volume differ primarily due to quenched defects (vacancies) in the solid. Above the transition region note that the curves all lie on top of each other. The temperature derivative of the liquid energy (specific heat) and volume (thermal expansion) are independent of quench rate. Hence the local structure in the liquid is not affected by the quench rate.

To determine whether there are any important finite size effects we compare the results in systems containing 12000, 2048 and 1500 atoms. The resulting energy curves are shown in Fig. 2. The point where the slope begins to increase as the temperature decreases indicates the onset of crystallization. Nucleation occurs at a slightly higher temperature for a given quench rate as the system size increases. This implies that the system size dependence for nucleation times is as expectedgunt  and finite size effects are not evident.

Analysis of the resulting microstructure showed that for some conditions the small systems are polycrystalline and others were single crystals. Similarly the large systems exhibit both polycrystalline and single crystalline states. The distribution of the polycrystalline domains is similar in all sizes examined. In addition the CNA as well as visual inspection(to be discussed below) show similar patterns of nucleation. We thus conclude that system size plays no important role in the resulting microstructure.

Using CNA and visual inspection we analyze the nucleation process in supercooled liquid Ni. Upon cooling for both the large system (12,000 atoms) and the small system (1500 atoms), a cluster of atoms having the fcc structure and the hcp structure appear to nucleate and grow at approximately $`t_o=1.36`$ nanoseconds into the quench. We use appear since we are looking for the saddle point object calculated in the theory.kle01  After some time into the growth of this cluster the fcc particles dominate the structure. One implication of this result is that the argument of Alexander and McTagueale78 , does not apply to this particular systemkle01 . This argument states that in systems where the densities of the liquid and solid are similar the solid phase has a bcc structure. In addition the prediction of icosahedral structures in supercooled liquids as suggested by Frankfra52  (Lennard-Jones), Jónsson et al.jon88  (Lennard-Jones, LJ), and Ramprasad et al.ram93  (Embedded Atom Method, EAM) is not observed. We attribute these differences to the spherical nature of the LJ and EAM potentials, which do not include angular dependence included in the MEAM potential. This characteristic appears to prevent the formation of icosahedral ordering in pure nickel.

As discussed above, the KL theory kle86 ; kle01  predicts stacked triangular planes or bcc for the droplet structure. The KL theory also predicts the development of a classical fcc/hcp core during the growth phase. Therefore we look visually for a halo surrounding the fcc/hcp structure described above. The possibility of a bcc halo was excluded by the CNA analysis so we search for stacked planes with a triangular structure. These were found as illustrated in Fig. 3. We then visually examine the system at 18 ps prior to the configuration at $`t_o`$ to investigate the possibility that there is a droplet consisting of only stacked triangular planes with no fcc/hcp core. This configuration is also found and is illustrated in Fig. 4. These droplets contain $`60`$ atoms. Finally, we investigate which of these structures; stacked triangular planes with, or without, the core correspond to the saddle point droplet predicted by the theory.kle86 ; kle01  This is done via an intervention technique. We take the system at a time which we consider it possible that a nucleation droplet exists in the sense of a saddle point solution. We then make $`N`$ copies of the system, scramble the velocities keeping the temperature fixed, and then restart each copy. In the $`N\mathrm{}`$ limit a saddle point implies that half of the copies will continue to grow and half will disappear.mmr ; mon  If the droplet is in the growth phase then $`N^{}>>N/2`$ droplets will grow. If the copies are made prior to the saddle point droplet then $`N^{\prime \prime }>>N/2`$ droplets will disappear. Using the configuration at $`t_o=1.36`$ nanoseconds we find that in all of the copies the droplet grows. We also find that the same result is obtained for configurations chosen at 1.5 and 11.5 ps after $`t_o`$. Copies made at 28.5 ps prior to $`t_o`$ all contain droplets that disappear. However, for configurations chosen at 8.5 and 18.5 ps prior to $`t_o`$ half of the droplets in the copies disappear and half grow. Since we used $`N=10`$ some spread in the location of the saddle point droplet was expected. In addition, the saddle point flattens at the top near a spinodal.uk85  Even though there was some spread in the location of the saddle point droplet it is important to note that the droplets at these two times consist only of randomly stacked planes with an in plane triangular structure and do not contain an fcc/hcp core.

In conclusion, our study shows that the local structure of the initial (saddle point) nucleation droplet consists of randomly stacked planes with a triangular ordering. During the initial growth phase a dense core develops which has both fcc and hcp atoms. Further into the growth phase the core becomes dominated by fcc atoms. This nucleation sequence is in direct contradiction of the argument of Alexander and McTagueale78 ; kle01  and appears to exhibit the non-classical behavior that Klein and Leyvrazkle86 ; kle01  predicted. This result has important implications for the calculation of nucleation rates and the understanding of the morphology in either the stable or metastable state that has been formed via nucleation. It also gives additional evidence that the spinodal, or limit of metastability of super-cooled liquids may play an important role in determining the local structure of a critical droplet. Finally, non-classical droplets have been found in Lennard-Jones fluidsyang88 ; yang90  and colloids,gasser01  however this is the first evidence for non-classical nucleation influenced by a spinodal in a model of a molecular material such as a metal.