Finding Halos with Bound Density Maxima code

Finding halos in dense environments is a challenge. The most widely used halo-finding algorithms - the friends-of-friends (e.g., Efstathiou et al.1985) and the spherical overdensity algorithm (e.g., Lacey & Cole 1994; Klypin 1996) - are not acceptable (Gelb & Bertschinger 1994, Summers et al.1995). The friends-of-friends (FOF) algorithm merges together apparently distinct halos if the linking radius is too large or misses some of the halos if the radius is small. Adaptive FOF (van Kampen 1995) seems to work better. But we find that it is difficult to find an optimal scaling of the linking radius with the density. We have developed a related algorithm (Klypin, Gottlober, Kravtsov 1997), which we call ``hierarchical friends-of-friends''. Because it uses all linking radii, it does not have the problem that the adaptive FOF algorithm has. The algorithms, either adaptive or hierarchical, can not work because they pick up many fake halos in very dense environments. Klypin et al.(1997) supplement the hierarchical FOF algorithm with an algorithm which checks if halos existed at previous moments. The algorithm which finds halos as maxima of mass inside spheres of a given overdensity works better than plain FOF, but no fixed overdensity limit can find halos in both low and high density environment. The DENMAX algorithm (Gelb & Bertschinger 1994) and its offspring, SKID (Governato et al.1997), make significant progress - they remove unbound particles, which is important for halos in groups and clusters. Recently, Summers et al. (1995) tried to perfect the idea of Couchman & Carlberg (1992) to trace the history of halo merging and to use it for halo identification. Starting at an early epoch, Summers et al.identify halos using the FOF algorithm with a linking radius corresponding to the ``virial overdensity'' of 200 and then trace particles belonging to halos at later times. It appears that it is impossible to make a working algorithm because halos interact too violently. A large fraction of mass is tidally stripped from some halos and a large fraction of mass is accreted. Some of the problems that any halo finding algorithm faces are not numerical. They exist in the real Universe. We select a few typical difficult situations.

1. A large galaxy with a small satellite. Examples: LMC and the Milky Way or the M51 system. Assuming that the satellite is bound, do we have to include the mass of the satellite in the mass of the large galaxy? If we do, then we count the mass of the satellite twice: once when we find the satellite and then when we find the large galaxy. This does not seem reasonable. If we do not include the satellite, then the mass of the large galaxy is underestimated. For example, the binding energy of a particle at the distance of the satellite will be wrong. The problem arises when we try to assign particles to different halos in an effort to find masses of halos. This is very difficult to do for particles moving between halos. Even if a particle at some moment has negative energy relative to one of the halos, it is not guaranteed that it belongs to the halo. The gravitational potential changes with time, and the particle may end up falling onto another halo. This is not just a precaution. This actually was found very often in real halos when we compared contents of halos at different redshifts. Interacting halos exchange mass and lose mass. We try to avoid the situation: instead of assigning mass to halos, we find the maximum of the ``rotational velocity'', $\sqrt{GM/R}$ , which is observationally a more meaningful quantity.

2. A satellite of a large galaxy. The previous situation is now viewed from a different angle. How can we estimate the mass or the rotational velocity of the satellite? The formal virial radius of the satellite is large: the big galaxy is within the radius. The rotational velocity may rise all the way to the center of the large galaxy. In order to find the outer radius of the satellite, we analyze the density profile. At small distances from the center of the satellite the density steeply declines, but then it flattens out and may even increase. This means that we reached the outer border of the satellite. We use the radius at which the density starts to flatten out as the first approximation for the radius of the halo. This approximation can be improved by removing unbound particles and checking the steepness of the density profile in the outer part.

3. Tidal stripping. Peripheral parts of galaxies, responsible for extended flat rotation curves outside of clusters, are very likely tidally stripped and lost when the galaxies fall into a cluster. The same happens with halos: a large fraction of halo mass may be lost due to stripping in dense cluster environments. Thus, if an algorithm finds that 90% of mass of a halo identified at early epoch is lost, it does not mean that the halo was destroyed. This is not a numerical effect and is not due to ``lack of physics''. This is a normal situation. What is left of the halo, given that it still has a large enough mass and radius, is a ``galaxy''.

We have developed our halo-finding algorithm (Klypin et al.1997) having in mind all these problems. The bound-density-maxima (BDM) algorithm stems from the DENMAX algorithm (Gelb & Bertschinger 1994). Just as in DENMAX, the BDM algorithm first finds positions of the density maxima on some scale and then removes unbound particles inside the halo radius. However, the algorithm finds maxima and removes unbound particles in a way different from DENMAX. The algorithm can work by itself or in conjunction with the hierarchical FOF. In the latter case, it takes positions of halos from the hierarchical friends-of-friends, and then removes unbound particles and finds parameters of halos.

In order to find positions of halos we choose a radius r_sp of a sphere for which we find maxima of mass. This defines the scale of objects we are looking for, but not exact radii or masses of halos. The radius of a halo can be either larger or smaller than r_sp , but distances between halos cannot be smaller than r_sp . We place a large number of spheres in the simulation box. The number of the spheres is typically an order of magnitude or more larger than the expected number of halos. For each sphere we find the mass inside the sphere and the center of the mass. The center of the sphere is displaced to the center of mass and the new mass and the center of mass is found. The process is iterated until convergence. Depending on specific parameters of the simulations, the number of iterations ranges from 10 to 100. This process finds local maxima of mass within sphere of radius r_sp .

The efficiency of finding local maxima of mass depends on how the spheres are chosen. In the present version of the code two algorithms were implemented. (1) A small fraction of all particles are chosen as centers of the spheres. The code asks you to enter N_seed , the ``Number of particles for initial seeds''. Then it will select every N<sub>particles</sub>/N<sub>seed</sub> particle as a center of a sphere. (2) Additional spheres will be placed in regions with relatively low density. The whole simulation box is divided into a mesh of large cells. The size of the cells is defined by the variable ``Cell'' in PMhalos.f. The ``Cell'' is typically equal to one or two PM cells. If a cell has many particles in it (``neighbors''), than some of them (not more than 3) will be chosen as centers of spheres. The code asks you to enter the minimum number of neighbors: ``Number of neighbors for a seed''.

In some cases one would need to improve the location of the halo. An example is if one is looking for groups of ``galaxies'' but also would like to have the groups always centered on a galaxy. The search radius for the groups may be chosen to be, say 500 kpc. Additional iterations with a smaller radius of the sphere will find the galaxy-size halo closest to the center of mass of the group and place the center of the group at the ``galaxy''. In the BDM code this option is realized in the following way. The code asks you to enter the ``smaller radius for final halos''. If this radius r_small is not equal to the search radius r_sp , the code will do additional iterations by gradually changing the search radius from r_sp to r_small . If r_small = r_sp , no additional iterations are made.

Some of the density maxima will be found many times because in the process of maximizing the mass some of spheres converge on the same local maximum. Spheres which find the same maximum are called ``duplicates''. We remove duplicates and keep only one halo for each maximum. Halos with too small number of particles (typically 5-10) and halos with too low central overdensity are removed from the final list. Parameters which control the removal are supplied by the user.

Once centers of potential halos are found, we start the procedure of removing unbound particles and finding the structure of halos. We place concentric spherical shells around each center. For each shell we find the mass of the dark matter particles, the mean velocity, the velocity dispersion relative to the mean, and the maximum of the rotational velocity V_max = $\sqrt{GM(r)/r}$$\mid_{max}^{}$ . In order to determine whether a particle is bound or not, we estimate the escape velocity at the distance r of the particle from the halo center:

$$\approx {\ln(1+2r/r_{max}) \over (r/r_{max})}$$ (25)

where r_max is the radius of the maximum of the rotational velocity. This expression for the escape velocity is valid for a halo with the Navarro-Frenk-White density profile. If the velocity of a particle is larger than the escape velocity, it is assumed to be unbound. We estimate the maximum rotational velocity V_max and radius of the maximum r_max using the density profile for the halo. Because V_max and r_max must be found before the unbound particles are removed and because the mean velocity is also found using all particles (bound and unbound), the whole procedure cannot be done in one step. We start by artificially increasing the value of the escape velocity by a factor of three. Only particles above the limit are removed. We find a new density profile, new mean velocities, and new V_max and r_max . The escape velocity is again increased, but this time by a smaller factor. The procedure is repeated 6 times. The last iteration does not have any extra factors for the escape velocity: all unbound particles are removed. Examples of halos identified by the code are presented in Klypin et al.(1997).

Finding a halo radius is straightforward for isolated halos: increase the radius of sphere untill the overdensity inside the sphere is equal to some limiting value provided by the user. For halos inside groups or clusters (halos inside halos) the algorithm is more complicated. It consists of three steps: (i) It starts with finding the radius of given overdensity limit (as for an isolated halo). (ii) Then, the algorithm checks how the mean overdensity inside given radius changes with the radius. It starts going from the very central shell outwards. If the mean overdensity stops declining, the algorithm assigns the radius to the halo radius. (iii) Now the algorithm goes from this radius inwards and checks the slope n of the overdensity profile: (M(r)/$\langle$M(r)$\rangle$) $\propto$ r ^{- n}. If the slope is too shallow (mass increases too fast with the radius), the radius of the halo is decreased because most of the mass at this radius does not belong to the halo. The radius is decreased untill the slope is steep enough. The limit for the slope is equal to n = 1 . This is a rather mild slope: for an isothermal sphere one expects n = 2 ; for the Navarro-Frenk-White profile n = 1.7 - 2.5 .

The BDM code will ask you to enter many parameters. A typical dialog may look as follows:

Enter Min. Center Overdensity for Halos = > 340. (1)

Enter Overdensity Threshold for Halos    = > 340. (2)

Enter Minimum halo mass in Msun/h     = > 5.e+9 (4)

Enter comoving search radius(Mpc/h)     = > 0.050 (5)

Enter smaller radius(Mpc/h) of final halos = > 0.030 (6)

Enter min.radius for halos(Mpc/h)         = > 0.030 (7)

Enter fraction of DM particles (1/4,1/2,1) = > 1 (8)

Enter rejection velocity limit (V/Vescape) = > 1.0 (9)

Distance to check for Velocity duplicates = > 0.300 (10)

Define duplicates if (v1-v2)/Vrms <           0.10 (11)

Enter Comoving Box size(Mpc/h)          = > 20. (12)

In lines (1-2), enter the minimum overdensity for halos. If both values are equal, the code will find halos with the overdensity above the limits you provided. You may choose only those halos which have higher central density, if you enter a larger number in the first line. Only halos with mass larger than the value entered in the third line and radius in the line 7 will be kept. For debugging of the code, you may read only 0.25, or 0.5 of all particles (line 8). If you would like to remove unbound particles as described above, enter 1 in line 9. In order to ignore the option (no removing of unbound particles) enter a number larger than 5. Lines 10 and 11 are used for two parameters needed for additional screening of duplicates. In case of very large halos, when the density profile in the center of the halo is rather flat, the iteration of spheres stops when the spheres are still far one from the other. The spheres have found the same object: masses, radii, velocities of the `` halos'' are very close, but their positions are slightly different. Because the fake halos have very close velocities, they can be identified. The parameter in the line 10 defines the maximum distance (in units of h ^{- 1}Mpc ) within which the code will look for the duplicates. If the difference in velocities $\sqrt{(v_{x1}-v_{x2})^2+(._y)+(._z)}$ is less than the parameter in line (11) as measured in units of maximum of the rms velocities for particles in the halos, only largest will be kept. You may ignore the option by entering two zeros.