Introduction

The main purpose of this book is to provide algorithms for direct N-body simulations, based on personal experience over many years. A brief description of the early history is included for general interest. We concentrate on developments relating to collisional direct integration methods but exclude three- and four-body scattering, which will be discussed in a separate chapter. In the subsequent section, we introduce some basic concepts which help to understand the behaviour of self-gravitating systems. The topics covered include two-body relaxation, violent relaxation, equipartition of kinetic energy and escape. Although the emphasis is on collisional dynamics, some of the theory applies in the large-N limit that is now being approached with modern hardware and improved numerical techniques. After these theoretical considerations, we turn to the problem at hand and introduce the general principles of direct integration as a beginner's exercise and also describe the first N-body method.

Historical developments

Numerical investigations of the classical N-body problem in the modern spirit can be said to have started with the pioneering effort of von Hoerner [1960]. Computational facilities at that time were quite primitive and it needed an act of faith to undertake such an uncertain enterprise. Looking back at these early results through eyes of experience, one can see that the characteristic features of binary formation and escape are already present for particle numbers as small as N = 16, later increased to 25 [von Hoerner, 1963].

Introduction

N-body simulations involve a large number of decisions and the situation becomes even more complex when astrophysical processes are added. The guiding principle of efficient code design must be to provide a framework for decision-making that is sufficiently flexible to deal with a variety of special conditions at the appropriate time. Since the direct approach is based on a star-by-star treatment at frequent intervals, this prerequisite is usually satisfied. However, we need to ensure that the relevant tests are not performed unnecessarily. The development of suitable criteria for changing the integration method or identifying procedures to be carried out does in fact require a deep understanding of the interplay between many different modes of interactions. Hence building up the network for decision-making is a boot-strapping operation needing much patience and experience. The aim of a good scheme should be that this part of the calculation represents only a small proportion of the total effort.

This chapter discusses several distinct types of decisions necessary for a smooth performance. First we deal with the important task of selecting the next particle, or block of particles, to be advanced in time. The challenge is to devise an optimized strategy in order to reduce the overheads. Another aspect concerns close encounters, either between single particles or where one or more subsystems already consist of binaries.

Introduction

In the last few years, the subject of dynamical planetary formation has undergone a remarkable transformation. Thus we now have an increasing database of actual observed systems which provides much material for theoretical and numerical work. It therefore seems appropriate to devote a small chapter to direct N-body simulations of planetary systems. The emphasis in the early days was on performing idealized calculations to ascertain whether the eccentricity growth due to weak perturbations could lead to significant accretion by collisions. With the increase of computing power, more realistic modelling has become feasible by direct methods, but powerful tree codes are an attractive alternative. The latter technique has proved effective for studying both planetesimal systems and planetary rings. Following the discovery of many extra-solar systems, the question of stability has become topical. Stability is often addressed using symplectic integrators which are outside the main scope of this book. In the following we distinguish between planetary formation and planetesimal dynamics. This division is somewhat arbitrary but planetesimal simulations are usually concerned with particles distributed in a thin annulus which therefore represents larger systems.

Planetary formation

Although there are early integrations of planetary dynamics relating to Bode's Law [Hills, 1970], it took another decade for the subject proper to get under way.

Introduction

The basic theory of chain regularization [Mikkola & Aarseth, 1990, 1993] is described in chapter 5, while algorithms that deal with different treatments of physical collisions are detailed in chapter 9. Here we are concerned with a number of additional features that deal with aspects relating to what might be termed the N-body interface, namely how to combine two different solution methods in a consistent way. Strong interactions in compact subsystems are usually of short duration, with the ejection of energetic particles a characteristic feature. First we give some algorithms for unperturbed triple and quadruple systems while the more extensive treatment of perturbed chain regularization is discussed in the subsequent sections. As far as the internal chain subsystem is concerned, this requires extra procedures that add to the program complexity and cost. Having selected a suitable subsystem for special treatment, we also need to consider the change of membership and possible astrophysical processes. Moreover, the question of when to terminate a given configuration requires suitable decision-making for the switch to alternative methods, as well as identification of hierarchical stability in order to prevent inefficiency. Finally, since the implementation of the time-transformed leapfrog scheme has many similarities with chain regularization, we include some relevant algorithms here.

Introduction

The question of numerical accuracy has a long history and is a difficult one. We are mainly concerned with the practical matter of employing convergent Taylor series in order to obtain statistically viable results. At the simplest level, the basic integration schemes can be tested for the two-body problem, whereas trajectories in larger systems exhibit error growth on short time-scales. However, it is possible to achieve solutions of high accuracy for certain small systems when using regularization methods. There are no generally agreed test problems at present but we suggest some desirable objectives, including comparison with Monte Carlo methods. Since large simulations inevitably require the maximum available resources, due attention must be paid to the formulation of efficient procedures. The availability of different types of hardware adds another dimension to programming design, which therefore becomes very specialized. Aspects of optimization and alternative hardware are also discussed, together with some performance comparisons.

Error analysis

It is a fact of computer applications that an error is made every time two arbitrary real numbers are added. Hence the task is to control the propagation of numerical errors and if possible keep them below an acceptable level. Since the N-body problem constitutes a system of non-linear differential equations the error growth tends to be exponential, as was demonstrated right at the outset of such investigations [Miller, 1964] and emphasized in a subsequent study [Miller, 1974].

Introduction

A large number of algorithms are connected with regularization. Many of these concern the KS treatment which plays a key role in the N-body simulation codes. In this chapter, we derive some expressions relating to the conversion of regularized time, followed by other considerations of a practical nature. A separate section provides essential details of the Stumpff KS method as employed in an N-body code. This is followed by an algorithmic discussion of KS termination. Next we describe decision-making procedures for unperturbed two-body motion which speed up the calculation by a large factor. Another important feature with the same objective is the so-called ‘slow-down device’, where the principle of adiabatic invariance is exploited. The theory was given previously in connection with chain regularization and here we discuss the KS implementation. Special treatments of stable hierarchies also contribute significantly to enhanced efficiency while retaining the essential dynamics. Finally, the last sections deal with several processes relating to tidal interactions in close binaries that are connected through an evolutionary sequence. We discuss tidal circularization and two-body capture, as well as Roche-lobe mass transfer which all contribute to making star cluster modelling such an exciting and challenging project.

General KS considerations

We first discuss various general features that are applicable to all the KS methods and also include some aspects of the divided difference scheme, while the next section deals specifically with the Stumpff version.

Introduction

We describe some relevant algorithms for initializing and terminating higher-order systems which are selected for the merger treatment discussed in chapter 11. Such configurations may consist of a single particle or binary in bound orbit around an existing triple or quadruple, or be composed of two stable triples. There are some significant differences from the standard case, the main one being that the KS solution of the binary is not terminated at the initialization. In other words, the complexity of the structure is increased, and the whole process is reminiscent of molecular chemistry. However, once formed as a KS solution, the new hierarchy needs to be restored to its original constituents at the termination which is usually triggered by large perturbations or mass loss. Special procedures are also required for removing all the relevant components of escaping hierarchies, and here we include merged triples and quadruples since the treatment is similar to that for higher-order systems.

Initialization

Consider an existing hierarchy of arbitrary multiplicity and mass, mi, which is to be merged with the mass-point mj, representing any object in the form of a single particle or even another hierarchy. Again we adopt the convention of denoting the component masses of a KS pair by mk and ml, respectively. Some of the essential steps are listed in Algorithm C.1.

Nearly all of these steps also appear in the standard case and therefore do not require comment. We note one important difference here, in that there is no termination of the KS solution.

Introduction

In the preceding chapter, we have considered several methods for dealing with the perturbed two-body problem. Formally all these methods work or quite large perturbations, provided the regularized time-step is chosen sufficiently small. However, the selection of the dominant pair in a triple resonance interaction frequently calls for new initializations where the intruder is combined with one of the components. Likewise, one may have situations in which two hard binaries approach each other with small impact parameter. Hence a description in terms of one dominant two-body motion tends to break down during strong interactions, precisely at times when interesting outcomes are likely to occur. Since the switching of dominant components reduces the efficiency and also degrades the quality of the results, it it highly desirable to seek alternative methods for improved treatment.

In this chapter, we discuss several multiple regularization methods that have turned out to be very beneficial in practical simulations. By multiple regularization it is understood that at least two separations in a compact subsystem are subject to special treatment where the two-body singularities are removed. We begin by describing a three-body formulation that may be considered the Rosetta Stone for later developments. The generalization to more members with just one reference body is also included for completeness. A subsequent section outlines the elegant global formulation and is followed by a detailed discussion of the powerful chain regularization.

Introduction

This book spans my entire life as a research worker at Cambridge. The circumstances that created this opportunity were based entirely on luck and this aspect played a vital part during subsequent developments. In the following chapters, I have tried to give details of the most relevant methods used in so-called ‘direct integration’ of the classical N-body problem, a method of attack somewhat analogous to scaling a mountain the hard way. This has been enhanced by an extensive discussion of the main algorithms implemented in the associated computer codes. A comprehensive review of related work in the field over the last 40 years is also presented. Throughout the term N-body simulations is used exclusively for methods based on direct summation, in keeping with tradition.

Although a wide range of problems is covered, the emphasis is on the dynamics of star clusters. This involves many aspects of stellar evolution. It is fortuitous that the University of Cambridge has a long tradition in this field that dates back to Eddington and Jeans. Fred Hoyle continued this school, which eventually gave rise to the application of synthetic stellar evolution. This subject was pioneered entirely at the Institute, mainly by the sequential efforts of Peter Eggleton, Christopher Tout and Jarrod Hurley, whose work has been vital for realistic star cluster simulations.