format: text

text:
# Galaxy Clustering and Galaxy Bias in a ΛCDM Universe

## 1 Introduction

The clustering of galaxies has long been an essential testing ground for cosmological models and the theory of galaxy formation, with comparisons between predicted and observed clustering driving much of the progress in cosmology during the 1970s and 1980s. Such comparisons provided early support for gravitational instability as the central mechanism of structure formation (e.g., Davis, Groth, & Peebles, 1977), for initial conditions with a power spectrum redder than white noise (e.g., Gott & Rees, 1975; Gott & Turner, 1977), for a “bottom up” rather than “top down” sequence of gravitational clustering (e.g., White, Frenk, & Davis, 1983; Davis et al., 1985; Fry & Melott, 1985), for an approximately Gaussian distribution of primordial fluctuations (e.g., Gott et al., 1989; Weinberg & Cole, 1992; Bouchet et al., 1993), and for a power spectrum significantly redder than that predicted by the “standard” ($`\mathrm{\Omega }_m=1`$, $`hH_0/100\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1=0.5`$) cold dark matter scenario (e.g., Efstathiou et al., 1990; Maddox et al., 1990; Park et al., 1994). The growth of galaxy redshift surveys has led to measurements of increasing precision and detail, and cosmological N-body simulations have developed into a powerful tool for calculating the gravitational clustering of collisionless dark matter from specified initial conditions. The main obstacle to drawing stronger inferences from the data is the dependence of theoretical predictions on both the relatively straightforward physics of gravitational clustering, which largely determines the distribution of dark matter, and the more complex physics of galaxy formation, which determines the relation between galaxies and mass, often referred to generically as “bias.”

In this paper, we study galaxy clustering and galaxy bias in a large ($`2\times 144^3`$ particles, $`50h^1\mathrm{Mpc}`$ cube), smoothed particle hydrodynamics (SPH) simulation of a low density, inflationary cold dark matter universe with a cosmological constant ($`\mathrm{\Lambda }`$CDM). Our goals are, first, to see whether the currently leading cosmological model and “standard” galaxy formation physics as incarnated in our SPH simulation can reproduce existing observations, and, second, to give guidance for the physical interpretation of those observations and for the design of more demanding tests of the theoretical predictions. The N-body+hydro approach to ab initio predictions of galaxy bias has a substantial history, beginning with dissipative “sticky particle” simulations by Carlberg & Couchman (1989) and continuing in the early 1990s with simulations using SPH (Katz, Hernquist, & Weinberg, 1992; Evrard, Summers, & Davis, 1994) or Eulerian grid hydrodynamics (Cen & Ostriker, 1992) to model galaxy formation in a representative (but small) cosmological volume. The present paper is a direct descendant of Katz et al. (1992), but we study a different cosmological model, we incorporate improvements to the physical treatment of radiative cooling and star formation (Katz, Weinberg & Hernquist, 1996, hereafter KWH), and, above all, we take advantage of the parallel implementation of TreeSPH (Davé, Dubinski, & Hernquist, 1997) and advances in computer technology to simulate a volume more than 90 times larger than that of Katz et al. (1992), at similar resolution. Our analysis overlaps significantly with other recent studies using large volume SPH (Pearce et al., 1999, 2001; Yoshikawa et al., 2001) or Eulerian hydrodynamics (Cen & Ostriker, 2000) simulations. Relative to these investigations, we have higher mass and/or spatial resolution and a somewhat smaller simulation volume, as discussed in §2 below.

Two other approaches to ab initio predictions of bias have gained currency in recent years: high-resolution, collisionless N-body simulations that identify galaxies with “subhalos” in the dark matter distribution (e.g., Colín et al., 1999; Kravtsov & Klypin, 1999), and a hybrid method that combines N-body simulations of the dark matter component with semi-analytic treatments of the galaxy formation physics (e.g., Kauffmann, Nusser, & Steinmetz, 1997; Governato et al., 1998; Kauffman et al., 1999a, b; Benson et al., 2000a, b, 2001; Hatton et al., 2002). We will discuss comparisons of our results to those from other hydrodynamic simulations and from the high-resolution N-body and hybrid approaches as they arise. We carry out a detailed comparison between our simulation and the semi-analytic model of Benson et al. (2000a) in a separate paper that focuses on the “halo occupation distribution” (HOD) predicted by the two methods (Berlind et al., 2002). The HOD description can be used to calculate many different galaxy clustering statistics, and it helps explain the origin of bias in a physically intuitive manner (see Berlind & Weinberg 2002 and references therein). Based on this comparison, we conclude that our SPH approach and Benson et al.’s (2000a) semi-analytic method should yield similar predictions for most galaxy clustering statistics.

Here we focus mainly on the “classic” measurements of galaxy clustering — the two-point correlation function, variance of galaxy counts, and pairwise velocity moments — and on one of the simplest measures of higher-order clustering, the third moment of counts-in-cells. Large redshift surveys like the 2dF Galaxy Redshift Survey (2dFGRS) and the Sloan Digital Sky Survey (SDSS) now allow precise measurements of these quantities for multiple classes of galaxies, defined by luminosity, color, morphology, or spectral type (see, e.g., Norberg et al. 2001, 2002b; Zehavi et al. 2002, and references therein). While we cannot model these variations with galaxy type in detail — our simulation volume is far smaller than these surveys, and we do not resolve the morphology of the simulated galaxies — we can examine the predicted trends of clustering with baryon mass, stellar population age, and local environment. We concentrate on present-day clustering, but we also compute the evolution of the two-point function, which can be compared to results from deep redshift surveys and to studies of Lyman-break galaxies at $`z3`$ (e.g., Adelberger et al., 1998, 2002). The discovery that the clustering of Lyman-break galaxies is similar to that of $`L_{}`$ galaxies at $`z=0`$, despite the weaker expected clustering of the underlying mass distribution, provides strong evidence that the bright galaxy population was highly biased at $`z3`$, even if galaxies roughly trace mass today.

Galaxy-galaxy weak lensing is an important new observational probe of the relation between galaxies and dark matter, measuring the galaxy-mass cross-correlation function and the extended mass distributions around galaxies of different types and in different environments. We will devote considerable attention to modeling recent observations of this phenomenon. The closest similar efforts are those of White, Hernquist, & Springel (2001), who present predictions from a hydrodynamic simulation at $`z=1`$ and $`z=0.5`$, and Guzik & Seljak (2001) and Yang et al. (2002), who analyze the N-body+semi-analytic galaxy distribution of Kauffman et al. (1999a).

One strength of the hydrodynamic simulation approach is that it predicts properties of the intergalactic medium (IGM) in addition to the galaxy population. We have used the simulation analyzed here to model the phase distribution of baryons in the present day universe (Davé et al., 2001), the X-ray background (Croft et al., 2001), X-ray emission from galaxy groups (Davé, Katz, & Weinberg, 2002), and X-ray absorption by the diffuse IGM (Chen et al., 2002). In combination with other simulations, we have used it to study cooling radiation and sub-millimeter emission from young galaxies (Fardal et al. 2001, 2002), the relative importance of mergers and smooth accretion in galaxy assembly (Murali et al., 2002), and the correlations between galaxies and the Ly$`\alpha `$ forest at high redshift (Kollmeier et al., 2002). The observed distributions and correlations of galaxy properties, in particular the luminosity function and the Tully-Fisher (1977) relation, are also essential tests of the galaxy formation model, and any discrepancies with these observations can provide insight into the model’s failings. We will consider these characteristics of the galaxy population in a separate paper (Katz et al., in preparation), though the comparison to observed galaxy-galaxy lensing measurements will bring up some of the same issues here. We proceed to a description of our numerical methods in §2, to predictions of galaxy clustering in §3, and to the galaxy-mass correlation function and comparison to lensing measurements in §4. We recap our findings and discuss their implications in §5.

## 2 Simulating Galaxy Formation

### 2.1 Simulation

We employ a hydrodynamic simulation of a $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_m=0.4`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.6`$, $`h=0.65`$, $`\mathrm{\Omega }_b=0.02h^2=0.0473`$, inflationary spectral index $`n=0.95`$, and a power spectrum normalization $`\sigma _8=0.80`$. This model is consistent with most available observational constraints (Bahcall et al., 1999; Jaffe et al., 2001), including COBE normalization (Bennett et al., 1996), cluster masses (Eke, Cole, & Frenk, 1996), cosmic shear (Hoekstra, Yee, & Gladders 2002, and references therein), the deuterium abundance in high-redshift Lyman limit systems (Tytler, Fan, & Burles, 1996; Burles & Tytler, 1997, 1998), the Hubble diagram of Type Ia supernovae (Riess et al., 1998; Perlmutter et al., 1999), and the flux power spectrum of the Ly$`\alpha `$ forest (Croft et al., 1999, 2002; McDonald et al., 2000).

Our initial conditions represent a random $`50h^1\mathrm{Mpc}`$ (comoving) cube, realized with $`144^3`$ gas and $`144^3`$ dark matter particles in a triply periodic volume, yielding particle masses of $`m_{\mathrm{gas}}=8.5\times 10^8M_{}`$ and $`m_{\mathrm{dark}}=6.3\times 10^9M_{}`$, respectively. These initial conditions are evolved from redshift $`z=490`$ using Parallel TreeSPH (Davé, Dubinski, & Hernquist, 1997). We soften the gravitational force via a $`10h^1\mathrm{kpc}`$ (comoving) spline kernel, equivalent to a Plummer softening of $`7h^1\mathrm{kpc}`$. For particle timesteps we use the criteria described by KWH and Quinn et al. (1997), setting the tolerance parameter $`\eta `$ (defined in these papers) to 0.4. The simulation was evolved to $`z=0.5`$ on the Cray T3E at SDSC and finished on the SGI Origin 2000 at NCSA.

We include radiative cooling for primordial composition gas as described by KWH. The tests in Weinberg et al. (1997) show that photoheating by an ionizing background artificially suppresses the formation of low mass galaxies in a simulation of this resolution, but that a simulation evolved without an ionizing background produces much the same galaxy population (above its resolution limit) as a higher resolution simulation with an ionizing background. Though Weinberg et al. (1997) considered an $`\mathrm{\Omega }_m=1`$ CDM model rather than $`\mathrm{\Lambda }`$CDM, we expect their conclusions to carry over to this model, and we therefore evolve the simulation with no ionizing background.

We incorporate star formation using the algorithm described and tested by KWH. This algorithm converts cold, dense gas into collisionless stars at a rate governed by the local dynamical and cooling timescales, returning supernova feedback as thermal energy to the surrounding medium. The tests in KWH and Weinberg et al. (1997) show that the resulting galaxy population is insensitive to the one free parameter of the algorithm and is similar to that identified using cold gas in a simulation without star formation; it is gravity, gas dynamics, and radiative cooling that determine where galaxies form and how massive they will be. A radical change to the star formation or feedback algorithm might alter the galaxy population in a significant way, but provided that it gave the same relative masses and ages of the galaxies (more precisely, that it preserved their rank order), it would not affect the clustering results presented here, except for the absolute values of mass-to-light ratios discussed in §4.

Our numerical approach is similar to those used by Pearce et al. (1999, 2001) and Yoshikawa et al. (2001), but we have opted for a somewhat smaller volume and larger particle number, and consequently higher resolution. Pearce et al. use $`128^3`$ particles in a $`70h^1\mathrm{Mpc}`$ cube, so our mass resolution is higher by $`(70/50)^3(144/128)^3=3.9`$ (though the ratio of SPH particle masses is only 2.9 because Pearce et al. use a lower value of $`\mathrm{\Omega }_b`$). The gravitational softening of the Pearce et al. simulation is similar to ours at $`z=0`$ ($`10h^1\mathrm{kpc}`$ vs. our $`7h^1\mathrm{kpc}`$), but larger at higher redshift, since they keep the softening fixed in physical coordinates at $`z<2.5`$. The mass resolution of the Yoshikawa et al. (2001) simulation is similar to that of Pearce et al. ($`128^3`$ particles in a $`75h^1\mathrm{Mpc}`$ cube), and the (Plummer equivalent) gravitational softening is $`41h^1\mathrm{kpc}`$, constant in comoving coordinates. The most important methodological difference is that we convert cold gas into collisionless stars, while Pearce et al. and Yoshikawa et al. “decouple” cold gas from hot gas when computing SPH densities but leave the cold gas subject to SPH forces and dissipation. Despite the differences in numerical resolution, input physics, and implementation details, we find quite good agreement with these previous investigations where our results overlap, as discussed in §3. The cosmological parameters of the three $`\mathrm{\Lambda }`$CDM simulations are similar but not identical. (Pearce et al. investigate an $`\mathrm{\Omega }_m=1`$ model in addition to $`\mathrm{\Lambda }`$CDM.)

Cen & Ostriker (2000) have investigated galaxy clustering and bias in the $`\mathrm{\Lambda }`$CDM model with a $`512^3`$ Eulerian mesh simulation of a $`100h^1\mathrm{Mpc}`$ volume. The initial mass resolution of this simulation is higher than ours, but the force resolution is much lower, with individual mesh cells $`200h^1\mathrm{kpc}`$. Because the simulation does not track individual galaxies within high density regions, Cen & Ostriker examine stellar mass weighted clustering statistics rather than the galaxy number weighted statistics that we investigate here. Nonetheless, our conclusions about galaxy bias on the scales $`r1h^1\mathrm{Mpc}`$ resolved by the Cen & Ostriker simulation are qualitatively similar, as discussed in §3.

### 2.2 Galaxy Identification

Figure 1 illustrates the state of the simulation at $`z=0`$, showing dark matter, gas, and star particles projected through the entire volume in a quarter of the box, and in expanded views 6, 1.5, and 0.375$`h^1\mathrm{Mpc}`$ on a side centered on the most massive system in our volume. The star particles are clumped into groups of extremely high overdensity, from a few kpc to a few tens of kpc across. Dark halos with the mass scale of galaxy groups and clusters contain many such dense clumps in addition to a smooth distribution of hot gas. Figure 1 shows that there is virtually no ambiguity in identifying the galaxies in this simulation; one simply needs an automated algorithm that picks out the distinct clumps of stars and cold gas. The algorithm that we use is Spline Kernel Interpolative DENMAX 111http://www-hpcc.astro.washington.edu/tools/SKID/ (SKID), which selects gravitationally bound groups of particles that are associated with a common local density maximum (see KWH). SKID selects essentially the same population of objects that one would pick out by eye from a representation like Figure 1.

Figure 2 shows the baryonic mass function of the simulated galaxy population at $`z=3`$, 1, and 0, where $`M_b`$ represents the total mass of a galaxy’s stars and associated cold ($`T<30,000`$K) gas. The number of galaxies steadily increases with time, though the rate of new galaxies forming is dramatically slower for $`z=10`$ than for $`z=31`$. We will discuss the $`z=0`$ luminosity function and Tully-Fisher relation in a separate paper (Katz et al., in preparation), and we have elsewhere compared predictions of the $`z=3`$ luminosity function from higher resolution simulations (of smaller volumes) to measurements of the Lyman-break galaxy population (Weinberg et al. 1999, 2002). The dashed vertical line in Figure 2 marks the mass of 60 SPH particles, and all of the simulated mass functions turn over shortly below this threshold. Our tests using a suite of simulations of different resolution show that galaxies below this threshold have underestimated masses or may be missed entirely, while the masses and locations of galaxies above this threshold are generally robust (see Weinberg et al. 1999). The dotted line shows the mass function at $`z=0`$ from a higher resolution ($`ϵ=3h^1\mathrm{kpc}`$, $`m_{\mathrm{gas}}=1.04\times 10^8M_{}`$), smaller volume ($`22.22h^1\mathrm{Mpc}`$ cube) simulation having $`2\times 128^3`$ particles and the same cosmological parameters, analyzed in the same way; a comparison with the solid line clearly shows that the turnover below the demarcated mass limit is due to numerical effects. We therefore adopt a baryon mass of $`M_b=60m_{\mathrm{gas}}=5\times 10^{10}M_{}`$ as our completeness limit, and we analyze the clustering only of galaxies above this threshold. Note that the mass function comparison also suggests that the lower resolution simulation may overestimate the baryon masses of the largest galaxies, a point that we will return to in §4.

Out of 4632 total SKID-identified galaxies at $`z=0`$, 2571 satisfy our $`M_b60m_{\mathrm{gas}}`$ completeness criterion. At $`z=1,2,3`$ and 4, the corresponding fractions are 2011/4012, 1034/2488, 420/1177, and 104/384. The sample of all galaxies with $`M_b60m_{\mathrm{gas}}`$ will be referred to as our “complete sample” of galaxies. At $`z=0`$, the space density of the complete sample, $`2571/(50h^1\mathrm{Mpc})^3=0.02h^3\mathrm{Mpc}^3`$, corresponds to that of observed galaxies brighter than about $`L_{}/4`$ (Blanton et al., 2001), where $`L_{}`$ is the characteristic luminosity of the Schechter (1976) function fit.

## 3 Galaxy Clustering

Figure 3 illustrates the evolution of structure in a $`10h^1\mathrm{Mpc}`$ thick slice of the simulation volume, running from $`z=3`$ (top) to $`z=1`$ (middle) to $`z=0`$ (bottom). The left hand panels show the dark matter distribution, which exhibits the usual behavior seen in N-body simulations: a steady increase of contrast between growing filaments and emptying voids, and fragmentation of filaments into increasingly massive clumps. The middle column shows the clustering of a subset of galaxies with nearly constant comoving space density, namely the 420 resolved galaxies at $`z=3`$ and the 500 most massive galaxies at $`z=1`$ and $`z=0`$. We will refer to this subset below as our “$`L_{}`$ sample.” The evolution of galaxy clustering is markedly different from that of the dark matter. The network of filaments and voids is already present at $`z=3`$, with nearly the full contrast that it achieves by $`z=0`$. The principal change in the clustering of this sample is on small scales, with the collapse of some extended structures into tight clumps.

The right column of Figure 3 shows the full resolved galaxy population at the three redshifts, with a point size that reflects the galaxy baryon mass. Here one can see the increase in the number of resolved galaxies (from 420 to 2011 to 2571), the growth in the average mass of these galaxies, and the preferential location of the most massive galaxies in dense systems at $`z=0`$. The changing number density and size of points makes it difficult to gauge the evolution of clustering strength visually, but we will see below that the correlation function of this complete galaxy sample stays nearly constant from $`z=3`$ to $`z=0`$.

### 3.1 The Two-Point Correlation Function

The two-point correlation function is the most well studied measure of galaxy clustering, in part because it is relatively straightforward to estimate from angular or redshift survey data, in part because it contains complete statistical information on fluctuations in the linear regime if these fluctuations are Gaussian, and, perhaps most of all, because its observed form is remarkably simple. For most galaxy samples, the correlation function is well described by a power law,

$$\xi (r)=\left(\frac{r}{r_0}\right)^\gamma ,$$
(1)
on scales $`r0.0110h^1\mathrm{Mpc}`$ (Totsuji & Kihara, 1969; Peebles, 1974; Gott & Turner, 1979). However, the values of $`r_0`$ and, to a lesser extent, $`\gamma `$, depend on the luminosity, color, and type of the sample galaxies (e.g., Guzzo et al. 1997; Norberg et al. 2001, 2002b; Zehavi et al. 2002, and numerous references therein). We list determinations of $`r_0`$ and $`\gamma `$ from a number of recent redshift surveys in Table 1. While the varying values of $`r_0`$ and $`\gamma `$ partly reflect statistical uncertainties, they arise mainly from the different selection criteria of the samples from which they are measured. In most cases, the values quoted in Table 1 are derived from flux-limited galaxy samples, and the weighting of galaxies of different luminosities depends on the specific procedure adopted for estimating $`\xi (r)`$. The quoted values for the 2dFGRS are derived from a volume-limited sample of galaxies with absolute magnitudes close to the characteristic magnitude $`M_{b_\mathrm{J}}^{}`$ obtained for 2dFGRS galaxies. The quoted values for the SDSS are from a flux-limited sample, but they are similar to those obtained for a volume-limited sample with absolute magnitudes close to $`M_r^{}`$. The difference between the 2dFGRS and SDSS correlation lengths, $`4.9\pm 0.3h^1\mathrm{Mpc}`$ vs. $`6.1\pm 0.2h^1\mathrm{Mpc}`$ (or $`6.3\pm 0.8h^1\mathrm{Mpc}`$ for the volume-limited sample) presumably reflects the difference between blue and red selection.

Early N-body simulations showed that gravitational clustering produces an approximate power-law $`\xi (r)`$, especially if the initial power spectrum has a slope $`n1`$ (Gott, Turner, & Aarseth, 1979). However, the much higher precision correlation functions calculated with modern simulations show significant departures from a power-law $`\xi (r)`$ for CDM or power-law initial spectra (e.g., Jenkins et al. 1998 and references therein). The origin of these departures can be understood in analytic models of matter clustering (Hamilton et al., 1991; Peacock & Dodds, 1996; Ma & Fry, 2000; Seljak, 2000; Scoccimarro et al., 2001). Thus, the observed power-law form of $`\xi (r)`$ requires either special features in the primordial power spectrum or scale-dependent bias between the galaxy and dark matter correlation functions at the present day.

Figure 4 illustrates the evolution of $`\xi (r)`$ in our simulation from $`z=3`$ to $`z=0`$. The top panel shows results for the complete galaxy sample, i.e., all galaxies above our $`60m_{\mathrm{gas}}=5\times 10^{10}M_{}`$ resolution threshold at the corresponding redshift, corresponding to the right hand column of Figure 3. As noted in §2.2, there are 420, 1034, 2011, and 2571 galaxies satisfying this criterion at $`z=3`$, 2, 1, and 0, respectively. Since the strength of galaxy clustering depends on mass (see §3.2 below), it is useful to identify samples of specified space density for comparison to observations. The middle panel shows results for the 500 most massive galaxies present at each redshift (except at $`z=3`$, where it contains only the 420 galaxies above the resolution threshold). The ranking of each galaxy is based on its total mass of stars and cold ($`T<30,000`$K) gas. The comoving space density of this sample, $`500/(50h^1\mathrm{Mpc})^3=0.004h^3\mathrm{Mpc}^3`$, is similar to that of $`L_{}`$ galaxies today (Blanton et al., 2001) and similar to (but slightly higher than) the space density of Lyman-break galaxies in the spectroscopic samples of Adelberger et al. (1998). We therefore refer to this as our “$`L_{}`$-sample,” illustrated by the middle column of Figure 3. At $`z=0`$, the baryon mass threshold for this sample is $`M_b=2.9\times 10^{11}M_{}`$. The bottom panel of Figure 4 shows the dark matter $`\xi (r)`$.

In each panel of Figure 4, we show as a reference the power-law fit to the SDSS correlation function found by Zehavi et al. (2002). Comparison of our predictions to results from a flux-selected observational sample rests on the implicit assumption that the luminosity of a galaxy is an approximately monotonic function of its mass, and we use the SDSS result in preference to the 2dFGRS result of Norberg et al. (2001) because this assumption seems more applicable to the $`r`$-band selection used in the SDSS. Lyman-break galaxy samples at $`z3`$ are selected in the rest-frame UV, where the instantaneous star formation rate is more important than the stellar mass in determining galaxy luminosity; however, in our simulations, these star formation rates are fairly well correlated (though not perfectly so) with baryon mass (Weinberg et al., 2002).

The main source of statistical error in our predictions is the finite number of coherent structures present in our $`(50h^1\mathrm{Mpc})^3`$ simulation volume. We estimate the error bars on the correlation function and on other statistical quantities that we compute later by applying jackknife resampling to the eight octants of the simulation cube. Specifically, we create eight jackknife subsamples by deleting each of the eight octants of the cube in turn, and we estimate $`\xi (r)`$ for each subsample; to avoid boundary effects that would complicate $`\xi (r)`$ estimation and reduce our effective volume, we identify pairs using the full periodic cube and assign each pair to an octant based on the position of one of its particles. The $`1\sigma `$ error bar $`\sigma _i`$ assigned to the value of $`\xi (r)`$ at $`r=r_i`$ is

$$\sigma _i^2=\frac{N}{N1}\underset{j=1}{\overset{N}{}}[\xi ^j(r_i)\widehat{\xi }(r_i)]^2,$$
(2)
where $`\xi ^j`$ is the estimate in subsample $`j`$, $`\widehat{\xi }`$ is the estimate from the full cube, and $`N=8`$. This approach is similar to using the dispersion among the values of $`\xi (r)`$ from the eight octants (in which case one would divide by $`8^{1/2}`$ to get the error on the mean), but it is more robust because each of the subsamples is closer in size to the full sample. Zehavi et al. (2002) use this method to estimate statistical errors on the observed correlation functions of SDSS galaxies. The jackknife error estimates automatically incorporate the contribution from galaxy shot noise in addition to the finite number of structures, but the latter usually dominates by a large factor over the former (which would be represented by “Poisson error bars”).

The striking difference between the evolution of galaxy and dark matter correlation functions in Figure 4 is similar to that found in previous studies based on hydrodynamic simulations, semi-analytic calculations of galaxy formation, and high-resolution N-body simulations that identify galaxies as “sub-halos” within larger virialized objects (Katz, Hernquist, & Weinberg, 1999; Colín et al., 1999; Kauffman et al., 1999b; Cen & Ostriker, 2000; Benson et al., 2001; Pearce et al., 2001; Somerville et al., 2001; Yoshikawa et al., 2001). As Figure 3 shows, galaxies (at least those massive enough to be resolved by our simulation) form at special locations in the density field, and at high redshift they already trace out the network of filaments nascent in the dark matter distribution. The dark matter correlation function grows in time as mass moves into this network from the surrounding regions, but the structure traced by galaxies stays relatively unchanged, and the galaxy correlation function is only weakly dependent on redshift.

Figure 5 quantifies these points in terms of the parameters $`r_0`$ and $`\gamma `$, which we determine from power-law fits to the correlation functions over the full range plotted in Figure 4, $`0.02h^1\mathrm{Mpc}<r<12h^1\mathrm{Mpc}`$. (Note that our figures plot the point for each radial separation bin at the computed mean separation of pairs in that bin.) For some sparse, high-redshift samples, the smallest radial bins have no pairs, in which we case we start the fit at slightly larger $`r`$. We use the jackknife error bars on the $`\xi (r)`$ points to weight the fit and determine its statistical errors. We do not incorporate the covariance of the $`\xi (r)`$ errors; in principle, we could estimate the full error covariance matrix by the jackknife method, but in practice such estimates are noisy and not easy to use. Treating the $`\xi (r)`$ errors as independent causes us to underestimate the statistical uncertainty in $`r_0`$. The dark matter correlation functions are not well described by a power law at $`z1`$, where the fits have $`\chi ^2/`$d.o.f. of $`1.52`$, so the derived values of $`r_0`$ and $`\gamma `$ are sensitive to the radial range selected for the fit. Power-law fits to the galaxy correlation functions are nearly always adequate, and they usually have $`\chi ^2/`$d.o.f. less than one ($`0.30.7`$) because we have ignored the covariance of the errors in $`\xi (r)`$.

Figure 5 shows steady growth of the comoving correlation length of the dark matter, as expected. The correlation length of the complete galaxy sample, on the other hand, starts at $`r_0=5.2h^1\mathrm{Mpc}`$ at $`z=4`$, declines to a minimum of $`3.0h^1\mathrm{Mpc}`$ at $`z=1.5`$, and climbs again to $`r_0=4h^1\mathrm{Mpc}`$ by $`z=0`$. The $`L_{}`$ sample has a higher correlation length than the complete sample at $`z=0`$, $`r_04.5h^1\mathrm{Mpc}`$, though the correlation lengths of the two samples converge at higher redshift as the $`L_{}`$ threshold approaches our completeness threshold. Correlation lengths for a sample with a higher mass threshold, the 200 most massive galaxies at each redshift, are higher still, $`r_0=5.2h^1\mathrm{Mpc}`$ at $`z=0`$, with only a slight decline towards higher redshift. The values of $`\gamma `$ for the galaxy correlation functions range from 1.75 to 2.05, with a tendency for larger $`\gamma `$ for the more massive galaxy samples or at higher redshifts. At high redshift, the dark matter correlation function is close to a power law, but a quite shallow one. As the scale of non-linear structure grows, it produces a “bulge” in the dark matter $`\xi (r)`$ (see Fig. 4), which drives the correlation function away from a power law and steepens the effective value of $`\gamma `$.

The slow evolution of the galaxy correlation length in Figure 5 agrees with results from the KPNO wide-area survey (Postman et al., 1998), which yields $`r_04.5h^1\mathrm{Mpc}`$ (comoving) at $`z0.5`$, and from the CNOC-2 redshift survey, which shows little evolution of the comoving $`r_0`$ out to $`z1`$ (Carlberg et al., 1999). However, the value of $`r_0`$ depends significantly on galaxy mass and age (see Figure 7 below), so caution is required when interpreting results of surveys that may select different types of galaxies at different redshifts. The predicted correlation length at $`z=3`$, $`r_0=4.2\pm 0.3h^1\mathrm{Mpc}`$, agrees very well with the value $`r_0=4.3\pm 0.3h^1\mathrm{Mpc}`$ measured from spectroscopic samples of Lyman-break galaxies (Adelberger et al. 2002; we have transformed from their $`\mathrm{\Omega }_m=0.3`$ cosmology to our $`\mathrm{\Omega }_m=0.4`$ cosmology using the ratio of angular diameter distances at $`z=3`$), though the slope $`\gamma =1.90\pm 0.07`$ is steeper than the observed value of $`\gamma =1.55\pm 0.15`$. Similar values of $`r_0`$ for Lyman-break galaxies have been obtained using semi-analytic galaxy formation models (Governato et al., 1998; Kauffman et al., 1999b; Benson et al., 2001; Somerville et al., 2001) and other hydrodynamic simulations (Cen & Ostriker, 2000; Pearce et al., 2001; Yoshikawa et al., 2001). The value of $`r_0`$ from our $`50h^1\mathrm{Mpc}`$ simulation is somewhat larger than the estimate $`r_0=3.1h^1\mathrm{Mpc}`$ (including a factor 2.2 correction for finite box size effects) derived by Katz, Hernquist, & Weinberg (1999) from a simulation of an $`11.111h^1\mathrm{Mpc}`$ cube, but a difference of this sort is expected because the mass resolution threshold of our present simulation is a factor of eight higher, and correlation strength increases with galaxy mass.

Perhaps the most significant result of this section is the power-law form of the galaxy correlation function at $`z=0`$. The dark matter correlation function shows the curvature found in previous N-body studies (e.g., Jenkins et al. 1998), and it is inconsistent with a pure power law. We demonstrate this point explicitly in Figure 6, where we plot the ratio of the dark matter, complete galaxy, and $`L_{}`$-galaxy correlation functions to a power law $`(r/4.0h^1\mathrm{Mpc})^{1.75}`$ that has the observed logarithmic slope. The dark matter correlation function rises above this power law at $`r0.5h^1\mathrm{Mpc}`$ and falls well below it at $`r0.05h^1\mathrm{Mpc}`$. (The departure from a power law is somewhat less prominent than that in Jenkins et al.’s  $`\mathrm{\Lambda }`$CDM simulation because we adopt $`\sigma _8=0.8`$ and $`n=0.95`$ instead of $`\sigma _8=0.9`$ and $`n=1`$.) The correlation function of the complete sample, on the other hand, follows the power law almost perfectly. The situation for the $`L_{}`$ sample is less clear; here $`\xi (r)`$ shows a rise at $`r1h^1\mathrm{Mpc}`$ reminiscent of that shown by the dark matter, but within our statistical errors, the $`L_{}`$ $`\xi (r)`$ is adequately described by a power law with $`\gamma =2.00\pm 0.08`$.

The correlation length of the $`L_{}`$ sample, $`r_0=4.5\pm 0.4h^1\mathrm{Mpc}`$, is smaller than the SDSS value of $`r_0=6.1h^1\mathrm{Mpc}`$, and slightly below the value $`r_0=4.9h^1\mathrm{Mpc}`$ found for blue-selected $`L_{}`$ galaxies in the 2dFGRS. The slope $`\gamma =2.00\pm 0.08`$ is steeper than the observed value of $`\gamma 1.75`$. The low amplitude appears to be partly an artifact of the particular realization of $`\mathrm{\Lambda }`$CDM initial conditions in our $`50h^1\mathrm{Mpc}`$ simulation volume. We have run particle-mesh N-body simulations (using the code of Park ) of these initial conditions and of four equivalent volumes with different random realizations, and we find that the dark matter correlation length of this realization is a factor of 1.2 lower than that derived from the average correlation function of the other four realizations ($`3.79h^1\mathrm{Mpc}`$ vs. $`4.60h^1\mathrm{Mpc}`$), simply by chance. There is also a systematic effect of missing power on scales larger than $`50h^1\mathrm{Mpc}`$. With PM simulations of $`100h^1\mathrm{Mpc}`$ volumes, we find a further increase in the mean $`r_0`$ by a factor of 1.1, and a slightly shallower slope ($`\gamma =1.92`$ over the range $`110h^1\mathrm{Mpc}`$, vs. $`\gamma =1.97`$ for the $`50h^1\mathrm{Mpc}`$ realization of the SPH simulation). Unfortunately, we do not know just how these differences in the dark matter $`\xi (r)`$ translate to the galaxy $`\xi (r)`$, but the statistical and systematic effects in our $`50h^1\mathrm{Mpc}`$ volume are large enough that we do not presently regard the quantitative discrepancies between the predicted and observed $`(r_0,\gamma )`$ for $`L_{}`$ galaxies as significant. The statistical uncertainties can be reduced by using the halo occupation distribution derived from the SPH simulation to populate a larger volume N-body simulation, as discussed by Berlind et al. (2002).

Berlind & Weinberg (2002) show that the key requirements for obtaining a power-law $`\xi (r)`$ in a $`\mathrm{\Lambda }`$CDM model are reducing the efficiency of galaxy formation in the most massive (group and cluster scale) halos and suppressing pair counts in the lowest mass halos by keeping fluctuations about the mean occupation number well below the Poisson level (see also Jing, Mo & Börner 1998; Seljak 2000; Peacock & Smith 2000; Scoccimarro et al. 2001). The galaxy population in our simulation satisfies both of these requirements, a point that we discuss in detail elsewhere (Berlind et al., 2002). Our conclusion that including galaxy formation physics removes most of the discrepancy between the predicted correlation function of dark matter and the observed correlation function of galaxies agrees with conclusions derived from other studies based on semi-analytic modeling (Kauffman et al., 1999a; Benson et al., 2000a, b; Somerville et al., 2001), high resolution N-body simulations (Colín et al., 1999), and numerical hydrodynamics (Pearce et al., 1999, 2001; Cen & Ostriker, 2000; Yoshikawa et al., 2001).

### 3.2 Correlation Functions of Sub-populations

As mentioned in §3.1, the galaxy correlation function is known to vary significantly with galaxy luminosity, color, and morphology (e.g., Davis & Geller 1976; Hamilton 1988; Loveday et al. 1995; Guzzo et al. 1997; Willmer, da Costa & Pellegrini 1998). The most comprehensive recent examinations are the studies of luminosity and color dependence for SDSS galaxies by Zehavi et al. (2002) and the studies of luminosity and spectral type dependence for 2dFGRS galaxies by Norberg et al. (2001, 2002b). Zehavi et al. (2002) find a steady increase of correlation amplitude with luminosity, going from $`1.5`$ magnitudes fainter than $`M_{}`$ to $`1.5`$ magnitudes brighter, and they find that red galaxies have a higher amplitude and steeper correlation function than blue galaxies, in agreement with earlier studies of the difference between clustering of early type (red) and late type (blue) samples. Norberg et al. (2001) find a somewhat different luminosity dependence for the $`b_\mathrm{J}`$-selected 2dFGRS galaxies, with little dependence of clustering amplitude on luminosity below $`L_{}`$ but a strong dependence at higher luminosities. Norberg et al. (2002b) find that the correlation amplitude varies strongly with galaxy spectral type, analogous to the usual dependence on morphological type, and that the luminosity dependence is present within each spectral type class. Galaxies with early spectral types have a somewhat steeper correlation function than those with late spectral types, though the difference in slope is not as striking as that for the SDSS red and blue galaxy samples.

Our simulation volume is too small and our computation of galaxy properties too simple to allow detailed comparisons to these observations, but we can examine qualitative trends. The color and spectral type of a galaxy depends largely on the age of its stellar population, and morphology is also known to correlate strongly with population age. For our simulated galaxy population, we define a galaxy’s age by the epoch at which it formed half of its stars. While this approach necessarily brings our heuristic star formation algorithm into play, the details of this algorithm should have little effect on the conclusions presented here, since we will use only the relative ages of the simulated galaxies.

Figure 7a shows the correlation functions of the older (dashed line) and younger (dotted line) halves of the complete galaxy sample, compared to the correlation function of the complete sample (solid line). The correlation amplitude exhibits a clear dependence on galaxy age; the older galaxies have a correlation length of $`r_0=5.2h^1\mathrm{Mpc}`$ and the younger galaxies have $`r_0=2.8h^1\mathrm{Mpc}`$. The fitted slopes of the correlation functions are nearly identical, $`\gamma =1.84`$ and 1.82, respectively. (Values of $`r_0`$ and $`\gamma `$ for all sub-populations appear in Table 2.) These results are in good qualitative agreement with the spectral type dependence found by 2dFGRS; relative to the SDSS results, they reproduce the observed color dependence of correlation amplitude but not the observed change of slope.

Figure 7b compares the correlation functions of the more massive (dashed) and less massive (dotted) halves of the complete sample. While the more massive galaxies have a higher $`\xi (r)`$ at all separations, the difference is marked only below $`0.5h^1\mathrm{Mpc}`$. Power-law fits yield nearly identical correlation lengths, $`r_04.0h^1\mathrm{Mpc}`$, but significantly different slopes ($`\gamma =1.90`$ vs. $`\gamma =1.57`$). Table 1 shows that $`r_0`$ does rise significantly if one takes the 500 most massive ($`r_0=4.5h^1\mathrm{Mpc}`$) or 200 most massive ($`r_0=5.2h^1\mathrm{Mpc}`$) galaxies, and in combination with Figure 7b it suggests that this mass dependence sets in at roughly the space density of $`L_{}`$ galaxies. The 2dFGRS results show a transition of this sort (Norberg et al., 2001), though the SDSS results show a steadier dependence of clustering amplitude on luminosity (Zehavi et al., 2002).

Figure 7c shows the effect of dividing the galaxy sample based on local density, computed by smoothing with a spline kernel (the same form used in the SPH calculation) whose radius is adjusted to always enclose 16 galaxies within the smoothing volume. Since $`\xi (r)`$ is itself a measure of overdensity at radius $`r`$, it is no surprise that this division produces the most marked change in the correlation function. The enhancement of $`\xi (r)`$ persists to scales much larger than the typical smoothing radius because of the bias of locally dense regions towards overdense large scale environments (Kaiser 1984; see Mann, Peacock, & Heavens 1998 and Narayanan, Berlind, & Weinberg 2000 for discussions in the context of N-body models). The observed dependence of galaxy morphology on local density (e.g., Dressler 1980; Postman & Geller 1984) should therefore lead to a large scale bias in the clustering of early-type galaxies.

To summarize the results of this and the previous section, we find that our simulated galaxy population reproduces many of the features of observed galaxy correlation functions: a power-law form with $`\gamma 1.75`$ down to separations $`r0.02h^1\mathrm{Mpc}`$, an increase of correlation strength with mass that becomes more pronounced at high masses, and a higher correlation amplitude for galaxies with older stellar populations. The trends of clustering strength with galaxy mass and age agree with the predictions from semi-analytic models (Kauffman et al., 1999a; Benson et al., 2000a; Somerville et al., 2001) and other hydrodynamic simulations (Cen & Ostriker, 2000; Pearce et al., 2001; Yoshikawa et al., 2001). The predicted correlation lengths are somewhat low and the $`\xi (r)`$ slope for $`L_{}`$ galaxies slightly too steep in comparison to the SDSS results, but our N-body experiments suggest that these discrepancies are largely explained by the truncation of power at the scale of our $`50h^1\mathrm{Mpc}`$ box and by a lower than average clustering amplitude in the particular $`\mathrm{\Lambda }`$CDM initial conditions used in our simulation.

### 3.3 Pairwise Velocities

Peebles (1980), Davis & Peebles (1983), and Bean et al. (1983) describe the techniques now commonly used in correlation function analyses of galaxy redshift surveys. Starting from $`\xi (r_p,\pi )`$, the correlation as a function of projected and line-of-sight separations, one integrates over $`\pi `$ to obtain the projected correlation function $`w(r_p)`$, which can be easily inverted to obtain the real space 3-d correlation function $`\xi (r)`$. The full $`\xi (r_p,\pi )`$ can be calculated, approximately, by a convolution of $`\xi (r)`$ with the galaxy pairwise velocity distribution (see Fisher  for a careful discussion of this approximation and its relation to the linear theory approximation of Kaiser  and Hamilton ). This distribution is specified by its first and second moments $`v_{12}(r)`$ and $`\sigma _{12}^2(r)`$, and by a functional form, usually assumed (on both theoretical and observational grounds) to be exponential. The pairwise dispersion can be inferred from $`\xi (r_p,\pi )`$ in a fairly robust fashion, though it is sensitive to the presence of rare, rich clusters in the survey volume (Mo, Jing, & Börner, 1993). The mean pairwise velocity can in principle be deduced as well, given sufficiently good data (see, for example, the analysis of Fisher et al. 1994); at larger separations, it can also be measured directly in galaxy peculiar velocity surveys (Ferreira et al., 1999).

Figure 8 shows the mean pairwise velocity (top) and pairwise velocity dispersion (bottom) at $`z=0`$ (left) and $`z=1`$ (right). At $`z=1`$, the mean pairwise velocity of the galaxies is higher than that of the dark matter, at all separations $`0.110h^1\mathrm{Mpc}`$. Although galaxies and dark matter feel the same gravitational forces and therefore have similar large scale velocity fields, pair-weighted statistics like $`v_{12}(r)`$ and $`\sigma _{12}(r)`$ are also sensitive to spatial bias; in this case, the preferential location of galaxies in overdense regions with high inflow velocities boosts $`v_{12}(r)`$. The mean pairwise velocities of the $`L_{}`$ sample are higher than those of the complete sample, but the difference is of marginal statistical significance.

By $`z=0`$ (top left), the bias of the galaxy mean pairwise velocities has largely disappeared, in correspondence with the small spatial bias at $`z=0`$ (see Figure 4, and Figure 9 below). The complete and $`L_{}`$ galaxy samples again exhibit similar $`v_{12}(r)`$, with velocities of the $`L_{}`$ sample being slightly higher. There is a much sharper difference, however, between the older and younger halves of the complete sample (circles and squares in Figure 8), with the older galaxies having nearly twice the value of $`v_{12}(r)`$ on scales of one to several $`h^1\mathrm{Mpc}`$. This difference reflects the preference of older galaxies for denser environments with larger inflow velocities.

The pairwise dispersion of the complete galaxy sample is below that of the dark matter both at $`z=1`$ and at $`z=0`$. The $`L_{}`$ sample’s pairwise dispersion is similar to that of the dark matter at $`z=1`$ and below it at $`z=0`$. All of the galaxy samples show a peak in $`\sigma _{12}(r)`$ at $`r0.3h^1\mathrm{Mpc}`$, roughly the half-mass radius of the richest groups. While one might expect suppression of the galaxy pairwise dispersion to go hand-in-hand with spatial anti-bias, our results show that the story is not so simple. For example, the correlation function of the $`L_{}`$ sample at $`z=0`$ is positively biased at all separations (Figure 6), but its pairwise dispersion is lower than that of the dark matter and nearly identical to that of the (more weakly clustered) complete galaxy sample. As emphasized by Berlind & Weinberg (2002) and Sheth et al. (2001), the pairwise dispersion is influenced by factors that do not affect $`\xi (r)`$ (or even affect it with opposite sign), making it sensitive to details of the relation between galaxies and their surrounding dark matter halos. While $`\sigma _{12}(r)`$ is not sensitive to galaxy baryon mass (at least in the comparison presented here), it is highly sensitive to stellar population age: the older half of the complete sample has a substantially higher pairwise dispersion than the younger half because of the location of older galaxies in more massive virialized structures.

The pairwise dispersion of the complete and $`L_{}`$ galaxy samples agrees reasonably well with the values $`\sigma _{12}550650\mathrm{km}\mathrm{s}^1`$ at $`r1h^1\mathrm{Mpc}`$ found in the CfA2, LCRS, and SDSS redshift surveys (Marzke et al., 1995; Jing, Mo & Börner, 1998; Zehavi et al., 2002). However, the sensitivity of $`\sigma _{12}(r)`$ to rare, rich clusters is reflected in the large jackknife error bars of Figure 8, so we cannot take this agreement as a major triumph of the simulation. The scale-dependence in Figure 8 contrasts with the rather flat $`\sigma _{12}(r)`$ inferred from the observations, but for the same reason we do not view this as a serious discrepancy. We can conclude from Figure 8 that the predicted difference between dark matter and galaxy pairwise dispersions makes it easier to understand the relatively low observed value of $`\sigma _{12}`$, and that the predicted difference between old and young galaxies naturally accounts for the factor of $`2`$ difference in $`\sigma _{12}`$ derived for red and blue galaxies in the SDSS (Zehavi et al., 2002) and the factor $`1.7`$ difference in $`\sigma _{12}`$ derived from optical redshift surveys (as cited above) and from the IRAS PSCz catalog (as analyzed by Jing, Börner, & Suto 2002).

Our conclusion that the galaxy pairwise dispersion at $`z=0`$ is significantly lower than that of the dark matter agrees with the results of Pearce et al. (2001), based on a hydrodynamic simulation, and Benson et al. (2000a), based on populating N-body halos according to semi-analytic prescriptions. We also agree fairly well with these calculations in terms of the amplitude of the pairwise dispersion, $`\sigma _{12}500600\mathrm{km}\mathrm{s}^1`$ at $`r1h^1\mathrm{Mpc}`$, and the shape of $`\sigma _{12}(r)`$, though our pairwise dispersion drops somewhat more rapidly at $`r<0.5h^1\mathrm{Mpc}`$. Kauffman et al. (1999a), also using N-body halos and semi-analytic galaxy models, find only a slight difference between the dark matter and galaxy pairwise dispersions, and their predicted amplitude at $`r1h^1\mathrm{Mpc}`$ is $`\sigma _{12}800\mathrm{km}\mathrm{s}^1`$. Benson et al. (2000a) argue that the difference between their results and Kauffmann et al.’s arises because their (Benson et al.’s) semi-analytic model places fewer galaxies in high dispersion halos. Since the halo occupation of our simulation agrees well with that of the Benson et al. (2000a) models (Berlind et al., 2002), we think that the same explanation applies to our results. We agree with Kauffman et al. (1999a) with respect to the mean pairwise velocity (not examined by the other groups), with both calculations showing little difference between the galaxies and the dark matter at $`z=0`$. While age divisions (or color divisions) were not examined in these other studies, the tendency to find older, redder galaxies in more massive halos appears to be a generic prediction of both hydrodynamic and semi-analytic calculations, so we expect the strong age-dependence of $`\sigma _{12}(r)`$ and $`v_{12}(r)`$ to be a fairly generic result.

### 3.4 Moments of Galaxy Counts

Moments of the galaxy count distribution in spherical (or other) cells provide an alternative to correlation functions for characterizing galaxy clustering. The $`n`$-th central moment of the galaxy counts is closely related to the $`n`$-point correlation function averaged over the cell volume. As pointed out by Szapudi & Szalay (1993), for a discrete galaxy distribution it is convenient to work with factorial moments such as $`N(N1)`$ and $`N(N1)(N2)`$ (where $`N`$ is the number of galaxies within a cell or sphere), since these automatically remove Poisson shot noise contributions and can therefore be related to the moments of an underlying continuum distribution in a straightforward way.

The ratio $`b_\sigma `$ of rms galaxy count fluctuations to rms mass fluctuations provides a simple and robust measure of galaxy bias, and the value of $`b_\sigma `$ for spheres of radius $`8h^1\mathrm{Mpc}`$ is the number most often quoted as the “bias factor” of the galaxy distribution. Figure 9 shows $`b_\sigma (8h^1\mathrm{Mpc})`$ and $`b_\sigma (2h^1\mathrm{Mpc})`$ as a function of redshift, with

$$b_\sigma \left(\frac{N(N1)}{N^2}1\right)^{1/2}\times \left(\frac{\delta \rho }{\rho }\right)^2^{1/2}$$
(3)
and the moments evaluated using $`10^6`$ spheres placed randomly in the simulation volume. At $`z=0`$ the bias factors are similar for $`r=2h^1\mathrm{Mpc}`$ and $`8h^1\mathrm{Mpc}`$, with $`b_\sigma (8h^1\mathrm{Mpc})=1.02`$ for the complete galaxy sample, 1.22 for the $`L_{}`$ sample, and 1.40 for the 200 most massive galaxies. Figures 58 show that the simulated galaxies do not trace the mass in detail and that galaxy clustering depends on age and luminosity, but our complete galaxy sample is, in an rms sense, nearly unbiased at $`z=0`$.

As shown earlier in terms of $`\xi (r)`$, the galaxy bias starts high and declines with redshift as the dark matter clustering catches up. For the complete sample, $`b_\sigma (8h^1\mathrm{Mpc})=3.3`$ at $`z=3`$. At redshifts $`z1`$, the bias factor is significantly larger at $`2h^1\mathrm{Mpc}`$ than at $`8h^1\mathrm{Mpc}`$, presumably reflecting the different behavior of $`b_\sigma `$ in the non-linear and linear regimes. Fry (1996) discusses an analytic model in which the galaxy population is born at some redshift $`z_{}`$ with a bias factor $`b_{}`$ but thereafter moves with the same velocity as the dark matter and conserves number (see also Tegmark & Peebles 1998). In this case the bias factor at redshift $`z<z_{}`$ is given by $`b(z)1=(b_{}1)D(z_{})/D(z)`$, where $`D(z)`$ is the linear growth factor at redshift $`z`$. Open circles in Figure 9 show the prediction of this model, normalized to the bias factor $`b_{}=3.1`$ of the 500-galaxy sample at $`z_{}=2.5`$. The actual bias factor of this sample declines significantly faster than the number-conserving model predicts, even though the sample contains 500 galaxies at all redshifts, presumably because mergers drive down the number of pairs in dense regions and galaxies that form at $`z<2.5`$ are born in less biased environments. Somerville et al. (2001) obtain a similar result with semi-analytic modeling.

Discussions of the third moment of the count distribution often focus on the hierarchical ratio $`S_3\delta ^3/\delta ^2^2`$. For spherical cells and a linear power spectrum $`P(k)k^n`$, the dark matter has $`S_3=34/7(3+n)`$ to second order in perturbation theory, independent of sphere radius and power spectrum amplitude and nearly independent of $`\mathrm{\Omega }_m`$ (Juszkiewicz, Bouchet, & Colombi, 1993). If the galaxy density contrast is a local function of the matter density contrast, $`\delta _g=f(\delta )`$, then the relation between $`S_{3g}`$ and $`S_3`$ (representing the galaxy and mass moments, respectively) depends on both the linear bias factor $`b_\sigma `$ and a second bias factor $`b_2`$ that characterizes the second derivative $`f^{\prime \prime }(\delta )`$ at $`\delta =0`$ (Fry & Gaztañaga, 1993; Juszkiewicz et al., 1995). The behavior of $`S_3`$ and the relation between $`S_{3g}`$ and $`S_3`$ becomes more complicated in the fully non-linear regime where second-order perturbation theory breaks down, but it remains the case that $`S_3`$ is only weakly dependent on scale and that the ratio $`S_{3g}/S_3`$ depends on the non-linearity of the relation between galaxy and mass densities.

Figure 10 shows the hierarchical ratio $`S_{3\mathrm{f}}`$ defined in terms of factorial moments,

$$S_{3\mathrm{f}}\left(\frac{N(N1)(N2)}{N^3}1\right)\left(\frac{N(N1)}{N^2}1\right)^2,$$
(4)
at $`z=0`$. For the complete galaxy sample, $`S_{3\mathrm{f}}3.54`$ for $`r=14h^1\mathrm{Mpc}`$, rising to $`S_{3\mathrm{f}}=6.5`$ at $`r=8h^1\mathrm{Mpc}`$. These values are slightly below those of the dark matter distribution. Ratios for the more massive, $`L_{}`$ sample are lower still; these galaxies are more clustered in an rms sense, but they are less “skewed” in the sense measured by $`S_{3\mathrm{f}}`$. As emphasized by Colombi et al. (2000), the statistical and systematic errors in measurements of $`S_3`$ from a volume this size are substantial, and our jackknife error bars are large. Thus, while the predicted value of $`S_{3\mathrm{f}}`$ for the complete galaxy sample agrees well with the value $`S_33.5`$ that Gaztañaga & Frieman (1994) infer by de-projecting the APM angular count moments, we cannot claim this as a major success. (The recent analyses of angular clustering in the SDSS by Gaztañaga  and Szapudi et al.  suggest that the APM numbers should be adjusted upward by $`20\%`$.) However, the relative values of $`S_{3\mathrm{f}}`$ for the different populations should be more robustly predicted, and our simulation indicates that galaxies comparable to those in the complete sample should have a hierarchical ratio slightly below that of the dark matter, and that $`S_{3\mathrm{f}}`$ should be smaller for more massive galaxies. The most striking result in Figure 10 is the strong scale-dependence of $`S_{3\mathrm{f}}`$ for young galaxies, a prediction that should be testable with the SDSS and 2dFGRS surveys in the near future.

## 4 Galaxy-Mass Correlations

Advances in wide-field CCD imaging have opened a new window on the relation between galaxies and mass: weak lensing measurements of the galaxy-mass cross-correlation function. This statistic can be derived from galaxy-galaxy lensing analyses, which measure the mean shear profile around a sample of foreground galaxies using the tangential distortion of background galaxies (Tyson et al., 1984; Brainerd, Blandford, & Smail, 1996; dell’Antonio & Tyson, 1996; Griffiths et al., 1996; Hudson et al., 1998; Fischer et al., 2000; Hoekstra, Yee, & Gladders, 2001; McKay et al., 2001; Smith et al., 2001; Wilson et al., 2001b). Alternatively, one can correlate the mass map derived from the ellipticity correlation function of faint galaxies with the light map of brighter, foreground galaxies in the same images (Wilson, Kaiser, & Luppino, 2001a). In addition, cosmic shear measurements are beginning to yield measurements of the matter correlation function itself (for specified $`\mathrm{\Omega }_m`$), with impressively good agreement among independent surveys (see Mellier et al. 2001; Hoekstra, Yee, & Gladders 2002; and references therein), allowing for the comparison of galaxy and mass fluctuation amplitudes.

Figure 11 shows the galaxy-mass cross correlation function for the complete galaxy sample and for the subsamples divided by age (top), mass (middle), and local density (bottom). We define $`\xi _{\mathrm{gm}}(r)=\delta _g(𝐱)\delta _m(𝐱+𝐫)`$, where $`r=|𝐫|`$ and $`\delta _g=n_g/\overline{n}_g1`$ and $`\delta =\rho /\overline{\rho }1`$ are the galaxy and mass (dark + baryonic) overdensities; we compute it by counting galaxy-particle pairs, with particles weighted by mass. While the form of $`\xi _{\mathrm{gm}}(r)`$ at small $`r`$ is connected to the typical galaxy halo profile, the connection is rather loose because the galaxies are not isolated and the samples contain objects with a range of halo virial radii (points emphasized by Guzik & Seljak (2001), White, Hernquist, & Springel (2001), and Yang et al. (2002)). We use the same format in Figure 11 and Figure 7, and the two figures bear a striking resemblance to each other. To a first approximation, the galaxy-mass correlation function has the same shape as the galaxy-galaxy correlation function, similar amplitude, and similar dependence on galaxy age, baryon mass, and environment. However, the variations of $`\xi _{\mathrm{gm}}(r)`$ for different subsamples are smaller than the variations of $`\xi _{\mathrm{gg}}(r)`$, as one would expect, given that the mass distribution stays the same regardless of the galaxy subsample.

With measurements of $`\xi _{\mathrm{gg}}(r)`$, $`\xi _{\mathrm{gm}}(r)`$, and the mass correlation function $`\xi _{\mathrm{mm}}(r)`$, one can define two “bias functions,” $`b_\xi (r)[\xi _{\mathrm{gg}}(r)/\xi _{\mathrm{mm}}(r)]^{1/2}`$ and $`b_{\xi \times }(r)\xi _{\mathrm{gg}}(r)/\xi _{\mathrm{gm}}(r).`$ The ratio $`r_\times (r)b_\xi (r)/b_{\xi \times }(r)=\xi _{\mathrm{gm}}(r)/[\xi _{\mathrm{gg}}(r)\xi _{\mathrm{mm}}(r)]^{1/2}`$ corresponds to a “correlation coefficient” as defined by, e.g., Tegmark & Bromley (1999), though here we work in terms of correlation functions rather than power spectra. In the case of linear, deterministic bias, where $`\delta _g(𝐱)=b\delta (𝐱)`$ at all positions $`𝐱`$, $`b_\xi `$ and $`b_{\xi \times }`$ would be equal to $`b`$ and independent of scale, and $`r_\times `$ would equal unity. Figure 12 plots $`b_\xi (r)`$, $`b_{\xi \times }(r)`$, and $`r_\times (r)`$ for the complete galaxy sample, the $`L_{}`$ sample, and the older and younger halves of the complete sample. We will focus on the $`z=0`$ results, but because measurements of these quantities come from deep imaging surveys, we also show the simulation predictions at $`z=0.5`$, which are qualitatively similar.

The complete sample has $`b_\xi 1`$, but there is a broad minimum around $`0.5h^1\mathrm{Mpc}`$, with $`b_\xi (0.5h^1\mathrm{Mpc})0.75b_\xi (10h^1\mathrm{Mpc})`$. It is precisely this mild scale-dependence of $`b_\xi (r)`$ that transforms the curved dark matter correlation function into a power-law galaxy correlation function (see Figure 6). The $`L_{}`$ sample has higher $`b_\xi `$, roughly flat (though noisy) at $`r0.5h^1\mathrm{Mpc}`$, and climbing towards small $`r`$. The old galaxies have $`b_\xi (r)`$ similar to that of the $`L_{}`$ sample, while the young galaxies have low $`b_\xi (r)`$ with a scale-dependence similar to that of the complete sample.

The scale-dependence of $`b_{\xi \times }(r)`$ is generally weaker than that of $`b_\xi (r)`$. At $`r1h^1\mathrm{Mpc}`$, values of $`b_{\xi \times }`$ are similar to those of $`b_\xi `$ for the corresponding sample, except in the case of the young galaxies, which have $`b_{\xi \times }>b_\xi `$ and thus $`r_\times (r)<1`$. This “anti-correlation” arises because the young galaxies tend to avoid the higher mass halos, which contribute significantly to $`\xi _{\mathrm{mm}}(r)`$. The complete, $`L_{}`$, and old galaxy samples have $`r_\times (r)1`$ at large scales and $`r_\times (r)>1`$ at small scales. While a correlation coefficient defined in terms of smoothed density fields, e.g. $`\delta _g(𝐱)\delta _m(𝐱)/[\delta _g^2(𝐱)\delta _m^2(𝐱)]^{1/2}`$, must lie in the range $`[1,1]`$, the same is not true for $`r_\times (r)`$, which is defined by ratios of correlation functions. The high values of $`r_\times (r)`$ at small $`r`$ reflect the location of a significant fraction of galaxies (especially in the $`L_{}`$ sample) near the central density maxima of dark matter halos. Galaxy-mass and galaxy-galaxy pairs involving these central galaxies boost $`\xi _{\mathrm{gm}}(r)`$ and $`\xi _{\mathrm{gg}}(r)`$, respectively, relative to the values expected if galaxies traced mass within halos. There is no corresponding boost to $`\xi _{\mathrm{mm}}(r)`$, and since $`r_\times (r)\xi _{\mathrm{gm}}(r)/[\xi _{\mathrm{gg}}(r)]^{1/2}`$, it tends to rise above unity on scales where central galaxies make major contributions to the correlation functions (for related discussions, see Seljak 2000 and Berlind & Weinberg 2002).

Hoekstra, Yee, & Gladders (2001) measure a ratio of galaxy-mass and galaxy-galaxy fluctuations that is directly analogous to $`b_{\xi \times }(r)`$, though it is defined in compensated apertures rather than annuli $`rr+\mathrm{\Delta }r`$. In our terminology, their results imply $`b_{\xi \times }=1.26_{0.12}^{+0.15}`$ for our flat, $`\mathrm{\Omega }_m=0.4`$ cosmology (converted from their $`b/r_\times =1.05`$ for $`\mathrm{\Omega }_m=0.3`$ using their equation 14), with no detectable scale-dependence in the range $`r0.11.5h^1\mathrm{Mpc}`$. This result is intermediate between that of our complete and $`L_{}`$ galaxy samples at $`z=0.5`$. (The median redshift of the Hoekstra et al. lens sample is $`z0.35`$, where the bias should be intermediate between our $`z=0`$ and $`z=0.5`$ values.) While we have not modeled the Hoekstra et al. procedure in detail, our prediction of a mild positive bias with no substantial scale dependence and a bias factor similar to that applying on larger scales appears to be fully compatible with their results. The approximate scale-independence of $`b_{\xi \times }(r)`$ also agrees with Wilson et al.’s (2001a) finding of similar shapes for the galaxy-galaxy and galaxy-mass angular correlation functions. It is difficult for us to compare the amplitude of $`b_{\xi \times }`$ to Wilson et al.’s result, in part because they weight galaxies by flux, and in part because they interpret their measurement in terms of a value of $`\mathrm{\Omega }_m`$ assuming that their (red, luminous) lens galaxy population traces the mass. The existence of positively and negatively biased samples with approximately scale-independent $`b_{\xi \times }(r)`$ in Figure 11 demonstrates that a constant ratio $`\xi _{\mathrm{gg}}(r)/\xi _{\mathrm{gm}}(r)`$ does not, on its own, justify this assumption.

One can also use galaxy-galaxy lensing to infer the average mass within a sphere of radius $`r`$ centered on a sample galaxy, an approach emphasized by McKay et al. (2001). (Note that this measurement makes no attempt to distinguish mass that “belongs” to the central galaxy from mass that “belongs” to its neighbors or to the common halo of its group or cluster; all mass counts, and it may be counted multiple times for different galaxies.) Figure 13a shows the average mass profile from $`25h^1\mathrm{kpc}`$ to $`1h^1\mathrm{Mpc}`$ in five different bins of galaxy baryon mass, each a factor of two in width. We express the mass in terms of an equivalent circular velocity $`v_c=\sqrt{G\mathrm{\Delta }M(r)/r}`$, which would be constant in the case of a singular isothermal sphere, though of course few galaxies would have any tracers in circular rotation at these large radii. We use the excess mass $`\mathrm{\Delta }M`$, subtracting the mean value $`4\pi r^3\overline{\rho }/3`$, because it is the excess mass that is constrained by lensing, but the difference is negligible on the scale of these plots.

For the lowest mass galaxies, $`M_b=(0.51)\times 10^{11}M_{}`$, the mean circular velocity curve is fairly flat at $`v_c200\mathrm{km}\mathrm{s}^1`$ out to $`r100h^1\mathrm{kpc}`$, but it rises at larger radii because these galaxies typically reside in groups whose characteristic circular velocities are larger than those of the galaxies themselves. Conversely, the highest mass galaxies have a mean circular velocity curve that falls all the way out to $`r=1h^1\mathrm{Mpc}`$, because these massive galaxies typically reside at the centers of potential wells whose mass profiles (baryon + dark) are steeper than isothermal. Intermediate mass ranges have intermediate behavior, with falling $`v_c`$ profiles at small $`r`$ changing to flat or rising profiles at $`r100200h^1\mathrm{kpc}`$. Dashed and dotted curves show $`v_c(r)`$ for the older and younger halves, respectively, of the $`(12)\times 10^{11}M_{}`$ sample. At small radii ($`r50h^1\mathrm{kpc}`$), the circular velocity is the same for the two subsamples, but the preferential location of older galaxies in more massive groups manifests itself as a higher $`v_c`$ at larger radii.

Figure 13b plots the ratio of the aperture mass $`\mathrm{\Delta }M(r)`$ to galaxy baryon mass $`M_b`$ (which is approximately equal to the stellar mass), as a function of $`M_b`$. From bottom to top, solid points connected by solid lines show results for $`r=50`$, 100, 260, 500, and $`1000h^1\mathrm{kpc}`$. We use $`260h^1\mathrm{kpc}`$ rather than $`250h^1\mathrm{kpc}`$ to facilitate comparison to McKay et al.’s (2001) analysis of galaxy-galaxy lensing in the SDSS. We average in the order appropriate for such a study, computing the mean value of the excess mass in bins of $`M_b`$, dividing by the mean $`M_b`$ of the bin, and plotting the result for each bin at the value of the mean $`M_b`$. Note that this ratio counts all excess mass within radius $`r`$, including that associated with neighboring galaxies, but divides by the baryon mass of the central galaxy only. For $`r=50h^1\mathrm{kpc}`$ the $`\mathrm{\Delta }M(r)/M_b`$ curve is essentially flat — on average, a galaxy’s baryon mass is proportional to the total mass in a $`50h^1\mathrm{kpc}`$ sphere around it. For larger $`r`$, the ratio rises towards smaller $`M_b`$, which makes sense if the apertures around these galaxies start to include more mass that “belongs” to their richer neighbors. Dotted and dashed lines are computed for the older and younger halves of the galaxies in each $`M_b`$ bin, at $`r=260h^1\mathrm{kpc}`$. For $`M_b5\times 10^{11}M_{}`$, the older galaxies have substantially higher $`\mathrm{\Delta }M(260h^1\mathrm{kpc})`$, by a factor $`1.52`$, analogous to the factor $`2.7`$ difference that McKay et al. (2001) find for elliptical and spiral subsamples. At higher masses, there is little difference in age between the older and younger halves of the sample (essentially all massive galaxies are old), and the systematic difference in $`\mathrm{\Delta }M(260h^1\mathrm{kpc})`$ narrows.

Figure 14 compares our predictions to the results of McKay et al. (2001), who estimate $`\mathrm{\Delta }M(260h^1\mathrm{kpc})`$ by fitting isothermal spheres to their measured shear profiles in four wide bins of galaxy luminosity. They find that $`\mathrm{\Delta }M(260h^1\mathrm{kpc})/L`$ is roughly independent of galaxy luminosity $`L`$ in the Sloan $`g`$, $`r`$, $`i`$, and $`z`$ bands. To convert $`\mathrm{\Delta }M(260h^1\mathrm{kpc})/L`$ to $`\mathrm{\Delta }M(260h^1\mathrm{kpc})/M_b`$ for comparison with our predictions, we draw on the results of Kauffman et al. (2002), who estimate ratios of stellar mass to observed luminosity by fitting population synthesis and dust extinction models to SDSS galaxy spectra and colors. The trend of mean stellar mass-to-light ratio with galaxy luminosity is well described by the relations $`M_{\mathrm{stel}}/L=\mathrm{{\rm Y}}_{}(L/L_{})^\beta `$, with $`(\mathrm{{\rm Y}}_{},\beta )=`$ (2.65,0.375) in $`g`$, (2.29,0.375) in $`r`$, (2.06,0.35) in $`i`$, and (1.57,0.275) in $`z`$, where $`L_{}=1.11`$, 1.51, 2.05, and $`2.58\times 10^{10}h^2L_{}`$ in the four bands, respectively (Blanton et al. 2001, for $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$).222Kauffmann et al.’s Figure 14 plots median $`M_{\mathrm{stel}}/L`$ values and interquartile ranges in bins of absolute magnitude. Since it is the mean $`M_{\mathrm{stel}}/L`$ that matters for our purposes, we have computed these means using the tabulated $`M_{\mathrm{stel}}/L`$ distributions kindly provided by G. Kauffmann, excluding the 10% tails of the distributions so that the mean is not distorted by extreme $`M_{\mathrm{stel}}/L`$ values that may be a result of poor model fits. Our conclusions would not be substantially different if we instead used the median $`M_{\mathrm{stel}}/L`$ values plotted by Kauffman et al. (2002). These relations fit the Kauffman et al. (2002) results to within $`10\%`$ for galaxies within two magnitudes of $`L_{}`$, though they break down outside of this range. We take the values of $`\mathrm{\Delta }M(260h^1\mathrm{kpc})/L`$ from McKay et al.’s Table 2, convert to $`h=0.65`$, and divide by the value of $`M_{\mathrm{stel}}/L`$ at the central luminosity of each bin (which McKay et al. compute using the same galaxy weights that enter the lensing analysis). In Figure 13b, triangles, squares, hexagons, and circles show results derived from the $`g`$, $`r`$, $`i`$, and $`z`$ band mass-to-light ratios, respectively. Results in the four bands are not independent because they rely on the same lensing data, but the generally good agreement among them suggests that there are no major systematic errors in the population synthesis and dust extinction modeling. The data points come in four groups reflecting the four luminosity bins, though the $`g`$ band points do not reach the highest $`M_b`$ values because of the higher $`M_{\mathrm{stel}}/L`$ in $`g`$ band. The solid line in Figure 14 shows the simulation result for all galaxies (same as the middle solid line of Fig. 13b).

The predicted and observationally inferred values of $`\mathrm{\Delta }M(260h^1\mathrm{kpc})/M_b`$ agree well in the lowest mass bin, but at higher masses the observational points lie above the simulation predictions, by a factor $`1.52`$ on average. We have treated our galaxies as entirely stellar, and we have ignored the light “missed” by SDSS Petrosian luminosities (see Strauss et al. 2002), but accounting for these effects would reduce the discrepancy by only $`1020\%`$. We conclude that the baryon masses of our simulated galaxies are on average too high by a factor $`1.52`$, for $`M_b2\times 10^{11}M_{}`$, provided that the Kroupa (2002) stellar IMF assumed in Kauffmann et al.’s (2002) population synthesis modeling is indeed universal across galaxy types and evolutionary stages. Other hints of a similar discrepancy include the high overall fraction of stellar mass in the simulation (Davé et al., 2001) and the relatively high stellar mass-to-light ratios required to match the observed galaxy luminosity function (Katz et al., in preparation).

We believe that these excessive galaxy baryon masses are largely numerical in origin. A comparison of two simulations of the same $`(22.22h^1\mathrm{Mpc})^3`$ volume, one with the same resolution as the $`(50h^1\mathrm{Mpc})^3`$ simulation analyzed here and one with a mass resolution eight times higher, shows that the same galaxies are on average more massive in the lower resolution simulation, by $`10\%`$ at the $`5\times 10^{10}M_{}`$ completeness limit rising to a factor $`2`$ at masses of several$`\times \mathrm{\hspace{0.17em}10}^{11}M_{}`$ (Fardal et al., in preparation; one can see hints of this effect in the mass function comparison of Figure 2). The higher baryon masses at low resolution probably arise from the overestimated cooling rate of hot gas particles in the vicinity of cold, dense clumps (Pearce et al., 1999; Croft et al., 2001). An alternative formulation of SPH that performs better in the presence of contact discontinuities between cold and hot phases indeed leads to lower galaxy baryon masses, typically by a factor $`2`$ (Springel & Hernquist, 2002). Empirically, we find that the relation $`M_{\mathrm{hr}}=5\times 10^{10}M_{}(M_{\mathrm{lr}}/5.4\times 10^{10}M_{})^{0.82}`$ provides a fairly accurate scaling between masses in the “low resolution” and “high resolution” simulations. The dashed curve in Figure 14 shows the effect of applying this empirical scaling. It represents our best guess at the results we would obtain from a $`2\times 288^3`$ particle simulation of the same $`(50h^1\mathrm{Mpc})^3`$ volume. Comparison to still higher resolution simulations at $`z=3`$ suggests that this is sufficient to yield numerically converged baryon masses, but we have not yet been able to carry out such a test at $`z=0`$.

Clearly, there are still systematic uncertainties in the observational results (the weak lensing measurements themselves and the $`M_{\mathrm{stel}}/L`$ scaling) and in the numerical convergence of the simulation predictions. However, the resolution-corrected predictions and the observationally inferred mass ratios agree well at log$`M_b/M_{}10.7`$ and 11.7, with a factor $`2`$ gap at the intermediate mass scale log$`M_b/M_{}11.4`$. The sharp drop between the second and third groups of observational data points is somewhat surprising, at least given the smooth predicted trend, so a full assessment of the comparison should await improved observational estimates, which should be available soon. Nonetheless, the current level of overall agreement does not leave much room for additional astrophysical processes, such as AGN feedback or stronger stellar feedback, to substantially reduce the baryon masses of galaxies with $`M_b5\times 10^{10}M_{}`$. We would probably predict lower values of $`\mathrm{\Delta }M(260h^1\mathrm{kpc})/M_b`$ (and thus require some suppression of baryon mass) if we adopted lower $`\mathrm{\Omega }_m`$ or higher $`\mathrm{\Omega }_b`$, without changing the assumed $`\sigma _8`$ or $`P(k)`$ shape. However, for our adopted cosmological parameters, the lensing comparison suggests that this simulation already incorporates the astrophysics needed to understand the baryon masses of luminous galaxies, within the current numerical and observational uncertainties. The constraints from galaxy-galaxy lensing complement those from the Tully-Fisher relation and the luminosity function, which we will discuss elsewhere (Katz et al., in preparation).

## 5 Summary and Discussion

For the full population of galaxies with $`M_b>5\times 10^{10}M_{}`$, our simulation predicts a bias factor $`b_\sigma (8h^1\mathrm{Mpc})`$ of unity at $`z=0`$. Beneath this rms similarity of fluctuation amplitudes, however, lie numerous differences between the galaxy and dark matter distributions, in terms of present-day structure and the evolution of that structure.

1. The bias defined by rms fluctuations or by correlation functions is strongly redshift dependent. If we choose a galaxy population of fixed comoving space density (above a redshift-dependent baryon mass threshold), then the galaxy correlation function stays approximately constant from $`z=3`$ to $`z=0`$, so that bias is initially large and declines towards unity as the dark matter catches up to the galaxies (see Figure 4).

2. The bias is scale-dependent in the non-linear regime, and this scale dependence itself depends on redshift. At $`z=0`$, the galaxy correlation function is depressed below the dark matter correlation function at $`r0.5h^1\mathrm{Mpc}`$ and boosted above it at $`r0.1h^1\mathrm{Mpc}`$. As a consequence, the galaxy correlation function is well described by a power law (of slope $`\gamma =1.78\pm 0.05`$), while the dark matter correlation function is not (see Figure 6). The dark matter correlation function is shallower at high redshift while the galaxy correlation function slope is roughly constant, so the rms bias factors at $`2h^1\mathrm{Mpc}`$ and $`8h^1\mathrm{Mpc}`$ start to diverge at $`z1`$ (see Figure 9).

3. The correlation amplitude depends on galaxy baryon mass, with only a slight difference between the more and less massive halves of the complete sample (Figure 7b) but significant increases of $`r_0`$ and $`b_\sigma `$ result if we select the 500 or 200 most massive galaxies in the simulation volume (Figures 5 and 9, Table 2).

4. The dependence of the correlation amplitude on stellar population age is more striking than the dependence on baryon mass. The correlation functions of the older and younger halves of the complete sample have almost identical slopes, but their correlation lengths differ by nearly a factor of two (Figure 7, Table 2).

5. At $`z=0`$, the pairwise velocity dispersion $`\sigma _{12}(r)`$ of galaxies is lower than that of the dark matter, by $`20\%`$ (Figure 8). There is a large difference between the pairwise dispersions of old and young galaxies, with the preferential location of old galaxies in dense environments making $`\sigma _{12}(r)`$ about twice as high at $`r2h^1\mathrm{Mpc}`$. The mean pairwise velocity $`v_{12}(r)`$ is similar for galaxies and dark matter at $`z=0`$, but the denser environments of older galaxies again cause a higher $`v_{12}(r)`$. At $`z=1`$, the positive spatial bias of the galaxy population makes the mean pairwise velocity higher than that of the dark matter, though the pairwise dispersion remains somewhat lower.

6. The skewness of galaxy count fluctuations, as quantified by the hierarchical ratio $`S_{3\mathrm{f}}`$, is lower than that of the dark matter by $`20\%`$ (Figure 10). The $`S_{3\mathrm{f}}`$ ratio and its scale dependence vary with galaxy baryon mass and age, in a fairly complex way.

7. The galaxy-mass cross-correlation function depends on galaxy mass and age in much the same way that the galaxy-galaxy correlation function itself does, though the trends are weaker because the mass distribution is the same in all cross-correlations (Figure 11). For each class of galaxies, the bias function defined by the ratio $`\xi _{\mathrm{gg}}(r)/\xi _{\mathrm{gm}}(r)`$ is roughly scale-independent, even when it is significantly different from unity, though it can be somewhat different in amplitude and shape from the analogous bias function $`[\xi _{\mathrm{gg}}(r)/\xi _{\mathrm{mm}}(r)]^{1/2}`$ (Figure 12).

8. The average extended mass distributions around galaxies depend significantly on baryon mass and age (Figure 13). The distribution around a typical high mass galaxy falls more steeply than an isothermal distribution, so that circular velocity curves $`[G\mathrm{\Delta }M(r)/r]^{1/2}`$ decline from $`r=25h^1\mathrm{kpc}`$ out to $`r=1h^1\mathrm{Mpc}`$. For low mass galaxies, on the other hand, the average mass profile is shallower than isothermal at $`r100h^1\mathrm{kpc}`$, where the individual galaxy halos begin to encounter the environment of a typical surrounding group. This effect is much more pronounced for old galaxies than for young galaxies. The ratio of mean aperture mass $`\mathrm{\Delta }M(r)`$ to galaxy baryon mass $`M_b`$ is roughly independent of baryon mass for apertures ranging from $`r=50h^1\mathrm{kpc}`$ to $`r=1h^1\mathrm{Mpc}`$, though the ratio is higher for the galaxies with the lowest $`M_b`$. The ratio is higher for old galaxies than for young galaxies, again demonstrating the preferential location of old galaxies in dense environments.

For the most part, the dark matter clustering of the $`\mathrm{\Lambda }`$CDM model and these biases of the galaxy population lead to good qualitative and, where we have adequate statistics for comparison, quantitative agreement with observed galaxy clustering. The predicted galaxy correlation function is a power law of the observed slope and approximately the observed amplitude (see further discussion below). The redshift-independence of the comoving correlation length agrees with measurements from Lyman-break galaxy surveys at $`z3`$ (Adelberger et al., 1998, 2002) and from deep redshift surveys extending to $`z1`$ (Postman et al., 1998; Carlberg et al., 1999). The predicted dependence of clustering strength on baryon mass and population age agrees with the trends as a function of luminosity, color, and spectral type found in the 2dFGRS and SDSS (Norberg et al., 2001, 2002b; Zehavi et al., 2002). The predicted value of the hierarchical ratio $`S_{3\mathrm{f}}`$ agrees with measurements from angular clustering catalogs (Gaztañaga & Frieman, 1994; Gaztañaga, 2002; Szapudi et al., 2002). The value and scale-independence of the bias inferred from the ratio of galaxy-galaxy and galaxy-mass correlation functions agrees with recent observational estimates (Hoekstra, Yee, & Gladders, 2001). The dependence of the galaxy-mass correlation function on population age parallels the dependence on galaxy type found by McKay et al. (2001), and the aperture mass $`\mathrm{\Delta }M(260h^1\mathrm{kpc})`$ is approximately proportional to galaxy baryon mass, as McKay et al. (2001) also find.

There are three quantitative discrepancies between our predictions and existing measurements, two fairly subtle and one more substantial. First, the power law slope of $`\xi (r)`$ for the $`L_{}`$-galaxy sample (the 500 most massive galaxies in the box at $`z=0`$) is slightly too steep, $`\gamma =2.00\pm 0.08`$ compared to observed values of $`1.71.8`$, and there are hints of a departure from a power law at $`r1h^1\mathrm{Mpc}`$. Second, the correlation length of this sample, $`r_0=4.5\pm 0.4h^1\mathrm{Mpc}`$, is low compared to the value $`r_0=6.3\pm 0.8h^1\mathrm{Mpc}`$ measured for $`r`$-selected galaxies with $`LL_{}`$ in the SDSS (Zehavi et al., 2002), though it is consistent with the value $`r_0=4.9\pm 0.3h^1\mathrm{Mpc}`$ found for $`L_{}`$ $`b_\mathrm{J}`$-selected galaxies in the 2dFGRS (Norberg et al., 2001). Because of the still-limited size of our simulation volume, these discrepancies are of marginal statistical significance, and the N-body tests discussed in §3.1 imply that the low correlation length, at least, is a consequence of the lower than average fluctuation amplitude in the particular realization of $`\mathrm{\Lambda }`$CDM initial conditions used here.

The more serious quantitative discrepancy is the ratio of aperture mass $`\mathrm{\Delta }M(260h^1\mathrm{kpc})`$ to galaxy baryon mass $`M_b`$. For $`M_b(0.51)\times 10^{11}M_{}`$, the simulation prediction agrees with the result derived by combining McKay et al.’s (2001) values of $`\mathrm{\Delta }M/L`$ with Kauffmann et al.’s (2002) stellar mass-to-light ratios, but at higher $`M_b`$ the predicted ratios are lower than the observationally inferred values by factors of $`1.52`$, implying that the galaxy baryon masses in the simulation are too high (Figure 14). As discussed in §4, we believe that this discrepancy is mostly numerical in origin, since our internal resolution tests (Fardal et al., in preparation) and comparisons between different SPH implementations (Springel & Hernquist, 2002) suggest that galaxy baryon masses are indeed overestimated by about this factor in this mass range, relative to a simulation with higher resolution or more accurate treatment of the interface between hot and cold gas phases. A comparison between our best guess at resolution-corrected results and the McKay et al. (2001) data points yields fairly good agreement in the average values of $`\mathrm{\Delta }M(260h^1\mathrm{kpc})/M_b`$, which leaves little room for additional astrophysical processes, such as more aggressive feedback, to substantially suppress gas cooling and star formation in this galaxy mass range. Consistent with this conclusion, Yang et al. (2002) find that the N-body + semi-analytic models of Kauffman et al. (1999a), which do include much stronger stellar feedback, predict $`\mathrm{\Delta }M(260h^1\mathrm{kpc})/L`$ ratios that are too low (by a factor $`2`$) relative to those of McKay et al. (2001). As emphasized by Guzik & Seljak (2002) and Seljak (2002a, b), the constraints from galaxy-galaxy lensing measurements complement those from the galaxy luminosity function and the Tully-Fisher and fundamental plane relations because they suffer different systematic uncertainties and have different sensitivity to astrophysical and cosmological parameters. For example, lowering $`\mathrm{\Omega }_m`$ while keeping other parameters (including the linear power spectrum) fixed would tend to lower $`\mathrm{\Delta }M(r)/L`$ while having relatively little impact on the luminosity function (though galaxy formation, and hence galaxy luminosities, would still be affected by the different timescales and baryon-to-dark matter ratios of the lower $`\mathrm{\Omega }_m`$ model). The weak lensing constraints will become considerably more powerful as the precision and detail of the measurements improves.

Where we have been able to compare, our results generally agree well with those of other hydrodynamic simulations (Pearce et al., 1999, 2001; Cen & Ostriker, 2000; Yoshikawa et al., 2001), high resolution N-body simulations that identify galaxies with sub-halos (Colín et al., 1999; Kravtsov & Klypin, 1999), and combined N-body + semi-analytic models (Kauffmann, Nusser, & Steinmetz, 1997; Governato et al., 1998; Kauffman et al., 1999a, b; Benson et al., 2000a, b, 2001). All of these methods make qualitatively similar predictions for the evolution of the galaxy correlation function, the bias between galaxies and mass at the present day, the dependence of that bias on galaxy type, and the relative amplitude of dark matter and galaxy pairwise dispersions. This agreement suggests that these aspects of galaxy clustering and bias arise from robust aspects of the physics of galaxy formation that all of these methods treat correctly, or at least similarly. Our conclusions about $`S_{3\mathrm{f}}`$ and galaxy-mass correlations are for the most part new, but we do not expect them to be fundamentally different from the predictions of other methods.

Returning to the main goals of the paper, we find that the combination of a $`\mathrm{\Lambda }`$CDM cosmology with standard ideas about galaxy formation is, on the whole, remarkably successful at reproducing observed galaxy clustering and galaxy-galaxy lensing measurements. In terms of physical interpretation, the simulation analysis offers a number of insights. It shows that the detailed relation between the galaxy and dark halo populations can account for the difference between the predicted correlation function of dark matter and the observed power law form of the galaxy correlation function, and between the predicted pairwise velocity dispersion of the dark matter and the observed, lower dispersion of galaxies. It shows that the observed dependence of galaxy clustering on galaxy properties emerges naturally from the dependence of galaxy mass and population age on environment — basically, older and more massive galaxies form in regions that collapse early and are today biased with respect to the overall mass distribution. It shows that values of $`b_\sigma 1`$ and $`b_{\xi \times }1`$ can emerge even when galaxies do not trace mass in detail. The quantitative discrepancies discussed above may also prove instructive, as the observations improve and the uncertainties in the simulations themselves are more thoroughly understood.

Our analysis makes a number of predictions that can be tested with future data, or with future analyses of existing data. Our quantitative comparison to the observed type dependence of galaxy clustering has been limited by our finite simulation volume, which prevents us from predicting the clustering of rare classes of galaxies. However, it would be straightforward to test the quantitative predictions presented in, for example, Figure 7 by taking an observed sample above a luminosity threshold that yields the same space density as our complete sample and dividing it in two based on luminosity or color (or, better still, on the basis of stellar mass and population age, using techniques like those of Kauffman et al. 2002). The predicted evolution of the mean pairwise velocity and pairwise dispersion (Figure 8) can be tested with upcoming large redshift surveys like the DEEP and VIRMOS-VLT programs. The systematic variations of $`S_{3\mathrm{f}}`$ with galaxy luminosity and age (Figure 10) can be tested with the 2dFGRS and SDSS. The most informative comparisons will probably come from improved measurements of galaxy-galaxy lensing and cosmic shear, which will allow detailed tests of the trends predicted in Figures 12 and 13.

On all of these fronts, observational efforts are advancing at a staggering pace. Comparisons between more comprehensive data and improved theoretical predictions over the next few years should tell us whether our understanding of dark matter clustering and the physics of galaxy formation is indeed complete, at least in terms of features that have major quantitative impact, or whether there are still important ingredients missing from our theoretical recipe.

We thank Andreas Berlind, Stephane Colombi, Mark Fardal, Henk Hoekstra, Guinevere Kauffmann, Tim McKay, Yannick Mellier, and Ludo van Waerbeke for valuable discussions on numerous aspects of this study. This work was supported by the NSF and by the NASA Astrophysical Theory Program and Long-Term Space Astrophysics Program. RD acknowledges the support of Hubble Fellowship grant number HST-HF-01128.01-A from the Space Telescope Science Institute, which is operated by AURA, Inc., under NASA contract NAS5-26555. DHW acknowledges the hospitality of the Institute for Advanced Study and the Institut d’Astrophysique de Paris, and the support of the Ambrose Monell Foundation and the French C.N.R.S., during the completion of this work.