Introduction

Direct N-body simulations are of necessity expensive because of the need to evaluate all the N(N –1)/2 force interaction terms. We have seen that the Ahmad–Cohen [1973] neighbour scheme only alleviates the problem to some extent. However, once the particle number becomes sufficiently large, the dynamical behaviour begins to change because close encounters are less important. This behaviour has inspired methods for collisionless systems to be developed, such as tree codes or multipole expansions. In this chapter, we are concerned with tree codes since some relevant aspects of the latter have already been discussed. First we review the basic features of the pioneering Barnes & Hut [1986] scheme which is widely used in a variety of applications. Since the emphasis in this book is on collisional stellar dynamics, we devote a section to describing a tree code for point-mass interactions [McMillan & Aarseth, 1993] in the hope that it might be revived. The final section deals with an independent development for flattened systems [Richardson, 1993a,b] that has been used to study different stages of planetary formation as well as ring dynamics, where collisions play an important role.

Basic formulation

In view of the rapid growth in the computational requirements for increasing particle numbers when using direct summation, it is not surprising that several tree-based approaches have been made to speed up the expensive force calculation.

Introduction

Recent years have seen impressive advances in the simulation of collisionless systems. However, a wide range of problems in galactic dynamics are still amenable to direct N-body integration and will therefore be discussed for completeness. This chapter also serves as a historical review of several topics that blossomed from a primitive beginning. All such problems are characterized by the employment of a small softening of the potential which reduces the effect of close encounters. We begin by describing an application of the grid perturbation method to a ring of interacting molecular clouds in the inner part of the Galaxy [Aarseth, 1988b]. This enables a realistic number of clouds to be considered. The tidal disruption of dwarf spheroidal galaxies orbiting the Milky Way has also been studied by direct means [Oh, Lin & Aarseth, 1995].

More extreme interactions of galaxies often result in the formation of one system, subsequently denoted as a remnant to distinguish this process from mergers used for hierarchical stellar configurations. Studies of galaxy interactions involving black holes have also become topical [Makino & Ebisuzaki, 1996; Makino, 1997]. This problem is particularly challenging because of the difficulty of scaling to realistic conditions. Small galaxy groups and clusters are ideal for N-body simulations if close encounters can be treated using a softened potential.

Introduction

A variety of procedures need to be carried out before the calculation proper can begin. The prescriptions for input parameters and options are discussed in chapter 7. Here we concentrate on different types of initial conditions for star cluster simulations, whereas planetary systems are described elsewhere. The cluster models are first generated for single stars with a specified initial mass function (hereafter IMF) and scaled to internal units. Since a variety of distributions may be considered, we provide several detailed algorithms. Next we present some procedures for including a realistic distribution of primordial binaries. Modelling of star clusters also requires external effects to be added. We distinguish between the motion of open clusters in circular orbits and globular clusters in 3D, with the galactic tidal force truncating the outer parts. Interstellar clouds form another perturbing agent which may be taken into account. Finally, with these procedures completed, the force polynomials for direct solutions as well as for any dominant two-body motions can be initialized.

Initial conditions for clusters

Although the choice of starting configurations for star cluster simulations is extremely wide, we may be guided by certain principles and observational constraints. Such models are usually represented by a smooth IMF and centrally concentrated density distribution. Depending on the objectives, the velocities may represent approximate equilibrium or initial collapse, whereas cosmological models are characterized by expansion.

Introduction

In the preceding chapters, we have described a variety of numerical methods and their implementation. Given these tools and algorithms, it should be possible in principle to construct direct N-body codes for dynamical studies. However, the combination of methods introduces considerable complications which have taken many years to master and it is therefore much easier to take over one of the large codes. On the other hand, it would be good programming practice to implement a stand-alone code based on direct summation with softened potential or three-body regularization, where decision-making is considerably simplified.

In the following we provide some hints to facilitate the use of one of the main codes, as well as guidelines for implementing movies. These range from practical comments on getting started to producing a meaningful description of the results. Special counters that record events of interest provide further information. Graphics packages are generally avoided in the codes for reasons of compatibility. However, two versions of stand-alone regularization codes are available with software for movie making which is a good way to study dynamical interactions. We also discuss some diagnostic aspects that may assist in the dissemination of the results and outline various strategies for identifying numerical problems.

Introduction

We now make an abrupt transition to a presentation of various algorithms utilized by the direct summation codes. Before proceeding further, it will be useful to include some practical aspects in order to have a proper setting for the subsequent more technical procedures. First we introduce the main codes that have been developed for studying different gravitational N-body problems. Where possible, the same data structure has been employed, except that the most recent versions are formulated in terms of the Hermite integration scheme. Since the largest codes are quite complicated, we attempt to describe the overall organization by tables and a flowchart to provide some enlightenment. Later sections give further details concerning input parameters, variables and data structure; each of these elements play an important role for understanding the general construction. We also discuss a variety of optional features which provide enhanced flexibility for examining different processes.

N-body codes

Before describing the characteristics of the codes, we introduce some short-hand notation to illustrate the different solution methods employed [cf. Makino & Aarseth, 1992]. Thus by ITS we denote the basic individual time-step scheme, whereas ACS defines the Ahmad–Cohen [1973] neighbour scheme. Likewise, HITS and HACS are used for the corresponding Hermite integration methods. Finally, MREG refers to the implementations of unperturbed three-body [Aarseth & Zare, 1974] and four-body chain regularization [Mikkola & Aarseth, 1990], as well as perturbed chain regularization [Mikkola & Aarseth, 1993].

Introduction

In this last chapter, we return to the subject of numerical experiments in the classical sense. Using the computer as a laboratory, we consider a large number of interactions where the initial conditions are selected from a nearly continuous range of parameters. Such calculations are usually referred to as scattering experiments, in direct analogy with atomic physics. Moreover, the results are often described in terms of the physicist's cross sections, with the results approximated by semi-analytical functions. The case of three or four interacting particles is of special interest. However, the parameters space is already so large that considerable simplifications are necessary. In addition to the intrinsic value, applications to stellar systems provide a strong motivation. Naturally, the conceptual simplicity of such problems has also attracted much attention.

In the following we discuss a number of selected investigations with emphasis on those that employ regularization methods. It is convenient to distinguish between simulations and scattering experiments, where in the former case all particles are bound. The aim is to obtain statistical information about average quantities such as escape times and binary elements, and determine their mass dependence. Small bound systems often display complex behaviour and therefore offer ample opportunities for testing numerical methods. On the other hand, scattering experiments are usually characterized by hyperbolic relative velocities, the simplest example being a single particle impacting a binary.

Introduction

In this chapter, we provide the tools needed for standard N-body integration. We first review the traditional polynomial method which leads to increased efficiency when used in connection with individual time-steps. This self-contained treatment follows closely an earlier description [Aarseth, 1985a, 1994]. Some alternative formulations are discussed briefly for completeness. We then introduce the simpler Hermite scheme [Makino, 1991a,b] that was originally developed for special-purpose computers but is equally suitable for workstations or laptops and is attractive by its simplicity. As discussed in a later section, the success of this scheme is based on the novel concept of using quantized time-steps (factor of 2 commensurate), which reduces overheads. Variants of the Hermite method were attempted in the past, such as the low-order scheme of categories [Hayli, 1967, 1974] and the full use of explicit Taylor series derivatives [Lecar, Loeser & Cherniack, 1974]. The former study actually introduced the idea of hierarchical time-steps with respect to individual force calculations using a low-order scheme, whereas the latter formulation is expensive (but accurate) even for modest particle numbers.

Force Polynomials

The force acting on a particle usually varies in a smooth manner throughout an orbit, provided the particle number is sufficiently large. Hence by fitting a polynomial through some past points, it is possible to extend the time interval for advancing the equations of motion and thereby reduce the number of force evaluations. In other words, we can use the past information to predict the future motion with greater confidence.

Introduction

Simulating star clusters by direct N-body integrations is the equivalent of scaling mountains the hard way. At any time the maximum particle number depends on hardware and is therefore limited by technology. Some of the methods that have been described in this book are ideally suited to studying the classical point-mass problem. In addition, a wide variety of astrophysical processes can be included for realistic modelling of actual clusters. Recently the simulations have been devoted to systems with up to N ≃ 104 particles which includes rich open clusters. However, with the construction of the GRAPE-6 special-purpose computer we are now able to investigate small globular clusters as observed in the Large Magellanic Cloud (LMC) [Elson et al., 1998].

In the following we concentrate on various aspects of star cluster simulations not covered in earlier chapters. We first describe algorithms for determining the core radius and density centre which are useful tools for data analysis. For historical reasons, idealized models (i.e. isolated systems) are also considered, particularly because of their relevance for more approximate methods. After further discussions of the IMF, we return to the subject of assigning primordial binaries and illustrate their importance by some general results. External effects due to the tidal field and interstellar clouds form an important ingredient in star cluster modelling even though the latter are rarely studied. Algorithms for combining stellar evolution with the dynamical treatment have been outlined previously.

Introduction

Direct N-body simulations on conventional computers benefit greatly from the use of the Ahmad–Cohen [1973] or AC neighbour scheme. Algorithms for both the divided difference method and Hermite formulation will therefore be discussed in the following sections. We also consider the implementation of the code N BODY6++ on a popular type of parallel computer [Spurzem, Baumgardt & Ibold, 2003], since it seems that the future of large-N calculations is evolving in this direction at least for those who do not use the special-purpose HARP or GRAPE machines described previously. The important problem of massive black hole binaries in galactic nuclei is very challenging and appears amenable to direct integration using parallel architecture and neighbour schemes. A direct solution method is described [Milosavljević & Merritt, 2001]. This treats the massive components by two-body regularization, whereas the formation process itself is studied by a tree code. Some of the drawbacks of this method inspired a new formulation where the massive binary is considered as part of a compact subsystem which is advanced by a time-transformed leapfrog method [Mikkola & Aarseth, 2002]. Over the years, the quest for larger particle numbers has also encouraged the construction of partially collisional methods.

Introduction

Most star clusters are characterized by large memberships that make direct N-body simulations very time-consuming. In order to study such systems, it is therefore necessary to design methods that speed up the calculations while retaining the collisional approach. One good way to achieve this is to employ a neighbour procedure that requires fewer total force summations. The AC neighbour scheme [Ahmad & Cohen, 1973, 1974] has proved very effective for a variety of collisional and collisionless problems. It is particularly suitable for combining with regularization treatments, where dominant particles as well as perturbers can be selected from the corresponding neighbour lists without having to search all the particles.

The AC method can also be used for cosmological simulations. The study of galaxy clustering usually employs initial conditions where the dominant motions are due to the universal expansion. By introducing comoving coordinates, we can integrate the deviations from the smooth Hubble flow, thereby enabling much larger time-steps to be used, at least until significant density perturbations have been formed. Naturally, such a direct scheme cannot compete with more approximate methods, such as the P3M algorithm [Efstathiou & Eastwood, 1981; Couchman, 1991] or the parallel tree code [Dubinski, 1996], which have been developed for large N, but in some problems the ability to perform detailed modelling may still offer considerable advantages.

Realistic N-body simulations of star cluster evolution require a substantial programming effort. Since it takes time to develop suitable software, published descriptions tend to lag behind or be non-existent. However, one large team effort has reached a degree of development that merits detailed comments, especially since many results have been described in this book. In the following we highlight some aspects relating to the integration method as well as the treatment of stellar evolution, based on one available source of information [Portegies Zwart et al., 2001].

N-body treatment

The kira integrator advances the particle motions according to the standard Hermite method [Makino, 1991a] using hierarchical (or quantized) time-steps [McMillan, 1986]. An efficient scheme was realized with the construction of the high-precision GRAPE computers which calculate the force and force derivative and also include predictions on the hardware.

One special feature here is the use of hierarchical Jacobi coordinates which is reminiscent of an earlier binary tree formulation [Jernigan & Porter, 1989]. This representation is equivalent to the data structure used in KS and chain regularization.

Introduction

Sooner or later during the integration of an N-body system close encounters create configurations that lead to difficulties or at best become very time-consuming if studied by direct methods. On further investigation one usually finds a binary of short period slowing down the calculation and introducing unacceptable systematic errors. Moreover, the eccentricity may attain a large value that necessitates small time-steps in the pericentre region unless special features are introduced. It can be seen that a relative criterion of the type (ηR/|F|)½ for a binary yields an approximate time-step Δt ∝ R3/2, where R is the two-body separation. From this it follows that eccentric orbits require more steps for the same period. Even a relatively isolated binary may therefore become quite expensive to integrate as well as cause a significant drift in the total system energy. It is convenient to characterize the systematic error of a binary integration by the relative change per Kepler orbit, α = Δa/a. For example, using the basic Hermite method we find α = –1.3×10–6 with 270 steps per orbit and an eccentricity e = 0.9. At this rate of inward spiralling, the binary energy would be significantly affected after 104 – 105 periods. Although better behaved, less eccentric systems are also time-consuming, giving α = –4×10–8 for e = 0.2 and 135 steps per orbit.

Clues to the chemical evolution of the Galaxy

The evolution of the composition of matter can be traced back through the various ages of the Galaxy by systematically examining surface abundances over a vary large population of stars by means of spectroscopic analysis (Table 8.1). One is particularly interested in elements observed in the spectra of ancient suns in the galactic halo. These little stars, the oldest we know of, are still shining valiantly today, boasting their exceptional longevity (Fig. 8.1).

Let us now describe the method used. The most accessible elements are those possessing clear lines in the optical spectra of these fossilised objects. In contrast, certain elements like neon and argon are not determined in these stars, whether they be dwarfs or giants. In their normal state, the noble gases produce no optical emission.

Families that lend themselves best to this evolutionary analysis are:

• the light nuclei Li, Be and B;

• the α nuclei, i.e. multiples of the helium nucleus, such as Mg, Si, S and Ca;

• nuclei around the iron peak, viz. Sc, Cr, Mn, Fe, Co, Ni, Cu and Zn;

• heavy s and r isotopes like Sr, Y, Ba and Eu.

Among these, iron is relatively easy to measure and serves as a reference, as a metallicity index, and thus as an indicator of the degree of evolution. Indeed, it is common practice in astronomy to treat the terms iron content (Fe/H) and metallicity (Z) as synonymous.

Does invisible matter exist?

Two arguments support the idea that some invisible substance exists in the Universe. The first is dynamic. It starts from observation of motions under the effect of gravity. The second is related to Big Bang nucleosynthesis, i.e. nuclear cosmology, which combines cosmology and nuclear physics.

Dynamical proof

The first observation is that, if we can go by what light is telling us, most of the matter in the Galaxy (and indeed any galaxy) is concentrated in the galactic bulge, a marked, reddened swelling at the centre of the star distribution. Therefore, if we assume that the mass distribution of luminous objects is representative of the total mass distribution in galaxies, every spiral galaxy should behave like a vast Solar System, with the stars and clouds playing the role of planets and the galactic bulge that of the Sun.

Nuclear statistical equilibrium

The most tightly bound nuclei, i.e. the most stable and robust, in the iron peak are not symmetric arrangements bringing together equal numbers of protons and neutrons (N = Z). Rather, they possess a neutron excess (N − Z) between 2 and 4. Close to iron, the most stable nucleus 56Fe has a number of neutrons which exceeds the number of protons by 4 units (N − Z = 4).

The isotopic and elemental abundance table shows that, in the Solar System, iron is more abundant than its neighbours. Analysis of stellar spectra confirms this result, giving it a universal character.

Theoretically, nuclear strength is enhanced by internal transmutations of protons into neutrons, under the mandate of the weak interaction, either by positron emission (p → n + e+ + ν) or by electron capture (p + e− → n + ν). However, the weak interaction is much slower than the strong interaction. The question remains as to whether it will happen inside the star, or outside, once the matter has been expelled, i.e. after the explosion. This is not just an academic question. The answer we give will determine whether or not we can corroborate explosive nucleosynthesis by observation.

Any attempt to understand the conditions in which iron and its kin were created, and identify the astrophysical site of their birth, must focus on the idea of nuclear statistical equilibrium. The situation is the exact nuclear analogy of the ionisation equilibrium occurring in hot gases.