Multi-level adaptive particle mesh (MLAPM): a c code for cosmological simulations

We present a computer code written in c that is designed to simulate structure formation from collisionless matter. The code is purely grid-based and uses a recursively refined Cartesian grid to solve Poisson's equation for the potential, rather than obtaining the potential from a Green's function. Refinements can have arbitrary shapes and in practice closely follow the complex morphology of the density field that evolves. The time-step shortens by a factor of 2 with each successive refinement.
Competing approaches to N -body simulation are discussed from the point of view of the basic theory of N -body simulation. It is argued that an appropriate choice of softening length ε is of great importance and that ε should be at all points an appropriate multiple of the local interparticle separation. Unlike tree and P³M codes, multigrid codes automatically satisfy this requirement. We show that at early times and low densities in cosmological simulations, ε needs to be significantly smaller relative to the interparticle separation than in virialized regions. Tests of the ability of the code's Poisson solver to recover the gravitational fields of both virialized haloes and Zel'dovich waves are presented, as are tests of the code's ability to reproduce analytic solutions for plane-wave evolution. The times required to conduct a ΛCDM cosmological simulation for various configurations are compared with the times required to complete the same simulation with the ART, AP³M and GADGET codes. The power spectra, halo mass functions and halo–halo correlation functions of simulations conducted with different codes are compared.
The code is available from http://www-thphys.physics.ox.ac.uk/users/MLAPM .     

On the effects of resolution in dissipationless cosmological simulations

We present a study of numerical effects in dissipationless cosmological simulations. The numerical effects are evaluated and studied by comparing the results of a series of 64³-particle simulations of varying force resolution and number of time-steps, performed using three of the N -body techniques currently used for cosmological simulations: the Particle–Mesh (PM), the Adaptive Particle–Particle–Particle–Mesh (AP³M) and the newer Adaptive Refinement Tree (ART) codes. This study can therefore be interesting both as an analysis of numerical effects and as a systematic comparison of different codes.
We find that the AP³M and the ART codes produce similar results given that convergence is reached within the code type. We also find that numerical effects may affect the high-resolution simulations in ways that have not been discussed before. In particular, our study revealed the presence of two-body scattering, the effects of which can be greatly amplified by inaccuracies in time integration. This process appears to affect the correlation function of matter, the mass function, the inner density of dark matter haloes and other statistics at scales much larger than the force resolution, although different statistics may be affected in a different fashion. We discuss the conditions for which strong two-body scattering is possible and discuss the choice of the force resolution and integration time-step. Furthermore, we discuss recent claims that simulations with force softening smaller than the mean interparticle separation are not trustworthy and argue that this claim is incorrect in general, and applies only to the phase-sensitive statistics. Our conclusion is that, depending on the choice of mass and force resolution and the integration time-step, a force resolution as small as 0.01 of the mean interparticle separation can be justified.     

TreePM: A code for cosmological N-body simulations

We describe the TreePM method for carrying out large N-Body simulations to study formation and evolution of the large scale structure in the Universe. This method is a combination of Barnes and Hut tree code and Particle-Mesh code. It combines the automatic inclusion of periodic boundary conditions of PM simulations with the high resolution of tree codes. This is done by splitting the gravitational force into a short range and a long range component. We describe the splitting of force between these two parts. We outline the key differences between TreePM and some other N-Body methods.     

Two-body relaxation in cosmological simulations

It is logically possible that early two-body relaxation in simulations of cosmological clustering influences the final structure of massive clusters. Convergence studies in which mass and spatial resolution are simultaneously increased cannot eliminate this possibility. We test the importance of two-body relaxation in cosmological simulations with simulations in which there are two species of particles. The cases of two mass ratios, √2:1 and 4:1, are investigated. Simulations are run with both a spatially fixed softening length and adaptive softening using the publicly available codes gadget and mlapm , respectively.
The effects of two-body relaxation are detected in both the density profiles of haloes and the mass function of haloes. The effects are more pronounced with a fixed softening length, but even in this case they are not so large as to suggest that results obtained with one mass species are significantly affected by two-body relaxation.
The simulations that use adaptive softening are less affected by two-body relaxation and produce slightly higher central densities in the largest haloes. They run about three times faster than the simulations that use a fixed softening length.     

An energy-conserving formalism for adaptive gravitational force softening in smoothed particle hydrodynamics and N-body codes

In this paper, we describe an adaptive softening length formalism for collisionless N -body and self-gravitating smoothed particle hydrodynamics (SPH) calculations which conserves momentum and energy exactly. This means that spatially variable softening lengths can be used in N -body calculations without secular increases in energy. The formalism requires the calculation of a small additional term to the gravitational force related to the gradient of the softening length. The extra term is similar in form to the usual SPH pressure force (although opposite in direction) and is therefore straightforward to implement in any SPH code at almost no extra cost. For N -body codes, some additional cost is involved as the formalism requires the computation of the density through a summation over neighbouring particles using the smoothing kernel. The results of numerical tests demonstrate that, for homogeneous mass distributions, the use of adaptive softening lengths gives a softening which is always close to the 'optimal' choice of fixed softening parameter, removing the need for fine-tuning. For a heterogeneous mass distribution (as may be found in any large-scale N -body simulation), we find that the errors on the least-dense component are lowered by an order of magnitude compared to the use of a fixed softening length tuned to the densest component. For SPH codes, our method presents a natural and an elegant choice of softening formalism which makes a small improvement to both the force resolution and the total energy conservation at almost zero additional cost.     

TreePM method for two-dimensional cosmological simulations

We describe the two-dimensional TreePM method in this paper. The 2d TreePM code is an accurate and efficient technique to carry out large two-dimensional N-body simulations in cosmology. This hybrid code combines the 2d Barnes and Hut Tree method and the 2d Particle-Mesh method. We describe the splitting of force between the PM and the Tree parts. We also estimate error in force for a realistic configuration. Finally, we discuss some tests of the code.     

Voids in a ΛCDM universe

We study the formation and evolution of voids in the dark matter distribution using various simulations of the popular Λ cold dark matter cosmogony. We identify voids by requiring them to be regions of space with a mean overdensity of −0.8 or less – roughly the equivalent of using a spherical overdensity group finder for haloes. Each of the simulations contains thousands of voids. The distribution of void sizes in the different simulations shows good agreement when differences in particle and grid resolution are accounted for. Voids very clearly correspond to minima in the smoothed initial density field. Apart from a very weak dependence on the mass resolution, the rescaled mass profiles of voids in the different simulations agree remarkably well. We find a universal void mass profile of the form  ρ(< r )/ρ( r _eff) ∝ exp[( r / r _eff)^α]  , where r _eff is the effective radius of a void and  α∼ 2  . The mass function of haloes in voids is steeper than that of haloes that populate denser regions. In addition, the abundances of void haloes seem to evolve somewhat more strongly between redshifts ∼1 and 0 than the global abundances of haloes.     

Performance characteristics of TreePM codes

We present a detailed analysis of the error budget for the TreePM method for doing cosmological N-body simulations. It is shown that the choice of filter for splitting the inverse square force into short and long range components suggested in Bagla [Bagla, J.S., 2002a. JApA 23, 185; Bagla, J.S., 2002b. Preprint. To appear in Proceedings of the Numerical Simulations in Astrophysics, Tokyo, July 2002] is close to optimum. We show that the error in the long range component of the force contributes very little to the total error in force. Errors introduced by the tree approximation for the short range force are different from those for the inverse square force, and these errors dominate the total error in force. We calculate the distribution function for error in force for clustered and unclustered particle distributions. This gives an idea of the error in realistic situations for different choices of parameters of the TreePM algorithm. We test the code by simulating a few power law models and checking for scale invariance.     

Dark matter halo structure in CDM hydrodynamical simulations

We have carried out a comparative analysis of the properties of dark matter haloes in N -body and hydrodynamical simulations. We analyse their density profiles, shapes and kinematical properties with the aim of assessing the effects that hydrodynamical processes might produce on the evolution of the dark matter component. The simulations performed allow us to reproduce dark matter haloes with high resolution, although the range of circular velocities is limited. We find that for haloes with circular velocities of [150–200] km s⁻¹ at the virial radius, the presence of baryons affects the evolution of the dark matter component in the central region, modifying the density profiles, shapes and velocity dispersions. We also analyse the rotation velocity curves of disc-like structures and compare them with observational results.     

The haloes of merger remnants

We perform collisionless N -body simulations of 1:1 galaxy mergers, using models which include a galaxy halo, disc and bulge, focusing on the behaviour of the halo component. The galaxy models are constructed without recourse to a Maxwellian approximation. We investigate the effect of varying the galaxies' orientation, their mutual orbit and the initial velocity anisotropy or cusp strength of the haloes upon the remnant halo density profiles and shape, as well as on the kinematics. We observe that the halo density profile (determined as a spherical average, an approximation we find appropriate) is exceptionally robust in mergers, and that the velocity anisotropy of our remnant haloes is nearly independent of the orbits or initial anisotropy of the haloes. The remnants follow the halo anisotropy – local density slope (β–γ) relation suggested by Hansen & Moore in the inner parts of the halo, but β is systematically lower than this relation predicts in the outer parts. Remnant halo axis ratios are strongly dependent on the initial parameters of the haloes and on their orbits. We also find that the remnant haloes are significantly less spherical than those described in studies of simulations which include gas cooling.     

Gravitational collapse in an expanding background and the role of substructure – II. Excess power at small scales and its effect on collapse of structures at larger scales

We study the interplay of clumping at small scales with the collapse and relaxation of perturbations at larger scales using N -body simulations. We quantify the effect of collapsed haloes on perturbations at larger scales using the two-point correlation function, moments of counts in cells and the mass function. The purpose of the study is twofold and the primary aim is to quantify the role played by collapsed low-mass haloes in the evolution of perturbations at large scales; this is in view of the strong effect seen when the large scale perturbation is highly symmetric. Another reason for this study is to ask whether features or a cut-off in the initial power spectrum can be detected using measures of clustering at scales that are already non-linear. The final aim is to understand the effect of ignoring perturbations at scales smaller than the resolution of N -body simulations. We find that these effects are ignorable if the scale of non-linearity is larger than the average interparticle separation in simulations. Features in the initial power spectrum can be detected easily if the scale of these features is in the linear regime; detecting such features becomes difficult as the relevant scales become non-linear. We find no effect of features in initial power spectra at small scales on the evolved power spectra at large scales. We may conclude that, in general, the effect on the evolution of perturbations at large scales of clumping on small scales is very small and may be ignored in most situations.     

Resolving cosmic structure formation with the Millennium-II Simulation

We present the Millennium-II Simulation (MS-II), a very large N -body simulation of dark matter evolution in the concordance Λ cold dark matter (ΛCDM) cosmology. The MS-II assumes the same cosmological parameters and uses the same particle number and output data structure as the original Millennium Simulation (MS), but was carried out in a periodic cube one-fifth the size  (100  h ⁻¹ Mpc)  with five times better spatial resolution (a Plummer equivalent softening of  1.0  h ⁻¹ kpc  ) and with 125 times better mass resolution (a particle mass of  6.9 × 10⁶  h ⁻¹ M_⊙  ). By comparing results at MS and MS-II resolution, we demonstrate excellent convergence in dark matter statistics such as the halo mass function, the subhalo abundance distribution, the mass dependence of halo formation times, the linear and non-linear autocorrelations and power spectra, and halo assembly bias. Together, the two simulations provide precise results for such statistics over an unprecedented range of scales, from haloes similar to those hosting Local Group dwarf spheroidal galaxies to haloes corresponding to the richest galaxy clusters. The 'Milky Way' haloes of the Aquarius Project were selected from a lower resolution version of the MS-II and were then resimulated at much higher resolution. As a result, they are present in the MS-II along with thousands of other similar mass haloes. A comparison of their assembly histories in the MS-II and in resimulations of 1000 times better resolution shows detailed agreement over a factor of 100 in mass growth. We publicly release halo catalogues and assembly trees for the MS-II in the same format within the same archive as those already released for the MS.     

Multimass spherical structure models for N-body simulations

We present a simple and efficient method to set up spherical structure models for N -body simulations with a multimass technique. This technique reduces by a substantial factor the computer run time needed in order to resolve a given scale as compared to single-mass models. It therefore allows to resolve smaller scales in N -body simulations for a given computer run time. Here, we present several models with an effective resolution of up to  1.68 × 10⁹  particles within their virial radius which are stable over cosmologically relevant time-scales. As an application, we confirm the theoretical prediction by Dehnen that in mergers of collisionless structures like dark matter haloes always the cusp of the steepest progenitor is preserved. We model each merger progenitor with an effective number of particles of approximately 10⁸ particles. We also find that in a core–core merger the central density approximately doubles whereas in the cusp–cusp case the central density only increases by approximately 50 per cent. This may suggest that the central regions of flat structures are better protected and get less energy input through the merger process.     

Simulating subhaloes at high redshift: merger rates, counts and types

Galaxies are believed to be in one-to-one correspondence with simulated dark matter subhaloes. We use high-resolution N -body simulations of cosmological volumes to calculate the statistical properties of subhalo (galaxy) major mergers at high redshift ( z = 0.6–5). We measure the evolution of the galaxy merger rate, finding that it is much shallower than the merger rate of dark matter host haloes at   z > 2.5  , but roughly parallels that of haloes at   z < 1.6  . We also track the detailed merger histories of individual galaxies and measure the likelihood of multiple mergers per halo or subhalo. We examine satellite merger statistics in detail: 15–35 per cent of all recently merged galaxies are satellites, and satellites are twice as likely as centrals to have had a recent major merger. Finally, we show how the differing evolution of the merger rates of haloes and galaxies leads to the evolution of the average satellite occupation per halo, noting that for a fixed halo mass, the satellite halo occupation peaks at   z ∼ 2.5  .     

Time evolution of galactic warps in prolate haloes

A recent observation with the Hipparcos satellite and some numerical simulations imply that the interaction between an oblate halo and a disc is inappropriate for the persistence of galactic warps. Following on from this , we have compared the time evolution of galactic warps in a prolate halo with that in an oblate halo. The haloes were approximated as fixed potentials, while the discs were represented by N -body particles. We have found that the warping in the oblate halo continues to wind up, and finally disappears. On the other hand, for the prolate halo model, the precession rate of the outer disc increases when the precession of the outer disc recedes from that of the inner disc, and vice versa. Consequently, the warping in the prolate halo persisted to the end of the simulation by retaining the alignment of the line of nodes of the warped disc. Therefore, our results suggest that prolate haloes could sustain galactic warps. The physical mechanism of the persistence of warp is discussed on the basis of the torque between a halo and a disc and that between the inner and outer regions of the disc.     

A new parallel PM code for very large-scale cosmological simulations

We have developed a parallel Particle–Particle, Particle–Mesh (P³M) simulation code for the Cray T3E parallel supercomputer that is well suited to studying the time evolution of systems of particles interacting via gravity and gas forces in cosmological contexts. The parallel code is based upon the public-domain serial Adaptive P³M-SPH (http://coho.astro.uwo.ca/pub/hydra/hydra.html) code of Couchman et al. (1995)[ApJ, 452, 797]. The algorithm resolves gravitational forces into a long-range component computed by discretizing the mass distribution and solving Poisson's equation on a grid using an FFT convolution method, and a short-range component computed by direct force summation for sufficiently close particle pairs. The code consists primarily of a particle–particle computation parallelized by domain decomposition over blocks of neighbour-cells, a more regular mesh calculation distributed in planes along one dimension, and several transformations between the two distributions. The load balancing of the P³M code is static, since this greatly aids the ongoing implementation of parallel adaptive refinements of the particle and mesh systems. Great care was taken throughout to make optimal use of the available memory, so that a version of the current implementation has been used to simulate systems of up to 10⁹ particles with a 1024³ mesh for the long-range force computation. These are the largest Cosmological N-body simulations of which we are aware. We discuss these memory optimizations as well as those motivated by computational performance. Performance results are very encouraging, and, even without refinements, the code has been used effectively for simulations in which the particle distribution becomes highly clustered as well as for other non-uniform systems of astrophysical interest.     

A cosmological hydrodynamic code based on the piecewise parabolic method

We present a hydrodynamical code for cosmological simulations which uses the piecewise parabolic method (PPM) to follow the dynamics of the gas component and an N -body particle–mesh algorithm for the evolution of the collisionless component. The gravitational interaction between the two components is regulated by the Poisson equation which is solved by a standard fast Fourier transform (FFT) procedure. In order to simulate cosmological flows we have introduced several modifications to the original PPM scheme which we describe in detail. Various tests of the code are presented including adiabatic expansion, single and multiple pancake formation and three-dimensional cosmological simulations with initial conditions based on the cold dark matter scenario.     

grommet: an N-body code for high-resolution simulations of individual galaxies

This paper presents a fast, economical particle-multiple-mesh N -body code optimized for large- N modelling of collisionless dynamical processes, such as black hole wandering or bar–halo interactions, occurring within isolated galaxies. The code has been specially designed to conserve linear momentum. Despite this, it also has variable softening and an efficient block-time-step scheme: the force between any pair of particles is calculated using the finest mesh that encloses them both (respecting Newton's third law) and is updated only on the longest time-step of the two (which conserves momentum). For realistic galaxy models with   N ≳ 10⁶  , it is faster than the fastest comparable tree code by factors ranging from ∼2 (using single time-steps) to ∼10 (multiple time-steps in a concentrated galaxy).