text:
# The Role of Heating and Enrichment in Galaxy Formation

## 1. Introduction

It has long been recognized that the X-ray luminosity-temperature ($`L_xT`$) relation of clusters does not obey the simple scaling laws that would hold if clusters were formed from the collapse of unheated primordial gas, and thus the gas within clusters is likely to have been heated before the formation of the clusters themselves (Kaiser 1991). Subsequent investigations have determined that preheating is also necessary to explain the $`L_xT`$ relation as a function of cluster mass (Cavaliere, Menci, & Tozzi 1999).

The high metallicity of cluster gas and the claimed over-abundance of alpha elements (Gibson and Matteucci 1996; Lowenstein & Mushotzky 1996) points to preheating by SNII driven winds (Renzini et al. 1993; Trentham 1994; Nath & Chiba 1995). High enrichment is not restricted to the most massive clusters, but appears to be widespread, extending to groups of galaxies (Buote 2000). Empirically the epoch of this enrichment is now known to occur at $`z>0.5`$, and approximately solar enriched cluster gas has been detected at redshifts as large as $`z1`$ (Hattori et al. 1997). Hence it is natural to suppose that pre-enrichment and preheating are the consequence of very early and vigorous massive star-formation.

Cluster gas is evidently a sink for enriched and heated material, and the most likely culprit for this enrichment is dwarf galaxies. Theoretical work has shown that supernovae and OB winds in these low-mass objects should lead to the production of energetic outflows with temperatures on the order of $`10^6`$ K (Larson 1974; Dekel & Silk 1986; Vader 1986). This behavior has been clearly identified in studies of both local starbursting galaxies (Axon & Taylor 1978; Marlowe et al. 1995; Heckman 1997; Hunter et al. 1998; Martin 1998) and spectroscopy of high-$`z`$ galaxies ( Franx et al. 1997; Pettini et al. 1998; Frye & Broadhurst 1998; Warren et al. 1998). Whether these outflows lead to a catastrophic loss of the interstellar gas however, is likely to depend on a number of factors (De Young & Heckman 1994), and is a subject of current investigations (Murakami & Babul 1999; Mac Low & Ferrara 1999; Strickland & Stevens 1999).

In the generally investigated hierarchical models of structure formation such as the Cold Dark Matter model (CDM), the existence of an era of widespread enrichment by outflows from dwarf galaxies is in fact quite natural, because low-mass galaxies are expected to form in large numbers and at early times (e.g., White & Frenk 1991). Since these early galaxies will be found preferentially in the large-scale overdense regions that later form clusters, dwarf outflows are the obvious candidates for pre-cluster heating.

The existence of large numbers of small galaxies at high redshifts is also favored by observations, to help understand the steep number counts and low luminosities of faint galaxies (Broadhurst, Ellis, & Glazebrook 1992), the sizes of faint galaxies in Hubble Deep Field images (Bouwens, Broadhurst, & Silk 1998a,b), and the small sizes of the distant Lyman-break galaxies (Steidel et al. 1999). Locally, however, the space density of dwarf galaxies relative to massive galaxies is far less than predicted on the basis of the steep Press-Schechter slope for the faint end of the mass function (Ferguson & Binggeli 1994), prompting theoretical studies of the disruption of dwarves by tidal forces from neighboring objects (Moore et al. 1998), external UV radiation (Kepner, Babul, & Spergel 1997; Norman & Spaans 1997; Corbelli, Galli, & Palla 1997), and catastrophic mass loss during outflows (Larson 1974; Dekel & Silk 1986; Vader 1986).

While many of these mechanisms may have had an impact on the formation of dwarf galaxies, relatively little attention has been directed towards the influence of outflows on neighboring galaxies. As the earliest galaxies to form were highly clustered (Kaiser 1984) and typical outflow temperatures and velocities were much larger than the virial temperatures and velocities of these galaxies, it is likely that dwarf galaxies were strongly influenced by their neighbors (Scannapieco, Ferrara, & Broadhurst 2000). Similarly, the sources responsible for pre-enrichment of the intracluster medium may well have enriched larger protogalaxies, with important consequences for their metallicity histories and cooling times.

In this work, we conduct an idealized investigation of the impact of dwarf outflows on the history of galaxy formation. While the consequences of homogeneous heating on galaxy formation have been examined in the past (Blanchard, Valls-Gabaud, & Mamon 1992), the inhomogeneous nature of this process and the associated pre-enrichment of galaxies have not been addressed. Our aim is not, however, to construct a complete model of galaxy formation, and thus we do not track processes such as the formation of second-generation stars in galaxies, production of dust, transfer of angular momentum, or the structure of the interstellar medium. Rather we focus on the properties of the inhomogeneously heated and enriched intergalactic medium (IGM) out of which galaxies coalesced, the likely consequences of galaxy formation in this environment, and to what degree these issues must be accounted for within more detailed simulations of galaxy formation.

The structure of this work is as follows. In §2 we review the observational evidence that leads us to consider a model in which widespread dwarf outflows shocked and enriched the medium out of which larger galaxies formed. In §3 we describe a simple numerical Monte Carlo code that we use to asses the overall features of such a model. The results of our simulations are given in §4 in which we examine which aspects of galaxy formation are most sensitive to the presence of outflows. In §5 we discuss the limitations of our modeling, and conclusions are listed in §6.

## 2. Observational Evidence

In this section we review the observational evidence that points towards widespread early enrichment by dwarf galaxies. We show that the chemical and thermal properties of nearby galaxy clusters are both suggestive of such pre-enrichment, as is the weak evolution of these properties with redshift. We show that such a picture is consistent with our knowledge of the intergalactic medium as well, and helps to explain the relatively high metallicity seen in many Ly absorption systems. The presence of outflows has been identified in emission line studies of high-$`z`$ galaxies, and several nearby dwarf galaxies have been caught “red-handed”, surrounded by clouds of enriched gas with temperatures greatly exceeding virial. Finally, we review some unsolved questions in galaxy formation that may depend on preheating and enrichment and serve to motivate our exploratory numerical studies.

### 2.1. Galactic Outflows and the Intra Cluster Medium

Many cluster properties can be understood in the context of self-similar models (Kaiser 1986, 1990) which are exact for power-law initial fluctuation spectra and a reasonable approximation for more realistic spectra such as that of the Cold Dark Matter (CDM) model. While the optical properties of galaxies in clusters are consistent with these predictions, the slope of the X-ray luminosity-temperature relationship of the hot intra-cluster gas is too steep and slowly evolving to be understood in this context.

A convincing explanation of $`L_xT`$ evolution has been made by Kaiser (1991) who speculated that at an epoch before the formation of present-day clusters, the gas which formed the intracluster medium (ICM) was preheated, injecting sufficient energy to expel gas from the small potential wells that were non-linear at early times. By assuming that this heated gas later cooled adiabatically onto the large cluster-size potential wells, he was able to simply reproduce the observed X-ray properties of clusters. Numerous investigations of preheating and $`L_xT`$ evolution have also reached similar conclusions (David, Jones, & Forman, 1996; Mushotzky & Scharf 1997; Eke, Navarro, Frenk 1998) as have studies of the $`L_xT`$ relation in groups of galaxies (Cavaliere, Menci, & Tozzi 1999).

The origin of this preheating is most naturally due to Supernovae Type-II (SNeII) activity in the dwarf starburst population. Several studies of expanding HI gas in nearby dwarf galaxies show clear evidence of dense expanding shells with velocities above $`15`$ km/s (Marlowe et al. 1995; Heckman 1997; Hunter et al. 1998; Martin 1998). $`\mathrm{Ly}\alpha `$ studies confirm that similar outflows are present in higher $`z`$ galaxies at redshifts of order $`3`$ (Pettini et al. 1998) and even higher (Frye & Broadhurst 1998; Warren et al. 1998 ), where the overall space density of small starbursting galaxies is larger.

A natural prediction of such a picture would not only be the heating of the ICM by the dwarf population, but the enrichment of this gas by metals expelled by the SNeII powering the galactic winds. In fact, the metallicity of the intracluster medium is observed to be quite high and constant ($`0.3Z_{}`$) over a large range of cluster masses (Renzini 1997). This roughly constant value is evidence that the clusters have neither lost or gained a large amount of baryons during their evolution, as massive cluster outflows or accretion of pristine gas would lead to a large scatter. Furthermore, this metallicity is roughly the same today as observed at $`z0.3`$ and perhaps even at $`z1`$ (Hattori et al. 1997), suggesting that IGM enrichment took place early on in the lifetime of the clusters (Renzini 1999).

Several authors have investigated the possibility that this enrichment is due to the present cluster members and that most of the X-ray emitting hot gas in clusters is due to outflows from galaxies that populate them today. Okazaki et al. (1993) considered gas ejection from bright, elliptical galaxies in clusters, concluding that these galaxies could contribute no more than 10% of the ICM. Trentham (1994) suggested that the ICM instead may be primarily due outflows from precursors to the observed populations of dwarf galaxies in clusters, but Nath and Chiba (1995) concluded that the metallicities generated in this type of scenario were sufficient to explain only clusters with low-metallicity gas.

Perhaps most convincingly Gibson and Matteucci (1996) have studied the possibility that the elliptical galaxies observed in clusters today were responsible for the majority of ICM enrichment. By developing models of galactic winds consistent with the observed properties of cluster ellipticals they concluded that even in their “maximal models” in which all the gas returned by dying stars is ejected into the ICM, neither the giant elliptical, the dwarf spheroidals, or both these populations combined can be responsible for more that 33/38% of the ICM. And yet this gas is hot and metal rich.

A final piece of evidence as to the origin of ICM metals comes from studies of relative abundances. While for many years, studies of metals in clusters were limited to measurements of iron abundances, the X-ray Advanced Satellite for Cosmology and Astrophysics (ASCA) (Tanaka, Inoue,& Holt 1994) provided the first opportunity to study other metals. This allowed Loewenstein & Mushotzky (1996) to observe excesses of alpha-elements such as O, Mg, and Si, which indicate that most of the ICM enrichment was due to SNeII rather than Type Ia (SNeIa). This result was later contested by Ishimaru & Arimoto (1997), who claimed that SNeIa enrichment could be responsible for 50% or more of ICM metals, based on the solar ‘meteoritic’ metallicity rather than the solar coronal gas. Gibson, Loewenstein, & Mushotzky (1997) pointed out, however, that this result was specifically linked to the SNeII yield used in their analysis. While they were unable to rule out Ishimaru & Arimoto’s claim within theoretical uncertainties, they showed that more recent SN models that treat convection and mass-loss, reduce the fractional contribution ICM by SNeIa to less than 5%. Renzini (1997) on the other hand while confirming the importance of early SneII, has argued that later enrichment by significant numbers of SNeIa are required to produce the observed high level of Fe relative to alpha elements.

In summary then, metallicities studies of clusters indicate that the ICM was enriched early on, that this enrichment is too high to be due to purely to the observed cluster members, and that the relative abundances of metals are suggestive of widespread SNeII enrichment; all indicative of enrichment by high-redshift dwarf galaxies.

### 2.2. Galactic Outflows and the Intergalactic Medium

The first galaxies to form would have necessarily had a huge impact not only in regions that would later form clusters, but on the Intergalactic Medium (IGM) as a whole. In standard schemes for hierarchical structure formation, the first baryonic objects formed at redshifts $`40`$ and with masses of order $`10^5M_{}`$, somewhat smaller than the dwarf galaxies observed today (Haiman, Thoul, & Loeb 1996). While several authors have shown that the total flux of UV photons from stars within these objects is sufficient to reionize the universe (Couchman & Rees 1986; Fukugita & Kawakasi 1994; Shapiro, Giroux, & Babul 1994), these objects are likely to suppress their own formation long before this occurs. As molecular hydrogen is easily photo-dissociated by 11.2-13.6 eV photons, to which the universe is otherwise transparent, the emission from the first stars quickly destroys all avenues for cooling by molecular line emission. This quickly raises the minimum virial temperature necessary to cool effectively to approximately 10,000 K, suppressing the further formation of objects with masses $`10^8M_{}`$ (Haiman Rees, & Loeb 1997; Ciardi et al. 1999).

Thus reionization, if achieved by galaxies, was relatively late ($`z15`$) and associated with objects of similar size or larger than dwarf galaxies. This high mass scale, implies that the IGM phase transition must have occurred in a clumpy and inhomogeneous manner, with important implications for small-scale microwave background anisotropies (see eg. Aghanim et al. 1996; Miralda-Escude, Haehnelt, & Rees 2000; Scannapieco 2000).

Early enrichment and heating is similarly expected to be inhomogeneous, following the spatial distribution of the first galaxies which are restricted to rare overdense regions. Hence measurements of a low average metallicity ($`0.01Z_{}`$) of the Ly-forest at high redshift (Songalia & Cowie 1996; Savaglio 1997; Songalia 1997) may be viewed not so much as evidence for a low metallicity intergalactic medium, as an inhomogeneous one. High-redshift regions of higher density that later develop into clusters would therefore represent the most polluted volumes, where one would naturally expect a significantly higher mean metallicity.

Direct evidence that the IGM is enriched by galaxy outflows may be inferred from detailed observations of local dwarf spheroidal galaxies. For example, ASCA observations (della Ceca et al. 1996) show that the star-forming dwarf galaxy NGC 1569 is surrounded by a hot ($`8\times 10^6`$ K) halo of gas whose temperature greatly exceeds the virial temperature of the galaxy and whose line strengths are consistent with $`0.25Z_{}`$. Similarly, X-ray observations of the dwarf irregular NGC 4449 by Bomans, Chu, & Hopp (1997) show that it is embedded in a supergiant shell of $`2\times 10^6`$K gas with a metallicity of $`0.3Z_{}.`$

### 2.3. Galactic Outflows and Galaxy Formation

The presence of an inhomogeneous IGM with strong temperature and metallicity fluctuations creates quite a different environment for galaxy formation than the primordial conditions often assumed. This becomes clear when one contrasts the virialization and cooling of primordial gas to that of gas pre-enriched and heated by dwarf outflows. In Figure 1 we replot a classic comparison first made in Dekel and Silk (1986). Collapse redshifts of $`1\sigma `$ to $`2\sigma `$ dark matter halos (dotted lines) are calculated from linear theory for both a flat CDM ($`\mathrm{\Omega }_0=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$, $`\mathrm{\Omega }_b=0.07`$, $`\sigma _8=0.6`$, and $`\mathrm{\Gamma }=0.44`$) and a $`\mathrm{\Lambda }`$CDM model ($`\mathrm{\Omega }_0=0.35`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$, $`\mathrm{\Omega }_b=0.06`$, $`\sigma _8=0.87`$ and $`\mathrm{\Gamma }=0.18`$), assuming spherical collapse. The solid lines show the lowest redshift that a sphere of gas can virialize and still have time to cool and form a galaxy by $`z=0`$, as calculated by the cooling models described in §3.2. In each pair the upper line corresponds to primordial gas, and the lower line to gas that has been enriched to $`0.1Z_{}`$. Here we see that even modest enrichment by outflows has the potential to greatly accelerate galaxy formation on the $`10^{12}M_{}/h`$ scale due to the long cooling times of large clouds without the additional avenue for cooling afforded by line emission by metals.

Also on this plot, we show the mass of halos associated with a virial temperature of $`5\times 10^5`$K, typical of galaxy outflows (dashed line). This serves as an estimate of the smallest mass of galaxies that can form in areas impacted by outflows without being disrupted by shocks. This suggests that while high-mass galaxy formation may be enhanced by dwarf outflows, the formation of dwarf galaxies has the potential to be suppressed.

The features shown on this plot invite comparison with the properties observed in elliptical galaxies, which are preferentially found in the most enriched regions in the universe, and whose mass functions are biased to large values relative to the rest of the galaxy population. It is also interesting to note that various analytical and N-body studies of cold dark matter models (Kauffmann, White, & Guiderdoni 1993; Klypin et al. 1999; Moore et al. 1999) have shown that $`50`$ satellites with circular velocities $`20`$ km/s should be found within 600 kpc of the Galaxy, while only 11 are observed. Again, this lack of local dwarf galaxies is suggestive of galaxy formation in an intergalactic medium that has been impacted by outflows.

Closer to home, the abundance distribution of long-lived stars in the solar neighborhood shows far too few low-metallicity stars compared with a simple “closed box” model of galactic chemical evolution. This is the long standing G-dwarf problem, first pointed out by van de Bergh (1962) and Schmidt (1963). The sudden drop in the number of G-dwarfs with metallicities below $`0.1Z_{}`$ has lead to a number of models in which an initial production spike or “prompt initial enrichment” adds metals to the gas out of which the majority of stars form (e.g., Truran & Cameron 1971; Ostriker & Thuan 1975; Köppen & Arimoto 1990). This ad hoc floor to the metallicity then allows good fits to the Galactic stellar data.

Similar pre-enrichment may also be necessary to explain the metallicity distribution of stars outside our own galaxy. Integrated spectra of elliptical galaxies have been studied with care (Worthey, Dorman, & Jones 1996) concluding rather puzzlingly that such galaxies do not have a closed-box history but again show the need for a “floor” of around $`0.1Z_{}`$ (Thomas, Greggio, & Bender 1999).

## 3. Cosmological Monte Carlo

The great preponderance of observational clues and pointers suggests that we should seriously consider a scenario in which the formation of modern-day galaxies was preceded by an era of IGM enrichment by a high-redshift population of dwarf starbursting galaxies. Note that from a theoretical point a view, such a model is not ad hoc, but rather a consequence of hierarchal structure formation. As smaller objects form early and are able to suppresses further small-scale formation, it is only natural that larger galaxies would have formed in a second wave of collapse within this pre-enriched medium.

The density dependence of the effect of outflows on galaxy formation requires a 3-D calculation. To date only the average properties of the IGM as a function of time have been explored in one dimensional calculations by Blanchard, Valls-Gabaud, & Mamon (1992) and by Nath and Trentham (1997). Such averages are appropriate to the IGM as a whole, but underestimate the effect on galaxy formation. Thus for example, the outflow model proposed by Nath and Trentham (1997) produces metallicities $`.01Z_{}`$, while the observed metallicities of cluster gas and elliptical galaxies are comparable to solar abundance (Renzini 1993).

In order to determine the impact of outflows on galaxy formation, we have developed a simple cosmological Monte Carlo code, in which we realize a small volume of the universe and study the linear evolution of objects in the mass range $`2\times 10^8\mathrm{\Omega }_0/h\mathrm{M}_{}M1\times 10^{13}\mathrm{\Omega }_0/h\mathrm{M}_{}`$ within it, with reasonable extensions to deal with collapse. Our philosophy is not to pretend that we can possibly reproduce in detail the complicated processes of dark matter halo collapse, galaxy formation, and gas infall and expulsion, but rather to explore a model that captures the essential features and predicts with some confidence rough magnitudes and trends that can be compared with observations.

### 3.1. Collapse of Dark Matter Halos

We choose for our simulation a cubic comoving volume of $`(12\mathrm{M}\mathrm{p}\mathrm{c}/h^1)^3`$ divided into $`512^3`$ cells and with periodic boundary conditions. On this mesh we construct a cosmological linear overdensity field $`\delta (𝐱,z)=\rho (𝐱,z)/\overline{\rho }(z)`$ where $`\overline{\rho }(z)`$ is global average density of the universe, $`𝐱`$ is a position in comoving coordinates, and $`z`$ is the redshift. Transforming into Fourier space the density field can be reexpressed as $`\stackrel{~}{\delta }(𝐤,z)d^3𝐱\mathrm{exp}(i𝐤𝐱)\delta (𝐱,z)`$. From linear theory, the evolution of the density field as a function of time is given simply by $`\stackrel{~}{\delta }(𝐤,z)=\stackrel{~}{\delta }_0(𝐤)D(z)/D_0`$ where $`D(z)`$ is the dimensionless growth factor, $`\stackrel{~}{\delta }_0(𝐤)\stackrel{~}{\delta }(𝐤,0)`$, and $`D_0D(0).`$

In the CDM model, the Gaussian random field $`\stackrel{~}{\delta }(𝐤)`$ can by constructed by randomly choosing modes such that the variance is given by $`\stackrel{~}{\delta }(𝐤)\stackrel{~}{\delta }(𝐤^{})=\stackrel{~}{\delta }(𝐤)\stackrel{~}{\delta }^{}(𝐤^{})=(2\pi )^3\stackrel{~}{\delta }(𝐤+𝐤^{})P(𝐤)`$. The power spectrum is then

$$P(k)=2\pi ^2\delta _H^2k^nT^2(k_p\mathrm{Mpc}/h\mathrm{\Gamma }),$$
(1)
where $`T(q)`$ is the CDM transfer function given by Eq. (G3) of Bardeen et al. (1996), and we choose $`a_0H_0=c`$, such that $`k_p=k/a_0=kH_0/c`$ is the physical wave number. The shape parameter $`\mathrm{\Gamma }`$ is defined as $`\mathrm{\Gamma }\mathrm{\Omega }_0h\mathrm{exp}(\mathrm{\Omega }_b\mathrm{\Omega }_b/\mathrm{\Omega }_0)`$ (Sugiyama 1995) and is observed be to $`\mathrm{\Gamma }=0.23_{0.034}^{+0.042}`$ given a spectral index of $`n=1`$ (Vianna & Liddle 1996). We fix the normalization factor, $`\delta _H`$, against the number abundance of clusters. The variance of the mass inclosed in a sphere of radius $`R`$ is

$$\sigma ^2(R)=\frac{1}{2\pi ^2}_0^{\mathrm{}}k^2𝑑kP(k)W^2(kR).$$
(2)
If we choose a spherical top-hat window function defined by $`W(x)3\left[\frac{\mathrm{sin}(x)}{x^3}\frac{\mathrm{cos}(x)}{x^2}\right]`$ then we can normalize our fluctuation spectrum by setting $`\sigma _8\sigma (8\mathrm{Mpc}h^1)=(0.6\pm 0.1)\mathrm{\Omega }_0^{C(\mathrm{\Omega }_0)},`$ where $`C(\mathrm{\Omega }_0)=.36+0.31\mathrm{\Omega }_00.23\mathrm{\Omega }_0^2`$ in the open case and $`C(\mathrm{\Omega }_0)=.59+0.16\mathrm{\Omega }_00.06\mathrm{\Omega }_0^2`$ if we take ($`\mathrm{\Omega }_0+\mathrm{\Omega }_\mathrm{\Lambda }=1`$) (Vianna & Liddle 1996).

Having constructed $`\stackrel{~}{\delta }_0(𝐤)`$ in this manner, we then convolve it with window functions with lengths scales corresponding to ten masses arranged logarithmically from $`M_12\times 10^8\mathrm{\Omega }_0/hM_{}`$ to $`M_{10}10^{13}\mathrm{\Omega }_0/hM_{}`$, to obtain $`\delta _0^{M_1}(𝐱)\delta _0(𝐱,R_{2\times 10^8\mathrm{\Omega }_0/hM_{}})=\frac{d^3𝐤}{(2\pi )^3}\mathrm{exp}(i𝐤𝐱)\stackrel{~}{\delta }(𝐤)W(R_{2\times 10^8\mathrm{\Omega }_0/hM_{}}k)`$, $`\delta _0^{M_2}(𝐱)\delta _0(𝐱,R_{6.6\times 10^8\mathrm{\Omega }_0/hM_{}})`$, … $`\delta _0^{M_{10}}(𝐱)`$, the linear density field at $`z=0`$ smoothed at each mass scale. Here the minimum mass of our simulations is set by molecular cooling constraints discussed further below.

From the spherical collapse model we can identify the linear overdensity $`\delta _{\mathrm{sc}}(z)`$ at which the true density field has virialized. For flat models this value is $`1.69`$ at all times while in open cases it is a weakly decreasing function of redshift as fitted in Appendix A of Kitayama & Suto (1996). For each peak in the linear overdensity field such that $`\delta _0^{M_i}(𝐱,R)>\delta _{\mathrm{sc}}(z=0)`$, we identify a collapse redshift $`z_{\mathrm{sc}}`$ such that

$$\delta _0^{M_i}(𝐱)=\delta _c(z_{\mathrm{sc}})\frac{D_0}{D(z_{\mathrm{sc}})}.$$
(3)
This approach has the limitation that it becomes inaccurate as $`\sigma (R)D(z)/D_0`$ approaches $`\delta _{\mathrm{sc}}`$. Thus while the spherical collapse model works well at determining the number and spatial clustering of rare objects (Lacy & Cole 1993; Mo, Jing, & White 1996, 1997) it fails to agree with numerical simulations for more common low-mass halos (Lacey & Cole 1994; Sheth & Tormen 1999). Furthermore, the object-by-object identification of linear peaks with collapsed objects has been shown to be unreliable, even in cases in which the statistical properties are in good agreement with N-body simulations (Bond et al. 1991).

Sheth, Mo, & Tormen (1999) have shown that these predictions can be improved by accounting for the ellipticity of collapsing clouds caused by tidal forces from nearby collapsing peaks. By including a correction factor for the critical linear overdensity, they have been able not only to improve statistical predictions over the spherical model, but do well on an object-by-object basis. As our Monte-Carlo approach depends on this object-by-object identification, we can greatly improve the accuracy of our simulations by comparing each peak with $`\delta _{\mathrm{ec}}(z)`$, such that

$`\delta _0^{M_i}=`$ $`\delta _{\mathrm{ec}}(z_{\mathrm{ec}}){\displaystyle \frac{D_0}{D(z_{\mathrm{ec}})}}`$
$`=`$ $`\delta _{\mathrm{sc}}(z_{\mathrm{ec}}){\displaystyle \frac{D_0}{D(z_{\mathrm{ec}})}}\left[1+\beta \left({\displaystyle \frac{\sigma ^2(R)}{\delta _{\mathrm{sc}}^2(z_{\mathrm{ec}})}}{\displaystyle \frac{D(z_{\mathrm{ec}})^2}{D_0^2}}\right)^\gamma \right],`$ (4)
where $`\beta =0.47`$ and $`\gamma =0.615`$.

Having identified the collapsed peaks by either algorithm, we then arrange them in order of decreasing collapse redshift, $`z_\mathrm{c}=z_{\mathrm{sc}}`$ or $`z_{\mathrm{ec}}`$, to obtain a list of candidate points and redshifts for collapsed halos at each of the different mass scales: $`𝒳_{\mathrm{candidate}}^{M_i}\{𝐱_1^{M_i},𝐱_2^{M_i},\mathrm{}\}\mathrm{and}𝒵_{\mathrm{c},\mathrm{candidate}}^{M_i}\{z_{\mathrm{c},1}^{M_i},z_{\mathrm{c},2}^{M_i},\mathrm{}\}.`$ From this list we exclude unphysical points corresponding to all halos that collapse within a previously collapsed halo at the same or greater mass scale. That is we remove all points such that

$$M_iM_j,𝐱_k^{M_i}𝐱_l^{M_j}<R_{M_j},\mathrm{and}z_{\mathrm{c},k}^{M_i}<z_{\mathrm{c},l}^{M_j},$$
(5)
where the distance is calculated accounting for the periodicity of the simulation volume.

With these points excluded, we then have the full history of the collapsing dark matter halos at each of the various mass scales: $`𝒳^{M_i}\{𝐱_1^{M_i},𝐱_2^{M_i},\mathrm{}\}`$, and $`𝒵_\mathrm{c}^{M_i}\{z_{\mathrm{c},1}^{M_i},z_{\mathrm{c},2}^{M_i},\mathrm{}\}`$. The merger history of these objects is also simply calculated by searching for all points such that

$$M_i<M_j,𝐱_k^{M_i}𝐱_l^{M_j}<R_{M_j},\mathrm{and}z_{\mathrm{c},k}^{M_i}>z_{\mathrm{c},l}^{M_j}.$$
(6)
Each such case indicates that the halo centered at $`𝐱_k^{M_i}`$ was absorbed into the larger halo centered at $`𝐱_l^{M_j}`$ at a redshift of $`z_{\mathrm{c},l}^{M_j}.`$

### 3.2. Comparison with Analytical Results

At each redshift, the number density of collapsed halos in the simulation can be compared to analytical predictions. For the spherical collapse model this gives (Press & Schechter 1974)

$`{\displaystyle \frac{dn_{\mathrm{sc}}(M,z)}{dM}}`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\rho (z)}{M}}{\displaystyle \frac{\delta _{\mathrm{sc}}(z)D_0}{\sigma (M)^2D(z)}}`$ (7)
$`\mathrm{exp}\left({\displaystyle \frac{\delta _{\mathrm{sc}}^2(z)D_0^2}{2\sigma ^2(M)D(z)^2}}\right){\displaystyle \frac{d\sigma (M)}{dM}}.`$
Strictly speaking, Eq. (7) corresponds to the number of peaks as determined by a sharp $`k`$-space filter, rather than the number determined by the top-hot filter used in our code. Equating these quantities has become somewhat of a common practice however, and has been shown to be a good estimate.

In Figure 2 we compare the number density of objects in our simulation as function of redshift with the number predicted by Eq. (7). As both the Press-Schechter integral, and our peak finding algorithm become inaccurate at late times, we exclude all peaks that collapse with a redshift lower that a minimum value, $`z_{\mathrm{min}}(M_i)`$ such that $`\sigma (M_i)D(z_{\mathrm{min}}(M_i))/D_0=1.3\delta _{\mathrm{sc}}.`$ We consider the $`\mathrm{\Lambda }`$CDM model shown in Figure 1, ($`\mathrm{\Omega }_0=0.35`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.65`$, $`\mathrm{\Omega }_b=0.06`$, $`\sigma _8=0.87`$ and $`\mathrm{\Gamma }=0.18`$, $`h=0.65`$).

Here we see that the number of identified objects agrees well at most mass scales and redshifts. At the smaller scales, there is some discrepancy in the collapse time of objects, with the code slightly over-predicting the number of objects formed at very early times. At larger-range scales, the Monte-Carlo approach does the best, matching the Press-Schechter predictions at all redshifts. Finally, at the highest mass scales, the low number of objects introduces significant statistical noise in our comparison.

The spatial distribution of collapsed halos can also be compared to analytical predictions. If we define the correlation function at a smoothing scale $`R`$ as the Fourier transform of the power spectrum smoothed by an appropriate top-hat window function, $`\xi _R(r)\frac{1}{2\pi ^2}_0^{\mathrm{}}k^2𝑑kP(k)W^2(kR)\frac{\mathrm{sin}(kr)}{kr}`$, then the correlation of collapsed peaks of scale $`R`$ can be approximated by (Kaiser 1984)

$`1+\xi _{R,\nu }(r)=`$ $`\left({\displaystyle \frac{2}{\pi }}\right)^{1/2}\left[\mathrm{erfc}(\nu /2^{1/2})\right]^2\times `$ (8)
$`{\displaystyle _\nu ^{\mathrm{}}}𝑑ye^{y^2/2}\mathrm{erfc}\left[{\displaystyle \frac{\nu y\xi _R(r)/\xi _R(0)}{\sqrt{22\xi _R^2(r)/\xi _R^2(0)}}}\right],`$
where erfc is the complimentary error function and $`\nu \delta _cD_0/\sigma (R)D(z).`$ Note that this expression arises from considering the fraction of Gaussian-distributed objects above a certain threshold, rather than from the more sophisticated excursion set approach (Bond, Efstathiou, & Kaiser 1991; Lacy & Cole 1993) used to derive Eq. (7).

In Figure 3 we compare the correlation function as given in Eq. (8) to that derived from the distances between the collapsed peaks at the five lowest mass-scales in our simulation, with redshifts chosen such that $`\nu =1.5`$ at each length scale. Eq. (8) is intended to represent the correlations between collapsed objects at or above a particular mass scale corresponding to the smoothing length scale. In practice, however, it makes little difference whether we compute $`\xi _{R_{M_i}}(r)`$ considering distances between all objects with masses above $`M_i`$ or only between those in the mass bin corresponding to $`M_i`$ itself. As comparisons only within a particular bin are both easier to calculate and understand physically, however, it is this quantity that we plot in Figure 3.

Here we see that the numerical and analytical expressions agree at lengths scales $`R^{M_i}`$, but diverge at smaller distances. As the simulation simply excludes all points closer than $`R_{M_i}`$ there is a sharp fall of in $`\xi _{R_{\mathrm{M}_\mathrm{i}}}(r)`$ below this distance. The analytical expression, however, does not exclude pairs of objects that are contained within a larger virialized object (as would be excluded in an excursion set calculation) and increases dramatically at small $`r`$. Thus these quantities are derived in a manner such that they must be discrepant at small distances, and we can have some confidence that the overall distribution of objects in our simulations is reasonable.

While the spherical collapse model is a logical testing ground for our code and can be simply compared to analytical predictions, a more accurate simulation results by applying the elliptical collapse criteria for halo formation, Eq. (4). This allows not only for a more accurate identification of collapsed peaks on an object-by-object basis, but delays the collapse of the more common peaks, and thus the minimum accurate redshift at each mass scale is decreased to $`z_{\mathrm{min}}(M_i)`$ such that $`\sigma (M_i)D(z_{\mathrm{min}}(M_i))/D_0=1.3\delta _{\mathrm{ec}}`$.

In this model, the Press Schechter relation, Eq. (7), is modified to

$`{\displaystyle \frac{dn_{\mathrm{ec}}(M,z)}{dM}}=`$ $`A\left[1+\left({\displaystyle \frac{\sigma (M)D(z)}{\delta _{\mathrm{sc}}D_0}}\right)^{2q}\right]\sqrt{{\displaystyle \frac{2}{\pi }}}{\displaystyle \frac{\rho (z)}{M}}`$ (9)
$`{\displaystyle \frac{\delta _{\mathrm{sc}}(z)}{\sigma (M)^2D(z)}}\mathrm{exp}\left({\displaystyle \frac{\delta _{\mathrm{sc}}^2(z)D_0^2}{2\sigma ^2(M)D(z)^2}}\right){\displaystyle \frac{d\sigma (M)}{dM}},`$
where $`A=0.322`$ and $`q=0.3.`$ Note that this function is essentially equivalent to the spherical collapse predictions for rare, high-mass peaks, and thus the $`\sigma _8`$ normalization from the number density of galaxy clusters is unchanged in this model. This expression has been checked against the number density of objects in our simulation and again is in good agreement.

### 3.3. Initial Gas Infall

While our simple nonlinear collapse model captures the overall history of the dark matter distribution, the evolution of the gas is more involved. While gas traces dark matter at early times, the formation of bound objects and subsequent star formation requires that the gas be able to cool to $`T0`$.

If gas collapses and virializes along with a dark matter perturbation, and we assume an isothermal distribution, then it will be heated to a temperature of (Eke, Cole, & Frenk 1996)

$$T_{\mathrm{vir}}=\frac{90\mathrm{K}}{\beta }\left(\frac{6.8}{5X+3}\right)M_6^{2/3}(1+z_c)\left(\frac{\mathrm{\Omega }_0\mathrm{\Delta }_c(z_c)}{\mathrm{\Omega }(z_c)18\pi ^2}\right)^{1/3},$$
(10)
where $`\beta `$ is the ratio of specific galaxy kinetic energy to specific gas thermal energy, $`X`$ is the hydrogen mass fraction which we take to be $`0.76`$, $`M_6M/(10^6M_{}/h)`$, and $`\mathrm{\Delta }_c`$ is the ratio of the mean halo density to the critical density at the redshift of collapse, a constant ($`18\pi ^2`$) for the $`\mathrm{\Omega }=1`$ case and a otherwise a weak function of $`z`$ fitted in Kitayama & Suto (1996). Navarro, Frenk, & White (1995) have shown that this relation with $`\beta =1.07\pm 0.05`$ is an accurate approximation to the results of N-body/hydrodynamic simulations in the case of an $`\mathrm{\Omega }=1`$ universe. Here we follow Eke, Cole & Frenk (1996) and fix $`\beta =1`$ for all cosmological models.

Once the gas has fallen into a potential well, we adopt a simple model to calculate the time scale at which the gas will cool and form stars. Neglecting the gravitational energy associated with further gas infall we have

$$n_e^2(r)\mathrm{\Lambda }(T(r),Z(r))=\frac{3}{2}n(r)k_B\frac{dT(r)}{dt},$$
(11)
where $`n_e(r)`$ is the number density of electrons at a radius $`r`$, $`n_{\mathrm{tot}}(r)`$ is the total number density, $`T(r)`$ is the temperature of the gas, $`k_B`$ is the Boltzmann constant, and $`\mathrm{\Lambda }`$ is the radiative cooling function, which is strongly dependent on the metallicity $`Z(r)`$ of the gas. Here we take $`\mathrm{\Lambda }`$ as tabulated by Sutherland and Dopita (1993) for equilibrium configurations, although strictly speaking these estimates are not exact due to the difference in abundance ratios for alpha elements in SNeII outflows as compared to solar proportions. Considering only helium and hydrogen the mean molecular weight is $`\mu m_p=\rho _g/n=m_p4/(85Y)`$ and the ratio of the electron and total number density are related by $`\eta \frac{n_e}{n_{\mathrm{tot}}}=\frac{4(Y^11)+2}{8(Y^11)+3}`$, where $`Y`$ is the helium fraction by mass, here taken to be $`0.25.`$

For a general configuration of gas, Eq. (11) is difficult to solve, and we therefore adopt a simple model for each collapsing halo as an isothermal sphere ($`\rho (r)r^2`$) with constant metallicity. We then follow the simple heuristic model of White & Frenk (1991) in which all the gas within some “cooling radius” $`r_{\mathrm{cool}}`$ cools instantaneously and all the gas outside this radius stays at the virial temperature of the halo, with $`r_{\mathrm{cool}}`$ moving outward with time. For an isothermal sphere this gives

$`{\displaystyle \frac{dM_{\mathrm{cool}}}{dt}}=`$ $`4\pi \rho _g(r_{\mathrm{cool}})r_{\mathrm{cool}}^2{\displaystyle \frac{dr_{\mathrm{cool}}}{dt}}`$
$`=`$ $`12f_{\mathrm{hot}}^{3/2}\left({\displaystyle \frac{T_{\mathrm{vir}}}{\mathrm{K}}}\right)\left({\displaystyle \frac{\mathrm{\Lambda }(T_{\mathrm{vir}},Z)}{10^{23}\mathrm{ergs}\mathrm{s}^1\mathrm{cm}^3}}{\displaystyle \frac{\mathrm{yr}}{t}}\right)^{1/2}M_{}\mathrm{yr}^1,`$ (12)
where $`f_{\mathrm{hot}}`$ is the fraction of the halo mass in the form of hot gas. While derived from a simple model, this expression is in good agreement with similarity solutions for cooling flows given by Bertschinger (1989) (which are $`28\%`$ smaller than this expression), and the one-dimensional simulations of Forcada-Miró & White (1996) (which are $`15\%`$ smaller). We therefore follow Somerville (1997) in multiplying the right hand side of this equation by an overall factor $`f_0`$ which we take to be $`0.8`$. Finally, we approximate $`f_{\mathrm{hot}}M_{\mathrm{baryon}}/M_{\mathrm{halo}}=\mathrm{\Omega }_b/\mathrm{\Omega }_0`$ at all times, and take the time at which $`M_{\mathrm{cool}}=M_{\mathrm{baryon}}`$ to be time at which a galaxy is formed.

We show this formation time as a function of mass for four different metallicities at $`z=1`$ in Figure 4. From this diagram we see that while cooling times are negligible for dwarf galaxies, cooling times for primordial-abundance halos with masses above $`10^{12}M_{}`$ are sufficiently long to prevent galaxy formation at low redshifts. This is dependent on metallicity, as we saw in Figure 1. In our simulation we therefore neglect cooling times for halos with masses $`10^{10}M_{}`$, while using the appropriate pre-enriched cooling times for more massive objects.

Note also the sharp rise in cooling times at temperatures below $`10,000`$ K, due to the lack of atomic transitions with energies below a few eV. Cooling at lower temperatures can only occur through molecular transitions, and thus the collapse of $`T_{\mathrm{vir}}10^4`$ primordial gas clouds depends strongly on the existence of molecular hydrogen. A detailed investigation of the production and dissociation of $`\mathrm{H}_2`$ in the early universe has been conducted by Haiman, Rees, & Loeb (1997). As $`\mathrm{H}_2`$ is easily photo-dissociated by 11.2-13.6 eV photons, to which the universe is otherwise transparent even before reionization, they find that a UV flux of $`10^{22}\mathrm{erg}\mathrm{cm}^2\mathrm{s}^1\mathrm{Hz}^1\mathrm{sr}^1`$, is capable of dissociating all $`\mathrm{H}_2`$ in collapsing halos (see also Ciardi et al. 1999). As this is more than two orders of magnitudes smaller that the reionizing flux, we assume in our simulations that a small population of early stars quickly depleted the primordial gas of molecular hydrogen. Thus we exclude the formation of all galaxies with virial temperatures below $`10,000`$ K, which sets the lower halo mass value to $`5\times 10^7\mathrm{M}_{}h^1`$, justifying the choice of spatial resolution for our simulations.

### 3.4. Heating and Enrichment of the IGM

Once a dwarf galaxy has been formed, we assume that a galactic outflow is immediately generated due to SNeII from the short-lived massive stellar population. In order to get a rough estimate of the mass of metals injected into the ICM, the radius over which they are mixed into the surrounding gas, and the heating of the ICM, we model them according to the simple analytical model developed in Tegmark, Silk, & Evrard (1993, hereafter TSE). Here the authors treat the outflow as an expanding shell that sweeps up most of the baryonic IGM and loses only a small fraction $`f_m=0.1`$ to the interior, such that the mass of the shell is

$$m(t)\frac{4\pi }{3}R(t)^3\rho _b(1f_m)$$
(13)
where $`R(t)`$ is the radius of the outflow. Acceleration of the shell is due to the internal pressure and deceleration from gravitational breaking, both estimated in the thin shell approximation (Ostriker & McKee 1988). With these approximations $`\ddot{R}`$ becomes

$$\ddot{R}=\frac{8\pi pG}{\mathrm{\Omega }_bH^2R}\frac{3}{R}(\dot{R}HR)^2\mathrm{\Omega }_0\frac{H^2R}{2}\frac{GM}{R^2}.$$
(14)
Note that the gravitational term here includes a correction for the presence of a halo of mass $`M`$ out of which the outflow expands. The pressure of the hot gas within the shell is given by $`p=\frac{2E_T}{3V}`$ where $`E_T`$ is the thermal energy of the interior, and energy conservation gives $`\frac{dE_T}{dt}=Lp\frac{dV}{dt}`$, where $`L`$ is the luminosity, incorporating all heating and cooling of the interior plasma.

The evolution of the blast wave is completely determined by the assumed luminosity history of the blast wave. Following TSE we consider five contributions to this luminosity: the cooling by Compton drag against the microwave background ($`L_{\mathrm{comp}}`$), the cooling due to bremsstrahlung and other two body interactions ($`L_{\mathrm{ne}^2}`$) as in Eq. (11), energy injection from supernovae ($`L_{\mathrm{sn}}`$), cooling by ionization of the IGM ($`L_{\mathrm{ion}}`$), and heating from collisions between the shell and the IGM ($`L_{\mathrm{diss}}`$); such that,

$$L=L_{\mathrm{comp}}L_{\mathrm{ne}^2}+L_{\mathrm{sn}}L_{\mathrm{ion}}+L_{\mathrm{diss}}.$$
(15)
The first of these we fix according to Eq. (4) of TSE, while $`L_{\mathrm{ne}^2}`$ is negligible in this regime. The other parameters require some calibration to the observed properties of outflows from dwarf ellipticals and our overall model of the IGM at the time of dwarf formation. Let us consider these individually.

Following TSE we assume that the energy in the blast wave is proportional to the number of SNeII and winds for O and B stars and that this in turn is proportional to the baryonic mass which has collapsed. Energy injection occurs over a relatively short time during which the dwarf galaxy is undergoing a starburst phase. We assume that a fraction of $`ϵ_{\mathrm{SF}}`$ of the baryons in each galaxy collapse to form stars, that one supernova is formed for every 100 $`M_{}`$ of stars in the halo, and that each supernova provides a typical energy output of $`10^{51}`$ ergs, with an equal contribution coming from winds from massive stars. If a fraction $`ϵ_{\mathrm{wind}}`$ of this input goes into the galactic outflow then the total energy from the supernovae is $`E=ϵ_{\mathrm{wind}}ϵ_{\mathrm{sf}}2.0\times 10^{55}\frac{\mathrm{\Omega }_b}{\mathrm{\Omega }_0}M_6.`$ Assuming that this energy is released at a constant rate during a period of $`t_{\mathrm{burn}}=1.0\times 10^8`$ years, we find

$$L_{\mathrm{sn}}=ϵ_{\mathrm{wind}}ϵ_{\mathrm{sf}}1.6\times 10^6L_{}/hM_6\times \frac{\mathrm{\Omega }_b}{\mathrm{\Omega }_0}$$
(16)
where the factor of $`h`$ is a result of the units we have chosen for $`M_6`$. Dekel and Silk (1986) have proposed that the observed differences in surface brightness and metallicity between the observed classes of diffuse dwarfs and normal galaxies can be understood in terms of outflows occurring in halos with virial velocities below critical value on of the order of 100 km/s, which corresponds to approximately $`2.5\times 10^9ϵ_{\mathrm{SF}}(1+z)^{3/2}\mathrm{\Omega }_b/\mathrm{\Omega }_0`$ supernovae. Here we take a conservative approach and exclude all bursts greater than $`2\times 10^7ϵ_{\mathrm{SF}}\mathrm{\Omega }_b/\mathrm{\Omega }_0,`$ which is typically a few hundred thousand supernovae and corresponds to a maximum luminosity of

$$L_{\mathrm{sn},\mathrm{max}}=ϵ_{\mathrm{wind}}ϵ_{\mathrm{sf}}3.2\times 10^9L_{}\times \frac{\mathrm{\Omega }_b}{\mathrm{\Omega }_0}$$
(17)
achieved by objects above $`2\times 10^9M_{}.`$

Another potential source of drag on the expanding shell is energy losses to ionizing the IGM, given by $`L_{\mathrm{ion}}=f_{\mathrm{neutral}}f_\mathrm{m}n_bE_o\mathrm{\hspace{0.17em}4}\pi R^2[\dot{R}HR],`$ where $`n_b`$ is the number density of baryons, $`E_013.6`$ eV, and $`f_{\mathrm{ion}}`$ is the ionization fraction of the IGM into which the shell is expanding. As they were trying to develop a model for IGM reionization by winds, TSE assume $`f_{\mathrm{neutral}}=1`$ at all times. Here we choose for our fiducial model the case in which blast waves expand into a medium that has completely been pre-ionized ($`f_{\mathrm{neutral}}=0`$) possibly by ionizing photons from active galactic nuclei, or from ionization fronts from the massive stars that precede the winds.

Finally $`L_{\mathrm{diss}}`$ accounts for the heating of the interior plasma due to collisions between the shell and the IGM. Taking this to be some fixed fraction $`f_{\mathrm{diss}}`$ of the kinetic energy lost by the shell gives

$$L_{\mathrm{diss}}=f_{\mathrm{diss}}\frac{3m}{2R}(\dot{R}HR)^3.$$
(18)
In TSE the authors consider the extreme cases in which $`f_{\mathrm{diss}}=0`$ and $`f_{\mathrm{diss}}=1`$, showing that their results are relatively insensitive to this parameter. Here we simply fix $`f_{\mathrm{diss}}=0.5`$ throughout.

With these approximations we can use the solution given in TSE to calculate the full expansion history of the galactic outflow, with their $`L_{}`$ and $`E_{}`$ appropriately modified to account for the gravitational potential of the collapsed halo as given by our Eq. (16). We step forward in time (backwards in redshift) in increments of $`\mathrm{\Delta }z=0.01`$ starting at a maximum redshift of $`z_{\mathrm{max}}=25`$, using a fourth-order Runge-Kutta technique. Figure 5 shows the first galaxy to undergo an outflow, a $`1.1\times 10^8M_{}`$ dwarf expanding into an ionized medium at $`z=25`$ in one such simulation.

As our simulations make no attempt to include ionization and do not account for the presence of ionizing radiation from stars within the dwarf galaxy or an external source of UV radiation, we expect our solution to only be accurate while the temperature of the material is greater than 10, 000 K ($`k_bT=13.6`$ eV). After this point, we assume as a reasonable approximation that the halo simply expands with the Hubble flow and remains at a fixed temperature of 10,000 K.

In this figure, we also include for comparison the results of one such blast wave expanding into a neutral medium. This closely traces the ionized solution up until the temperature of the outflow drops to $`10,000`$ K.

Finally, we model the metallicity of the expanding halo by assuming that each supernova typically ejects 2$`M_{}`$ of metals into the intergalactic medium. This number is consistent with the average stellar yields in SNeII simulations as compiled in Nagataki & Sato (1997), although there are significant theoretical uncertainties between the various simulations. The total mass of metal injected into the interstellar medium of the star forming galaxies is then

$$M_{6,z}=0.02ϵ_{\mathrm{sf}}(\mathrm{\Omega }_b/\mathrm{\Omega }_0)M_6.$$
(19)
Our model has two important free parameters, $`ϵ_{\mathrm{sf}}`$ and $`ϵ_{\mathrm{wind}}`$. While the amount of energy in the winds is directly dependent on the product of these two parameters, their ratio is somewhat more unconstrained. Our approach here then will be to fix the star formation efficiency at a somewhat standard value of $`ϵ_{\mathrm{sf}}=0.1`$ and allow the fraction of the SN energy channeled into the wind to vary between 5 and 20 percent to quantify model uncertainties.

With this simple model we can construct from $`𝒳^{M_i}`$ and $`𝒵_\mathrm{c}^{M_i}`$ the full histories of the blast waves expanding into the intergalactic medium: the center of each expanding bubble ($`𝒳^{\mathrm{blast}}\{𝐱_1^{\mathrm{blast}},𝐱_2^{\mathrm{blast}},\mathrm{}\}`$), its comoving radius as a function of redshift ($`^{\mathrm{blast}}(z)\{R_1^{\mathrm{blast}}(z),R_2^{\mathrm{blast}}(z),\mathrm{}\}`$), temperature as a function of redshift ($`𝒯^{\mathrm{blast}}(z)\{T_1^{\mathrm{blast}}(z),T_2^{\mathrm{blast}}(z),\mathrm{}\}`$), and mass of ejected metals ($`_Z^{\mathrm{blast}}\{M_{Z,1}^{\mathrm{blast}},M_{Z,2}^{\mathrm{blast}},\mathrm{}\})`$. While $`^{\mathrm{blast}}(z)`$ and $`𝒯^{\mathrm{blast}}(z)`$ are updated in redshifts intervals of $`\mathrm{\Delta }z=0.01`$, the history is only stored in intervals of $`\mathrm{\Delta }z=0.05`$ to save memory. Note that the jagged appearance at late times of the comoving size of the outflow in Figure 5 is caused by this discrepancy in times steps, rather than by a numerical instability. Note also that as outflows only occur in objects with masses $`10^{10}`$ for which cooling times are negligible, we can compute the entire outflow history of our Monte Carlo volume without considering the impact of metal enrichment on the gas cooling times. Suppression of galaxy formation due to outflows from nearby objects, however, must be considered for outflow-scale second generation objects.

### 3.5. Formation of ‘Second Generation’ Objects

As we saw in §2, galactic outflows affect subsequent galaxy formation both by shocking the IGM as well as enriching it with metals. Shocks can suppress the formation of nearby galaxies by two main mechanisms, heating the halo gas or striping it from the dark matter. In the first case the gas associated with a collapsing halo is heated to above the viral temperature of the halo but remains bound to the collapsing dark matter. The thermal pressure of the gas then overcomes the dark matter potential and the gas expands out of the halo, preventing galaxy formation. In the second case suppression is caused by stripping of the baryons by a shock from a nearby source. The momentum of the moving shock is sufficient to carry with it the gas associated with the halo of an unvirialized perturbation, thus emptying the dark matter halo of its gas and preventing a galaxy from forming.

In a companion letter (Scannapieco, Ferrara, & Broadhurst 2000), we evaluate both these scenarios of halo suppression and find the stripping mechanism to be dominant, because most shock-heated clouds can radiatively cool within a sound crossing time. As in the case of outflow generation, momentum transfer between the outflows and the collapsing halos is a complicated process, dependent on the overall density profile of the collapsing halo and the possibility of smaller collapsed objects within it. Again, our approach in this exploratory study will be to adopt simple criteria to estimate on average when the momentum in a shell is sufficient to remove material from a collapsing halo.

We allow for baryonic stripping to occur only if a halo has not yet virialized, and assume that stripping occurs in all cases in which a shock moves through the center of the pre-virialized halo with sufficient momentum to accelerate the gas to the escape velocity. Thus we exclude all objects for which

$$fM_sv_sM_cv_e,$$
(20)
where $`f=\mathrm{}^2/4r_s^2`$ is the solid angle of the shell that is subtended by the collapsing halo, $`\mathrm{}`$ is the radius of the collapsing region when the shock moves through its center, $`M_s=\frac{4\pi }{3}\mathrm{\Omega }_b\mathrm{\Omega }(z)/\mathrm{\Omega }_0\rho _cr_s^3`$ is the mass of the material swept up by the shock and $`M_c`$ is the baryonic mass of the collapsing halo. We can estimate the comoving radius of the collapsing halo as

$$\mathrm{}=R_{M_j}(1+\delta _{\mathrm{NL}})^{1/3},$$
(21)
where $`\delta _{\mathrm{NL}}`$ is the nonlinear overdensity of the region, which is given in the spherical collapse model as

$$1+\delta _{\mathrm{NL}}=\frac{9}{2}\frac{(\theta \mathrm{sin}\theta )^2}{(1\mathrm{cos}\theta )^3},$$
(22)
where the collapse parameter $`\theta `$ is given by $`(\theta \mathrm{sin}\theta )^{2/3}\pi ^{2/3}=D(z_{\mathrm{cross}})/D(z_c)`$ where $`z_{\mathrm{cross}}`$ is the redshift at which the shock moves through the center of the halo. Note that in the elliptical collapse case, $`z_c`$ is computed according to Eq. (4), slowing down the condensation of the cloud somewhat. Note also that we treat halos with subcondensations that have already virialized in the same way. In this case the shock would lose energy impinging on the condensed sub-objects but would also have less uncondensed gas to accelerate to escape velocity. We assume for our purposes here that these effects roughly compensate, although more detailed hydrodynamical studies would help to refine this comparison.

Finally, we calculate the metallicity of each collapsing sphere by volume averaging the contributions of each of the expanding blast waves that pass within the collapse radius $`r_{\mathrm{col}}`$. Each halo is assigned a collapse mass in metals $`M_{Z,k}^{M_i}`$, which is taken to be zero initialy and modified by each blast front that passes within the collapse radius. For each such occurrence, the mass in metals is updated to

$$M_{Z,k}^{M_i}M_{Z,k}^{M_i}+\frac{V_{\mathrm{overlap}}}{\frac{4\pi }{3}R_{M_i}^3}M_Z^{\mathrm{blast}}$$
(23)
where $`V_{\mathrm{overlap}}`$ is the volume of intersection of the two spheres. By dividing this mass by the total baryonic mass of the galaxy we can compute the initial metallicity of the object. This value is then used to compute the collapse time of the object by using the appropriate cooling function, $`\mathrm{\Lambda }(T,Z)`$.

## 4. Results

In this section we summarize the results of our simulations, both in relation to the suppression of low-mass galaxy formation by baryonic stripping of pre-virialized halos, as well as the pre-enrichment of high mass galaxies by metals carried in galactic outflows. We show that the tendency for small galaxies to disrupt the formation of their neighbors helps to explain the small number of Milky Way satellites relative to Cold Dark Matter model predictions, and leads to a lack of small-scale galaxy mergers, suggesting a bell-shaped luminosity function for elliptical galaxies. Pre-enrichment helps to explain the high G-dwarf metallicities found in our own and other galaxies, but the shorter cooling times of metal enriched clouds seems to have little effect on the overall number of large galaxies formed. We examine the effect of varying model parameters on these features and discuss the compatibility between the star formations rates and IGM heating in our models and constraints from optical observations and measurements of thermal distortions in the cosmic microwave background.

### 4.1. Shocks and Low-Mass Galaxy Formation

Having constructed a distribution of halos with collapse times and spatial orientation in good agreement with analytical expressions, we now consider the impact of outflows and enrichment on galaxy formation within this distribution. We first restrict our attention to the $`\mathrm{\Lambda }`$CDM model of Figure 1, as the currently favored cosmological model, and consider the effects of varying the cosmology and other model parameters in §4.4. In Figure 6 we show the collapse and formation of outflows in our fiducial model with $`ϵ_{\mathrm{wind}}=0.1`$. Here we see that regions effected by galactic outflows are highly correlated with regions in which new objects are collapsing, illustrating the necessity of an inhomogeneous approach.

In Figure 7 we show the volume fraction in the expanding shells as a function of $`z`$. This quantity is not estimated on a cell-by-cell basis, but rather approximated by its value in the case in which the positions of the bubbles are uncorrelated,

$$\mathrm{F}(z)1\mathrm{exp}\left(\frac{\frac{4\pi }{3}_{i=1}^NR_i^{\mathrm{blast}}(z)^3}{(12\mathrm{M}\mathrm{p}\mathrm{c}/h)^3}\right).$$
(24)
Note that this quantity overestimates the volume within outflows as the correlations between outflow positions can be significant, and are even greater in a full N-body approach than in our simple Monte Carlo.

Also shown on this plot is the volume-averaged temperature within the shells and the overall volume-averaged temperature of the simulation. Here we see that at redshifts much greater than $``$ 5, the areas near dwarf starburst galaxies were significantly hotter than the IGM as a whole. This discrepancy is somewhat reduced at later times, when the outflows begin to fill a significant fraction of the total volume of the universe. Even at these late times, however, the universe remains inhomogeneous, with new starbursting dwarf galaxies forming within regions only mildly impacted by the earliest outflowing objects. Thus the root mean squared (RMS) fluctuation in the temperature remains significant at redshifts $``$ 1.

This high degree of inhomogeneous shocking naturally has a large impact on galaxy formation at masses below $`10^{10}M_{}`$. In Figure 8 we plot the number of collapsed halos as a function of redshift, along with the number of unsuppressed objects as described in §3.4. While the earliest dwarf galaxies form along with the collapsing halos, by $`z12`$ outflows become important. As $`M10^9M_{}`$ halos collapse after this redshift, galaxy formation at this scale is highly suppressed, leading to a “mass desert” between the major population of polluting dwarf galaxies and the larger galaxies that are seen today. It is also clear from these figures that galaxy formation is likely to have occurred in two stages, with the dwarf outflow population forming mostly at the highest redshifts, and larger galaxies forming only later when halos of sufficient mass began to collapse.

Also in this figure, we plot the total number density of collapsed halos with masses between $`1.1\times 10^8M_{}`$ and $`1.4\times 10^{11}M_{}`$, along with the total number of galaxies in this mass range. Notice that only $`30\%`$ of the halos are populated, consistent with the factor of $`4`$ suppression needed to reconcile the number of predicted and observed Milky-Way satellites (Kauffmann, White, & Guiderdoni 1993; Klypin et al. 1999; Moore et al. 1999). Our scenario thus provides a natural mechanism for the formation of “dark halo” satellites around our galaxy, which may be associated with the abundant High-Velocity Clouds as discussed by Blitz et al. (1999).

### 4.2. Enrichment and High-Mass Galaxy Formation

In Figure 9, we plot the volume-averaged metallicity within the outflows and the simulation overall, as well as the RMS variation between the shells. Here we see that like the temperature, the mean metallicity within the outflows is the greatest at early times, and decreases to a few times the mean metallicity of the universe at a redshift of $`5`$. The mean IGM metallicity from $`z=1`$ to $`z=5`$ ranges from about $`\mathrm{log}_{10}[Z/Z_{}]1.8`$ to $`\mathrm{log}_{10}[Z/Z_{}]2.6`$ to in good agreement with observed metallicities in Ly absorption systems (e.g. Cowie & Songaila 1998, McDonald et al 1999, Ellison et al. 1999), as well as previous studies of enrichment by dwarfs (Nath & Trentham 1997). The RMS variations in the metallicities within the outflows are even more severe than the temperature variations however, and remain several times larger than the mean metallicity at all times.

This highly inhomogeneous distribution of metals is consistent with observed metallicity inhomogeneities as discussed in §2.2 and suggests that the mean metallicity of the gas out of which most galaxies formed may be significantly higher than that of the IGM. Figure 9 also shows the initial metallicity of the IGM newly formed galaxies. Note that this does not include the contribution from the metals in the progenitor galaxies themselves. Here we see that indeed the mean metallicity of collapsed objects is an increasing function of mass which is expected to be $`\mathrm{log}_{10}[Z/Z_{}]1`$ for Milky-Way sized objects, and higher than the mean IGM metallicity at for all mass scales. This value is comparable to the $`0.1Z_{}`$ pre-enrichment necessary to explain the lack of low-metallicity G-dwarfs in the Solar neighborhood. The overall relatively high initial metallicity is also consistent with observations that suggest the need for a “floor” of around $`0.1Z_{}`$ in the metallicity distribution of stars in other massive galaxies. Note that we make no attempt to trace subsequent star formation and interstellar medium (ISM) enrichment, and our values only provide a lower limit on the lowest metallicity stars in disk galaxies that formed largely from the IGM gas. Note also that these metallicities are slowly increasing with mass, suggesting an alternative explanation for the mean metallicity luminosity relationship, which is usually understood as ISM enrichment within galaxies and their progenitors (see eg. Vader 1986; Ferguson & Binggeli 1994; Kauffmann & Charlot 1998).

As collapsing halos in which galaxy formation is suppressed are found closest to outflowing galaxies, where both shock velocities and metallicities are high, it is likely that whatever residual gas remains in these objects would be of similar of even higher metallicity than that of unsupressed objects. While the calculation of these metallicities is beyond the scope of our simulations, it is interesting to note that some of the high velocity clouds may also have metallicities that actually exceed those of many dwarf galaxies although may still be somewhat less than suprasolar values predicted by a “galactic-fountain” model (Sembach et. al 1999; Wakker et al. 1999).

In Figure 10 we show the effect of cooling on galaxy formation in large-mass halos. Here we see that while halo collapse and galaxy formation are almost simultaneous at masses $`5\times 10^{10}M_{}`$, long cooling times suppress all galaxy formation on scales $`5\times 10^{12}M_{}`$. The increased initial metallicities help to accelerate galaxy formation in objects between these mass limits somewhat, although the overall final number densities remain unchanged.

### 4.3. Properties of Elliptical Galaxies

As elliptical galaxies tend to be found in the most dense and enriched regions of space, it is natural to expect that outflows would have had the largest impact on these objects. In order to study this connection, we follow the conventional wisdom that ellipticals correspond to mergers of (see e.g., Barnes 1992; Hernquist 1993) large progenitors.

We therefore define $`M_{\mathrm{merger},k}^{M_i}`$ as the total contribution to the mass of halo $`k`$ of mass $`M_i`$ due to mergers of unsupressed galaxies of mass scales $`M_{i1}`$ and $`M_{i2}`$, that is

$`M_{\mathrm{merger},k}^{M_i}`$ $`{\displaystyle \underset{l}{}}\{\begin{array}{cc}M_{i1}\hfill & \text{if }𝐱_k^{M_i}𝐱_l^{M_{i1}}R_{M_i}\hfill \\ 0\hfill & \text{ otherwise}\hfill \end{array}`$ (25)
$`+{\displaystyle \underset{l}{}}\{\begin{array}{cc}M_{i2}\hfill & \text{if }𝐱_k^{M_i}𝐱_l^{M_{i2}}R_{M_i}\hfill \\ 0\hfill & \text{ otherwise}\hfill \end{array}`$
where we consider only objects that have cooled and have not been swept away by shocks while forming. We then identify as ellipticals all collapsed objects with $`M_{\mathrm{merger},k}^{M_i}0.5M_i`$, in which half of the total mass comes from a merger of large progenitors.

We fix this threshold in order to approximate the observed field elliptical fraction of $`15\%`$ (Baugh, Cole, & Frenk 1996), and ask the question of what mass scales correspond to these objects.

In Figure 11 we plot the total space density of galaxies in our simulations, along with the space density of galaxies identified as ellipticals. In this figure we see that the mass distribution of ellipticals is quite different than that of the overall galaxy population. As the outflows from dwarf galaxies suppress the formation of nearby dwarfs, there are very few mergers at small masses, with virial temperatures much smaller than that typical of outflows. Thus while the number density of halos continues to increase with decreasing mass, the number density of ellipticals is the same at $`1.4\times 10^{11}M_{}`$ and $`4.4\times 10^{10}M_{}`$, and decreases dramatically for even lower masses.

Using the observed Faber-Jackson relation between the luminosity of elliptical galaxies and the velocity dispersion, we can recast this as an absolute r-band magnitude of $`M_r=1810\mathrm{log}_{10}\sigma _{100}`$, (Oegerle & Hoessel 1991) where $`\sigma _{100}`$ is the velocity dispersion in units of $`100`$ km/s. This gives the luminosity function of elliptical galaxies shown in the bottom of Figure 11. While these results are necessarily crude, and are bounded at the high-luminosity end by the finite size of out simulations, they nevertheless naturally reproduce the observed bell-shaped luminosity function of elliptical galaxies (see eg., Bromley et al. 1998).

### 4.4. Effects of Model Uncertainties

We have only two free parameters in our simplified modeling, both relating to the gas outflow, and all other modeling based on a simple spherical or elliptical collapse model, involving the usual assumptions. Our conclusions are most strongly dependent on the product of the two free parameters $`ϵ_{\mathrm{wind}}`$ and $`ϵ_{\mathrm{sf}}`$, which together determine the total energy being channeled into galactic outflows. In our fiducial model, we took $`ϵ_{\mathrm{sf}}`$ = 0.1 and $`ϵ_{\mathrm{wind}}ϵ_{\mathrm{sf}}=.01`$. These canonical values are motivated by the current observational work on dwarf galaxy outflow energetics, but a large spread is naturally expected and thus we have some freedom in changing this efficiency. We therefore examined two extreme cases with $`ϵ_{\mathrm{wind}}=.05`$ and $`ϵ_{\mathrm{wind}}=.2`$ to asses the robustness of the features described above. These cases bound the shaded regions of the results plotted in Figures, 7, 8, 11, 12.

In Figure 7 we see the impact of model uncertainties on the total volume fraction and temperatures of outflows in our simulations. The impact of these differences on the suppression of dwarf galaxies is shown in the shaded regions in Figure 8. Here we see that the existence of two-stage evolution with a “mass desert” in the $`10^910^{10}M_{}`$ range persists for all values of these parameters. The total number of unsuppressed dwarf galaxies ranges from $`10\%30\%`$ reinforcing the idea that outflow suppression may be important in reproducing the small number of Milky-Way satellites.

The effect of varying $`ϵ_{\mathrm{sf}}ϵ_{\mathrm{wind}}`$ on the distribution of elliptical galaxies is shown in Figure 11. While the uncertainties in this case are somewhat more severe than in Figure 8, a large fall-off in the number of objects with masses below $`5\times 10^{10}M_{}`$ is present for all cases.

We show the metallicities of the IGM component of the forming galaxies in each of these cases in the upper panel of Figure 12. In both cases the mean metallicity of collapsed objects is an increasing function of mass which is higher than the mean IGM metallicity. The predicted Milky-Way scale IGM metallicity is again always $`\mathrm{log}_{10}[Z/Z_{}]1`$ and thus the G-dwarf problem is easily understood.

In the lower panel of this figure, we compute the star formation rate in dwarf galaxies to assuage any concern that our code might be pumping an inordinate amount of energy into the IGM. Note that our models do not include subsequent star-formation in disk galaxies and thus this figure represents a lower bound to the number of stars formed at each redshift. Nevertheless, our values are well within observational constraints (see eg., Adelberger & Steidel 2000).

We can also compare our models to the COBE constraints on the optical depth to the surface of last scatter as well as the the distortion of the cosmic microwave background (CMB) spectrum due to the Sunyaev-Zel’dovich effect. We can estimate the optical depth in our model as

$$\tau =5.9\times 10^3\frac{\mathrm{\Omega }_b}{\mathrm{\Omega }_0}_0^{z_{\mathrm{max}}}𝑑z\frac{dx}{dz}\mathrm{\Omega }(z)(1+z)^2\mathrm{F}(z),$$
(26)
where $`x`$ is again the comoving distance defined such that $`a_0H_0=c`$. This analysis gives $`\tau =0.05`$ for the fiducial model, and varies from $`\tau =0.04`$ to $`\tau =0.06`$ in the low and high energy outflow cases, which are all well within the observational limit of $`\tau 0.5`$ (Griffiths, Barbosa, & Liddle 1999).

The degree of CMB spectral distortions is given by the Compton-$`y`$ parameter which is the convolution of the optical depth with the electron temperature along the line of sight (Zel’dovich & Sunyaev 1969; Sunyaev & Zel’dovich 1972). Thus

$$y=1.0\times 10^8\frac{\mathrm{\Omega }_b}{\mathrm{\Omega }_0}_0^{z_{\mathrm{max}}}𝑑z\frac{dx}{dz}T_5(z)\mathrm{\Omega }(z)(1+z)^2\mathrm{F}(z),$$
(27)
where $`T_5(z)`$ is the mean temperature within the outflows as a function of redshift in units of $`10^5`$K. In this case $`y`$ is $`1.6\times 10^6`$ in the fiducial model, and varies from $`1.0\times 10^6`$ to $`2.3\times 10^6`$ for the extreme cases. These values are again within the observational constraint of $`y1.5\times 10^5`$ (Fixsen et al. 1996). Nevertheless, our modeling suggests that the spectral distortions due to galactic outflows should be observable with the next generation of cosmic microwave background experiments.

Varying the second parameter in our model, $`ϵ_{\mathrm{sf}}`$, while keeping $`ϵ_{\mathrm{sf}}ϵ_{\mathrm{wind}}`$ fixed acts as an overall shift in the star formation rate and metallicities in Figures 9 and 12 and can be estimated simply by eye. Including the ionization drag term in Eq. (11) and allowing the shocks to act as the ionizing source of the IGM has little effect on our results as suggested by Figure 5 and verified by a full simulation. Thus we see that all the major features of the fiducial model persist over a wide range of parameters.

We have also examined the cosmological dependence of our results by examining a flat model with parameters as described in §2.3 and $`ϵ_{\mathrm{wind}}=0.1`$. In this case only $`50\%`$ suppression of galaxies occurs, again with the most impact on galaxies in the few times $`10^9M_{}`$ range. In this model elliptical formation falls off below the $`10^{11}M_{}`$ scale but more gradually than in the open case, and initial Milky-Way metallicities again about $`\mathrm{log}_{10}[Z/Z_{}]=1.0`$. Thus even in a model in which structure forms much more quickly than suggested by observations (see eg., Bahcall, Fan, & Cen 1997), the major features in our model persist.

## 5. Discussion

Our treatment of galactic outflows is intended as an exploratory study of these processes, and makes a number of approximations that should be made explicit. While a Monte Carlo approach uncovers many of the issues of inhomogeneity that can not properly be studied analytically, it fails to capture the nonlinear structure and clustering present in a full N-body treatment. At late times peaks will be grouped along sheets and filaments, and thus clustering is likely to be more severe than in our simulations and of a more complicated nature.

A number of uncertainties also arise from our simple treatment of IGM enrichment and shocking. One concern is the presence of radiative feedback from ionization fronts around dwarf galaxies. Our assumption that fronts precede the formation of outflows, and thus each shell can be treated as expanding into an ionized medium is probably a reasonable one, while our placement of a “temperature floor” at the ionization temperature of hydrogen at all times in our simulations is stronger approximation. Fortunately the extreme fragility of molecular hydrogen suggests that the formation of objects with virial temperatures below $``$10,000 K may have not been largely affected by reionization. While the inhomogeneous structure of reionization is likely to have had a great impact on the CMB fluctuations (Aghanim et al. 1996; Miralda-Escude, Haehnelt, & Rees 2000; Scannapieco 2000), the formation of $`10^7M_{}`$ galaxies is likely to have been halted by the dissociation of molecular hydrogen by the first stars, long before reionization took place.

A bigger concern may be the structures of outflows and the ejection fractions of interstellar gas and metals. This has been studied in detail by Mac Low & Ferrara (1999) and Ferrara & Tolstoy (2000), who conclude that efficient ejection of the interstellar medium or “blowaway” occurs only in halos with masses $`10^7M_{}`$. The issue of whether dwarfs retain a sizeable fraction of their gas, however, is to some degree decoupled from the formation of outflows. In halos in the mass range $`10^7M_{}M10^9M_{}`$, a “blowout” occurs in which the super-bubbles around groups of SNeII punch out of the galaxy, shocking the surrounding IGM and efficiently ejecting metals while failing to excavate the interstellar medium of the galaxy as a whole. Such a scheme in which outflows are formed primarily of the enriched IGM surrounding dwarf galaxies is equivalent to outflows of galactic gas for the purposes of our simulations. Additionally such a picture would help to reconcile observations of expanding high metallicity shells around dwarf galaxies as described in §2 with observations of multiple episodes of and HI gas in dwarf spheroidal galaxies (Smecker-Hane et al. 1994; Grebel 1998) some of which even suggest that many of these objects are gas-rich, but with extended HI envelopes (Blitz & Robishaw 2000).

That being said, the TSE model adopted for outflow evolution in our simulations is almost certainly oversimplified. The “blowout” scenario described by Mac Low & Ferrara (1999) ejects matter perpendicular to the thick disk, and is thus quite different from the spherical outflow modeled here. Even if the winds from dwarf galaxies can be reasonably approximated as spherical shells at early times, such outflows of low-density heated gas would necessarily become Rayleigh-Taylor unstable as they expanded into the denser IGM. Thus material is most likely to flow out in a number of heated clumps at large distances.

Also ignored in our model are shell-shell interactions which would contribute additional heating but slow expansion in regions which two shocks meet. Similarly, an additional source of IGM enrichment is the merger mechanism described in Gnedin & Ostriker (1997) and Gnedin (1998), in which a significant fraction of the interstellar medium of merging galaxies is ejected.

Finally, the criteria for halo suppression used in our simulation is over-simplified. While regions of space in our simulations that meet the stripping criterion, Eq. (20), are likely to have delayed or suppressed galaxy formation, the one-to-one correspondence assumed in our simulations is unlikely. In reality the formation of galaxies is likely to be a complicated function of the number and strength of the shocks moving through the regions, and the structure of these regions during and after halo collapse. This subject merits further investigation and would help to sharpen our conclusions.

## 6. Conclusions

While the impact of preheating and enrichment on the observed properties of galaxy clusters has long been recognized, the impact on the observed properties of galaxies themselves has been little explored. By accounting for the observed outflows from dwarf galaxies, we have been able to show how many of the unexplained properties of galaxies and the IGM can be naturally understood.

Firstly, the suppression of low-mass galaxy formation by outflows provides a natural explanation for the factor of $`4`$ discrepancy between the number of observed Milky-Way satellites and predictions from standard CDM models that do not include outflows. Suppression also provides a natural mechanism for the formation of “dark halos” which may be associated with the High-Velocity clouds.

Secondly, baryonic stripping results in a bell-shaped luminosity function of ellipticals. The lower the mass of a halo, the more likely it is to generate an outflow that strips material from a similar mass neighboring pre-virialized halo that would otherwise later form into a galaxy. This results in very few pairs of neighboring low-mass galaxies, and hence a relative deficit of major low mass-mergers.

Finally, our models of enrichment of protogalactic gas predict a trend of increasing metallicity with galaxy mass in good agreement with inferences from observations. The initial metallicity predicted for a Milky-Way mass galaxy is $`0.1Z_{}`$ providing a natural initial floor at the level required to solve the G-dwarf problem and the more general lack of low metallicity stars in well studied massive elliptical galaxies relative to “closed-box” models of chemical enrichment.

These results are persistent over a wide range of model parameters and cosmologies, and are not a result of fine-tuning parameters or invoking additional physics.

While galaxy outflows are already incorporated into modern studies of galaxy formation, this is done only as an internal modification, ignoring pre-enrichment. Keeping track of the effect of outflows on neighboring halos is essential in understanding the properties of galaxy clusters, and hence it is not surprising that these effects would play a major role in the formation of galaxies as well. While the details await further investigation, it is clear that any complete picture of galaxy formation must account for heating and enrichment.

We would like to thank Rychard Bouwens, Andrew Cumming, Julianne Dalcanton, Marc Davis, Andrea Ferrara, Ignacio Ferreras, Brenda Frye, John Huchra, Siang Peng Oh, Alvio Renzini, David M. Sherfesee, Joseph Silk, Jonathan C. Tan, Robert Thacker, and Simon White for helpful comments and discussions. This research was supported in part by the National Science Foundation under Grant PHY94-07194. TJB acknowledges NASA grant AR07522.01-96A.


format: text

text:
# Symmetry Dependence of Localization in Quasi- 1- dimensional Disordered Wires

## Abstract

The crossover in energy level statistics of a quasi-1-dimensional disordered wire as a function of its length $`L`$ is used, in order to derive its averaged localization length, without magnetic field, in a magnetic field and for moderate spin orbit scattering strength. An analytical function of the magnetic field for the local level spacing is obtained, and found to be in excellent agreement with the magnetic field dependent activation energy, recently measured in low-mobility quasi-one-dimensional wires. This formula can be used to extract directly and accurately the localization length from magnetoresistance experiments. In general, the local level spacing is shown to be proportional to the excitation gap of a virtual particle, moving on a compact symmetric space. Pacs- numbers: 72.15.Rn,73.20.Fz,02.20.Qs

In disordered wires quantum interference results in localization of all states for arbitrary disorder strength, if the wire is infinite. It has been discovered that the localization length depends on the global symmetry of the wire : $`L_c=\beta \pi \mathrm{}\nu SD_0`$, where $`\beta =1,2,4`$, corresponding to no magnetic field, finite magnetic field, and strong spin- orbit scattering or magnetic impurities, respectively. $`\nu (E)`$ is the electronic density of states in the wire. $`D_0`$ is the classical diffusion constant of the electrons in the wire, and $`S`$ its crossection. This result was obtained by calculating the spatial decay of the density correlation function for wires whose thickness exceeds the mean free path $`l`$. Independently, it was obtained by calculating the transmission probability through thin, few channel wires. Recently, the doubling of the localization length was observed in sub-micron thin wires of doped GaAs by Khavin, Gershenson and Bogdanov, who found a continously decreasing activation energy when the magnetic field is increased, saturating indeed at one half of its field free value . However, the crossover as function of moderate symmetry breaking fields defied any attempt to study it with nonperturbative methods. The only results were based on a heuristic approach by Bouchaud, a kind of semiclassical analysis by Imry and Lerner and numerical studies. Recently, however, Kolesnikov and Efetov succeeded to tackle that complex problem and derived the density- density correlation function with the supersymmetry method in the crossover regime. The result does not scale with a single localization length.

Here, the problem of the crossover is addressed from a different perspective. The statistics of discrete energy levels of a finite coherent, disordered metal particle is an efficient way to characterize its properties . This can be studied by calculating a disorder averaged autocorrelation function between two energies at a distance $`\omega `$ in the energy level spectrum. It is an oscillatory function whose amplitude decays with a power law when the energy levels in the vicinity of the central energy $`E`$ are extended. A Gaussian decaying function is a strong indication that all states are localized. The autocorrelation function of spectral determinants (ASD) is defined by $`C(\omega )=\overline{C}(\omega )/\overline{C}(0),\overline{C}(\omega )=\text{det}(E+\omega /2H)\text{det}(E\omega /2H).`$ $`E`$ is a central energy. We consider the hamiltonian of disordered electrons

$$H=\left[𝐩q𝐀\right]^2/2m+V(𝐱),$$
(1)
where $`q`$ is the electron charge, $`c`$ the velocity of light. $`V(𝐱)`$ is taken to be a Gaussian distributed random function $`V(𝐱)=0`$, and $`V(𝐱)V(𝐱^{})=(\mathrm{\Delta }SL/2\pi \tau )\delta (𝐱𝐱^{}),`$ which models randomly distributed, uncorrelated impurities in the sample. $`\mathrm{\Delta }`$ is the mean level spacing, $`1/\tau `$ the elastic scattering rate. The vector potential is used in the gauge $`𝐀=(By,0,0)`$, where $`x`$ is the coordinate along the wire of length L, $`y`$ the one in the direction perpendicular both to the wire and the magnetic field $`𝐁`$, which is directed perpendicular to the wire. The angular brackets denote averaging over impurities.

The ASD appears as an intermediate step in the Grassman replica trick and as the compact sector of the supersymmetric field theory of disordered systems . It is a non-self-averaging quantity. It has been shown recently, that the ASD can distinguish between an uncorrelated spectrum of localized states and a correlated spectrum of extended states in the vicinity of $`E`$. This way, the metal- insulator crossover as a function of the length of a disordered wire has been explored with the ASD. A localization length $`\xi _c`$ in a moderate magnetic field has been derived as the crossover length scale and shown to coincide with $`L_c(\beta =2)`$. The ASD also yields qualitative information on the location of localized states in a quantum- Hall- system .

The impurity averaged ASD can be written as a partition function

$$\overline{C}(\omega )=\text{Tr}\mathrm{exp}(L\overline{H}\left[Q\right]),$$
(2)
where $`\overline{H}`$ is an effective Hamiltonian of matrices Q on a compact manifold, determined by the symmetries of the Hamiltonian $`H`$ of disordered electrons. Thus, the problem reduces to the one of finding the spectrum of the effective Hamiltonian $`\overline{H}`$.

There is a finite gap $`E_G`$ between the ground state energy and the energy of the next excited state of $`\overline{H}(\omega =0)`$. For a long wire, $`LE_G1`$, the ASD becomes, according to Eq. (2), $`C(\omega )=\mathrm{exp}(const.L\omega ^2/E_G)`$. This is typical for a a spectrum of localized states. In the other limit $`LE_G1`$, all modes of $`\overline{H}`$ do contribute to the trace in the partition function Eq. (2) with equal weight, yielding the correlation function of a spectrum of extended states. Thus, the crossover length $`\xi _c1/E_G`$ can be identified with an averaged localization length.

Being a product of two spectral determinants, $`4\alpha `$ -component Grassman fields are needed to get the functional integral representation of the ASD. Here, $`\alpha =1`$, when the Hamiltonian is independent of the spin of the electrons, and each level is doubly spin degenerate. There is one pair of Grassman fields for each determinant in the ASD and each pair is composed of a Grassman field and its time reversed one, as obtained by complex conjugation. $`\alpha =2`$ has to be considered, when the Hamiltonian does depend on spin, as for the case with moderately strong magnetic impurity or spin- orbit scattering. This necessitates the use of a vector of a spinor and the corresponding time reversed one. One notes a global invariance under rotations of these vectors.

Averaging over impurities according to Eq. (1), one gets an interacting theory of Grassman fields with interaction strength $`1/\tau `$. It can be decoupled with a Gaussian integral over $`4\alpha \times 4\alpha `$ matrices which preserve the global rotational invariance. The Grassman fields can now be integrated out exactly. Finding the saddle point of the integral over matrices and integrating out longitudinal (massive) modes, for $`\omega 1/\tau `$ and $`1/\tau \mathrm{\Delta }`$, the ASD reduces to a functional integral over transverse modes Q, being elements of respective symmetric spaces. For disordered wires whose thickness $`W`$ is larger than the mean free path $`l`$ but smaller than the magnetic length $`W<l_B`$, the action reduces to the one of the one-dimensional compact nonlinear sigma model ,

$$\overline{C}(\omega )=\underset{x}{}\text{d}Q(x)\mathrm{exp}(F[Q]),$$
(3)
$`F[Q]`$ $`=`$ $`\alpha {\displaystyle \frac{1}{16}}L_{CU}{\displaystyle _0^L}\text{d}x\text{Tr}\left[(_xQ(x))^2{\displaystyle \frac{q^2}{\mathrm{}^2}}\overline{y^2}B^2[Q,\tau _3]^2\right]`$ (4)
$`+`$ $`i\alpha {\displaystyle \frac{\pi }{4}}{\displaystyle \frac{\omega }{\mathrm{\Delta }}}{\displaystyle \frac{\text{d}x}{L}\text{Tr}\mathrm{\Lambda }_3Q(x)}.`$ (5)
where $`L_{CU}=L_C(\beta =2)=2\pi \mathrm{}\nu SD_0`$ is the localization length in the wire in a moderately strong magnetic field . The integral over the crossection $`S`$ of the wire is done, giving $`\overline{y^2}`$. Here, and in the following, $`\mathrm{\Lambda }_i`$ are the Pauli matrices in the subbasis of the left and the right spectral determinant, $`\tau _i`$ the ones in the subbasis spanned by time reversal and $`\sigma _i`$ the ones in the space spanned by the spinor, for $`i=1,2,3`$. The fluctuations of the matrices $`Q`$ are transverse, $`Q^2=1`$, and restricted by global rotational invariance to $`Q^+=Q`$. For $`\alpha =1`$ in addition $`Q^TC=CQ`$, where $`C=i\tau _2`$. This constrains $`Q`$ to be an element of the group defined on the compact symplectic symmetric space Sp(2)$`/`$( Sp(1) $`\times `$ Sp(1)). For a moderate magnetic field, $`Q`$ is reduced to a $`2\times 2`$\- matrix by the broken time reversal symmetry. This reduces the space of Q to $`U(2)/(U(1)\times U(1))`$. For $`\alpha =2`$ the matrix $`C`$ is, due to the time reversal of the spinor, substituted by $`i\sigma _2\tau _1`$. Both magnetic impurities and spin-orbit scattering reduce the Q matrix to unity in spin space. Thus, C has effectively the form $`\tau _1`$. The condition $`Q^TC=CQ`$ leads therefore to a new symmetry class, when the spin symmetry is broken but the time reversal symmetry remains intact. This is the case for moderately strong spin-orbit scattering. Then, $`Q`$ are $`4\times 4`$\- matrices on the orthogonal symmetric space $`O(4)/(O(2)\times O(2))`$ . With magnetic impurities both the spin and time reversal symmetry is broken, and the Q- matrices are in the unitary symmetric space $`U(2)/(U(1)\times U(1))`$ as for a moderate magnetic field and spin degenerate levels. The difference in the prefactor $`\alpha `$ in the action, Eq.(4), remains. One can extend this approach to other compact symmetric spaces with physical realizations, see Ref.  for a complete classification. For $`\omega /\mathrm{\Delta }<L_{CU}/L`$, the spatial variation of $`Q`$ can be neglected and one retains the same ASD as for random matrices of orthogonal, unitary and symplectic symmetry, characterizing the energy levels of an ergodic particle without magnetic field, with magnetic field, and with spin- orbit- scattering, respectively. The confusing fact, that random orthogonal matrices result for the ASD in a functional integral over compact symplectic Q- matrices and vice versa, results from the sign change of the Grassman variables under time reversal.

We can derive the corresponding Hamiltonian $`\overline{H}`$ by means of the transfer matrix method, reducing the one-dimensional integral in Eq. (3) to a single functional integral. Thus, the ASD is obtained in the simple form of Eq. (2), with the effective Hamiltonian

$$\overline{H}(\omega =0)=\frac{1}{\alpha L_{CU}}(4\mathrm{\Delta }_Q\frac{1}{16}X^2Tr_Q[Q,\tau _3]^2).$$
(6)
$`\mathrm{\Delta }_Q`$ is that part of the Laplacian on the symmetric space, which couples to $`Tr[\mathrm{\Lambda }_3Q]`$. $`X=2\pi \alpha L_{CU}(\overline{y^2})^{(1/2)}B/\varphi _0`$, where $`\varphi _0=q/h`$ is the flux quantum.

The problem is now equivalent to a particle with “mass” $`(\alpha /8)L_{CU}(E)`$ moving on the symmetric space of $`Q`$ in a harmonic potential with “frequency” $`X/(\alpha L_{CU})`$, and in an external field $`i\alpha (\pi /4)\omega /(L\mathrm{\Delta })`$, in “time” $`x`$, the coordinate along the wire. To find the ASD as a function of $`\omega `$ and the length of the wire $`L`$, one can do a Fourier analysis in terms of the spectrum and eigenfunctions of the effective Hamiltonian at zero frequency, $`\overline{H}(\omega =0)`$ . For $`L\xi _c`$, only the gap between the ground and the first excited state of $`\overline{H}(\omega =0)`$, determines the ASD.

Without magnetic field, $`B=0`$, the Laplacian is

$`\mathrm{\Delta }_Q`$ $`=`$ $`_{\lambda _C}(1\lambda _C^2)_{\lambda _C}+2{\displaystyle \frac{1\lambda _C^2}{\lambda _C}}_{\lambda _C}`$ (7)
$`+`$ $`{\displaystyle \frac{1}{\lambda _C^2}}_{\lambda _D}(1\lambda _D^2)_{\lambda _D},`$ (8)
where $`\lambda _{C,D}[1,1]`$. Its ground state is $`1`$ and the first excited state is $`\lambda _C\lambda _D`$. Thus, the gap is

$$E_G(B=0)=16/L_{CU}.$$
(9)
For moderate magnetic field, with the condition $`L_{CU}(\overline{y^2})^{(1/2)}B\varphi _0=h/q`$, all degrees of freedom arising from time reversal invariance are frozen out, due to the term $`Tr_Q[Q,\tau _3]^2=16(\lambda _C^21)`$ which fixes $`\lambda _C^2=1`$. Then, the Laplacian reduces to

$`\mathrm{\Delta }_Q=_{\lambda _D}(1\lambda _D^2)_{\lambda _D}.`$ (10)
Its eigenfunctions are the Legendre polynomials. There is a gap above the isotropic ground state of magnitude

$$E_G(X1)=8/L_{CU}.$$
(11)
For moderate magnetic impurity scattering exceeding the local level spacing, $`1/\tau _s>\mathrm{\Delta }_C`$, the Laplacian is given by Eq.(8). Thus, due to $`\alpha =2`$, the gap is reduced to $`E_G(1/\tau _S>\mathrm{\Delta }_C)=4/L_{CU}`$. For moderately strong spin- orbit scattering $`1/\tau _{SO}>\mathrm{\Delta }_C`$, the Laplace operator is

$$\mathrm{\Delta }_Q=\underset{l=1,2}{}_{\lambda _l}(1\lambda _l^2)_{\lambda _l},$$
(12)
where $`\lambda _{1,2}[1,1]`$. The ground state is $`\psi _0=1`$, the first excited state is doubly degenerate, $`\psi _{11}=\lambda _1`$, $`\psi _{12}=\lambda _2`$. Thus, the gap is the same as for magnetic impurities,

$$E_G(1/\tau _{SO}>\mathrm{\Delta }_C)=4/L_{CU}.$$
(13)
An external magnetic field lifts this degeneracy but does not change the gap.

Using the crossover in energy level statistics as the definition of a localization length as above, we get in a quasi- 1 -dim. wire,

$$\xi _c=1/E_G(\beta )=(1/16)\beta L_{CU},$$
(14)
where $`\beta =1,2,4`$ corresponding to no magnetic field, finite magnetic field, and strong spin- orbit scattering or magnetic impurities, respectively. Comparing with the known equation for the localization length, $`L_c`$, we find that the dependence of the ratios $`\beta `$ on the symmetry are in perfect agreement with the result as obtained from the spatial decay of the density- density- correlation function, while it defers by the overall constant $`1/8`$.

This relation can be proven directly. The ASD at zero frequency $`\overline{C}(0)_L`$ of the wire of length $`L`$, becomes, when the wire is divided into two parts, $`\overline{C}(0)_{L/2}^2`$. For $`L\mathrm{}`$, we find that the relative difference is:

$$f(L)=\frac{\overline{C}(0)_{L/2}^2}{\overline{C}(0)_L}1=2\mathrm{exp}(LE_G/2),$$
(15)
exponentially decaying with the length $`L`$. $`f(L)`$ can be estimated, following an argument by Mott: When the two halves of the wire get connected, the Eigenstates of the two separate halves become hybridized and the Eigenenergy of a state $`\psi _n`$ is changed by $`\pm \mathrm{\Delta }_C\mathrm{exp}(2x_n/L_C)`$. $`x_n`$ is random, depending on the position of an eigenstate with closest energy in the other half of the wire. Thus, averaging over $`x_n`$ gives:

$$f(L)+\mathrm{exp}(4L/L_C).$$
(16)
Comparison with Eq. (15) yields indeed $`1/L_C=8E_G`$.

The Crossover Behaviour of the Localization Length

This direct relation of the ASD to the spatial overlap between wavefunctions allows us to extend this approach to get an analytical solution for the crossover behaviour of the localization length and the local level spacing as a magnetic field is turned on, and without strong spin- orbit scattering. While a heuristic approach and numerical studies seemed to indicate a continous increase of the localization length, the analytical result  does show that both limiting localization lengths $`L_c(\beta =1)`$ and $`L_c(\beta =2)`$ are present in the crossover regime and that there is no single parameter scaling. This is explained by arguing that the far tails of the wavefunctions do cover a large enough area to have fully broken time reversal symmetry, decaying with the length scale $`L_c(\beta =2)`$ even if the magnetic field is too weak to affect the properties of the bulk of the wavefunction, which does decay at smaller length scales with the shorter localization length $`L_c(\beta =1)`$, corresponding to the time reversal symmetric case. The quantity studied there is the imprurity averaged correlation function of local wavefunction amplitudes and its momenta at a fixed energy $`ϵ`$: $`Y(ϵ)=<_\alpha \psi _\alpha (0)^2\psi _\alpha (r)^2\delta (ϵϵ_\alpha )>`$. It is averaged over a distribution of eigenfunctions in different impurity representations. Thus, each eigenfunction could decay exponentially with a single localization length, but with a distribution which has two maxima, at $`L_c(\beta =1)`$ and $`L_c(\beta =2)`$, whose weight is a function of the magnetic field in the crossover regime. It is interesting to ask if this property of the eigenfunctions does have a consequence on the energy level statistics as well, or, if that is still governed by a single parameter as the magnetic field is varied.

The effective Hamiltonian for moderate magnetic fields is given, without spin dependent scattering, $`\alpha =1`$, by:

$$\overline{H}=\frac{1}{L_{CU}}(4\mathrm{\Delta }_Q+X^2(1\lambda _C^2)),$$
(17)
where the Laplacian is Eq. and $`X=2\pi \varphi /\varphi _0`$, where $`\varphi `$ is the magnetic flux through an area $`A=L_{CU}(\overline{y^2})^(1/2)`$.

In the limit $`X0`$ the ground state and first excited state approach $`1,\lambda _C\lambda _D`$, respectively. In the limit $`X1`$, $`\lambda _C^2`$ becomes fixed to 1. Thus, the Ansatz $`\psi _0(\lambda _C)=\mathrm{exp}(A_0(1\lambda _C^2))`$, and $`\psi _1(\lambda _C,\lambda _D)=\lambda _C\lambda _D\mathrm{exp}(A_1(1\lambda _C^2))`$, solves $`\overline{H}\psi =\overline{E}\psi `$ to first order in $`z=X^2(1\lambda _C^2)`$ and yields the gap:

$$E_G(X)=4(2+\sqrt{49+X^2}\sqrt{25+X^2})/L_{CU}.$$
(18)
This solution is valid in both the limits $`X1`$ and $`X1`$, interpolating the region $`X1`$. The resulting ratio of local energy level spacings $`\mathrm{\Delta }_C(B)/\mathrm{\Delta }_C(0)=E_G(B)/E_G(0)`$, is shown in Fig. 5 to be in excellent qualitative agreement with experimental data for the magnetic field dependent activation energy, measured recently in transport experiments.

There is a quantitative discrepancy by a factor $`.6`$ between the best fit $`X=.036H/Oe`$, and $`X=2\pi \varphi /\varphi _0`$, $`\varphi =\mu _0HL_{CU}(\overline{y^2})^{(1/2)}`$. With the experimental parameters $`\alpha =1,L_{CO}=.61\mu m`$, width $`W=.2\mu m`$ of sample 5 in Ref.  and $`\overline{y^2}=W^2/12`$ for a 2- dimensional wire, it yields $`X=.021H/Oe`$. We note that smooth confinement can give $`\overline{y^2}>W^2/12`$, thus explaining this discrepancy and the observed difference between $`W`$ as obtained from the sample resistance and estimated from the analysis of the weak localization magnetoresistance, which also depends on $`\overline{y^2}`$. When $`W^2B>\varphi _0`$ or $`H>H^{}=2.110^{11}Oe[m^2/W^2]=525Oe`$, the derivation has to be extended to the 2- dimensional nonlinear sigma- model, and Eq. 18 seizes to be valid. Then, $`L_C`$ increases as function of magnetic field more strongly.

The asymptotic behaviour, of $`\delta L_C(B)B^2`$ for small and $`1/B`$ at large magnetic fields of Eq. (18) does agree with the suggestion by Bouchaud.

The two length scale physics of Ref.  has thus no consequence for the energy level statistics as studied with the ASD. Since we could also show a direct relation between this spectral statistics and the exponential decay of localized eigenfunctions, it is suggestive to be the non self averaging property of the ASD which washes out the two scale physics.

The author gratefully acknowledges, usefull comments by Isa Zarekeshev, discussions with Bernhard Kramer and support from DFG. He thanks Yuri Khavin for usefull communications and Martin Zirnbauer for pointing out a serious confusion in the notation.