arXiv Vanity renders academic papers from arXiv as responsive web pages so you don’t have to squint at a PDF. Read this paper on arXiv.org.

# Halo Abundance Matching: accuracy and conditions for numerical convergence

Anatoly Klypin, Francisco Prada, Gustavo Yepes,Steffen Heß, Stefan Gottlöber
, Astronomy Department, New Mexico State University, MSC 4500, P.O.Box 30001, Las Cruces, NM, 880003-8001, USA
Campus of International Excellence UAM+CSIC, Cantoblanco, E-28049 Madrid, Spain
Instituto de Astrof´ısica de Andaluc´ıa (CSIC), Glorieta de la Astronom´ıa, E-18080 Granada, Spain
Leibniz-Institut f¨ur Astrophysik (AIP), An der Sternwarte 16, D-14482 Potsdam, Germany
E-mail:
###### Abstract

Accurate predictions of the abundance and clustering of dark matter haloes play a key role in testing the standard cosmological model. Here, we investigate the accuracy of one of the leading methods of connecting the simulated dark matter haloes with observed galaxies – the Halo Abundance Matching (HAM) technique. We show how to choose the optimal values of the mass and force resolution in large-volume -body simulations so that they provide accurate estimates for correlation functions and circular velocities for haloes and their subhaloes – crucial ingredients of the HAM method. At the 10% accuracy, results converge for 50 particles for haloes and 150 particles for progenitors of subhaloes. In order to achieve this level of accuracy a number of conditions should be satisfied. The force resolution for the smallest resolved (sub)haloes should be in the range , where is the scale radius of (sub)haloes. The number of particles for progenitors of subhaloes should be . We also demonstrate that the two-body scattering plays a minor role for the accuracy of -body simulations thanks to the relatively small number of crossing-times of dark matter in haloes, and the limited force resolution of cosmological simulations.

###### keywords:
cosmology: large-scale structure of the Universe – cosmology: theory – methods: numerical

## 1 Introduction

The standard CDM  cosmological model phases significant challenges when it comes to testing its predictions for the distribution and properties of galaxies. The model is able to make detailed predictions on the distribution of dark matter, but connecting the luminous galaxies with their dark matter haloes is a difficult task. There are different possibilities to make this galaxy-halo connection. Halo Abundance Matching (HAM) is a simple and yet realistic way to bridge the gap between dark matter haloes and galaxies (Kravtsov et al., 2004a; Tasitsiomi et al., 2004; Vale & Ostriker, 2004; Conroy et al., 2006; Guo et al., 2010; Trujillo-Gomez et al., 2011; Reddick et al., 2012; Kravtsov, 2013). HAM resolves the issue of connecting observed galaxies to simulated haloes and subhaloes by setting a correspondence between the stellar and halo masses: more luminous galaxies are assigned to more massive haloes. By construction, the method reproduces the observed stellar mass and luminosity functions. It also reproduces the galaxy clustering over a large range of scales and redshifts (Conroy et al., 2006; Guo et al., 2010; Trujillo-Gomez et al., 2011; Reddick et al., 2012; Nuza et al., 2013). When HAM is applied using, for example, the observed SDSS stellar mass function (Li & White, 2009), it gives also a reasonable fit to lensing results (Mandelbaum et al., 2006), to the galaxy clustering and the relation between stellar and halo virial masses (Guo et al., 2010). Trujillo-Gomez et al. (2011) show that accounting for baryons drastically improves the match of predicted and observed Tully-Fisher and Baryonic Tully-Fisher relations. HAM modeling was also successful in reproducing other properties of galaxies (Behroozi et al., 2010; Leauthaud et al., 2011a; Reddick et al., 2012; Hearin et al., 2012; Kravtsov, 2013).

Yet, using the HAM technique to connect galaxies to (sub)haloes requires substantially more accurate -body simulations as compared with those used for the Halo Occupation Distribution (HOD) model (Kravtsov et al., 2004a; Tasitsiomi et al., 2004; Vale & Ostriker, 2004). Because of the demanding numerical requirements, there are relatively few studies of abundance and clustering of haloes based on circular velocities, which are often used in conjunction with HAM models (Klypin et al., 2011; Trujillo-Gomez et al., 2011; Reddick et al., 2012; Guo & White, 2013). Accuracy of these statistics was challenged by Guo & White (2013), who found a very poor numerical convergence of results for subhalo population in the Millennium I (Springel et al., 2005) and Millennium II (Boylan-Kolchin et al., 2009) simulations. Though we find much better convergence in our Bolshoi and suite of MultiDark simulations, we agree with Guo & White (2013) that it is important to investigate the accuracy and limits of HAM technique to connect galaxies with (sub)haloes using only -body simulations. In particular, physically robust results demand that statistics such as the abundance and clustering of subhaloes to be unaffected by numerical resolution.

Many large volume -body cosmological simulations are not suited and should not be used for HAM models. The numerical and physical processes which affect the accuracy of results based on cosmological simulations were discussed in many publications (e.g., Knebe et al., 2000; Klypin et al., 2001; Power et al., 2003; Hayashi et al., 2003). It is one of the goals of this paper to find conditions, that a simulation should pass in order to be used for HAM models.

The key ingredient of HAM models are subhaloes: satellites (subhaloes) of distinct haloes must be part of the abundance matching prescription because each lump of dark matter with enough mass and concentration should host a galaxy regardless whether that is the central object or a satellite. We call a halo distinct if its center is not inside the virial sphere of even a larger halo. By definition, a subhalo is always within the virial radius of another halo, which in this case is called parent halo. In the sense of dynamical evolution subhaloes are different objects because their physical properties can be significantly affected (mostly through tidal stripping) by their parent haloes (Tasitsiomi et al., 2004; Kravtsov et al., 2004b). In reality, the boundary separating distinct haloes and subhaloes is blurry. First, there are different definitions of virial radius. As a result, the same object can be called distinct or subhalo depending on the virial ghalo definition. Second, distinct haloes may experience strong interaction with their environment long before they formally cross the virial radius of their parent (Kravtsov et al., 2004b; Behroozi et al., 2013a). In that respect some distinct haloes evolve as subhaloes.

There are different flavors of HAMs. Generally, one does not expect a pure monotonic relation between stellar and halo masses. There should be some degree of stochasticity in this relation due to deviations in the halo merger history, angular momentum, and concentration. Even for haloes (or subhaloes) with the same mass, these properties should be different for different systems, which would lead to deviations in stellar mass. Observational errors are also responsible in part for the non-monotonic relation between halo and stellar masses. Most of modern HAM models already implement prescriptions to account for the stochasticity (e.g., Behroozi et al., 2010; Trujillo-Gomez et al., 2011; Leauthaud et al., 2011a). The difference between monotonic and stochastic models depends on the magnitude of the scatter and on the stellar mass. The typical value of the scatter in the -band is expected to be mag (e.g., Trujillo-Gomez et al., 2011). For the Milky Way-size galaxies the differences are practically negligible (Behroozi et al., 2010).

Because haloes may experience tidal stripping, their dark matter mass may be significantly reduced (e.g. Klypin et al., 1999a; Hayashi et al., 2003; Kravtsov et al., 2004b; Arraki et al., 2012). Galaxies hosted by these stripped haloes are not expected to lose their stellar mass to the same degree as the dark matter mass because the stellar component is much more concentrated and thus less susceptible to tidal forces. This is the reason why the halo dark matter mass is expected to be a poor indicator of the stellar mass of galaxy hosted by the halo. There are quantities that should work better for HAM and are often used as proxies for stellar mass: dark matter mass before the stripping started, maximum circular velocity of the dark matter or the maximum circular velocity of the dark matter before the stripping (e.g., Conroy et al., 2006; Trujillo-Gomez et al., 2011; Reddick et al., 2012).

Here we use the maximum circular velocity, not mass. The maximum is reached in the central halo region, which is expected to correlate better with the stellar or luminous component and it is less sensitive to tidal stripping. This is the main reason why in this paper we focus our numerical convergence study on the abundance and clustering statistics related with the (sub)halo circular velocities.

While the maximum circular velocity is a better quantity to characterize (sub)haloes (Conroy et al., 2006; Klypin et al., 2011; Trujillo-Gomez et al., 2011), it is more difficult to accurately measure it in numerical simulations. For example, for galaxy-size haloes the maximum of the circular velocity happens at of the virial radius. At this radius, the resolution should be sufficient to estimate the maximum circular velocity with better than, say, an accuracy better than 10 percent. In addition, subhaloes also should be resolved and their circular velocities are estimated. Thus, the resolution of simulations intended for HAM models should be at least an order of magnitude better than those intended, for example, for HOD models.

Our paper is organized as follows. In Section 2 we present simulations used for our analysis and discuss halo identification. Section 3 gives results on the abundance of haloes and subhaloes. The numerical convergence study of the correlation function for (sub)haloes selected by the maximum circular velocity is discussed in Section 4. In Section 5 we investigate numerical and physical processes, which affect accuracy of cosmological -body simultions. Here we also provide the conditions for the numerical convergence of HAM results. Our results are summarized in Section 6

## 2 Simulations and dark matter Haloes

In order to study the numerical convergence of cosmological -body simulations we use three large simulations, which cover a mass range of five orders of magnitude. The same cosmological parameters were used for the simulations. The simulations were performed using two codes: the Adaptive Refinement Tree (ART Kravtsov et al., 1997; Gottlöber & Klypin, 2008) and Gadget (Springel, 2005). Details of the ART simulations Bolshoi and MultiDark were given in Klypin et al. (2011) and Prada et al. (2012). The simulation done with the Gadget is from a large suit of BigMultiDark simulations(Hess et al., 2013). Main parameters of our simulations are presented in Table 1, where we also present parameters of the two Millennium simulations (Springel et al., 2005; Boylan-Kolchin et al., 2009). The table highlights one substantial difference between our simulations and Millennium: our force resolution is substantially better for the same mass resolution. Another difference, which is not shown in the table is that we systematically use significantly smaller time-steps.

The choice of numerical parameters (combination of mass, force, and time resolutions) in our simulations is not accidental. We made many tests (some of them discussed below) studying the convergence of results.

We use the Bound Density Maxima (BDM) spherical overdensity code (Klypin & Holtzman, 1997; Riebe et al., 2011) to identify and find the properties of haloes and subhaloes in all our simulations. The code was extensively tested and compared with other halofinders (Knebe et al., 2011; Behroozi et al., 2013b; Knebe et al., 2013). Halo tracking was made for two simulations (Bolshoi and MultiDark) using the ROCKSTAR code (Behroozi et al., 2013b). No halo tracking so far is available for the bigMultiDark simulation.

The BDM halo finder provides many parameters for each (sub)halo including the virial mass, concentration, and spin parameter. It also includes the maximum circular velocity . For Bolshoi and MultiDark we also find the peak of over the history of evolution of each halo and subhalo. This is done using the information over the history of the major progenitor of each halo. The peak velocity serves as proxy for the rotation velocity of the galaxy hosted by its halo. This is why plays a key role for HAM models. However, from the point of view of the accuracy of simulations, is a more important quantity with simply considered to be mapping .

The virial radius for distinct haloes is defined as radius within which the mean matter density is equal to the virial overdensity found from the solution of the top-hat model for cosmological model with a cosmological constant. For our set of cosmological parameters this corresponds to a matter overdensity of 360 times the average matter density in the universe.

When studying convergence of numerical results, we quote the number of particles within the virial radius. Subhalos normally do not extend to virial radius because their radius is truncated by tidal stripping. In order to make the number of particles in haloes and subhalos compatible, for subhalos we quote the number of particles inside the virial radius of a distinct halo with the same maximum circular velocity as that of the subhalo. The extrapolation to the virial radius depends on halo concentration. We use the average virial mass - relation for distinct halos to make the extrapolation to the virial radius. The actual number of particles in a subhalo can only be smaller than the extrapolated value. We find that the extrapolated virial masses of subhalos are within 10–20% of virial masses of 90% of progenitors of the subhalos defined at the moment of the peak of . This is why we call these estimates of the subhalo masses and number of particles the masses and particles numbers of subhalo progenitors.

## 3 Halo velocity function

We start with the analysis of the velocity functions using the peak values of in Bolshoi and MultiDark. Figure 1 shows the velocity functions separately for subhaloes and the total populations of haloes. Overall there is an excellent agreement between both simulations for a wide range of circular velocities , where (sub)haloes are resolved with enough particles and the cosmic variance does not effect the results. At large velocities ( there are clearly some effects due to the lack of long waves and cosmic variance in the Bolshoi simulation due to its smaller box size. Deviations of the number of haloes at small velocities are indication of numerical resolution effects in the MultiDark simulation. The smaller number of haloes in Multidark at smaller velocities are due to the lack of force and mass resolution as compared to Bolshoi (see Table 1).

At the 10% level of accuracy the number-density of haloes converge for peak velocities , which corresponds to the average number of 52 particles inside the virial radius in the Multidark simulation. As expected, convergence is worse for subhaloes. However, it is still very good: at 10% it is . For isolated haloes this corresponds to an average of 150 particles. This is a remarkably small number of particles. For comparison, convergence of the velocity function of subhaloes in the Millennium simulations, at the same accuracy level, was achieved only with particles (see Guo & White, 2013)

Accuracy of HAM results depends on a number of numerical effects. Each of them should produce small enough errors so that the final result can be trustful. The effects include the tracking of subhaloes as they fall into larger haloes and become subhaloes (Behroozi et al., 2010). Other effects include mass and force resolution, time-integration of equations of motion, and two-body scattering. Accuracy of halo tracking was discussed in Behroozi et al. (2010). Here we focus on the other effects. In order to study those, we investigate the statistics based on the instantaneous maximum of the circular velocity .

Figure 2 presents the results of three simulations: Bolshoi, MultiDark, and BigMultiDark. To large degree the results are very similar to those based on . Again, the comparison between Bolshoi and MultiDark shows convergence at 10% level for distinct haloes with and for subhaloes with . This implies that subhaloes need about 150 particles to have this accuracy. The comparison between BigMultiDark and MultiDark shows the same, i.e. in order to reach the 10% accuracy, one needs particles for subhaloes.

## 4 Convergence of the correlation function

Convergence of the two-point correlation function for haloes and subhaloes, selected by either or , is demonstrated in Figures 3 and 4, where we compare results for Bolshoi and MultiDark. Just as for the circular velocity functions, the agreement of the results on small scales () indicates convergence of the results regaring the numerical parameters such as mass and force resolution. The differences in correlation functions on larger scales are due to finite box size and cosmic variance.

On scales the agreement between Bolshoi and MultiDark is about 5%, with Bolshoi exhibiting slightly larger clustering signal as expected due to it better force and mass resolution (see Table 1). The situation is unclear on even smaller scales mostly because of very small number of pairs within these separations (and thus large noise) in Bolshoi.

Because of much larger statistics of halos in MultiDark and BigMultiDark simulations, the level of shot noise on scales is smaller, which allows us to test those simulations on these scales. Indeed, in this case the converge of results is substantially better as illustrated by Figure 5. It is also much better on large scales due to sugnificantly larger box sizes. Both simulations have nearly the same (within few percent) correlation functions from   to . One starts to run on substantially larger differences on scales below . This is due to a combination of force resolution and real-space halofinder, which misses subhalos, when their centers are too close to the center of parent halo.

Another interesting statistics, which can be used to probe the numerical convergence, is the halo-dark matter cross correlation function . It is an important quantity on its own because it is required for theoretical predictions of the weak galaxy lensing signal (e.g., Leauthaud et al., 2011b). Figure 6 shows for MultiDark and BigMultiDark with halos selected at the low limit of resolution: . The difference between the simulations is less than % for scales larger than . There is a 5–10% error on scales with the lower-resolution BigMultiDark simulation predicting stronger clustering signal. By using simple -body toy models we find that the reason for this extra clustering is due to a side-effect of force resolution. Because of lack of resolution, the dark matter, which should have been at small distances, ends up at larger distances, where it slightly increases the local density and, thus, increases the cross correlation signal.

On scales above the correlation function in Bolshoi is systematically lower than in MultiDark. The large-scale differences between Bolshoi and Multidark simulations are likely due to both the cosmic variance and the finite box size. It is easy to test the effect of box size. For that, we estimate the correlation function using linear theory and truncate the power spectrum at the lower wavenumber limit given by the fundamental mode of the simulation box: , where is the length of computational box. Figure 7 shows correlation functions of the dark matter, which are obtained for different box sizes.

As expected, the finite size of the simulation box affects the large scales in a profound way. For example, for the box used in the Millinium-II simulation (Boylan-Kolchin et al., 2009), the correlation function becomes negative for while it should be positive. Even smaller scales are affected: at the correlation function is 30% below the true one (solid cirve in Figure 7. Note that often used a rule of thumb that scales smaller than 1/10 of the box length are not affected by the box size, does not really hold. Indeed, the deviations of the linear correlation function at strongly depend on the box length . One may expect that nonlinear effects exacerbate the differences.

Because of the box size effects, the comparison between simulations on large scales should be done using a large box size. Indeed, the converge of results is substantially better for large-box simultions. Figure 5 illustrates this point by comparing MultiDark and BigMultiDark simulations with Gpc and Gpc. Both simulations have nearly the same (within few percent) correlation functions from   to for haloes and progenitors of subhaloes with particles.

Convergence and accuracy of halo clustering and abundances in our set of simulations are significantly better – by a factor of 20 in mass – than for the Millennium simulations (Guo & White, 2013). As a way to compensate the lack of resolution in the Millennium simulations, a trick is used: if a subhalo is lost, it is assigned to the most bound dark matter particle of this subhalo, and called an “orphan”. A simplified dynamical friction estimate is then used to set a clock when the subhalo is assumed to merge with its parent (Guo et al., 2011). It is not clear to us what causes such a poor performance of the Millennium simulations. This is not related with the -body code: we achieve satisfactory results aslo with the Gadget code. There are some concerns regarding the Millennium simulations, but none of them sepearately seems to be enough to explain the poor convergence in the Millennium simulations. For example, the number of time-steps is low. For the Millennium-I run the parameter ErrTolIntAccuracy, which defines the number of steps, is twice larger than what we use in our Gadget simulation bigMultidark. The force resolution is also low in the Millennium simulations. In addition, the small size of the computational box of the Millennium-II run raises doubts that it was large enough to provide the true solution when compared with the Millennium-I. It is possible that a combination of all these small defects resulted in much worse results for subhaloes. We highlight that the better accuracy of our simulations eliminates the need for “orphans” in our simulations.

## 5 Conditions for convergence

Usually, numerical convergence of cosmological simulations is discussed in the context of simulations of individual dark matter haloes with large number of particles. For example, convergence of dark matter density profiles is studied by running simulations with increasing very large number of particles (e.g. Klypin et al., 2001; Springel et al., 2008; Stadel et al., 2009). Here, we are interested in convergence on the opposite side of the spectrum: how reliable are results with very small number of particles? The small haloes are very important. After all, most of haloes in each simulation are tiny, and numerous statistics (e.g. the halo mass functions) use those small halos. Results presented in Section 3 indicate convergence of the velocity function of distinct haloes with particles. How can we understand this?

The accuracy and convergence of results used for halo abundance matching depend on a number of physical and numerical processes. We split those into two categories: (a) accuracy of the maximum circular velocity for isolated haloes and (b) physics and numerics of subhaloes. We focus on some rather specific problems, which to large degree define convergence of numerical results. The goal is to study these processes separately by isolating each one and investigating it by using a simple realistic model. In reality, all processes are interconnected and only a full-scale cosmological simulation have all all of them included. However, understanding the conditions for convergence of large cosmological simulations is difficult because it is not clear what are the processes and how they potentially affect the results. We get much better insights by studying separate processes.

### 5.1 Accuracy of circular velocities of isolated haloes

Here we consider three effects: (a) force resolution, (b) mass resolution (number of particles) and (c) two-body scattering. We investigate a series of simulations of isolated NFW halos. In all cases we start with a NFW halo, which is set in equilibrium. The equilibrium initial conditions are constructed assuming an ideal phase-space distribution function for a NFW halo with isotropic velocities and with no correction for force softening or discreteness. Because of different numerical effects – the force softening, discreteness effects, and time-stepping – the system of particles does not stay in the equilibruim. It evolves. In most of the cases, the system reaches another equilibrium after few dynamical times. By comparing this new state with the initial one, we estimate the magnitude of the numerical effect responsible for the evolution.

When dealing with the NFW halo, it is convenient to express all physical quantities using the scale radius and the maximum circular velocity . Those quantities are related to the virial mass and concentration through the following relations:

 V2max = GMvirRvirf(xmax)f(c)cxmax, (1) c = Rvirrs,x≡rrs,xmax=2.163 (2) M(x) = Mvirf(x)f(c),f(x)=ln(1+x)−x1+x, (3) V2(x) = V2maxxmaxxf(x)f(xmax) (4) ρ(x) = V2max4πGr2sxmaxf(xmax)1x(1+x)2. (5)

Written in this way, the halo structure does not depend on and . This allows us to simulate just one system and then to rescale it to any particular values of and .

It is convenient to use some fiducial virial radius and virial mass. For this, we use concentration and fiducial virial radius which are typical parameters for galaxy-size haloes with mass (Prada et al., 2012). A useful scale for time is the crossing-time at the radius , at which the circular velocity reaches the maximum :

 tcross≡rmaxσv(rmax)≈xmaxrsVmax, (6)

where is the 3D r.m.s. velocity at . The crossing-time in physical units depends on halo concentration. For virial mass defined at overdensity , where being the critical density of the Universe, the crossing-time is:

 tcross=9.8h−1f(c)c3/2Gyrs. (7)

He have  Gyr for a Milky Way-type galaxy with . The crossing-time depends only very weakly on halo mass: it is twice smaller for haloes of dwarf galaxies with and twice larger for clusters of galaxies with . Note that haloes do not live long: during 10 Gyrs of evolution they have only crossing-times.

Initial density profiles extend to , which is significantly larger than any realistic halo in cosmological simulations. All simulations use a simple direct summation code with a leap-frog integration scheme. The code adopts the Plummer softening (force resolution) . All simulations were done with time-steps small enough to have energy conservation .

As a test, we run a simulation with large number of particles and small force softening . With the exception of very small non-equilibrium effects at the center, the halo did not change for many dynamical time-scales. However, when simulations are run with either small or large , the haloes do get modified. The comparison of the evolved halo with its initial equilibrium distribution gives as a way to measure effects of mass and force resolution separately. In cosmological simulations haloes grow in a complex fashion through merging and accretion. However, at late stages of evolution the accretion slows down and haloes spend most of their time in the regime of slowly growing quasi equilibrium NFW distribution. Ability of code to maintain this equilibrium distribution is an important indication of its accuracy.

#### 5.1.1 Force resolution

We first study the effects of force resolution by studying relatively large numerical experiments with 50,000 particles within a radius of and 100,000 particles in total. Within the radius the halo has about 15,000 particles. This configuration was simulated with different force resolutions ranging from to . Each simulation started with exactly the same initial conditions of equilibrium NFW halo without any consideration for the force resolution. For simulations with small there was very little evolution in the distribution of particles. However, in the cases of large the force was too weak in the central region to keep particles in equilibrium. As the result, each simulation with large evolved starting with the center. All simulations with settled into new equilibrium after few crossing-times. The simulation with was changing even after 5 crossing-times. At this moment the simulation was stopped because its very peripheral regions with were affected.

Figure 8 shows circular velocity profiles for runs with different force resolutions . Qualitatively the results are similar to those found in convergence studies of individual haloes in cosmological simulations: as the resolution increases, halo circular velocity profiles converge.

We can use these results to estimate errors in . The top panel in Figure 8 shows the dependance of on . At small the error in is relatively small. For example, we have for . However, the error of significantly increases with increasing : for . One must keep in mind that should be quite accurate because the relevant statistics, which use , are very sensitive to it. For example, the abundance of haloes and subhaloes scales with the maximum circular velocity as .

It is interesting to compare the convergence results in these experiments with those in our cosmological simulations. Our results indicated the convergence of the halo velocity function at the 5% accuracy for haloes with in the BigMultiDark simulation. This corresponds to haloes with 120 particles and typical scale radius . The force resolution is for the BigMultiDark. This gives . According to Figure 8 the error in is less than 2%. In turn, the error in the velocity function should be less than 5%, which is close to what we find in our cosmological simulations.

The following model can be used to make an approximation for the errors in . We note that profiles of shown in the bottom panel of Figure 8 look similar to the NFW profile, which scale radius is larger and is smaller than in the initial NFW halo. At large radii all the velocity curves converge to the same values because the mass distribution at large distances is not affected by the resolution. The increase in the scale radius due to the force resolution can be approximated as

 rs,ϵ=β(ϵ)rs,β(ϵ)≈√1+(1.7ϵ/rs)2. (8)

Using this expression and normalizing the modified NFW profile so that at large radii it has the same mass as the initial NFW, one obtains an approximation for the decline of due to the force resolution:

 (Vmax,ϵVmax)≈β(ϵ)−1/4=(1+(1.7ϵ/rs)2)−1/4. (9)

This approximation is shown in the top panel of Figure 8.

The increase in the scale radius due to poor force resolution leads to a decrease in halo concentration. Equation (8) implies that concentration measured for a halo, which was simulated with the force resolution , declines as . For this gives 8% error in concentration. If equation (1) is used to find halo concentration from measured and virial velocity (e.g. Prada et al., 2012), then reqirements to the force resolution are even more stringent. Error analysis shows that in order to have an error in concentration less than 5%, the force softening must be smaller than , which in turn implies not more than 1% error in .

Equation (8) can be used to recover the true concentration . We assume that a simulation, which is performed with the force resolution , provides the virial radius and the concentration for a halo. The true concentration for the halo can be estimaed as:

 c=cϵ√1−(1.7ϵcϵ/Rvir)2. (10)

These approximations can be used to set the upper limit on the force resolution needed to achieve a given accuracy of the maximum of the circular velocity . For example, in order to have no more than 2(5)% error in , the force resolution should be smaller than .

#### 5.1.2 Mass resolution

Here we address two issues: (a) are there any systematic effects that force the value of , measured in simulations with small number of particles, to deviate from its true value and (b) what is the level of statistical fluctuations of .

In order to address these issues, we made three simulations of NFW halos, which were initially in equilibrium, each with 300 particles inside radius of . These simulations initially had only 150 particles within the fiducial virial radius and 25 particles inside . This is the smallest number of particles, for which our cosmological simulations show convergence. Results presented below are averages of the three realizations. The force resolution is , which is typical for our cosmological simulations for the smallest resolved haloes. Results presented in Sec. 5.1.1 indicate that this force resolution is sufficient for accurate estimates of .

Figure 10 shows the circular velocity profiles in these simulations with very small number of particles. We can only estimate the velocity profile for distances larger than : there is not enough particles to probe smaller radii. However, on average the models stay close to the initial NFW profile: deviations are less than 5% for radii and for time .

The simulations indicate that the value of does not show any systematic changes over the same period of time and that even with 150 particles one can accurately estimate . At later moments gradually increases presumably due to accumulation of the two-body scattering effects. However, the change is small: 5% after 200 crossing times.

In order to measure the level of statistical fluctuations of , the value of is estimated for each snapshot and the statistics is accumulated for many snapshots. Results are presented in the top panel of Figure 10. The level of fluctuations of is just 3%, which is remarkably small considering that there are only 50 particles inside the radius , which typically defines the maximum of the circular velocity. One would naively expect larger fluctuations of . The real fluctuations are about five times smaller.

Assuming that the number of particles inside a given radius is distributed according to the Poissonian statistics, we expect that the r.m.s. fluctuations of circular velocity are , which gives . Additional reduction in the fluctuations is related with deviations from the Poissonian noise. Fluctuations of the number of particles are actually not expected to be Poissonian. One can design few examples to illustrate this. Orbital eccentricity should be a factor: there are no fluctuations, if the orbits are circular. For a given radius there are always some number of particles, which do not have energy to leave the radius. These particles always stay inside the radius and, thus, their number does not fluctuate. The fraction of these particles depends on the radius. In the central region, , a large fraction of particles is not bound to the center. They travel to large distances and only occasionally are found within the radius . At large distances the fraction of bound particles increases resulting in smaller fluctuations.

To measure the deviations from the Poissonian statistics, we run a number of different numerical experiments trying to estimate the r.m.s. fluctuations at different radii. The results, when expressed in units of the Poissonian noise, do not depend on how many particles are used for a particular simulation. Figure 9 shows the results for a simulation with 2,000 particles, which was run for 200 crossing-times. The fluctuations in the number of particles are close to Poissonian for small radii. At radius the r.m.s. fluctuations are about 1/2 of those expected for the Poissonian distribution. Combining all the effects, we estimate that the r.m.s. fluctuations of the maximum circular velocity due to uncorrelated statistical fluctuations of the number of particles can be estimated as:

 ΔVmaxVmax=14√N, (11)

where is the number of particles within a sphere of radius . For this approximation gives , which is what we found in our simulations.

#### 5.1.3 Two-body scattering

Effects of two-body scattering for the NFW haloes were extensively studied for haloes with relatively large number of particles (Hayashi et al., 2003; Power et al., 2003; Valenzuela & Klypin, 2003; Diemand et al., 2004). The smallest haloes studied had particles (Hayashi et al., 2003; Diemand et al., 2004). The main focus of these studies was the inner cusp of the haloes: to what degree the scattering can affect the inner structure of haloes. Instead, our analysis is focused on very small haloes with particles. Our main interest is the stability of the region which defines . All haloes – even the largest ones – go through a stage when they have very few particles. Hence, it is important that two-body scattering does not affect the density of these small haloes.

We start with analytical estimates of the two-body scattering, which we write using the notations introduced before. The time-scale for the scattering is (Binney & Tremaine, 2008, equation (7.106)):

 trelax=0.34σ3G2mρlnΛ, (12)

where is the 1d r.m.s. velocity of particles, is the particle mass, and is the Coulomb logarithm. For the NFW haloes with isotropic velocities one can use the Jeans equation to find :

 [σ(r)Vmax]2=xmaxf(xmax)x(1+x)2∫∞xf(x)dxx3(1+x)2. (13)

The particle mass can be parametrized using the number of particles inside the radius . Combining these expressions and writing the two-body scattering time in units of the crossing time (eqs.6 and 7), we get:

 trelaxtcross=0.344πNmaxlnΛ[x(1+x)2]5/2x3/2maxf1/2(xmax)[∫∞xf(x)dxx3(1+x)2]3/2. (14)

The integral in equation (14) may be taken numerically. However, it can be approximated in a simple way. We note that the relaxation time depends on the course-grained phase-space density , which is know to be nearly a power law. Indeed, the following approximation is accurate within 15% for a wide range of radii :

 [x(1+x)2]5/2[∫∞xf(x)dxx3(1+x)2]3/2≈x1.927.5. (15)

Substituting this into equation (14) and estimating numerical factors, we get:

 (16)

The term is the ratio of the maximum to the minimum impact parameters for the two-body scattering problem. In realistic cosmological simulations and in our simple numerical models the minimum impact parameter is defined by the force resolution. For we use the distance at which the force becomes Newtonian: . For one may use the whole region, which is relevant for the inner structure of the NFW distribution . This leads to .

We can get another useful approximation for the relaxation time, if in equation 16 we replace with its value given by equation (7) and also use the number of particles inside a given radius: . The following aprroximation gives for a given halo concentration and a number of particles :

 trelax≈1.4h−1GyrsN(r)cln(3rs/ϵ)[1+(rrs)5/4]. (17)

We can use equations (16) or (17) to estimate for our simulations. For simulations presented in Figure 10, and one finds , which is  Gyrs when scaled to Milky Way-size haloes ().

Note that during the halo density profile does not change much: within % the circular velocity profile is the same as the initial NFW (see Figure (reffig:TwoBody2)). The halo does suffer from the scattering, but on much longer time-scale of . To some degree, the scattering in these simulations was suppressed by choosing a reasonable force resolution, i.e. small enough to resolve well the halo and not too small to avoid excessive scattering. Indeed, we run the same simulations with much smaller force softening and find that scattering significantly increases. Figure 11 shows the results for . It is clear that in this case the two-body scattering affects the whole halo. However, even for this unrealistic case the scattering changes by 5% for the life-span of of these haloes in cosmological runs.

As equation (16) shows, the relaxation time-scale very strongly depends on radius . This makes the averaging of over the whole halo not very useful for haloes with steep cusps (Quinlan, 1996; Hayashi et al., 2003). E.g., a long time-scale of scattering estimated at half-mass radius does not imply that the scattering is negligable because it still may be important close to the halo center. The opposite is also not true: the fact that the scattering time is short at some distance does not mean that the other parts of the halo are affected.

All the results present above indicate that the effects of two-body scattering are localized, with the outer halo regions being basically collisionless. This distance, at which the scattering starts to make an effect may or may not be important for a particular situation. To illustrate the point we made numerical experiments with larger number of particles, which allow us to probe small radii. We simulated three realizations of haloes with 50,000 particles inside the fiducial virial radius of and particles in total. The force resolution was . Figure 12 shows that in these haloes the central region is affected by the two-body scattering while the outer regions are not. As time goes on, the radius separating the inner affected region from our collisionless halo gradually increases.

Note that at the radius containing the first 50 particles the decline in circular velocity is worse for these haloes with larger number of particles as compared with the previously studied haloes with only 150 particles inside the fiducial virial radius. This increase in the scattering at small radii can be shown analytically using equation (17). Indeed, the number of particles inside radius with a given value of must be times larger for than for . In that sense, the two-body scattering has the smallest impact for small haloes. For example, to have a halo should have particles inside and particles inside .

Our results are in overall agreement with those of Hayashi et al. (2003) when we compare their simulation of a NFW halo with 3,000 particles and take into account differences in the definitions of the crossing-time and the Coulomb logarithm. Just as Hayashi et al. (2003), we also find that two-body scattering starts to modify the resolved density profile on a time-scale, which is many times longer than the formal analytical estimates. Our conclusions are in strong disagreement with those of Diemand et al. (2004), who argue that the two-body scattering so much affects small haloes that the whole hierarchical growth of haloes is collisional. This strictly contradicts our results: the two-body scattering does not affect much the value of and the resolved parts of small haloes. A combination of two factors significantly reduces the role of the scattering: (a) small life-time of cosmological haloes (), and (b) low force resolution of cosmological simulations.

Our numerical results are also in a very good agreement with those of Power et al. (2003), but we disagree on the interpretation and final conclusions. Results presented in Figure 12 show that the circular velocity deviates from the analytical solution by less than % for central region containing 1,000 particles. This is in agreement with Figure 14 in Power et al. (2003), which shows the mean halo density contrast as a function of the enclosed number of particles. However, this is somewhat misleading because both statistics (circular velocity and mean enclosed density) are integral characteristics. The plot of the density at a given radius (a differential statistics) presented in Figure 13 shows that the density is not affected by the scattering already at radius containing particles. This agrees with the analysis of convergence of density profiles in cosmological simulations presented by Klypin et al. (2001).

It is difficult to justify an analytical approximation for the two-body scattering time used by Power et al. (2003, equation (20)). It gives the impression that it is a slightly modified form of the Chandrasekhar formula, but it is actually not. At the end the main difference is related with the r.m.s. velocities. In the Chandrasekhar approximation for two-body scattering (see equation (14)) the term is the r.m.s. velocity of the particles, which in the inner halo regions is substantially larger than the circular velocity given by equation (4). Instead, Power et al. (2003, equation (20)) use the circular velocity. This results in a substantial overestimate of the two-body scattering time and in a different scaling relation. One may treat Power et al. (2003, equation (20)) as a pure numerical fit to simulation data, which is not related to dynamical arguments. However, after correcting for the integral nature of the data presented in Power et al. (2003, equation (20)), a better approximation to the data is a simpler form of the radius contaning   particles, which we find in our simulations.

To summarize, we find that two-body scattering has less than 2% effect on the value of for haloes with as little as particles inside the radius defining . In order to achieve the same accuracy for smaller radii , the number of particles should be .

### 5.2 Subhaloes: tidal stripping and limits of survival

When a dark matter halo falls into another halo and becomes a subhalo, the tidal forces of the “parent” start to remove the least bound particles from the subhalo, thus reducing its mass. Effects of the tidal stripping experienced by an orbiting dark matter halo with a NFW profile have been extensively studied in the past (e.g., Klypin et al., 1999b; Hayashi et al., 2003; Peñarrubia et al., 2008; Arraki et al., 2012). Tidal stripping depends on the pericenter of the subhalo orbit, on the concentration of the subhalo, and on the number of orbits which the subhalo makes after falling into the parent halo. So, the situation is complex and it is difficult to find a realistic description of the whole problem without running complicated simulations. Nevertheless, we can derive simple relations, which can be used to understand the survival of subhaloes in our cosmological simulations. The main question, which we try to answer here, is how can a small halo with only particles survive when it falls into a larger halo with a strong tidal force?

To be more specific, we provide here simple estimates of the minimum distance from the center of the parent at which the subhalo, with a given number of particles, can be detected in cosmological simulations. Our analysis is simplified by the fact that the decrease in the mass of the central region of the halo with is very sensitive to the ratio between the tidal radius and the scale radius of the NFW density profile. For the central density and do not change much. For example, according to Figure 8 in Arraki et al. (2012) for the value of declines just by 20-30% as compared with the initial NFW value. This subhalo potentially can be detected by a halofinder. The situation is quite different for a smaller tidal radius. For example, for the maximum circular velocity decreases to 1/3-1/2 of its initial value. This means that the mass of the object becomes less than particles and subhaloes with this few particles cannot be detected.

The reason for this steep decline of bound mass is related to the fact that the central region is not self bound: kinetic energy of all particles inside this region is larger than the potential energy of these particles. Once the tidal radius is close to this unbound region, tidal stripping dramatically increases. It does not mean that the whole subhalo is totally disrupted because a small central region still may survive (Peñarrubia et al., 2008). However, in simulations with relatively small number of particles, the subhalo is lost. This is the reason why we assume that a subhalo is destroyed and cannot be detected in simulations once the tidal radius becomes less than . Our goal is to find the distance to center of the parent halo at which this happens.

There are two conditions which are used to find the tidal radius (e.g., Klypin et al., 1999b). The first is found by equating the external tidal force to the force due to the interior of the subhalo. The second is related with the resonance between the external force and the internal orbital motion of particles. We use the resonant condition for orbital frequencies, which gives slightly smaller tidal radius in central halo regions. If is the frequency of a particle on a circular orbit with radius from the center of the subhaloes and is the frequency of the perturbing external force, then the tidal radius is defined by condition .

The ratio of orbital frequencies can be written in the following form:

 ω2(r)Ω2(R)=(VmaxVp,max)2(Rsrs)2(yx)3f(x)f(y),y≡RRs, (18)

where and are the maximum circular velocity and the characteristic radius of the parent halo. Equation (18) can be simplified if we use halo concentrations and for the parent and the subhalo correspondingly. In this case, the equation for the tidal radius takes a simple form:

 u(xtide)=u(cs)u(cp)u(y),u(x)=x3f(x). (19)

Note that function has a simple meaning. It is inverse of the mean density inside radius . Equation (19) can be solved numerically to find the tidal radius for given halo concentration pair and for the distance from the parent center . It also can be used to find the distance from the parent at which the subhalo is destroyed. To find we assume that the subhalo is destroyed when its tidal radius becomes too small. If is the radius, than can be found by solving the following equation

 u(ylim)=u(cp)u(cs)u(xd). (20)

The following approximations give 15% accuracy for the tidal radius and for the destruction distance:

 rtide ≈ R1+0.04(R/Rs)3/2(rsRs)(cscp)3/2, (21) Rlim ≈ Rs0.9xd1+0.04x3/2d(cpcs)3/2,xd=0.1−3. (22)

One interesting consequence of equations (21-22) is that in a typical large-scale cosmological simulation subhaloes should be destroyed if their distance to the parent becomes less than the scale radius of the parent. Indeed, unless an exuberant number of particles is used for a subhalo, the subhalo will be destroyed once the tidal radius becomes less than . Denser subhaloes survive better against the tidal field: the tidal radius . However, the average concentration depends very weakly on mass . For subhaloes which are 100-1000 times less massive than the parent halo, . Putting these estimates into equation (22), we find that a combination of tidal stripping and small number of particles inside subhalos make survival of subhalos difficult for distances from the center of parent halo smaller than .

## 6 Conclusions

Large cosmological -body simulations, which use billions of particles and provide properties of millions of dark matter haloes, play a very important role in testing predictions of the standard cosmological model. It takes significant amount of computer resources and manpower to run and analyze these simulations. In this respect, it is crucial to optimize numerical parameters used to make the simulations and to understand their limitations. There is no unique set of requirements for the simulations: it all depends how the simulations will be used for particular science applications. For numerous purposes resolving the interior structure of haloes and identifying subhaloes are needed. Using state-of-the-art cosmological simulations and simplified models of individual dark matter haloes, we investigate the numerical accuracy and convergence of (sub)halo properties used for the Halo Abundance Matching method.

In the large volume cosmological simulations with billions of particles, which use the -body codes ART and Gadget, we find that convergence at the 10% accuracy for the abundance of haloes and subhaloes, and the correlation functions can be achieved with particles when parameters of the simulations – the force resolution and time-stepping – are chosen to satisfy a number of conditions.

Using simplified models of dark matter haloes, we study different effects, which may play a key role in the accuracy of results. Force resolution is the key parameter. Figure 8 shows how the circular velocity profile is affected by the lack of force resolution. Equation (9) can be used to estimate the error of for a given force resolution . Because the error in produces much larger error in the velocity function (), we recommend that the force resolution (equivalent of the Plummer softening) should be

 0.1rs<ϵ<0.3rs, (23)

where is the scale radius of smallest resolved haloes. Here the upper limit is defined by condition that the error in is less than 5%. At the lower limit the error is 1%. The lower limit is set to avoid too fine resolution, which may enhance the two-body scattering and may require too small time-step.

These constraints on the force resolution should be contrasted with the “optimal gravitational softening” recommended by Power et al. (2003): , where is the halo concentration and is the number of particles inside virial radius. Assuming and , we get . This force softening would result in errors in so large that it would render the simulation useless. Note that our recommended upper limit on is ten times smaller than the lower limit in Power et al. (2003).

We also find that the two-body scattering, though clearly detected, plays relatively minor role for small haloes with as little as particles. Estimates based on standard analytical approximations such as equation (12) and equation (16) are still useful because they provide insights on scaling of the relaxation time with distance and number of particles. However, analytical approximations overestimate the effect of the scattering. For the region which defines the value of the maximum circular velocity the analytical estimates give relaxation times, which are 5–6 times too short.

Two effects seem to contribute to the reduced impact of the two-body scattering in cosmological haloes: small number of crossing times for particles and reduced force resolution.

Time-stepping plays an important role in defining the accuracy of simulations (e.g., Power et al., 2003; Klypin et al., 2009). However, because of different implementations of time-stepping schemes and conditions for time-step refinement, it is difficult to give unique recommendations. For our Gadget simulations we use parameter ErrTolIntAccuracy, which is twice smaller than what was used for the Millennium-I simulation. The time-step was even smaller in the case of MultiDark, which used steps at the highest level of resolution.

Tidal stripping and numerical destruction of subhaloes also set limitations on the HAM models. As equation (22) shows, in large cosmological simulations subhaloes are destroyed once they come too close to the parent’s center. For realistic combinations of halo concentrations, this subhalo destruction occurs at the distance from the parent halo . Subhaloes can survive at smaller distances, if they have a very large number of particles. They also can be found there temporarily because it takes few orbits to get severely stripped, which can take few billion years.

We recommend the following steps to estimate parameters of -body simulations, which can be used to resolve central halo regions that define (sub)halo maximum circular velocity:

1. Using particle mass find virial mass containing particles. Find the halo virial radius. Use a concentration - mass relation to find halo concentration for smallest resolved halos.

2. Using the halo concentration and the virial radius, find scale radius for the smallest halo. Plummer force softening is defined by required accuracy of the maximum circular velocity . Use equation (9) to find , which gives this accuracy.

3. Run small tests with these parameters to find the time-step (or parameters, which define it) that give converge of halo profiles and (sub)halo abundances. The smaller is the force softening, the smaller is the time step.

## Acknowledgments

AK acknowledges the support of NSF and NASA grants to NMSU. GY acknowledges support from the Spanish MINECO under research grants AYA2012-31101, FPA2012-34694 and Consolider Ingenio SyeC CSD2007-0050 and from Comunidad de Madrid under ASTROMADRID project (S2009/ESP-1496). The BigMultiDark simulation suite have been performed in the Supermuc supercomputer at LRZ thanks to the cpu time awarded by PRACE (proposal number 2012060963). S.H. acknowledges support by the Deutsche Forschungsgemeinschaft under the grant .

## References

• Arraki et al. (2012) Arraki, K. S., Klypin, A., More, S., & Trujillo-Gomez, S. 2012, ArXiv e-prints
• Behroozi et al. (2010) Behroozi, P. S., Conroy, C., & Wechsler, R. H. 2010, ApJ, 717, 379
• Behroozi et al. (2013a) Behroozi, P. S., Wechsler, R. H., Lu, Y., Hahn, O., Busha, M. T., Klypin, A., & Primack, J. R. 2013a, ArXiv e-prints
• Behroozi et al. (2013b) Behroozi, P. S., Wechsler, R. H., Wu, H.-Y., Busha, M. T., Klypin, A. A., & Primack, J. R. 2013b, ApJ, 763, 18
• Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition (Princeton University Press)
• Boylan-Kolchin et al. (2009) Boylan-Kolchin, M., Springel, V., White, S. D. M., Jenkins, A., & Lemson, G. 2009, MNRAS, 398, 1150
• Conroy et al. (2006) Conroy, C., Wechsler, R. H., & Kravtsov, A. V. 2006, ApJ, 647, 201
• Diemand et al. (2004) Diemand, J., Moore, B., Stadel, J., & Kazantzidis, S. 2004, MNRAS, 348, 977
• Gottlöber & Klypin (2008) Gottlöber, S., & Klypin, A. 2008, ArXiv e-prints
• Guo & White (2013) Guo, Q., & White, S. 2013, ArXiv e-prints
• Guo et al. (2011) Guo, Q., White, S., Boylan-Kolchin, M., De Lucia, G., Kauffmann, G., Lemson, G., Li, C., Springel, V., & Weinmann, S. 2011, MNRAS, 413, 101
• Guo et al. (2010) Guo, Q., White, S., Li, C., & Boylan-Kolchin, M. 2010, MNRAS, 404, 1111
• Hayashi et al. (2003) Hayashi, E., Navarro, J. F., Taylor, J. E., Stadel, J., & Quinn, T. 2003, ApJ, 584, 541
• Hearin et al. (2012) Hearin, A. P., Zentner, A. R., Berlind, A. A., & Newman, J. A. 2012, ArXiv e-prints
• Hess et al. (2013) Hess, S., Yepes, G., Prada, F., Klypin, A., & Gottlöber, S. 2013, ArXiv e-prints
• Klypin et al. (1999a) Klypin, A., Gottlöber, S., Kravtsov, A. V., & Khokhlov, A. M. 1999a, ApJ, 516, 530
• Klypin et al. (1999b) —. 1999b, ApJ, 516, 530
• Klypin & Holtzman (1997) Klypin, A., & Holtzman, J. 1997, ArXiv Astrophysics e-prints
• Klypin et al. (2001) Klypin, A., Kravtsov, A. V., Bullock, J. S., & Primack, J. R. 2001, ApJ, 554, 903
• Klypin et al. (2009) Klypin, A., Valenzuela, O., Colín, P., & Quinn, T. 2009, MNRAS, 398, 1027
• Klypin et al. (2011) Klypin, A. A., Trujillo-Gomez, S., & Primack, J. 2011, ApJ, 740, 102
• Knebe et al. (2000) Knebe, A., Kravtsov, A. V., Gottlöber, S., & Klypin, A. A. 2000, MNRAS, 317, 630
• Knebe et al. (2013) Knebe, A., Pearce, F. R., Lux, H., Ascasibar, Y., Behroozi, P., Casado, J., Corbett Moran, C., Diemand, J., Dolag, K., Dominguez-Tenreiro, R., Elahi, P., Falck, B., Gottloeber, S., Han, J., Klypin, A., Lukic, Z., Maciejewski, M., McBride, C. K., Merchan, M. E., Muldrew, S. I., Neyrinck, M., Onions, J., Planelles, S., Potter, D., Quilis, V., Rasera, Y., Ricker, P. M., Roy, F., Ruiz, A. N., Sgro, M. A., Springel, V., Stadel, J., Sutter, P. M., Tweed, D., & Zemp, M. 2013, ArXiv e-prints
• Knebe et al. (2011) Knebe, A., et al. 2011, MNRAS, 415, 2293
• Kravtsov (2013) Kravtsov, A. V. 2013, ApJL, 764, L31
• Kravtsov et al. (2004a) Kravtsov, A. V., Berlind, A. A., Wechsler, R. H., Klypin, A. A., Gottlöber, S., Allgood, B., & Primack, J. R. 2004a, ApJ, 609, 35
• Kravtsov et al. (2004b) Kravtsov, A. V., Gnedin, O. Y., & Klypin, A. A. 2004b, ApJ, 609, 482
• Kravtsov et al. (1997) Kravtsov, A. V., Klypin, A. A., & Khokhlov, A. M. 1997, ApJS, 111, 73
• Leauthaud et al. (2011a) Leauthaud, A., Tinker, J., Behroozi, P. S., Busha, M. T., & Wechsler, R. H. 2011a, ApJ, 738, 45
• Leauthaud et al. (2011b) —. 2011b, ApJ, 738, 45
• Li & White (2009) Li, C., & White, S. D. M. 2009, MNRAS, 398, 2177
• Mandelbaum et al. (2006) Mandelbaum, R., Seljak, U., Kauffmann, G., Hirata, C. M., & Brinkmann, J. 2006, MNRAS, 368, 715
• Nuza et al. (2013) Nuza, S. E., Sánchez, A. G., Prada, F., Klypin, A., Schlegel, D. J., Gottlöber, S., Montero-Dorta, A. D., Manera, M., McBride, C. K., Ross, A. J., Angulo, R., Blanton, M., Bolton, A., Favole, G., Samushia, L., Montesano, F., Percival, W. J., Padmanabhan, N., Steinmetz, M., Tinker, J., Skibba, R., Schneider, D. P., Guo, H., Zehavi, I., Zheng, Z., Bizyaev, D., Malanushenko, O., Malanushenko, V., Oravetz, A. E., Oravetz, D. J., & Shelden, A. C. 2013, MNRAS
• Peñarrubia et al. (2008) Peñarrubia, J., Navarro, J. F., & McConnachie, A. W. 2008, ApJ, 673, 226
• Power et al. (2003) Power, C., Navarro, J. F., Jenkins, A., Frenk, C. S., White, S. D. M., Springel, V., Stadel, J., & Quinn, T. 2003, MNRAS, 338, 14
• Prada et al. (2012) Prada, F., Klypin, A. A., Cuesta, A. J., Betancort-Rijo, J. E., & Primack, J. 2012, MNRAS, 423, 3018
• Quinlan (1996) Quinlan, G. D. 1996, NewA, 1, 255
• Reddick et al. (2012) Reddick, R. M., Wechsler, R. H., Tinker, J. L., & Behroozi, P. S. 2012, ArXiv e-prints
• Riebe et al. (2011) Riebe, K., Partl, A. M., Enke, H., Forero-Romero, J., Gottloeber, S., Klypin, A., Lemson, G., Prada, F., et al. 2011, ArXiv e-prints
• Springel (2005) Springel, V. 2005, MNRAS, 364, 1105
• Springel et al. (2008) Springel, V., Wang, J., Vogelsberger, M., Ludlow, A., Jenkins, A., Helmi, A., Navarro, J. F., Frenk, C. S., & White, S. D. M. 2008, MNRAS, 391, 1685
• Springel et al. (2005) Springel, V., White, S. D. M., Jenkins, A., Frenk, C. S., Yoshida, N., Gao, L., Navarro, J., Thacker, R., Croton, D., Helly, J., Peacock, J. A., Cole, S., Thomas, P., Couchman, H., Evrard, A., Colberg, J., & Pearce, F. 2005, Nature, 435, 629
• Stadel et al. (2009) Stadel, J., Potter, D., Moore, B., Diemand, J., Madau, P., Zemp, M., Kuhlen, M., & Quilis, V. 2009, MNRAS, 398, L21
• Tasitsiomi et al. (2004) Tasitsiomi, A., Kravtsov, A. V., Wechsler, R. H., & Primack, J. R. 2004, ApJ, 614, 533
• Trujillo-Gomez et al. (2011) Trujillo-Gomez, S., Klypin, A., Primack, J., & Romanowsky, A. J. 2011, ApJ, 742, 16
• Vale & Ostriker (2004) Vale, A., & Ostriker, J. P. 2004, MNRAS, 353, 189
• Valenzuela & Klypin (2003) Valenzuela, O., & Klypin, A. 2003, MNRAS, 345, 406
{"extraction_info": {"found_math": false, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 0, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.9437059760093689, "perplexity": 1747.2308651706273}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2020-45/segments/1603107922746.99/warc/CC-MAIN-20201101001251-20201101031251-00713.warc.gz"}