# The Adaptive TreePM: An Adaptive Resolution Code for Cosmological N-body Simulations

J. S. Bagla and Nishikanta Khandai
Harish-Chandra Research Institute, Chhatnag Road, Jhunsi,
E-mail: jasjeet@hri.res.in, nishi@hri.res.in
###### Abstract

Cosmological -Body simulations are used for a variety of applications. Indeed progress in the study of large scale structures and galaxy formation would have been very limited without this tool. For nearly twenty years the limitations imposed by computing power forced simulators to ignore some of the basic requirements for modeling gravitational instability. One of the limitations of most cosmological codes has been the use of a force softening length that is much smaller than the typical inter-particle separation. This leads to departures from collisionless evolution that is desired in these simulations. We propose a particle based method with an adaptive resolution where the force softening length is reduced in high density regions while ensuring that it remains well above the local inter-particle separation. The method, called the Adaptive TreePM, is based on the TreePM code. We present the mathematical model and an implementation of this code, and demonstrate that the results converge over a range of options for parameters introduced in generalizing the code from the TreePM code. We explicitly demonstrate collisionless evolution in collapse of an oblique plane wave. We compare the code with the fixed resolution TreePM code and also an implementation that mimics adaptive mesh refinement methods and comment on the agreement, and disagreements in the results. We find that in most respects the ATreePM code performs at least as well as the fixed resolution TreePM in highly over-dense regions, from clustering and number density of haloes, to internal dynamics of haloes. We also show that the adaptive code is faster than the corresponding high resolution TreePM code.

###### keywords:
gravitation, methods: N-Body simulations, cosmology: large scale structure of the universe
pubyear: 2008pubyear: 2008pagerange: LABEL:firstpageA

## 1 Introduction

Large scale structures traced by galaxies are believed to have formed by amplification of small perturbations (Peebles, 1980; Peacock, 1999; Padmanabhan, 2002; Bernardeau et al., 2002). Galaxies are highly over-dense systems, matter density in galaxies is thousands of times larger than the average density in the universe. Typical density contrast () in matter at these scales in the early universe was much smaller than unity. Thus the problem of galaxy formation and the large scale distribution of galaxies requires an understanding of the evolution of density perturbations from small initial values to the large values we encounter today.

Initial density perturbations were present at all scales that have been observed (Spergel et al., 2007; Percival et al., 2007). The equations that describe the evolution of density perturbations in an expanding universe have been known for several decades (Peebles, 1974) and these are easy to solve when the amplitude of perturbations is small. Once density contrast at relevant scales becomes comparable to unity, perturbations becomes non-linear and coupling with perturbations at other scales cannot be ignored. The equation for evolution of density perturbations cannot be solved for generic initial conditions in this regime. N-Body simulations (e.g., see Efstathiou et al. (1985); Bertschinger (1998); Bagla & Padmanabhan (1997); Bagla (2005); Dolag et al. (2008)) are often used to study the evolution in this regime. Alternative approaches can be used if one requires only a limited amount of information and in such a case either quasi-linear approximation schemes (Bernardeau et al., 2002; Zel’Dovich, 1970; Gurbatov, Saichev, & Shandarin, 1989; Matarrese et al., 1992; Brainerd, Scherrer, & Villumsen, 1993; Bagla & Padmanabhan, 1994; Sahni & Coles, 1995; Hui & Bertschinger, 1996) or scaling relations (Davis & Peebles, 1977; Hamilton et al., 1991; Jain, Mo, & White, 1995; Kanekar, 2000; Ma, 1998; Nityananda & Padmanabhan, 1994; Padmanabhan et al., 1996; Peacock & Dodds, 1994; Padmanabhan, 1996; Peacock & Dodds, 1996; Smith et al., 2003) suffice. However, even the approximation schemes and scaling relations must be compared and calibrated with simulations before these can be used with confidence.

Last three decades have seen a rapid development of techniques and computing power for cosmological simulations and the results of these simulations have provided valuable insight into the study of structure formation. The state of the art simulations used less than particles two decades ago (Efstathiou et al., 1988) and if the improvement had been due only to computing power then the largest simulation possible today should have been around particles, whereas the largest simulations done till date used more than particles (Springel et al., 2005). Evidently, development of new methods and optimizations has also played a significant role in the evolution of simulation studies (Efstathiou et al., 1985; Barnes & Hut, 1986; Greengard & Rokhlin, 1987; Bouchet & Hernquist, 1988; Jernigan & Porter, 1989; Hernquist, 1990; Makino, 1990, 1991; Hernquist, Bouchet, & Suto, 1991; Couchman, 1991; Ebisuzaki et al., 1993; Theuns, 1994; Brieu, Summers, & Ostriker, 1995; Suisalu & Saar, 1995; Xu, 1995; Dubinski, 1996; Kravtsov, Klypin, & Khokhlov, 1997; Macfarland et al., 1998; Bode, Ostriker, & Xu, 2000; Brieu & Evrard, 2000; Dehnen, 2000; Knebe, Green, & Binney, 2001; Springel, Yoshida, & White, 2001; Kawai & Makino, 2001; Makino, 2002; Dehnen, 2002; Bagla, 2002; Bagla & Ray, 2003; Makino et al., 2003; Bode & Ostriker, 2003; Ray & Bagla, 2004; Dubinski et al., 2004; Makino, 2004; Springel, 2005; Merz, Pen, & Trac, 2005; Yoshikawa & Fukushige, 2005; Wadsley, Stadel, & Quinn, 2004; Thacker & Couchman, 2006; Khandai & Bagla, 2008) Along the way, code developers have also successfully met the challenge posed by the emergence of distributed parallel programming.

In modeling gravitational clustering, cosmological N-Body codes should ensure the following:

• The universe does not have a boundary. Therefore cosmological simulations need to be run with periodic boundary conditions111An exception are simulations of large spherical volumes that do not suffer significant deformation during the course of evolution.. The simulation volume should be large enough for the effects of missing modes to be small, (Bagla & Ray, 2005; Bagla & Prasad, 2006; Power & Knebe, 2006; Prasad, 2007; Bagla, Prasad, & Khandai, 2008).

• The mass of each particle in simulations should be much smaller than mass scales of interest in the simulation output. This is to ensure adequate mass resolution.

• Each particle in an N-Body simulation represents a very large number of particles/objects in the universe. Thus we must ensure that pair-wise interaction of N-Body particles is softened at scales comparable with the local inter-particle separation. If this is not ensured then the resulting two body collisions introduce errors in the resulting distribution of particles, (Splinter et al., 1998; Binney & Knebe, 2002).

In spite of the vast improvement in computing power, simulators have often had to compromise on one or more of these points. Often, errors also creep in due to the approximate methods used for solving for force in simulations. A large fraction of current cosmological N-Body codes suffer from collisionality or force biasing. Force is biased when softening lengths are much larger than the local mean inter-particle separation, . Whereas a complementary effect, collisions, occur whenever . The reader is referred to Dehnen (2001) for a detailed discussion on these two effects. Codes that adapt their softening lengths to local densities generally are mostly of the adaptive mesh refinement type. These use a grid for solving the Poisson equation and often have anisotropies in force at the scales comparable to the size of a grid cell. Codes which use fixed softening lengths are not entirely collisionless, and, in highly over-dense regions biasing of force also exists. The TPM, (Xu, 1995; Bode, Ostriker, & Xu, 2000; Bode & Ostriker, 2003) is a particle based code with a one step adaptive resolution. However, in this case the use of the unmodified PM approach for computing the long range force introduces significant errors at scales comparable to the grid. In this work we describe a code which addresses all three issues of force anisotropy, collisionality and force biasing by employing an adaptive softening formalism with the hybrid TreePM code.

The choice of the optimal softening length has been discussed at length (Merritt, 1996; Athanassoula et al., 2000; Dehnen, 2001; Price & Monaghan, 2007). These studies were carried out in the context of isolated haloes in dynamical equilibrium, e.g., Plummer and Hernquist profiles, therefore errors could be defined clearly. It is also possible to compute and compare various physical quantities with analytical expressions derived from the distribution function. Dehnen (2001) derived analytical expressions for errors in the context of these profiles. This work suggested that the optimal softening length must adapt to the local inter-particle separation as a function of space and time. Price & Monaghan (2007) have developed an energy-momentum conserving formalism with adaptive softening and demonstrated that it was superior to fixed softening.

The evolution of perturbations at small scales depends strongly on the mass and force resolution. High force resolutions can lead to better modeling of dense haloes, but gives rise to two body collisions in regions where the softening length is smaller than the local inter particle separation (Splinter et al., 1998; Binney & Knebe, 2002; Diemand et al., 2004; Binney, 2004; El-Zant, 2006; Romeo et al., 2008). As all particles in very high density regions go through such a phase during evolution, any errors arising due to two body collisions can potentially effect the structure of high density haloes that form. A high force resolution without a corresponding mass resolution can also give misleading results as we cannot probe shapes of collapsed objects Kuhlman, Melott, & Shandarin (1996). In addition, discreteness and stochasticity also limit our ability to measure physical quantities in simulations, and these too need to be understood properly (Thiébaut et al., 2008; Romeo et al., 2008). In all such cases, the errors in modeling is large at small scales. It is important to understand how such errors may spread to larger scales and affect physical quantities.

This work is organized as follows. In §2 we describe the formalism for adaptive softening in a cosmological -body code. In §3 we review the TreePM method, which forms the backbone of our gravity solver and describe how adaptive softening is implemented with it. We briefly discuss the performance characteristics of ATreePM in §4. We discuss validation of the ATreePM code in §5 and we conclude in §6.

### 2.1 Formalism

The aim of any collisionless -body code is to self consistently evolve the phase-space distribution function (DF) under its own gravitational force field :

 ∂tf + (∂rf).v+(∂vf).F=0 (1) F(r,t) = −G∫∫d3r′d3vr−r′|r−r′|3f(r′,v,t) (2)

The approach that one takes is to sample the DF, by  phase-space points, , at initial time . Liouville’s theorem then states that evolving the trajectories of these points to any time will be a representation of the DF at that time. Since the system is a collisionless one, one has to suppress artificial two-body collisions arising out of interactions between particles which were generated for sampling the density field. One therefore assigns a finite size to N-Body particles which ensures softening of force at small scales, instead of assuming these to be point particles. The density field when sampled by point particles,

 ρ(r)=∫d3vf(r,v,t)≡N∑j=1mjδ3D(r−rj) (3)

is now smoothed at small scales if we assign a finite size to every particle:

 ρ(r)=∫d3vf(r,v,t)≡N∑j=1mjW(|r−rj|,ϵ) (4)

Where is known as the smoothing kernel and we have assumed that particles are spherical in shape. Here is the mass of the particle. We can now integrate the Poisson equation to obtain the expression for the kernel for computing force and potential. Both the quantities are softened at scales below the softening length . We choose to work with the cubic spline kernel (Monaghan & Lattanzio, 1985) whose expression is given below. Complete expressions for the potential and force are given in the Appendix (See Eqn.(A, 28)).

 W(u,ϵ)=8πϵ3⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩1−6(uϵ)2+6(uϵ)3,0≤uϵ<0.52(1−uϵ)3,0.5≤uϵ<1.00,1.0≤uϵ (5)

In the context of individual softening lengths the density at the location of the particle is given by:

 ρ(ri)=N∑j=1mjW(rij,ϵi)=nn∑j=1mjW(rij,ϵi) (6)

where i,j indicate the particle indices, and . The summation in principle can be extended upto infinity if the kernel is infinite in extent (e.g. plummer, gaussian kernels). But since such kernels tend to bias the force (Dehnen, 2001; Price & Monaghan, 2007) we will work only with those kernels which have compact support, in particular the cubic spline kernel. Such kernels ensure that the force is Newtonian beyond the softening scale. is the number of the nearest nearbors within for the particle . In the discussion that follows we assume that this number is fixed for every particle and sets the value of the softening length . We implicitly assume it in our summation. Integrating the Poisson equation we obtain the Green’s function for the potential , where the functional form for is given in the Appendix (See Eqn.(A, 28)).

With the introduction of individual softening lengths for particles, the symmetry of the potential is lost and momentum conservation is violated. Since is now a local quantity, the EOM is incomplete if one takes the expression of force with the fixed softening length replacing it by a local softening length . Hence energy conservation also gets violated. This is because the force is derived from the potential and with the introduction of a local softening length in the potential, the gradient must also act on giving us an extra term. Traditionally this grad- () term has been ignored since in typical applications these were found to be subdominant when compared to the usual force (Gingold & Monaghan, 1982; Evrard, 1988; Hernquist & Barnes, 1990; Monaghan, 1992). It has been shown recently (Price & Monaghan, 2007), that this term plays an important role in N-Body simulations. We study the impact of ignoring this term in Cosmological simulations.

A remedy for momentum non-conservation is to use a symmetrized softening length

 ϵij=12(ϵi+ϵj) (7)

and plug it into the expression for density to re-derive a symmetrized expression for the potential as . This prescription changes the softening length and hence the neighborlist, which one has to recompute. Another disadvantage is that with this prescription for symmetrization, the number of neighbors is not fixed for every particle and hence errors in all smoothed estimates are not the same for every particle. An alternate method (Hernquist & Katz, 1989) is to symmetrize the kernel itself.

 ˜Wij = ˜Wji=12[W(rij,ϵj)+W(rij,ϵi)] ˜ϕij = 12[ϕ(rij,ϵj)+ϕ(rij,ϵi)] (8)

The total potential is thus

 Φtot=12N∑i,j˜ϕij (9)

We can use this to write a Lagrangian which is manifestly symmetric and this ensures momentum conservation. Energy is conserved only if the term is retained in the EOM.

 L=N∑i12miv2i−G2N∑i,jmimj˜ϕij (10)

The EOM of motion can be derived with this Lagrangian (the reader is referred to Price & Monaghan (2007) for details)

 dvidt= − G∑jmj˜ϕ′ijri−rj|ri−rj| (11) − G2∑jmj[ζiΩi∂Wij(ϵi)∂ri+ζjΩj∂Wij(ϵj)∂ri]

Where the first term is the standard Newtonian force term (we refer to it as the term). The second is the energy conserving term which would be zero for fixed softening. Notice that all terms are antisymmetric in and hence the total momentum is conserved. Here and are:

 ζi≡∂ϵi∂ρi∑jmj∂ϕ(rij,ϵi)∂ϵi (12) Ωi=[1−∂ϵi∂ρi∑jmj∂Wij(ϵi)∂ϵi] (13)

Expressions for , and are given in the Appendix (See 29,30 and31). As we have . The term term ensures that the EOM is accurate to all orders in (Springel & Hernquist, 2002).

## 3 The Adaptive TreePM Method

### 3.1 The TreePM Method

The TreePM algorithm (Bagla, 2002; Bagla & Ray, 2003) is a hybrid N-Body method which improves the accuracy and performance of the Barnes-Hut (BH) Tree method (Barnes & Hut, 1986) by combining it with the PM method (Bagla & Padmanabhan, 1997; Merz, Pen, & Trac, 2005; Klypin & Shandarin, 1983; Miller, 1983; Bouchet & Kandrup, 1985; Bouchet, Adam, & Pellat, 1985; Hockney & Eastwood, 1988). The TreePM method explicitly breaks the potential into a short-range and a longe-range component at a scale . The PM method is used to calculate long-range force and the short-range force is computed using the BH Tree method. Use of the BH Tree for short-range force calculation enhances the force resolution as compared to the PM method.

The gravitational force is divided into a long range and a short range part using partitioning of unity in the Poisson equation.

 φk = −4πGρkk2 φk = φlrk+φsrk φlrk = −4πGρkk2exp(−k2r2s) (14) φsrk = −4πGρkk2[1−exp(−k2r2s)] (15)

Here and are the short-range and long-range potentials in Fourier space. is the density, G is the gravitational coupling constant and is the scale at which the splitting of the potential is done. It has been shown that this particular split between the long range and the short range force is optimal amongst a large class of suitable functional forms (Bagla & Ray, 2003). The long-range force is computed in Fourier space with the PM method and the short-range force is computed in real space with the Tree method. The short range potential and force in real space are:

 φsr(r,ϵ) = Gmϕ(r,ϵ)erfc(r2rs) (16) Fsr(r,ϵ) = Gmf(r,ϵ)C(rrs) (17) C(rrs) = [erfc(r2rs)+rrs√πexp(−r24r2s)] (18)

Here erfc is the complementary error function. and are the usual potential and force kernels, respectively. modifies the softened Newtonian force kernel to the short-range force . The expression for and depend on the kernel used for smoothing and are given for cubic spline in the Appendix. We find that tabulating and and using interpolation to compute these functions is much more effective than calculating them every time. The short range force is below of the total force at . The short range force is therefore computed within a sphere of radius using the BH tree method.

The long range force falls below of the total force for . Thus the choice of softening length does not affect the long range force in a significant manner as long as the force softening is done with a kernel that has a compact support and the softening length is below .

The BH tree structure is built out of cells and particles. Cells may contain smaller cells (sub-cells) within them. Sub-cells can have even smaller cells within them, or they can contain a particle. In three dimensions, each cubic cell is divided into eight cubic sub-cells. Cells, as structures, have attributes like total mass, location of center of mass and pointers to sub-cells. Particles, on the other hand have the usual attributes: position, velocity and mass.

Force on a particle is computed by adding contribution of other particles or of cells. A cell that is sufficiently far away can be considered as a single entity and we can add the force due to the total mass contained in the cell from its center of mass. If the cell is not sufficiently far away then we must consider its constituents, sub-cells and particles. Whether a cell can be accepted as a single entity for force calculation is decided by the cell acceptance criterion (CAC). We compute the ratio of the size of the cell and the distance from the particle in question to its center of mass and compare it with a threshold value

 θ=Lcellr≤θc (19)

The error in force increases with . Poor choice of can lead to significant errors (Salmon & Warren, 1994). Many different approaches have been tried for the CAC in order to minimize error as well as CPU time usage (Salmon & Warren, 1994; Springel, Yoshida, & White, 2001). The tree code gains over direct summation as the number of contributions to the force is much smaller than the number of particles.

The TreePM method is characterized therefore by three parameters, , and . For a discussion on the optimum choice of these parameters the reader is referred to Bagla & Ray (2003). For all our tests we choose conservative values , and which give errors below in force. All lengths are specified in units of the PM grid.

### 3.2 The Adaptive TreePM Method

We choose to implement adaptive softening with a modified TreePM code (Khandai & Bagla, 2008), which incorporates Barnes optimization using groups (Barnes, 1990; Makino, 1991; Yoshikawa & Fukushige, 2005), into the TreePM code. In principle one can also incorporate a similar formalism for treecodes (Springel, Yoshida, & White, 2001), PM codes (Couchman, 1991; Couchman, Thomas, & Pearce, 1995) and other variants like TPM (Xu, 1995), TreePM (Bagla, 2002; Springel, 2005) and GOTPM (Dubinski et al., 2004). We shall discuss one advantage of using the TreePM code with the group optimizations below.

#### 3.2.1 Estimating the Softening Length

Our first task is to get an estimate of the softening length for each particle. A natural way to extract a local length scale uses the numerical value of local density. The local number density is related to the softening length as:

 ϵ(r)∼(nnn(r)))13 (20)

is the number density at the location of the particle and we assume all particles have the same mass. Here, is a reference number and we take it to be the number of neighbors used for estimation of the number density. The above equation is implicit and can be solved iteratively, see, e.g., Springel & Hernquist (2002); Springel (2005); Price & Monaghan (2007) for details. Price & Monaghan (2007) have shown that errors are not very sensitive to the exact value of . We choose in our simulations, and also comment on variation in results with this choice.

We are using the formalism developed by Price & Monaghan (2007) for achieving an adaptive resolution in gravitational interactions of particles. As the formalism was developed in the context of SPH codes, and some of the quantities required can be computed naturally using the SPH method, the same was used even though the gravitational interaction is completely collisionless. We follow a similar implementation and use methods commonly used in SPH simulations, even though there are no hydrodynamical effects present in the gravitational interaction being studied here. For an overview of SPH methods, please see Lucy (1977); Gingold &Monaghan (1977). The SPH methods assign values for functions like the density to particles by averaging over nearest neighbors. The list of nearest neighbors is an essential requirement for computing anything using these methods. All quantities (, , , etc.) can be computed at runtime with this neighborlist once we have converged to a value for by solving Eqn.(20). We compute the neighborlist using linked lists (Hernquist & Katz, 1989).

We put bounds on the maximum and minimum softening lengths, and . A maximum bound is required so that force softening is restricted to the short-range force only, we choose in order to ensure that the long range force is of order of the total force (or smaller) at scales where force softening is important. This ensures that any errors arising from non-modification of the long range force are smaller than , if one puts a lower threshold on the maximum allowed error then the scale has to be lowered correspondingly. Alternatively, one can work with a larger and then it becomes possible to allow a larger . A lower bound is also required for the reason that we do not want a few isolated highly over-dense regions to dominate the CPU time requirements. The value of the lower bound must correspond to densities that are much higher than highest over-densities of interest in the simulation. We set this lower bound , which corresponds to . The lower bound however is not critical to the structure of the code and may be omitted.

In our implementation a neighbor search is carried out only upto . Particles which do not have neighbors within are assigned and the spherical top hat (STH) density is assigned with the number of neighbors within . For these particles we assign and , which makes their . A similar assignment is carried out for particles which have . The term is calculated before the short-range force so that the individual softening lengths needed for short-range force calculation are also assigned in the process.

#### 3.2.2 Memory Requirements

Even though we use two separate data structures, namely linked lists for and tree for in order to compute the total short-range force, additional memory requirements compared to TreePM are minimal: we require one additional array for storing the softening lengths. The term does not require an additional array since it is a component of the short-range force and it can be computed at run time. This is because the two data structures are never required at the same time. We require specification of the largest force softening length in a given cell (See the subsection on the cell acceptance criterion below.). This amounts to a single precision array of the same size as the number of cells in the tree.

An advantage of using an analytical splitting of force, in the manner TreePM does, is that computation of short-range force does not need global data structures. For example one can geometrically divide regions into smaller regions and construct local trees and linked lists in them (just like one would go about doing it in a distributed code) and iterate through these regions for computing short-range force instead of constructing one global data structure for the entire simulation volume for computing the short-range force (Dubinski et al., 2004). This reduces memory usage significantly and the dominant part is taken up by the arrays required for computing the long range PM force.

#### 3.2.3 Timestepping Criterion and Integration

We have implemented a hierarchical time integrator similar to that used in GADGET-2 Springel (2005), in which particle trajectories are integrated with individual timesteps and synchronized with the largest timestep. As we allow the block time step222Same as the largest time step. to vary with time, we work with the so called KDK approach (Kick-Drift-Kick) in which velocities are updated in two half steps whereas position is updated in a full step. It can be shown that with a variable time step, KDK performs better than DKD (Drift-Kick-Drift) (see Springel (2005) for details.). We give separate PM (long range, global) kicks and Tree (short range, individual) kicks333Tree kick includes the contribution of the term.. The block timestep is determined by the particle which has the maximum PM acceleration :

 ΔtPM=δt(ϵmaxaPMmax)1/2 (21)

Here is the dimensionless accuracy parameter. In our implementation of the hierarchy of time steps, the smaller time steps differ by an integer power () of from the largest, block time step. An array is then used to store the value which determines the timestep of the particle. Individual timesteps are first calculated:

 Δtsri=δt(ϵiasri)1/2 (22)

and then the appropriate hierarchy is chosen depending on this value. Here and are the modulus of the individual short-range acceleration (sum of the and terms) and the softening lengths respectively. TreePM has a similar time-stepping criterion with replace by .

The code drifts all the particles with the smallest timestep to the next time, where a force computation is done for particles that require an updation of velocity (Kick). However the neighborlist and individual softening lengths are computed for all particles at every small timestep. This is because, even though some particles do not require a velocity update, their neighbors might require one for which they would contribute through their updated softening lengths. The term however can be computed only for those particles requiring a velocity update.

Within a given block time step, the smaller time steps are constant for a given particle. The time step changes across block time step and this brings in inaccuracies in evolution of trajectories. It is possible, in principle, to ensure that the second order accuracy is maintained here. However, we find that the time steps for particles change very slowly and this change does not affect trajectories in a significant manner.

The Courant condition is satisfied for the choice of we use. Indeed, we chose by requiring that two particles in a highly eccentric orbit around each other maintain the trajectory correctly for tens of orbits.

#### 3.2.4 Cell Acceptance Criterion For ATreePM

The adaptive force resolution formalism requires us to symmetrize force between particles that are separated by a distance smaller than the larger of the two softening lengths. Without this, the momentum conservation cannot be ensured. For pairs of particles separated by larger distances, there is no need to explicitly symmetrize force as there is no dependence on the softening length at these scales444As an aside we would like to note that the Tree method does not conserve momentum explicitly. This is because the tree traversal approximates the force due to pairwise interactions, and in the process the pairwise symmetry is lost. In the modified Tree method (Barnes, 1990; Makino, 1991; Yoshikawa & Fukushige, 2005; Khandai & Bagla, 2008) explicit pairwise force is computed for particles within each group. Particles within a group have a common interaction list for force due to particles outside the group. Exact pairwise PP force is computed explicitly for intra-group particles and we have a pairwise symmetry for this component, but for interaction with particles outside the group there is no explicit pairwise symmetry and hence no explicit momentum conservation.. Thus the cell acceptance criterion needs to be changed within the tree part of the code to ensure that for pairs of particles separated by the critical distance, forces are computed through pairwise particle-particle (PP) interaction. In our implementation, along with particles, cells (and hence groups) are assigned softening lengths corresponding to the particle contained within the cell that has the largest softening length. The CAC, that is modified when we go from the TreePM to the modified TreePM (Khandai & Bagla, 2008) has to be changed again to take into account the largest softening lengths for particles within cells and groups whose interaction is to be computed.

 (Lcell+ϵmaxcell)(rgcm−Lgroup−ϵmaxgroup)≤θc (23)

are the softening lengths of the cell and group (in the sense explained above) respectively. is the distance separating the centers of mass of the group and the cell. is the size of the cell and is the distance to the furthest particle from the center of mass of the group. This CAC ensures that the interaction of particles separated by less than the softening length is computed in a direct pairwise manner and hence can be explicitly symmetrized. These are also the pairs for which the term needs to be computed.

## 4 Performance characteristics

### 4.1 Timing

We now look at the wall clock time as a measure of performance between different codes. We studied evolution of a power law model with using particles up to the stage where the scale of non-linearity is grid lengths. More details of the run are given in the section on validation of the ATreePM code. Figure 1 shows wall clock time as a function of softening length for TreePM (left panel) and the wall clock time as a function of for the ATreePM codes.

We can qualitatively understand the slope for the TreePM curve by looking at the equivalent time-stepping criterion as Eqn.(22) for TreePM.

 Δtsri=δt(ϵasri)1/2 (24)

From here the naive expectation is that the time taken should scale as . The slope of the curve is in the range to . The reason for this small deviation lies in our use of a hierarchy of time steps, where trajectories of all the particles are not updated at every time step. However, positions of particles are updated at every time step and this operation as well as those related to creating the tree structure at every step add an overhead. This overhead becomes more and more important at small where we have many more levels of hierarchy realized in a simulation. This leads to steepening of the curve from the simple expectation given above.

For the ATreePM, we expect softening lengths to be larger for larger . On the other hand, a larger implies a larger neighborlist and the time taken for setting up the neighborlist increases. The second effect is the dominant one and we see that the time taken for Adaptive TreePM increases with . We see that time taken by both variants of ATreePM is similar for and . This indicates that the time taken for calculation of the term is negligible. There is a difference between the timing for , as the code with the does not evolve the system correctly: this can be seen in all the indicators like the amplitude of clustering, mass function, etc., presented in the next section. We see that TreePM with takes more time than ATreePM with , whereas the time taken are comparable with ATreePM with .

Since is assigned by hand in TreePM, the use of a small softening length means a considerable number of particles obtain a small timestep even though the acceleration for these particles is small and does not vary rapidly with time. In a hierarchical integrator that we use, this means that force computation is done more often within a block timestep. In ATreePM on the other hand a small softening is only assigned to particles in over-dense regions. This saves time in the under-dense and not so over-dense regions and the ATreePM code devotes more time to evolving trajectories in highly over-dense regions. This is illustrated in Figure 2 where we have plotted the cumulative distribution of softening lengths at the final epoch in one of the simulations used here (see lower-left panel of Figure 8). The softening length was computed here for . We find that only of the particles have a softening length smaller than the smallest softening length used for the fixed resolution simulations. Even at this epoch where highly non-linear clustering is seen, nearly half the particles have a softening length corresponding to the maximum value of . This shows how we are able to evolve the system in an ATreePM simulation with a lower computational cost while resolving highly over-dense regions.

### 4.2 Errors

In this section we discuss the dependence of errors on . In cosmological simulations it is difficult to define errors when softening lengths are varied due to the lack of a reference setup. Nevertheless we can choose the optimal softening length such that one minimizes the globally averaged fluctuation in force as the softening lengths are varied. Following Price & Monaghan (2007), we define average square error:

 ASE(nn)=CNN∑i|fi(nn)−fi(nn+Δnn)|2 (25)

is the total number of particles. is a normalization constant which is taken as , with being the largest value of force in either runs. is the change in and we choose this to be . With other values of the overall behaviour of the error plot remains qualitatively the same. Figure 2 shows as function of for ATreePM with (solid line) and ATreePM without (dashed line). These were computed for the same clustered distribution of particles. The qualitative behavior of here is same as seen in (Price & Monaghan, 2007). At small ATreePM with has larger errors than ATreePM without . This is also reflected in the poor evolution of density fluctuations with for the ATreePM with , discussed in the next section. Both variants have a minima, the value of error at the minima being lower for ATreePM with the term. Another interesting feature is that the region where errors are small is fairly broad for the ATreePM with the term. For larger the error increases sharply for ATreePM without the term. Thus the optimal configuration is the one with the term and .

## 5 Validation of the Adaptive TreePM code

Cosmological N-Body simulations lack the equivalent of equilibrium distribution functions for haloes, e.g., Plummer halo that may be used to validate a new code. Use of such equilibrium distributions allows one to quantify errors in a clean manner. However, given that the formalism we use is the same as that presented by Price & Monaghan (2007), and that their implementation works well with such tests gives us confidence in the formalism. In the following discussion, we will test the Adaptive TreePM code in a variety of ways and look for numerical convergence. We compare the performance of the Adaptive TreePM with the fixed resolution TreePM. We also study the role of the term and check if dropping this term leads to any degradation in performance of the code. This last point is important as most AMR codes do not have an equivalent term in the equation of motion.

### 5.1 Two-body collisions: plane wave collapse

We start by checking whether the adaptive force softening suppresses two body collisions in the Adaptive TreePM code. We repeat a test recommended by Melott et al. (1997) where they study the collapse of an oblique plane wave. In this test, the collapse should not lead to any transverse motions if the evolution is collisionless. Melott et al. (1997) had shown that codes where the force softening length is much smaller than the inter-particle separation are collisional and lead to generation of significant transverse motions. Some authors (Heitmann et al., 2005) have stated that the failure to ensure planar symmetry may not have a one to one correspondence with collisionality. However, the test is nevertheless an important one and we present the results with the ATreePM code here.

We use exactly the same initial perturbations as used by Melott et al. (1997). The simulations are done with particles and a grid. We choose grid length. Other simulation parameters are as described in the subsection §3.1 on the TreePM code. The output is studied at , following Melott et al. (1997).

We first conduct the test with the fixed resolution TreePM code. Top row of the Figure 3 shows the phase portrait along the direction of collapse for different choices of the force softening length. The softening length varies between in units of the mean inter-particle separation. We see that the phase portrait in the multi-stream region is heavily distorted for the smallest force softening length but is correct for the largest force softening length used here. This reinforces the conclusions of Melott et al. (1997) that using a force softening length that is much smaller that the mean inter-particle separation leads to two body collisions. This point is presented again in the lower row of Figure 3 where we show a scatter plot of modulus of the transverse (to the direction of collapse) velocities and radial velocities. This is shown for a random sub-set of all particles. We see that transverse motions are significant for the simulation with the smallest force softening length, but are under control for the largest force softening length. Thick line in the lower panels connects the average modulus of the transverse velocities in bins of magnitude of radial velocity. The visual impression gathered from the scatter plot is reinforced in that with decreasing force softening length, we get larger transverse motions.

Figure 4 presents results of simulations with the same initial conditions carried out with the Adaptive TreePM code. We show the results for Adaptive TreePM without the term, (left column); with the term, (middle column), and, with the term, (middle column). In general we do not recommend use of due to reasons discussed in the preceeding subsection on errors, discussion in the following subsections, and in Price & Monaghan (2007), but we nevertheless use it in order to look for early signs of two body collisions in a simulation with a relatively small number of particles. The phase portrait for all the three Adaptive TreePM runs is a faithful representation of the expectations. The transverse motions are suppressed strongly, almost to the same level as the TreePM simulation with a force softening of . One of the reasons for this is that the highest overdensities reached in this experiment do not lead to a considerable reduction in the force softening length. The other reason is that the adaptive nature of the code leads to presence of a sufficient number of neighbours within a force softening length and hence the anisotropy in the transverse force is reduced substantially.

Thick lines in the lower panels show the modulus of the average transverse velocities in a few bins of the modulus of the longitudinal velocity. We see that for TreePM, the magnitude of transverse motions drops rapidly as we increase the force softening length. We also note that for the adaptive TreePM, the magnitude of transverse motions is similar to that seen with the TreePM when for the TreePM coincides with the for the ATreePM.

We conclude the discussion of this test by noting that our requirement of , where is the local inter-particle separation, ensures collisionless behavior in evolution of the system.

### 5.2 Convergence with nn and relevance of the ∇ϵ term

For the following discussions we run a power-law model with index . The spectrum is normalized at a an epoch when the scale of non-linearity in grid units. The scale of non-linearity is defined as the scale at which the linearly extrapolated mass variance, defined using a top hat filter, is unity . The simulations are done with particles and a grid. We choose grid length. Other simulation parameters are as described in the subsection §3.1 on the TreePM code. We assume that the Einstein-de Sitter model describes the background universe. In this case self-similar evolution of quantities for power law models provides an additional test for simulation results.

#### 5.2.1 Clustering Properties

We compute the volume averaged point correlation function .

 ¯ξ(r)=3r3∫r0ξ(x)x2dx (26)

is the two point correlation function (Peebles, 1980). To compute from simulation output, we take independent random subsets of particles each. We then estimate the average value of over these subsets. The maximum and the minimum values of in these subsets are our estimate of the errors. Figure 5 shows at an epoch when . TreePM (left column) with , , , , are shown in black, red, blue, ochre and magenta respectively. ATreePM with (middle column) and ATreePM without (right column) with , , are shown in blue, red, black lines respectively. The lower row is a zoom-in of so as to highlight the differences at small scales, due to the variation in and in different runs.

One can show that we require a minimum of neighbors to solve for Eqn.(20). The term is important in over-dense regions where , and is a rapidly varying function of density, one therefore requires a reasonably large to compute it accurately. Use of a small leads to a noisy estimate of and Figure 5 illustrates this point555Also see Figure 2.. With , deviates significantly from our expectations. Disagreement is worse for ATreePM with , as this is the case where a poor estimate of the extra term leads to larger errors. If we use a larger , we expect to see convergence at some point. Curves for and agree with each other (within error-bars) for ATreePM with indicating that evaluation of the extra term is stable. It also indicates that this extra term compensates for the use of a larger number of particles, otherwise typical softening length is expected to be larger for higher and this should lead to lowering of the clustering amplitude at small scales. This is very clearly illustrated for TreePM (left column in Figure 5) where increasing reduces monotonically. We also see this effect for ATreePM without the term where the amplitude of clustering is much smaller for as compared to and the difference between the two is more than twice that seen for ATreePM with the term. Such a behavior is expected, and was also seen by Price & Monaghan (2007). They noted that the increase in force softening length with leads to a bias in the force at small scales. The term corrects for this and biasing of force is less important. Since the distribution of particles gets more strongly clustered with time, we expect over-softening of the force field to degrade further evolution for the TreePM and the ATreePM without , beyond a certain epoch.

We have plotted the scaled 2-point volume averaged correlation at small scales in Figure 6 for TreePM with (left panel), ATreePM with the term (middle panel) and ATreePM without the term (right panel). We used simulations with for both the ATreePM runs shown here. has been computed at epochs when the scale of non-linearity , and (black, red and blue lines respectively). Since the only scale in power-law models is the scale of non-linearity , introduced by gravity, one expects to evolve in a self-similar manner. The scale of self-similarity, , for a given epoch is is determined by finding the scale at which matches with the of the earlier epoch within the error bars computed in the manner described above.

We illustrate numerical convergence properties of the three codes in table 1. Here, we list the scale and at different epochs. For ATreePM is the scale beyond which with the given matches with a reference with . For TreePM this is the scale where with a given force softening length matches with the reference with .

We find that both variants of ATreePM converge with time, the convergence being more rapid for ATreePM with as compared to ATreePM without . The scale of convergence is also smaller for ATreePM with . On the other hand self-similar behavior is degraded with evolution for ATreePM without as increases with time, whereas for ATreePM with the term the evolution is self-similar over a larger range of scales.

Thus we may conclude that the ATreePM with the has well defined numerical convergence as we vary , and the evolution of power law models for this code is self-similar over a wide range of scales. This sets it apart from the ATreePM without the term, where the convergence with is not well defined and the evolution of a power law model is self-similar over a smaller range of scales. Better match with the fixed resolution TreePM code is an added positive feature for the code with the term. This raises obvious questions about AMR codes, where no such term is taken into account as we increase the resolution of the code. We comment on this issue in the Discussion section.

We briefly comment on the convergence and self-similar evolution in the fixed resolution TreePM code. As seen in Table 1, we find that at early times the scale above which we have self-similar evolution is almost the same for . This may indicate discreteness noise, or it may be due to two body collisions during early collapse (Splinter et al., 1998; Melott et al., 1997; Joyce, Marcos, & Baertschiger, 2008; Romeo et al., 2008). At late times, the scale above which evolution is self-similar is . This indicates that any transients introduced at early times by discreteness, etc. have been washed out by the transfer of power from large scales to small scales (Little, Weinberg, & Park, 1991; Evrard & Crone, 1992; Bagla & Padmanabhan, 1997; Bagla & Prasad, 2008). Thus one gets an improved self-similar evolution at late times in the fixed resolution TreePM with the use of a smaller softening length.

On the other hand we see that for TreePM, increases with time. As defined above, this is the scale at which obtained with a given value of matches with the value obtained using . We find that at the last epoch , whereas at early times the limit is . This trend is contrary to that seen with the ATreePM for which decreases with time, the rate at which it comes down being faster for the ATreePM with the . We would like to point out that we do not probe for scales as the number of pairs with a smaller separation is often too small to get reliable estimates of .

Another remarkable fact is that for the TreePM code, . Thus we have the desired behavior in terms of evolution even though we are yet to achieve numerical convergence at the relevant scales. Such a problem does not exist for the ATreePM code.

Lastly, as we show below, the ATreePM code is very good at resolving highly over-dense regions very well. This would appear to be in contradiction with the slightly lower two point correlation function. The reason for the slightly lower correlation function is that the ATreePM code does not resolve small haloes, in particular those with fewer than particles (see the discussion of mass function of collapsed haloes). Indeed, the number density of small haloes is severely under-estimated in the ATreePM simulations and we believe that this is the main reason for a weaker two point correlation function.

#### 5.2.2 Collapsed Haloes

We now look at another global indicator, namely the mass function of collapsed halos, to compare the performance of the TreePM and the ATreePM codes. The number density of halos in the mass range is plotted in Figure 7. Halos were identified using the Friends-Of-Friends algorithm with a linking length in grid units and only haloes with at least particles were considered.

For TreePM we see that as we increase the force softening length there is a clear change in the mass function. The mass function converges for all masses with . However, the mass function for shows fewer haloes of up to about particles. Thus the effect of a large softening length is seen in mass function up to fairly large mass scales. This qualitative behavior is insensitive to the choice of linking length, though the scale of convergence is smaller if we choose a larger linking length.

The ATreePM does not sample the haloes with fewer than members properly. This is expected as such over-densities are not sufficient to reduce the softening length from . At larger scales we have a convergence for and , where the mass function agrees with that for TreePM with . Mass function for deviates significantly from the other two for ATreePM, deviation being larger for the ATreePM with the term. We believe that this is due to the noisy estimation of the force softening length and the term, and resulting errors in force.

As a reference the Press-Schecter mass function is also plotted in all three panels (black dot-dashed line). However it has been shifted up vertically by multiplying by a factor of for ease of comparison.

Distribution of particles in a simulation is a useful representation for comparing large scale features. The distribution of particles in a thin slice is shown in Figure 8. We have shown the distribution for four simulations: TreePM with (top-left), TreePM with (top-right), ATreePM with the term and (lower-left) and ATreePM without the term and (lower-right). The distribution of particles is for the last epoch, corresponding to grid lengths. This slice contains the largest halo in our simulation on the top left corner. We note that the large scale distribution is the same in all cases, indicating that the Adaptive TreePM is not changing anything at large scales.

We zoom in on the region around the largest halo in Figure 9, where the distribution of particles is shown from the set of simulations used in Figure 8. Distribution of mass at scales larger than a gid length is the same in all simulations but at small scales we begin to see some differences between the distribution of particles in different simulations. TreePM with differs most from the other three, in that it does not resolve small scale structures. This happens due to the relatively large force softening length that inhibits collapse if the expected size of the collapsed halo is smaller than . The notable differences between the TreePM () and the ATreePM slices are as follows:

• Small haloes in the region away from the large haloes are more compact for TreePM than for ATreePM. This is perhaps caused by the ATreePM not having a constant and it is likely that in these clumps the value of is larger than that is used in the TreePM.

• Location of sub-structure in the central halo differs somewhat between the different simulations.

Figure 10 shows the inner regions of the halo seen in Figure 9. We show the inner parts of the halo in the four cases (TreePM (), ATreePM with the term (), ATreePM without the term (). The top row shows all the particles in the inner parts of the halo. The second row shows all the particles with an SPH over-density estimate (computed with and using the kernel specified in Eqn. 5) of greater than , the third row shows the particles with an over-density of greater than , and, the last row shows particles with an overdensity larger than . We see that the ATreePM with the term shows the most pronounced central part of the halo in the top panel, and this impression gains strength as we move to lower panels. Table 2 gives us the number of particles in the various panels of Fig. 10. TreePM with manages to trace some of the highly over-dense substructure seen in the ATreePM simulations, though with a much smaller number of particles in these clumps. The ATreePM with the term does better than the one without, particularly at the highest over-densities used here. Indeed ATreePM with the term retains nearly times as many particles in the halo core as compared to TreePM and ATreePM without the term. The same thing is also reflected in Figure 5 where the ATreePM without the term is seen to underestimate clustering at small scales. It is noteworthy that the size of highly over-dense structures in the core of this halo are much bigger than the force softening length for the TreePM.

It is notable that the ATreePM is able to resolve highly over-dense regions while taking much less time than the fixed resolution TreePM. Thus we have the added performance at the cost of fewer resources.

## 6 Discussion

In this paper we have introduced the first cosmological N-Body code that has a continuously varying force resolution. This is based on a well defined formalism that ensures energy and momentum conservation. An important aspect of this formalism is that it modifies the equation of motion at scales below the force softening length. Our implementation of the formalism uses methods developed for SPH codes, this is required as we need to assign quantities like density to particles. We have described our implementation of the code in detail here. We have given a brief summary of errors in force for the Adaptive TreePM code with the new parameter . Given that this code is based on the TreePM code, it inherits the errors in force with respect to all the other parameters. As we have seen in earlier work (Bagla, 2002; Bagla & Ray, 2003; Khandai & Bagla, 2008), errors in force can be minimised to a fairly low level for the TreePM with a judicious choice of parameters.

We find that the time taken by the Adaptive TreePM (for ) is comparable to that taken by the fixed resolution TreePM code if the latter uses a force softening length that is about of the mean inter-particle separation.

One of the key reasons for considering adaptive resolution is to avoid two body collisions. Cosmological simulations require collisionless evolution. We test the Adaptive TreePM by simulating collapse of an oblique plane-wave. We find that unlike the fixed resolution TreePM code that leads to large transverse motions, a clear sign of two body collisions, the Adaptive TreePM code does not have any significant transverse motions. (See Figures 3, 4)

We have compared the performance of the Adaptive TreePM with that of the fixed resolution TreePM code for power law initial conditions. We used a power law model and the Einstein-de Sitter background as it allows us to test codes by requiring self-similar evolution of the correlation function. We find that for , the resulting density distribution in the Adaptive TreePM is comparable to that seen in the highest resolution TreePM in many respects. The correlation function for the two match at all scales and the mass function of collapsed haloes matches for haloes with more than particles. The Adaptive TreePM with the term in the equation of motion does much better than the TreePM as well as the Adaptive TreePM without the extra term in resolving highly over-dense cores of massive haloes. It is noteworthy that the Adaptive code takes less time than the TreePM it has been compared with (see fig. 1).

We have tested the codes by looking for convergence in results as we modify the key parameters that describe the force resolution. We find that the fixed resolution TreePM code converges slower than expected from, for example, self-similar evolution of clustering. The Adaptive TreePM code with the term converges fairly quickly and offers an effective dynamical range that is slightly smaller than the fixed resolution TreePM that takes almost the same time to run.

Given our analysis of errors, and the comparison of the performance of the ATreePM code with and without the extra term in the equation of motion, the natural conclusion is that we require the extra term in order to obtain low errors and numerical convergence. This raises the obvious question about the AMR codes where no such extra term is used. Further, in most AMR codes, resolution is increased when there are of order particles in the lower resolution elements. We find that the errors are minimized when this number is around even when the extra term is not taken into account. It is not clear how serious these issues are, given that AMR codes have been developed and tested in a variety of ways over the last three decades, but it is a concern666One aspect where AMR codes do relatively poorly is in getting the halo mass function at the low mass end (O’Shea et al., 2005; Heitmann et al., 2008). As one can see in figures in the papers cited here, the shortfall is often spread over a decade in mass of haloes. We do not see any gradual decline in the number density of haloes in ATreePM up to the cutoff set by ..

We have established in this paper that an adaptive resolution code for evolution of perturbations in collisionless dark matter can give reliable results for a range of indicators from clustering properties to mass function of collapsed haloes and even get the internal dynamics of collapsed haloes tight. Further, we show that such a code is efficient in that it is faster than fixed resolution codes that give us comparable force resolution in highly over dense regions. Such a code is very useful as it allows us to probe clustering at small scales in a reliable manner. Studies of box-size effects have shown that large simulation boxes are required in order to limit the effect of perturbations at larger scales that are not taken into account (Bagla & Prasad, 2006; Bagla, Prasad, & Khandai, 2008). In such a situation an adaptive code provides us with a reasonable range of scales over which the results can be trusted. We expect to use this code to address several issues related to gravitational clustering in an expanding universe. We also plan to revisit issues related to density profiles of collapsed haloes.

Given the conclusions listed above, we feel that the Adaptive TreePM methods represents an exciting development where we can set aside worries about the impact of collisionality. The relative speed of the adaptive code also makes this a more pragmatic option.

## Acknowledgments

The authors are indebted to Daniel J. Price for clarifying several issues regarding their implementation of the adaptive code. The authors thank Martin White, Volker Springel and Alessandro Romeo for useful discussions and comments. Numerical experiments for this study were carried out at cluster computing facility in the Harish-Chandra Research Institute (http://cluster.hri.res.in). This research has made use of NASA’s Astrophysics Data System.

## References

• Athanassoula et al. (2000) Athanassoula E., Fady E., Lambert J. C., Bosma A., 2000, MNRAS, 314, 475
• Bagla & Padmanabhan (1994) Bagla J. S., Padmanabhan T., 1994, MNRAS, 266, 227
• Bagla & Padmanabhan (1997) Bagla J. S., Padmanabhan T., 1997, Pramana, 49, 161
• Bagla & Padmanabhan (1997) Bagla J. S., Padmanabhan T., 1997, MNRAS, 286, 1023
• Bagla (2002) Bagla J. S., 2002, JApA, 23, 185
• Bagla & Ray (2003) Bagla J. S., Ray S., 2003, NewA, 8, 665
• Bagla & Ray (2005) Bagla J. S., Ray S., 2005, MNRAS, 358, 1076
• Bagla (2005) Bagla J. S., 2005, CSci, 88, 1088
• Bagla & Prasad (2006) Bagla J. S., Prasad J., 2006, MNRAS, 370, 993
• Bagla & Prasad (2008) Bagla J. S., Prasad J., 2008, arXiv, arXiv:0802.2796
• Bagla, Prasad, & Khandai (2008) Bagla J. S., Prasad J., Khandai N., 2008, arXiv, arXiv:0804.1197
• Barnes & Hut (1986) Barnes J., Hut P., 1986, Natur, 324, 446
• Barnes (1990) Barnes J. E., 1990, JCoPh, 87, 161
• Bernardeau et al. (2002) Bernardeau F., Colombi S., Gaztañaga E., Scoccimarro R., 2002, PhR, 367, 1
• Bertschinger (1998) Bertschinger E., 1998, ARA&A, 36, 599
• Binney & Knebe (2002) Binney J., Knebe A., 2002, MNRAS, 333, 378
• Binney (2004) Binney J., 2004, MNRAS, 350, 939
• Bode & Ostriker (2003) Bode P., Ostriker J. P., 2003, ApJS, 145, 1
• Bode, Ostriker, & Xu (2000) Bode P., Ostriker J. P., Xu G., 2000, ApJS, 128, 561
• Bouchet, Adam, & Pellat (1985) Bouchet F. R., Adam J.-C., Pellat R., 1985, A&A, 144, 413
• Bouchet & Kandrup (1985) Bouchet F. R., Kandrup H. E., 1985, ApJ, 299, 1
• Bouchet & Hernquist (1988) Bouchet F. R., Hernquist L., 1988, ApJS, 68, 521
• Brainerd, Scherrer, & Villumsen (1993) Brainerd T. G., Scherrer R. J., Villumsen J. V., 1993, ApJ, 418, 570
• Brieu & Evrard (2000) Brieu P. P., Evrard A. E., 2000, NewA, 5, 163
• Brieu, Summers, & Ostriker (1995) Brieu P. P., Summers F. J., Ostriker J. P., 1995, ApJ, 453, 566
• Couchman (1991) Couchman H. M. P., 1991, ApJ, 368, L23
• Couchman, Thomas, & Pearce (1995) Couchman H. M. P., Thomas P. A., Pearce F. R., 1995, ApJ, 452, 797
• Davis & Peebles (1977) Davis M., Peebles P. J. E., 1977, ApJS, 34, 425
• Dehnen (2000) Dehnen W., 2000, ApJ, 536, L39
• Dehnen (2001) Dehnen W., 2001, MNRAS, 324, 273
• Dehnen (2002) Dehnen W., 2002, JCoPh, 179, 27
• Diemand et al. (2004) Diemand J., Moore B., Stadel J., Kazantzidis S., 2004, MNRAS, 348, 977
• Dolag et al. (2008) Dolag K., Borgani S., Schindler S., Diaferio A., Bykov A. M., 2008, SSRv, 134, 229
• Dubinski (1996) Dubinski J., 1996, NewA, 1, 133
• Dubinski et al. (2004) Dubinski J., Kim J., Park C., Humble R., 2004, NewA, 9, 111
• Dyer & Ip (1993) Dyer C. C., Ip P. S. S., 1993, ApJ, 409, 60
• Ebisuzaki et al. (1993) Ebisuzaki T., Makino J., Fukushige T., Taiji M., Sugimoto D., Ito T., Okumura S. K., 1993, PASJ, 45, 269
• Efstathiou et al. (1985) Efstathiou G., Davis M., White S. D. M., Frenk C. S., 1985, ApJS, 57, 241
• El-Zant (2006) El-Zant A. A., 2006, MNRAS, 370, 1247
• Evrard (1988) Evrard A. E., 1988, MNRAS, 235, 911
• Evrard & Crone (1992) Evrard A. E., Crone M. M., 1992, ApJ, 394, L1
• Efstathiou et al. (1988) Efstathiou G., Frenk C. S., White S. D. M., Davis M., 1988, MNRAS, 235, 715
• Gingold &Monaghan (1977) Gingold R. A., Monaghan J. J., 1977, MNRAS, 181, 375
• Gingold & Monaghan (1982) Gingold R. A., Monaghan J. J., 1982, JCoPh, 46, 429
• Greengard & Rokhlin (1987) Greengard L., Rokhlin V., 1987, JCoPh, 73, 325
• Gurbatov, Saichev, & Shandarin (1989) Gurbatov S. N., Saichev A. I., Shandarin S. F., 1989, MNRAS, 236, 385
• Hamana, Yoshida, & Suto (2002) Hamana T., Yoshida N., Suto Y., 2002, ApJ, 568, 455
• Hamilton et al. (1991) Hamilton A. J. S., Kumar P., Lu E., Matthews A., 1991, ApJ, 374, L1
• Heitmann et al. (2005) Heitmann K., Ricker P. M., Warren M. S., Habib S., 2005, ApJS, 160, 28
• Heitmann et al. (2008) Heitmann K., et al., 2008, CS&D, 1, 015003
• Hernquist (1987) Hernquist L., 1987, ApJS, 64, 715
• Hernquist & Katz (1989) Hernquist L., Katz N., 1989, ApJS, 70, 419
• Hernquist (1990) Hernquist L., 1990, JCoPh, 87, 137
• Hernquist & Barnes (1990) Hernquist L., Barnes J. E., 1990, ApJ, 349, 562
• Hernquist, Bouchet, & Suto (1991) Hernquist L., Bouchet F. R., Suto Y., 1991, ApJS, 75, 231
• Hockney & Eastwood (1988) Hockney R. W., Eastwood J. W., 1988, Computer Simulation using Particles, McGraw-Hill
• Hui & Bertschinger (1996) Hui L., Bertschinger E., 1996, ApJ, 471, 1
• Jain, Mo, & White (1995) Jain B., Mo H. J., White S. D. M., 1995, MNRAS, 276, L25
• Jernigan & Porter (1989) Jernigan J. G., Porter D. H., 1989, ApJS, 71, 871
• Joyce, Marcos, & Baertschiger (2008) Joyce M., Marcos B., Baertschiger T., 2008, arXiv, arXiv:0805.1357
• Kanekar (2000) Kanekar N., 2000, ApJ, 531, 17
• Kawai & Makino (2001) Kawai A., Makino J., 2001, ApJ, 550, L143
• Klypin & Shandarin (1983) Klypin A. A., Shandarin S. F., 1983, MNRAS, 204, 891
• Knebe, Green, & Binney (2001) Knebe A., Green A., Binney J., 2001, MNRAS, 325, 845
• Knebe et al. (2000) Knebe A., Kravtsov A. V., Gottlöber S., Klypin A. A., 2000, MNRAS, 317, 630
• Khandai & Bagla (2008) Khandai N., Bagla J. S., 2008, arXiv:0802.3215 [astro-ph]
• Kravtsov, Klypin, & Khokhlov (1997) Kravtsov A. V., Klypin A. A., Khokhlov A. M., 1997, ApJS, 111, 73
• Kuhlman, Melott, & Shandarin (1996) Kuhlman B., Melott A. L., Shandarin S. F., 1996, ApJ, 470, L41
• Little, Weinberg, & Park (1991) Little B., Weinberg D. H., Park C., 1991, MNRAS, 253, 295
• Lucy (1977) Lucy L. B., 1977, AJ, 82, 1013
• Ma (1998) Ma C.-P., 1998, ApJ, 508, L5
• Macfarland et al. (1998) Macfarland T., Couchman H. M. P., Pearce F. R., Pichlmeier J., 1998, NewA, 3, 687
• Makino (1991) Makino J., 1991, PASJ, 43, 621
• Makino (2004) Makino J., 2004, PASJ, 56, 521
• Makino (2002) Makino J., 2002, NewA, 7, 373
• Makino (1990) Makino J., 1990, JCoPh, 87, 148
• Makino et al. (2003) Makino J., Fukushige T., Koga M., Namura K., 2003, PASJ, 55, 1163
• Matarrese et al. (1992) Matarrese S., Lucchin F., Moscardini L., Saez D., 1992, MNRAS, 259, 437
• Melott et al. (1997) Melott A. L., Shandarin S. F., Splinter R. J., Suto Y., 1997, ApJ, 479, L79
• Merritt (1996) Merritt D., 1996, AJ, 111, 2462
• Merz, Pen, & Trac (2005) Merz H., Pen U.-L., Trac H., 2005, NewA, 10, 393
• Miller (1983) Miller R. H., 1983, ApJ, 270, 390
• Monaghan & Lattanzio (1985) Monaghan J. J., Lattanzio J. C., 1985, A&A, 149, 135
• Monaghan (1992) Monaghan J. J., 1992, ARA&A, 30, 543
• Nityananda & Padmanabhan (1994) Nityananda R., Padmanabhan T., 1994, MNRAS, 271, 976
• O’Shea et al. (2005) O’Shea B. W., Nagamine K., Springel V., Hernquist L., Norman M. L., 2005, ApJS, 160, 1
• Padmanabhan (2002) Padmanabhan T., 2002, Theoretical Astrophysics, Volume III: Galaxies and Cosmology. Cambridge University Press.
• Padmanabhan et al. (1996) Padmanabhan T., Cen R., Ostriker J. P., Summers F. J., 1996, ApJ, 466, 604
• Peacock & Dodds (1996) Peacock J. A., Dodds S. J., 1996, MNRAS, 280, L19
• Peacock & Dodds (1994) Peacock J. A., Dodds S. J., 1994, MNRAS, 267, 1020
• Peacock (1999) Peacock J. A., 1999, Cosmological Physics, Cambridge University Press
• Peebles (1980) Peebles P. J. E., 1980, The Large-Scale Structure of the Universe, Princeton University Press
• Peebles (1974) Peebles P. J. E., 1974, A&A, 32, 391
• Percival et al. (2007) Percival W. J., et al., 2007, ApJ, 657, 645
• Power & Knebe (2006) Power C., Knebe A., 2006, MNRAS, 370, 691
• Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
• Ray & Bagla (2004) Ray S., Bagla J. S., 2004, astro, arXiv:astro-ph/0405220
• Romeo et al. (2008) Romeo A. B., Agertz O., Moore B., Stadel J., 2008, ApJ, 686, 1
• Sahni & Coles (1995) Sahni V., Coles P., 1995, PhR, 262, 1
• Salmon & Warren (1994) Salmon J. K., Warren M. S., 1994, JCoPh, 111, 136
• Smith et al. (2003) Smith R. E., et al., 2003, MNRAS, 341, 1311
• Spergel et al. (2007) Spergel D. N., et al., 2007, ApJS, 170, 377
• Splinter et al. (1998) Splinter R. J., Melott A. L., Shandarin S. F., Suto Y., 1998, ApJ, 497, 38
• Springel, Yoshida, & White (2001) Springel V., Yoshida N., White S. D. M., 2001, NewA, 6, 79
• Springel & Hernquist (2002) Springel V., Hernquist L., 2002, MNRAS, 333, 649
• Springel (2005) Springel V., 2005, MNRAS, 364, 1105
• Springel et al. (2005) Springel V., et al., 2005, Natur, 435, 629
• Suisalu & Saar (1995) Suisalu I., Saar E., 1995, MNRAS, 274, 287
• Szebehely & Peters (1967) Szebehely V., Peters C. F., 1967, AJ, 72, 1187
• Thacker & Couchman (2006) Thacker R. J., Couchman H. M. P., 2006, CoPhC, 174, 540
• Theis (1998) Theis C., 1998, A&A, 330, 1180
• Theuns (1994) Theuns T., 1994, CoPhC, 78, 238
• Thiébaut et al. (2008) Thiébaut J., Pichon C., Sousbie T., Prunet S., Pogosyan D., 2008, MNRAS, 387, 397
• Xu (1995) Xu G., 1995, ApJS, 98, 355
• Yoshikawa & Fukushige (2005) Yoshikawa K., Fukushige T., 2005, PASJ, 57, 849
• Zel’Dovich (1970) Zel’Dovich Y. B., 1970, A&A, 5, 84
• Zhan (2006) Zhan H., 2006, ApJ, 639, 617

## Appendix A Expressions for Cubic Spline Softened Potential and Force Kernels

One can integrate the Poisson equation for a given kernel to obtain the softened two-body potential kernel. For the case of the cubic spline kernel we have the expression for the two-body potential kernel and the regular two-body force :

 ϕ(u,ϵ) = ⎧⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪⎩16ϵ[13(uϵ)2−35(uϵ)4+25(uϵ)5]−145ϵ,0≤uϵ<0.51ϵ[115(ϵu)+323(uϵ)2−161(uϵ)3+485(uϵ)4−3215(uϵ)5]−165ϵ,0.5≤uϵ<1.0−1u,1.0≤uϵ f(u,ϵ) = (28)

For variable or local smoothing length there will be an additional term given in eqs.(10-12) for which one needs to compute , and For the cubic spline kernel these expressions are:

 ∂W∂ϵ = 8πϵ4⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩3[−1+10(uϵ)2−12(uϵ)3],0≤uϵ<0.56[−1+4(uϵ)−5(uϵ)2+2(uϵ)3],0.5≤uϵ<1.00,1.0≤uϵ (29) ∂W∂u = (30) ∂ϕ∂ϵ = ⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩16ϵ2[−(uϵ)2+3(uϵ)4−125(uϵ)5]+145ϵ20≤uϵ<0.51ϵ2[−32(uϵ)2+64(uϵ)3−48(uϵ)4+645(uϵ)5]+165ϵ2% 0.5≤uϵ<1.00,1.0≤uϵ (31)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters