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.     

Dark matter halo response to the disc growth

We consider the sensitivity of the circular-orbit adiabatic contraction approximation to the baryon condensation rate and the orbital structure of dark matter haloes in the Λ cold dark matter (ΛCDM) paradigm. Using one-dimensional hydrodynamic simulations including the dark matter halo mass accretion history and gas cooling, we demonstrate that the adiabatic approximation is approximately valid even though haloes and discs may assemble simultaneously. We further demonstrate the validity of the simple approximation for ΛCDM haloes with isotropic velocity distributions using three-dimensional N -body simulations. This result is easily understood: an isotropic velocity distribution in a cuspy halo requires more circular orbits than radial orbits. Conversely, the approximation is poor in the extreme case of a radial orbit halo. It overestimates the response of a core dark matter halo, where radial orbit fraction is larger. Because no astronomically relevant models are dominated by low angular momentum orbits in the vicinity of the disc and the growth time-scale is never shorter than a dynamical time, we conclude that the adiabatic contraction approximation is useful in modelling the response of dark matter haloes to the growth of a disc.     

Galaxies in clusters: the observational characteristics of bow shocks, wakes and tails

The dynamical signatures of the interaction between galaxies in clusters and the intracluster medium (ICM) can potentially yield significant information about the structure and dynamical history of clusters. To develop our understanding of this phenomenon we present results from numerical modelling of the galaxy–ICM interaction, as the galaxy moves through the cluster. The simulations have been performed for a broad range of ICM temperatures ( kT _cl=1, 4 and 8 keV), representative of poor clusters or groups through to rich clusters.
There are several dynamical features that can be identified in these simulations. For supersonic galaxy motion, a leading bow shock is present, and also a weak gravitationally focused wake or tail behind the galaxy (analogous to Bondi–Hoyle accretion). For galaxies with higher mass replenishment rates and a denser interstellar medium (ISM), the dominant feature is a dense ram-pressure stripped tail. In line with other simulations, we find that the ICM/galaxy–ISM interaction can result in complex time-dependent dynamics, with ram-pressure stripping occurring in an episodic manner.
In order to facilitate this comparison between the observational consequences of dynamical studies and X-ray observations we have calculated synthetic X-ray flux and hardness maps from these simulations. These calculations predict that the ram-pressure stripped tail will usually be the most visible feature, though in nearby galaxies the bow shock preceding the galaxy should also be apparent in deeper X-ray observations. We briefly discuss these results and compare them with X-ray observations of galaxies where there is evidence of such interactions.     

Dual black holes in merger remnants – I. Linking accretion to dynamics

We study the orbital evolution and accretion history of massive black hole (MBH) pairs in rotationally supported circumnuclear discs up to the point where MBHs form binary systems. Our simulations have high resolution in mass and space which, for the first time, makes it feasible to follow the orbital decay of a MBH either counter- or corotating with respect to the circumnuclear disc. We show that a moving MBH on an initially counter-rotating orbit experiences an 'orbital angular momentum flip' due to the gas-dynamical friction, i.e. it starts to corotate with the disc before a MBH binary forms. We stress that this effect can only be captured in very high resolution simulations. Given the extremely large number of gas particles used, the dynamical range is sufficiently large to resolve the Bondi–Hoyle–Lyttleton radii of individual MBHs. As a consequence, we are able to link the accretion processes to the orbital evolution of the MBH pairs. We predict that the accretion rate is significantly suppressed and extremely variable when the MBH is moving on a retrograde orbit. It is only after the orbital angular momentum flip has taken place that the secondary rapidly 'lights up' at which point both MBHs can accrete near the Eddington rate for a few Myr. The separation of the double nucleus is expected to be around ≲10 pc at this stage. We show that the accretion rate can be highly variable also when the MBH is corotating with the disc (albeit to a lesser extent) provided that its orbit is eccentric. Our results have significant consequences for the expected number of observable double active galactic nuclei at separations of ≲100 pc.     

The tidal stripping of satellites

We present an improved analytic calculation for the tidal radius of satellites and test our results against N -body simulations.
The tidal radius in general depends upon four factors: the potential of the host galaxy, the potential of the satellite, the orbit of the satellite and the orbit of the star within the satellite . We demonstrate that this last point is critical and suggest using three tidal radii to cover the range of orbits of stars within the satellite. In this way we show explicitly that prograde star orbits will be more easily stripped than radial orbits; while radial orbits are more easily stripped than retrograde ones. This result has previously been established by several authors numerically, but can now be understood analytically. For point mass, power-law (which includes the isothermal sphere), and a restricted class of split power-law potentials our solution is fully analytic. For more general potentials, we provide an equation which may be rapidly solved numerically.
Over short times (≲1–2 Gyr ∼1 satellite orbit), we find excellent agreement between our analytic and numerical models. Over longer times, star orbits within the satellite are transformed by the tidal field of the host galaxy. In a Hubble time, this causes a convergence of the three limiting tidal radii towards the prograde stripping radius. Beyond the prograde stripping radius, the velocity dispersion will be tangentially anisotropic.     

Variable accretion and emission from the stellar winds in the Galactic Centre

We present numerical simulations of stellar wind dynamics in the central parsec of the Galactic Centre, studying in particular the accretion of gas on to Sgr A*, the supermassive black hole. Unlike our previous work, here we use state-of-the-art observational data on orbits and wind properties of individual wind-producing stars. Since wind velocities were revised upwards and non-zero eccentricities were considered, our new simulations show fewer clumps of cold gas and no conspicuous disc-like structure. The accretion rate is dominated by a few close 'slow-wind stars' ( v _w≤ 750 km s⁻¹), and is consistent with the Bondi estimate, but variable on time-scales of tens to hundreds of years. This variability is due to the stochastic infall of cold clumps of gas, as in earlier simulations, and to the eccentric orbits of stars. The present models fail to explain the high luminosity of Sgr A* a few hundred years ago implied by Integral observations, but we argue that the accretion of a cold clump with a small impact parameter could have caused it. Finally, we show the possibility of constraining the total mass-loss rate of the 'slow-wind stars' using near infrared observations of gas in the central few arcseconds.     

A hybrid N-body code incorporating algorithmic regularization and post-Newtonian forces

We describe a novel N -body code designed for simulations of the central regions of galaxies containing massive black holes. The code incorporates Mikkola's 'algorithmic' chain regularization scheme including post-Newtonian terms up to PN2.5 order. Stars moving beyond the chain are advanced using a fourth-order integrator with forces computed on a GRAPE board. Performance tests confirm that the hybrid code achieves better energy conservation, in less elapsed time, than the standard scheme and that it reproduces the orbits of stars tightly bound to the black hole with high precision. The hybrid code is applied to two sample problems: the effect of finite- N gravitational fluctuations on the orbits of the S-stars, and inspiral of an intermediate-mass black hole into the Galactic Centre.     

Mass distribution in nearby Abell clusters

We study the mass distribution in six nearby  ( z < 0.06)  relaxed Abell clusters of galaxies A0262, A0496, A1060, A2199, A3158 and A3558. Given the dominance of dark matter in galaxy clusters, we approximate their total density distribution by the Navarro, Frenk & White (NFW) formula characterized by virial mass and concentration. We also assume that the anisotropy of galactic orbits is reasonably well described by a constant and that galaxy distribution traces that of the total density. Using the velocity and position data for 120–420 galaxies per cluster we calculate, after removal of interlopers, the profiles of the lowest order even velocity moments, dispersion and kurtosis. We then reproduce the velocity moments by jointly fitting the moments to the solutions of the Jeans equations. Including the kurtosis in the analysis allows us to break the degeneracy between the mass distribution and anisotropy and constrain the anisotropy as well as the virial mass and concentration. The method is tested in detail on mock data extracted from the N -body simulations of dark matter haloes. We find that the best-fitting Galactic orbits are remarkably close to isotropic in most clusters. Using the fitted pairs of mass and concentration parameters for the six clusters, we conclude that the trend of decreasing concentration for higher masses found in the cosmological N -body simulations is consistent with the data. By scaling the individual cluster data by mass, we combine them to create a composite cluster with 1465 galaxies and perform a similar analysis on such sample. The estimated concentration parameter then lies in the range  1.5 < c < 14  and the anisotropy parameter in the range  −1.1 < β < 0.5  at the 95 per cent confidence level.     

ω Cen – An ultra compact dwarf galaxy?

We study the merging of star clusters out of cluster aggregates similar to Knot S in the Antennae on orbits close to the one of ω Cen by carrying out high resolution numerical N-body simulations. We want to constrain the parameter space which is able to produce merger objects with similar properties as ω Cen. This revised version was published online in August 2006 with corrections to the Cover Date.     

Dynamical friction of bodies orbiting in a gaseous sphere

The dynamical friction experienced by a body moving in a gaseous medium is different from the friction in the case of a collisionless stellar system. Here we consider the orbital evolution of a gravitational perturber inside a gaseous sphere using three-dimensional simulations, ignoring however self-gravity. The results are analysed in terms of a 'local' formula with the associated Coulomb logarithm taken as a free parameter. For forced circular orbits, the asymptotic value of the component of the drag force in the direction of the velocity is a slowly varying function of the Mach number in the range 1.0–1.6. The dynamical friction time-scale for free decay orbits is typically only half as long as in the case of a collisionless background, which is in agreement with E. C. Ostriker's recent analytic result. The orbital decay rate is rather insensitive to the past history of the perturber. It is shown that, similarly to the case of stellar systems, orbits are not subject to any significant circularization. However, the dynamical friction time-scales are found to increase with increasing orbital eccentricity for the Plummer model, whilst no strong dependence on the initial eccentricity is found for the isothermal sphere.     

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 .     

Resolving the timing problem of the globular clusters orbiting the Fornax dwarf galaxy

We re-investigate the old problem of the survival of the five globular clusters (GCs) orbiting the Fornax dwarf galaxy in both standard and modified Newtonian dynamics (MOND). For the first time in the history of the topic, we use accurate mass models for the Fornax dwarf, obtained through Jeans modelling of the recently published line-of-sight (LOS) velocity dispersion data, and we are also not resigned to circular orbits for the GCs. Previously conceived problems stem from fixing the starting distances of the globulars to be less than half the tidal radius. We relax this constraint since there is absolutely no evidence for it and show that the dark matter (DM) paradigm, with either cusped or cored DM profiles, has no trouble sustaining the orbits of the two least massive GCs for a Hubble time almost regardless of their initial distance from Fornax. The three most massive globulars can remain in orbit as long as their starting distances are marginally outside the tidal radius. The outlook for MOND is also not nearly as bleak as previously reported. Although dynamical friction (DF) inside the tidal radius is far stronger in MOND, outside DF is negligible due to the absence of stars. This allows highly radial orbits to survive, but more importantly circular orbits at distances more than 85 per cent of Fornax's tidal radius to survive indefinitely. The probability of the GCs being on circular orbits at this distance compared with their current projected distances is discussed and shown to be plausible. Finally, if we ignore the presence of the most massive globular (giving it a large LOS distance), we demonstrate that the remaining four globulars can survive within the tidal radius for the Hubble time with perfectly sensible orbits.     

Gas dynamic stripping and X-ray emission of cluster elliptical galaxies

Detailed three-dimensional numerical simulations of an elliptical galaxy orbiting in a gas-rich cluster of galaxies indicate that gas dynamic stripping is less efficient than the results from previous, simpler calculations by Takeda et al. and Gaetz et al. implied. This result is consistent with X-ray data for cluster elliptical galaxies. Hydrodynamic torques and direct accretion of orbital angular momentum can result in the formation of a cold gaseous disc, even in a non-rotating galaxy. The gas lost by cluster galaxies via the process of gas dynamic stripping tends to produce a colder, chemically enriched cluster gas core. A comparison of the models with the available X-ray data of cluster galaxies shows that the X-ray luminosity distribution of cluster galaxies may reflect hydrodynamic stripping, but also that a purely hydrodynamic treatment is inadequate for the cooler interstellar medium near the centre of the galaxy.     

The globular cluster system of the young elliptical NGC 6702

We study the globular cluster (GC) system of the dust-lane elliptical galaxy NGC 6702, using B -, V - and I -band imaging observations carried out at the Keck telescope. This galaxy has a spectroscopic age of ≈2 Gyr suggesting recent star formation. We find strong evidence for a bimodal GC colour distribution, with the blue peak having a colour similar to that of the Galactic halo GCs. Assuming that the blue GCs are indeed old and metal-poor, we estimate an age of 2–5 Gyr and supersolar metallicity for the red GC subpopulation. Despite the large uncertainties, this is in reasonable agreement with the spectroscopic galaxy age. Additionally, we estimate a specific frequency of S _N =2.3±1.1 for NGC 6702. We predict that passive evolution of NGC 6702 will further increase its specific frequency to S _N ≈2.7 within 10 Gyr, in closer agreement to that of typical present-day ellipticals. We also discuss evidence that the merger/accretion event that took place a few Gyr ago involved a high gas fraction.     

The formation and evolution of clusters of galaxies in different cosmogonies

The formation of galaxy clusters in hierarchically clustering universes is investigated by means of high-resolution N -body simulations. The simulations are performed using a newly developed multimass scheme which combines a PM code with a high-resolution N -body code. Numerical effects resulting from time-stepping and gravitational softening are investigated, as well as the influence of the simulation box size and of the assumed boundary conditions. Special emphasis is laid on the formation process and the influence of various cosmological parameters. Cosmogonies with massive neutrinos are also considered. Differences between clusters in the same cosmological model seem to dominate over differences caused by differing background cosmogony. The cosmological model can alter the time evolution of cluster collapse, but the merging pattern remains fairly similar, e.g. the number of mergers and the mass ratio of mergers. The gross properties of a halo, such as its size and total angular momentum, also evolve in a similar manner for all cosmogonies, and can be described using analytical models. It is shown that the density distribution of a halo shows a characteristic radial dependence which follows a power law with a slope of =1 at small radii and =3 at large radii, independent of the background cosmogony or the considered redshift. The shape of the density profiles follows the generic form proposed by Navarro et al. for all hierarchically clustering scenarios, and retains very little information about the formation process or the cosmological model. Only the central matter concentration of a halo is correlated with the formation time and therefore the corresponding cosmogony. We emphasize the role of non-radial motions of the halo particles in the evolution of the density profile.     

Ram pressure stripping the hot gaseous haloes of galaxies in groups and clusters

We use a large suite of carefully controlled full hydrodynamic simulations to study the ram pressure stripping of the hot gaseous haloes of galaxies as they fall into massive groups and clusters. The sensitivity of the results to the orbit, total galaxy mass, and galaxy structural properties is explored. For typical structural and orbital parameters, we find that ∼30 per cent of the initial hot galactic halo gas can remain in place after 10 Gyr. We propose a physically simple analytic model that describes the stripping seen in the simulations remarkably well. The model is analogous to the original formulation of Gunn & Gott, except that it is appropriate for the case of a spherical (hot) gas distribution (as opposed to a face-on cold disc) and takes into account that stripping is not instantaneous but occurs on a characteristic time-scale. The model reproduces the results of the simulations to within ≈10 per cent at almost all times for all the orbits, mass ratios, and galaxy structural properties we have explored. The one exception involves unlikely systems where the orbit of the galaxy is highly non-radial and its mass exceeds about 10 per cent of the group or cluster into which it is falling (in which case the model underpredicts the stripping following pericentric passage). The proposed model has several interesting applications, including modelling the ram pressure stripping of both observed and cosmologically simulated galaxies and as a way to improve present semi-analytic models of galaxy formation. One immediate consequence is that the colours and morphologies of satellite galaxies in groups and clusters will differ significantly from those predicted with the standard assumption of complete stripping of the hot coronae.     

An MHD gadget for cosmological simulations

Various radio observations have shown that the hot atmospheres of galaxy clusters are magnetized. However, our understanding of the origin of these magnetic fields, their implications on structure formation and their interplay with the dynamics of the cluster atmosphere, especially in the centres of galaxy clusters, is still very limited. In preparation for the upcoming new generation of radio telescopes (like Expanded Very Large Array, Low Wavelength Array, Low Frequency Array and Square Kilometer Array), a huge effort is being made to learn more about cosmological magnetic fields from the observational perspective. Here we present the implementation of magnetohydrodynamics (MHD) in the cosmological smoothed particle hydrodynamics (SPH) code gadget . We discuss the details of the implementation and various schemes to suppress numerical instabilities as well as regularization schemes, in the context of cosmological simulations. The performance of the SPH–MHD code is demonstrated in various one- and two-dimensional test problems, which we performed with a fully, three-dimensional set-up to test the code under realistic circumstances. Comparing solutions obtained using athena , we find excellent agreement with our SPH–MHD implementation. Finally, we apply our SPH–MHD implementation to galaxy cluster formation within a large, cosmological box. Performing a resolution study we demonstrate the robustness of the predicted shape of the magnetic field profiles in galaxy clusters, which is in good agreement with previous studies.     

Evaluating approximations for halo merging histories

We study the merging history of dark matter haloes in N -body simulations and semi-analytical 'merger trees' based on the extended Press–Schechter (EPS) formalism. The main focus of our study is the joint distribution of progenitor number and mass as a function of redshift and parent halo mass. We begin by investigating the mean quantities predicted directly by the Press–Schechter (PS) and EPS formalism, such as the halo mass and conditional mass functions, and compare these predictions with the results of the simulations. The higher moments of this distribution are not predicted by the EPS formalism alone and must be obtained from the merger trees. We find that the Press–Schechter model deviates from the simulations at the level of 30–50 per cent on certain mass scales, and that the sense of the discrepancy changes as a function of redshift. We show that this discrepancy is reflected in the higher moments of the distribution of progenitor mass and number. We investigate some related statistics such as the accretion rate and the mass ratio of the largest two progenitors. For galaxy sized haloes ( M ∼10¹² M_⊙), we find that the merging history of haloes, as represented by these statistics, is well reproduced in the merger trees compared with the simulations. The agreement deteriorates for larger mass haloes. We conclude that merger trees based on the extended Press–Schechter formalism provide a reasonably reliable framework for semi-analytical models of galaxy formation.     

The dynamical evolution of substructure

The evolution of substructure embedded in non-dissipative dark haloes is studied through N -body simulations of isolated systems, both in and out of initial equilibrium, complementing cosmological simulations of the growth of structure. We determine by both analytic calculations and direct analysis of the N -body simulations the relative importance of various dynamical processes acting on the clumps, such as the removal of material by global tides, clump–clump heating, clump–clump merging and dynamical friction. The ratio of the internal clump velocity dispersion to that of the dark halo is an important parameter; as this ratio approaches a value of unity, heating by close encounters between clumps becomes less important, while the other dynamical processes continue to increase in importance. Our comparison between merging and disruption processes implies that spiral galaxies cannot be formed in a protosystem that contains a few large clumps, but can be formed through the accretion of many small clumps; elliptical galaxies form in a more clumpy environment than do spiral galaxies. Our results support the idea that the central cusp in the density profiles of dark haloes is the consequence of self-limiting merging of small, dense haloes. This implies that the collapse of a system of clumps/substructure is not sufficient to form a cD galaxy, with an extended envelope; plausibly, subsequent accretion of large galaxies is required. The post-collapse system is in general triaxial, with rounder systems resulting from fewer, but more massive, clumps. Persistent streams of material from disrupted clumps can be found in the outer regions of the final system, and at an overdensity of around 0.75, can cover 10 to 30 per cent of the sky.