The shortest length and time scales of interest in globular cluster
evolution are posed by the closest encounters of the most compact
types of stars. This leads us to a near-grazing encounter of two
neutron stars, an event with a duration of order of a few
milliseconds. Compared with the age of a globular cluster,
years, we have a discrepancy in time scales of no less than 20 orders
of magnitude!
A similar problem shows up when we compare length scales. The
diameter of a neutron star, expressed in units of the tidal radius of
a globular cluster, is of order . At first sight, modeling
a globular cluster on a star-by-star basis would seem pretty hopeless.
The fact that we are able to contemplate such an enterprise at all is
largely due to the ingenuity and persistence of Sverre Aarseth, who
has provided a framework in which to tackle these seemingly
unsurmountable problems (together with efficient implementations).
The key words for facing up to the length scale and time scale
problems are individual timesteps and local coordinate
transformations.
Time Scale Problem
In the course of the sixties, it was gradually realized how important binaries are in determining the dynamics of star systems large or small - and at the same time how devastating their presence was for one's computer budget. Even in a modest 10-body simulation, the formation of a single hard binary is guaranteed to slow down the speed of the whole simulation by orders of magnitude, when using one of the standard integration schemes for coupled differential equations.
The reason is that even with variable time steps, these standard schemes force all stars to share the same time step size. As a result, a typical timestep size of a few hundred years in a loose association of ten stars can be reduced to a month or less with the first moderately tight binary appearing on the scene.
The answer to this problem, provided by Aarseth (cf. Aarseth 1985), was the introduction of individual time steps. Instead of viewing a star cluster with three-dimensional eyes, as made up of a set of mass points in space, Aarseth took a four-dimensional view, in which each star was replaced by an orbit segment. Whenever the time had come for a particular star to move, it could determine the gravitational acceleration by all its neighbors by asking them to slide along their orbit segments (implemented as polynomial approximations) to the desired point in time, even though each of them in turn had computed its own positions, velocities, and accelerations only at earlier times.
This single algorithmic improvement did more to speed up star cluster calculations than decades worth of hardware speed improvement. Even on today's fastest machines, it would be difficult to reproduce Aarseth's calculations of twenty years ago without using individual timesteps.
Length Scale Problem
The large range in length scales, although somewhat smaller than the range in time scales, has turned out to pose a more significant problem for simulations of star cluster evolution. In contrast to the time scale problem, which simply slows things down to a crawl, the length scale problem can easily make a whole calculation meaningless.
The main problem lies in round-off errors. When we take a binary consisting of two neutron stars, moving in the outskirts of a globular cluster, there is no need to follow its internal orbit as long as other perturbers are far away (another significant software improvement provided by Aarseth). However, in the rare case that a third star would interact with that type of binary, we would have to follow all three point-particles numerically, during the time of the interaction.
And here the problem shows up most clearly: for the position vectors of the two neutron stars, with respect to the center of mass of the cluster, the first fifteen digits could turn out to be the same, as we saw above. Alas, even a double precision (64-bit) representation of floating point numbers typically carries a mantissa of only fifteen significant digits. So much for integrating equations of motion.
Admittedly, this is an extreme case, but it makes a point: loss of
accuracy due to round-off when subtracting large numbers is a very
serious worry in -body calculations. Aarseth's answer was to
implement a series of ever more intricate treatments of two-body,
three-body and more complex multiple encounters. These treatments are
largely based on Kustaanheimo-Stiefel regularization procedures, in
which the three-dimensional Kepler singularity is `unfolded' by a
coordinate transformation to four dimensions, mapping the
three-dimensional Kepler problem into that of a four-dimensional
harmonic oscillator through the inverse of one of the Hopf maps from
to
(cf. Stiefel &Scheifele 1971), introducing a U(1)
gauge symmetry in the process.
Code Development
Aarseth has implemented various other algorithmic improvements as
well, most notably the Ahmad-Cohen neighbor scheme, in which the
individual time step scheme is further refined to a two-time-scale
approach (by a more frequent calculation of nearby interactions,
compared to remote interactions). The resulting code, NBODY5, that
includes all these improvements, has been the tool of choice for any
type of detailed star cluster -body simulation, for well over a
decade. The community owes a debt of gratitude to the generosity and
dedication of Sverre Aarseth, who not only has made his codes widely
available, but has made himself available as well - for friendly
advice as much as for code repair and maintenance.
All these developments notwithstanding, there are still some major
problems facing us, if we want to carry out a star-by-star -body
simulation of globular cluster evolution. First of all, if the past
can offer any guidance, increasing the number of particles by a
significant factor is guaranteed to uncover new algorithmic
bottlenecks, requiring new solutions and fine-tuning. The present
state of Aarseth's codes is the outcome of an evolutionary process of
successive attempts to deal with increasingly harder problems, posed
by the increase in complexity of the systems to which it has been
applied. This process is likely to continue for quite a while.
Secondly, the large set of heuristic improvements, valuable as they have been in making calculations possible in the first place, pose a formidable problem to the implementation of stellar evolution recipes. Even simple questions such as the value of the distances between various particles at a given time becomes less straightforward as it may seem at first, when one realizes that one has to trace these distances in the guise of their four-dimensional transformations in which time itself is an extra coordinate, given as a function of the independent integration parameter. Add to this the veritable complexities in the Ahmad-Cohen bookkeeping of the two time scales used (separately for each individual particle, introducing a multiplicity of orbit segment representations), and the outline of the daunting task one is facing begins to appear.
Thirdly, and equally important, it would seem less than ideal if all
simulations in an entire field of astrophysics would be continued with
a single code, without any independent form of comparison or
calibration. In a sense, winning the star cluster -body space
race hands-down, twenty years ago, has been a mixed blessing for
Sverre Aarseth. He has gained a lot of friends, but at the same time
has lacked any significant form of competition. It would seem high
time for someone to explore alternative approaches to
-body code
writing.
A Recursive Approach
For all these reasons, it has been clear for many years that an independent approach was called for. With the prospect of teraflops speeds becoming available in the near future, Jun Makino, Steve McMillan and I decided that the time had come to take on the somewhat daunting task of engaging Sverre Aarseth in an amicable competition. Apart from the adviced to `stay tuned', I briefly sketch what lies at the core of our approach (Hut et al. 1993).
The central idea we have introduced is that of recursive coordinate
transformations. The underlying notion is that of hierarchical
simplicity. Rather than giving a special treatment to each of a
number of different closely interacting groups (binaries, triples,
binary-binary encounters, etc.), we use a recursive approach to split
up an interacting group of stars in subgroups, down to the level of
individual stars. The advantage is twofold: we can now handle
arbitrary -body subsystems, with
as well; and we can treat
especially tight subgroups that may appear deep inside an already
tight group (or even tighter sub-subgroups).
Most stars are unaffected by these special treatments, and are
simply represented as leaves of a flat top-level tree. The stars
that take part in close encounters or are members of multiple star
systems, however, play a crucial role in cluster simulations. The
bulk of NBODY5, for example, is concerned with special treatments of
closely interacting groups of stars. In addition, in most
simulations so far most of the computer time has been taken up by the
integration of the equations of motion of this stellar minority.
In our treatment of such groups of closely interacting stars, the center of mass of the group is represented by a node in the top-level tree. This node in turn serves as the local root node for a binary tree (a tree with two branches per node), in which each star forms a leaf. The structure of this tree is changed dynamically, to guarantee that the tree structure closely reflects, at any time, the configuration of the (sub)groups of the stars. Our approach somewhat resembles that taken by Jernigan and Porter (1989), the major differences being that we limit ourselves to small subsets of stars, and at least for the time being forgo any fancy regularization technique.
Whether we can get away with a bare-bones set of recursive coordinate transformations, rather than four-dimensional regularizations, remains to be seen. But in any case, our recursive approach to a dynamical maintenance of a tree configuration is likely to alleviate our task of providing a clean interface between the stellar dynamics and the stellar evolution parts of our code.
Clearly, the art of -body modeling is a vast subject that still
leaves room for many novel approaches. Extending its realm of
application to the dynamics of
stars with, say,
primordial binaries, is a sure-fire way to stimulate further
exploration.