text:
# Clusters of Galaxies in the Local Universe

## 1. Introduction

Clusters of galaxies have become one of our most important cosmological probes because they are relatively easy to discover yet have physical properties and abundances that are very sensitive to our model for the formation and evolution of structure in the universe. Of particular interest to us, a cluster sample provides the means to study the high-mass end of the halo multiplicity function, the average number of galaxies in a halo of mass $`M`$, which provides important insight into the process of galaxy formation.

In this paper we have two objectives. First, we will demonstrate the use of a matched filter approach to finding clusters in a redshift survey using both synthetic catalogs and a large sample of galaxies from the 2MASS survey (Skrutskie et al. Skrutskie97 (1997), Jarrett et al. Jarrett00 (2000)). With the synthetic data we can test the algorithm and our ability to extract the input halo multiplicity function from the output cluster catalog. Second, we will determine the halo multiplicity function, $`N(M)`$, from the 2MASS survey. This study is the first phase in a bootstrapping process – based on a synthetic catalog known to have problems matching reality in detail we can calibrate an algorithm which when applied to the real data can supply the parameters for an improved model of the data. Then the improved model can be used to improve the algorithm and so on.

In §2 we describe the synthetic and real 2MASS data. In §3 we describe our version of the matched filter algorithm. In §4 we test the algorithm on the synthetic catalog, focusing on our ability to determine the halo multiplicity function. In §5 we apply the algorithm to the 2MASS sample. Finally, in §6 we discuss the steps which can improve both the synthetic catalog and the algorithm.

## 2. Data: Real and Synthetic

We searched for clusters using galaxies in the 2MASS survey (Skrutskie et al. Skrutskie97 (1997)) and in simulations of the 2MASS survey. We have chosen the 2MASS survey because of its uniform photometry and large areal coverage, which makes it ideal for finding the relatively rare rich clusters in the local neighborhood where a wealth of auxiliary information is available.

### 2.1. 2MASS Galaxies

The 2MASS survey provides J, H, K photometry of galaxies over the full sky to a limiting magnitude of $`K<13.75`$ mag (Jarrett et al. Jarrett00 (2000)). Based on the Schlegel, Finkbeiner & Davis (Schlegel98 (1998)) Galactic extinction model and the available 2MASS catalogs, we selected all galaxies with Galactic extinction-corrected magnitudes of $`\text{K}12.25`$ mag111We use the 20 mag/arcsec2 circular isophotal magnitudes throughout the paper. and Galactic latitude $`|b|>5^{}`$. The final sample contains a total of 90989 galaxies distributed over 91% of the sky. The redshift measurements are 89% complete for $`\text{K}11.25`$ mag but only 36% complete between $`11.25<\text{K}12.25`$. The galaxies with redshifts are not a random sample of all of the galaxies – redshifts are more likely to be of cluster than field galaxies.

With the data corrected for Galactic extinction, the relation between apparent and absolute magnitudes is

$$M_K=K5\mathrm{log}(D_L(z)/r_0)k(z)=K𝒟(z)$$
(1)
where $`r_0=10`$ pc, $`D_L(z)`$ is the luminosity distance to redshift $`z`$ and $`k(z)=6\mathrm{log}(1+z)`$ is the $`k`$-correction. The $`k`$-correction is negative, independent of galaxy type, and valid for $`z<0.25`$ (based on the Worthey (Worthey94 (1994)) models). To simplify our later notation, we introduce $`𝒟(z)=5\mathrm{log}(D_L(z)/r_0)+k(z)`$ as the effective distance modulus to the galaxy. We use an $`\mathrm{\Omega }_{\mathrm{mat}}=1`$ cosmological model for the distances and assume a Hubble constant of $`H_0=100h`$ km/s Mpc. For our local sample the particular cosmological model is unimportant.

Kochanek et al. (Kochanek01 (2001)) derived the K-band luminosity function for a subset of the 2MASS sample at comparable magnitudes. The luminosity function is described by a Schechter function,

$$\varphi (M)=\frac{dn}{dM}=0.4\mathrm{ln}10n_{}\left(\frac{L}{L_{}}\right)^{1+\alpha }\mathrm{exp}(L/L_{})$$
(2)
where $`M_K=M_K2.5\mathrm{log}(L/L_{})`$ and the integrated luminosity function is

$$\mathrm{\Phi }(M)=_{\mathrm{}}^M\varphi (M)𝑑M=n_{}\mathrm{\Gamma }[1+\alpha ,L/L_{}].$$
(3)
The parameters for the global luminosity function are $`n_{}=(1.16\pm 0.10)\times 10^2h^3`$/Mpc3, $`\alpha =1.09\pm 0.06`$ and $`M_K=23.39\pm 0.05`$ using the 2MASS 20 mag/arcsec2 isophotal magnitudes. This luminosity function is consistent with that derived by Cole et al. (Cole01 (2001)). The total luminosity of a galaxy is approximately $`1.20\pm 0.04`$ times the isophotal luminosity.

We should keep in mind that the early-type and late-type galaxies have different luminosity functions, with the early-type galaxies being brighter but less common than the late-type galaxies. The algorithm for finding clusters requires both a field luminosity function $`\varphi _f(M)`$ and a cluster luminosity function $`\varphi _c(M)`$ which will be different because of the morphology-density relation (e.g. Dressler Dressler80 (1980)). In our present study we assume that we can model both using the global luminosity function, $`\varphi _f(M)=\varphi _c(M)=\varphi (M)`$. Once we have catalogs of 2MASS galaxies in different environments we can go back and estimate the environment-dependent luminosity functions. Note that we only need the shape of the luminosity function for the clusters – all expressions using $`\varphi _c(M)`$ will appear in the ratio $`\varphi _c(M)/\mathrm{\Phi }_c(M_K)`$, which eliminates any dependence on $`n_{}`$.

### 2.2. Simulations of 2MASS

We generated a simulated 2MASS catalog based on the approach described in detail in White & Kochanek (WhiKoc (2001); hereafter WK). We used as a starting point the $`z=0`$ output of a $`256^3`$ particle N-body simulation of a $`\mathrm{\Lambda }`$CDM ‘concordance’ cosmology whose predicted mass function agrees well with recent observational estimates (e.g. Pierpaoli, Scott & White PSW (2001)). Halos were identified using the friends-of-friends (FoF) algorithm with a linking length $`b=0.2`$ times the mean inter-particle spacing, and galaxies were ‘assigned’ to halos based on the multiplicity function. The average number of galaxies assigned to a cluster halo with a FoF mass $`M_{FoF}`$ is

$$N_{FoF}(M)=\left(\frac{M_{FoF}}{M_R}\right)+0.7\text{max}[1,\left(\frac{M_{FoF}}{M_B}\right)^{0.9}]$$
(4)
galaxies where $`M_R=2.5\times 10^{12}h^1M_{}`$ and $`M_B=4.0\times 10^{12}h^1M_{}`$ based on the halo occupancy models of Scoccimarro et al. (SSHJ (2001)) which we have modified to better match the observed number densities of galaxies and the expectation that 2MASS galaxies will trace the mass more faithfully than a blue selected sample. We sometimes used the virial mass $`M_{200}`$ defined as the mass interior to a radius $`r_{200}`$ within which the mean density exceeds 200 times the critical density. The virial radius, which scales as $`r_{200}=1.6(M_{200}h/10^{15}M_{})^{1/3}h^1\mathrm{Mpc}`$, is slightly smaller than the average radius of the FoF cluster with $`b=0.2`$, so that on average $`M_{200}=0.66M_{FoF}`$ (median $`M_{200}=0.71M_{FoF}`$). For massive halos the number of galaxies actually assigned to a halo is Poisson distributed for an expectation value of $`N_{FoF}`$. Each halo hosts a ‘central’ galaxy which inherits the halo center of mass position and velocity. Any additional galaxies are randomly assigned the positions and velocities of the dark matter particles identified with the cluster by the FoF algorithm. In this way the satellite galaxies trace the density and velocity field of the dark matter.

One problem with standard halo multiplicity functions is that they are one-dimensional functions, $`N(M)`$, normalized to match a specific sample (but see Yang, Mo & van den Bosch YanMoBos (2002) for improvements in this regard). Even for galaxy mass scales, current models of the halo multiplicity function lead to too few galaxies to match the luminosity or velocity functions of galaxies over broad ranges (see Kochanek Koc (2001); Scoccimarro et al. SSHJ (2001); Chiu, Gnedin & Ostriker Chiu01 (2001)). For general theoretical use, we need the luminosity-dependent multiplicity function, $`N(>L|M)`$, for the average number of galaxies more luminous than $`L`$ in a halo of mass $`M`$. For clusters, we might expect a halo multiplicity function to factor as $`N_{}(M)\mathrm{\Phi }_c(>L)/\mathrm{\Phi }_c(>L_{})`$ where $`\mathrm{\Phi }_c(>L)`$ is the integrated cluster luminosity function and $`N_{}(M)`$ is the mass-dependent number of galaxies more luminous than $`L_{}`$. In the language of Yang et al. (YanMoBos (2002)) this assumption corresponds to assuming a weak dependence of their $`\stackrel{~}{\mathrm{\Phi }}_{}`$, $`\stackrel{~}{L}_{}`$ and $`\stackrel{~}{\alpha }`$ on $`M`$ in this range.

With our modified multiplicity function, Eq. (4), we can match the number density of galaxies in the global Kochanek et al. (Kochanek01 (2001)) K-band luminosity function or the K-band galaxy number counts if we view these galaxies as a complete sample truncated at $`L_{\mathrm{min}}=0.01L_{}`$. The average number of galaxies assigned to a cluster of mass $`M`$ at redshift $`z`$ in a catalog with magnitude limit $`K_{\mathrm{lim}}`$ is $`N(M,z)=N_{FoF}(M)C(z)`$. The completeness function is

$$C(z)=\frac{\mathrm{\Gamma }[1+\alpha ,L(z)/L_{}]}{\mathrm{\Gamma }[1+\alpha ,L_{\mathrm{min}}/L_{}]}$$
(5)
where $`L(z)=\text{max}(L_{\mathrm{min}},L_{\mathrm{lim}}(z))`$ for the limiting luminosity of $`2.5\mathrm{log}(L_{\mathrm{lim}}(z)/L_{})=M_K(K_{\mathrm{lim}}𝒟(z))`$. For our parameters, $`L_{\mathrm{lim}}(z)=L_{\mathrm{min}}`$ at redshift $`cz=1350`$ km/s and $`L_{\mathrm{lim}}(z)=L_{}`$ at $`cz=14000`$ km/s. The average number of galaxies with $`LL_{}`$ is

$$N_{FoF}(M)=\frac{N_{FoF}(M)\mathrm{\Gamma }[1+\alpha ,1]}{\mathrm{\Gamma }[1+\alpha ,L_{\mathrm{min}}/L_{}]}=0.042N_{FoF}(M)$$
(6)
for our standard parameters. On average, clusters with $`Mh/M_{}=10^{15}`$, $`10^{14}`$ and $`10^{13}`$ are assigned $`500`$, $`53`$ and $`6`$ galaxies of which $`21`$, $`2`$ and $`0.2`$ are brighter than $`L_{}`$. Over the range we actually find clusters, this relation is well fit by a power-law of the form $`\mathrm{log}N_{666}=A+B\mathrm{log}(M_{200}h/10^{15}M_{})`$ with $`A=1.32`$ and $`B=0.98`$. We adopt this form to describe the halo occupancy function and summarize all of our subsequent estimates for it in Table 2.

The model catalog contained 100706 galaxies, whose clustering properties are close to those of the observed population. Unlike WK, where we characterized the observed properties of the galaxy simply by the accuracy with which its redshift could be determined, we assigned K-band magnitudes to each galaxy by drawing randomly from the luminosity function. This process oversimplifies several aspects of the real survey. First, we did not include the environment dependence of the luminosity function. Second, for lower mass halos where there is essentially one galaxy per halo, we did not correlate the luminosity with the mass of the halo. Third, to mimic the redshift measurements in the 2MASS sample we assumed that all redshifts were measured for K$`11.25`$ mag and that one-third, randomly selected, were measured between $`11.25<\text{K}12.25`$. This mimics the statistics of the real data, but the redshift measurements in the real data are more likely to be of cluster galaxies than of field galaxies. Finally we ignored the $`|b|>5^{}`$ latitude cut made for the 2MASS data and simply analyzed the synthetic catalogs for the whole sky. Most of these simplifications should make it more difficult to identify clusters in the synthetic catalog than in the real data. Indeed, the visual impression when comparing slices of the model survey and 2MASS data is that the “fingers of god” are weaker in the model. We intend to significantly improve our mock galaxy catalogs in the next iteration, using the relations we derive from the data in later sections (see §6) combined with a better underlying simulation.

## 3. Finding and analyzing the clusters

Large cluster samples have been constructed using five general approaches. In the first method, used to construct the original large catalogs and their successors (e.g. Abell Abe58 (1958); Shectman She85 (1985); Dalton et al. APM (1992); Lumsden et al. EDCC (1996); Ostrander et al. ONRG (1998); Scoddeggio et al. Sco99 (1999); Gladders & Yee GlaYee (2000); White et al. WBL (1999)), clusters were selected as overdensities in the projected density of galaxies on the sky. It was quickly realized such surveys suffer from projection effects and there is a large scatter between optical richness and cluster mass (for recent theoretical studies see e.g. van Haarlem, Frenk & White vHaFreWhi (1997); Reblinsky & Bartelmann RebBar (1999); WK).

Recently, more quantitative versions of these methods based on matched filter algorithms to find either the discrete cluster galaxies (Postman et al. MF (1996); Kepner et al. Kepner99 (1999); Kim et al. Kim01 (2001)) or the excess luminosity from unresolved galaxies (Dalcanton Dal96 (1996); Zaritsky et al. Zaritsky97 (1997); Gonzalez et al. LCDCS (2001)) have been developed to find clusters at intermediate redshifts. A host of approaches have recently been applied to the data from the SDSS (see Bahcall et al. SDSS (2002) for recent discussion and Goto et al. Goto01 (2001) and Nichol et al. C4 (2000) for early work).

Second, with the advent of large redshift surveys, cluster catalogs were constructed by finding three-dimensional overdensities in the galaxy distribution using the friends-of-friends (FoF) algorithm (e.g. Huchra & Geller HG82 (1982); Geller & Huchra GH83 (1983); Ramella et al. Ramella94 (1994); Ramella, Pisani & Geller Ramella97 (1997); Christlein Christlein00 (2000), Ramella et al. Ramella02 (2002)). Studies of the FoF algorithm using N-body simulations (Nolthenius & White Nolthenius87 (1987); Frederic Frederic95 (1995); Diaferio et al. Diaferio99 (1999)) demonstrate the FoF algorithm works well for appropriate choices of the linking parameters. The tree or “dendogram” methods (e.g. Tully Tully87 (1987)) are related to FoF methods.

The third approach to building cluster catalogs is to use X-ray surveys (Gioia et al. Gio90 (1990); Edge et al. Edg90 (1990); Henry & Arnaud HenArn (1991); Rosati et al. Ros95 (1995); Jones et al. Jon98 (1998); Ebeling et al. Ebe98 (1998); Vikhlinin et al. Vik98 (1998); Romer et al. Rom00 (2000); Henry Hen00 (2000); Blanchard et al. BSBD (2000); Scharf et al. ROXS (2000); Böhringer et al. REFLEX (2001); Ebeling, Edge & Henry EbeEdgHen (2001); Gioia et al. Gio01 (2001)). X-ray surveys avoid222But see Lewis et al. (LSEG (2001)). many of the problems in selection and characterization which have made optical catalogs difficult to use, but frequently lack the sensitivity to detect and characterize groups or clusters at $`z0.5`$. Finally, two rapidly developing approaches are to use surveys based on the Sunyaev-Zel’dovich effect (Carlstrom et al. JC (2000)) or on weak gravitational lensing (Schneider & Kneib SchKne (1998); Wu et al. WCFX (1998); Hattori et al. HKM (1999); Wittman et al. Witetal01 (2001); Möller et al. MPKB (2001); but see White, van Waerbeke & Mackey WvWM (2002)).

### 3.1. The matched filter algorithm

We identify groups and clusters in the 2MASS catalog using the modified optimal filter method we developed in WK, which is itself derived from the similar algorithms of Kepner et al. (Kepner99 (1999)) and Postman et al. (MF (1996)). These methods identify clusters as 2- or 2.5-dimensional overdensities in three-dimensional redshift surveys. This is a departure from the essentially uniform use of the FoF algorithm to find clusters in redshift surveys. While the FoF algorithm clearly works to identify groups and clusters with very small numbers of galaxies, a matched filter approach provides several advantages. Matched filters provide likelihood estimates for the detection, membership probabilities for individual galaxies, a range of cluster properties and their uncertainties and (importantly for us) can naturally incorporate galaxies both with and without redshift measurements.

In WK we used N-body simulations to explore the selection of clusters in a series of model surveys for clusters at intermediate redshift. In particular, we developed strategies for selecting samples above a fixed mass threshold over a range of redshifts, showing that it was relatively easy to obtain highly complete samples at the price of significant contamination from real but somewhat less massive clusters. Contamination by detections of non-existent clusters was rare.

The algorithm works as follows. We model the universe as a uniform background and a distribution of $`k=1\mathrm{}n_c`$ clusters. Cluster $`k`$ is described by its angular position $`\stackrel{}{\theta }_k`$, (proper) scale length, $`r_{ck}`$, a fixed concentration, $`c`$, redshift $`z_k`$ and galaxy number $`N_k`$. To ease comparisons to theoretical models we describe the galaxy density profile by an NFW (Navarro, Frenk & White NFW (1996)) profile. In three dimensions, the galaxy distribution normalized by the number of galaxies inside an outer radius $`r_{\mathrm{out}}=cr_c`$ is

$$\rho (r)=\frac{1}{4\pi r_c^3F(c)}\frac{1}{x(1+x)^2}$$
(7)
where $`x=r/r_c`$, $`F(x)=\mathrm{ln}(1+x)x/(1+x)`$ and $`4\pi _0^{cr_c}r^2𝑑r\rho 1`$. The relative number of galaxies inside radius $`r`$ is $`N(<r)=F(r/r_c)/F(c)`$. The projected surface density is more complicated,

$$\mathrm{\Sigma }(R)=\frac{f(R/r_c)}{2\pi r_c^2F(c)}$$
(8)
where

$$f(x)=\frac{1}{x^21}\left[1\frac{2}{|x^21|^{1/2}}\text{tann}^1\left|\frac{x1}{x+1}\right|^{1/2}\right]$$
(9)
and $`\text{tann}^1(x)=\mathrm{tan}^1(x)`$ ($`\text{tanh}^1(x)`$) for $`x>1`$ ($`x<1`$). The relative number of galaxies projected inside cylindrical radius $`R`$ is $`N(<R)=g(R/r_c)/F(c)`$ where

$$g(x)=\mathrm{ln}(x/2)+\frac{2}{|x^21|^{1/2}}\text{tann}^1\left|\frac{x1}{x+1}\right|^{1/2}$$
(10)
(Bartelmann Bart96 (1996)). The projected quantities are considerably more complicated than the simple model ($`\mathrm{\Sigma }(1+x)^2`$) we used in WK, but the earlier model corresponds to a three-dimensional density which is too complicated for convenient use. We fixed the halo concentration to $`c=4`$ and the break radius to $`r_c=200h^1`$ kpc, as typical parameters for cluster-mass halos. We estimate the probability a survey galaxy has its observed properties as either a field or a cluster member galaxy. These properties can include the flux, redshift, angular position and color, with a basic distinction being made between galaxies with redshifts and those without.

### 3.2. Probabilities for Field Galaxies

The probability of finding a field galaxy of magnitude $`K_i`$ and redshift $`z_i`$ is

$$P_f(i)=0.4\mathrm{ln}10D_C^2(z_i)\frac{dD_C}{dz}\varphi _f(K_i𝒟(z_i))$$
(11)
where $`D_C`$ is the comoving distance. The probability of finding a field galaxy of magnitude Ki is

$$P_f(i)=\frac{dN}{dm}(K_i)$$
(12)
where $`dN/dm=_0^{\mathrm{}}𝑑zP_f(K,z)`$ are the differential number counts.

### 3.3. Probabilities for Cluster Galaxies

The probability of finding a cluster galaxy of magnitude $`K_i`$ and redshift $`z_i`$ at an angular separation $`\theta _{ik}`$ from a cluster at redshift $`z_{ck}`$ with (proper) scale length $`r_{ck}`$, galaxy number $`N_{ck}`$, and rest frame velocity dispersion $`\sigma _{ck}`$ is

$`P_c(i,k)=N_{ck}{\displaystyle \frac{\varphi _c(K_i𝒟(z_{ck}))}{\mathrm{\Phi }_c(M_K)}}\mathrm{\Sigma }\left({\displaystyle \frac{D_A(z_{ck})\theta _{ik}}{r_{ck}}}\right)`$
$`{\displaystyle \frac{\mathrm{exp}(c^2(z_iz_{ck})^2/2\sigma _{ck}^2(1+z_{ck})^2)}{\sqrt{2\pi }\sigma _{ck}(1+z_{ck})}}.`$ (13)
The number of galaxies $`N_{ck}`$ is normalized to represent the number of galaxies brighter than $`L_{}`$ inside the radius $`r_{\mathrm{out}}cr_c`$. This normalization is convenient for calculation, and we discuss other normalizations in §3.6. The performance of the algorithm in correctly determining galaxy membership and cluster velocity dispersions is significantly enhanced if we optimize the width of the velocity filter for each cluster candidate. We used the range $`150\text{km s}^1\sigma _{ck}1200\text{km s}^1`$. For a galaxy with an unknown redshift, the likelihood is

$$P_c(i,k)=N_{ck}\frac{\varphi _c(K_i𝒟(z_{ck}))}{\mathrm{\Phi }_c(M_K)}\mathrm{\Sigma }\left(\frac{D_A(z_{ck})\theta _{ik}}{r_{ck}}\right).$$
(14)
We calculate the likelihood in $`1h^1`$ Mpc regions around each galaxy based on either the measured redshift or the average redshift for galaxies with the observed magnitude. The search area was limited to $`\mathrm{\Delta }\theta =4^{}`$. The fraction of cluster galaxies inside search radius $`R=D_A(z_{ck})\mathrm{\Delta }\theta `$ of a cluster at redshift $`z_{ck}`$ is

$$A_k(\mathrm{\Delta }\theta )=\frac{g(R/r_c)}{F(c)}\frac{\mathrm{\Phi }_c(K_{\mathrm{lim}}𝒟(z_k))}{\mathrm{\Phi }_c(M_K)}.$$
(15)
### 3.4. Likelihood Function

The likelihood function for finding the first cluster $`n_c=k=1`$ is

$$\mathrm{\Delta }\mathrm{ln}(k)=N_{ck}A_k+\underset{j=1}{\overset{N_g}{}}\mathrm{ln}\left(\frac{P_f(j)+N_{ck}P_{ck}(j,k)}{P_f(j)}\right)$$
(16)
where the sum extends over all galaxies within angle $`\mathrm{\Delta }\theta `$ of the cluster. We compute the likelihood at the positions of each survey galaxy over a range of trial redshifts, optimizing the number of galaxies $`N_{ck}`$ for each trial. We then add the trial cluster which produced the largest change in the likelihood to the cluster sample. We automate the process and avoid finding the same cluster by including each new cluster in the density model used to find the next cluster,

$$\begin{array}{cc}\mathrm{\Delta }\mathrm{ln}(n_c)\hfill & =N_{cn_c}A_{n_c}\hfill \\ & +_{j=1}^{N_g}\mathrm{ln}\left(\frac{P_f(j)+_{k=1}^{n_c}N_{ck}P_{ck}(j,k)}{P_f(j)+_{k=1}^{n_c1}N_{ck}P_{ck}(j,k)}\right).\hfill \end{array}$$
(17)
The clusters $`k=1\mathrm{}n_c1`$ are our current catalog and have fixed properties, while $`k=n_c`$ is our new trial cluster. By this method we “clean” the clusters out of the sample in order of their likelihood (for further discussion see WK).

We found that several small modifications improved the performance of the algorithm. In most cases these are priors on the filter variables $`N_c`$ and $`\sigma _{ck}`$ designed to stabilize parameter estimation for systems with very few galaxies. First, we added a prior to the likelihood, Eq. (17), of the form $`\mathrm{ln}(N_0^2+N_c^2)`$ with $`N_0=0.1`$ comparable to the value of $`N_c`$ expected for a very poor group. This encodes the information that massive clusters are rare, with a slope roughly matching the power-law slope of the halo mass function. In some senses the prior acts to control the problem of Malmquist biases created by combining a steep number function $`dn/dN_c`$ with the uncertainties in estimates of $`N_c`$. Second, when optimizing the likelihood over the width of the velocity filter $`\sigma _c`$, we included a prior on the likelihood of the form $`\mathrm{ln}(\sigma _c/\sigma _0)`$ with the arbitrary normalization of $`\sigma _0=10^3`$ km/s. This prior makes the likelihood independent of $`\sigma _c`$ when a candidate cluster contains only one galaxy with a redshift measurement333For a maximum likelihood estimate of the velocity dispersion from a set of velocity measurements $`v_i`$, adding this prior leads to an unbiased estimator of the velocity dispersion. If $`\mathrm{ln}=N\mathrm{ln}\sigma +_i(v_i\overline{v})^2/(2\sigma ^2)`$ adding a $`\mathrm{ln}\sigma `$ prior and solving for $`\sigma `$ leads to the unbiased estimator $`\sigma ^2=(N1)^1_i(v_i\overline{v})^2`$. Without the prior the estimator is biased by $`N/(N1)`$.. Third, we added a weak prior on the relation between $`N_c`$ and $`\sigma _c`$ based on our initial results. Preliminary, consistent models found that $`\mathrm{log}N_c=1.13+1.90\mathrm{log}(\sigma _c/1000\text{km/s})`$ with a dispersion of 0.39 dex. After determining $`N_c`$ as a function of $`\sigma _c`$ on a coarse grid, the likelihoods were weighted by this log-normal prior. The prior is only important for setting the filter width $`\sigma _c`$ in systems with too few galaxies for a direct estimate of the velocity dispersion. It also helps to suppress the superposition of clusters into a single larger cluster because superpositions have large apparent velocity dispersions for the number of galaxies. Fourth, for the trial galaxy used to compute the likelihood we used the mean surface density inside the break radius $`r_c`$ rather than the central density of the profile in computing the likelihood. This reduces a bias towards finding clusters centered on isolated galaxies given our cusped density profile.

### 3.5. Matching to Existing Catalogs

Once we have identified a cluster candidate we want to compare its properties either to the known properties of our synthetic catalog or to other properties of the real clusters. We cross reference our cluster catalog to either published or the synthetic cluster catalogs both through the member galaxies and by matching the cluster coordinates and redshifts, separately tracking both identifications. Each galaxy is associated with its most probable cluster, the cluster which maximizes its membership likelihood, and has a likelihood ratio $`\delta _i`$ between the probability that it is a member of that most probable cluster and the probability that it is a field galaxy. This likelihood ratio defines the probability that the galaxy is a cluster member, $`p_i=\delta _i/(1+\delta _i)`$.

In matching our new catalog to existing catalogs we require algorithms which are fully automated and robust. They will not be perfect, but we cannot afford to conduct laborious individual case studies in a broad survey. For our synthetic data, each galaxy is labeled by its parent cluster, allowing us to match the output and input catalogs by identifying the most common parent cluster for all the galaxies with membership probability $`p_i>0.5`$ ($`\delta _i>1`$) in each output catalog. This means of identification is robust, with relatively few multiple identifications or false identifications.

For the real data, the inhomogeneity of local cluster identification enormously complicates cross references between our cluster identifications and previous studies. The primary problem is that cross referencing is best done by cross referencing the proposed member galaxies, but few local group/cluster surveys supply such data. We settled on an ecumenical approach using the NASA Extragalactic Database (NED) to identify the galaxies likely to be associated with each known cluster and then matching our cluster catalog to the existing catalogs based on the member galaxies. We searched NED for all “clusters” (object types GClstr, GGroup, GTrpl & GPair). From this list of 31130 clusters we looked for overlapping identifications which were not recognized by NED. This required some iteration, and was a hopeless task at low redshift (near Virgo) where a nearly infinite number of “groups” have been found, many of which appear to be substructures of a larger object. Based on these identifications we assigned cluster identifications to galaxies within a projected radius $`0.5h^1`$ Mpc and a redshift difference less than $`1000`$ km/s. We then attempted to match our cluster detections to these galaxies using both the match to the cataloged cluster position and redshift and the most common cluster names associated with the galaxies we assign to our detection. The final cross-matching system works reasonably well but is by no means perfect.

### 3.6. Mass Estimates

In addition to identifying the parent cluster and the member galaxies we would also like to estimate the mass of the cluster. We have four different ways of making the comparison. First, for the synthetic catalogs we have direct estimates of the masses. Second, for the real data we have measurements of X-ray luminosities and temperatures which may serve as surrogate estimates of the mass. Third, we have the characteristic number of galaxies in the clusters. Fourth, we have estimates of the galaxy velocity dispersion in the clusters.

Because clusters lack sharp edges, we must exercise some care in defining and comparing our mass estimates (White HaloMass (2001),MassFunction (2002)). In fact, we found that failure to properly match mass definitions affects both the normalization and the slope of our estimates of the halo multiplicity function. For this reason we give explicit details of our procedure below.

We use two standard mass estimates for the clusters in our synthetic catalog. The friends-of-friends mass, $`M_{FoF}`$, is the total mass of the N-body particles linked into a group using the FoF algorithm with a linking length of $`b=0.2`$ in units of the mean inter-particle spacing. We assign galaxies to the clusters based on the FoF mass (Eq. 4) because it can be computed quickly while we are building the synthetic galaxy catalogs. We also calculate separately the mass $`M_{200}`$ inside the radius $`r_{M200}`$ for which the enclosed density is $`\mathrm{\Delta }_M=200`$ times the critical density. When we compare these two estimates we find that $`M_{200}0.66M_{FoF}`$ because the FoF mass corresponds to a lower density threshold of $`\mathrm{\Delta }_M\mathrm{\Omega }_M/b^340`$.

For the real data, where we do not measure the masses directly, we will compare our estimates from the galaxies to X-ray luminosities and temperatures. We used compilations of X-ray luminosities from Böhringer et al. (Bohringer00 (2000)), Cruddace et al. (Cruddace01 (2001)), Ebeling et al. (Ebe96 (1996), Ebe98 (1998), Ebeling00 (2000)), De Grandi et al. (deGrandi99 (2000)), Mahdavi et al. (Mahdavi00 (2000)), and Reiprich & Böhringer (Reiprich02 (2002)), and compilations of X-ray temperatures from David et al. (David93 (1993)), Ponman et al. (Ponman96 (1996)), Markevitch (Mar98 (1998)), Helsdon & Ponman (Helsdon00 (2000)), Ikebe et al. (Ikebe02 (2002)) and White (dwhite00 (2000)). Most of the X-ray luminosity data is drawn from the ROSAT All-Sky Survey (RASS, Voges et al. Voges99 (1999)). We rescaled the X-ray luminosities and fluxes to a common scale by allowing for logarithmic offsets relative to Reiprich & Böhringer (Reiprich02 (2002)) for each survey, with the offsets determined from clusters with measurements in multiple surveys. The offsets we applied are listed in Table 1. Measured scaling relations for clusters of galaxies are (Markevitch Mar98 (1998))

$$\mathrm{log}\left(\frac{L_Xh^2}{10^{44}\text{ergs/s}}\right)=(0.15\pm 0.04)+(2.10\pm 0.24)\mathrm{log}\left(\frac{T}{6\text{keV}}\right)$$
(18)
and (Horner, Mushotsky & Scharf HorMusSch (1999); Nevalainen, Markevitch & Forman NevMarFor (2000); Finoguenov, Reiprich & Böhringer FinReiBoh (2001); Xu, Jin & Wu XJW01 (2001))

$$M_\mathrm{\Delta }10^{15}h^1M_{}\left(\frac{T_X}{1.3\mathrm{keV}}\right)^{3/2}\left(\mathrm{\Delta }_ME^2\right)^{1/2}$$
(19)
where $`\mathrm{\Delta }_M`$ is the density threshold used to estimate the mass $`M_\mathrm{\Delta }`$, and $`E(z)H(z)/H_0=1`$ for our low redshift sample. The normalization of this last relation is currently quite uncertain, at the 30-50% level, due to the notorious difficulties inherent in measuring the mass of a cluster (e.g. Evrard, Metzler & Navarro EMN96 (1996)) and differences in methods for fitting temperatures to the X-ray emitting ICM. (See Fig. 2 of Huterer & White HutWhi (2002) for a compilation of theoretical and observational results.) Different authors use different methods to model the ‘temperature’ of the plasma, different methods of estimating the mass and indeed even different definitions of the mass! To cloud the picture even further, numerical simulations based on purely adiabatic gas physics find a different normalization of this relation than most of the observations (although it must be noted that significant progress in this regard has occurred recently, e.g. Muanwong et al. MTKP (2002), Voit et al. Voit02 (2002)).

The number of member galaxies or their aggregate luminosity must be carefully defined before we can compare it to other results or theoretical models. Our fit parameter $`N_c`$ is defined to be the number of $`LL_{}`$ galaxies inside the fixed, spherical outer radius $`r_{\mathrm{out}}=cr_c=0.8h^1`$ Mpc. Postman et al. (MF (1996), Postman02 (2002)) and Donahue et al. (Donahue01 (2001)) use the luminosity, measured in units of $`L_{}`$, rather than the number. We prefer the number of galaxies, $`N_{}`$, to this equivalent luminosity, $`\mathrm{\Lambda }_{}`$, because it has a closer relation to the Poisson statistics which determine the detectability of clusters and may depend less on the defining luminosity function or wavelength. For the same normalizing aperture, the two quantities are related by

$$\mathrm{\Lambda }_{}L_{}=L_{}N_{}\mathrm{\Gamma }[2+\alpha ]/\mathrm{\Gamma }[1+\alpha ,1]=5.0N_{}L_{}$$
(20)
for $`\alpha =1.09`$. We cannot precisely match the Postman et al. (MF (1996)) definition of $`\mathrm{\Lambda }_{cl}`$ because we use different filters and luminosity functions. To the extent that $`\mathrm{\Lambda }_{cl}`$ is defined by a $`r_{cl}=1.1h^1`$ Mpc ($`1.5h_{75}^1`$ Mpc) radius aperture, we expect $`N_{cl}(g(r_{cl}/r_c)/F(c))N_c=1.6N_c`$.

While $`N_c`$ and $`\mathrm{\Lambda }_{cl}`$ are natural variables for finding clusters, they are poor choices for comparing to theoretical models. For comparisons to theoretical models we will use the number of $`L>L_{}`$ galaxies $`N_\mathrm{\Delta }`$ inside the spherical radius $`r_{N\mathrm{\Delta }}`$ such that the enclosed galaxy density is $`\mathrm{\Delta }_N`$ times the average density of galaxies. Given a cluster with $`N_c`$ galaxies brighter than $`L_{}`$ inside $`r_{\mathrm{out}}`$ we determine $`r_\mathrm{\Delta }`$ and $`N_\mathrm{\Delta }`$ by solving the defining relations

$$N_\mathrm{\Delta }N_cF(r_{N\mathrm{\Delta }}/r_c)/F(c)=\frac{4\pi }{3}n_{}\mathrm{\Delta }r_{N\mathrm{\Delta }}^3\mathrm{\Gamma }[1+\alpha ,1].$$
(21)
For reasons of computational efficiency and stability we keep $`c`$ fixed when we determine $`r_{N\mathrm{\Delta }}`$. In most theoretical models, the distribution of galaxies basically traces the distribution of mass (Katz, Hernquist & Weinberg KatHerWei (1999); Gardner et al. GKHW (1999); Pearce et al. Pearce (1999); White, Hernquist & Springel WhiHerSpr (2001)). However, if we assume the galaxies are unbiased444The ‘galaxies’ in the simulated catalog are very nearly unbiased. The bias of the 2MASS galaxies is currently unknown, but we shall argue later it is likely they are relatively unbiased., to isolate the same physical region we must set $`\mathrm{\Delta }_N=\mathrm{\Delta }_M/\mathrm{\Omega }_M`$ because the galaxy overdensity is defined relative to the average galaxy density while the mass overdensity is defined relative to the critical density rather than the mean mass density. Thus, taking $`\mathrm{\Omega }_M=0.3`$, the number of galaxies brighter than $`L_{}`$ inside the region with the standard mass overdensity of $`\mathrm{\Delta }_M=200`$ is $`N_{666}=0.66N_{FoF}`$ with $`r_{M200}=r_{N666}`$.

Our best theoretical estimate of the conversion to Postman et al. (MF (1996), Postman02 (2002)) and Donahue et al. (Donahue01 (2001)) is

$$\mathrm{\Lambda }_{cl}11N_{666}^{3/4}$$
(22)
for $`\mathrm{\Lambda }_{cl}>10`$. The slope difference is created by the difference between a fixed $`r_{cl}=1.5h_{75}^1`$ Mpc aperture for $`\mathrm{\Lambda }_{cl}`$ and a variable aperture $`r_{N666}`$ for $`N_{666}`$. In practice, both the normalization and slope could differ from this estimate because of differences in the luminosity functions, wavelengths (K band versus I band) and matched filter structures. A simple check is to compare $`\mathrm{\Lambda }_{cl}`$ and $`N_{666}`$ values for clusters of different Abell (Abe58 (1958)) cluster richness classes. Our overlap with the Abell (Abe58 (1958)) catalog is limited because the average redshift of the Abell clusters is significantly higher than our 2MASS catalog. A quick match based on NED cluster names found 31, 39, and 15 Abell richness class 0, 1 and 2 clusters with fractional errors in $`N_{666}`$ smaller than 25% in the catalogs of §5. We find an average richness and dispersion of $`N_{666}=5.3\pm 2.7`$, $`8.9\pm 4.3`$ and $`13.4\pm 6.4`$ for Abell richness classes of 0, 1 and 2 respectively. The ratios of the average $`N_{666}`$ values roughly track the ratios of the number of galaxies associated with each richness class. Postman et al. (Postman02 (2002)) estimate that $`\mathrm{\Lambda }_{cl}=50`$, $`80`$ and $`130`$ should correspond to the typical Abell richness classes 0, 1 and 2 clusters which differs significantly from our estimated conversion to $`\mathrm{\Lambda }_{cl}`$ (Eq. 22) of $`\mathrm{\Lambda }_{cl}=41\pm 16`$, $`60\pm 22`$, and $`84\pm 28`$ for the same richness classes. While using Abell richness classes to compare modern algorithms can hardly be recommended, the comparison suggests that we should explore empirical conversions between $`\mathrm{\Lambda }_{cl}`$ and $`N_{666}`$.

We find clusters with richnesses as low as $`N_{666}0.1`$ in our catalogs, and this requires some explanation. First, low values of $`N_{666}`$ only imply that it is unlikely the cluster contains any bright ($`L>L_{}`$) galaxies. It does not imply the cluster is unlikely to contain galaxies. For example, an $`N_{666}=0.3`$ cluster should contain 3 (7) galaxies with $`L>L_{}/10`$ ($`L>L_{}/100`$) inside its virial radius. We can routinely detect such systems provided the redshift is low enough for us to detect the lower luminosity galaxies. Second, we detect the systems in projection, and the number of galaxies projected inside the virial radius is larger than the number inside the spherical virial radius. This is a modest effect for the NFW profile, but real clusters are embedded in more extended infall regions which boost projection effects. Third, there are significant errors in these low estimates of $`N_{666}`$ because there are so few galaxies. Thus, these low richness clusters can be real virialized systems, but they must be treated with care.

Finally, we have estimates of the velocity dispersion of the clusters. Although we optimized the cluster velocity width $`\sigma _{ck}`$ as part of our matched filter, we do so at low resolution (150 km/s) and it includes contributions from galaxies at large distances from the cluster center. To better mimic a virial velocity for a cluster we use a membership probability weighted estimate of the velocity dispersion. We take all galaxies with measured redshifts, membership probabilities $`p_i0.5`$ and projected separations less than the model outer radius $`Rcr_c=0.8h^1`$ Mpc and estimate a mean redshift

$$cz=cz_ip_i/p_i$$
(23)
and a velocity dispersion

$$\sigma _0^2=\frac{(p_i)^2}{(p_i)^2p_i^2}\frac{(cz_icz)^2p_i}{p_i}$$
(24)
which we then correct to the rest frame $`\sigma =\sigma _0/(1+z)`$. The leading term plays the role of the familiar $`N/(N1)`$ factor needed to make the variance of the velocities about the mean an unbiased estimator of the velocity dispersion. The weighting by membership probability plays an equivalent role to the clipping of outliers in standard approaches to estimating cluster velocity dispersions (see discussion in Borgani et al. BGCYE (1999)). and the method is similar to that used by Barmby & Huchra (Barmby98 (1998)).

### 3.7. The Cluster Number Function

We can also estimate the cluster number function, $`dn/dN`$, which is the distribution of clusters with respect to the number of member galaxies. If the number of member galaxies, $`N(M)`$, is a simple function of the cluster mass then the number function is simply related to the cluster mass function $`dn/dM`$ by

$$\frac{dn}{dN}=\frac{dn}{dM}\left|\frac{dM}{dN}\right|.$$
(25)
If we assume that the cluster mass function is known, then we can infer the halo occupancy function by determining the function $`N(M)`$ needed to transform the assumed mass function into the observed number function. We use the Jenkins et al. (JFWCCEY (2001)) CDM mass function, which should closely match that of our N-body simulation since it is based on numerical simulations with similar cosmological parameters (see White MassFunction (2002) for further discussion). The Jenkins et al. (JFWCCEY (2001)) mass function was derived using $`M_{FoF}`$ masses for the clusters, so we transformed it from $`dn/dM_{FoF}`$ to $`dn/dM_{200}`$ assuming the average relation that $`M_{200}=0.66M_{FoF}`$ we observe in our simulations for the same linking length. This method is similar to the approach used by Berlind & Weinberg (Berlind02 (2002)) and Berlind et al. (BerSDSS (2002)).

We can estimate the number function using the $`V/V_{\mathrm{max}}`$ method (e.g. Moore et al. MooFreWhi (1993)) as follows. At redshift $`z`$, a cluster normalized by $`N_{666}`$ contains

$$N_g(z)=N_{666}\mathrm{\Gamma }[1+\alpha ,L_{\mathrm{lim}}(z)/L_{}]/\mathrm{\Gamma }[1+\alpha ,1]$$
(26)
galaxies, where $`L_{\mathrm{lim}}(z)`$ is the luminosity corresponding to the survey magnitude limit at redshift $`z`$. Suppose we detect all clusters with at least $`N_gN_{\mathrm{thresh}}`$ galaxies. Then we detect clusters in the volume $`V_{\mathrm{max}}(N_{666})`$ defined by the comoving survey volume out to the redshift $`z_{\mathrm{lim}}`$ which solves $`N_{\mathrm{thresh}}=N_g(z_{\mathrm{lim}})`$. The integral number function is then

$$n(>N)=\underset{𝒩N}{}V_{\mathrm{max}}(𝒩)^1$$
(27)
where the sum includes only clusters with more than the threshold galaxy number at their observed redshift. For large values of $`N_{\mathrm{thresh}}`$ we are less sensitive to the Poisson errors in estimating $`N_{666}`$ and minimize the risk of including false positives as clusters. This comes at the price of increased Poisson errors in the number function, because we include fewer clusters, and increasing sample variance because we include less volume. We limited the sample used and $`V_{\mathrm{max}}`$ to $`cz<15000`$ km/s and used bootstrap resampling of the catalog to estimate the uncertainties in the number function.

## 4. Tests With the Synthetic Catalog

We start by testing the algorithm and our recovery of the halo occupancy function in the synthetic data designed to mimic the real 2MASS galaxy sample.

### 4.1. Completeness and contamination

We first consider the completeness of the resulting catalogs. As we discuss in WK, the likelihood associated with a cluster roughly scales with the number of member galaxies, so the distance at which we can detect clusters can be roughly estimated from the number of galaxies. In our models, the average number of cluster galaxies drops to $`N_{FoF}=10`$ ($`3`$) at redshifts of $`cz=12000`$ ($`17000`$), $`7400`$ ($`13000`$), and $`2700`$ ($`7500`$) km/s for clusters of mass $`Mh/M_{}=3\times 10^{14}`$, $`10^{14}`$ and $`3\times 10^{13}`$ respectively. The spread in the likelihoods for a fixed mass and redshift will be considerable because of the Poisson variations in the number of galaxies bright enough to be visible and because of the random variations in the numbers of measured redshifts for the fainter galaxies. In WK we showed that the cluster likelihood essentially scales with the number of member galaxies, so that by scaling the selection likelihood as $`\mathrm{\Delta }\mathrm{ln}(z)=C(z)\mathrm{\Delta }\mathrm{ln}(0)`$ we can extract clusters with a nearly constant mass threshold as a function of redshift. The choice of the zero-redshift likelihood threshold, $`\mathrm{\Delta }\mathrm{ln}(0)`$, is related to the mass cut and its completeness.

Fig. 1 illustrates this for model clusters with $`M10^{14}h^1M_{}`$. Of the 101 clusters with $`M10^{14}h^1M_{}`$ and $`cz<10^4`$ km/s, we miss only 5 at our likelihood cutoff. Of these, one was merged into another cluster, 4 were systems with $`M10^{14}h^1M_{}`$ and fewer galaxies than expected given their mass. We can achieve very high completeness out to the redshift where the typical cluster contains only 3 galaxies. However, all samples with high completeness also have high false positive fractions. If the likelihood threshold is kept high enough, the vast majority of these false positives are real clusters with masses just below the cutoff mass scale – very few of these false positives lack identification with any $`M>3\times 10^{13}h^1M_{}`$ halo. For lower likelihood thresholds the mass range of the false positives broadens and then slowly becomes dominated by unmatched systems where we could not identify the detection with a cluster in the input catalog.

Fig. 2 shows the redshift dependence of the completeness and false positive rates after choosing $`\mathrm{\Delta }\mathrm{ln}(0)`$ to have a 50% false positive rate. The completeness declines slowly with redshift, but remains close to 90% complete from $`0cz10000`$ km/s. The false positive rate averages 50%, and the false positives consist almost entirely of real clusters with masses between $`3\times 10^{13}h^1M_{}`$ and $`10^{14}h^1M_{}`$. False positives that correspond to no cluster become important only for redshifts beyond which a $`10^{14}h^1M_{}`$ cluster contains fewer than 3 galaxies. Within our $`cz10000`$ km/s selection window there is a total volume of $`4.1V_0`$ where $`V_0=(100/h)^3`$ Mpc3, containing 101 clusters more massive than $`10^{14}h^1M_{}`$. Using these clusters to estimate the mass function would be limited by the correction for the contamination by lower mass clusters rather than the correction for completeness or statistical uncertainties including sample variance.

The true completeness of the recovered cluster sample is higher than we have quoted because of the way we have generated the synthetic catalog. In regions of high density, where we expect the most massive clusters, the FoF algorithm with a linking length of 0.2 is prone to merging sub-clumps into larger objects. Our search algorithm instead detects the individual sub-clumps with lower likelihoods than would be expected for an object as massive as the FoF cluster. The sub-clumps may have too low a likelihood to be included (creating a false impression of incompleteness) or multiple clumps may be included (adding to the false positives). Most peculiar high FoF mass but low detection likelihood systems could be traced to this feature of the FoF algorithm.

### 4.2. Halo occupancy function

Given a well-defined cluster catalog, the next step is to examine how well the cluster parameters agree with the true properties of the clusters. We first test our ability to recover the true halo occupancy function $`N(M)`$ given our knowledge of the halo masses for the synthetic catalogs. For the real 2MASS data we will consider cluster X-ray luminosities and temperatures or galaxy velocity dispersions as surrogate, but independent, estimates for the mass.

We test our ability to recover the halo occupancy function by comparing the number of galaxies $`N_{666}`$ inside the $`\mathrm{\Delta }_N=666`$ contrast level. This corresponds to the number of galaxies inside a mass contrast level of $`\mathrm{\Delta }_M=200`$ and should, on average, match $`0.66N_{FoF}`$. There are two critical issues to estimates of the halo occupancy function.

First, definitions are critical! The “number of galaxies in a cluster” must be precisely and similarly defined both for the algorithm used to find the clusters and the theoretical model used to interpret the results. Failure to match the definitions, for example by using $`N_c`$ instead of $`N_{666}`$, makes it impossible test the algorithm (for our synthetic catalog) or to interpret the results (for the 2MASS catalog). In particular, $`N_c>N_{666}`$ for low mass clusters where $`r_{M200}<cr_c=0.8h^1`$ Mpc and $`N_c<N_{666}`$ for high mass clusters where $`r_{M200}>cr_c=0.8h^1`$ Mpc. Thus, the slope of $`N_c(M)M^{0.75}`$ is significantly flatter than the slope of $`N_{666}M`$. We expect the same flattening of the slope will be found for $`\mathrm{\Lambda }_{cl}`$ as used by Postman et al. (MF (1996), Postman02 (2002)) and Donahue et al. (Donahue01 (2001)) for similar reasons.

Second, in estimating relations like $`N(M)`$ it is important to select subsamples which minimize the introduction of mass and number-dependent biases. For example, estimating $`N(M)`$ using a sub-sample selected such that the standard error in $`N`$ is less than a threshold or that the number of galaxies in the cluster exceeds some threshold will preferentially include low mass clusters with high values of $`N`$ over those with low values of $`N`$. This biases estimates of $`N(M)`$ to have flatter slopes than the true relation.

With the logarithmic prior for $`N_c`$ included in the likelihood, our ability to estimate the halo occupancy becomes very robust. Simple redshift cuts work well, but we will define our standard sample by systems which contained at least $`N_v5`$ galaxies with redshifts. After selecting a cluster sample we fit the distribution as a power law including the formal uncertainties in the estimate of $`N_{666}`$ for each cluster. The results, illustrated in Fig. 3, are remarkably good. Above $`10^{13}h^1M_{}`$ the input relation is well approximated by a power law (see Table 2), $`\mathrm{log}N_{666}=A+B\mathrm{log}(M_{200}h/10^{15}M_{})`$ with $`A=1.32`$ and $`B=0.98`$. With our standard sample we find $`A=1.25\pm 0.03`$ and $`B=0.89\pm 0.02`$. The error bars are small because of the large number of clusters (484), so the fits formally disagree with the input relation by 3–4 standard deviations.

The uncertainties in the halo occupancy function are dominated by systematic errors arising from sample selection rather than statistical errors. For example, if we simply include all systems with $`cz15000`$ km/s, roughly the limit for detecting $`L_{}`$ galaxies, we find parameters ($`A=1.20\pm 0.02`$ and $`B=0.92\pm 0.02`$) which has a lower normalization, a steeper slope and more scatter about the mean relation. For systems with $`cz10000`$ km/s, we find $`A=1.23\pm 0.02`$ and $`B=0.95\pm 0.01`$. If we use a likelihood-dependent redshift cut such as $`cz<10000\mathrm{log}(\mathrm{ln})`$ km/s, we find $`A=1.26\pm 0.02`$ and $`B=0.95\pm 0.01`$. We will adopt the $`N_v5`$ sample as our standard from this point, and summarize the estimates of the occupancy function in Table 2.

In summary, we can recover the true halo occupancy function of the synthetic catalog reasonably well. The uncertainties are dominated by systematic errors at the level of 0.05 in the slope and 10% in the normalization of the relation. Some of these problems are due to the generation of the synthetic catalog rather than the search for clusters in the synthetic catalog. The FoF algorithm used to find clusters in the N-body simulation is prone to linking objects which both the eye and our matched filter will regard as neighboring but separate objects. These then appear in our test as objects with high ‘true’ mass but low ‘recovered’ mass, possibly creating a modest bias towards flatter slopes and lower normalizations. While this aspect could be improved upon, we will leave exploration of ‘unlinking’ algorithms to future work.

Unfortunately in the real data we do not know the intrinsic masses of the halos and must use a surrogate estimate such as the cluster velocity dispersion, X-ray luminosity or X-ray temperature. The surrogate we estimate as part of our algorithm is the cluster velocity dispersion $`\sigma `$, and in Fig. 4 we compare our estimated velocity dispersions $`\sigma `$ based on Eq. (24) to the velocity dispersion of the dark matter particles $`\sigma _{200,DM}`$ inside $`r_{M200}`$ of the cluster center. We select the sample for comparison with both our standard limits on the redshift ($`cz<15000`$ km/s) and a minimum number $`N_v5`$ of velocities to be included in the redshift estimate. Fit as a power law with $`(\sigma /10^3\text{km/s})=a(\sigma _{200,DM}/10^3\text{km/s})^b`$ we find little offset, $`a=0.93\pm 0.03`$, or slope difference, $`b=1.04\pm 0.04`$. Restricting the analysis to systems with $`N_v10`$ produces similar results, with $`a=0.89\pm 0.05`$ and $`b=0.94\pm 0.06`$. The scatter is consistent with the uncertainties in the velocity dispersion estimates.

We can also estimate $`N_{}(\sigma )`$, as shown in Fig. 5. In the simulations the velocity dispersion is related to the mass by

$$\sigma _{200}925\left(\frac{M_{200}h}{10^{15}M_{}}\right)^{0.34}\text{km/s}$$
(28)
with a 10% dispersion about the power-law. The power law describing the input $`N_{666}(\sigma _{200})`$ relation has $`N_0=25.6`$ and $`x=2.75`$ for $`N_{666}=N_0(\sigma /10^3\text{km/s})^x`$. The observed relation, with our standard selection procedures, is nearly identical, with $`x=2.76\pm 0.08`$ and $`N_0=21.1\pm 2.0`$, to the input relation. The uncertainties are large because both quantities have relatively large statistical uncertainties. These parameters change little as we adjust the sample selection criteria. Alternatively, we can use Eq. (28) to transform from $`\sigma `$ to $`M_{200}`$ to estimate the occupation number directly, with results very close to the input relation (see Table 2).

### 4.3. The Cluster Number Function

Finally we examine how well we can recover the cluster number function $`dn/dN_{666}`$ of the synthetic catalog. Fig. 6 shows the distributions estimated using thresholds of $`N_{\mathrm{thresh}}=3`$, $`5`$ and $`10`$ in defining $`V_{\mathrm{max}}`$ (see §3.7). The three estimates are mutually consistent, extending over two decades in $`N_{666}`$. We can also compute the distribution for the $`(200h^1\mathrm{Mpc})^3`$ N-body simulation volume from which our synthetic catalog is drawn. The characteristic volume for fair samples of the universe is roughly $`V_0=(100h^1\mathrm{Mpc})^3`$, which we shall use as a unit of volume. In these units the simulation has a volume of $`8V_0`$, and to generate the synthetic catalogs we periodically replicate the simulation. The number function of low $`N_{666}`$ clusters is derived from a small region of the simulation, while the number function of high $`N_{666}`$ clusters is drawn from several of the periodic repetitions of the data. There is also good agreement with the semi-analytic number function found by combining the the synthetic catalog’s multiplicity function with the Jenkins et al. (JFWCCEY (2001)) CDM mass function following the procedures in §3.7 (Eq. 25). The only significant difference is the absence of high $`N_{666}`$ clusters, but this is a real feature of the N-body simulation used to produce the synthetic catalog. If we combine the mass function and the number function to estimate the occupancy function, we obtain mutually consistent estimates of $`N_{666}(M_{200})`$ for $`N_{\mathrm{thresh}}=3`$, $`5`$ and $`10`$ which have slopes and normalizations slightly above (1–2$`\sigma `$) the true relation (see Table 2).

## 5. 2MASS Clusters

With this understanding of the algorithm, we can now search for clusters in the real 2MASS data. We present no catalog at present pending a reanalysis of the full 2MASS data to a fainter limiting magnitude than is presently available. Instead, we focus on the properties of the clusters. In any fitted relation, we have selected the cluster sample used to perform the fit in exactly the same manner as was used for the synthetic catalog, and the order of the analyses parallels that used to test the algorithm on the synthetic catalog. Our standard sample consists of all galaxies identified from the 2MASS extended source catalog (see Jarrett et al. Jarrett00 (2000)) with $`|b|>5^{}`$ to an extinction corrected magnitude limit of K$`12.25`$ mag.

### 5.1. Tests of Completeness

We start with two tests of the completeness of the algorithm for the real data. The first test compares the “shallow” cluster catalogs generated using a bright subsample of the galaxy catalog, truncated at K$`11.25`$ mag, to the “deep” cluster catalogs generated using the complete sample to K$`12.25`$ mag. The second test is to compare our cluster catalog to local X-ray surveys.

#### 5.1.1 Comparison of shallow and deep surveys

We first discuss the results from comparing the shallow and deep cluster catalogs. With a one magnitude change in the survey depth, individual clusters will have significant changes in their membership, as a cluster which contained only $`L>L_{}`$ ($`L>0.1L_{}`$) galaxies in the bright subsample, would on average contain 3.4 (1.5) times as many galaxies in the complete sample. Using the same termination criterion for the search, the shallow catalog contains only 375 clusters, as compared to 1743 in the deep catalog.

We first consider the reliability of the catalog by trying to identify shallow catalog clusters which have no counterpart in the deeper catalog. By cross-referencing the cluster assignments of the galaxies in the two catalogs we estimated a matching probability for identifying a cluster in the shallow catalog with each cluster in the deep catalog. We included all galaxies with a 50% or higher cluster membership probability. The cluster matching probability can range from 0, if no galaxy assigned to a shallow cluster has been assigned to a deep cluster, to unity if every galaxy assigned to the shallow cluster is assigned to the same deep cluster. Only 4 of the 375 clusters lack a firm counterpart in the deeper catalog (defined by a matching probability $`<50\%`$), and an additional 8 may have ambiguous assignments (matching probability of 50–70%). Most shallow clusters have unambiguous counterparts, with 78%, 85% and 91% of the shallow clusters having match probabilities exceeding 95%, 90% and 85% respectively. The match probabilities are not unity because galaxies with modest membership probabilities can be dropped or assigned to a different cluster based on the deeper catalogs, reducing the match probability, even if the galaxies in the cluster cores remain perfectly matched. Thus, if we define a false positive as a cluster which will lack a clear counter part in a deeper survey, our method appears to have a false positive rate of 1–3% depending on the desired threshold. Be warned, however, that this false positive rate is not the same as the true false positive rates we considered for the synthetic data – a real but unvirialized overdensity of galaxies can be identified as a false positive in the synthetic data, but our deeper survey may still identify the clump as a cluster.

Next we examine the relative completeness of the two catalogs. We divided the deep catalog into bins of $`N_{666}`$ and computed the fraction of the deep catalog clusters found by the shallow survey as a function of redshift. For each richness bin, the shallow survey is at least 80% complete until the redshift at which the least rich clusters in the bin are expected to have fewer than 3 galaxies. For example, for the clusters in the deep catalog with $`1<N_{666}<3`$, corresponding to masses near $`10^{14}h^1M_{}`$, the completeness of the shallow catalog is 92% (61%) for redshifts below that where an $`N_{666}=1`$ ($`3`$) cluster is expected to contain 3 galaxies. The deep catalog presumably has similar properties.

The final comparison we make is to compare the estimates of $`N_{666}`$ for the 363 well-matched clusters from the two catalogs. We find that

$$\mathrm{log}N_{666}^{\mathrm{shallow}}=(0.07\pm 0.01)+(1.00\pm 0.01)\mathrm{log}N_{666}^{\mathrm{deep}}$$
(29)
with a variance of about the relation of $`0.21`$ dex. Formally, the error estimates for $`N_{666}`$ would have to be raised by 14% to make the variance consistent with the error estimates. The shallow catalog estimates of $`N_{666}`$ are higher by 17% on average (also found as the median difference). The offset is probably a real systematic effect – the shallow catalog will preferentially include clusters with excess numbers of $`K<11.25`$ mag galaxies because an upward Poisson fluctuation increases the likelihood of including the cluster in the catalog. The same cluster is unlikely to show the same upward fluctuation in the $`K<12.25`$ mag catalog, producing a net bias. This interpretation predicts that a comparison restricted to clusters in the shallow survey with more galaxies, corresponding to clusters with smaller uncertainties in $`N_{666}`$, should have less of a bias. For example, if we restrict the comparison to the 41 clusters in the shallow catalog with errors in $`\mathrm{log}N_{666}`$ smaller than $`0.10`$, the average offset is reduced to 7% ($`0.03\pm 0.02`$) from 17%.

#### 5.1.2 Comparison with X-ray catalogs

The second test compares our complete catalog to X-ray selected samples. It is difficult to do so in a completely satisfactory manner because the selection methods are so different. We first compare our results to the 63 clusters in the HIFLUGCS X-ray flux limited sample ($`f_X>2\times 10^{11}`$ ergs/cm2 s, Reiprich & Böhringer Reiprich02 (2002)). The HIFLUGCS survey covers most of the sky with $`|b|>20^{}`$ plus small holes near the LMC, the SMC and the Virgo cluster. Fig. 7 shows that we find all HIFLUGCS clusters out to the redshifts where we expect even rich clusters to have too few galaxies to be detected by our algorithm ($`cz15000`$ km/s). Our automated matching system identified counterparts to 97% (35 of 36) of the HIFLUGCS clusters with $`cz<15000`$ km/s. We missed one very low redshift system, the NGC 4636 group at $`cz=1320`$ km/s, but closer examination strongly suggests that this is a matching and coordinate problem rather than a genuine failure. Between 15000 km/s and 20000 km/s we find 67% (8 of 12) of the HIFLUGCS clusters. We miss the lower temperature clusters ($`T_X=4.3`$ keV) and find the higher temperature clusters ($`T_X=6.6`$ keV). At higher redshifts, an increasing fraction of X-ray clusters are genuinely missed. Spot checks of these regions show only 0–2 galaxies at the position and redshift of the X-ray cluster.

We can also use our $`N_{666}`$ estimates to estimate an X-ray flux (based on the correlation between $`N_{666}`$ and the X-ray luminosity in Eq. 36). This leads to 18 (14 below $`cz=15000`$ km/s) clusters whose $`N_{666}`$ values and coordinates suggest the cluster should have an X-ray flux large enough to be included in HIFLUGCS. These flux estimates are not very accurate, but we show their distribution, in Fig. 7. They can be divided into three loose categories. First, there are 7 very low redshift systems with $`cz<2000`$ km/s. Some have very high likelihoods in our catalog (the two most prominent are the Eridanus group and the Ursa Major cluster which are the 26th and 38th clusters found), but no published X-ray studies. These systems have virial sizes significantly larger than the ROSAT field of view. Second, there are 5 at intermediate redshifts $`cz5000`$ km/s. The two most prominent systems here, Abell S0805 and Abell 3574, are the 9th and 17th clusters found. Finally, there are 6 systems above $`cz>10000`$ km/s of which the most prominent are Abell 1913 and Abell S0740. After considerable and time consuming effort, we gave up on tracking down the origins of these objects. Many of the systems are well studied clusters with X-ray emission (e.g. Abell S0805) that do not appear in any of the ROSAT catalogs we have used for our comparisons. Some are systems contaminated by X-ray emission from AGN (e.g. Abell 3574). The very low redshift systems have angular extents comparable (or larger) than the ROSAT aperture. In short, understanding these systems in detail requires a careful analysis of the ROSAT data which is beyond the scope of our present effort. We note, however, that Donahue et al. (Donahue02 (2002)) also find examples of optically detected clusters which appear to have fainter than expected X-ray counterparts.

Fig. 7 also shows all other clusters for which we have found a ROSAT X-ray flux measurement but are not derived from complete, flux-limited samples. As these are largely derived from the RASS survey as well, they tend to have X-ray fluxes $`f_X>2\times 10^{12}`$ ergs/cm2 s. The automatic matches find 84% (178 of 213) of the clusters with X-ray flux measurements in our input local cluster catalog and $`cz<15000`$ km/s. As with the HIFLUGCS sample, some of the mismatches at low redshift are due to cross-matching and coordinate problems.

In summary, by comparing our standard cluster catalog to either a shallow catalog or X-ray cluster catalogs, we can show that we achieve very high completeness up to the redshifts where clusters contain too few galaxies for robust detection (about 3 galaxies). It is difficult, however, to use either comparison to robustly estimate a false positive rate. The comparison of the shallow and deep catalogs can only set a lower bound (of a few percent) on the false positive rate. The HIFLUGCS sample sets an upper bound of 40%, but this mainly reflects the limitations of the RASS and our ability to estimate X-ray fluxes from $`N_{666}`$ rather than a true estimate of the false positive rate.

### 5.2. Further Checks of the Velocity Dispersion

Before proceeding to estimates of the halo occupancy function, we made three additional checks of our velocity dispersion estimates, comparing directly with other published values and through comparisons to X-ray observations of the clusters. All fits of variable correlations are performed as $`\chi ^2`$ fits including the errors in both variables and including all clusters with at least $`N_v5`$ associated redshifts that we could match to the other data.

First we checked our velocity dispersions against the velocity dispersions $`\sigma _{\mathrm{pub}}`$ from Girardi et al. (Gir98 (1998); Table 2) and Wu et al. (WuXuFa (1999)). We matched 98 clusters with at least $`N_v5`$ redshift measurements to these previous estimates, and we found good agreement within the rather large errors (Fig. 8). If we fit a power-law, our $`\sigma `$ is related to $`\sigma _{\mathrm{pub}}`$ by

$$\frac{\sigma }{1000\text{km/s}}=(1.09\pm 0.05)\left(\frac{\sigma _{\mathrm{pub}}}{1000\text{km/s}}\right)^{1.06\pm 0.08}$$
(30)
with a dispersion of $`0.16`$ dex that is consistent with the uncertainties. If we restrict the fit to the 69 systems with $`N_v10`$ redshift measurements, we find essentially the same parameters. Hence our membership probability weighting method for determining velocity dispersions is consistent with the more standard methods based on rejecting outliers.

We also find 173 (108) clusters with X-ray luminosities (temperatures) and at least $`N_v5`$ galaxies with redshifts from from the sources summarized in §3.6. Figs. 9 and 10 show our estimates of the relations between the X-ray luminosity $`L_X`$ and the X-ray gas temperature $`T_X`$ to the velocity dispersion estimates. The correlation between the velocity dispersion and the X-ray luminosity, including errors in both quantities, is

$$\mathrm{log}L_{44}=(0.20\pm 0.07)+(4.47\pm 0.24)\mathrm{log}\left(\frac{\sigma }{10^3\text{km/s}}\right)$$
(31)
where $`L_{44}=L_Xh^2/10^{44}\text{ergs/s}`$. The variance in the X-ray luminosity at fixed velocity dispersion is a (typically) large 0.76 dex. If we restrict the fit to systems with $`N_v10`$ the normalization of $`0.18\pm 0.08`$ is little changed, but the slope steepens to $`4.69\pm 0.28`$ with a reduced variance of 0.55 dex. The slope and normalization are close to those found by Mulchaey & Zabludoff (Mulchaey98 (1998), zero point $`0.48\pm 1.09`$, slope $`4.29\pm 0.37`$) or Mahdavi & Geller (Mahdavi01 (2001), zero point $`0.40_{2.0}^{+0.9}`$, slope $`4.4_{0.3}^{+0.7}`$). Wu et al. (WuXuFa (1999)) found a range of slopes from $`2.56\pm 0.21`$ to $`5.24\pm 0.29`$ depending on their fitting methodology. The zero-point uncertainties we cite for these published relations are overestimates because we had to correct the relations to a velocity scale $`10^3`$ km/s from one of $`1`$ km/s without knowing the full covariance matrix of the fit. For the 108 clusters with X-ray temperatures, we find

$$\mathrm{log}\left(\frac{\sigma }{10^3\text{km/s}}\right)=(0.01\pm 0.02)+(0.63\pm 0.04)\mathrm{log}\left(\frac{T_X}{6\text{keV}}\right)$$
(32)
with a spread of 0.16 dex, which is close to the relations found by Mulchaey & Zabludoff (Mulchaey98 (1998), zero point $`0.01\pm 0.04`$, slope $`0.51\pm 0.05`$), Wu et al. (WuXuFa (1999), zero point $`0.00\pm 0.03`$, slope $`0.65\pm 0.02`$), and Girardi et al. (Gir98 (1998), zero point $`0.01\pm 0.03`$, slope $`0.62\pm 0.04`$). Fitting to the 79 systems with $`N_v10`$ gives the same parameters ($`0.01\pm 0.02`$ and $`0.62\pm 0.04`$) with a reduced scatter of 0.12 dex. The similarity with these published analyses suggests that our method for determining velocity dispersions is comparable to those used for other studies despite the small numbers of galaxies used in our dispersion estimates compared to these more extensive surveys of individual clusters.

### 5.3. The Halo Multiplicity Function

We now derive the halo multiplicity function based on the three observables, velocity dispersion, X-ray luminosity and X-ray temperature, which we can later relate to the cluster mass. We first derive the observed correlations and then estimate the conversion to the desired correlation with cluster mass.

Fig. 11 shows that the observed relationship between velocity dispersion and galaxy number is very similar to that in the synthetic data (§4, Fig. 5). The best fit relation for the 939 systems with $`N_v5`$, including the errors in both quantities, is

$$\mathrm{log}N_{666}=(1.37\pm 0.03)+(2.63\pm 0.06)\mathrm{log}(\sigma /10^3\text{km/s})$$
(33)
compared to $`1.32\pm 0.04`$ and $`2.76\pm 0.08`$ for the zero point and slope of the fits to the synthetic catalog ($`1.41`$ and $`2.75`$ for the input relation). Only one system is dropped from the fits. The scatter in $`N_{666}`$ at fixed dispersion is slightly larger in the real data ($`0.53`$ dex versus $`0.50`$ dex). If we restrict the fit to the 238 systems with $`N_v10`$, we obtain consistent parameters ($`1.31\pm 0.04`$ and $`2.60\pm 0.10`$) with a reduced scatter of $`0.41`$ dex. While the significance of the fit in Eq. (33) is quite high, the very large scatter makes it hard to see by eye. Therefore as a check we also used the dispersion in the distribution of galaxy-group velocities in ‘stacked’ clusters, binned in 5 logarithmic bins of $`N_{666}`$, to make the trend more clearly visible. We obtained a fit very close to Eq. (33).

Next we compare the halo occupation number to the X-ray properties of the clusters, as shown in Figs. 12 and 13. We must face two problems in correlating $`N_{666}`$ with the X-ray data. First, we find that the scatter in the correlations is significantly larger than is consistent with the uncertainties in either quantity. In part this is due to the presence of significant intrinsic scatter in the X-ray properties of clusters. Second, there is considerable evidence for a break in the X-ray properties between groups and clusters, generally believed to be due to significant cooling and heating from star formation processes in the lower mass systems (David et al. David93 (1993); Ponman et al. Ponman96 (1996); White, Jones & Forman WJF97 (1997); Allen & Fabian AF98 (1998); Markevitch Mar98 (1998); Arnaud & Evrard ArnEvr (1999); Helsdon & Ponman Helsdon00 (2000); Finoguenov, Reiprich & Böhringer FinReiBoh (2001)). Between the intrinsic scatter and the break in the properties, the parameters for power-law fits to X-ray correlations can depend on the fitting methods and the data used.

In order to include the effects of these problems we made two modifications to our fitting procedures. First, we added a logarithmic systematic error (or intrinsic scatter) $`\sigma _{\mathrm{sys}}`$ in quadrature with the measurement errors for the fitted variables. We then adjusted $`\sigma _{\mathrm{sys}}`$ so that the fit has $`\chi ^2=N_{\mathrm{dof}}`$. Second, we fit power-law relations including the X-ray data only for the more massive clusters with $`T_X1`$ keV and $`L_X10^{42}h^2`$ ergs/s. This reduces the biases in the correlations from attempting to fit the lower mass systems where non-adiabatic effects begin to modify the X-ray properties while leaving enough systems to obtain a statistically significant fit. As in all our standard fits we include only the systems with $`N_v5`$ associated redshifts.

If we now find the best power-law relation between $`L_X`$ and $`T_X`$, we find

$$\mathrm{log}L_{44}=(0.21\pm 0.04)+(2.66\pm 0.12)\mathrm{log}\left(\frac{T}{6\text{keV}}\right)$$
(34)
based on 84 clusters with $`\sigma _{\mathrm{sys}}=0.08`$ and a mean scatter of $`0.30`$ dex. The normalization is consistent with our adopted standard from Markevitch (Mar98 (1998)) in Eq. (18), but it has a significantly steeper slope. If we include no limits on the X-ray properties we find a similar zero-point of $`0.28\pm 0.05`$ but a still steeper slope of $`3.04\pm 0.10`$ and a larger scatter of $`0.35`$ dex. This fit matches that of Xue & Wu (Xue00 (2000)) under similar assumptions. There is clearly a break in the slope of the $`L_X`$$`T_X`$ relation near $`T_X=1`$ keV which is incompatible with a single power law providing a good fit.

We first fit the number-temperature relation as a power law using 84 clusters, as shown in Fig. 12, to find that

$$\mathrm{log}N_{666}=(1.38\pm 0.05)+(2.09\pm 0.17)\mathrm{log}(T/6\text{keV}).$$
(35)
with $`\sigma _{\mathrm{sys}}=0.12`$ dex and a mean scatter of $`0.35`$ dex. This relationship changes little if we drop the restrictions on the X-ray properties. The fit for all 102 clusters with X-ray temperatures, luminosities and $`N_v5`$, has a zero-point of $`1.39\pm 0.05`$, a slope of $`1.97\pm 0.12`$, and a negligible increase in $`\sigma _{\mathrm{sys}}`$ and the scatter.

Fig. 13 shows the results for fitting the 153 clusters with $`L_X10^{42}h^2`$ ergs/s and $`N_v5`$. We find

$$\mathrm{log}N_{666}=(1.29\pm 0.05)+(0.75\pm 0.05)\mathrm{log}L_{44}.$$
(36)
for $`\sigma _{\mathrm{sys}}=0.24`$ and with a scatter of $`0.33`$ dex. If we include all 173 clusters with X-ray luminosity measurements and $`N_v5`$, the zero-point becomes $`1.23\pm 0.04`$ and the slope of $`0.64\pm 0.03`$ is significantly flatter. The distribution of clusters in Fig. 13 is clearly inconsistent with a simple power law if we extend the fit to low X-ray luminosities, so we will not consider this case further. As a consistency check, we can combine the $`N_{666}`$$`T_X`$ and $`N_{666}`$$`L_X`$ relations (Eqs. 36 and 35) to infer a zero-point of $`0.12\pm 0.09`$ and a slope of $`2.79\pm 0.29`$ for the $`L_X`$$`T_X`$ relation which are consistent with the direct fit in Eq. (34). Similarly, if we combine the $`N_{666}`$$`\sigma `$ and $`N_{666}`$$`T_X`$ relations (Eqs. 33 and 35), we infer a zero-point of $`0.00\pm 0.02`$ and a slope of $`0.79\pm 0.07`$ for the $`\sigma `$$`T_X`$ relation which are reasonably consistent with the direct fit in Eq. (32).

We can also compare our results to those of Donahue et al. (Donahue01 (2001)), who compared the properties of clusters found in X-ray surveys with those found using the Postman et al. (MF (1996)) matched filter algorithm. In this algorithm, clusters are characterized by their total luminosity in units of $`L_{}`$, $`\mathrm{\Lambda }_{cl}`$, where we estimated a conversion of $`\mathrm{\Lambda }_{cl}=11N_{666}^{0.75}`$ (see Eq. 22) between our cluster normalizations. Given this relation, we can convert their estimated correlation with X-ray luminosity into our units,

$$\mathrm{log}N_{666}=(1.8\pm 1.0)+(0.8\pm 0.2)\mathrm{log}L_{44}.$$
(37)
While the slope and the normalization are consistent with our estimate given the uncertainties in their relation, the normalization is systematically higher. If we estimate the conversion between $`\mathrm{\Lambda }_{cl}`$ and $`N_{666}`$ by matching the X-ray correlations for the two variables, we find that $`\mathrm{\Lambda }_{cl}=21N_{666}^{0.83}`$, where the slope is consistent with our theoretical estimate (Eq. 22) but the normalization is nearly doubled. Since the sense of the change is the same as that needed to better match $`\mathrm{\Lambda }_{cl}`$ and $`N_{666}`$ for clusters of the same Abell richness class (see §3.6), we will adopt an empirical relation for the conversion of

$$\mathrm{\Lambda }_{cl}21N_{666}^{3/4}$$
(38)
for subsequent comparisons. Using the original $`L_X`$ and $`\mathrm{\Lambda }_{cl}`$ data points, we confirmed that this empirical conversion was correct, while the theoretical conversion from Eq. (22) showed an offset whose origin we do not fully understand.

The final step in the analysis is to convert the observed occupancy functions ($`N_{666}(\sigma )`$, $`N_{666}(T_X)`$ and $`N_{666}(L_X)`$) into an estimate of the occupancy function $`N_{666}(M_{200})`$ as a function of the virial mass. The results, summarized in Table 2, are dominated by systematic errors associated with the choice of the mass scale. The average normalization $`A=1.25\pm 0.03`$ and slope $`B=1.02\pm 0.10`$ found when we set the mass scale using the velocity dispersion ($`\sigma (M_{200})`$, Eq. 28) are lower than the average normalization $`A=1.58\pm 0.02`$ and slope $`B=1.21\pm 0.14`$ found when we set the mass scale using the X-ray temperatures ($`M_{200}(T_X)`$, Eq. 19). The results based on the same mass scale are largely independent of the observed occupancy function used to make the estimate. This is particularly true of the normalization $`A`$, where the scatter between the individual estimates is less than the formal errors, and less true of the slope $`B`$ where the scatter between the individual estimates is somewhat larger than the formal errors. A shallower mass-temperature relation with $`TM^{1/2}`$ rather than $`TM^{2/3}`$, perhaps better suited to the mix of high and low temperature clusters we use, would bring the two estimates of the slopes into agreement. In addition to these systematic differences, recall that in the synthetic data our fits to the occupancy function underestimated the normalization by $`(19\pm 7)\%`$ and the slope by $`0.06\pm 0.03`$ (see §4.2 and Table 2).

### 5.4. The Cluster Number Function

Next we estimate the cluster number function, $`dn/dN_{666}`$, for the 2MASS sample as shown in Fig. 14. The number function for the real sample extends to higher values of $`N_{666}`$ than that for the synthetic sample (see Fig. 6) but otherwise looks very similar. The results for the three values of $`N_{\mathrm{thresh}}`$ are again mutually consistent. Note that for large $`N_{666}`$ the number function is systematically above that of the synthetic catalog.

Comparing our results with previous estimates of the same quantity is quite difficult, due to different conventions for estimating ‘richness’. The best we can do is to apply reasonable transformations based on scaling relations between the widely differing definitions, realizing that this method is far from perfect.

Donahue et al. (Donahue01 (2001)) and Postman et al. (Postman02 (2002)) have estimated the density of clusters at intermediate redshift in terms of $`\mathrm{\Lambda }_{cl}`$ (see §3.6 and Eq. 20) over the ranges $`20<\mathrm{\Lambda }_{cl}<100`$ and $`\mathrm{\Lambda }_{cl}>30`$ respectively. They fit the number density as a power law of the form

$$\frac{dn}{d\mathrm{\Lambda }_{cl}}=N_0\left(\frac{\mathrm{\Lambda }_{cl}}{40}\right)^\alpha $$
(39)
finding $`N_0=6_1^{+3}\times 10^6h_{75}^3\mathrm{Mpc}^3`$ with $`\alpha =5.3\pm 0.5`$ and $`N_0=(1.55\pm 0.40)\times 10^6h_{75}^3\mathrm{Mpc}^3`$ with $`\alpha =4.40\pm 0.30`$ respectively. If we use the empirical conversion between $`\mathrm{\Lambda }_{cl}`$ and $`N_{666}`$ derived from comparing the two richness estimates at fixed X-ray luminosity (Eq. 38), then we find good agreement for the number density of $`N_{666}10`$ clusters (see Fig. 6). If, however, we use our theoretical estimate of the conversion (see §3.6, Eq. 22), then the converted Donahue et al. (Donahue01 (2001)) and Postman et al. (Postman02 (2002)) number densities are nearly an order of magnitude higher than our own estimates. Since both our catalogs (§5.1) and Donahue et al. (Donahue01 (2001), Donahue02 (2002)) have similar completeness compared to X-ray surveys, the need to use the empirical conversion in order to obtain similar number densities suggests that there is a problem with the theoretical conversion between the two richness estimates.

Marinoni, Hudson & Giuricin (Marinoni01 (2000)) used the Nearby Optical Galaxy (NOG, Marinoni NOG (2001)) sample to make a similar estimate based on the group catalogs of Giuricin et al. (Giuricin00 (2001)). If $`\mathrm{\Lambda }_G=L_{\mathrm{tot}}/L_{}`$ is the estimate for the total group luminosity $`L_{\mathrm{tot}}`$ in units of the luminosity of an $`L_{}`$ galaxy, they derive a number function of

$$\frac{dn}{d\mathrm{\Lambda }_G}=\frac{n_0}{\mathrm{\Lambda }_{G0}}\left(\frac{\mathrm{\Lambda }_G}{\mathrm{\Lambda }_{G0}}\right)^{\alpha _G}\mathrm{exp}(\mathrm{\Lambda }_G/\mathrm{\Lambda }_{G0})$$
(40)
where $`\alpha _G1.4`$, $`n_0=4.8h_{75}^3\times 10^4`$ Mpc-3, and $`\mathrm{\Lambda }_{G0}9.9`$. For a galaxy luminosity function with their estimated Schechter slope of $`\alpha =1.1`$, the equivalent number of $`L>L_{}`$ galaxies is $`N_{G0}\mathrm{\Lambda }_{G0}/5.1=1.9`$. For groups selected with a standard FoF algorithm, we expect $`N_{666}0.66N_G`$, but the applicability of this conversion to the Marinoni et al. (Marinoni01 (2000)) estimates is just a rough estimate. The number function $`dn/dN_{666}`$ is identical to the Schechter function $`dn/d\mathrm{\Lambda }_G`$ with $`\mathrm{\Lambda }_G`$ replaced by $`N_{666}`$ and $`\mathrm{\Lambda }_{G0}=1.3`$. This estimate, which is also shown in Fig. 14, is a factor of 5 higher than ours for low $`N_{666}`$ and cuts off more sharply. The sharper cutoff is simply a consequence of the $`cz<6000`$ km/s redshift cutoff of the NOG Sample, which eliminates all the nearby rich clusters. The normalization could be brought into agreement with our estimate if we have overestimated the coefficient in the relation $`N_{666}0.66N_G`$ used to convert the FoF membership into our standard overdensity. Finally we have converted the luminosity function of the CfA sample, as derived by Moore, Frenk & White (MooFreWhi (1993)), to $`dn/dN_{666}`$ in the same manner as described above and with the same difficulties. This is plotted as the other dashed curve which agrees with the NOG estimate at small $`N_{666}`$ but has significantly more objects at high $`N_{666}`$. As with the comparison to Donahue et al. (Donahue01 (2001)) and Postman et al. (Postman02 (2002)), differences in the absolute number density are very sensitive to errors in the conversion between richness estimates, and we lack any means of doing this more precisely given the information available in Marinoni et al. (Marinoni01 (2000)) and Moore et al. (MooFreWhi (1993)).

Table 2 also gives the estimates for the halo occupancy function needed to transform the Jenkins et al. (JFWCCEY (2001)) mass function for our assumed cosmology into the observed number function (as in section 4.3). The results for $`N_{\mathrm{thresh}}=3`$, $`5`$ and $`10`$ are mutually consistent and intermediate to the results from fitting the observed occupancy functions in §5.3.

### 5.5. Derived Quantities

Finally we discuss the implied infrared mass-to-light ratio of clusters and the fraction of the local luminosity density associated with the virialized regions of clusters. We start by adopting a standard model for the halo occupation function. Our uncertainties in the occupation function are dominated by systematic errors in how to relate the observed quantities to cluster masses. Depending on whether we set the mass scale using the velocity dispersion, the X-ray temperature, or the cluster mass function, we find significant differences in both the normalization and the slope of the occupancy function (see Table 2). We decided to simply average the results for the three possible mass scales and use their dispersions as the standard errors:

$$\mathrm{log}N_{666}=(1.44\pm 0.17)+(1.10\pm 0.09)\mathrm{log}(M_{200}h/10^{15}M_{}).$$
(41)
The resulting errors are larger than the statistical uncertainties for any given mass normalization method and somewhat larger than the systematic offsets we found when estimating the occupancy function of the synthetic catalog. The normalization is higher than we used in the synthetic catalog, which seems reasonable given the visibly weaker fingers of god in the synthetic catalogs.

From Eq. (41) we can estimate more traditional quantities like the K-band mass-to-light ratios of the clusters. Including the uncertainties in the luminosity function, the correction from aperture to total magnitudes, and a zero point of $`M_{K_s}=3.39`$ mag (see Kochanek et al. Kochanek01 (2001) and Cole et al. Cole01 (2001)) we find that

$$\left(\frac{M_{200}}{L_{666}}\right)_K=(116\pm 46)h\left(\frac{M_{200}h}{10^{15}M_{}}\right)^{0.10\pm 0.09}$$
(42)
where almost all the uncertainty comes from the uncertainty in the normalization $`A=1.44\pm 0.17`$ of the occupancy function (it would be $`(116\pm 10)h`$ for $`A1.44`$). This is nearly identical to the baseline semi-analytic model of Cole et al. (Cole00 (2000)) where the K-band mass-to-light ratio is $`118(M_{200}h/10^{15}M_{})^{0.06}`$ for clusters with masses above $`10^{13.5}h^1M_{}`$. The mass-to-light ratio is essentially independent of mass, which differs somewhat from estimates in the B-band (e.g. Girardi et al. girardi00 (2000), Marinoni & Hudson Marinoni01b (2001)) which suggest a steeper slope of $`(M/L)_BM^{0.10.3}`$, albeit with significant uncertainties in the slopes. If there is a real slope difference, with $`(M/L)_B/(M/L)_KM^x`$ then we can explain this as a trend in the colors of the galaxy population with $`\mathrm{\Delta }(BK)=5x0.5`$ to $`1.5`$ mag over the range from $`10^{13}M_{}`$ to $`10^{15}M_{}`$. All general trends in cluster properties have the sense required to produce the slope difference. Starting from a constant K-band mass-to-light ratio, the more massive clusters could be dimmer in the B-band because of the increasing early-type galaxy fraction (a B–K color shift of roughly 0.5 mag from Sa to E), increasing metallicity ($`\mathrm{\Delta }(`$B–K$`)/[Fe/H]1.2`$ mag/dex) or increasing age ($`\mathrm{\Delta }(`$B–K$`)/\mathrm{log}t0.5`$ to $`1.0`$ mag/dex). Like the halo occupation number itself, it is important to match definitions when comparing our mass-to-light ratios to other estimates. Our estimate is normalized to be the ratios of the two quantities in spheres of radius $`r_{N666}r_{M200}`$, leading to higher mass to light ratios than estimates comparing the light in cylinders to the mass in spheres. For example, Rines et al. (Rines01 (2001)) found $`(M/L)_K=(75\pm 23)h`$ for Coma using the 2MASS survey, and the Kochanek et al. (Kochanek01 (2001)) luminosity function. At the $`r_{N666}=1.5h^1`$ Mpc, which also matches the Rines et al. (Rines01 (2001)) estimate of the virial radius, we must raise their estimate of the mass to light ratio by 25% to $`(M/L)_K=(92\pm 28)h`$ to compare it to our estimates of the mass-to-light in spheres. This is quite close to our estimate of $`(M/L)_K=(116\pm 46)h`$ for the mass-to-light ratio of $`10^{15}h^1M_{}`$ cluster like Coma in Eq. (42).

Since the K-band luminosity is well correlated with the total stellar mass, we can use the cluster number function to obtain a rough estimate of the fraction of stars in systems above a given mass or number threshold. The total luminosity density of all galaxies is $`j_{\mathrm{tot}}=n_{}L_{}\mathrm{\Gamma }[2+\alpha ]`$ while the equivalent luminosity density contained clusters with $`N_{666}>N`$ (and in the region with $`r<r_{N666}`$) is

$$j_{\mathrm{clust}}(>N)=\frac{L_{}\mathrm{\Gamma }[2+\alpha ]}{\mathrm{\Gamma }[1+\alpha ,1]}_N^{\mathrm{}}𝑑N_{666}N_{666}\frac{dn}{dN_{666}}.$$
(43)
The ratio $`j_{\mathrm{clust}}/j_{\mathrm{tot}}`$ is shown in Fig. 15 as a function of $`N_{666}`$. We find that the virialized regions of clusters more massive than $`10^{14}h^1M_{}`$ contain $`7\%`$ of the local luminosity. This is quite comparable to the fraction of the mass contained in such clusters in popular $`\mathrm{\Lambda }`$CDM models with $`\mathrm{\Omega }_{\mathrm{mat}}0.3`$ and $`\sigma _81`$, though the number is particularly sensitive to the latter assumption. The enclosed fraction also depends on the density contrast to which we extend the cluster edges. For example, if we include galaxies out to $`\mathrm{\Delta }_N=100`$, the luminosity fraction in clusters doubles.

Since our data is consistent with $`N(M)M`$ we would predict that 2MASS galaxies provide good tracers of the mass with a small and relatively scale-independent bias on Mpc scales. This would suggest that the correlation function (or power spectrum) of 2MASS galaxies should show the virial ‘inflection’ of the dark matter power spectrum on Mpc scales, where the growth of fluctuations is first enhanced by non-linear infall and then suppressed by virial motions within halos. There is weak evidence for this in Fig. 3 of Allgood, Blumenthal & Primack (ABP (2001)), though further work with the full catalog needs to be done to really test this prediction. That $`N(M)M`$ is perhaps not unexpected, as both simulations and observations show that galaxies selected in the “red” trace the matter distribution better than those selected in the “blue”. The original semi-analytic models (Kauffman et al. KCDW (1999)) predicted a steeper slope for $`N(M)`$ for their red galaxies, with $`N(M)M^{0.9}`$, which is consistent with our results.

Theory predicts, and our simulations have assumed, that the number of galaxies in a cluster mass halo should follow a Poisson distribution (at fixed mass). It is difficult to test this with our current sample due to the difficulty of estimating the mass of our clusters. The ratio of $`N(N1)^{1/2}/N`$, where $`N`$ is the number of member galaxies determined from the matched filter probabilities $`p_i`$, is very close to unity for clusters binned by $`N_{666}`$. (The Poisson prediction is $`N(N1)^{1/2}=N`$.) However $`N_{666}`$ is not a good proxy for $`M_{200}`$, having significant scatter. For the small number (110) of clusters for which we can detect all $`L_{}`$ galaxies and we have an X-ray temperature, the ratio $`N(N1)^{1/2}/N`$ is also consistent with unity, though we need to choose broad bins in $`T_X`$ and the statistics are poor.

## 6. Conclusions

We have applied the matched filter technique to both simulated galaxy catalogs and the 2MASS galaxy catalogue to search for clusters over approximately 90% of the sky to a redshift limit of $`z0.05`$. We have matched our 2MASS derived catalog to existing catalogs in an automated way using the NASA Extragalactic Database. Our algorithm appears to be both robust and efficient, returning quite complete samples of clusters out to distances where the typical cluster contains as few as 3 galaxies. The algorithm finds almost no ‘false groups’, the main source of contamination being the inclusion of objects below the mass cut in the catalog.

The matched filter algorithm gives us a way of estimating cluster membership and hence cluster properties like the multiplicity function, the number function and the velocity dispersion in a new way. While our estimates for the velocity dispersion agree well with earlier work, we typically have far fewer galaxies per cluster than dedicated surveys and hence larger errors. This is offset by the fact that we have a large number of clusters, over much of the sky, with which to search for correlations between cluster properties such as velocity dispersion and X-ray temperature. Where there is overlap we find quite good agreement with earlier work, though often with improved statistics.

Although it is necessary to be very careful about the definition of $`N`$ in order to compare with theory, we find that with sufficient care we have been able to estimate the cluster number function, $`dn/dN`$, over more than 2 orders of magnitude in $`N`$.

By using the velocity dispersion, X-ray luminosity or X-ray temperature as a surrogate for mass we are able to estimate the multiplicity function $`N(M)`$. Although there are serious issues in converting observables into cluster masses, all of our results are fairly consistent and suggest that $`N(M)M`$ or perhaps slightly steeper, in reasonable agreement with earlier estimates. With all quantities normalized by the spherical radius corresponding to a mass overdensity of $`\mathrm{\Delta }_M=200`$ or the equivalent galaxy number overdensity of $`\mathrm{\Delta }_N=200\mathrm{\Omega }_M^1=666`$, we find that the number of $`L>L_{}`$ galaxies in a cluster of mass $`M_{200}`$ is

$$\mathrm{log}N_{666}=(1.44\pm 0.17)+(1.10\pm 0.09)\mathrm{log}(M_{200}h/10^{15}M_{}).$$
(44)
The uncertainties in this relation are largely due to the choice made for relating the observed quantities to the cluster mass scale. For a fixed mass scale the scatter resulting from the different observed correlations is considerably smaller. Correlations of $`N`$ with other cluster properties, X-ray luminosity, temperature or galaxy velocity dispersion, are given in §5.3. The region has a K-band cluster mass-to-light ratio of $`(M/L)_K=(116\pm 46)h`$ which is essentially independent of cluster mass. The uncertainties are again dominated by the choice of the mass scale. This scaling is consistent with $`N(M)M`$, though if we take our best fit seriously and $`N(M)`$ is steeper than $`M`$ we expect $`M/L`$ would fall slowly with increasing mass. Integrating over all clusters more massive than $`M_{200}=10^{14}h^1M_{}`$, the virialized regions of clusters contain 7% of the local stellar luminosity, quite comparable to the (somewhat theory dependent) mass fraction in such objects in currently popular $`\mathrm{\Lambda }`$CDM models.

The cluster likelihoods tend to be larger in the real data than our mock catalogs. This difference could have several sources. It could be due to the concentration of redshift measurements in the real catalog towards groups and clusters rather than having the random distribution of the synthetic catalog. It may also indicate that the cosmology adopted in the underlying simulation is incorrect (the number of clusters is very sensitive to $`\sigma _8`$ for example), the normalization of the input $`N(M)`$ may be too low, our assumed $`N(M)`$ may have too many galaxies in low-mass halos compared to high-mass halos, or the luminosity function could vary systematically with the parent halo mass rather than being fixed (as we have assumed). Along these lines, both the expectation that light closely trace mass in $`K`$-band and that the luminosity function change with mass suggest that clusters contain a larger fraction of galaxies than we have assumed in our current generation of mock catalogs. These avenues will be explored in the future with larger simulations and better data.

We view this work as but the first step in an iterative sequence. Based on simulations which we know describe reality imperfectly, we have calibrated our cluster finding algorithm. When applied to the real data this algorithm allows us to estimate correlations between different properties of a halo which can be used in the next stage to improve the simulated catalogs. As the 2MASS catalog becomes increasingly complete, and more redshifts become available, we will estimate the global relations between different halo parameters by providing them as priors to the fitting, and varying the parameters to optimize the global likelihood. These relations can then be used in the construction of the galaxy catalogs from improved simulations which will allow us to further optimize and understand the cluster finding algorithm.

## Acknowledgments

We would like to thank the 2MASS collaboration for their cooperation and comments on this paper. We would also like to thank N. Caldwell, M. Donahue, D. Eisenstein, C. Jones, M. Pahre, M. Postman and K. Rines for helpful conversations and C. Baugh & S. Cole for providing additional IR properties for clusters in their semi-analytic models. M.W. would like to thank H. Mo and R. Yan for enlightening conversations on the halo model of galaxy clustering. This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. C.S.K. and J.P.H. were supported by the Smithsonian Institution. M.W. was supported by a Sloan Fellowship, the NSF and NASA. Simulations were carried out at CPAC, through grants PHY-9507695, and at the National Energy Research Scientific Computing Center.

## Appendix A Including Color Information

In this appendix we outline the inclusion of color information in our algorithm. To include color information in our algorithm we require a model for the multi-variate luminosity function in luminosity and colors, $`\varphi (M,C_0)`$, and its evolution with redshift. The multi-variate and standard luminosity functions are related by $`\varphi (M)=\varphi (M,C_0)𝑑C_0`$. Assuming that the numbers of galaxies are not evolving, the measured magnitudes and colors, $`m`$ and $`c`$ are related to the local values by $`M=m𝒟(z)`$ (Eq. 1) and $`C_0=c𝒞(M,z)`$ to remove the effects of distance, K-corrections and evolution.

For the 2MASS survey the description simplifies further because only the H–K color has a significant correlation with redshift and neither the H–K nor the J–H color has a significant correlation with luminosity. This means that we can factor the bivariate luminosity function, $`\varphi (M,C_0)=\varphi (M)\xi (C_0)`$, with $`\xi (C_0)𝑑C_01`$. The color distribution, which consists of a narrow peak with a red tail, is well modeled by the sum of two Gaussians,

$$\xi (HK)=\frac{f}{\sqrt{2\pi }\sigma _0}e^{(HKa_0)/2\sigma _0^2}+\frac{1f}{\sqrt{2\pi }\sigma _1}e^{(HKa_1)/2\sigma _1^2}$$
(A1)
where $`f=0.698`$, $`a_0=0.226`$, $`a_1=0.268`$, $`\sigma _0=0.025`$ and $`\sigma _1=0.081`$. The galaxies become steadily redder at higher redshift, with

$$𝒞(z)=2.333z0.770z^2,$$
(A2)
and it appears broader at fainter magnitudes due to measurement errors that can be modeled by convolving the distribution with a Gaussian of width $`\sigma =0.092\times 10^{0.2(K13)}`$. While the errors in the flux are unimportant to our analysis given the sample magnitude limit, the shallow slope of $`𝒞(z)`$ means that we must take into account the errors in the H–K colors.

Given the bivariate luminosity function and the dependence of the colors on redshift, we can modify the terms of the likelihood function (see §3.2 and §3.3) to include the additional information. The probability of finding a field galaxy with magnitude $`m_i`$, color $`c_i`$ and redshift $`z_i`$ is

$$P_f(m_i,c_i,z_i)=0.4\mathrm{ln}10D_C^2(z_i)\frac{dD_c}{dz}\varphi _f(m𝒟(z_i))\xi (c_i𝒞(z_i)).$$
(A3)
Integrating over redshift, we obtain the number counts of galaxies with a given color, which is equivalent to the probability of finding a galaxy with magnitude $`m_i`$ and color $`c_i`$,

$$P_f(m_i,c_i)=_0^{\mathrm{}}P_f(m_i,c_i,z)𝑑z.$$
(A4)
The probability distribution for an unknown galaxy redshift is simply the ratio $`P_f(m_i,c_i,z_i)/P_f(m_i,c_i)`$. Because the slope of $`𝒞`$ is shallow, the addition of the H–K color information only modestly improves the accuracy of redshift estimates. For the 50000 galaxies with known redshifts, the observed errors in the predicted redshifts are modestly smaller than their theoretical estimates of $`\sigma _z=0.028`$ (using only fluxes), $`0.036`$ (using only colors) and $`0.024`$ (using both). The errors in redshifts estimated using colors (fluxes) grow slowly (linearly) with redshift, so at low redshifts fluxes always become better redshift estimators than colors. We must also modify the cluster membership probabilities, replacing $`\varphi _c(M)`$ by $`\varphi _c(M,C_0)`$, while leaving the normalization factor in the denominator, $`\mathrm{\Phi }_c(M)`$, unchanged.

This formalism for including color information in a matched filter algorithm can be generalized to additional colors and more complicated descriptions for the evolution of colors with redshift and luminosity, although tabulated models will probably be required. Matched filters can also exploit the redder optical colors of cluster galaxies (as used in other algorithms, e.g. Goto et al. Goto01 (2001)) by using a multi-variate cluster luminosity function with a different color distribution from that for the field.