The Adaptive TreePM: An Adaptive Resolution Code for Cosmological body Simulations
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 interparticle 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 interparticle 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 overdense 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: NBody simulations, cosmology: large scale structure of the universe1 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 overdense 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 nonlinear 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. NBody 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 quasilinear 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 NBody codes should ensure the following:

The universe does not have a boundary. Therefore cosmological simulations need to be run with periodic boundary conditions^{1}^{1}1An 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 NBody simulation represents a very large number of particles/objects in the universe. Thus we must ensure that pairwise interaction of NBody particles is softened at scales comparable with the local interparticle 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 NBody codes suffer from collisionality or force biasing. Force is biased when softening lengths are much larger than the local mean interparticle 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 overdense 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 interparticle separation as a function of space and time. Price & Monaghan (2007) have developed an energymomentum 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; ElZant, 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 Adaptive Force Softening
2.1 Formalism
The aim of any collisionless body code is to self consistently evolve the phasespace distribution function (DF) under its own gravitational force field :
(1)  
(2) 
The approach that one takes is to sample the DF, by phasespace 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 twobody collisions arising out of interactions between particles which were generated for sampling the density field. One therefore assigns a finite size to NBody 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,
(3) 
is now smoothed at small scales if we assign a finite size to every particle:
(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)).
(5) 
In the context of individual softening lengths the density at the location of the particle is given by:
(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 NBody simulations. We study the impact of ignoring this term in Cosmological simulations.
A remedy for momentum nonconservation is to use a symmetrized softening length
(7) 
and plug it into the expression for density to rederive 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.
(8) 
The total potential is thus
(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.
(10) 
The EOM of motion can be derived with this Lagrangian (the reader is referred to Price & Monaghan (2007) for details)
(11)  
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:
(12)  
(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 NBody method which improves the accuracy and performance of the BarnesHut (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 shortrange and a longerange component at a scale . The PM method is used to calculate longrange force and the shortrange force is computed using the BH Tree method. Use of the BH Tree for shortrange 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.
(14)  
(15) 
Here and are the shortrange and longrange 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 longrange force is computed in Fourier space with the PM method and the shortrange force is computed in real space with the Tree method. The short range potential and force in real space are:
(16)  
(17)  
(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 shortrange 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 (subcells) within them. Subcells can have even smaller cells within them, or they can contain a particle. In three dimensions, each cubic cell is divided into eight cubic subcells. Cells, as structures, have attributes like total mass, location of center of mass and pointers to subcells. 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, subcells 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
(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:
(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 shortrange 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 nonmodification 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 overdense regions to dominate the CPU time requirements. The value of the lower bound must correspond to densities that are much higher than highest overdensities 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 shortrange force so that the individual softening lengths needed for shortrange 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 shortrange 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 shortrange 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 shortrange 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 shortrange force instead of constructing one global data structure for the entire simulation volume for computing the shortrange 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 GADGET2 Springel (2005), in which particle trajectories are integrated with individual timesteps and synchronized with the largest timestep. As we allow the block time step^{2}^{2}2Same as the largest time step. to vary with time, we work with the so called KDK approach (KickDriftKick) 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 (DriftKickDrift) (see Springel (2005) for details.). We give separate PM (long range, global) kicks and Tree (short range, individual) kicks^{3}^{3}3Tree kick includes the contribution of the term.. The block timestep is determined by the particle which has the maximum PM acceleration :
(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:
(22) 
and then the appropriate hierarchy is chosen depending on this value. Here and are the modulus of the individual shortrange acceleration (sum of the and terms) and the softening lengths respectively. TreePM has a similar timestepping 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 scales^{4}^{4}4As 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 intragroup 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 particleparticle (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.
(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 nonlinearity 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 timestepping criterion as Eqn.(22) for TreePM.
(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 overdense regions. This saves time in the underdense and not so overdense regions and the ATreePM code devotes more time to evolving trajectories in highly overdense 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 lowerleft 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 nonlinear 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 overdense 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:
(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 NBody 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 Twobody 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 interparticle 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 interparticle separation. We see that the phase portrait in the multistream 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 interparticle 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 subset 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 interparticle separation, ensures collisionless behavior in evolution of the system.
5.2 Convergence with and relevance of the term
For the following discussions we run a powerlaw model with index . The spectrum is normalized at a an epoch when the scale of nonlinearity in grid units. The scale of nonlinearity 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 Einsteinde Sitter model describes the background universe. In this case selfsimilar 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 .
(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 zoomin 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 overdense 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 point^{5}^{5}5Also 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 errorbars) 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 oversoftening of the force field to degrade further evolution for the TreePM and the ATreePM without , beyond a certain epoch.
We have plotted the scaled 2point 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 nonlinearity , and (black, red and blue lines respectively). Since the only scale in powerlaw models is the scale of nonlinearity , introduced by gravity, one expects to evolve in a selfsimilar manner. The scale of selfsimilarity, , 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.
Runs  

TreePM ()  1.2  0.96  0.61  1.1  0.67 
TreePM ()  0.61  0.44  0.36  0.42  0.15 
TreePM ()  0.25  0.20  0.10  0.18  0.16 
TreePM ()  0.10  0.10  0.10  0.10  0.17 
TreePM ()        0.10  0.18 
ATreePM ()  0.10  0.15  0.30  0.27  0.31 
ATreePM ()        0.42  0.47 
ATreePM ()  0.27  0.28  0.38  0.54  0.42 
ATreePM ()        0.59  0.50 
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 selfsimilar behavior is degraded with evolution for ATreePM without as increases with time, whereas for ATreePM with the term the evolution is selfsimilar 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 selfsimilar 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 selfsimilar 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 selfsimilar evolution in the fixed resolution TreePM code. As seen in Table 1, we find that at early times the scale above which we have selfsimilar 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 selfsimilar 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 selfsimilar 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 overdense 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 underestimated 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 FriendsOfFriends 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 overdensities 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 PressSchecter mass function is also plotted in all three panels (black dotdashed 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 (topleft), TreePM with (topright), ATreePM with the term and (lowerleft) and ATreePM without the term and (lowerright). 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 substructure 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 overdensity estimate (computed with and using the kernel specified in Eqn. 5) of greater than , the third row shows the particles with an overdensity 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 overdense 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 overdensities 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 overdense 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 overdense regions while taking much less time than the fixed resolution TreePM. Thus we have the added performance at the cost of fewer resources.
Runs  
TreePM ()  320  106  18 
ATreePM ()  568  297  125 
ATreePM ()  374  169  18 
5.2.3 Dynamics within Collapsed Haloes
The equation of motion for the adaptive code has the additional term, and it is important to test whether this term leads to a change in dynamics in highly overdense regions. We test this by studying the dynamics in central cores of the largest haloes in the simulation. We study the ratio at a scale corresponding to , with the averaging done around the centre of mass of the halo. We chose the reference scale in this manner as the density profiles for haloes with different softening length differ considerably and using a high density contour to define the radius leads to very different physical scales. Figure 11 shows as a function of softening length for the fixed resolution TreePM simulations. This has been done for the five largest haloes in the simulation. We see that as we go to a larger softening length, the ratio becomes larger than unity, indeed for the ratio is closer to two for one of the haloes. This is likely to happen if the size of the overdense core is comparable to the softening length, and indeed this is the case for some of the haloes used for this plot. The same plot shows the ratio as seen in the adaptive code (with the term) with empty squares, and as filled triangles for the adaptive code without the term. These circles are plotted at small values of , not used for TreePM runs. The values of used for positioning these symbols in the plot has no relevance to the simulation, and this has been done purely for the purpose of plotting the values. We see that the values for the ratio in both the adaptive codes correspond closely to the values of the ratio seen in fixed resolution TreePM simulations with a small . Thus we may conclude that the additional term in the Adaptive code is not leading to any significant changes in the dynamics, as compared to the fixed resolution codes.
6 Discussion
In this paper we have introduced the first cosmological NBody 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 interparticle 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 planewave. 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 Einsteinde Sitter background as it allows us to test codes by requiring selfsimilar 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 overdense 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, selfsimilar 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 concern^{6}^{6}6One 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 boxsize 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 HarishChandra 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
 ElZant (2006) ElZant 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, McGrawHill
 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 [astroph]
 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 (1996) Padmanabhan T., 1996, MNRAS, 278, L29
 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 LargeScale 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
 Prasad (2007) Prasad J., 2007, JApA, 28, 117
 Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
 Ray & Bagla (2004) Ray S., Bagla J. S., 2004, astro, arXiv:astroph/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
 Wadsley, Stadel, & Quinn (2004) Wadsley J. W., Stadel J., Quinn T., 2004, NewA, 9, 137
 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 twobody potential kernel. For the case of the cubic spline kernel we have the expression for the twobody potential kernel and the regular twobody force :
(28) 
For variable or local smoothing length there will be an additional term given in eqs.(1012) for which one needs to compute , and For the cubic spline kernel these expressions are:
(29)  
(30)  
(31) 