Contact us Heritage collections Image license terms
HOME ACL Associates Technology Literature Applications Society Software revisited
Further reading □ OverviewRadiative transfer in planetary atmospheresRelativistic Hartree-FochN-body ProblemTransport TheoryComputational modelling of radiative transferConference, 1970
ACD C&A INF CCD CISD Archives Contact us Heritage archives Image license terms

Search

   
ACLApplicationsAstrophysics :: Astrophysics at ACL
ACLApplicationsAstrophysics :: Astrophysics at ACL
ACL ACD C&A INF CCD CISD Archives
Further reading

Overview
Radiative transfer in planetary atmospheres
Relativistic Hartree-Foch
N-body Problem
Transport Theory
Computational modelling of radiative transfer
Conference, 1970

The N-Body Problem

S. J. Aarseth

1968

Expanding Universe

Ever since the introduction of Newton's law of gravitation astronomers and mathematicians have been interested in the star cluster problem. It is paradoxical that the law governing the motion of each member is so simple and yet only a general description can be given for the dynamical behaviour of more than a few particles. In fact exact solutions are known only for the two-body problem and a few special cases for N = 3.

The simple formulation of the problem of interacting mass-points lends itself readily to a numerical attack. Given the equation of motion for a particle with index i, co-ordinate ri and mass mi, viz.

d 2 r 1 d t 2 = - G j = 1 , j i N m j r i - r j | r i - r j | 3 d 2 r 1 d t 2 = - G j = 1 , j i N m j r i - r j r i - r j 3 (1)

where G is the universal constant of gravitation, the problem is to find the future positions of all members as a function of the time t. A unique solution of (1) exists for any given set of initial conditions since the 6N required constants are provided by the initial position and velocity vectors of each particle. For convenience of numerical solution we replace the 3N second-order differential equations by an equivalent set of 6N equations of first order. Denoting by Fi the force per unit mass obtained from the summation in (1) we can then write a Taylor series solution for positions and velocities, valid for a small interval of time Δti.

d r i t d t = d r i 0 d t + F i 0 Δ t i + 1 2 d F i 0 d t Δ t i 2 + ... d r i ( t ) d t = d r i ( 0 ) d t + F i ( 0 ) Δ t i + 1 2 d F i ( 0 ) d t Δ t i 2 + ...

r i t = r i 0 + d r i 0 d t Δ t i + 1 2 F i 0 Δ t i 2 + 1 6 d F i 0 d t Δ t i 3 + ... r i ( t ) = r i ( 0 ) + d r i ( 0 ) d t Δ t i + 1 2 F i ( 0 ) Δ t i 2 + 1 6 d F i ( 0 ) d t Δ t i 3 + ... (2)

where all quantities on the right-hand side are known at a given time, say t = 0. The size of the integration step Δti depends on how many orders are retained in the expansion and is closely related to the desired accuracy. Thus in principle any accuracy can be obtained, provided that the step is small enough, but a practical limit is set by the finite size of the computer word. New values of the accelerations are required after advancing the velocities and positions by one step. It is evident that most of the computer time will be spent on the summation (1) for large particle numbers; e.g., αN2 operations are needed to obtain all accelerations whereas one integration step only requires βN operations, α and β being constants with similar values. However, a considerable gain in efficiency can be achieved by choosing different time-steps Δti for each particle, since the expansion (2) will not be equally convergent for all members. Dynamically this amounts to saying that a longer step suffices to integrate the smoother parts of an orbit so that the relative error tends to be preserved. An additional 25% of the computer time is required to obtain the interpolated positions of all particles at the time of a force calculation but this is a small price to pay for an individual time-step method where typically the steps at a given time may vary by factors of 103 or more.

The storage requirement for the scheme is 28N when terms including third derivatives of the accelerations are retained explicitly. The machine time requirement is, however, more serious since it takes about one Atlas hour to follow the motion of 100 particles over one mean diameter and the crossing time increases as N2 for larger systems. For this reason a small but time-consuming part of the program was written in basic machine language, thereby gaining about 40% in efficiency over Hartran. Also a somewhat more sophisticated scheme than indicated by eqn. (2) is used in the actual computation and a detailed discussion is given in ref. (1). Normally one is interested in cluster simulation over many crossing times in order to obtain meaningful statistics about the dynamical evolution. Luckily there are many real clusters, both of stars and individual galaxies which have less than a few hundred members and are therefore within reach of calculation. Thus in one hour the machine performs a star cluster model calculation for 100 members which in real time would take about 107 yrs.

Before describing some results it should be emphasized that the numerical approach cannot yield the unique solution for all time specified by a set of initial conditions. Such an aim must be abandoned because of the finite word-length as well as the need to approximate a continuous solution by a sequence of integration steps. Spurious errors are therefore introduced into the solutions which diverge because of the strong non-linearity of gravitational interactions. A deterministic description of individual positions and velocities must therefore give way to a probabilistic distribution of orbits. The results may then be used for statistical purposes or a description of global properties, provided that the evolution rate depends mainly on the overall distribution and systematic errors do not enter into the calculation of individual encounters. Even though nothing is known about the detailed behaviour of an N-body system one can use the ten general constants of the motion as integration tests. The conservation of the total energy, angular momentum and centre of mass motion is no guarantee for the accuracy of a solution since errors may cancel, but this is not very likely. An additional and more detailed check, however, is provided by the reversibility of the equations of motion and such tests show that significant dynamical interactions can be retraced to a sufficient accuracy and are therefore physically meaningful.

A case with N = 100 simulating a rotating cluster of galaxies has been studied for about 100 initial crossing times. During this period substantial evolution takes place and 20 of the cluster members escape from the bound system with enough kinetic energy to reach infinity. Clearly such members do not play an important role in the further evolution of the remaining cluster, and by eliminating distant escapers one also saves computer time. The escape process itself can be understood in terms of an energetic interaction involving one or more members in which there is an exchange of kinetic energy. This process is particularly efficient when as in the above example different masses are included since it represents the tendency for a sharing of the kinetic energy of the particles involved. One therefore has a qualitative explanation of the preferential evaporation of the lighter members.

As a consequence of the conservation of total energy it follows that the more massive bodies tend to increase their central concentration. The formation of a dense nucleus leads to strong interactions between a few heavy members which may result in the formation of binaries by dynamical means. A binary represents an energy sink and further interactions may increase its share of the total binding energy. Such encounters often produce highly energetic escapers, especially if a third body has been temporarily captured to form a triple system. It is possible for a substantial proportion of the total energy to be absorbed in a close binary and this process is very suggestive of a final state where only a central binary is left. However, the numerical difficulties of studying such cases of advanced evolution are formidable.

An encounter which is insufficient to produce an escaper leads to one of the particles moving out into the halo where it spends a considerable time following a smooth orbit before returning for further interactions. Hence an initially bound cluster model selected at random from a spherical distribution of constant density tends to increase its volume considerably by evolving towards a more centrally concentrated state. The dynamical importance of the rotation decreases as a result of the expansion and the evolved system does not show a flattening as expected on theoretical grounds. A further 61 distant halo orbits are also removed from the main cluster because they would become part of the general field of galaxies and therefore not be recognized as members by an observer. This brings out one difficulty of comparison with observations since isolated systems do not exist in nature. However, it is relatively easy to modify the force law to include external effects such as a galactic tidal force acting on star clusters and such models have already been considered. A further complication as far as star cluster models are concerned is the assumption of keeping individual masses constant; e.g., one would expect modifications of the mass spectrum by physical evolution of the more massive members which are also most important dynamically.

A simulated cluster tries to set up an equilibrium distribution but particles are continually lost by evaporation and this leads to further adjustments in the structure. No general equilibrium configuration is known for more than two bodies and the evolution therefore proceeds in this direction although the last stage must take a very long time to reach. However, it is meaningful to discuss the overall properties of a system in quasi-equilibrium which is approached slowly in the centrally concentrated state. The distributions of density and velocity may then be compared with theoretical values which have been obtained by analytical considerations. In fact, the direct calculations are made with a view to test theories using different types of assumptions. A relaxation time for the gravitational interactions may also be computed from the change in binding energy of a particle along the orbit and the agreement with theory is quite satisfactory for homogeneous systems. The resulting escape rates, however, are rather more uncertain, even for the simpler case of equal masses. It is found that the mass loss is increased considerably when distributions of different mass groups are studied. Such cases are also more realistic since it is known that actual star clusters contain a large number of light members.

No firm conclusions about the disintegration rate of small star clusters exist as yet but numerical work has shown that the results depend upon the adopted mass spectrum, rather than the selected values of initial positions and velocities. However, there can be little doubt that the half-lives of such systems are considerably shorter than the age of the Galaxy. So far the numerical effort has been relatively modest but it has increased our general understanding of the dynamics of stellar systems and the direct method of cluster simulation on computers therefore offers great promise for the future.

Acknowledgements

It is a great pleasure to thank the Director and staff of the Atlas Computer Laboratory for the generous provision of computer time as well as the excellent service.

References

1. S. J. Aarseth, Bulletin Astronomique, 3, 1968.

⇑ Top of page
© Chilton Computing and UKRI Science and Technology Facilities Council webmaster@chilton-computing.org.uk
Our thanks to UKRI Science and Technology Facilities Council for hosting this site