The theory of collisional relaxation in stellar systems is discussed in terms of an expansion in powers of 1/N, the inverse of the total number of stars. The results are expressed in terms of the concept of gravitational polarization.

The basic assumption in Chandrasekhar’s approach of statistical stellar dynamics (Chandrasekhar, 1942) is the postulate that a test star within a stellar system being stationary in the sense of collisionless continuum theory suffers random displacements in velocity space generated by the fluctuating part of the gravitational field in a manner that can be described in terms of a random walk. This is equivalent to the assertion that the increments of velocity are regarded as stochastically independent in disjoint time intervals. From this Chandrasekhar derived a diffusion process in velocity space.

It is shown that the time of relaxation by particle encounters of self-gravitating systems in the plane interacting by 1/r2 forces is of the same order of magnitude as the mean orbit time. Therefore such a system does not have a Vlasov limit for large numbers of particles, unless appeal is made to some non-zero thickness of the disk. The relevance of this result to numerical experiments on galactic structure is discussed.

This review first discusses the different types of numerical methods available for integrating the equations of motion of N-body systems. It is desirable to supplement ordinary integration schemes with special treatments of close encounters using a two-body perturbation description or introducing regularizing transformations of the co-ordinates and time. Direct methods are at present limited to the study of a few hundred particles but larger systems may be investigated using Monte Carlo techniques or the Boltzmann moment equations.

N-body computations have been performed for a whole range of initial conditions and the general results are summarized. Numerical investigations have already clarified a number of important aspects of cluster evolution and the qualitative behaviour of small stellar systems is now quite well understood. Recent theoretical modifications have reduced the disagreement with experiments but further improvements are still needed.

Certain features of a practical regularisation method, in which the potential or kinetic energy is used as a time-regularising function, are described. For two-body encounters the method is less powerful than Kustaanheimo-Stiefel regularisation, but it has wider applicability.

The numerical integration of systems of differential equations that possess integrals is often approached by using the integrals to reduce the number of degrees of freedom or by using the integrals as a partial check on the resulting solution, retaining the original number of degrees of freedom.

Another use of the integrals is presented here. If the integrals have not been used to reduce the system, the solution of a numerical integration may be constrained to remain on the integral surfaces by a method that applies corrections to the solution at each integration step. The corrections are determined by using linearized forms of the integrals in a least-squares procedure.

The results of an application of the method to numerical integrations of a gravitational system of 25-bodies are given. It is shown that by using the method to satisfy exactly the integrals of energy, angular momentum, and center of mass, a solution is obtained that is more accurate while using less time of calculation than if the integrals are not satisfied exactly. The relative accuracy is ascertained by forward and backward integrations of both the corrected and uncorrected solutions and by comparison with more accurate integrations using reduced step-sizes.

Indications are that instabilities play an important role in many of the phenomena of stellar dynamics. Examples of such phenomena are: the formation of spiral arms, and the evolution of stellar clusters at a rate faster than one would expect from normal two-body collisions. The analogy with the situation in plasma physics, where similar phenomena are known to be dominated by instabilities, is very suggestive that one might seek for instabilities in stellar dynamics that correspond to similar ones in plasma physics.

The statistical mechanics of an isolated self-gravitating system consisting of N uniform mass sheets is considered using both canonical and microcanonical ensembles. The one-particle distribution function is found in closed form. The limit for large numbers of sheets with fixed total mass and energy is taken and is shown to yield the isothermal solution of the Vlasov equation. The order of magnitude of the approach to Vlasov theory is found to be 0(1/N). Numerical results for spatial density and velocity distributions are given.

Gravitational n-body calculations handling 105 particles have become routine in the last three or four years. The plasma physicists have also run calculations with comparable numbers of particles – you will hear of these from Dawson this afternoon. In the stellar dynamics area, two groups have been making these calculations – our own group with Kevin Prendergast and William Quirk, and the competition, the group that formed when Roger Hockney and Frank Hohl merged (and later parted ways). Our calculation is described in Miller and Prendergast (1968); theirs is best described by Hohl and Hockney (1969), and some very thorough discussions of computational speed of force calculations by Hockney (1970).

The evolution of initially balanced rotating disks of stars is investigated with a computer model for isolated disks of stars. An isolated, initially cold balanced disk is found to be violently unstable. Adding a sufficient amount of velocity dispersion will stabilize all small-scale disturbances. However, the disks are still unstable against slowly growing long wave-length modes and after about two rotations most disks tend to assume a bar-shaped structure. It is found that the final mass distribution over most of the disk can be closely approximated by an exponential variation, irrespective of the initial mass distribution. The gravitational two-stream instability is investigated by means of a modified computer model for infinite doubly periodic stellar systems.

Dynamical systems with three degrees of freedom can be reduced to the study of a four-dimensional mapping. We consider here, as a model problem, the mapping given by the following equations:

We have found that as soon as b ≠ 0, i.e. even for a very weak coupling, a dynamical system with three degrees of freedom has in general either two or zero isolating integrals (besides the usual energy integral).

A fourth-order polynomial method for the integration of N-body systems is described in detail together with the computational algorithm. Most particles are treated efficiently by an individual time-step scheme but the calculation of close encounters and persistent binary orbits is rather time-consuming and is best performed by special techniques. A discussion is given of the Kustaanheimo-Stiefel regularization procedure which is used to integrate dominant two-body encounters as well as close binaries. Suitable decision-making parameters are introduced and a simple method is developed for regularizing an arbitrary number of simultaneous two-body encounters.

One of the main difficulties encountered in the numerical integration of the gravitational n-body problem is associated with close approaches. The singularities of the differential equations of motion result in losses of accuracy and in considerable increase in computer time when any of the distances between the participating bodies decreases below a certain value. This value is larger than the distance when tidal effects become important, consequently, numerical problems are encountered before the physical picture is changed. Elimination of these singularities by transformations is known as the process of regularization. This paper discusses such transformations and describes in considerable detail the numerical approaches to more accurate and faster integration. The basic ideas of smoothing and regularization are explained and applications are given.

We give here a detailed technical description of a Monte Carlo scheme for the dynamical evolution of spherical stellar systems. The philosophy of the method, as well as a few illustrative results, are given elsewhere (Hénon, 1971, hereafter called I).

Many problems in stellar dynamics involve phenomena occurring in inhomogeneous systems in which the interaction between the particles is fully described by a self-consistent field operating in phase space. Because the particles interact by means of the long-range Coulomb force, each particle is under the simultaneous influence of a large number of other particles. Therefore, stellar systems will respond to any perturbation in a collective manner, and a study of such systems is concerned essentially with the N-body problem.

A model of the galaxy is constructed and evolved in which the integrated influence of stellar and supernova nucleosynthesis on the composition of the interstellar gas is traced numerically. Our detailed assumptions concerning the character of the matter released from evolving stars and supernovae are guided by the results of recent stellar evolutionary calculations and hydrodynamic studies of supernova events. Stars of main sequence mass in the range 4 ≤ M ≤ 8 M⊙ are assumed to give rise to supernova events, leaving remnants we identify with neutron stars and pulsars and forming both the carbon-to-iron nuclei and the r-process heavy elements in the explosive ejection of the core material. For more massive stars, we assume the core implosion will result in the formation of a Schwarzschild singularity, that is, a black hole or ‘collapsar’. The straightforward assumptions (1) that the gas content of the galaxy decreases exponentially with time to its present level of ~ 5 % and (2) that the luminosity function characteristic of young clusters and the solar neighborhood is appropriate throughout galactic history, lead to the prediction that ≈ 20% of the unevolved stars of approximately one solar mass (M⊙) in the galaxy today should have metal compositions Z ≲ 0.1 Z⊙. As Schmidt has argued from similar reasoning, this is quite inconsistent with current observations; an early generation dominated by more massive stars – which would by now have evolved – is suggested by this difficulty. Many of these massive stars, according to our assumptions, will end their lives as collapsed black hole remnants. It is difficult to visualize an epoch of massive star formation in the collapsing gas cloud which formed our galaxy which would enrich the gas rapidly enough to account for the level of heavy element abundances in halo population stars; we have therefore proposed a stage of star formation which is entirely pregalactic in character. We suggest that the Jeans’ length-sized initial condensations in the expanding universe discussed by Peebles and Dicke may provide the appropriate setting for this first generation of stars. Guided by these considerations, and by the need for a substantial quantity of ‘unseen’ mass to bind our local group of galaxies, we have constructed a model of the galaxy in which this violent early phase of massive star formation produces both (1) approximately 25 % of the level of heavy elements observed in the solar system and (2) an enormous unseen mass in the form of black holes. The implications of our model for other features of the galaxy, including supernova nucleosynthesis, the cosmic ray production of the light elements, and cosmochronology, are discussed in detail.

In this paper, the Combined Operational Method developed by Busbridge (1961) in connection with the radiative transfer problems in plane-parallel atmospheres has been extended to similar problems in isotropic scattering, homogeneous spherical media. The relevant auxiliary equation has been formulated, the scattering function defined and the integro-differential equation for such function deduced. For a medium having radial distribution of source in addition to the incident flux at the outer surface, the integro-differential equations for source function and emergent intensity have been established.

The energy spectra observed in a CsI crystal in the 20 keV–2 MeV range, due to the decay of radioactive isotopes produced in the crystal by bombardment with 155 MeV protons, are presented as a function of time after irradiation.

It is shown that the large number of decay products produced by spallation can account for these spectra and that these spectra are in quantitative agreement with the predictions of a semi-empirical formula due to Rudstam, which gives the numbers of different isotopes produced. This formula is used to predict the spallation that would occur in such a crystal on board a satellite due to cosmic rays and passages through the South Atlantic Anomaly. Inspection shows that the spallation produced in the latter case is well approximated by that at 155 MeV. Hence the experimental results are used to explain previously observed background rates and to predict the background rates that would occur in the U.K.5 X-ray telescope of Imperial College. Using the Rudstam formula an estimate of cosmic ray induced background is also made.

The relative importance of activity resulting from neutron interactions (atmospheric albedo and spacecraft secondaries) is considered.

It is suggested that the Rudstam formula provides a general method of predicting induced radioactivity in satellite materials and that observed breaks in the diffuse cosmic X-ray spectrum could be due to inadequate allowance for this source of background.