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.

Think of a scene rare on home-grown European TV, but common in the USA, which abounds with images of games of skill. Tenpin bowling is a highly refined art to those who are experts. In this sport, there's science in every step. From the calculation of the trajectory to the study of the resistance of the lanes, which are not uniformly rough. You can empathise with the fact that many years of skill are required to add the correct spin to the ball in order to attain, at the end of the lane, just the right curve which lets it hit the pins just slightly to the side in order that they fall like dominoes.

For beginners, the feat is infinitely more difficult. Either the ball, clearly wanting to be uncooperative, rushes off into the gulley, or else in spite of a nice, straight and apparently effective trajectory, it only removes the middle pins, leaving two separate groups of survivors, thereby removing all hope of cleaning out the set on the second throw.

So, frustration is often the lot of the novice, who is left with no option other than to persist in the hope that maybe one day. … Luckily, novices can count on an irregular ally: chance. Nothing seems to distinguish one bowl from the previous ones, yet miraculously the ball starts out in the lane, glides smoothly along before arriving at the rough bit, which it rapidly grasps, and perfectly curves its trajectory in just the way needed to leave no pin standing.

The discovery of the exoplanets is undoubtedly a technological feat. But without exceptional people, there would not be any technological feats. So it's equally – and maybe even primarily – a human adventure, a personal experience, and, for this reason, is better narrated in the first person singular, in the voice of Michel Mayor. This will essentially be the case in the first chapter, but also every now and then in later chapters (especially Chapters 6, 7 and 8). These uses of ‘I’ are like memories which suddenly bubble up to the surface during particular scenes. And during the numerous interviews that took place between the scientist and the journalist, there were many of these bubblings. So, henceforth, ‘I’ and Michel Mayor will be indistinguishable.

Having made this comment, all that remains to be done is to set the scene by beginning with a trip through time and space. It's October 1995, in Florence, the seductive Tuscan town where art and science live hand in hand in idyllic happiness. It's there that it all starts, where the discovery would see the light of day.

For Didier Queloz, my young collaborator, it's his first trip ever to Italy. Given the (happy) circumstances which brought us here, he decided to celebrate the event, sharing a room in a beautiful hotel, combining luxury and quaint charm, with his wife, Valérie.

From the most ancient times right up to the sixteenth century, astronomers knew of only five planets (Mercury, Venus, Mars, Jupiter, Saturn) as well as the Earth, which was delicately nested at the heart of the Universe. For a very long time, this geocentrism constituted the dominant vision of the Western world, right up to the day when Nicolas Copernicus put the church back in the middle of the village and the Sun in the centre of the Solar System.

The Earth no longer reigned at the centre of everything. It became an appendage of the day star, a planet like any other. The Universe was turned topsy-turvy. And in addition to this, there were Galileo's Medicean planets. The Solar System was taking shape, and there was nothing to stop the telescope from finding new worlds. Anything was possible, except, perhaps, finding a planet beyond Saturn. It was still thought that the ‘lord of the rings’ ended the world of the planets and that after it there was nothing but stars. It took nearly two centuries for this model to be laid to rest.

HINTS OF URANUS

The German William Herschel was born in 1738 in Hanover, Prussia. His family was a big one. His parents, Isaac and Anne, had ten children. Four died young. The other six were raised to the regular rhythm of scales played by their father, an oboe player in the military band of the Guards of Hanover.

It had to happen at some time or another that someone would look up at the sky and wonder about the nature of the stars. When did this first happen? Undoubtedly, long ago. The first explicit clues linked with astronomical activities date to several millenia before Christ. Just think of Stonehenge, the famous site in England, or of some of the ancient ruins inherited from the Sumerian and Babylonian civilisations.

It was with Greeks that astronomy started to distance itself from the influence of myth and religion. The sky, as well as the Earth, became an object of study, an object of observation, an object of science. Nature became less and less spiritual, and more and more material. However, the arrival of Greek thought did not stop speculation.

In the fourth century BC, the Greek philosopher Epicurus (341–270 BC) asked the fundamental and dizzying question: are we alone in the Universe? Nowadays we know that this question has a real scientific relevance. At the time, it was much less obvious. For the immense majority of Epicurus' contemporaries, at least for those interested in the question, the Universe was closed, bounded by a sphere on which the stars were fixed.

Will we find extraterrestrial life? Has another planet in the Universe succeeded in assembling the extraordinary rainbow of conditions that life seems to require in order to appear? This is the ultimate question that lies behind the quest for exoplanets. But even if you're a member of the club of those who believe that life is not a terrestrial privilege and that it has undoubtedly developed on other planets, in other solar systems, there's still a challenging question: how will we find extraterrestrial life?

If one day humanity succeeded in finding such proof, we would confront the fourth cultural shock of our history. After having learnt from Copernicus that we're not at the centre of the Universe, from Darwin that we're the ‘descendants’ of an ape who herself is the very distant grandchild of a simple cell, and from Freud that we're subject to the whims of our subconscious, we would also have to cope with the idea that we're not the only living beings in the Universe.

Based on our present knowledge, it's becoming more and more difficult to imagine that the Earth is the only host of life in the Cosmos. In our Galaxy alone, there are more than 100 billion stars, while there are billions of galaxies in the observable Universe. Why would life have contented itself with appearing on a single planet, as beautiful and blue as it is?

In May 1998 the press announcement that followed the observation of a young binary in the Taurus constellation took the whole astronomical community aback. However, since it was based on observations by the Hubble Space Telescope (which was then the most sharp-sighted telescope available, able to resolve fine details of distant objects more than 10 billion years old) there was no reason to doubt it. Susan Tereby of the Extrasolar Research Corporation in Pasadena had aimed the Near Infrared Camera and Multi-Object Spectrometer at the Taurus constellation which is known to host numerous young stars and by chance observed a binary just at the moment of formation (it was a few thousand years old at most) at about 450 light-years from the Sun. And this binary exhibited some unexpected properties.

Streaming off this stellar couple is a long filament of luminous gas extending nearly 200 billion kilometres and ending at a bright, point-like object. Together they look like a sort of cosmic exclamation mark which perfectly illustrates the American team's puzzlement, though they searched for an explanation. It was suggested that the luminous point was a giant planet of 2–3 jovian masses, which, victim of the gravitational interactions between the two stars being born, could have simply been catapulted out of the system. The light trail could be matter from the protoplanetary disc trailed by the planet.