# Group Identification in N-Body Simulations: SKID and DENMAX Versus Friends-of-Friends**footnote *Preprint BROWN HET-1140

## 1. Introduction

N-body simulations are a frequently used tool in cosmology to study the origin of large-scale structure in the universe or the formation of galaxies. For these applications, it is necessary to identify the dark matter halos in the models, as those are interpreted to be the locations of galaxies or galaxy clusters. Ideally, the way in which the dark matter halos are extracted should not influence the results sought in the simulations, e.g. the mass distribution, correlation function, formation and merging rates, and many other properties of the halos. Fortunately or unfortunately, many different algorithms to identify groups of particles in N-body simulations have been proposed. The question we address here is how the choice of algorithm affects the properties of halos.

The attention of the authors was drawn to this question by the results shown in Fig. 1. This shows the mass distribution of halos, in the form of the differential mass function, at redshift $`z=0`$ for a series of N-body simulations, run with a particle-particle/particle-mesh (P<sup>3</sup>M) code by E. Bertschinger and J. M. Gelb . The simulations, originally carried out for a different purpose, are of a standard cold dark matter (SCDM) model ($`\mathrm{\Omega }_0=1`$ and $`H_0=50\text{ km/s/Mpc}`$) with $`64^3`$ particles in boxes of sizes varying from 50 Mpc to 1600 Mpc. The mass distribution of halos was determined in two different ways, using two popular algorithms: friends-of-friends (FOF), a percolation algorithm, and DENMAX. (See below for a description of how these algorithms work.) Fig. 1 reveals that the two algorithms give very different results for the number of massive halos in large simulations. In particular, DENMAX overpredicts their number when compared to FOF or the theoretical prediction based upon Press-Schechter theory. (Again, see below for an explanation of Press-Schechter theory.) Thus the choice of group identification algorithm indeed has an effect on the properties of halos, and this article attempts to explore this by comparing three of the more popular algorithms.

Section 2 gives a brief overview of the N-body simulations used to provide the ‘raw material’ for the comparison. Section 3 introduces the three different algorithms which had been chosen for that purpose. Section 4 describes Press-Schechter theory, which is used as a baseline to compare the mass distribution of halos as identified by the different algorithms, and how the agreement (or disagreement) with Press-Schechter theory can be quantified. The next section contains the actual results of the comparisons, how these depend upon the choice of the free parameters in the algorithms, the box sizes of the simulation, and the cosmology. This is then followed by the conclusions.

## 2. The Simulations

To test the reliability of group identification algorithms we use representative results of N-body simulations. The first indications of the effects of different group identification algorithms were found in N-body simulations with the P<sup>3</sup>M code by E. Bertschinger and J. M. Gelb . We decided to use a different code for this systematic study. The simulations were run with the Hydra code by H. M. P. Couchman et al. . Hydra incorporates the evolution of both, dark matter via an adaptive particle-particle/particle-mesh method (AP<sup>3</sup>M), and of baryonic gas with smoothed particle hydrodynamics (SPH). It was run with dark matter particles only since the goal was to provide typical cosmological settings for halo identification. Hydra was chosen because of its higher computational efficiency in high density regions, which are the regions we are particularly interested in, and because software to identify groups in the output files produced by Hydra is readily available. Two “favorite” cosmological scenarios were chosen: an open cold dark matter model (OCDM) with $`\mathrm{\Omega }_0=0.4`$, a Hubble constant of $`H_0=60\text{ km/s/Mpc}`$, normalized in such a way that $`\sigma _8=0.8`$ at present time, consistent with the observed abundance of galaxy clusters . The second model was standard cold dark matter (SCDM) with $`\mathrm{\Omega }_0=1`$, $`H_0=50\text{ km/s/Mpc}`$, and $`\sigma _8=1`$. For each model, simulations were run with $`64^3`$ particles with box sizes varying from 50 Mpc to 1600 Mpc in steps of a factor of two both to insure a wide dynamic range in the masses of the dark matter halos and to see at which resolution the group identification algorithms break down. Table 1 lists all the parameters of the simulations, including the gravitational softening employed, initial redshift, and the number of timesteps.

## 3. Identifying Groups

At redshift zero, groups in the distribution of dark matter particles were identified using three common algorithms: the friends-of-friends (FOF) method, DENMAX, and SKID. With FOF , particles are joined into groups if the separation to the nearest neighbor is less than a given threshold, called the linking length $`b`$, which is the only free parameter for this algorithm. We will express $`b`$ in units of the mean particle separation throughout this paper. Then $`1/b^3`$ corresponds to an overdensity, and FOF approximately groups together particles which lie inside the corresponding level surfaces.

Group identification in DENMAX  and SKID  is a three-step process. First, particles are moved in small steps in the direction of the local density gradient, until their locations stay within a small distance, the convergence radius, for a certain number of steps, which means that they oscillate around a local density maximum. In the second step, groups among the particles at their moved positions are identified with a FOF algorithm with a linking length $`b`$ twice the convergence radius. The third and last step then removes particles which are not gravitationally bound to the groups identified in the previous step. This is done by calculating the potential energies of all of the particles in the group (now at their original, unmoved locations) and their kinetic energies (with peculiar velocities with respect to the center of mass of the group and an additional velocity due to the Hubble expansion of the halo). If there are unbound particles with positive total energy, the one with the highest total energy is removed from the group, and the process of calculating potential and kinetic energies is repeated, until there are no unbound particles left. Note that this implementation of DENMAX actually follows SKID and is slightly different from the originally proposed . The difference between DENMAX and SKID lies in the way densities, and thus their gradients, are determined in the first step. For DENMAX, they were calculated on a fixed grid with a triangular shaped cloud (TSC) scheme \[10, Ch. 5-3-2\], the same as in the P<sup>3</sup>M code. SKID instead uses a symmetric SPH smoothing kernel  considering a certain number of neighbors around each particle to smooth over. Hence the main difference between the two algorithms is that DENMAX computes density gradients on a fixed Eulerian scale, whereas SKID computes them on a fixed Lagrangian scale, adjusting the resolution of this calculation to the local density. Different to FOF, DENMAX and SKID have two free parameters each: the linking length $`b`$ and the size of the grid (for DENMAX) or the number of nearest neighbors (for SKID) used in calculating the density gradients.

## 4. Press-Schechter Theory and How to Compare With It

One of the most basic quantities derived from identifying halos is the distribution of their masses, which should be insensitive to different algorithms. To this end, the differential mass function $`n(M)`$, i.e. the number density of halos with masses between $`M`$ and $`M+\mathrm{d}M`$ as obtained in the simulations and by the different algorithms, was compared with the prediction from Press-Schechter theory  based upon the power spectrum $`P(k)`$ used for the initial conditions in the simulations

$$\begin{array}{c}n_{\mathrm{PS}}(M)\mathrm{d}M=\\ \sqrt{\frac{2}{\pi }}\frac{\overline{\rho }}{M^2}\frac{\delta _\mathrm{c}}{D\sigma (M)}\frac{\mathrm{d}\mathrm{ln}\sigma }{\mathrm{d}\mathrm{ln}M}\mathrm{exp}\left[\frac{\delta _\mathrm{c}^2}{2D^2\sigma ^2(M)}\right]\mathrm{d}M.\end{array}$$
(1)   
We use the form of the differential mass function as given by  for redshift $`z=0`$, the case we are interested in. $`\overline{\rho }`$ is the present mean mass density of the universe, $`\delta _\mathrm{c}`$ the linearly extrapolated overdensity at which a top-hat shaped overdensity has collapsed ($`\delta _\mathrm{c}=1.686`$ for SCDM and $`\delta _\mathrm{c}=1.660`$ for OCDM at $`z=0`$ — see  for a discussion of this quantity and its dependence upon cosmological parameters and redshift). $`\sigma `$ is the fractional r.m.s. mass fluctuation within a top-hat of radius $`R`$ in the initial conditions of the simulation

$$\sigma ^2(R)=4\pi _0^{\mathrm{}}k^2dkP(k)W_R^2(k),$$

where $`W_R(k)`$ is the Fourier transform of the top-hat window function

$$W_R(k)=3\left[\frac{\mathrm{sin}kR}{(kR)^3}\frac{\mathrm{cos}kR}{(kR)^2}\right].$$

The connection between $`\sigma `$ as a function of radius and sigma as a function of mass, which is needed in evaluating eq. (1), is then

$$M=\frac{4\pi }{3}R^3\overline{\rho },$$

via the average mass contained within the top hat. Finally, $`D`$ in eq. (1) is the ratio of the growing modes of linear perturbation theory $`D(z)`$ at present time and at the redshift of the initial conditions $`D=D(z=0)/D(z=z_{\mathrm{initial}})`$.

In the interpretation and derivation of , the Press-Schechter prediction (1) gives the mass distribution of regions which lie inside contours with density contrast $`\delta _\mathrm{c}`$, which is the linearly extrapolated density contrast at which a full calculation, which is possible for a top-hat shaped overdensity, shows that the region has already collapsed.

Since one of the group identification algorithms (FOF) showed a good agreement with Press-Schechter theory in the simulations which drew our attention to this problem (Fig. 1), we decided to use Press-Schechter theory as a baseline for comparing the different algorithms. To quantify the differences in the differential mass functions obtained in the simulations, $`n(M)`$, and from Press-Schechter theory, $`n_{\mathrm{PS}}(M)`$, we simply integrate the difference squared between appropriate masses $`M_0`$ and $`M_1`$ to obtain a goodness of agreement $`\mathrm{\Delta }`$

$$\mathrm{\Delta }=_{\mathrm{lg}M_0}^{\mathrm{lg}M_1}\left[\mathrm{lg}n(M)\mathrm{lg}n_{\mathrm{PS}}(M)\right]^2\mathrm{d}\mathrm{lg}M,$$
(2)   
but with the logarithms taken of all relevant quantities. (Here, we use logarithms to base 10 and express the differential mass functions in units of Mpc$`{}_{}{}^{3}M_{}^{1}`$ and masses in $`M_{}`$.) Since the differential mass functions approximately drop as a power law with increasing mass, introducing the logarithms in (2) assures that differences at the low-mass end contribute as much to $`\mathrm{\Delta }`$ as differences at the high-mass end.

## 5. Results

The three different group finding algorithms, with varying parameters, were applied to the results of the N-body simulations, and the resulting differential mass functions of halos were compared with Press-Schechter theory. First, we looked at how the parameters of the algorithms influence the mass function, and then how varying the box size of the simulation (and thus effectively changing the resolution of the halos) affected the results.

### 5.1. Varying the Parameters

Fig. 2 represents the differential mass functions in a typical case, here a 400 Mpc box with $`64^3`$ particles from the series of OCDM Hydra simulations. The linking lengths used are $`b=0.25`$, $`b=0.2`$, and $`b=0.1`$, in units of the mean particle separation. These values cover the range used in the literature. We see that the differential mass function obtained with FOF depends strongly upon $`b`$, and that the value $`b=0.2`$ gives good agreement with Press-Schechter theory. (This has been noticed earlier in  and  — see  for a comparison of FOF with virialized halos.) For FOF, the different values of the linking length effectively correspond to different values of the density contours, inside which particles are grouped together, and therefore we expect a strong dependence upon that parameter. Since FOF does not require the calculation of particle densities, $`b`$ is the only free parameter for this algorithm. In contrast, DENMAX and SKID show hardly any variation with the linking length. This is not surprising. During the first phase of these two algorithms the particles are moved towards the closest density maximum until their locations vary by less than the convergence radius. But the linking length of the FOF step during the second phase is twice the convergence radius. Hence changing $`b`$ also changes how tight the moved particles cluster around the density maxima. But the location of the density maxima is independent of the linking length, and so the FOF step scoops up the same particles into the same groups before the removal of gravitationally unbound particles takes place, which is based upon the original, unmoved locations.

Instead, DENMAX and SKID show a strong dependence upon the resolution at which the density gradients are calculated. (The resolution of the simulation stays the same throughout this discussion as the number of particles is kept the same.) For DENMAX, the grid size was increased from $`64^3`$ (the same as the number of particles) to $`192^3`$ (the largest size handled by the workstation on which the calculations were performed). For SKID, the number of nearest neighbors was reduced from 256 to 16. It can be seen from Fig. 2 that a low resolution of the density field (corresponding to a small grid or a large number of neighbors) produces preferably large halos, since there are few density maxima present inside the volume towards which all of the particles are moved. Each of the density maxima will then end up with a large number of particles, and hence larger groups are produced. On the other hand, with high resolution of the density field (i.e. large grids or small number of neighbors), it becomes bumpy, and, except for the densest regions, each particle will correspond to a density maximum at its location. Thus the particles will not move much during the first phase of DENMAX and SKID, and these algorithms become close to plain FOF. This can be seen clearly in Fig. 2 where the curves for DENMAX and SKID come closer to the Press-Schechter prediction (to which FOF with $`b=0.2`$ is very close as seen above) as the resolution of the density field increases. We also see an emergence of variation in the differential mass function with the linking length for the high resolution version of SKID with 16 neighbors, as we would expect if this algorithm became more similar to FOF.

Table 2 summarizes these results in numeric form. FOF with $`b=0.2`$ has the overall smallest value of the goodness of agreement $`\mathrm{\Delta }`$. This is the case even when compared to those values for SKID, where a smaller mass range for $`\mathrm{\Delta }`$ had to be used because of the lack of groups at the low-mass or high-mass ends. The values of $`\mathrm{\Delta }`$ in Table 2 for SKID with 256, 128, and 16 neighbors (there for $`b=0.1`$ only) are only lower limits because of the restricted mass ranges. For DENMAX and SKID we see, at the same grid size or number of neighbors, very little variation of $`\mathrm{\Delta }`$ with the linking length, but an improvement towards the values obtained with FOF as the resolution is increased.

To see more clearly where the large halos in low-resolution DENMAX and SKID originate we have isolated three typical halos in the 400 Mpc OCDM simulation, confining ourselves to a linking length of $`b=0.2`$ for the moment. Figs. 3 and 4 compare the group identifications of FOF with those of DENMAX and SKID, respectively. We see that low-resolution DENMAX and SKID indeed produce very extended halos, as argued above from the smaller number of density maxima in the simulation volume. The gravitational unbinding step in both algorithms does not seem to be able get rid of all of these outlying particles, which based upon the visual appearance should not be regarded as members of the halo. As can be seen from Figs. 5 and 6, it retains particles which have a low peculiar velocity with respect to the center of mass of the halo in these outlying regions. There is no shortage of such particles, as we expect the particles in these regions to have low peculiar velocities due to the shallow potential there. The situation improves for the higher resolution versions of DENMAX and SKID when these outlying particles are not considered to be group members in the first place because these two algorithms become closer to plain FOF.

The inclusion of outlying particles with low peculiar velocities in low-resolution DENMAX and SKID biases the velocity dispersion to lower values in these cases, as can be seen in Table 3. Thus the choice of group identification not only affects the masses of the halos, but also actual physical observables like the 1D velocity dispersion of a halo.

The same results are observed in the 400 Mpc box of the SCDM simulation.

### 5.2. Varying the Box Size

To compare the group finding algorithms at different resolutions of the simulations we looked at different box sizes (since the number of dark matter particles is kept the same). The best values for the parameters found in the discussion above are used, i.e. $`b=0.2`$ for the linking length of FOF and the highest resolution in calculating the density field for DENMAX and SKID. Since the choice of $`b`$ for DENMAX and SKID does not make much of a difference, the same value $`b=0.2`$ was chosen for these two algorithms.

The resulting differential mass functions are plotted in Fig. 7 for the OCDM Hydra simulations, and in Fig. 8 for the SCDM Hydra simulations. At visual inspection, the three algorithms seem to perform equally up to box sizes of 200 Mpc (for OCDM) and 400 Mpc (for SCDM), above which DENMAX and SKID start again to produce a larger number of massive halos as compared with FOF or Press-Schechter theory. At the low-mass ends, FOF and DENMAX show an overabundance with respect to Press-Schechter theory. But the smallest groups plotted in Figs. 7 and 8 consist of only two particles, and hence this overabundance is not worrisome as such small groups would be neglected in real applications anyhow.

The goodness of agreement $`\mathrm{\Delta }`$ with Press-Schechter theory, calculated over a mass range which excludes these smallest groups, confirms the trend in the plots. For the OCDM Hydra simulations, Table 4 shows that actually DENMAX has the best agreement up to 200 Mpc, after which FOF fares the best. SKID is always worse than FOF and, except for one case, even DENMAX. Similarily, Table 5 for the SCDM Hydra simulations prefers DENMAX up to 200 Mpc, SKID slightly for 400 Mpc, and FOF for 800 Mpc, when compared with Press-Schechter theory. For the 1600 Mpc simulations, FOF produced a very small number of halos, such that it was not possible to construct a reliable mass function in the OCDM model, and only a very unreliable one in the SCDM model. Therefore, the 1600 Mpc box has been omitted from the comparison in Tables 4 and 5.

## 6. Conclusions

The most basic property of a dark matter halo, its mass, depends upon the choice of group finding algorithm and its free parameters. Many of the other properties are likely to do so as well, and we have shown circumstantial evidence of this for the velocity dispersion. Much care has to be taken in choosing the right algorithm for the purpose in question. If the task is to reproduce the mass distribution as predicted by Press-Schechter theory (which might not be correct, though), the popular DENMAX and SKID algorithms have to be run with a very high resolution of the density field. For DENMAX this means a large grid (larger than the number of particles), and for SKID a small number of neighbors (smaller than 64 neighbors, which is the default number recommended by the authors of SKID ). Under these circumstances, DENMAX and SKID become close to a simple FOF algorithm, which, for the right linking length, generally gives the best agreement with Press-Schechter theory and at a much smaller computational cost. Similar effects are seen when instead of the resolution of the algorithm the resolution of the simulation is changed. With fixed number of particles, DENMAX and SKID perform worse, with respect to Press-Schechter theory, as the box size is increased, with FOF faring better. These results hold up for both cosmological scenarios studied.

In cases with low resolution, DENMAX and SKID produce more massive halos as compared with Press-Schechter theory and FOF. These halos are larger because DENMAX and SKID have the tendency of including distant particles which have a low peculiar velocity, a problem from which FOF does not seem to suffer. This would make FOF the preferable algorithm, and imply with its close agreement with Press-Schechter theory, that its predictions are doing quite well, too.

## Acknowledgments

We wish to thank Edmund Bertschinger for providing the P<sup>3</sup>M code, and Hugh Couchman and the Hydra Consortium and the University of Washington group of the High Performance Computing and Communication Program for making the N-body simulation and group identification codes available to the cosmology community on the World Wide Web. Further thanks go to Lars Hernquist, Ue-Li Pen, and David Weinberg for discussions and comments on a first draft of this paper. Part of the computational work in support of this research was performed at the Theoretical Physics Computing Facility at Brown University. MG is supported by an SAO Predoctoral Fellowship and a Manning Fellowship from Brown University.