Nbodie Problem

Download as pdf or txt
Download as pdf or txt
You are on page 1of 18

The Classical Gravitational N -Body Problem

arXiv:astro-ph/0503600v2 11 Aug 2005

Douglas C. Heggie
University of Edinburgh,
School of Mathematics,
King’s Buildings,
Edinburgh EH9 3JZ,
UK

Keywords: many-body problem - classical gravitation - Kepler problem -


three-body problem - few-body problem - periodic orbits - perturbation meth-
ods - numerical methods - Vlasov-Poisson system - Fokker-Planck equation
- relaxation - gravothermodynamics

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.

2.2 The restricted three-body problem


The simplest three-body problem is given by the motion of a test particle in
the gravitational field of two particles, of positive mass m1 , m2 , in circular
Keplerian motion. This is called the circular restricted three-body problem,
and the two massive particles are referred to as primaries. In a rotating frame
of reference, with origin at the centre of mass of these two particles, which
are at rest at positions r1 , r2, the equation of motion is
!
m1 m2
r̈ + 2Ω × ṙ + Ω × (Ω × r) = G∇ + , (5)
|r − r1 | |r − r2 |

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

1. the elliptic restricted three-body problem, in which the primaries move


on an elliptic Keplerian orbit; in suitable coordinates the equation
of motion closely resembles eq.(5), except for a factor on the right
side which depends explicitly on the independent variable (transformed
time); this system has no first integral.

2. Sitnikov’s problem, which is a special case of the elliptic problem, in


which m1 = m2 , and the motion of the massless particle is confined
to the axis of symmetry of the Keplerian motion; this is still non-
integrable, but simple enough to allow extensive analysis of such fun-
damental issues as integrability and stochasticity;

3. Hill’s problem, which is a scaled version suitable for examining motions


close to one primary; its importance in applications began with studies
of the motion of the moon, and it remains vital for understanding the
motion of asteroids.

2.3 The general three-body problem


2.3.1 Exact solutions
When all three particles have non-zero masses, the equations of motion be-
come
mi r̈i = −∇i W,
where the potential energy is
X mi mj
W = −G .
1≤i<j≤3 |ri − rj |

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

The existence of other important classes of periodic solutions can be


proved analytically, even though it is not possible to express the solution
in closed form. Examples include hierarchical three-body systems, in which
two masses m1 , m2 exhibit nearly elliptic relative motion, while a third mass
orbits the barycentre of m1 and m2 in another nearly elliptic orbit. In the
mathematical literature this is referred to as motion of elliptic-elliptic type.
More surprisingly, the existence of a periodic solution in which the three
bodies travel in succession along the same path, shaped like a figure 8 (cf.
fig2), was established by Chenciner & Montgomery (2000), following its inde-
pendent discovery by Moore using numerical methods. Another interesting
periodic motion that was discovered numerically, by Schubart, is a solution
of the collinear three-body problem, and so collisions are inevitable. In this
motion the body in the middle alternately encounters the other two bodies.

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.

2.3.3 Qualitative ideas


Reference has already been made to motion of elliptic-elliptic type. In motion
of elliptic-hyperbolic type there is again an “inner” pair of bodies describing
nearly Keplerian motion, while the relative motion of the third body is nearly
hyperbolic. In applications this is referred to as a kind of scattering encounter
between a binary and a third body. When the encounter is sufficiently close,
it is possible for one member of the binary to be exchanged with the third
body. One of the major historical themes of the general three-body problem is
the classification of connections between these different types of asymptotic
motion. It is possible to show, for instance, that the measure of initial
conditions of hyperbolic-elliptic type leading asymptotically to elliptic-elliptic
motion (or any other type of permanently bound motion) is zero. Much of
the study of such problems has been carried out numerically.
There are many ways in which the stability of three-body motions may be
approached. One example is furnished by the central configurations already
referred to. They can be used to establish sufficient conditions for ensuring
that exchange is impossible, and similar conclusions.
A powerful tool for qualitative study of three-body motions is Lagrange’s
identity, which is now thought of as the reduction to three bodies of the
virial theorem. Let the size of the system be characterised by the “moment
of inertia” I. Then it is easy to show that

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.

2.3.5 Numerical methods


Numerical integrations of the three-body problem were first carried out near
the beginning of the 20th century, and are now commonplace. For typical
scattering events, or other short-lived solutions, there is usually little need to
go beyond common Runge-Kutta methods, provided that automatic step-size
control is adopted. When close two-body approaches occur, some regularisa-
tion based on the KS transformation is often exploited. In cases of prolonged
elliptic-elliptic motion, an analytic approximation based on Keplerian motion
may be adequate. Otherwise (as in problems of planetary motion, where the
evolution takes place on an extremely long time scale), methods of very high
order are often used. Symplectic methods, which have been developed in the
context of Hamiltonian problems, are increasingly adopted for problems of
this kind, as their long-term error behaviour is generally much superior to
that of methods which ignore the geometrical properties of the equations of
motion.

2.4 Four- and five-body problems


Many of the foregoing remarks, on central configurations, numerical methods,
KAM theory, etc, apply equally to few-body problems with N > 3. Of special
interest from a theoretical point of view is the occurrence of a new kind of
singularity, in which motions become unbounded in finite time. For N = 4
the known examples also require two-body collisions, but non-collision orbits
exhibiting finite-time blow-up are known for N = 5.
One of the practical (or, at least, astronomical) applications is again to
scattering encounters, this time involving the approach of two binaries on
a hyperbolic relative orbit. Numerical results show that a wide variety of
outcomes is possible, including even the creation of the figure-8 periodic
orbit of the three-body problem, while a fourth body escapes (Fig.2).

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 .

3.1 Evolution of the distribution function


The most useful statistical description is obtained if we neglect correlations
and focus on the one-particle distribution function f (r, v, t), which can be
interpreted as the number-density at time t at the point in phase space
corresponding to position r and velocity v. Several processes contribute to
the evolution of f .

3.1.1 Collective effects


When the effects of near neighbours are neglected, the dynamics is described
by the Vlasov-Poisson system

∂f ∂f ∂φ(r, t) ∂f
+ v. − . = 0 (6)
∂t ∂r ∂r ∂v Z
∇2 φ = 4πGm f (r, v, t)d3v, (7)

in which φ is the gravitational potential, and m is the mass of each body.


Obvious extensions are necessary if not all bodies have the same mass.
Solutions of eq.(6) may be found by the method of characteristics, which
is most useful in cases where the equation of motion r̈ = −∇φ is integrable,
e.g. in stationary, spherical potentials. An example is the solution

f = |E|7/2 , (8)

where E is the specific energy of a body, i.e. E = v 2 /2 + φ. This satisfies


eq(6) provided that φ is static. Eq.(7) is satisfied provided that φ satisfies a
case of the Lane-Emden equation, which is easy to solve in this case.

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:

1. Relaxation is possible on each energy hypersurface, even in a static


potential, if the potential is non-integrable.

2. The range of collective phenomena becomes remarkably rich if the sys-


tem exhibits ordered motions, as in rotating systems. Then an impor-
tant role is played by resonant motions, especially resonances of low
order. The corresponding theory lies at the basis of the theory of spiral
structure in galaxies, for instance.

3.1.2 Collisional effects


The approximations of collisionless stellar dynamics suppress two important
processes:

1. the exponential divergence of stellar orbits, which takes place on a time


scale of order tcr . Even in an integrable potential, therefore, f evolves
on each energy hypersurface.
N
2. two-body relaxation. It operates on a time scale of order tcr ,
ln N
where N is the number of particles. Though this two-body relaxation

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.

3.2 Collisional equilibrium


We consider the collisional evolution of an N-body system further in Sec.3.4,
and here develop the fundamental ideas about the isothermal model. The
isothermal model has infinite mass, and much has been learned by considering
a model confined within an adiabatic boundary or enclosure. There is a series
of such models, characterised by a single dimensionless parameter, which can
be taken to be the ratio between the central density and the density at the
boundary, ρ0 /ρe (Fig.3).
These models are extrema of the Boltzmann entropy S = −k f ln f d6 τ ,
R

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.

3.3 Numerical methods


The foregoing considerations are difficult to extend to systems without a
boundary, though they are a vital guide to the behaviour even in this case.
Our knowledge of such systems is due largely to numerical experiments, which
fall into several classes:

1. Direct N-body calculations. These minimise the number of simplifying


assumptions, but are expensive. Special-purpose hardware is readily
available which greatly accelerates the necessary calculations. Great
care has to be taken in the treatment of few-body configurations, which
otherwise consume almost all resources.

2. Hierarchical methods, including tree methods, which shorten the cal-


culation of forces by grouping distant masses. They are mostly used
for collisionless problems.

3. Grid-based methods, which are used for collisionless problems

4. Fokker-Planck methods, which usually require a theoretical knowledge


of the statistical effects of two-, three- and four-body interactions. Oth-
erwise they can be very flexible, especially in the form of Monte Carlo
codes.

5. Gas codes. The behaviour of a self-gravitating system is simulated


surprisingly well by modelling it as a self-gravitating perfect gas, rather
like a star.

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

Figure 3: The density profile of the non-singular isothermal model, with


conventional scalings.

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

You might also like