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 .     

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.     

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.     

The Adaptive TreePM: an adaptive resolution code for cosmological N-body simulations

Cosmological N -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 20 yr the limitations imposed by computing power forced simulators to ignore some of the basic requirements for modelling 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 (ATreePM), 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.     

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.     

Numerical simulations of the Lyman α forest – a comparison of gadget-2 and enzo

We compare simulations of the Lyman α forest performed with two different hydrodynamical codes, gadget-2 and enzo . A comparison of the dark matter power spectrum for simulations run with identical initial conditions show differences of 1–3 per cent at the scales relevant for quantitative studies of the Lyman α forest. This allows a meaningful comparison of the effect of the different implementations of the hydrodynamic part of the two codes. Using the same cooling and heating algorithm in both codes, the differences in the temperature and the density probability distribution function are of the order of 10 per cent. The differences are comparable to the effects of box size and resolution on these statistics. When self-converged results for each code are taken into account, the differences in the flux power spectrum – the statistics most widely used for estimating the matter power spectrum and cosmological parameters from Lyman α forest data – are about 5 per cent. This is again comparable to the effects of box size and resolution. Numerical uncertainties due to a particular implementation of solving the hydrodynamic or gravitational equations appear therefore to contribute only moderately to the error budget in estimates of the flux power spectrum from numerical simulations. We further find that the differences in the flux power spectrum for enzo simulations run with and without adaptive mesh refinement are also of the order of 5 per cent or smaller. The latter require 10 times less CPU time making the CPU time requirement similar to that of a version of gadget-2 that is optimized for Lyman α forest simulations.     

Realistic triaxial density–potential–force profiles for stellar systems and dark matter haloes

Popular models for describing the luminosity–density profiles of dynamically hot stellar systems (e.g. Jaffe, Hernquist, Dehnen) were constructed to match the deprojected form of de Vaucouleurs' R ^1/4 light-profile. However, we now know that elliptical galaxies and bulges display a mass-dependent range of structural profiles. To compensate this, the model in Terzić & Graham was designed to closely match the deprojected form of Sérsic R ^{1/ n} light-profiles, including deprojected exponential light-profiles and galaxies with partially depleted cores. It is thus applicable for describing bulges in spiral galaxies, dwarf elliptical galaxies, both 'power-law' and 'core' elliptical galaxies, also dark matter haloes formed from Λ cold dark matter cosmological simulations. In this paper, we present a new family of triaxial density–potential–force triplets, which generalizes the spherical model reported in Terzić & Graham to three dimensions. If the (optional) power-law core is present, it is a five-parameter family, while in the absence of the core it reduces to three parameters. The isodensity contours in the new family are stratified on confocal ellipsoids and the potential and forces are expressed in terms of integrals which are easy to evaluate numerically. We provide the community with a suite of numerical routines for orbit integration, which feature: optimized computations of potential and forces for this family; the ability to run simulations on parallel platforms; and modular and easily editable design.     

Chemical evolution of galaxies – I. A composition-dependent SPH model for chemical evolution and cooling

We describe an smooth particle hydrodynamics (SPH) model for chemical enrichment and radiative cooling in cosmological simulations of structure formation. This model includes: (i) the delayed gas restitution from stars by means of a probabilistic approach designed to reduce the statistical noise and, hence, to allow for the study of the inner chemical structure of objects with moderately high numbers of particles; (ii) the full dependence of metal production on the detailed chemical composition of stellar particles by using, for the first time in SPH codes, the   Q _ij   matrix formalism that relates each nucleosynthetic product to its sources and (iii) the full dependence of radiative cooling on the detailed chemical composition of gas particles, achieved through a fast algorithm using a new metallicity parameter ζ( T ) that gives the weight of each element on the total cooling function. The resolution effects and the results obtained from this SPH chemical model have been tested by comparing its predictions in different problems with known theoretical solutions. We also present some preliminary results on the chemical properties of elliptical galaxies found in self-consistent cosmological simulations. Such simulations show that the above ζ-cooling method is important to prevent an overestimation of the metallicity-dependent cooling rate, whereas the   Q _ij   formalism is important to prevent a significant underestimation of the [α/Fe] ratio in simulated galaxy-like objects.     

Simulating galactic outflows with kinetic supernova feedback

Feedback from star formation is thought to play a key role in the formation and evolution of galaxies, but its implementation in cosmological simulations is currently hampered by a lack of numerical resolution. We present and test a subgrid recipe to model feedback from massive stars in cosmological smoothed particle hydrodynamics simulations. The energy is distributed in kinetic form among the gas particles surrounding recently formed stars. The impact of the feedback is studied using a suite of high-resolution simulations of isolated disc galaxies embedded in dark haloes with total mass 10¹⁰ and  10¹²  h ⁻¹ M_⊙  . We focus, in particular, on the effect of pressure forces on wind particles within the disc, which we turn off temporarily in some of our runs to mimic a recipe that has been widely used in the literature. We find that this popular recipe gives dramatically different results because (ram) pressure forces on expanding superbubbles determine both the structure of the disc and the development of large-scale outflows. Pressure forces exerted by expanding superbubbles puff up the disc, giving the dwarf galaxy an irregular morphology and creating a galactic fountain in the massive galaxy. Hydrodynamic drag within the disc results in a strong increase in the effective mass loading of the wind for the dwarf galaxy, but quenches much of the outflow in the case of the high-mass galaxy.     

Sunyaev–Zel'dovich cluster survey simulations for Planck

We examine the ability of the future Planck mission to provide a catalogue of galaxy clusters observed via their Sunyaev–Zel'dovich (SZ) distortion in the cosmic microwave background (CMB). For this purpose we produce full-sky SZ maps based on N -body simulations and scaling relations between cluster properties for several cosmological models. We extrapolate the N -body simulations by a mass function to high redshifts in order to obtain a realistic SZ background. The simulated Planck observations include, besides the thermal and kinematic SZ effects, contributions from the primordial CMB, extragalactic point sources as well as Galactic dust, free–free and synchrotron emission. A harmonic-space maximum-entropy method is used to separate the SZ signal from contaminating components in combination with a cluster detection algorithm based on thresholding and flux integration to identify clusters and to obtain their fluxes. We estimate a survey sensitivity limit (depending on the quality of the recovered cluster flux) and provide cluster survey completeness and purity estimates. We find that, given our modelling and detection algorithm, Planck will reliably detect at least several thousands of clusters over the full sky. The exact number depends on the particular cosmological model (up to 10 000 cluster detections in a concordance ΛCDM model with  σ₈= 0.9  ). We show that the Galaxy does not significantly affect the cluster detection. Furthermore, the dependence of the thermal SZ power spectrum on the matter variance on scales of  8 h ⁻¹  Mpc and the quality of its reconstruction by the employed method are investigated. Our simulations suggest that the Planck cluster sample will not only be useful as a basis for follow-up observations, but also will have the ability to provide constraints on cosmological parameters.     

In-shock cooling in numerical simulations

We model a one-dimensional shock-tube using smoothed particle hydrodynamics and investigate the consequences of having finite shock-width in numerical simulations caused by finite resolution of the codes. We investigate the cooling of gas during passage through the shock for three different cooling regimes.
For a theoretical shock temperature of 10⁵ K, the maximum temperature of the gas is much reduced. When the ratio of the cooling time to shock-crossing time was 8, we found a reduction of 25 per cent in the maximum temperature reached by the gas. When the ratio was reduced to 1.2, the maximum temperature reached dropped to 50 per cent of the theoretical value. In both cases the cooling time was reduced by a factor of 2.
At lower temperatures, we are especially interested in the production of molecular hydrogen, and so we follow the ionization level and H₂ abundance across the shock. The effect of in-shock cooling is substantial: the maximum temperature the gas reaches compared with the theoretical temperature is found to vary between 0.15 and 0.81, depending upon the shock strength and mass resolution. The downstream ionization level is reduced from the theoretical level by a factor of between 2.4 and 12.5, and the resulting H₂ abundance by a factor of 1.35 to 2.22.
At temperatures above 10⁵ K, radiative shocks are unstable and will oscillate. We find that the shock jump temperature varies by a factor of 20 because of these oscillations.
We conclude that extreme caution must be exercised when interpreting the results of simulations of galaxy formation.     

The alignment of clusters using large-scale simulations

The alignment of clusters of galaxies with their nearest neighbours and between clusters within a supercluster is investigated using simulations of 512³ dark matter particles for ΛCDM and τ CDM cosmological models. Strongly significant alignments are found for separations of up to 15  h ⁻¹ Mpc in both cosmologies, but for the ΛCDM model the alignments extend up to separations of 30  h ⁻¹ Mpc. The effect is strongest for nearest neighbours, but is not significant enough to be useful as an observational discriminant between cosmologies. As a check of whether this difference in alignments is present in other cosmologies, smaller simulations with 256³ particles are investigated for four different cosmological models. Because of poor number statistics, only the standard CDM model shows indications of having different alignments from the other models.     

Lyman break galaxies at z= 4–6 in cosmological smoothed particle hydrodynamics simulations

We perform a spectrophotometric analysis of galaxies at redshifts z = 4–6 in cosmological smoothed particle hydrodynamics simulations of a Λ cold dark matter universe. Our models include radiative cooling and heating by a uniform ultraviolet (UV) background, star formation, supernova feedback, and a phenomenological model for galactic winds. Analysing a series of simulations of varying box size and particle number allows us to isolate the impact of numerical resolution on our results. Specifically, we determine the luminosity functions in B , V , R , i ' and z ' filters, and compare the results with observational surveys of Lyman break galaxies (LBGs) performed with the Subaru telescope and the Hubble Space Telescope . We find that the simulated galaxies have UV colours consistent with observations and fall in the expected region of the colour–colour diagrams used by the Subaru group. The stellar masses of the most massive galaxies in our largest simulation increase their stellar mass from   M _★∼ 10¹¹ M_⊙  at z = 6 to   M _★∼ 10^11.7 M_⊙  at z = 3. Assuming a uniform extinction of E ( B − V ) = 0.15, we also find reasonable agreement between simulations and observations in the space density of UV bright galaxies at z = 3–6, down to the magnitude limit of each survey. For the same moderate extinction level of E ( B − V ) ∼ 0.15, the simulated luminosity functions match observational data, but have a steep faint-end slope with α∼−2.0. We discuss the implications of the steep faint-end slope found in the simulations. Our results confirm the generic conclusion from earlier numerical studies that UV bright LBGs at z ≥ 3 are the most massive galaxies with E ( B − V ) ∼ 0.15 at each epoch.     

Constraining quasar host halo masses with the strength of nearby Lyα forest absorption

Using cosmological hydrodynamic simulations, we measure the mean transmitted flux in the Lyα forest for quasar sightlines that pass near a foreground quasar. We find that the trend of absorption with pixel quasar separation distance can be fitted using a simple power-law form including the usual correlation function parameters r ₀ and γ, so that     . From the simulations, we find the relation between r ₀ and quasar host mass, and formulate this as a way to estimate quasar host dark matter halo masses, quantifying uncertainties due to cosmological and IGM parameters, and redshift errors. With this method, we examine data for ∼9000 quasars from the Sloan Digital Sky Survey (SDSS) Data Release 5, assuming that the effect of ionizing radiation from quasars (the so-called transverse proximity effect) is unimportant (no evidence for it is seen in the data). We find that the best-fitting host halo mass for SDSS quasars with mean redshift z = 3 and absolute G -band magnitude −27.5 is  log  M /M_⊙= 12.68^+0.81_−0.67  . We also use the Lyman-Break Galaxy (LBG) and Lyα forest data of Adelberger et al. in a similar fashion to constrain the halo mass of LBGs to be  log₁₀  M /M_⊙= 11.41^+0.54_−0.59  , a factor of ∼20 lower than the bright quasars. In addition, we study the redshift distortions of the Lyα forest around quasars, using the simulations. We use the quadrupole to monopole ratio of the quasar Lyα forest correlation function as a measure of the squashing effect. We find its dependence on halo mass difficult to measure, but find that it may be useful for constraining cosmic geometry.     

Parameter tests within cosmological simulations of galaxy formation

Numerical simulations of galaxy formation require a number of parameters. Some of these are intrinsic to the numerical integration scheme (e.g., the time-step), while others describe the physical model (e.g., the gas metallicity). In this paper we present results of a systematic exploration of the effects of varying a subset of these parameters on simulations of galaxy formation. We use N -body and 'Smoothed Particle Hydrodynamics' techniques to follow the evolution of cold dark matter and gas in a small volume. We compare a fiducial model with 24 different simulations, in which one parameter at a time is varied, focusing on properties such as the relative fraction of hot and cold gas, and the abundance and masses of galaxies. We find that for reasonable choices of numerical values, many parameters have relatively little effect on the galaxies, with the notable exception of the parameters that control the resolution of the simulation and the efficiency with which gas cools.     

The clustering and host haloes of galaxy mergers at high redshift

High-resolution simulations of cosmological structure formation indicate that dark matter substructure in dense environments, such as groups and clusters, may survive for a long time. These dark matter subhaloes are the likely hosts of galaxies. We examine the small-scale spatial clustering of subhalo major mergers at high redshift using high-resolution N -body simulations of cosmological volumes. Recently merged, massive subhaloes exhibit enhanced clustering on scales  ∼100–300  h ⁻¹ kpc  , relative to all subhaloes of the same infall mass, for a short time after a major merger (<500 Myr). The small-scale clustering enhancement is smaller for lower mass subhaloes, which also show a deficit on scales just beyond the excess. Haloes hosting recent subhalo mergers tend to have more subhaloes; for massive subhaloes, the excess is stronger and it tends to increase for the most massive host haloes. The subhalo merger fraction is independent of halo mass for the scales we probe. In terms of satellite and central subhaloes, the merger increase in small-scale clustering for massive subhaloes arises from recently merged massive central subhaloes having an enhanced satellite population. Our mergers are defined via their parent infall mass ratios. Subhaloes experiencing major mass gains also exhibit a small-scale clustering enhancement, but these correspond to two-body interactions leading to two final subhaloes, rather than subhalo coalescence.     

Investigating cosmological weak lensing with the ray-bundle method

Using the ray-bundle method for calculating gravitational lens magnifications, we outline a method by which the magnification probability may be determined specifically in the weak lensing limit for cosmological models obtained from N -body simulations.
16 different models are investigated, which are variations on three broad classes of cold dark matter model: the standard model with  (Ω₀, λ ₀)=(1.0,0.0)  , the open model with  (Ω₀, λ ₀)=(0.3,0.0)  and the lambda model, which is a flat model with a cosmological constant  (Ω₀, λ ₀)=(0.3,0.7)  .
The effects of varying the Hubble parameter, H ₀, the power spectrum shape parameter, Γ, and the cluster mass normalization, σ ₈, are studied. It is shown that there is no signature of these parameters in the weak lensing magnification distributions. The magnification probability distributions are also shown to be independent of the numerical parameters such as the lens mass and simulation box size in the N -body simulations.     

Constraints on cosmological and biasing models using active galactic nucleus clustering

We attempt to put constraints on different cosmological and biasing models by combining the recent clustering results of X-ray sources in the local ( z ≤0.1) and distant Universe ( z ∼1) . To this end we compare the measured angular correlation function for bright (Akylas et al.) and faint (Vikhlinin & Forman) ROSAT X-ray sources respectively with those expected in three spatially flat cosmological models. Taking into account the different functional forms of the bias evolution, we find that there are two cosmological models which match the data well. In particular, low-Ω_○ cosmological models (Ω_Λ=1−Ω_○=0.7) that contain either (i) high σ ₈^mass=1.13 value with galaxy merging bias, b ( z )∝(1+ z )^1.8 or (ii) low σ ₈^mass=0.9 with non-bias, b ( z ) ≡ 1 best reproduce the AGN clustering results, while τ CDM models with different bias behaviour are ruled out at a high significance level.     

Gravitational lensing effects on high-redshift Type II supernova studies with the Next Generation Space Telescope

We derive the expected Type II supernova (SN) differential number counts, N ( m ), and Hubble diagram for SCDM and LCDM cosmological models, taking into account the effects of gravitational lensing (GL) produced by the intervening cosmological mass. The mass distribution of dark matter haloes (i.e. the lenses) is obtained by means of a Monte Carlo method applied to the PressSchechter mass function. The haloes are assumed to have a Navarro, Frenk & White (NFW) density profile, in agreement with recent simulations of hierarchical cosmological models. Up to z =15, the (SCDM, LCDM) models predict a total number of (857, 3656) SNII yr¹ in 100 surveyed 44 arcmin² fields of the Next Generation Space Telescope ( NGST ). NGST will be able to reach the peak of the N ( m ) curve, located at AB 30(31) for SCDM (LCDM) in J and K wavelength bands, and detect (75 per cent, 51 per cent) of the above SN events. This will allow a detailed study of the early cosmic star formation history, as traced by SNIIe. N ( m ) is only very mildly affected by the inclusion of lensing effects. In addition, GL introduces a moderate uncertainty in the determination of cosmological parameters from Hubble diagrams, when these are pushed to higher z . For example, for a 'true' LCDM with (_M=0.4, =0.6), without proper account of GL, one would instead derive We briefly compare our results with previous similar work and discuss the limitations of the model.     

Collapse and fragmentation of rotating magnetized clouds – I. Magnetic flux–spin relation

We discuss the evolution of the magnetic flux density and angular velocity in a molecular cloud core, on the basis of three-dimensional numerical simulations, in which a rotating magnetized cloud fragments and collapses to form a very dense optically thick core of  >5 × 10¹⁰ cm⁻³  . As the density increases towards the formation of the optically thick core, the magnetic flux density and angular velocity converge towards a single relationship between the two quantities. If the core is magnetically dominated its magnetic flux density approaches  1.5( n /5 × 10¹⁰ cm⁻³)^1/2 mG  , while if the core is rotationally dominated the angular velocity approaches  2.57 × 10⁻³ ( n /5 × 10¹⁰ cm⁻³)^1/2 yr⁻¹  , where n is the density of the gas. We also find that the ratio of the angular velocity to the magnetic flux density remains nearly constant until the density exceeds  5 × 10¹⁰ cm⁻³  . Fragmentation of the very dense core and emergence of outflows from fragments will be shown in the subsequent paper.