In this chapter we will present examples of so-called open boundary problems. By this we mean that the simulation contains all particles relevant to the problem, and the size of the simulation region is adjusted accordingly at each timestep. This type of problem is by far the easiest one to which one can apply hierarchical data structures. The additional difficulties posed by periodic boundaries will be considered in Chapters 5 and 6.

Gravitational Problems in Astrophysics

Hierarchical tree codes were first developed in the context of astrophysics. This is not surprising because there is a big discrepancy between the number of bodies one would like to study – for example, O (1011) for a galaxy – and the number one can afford to model with a standard N-body code – at present O(105). PIC codes, which employ a grid structure to represent the fields in space, usually cannot handle these problems for two reasons: the complex structure of the investigated object and the large density contrasts such as those found in galaxies. N-body codes are able to avoid the first of these difficulties, because they are gridless and can therefore cope with arbitrarily complicated structures. The second difficulty remains, however, because of the N2 scaling of computation time. This can have two consequences. If the number of simulation particles is too small, the spatial resolution of the simulation might not be good enough to reveal the real dynamic behaviour of the system.

In this chapter we will consider a variety of fields where the tree algorithm in combination with periodic boundaries can be applied and where the speedup in comparison to standard MD and MC codes enables previously inaccessible problems to be investigated. This set of applications is in no way exhaustive, but is intended to indicate the types of problems where the algorithm might best be put to use.

Practically every N-body MD or MC code could incorporate the tree algorithm. However, for systems with short-range potentials, like the Lennard–Jones potentials for the description of solids, this does not bring much of an advantage. Because the potential falls off very rapidly (see Fig. 6.1), a sharp cutoff can be used to exclude interactions of more distant particles whose contribution is negligible. The tree algorithm is mainly suited to systems with long-range forces such as Coulomb, where the summed effect of distant particles is important.

An ideal application is dense, fully-ionized plasmas. Here the particles are so closely packed that an analytical treatment is difficult, and many MD and MC calculations have been carried out to investigate their properties. Limitations in the number of simulation particles make some problems difficult to address due to a combination of poor statistics and small system size. The tree algorithm could be successfully applied here, because the particles interact purely through Coulomb forces.

In previous chapters we occasionally referred to an alternative type of tree code, namely the Fast Multipole Method (FMM). This technique, an elegant refinement of the basic Barnes–Hut algorithm, appears to be best suited to ‘static’ problems, where the particle distribution is more or less uniform. Although it has not been as widely used as the Barnes–Hut (BH) method for dynamic problems – because of either its increased mathematical complexity or the additional computational overhead – it may well become the basis of ‘multimillion’ N-body problems in the near future. We therefore include an introduction to FMM here, based primarily on works by Greengard (1987, 1988, 1990) and Schmidt and Lee (1991). At the same time, we will try to maintain a consistency of notation with the Barnes–Hut algorithm (hereafter referred to as the ‘tree method’ or ‘tree algorithm’), as described in Chapter 2.

Outline of the Fast Multipole Algorithm

The Fast Multipole Method makes use of the fact that a multipole expansion to infinite order contains the total information of a particle distribution. As in the BH algorithm, the interaction between near neighbours is calculated by direct particle–particle force summation, and more distant particles are treated separately. However, the distinction between these two contributions is obtained in a different way. In FMM the distant region is treated as a single ‘far-field’ contribution, which is calculated by a high-order multipole expansion.

The FMM was first formulated by Greengard and Rokhlin (1987).

The difficulty in writing a ‘how-to’ book on numerical methods is to find a form which is accessible to people from various scientific backgrounds. When we started this project, hierarchical N-body techniques were deemed to be ‘too new’ for a book. On the other hand, a few minutes browsing in the References will reveal that the scientific output arising from the original papers of Barnes and Hut (1986) and Greengard and Rohklin (1987) is impressive but largely confined to two or three specialist fields. To us, this suggests that it is about time these techniques became better known in other fields where N-body problems thrive, not least in our own field of computational plasma physics. This book is therefore an attempt to gather everything hierarchical under one roof, and then to indicate how and where tree methods might be used in the reader's own research field. Inevitably, this has resulted in something of a pot-pourri of techniques and applications, but we hope there is enough here to satisfy the beginners and connoisseurs alike.

Tree Construction

We have seen in the preceding chapter that in grid-based codes the particles interact via some averaged density distribution. This enables one to calculate the influence of a number of particles represented by a cell on its neighbouring cells. Problems occur if the density contrast in the simulation becomes very large or the geometry of the problem is very complex.

So why does one bother with a grid at all and not just calculate the interparticle forces? The answer is simply that the computational effort involved quite dramatically limits the number of particles that can be simulated. Particularly with 1/r-type potentials, calculating each particle–particle interaction requires an unnecessary amount of work because the individual contributions of distant particles is small. On the other hand, gridless codes cannot distinguish between near-neighbours and more distant particles; each particle is given the same weighting.

Ideally, the calculation would be performed without a grid in the usual sense, but with some division of the physical space that maintains a relationship between each particle and its neighbours. The force could then be calculated by direct integration while combining increasingly large groups of particles at larger distances. Barnes and Hut (1986) observed that this works in the same way that humans interact with neighbouring individuals, more distant villages, and larger states and countries. A resident of Lower-Wobbleton, Kent, England, is unlikely to undertake a trip to Oberfriedrichsheim, Bavaria, Germany, for a beer and to catch up on the local gossip.

In Chapter 2 we saw the basic workings of the tree algorithm. Now we will discuss some methods that can be used to optimise the performance of this type of code. Although most of these techniques are not specific to tree codes, they are not always straightforward to implement within a hierarchical data structure. It therefore seems worthwhile to reconsider some of the common tricks of the N-body trade, in order to make sure that the tree code is optimised in every sense – not just in its N log N scaling.

There are basically two points of possible improvement:

• Improvement of the accuracy of the particle trajectory calculation by means of higher order integration schemes and individual timesteps. This is especially important for problems involving many close encounters of the particles, that is, ‘collisional’ problems.

• Speedup of the computation time needed to evaluate a single timestep by use of modern software and hardware combinations, such as vectorisation, and special-purpose hierarchical or parallel computer architecture.

Individual Timesteps

For most many-body simulations one would like the total simulated time T = ntΔt (where nt is the number of timesteps) to be as large as possible to approach the hydrodynamic limit. However, the choice of the timestep Δt has to be a compromise between this aim and the fact that as Δt increases, the accuracy of the numerical integration gets rapidly worse.

The hierarchical tree method can not be adapted only for Monte Carlo applications: It can also be modified to perform near-neighbour searches efficiently. This means that the tree algorithm could also have applications for systems with short-range or contact forces. Hernquist and Katz (1989) first showed how the tree structure can be used to find near neighbours through range searching. Following their method, the near-neighbour search is performed the following way.

Consider a system in which only neighbours lying within a distance h will interact with the particle i in question. For the near-neighbour search this sphere is enclosed in a cube whose sides are of length 2h. The tree is built the usual way and the tree search starts at the root. The tree search is performed in a very similar way to the normal force calculation of Section 2.2 by constructing an interaction list. The main difference is that the s/d criterion is substituted by the question: ‘Does the volume of the search cube overlap with the volume of the pseudoparticle presently under consideration?’

If there is no overlap, this branch of the tree contains no near neighbours and is not searched any further. However, if there is an overlap, the cell is subdivided into its daughter cells and the search continues on the next highest level. If the cell is a leaf – meaning there is only one particle in the cell – one has to check whether it actually lies within the radius h of particle i.

THE RESTLESS UNIVERSE

From ancient Hindu mythology comes this story about the Pole Star: King Uttanapada had two wives. The favourite, Suruchi, was haughty and proud, while the neglected Suniti was gentle and modest. One day Suniti's son Dhruva saw his co-brother Uttama playing on their father's lap. Dhruva also wanted to join him there but was summarily repulsed by Suruchi, who happened to come by. Feeling insulted, the five-year-old Dhruva went in search of a place from where he would not have to move. The wise sages advised him to propitiate the god Vishnu, which Dhruva proceeded to do with a long penance. Finally Vishnu appeared and offered a boon. When Dhruva asked for a place from where he would not have to move, Vishsnu placed him in the location now known as the Pole Star – a position forever fixed.

Unlike other stars and planets, the Pole Star does not rise and set; it is always seen in the same part of the sky. This immovability of the Pole Star has proved to be a useful navigational aid to mariners from ancient to modern times. Yet, a modern-day Dhruva could not be satisfied with the Pole Star as the ultimate position of rest. Let us try to find out why.

The Pole Star does not appear to change its direction in the sky because it happens to lie more or less along the Earth's axis of rotation. As the Earth rotates about its axis, other stars rise over the eastern horizon and set over the western horizon.

FUSION

It is often argued that man's growing energy needs will be met if he succeeds in making fusion reactors. In a fusion reactor, energy is generated by fusing together light atomic nuclei and converting them into heavier ones. The primary fuel for such a fusion reactor on the Earth would be heavy hydrogen, whose technical name is deuterium. Through nuclear fusion, two nuclei of deuterium are brought together and converted to the heavier nucleus of helium, and in this process nuclear energy is released.

The following is the recipe for a fusion reactor. First, heat a small quantity of the fusion fuel, deuterium, above its ignition point – to a temperature of some 100 million degrees Celsius. Second, maintain this fuel in a heated condition long enough for fusion to occur. When this happens, the energy that is released exceeds the heat input, and the reactor can start functioning on its own. The third and final part of the operation involves the conversion of the excess energy to a useful form, such as electricity.

The primary fuel for this process, the heavy hydrogen, is chemically similar to but a rarer version of the commonly known hydrogen. An atom of ordinary hydrogen is made up of a charged electrical particle called the proton at the nucleus with a negatively charged particle, the electron, going round it. The nucleus of heavy hydrogen carries an additional particle called the neutron in its nucleus. The neutron has no electric charge so the total charge of the nucleus of heavy hydrogen is the same as that of ordinary hydrogen.

It is often said that modern theoretical physics began with Newton's law of gravitation. There is a good measure of truth in this remark, especially when we take into account the aims and methods of modern physics – to describe and explain the diverse and complex phenomena of nature in terms of a few basic laws.

Gravity is a basic force of the Universe. From the motions of ocean tides to the expansion of the Universe, a wide range of astronomical phenomena are controlled by gravity. Three centuries ago Newton summed up gravity in his simple inverse-square law. Yet, when asked to say why gravity follows such a law, he declined to hazard an opinion, saying ‘Non fingo hypotheses’ (I do not feign hypotheses). A radically new attempt to understand gravity was made in the early part of this century by Einstein, who saw in it something of deeper significance that linked it to space and time. The modern theoretical physicist is trying to accommodate it within a unified theory of all basic forces. Yet, gravity remains an enigma today.

In this book I have attempted to describe the diversity, pervasiveness, and importance of this enigmatic force. It is fitting that I have focused on astronomical phenomena, because astronomy is the subject that first provided and continues to provide a testing ground for the study of gravity. These phenomena include the motions of planets, comets, and satellites; the structure and evolution of stars; tidal effects on the Earth and in binary star systems; gigantic lenses in spaced highly dense objects, such as neutron stars, black holes, and white holes; and the origin and evolution of the Universe itself.

THE BLACK HOLE IN HISTORY AND ASTRONOMY

The oldest mention of a black hole is found not in books of physics or astronomy but in books of history. In the summer of the year 1757, Nawab Siraj-Uddaula, the ruler of Bengal in eastern India, marched on Calcutta to settle a feud with the British East India Company. The small garrison stationed in Fort William at Calcutta was hardly a match for the Nawab's army of 50 000. In the four-day battle that ensued, the East India Company lost many lives, and a good many, including the company's governor, simply deserted. The survivors had to face the macabre incident now known as the Black Hole of Calcutta.

The infuriated Nawab, whose army had lost thousands of lives in the battle, ordered the survivors to be imprisoned in what came to be known as the Black Hole, a prison cell in Fort William. In a room 18 feet by 14 feet, normally used for housing three or four drunken soldiers, the 146 unfortunate survivors were imprisoned. The room had only two small windows (see Figure 7–1). During the 10 hours of imprisonment, from 8 p.m. on 20 June to 6 a.m. on 21 June in the hottest part of the year, 123 prisoners died. Only 22 men and 1 woman lived to tell the tale.

Apart from its macabre aspect, the Black Hole of Calcutta did bear some similarity to its astronomical counterpart, involving as it did a large concentration of matter in a small space from which there was no escape.