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).     

A tree code for planetesimal dynamics: comparison with a hybrid direct code

We present a tree code for simulations of collisional systems dominated by a central mass. We describe the implementation of the code and the results of some test runs with which the performance of the code was tested. A comparison between the behaviour of the tree code and a direct hybrid integrator is also presented. The main result is that tree codes can be useful in numerical simulations of planetary accretion, especially during intermediate stages, where possible runaway accretion and dynamical friction lead to a population with a few large bodies in low-eccentricity and low-inclination orbits embedded in a large swarm of small planetesimals in rather excited orbits. Some strategies to improve the performance of the code are also discussed.     

A momentum-conserving N-body scheme with individual time steps

We present an N-body code called Taichi for galactic dynamics and controlled numerical experiments. The code includes two high-order hierarchical multipole expansion methods: the Barnes-Hut (BH) tree and the fast multipole method (FMM). For the time integration, the code can use either a conventional adaptive KDK or a Hamiltonian splitting integrator. The combination of FMM and the Hamiltonian splitting integrator leads to a momentum-conserving N-body scheme with individual time steps. We find Taichi performs well in the typical applications in galactic dynamics. In the isolated and interacting galaxies tests, the momentum conserving scheme produces the same result as a conventional BH tree code. But for similar force accuracies, FMM significantly speeds up the simulations compared to the monopole BH tree. In the cold collapse test, we find the inner structure after relaxation can be sensitive to the force accuracies. Taichi is ready to incorporate special treatment of close encounters thanks to the Hamiltonian splitting integrator, suitable for studying dynamics around central massive bodies.     

Galaxies in the connection machine

Many problems in galactic dynamics require computer speeds orders of magnitude larger than what can be achieved on current single-processor machines. In the near future such speeds are likely to become available through computer architectures based on large-scale, fine-grained parallelism. An example of a highly parallel computer is the Connection Machine, with up to 65,636 processors. We have benchmarked gravitationalN-body algorithms on the Connection Machine, and compared those with similar benchmarks which we have obtained on more traditional vector supercomputers. Our conclusions are: (1) The direct summation algorithm, with of orderN ² interactions forN particles, can be made to run with high efficiency on either type of computer. As a result, the Connection Machine clearly wins in speed over all supercomuters tested, with the exception of an 8-processor ETA, which shows a comparable performance. (2) A more efficient tree algorithm reduces the growth of the number of interactions fromN ² toN logN. However, the greater complexity of this algorithm causes a considerable degradation of efficiency, by a factor which is larger on the Connection Machine than on vector supercomputers. As a result, out tree code runs at comparable speeds on both types of machines, with the notable exception of the 8-processor ETA, which has an extrapolated speed for running our tree code which is higher than any of the other machines we have tested.     

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 .     

Computational advances in gravitational microlensing: A comparison of CPU,GPU, and parallel,large data codes

To assess how future progress in gravitational microlensing computation at high optical depth will rely on both hardware and software solutions, we compare a direct inverse ray-shooting code implemented on a graphics processing unit (GPU) with both a widely-used hierarchical tree code on a single-core CPU, and a recent implementation of a parallel tree code suitable for a CPU-based cluster supercomputer. We examine the accuracy of the tree codes through comparison with a direct code over a much wider range of parameter space than has been feasible before. We demonstrate that all three codes present comparable accuracy, and choice of approach depends on considerations relating to the scale and nature of the microlensing problem under investigation. On current hardware, there is little difference in the processing speed of the single-core CPU tree code and the GPU direct code, however the recent plateau in single-core CPU speeds means the existing tree code is no longer able to take advantage of Moore’s law-like increases in processing speed. Instead, we anticipate a rapid increase in GPU capabilities in the next few years, which is advantageous to the direct code. We suggest that progress in other areas of astrophysical computation may benefit from a transition to GPUs through the use of “brute force” algorithms, rather than attempting to port the current best solution directly to a GPU language – for certain classes of problems, the simple implementation on GPUs may already be no worse than an optimised single-core CPU version.     

Application of a massively parallel computer to the N-body problem

Parallel processor computers represent a new technology that has recently become available for astronomical applications. We have implemented an N-body code on a TMC Connection Machine CM-2 in order to investigate the advantages of a massively parallel computer over serial machines, including conventional supercomputers. For collisionless problems following N stars, a direct integration code scales as O(N²) on serial machines and on the CM-2 as O(log(N)) for small N and O(N log(N)) for large N. The CM-2 outperforms workstations for N>50 and supercomputers for N>4000.     

Multi-fluid model of comet 1P/Halley

A new three-dimensional magnetohydrodynamic model of the coma of a comet has been developed and applied to simulations of a Halley-class coma using the solar-wind conditions of the Giotto flyby of Halley in 1986. The code developed for high-performance parallel processing computers, combines the high spatial resolution of smaller than 1 km grid spacing near the nucleus, with a large computational domain that enables structures nearly 10 million km down the comet tail to be modeled. Ions, neutrals, and electrons are considered as separate interacting fluids. Significant physical processes treated by the model include both photo and electron impact ionization of neutrals, recombination of ions, charge exchange between solar-wind ions and cometary neutrals, and frictional interactions between the three fluids considered in the model. A variety of plasma structures and physical parameters that are the output of this model are compared with relevant Giotto data from the 1986 Halley flyby.     

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.     

Monte Carlo simulations of star clusters — I. First Results

A revision of Stodoíkiewicz's Monte Carlo code is used to simulate evolution of star clusters. The new method treats each superstar as a single star and follows the evolution and motion of all individual stellar objects. The first calculations for isolated, equal-mass N -body systems with three-body energy generation according to Spitzer's formulae show good agreement with direct N -body calculations for N  = 2000, 4096 and 10 000 particles. The density, velocity, mass distributions, energy generation, number of binaries, etc., follow the N -body results. Only the number of escapers is slightly too high compared with N -body results, and there is no level-off anisotropy for advanced post-collapse evolution of Monte Carlo models as is seen in N -body simulations for N  ≤ 2000. For simulations with N  > 10 000 gravothermal oscillations are clearly visible. The calculations of N   2000, 4096, 10 000, 32 000 and 100 000 models take about 2, 6, 20, 130 and 2500 h, respectively. The Monte Carlo code is at least 10⁵ times faster than the N -body one for N  = 32 768 with special-purpose hardware. Thus it becomes possible to run several different models to improve statistical quality of the data and run individual models with N as large as 100 000. The Monte Carlo scheme can be regarded as a method which lies in the middle between direct N -body and Fokker–Planck models and combines most advantages of both methods.     

Enigmatic aspects of entropy inside the black hole: what do falling comoving observers see?

In a recent paper it was suggested that inclusion of mutual gravitational interactions can give a possible scenario for reversing gravitation collapse and averting a singular phase. We extend this idea to the still unsolved problem of matter collapsing beyond black hole event horizons. For a comoving observer there is no change in entropy as he goes through the horizon. Matter collapses to a minimum radius, and then can re-expand with the same entropy. It is shown that phase space inside a collapsing black hole is also invariant.     

Numerical simulations of protostellar encounters — III. Non-coplanar disc–disc encounters

It is expected that an average protostar will undergo at least one impulsive interaction with a neighbouring protostar whilst a large fraction of its mass is still in a massive, extended disc. If protostars are formed individually within a cluster before falling together and interacting, there should be no preferred orientation for such interactions. As star formation within clusters is believed to be coeval, it is probable that, during interactions, both protostars possess massive, extended discs.   We have used an SPH code to carry out a series of simulations of non-coplanar disc–disc interactions. We find that non-coplanar interactions trigger gravitational instabilities in the discs, which may then fragment to form new companions to the existing stars. (This is different from coplanar interactions, in which most of the new companion stars form after material in the discs has been swept up into a shock layer, and this then fragments.) The original stars may also capture each other, leading to the formation of a small- N cluster. If every star undergoes a randomly oriented disc–disc interaction, then the outcome will be the birth of many new stars and substellar objects. Approximately two-thirds of the stars will end up in multiple systems.     

Simulations of spheroidal systems with substructure: trees in fields

We present a hybrid technique of N -body simulation to deal with collisionless stellar systems having an inhomogeneous global structure. We combine a treecode and a self-consistent field code such that each of the codes models a different component of the system being investigated. The treecode is suited to treatment of dynamically cold or clumpy components, which may undergo significant evolution within a dynamically hot system. The hot system is appropriately evolved by the self-consistent field code. This combined code is particularly suited to a number of problems in galactic dynamics. Applications of the code to these problems are briefly discussed.     

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.     

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.     

Regular and chaotic motion in softened gravitational systems

The stability of the dynamical trajectories of softened spherical gravitational systems is examined, both in the case of the full N -body problem and that of trajectories moving in the gravitational field of non-interacting background particles. In the latter case, for   N 10 000  , some trajectories, even if unstable, had exceedingly long diffusion times, which correlated with the characteristic e-folding time-scale of the instability. For trajectories of   N ≈100 000  systems this time-scale could be arbitrarily large – and thus appear to correspond to regular orbits. For centrally concentrated systems, low angular momentum trajectories were found to be systematically more unstable. This phenomenon is analogous to the well-known case of trajectories in generic centrally concentrated non-spherical smooth systems, where eccentric trajectories are found to be chaotic. The exponentiation times also correlate with the conservation of the angular momenta along the trajectories. For N up to a few hundred, the instability time-scales of N -body systems and their variation with particle number are similar to those of the most chaotic trajectories in inhomogeneous non-interacting systems. For larger N (up to a few thousand) the values of the these time-scales were found to saturate, increasing significantly more slowly with N . We attribute this to collective effects in the fully self-gravitating problem, which are apparent in the time variations of the time-dependent Liapunov exponents. The results presented here go some way towards resolving the long-standing apparent paradoxes concerning the local instability of trajectories. This now appears to be a manifestation of mechanisms driving evolution in gravitational systems and their interactions – and may thus be a useful diagnostic of such processes.     

A GPU accelerated Barnes–Hut tree code for FLASH4

We present a GPU accelerated CUDA-C implementation of the Barnes Hut (BH) tree code for calculating the gravitational potential on octree adaptive meshes. The tree code algorithm is implemented within the FLASH4 adaptive mesh refinement (AMR) code framework and therefore fully MPI parallel. We describe the algorithm and present test results that demonstrate its accuracy and performance in comparison to the algorithms available in the current FLASH4 version. We use a MacLaurin spheroid to test the accuracy of our new implementation and use spherical, collapsing cloud cores with effective AMR to carry out performance tests also in comparison with previous gravity solvers. Depending on the setup and the GPU/CPU ratio, we find a speedup for the gravity unit of at least a factor of 3 and up to 60 in comparison to the gravity solvers implemented in the FLASH4 code. We find an overall speedup factor for full simulations of at least factor 1.6 up to a factor of 10.     

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.     

Problems on the Galilean satellites of Jupiter

The purpose of this paper is to present a critical review of some problems concerning the dynamics of Jupiter's Galilean system of satellites. Theory, ephemeris and observation are considered.Two theories were proposed by Ferraz-Mello and by Sagnier. The main characteristics of these theories are that the frequencies are allowed to be kept fixed for all times from the earlier stages, and so to have a purely trigonometric solution.For a completely satisfactory work we need many more observations than actually exist. Two kinds of observations seem to be the best suitable: long-focus photographic plates and photometric records of mutual events.The most recent photographic observations are discussed in order to state guidelines for future work. The problem of the precision of Sampson's tables is discussed on the grounds of the recent observations.Paper presented at IAU Colloquium, No. 28, Ithaca, N.Y., August, 1974.     

椭圆星系的模拟试验

本文用树形法,模拟试验由一种星系相互作用模式产生椭圆星系的可能性．并用惯量张量分析方法和曲线拟合对椭球体进行结构分析,得到椭球体的半长轴相对长度．