Nbodie Problem
Nbodie Problem
Nbodie Problem
Douglas C. Heggie
University of Edinburgh,
School of Mathematics,
King’s Buildings,
Edinburgh EH9 3JZ,
UK
1 Introduction
Let a number, N, of particles interact classically through Newton’s Laws of
Motion and Newton’s inverse square Law of Gravitation. Then the equations
of motion are
j=N
X ri − rj
r̈i = −G mj . (1)
j=1,j6=i |ri − rj |3
where ri is the position vector of the ith particle relative to some inertial
frame, G is the universal constant of gravitation, and mi is the mass of the
ith particle. These equations provide an approximate mathematical model
with numerous applications in astrophysics, including the motion of the moon
and other bodies in the Solar System (planets, asteroids, comets and meteor
particles); stars in stellar systems ranging from binary and other multiple
stars to star clusters and galaxies; and the motion of dark matter particles
in cosmology. For N = 1 and N = 2 the equations can be solved analyti-
cally. The case N = 3 provides one of the richest of all unsolved dynamical
problems – the general three-body problem. For problems dominated by one
1
massive body, as in many planetary problems, approximate methods based
on perturbation expansions have been developed. In stellar dynamics, as-
trophysicists have developed numerous numerical and theoretical approaches
to the problem for larger values of N, including treatments based on the
Boltzmann equation and the Fokker-Planck equation; such N-body systems
can also be modelled as self-gravitating gases, and thermodynamic insights
underpin much of our qualitative understanding.
2 Few-Body Problems
2.1 The two-body problem
For N = 2 the relative motion of the two bodies can be reduced to the force-
free motion of the centre of mass and the problem of the relative motion. If
r = r1 − r2 , then
r
r̈ = −G(m1 + m2 ) 3 , (2)
|r|
often called the Kepler Problem. It represents motion of a particle of unit
mass under a central inverse-square force of attraction. Energy and angu-
lar momentum are constant, and the motion takes place in a plane passing
through the origin. Using plane polar coordinates (r, θ) in this plane, the
equations for the energy and angular momentum reduce to
!
1 2 L2 G(m1 + m2 )
E = ṙ + 2 − (3)
2 r r
L = r 2 θ̇. (4)
(Note that these are not the energy and angular momentum of the two-body
problem, even in the barycentric frame of the centre of mass; E and L must
be multiplied by the reduced mass m1 m2 /(m1 + m2 ).) Using eqs.(3) and (4)
the problem is reduced to quadratures. The solution shows that the motion
is on a conic section (ellipse, circle, straight line, parabola or hyperbola),
with the origin at one focus.
This reduction depends on the existence of integrals of the equations
of motion, and these in turn depend on symmetries of the underlying La-
grangian or Hamiltonian. Indeed eqs.(1) yield ten first integrals: six yield
the rectilinear motion of the centre of mass, three the total angular momen-
tum and one the energy. Furthermore, eq.(2) may be transformed, via the
2
Kustaanheimo-Stiefel transformation, to a four-dimensional simple harmonic
oscillator. This reveals further symmetries, corresponding to further invari-
ants: the three components of the Lenz vector. Another manifestation of the
abundance of symmetries of the Kepler problem is the fact that there exist
action-angle variables in which the Hamiltonian depends on only one action,
i.e. H = H(L). Another application of the KS transformation is one that
has practical importance: it removes the singularity of (i.e. regularises) the
Kepler problem at r = 0, which is troublesome numerically.
To illustrate the character of the KS transformation, we consider briefly
the planar case, which can be handled with a complex variable obeying the
equation of motion z̈ = −z/|z|3 (after scaling eq.(2)). By introducing the
Levi-Civita transformation z = Z 2 and Sundman’s transformation of the
time, i.e. dt/dτ = |z|, the equation of motion transforms to Z ′′ = hZ/2,
where h = |ż|2 /2 − 1/|z| is the constant of energy. The KS transformation is
a very similar exercise using quaternions.
where r is the position of the massless particle and Ω is the angular velocity
of the frame.
This problem has three degrees of freedom but only one known integral:
it is the Hamiltonian in the rotating frame, and equivalent to the Jacobi
integral, J. One consequence is that Liouville’s theorem is not applicable, and
more elaborate arguments are required to decide its integrability. Certainly,
no general analytical solution is known.
There are five equilibrium solutions, discovered by Euler and Lagrange
(see Fig.1). They lie at critical points of the effective potential in the rotating
frame, and demarcate possible regions of motion.
3
Throughout the twentieth century, much numerical effort was used in
finding and classifying periodic orbits, and in determining their stability and
bifurcations. For example there are families of periodic orbits close to each
primary; these are perturbed Kepler orbits, and are referred to as satellite
motions. Other important families are the series of Liapounov orbits starting
at the equilibrium points.
Some variants of the restricted three-body problem include
Then the exact solutions of Euler and Lagrange survive in the form of homo-
graphic solutions. In these solutions the configuration remains geometrically
4
similar, but may rotate and/or pulsate in the same way as in the two-body
problem.
Let us represent the position vector ri in the planar three-body problem
by the complex number zi . Then it is easy to see that we have a solution of
the form zi (t) = z(t)z0i , provided that
z
z̈ = −C
|z|3
and
mi Cz0i = ∇i W (z01 , z02 , z03 ),
for some constant C. Thus z(t) may take the form of any solution of the
Kepler problem, while the complex numbers z0i must correspond to what is
called a central √configuration. These are in fact critical points of the scale-
free function W I, where I (the “moment of inertia of the system”) is given
by I = 31 mi ri2 ; and C = −W/I.
P
2.3.2 Singularities
As Schubart’s solution illustrates, two-body encounters can occur in the
three-body problem. Such singularities can be regularised just as in the
pure two-body problem. Triple collisions can not be regularised in general,
and this singularity has been studied by the technique of “blow-up”. This
has been worked out most thoroughly in the collinear three-body problem,
which has only two degrees of freedom. The general idea is to transform
to two variables, of which one (denoted by r, say) determines the scale of
5
the system, while the other (s) determines the configuration (e.g. the ratio
of separations of the three masses). By scaling the corresponding velocities
and the time, one obtains a system of three equations of motion for s and
the two velocities which are perfectly regular in the limit r → 0. In this
limit, the energy integral restricts the solutions of the system to a manifold
(called the collision manifold). Exactly the same manifold results for zero-
energy solutions, which permits a simple visualisation. Equilibria on the
collision manifold correspond to the Lagrangian collinear solutions in which
the system either expands to infinity or contracts to a three-body collision.
d2 I
= 4T + 2W,
dt2
where T, W are, respectively, the kinetic and potential energies of the system.
Usually the barycentric frame is adopted. Since E = T + V is constant and
T ≥ 0, it follows that the system is not bounded for all t > 0 unless E < 0.
6
2.3.4 Perturbation theory
The question of the integrability of the general three body problem has stim-
ulated much research, including the famous study by Poincaré which estab-
lished the non-existence of integrals beyond the ten classical ones. Poincaré’s
work was an important landmark in the application to the three-body prob-
lem of perturbation methods. If one mass dominates, i.e. m1 ≫ m2 and
m1 ≫ m3 , then the motion of m2 and m3 relative to m1 is mildly perturbed
two-body motion, unless m2 and m3 are close together. Then it is beneficial
to describe the motion of m2 relative to m1 by the parameters of Keplerian
motion. These would be constant in the absence of m3 , and vary slowly be-
cause of the perturbation by m3 . This was the idea behind Lagrange’s very
general method of variation of parameters for solving systems of differential
equations. Numerous methods were developed for the iterative solution of the
resulting equations. In this way the solution of such a three-body problem
could be represented as a type of trigonometric series in which the arguments
are the angle variables describing the two approximate Keplerian motions.
These were of immense value in solving problems of celestial mechanics, i.e.
the study of the motions of planets, their satellites, comets and asteroids.
A major step forward was the introduction of Hamiltonian methods. A
three-body problem of the type we are considering has a Hamiltonian of the
form
H = H1 (L1 ) + H2 (L2 ) + R,
where Hi , i = 1, 2 are the Hamiltonians describing the interaction between mi
and m1 , and R is the “disturbing function”. It depends on all the variables,
but is small compared with the Hi . Now perturbation theory reduces to the
task of performing canonical transformations which simplify R as much as
possible.
Poincaré’s major contribution in this area was to show that the series
solutions produced by perturbation methods are not, in general, convergent,
but asymptotic. Thus they were of practical rather than theoretical value.
For example, nothing could be proved about the stability of the solar system
using perturbation methods. It took the further analytic development of
KAM theory to rescue this aspect of perturbation theory. This theory can be
used to show that, provided that two of the three masses are sufficiently small,
then for almost all initial conditions the motions remain close to Keplerian
for all time. Unfortunately now it is the practical aspect of the theory which
is missing; though we have introduced this topic in the context of the three-
7
body problem, it is extensible to any N-body system with N −1 small masses
in nearly-Keplerian motion about m1 , but to be applicable to the solar system
the masses of the planets would have to be ridiculously small.
8
3 Many-Body Problems
Many of the concepts already introduced, such as the virial theorem, apply
equally well to the many-body classical gravitational problem. In this section
we refer mainly to the new features which arise when N is not small. In
particular, statistical descriptions become central. The applications also have
a different emphasis, moving from problems of planetary dynamics (celestial
mechanics) to those of stellar dynamics. Typically, N lies in the range 102 –
1012 .
∂f ∂f ∂φ(r, t) ∂f
+ v. − . = 0 (6)
∂t ∂r ∂r ∂v Z
∇2 φ = 4πGm f (r, v, t)d3v, (7)
f = |E|7/2 , (8)
9
The solution just referred to is an example of an equilibrium solution.
In an equilibrium solution the virial theorem takes the form 4T + 2W = 0,
where T, W are appropriate mean-field approximations for the kinetic and
potential energy, respectively. It follows that E = −T , where E = T + V is
the total energy. An increase in E causes a decrease in T , which implies that
a self-gravitating N-body system exhibits a negative specific heat.
There is little to choose between one equilibrium solution and another,
except for their stability. In such an equilibrium, the bodies orbit within the
potential on a timescale of the crossing time, which is conventionally defined
GM 5/2
to be tcr = .
(2|E|)3/2
The most important evolutionary phenomenon of collisionless dynamics
is violent relaxation. If f is not time-independent then φ is time-dependent
in general. Also, from the equation of motion of one body, E varies according
to dE/dt = ∂φ/∂t, and so energy is exchanged between bodies, which leads
to an evolution of the distribution of energies. This process is known as
violent relaxation.
Two other relaxation processes are of importance:
10
timescale, tr , is much longer than any other timescale we have con-
sidered, this process leads to evolution of f (E), and it dominates the
long-term evolution of large N-body systems. It is usually modelled
by adding a collision term of Fokker-Planck type on the right side of
eq.(6).
In this case the only equilibrium solutions in a steady potential are
those in which f (E) ∝ exp(−βE), where β is a constant. Then eq.(7)
becomes Liouville’s equation, and for the case of spherical symmetry
the relevant solutions are those corresponding to the isothermal sphere.
where k is Boltzmann’s constant, and the integration is taken over all avail-
able phase-space. Their stability may be determined by evaluating the second
variation of S. It is found that it is negative definite, so that S is a local
maximum and the configuration is stable, only if ρ0 /ρe < 709 approximately.
A physical explanation for this is the following. In the limit when ρ0 /ρe ≃ 1
the self-gravity (which causes the spatial inhomogeneity) is weak, and the
system behaves like an ordinary perfect gas. When ρ0 /ρe ≫ 1, however,
the system is highly inhomogeneous, consisting of a core of low mass and
high density surrounded by an extensive halo of high mass and low density.
Consider a transfer of energy from the deep interior to the envelope. In the
envelope, which is restrained by the enclosure, the additional energy causes
a rise in temperature, but this is small, because of the huge mass of the
halo. Extraction of energy from around the core, however, causes the bodies
there to sink and accelerate; and, because of the negative specific heat of a
self-gravitating system, they gain more kinetic energy than they lost in the
original transfer. Now the system is hotter in the core than in the halo, and
the transfer of energy from the interior to the exterior is self-sustaining, in a
11
gravothermal runaway. The isothermal model with large density contrast is
therefore unstable.
The negative specific heat, and the lack of an equilibrium which maximises
the entropy, are two examples of the anomalous thermodynamic behaviour of
the self-gravitating N-body problem. They are related to the long-range na-
ture of the gravitational interaction, the importance of boundary terms, and
the non-extensivity of the energy. Another consequence is the inequivalence
of canonical and microcanonical ensembles.
12
3.4 Collisional evolution
Consider an isolated N-body system, which we suppose initially to be given
by a spherically symmetric equilibrium solution of eqs.(6-7), such as eq.(8).
The temperature decreases with increasing radius, and a gravothermal run-
away causes the “collapse” of the core, which reaches extremely high density
in finite time. (This collapse takes place on the long two-body relaxation
time scale, and so it is not the rapid collapse, on a free-fall time scale, which
the name rather suggests.)
At sufficiently high densities the time scale of three-body reactions be-
comes competitive. These create bound pairs, the excess energy being re-
moved by a third body. From the point of view of the one-particle distribution
function, f , these reactions are exothermic, causing an expansion and cool-
ing of the high-density central regions. This temperature inversion drives the
gravothermal runaway in reverse, and the core expands, until contact with
the cool envelope of the system restores a normal temperature profile. Core
collapse resumes once more, and leads to a chaotic sequence of expansions
and contractions, called gravothermal oscillations.
The monotonic addition of energy during the collapsed phases causes a
secular expansion of the system, and a general increase in all time scales. In
each relaxation time a small fraction of the masses escape, and eventually (it
is thought) the system consists of a dispersing collection of mutually unbound
single masses, binaries and (presumably) stable higher-order systems.
It is very remarkable that the long-term fate of the largest self-gravitating
N-body systems appears to be intimately linked with the three-body prob-
lem.
4 Further Reading
Arnold, V.I., Kozlov, V.V., Neishtadt, A.I., 1993, Mathematical Aspects of
Classical and Celestial Mechanics, 2nd ed, (Berlin: Springer)
Barrow-Green, J., 1997, Poincaré and the three body problem, (Providence:
American Mathematical Society)
Binney, J., Tremaine, S., 1988, Galactic Dynamics (Princeton: Princeton
University Press)
Chenciner A., Montgomery R., 2000, A Remarkable Periodic Solution of the
Three-Body Problem in the Case of Equal Masses, Ann. Math., 152, 881-901
13
Dauxois, T., et al (eds), 2003, Dynamics and Thermodynamics of Systems
with Long Range Interactions (Berlin: Springer-Verlag)
Diacu, F., 2000, A Century-Long Loop, Math. Intelligencer, 22, 2, 19-25
Heggie, D.C., Hut, P., 2003, The Gravitational Million-Body Problem (Cam-
bridge: Cambridge University Press)
Marchal, C., 1990, The Three-body Problem (Amsterdam: Elsevier)
Murray, C.D., Dermott, S.F., 2000, Solar System Dynamics, (Cambridge:
Cambridge University Press)
Padmanabhan, T., 1990, Statistical mechanics of gravitating systems, Phys.
Rept. 188, 285-362
Siegel, C.L., Moser, J.K., Kalme, C.I., 1995, Lectures on Celestial Mechanics
(Berlin: Springer-Verlag)
Steves, B.J., Maciejewski, A.J. (eds), 2001, The Restless Universe. Applica-
tions of Gravitational N-Body Dynamics to Planetary, Stellar and Galactic
Systems (Bristol: IoP and SUSSP)
Szebehely, V., 1967, Theory of orbits : the restricted problem of three bodies
(New York : Academic Press)
5 Figure Captions
14
Figure 1: The equilibrium solutions of the circular restricted three-body
problem. We choose a rotating frame of reference in which two particles are
at rest on the x-axis. The massless particle is at equilibrium at each of the
five points shown. Five similar configurations exist for the general three-body
problem; these are the “central” configurations.
15
Figure 2: A rare example of a scattering encounter between two binaries
(which approach from upper right and lower left) which leads to a perma-
nently bound triple system describing the “figure 8” periodic orbit. A fourth
body escapes at the bottom. Note the differing scales on the two axes. Orig-
inally published in MNRAS
16
0
-1
-2
log density
-3
-4
-5
-6
-1 -0.5 0 0.5 1 1.5 2 2.5 3
log radius
17
Figure 4: Gravothermal oscillations in an N-body system with N = 65536.
The
√ central density is plotted as a function of time, in units such that tcr =
2 2. Source: H. Baumgardt, P. Hut, J. Makino, with permission
18