Publ. Natl. Astron. Obs. Japan, Vol. 9, 55–112 (2007)

Trends in 20th Century Celestial Mechanics
Takashi Ito and Kiyotaka Tanikawa(Received April 28, 2006; Accepted March 14, 2007)
Abstract
We review some of the major achievements of celestial mechanics in the twentieth century, and discuss someunsolved problems that are left to the twenty-first century. The four major research fields in celestial mechanics aretreated: dynamical systems, three-body problems, solar system dynamics, and numerical methods. Over the past
decades, celestial mechanics has become extremely specialized and categorized along many trajectories, from purely
mathematical to quite practical points of view, as well as in what has involved digital computer innovation. Thesevarious ends of celestial mechanics sometimes exist in isolation, losing their connections with each other. The purposeof this manuscript is to find clues for rediscovering connections between each of these fields by listing some of majorachievements in hoping to find new paths for celestial mechanics in the new century.
Key words: celestial mechanics, dynamical system, three-body problem, planetary dynamics, numerical method

Contents
1
Introduction
2
Dynamical systems—from Poincaré to chaos
2.1 The KAM theory
2.2 Chaos
2.3 Twist maps
2.4 Integrability
3
The three-body problem
3.1 Analytical and qualitative studies
3.2 Numerical studies
3.3 Hierarchical triple systems
3.4 New frontiers of the three-body problems
3.5 Applications to actual planetary and stellar systems
4
Solar system dynamics
4.1 Progress of perturbation theories
4.2 Long-term dynamical stability of planetary motion
4.3 Dynamics in planet formation study
4.4 Dynamics of small bodies
4.4.1 Asteroids and resonant dynamics
4.4.2 The Yarkovsky effect and asteroid families
4.4.3 Kuiper Belt objects
4.4.4 The Oort Cloud
4.5 Planetary rotation
4.5.1 Nutation and precession of the Earth
4.5.2 The Earth–Moon system
4.5.3 Evolution of planetary rotations
5
Numerical methods
5.1 Regularization
5.2 Symplectic integrators
5.3 Other promising numerical techniques
5.4 GRAPE: A special-purpose computer for N -bodydynamics
6
Concluding Remarks
Acknowledgments
References

55
56
58
60
63
65
68
68
69
70
71
73
73
74
76
78
79
79
83
84
85
86
86
87
87
88
88
89
91
93
93
95
95

1. Introduction
This manuscript represents our effort to summarize
some of the major trends in celestial mechanics in the
Corresponding author: ito.t@nao.ac.jp (Takashi Ito).

recently ended twentieth century. Not only celestial mechanics but the entire field of astronomy has evolved anddrastically changed over the past one hundred years. To
commemorate the beginning of the new century, we venture this review of the history of celestial mechanics in thepast century in order to identify new paths to explore in theforthcoming years. It goes without saying that we possessneither sufficient knowledge nor experience to thoroughly
review the celestial mechanics of the past century. Neither
can we predict its way in the new century. What is written
in this manuscript is our recognition and anticipation of
the past and the future of several specific areas of celestialmechanics in which we are interested. In this sense, ourdescription may be biased toward the authors’ favor. Long
before our attempt, V. Szebehely (1997) stated Open problems on the eve of the next millennium. Though the scopeof our manuscript is not as large as the next millennium,
we hope our attempt serves as a milestone for getting clues
and opening new frontiers in celestial mechanics.
Major flows in nineteenth-century celestial mechanics
converged at H. Poincaré. And major flows in twentieth century celestial mechanics diverged out from G. D. Birkhoff.
The stability analysis by A. M. Lyapunov, which he started
in the late nineteenth century, also grew into a large flowin the early twentieth century. These three great masters
established many of the fundamental concepts of modern
celestial mechanics—integrability, convergence of perturbative developments, normal forms, stability of dynamicalsystems, and so on. In the 1940s, C. L. Siegel worked on the
singularity of the three-body and small denominator problems. Based on what these researchers had established, theresearch of dynamical systems in the 1950s and the 1960s
produced the greatest achievement of celestial mechanics
in the twentieth century: the Kolmogorov–Arnold–Moser(KAM) theory. Later this research flow morphed into theinvestigation of chaos phenomena, thanks to the development of computer technology.

56

Takashi Ito and Kiyotaka Tanikawa

In the first half of the twentieth century, after Poincaréand Birkhoff, celestial mechanics flourished in a variety offorms. The understanding of planetary motion in our solar
system progressed, led by the development of accurate perturbation theories. The progress of planetary perturbation
theories stimulated the search for the ninth planet, resulting in the discovery of Pluto by C. W. Tombaugh in 1930.
Research on major and minor planets alike became one of
the central issues of celestial mechanics in this era as astrometric observation became increasingly precise. Motivatedby the identification of asteroid families by K. Hirayama in
1918, the quantitative research of asteroid dynamics began.
The Soviet Union’s launch of the world’s first artificial
satellite, Sputnik in 1957, turned out to be one of the mostsignificant stimuli for celestial mechanics. This epochmaking event marked the start of the space age, and owed muchto the development of advanced technology. From the viewpoint of celestial mechanics, Sputnik represented the triumph of Newtonian mechanics and perturbation theories.
Theoretical studies of the motion of artificial satellites madegreat advances after Sputnik’s launch.
As well as to many other fields of science and engineering, the development of fast computers has made an enormous contribution to celestial mechanics, replacing the numerical studies by hands (“human” computers) that had
already started in the early twentieth century. The use ofelectronic digital computers for numerical studies in celestial mechanics began in the 1950s. At first, the generalthree-body problem was targeted as an example study using electronic computers. Regularization methods that had
been initially devised for theoretical purposes have become
important since the 1960s. The techniques of regularization
partially overcame one of the biggest difficulties in celestialmechanics: collisions. Computers are also used as instruments for purely theoretical research. For example, developments of dynamical theories were enhanced a great deal
owing to the detailed numerical survey of phase space.
The progress of numerical study in celestial mechanics
owes a lot not only to the development of computer hardware, but also to sophisticated numerical algorithms. A variety of numerical integration schemes such as the Runge–
Kutta methods, various predictor–correctors, the extrapolation method, and the family of symmetric integration
schemes has contributed greatly to the progress of numerical celestial mechanics. In the 1990s, the knowledge of thesenumerical techniques was combined with that of Hamiltonian dynamics, resulting in a totally new type of numerical
scheme: symplectic integrators.
Digital computer technology yielded huge innovations
not only in the field of theoretical computation but also in
the area of advanced engineering of astronomical observation that has given celestial mechanics a substantial stimulation. For example, detailed observations of solar systemobjects such as asteroids, comets, and the trans-Neptunian
objects including the Kuiper Belt and the Oort Cloud ob-

jects have helped solar system dynamics achieve a hugeprogress. A number of spacecrafts with cutting-edge digitaltechnology have explored the solar system objects, giving
us precious information about the dynamical nature of ourplanetary system. In addition to this, the recent rush of discoveries of extrasolar planetary systems, owing to the development of new observing technologies, has significantly extended the field of celestial mechanics since the mid 1990s.
Bearing this background in mind, from the next section
we review some of the major achievements and remaining problems of twentieth century celestial mechanics 1 .
We have categorized the major fields of celestial mechanics
into four parts: dynamical systems (Section 2), the threebody problem (Section 3), solar system dynamics (Section4), and numerical methods (Section 5). Section 6 presents
some discussion and concluding remarks. We are fully awareof our inability to describe several important topics suchas relativistic celestial mechanics, boundary area to stellar
dynamics, current progress of artificial satellite theories,
observation-related research of near-Earth objects, detailed
description about ring-satellite dynamics, and many otherfields. This ignorance just reflects the limitation of our abilities and experience in these fields. We have noted major
missing topics and their importance at the beginning ofeach section. Throughout this manuscript, we would like toprovide some thoughts and questions for our readers: Whatis celestial mechanics? How should we define it? Since theprogress of celestial mechanics is quite rapid, the answer
to these questions might not be straightforward. But thesubtleties of such questions could help us to delineate newacademic frontiers in the new century.
2. Dynamical systems—from Poincaré to chaos
Towards the end of the nineteenth century, H. Poincarérevolutionized the human view of nature (Poincaré, 1892).
Nature is not so simple as the simple form of the Newtonian equations of motion suggests, something Poincarérealized through his own experience. His first experience
perhaps was his failure to show the solvability of the restricted three-body problem in the Grand prix paper for
the Swedish prize in 1889. The failure was said to be related
to the stable and unstable manifolds of periodic points. At
first, Poincaré implicitly assumed that these manifolds wereconnected smoothly. But later, he realized that they generally intersect transversely. His second experience came withhis success in proving the analytic non-integrability of therestricted three-body problem in 1890. His failures and successes might have made Poincaré a modest person. Afterthese experiences, Poincaré proposed to study the general
1

Note that our definition of the border between the twentieth andtwenty-first century is somewhat ambiguous. Sometimes we consider
the publications from 2001 to today as those of the “previous century”. Basically we regard research undertaken now and in the near
future as that of the “new century”.

Trends in 20th Century Celestial Mechanics
problem of dynamics with a Hamiltonian F , developed interms of the powers of a small parameter µ in the following
manner:F = F0 (p) + µF1 (p, q) + µ2 F2 (p, q) + · · · ,(1)where p and q are the action–angle variables, F0 is the unperturbed part of F depending only on p, and F1 , F2 , . . .
are periodic functions of period 2π with respect to q. This
problem became known as the fundamental problem of dynamics. The moral of the fundamental problem is as follows: Do not swim far from safe land (i.e. the integrable
system); if one ignores this warning, one may be swept
up in the retreating flow and drowned in the vast sea ofthe unknown. However, an infinite sequence of perturbative transformations and an infinite series of perturbative
developments do not converge due to the existence of small
denominators. Any perturbation theory would not yield exact solutions, no matter how high its order is. If an infinite
series developed in powers of small parameters converges,
it would mean the dynamical system under consideration
is integrable. A common understanding at the end of thenineteenth century was that the general problems of celestial mechanics are non-integrable. Integrable systems suchas a special spinning top are very rare.
Poincaré may have sensed the limitations of differential equations, or the difficulties of solving differential equations. Newton’s laws are differential, whereas Kepler’s laws
are integral. Differential laws describe the relations in nature locally; they connect the relations of physical quantities spatially and temporally in the close vicinity. On theother hand, Kepler’s laws describe the global relations between physical quantities. Natural laws might be describedusing procedures other than differential equations. Nevertheless, according to Daniel L. Goroff who is the editor of
the English translation of Poincaré (1892), Poincaré appreciated the great conceptual power of differential equations.
According to Poincaré, thanks to the postulate that theactual state of the world depends only on the most recent
past, without being directly influenced by the memory ofthe distant past, one can start from writing out differential
equations.
Poincaré kept his distance from pure differential equations; he undertook the qualitative and geometrical studyof dynamical systems. Poincaré surfaces and Poincaré maps
are typical examples of this direction in his research. He
was interested in periodic orbits as phenomena forming the
backbone of phase space of dynamical systems. He founddoubly asymptotic solutions. These are generally transverse intersections of singly asymptotic manifolds of solutions to periodic solutions, solutions which Poincaré called“homoclinic” or “heteroclinic”. These solutions now play
central roles in the study of dynamical systems. He wasalso interested in the stability of the solar system. Much of
Poincaré’s work other than in the field of celestial mechanics, such as on triangulation, homology, torsion, and the57

fundamental group in algebraic topology, was “motivatedby Poincaré’s desire to extend his ideas about dynamics tohigher dimensions,” again according to Daniel L. Goroff.
Virtually independently in Russia, motivated by thestability problem of planetary motion, A. M. Lyapunovdeveloped a stability theory of dynamical systems (Lyapunov, 1892). Lyapunov’s first method is related to the stability analysis of equilibrium points of non-linear systems,
in which equilibrium points are linearly asymptotically stable. But when the linear stability of the equilibrium point
is marginal, his first method is not applicable. So, Lyapunov proposed his second method for this case. One looks
for a Lyapunov function which is positive and whose time
derivative is zero or negative. The equilibrium point is stable if the time derivative of the Lyapunov function is zero,
while the point is asymptotically stable if the derivative isnegative. Lyapunov’s results had a profound influence on
the celestial mechanics of the twentieth century. We recognize this fact by referring to the notions named after Lyapunov, such as the Lyapunov characteristic exponent, theLyapunov function, Lyapunov’s first and second methods,
the Lyapunov stability, and the Lyapunov asymptotic stability. Actually, equilibrium points in Hamiltonian systems
are always marginal, and the time derivative of the desired
Lyapunov function is zero, i.e., the Lyapunov function is
the integral of motion. Thus, the search for Lyapunov functions is equivalent to the search for new integrals of motion.
G. D. Birkhoff began his substantial scientific career by
proving Poincaré’s last geometric theorem on the existenceof fixed points for twist maps. Then following Poincaré,
Birkhoff qualitatively studied the motion of the third bodyin the restricted three-body problem (Birkhoff, 1915). He
returned to this subject several times in the course of hislife. He introduced various concepts in general dynamicalsystems that are extensively used at present, such as α- andω-limit points, minimal motions, central motions, wandering and non-wandering points. Formally, he obtained thenormal form of great generality for Hamiltonian and Pfaffian systems (Morse, 1946). Though many of these notions
are valid for an arbitrary degree of freedom, Birkhoff’s
analysis was restricted to systems with two degrees of
freedom. The analysis of two-dimensional area-preserving
maps, which was started by Poincaré and deepened by
Birkhoff, is extremely popular at present. Even when wecannot follow the individual orbits of celestial bodies, weare able to follow a group of orbits, or a bunch of orbits.
This was Birkhoff’s basic idea which commemorated thebirth of the research into dynamical systems.
Many researchers have followed the legacy of Poincaré,
Birkhoff and Lyapunov. C. L. Siegel was one of them. His
most important contribution to the study of dynamicalsystems is twofold: one aspect is the solution of the small
denominator problem in a simpler form. The second is
the proof of general divergence of the normal form procedure proposed by Birkhoff (Siegel and Moser, 1971). The

58

Takashi Ito and Kiyotaka Tanikawa

KAM theory (see Section 2.1) started from a brief note by
A. N. Kolmogorov (1954) where the proof of the KAM theorem was only outlined with the use of Newton’s method
to obtain zeros of a function. V. I. Arnold (1963) proved
the theorem for the analytical Hamiltonian systems (seeSection 2.1), and J. K. Moser (1962) proved the theorem
for smooth 2-dimensional twist maps (see Section 2.3).
S. Smale (1967) devised the so-called horseshoe in an illustrative manner. M. Hénon (1976) proposed a map, now
called the “Hénon map” which may exhibit strange attractors. Benedicks and Carleson (1991) proved that the Hénon
map has a strange attractor for some parameters. Strange
attractors have been found numerically by Lorenz (1963)
and Ueda et al. (1973), eventually leading to the world of
chaos.
In the following subsections, we describe major four issues that derived from the work of Poincaré, Lyapunov, andBirkhoff. The first of these is the KAM theory including
the Arnold diffusion and the Aubry–Mather sets. The second is the world of chaos. The third is the area-preserving
twist maps derived from Hamiltonian systems. The final
issue is integrability. Major important topics of this fieldthat are not described in this section include integrability of partial differential equations and advance of canonical perturbation theories (though we will briefly mention
some perturbation theories in relation to the solar system
dynamics in Section 4.1). Readers can study more about
these topics from Brouwer and Clemence (1961), Hagihara(1970, 1972a,b), Novikov (1981), Soffel (1989), Brumberg(1991), Arnold et al. (1993), or Boccaletti and Pucacco(1996, 1998).

that the future is stochastic: the longevity of the solar system (see also Section 4.2). The lifetime of the solar systemthat we know has grown along with our understanding of
the universe. In the nineteenth century, our universe was
supposed to be only several tens of million years old. Evenin the 1920s, the age of the universe was considered to be
as young as two billion years. In the 1950s, the age of thesolar system turned out to be about 4.5 billion years oldfrom the results of radioisotope analysis of meteorites. This
finally indicated the possibility that planetary motion has
been very stable for billions of years since its formation,
maintaining a nearly regular (in the terminology of Hamiltonian dynamics, “quasi-periodic”) motion. On the other
hand, the equations of motion describing the solar system
are non-integrable, and planetary motions can be, and actually are, chaotic (see Section 2.2). If the motion of theplanets has been bounded in a small region of the phasespace for billions of years, we may expect that there exist
barriers that prevent planets from wandering out of theiraccustomed place. Is this expectation true?
Here we have the very theory that answers precisely tothis question: the KAM theory (Kolmogorov, 1954; Arnold,
1963; Moser, 1962). According to the KAM theory, when
non-integrability is weak, almost all of phase space is occupied by such barriers, i.e., invariant sets of points that arefilled with quasi-periodic orbits. As the system goes further from integrable state, these barriers become sparse.
These barriers are formed by quasi-periodic motions. Here,a quasi-periodic motion is a motion described by a quasiperiodic function of time:

2.1. The KAM theory

with frequencies ω1 , . . . , ωk that are rationally independent.
All other motions in phase space are sandwiched by thesequasi-periodic motions. This is topologically true in two
degrees of freedom case since quasi-periodic motions form
two-dimensional torus which divides three-dimensional energy surface into two parts. The solar system planetary
motion might not be described by (2). However, planetary
orbits will not go beyond the barriers whose existence is assured by the KAM theory. Again, the pendulum swung tothe opposite side. “Okay, our future is not so uncertain.”Not only to celestial mechanics, but to nearly all fields ofmodern mathematical sciences, the KAM theory was a relief.
Let us briefly explain the proof of the KAM theorem
along with Arnold (1963). In his paper, after showing the
preparatory results from celestial mechanics and mathematics, Arnold explains the classical perturbation theoryof celestial mechanics. The perturbation theory has two
essential difficulties: One is the so-called small denominator problem which comes from resonance and whose existence is dense in the frequency domain. The other is thedivergence of the infinite sequence of canonical transformations to drive off the perturbation to higher and higher

Pierre Simon de Laplace, who lived in the latter halfof the eighteenth century and in the first quarter of thenineteenth century, imagined the so-called “Daemon of
Laplace” in his book (Laplace, 1814). Literally, Laplace
wrote “Give me all the initial conditions. Then I will
predict as precisely as possible the future of the world,”
expressing the victory of Newtonian mechanics. Celestialmechanists, scientists, and the general audience of thenineteenth century anticipated or wished the world to beprecisely predictable. But at the same time, strange phenomena which could not be explained by simple mechanics
were also being discussed—thermal and electromagnetic
phenomena. These gradually became the subjects of modern sciences. Confidence in the mechanical world was
shaken like an unstable sky scraper. Eventually, Poincaré
proved that perturbative expansion does not converge in
general. Almost at the same time, came the emergence ofstatistical mechanics. The pendulum of the world seemedto swing toward extreme ambiguity, and suddenly, thefuture looked statistically uncertain.
However, we have strong counterevidence to the ideax(t) = f (ω1 t, . . . , ωk t),(2)

Trends in 20th Century Celestial Mechanics
orders. Then, Arnold introduces a program proposed byKolmogorov (1954) to circumvent these difficulties, as wellas to show the stability of fixed points and periodic orbitsof near integrable systems.
Kolmogorov’s (1954) program is twofold: In order toovercome the two difficulties that we mentioned above, atfirst Kolmogorov recommends that we be satisfied withlooking at particular frequencies, and suggests that we look
for tori where the quasi-periodic motions with these frequencies persist. Kolmogorov’s second recommendation isto restrict the transformation to the neighborhood of suchtori, avoiding resonant regions. As a result, new perturbations will be quadratically small compared with the starting perturbation. Repeating this transformation, we willobtain a smaller region of applicability of the transformation in phase space. Finally we will arrive at invariant tori.
The convergence will be very fast, like popular Newton’s
method for obtaining the zeros of functions.
Arnold followed Kolmogorov’s suggestion and carriedout the program. Actually, Arnold did more. He was interested in the stability of the planetary system. In this case,the situation is not so simple: He needed to treat a degenerate problem. The integrable Hamiltonian is “properly
degenerate”, that is, the number of angle variables or thenumber of independent frequencies is less than the numberof degrees of freedom. In addition, perturbed systems havefast and slow variables, i.e. two different time scales. The
variables representing the orbital period are the fast ones,
whereas the variables that do not appear in the integrable
case are the slow ones. A typical example of fast variablesin the Keplerian motion is mean anomaly, and a typical example of slow variables is argument of perihelion. In order
to treat them, Arnold made precise the notion of adiabatic
invariant. He divided the variables into fast and slow ones
and carried out averaging over the fast variables, which isa standard technique in modern celestial mechanics.
Arnold’s system has the Hamiltonian H with n degrees
of freedom:H(p, q) = H0 (p0 ) + µH1 (p, q) + O(µ2 ),(3)

where H0 represents the integrable part, µ a small parameter, p and q canonical momenta and coordinates, part ofwhich being action–angle variables. The last term on theright-hand side of (3) is the higher-order terms in µ. We
should notice that p0 is the action variables corresponding
to the fast angular variables, and its dimension is less thann. Hence the system is properly degenerate. However, as
an independent Hamiltonian, H0 (p0 ) is assumed to be nondegenerate. H1 (p, q) is assumed to consist of three termsof different characters representing perturbations. The analytic Hamiltonian H(p, q) is defined explicitly in a complex domain. This is essential in the following discussion.
The assertion is that the real part of the domain of the definition of the Hamiltonian is divided into two sets: One ofthem is invariant with respect to the canonical equations

59

with this Hamiltonian, and the other is small. The larger
set consists of invariant n-dimensional analytic tori which
are little deformed from the integrable ones. Motions onthese tori are quasi-periodic.
In order to prove this fundamental theorem, Arnold
first carried out the integration of the Hamiltonian overthe fast variables, changing slow variables into the polar
coordinates and reducing the original Hamiltonian to theform
H(p, q) = H0 (p0 ) + H1 (p0 , p1 , q1 ) + H2 (p0 , p1 , q0 , q1 ), (4)where dim p0 = n0 and dim p1 = n1 with n0 + n1 = n, andthe suffix of H indicates the order of the term. H is definedin a complex domain F as
F = {p0 ∈ G0 , |Im q0 | ≤ ρ, p21 + q12 ≤ R}.

(5)

where G0 is a domain, and ρ and R are positive numbers.
Then, Arnold constructed an inductive scheme with an inductive theorem, an inductive lemma, and a fundamental
lemma accompanied by a series of technical lemmas. He applied this scheme to the Hamiltonian (4). His purpose was,through an infinite sequence of transformations, to nullify
the H2 term. He showed in the fundamental lemma that
there is a canonical transformation such that in the new
Hamiltonian, the new H2 is quadratically small comparedwith the original H2 , though the domain of the Hamiltonian
becomes smaller with certain resonant domains excluded
from the original domain. Using the fundamental lemma,
he proved the inductive lemma that roughly states that the
Hamiltonian system of the form of (4), which is slightly different in Arnold’s original paper, can be transformed into
a new Hamiltonian system. Arnold also proved that the
new system possesses a form similar to that of the original Hamiltonian with a quadratically small H2 term andwith a reduced domain. Then in the inductive theorem, herepeated an arbitrary number of times the transformation
defined in the inductive lemma, and reduced the absolutevalue of the H2 term to an arbitrarily high order. And finally, Arnold went to the limit and proved that the phasespace is divided into two parts: One of them occupies almost the entire phase space and is filled with tori, and theother occupies a small set in the sense of measure.
Some serious questions remained open after the proofof the KAM theorem. A big one was that of the higher dimensional case. n-dimensional tori of quasi-periodic motion
do not divide the 2n-dimensional phase space. Therefore,
tori are not the barriers in the higher dimensional cases.
General motion between two tori is not confined, and thesolution gradually diffuses. This phenomenon is called theArnold diffusion (Arnold, 1964; Nekhoroshev, 1977). How
fast is the diffusion, then? Could the solar system disintegrate soon? Fortunately, the diffusion speed is slow. The diffusion proceeds along the web of stable and unstable manifolds from various periodic orbits. For the Hamiltonian

60

Takashi Ito and Kiyotaka Tanikawa

H0 (p) + εH1 (p, q),(6)

with action p and angle q with a small parameter ε, there
exist constants a, b such that forµ ¶1
1
0 < t < T = exp a ,(7)εεthe action satisfies|p(t) − p(0)| < εb .

(8)

The diffusion velocity should depend on the distanced from the invariant torus. More precisely (Morbidelli and
Giorgilli, 1995), the diffusion is super-exponentially slow(∼ 1/ exp[exp(1/d)]) until some threshold distance. The
diffusion is exponentially slow (∼ 1/ exp(1/d)) until thenext threshold distance. And, it is quadratically slow (∼ d2 )until the final threshold distance (Morbidelli and Guzzo,1997). The last threshold is the boundary to the sea ofstrong chaos.
An ideal concept like the Arnold diffusion and its confirmation are always two different things. Xia (1992) claimed
to have confirmed the existence of the Arnold diffusion in
the planar isosceles three-body problem. However, in a realistic dynamical system with a large number of degrees of
freedom, nobody can say for sure that he/she has observed
the Arnold diffusion yet. For us to “see” the Arnold diffusion, we need to take a special care to somehow shorten its
timescale. Minor bodies in our solar system might give us
a clue to “see” the Arnold diffusion, because the dynamical lifetime of some minor objects could be so short thatthe timescale of their Arnold diffusion might be observablein numerical experiments. Some such trials have alreadystarted. For example, Morbidelli (1997) proposed using asteroids with short dynamical lifetimes in order to observe
the Arnold diffusion in numerical experiments. Morbidelliand Guzzo (1997) explored the dynamical structure of the
2:3 mean motion resonance with Neptune in the Kuiper
Belt objects using numerical integrations. They focused on
the existence of slowly diffusing chaotic orbits that escapefrom the 2:3 resonance after billions of years, which might
be a part of the Arnold diffusion. In any case, we need avery large amount of computation with a very high accuracy (hopefully with a large-number of arithmetic digits)
to observe the Arnold diffusion within a realistic timescale.
Before the end of this subsection, let us give a remark
on an interesting application of Nekhoroshev’s spirit to thestability of asteroid orbits (Giorgilli and Skokos, 1997). This
work shows that the Nekhoroshev-type estimate is applicable not only to nearly integrable systems but also to systems near an equilibrium using the Birkhoff normal form.
The essence of the idea is contained in the trivial inequality|pj (t) − pj (0)| ≤ |t|

supdomain ∆|ṗj |(9)

where pj is the action variable under consideration, t is
time, ṗj is the time derivative of pj , and ∆ is a domain

where the diffusion takes place. We transform the expanded
Hamiltonian to a normal form up to an optimal order r as
H=

rX

Hk + Rr+1

(10)

k=2

where Hr is a function only of pj (j = 1, . . . , n) of order r/2,
and Rr+1 is a remainder of higher orders. Then, we obtainṗj = {pj , H} = {pj , Rr+1 }.

(11)

Giorgilli and Skokos (1997) succeeded in estimating thesupremum of the norm ||{pj , Rr+1 }|| in a neighborhood of
L4 in the restricted three-body problem, and discussed thestability of actual Jupiter’s Trojan asteroids.
What happens if quasi-periodic motions no longer constitute a surface? How would it look like when the torus
of quasi-periodic motion gets disintegrated? This question
had been open since the dawn of KAM theory research, andwas solved for two-dimensional twist maps by Aubry and
LeDaeron (1983) and Mather (1982). A KAM curve, when
it disintegrates like a sweater with its weft lost, becomesloose and wavy. When we look at the surface of section of
the KAM curve, a continuous curve is disintegrated into a
Cantor-like discontinuous set: an Aubry–Mather set. Now
we are in the world of chaos. We will come back to the
Aubry–Mather theory in the section on twist maps (Section 2.3).
2.2. Chaos
Chaos was recognized by the ancient Chinese philosopher, Zhuangzi, in the 4th century BC. Zhuangzi realizedthat chaos was essential in the world, and that the majorpart of the world could become meaningless, dark, or evendead without chaos. Mathematical chaos has its origin in
the work of Poincaré. Poincaré carefully observed special
dynamical flows going into and emanating from hyperbolic
fixed points of the restricted three-body problem. He also
observed complicated tangles of stable and unstable manifolds. Poincaré shuddered at the glimpse of infinity that
Newtonian dynamics let him see. Just like a baker’s dough,our world is folded and extended.
Later in the twentieth century, the existence of chaos
has been confirmed independently and concurrently bymany authors in various fields: Lorenz (1963) in meteorology, May (1976) in ecology, Ueda et al. (1973) and Ueda
(1979) in electric circuits, Li and Yorke (1975) and Feigenbaum (1978) in one-dimensional maps, and Ruelle andTakens (1971) in turbulent flows. The so-called “fractal”
of Mandelbrot (Mandelbrot, 1977) turned out to underlie
these phenomena, and their relation led to the popularization of chaos. Later it turned out that Poincaré, Birkhoff,and Smale were observing the same phenomenon years
before.
The first astonishment among physicists concerned withchaos came from the observation that phenomena of ran-

Trends in 20th Century Celestial Mechanics
dom nature were visible in systems with low degrees of freedom. In fact, Landau and Lifshitz (1959) proposed expressing the motion of fluid by a quasi-periodic function of time(2). Landau and Lifshitz (1959) presumed that the number of frequency k would increase as the fluid approachedturbulence; the degree of freedom of a dynamical system
would increase as the flow became laminar to turbulent.
They expected that unpredictable motion would appear
only in a system with large degrees of freedom. Full turbulence should have an infinite degree of freedom. The reality is, however, even a system described by three ordinary differential equations exhibits chaos. Chaos appears
in a Hamiltonian system with only two degrees of freedom.
Interestingly, larger astonishment for celestial mechanists
arose from the fact that chaos was not restricted to Hamiltonian systems. Chaos is also seen in dissipative systems.
The essence of chaos is folding and extending. Its simplest representation was given by Smale (1967), a horseshoe: Extend a square, fold it in the middle, and put it on
the original square. This model provides us a glimpse intochaos. A chaotic system is exponentially unstable (i.e. coming from extension), but it comes back close to its original
location (i.e. folding). Repetition of folding and extending
forces the system to forget the memory of its initial conditions. Hence a chaotic system has a very high sensitivity
to initial conditions. Now, even the Laplace daemon getsseriously annoyed, faced with infinite repetitions. Is there
any point in predicting the future of a chaotic system? The
intricate and complicated structure of dynamical systems
that Poincaré did not want to describe troubles scientists
again. The complexity of chaos is ubiquitous in every scaleof our universe.
As the system becomes chaotic and loses regularity,
KAM tori eventually disintegrate. The final motion (i.e.
the states when time t → ∞) of a chaotic system depends
sensitively on its initial conditions. The initial value spacethat is mapped according to the different types of final motions shows us fractals. This fractal figure bifurcates andchanges its features when the external parameters change.
Moreover, this bifurcation itself takes place chaotically. The
remnant of the KAM curve is a Cantor set. It takes a longtime for a phase point to move through the holes of the Cantor set. In a Hamiltonian system with two degrees of freedom, when the last KAM curve disintegrates, a phase point
so far confined within a bounded region is now set free tomove toward infinity. However, according to numerical integrations, it is a slow process just after the disintegration.
One evidence of the slowness is the difficulty of the determination of the parameter value when the last KAM curve
disintegrates. See Fig. 5.3 of Chirikov (1979). In relation to
the disintegration of KAM curves, Hénon and Heiles (1964)
numerically discovered Hamiltonian chaos by observing thedisintegration of KAM curves. Hénon (1965) discussed thesame phenomenon using the restricted three-body problem
as an example. He devised a planar quadratic map which

61

exhibits a strange attractor (Hénon, 1976). This map is now
called the Hénon map.
Chaos is the ubiquitous and normal state of the universe. But this was not common sense until the final decadeof the twentieth century when the efforts of scientists began
to reveal chaotic phenomena in a variety of fields. Examples of chaos in our solar system were extensively dug up
and exposed to the light of science by Jack Wisdom fromthe viewpoint of the Lyapunov exponent (Wisdom, 1982,
1983, 1987a,b; Wisdom et al., 1984). Wisdom dissertatedabout chaos in the motion of small bodies in the solar system such as asteroids and a saturnian satellite, Hyperion.
He also found chaos in the orbital motion of major planets. Solar system chaos is extensively reviewed in Lissauer(1999).
Let us briefly introduce the definition of the Lyapunov
exponent. First, consider a one-dimensional discrete dynamical system denoted by f . Then, Lyapunov exponent
λ(x) of the orbit of x is defined by¯¯
1 ¯¯ df k (x) ¯¯λ(x) = lim ln ¯
.
(12)
k→∞ k
dx ¯
λ(x) represents the growth of the infinitesimal differenceof the initial conditions in phase space as a function of thenumber of iteration, k. If the growth is slower than an exponential one, the Lyapunov exponent is zero. Once a point
x is given, then the Lyapunov exponent is unique, i.e., any
point f j (x) on the orbit has the same Lyapunov exponent.
So the notion of the Lyapunov exponent corresponds to a
trajectory in phase space. For an n-dimensional discrete
dynamical system, there are n Lyapunov exponents. We
usually deal with the largest Lyapunov exponent.
The characteristic time of a chaotic system is often measured in the Lyapunov time (cf. Benettin et al., 1976). The
Lyapunov time is defined as the inverse of the maximum
Lyapunov exponent. By definition, the Lyapunov time is anaverage e-folding time of the separation of two different trajectories. The series of work by J. Wisdom, just mentioned,
extensively used Lyapunov exponents and Lyapunov time
to render visible the degree of chaos in solar system dynamics. Sussman and Wisdom (1988) found that the orbital motion of outer planets including Pluto was chaotic,
judging from the largest Lyapunov exponent in their longterm numerical integration. Sussman and Wisdom (1992)
also discovered that the Lyapunov time of the motions of
the four inner planets is only a few million years, indicatinga typical chaotic system. Their conclusions were later confirmed by an accurate semi-analytical perturbation theoryof J. Laskar (Laskar, 1989, 1990, 1994, 1996; Laskar et al.,1992).
However, a short Lyapunov time does not mean thata chaotic system becomes unstable and disintegrates very
quickly; the Lyapunov time is different from the global instability time of a dynamical system. One of the possible
explanations for this fact is that the Lyapunov exponent is62Takashi Ito and Kiyotaka Tanikawa

a local notion. Exponential separation, when it ever happens, takes place in a linear regime: two trajectories inphase space separate from each other exponentially in time,
starting from their infinitesimally small distance. Strictly
speaking, we do not know what happens when the separation becomes macroscopically large. As we mentioned, theessence of chaos is folding and extending. The positive Lyapunov exponent implies extending. But to which extent?
The folding process may bring two trajectories back to theinitial small separation.
A lot of evidence indicates that the planets in our solar system have survived nearly 4.6 billion years, far longer
than their Lyapunov time, with a global stability (e.g. Ito
and Tanikawa, 2002). This typically shows that Lyapunovtime TL does not help us very much in predicting the real
instability time TI of the system. There have been, however, some efforts to understand the relationship between
TL and TI , utilizing the dynamical behavior of solar systemminor bodies. For example, some asteroids go into unstableorbits over only a few hundred to a few thousand years, sothese objects might be used for determining the relationship between TL and TI (Lecar and Franklin, 1992; Lecaret al., 1992). Lecar et al. (1992) proposed a relation
TI ∝ TLγ ,(13)through a variety of numerical experiments for asteroids
with unstable orbits. The exponent γ looked close to 1.8
in their experiments. But the dynamical meaning of either
relationship (13) or the value of γ ∼ 1.8 is not clear yet,though some qualitative explanations exist (e.g. Morbidelli
and Froeschlé, 1996). If relationship (13) holds true in manycases, a numerical integration of the length of Lyapunovtime TL of the system would suffice to tell us the real instability time TI , which could save us computer time.
A chaos that does not cause any real or global instability over a much longer time than the Lyapunov time is often
dubbed a “stable chaos” or “weak chaos” (e.g. Froeschlé
et al., 1997; S̆idlichovský, 1999; Tsiganis et al., 2002). For
example, many asteroids are dynamically stable for much
longer than their Lyapunov timescales. This kind of stablechaos in asteroid dynamics is relevant to high-order mean
motion resonances with Jupiter in combination with secular perturbations on the perihelia of the asteroids. These
perturbations move the asteroid orbit from one high-orderresonance to another in some irregular ways (e.g. Milani
and Nobili, 1997). Chaos in asteroidal orbital motion has
also been extensively studied by N. Murray and M. Holman (Murray and Holman, 1997, 1999; Murray et al., 1998).
Murray and Holman (1997), in particular, presented ananalytic theory of asteroid motion near resonances in theplanar elliptic restricted three-body problem. Their theory predicts the location and extent in the semimajor axis
and eccentricity space of the chaotic motion, the Lyapunov
time, and the time for objects on chaotic orbits to be removed from the system. Their theory predicts that aster-

oids in a number of high-order mean motion resonances
possess very short Lyapunov times, such as 105 years, butwith instability times as long as the lifetime of our solar
system. Incidentally, we should notice that apparent chaos
in solar system dynamics sometimes comes out as a resultof numerical error or instability (cf. Rauch and Holman,
1999; Ito and Kojima, 2005).
Before closing this subsection, we would like to mention
a bit about classification of the degrees of chaos in Hamiltonian systems. Hamiltonian systems have neither attractors(i.e. sinks) nor repellers (i.e. sources). Phase points straggle
out almost without purpose in the phase space. Everywhereis filled with complicated structure. When points approach
the local KAM tori, they are trapped for a long time ina stagnant region. It seems difficult to obtain an averagetime for a point to stay at a particular place (or region) inphase space. This structure becomes more and more entangled as the system goes further from the integrable state.
Strange as it might sound, there appears a random system
as a limit of this complexity, and the average staying timeof a point can be calculated.
The complexity is measured by topological entropy
and/or measure entropy. The topological entropy is a measure of maximum complexity of a system, whereas the measure entropy is the average complexity over phase space.
The definition of topological entropy is not so straightforward. For simplicity, consider a discrete dynamical system
f of a compact phase space X. In general, two orbits ina dynamical system separate with time. Let B(x, ε, n) bea set of points y such that the largest of d(f i (x), f i (y))(i = 0, . . . , n − 1) is less than ε. Here f i (·) is the i-th iterateof the map, and d(·, ·) is a distance function. Let N (ε, n)be the minimum number of B(x, ε, n) to cover X. Then,
topological entropy htop (f ) is defined as·¸
1
htop (f ) = lim lim sup ln N (ε, n) .
(14)
ε→0
n
n
Suppose that we have a finite resolution ε. We can notrecognize two different states as “different” when their separation is smaller than ε. We can roughly say that N (ε, n)is the number of states we can recognize after the n-th iteration. Then n1 ln N (ε, n) is the exponential increase (periteration) of the number of distinguishable states. Finally,
htop (f ) is the asymptotic value of information production
rate per iteration for infinitesimal resolution.
Both of the entropies, particularly measure entropy, arevery difficult to calculate. Some methods have been developed for estimating the lower bound of the topologicalentropy for a given system. The simplest way is to find ahorseshoe. For example, the lower bound is ln 2 if the Smale
horseshoe exists (cf. Section 15.2 of Katok and Hasselblatt(1995)). Another method is to find a periodic orbit fromwhich a non-trivial braid can be constructed. In the caseof the Hénon map, the lower bound of the topological entropy is measured by the increasing rate of the number ofTrends in 20th Century Celestial Mechanics

63

periodic orbits per period (cf. Sterling et al., 1999).
In general, the larger the structure, the more responsible shorter periodic solution is. Here, by the structure wemean the distribution of the so-called islands and sea in thephase space. Some of the bifurcated solutions are unstablefrom the start, and others are initially stable and become
unstable later. Stable and unstable manifolds emanatingfrom the unstable solutions intersect each other, as well asother unstable and stable manifolds from other unstable
solutions. They together make a complicated network of invariant manifolds. Phase points move along this network.
The qualitative nature of the network differs from system
to system. The network of the phase space of the solar system and that of a simple three-body problem are its examples. Our principal goal may be to clarify the structure of
the network of individual systems. A global network connected to infinity leads to the instability of the system. In
systems with degrees of freedom larger than two, any phase
point seems to be close to this global network. This is theorigin of the Arnold diffusion that leads to instability.

also express the dynamics of Hamiltonian systems with twodegrees of freedom.
Birkhoff later made some remarkable contributions tothe study of twist maps (Birkhoff, 1920). One of them wason the functional form of the invariant curves, now called
KAM curves. Another was on the property of points in
the instability zone sandwiched by two invariant curves.
The following theorem represents an important part of
Birkhoff’s contribution:

2.3. Twist maps

Theorem 2.3 (Birkhoff ) Let C 0 , C 00 be entirely distinct
invariant curves forming the boundary curves of a ring of
instability. Then, for any ε > 0 an integer N can be assigned
such that a point P 0 exists within a distance ε of any point
P of C 0 (or C 00 ) which goes into a point Q0 within ε of any
point Q of C 00 (or C 0 ) in n < N iteration of T (or T 0 ).

There is a direction of studies in which some twodimensional maps of a Hamiltonian system with two degrees of freedom, in particular twist maps, are used toinvestigate the start, development, and degree of chaos.
The study of twist maps again goes back to Poincaré.
He considered a theory of consequents: a theory of the
Poincaré surface and Poincaré map in the present terminology. The example object of the study was the restricted
three-body problem. Observing periodic orbits of the problem, Poincaré took the Poincaré map of its neighbor. Ifneighboring points move on a closed curve around the fixed
point, then the fixed point is stable. Behavior of points near
the fixed point can be viewed like that in a punctured disk,
i.e. a ring-like structure when the fixed point is removed
from its center and when the center is blown up. Poincaréarrived at his “last geometric theorem” (Poincaré, 1912)
starting from the surface map of the periodic orbits of therestricted three-body problem. Let us describe a famous
theorem known as Poincaré–Birkhoff:
Theorem 2.1 (Poincaré–Birkhoff ) Let us suppose thata continuous one-to-one transformation T takes a ring R,
formed by concentric circles Ca and Cb of radii a and b
respectively (a > b > 0), into itself in such a way as toadvance the points of Ca in a positive sense, and the pointsof Cb in the negative sense, and at the same time to preserveareas. Then there are at least two invariant points.
This kind of map is now called a twist map because
it twists a ring. Birkhoff (1913, 1925) started his substantial scientific career proving this theorem. Twist maps havebeen used in various contexts. One of the main reasons forthis is that we can treat the maps very easily. The maps

Theorem 2.2 (Birkhoff ) In a small neighborhood of anelliptic point with irrational rotation number where the twist
is monotone, any invariant curve enclosing the invariant
point meets every radius vector through the invariant pointin only one point. If the barred angle in the plane be drawn
at the corresponding point, the curve lies entirely within it
on either side in the vicinity of the point.
In other words, the invariant curve encircling the fixed
point is a graph of a Lipschitz function of polar angle. The
following theorem describes the transverse motion of points
in a zone of instability.

The next major step marked in the context of twist
maps is the KAM theorem (or the KAM theory). In Section
2.1 we explained Arnold’s proof (Arnold, 1963). In the context of twist maps, Moser (1962) discussed the persistenceof invariant curves of monotone twist maps of the annulus
θ1 := θ + α(r) + F (r, θ),(15)

r1 := r + G(r, θ),(16)which are derived by adding small perturbations to an integrable case
θ1 := θ + α(r),(17)

r1 := r,

(18)

where
dθ1
> 0,
dr

(19)

and F, G are small perturbations periodic in θ. Literally, given an integrable monotone twist map and smoothenough small perturbations, then invariant curves of rotation number ω satisfying Diophantine conditions
ε|2πp − qω| ≥ 3 ,(20)
q2
survive the perturbation.
It is to be noted here that recently monotone twist maps
are most frequently studied. So we usually abbreviate these64

Takashi Ito and Kiyotaka Tanikawa

maps simply as twist maps. Original twist maps considered by Poincaré and Birkhoff do not necessarily satisfy
the monotone condition (19). Hereafter in this manuscript,“twist maps” basically denotes “monotone twist maps” unless otherwise noted.
Chirikov (1979) discussed the Arnold diffusion mainly
using so-called standard maps defined on the cylinder as
θ0 := θ + y + a sin θ,(21)

y 0 := y + a sin θ,(22)

where a is a parameter. The Arnold diffusion cannot be
observed in two-dimensional maps. However, Chirikov proposed a concept called “resonance overlap” which is a lower
dimensional analogue of the heteroclinic tangle of stableand unstable manifolds. Resonance overlap is the main
pathway of the Arnold diffusion in higher dimensional systems (see also Section 4.4.1).
At nearly the same time as Chirikov, Greene (1979)
paid attention to the disintegration of the KAM curves. According to the KAM theory, invariant curves survive longer
when they are with the rotation numbers that are difficult
to be approximated by rationals. Greene tried to determine the critical parameter value of KAM curves with a
special type of irrational rotation number, hoping that the
KAM theory would be still valid for this parameter range.
Greene’s tool was the residue criterion which has been waiting for a rigorous proof until now. Here the residue R of aperiodic point isR=

2 − λ − λ−1
,
4

(23)

where λ is the eigenvalue of the periodic point. A regularsaddle has R < 0, an elliptic point has 0 < R < 1, and ainversion saddle has R > 1. His idea is that if the KAM
curve is disintegrated, the nearby periodic points will be
unstable. Let pn /qn be the convergent rationals, like a Fibonacci sequence, to irrational number ω. Then, according
to Mackay (1992), the residue criterion is described as follows:
Conjecture 2.1 Let Rn be the residue of a Birkhoff periodic point with rotation number pn /qn . Thenν(ω) = lim qn−1 ln |Rn |

(24)

n→∞

exists, and ν(ω) ≤ 0 implies there is a KAM curve withrotation number ω, whereas ν(ω) > 0 implies there is none.
√

The golden mean rotation number 1+2 5 has the continued fraction expression [1, 1, 1, . . .], and is the most distant number from rationals. The last KAM curve is the onethat has the golden mean rotation number. Greene (1979)
estimated the critical parameter value for the destructionof the last KAM curve to be 0.971635 · · ·, setting R = 0
at the critical situation. Even now, this value is not fully
explained theoretically.

Later, Mather (1982) as well as Aubry and LeDaeron(1983) succeeded in showing what happens if KAM curves
are destroyed. Katok (1982) gave elementary proof of this.
Motion on each KAM curve, if destroyed, remains as quasiperiodic, but the curve itself becomes a Cantor set. Thiskind of motion has already been investigated in a circle map
by Denjoy (1932). A so-called “Denjoy counter-example”
exists when the smoothness of the map is less. In this case,the motion with irrational rotation number on the circle
does not fill the circle itself, but fills only a Cantor set. In
twist maps, smoothness conditions of the maps are different from those of one-dimensional maps. KAM curves become Cantor sets even when the map is analytic. In twist
maps, after the destruction of KAM curves, transverse motion across the former KAM curves is possible. The hole of
the Aubry–Mather set (remnants of the KAM curve) becomes larger as the perturbation parameter grows. In deriving the Aubry–Mather theory, the classical variational
principle is used. Monotone twist maps are extensively investigated in the 1980s and the early 1990s (e.g. Mather,
1984, 1985, 1991; Aubry and Abramovici, 1990; Herman,
1983, 1986; Le Calvez, 1987, 2000; Mackay and Percival,
1985; Mackay, 1992, 1993). Bangert (1988) gave an excellent review. Mather (1991) demonstrated the existence oforbits with complicated behavior in a zone of instability.
In the meantime, efforts to estimate the topological entropy of low-dimensional dynamical systems by constructing braids has begun (e.g. Matsuoka, 1983). This method
has been applied to twist maps (e.g. Hall, 1984). A periodic orbit of a surface map can be represented as a braid.
In fact, suppose a discrete dynamical system f is given on
a disk, and f has a period-n orbit. We put the second disk
below the first one separated by a distance. Let pi and p0i(i = 1, 2, . . . , n) be orbital points of the same periodic orbit on the disks above and below. We connect these points
by strings so that p1 is connected to p02 , p2 to p03 , . . ., andfinally pn to p01 . This is exactly the braid with n strings(“n-braid”, in short). We call pi and p0i (i = 1, . . . , n) thestarting and end points. We regard any two braids equivalent if they transform each other under continuous movements of starting and end points together. An equivalentclass forms a braidtype. The set of braidtypes is the set of
inequivalent braids. The braid is an invariant characterizing topological behavior of orbits for 2-dimensional maps.
If the system becomes more chaotic, braids that are formed
with newly born periodic orbits may become complicated.
We have a convenient theory to estimate topological entropy of braids (Matsuoka, 1993). Thus, we can estimate
topological entropy of dynamical systems by constructing
braids from periodic orbits.
At the turn of the century, the study of twist maps seems
to be in an interlude: Researchers seem to have lost theirway. In which direction are we to go? In this opaque age,some people raise their hands and point in the direction of
non-Birkhoff periodic orbits, i.e., non-monotone periodic

Trends in 20th Century Celestial Mechanics
orbits. We hope that this direction yields a promised land.
Non-monotone periodic orbits or non-Birkhoff periodic orbits are the orbits that do not preserve the orbital orders.
Let us briefly explain these orbits. Consider a twist
map f defined in the infinite cylinder S 1 × R. We lift themap to the universal cover R × R and denote it by fb. Let{fbk (bz )}k∈Z be the lift of orbit {f k (z)}k∈Z where zb is a liftof point z. We define the extended orbit of zb byeo(bz ) = {fbk (bz ) + (l, 0)}k,l∈Z .

(25)Then, the orbit of z is monotone if the following condition
π1 (ξ) < π1 (ζ) =⇒ π1 (f (ξ)) < π1 (f (ζ))(26)

is satisfied for any ξ, ζ ∈ eo(bz ) where π1 is the projection to
the x-coordinate. If a periodic orbit is monotone, it is sometimes called “Birkhoff”. If a periodic orbit is non-monotone,
it is sometimes called “non-Birkhoff”. In a one-parameterfamily of twist maps containing the integrable map, nonBirkhoff periodic points appear (or bifurcate) when themap becomes non-integrable. There are many types of nonBirkhoff orbits that represent the chaosity of the system.
Reversibility is the property of a dynamical system
originated from the Hamiltonian (Zare and Tanikawa,
2002). In reversible area-preserving twist maps, we can findperiodic orbits without appealing to variational principles(e.g. Tanikawa and Yamaguchi, 1987, 1989). Dynamicalorders of periodic orbits in these maps are derived andlower bounds of topological entropy of the system are obtained by constructing braids from periodic orbits (e.g.
Yamaguchi and Tanikawa, 2000, 2001a,c, 2002a,b,c, 2003,
2005a,b,c; Tanikawa and Yamaguchi, 2001, 2002a,b, 2005).
Non-Birkhoff periodic orbits promise us a rich structure oftwist maps, such as the existence of oscillatory orbits inthe standard map (e.g. Yamaguchi and Tanikawa, 2004a).
Exponential splitting of stable and unstable manifolds isnumerically inferred from a sequence of homoclinic points
on the y axis accumulating at the saddle fixed point in
the standard map (e.g. Yamaguchi and Tanikawa, 2001b).
Quadruply reversible non-twist maps have also been discussed (e.g. Zare and Tanikawa, 2005; Yamaguchi andTanikawa, 2004b).
2.4. Integrability
The study of the integrability or non-integrability ofHamiltonian systems, as far as we understand, seems related to the formal side of celestial mechanics. As is wellknown, celestial mechanics started in the seventeenth century in Europe with the aim of mathematically explaining Kepler’s three laws of planetary motion that were discovered at the beginning of that century. The explanation
was first accomplished by I. Newton in his Principia, finding a force law that satisfied Kepler’s laws. For the one
hundred years after Principia, a lot of ingenious celestial

65

mechanists and mathematicians developed differential calculus, trying to establish the laws of force on metaphysicalprinciples: the principle of virtual velocity, the principle ofd’Alembert, the principle of Maupertuis, and so on. Finallynear the end of the eighteenth century, J. L. Lagrange formulated the equations named after him, which are equivalent to Newton’s equations of motion, under the principle of least action (Lagrange, 1788). Lagrange gave a unified interpretation of various principles. In the nineteenth
century, W. R. Hamilton, starting with the formulation of
the path of light rays in geometrical optics through variational calculus, gave his Hamiltonian formalism of equations of motion, and discussed the transformation theory ofdynamical systems. The existence of ignorable or cyclic coordinates motivates transformations from the original system to a system with fewer degrees of freedom. Finally, it
has been recognized that the ultimate goal of the study ofdynamical systems is to solve the Hamilton–Jacobi equation. This partial differential equation is for a generating
function which gives the transformation from the original
dynamical system to a system in which all conjugate variables are constant through motion. Thus, the formal sideof mechanics has been completed in the sense that the goalof the study was revealed.
How to attain the goal, then? We do not have a systematic or automatic procedure for solving the Hamilton–
Jacobi equation. We know only a few systems for which the
Hamilton–Jacobi equation can be solved. One of the possible procedures for attaining the goal is the perturbative
approach when the system is known to be a perturbation of
an integral system. We repeat canonical transformations,and move to systems which have terms of non-integrable
nature in higher and higher orders. We hope these terms tobe driven away to an infinitely higher order, eventually obtaining an integrable system. This approach is not applicable to a system that is not close to any of known integrable
systems. Another approach is to find additional integrals.
As is easily shown, once these integrals are incorporated
into the equations of motion, the degrees of freedom of thesystem can be reduced. At the end of the nineteenth century, Poincaré discussed various kinds of integral invariantsof Hamiltonian systems. He related them to the existenceof additional integrals of the system.
In the first part of Section 2 (p. 56), we already mentioned that Poincaré had succeeded in proving the analytic non-integrability of the restricted three-body problem in the late nineteenth century. More precisely, Poincaré
proved the non-existence of an integral that is independentof the Hamiltonian and that can be developed in a power series of a small parameter 2 . Shortly before Poincaré, Bruns
2This latter condition was too restrictive. So researchers later relaxed the condition, and discovered that there can be some integrable
systems for fixed perturbation parameters. For example, three-bodysystems are known to be integrable when they have particular mass
combinations. This fact was proved through the discussions on col-

66

Takashi Ito and Kiyotaka Tanikawa

proved the algebraic non-integrability of the general threebody problem. In what follows, let us briefly review the
proofs by Bruns (1887) and by Poincaré (1890), both ofwhich are quite suggestive to us even now. The ideas underlying the two proofs differ very much.
It is not quite easy to summarize Bruns’ proof of algebraic non-integrability in a few paragraphs, because it isvery deep and contains a lot of steps (cf. Whittaker (1904),
Chapter XIV, §164). Understanding process of this kindof proof is to dive into deep water and stay there, looking
for precious stones in the depths. Here we only try to describe a rough story. Bruns first expresses formally a supposed integral as a function of coordinates and momenta,
and proves that the integral must involve some of the momenta. Mutual distances are irrational functions of coordinates. Bruns introduces the sum s of mutual distancesri , adds s to dependent variables, and tries to express integrals as an algebraic function of coordinates, momenta,
and s. From this point begins the main stage of the proof.
Bruns reduces the problem of finding algebraic integrals
into that of finding an integral as a quotient of real polynomials of (qi , pi , s). Then, after a rather long argument,
he shows that he can make an integral from the numerator
or denominator. In fact, he expresses the numerator or denominator as φ = φ0 + φ1 + φ2 + . . ., where φi are homogeneous functions of momenta and φ0 is the highest term in
momenta. Bruns shows that if φ0 does not contain s, thenmultiplying a rational function to φ, he can obtain an integral. Then he shows that φ0 actually does not contain s.
Here, s = ±r1 ± r2 ± r3 has generally eight different values.
However, it may have the same value for special values ofri . Bruns uses this fact to derive a contradiction. Finally,
Bruns shows that φ0 is a function of the momenta and theintegrals of angular momentum, and that φ is expressible
as a function of the classical integrals. It seems that theparticular properties of the N -body problem used by Bruns
are the weighted homogeneousness of the equations of motion, which is later used by Yoshida (1983a,b).
On the other hand, Poincaré proved generally the nonexistence of an integral of the restricted three-body problem
depending analytically on coordinates and momenta andalso on the mass parameter. His proof proceeds as follows.
Poincaré expands the supposed integral Φ as a power series
in the mass parameter µlisional manifold using the McGehee variables. Readers can consult
McGehee (1974) for the case of one-dimensional three-body problem, as well as Devaney (1980) for the case of planar isosceles threebody problem. Another, simpler example is a system described by a
HamiltonianH=

p21 + p222

q3
+ q12 q2 + ε 2
3(27)

with a fixed perturbation parameter ε. A dynamical system dominated by H is algebraically integrable when ε = 0, 1, 6 or 16 (Yoshida,
1983a,b).

Φ = Φ0 + µΦ1 + µ2 Φ2 + . . . ,(28)

where Φi is analytic in (qi , pi ), and periodic in q1 and q2 . It
is to be noted that Φ itself is analytic both in (qi , pi ) and µ.
He expands the necessary and sufficient condition for the
integrability, i.e., the vanishing of the Poisson bracket{H, Φ} = 0,(29)

and gets two conditions{H0 , Φ0 } = 0,

{H1 , Φ0 } + {H0 , Φ1 } = 0,(30)

from the lowest two orders of the expansion.
Poincaré first proves that Φ0 is not a function of H0by a smart, iterative argument. Then he proves that Φ0
includes neither q1 nor q2 by a non-degeneracy conditionof the unperturbed Hamiltonian. This non-degeneracy waslater used as a basic hypothesis in the proof of the KAM theorem. With H0 and Φ0 not involving coordinates, Poincaré
proceeds to show that the existence of the uniform integral
is not compatible with the property that Φ0 is not a function of H0 . This problem is reduced to a commensurability
problem: the commensurability condition
m1

∂H0
∂H0
+ m2
=0
∂p1
∂p2

(31)

is satisfied in every small region of (p1 , p2 ). The vanishing of
the Poisson bracket brings the functional dependency of H0and Φ0 at infinite number of commensurable frequencies.
This leads to the analytical functional dependency of Φ0and H0 , which contradicts the independency of Φ0 and H0 .
Ever since the proofs of algebraic and analytic nonintegrability of the three-body problem by Bruns andPoincaré that we have summarized, there might have beenan acceptance that the study of integrable systems would
not contribute to the understanding of general systems.
This is because integrability requires a stringent condition, and also because integrable systems are special andvery rare. General systems must be treated perturbatively
or numerically. Of course, a deeper understanding of thestructure of Hamiltonian or Lagrangian equations of motion or mechanical systems is required as well. Looking at
the swarm of researchers attacking Hamiltonian systems
from various aspects, it seems to us that the study of integrability is one of the most orthodox paths leading to
the next reformulation or reinterpretation of Hamiltonianmechanics in general.
In the latter half of the twentieth century, the notionof integrability itself evolved. In the 1960s, geometrical interpretation was given to the concept of complete integrability. Naive integrability may just mean an integration ofequations of motion by quadrature. Hamiltonian systems
with n degrees of freedom are said to be integrable in thesense of Liouville, or simply Liouville integrable when there
are n independent first integrals in involution. Here those
first integrals are said to be in involution if Poisson brackets of each pair of them vanish. The first integral I(x) is aTrends in 20th Century Celestial Mechanics
function whose time derivative in the direction of the vector field vanishes. Arnold added geometrical interpretation
to this notion (Arnold, 1978). The Arnold–Liouville integrability (in the simplest setting) in addition to the statement of the Liouville integrability says that phase spaceis divided into n-dimensional invariant tori, and that the
motion takes place in each of these tori.
According to Goriely (2000), there are three approaches
to the study of integrability or non-integrability of dynamical systems: a dynamical system approach, an algebraic approach, and an analytic approach. A dynamicalsystem approach starts with a search for fixed points or periodic points. This approach proceeds to an analysis of theirstability, followed by various analyses such as topologicalconjugacy, bifurcation analysis, and so on. This approachseems to fit the study of non-integrable systems.
The algebraic approach has its origin in Bruns. One
looks for the first integrals of algebraic functions, in particular, polynomial functions. If one substitutes a polynomial of degree k into the equation satisfied by the supposed
first integral, then the problem of finding the first integral
is reduced to an algebraic problem: the problem of solving
a system of linear equations for its parameters (or coefficients). There is a notion of “second integral”. Second integral is an invariant relation for a subset of phase space,and it serves as the building blocks for the first integral.
In the analytic approach, one looks at the local behavior of the solutions of equations of motion around their singularities in complex time. The singularity may be a pole,a branch point, or an essential singularity. The coefficientmatrix of the variational equations around some particularsolution is called Kovalevskaya matrix, and its eigenvalues
are called Kovalevskaya exponents (Yoshida, 1983a). It has
been shown that the behavior of these exponents is strongly
related to the integrability of the system (Yoshida, 1983b).
In this approach, one looks for meromorphic integrals, i.e.
holomorphic (or analytic) integrals that do not have singularities other than poles.
The study of integrability was partly revived when solitons were found in partial differential equations. The Hirota method (Hirota and Satsuma, 1976) and the famous
GGKM paper (Gardner et al., 1967) should be consulted
here. The study of integrability before 1983 was extensively
summarized by Kozlov (1983).
In the mid 1970s, the integrability of the finite Toda
lattice was proved (Hénon, 1974; Flaschka, 1974). Ziglin(1983) gave a necessary condition for the integrability ofHamiltonian systems using variational equations aroundparticular solutions. This condition can be used as sufficientfor the nonintegrability of the system. Yoshida (1983a,b)
has made a pioneering contribution to the integrability
and nonintegrability of Hamiltonian systems. Yoshida proposed a procedure to obtain necessary conditions for theweighted-homogeneous systems of differential equations to
be integrable. Algebraic integrals do not exist if none of the67

Kovalevskaya exponents is rational, where Kovalevskaya
exponents characterize the stability of certain singular solutions. This property of Kovalevskaya exponents is related
to the single-valuedness of the solution of differential equations.
Yoshida considers a similarity invariant system of differential equations. Here, similarity invariance is the invariance under transformation
−1
 t → α t,





x1 → αg1 x1 ,



(32)
x2 → αg2 x2 ,



.

..






 x → αgn x ,nnwith rational numbers g1 , . . . , gn and a constant α. Thissystem of equations permits a particular solution originating in the similarity itself and corresponding to a collision(collapse) orbit in n-body systems. Solutions of the variational equations around this particular solution can be obtained from the eigenvalues and eigenvectors of an n × ncomplex constant matrix K = (Kij ) with
Kij =

∂Fi
(c1 , . . . , cn ) + δij gi
∂xj

(33)

formed from the coefficients of the variational equations.
This matrix is given the name “Kovalevskaya” after the celebrated contribution of S. V. Kovalevskaya on the motionof a rigid body around a fixed point in the late nineteenth
century. The characteristic equations for the Kovalevskayamatrix is called the “Kovalevskaya determinant”, and theroots of this equation are called the “Kovalevskaya exponents”. Later, these notions play important roles in the research in this area.
Now Yoshida introduces a weighted homogeneous polynomial φ of weighted degree M asφ(α−1 t, αg1 x1 , . . . , αgn xn ) = αM φ(t, x1 , . . . , xn ).

(34)

Repeating Bruns’ discussion on the non-existence of algebraic integrals, Yoshida reduces the existence of integrals
into the existence of the rational first integrals of weighted
homogeneous functions. Then he arrives at the main results. If the similarity invariant system of differential equations has a weighted homogeneous first integral of degree
M , then M is a Kovalevskaya exponent. If the system is
described by a Hamiltonian, and if the Hamiltonian has
a weighted degree h and admits a weighted homogeneous
first integral of weighted degree M , then M and h − 1 − M
are Kovalevskaya exponents.
Yoshida (1983b) shows that irrational or imaginary Kovalevskaya exponents are inconsistent with the existence ofa sufficient number of algebraic first integrals. Before showing this, he introduces the algebraic integrability of general

68

Takashi Ito and Kiyotaka Tanikawa

differential equations. Then he tries to extend his result tonon-similarity invariant systems.
The integrability of Hamiltonian systems and thesingle-valuedness of solutions near singularities have long
been suspected to have intimate connections. A system ofordinary differential equations is said to have the Painlevé
property when its general solution does not have movablecritical singularities. Here, a movable singularity meansthat the position of the singularity depends upon theinitial conditions. A “critical singularity” means that asingle-valued function, after analytic continuation, has different values at this singularity. The procedure for checking
whether or not a given system has the Painlevé property
is called the Painlevé test. Though it is widely believedthat the Painlevé property is incompatible with chaotic
motions, there is no rigorous proof of this simple statement
yet to date (Goriely, 2000). For reference, Morales-Ruizand Ramis (2001a,b) obtained a strong necessary condition
for the integrability of Hamiltonian systems, extending
the work of Ziglin (1983). The condition says that if theoriginal Hamiltonian system is Liouville-integrable, thenthe variational equation around a particular solution issolvable by a combination of quadratures, exponential ofquadratures, and algebraic functions.
Applying a theorem of Morales-Ruiz and Ramis(2001a,b) to the two-degrees-of-freedom Hamiltonian,
Yoshida (1999) solved the weak Painlevé conjecture positively. The conjecture in Yoshida (1999) says that when
the Hamiltonian H = 21 (p21 + . . . + p2n ) + V (x1 , . . . , xn ) has
a weak Painlevé property, then it is Liouville integrable. As
to the definition of the weak Painlevé property, see Goriely(2000) for example. Nakagawa and Yoshida (2001a,b) derived a necessary condition for the integrability of homogeneous Hamiltonian systems with two degrees of freedom.
They made a list of all integrable two-dimensional homogeneous polynomial potentials with a polynomial integralof the order of at most four in momenta.
It has been shown that the transversal intersections ofstable and unstable manifolds, i.e. the existence of a horseshoe, implies the non-integrability of dynamical systems(Moser, 1974). Conversely, the non-integrability implies theexistence of transversal intersections of stable and unstablemanifolds in some dynamical systems. One famous example is the standard map (Lazutkin et al., 1989). It seems
strange to the authors that the latter property is not yet
proved rigorously. Numerically, the latter property seems
obvious. In fact, “chaos comes from infinity,” i.e. the number of iterates necessary for completing the horseshoe becomes infinite as the system approaches integrable.
Suppose we take a look at an arbitrary small neighborhood of a point in phase space. How are the set of solutions
passing near the point geometrically disposed? The disposition is closely related to the integrability of the system.
The range of potential dispositions increases with the degrees of freedom. At one extreme, there is a possibility thatphase space is strictly stratified by hypersurface. At theother extreme, any solution curve wanders through all dimensions in the neighborhood of the point. This difference
should be reflected in the functional form of solutions.
3. The three-body problem
The three-body problem has a long history, and there
has been a great deal of previous research. This problem
has attracted so many astronomers, physicists, and mathematicians, in particular in the twentieth century. This isbecause this problem contains many aspects of general dynamical systems. Let us point out some of them: The threebody problem contains (collisional) singularities which giverise to rich dynamical phenomena. There are various kindsof periodic orbits and related homoclinic and heteroclinic
phenomena in the three-body problem. This problem canbe used as a testbed of newly devised numerical methods such as symmetric schemes and symplectic integrators.
New concept of integrability, such as integrability in thesense of C n , is introduced as a part of the three-body problem research. There are various versions of the three-bodyproblem: restricted, collinear, isosceles, planar, and three
dimensional, each of which asserts its own raison d’etre.
For a particular choice of mass, this problem is reduced toperturbed two-body problem that is close to be integrable.
Thus, the three-body problem covers from near-integrable
realm to chaotic realm of dynamical systems. The degreeof freedom is from two to four. The three-body problem
has many applications in astronomy and in physics. Tripleinteractions of particles are regarded as a basic process inthe dynamics of stars and galaxies. Solutions like oscillatory and non-collision singularity, or “figure-eight”, may
be particular to this problem. Finally but not the last, thethree-body problem is beautiful as it is.
In this section we categorize the study of the three-bodyproblem into a few parts: Analytical and qualitative studies, numerical studies, and hierarchical triple systems. We
also mention some current frontiers of three-body problem
research including some recent applications to actual stellar and planetary systems. Due to our inability, many important topics related to the three-body problem are notvery well described in this section, such as the restricted
three-body problem as an independent field of study.
3.1. Analytical and qualitative studies
After the rigorous proofs of algebraic and analytical
non-integrability by Bruns and Poincaré at the end of
the nineteenth century (see Section 2.4), studies of thethree-body problem moved toward qualitative directions(cf. Whittaker, 1904). Two famous example problems were
proposed. One is about the existence of non-collision singularity in N -body systems proposed by Painlevé (1897). The
other is about the existence of oscillatory solutions in theTrends in 20th Century Celestial Mechanics
three-body system (Chazy, 1922). The former was solved
by Xia (1992). The latter was first solved by Sitnikov (1960)
for the restricted three-body problem, then later solved
through a perturbative approach in the three-dimensional
isosceles three-body problem by Alekseev (1968a,b, 1969),
and finally, in the planar isosceles problem by Xia (1994).
C. Siegel questioned why the three-body problem was not
integrable. Siegel (1941) showed that triple collision is an
essential singularity, and that triple collision is the originof bad behaviors in the three-body problem.
Chazy (1922) described all possible types of final motion, i.e., behavior as t → ±∞, of the three-body problem.
After Chazy, from the 1940s to the 1960s, qualitative studies of the three-body problem flourished in the former Soviet Union (USSR). Final motions are classified into H, P,
E, B, and OS, corresponding to hyperbolic, parabolic, elliptic, bounded, and oscillatory motion. The initial motion
is suffixed by − and final motions are suffixed by +. For
example, HEj− → HEi+ means that as t → −∞, one particle escapes (H) and the remaining two particles form abinary (Ej ) where the escaping particle is the j–th particle.
As t → ∞, one particle escapes (H) and the remaining two
particles form a binary (Ei ) where the escaping particle isthe i–th particle. In this particular example, the motion iscalled “exchange” when i 6= j, and the measure of this kindof orbits in phase space has been proved to be positive.
Classifications of the final and initial motions of thethree-body problem were intensively carried out (e.g. Merman, 1958; Khilmi, 1961). This was a modest approachcompared with the conventional pursuit of individual orbits. In the new approach, we are content with understanding the behavior of a bundle of orbits. The possibility ofcapture, escape, and exchange was discussed, and some important quantitative results have been obtained. See thereviews by Hagihara (1971) and Alekseev (1981) for moredetail.
Saari (1971) studied the relation between collision
singularity and the N -body problem. McGehee (1974) introduced a remarkable set of coordinates in the three-bodyproblem. The coordinates are remarkable in the sense that
triple collision is blown up to a manifold and orbits approaching triple collision move along the fictitious orbits
on this collision manifold. Owing to these coordinates,
theoretical research have been revitalized. People such as
R. Moeckel, C. Simo, and R. L. Devaney have deployed the
McGehee variables and analyzed the behavior of solutionsof the three-body problem near triple collision. Recently,
computer-aided research has tried to reveal the structureof general phase space in various settings (e.g. Tanikawa
et al., 1995; Umehara, 1997; Zare and Chesley, 1998;
Umehara and Tanikawa, 2000; Tanikawa and Mikkola,
2000a,b; Tanikawa, 2000). There is a special field of study
in which escape criteria are the target (e.g. Marshal, 1990);
Yoshida’s criteria (Yoshida, 1972, 1974) are considered as
the best at present.

69

3.2. Numerical studies
Numerical studies of the three-body problem have a relatively long history, dating back to an era before the arrivalof digital computers. Numerical integrations of the periodic
solutions of the restricted three-body problem were carriedout by two groups in the first half of the twentieth century,
either by hand or using the Tiger accumulators. One ofthese groups was called the Copenhagen school, headed by
E. Strömgren (Strömgren, 1935). The other was a Japanesegroup headed by T. Matukuma (Matukuma, 1933). Afterthe introduction of electronic computers, Hénon (1965) systematically followed the scenario of the Copenhagen school,and started the quest for periodic orbits. Throughout thelatter half of the twentieth century, following the introduction of electronic computers, the restricted three-bodyproblem has been used as a mathematical model to express
real dynamical phenomena in the real world. Examples are
the studies on instability of asteroidal motion (e.g. Wisdom,1983), origin of planetary spins (e.g. Tanikawa et al., 1989,
1991), collision probability of planetesimals (e.g. Nakazawa
et al., 1989a,b), stability of planetary ring systems (e.g.
Hénon and Petit, 1986; Petit and Hénon, 1987), or impossibility of the capture of satellites (e.g. Tanikawa, 1983).
Hénon (1965) was one of the first group of people tohave started the numerical check of the KAM theory andthe numerical study of the stability of periodic orbits, using
the restricted three-body problem as a tool. It seems obvious that Hénon followed the tradition of Poincaré, treating
the restricted three-body problem as a representative example of dynamical systems. In addition, Hénon was among
the first generation of so-called “hackers”. Since Hénon, therestricted three-body problem has been the target of research, as well as a good testbed, for electronic computers.
V. Szebehely, W. H. Jefferys, V. V. Markellos, and manyother people joined these researches at the frontier.
In addition to the numerical study of the restricted
three-body problem, numerical integrations of the generalthree-body problem also started along with the popularization of electronic computers. Celestial mechanists allover the world began to take part, such as S. Aarseth (UK),J. P. Anosova (USSR), R. Broucke (USA), S. Mikkola(Finland), D. C. Heggie and P. Hut (UK), K. Zare
(USA), V. Szebehely (USA), G. Contopoulos (Greece),
S. Yabushita (Japan), or K. Tanikawa (Japan). There area lot of variations within the general three-body problem:the collinear problem, planar and spatial isosceles problem, general planar problem, hierarchical problem, free-fallproblem, and many others. For astronomical application,encounters of single and binary stars are important, andthe timescale for the disintegration of such systems has
been studied (e.g. Heggie, 1975; Heggie and Hut, 1993).
The research on stability of the hierarchical three-bodyproblem that was started by Harrington (1975) may beclosely related to the study of extrasolar planetary systems.

70

Takashi Ito and Kiyotaka Tanikawa

Numerical study of the three-body problems have also
demonstrated how difficult accurate numerical integrationsof gravitational N -body problems are, even when N is
not large. For example, Szebehely and Peters (1967) numerically solved the three-body Pythagorean problem inwhich the three bodies initially fixed on the vertices of a
Pythagorean triangle with edges of lengths 3, 4, and 5,
providing a movie of the results as well. People looking at
this movie understood the difficulty in predicting final motions of the three-body problem with many close encounters. Miller (1964) numerically integrated back and fourtha few-body problem and showed the deterioration of theaccuracy when the system experiences a close encounter.
He measured the accuracy by the phase distance of thesystem at t = 0 and the system returned to t = 0 after some excursions in the phase space. Miller pointed outthat Hamiltonian systems are exponentially unstable, andthat this causes unreliability of numerical integrations. At
nearly the same time, international comparison of the accuracy of numerical integrations using the same 25-bodyproblem has been carried out (Lecar, 1968; Hayli, 1970).
Nine of eleven participating groups integrated “the standard” 25-body problem beginning with the free-fall initialconditions until ∼ 2.5 crossing times. It has been shownthat depending on computers, different bodies escape fromthe system at different times. Further, the number of escapers differs from computer to computer. The problem raisedin this comparison is still open. Another comparison using
32 bodies was performed (Miller, 1971a,b), which raised asimilar question (Aarseth and Lecar, 1975). These difficulties in numerical integration of N -body systems have led us
to the development of more sophisticated numerical techniques, especially those of regularization (see Section 5.1).
3.3. Hierarchical triple systems
It might not be wrong to say that most of the threebody systems in the universe have hierarchical structureswith a high contrast of component masses. In these systems,
declaring a clear definition or a clear criterion of stability
and instability is essential.
So far there are three different stability criteria in theliterature relevant to hierarchical triple systems. First: ahierarchical system is stable if there is no escape of a body(e.g. Anosova, 1996). The stability in this sense is called theE-stability. Second: a hierarchical system is stable if there is
no change in configuration (e.g. Szebehely and Zare, 1977;
Kiseleva et al., 1994a,b). Here the configuration is defined
as being changed when the partner of the inner binary is replaced, or when the hierarchical configuration itself is lost.
This stability is called the CN-stability. Third: a hierarchical system is stable if there are neither secular changes norlarge variations in semimajor axis and/or eccentricity (e.g.
Harrington, 1975; Black, 1982; Donnison and Mikulskis,
1992). The stability in this sense is called the T-stability.

Among the three stability criteria, the T-stability is,
as Harrington himself admitted, subjective and rather apoorly specifying one. For example, Black (1982) considered that the upper limit of the variations of semimajor axis
should be 10%. He thought that systems with larger variations would be E-unstable. The value of 10% has no theoretical background. Hence the criterion may change from
author to author, as Black (1982) or Pendleton and Black(1983) indicated.
The CN-stability is a better-defined concept since twobody energies of pairs can be estimated from the orbital
information. However, the change of the three-body configuration that this criterion supposes might be slightly too
drastic, judging from that fact that there exist triple stellar systems such as the CH Cyg system (e.g. Mikkola andTanikawa, 1998a) that show cataclysmic variation presumably without a drastic change in three-body configuration.
We need another, milder criterion.
Here we propose yet another stability criterion: A hierarchical system is stable when there occurs neither collision of components nor escape of a body. Hereafter wecall the stability in this sense as C-stability. Suppose that acollision takes place between the inner binary components.
Then, the C-stability criterion distinguishes prograde andretrograde hierarchical systems.
Checking the C-stability numerically is easy, if we give
finite non-zero radii to the components, or if we use a numerical technique developed by Tanikawa et al. (1995). It
is conceivable that near-collision orbits develop changes in
configuration. Therefore, we cannot say in general which of
CN- and C-stabilities is stronger.
To organize the relationship between E-, C-, and Tstabilities, consider the sets of E-, C-, and T-stable trajectories in phase space for a given combination of masses.
Evidently, we have T ⊂ C ⊂ E. T is a proper subset of C.
On the other hand, C and E may have common boundaries,
i.e., some C-unstable systems may also be E-unstable.
There is a series of numerical studies for the stabilityof hierarchical triples initiated by Harrington (1975, 1977).
His criterion is that there is no secular change in a and e ofthe triple. Let q2 be the periastron of the outer binary anda1 be the semimajor axis of the inner binary. Three masses
are m1 , m2 , and m3 , among which m3 is of the outer body.
Harrington proposed a sufficient condition for stability ofthe planar three-body problem in the formµ¶
q2(q2 /a1 )0m3≥
log10 1 +,(35)
a1
log10 1.5
m1 + m2
where (q2 /a1 )0 is the least stable value for the equal-mass
triples. (q2 /a1 )0 has different values depending on the direction of rotation of the outer body. Following Harrington,
Szebehely and Zare (1977), Black (1982), Donnison and
Mikulskis (1992) derived better criteria. The newest criterion is due to Mardling and Aarseth (1999). Their criterion
for the planar prograde motion isTrends in 20th Century Celestial Mechanics
q2
≥C
a1

·µ
1+

m3
m1 + m2¶

1 + eout√
1 − eout

¸ 52
,

(36)

where eout is the eccentricity of the outer binary. The constant C ' 2.8 is determined empirically.
As we will mention later (such as in Section 4.2), recent
investigations into the long-term stability of planetary orbits have exhibited a large volume of chaotic orbits thatare actually stable in the sense that they survive billionsof years without resulting in collisions between planets.
These results suggest the existence a class of orbits that arechaotic as well as C-stable in N -body systems with N > 3.

3.4. New frontiers of the three-body problems
The three-body problem has long been an attractive but
dangerous subject for students. This is because the threebody problem has quite a simple setting that even a highschool student can understand, and it appears relatively
easy to find something new even for a beginner. However,
the three-body problem has been investigated in such agreat detail for so many years that it is actually very difficult for a newcomer to quickly obtain anything new in thisfield. Thus, experienced scientists often give an alert to newstudents, “Never be attracted to the three-body problem.
It is too dangerous.”
Although purely theoretical aspects of the three-bodyproblem are quite formidable for students, the advent offast computers allows us to observe and analyze the everincreasing complexity of the phase space created by thethree-body problem. There is an opinion that the generalthree-body problem can only survive as a toy for eccentric scientists. What is the purpose of solving the motion
only of three particles? Is there any relation to the otherfield of sciences even if one knows how three bodies move?
Fortunately, however, the three-body problem has rather
strong connections to other fields. From a practical pointof view, the three-body problem has yielded many new
applications to astrophysics. Some examples are given in
Section 3.5. From a theoretical point of view, connections
between the three-body problem and general Hamiltoniansystems are rather strong. Indeed, the equations of motionof the three-body problem are among the simplest and themost meaningful examples of differential equations appearing in astronomy, physics, and mathematics. Moreover, theequations are non-integrable: the set of solutions may havea very complicated and intriguing structure. Theoreticaland numerical techniques developed in general Hamiltonian systems and in general differential equations have beenchecked of their effectiveness in terms of the three-bodyproblem. Conversely, techniques developed in the threebody problem have been extended to a broader class of differential equations. This tendency will continue in the new
century.

71

One promising approach to the three- and several-bodyproblem is to apply symbolic dynamics. It might be called“numerical symbolic dynamics”. In symbolic dynamics, anorbit is replaced by a finite, infinite, or bi-infinite sequenceof symbols that represent special events along the orbit.
Typically, a special event corresponds to the occupation of
special place in phase space. The number of special places
is equal to the number of symbols necessary to construct
a symbol sequence. The origin of symbolic dynamics dates
back to Hadamard (1898) and the textbook of Birkhoff(1927); the method was soon developed in a paper entitledSymbolic Dynamics by Morse and Hedlund (1938). Alekseev (1968a,b) applied symbolic dynamics to the threebody problem (see also Alekseev, 1981). Alekseev used symbolic dynamics to prove the existence of oscillatory orbits
in the general three-body system. On an oscillatory orbitof the three-body system, one of the bodies is repeatedly
ejected from the remaining two bodies, and the amplitudeof the ejection distance becomes unbounded in time.
Numerical symbolic dynamics has recently been used
to unveil qualitative properties of three-body orbits. Examples of successful applications of numerical symbolic
dynamics to the few-body problem are Tanikawa andMikkola (2000a,b), Zare and Chesley (1998), Sekiguchi andTanikawa (2004), Sano (2004), and Saito (2005). Tanikawa
and Mikkola (2000a,b), in the one-dimensional three-bodyproblem, explicitly used two types of binary collision and a
triple collision along orbits as three symbols, and replaced
the orbits with sequences of these three symbols. Using
general techniques of symbolic dynamics such as making
the transition graph, and searching inadmissible symbol
sequences and periodic sequences, they succeeded in obtaining qualitative properties of orbits and phase spacestructure of the one-dimensional three-body problem withequal mass. These studies happened to clarify the relation
between the phase space structure and triple collisions.
Collision orbits are on the stable or unstable manifolds ofthe fixed points on the collision manifold in the extended
phase space. It has been clearly illustrated that stable andunstable manifolds form the backbone of the phase space.
In addition, triple collisions have been easily identified.
This method has been applied to the symmetric four-bodyproblem (Sekiguchi and Tanikawa, 2004) and the collinear
Helium nucleus (eZe) problem (Sano, 2004). The technique is expected to be applicable to the planar three-bodyproblem with zero angular momentum.
At the turn of the century, unexpected and surprising
news spread through the celestial mechanics community:
A proof of the existence of a new type of periodic orbit,
the figure-eight (“8”) orbit (Chenciner and Montgomery,
2000). This solution has been called “a choreographic solution” (Simó, 2002). This is because bodies dance along
a single track on the plane in this type of orbits, with all
the bodies having the same mass with time shift equal to
1/3 of the period. The figure-eight solution lives on the

72

Takashi Ito and Kiyotaka Tanikawa

hypersurface of zero angular momentum, and it is stable.
Chenciner and Montgomery (2000) noted that the space oforiented triangles in the plane formed by three bodies up totranslation and rotation, if normalized by the inertial moment, is a 2-sphere. They called this the “shape sphere”.
The north and south poles correspond to equilateral triangles (Lagrange configurations), the equator corresponds todegenerate triangles (collinear configurations), and threeequally separated meridians correspond to isosceles triangles. Cross points of these meridians and the equator are
either collinear Euler configurations Ei , or binary collisions
Ci . The shape sphere is divided into twelve similar regions
by three meridians and the equator. A path in one of thetwelve regions, starting at the Euler configuration and ending at an isosceles configuration, is copied to the remaining eleven regions with the aid of symmetry. Upon connection, these twelve paths form one closed path on the shapesphere. Chenciner and Montgomery (2000) first showedthat the action-minimizing path has zero angular momentum. Then they showed that the action-minimizing pathdoes not experience collision. Finally they showed that theconnected path actually corresponds to a closed path in theconfiguration space and to the figure-eight solution.
Simó (2002) studied numerically the figure-eight orbit
and its neighboring orbits. He took a surface of section neara phase space point where the figure-eight orbit becomes
collinear. He integrated 36 million orbits around there until one of collision, long ejection to some distance, and thecompletion of 5,000 revolutions of the figure-eight orbit occurs. He plotted the initial points of the surviving orbits,
and discussed various properties of these orbits. There wasno choreographic orbit other than the figure-eight orbit on
this local surface of section. Then Simo changed his strategy. He extended the figure-eight orbit problem in two different directions. One direction was to try different mass
combinations. He found periodic orbits, but they were nonchoreographic. The other direction was to give non-zero angular momentum. In this case, Simo found a lot of choreographic orbits: the figure-eight rotates and closes after somerevolutions.
Seven years before Chenciner and Montgomery (2000),
C. Moore (1993) had numerically found the figure-eight
periodic solution of the three-body problem. His ideas are
remarkable in three-fold aspects. First, Moore considered
the planar N -body problem and construct braids in 2+1
space-time. Second, Moore embedded the gravitational N body problem in the dynamics with potential
V =

N X
N
X
i

Vij ,

α
Vij = Ami mj rij(−2 ≤ α ≤ 2).

(37)

j

The third idea is that Moore put the problem in the extrema principle. He observed that the system is integrable
for α = 2, and that the bodies move around the center ofmass harmonically. So, all braids are trivial when α = 2,

whereas all braid types exist when α = −2. For intermediate values of α, some braids exist, while others do not: In
the word of the N -body problem, some periodic orbits exist, while others do not. Moore confirmed that the changes
in accordance with α is monotone. To obtain periodic orbits, Moore first specified the braid type, and started witha test orbit that realizes this braid type. His action reads¶2
Z µ 2
d x F (x)∂τ S = −m
−
dt.
(38)
dt2
m
This becomes maximal when the acceleration on thepath is derived from force. Moore moved from one path toanother in order that ∂τ S increases. If a collision takes place
before the action maximizing path is attained, then there isno periodic orbit corresponding to the given braid, because
the braid type changes through collision. In this way, Moore
found the figure-eight orbit for α < 2. He found other kindof orbits as well. However, the figure-eight orbit is the only
choreographic orbit (see Table I of Moore (1993)). We mightwant to look for other choreographic orbits using Moore’s
method.
The discovery of the figure-eight orbit triggered an explosion of the study of choreographic (or dancing) N -bodysolutions (note that these names were given by C. Simó).
Together with many other publications, works by T. Fujiwara and collaborators are worth being mentioned (Fujiwara et al., 2003a,b, 2004a; Fujiwara and Montgomery,
2003). They found very interesting and simple propertiesof the zero angular momentum solution of the three-bodyproblem. One is the three-tangents theorem which says thatin the center of gravity system, three tangents to the directions of motion of mass points always cross at a point or at
infinity. Similarly, in the center of gravity system, if the moment of inertia is constant, then normal lines to the velocity
vectors at bodies meet at a point or at infinity. Kuwabara
and Tanikawa (2006) extended the three-tangents theorem
and three-normals theorem to the planar three-body problem with non-zero angular momentum. Fujiwara and his
collaborators also showed that the form of the figure-eightsolution is a lemniscate when the force has an additional
repulsive term (Fujiwara et al., 2003a). One other property
they have shown is the convexity of the figure-eight orbit(Fujiwara et al., 2003b). Though the exact functional formof the figure-eight orbit is still unknown, the figure-eight
orbit turned out to have a non-constant moment of inertia(Fujiwara and Montgomery, 2003). Fujiwara et al. (2004a)suggested a direction of general study of the planar threebody problem with the aid of a three-tangents theorem.
Choreographic N -body solutions and related topics
will remain an important and exciting subject in the earlytwenty-first century. A lot of solutions for the N -bodyproblem of this kind have been found, and continuously
being found. Mathematical structures of the solutions andtheir positions in phase space are attracting many scientists, which may affect all other branches of mathematical

Trends in 20th Century Celestial Mechanics
sciences as the three-body problem has done. It has beenfound that embedding the gravitational N -body problem
in a sequence of different force laws characterizes the gravitational N -body problem itself. Works by Fujiwara andMoore are good examples. We will find further interesting issues of the choreographic N -body problem in thisdirection.
3.5. Applications to actual planetary and stellar systems
The three-body problem shows even larger importancewhen it is applied to actual planetary or stellar systems.
Since there have been so many applications of the restricted
three-body problem to astronomy, it is impossible to cover
the whole range of applications in this manuscript. Someof them are dispersed in various sections. As for the future direction of the application, the restricted three-bodyproblem will be applied to newly discovered multiple stellar systems, such as extrasolar planetary systems or stars
around binary blackholes. As for the general three-bodyproblem, the authors would like to point out two promising directions of study. The Kozai mechanism is involvedin both of them.
The first one is the study of long-periodic variable stars,
for example, a very strange variable, CH Cyg (Mikolajewski et al., 1990). In the nineteenth century, CH Cyg used
to be a standard star. Later it became a semi-regular star.
Then, it suddenly started its prominent activity in 1963.
Since then it has been almost constantly emitting radio,
optical, and even X-ray jets. Its brightness changes by a
few magnitude. Actually CH Cyg was recognized as a binary star (Yamashita and Maehara, 1979), but this kind of
sudden start of stellar activity cannot be simply explainedby a binary star model. Through infrared spectral observations, CH Cyg turned out to be a triple (Hinkle et al.,1993). As a binary, it has a period of around two years. But
as a triple, it has a period of 14 years. How can a triple likethis maintain this long-periodic activity, apparently longer
than a hundred years? The first trial to solve this problem
was done by Mikkola and Tanikawa (1998b). They proposed
the Kozai mechanism as a dynamical model to explain the
behavior of CH Cyg. The Kozai mechanism has increasing
importance in the various fields of astronomy. In the hierarchical three-body problem, the Kozai mechanism has along periodicity called the Kozai cycle. The period TKozai
of the Kozai cycle is numerically fitted to a formula as follows (Aarseth, 2003):
µ
¶
2
¢3
Tout
1 + qout ¡TKozai =1 − e2out 2 g(ein , ωin , ψ)
(39)
Tin
qout
where Tout and eout are the orbital period and eccentricity of the outer binary, Tin , ein , ωin are the orbital period,
eccentricity, and argument of pericenter of the inner binary, ψ is the orbital inclination, qout = m3 /(m1 + m2 ),and g is a function of order unity depending on its argu-

73

ments. As is seen in (39), TKozai is longer than the orbital
period of the outer component of a triple by a factor often or more. When the orbital planes of the three bodies coincide, the eccentricity of the inner binary increases.
Hence at the pericenter, the supergiant component of CH
Cyg may fill its Roche lobe. Overflowed material accretes
onto a companion—neutron star or white dwarf. During
this stage, CH Cyg repeats activity in two year periodicity.
On the other hand when the orbital planes of two binaries
differ, the inner binary becomes a detached system without
any activity. This is the basic explanation for the strange
behavior of CH Cyg.
Related to this study is the direct search for long periodic variables. For example, Fujiwara et al. (2004b) compared seven historical star catalogues from Almagest to
Uranometria, and tried to find long periodic variables. Fujiwara and her colleagues found 57 Peg and 19 Psc as candidates (Fujiwara and Hirai, 2006). Dynamical work (Mikkola
and Tanikawa, 1998b) and astrophysical work (Fujiwara
and Hirai, 2006) may give impetus to the study of variable
stars that have been considered irregular until now and left
unanalyzed. In this respect, a systematic dynamical studyof three-dimensional hierarchical triple systems is strongly
desired.
The second example of the promising directions of thegeneral three-body problem is the study of blackholes in
the center of galaxies which recently draw a large attention.
Merged galaxies may contain two or more super-massive
blackholes in their center. Then, for the evolution of theentire galaxy, the evolution of central super-massive blackholes is crucial. There arises a relativistic three-body problem. If mergings are numerous, the frequency of triple blackholes is high. There can be a considerable change in thefate of the central part of a galaxy whether it has a blackhole binary or a triple blackhole. Merging timescale of abinary blackhole is said to be longer than the age of theuniverse. Iwasawa et al. (2006) studied this problem. They
found that three-dimensional three-body problem of supermassive blackholes significantly shorten the timescale. Indeed, as in the case of stellar triples, there are configurations in which the Kozai mechanism is effective. In thesecases, when the orbital planes coincide, the inner binary
may emit gravitational waves at their pericenter and theirorbit shrinks. This substantially speeds up the orbital evolution of the inner binary. Then, the binary merges during
the timescale shorter than the Hubble time.
4. Solar system dynamics
In the latter half of the nineteenth century, discoveriesof major and minor planets were the most significant issue
in celestial mechanics. Also, earlier in the twentieth century, the existence of another planet was predicted, and theeffort to find it bore fruit when C. W. Tombaugh made his
discovery of Pluto in 1930. This kind of activity was directly

74

Takashi Ito and Kiyotaka Tanikawa

connected to the development of analytical perturbation
theories, now gradually being replaced by numerical integrations. In what follows, we describe several major aspectsof the research in solar system dynamics: progress of analytical perturbation theories, long-term dynamical stabilityof the solar system planetary motion, dynamical studies ofplanet formation, dynamical studies of small bodies, andplanetary rotation. Many important topics related to solar
system dynamics are not described here due to authors’ incompetence such as dynamics of near-Earth asteroids andthat of ring-satellite systems, motion and control of artificial satellites and spacecrafts, theory of planetary tides,
interaction between solar system bodies and protoplanetary nebula, and orbital determination. The description
about bodies at the outer edge of the solar system (transNeptunian objects including the Kuiper Belt, the scattereddisk, and the Oort Cloud objects) is extremely limited, andshould be expanded far more. The subsection 4.5 for the dynamics of planetary rotation should have been longer andmore extensive, incorporating the latest observation results
and their dynamical interpretation. Also, we do not havea specific section about extrasolar planetary systems thatoffer us an extremely large variety of testbeds for solar system dynamics. This is because we think the most methodology of extrasolar planetary systems research is alreadyembedded in the conventional studies of solar system bodies, many of which we have described in this manuscript.
4.1. Progress of perturbation theories
In the era of Laplace or Lagrange, solar system planetary motion and its stability were the central issues of astronomy. The perturbation equations of orbital elements
derived by the method of variation of constant are now
called Lagrange’s planetary equations. Also, a secular perturbation theory up to the lowest order in planetary eccentricity and inclination is called the Laplace–Lagrange
method. In the nineteenth century, urged on by the adventof accuracy in astronomical observation of planetary motion, greater precision and reliability were required of analytical perturbation theories in order to describe and predict planetary motion. Many famous astronomers such as
U. J. J. LeVerrier, S. Newcomb, and G. W. Hill devoted alarge part of their lives to construct accurate perturbation
solutions of planetary motion.
In addition to planetary motion, lunar motion had
drawn the attention of scientists working on analytical perturbation theories until the end of the nineteenth century.
Newton showed in Principia that the principal periodic
perturbation in lunar motion can be ascribed to the actionof the Sun. For constructing the accurate ephemeris of theMoon, the method of the variation of arbitrary constants
was studied in the eighteenth century by Euler, and waslater discussed in more detail by Poisson in the mid nineteenth century. This method, which had been an important

foundation of canonical perturbation theories in the twentieth century, was further developed by Delaunay in the1860s. Delaunay completed the application of his method
to the solution of the lunar “main problem” where only theMoon, the Earth, and the Sun are taken into account andconsidered as point masses. After Delaunay in the 1880s,
Hill provided the calculation of the terms arising fromthe non-sphericity of the Earth, and Radau calculated theplanetary perturbations in the 1890s. All this work createda strong basis for the analytical perturbation theories thatwere developed later in the twentieth century. Brouwerand Hori (1962) provide a concise summary of the historyof lunar perturbation theories.
Based on the legacy of the nineteenth century, aclassical planetary perturbation theory by Brouwer andvan Woerkom (1950) was developed as an extension
and a rebuild of Hill’s theory (Hill, 1897). Brouwer andvan Woerkom (1950) included the terms up to the sixth
degree in the eccentricities and inclinations of Jupiter andSaturn. They also considered terms related to the “great
inequality” (i.e. the near equality to the 5:2 mean motion
resonance) between Jupiter and Saturn up to the second
order in planetary masses.
Bretagnon (1974) and Bretagnon and Francou (1992)
extended the classical theory of Brouwer and van Woerkom(1950). Bretagnon’s analytical theory includes the terms
up to the fourth degree in planetary eccentricities andinclinations, and up to second order in planetary masses.
Based on the preceding research, Laskar developed a very
high-accuracy semi-analytical secular perturbation theory(Laskar, 1985, 1986, 1988, 1990). Laskar’s theory includes
terms up to the sixth degree in planetary eccentricities andinclinations, and up to second order in planetary masses.
Laskar’s secular perturbation theory is an extension ofDuriez’s (1977; 1982) preceding work concerning the four
jovian planets. Laskar solved the equations of motion of
secular planetary perturbation, decomposed the solutions
into Fourier components, and obtained the frequencies andthe amplitudes of major periodic oscillations. This theory
now serves as a standard in the area of classical secularplanetary theories.
After the epochal launch of Sputnik in 1957, the application of analytical perturbation theories suddenly expanded toward the motion of artificial satellites. Startingfrom the three famous papers (Garfinkel, 1959; Kozai, 1959;
Brouwer, 1959) published consecutively in the same issue ofthe same journal, a large number of theories have been published about the motion of artificial satellites around theEarth, composing the fundamental theoretical basis of themore than 5,000 artificial satellites that have been launchedsince Sputnik. Kozai (1962a) and Kinoshita (1977b) established very high-accuracy theories of the motion of artificial satellites. Analytical perturbation theories for higheccentricity and high-inclination satellite orbits have beendeveloped by E. Brumberg and V. A. Brumberg, exploit-

Trends in 20th Century Celestial Mechanics
ing elliptic functions (Brumberg et al., 1995; Brumberg andBrumberg, 1995, 2001). New theories are still being published in the field of satellite control and astronautics, withplenty of good and detailed textbooks dealing with this issue such as those by Taff (1985), Battin (1987), or Montenbruck and Gill (2000).
One class of the newer perturbation methods is based on
Hamiltonian dynamics: The canonical perturbation theory.
In a general problem with small perturbation, suppose thata simple, time-independent Hamiltonian H with n degrees
of freedom be written as
H(θ1 , θ2 , · · · , θn , J1 , J2 , · · · , Jn )= H0 (J1 , J2 , · · · , Jn )+Hpert (θ1 , θ2 , · · · , θn , J1 , J2 , · · · , Jn )(40)

with Hamilton’s canonical equations
dθi
∂H
=
,
dt
∂Ji

dJi
∂H
=−
dt
∂θi

(i = 1, · · · , n)(41)

where θi and Ji are the angles and actions of the unperturbed (and integrable) system, H0 . The perturbed part of
the Hamiltonian Hpert is assumed to be much smaller thanthe unperturbed part (i.e. Hpert /H0 = O(ε)).
In canonical perturbation theories, we apply canonicaltransformations to the system described by (40) and (41),trying to reduce a system with terms of non-integrable nature (Hpert ) to an integrable form. This procedure is equivalent to eliminating short-periodic terms (i.e. angle-like variables) from the Hamiltonian by averaging processes. Oureventual goal is to make the system depend only on actions
which are constant. More specifically, we want to convert
the Hamiltonian (40) by a canonical transformation into aform asH ∗ (−, θ2∗ , · · · , θn∗ , J1∗ , J2∗ , · · · , Jn∗ )(42)

where “−” denotes that a variable is eliminated by an averaging process (in this case, θ1 is supposed to be eliminatedfrom the system), and the superscript ∗ denotes that thissystem has been canonically transformed. Canonical equations using a time-like variable t∗ will be
dθi∗
∂H ∗
=
,
dt∗
∂Ji∗

dJi∗
∂H ∗
=
−
,
dt∗
∂θi∗

1, · · · , n. Similarly, we eliminate θ2∗

for i =
another averaging process as(43)∗from H by

H ∗∗ (−, −, θ3∗∗ , · · · , θn∗∗ , J1∗∗ , J2∗∗ , · · · , Jn∗∗ )= H0∗∗ (J1∗∗ , J2∗∗ , · · · , Jn∗∗ )∗∗+Hpert(−, −, θ3∗∗ , · · · , θn∗∗ , J1∗∗ , J2∗∗ , · · · , Jn∗∗ )

with
dθi∗∗
∂H ∗∗
=
,
dt∗∗
∂Ji∗∗

dJi∗∗
∂H ∗∗
=
−
,
dt∗∗
∂θi∗∗

for i = 1, · · · , n. If this series of canonical transformations
remains possible until we have eliminated all the angle-likevariables, we will end up with obtaining a Hamiltonian that
contains only actions as
H ∗···∗ (−, · · · , −, J1∗···∗ , J2∗···∗ , · · · , Jn∗···∗ )= H0∗···∗ (J1∗···∗ , J2∗···∗ , · · · , Jn∗···∗ )∗···∗+Hpert(−, · · · , −, J1∗···∗ , J2∗···∗ , · · · , Jn∗···∗ )(46)

with a set of canonical equations of motiondθi∗···∗∂H ∗···∗=,
dt∗···∗∂Ji∗···∗dJi∗···∗∂H ∗···∗=−,
dt∗···∗∂θi∗···∗(47)

for i = 1, · · · , n. Then, the system described by (40) and(41) would be integrable: The second equation of (47) tells
us that all Ji∗···∗ are constant; they are true actions of H ∗···∗ .
Also, from the first equation of (47), it turns out that theangles θi∗···∗ are just proportional to the time-like variablet∗···∗ .
In the field of general canonical perturbation methods,
we should mention the epochmaking transformation theoryof G. Hori (Hori, 1966, 1967, 1970, 1971) which utilizes the
Lie transformation and accomplishes an explicit canonicaltransformation. This explicitness is the most significant advantage of Hori’s method in its practical use, making it far
more advanced than the traditional von Zeipel method (cf.
von Zeipel, 1916). Later, Deprit (1969) proposed anotherform of the explicit canonical transformation theory equivalent to Hori’s.
The explicitness of Hori’s general perturbation theory may be summarized as follows (Hori, 1966). First, let
us assume that (ξj , ηj ) is a set of 2n canonical variables
with j = 1, 2, · · · , n. Let f (ξ1 , ξ2 , · · · , ξn , η1 , η2 , · · · , ηn ) andS(ξ1 , ξ2 , · · · , ξn , η1 , η2 , · · · , ηn ) be arbitrary functions of(ξj , ηj ). Here we introduce a differential operator DSn as:
DS0 f = f,

= H0∗ (J1∗ , J2∗ , · · · , Jn∗ )∗(−, θ2∗ , · · · , θn∗ , J1∗ , J2∗ , · · · , Jn∗ )+Hpert

75

(44)

(45)

DS1 f = {f, S} ,(48)

DSn f = DSn−1 (DS1 f ),(n ≥ 2)

with the Poisson brackets {, }. The theoretical basis thatHori (1966) relies upon is the fact that was discovered by Lie(1888): A set of 2n variables xj , yj (j = 1, · · · , n) defined byf (x1 , · · · , xn , y1 , · · · , yn )
∞
X
εn n
=
D f (ξ1 , · · · , ξn , η1 , · · · , ηn ),
n! S
n=0

(49)

is canonical if the series in the right-hand side of (49) converges, where ε is a small constant that is independent ofξj and ηj . When f (ξ1 , · · · , ηn ) = ξj or f (ξ1 , · · · , ηn ) = ηj ,from (49) we get
x j = ξj +yj = ηj −

∞
X
εn
n=1
∞
X

n!

DSn−1

∂S
,
∂ηj

εn n−1 ∂S
D.
n! S ∂ξj
n=1

(50)

76

Takashi Ito and Kiyotaka Tanikawa

What is most remarkable in (50) is that the new variables (xj , yj ) and the original variables (ξj , ηj ) are notmixed: xj and yj are only on the left-hand sides, while theright-hand sides of (50) are the functions only of ξj and ηj .
We can compare (50) with the conventional form of canonical transformation as
x j = ξj + ε

∂ S̃
,
∂yj

yj = ηj − ε

∂ S̃
,
∂ξj

(51)

where S̃(ξ1 , · · · , ξn , y1 , · · · , yn ) is a function of the mixed
set of original and new variables, ξi and yi . When ε is smallenough, S(ξ1 , · · · , ξn , η1 , · · · , ηn ) and S̃(ξ1 , · · · , ξn , y1 , · · · , yn )are related to each other by the following equation:nε X ∂ S̃ ∂ S̃
2 j=1 ∂ξj ∂ηjÃnn
ε2 X X ∂ S̃ ∂ S̃ ∂ 2 S̃+
12 j=1∂ξj ∂ξk ∂ηj ∂ηk

S = S̃ −k=1+4∂ S̃ ∂ S̃ ∂ 2 S̃∂ S̃ ∂ S̃ ∂ 2 S̃+∂ξj ∂ηk ∂ηj ∂ξk∂ηj ∂ηk ∂ξj ∂ξk+O(ε3 ).

even if this ultimate situation arises, someone will have tounderstand the meaning of the output from the package,
lest the package be a total black box. Thus, analytical
perturbation theories should be passed on from one toanother as an important traditional art in solar system
dynamics. It might not be necessary that a huge number ofpeople are involved, but there must be a certain number ofpeople (or groups) that dedicate themselves to this kind ofanalytical work. We must also keep in mind that work inthis field generally requires a much longer time to publish
than does work in the field of numerical studies. Hence, it
might not be very easy for us to encourage young students
to get involved with time-consuming work on analytical
perturbation theory, especially in the current “Publish or
perish” atmosphere of the academic world. However, wemust never let this field of research perish.
4.2. Long-term dynamical stability of planetary motion

!

(52)

In (51), the relationship between the new variables(xj , yj ) and the original variables (ξj , ηj ) is implicit, while
it is explicit in (50). Obviously, the explicit relationship(50) is far more useful in many of the realistic problemsof actual dynamical application of canonical perturbation
theories.
One of the major obstacles for traditional analytical
perturbation theories is that they are hard to use for treating the orbits that cross each other, such as those of nearEarth asteroids. This kind of orbit introduces a singularity in the Hamiltonian. The averaging of the equations of
motion that is traditionally used to compute secular perturbations is not definable. This singularity problem hasrecently been overcome in a series of efforts in a secular theory by G. F. Gronchi (Gronchi and Milani, 1999a,b, 2001;
Gronchi and Michel, 2001; Gronchi, 2002). He showed that
it is possible to define generalized averaged equations of
motion for crossed orbits. These are obtained using Kantorovich’s method of extraction of singularities. The modified distance used to approximate the singularity is thesame one used by Wetherill (1967).
Current perturbation theories, though they are still
called “analytical,” depend so strongly on computer algebra that we are better off calling them “computationallyanalytical”. Indeed in the near future, we may be able toanticipate the development of a “fully automatic analyticaldevelopment software” package, to be used as a development tool exclusively for perturbation theories in celestialmechanics. When using this kind of software, all we haveto do would be to specify orders and parameters for thedevelopment of perturbation. Then, computers would provide us with solutions as accurate as we desired. However,Long-term dynamical stability of solar system planetary motion has been a central problem of celestial mechanics for centuries (cf. Tremaine, 1995). This problem
has been debated over several hundred years since Newton’s era, and has attracted many famous mathematicians
over the years. This problem has also played a significant
role in the development of nonlinear dynamics and chaostheory. However, we do not yet have a definite answer tothe question of whether our solar system is stable or not.
This is partly due to the fact that the definition of “stability” is vague when it is used in relation to the problem ofplanetary motion in the solar system. Nobili et al. (1989)
fluently expressed this confusion in the introductory partof her paper:“When the stability of our solar system is discussed,
two objections often arise. Firstly, this problem has beenaround for too long, never getting to the point of stating clearly whether the system is stable or not; the few
definite results refer to mathematical abstractions suchas N -body models and do not really apply to the real
solar system. Secondly, the solar system is macroscopically stable—at least for a few 109 years—since it is still
there, and there is not much point in giving a rigorousargument for such an intuitive property.”
The traditional approach to dealing with this problem
was to resort to analytical perturbation theories. But the recent progress of computer technology pushes us toward numerical studies. Numerical research on this problem started
in the 1950s when electronic computing machines became
available. Because of the very limited computational resources, the object of numerical integrations then was focused on the five outer planets 3 (Jupiter, Saturn, Uranus,3

As is already well known, Pluto has been “demoted” from the category of planet to dwarf planet by International Astronomical Union

Trends in 20th Century Celestial Mechanics
Neptune, and Pluto) whose orbital timescales are much
longer than those of the four inner planets; Mercury, Venus,
the Earth, and Mars. The first famous numerical computation of planetary orbits was by Eckert et al. (1951) who did a
350-year simulation of the four outer planets (Jupiter, Saturn, Uranus, and Neptune) using a large mainframe computer. This was extended by Cohen and Hubbard (1965)
to 120,000 years, and by Cohen et al. (1973) to one million
years. These integrations agreed well with the perturbation
calculations of Brouwer and van Woerkom (1950), showing
quasiperiodic behavior for the four major outer planets.
Later, Kinoshita and Nakai (1984) extended the length of
the integration of the five outer planetary motion (Jupiter,
Saturn, Uranus, Neptune, and Pluto) to five million years.
One of the longest numerical integrations using conventional integration schemes on a general purpose computer
was the LONGSTOP 1B integration done by Nobili et al.
(1989) which ran for 100 Myr. Consult Roy et al. (1988)
or Milani et al. (1987, 1989) for details of the LONGSTOP
project. Their integrations followed the mutual interactionsof the five outer planets, but also included the secular effectof the four inner planets and the effect of general relativity.
Integrations of the five outer planets for periods up to845 Myr was done with a special-purpose computer, theDigital Orrery (Applegate et al., 1986; Sussman and Wisdom, 1988). Sussman and Wisdom (1988) found in theirlong-term integration that the orbital motions of five outerplanets including Pluto are chaotic. The Digital Orrery wasa specialized but programmable high-performance computer designed for the efficient numerical integration of theequations of motion for systems with a small number ofbodies that move in roughly circular orbits. This machineconsisted of one CPU per planet, with all the CPUs arranged in a ring. Wisdom and Holman (1991) extended thelength of the numerical integrations of five outer planetary
motion to 1.1 Gyr using their original high-speed symplectic map. The integration period was further extended to±5.5 Gyr for the five outer planets by Kinoshita and Nakai(1995, 1996) using a PC with a dedicated numerical processor.
When we include four inner planets in the orbital integrations, the amount of computation is multiplied about 150 times: Stepsize of integration must be
as of August 24, 2006 (see http://www.iau.org/ for more detail).
However in this manuscript, our treatment of Pluto goes along how
original publications dealt with this body; Pluto has been treated asa planet in most previous literature. In addition, we are aware of thefact that the definition of “inner planets” or “outer planets” sometimes differs by who uses these terms. Generally in this section, “thefive outer planets” denotes Jupiter, Saturn, Uranus, Neptune, andPluto, “the four outer planets” or “the four jovian planets” denotes
Jupiter, Saturn, Uranus, and Neptune, and “the four inner planets”
means Mercury, Venus, the Earth, and Mars. We also try to list theplanet names explicitly such as “the five outer planets (Jupiter, Saturn, Uranus, Neptune, and Pluto)” or “the four inner planets (Mercury, Venus, the Earth, and Mars)” whenever we think necessary.

77

PJupiter /PMercury ∼ 50 times smaller (where P is the orbital period of planets), and the number of force combinations becomes 9+1 C2 /5+1 C2 = 3 times larger. In addition,
until recently the length of direct integrations was limited
not by the CPU time, but by numerical errors. Milani and
Nobili (1988) concluded that it was impossible to reliably
integrate the orbits of the outer planets for a period of 109
years or more with the computer hardware and conventional algorithm current at that time. The main problem
they found was the limited machine precision, and empirical evidence that the longitude error after n steps is
proportional to n2 . Quinn and Tremaine (1990) proposedseveral corrections to conventional integration algorithms
which have considerably reduced the roundoff error. Quinn
et al. (1991) proposed a high-order multistep symmetricscheme. Quinn et al. (1991) utilized some of the numerical
techniques that were state-of-the-art at the time to makean accurate integration of all nine planets and Earth’s
spin axis for 3.05 Myr into the past and future. Previously,
Richardson and Walker (1989) performed 2 Myr integrations of all nine planets, but they neglected the effects ofgeneral relativity and the finite size of the Earth–Moon
system.
After the publication of Sussman and Wisdom’s (1992)
work that confirmed that the Lyapunov time of the four
inner planets is less than several million years, the boom
in long-term numerical integration seemed to slow somewhat. But thanks to the recent availability of low-cost andhigh-performance PCs, much longer-term numerical integrations are now possible. Two of the longest numericalintegrations currently published are Duncan and Lissauer(1998) and Ito and Tanikawa (2002). Duncan and Lissauer’s(1998) main target was the effect of post-main-sequence solar mass loss on the stability of our planetary system. They
performed many numerical integrations up to 1011 years offour jovian planets’ orbital motion, and several 109 yr integrations of seven planets (from Venus to Neptune). The
initial orbital elements of the planets are the same as thoseof our solar system planets, but the sun’s mass decreases
relative to those of the planets in their experiments. They
found that orbit crossing times of planets, which can be atypical indicator of an instability timescale, are quite sensitive to the degree of the sun’s mass decrease. The reduction of the sun’s mass induces reduction of the instability, probably increasing the mutual Hill radii (also known
as tidal radii; cf. Gladman (1993)) among planets. Whenthe sun’s mass is not so different from its present value,
jovian planets remain stable over 1010 years, or perhaps
more. Duncan and Lissauer (1998) also performed four similar experiments on the orbital motion of the seven planetsfrom Venus to Neptune spanning about 109 years. It seems
that the inner planets (Venus, the Earth, and Mars) also
remain stable during the integration period, maintaining
quasiperiodic oscillations similar to the present one.

78

Takashi Ito and Kiyotaka Tanikawa

Ito and Tanikawa (2002) presented the results of very
long-term numerical integrations of planetary orbital motions over ±5 × 109 -yr timespan, encompassing all nine
planets from Mercury to Pluto. Their numerical model
includes only the classical Newtonian gravitation among
planets and the Sun. Satellites, general relativity, solar massloss, or other non-gravitational forces were neglected. Theirnumerical data shows that the planetary motion seems
quite stable even for a 1010 -year timespan. A closer lookat the lowest-frequency oscillations in their result indicatesthe potentially diffusive character of the four inner planetary motion, especially that of Mercury. The behavior of
Mercury’s orbital eccentricity in their integrations is qualitatively similar to that found in the results of Laskar’s secular perturbation theory. However, there is no apparent secular increase in eccentricity or inclination in any planet’s
orbital elements. Ito and Tanikawa (2002) have also performed a couple of integrations including only the five outerplanets for ±5 × 1010 years. They found that the threemajor resonances in the Neptune–Pluto system have been
maintained over a 1011 -yr timespan, with rigorous stability
in the orbital motion of all five outer planets.
The experience and knowledge gained through the longterm numerical study of our solar system’s planets havebeen inherited by the new generation of researchers of extrasolar planetary systems. The overwhelming number of
discoveries of extrasolar planets (more than 200 as of April
2007. See Butler et al. (2006) for more detail) suggests that
one of our future tasks is to count and categorize the mutual interactions of planets in our solar system as well as inother planetary systems, for as many timescales as possible: a taxonomy of planetary stability and instability. This
will be an important theme of solar system dynamics in thenear future, and will form a basis for the dynamical study of
extrasolar planetary systems. In this sense, a set of numerical experiments by Lissauer et al. (2001) is quite interesting and intriguing. If a planetary mass body were present
in the asteroid belt, the orbits of the four inner planets andthose of the outer giant planets would be more closely coupled. A greater exchange of angular momentum could affectthe stability of the four inner planets. More extensive andmore comprehensive studies in this line should be done fora larger range of parameters of various planetary systems.
4.3. Dynamics in planet formation study
Equipped with the evidence that the current solar system planetary motion is stable for a very long time, research
into the origin and evolution of planetary systems has progressed remarkably with the help of celestial mechanics.
This field involves celestial mechanics mainly in terms ofthe dynamical motion of planetesimals (or protoplanets),
created by the condensation of dusts in protoplanetary nebula (e.g. Hayashi et al., 1985). Celestial mechanics becomes
more important in the latter stages of planet formation

when we can approximate the dynamics of particles (planetesimals or protoplanets) as that of point masses. One of
the biggest problems of the so-called “standard model” ofplanet formation was that the accretion timescales of jovian
planets, especially those of icy Uranus and Neptune could
be very long, even longer than the age of the solar systemitself. In the 1980s, efforts to solve this problem have led toinvestigation of the runaway growth of planetesimals (e.g.
Wetherill and Stewart, 1989).
As for the dynamical research in this field, we should
mention the enormous success of the special-purpose supercomputer for the gravitational N -body problem, GRAPE(see Section 5.4 for detail). Fully exploiting the ultra-highcomputational speed of the GRAPE system, S. Ida andE. Kokubo succeeded in performing a direct calculationof planetary accretion from planetesimals to protoplanets.
According to their series of papers (Kokubo and Ida, 1995,
1996, 1998, 2000, 2002), the mode of planetesimal accretionis inevitably “runaway” under the physical and dynamical
conditions of the protosolar nebula. Swarms of planetesimals are accumulated and grow oligarchically into much
fewer number of protoplanets (or planetary embryos). Now
the concept of oligarchic growth is also common to the formation of giant planets (e.g. Thommes et al., 2003).
There is an interesting dynamical feature of the final
stage of planet formation. When we look at the relationshipbetween the separation of protoplanets and their stability,
we can see a kind of an experimental scaling law. It is a relationship between the normalized distance between protoplanets and time until initial instability of the system. This
relationship was first confirmed by Chambers et al. (1996)
based on a previous research by Gladman (1993). Chambers et al. (1996) found that the time TE , until the firstclose encounter between two protoplanets, is exponentially
proportional to the protoplanetary distance ∆ normalized
by the mutual Hill radii; TE ∝ e∆ . They integrated systems of several protoplanets with equal-mass, initially onequally spaced coplanar and circular orbits. This relationship was later confirmed also by Yoshinaga et al. (1999)
when the orbits of the protoplanets are elliptic and inclined,
and by Ito and Tanikawa (1999, 2001) when Jupiter is outside the protoplanetary system as a perturber. In Ito andTanikawa’s (1999) numerical experiments, this exponential
law stands over five orders of magnitude in TE . However,the dynamical fundamentals of this relationship are still
unknown, though it is expected that they are based on acertain chaotic diffusion in phase space.
A similar but different scaling relationship is known
to exist between the mass and the instability timescale ofplanet-like bodies. Duncan and Lissauer (1997) integrated
systems with initial orbits that are identical to those of thesatellites of Uranus. In their experiments, the masses of thesatellites are all increased from actual values by the same
factor mf . They found that the orbital crossing time TC
obeyed a relationship as TC ∝ mα
f where α is a negative

Trends in 20th Century Celestial Mechanics
constant. These power law fits even extend over seven orders of magnitude in the crossing time TC . The theoreticalfundamental of the relationship TC ∝ mα
f is not known either. The two scaling relationships, TE ∝ e∆ and TC ∝ mαf ,may be two special cases of a more general relationship, andmay be involved with the relationship between Lyapunovtime TL and real instability time TI that we mentioned in
Section 2.2: TI ∝ TLγ in Eq. (13).
Let us observe the current solar system planets fromthe standpoint of normalized planetary separation, ∆. ∆between terrestrial planets is very large, such as ∆ > 26RH
where RH is the mutual Hill radius of planets. On the other
hand, those separations among jovian planets are muchsmaller, such as ∆ < 14RH . Although this difference could
be one of the essential factors in maintenance of the solar
system stable planetary motion over billions of years (cf.
Ito and Tanikawa, 1999, 2002), the reason of this differenceis yet to be precisely known. There have been some efforts
to explain the origin and difference of the planetary separations as being caused by dynamical evolution after theformation of the planetary system (e.g. Laskar, 1997), but
it would be more natural to ascribe these differences to theformation process of the planets. The results of N -bodynumerical simulations of the runaway growth of planetesimals tell us that the initial value of the normalized separation of protoplanetary systems could be ∆ = 5–10RH(Kokubo and Ida, 1998). So the current separations between the jovian planets are consistent with what N -body
integration predicts. The problem left is that of when andhow the normalized separation between terrestrial planets
became so large over the course of their growth from protoplanets. This question naturally leads us to the study ofthe late stage of terrestrial planet formation. A large number of papers have been, and will be, published along thisline of research, employing accurate and long-term N -body
integrations such as Chambers and Wetherill (1998), Levison et al. (1998), Agnor et al. (1999) Chambers (2001),
Kominami and Ida (2002, 2004), Kominami et al. (2005),
and Daisaka et al. (2006). Together with dynamical studiesof the diversity of extrasolar planetary systems, this field
will remain one of the major research areas of planetary
dynamics in the next decade.
As well as the late stage of planet formation where planetesimals and protoplanets play a central role, the very beginning stage of planet formation is also dominated by acertain kind of dynamics: the dynamics of dust. Dust is oneof the key ingredients of the solar system. As well as gas,
the dynamics of dust controls the early stage of planetary
formation toward the accretion of planetesimals. However,the dynamics of dust is somewhat different from that oflarger bodies, particularly in connection to the solar radiation. Dust particles are so small that they experience solar
radiation and stellar wind forces. These pressures change
the orbital elements of dust particles a great deal. Dust particles even escape from the system on hyperbolic orbits be-79

cause of the additional angular momentum and orbital energy from the radiation pressure. The Poynting–Robertson
effect and solar wind drag tend to circularize their orbits,
forcing the dust particles to slowly drift toward the central
star (Burns et al., 1979). As for a more detail of the dynamics of dust and its significance in the planet formation, consult publications such as Backman and Paresce (1993), Liouet al. (1995), Stern (1996), Yamamoto and Mukai (1998),
and Moro-Martı́n and Malhotra (2002, 2003).
4.4. Dynamics of small bodies
Dynamics of minor bodies in our solar system always
provides numerous interesting exercises for celestial mechanics. Now that we know of a huge number of small bodies
in our solar system, the number of dynamical phenomenathat requires celestial mechanics for explanation is soaring. Recent reviews such as Lazzaro et al. (2006) will be agood guide for readers to know how extensively and fast
the research of this field goes on. Among many kinds of minor body population research, particularly we focus on thefollowing topics in this subsection: asteroids and resonant
dynamics, the Yarkovsky effect and asteroid family, andtrans-Neptunian objects including the Kuiper Belt and theOort Cloud objects.
4.4.1. Asteroids and resonant dynamics
More than two hundred years have passed since the first
discovery of an asteroid (Ceres: discovered in 1801), andalmost ninety years have passed since the discovery of theHirayama asteroid families (Hirayama, 1918). More than
360,000 asteroids have been discovered as of April 2007, anda vast amount of knowledge has been accumulated about
asteroids; knowledge of spatial distribution, spectroscopic
information, taxonomy, shape, rotational motion, and orbital dynamics. In terms of long-term orbital dynamics, asteroids are particularly interesting from the standpoint of
their complicated resonant behavior in relation to the existence of numerous Kirkwood gaps in the main asteroid belt
(e.g. Murray and Dermott, 1999).
Mean motion and secular resonances
Two sorts of resonances that prevail in asteroid dynamics are mean motion resonances and secular resonances. A
mean motion resonance occurs when two bodies have orbital periods in a commensurability (a simple integer ratio,such as 1:2). In the main asteroid belt within 3.5 AU fromthe sun, the major mean motion resonances with Jupiterform gaps in the asteroid distribution such as at 3:1, 5:2,
7:3, or 2:1 resonances. They are the Kirkwood gaps (Kirkwood, 1867). Most asteroids have been ejected from thesezones by repeated encounters with Jupiter. On the other
hand in the outer main asteroid belt, 3:2, 4:3 and 1:1 resonances with Jupiter are populated by groups of asteroids:
the Hilda group (3:2), the Thule group (4:3), and the Tro-

80

Takashi Ito and Kiyotaka Tanikawa

jans (1:1). Dynamics of mean motion resonances is well illustrated and explained in Murray and Dermott (1999).
A secular resonance arises when the rate of change ofthe proper longitude of pericenter or of proper longitude ofascending node of a test body and one of the eigenfrequencies of the system of perturbing bodies are in a commensurability. The location of secular resonances in the solar
system (known under the names such as ν5 , ν6 , or ν16 ) israther complicated because of the coupling of terms related
to eccentricities and those related to inclinations in the disturbing function. The location corresponds to surfaces inthree-dimensional (a, e, I) space where a is the semimajoraxis, e is the eccentricity, and I is the inclination. Hencesecular resonances form the boundaries of the main asteroid belt in the (a, e, I) space. For more detail about thetheories of secular resonances, consult publications suchas Williams (1969, 1979), Williams and Faulkner (1981),
Scholl et al. (1989), Morbidelli and Henrard (1991a,b), Milani and Knez̆ević (1990), Froeschlé and Morbidelli (1994),
Knez̆ević and Milani (1994), Knez̆ević et al. (1995), andMurray and Dermott (1999).
In the early 1980s, J. Wisdom gave asteroid dynamists
a huge shock. He introduced a mapping method to simulate asteroid dynamics, exploiting Dirac’s δ-function (Wisdom, 1982, 1983). In his map, the continuous action of force
is replaced by impulses, resulting in a drastic shorteningof the computation time. Wisdom observed that asteroid
orbits in the 3:1 mean motion resonance attain a high eccentricity and may collide with Mars or the Earth. Thismethod and result were initially regarded with some doubt,
but soon confirmed and recognized as correct. The speedof his mapping method revolutionized the study of asteroid dynamics. Before Wisdom, a one-million-year integration of an asteroid orbit was a tremendous task. With themap derived from the planar elliptic restricted three-bodyproblem, Wisdom found that the orbital eccentricity of thesmall body can suddenly increase after a long time of quiescence. In addition, it turned out that the increase and decrease in eccentricity is repeated chaotically. This was thefirst work to explicitly deal with solar system chaos. Wisdom’s result was later confirmed by himself through direct
integrations of the equations of motion (Wisdom, 1987b).
After Wisdom, a large number of celestial mechanists
have struggled through an enormous number of publications to explain the origin of the Kirkwood gaps at
various resonances. For example, Yoshikawa (1987, 1989),
within the framework of the planar (and partially threedimensional) elliptic restricted three-body problem, semianalytically considered the motion of asteroids in 3:1, 5:2,
7:3, and 2:1 resonances. He also carried out numerical integrations and illustrated the orbital changes in mean motion
resonances. Morbidelli and Moons (1993) and Moons andMorbidelli (1995), in the restricted four-body problem of
the Sun, Jupiter, Saturn, and an asteroid, studied secularresonances in 3:2 and 2:1 mean motion resonances three-dimensionally and 4:1, 3:1, 5:2, and 7:3 mean motion resonances. They mapped effective places of secular resonancesin the region of the mean motion resonance. As an example, in the 3:2 resonance, ν5 and ν6 occupy large chaotic
areas, whereas ν16 resonance works in the 2:1 resonance.
The presence of Saturn introduces a number of additional frequencies into Jupiter’s orbit, which can significantly influence the orbits of asteroids. The resonant arguments of these three-body resonances contain the longitudes of Jupiter and asteroid together with either the secular frequency g6 , or the longitude of Saturn. Resonances
involving the longitude of Saturn are analogs of the Laplaceresonance in the Jovian satellite system (e.g. Greenberg,
1977). Murray et al. (1998) considered the effect of gravitational perturbations from Jupiter on the dynamics of asteroids, when Jupiter is itself perturbed by Saturn. Murray
et al. (1998) showed that many three-body resonances involving the longitude of Saturn are chaotic, and they gave
simple expressions for the width of the chaotic region andthe associated Lyapunov time.
Kozai mechanism
In addition to ordinary secular and mean motion resonances, Kozai (1962b) investigated the long-term motionof asteroids with high inclination and eccentricity. As we
mentioned in Section 3.5, Kozai found what we now callthe Kozai mechanism (or the Kozai behavior) that drives
the eccentricity and the inclination of asteroids very high,such as e ∼ 1 or I ∼ 90◦ under certain conditions. Kozai(1962b) is one of the pioneers in exploiting equi-energy
curves (equi-Hamiltonian maps) with one degree of freedom for the purpose of investigating the global dynamical
behavior of test particles undergoing the perturbation ofmassive bodies. The concept of the Kozai mechanism was
extended by Michel and Thomas (1996) to include the four
jovian planets. Michel and Thomas (1996) showed that the
Kozai mechanism provides a temporary protection mechanism for near-Earth asteroids from close planetary encounters. A near-Earth asteroid can be locked in the Kozai
mechanism in a stable state, where it can stay for a comparatively long time before a close approach with a planet
drastically changes its orbit. Now the concept of the Kozai
mechanism is ubiquitous throughout celestial mechanics,
not only in the asteroidal motion field, but in the studyof the motion of Kuiper Belt objects (e.g. Duncan et al.,
1995; Kuchner et al., 2002), comets (e.g. Bailey et al., 1992;
Bailey and Emelýanenko, 1996), Pluto (e.g. Williams andBenson, 1971; Wan et al., 2001), natural and artificial satellites (e.g. Giacaglia et al., 1970; Mignard, 1975; Nesvorný
et al., 2003), meteoroids (e.g. Wetherill, 1974; Valsecchi
et al., 1999), binaries and triple stars (e.g. Innanen et al.,
1997; Mikkola and Tanikawa, 1998a; Mardling and Aarseth,2001), forming planets (e.g. Chambers and Wetherill, 1998;
Levison et al., 1998; Levison and Agnor, 2003), extrasolar
planets (e.g. Laughlin and Adams, 1999; Ford et al., 2000;

Trends in 20th Century Celestial Mechanics
Drovak et al., 2003), and even supermassive binary black
holes (e.g. Blaes et al., 2002). Gronchi and Milani (1999b)
devised a semi-analytical method to extend the applicationof the Kozai mechanism to planet-crossing orbits.
Secondary resonances
The dynamics of asteroids is characterized not only bymean motion resonances with Jupiter or simple secular resonances, but by complicated secondary resonances and resonance overlaps. For example, when we consider the planar circular restricted three-body problem, the resonancecan be classified into primary and secondary resonances.
Primary resonances are due to commensurabilities between
the mean motions of the perturbing body and the test particle. Secondary resonances occur when the libration frequency of a primary resonance is commensurate with oneof the higher-order primary resonances (Murray and Dermott, 1999). Secondary resonances are not just a productof theory, but an actual phenomenon that can have a significant effect on the asteroid and satellite dynamics. We
can see a typical example of secondary resonances in theorbit of a Uranian satellite, Miranda. The series of works
by R. Malhotra, S. F. Dermott, and C. D. Murray (e.g.
Dermott et al., 1988; Malhotra and Dermott, 1990; Malhotra, 1990) on the role of secondary resonances in theorbital history of Miranda established an important basis
for the dynamical treatment of secondary resonances in solar system dynamics, and is concisely summarized in Kortenkamp et al. (2004). Malhotra and her colleagues showedthat the anomalously large orbital inclination of Miranda
is naturally explained as a consequence of the secondaryresonance sweeping due to the tidal evolution within a primary 3:1 inclination-type mean motion resonance with Umbriel, another satellite of Uranus. The primary resonance increased Miranda’s inclination from an initially small value,
but the resonance was temporary. As Miranda’s inclination
approached its current value, the satellites were captured
in a secondary resonance which amplified their primaryresonance libration amplitude. This eventually caused thesatellites to escape from the 3:1 mean-motion resonance,leaving Miranda with an inclination that is preserved to thepresent. The secondary resonance implicated in this casewas due to a 1:3 commensurability between the libration
frequency of the primary mean motion resonance angle andthe secular frequency of precession of the relative lines ofnodes of the two satellites.
Another famous example of secondary resonances in thesolar system is seen in the motion of Pluto, particularly
in its resonant relation with Neptune. Pluto’s orbital motion is characterized by three major resonances with Neptune (Kinoshita and Nakai, 1995, 1996) which ensure that
Pluto and Neptune do not encounter each other: (i) Pluto
and Neptune are in a 2:3 mean motion resonance with thecritical argument θ1 = 3λP − 2λN − $P . (ii) Pluto’s argument of perihelion ωP = θ2 = $P − ΩP librates around81

90◦ with a period of about 3.8 × 106 years. This is a typical example of the Kozai mechanism. (iii) The longitudeof Pluto’s node referred to the longitude of Neptune’s node(θ3 = ΩP − ΩN ) circulates, and the period of θ3 circulation
is equal to the period of θ2 libration. The third resonance,sometimes called the “superresonance” (Milani et al., 1989;
Wan et al., 2001), is a 1:1 secondary resonance between θ2and θ3 .
Let us give some examples of the study of secondary
resonances in the asteroidal motion. Moons and Morbidelli(1993) studied the motion of asteroids in the mean motioncommensurabilities through the planar restricted threebody problem. They provided global pictures of the dynamics in the region of secondary resonances. Ferraz-Mello(1994) studied the structure of the phase space of the 2:1
resonance through the planar averaged asteroidal threebody problem. The only chaotic regions he found were those
associated with secondary resonances, confined to low eccentricities, and those associated with large libration amplitudes. Nesvorný and Ferraz-Mello (1997) applied Laskar’s(1993) frequency map analysis to the dynamical models of
the 2:1 asteroidal mean motion resonance; the planar restricted three-body model and the planar restricted fourbody model with Saturn as the third primary. Their result reproduced the chaotic region formed by the overlapof secondary resonances in low eccentricities. Nesvorný andFerraz-Mello (1997) examined the chaos generated by highorder secondary resonances in moderate eccentricities.
Secondary resonances might have worked as a clearing
mechanism of trojan-type companions of Neptune during
primordial migration of four giant planets (Kortenkamp
et al., 2004). The loss of Neptune trojans (which we know
very few of now) can occur when trojan particles are
swept by secondary resonances associated with mean motion commensurabilities of Uranus with Neptune (2:1).
These secondary resonances arise when the circulation frequencies of critical arguments for Uranus–Neptune mean
motion near-resonances are commensurate with harmonics of the libration frequency of the critical argument for
the Neptune–trojan 1:1 mean motion resonance. Trojans
that are trapped in the secondary resonances typicallyhave their libration amplitudes amplified until they escapethe 1:1 resonance with Neptune. Trojans with large libration amplitudes are susceptible to loss during sweeping by
numerous high-order secondary resonances.
Resonance overlap
In general, a resonance can have its own finite width.
This finite width can cause an overlap of nearby resonances,which frequently results in the onset of chaos. As B. V.
Chirikov predicted nearly 30 years ago (Chirikov, 1979),
this resonance overlap plays a significant role in the dynamical behavior of asteroids. Right after Chirikov (1979), Wisdom (1980) featured the resonance overlap in relation to theonset of stochastic behavior of the planar circular restricted

82

Takashi Ito and Kiyotaka Tanikawa

three-body problem. After that, the concept of resonanceoverlap was gradually recognized and accepted by the society of solar system dynamics. For example, Froeschlé and
Scholl (1989) performed numerical experiments that indicated possible chaotic motion due to overlapping resonances in the asteroid belt. A secular resonance may overlap
with another secular resonance or with a mean motion resonance. Morbidelli and Moons (1995) carried out numericalintegrations of objects that are in the 3:1 mean motion resonance. Their result indicates that the dynamics in the 3:1
resonance zone is strongly chaotic, due to the overlappingof secular resonances. As a consequence, the eccentricity ofbodies in this zone can increase approximately to unity on atime scale of one million years. Holman and Murray (1996)
investigated the origin of chaos near high-order mean motion resonances in the outer asteroid belt. They surveyed
the variation of the Lyapunov time with semimajor axis
throughout the outer belt region with the elliptic restricted
three-body problem. They also developed an analytic theory based on resonance overlap that predicts the widths insemimajor axes as well as typical Lyapunov times of thechaotic zones associated with these resonances. Michel andFroeschlé (1997) presented the location of linear secularresonances in the region of semimajor axes a less than 2
AU where many near-Earth asteroids are found. The positions of the secular resonances in the plane (a, I) show thatthe inner solar system is dynamically very complex: a lotof resonances are present, and some of them overlap. An
extensive review has been done by Lecar et al. (2001) for
the relationship between asteroidal resonances and chaos,including resonance overlaps.
Trojans and migration of resonances
Talking about resonances, minor bodies around the 1:1
mean motion resonance with major planets, Trojan-typecompanions or co-orbital objects, still provide us with a big
challenge in terms of celestial mechanics, especially as a realistic application of the restricted three-body problem. As
of April 2007, more than 2000 asteroids are known to share
Jupiter’s orbit, five others are known as Neptune trojans,
but none is known around Saturn and Uranus. The originof this difference is not yet fully understood, though several hypotheses have been proposed. Consult the introduction of Kortenkamp et al. (2004) and the references therein.
Also, recent discoveries of co-orbital objects around theEarth (Wiegert et al., 1997, 1998, 2000; Connors et al.,2002) and Venus (Mikkola et al., 2004; Brasser et al., 2004;
Morais and Morbidelli, 2006) have prompted a number ofquestions about their origin and stability.
Among many theories and hypotheses as to the origin
of Jupiter’s trojans, a series of papers by A. Morbidelli,
R. Gomes, K. Tsiganis, and H. F. Levison are of particular interest. In previous literature on asteroid dynamics,
Jupiter’s trojans have most frequently been considered to
be planetesimals that formed near Jupiter and were cap-

tured onto their current orbits while Jupiter was growing (e.g. Marzari and Scholl, 1998; Fleming and Hamilton, 2000), possibly with the help of gas drag and/or collisions (e.g. Shoemaker et al., 1989). However, this hypothesis cannot explain some basic properties of the trojan population, in particular its broad orbital inclination distribution, which ranges up to ∼ 40◦ (Marzari et al., 2002). Morbidelli et al. (2005) showed that Jupiter’s trojans could haveformed in more distant regions and been subsequently captured into co-orbital motion with Jupiter when the giant
planets migrated by removing neighboring planetesimals.
The capture was possible during a short period of time, justafter Jupiter and Saturn crossed their mutual mean motion
2:1 resonance, when the dynamics of the trojan region werecompletely chaotic.
Morbidelli et al.’s (2005) simulations of the process involving the resonance passing of Jupiter and Saturn do not
only satisfactorily reproduce the orbital distribution of thetrojans and their total mass, but it can reasonably accountfor an intense period of planetesimal bombardment in the
inner solar system, collectively called the “late heavy bombardment” or the LHB. The petrology record on the Moonsuggests that a cataclysmic spike in the cratering rate occurred about 700 million years after the planets formed(Tera et al., 1974; Ryder, 1990; Hartmann et al., 2000).
Planetary formation theories cannot naturally account fora bombardment period so late in the solar system history(Morbidelli et al., 2001; Levison et al., 2001a). Using thesame dynamical model as Morbidelli et al. (2005), Gomes
et al. (2005) proposed that the LHB was triggered by the
rapid migration of the giant planets, which occurred after a long quiescent period. During this burst of migration, the planetesimal disk outside the orbits of the planets was destabilized, causing a sudden massive delivery ofplanetesimals to the inner solar system. The asteroid belt
was also strongly perturbed, with these objects supplying
a significant fraction of the LHB impactors in accordancewith recent observational evidence (e.g. Bogard, 1995; Cohen et al., 2000; Kring and Cohen, 2002; Strom et al., 2005;
Ito and Malhotra, 2006).
Incidentally, the dynamical model used by Morbidelli
et al. (2005) and Gomes et al. (2005) also explains well
the current values of eccentricities of major giant planets;
Jupiter, Saturn, and Uranus. Again using the same dynamical model, Tsiganis et al. (2005) showed that a planetarysystem with initial quasi-circular, coplanar orbits would
have evolved to the current solar system planetary orbital
configuration, provided that Jupiter and Saturn crossed
their 2:1 mean motion resonance. Tsiganis et al. (2005)
showed that this resonance crossing could have occurred as
the giant planets migrated owing to their interaction witha disk of planetesimals (Fernández and Ip, 1984; Malhotra, 1995). Tsiganis et al.’s dynamical model reproduces
all the important characteristics of the giant planets’ orbits; their final semimajor axes, eccentricities and mutualTrends in 20th Century Celestial Mechanics
inclinations. In addition, it turned out that this migration
scenario suggests that the present obliquities of the giant
planets were presumably achieved when Jupiter and Saturn crossed the 2:1 mean motion resonance (Brunini, 2006).
The existence of the regular satellites of the giant planets
does not represent a problem in this scenario because, although they formed soon after the planet formation, they
can follow the slow evolution of the equatorial plane of planets during the migration stage.
4.4.2. The Yarkovsky effect and asteroid families
A group of asteroids that have similar proper orbital
elements are called an asteroid family. The concept of theasteroid families, i.e. the Hirayama asteroid families, wasfirst advocated by K. Hirayama (1918), followed by numerous research efforts. Hirayama found three families in his
celebrated 1918 paper: Koronis, Eos, and Themis families.
The Hirayama asteroid families are considered remnantsof collisional disruption of Myr to a few Gyr ago (Zappalà
et al., 1994). The proper elements are analytically definedas constants of motion of a suitably simplified dynamicalsystem, obtained by classical transformation theory froma real dynamical system (Milani and Knez̆ević, 1990; Morbidelli, 2002; Knez̆ević et al., 2002). If the real dynamics
is regular, proper elements show only very small oscillations around a mean value that stays constant with time
(Knez̆ević et al., 1995). However, if the motion is chaotic,
the mean value of the proper elements can drift in time at
a rate that depends on the properties of the chaotic motion(Milani and Farinella, 1994). Now it has turned out that theproper elements could be altered, and asteroid families are
diffused a great deal over the long timescale of 100 Myr toGyr, owing to numerous weak resonances (Morbidelli and
Nesvorný, 1999), to mutual close encounters with massive
asteroids (Nesvorný et al., 2002), as well as to the Yarkovsky
thermal effect (Bottke et al., 2001; Nesvorný et al., 2002).
The Yarkovsky effect is caused by a thermal radiation
force that makes objects undergo slow but steady semimajor axis drift and spin-up/down as functions of theirrotation, orbit, and material properties (Rubincam, 1995,
1998). Numerical results suggest that this mechanism canbe used to: (i) deliver asteroids with diameter D . 20 km
from their original locations in the main belt to resonanceregions that are capable of transporting them to Earthcrossing orbits, (ii) disperse asteroid families, with drifting bodies jumping or becoming trapped in mean-motion
and secular resonances within the main belt, and (iii) modify the rotation rates of asteroids a few km in diameter orsmaller.
The Yarkovsky effect solves several definite inconsistencies that gravitational forces cannot explain in terms of asteroid dynamics. For example, the meteorite cosmic-ray exposure age of asteroids is an order of magnitude longer thanconventional model predictions. The mechanism for continuously supplying near-Earth asteroids to the mean motion

83

resonance zones in the main belt is not well explained, either (Bottke et al., 2002). Since the phenomena in the solar
system are not always “clean” (i.e. expressed only by thelaw of gravity), sometimes we have to extend the applicablearea of celestial mechanics from what ordinary gravity law
governs to something that simple gravity law cannot cope
with, especially in asteroid dynamics. The Yarkovsky effectis a typical and very important one of these “dirty” forces.
The Yarkovsky effect should now be considered as being as important as collisions and gravitational perturbations are to our overall understanding of asteroid dynamics,
especially that of the Hirayama asteroid families. The orbital distributions of prominent families are thought to be
by-products of catastrophic disruption events. But ejection
velocities derived from the orbital elements of families aremuch higher compared with the results of impact experiments and simulations. One way to resolve this contradiction is to assume that D . 20 km family members have undergone semimajor axis drift due to the Yarkovsky effect.
W. Bottke demonstrated how the Yarkovsky effect could
influence creation of some of the existing Hirayama asteroid families (Bottke et al., 2001). Their result explains why
families are sharply bounded by nearby Kirkwood gaps,
why some families have asymmetric shapes, and the curious presence of family members on short-lived orbits.
The Yarkovsky effect also influences the rotationalmotion of asteroids. The Yarkovsky–O’Keefe–Radzievskii–
Paddack (YORP) effect, a version of the Yarkovsky effect,can spin up or down asteroids with 10 km diameter on a
108 -year timescale (Rubincam, 2000). Smaller asteroids
spin up or down even faster due to the radius-squared dependence of the YORP timescale. The YORP effect may
explain the rapid rotation, slow tumbling, and spin-orbit
coupling appearing in some asteroid rotational motions.
Recently, the Yarkovsky and the YORP effects haveeven been directly detected by accurate observations.
Vokroulický et al. (2000) argued that a precise radar refinement of the orbits of near-Earth asteroids offers the possibility of detecting the Yarkovsky effect during the next few
decades. Chesley et al. (2003) reported the prospects forthe detection of the Yarkovsky effect from fine-precision
Doppler radar observations of the half-kilometer asteroid,(6489) Golevka. Lowry et al. (2007) presented precise optical photometric observations of a small near-Earth asteroid, (54509) 2000 PH5, acquired over 4 years. They foundthat the asteroid has been continuously increasing its rotation rate ω over this period by dω/dt = 2.0(±0.2) × 10−4
degrees per day squared. They simulated the asteroid’s
close Earth approaches from 2001 to 2005, showing thatgravitational torques cannot explain the observed spin
rate increase. Dynamical simulations suggest that 2000
PH5 may reach a rotation period of ∼ 20 seconds toward
the end of its expected lifetime. Also, Taylor et al. (2007)
showed that radar and optical observations revealed thatthe continuous increase in the spin rate of this asteroid

84

Takashi Ito and Kiyotaka Tanikawa

can be attributed to the YORP effect. The change in spin
rate is in reasonable agreement with theoretical predictions for the YORP acceleration of a body with the radardetermined size, shape, and spin state of 2000 PH5. The
detection of asteroid spin-up supports the YORP effectas an explanation for the anomalous distribution of spin
rates for asteroids under 10 kilometers in diameter and as
a binary formation mechanism. More accurate observation
results will reveal the true nature of the Yarkovsky effect,
which also gives us significant information about the dynamical history of asteroids as well as of other small bodies
in the solar system.
4.4.3. Kuiper Belt objects
Owing to the recent development of high-accuracy observations, physical and dynamical research into the outermost area of our solar system seems to be approaching its
zenith. Small bodies in the outermost area of the solar system that orbit the sun at distances farther than Neptune arecollectively called trans-Neptunian objects (TNOs). The
Kuiper Belt, the scattered disk (which is occasionally included in the Kuiper Belt), and the Oort Cloud are thenames for three divisions of this volume of space.
Among the trans-Neptunian orbital zones, the Kuiper
Belt is the innermost region: an area extending from theorbit of Neptune to r ∼ 50 AU (r denotes heliocentric distance). Although the existence of the objects in this region
had been predicted since the 1940s to the 1950s (Edgeworth, 1949; Kuiper, 1951), it took more than 40 years forthe first object to be discovered (Jewitt and Luu, 1993).
Since then, more than 1000 Kuiper Belt objects (KBOs)
have been found, as of April 2007. The recent status of the
KBO research is summarized in Weissman (1995), Jewitt(1999), Luu and Jewitt (2002), and Gladman (2005).
KBOs have been dynamically categorized into three
types (Trujillo, 2003). About half of the KBOs have nearlycircular orbits (e < 0.2) and form a kind of ring structureat around 42–48 AU. These are called “classical” objects.
About 10% of the KBOs are in 2:3 mean motion resonancewith Neptune, as Pluto is. These are called “plutinos”. The
rest of the KBOs, about 40% of them, have highly eccentric
orbits, and are therefore called “scattered” objects. The orbits of the scattered KBOs are distributed very widely fromthe inner part of the classical object zone to the furthest
area, at around 900 AU. Since the orbital characteristics ofthe scattered objects are so different from that of classicalobjects and plutinos, the space volume that scattered objects occupy is recently dubbed “the scattered disk”, andis regarded as a separate category among the TNOs.
In spite of the numerous research efforts devoted to
KBOs, detailed dynamical structure and origin of the
KBOs are not well known yet. First of all, there seems
to be a relatively sharp outer boundary for the KBOs,
especially the classical Kuiper Belt (cf. Allen et al., 2002).
Beyond a = 48AU, no KBO has a circular orbit. A few

explanations have been proposed to account for the formation of this edge, but none of them is widely accepted
yet. One theory claims that the outermost particles werepretty much dynamically excited and stripped away by aclose encounter with a passing star in the early stages ofsolar system history (e.g. Ida et al., 2000; Kobayashi and
Ida, 2001). Another says that an Earth-sized planet(s)may have passed through the primordial Kuiper Belt andpumped up the eccentricities and inclinations of the particles there (e.g. Petit et al., 1999). Yet another possibility
is that the primordial solar nebula and its dispersal caused
the migration of some strong secular resonances outside of
Neptune, which enhanced the random velocity of a certain
part of the KBOs and created the edge (e.g. Nagasawa andIda, 2000).
In addition to the existence of the outmost edge of theclassical KBOs, recent observations (e.g. Burns, 2002; Veillet et al., 2002; Noll et al., 2002; Stephens and Noll, 2006)
have shown that the fraction of binaries among the KBOs is
unexpectedly high, much higher than the fraction of asteroid binaries in the asteroid belt. Also, the orbital separationof KBO binaries is usually larger than that of asteroid binaries, and their mass ratio is closer to unity (Veillet et al.,2002). To theoretically understand these characteristics of
KBO binaries, Funato et al. (2004) showed that the fraction of the KBO binary population was initially as small
as that of the current asteroid binaries, but was gradually
increased through three-body dynamical interactions. Thiswas possible because the tidal force from the Sun is weakerin the Kuiper Belt area than in the asteroid belt. Different explanations have been attempted by Weidenschilling(2002) and Goldreich et al. (2002).
Some of the important dynamical characters of the
KBOs probably have their origins in their formation stage.
One of the most intriguing ideas that are derived fromthe dynamical characteristics of the KBOs is that thecurrent orbital distribution of the KBOs was producedby the so-called resonance capture mechanism, triggeredby the orbital migration of Neptune (and probably otherthree giant planets). This idea was first explored extensively by Malhotra (1993) to explain the origin of Pluto’s
orbit. Pluto’s orbit is quite unusual compared with otherplanetary orbits, being the most eccentric and the mostinclined. The orbits of Pluto and Neptune overlap. But
close approaches of these two planets are prevented by
the existence of a resonance condition such that Pluto’s
orbital period is exactly 3:2 that of Neptune. This ensuresthat the conjunctions always occur near Pluto’s aphelion.
Malhotra showed that Pluto could have acquired its current orbit during the late stages of planetary accretion andthe formation of the Oort Cloud comets, when the jovian
planets underwent significant orbital migration as a resultof encounters with residual planetesimals. As Neptune
moved outwards, a small body like Pluto in an initiallycircular orbit could have been captured into the 3:2 mean

Trends in 20th Century Celestial Mechanics
motion resonance. After this event, its orbital eccentricity
would rise rapidly to its current Neptune-crossing value. If
this mechanism occurred, the entire region from the orbit
of Neptune to r ∼ 50 AU would have been swept by strong
mean motion resonances. This means that the resonancecapture could have occurred not only for Pluto but forother KBOs. Malhotra’s (1995) classic paper dealing withthis theory explains some of the major consequences in the
Kuiper Belt, such as the fact that most objects in the region beyond Neptune and up to r ∼ 50 AU are confined in
the narrow resonance zone, especially in 3:2 and 2:1 mean
motion resonances.
A vast amount of literature has followed Malhotra’s pioneering work of planetary orbital migration and resonant
capture, not only in the area of KBOs formation research,
but covering almost all stages of planetary formation. Oneof the recent reports has come from Levison and Morbidelli(2003). They have shown that the objects currently observed in the dynamically cold Kuiper Belt were most probably formed within ∼ 35 AU, and were subsequently pushed
outward by Neptune’s 2:1 mean motion resonance during
its final phase of migration. They have concluded that theentire Kuiper Belt formed closer to the Sun and was transported outward during the final stages of planet formation.
All these questions, hypotheses, and enigmas on the dynamical characteristics of the KBOs as well as of the scattered objects should be, and will be confirmed or eliminatedby more accurate observations of distant regions such as
r > 50 AU (Brown et al., 2004). Accordingly, discoveriesof very large KBOs or scattered objects such as (90377)
Sedna or (136199) Eris will expand the outmost region theouter solar system to a great deal (e.g. Brown et al., 2005).
Actually, the discovery of Eris led astronomers and planetary scientists to question the definition of the term planet,as Eris turned out be larger than Pluto. This epochmaking finding eventually ended up with the demotion of Pluto
from a planet to an object among the new category of small
body population in the solar system: dwarf planet. The dynamical characteristics and origin of the outermost dwarf
planets, such as Eris or Sedna, are now the subjects of intense research, and will be one of the central topics of thesolar system dynamics in the new century.
Speaking of KBOs, though it seems digressive to talkabout astronomical engineering in this manuscript, the ideaof gradually expanding the Earth’s orbit artificially using a
KBO over a very long timespan (Korycansky et al., 2001) ismore than just an SF story, but it is an interesting thought
experiment from the standpoint of the KBO and planetary
dynamics.
4.4.4. The Oort Cloud
The existence of a small body population at the very
outer edge of the solar system, now known as the Oort
Cloud, has also been imagined for a long time since the1950s. In 1950, J. H. Oort advocated the idea that our so-

85

lar system is surrounded by a distant cloud of comet stuff(Oort, 1950). Although the existence of the Oort Cloud was
proposed almost at the same time as that of the Kuiper
Belt (Edgeworth, 1949; Kuiper, 1951), there is a significant
difference between them: many KBOs have already been
discovered and confirmed, while none of the Oort Cloudobjects has yet been directly detected. This is mainly because of their very large distance from the Earth and theSun. However, we can anticipate their existence by the fluxand distribution of long-period comets coming all the way
to the inner solar system. Comet clouds such as the Oort
cloud are expected to exist not only in our solar system butaround other stars. Consult Dones et al. (2004) for a comprehensive review of the cometary dynamics and the Oort
Cloud formation.
The pioneering work in comet cloud dynamics has been
carried out by M.J. Duncan, T. Quinn, and S. Tremaine(e.g. Duncan et al., 1987), based on a dynamical modelof galactic tidal fields (e.g. Heisler and Tremaine, 1986;
Heisler et al., 1987). Duncan et al. (1987) simulated theformation of the comet cloud and its subsequent evolution
over the age of the solar system. Their simulations showedthat the formation of the current comet cloud was driven
mainly by an interaction between planetary perturbations
and torquing due to galactic tides. The inner edge of thecloud is estimated at about 3000 AU, the radius where thetimescales for the two processes are comparable. Duncan
et al. (1987) also suggested that the flux of comets into theinner solar system initiated by the close passage of a star
may be up to 20 times higher than the steady state rate. See
also Duncan et al. (1988), Quinn et al. (1990), or Hogg et al.
(1991) for other publications about comet cloud dynamics
by S. Tremaine and his colleagues.
The dynamical source of long-period comets is often
assumed to be around the Uranus–Neptune zone, since
Jupiter and Saturn are considered to have ejected most of
their icy planetesimals to interstellar space. The expected
fraction of planetesimals ejected from each giant planet
zone to the Oort Cloud shows that the relative efficiency of
Uranus and Neptune in placing objects in the Oort Cloud
far exceeded that of Jupiter and Saturn, by factors up to 20
or more (cf. Fernández and Ip, 1981; Duncan et al., 1987;
Fernández, 1997). On the other hand, Higuchi et al. (2006)
investigated the first dynamical stage of comet cloud formation, the scattering of planetesimals by a planet. Applying their results to the solar system Oort Cloud, they foundthat Jupiter is the planet most responsible for producing
candidate elements of the Oort Cloud, as long as the inclination of planetesimals is constant or proportional to the
reduced Hill radius of each planet.
Another group that has been working on dynamical origin of the Oort Cloud is L. Dones, H. Levison, P. R. Weissman and their colleagues. In their series of numericalintegrations (e.g. Dones et al., 1998; Levison et al., 1999;
Dones et al., 2000) they calculated the orbits of thousands

86

Takashi Ito and Kiyotaka Tanikawa

of test particles from initially low-inclination and loweccentricity orbits with semi-major axes of 4–40 AU over
4 billion years. They assumed that the long-period cometsreached the Oort Cloud through planetary perturbations
and galactic tides, and tried to determine the birthplace of
the Oort Cloud comets. Levison et al. (2001b) studied theorigin of Halley-type comets using the orbital integrationof many test particles initially entering the planetary system from the Oort Cloud. They found that an isotropically
distributed Oort Cloud does not reproduce the observed
orbital element distribution of the Halley-type comets.
Levison et al. (2002) calculated the number of dormant,
nearly isotropic Oort Cloud comets using various orbital
distribution models. According to them, we should havediscovered ∼ 100 times more dormant, nearly isotropic
comets than are actually seen under the assumption thatcomets are never destroyed. This result indicates that the
majority of comets must physically disrupt as they evolve
inward from the Oort Cloud. The Oort Cloud is suspected,
incidentally, to be seamlessly connected to the inner partof the TNO region: the scattered disk. See Levison andDuncan (1997), Duncan and Levison (1997), Levison et al.
(2001b), or Levison et al. (2006) for the relation between
the Oort Cloud dynamics and the scattered disk.
The question on the dynamical source of the Oort Cloudobjects (or long-period comets) will be open until we haveclearly known their accurate spatial distribution. Overall,
the dynamical study of the Oort Cloud is not yet quite
strongly constrained by observations; we have not seen anyof the Oort Cloud object at their home region of several
thousand AU away. Although the direct detection of Oort
Cloud objects in their home region of 3000 AU or further
is not easy, it is not impossible either, using star-object occultation event sensing (cf. Kaplan et al., 2003) or spacecraft missions. Once the Oort Cloud objects are directly
detected, a significant amount of information on the dynamical characteristics of the small body population in
the outermost area of the solar system will be delivered
to us. Also, accurate observational data obtained from theground-based and satellite astrometry provide us with more
and more definitive information about the structure of ourgalaxy and the possible history of stellar passages near theedge of our solar system (e.g. Garcı́a-Sánchez et al., 1999,
2001). This will tell us more about dynamical history andevolution process of the Oort Cloud.
4.5. Planetary rotation
Dynamics of planetary rotation, especially those of the
Earth, has been examined in a wide range of research fields,
according to the timescales considered. The timescales
range from an hour or a day for the nutation with small
amplitude, to billions of years for the secular change of
the Earth’s spin state due to the tidal evolution of the
Earth–Moon system.

4.5.1. Nutation and precession of the Earth
Research on the Earth’s rotational motion with shorttimescale (days to months) is renowned for its high accuracy. A hundred years have passed since H. Kimura discovered the so-called z-term in the Earth’s latitude variations (Kimura, 1902). The z-term was finally explained
as an error of the coefficient of the semi-annual nutation
term (Wako, 1970). There is a long history for this solution. In the ILS (International Latitude Service) network,
before 1955, two groups of stars (six pairs of star in eachgroup in two hours) had been observed every night at eachof the latitude stations. The observational instrument was
the VZT (Visual Zenith Telescope). One group of stars was
replaced by another group in the next month. One year
later, the first group of stars came back, and the cycle wasclosed. This chain method was to equalize positional errors of the stars in the catalogue. Kimura somehow noticed
that the two-group observations are not enough to separate
annual meteorological variation and semi-annual nutation.
In the 1938 IAU (International Astronomical Union) general assembly, Kimura proposed, via Y. Hagihara, an observation program consisting of three groups of stars pernight where one group of stars is replaced by another every month. Though this proposal was not adopted by the
IAU in 1938, the ILS started the three-group observationsin 1955. Then, soon the z-term turned out to represent thedifference of the rigid Earth model and the non-rigid Earthmodel. Wako’s result gave substantial support to the theory of the Earth rotation with a liquid core and an inner
solid core developed by Jeffreys and Vicente (1957) andothers. From the age of Kimura until the 1970s, this kindof research on the Earth rotation was done using latitude
and polar motion observations by way of positional astronomy. Recently, these classical methods have given way to
new space technology with very high accuracy, such as SLR(Satellite Laser Ranging), VLBI (Very Long Baseline Interferometry), and GPS (Global Positioning System). Now
the accuracy of these observations is extremely high, requiring of the corresponding theories the highest accuracy
possible.
As for the nutation of the rigid Earth, which has theshortest timescales among these sorts of phenomena, H. Kinoshita has made a significant contribution. The IAU (International Astronomical Union) had officially adopted thetheory of E. M. Woolard (Woolard, 1953) as a standard
model of the Earth nutation. However, Woolard’s theory
was no longer sufficiently accurate to compare with observations having very high accuracy using VLBI and othermethods in the 1970s. Woolard’s theory also had the defectof calculating the nutation of the instantaneous rotational
axis of the Earth, not of the axis of the Earth’s shape. Kinoshita (1977a) had brought a Hamiltonian dynamics into
the field of the rigid body nutation of the Earth. He recalculated the nutation in terms of the axis of shape, the axis
of instantaneous rotation, and the axis of angular momen-

Trends in 20th Century Celestial Mechanics
tum of the rigid Earth. The number of periodic components
calculated in Kinoshita (1977a) is 106, and the order of error is better than 10−4 arc second. His theory was extendedin Kinoshita and Souchay (1990) as well as Souchay andKinoshita (1991).
At present, the difference between theoretical and observational values of the Earth nutation is recognized as
being caused by the non-rigid components of the Earth’s
interior. T. Sasao had made a large and significant contribution to this field with a very accurate nutation theory encompassing the fluid dynamics of the Earth’s core
(Sasao et al., 1980); the theory is succeeded by Shirai and
Fukushima (2000, 2001). J. Getino is trying to incorporate
the non-rigidity of the Earth rotation in Hamiltonian formalism (Getino et al., 2001).
4.5.2. The Earth–Moon system
When the timescale of the variation in the Earth rotation becomes as long as 108 to 109 years, the dynamics
inevitably involves the tidal evolution of the Earth–Moon
system (e.g. Goldreich, 1966; Mignard, 1982). The consequence of secular changes in the Earth–Moon system isrecorded in the geological and fossil record, such as in the
striped bands on shell fish (e.g. Williams, 1990; Ito et al.,1993). As is widely accepted now, the Moon was formed bya giant impact of a Mars-sized protoplanet on the protoEarth (e.g. Hartmann et al., 1986; Lee et al., 1997a; Ida
et al., 1997; Canup and Righter, 2000; Kokubo et al., 2000).
Since then, the rotational rate of the Earth has been decreasing and the Moon has been moving gradually awayfrom the Earth. It will soon be possible to numerically calculate the evolution of Earth–Moon dynamics fully including the interaction between the Earth’s rotation and theMoon’s orbital motion for billions of years. Some such trials
have already been done by J. Touma and J. Wisdom using a
sophisticated symplectic integrator (Touma and Wisdom,
1994, 1998, 2001). For this purpose, we need more accurate
information on the time variation of ocean-continent distribution on the Earth as well as that on the Earth’s andMoon’s interior structure so that we can make a precise
dynamical model of the lunar tidal torque.
As for relatively short timescales up to O(103 ) years,
we have to depend on ancient historical records concerning
astronomical phenomena: solar and lunar eclipses, occultations of planets/stars by the Moon. Here the research on the
Earth rotation becomes an interdisciplinary area involving
the study of the ancient history of mankind (e.g. Stephenson, 1997; Tanikawa and Sôma, 2001). Most of these data
come from Chinese chronicles and Babylonian astronomical
diaries written on clay tablets. An accurate determinationof ∆T = TT−UT has just recently begun (e.g. Stephenson
and Morrison, 1995; Stephenson, 1997). Japanese and Chinese contributions to this field are highly anticipated (e.g.
Han, 1997; Tanikawa and Sôma, 2001, 2004a,b; Kawabata
et al., 2002, 2004; Sôma et al., 2004).

87

When the timescale becomes longer, up to several 104
to 105 years, what we can rely on is geological records suchas the fossils of plankton in ocean drilling cores, data in ice
drilling cores from the polar regions, or the annual ringsof trees. These geological data constitute the fundamental
materials underpinning a theory of long-term solar insolation radiation called the “Milankovitch cycles”; these cycles are due to the precession of the Earth’s spin axis combined with the secular change in the Earth’s orbital elements (cf. Milankovitch, 1920, 1930, 1941; Sharaf and Budnikova, 1969a,b; Berger, 1976; Berger et al., 1984; Berger,
1988, 1989; Paillard, 2001). Now the Milankovitch cycles
are being considered for bodies other than the Earth. The
polar caps of Mars might retain some evidence of its climatechange through the severe variation in its obliquity, which
drives some theoretical work based on data obtained byplanet exploring spacecrafts (e.g. Ward, 1974; Ward et al.,
1974; Laskar et al., 2002; Head et al., 2003).
4.5.3. Evolution of planetary rotations
Understanding evolution of planetary spins is an interesting but formidable task for solar system dynamics. Since
the timescale of rotational motion is generally much shorter
than that of orbital motion, it is not easy to investigate rotational motion as accurately and for such long terms as
we do orbital dynamics over the timespan of solar systemevolution. Also in the long-term evolution of planetary rotational and orbital motion, we have to consider the effectof their complicated coupling, so called “spin-orbital coupling” (e.g. Goldreich and Peale, 1966; Tomasella et al.,
1997; Murray and Dermott, 1999; Ward and Canup, 2006).
Partly due to these difficulties, quite a few interesting problems still remain open in this line of research. For example,
the origin of the retrograde rotation of Venus, the originof the synchronization of Mercury’s and the Moon’s rotational motions with their orbital motions, and the originof the inclined rotational axis of Uranus. Readers can consult discussions on the origin and evolution of planetaryspin in the context of planet formation, for example, in
Tanikawa et al. (1991), Ohtsuki and Ida (1998), or Agnoret al. (1999). Agnor et al. (1999) determined the spin angular momentum states of the growing planets by summing
the contributions from each collisional encounter. Their results showed that the spin angular momentum states of the
final planets are generally the result of contributions madeby the last few large impacts. One of the latest treatmentof the evolution of planetary spin is presented as a series ofwork by A. C. M. Correia (Correia and Laskar, 2001, 2003;
Correia et al., 2003) that presented a detailed formulationof spin axis dynamics of Venus with dissipation.
One additional complication to this research field isthe fact that planetary rotations are considered strongly
chaotic (cf. Laskar and Robutel, 1993; Laskar et al.,1993a,b), and the detailed structure of the chaos is not yet
well known. In other words, however, this fact means that

88

Takashi Ito and Kiyotaka Tanikawa

the planetary rotation study can create a series of good
example problems for the research of general dynamicalsystems.
Planetary rotational dynamics is quite relevant to research on extrasolar planetary systems in terms of the question of the possible existence of habitable planets. It is generally agreed that a habitable planet would be in a circumstellar “habitable zone,” on a nearly circular and planar
orbit. The obliquity of the planet should not be so large,
so as to avoid an extreme change of seasons (e.g. Williams
and Kasting, 1997; Williams et al., 1997). Condition for
long-term dynamical stability and instability of planetary
obliquity and its relationship to planetary climate is one of
the key ingredients of the study of habitability in the future extrasolar planet research (e.g. Atobe et al., 2004; Abeet al., 2005; Atobe and Ida, 2007).

5. Numerical methods
Gordon E. Moore, one of the founders of Intel Corporation, made his famous observation in 1965, four years after the first planar integrated circuit was invented. It wascalled “Moore’s Law”. Moore (1965) observed an exponential growth in the number of transistors per integrated
circuit and predicted that this trend would continue. As
he predicted, the development of general-purpose microprocessors has continued, roughly doubling the number oftransistors every 1.5 years for the past 40 years. For example, the Intel 80386 processor introduced in 1985 had only
275,000 transistors per chip, whereas the Pentium 4 processor in 2000 had about 42,000,000 transistors. This trend is
expected to continue until around 2020 when atomic levelaccuracy will be required to produce integrated circuits.
Needless to say, celestial mechanics also owes a great dealof its progress to the progress of this digital technology.
Averaged calculation speed of digital computers has increased nearly a hundred times in the past decade. In addition to the development of integrated circuits describedby Moore’s Law, the architecture of computers has madegreat progress in recent years; from simple scalar computers
through vector computers, to massively-parallel computers.
However, the situation may change from this point on. The
cost of producing new LSIs (large-scale integrated circuits)
is getting higher, with much longer times involved. Designing a large LSI is no longer as easy a task as before. Thesedifficulties make us hesitate to develop new LSIs, leading tohigher production cost of LSIs. Then, the high cost makes
the difficulties of designing new LSIs even greater: This is abad positive feedback. Vector computers are almost dyingout because of their low performance/cost ratio as well as
because of their small market. Very few computer makers
manufacture vector processors now, and it is just a matter of time before vector computers just disappear like dinosaurs.

Under this situation, the future of high-performancecomputing in celestial mechanics will be multidirectional;
entailing sophisticated algorithms to accurately and efficiently solve equations of motion of celestial bodies,
massively-parallel computers with relatively small cost,and special-purpose computers such as GRAPE. Having
this background in mind, we describe in this section somepoints in terms of the development of numerical methods
in celestial mechanics: some issues on regularization, recent progress in symplectic integrators, the special-purpose
GRAPE supercomputer, and several related topics to thisfield. We are quite aware that many interesting and important aspects in this field are missing from this section due
to authors’ inability. For example, we should have mentioned a lot more about symmetric integrators and theirapplication to celestial mechanics as well as those of geometric integrators in a broader sense. Various techniques
to deal with a huge number of particles that interact witheach other through gravitational force, such as the tree
method, are not described. Also, perhaps a brief summary of general development of computer hardware thatis related to celestial mechanics should have been done.
However, we hope that readers are able to get some clues
about these issues from the descriptions and the references
cited in this section.
5.1. Regularization
One of the most difficult tasks in numerical integrationin celestial mechanics is the treatment of collisions. This is
not only a numerical but also a theoretical problem. Regularization is a method used to deal with dynamical collisions numerically as well as analytically.
The restricted three-body problem admits binary collisions. Sundman (1912) has shown that the singularity of theequations of motion corresponding to a binary collision is
not essential, but can be removed via a change of variables.
This process is called a regularization. Regularization of abinary collision has been carried out by Thiele (1896), LeviCivita (1906), Burrau (1906), and Birkhoff (1915). Thiele’s(1896) work concerned an equal-mass binary system, later
extended to arbitrary mass ratio by Burrau (1906). It is tobe noted that these regularizations are all for planar problems. Levi-Civita (1906) provided planes around the singularity. Whenever an orbit moves from one plane to another,
it experiences a collision. Birkhoff (1915) devised a changeof variables so that a binary collision either at the first or
the second primary can be treated through a single formula.
Thiele (1896) and Burrau (1906) transformed the variables
to hyperbolic–elliptic coordinates. The Lemaı̂tre transformation (Lemaı̂tre, 1952) looks like that of Birkhoff, but thecorrespondence of old and new variables is 1 to 2 in thelatter, whereas that is 1 to 4 in the former. Based on thesepioneering studies, Kustaanheimo and Stiefel (1965) havesucceeded in regularizing the three-dimensional problem,

Trends in 20th Century Celestial Mechanics
embedding three variables in the four-dimensional space.
This is the celebrated Kustaanheimo–Stiefel (K–S) regularization (or the K–S transformation). See also Section 5.2
for incorporation of the K–S regularization into symplectic
integrators by S. Mikkola 4 .
Currently, the most significant and promising flow ofregularization research goes on around the K–S transformation. For example, Arakida and Fukushima (2000, 2001)
confirmed that the positional error of a perturbed two-bodyproblem expressed in the K–S variables is proportional to
the fictitious time s, which is the independent variable in
the K-S transformation. This property does not depend on
the type of perturbation, on the integrator used, or on theinitial conditions, including the nominal eccentricity. The
error growth of the physical time evolution and the Kepler
energy is proportional to s when the perturbed harmonic
oscillator part of the equation of motion is integrated by atime-symmetric integration formula, such as the leapfrog
or the symmetric multistep method. The error growth is
proportional to s2 when using traditional integrators, suchas the Runge-Kutta, Adams, Störmer, and extrapolation
methods. Also, it turned out that the K–S regularization
avoids the step size resonance or instability of the symmetric multistep method that appears in the non-regularizedcases. Therefore, the K–S regularized equations of motion
are found quite useful for investigating the long-term behavior of perturbed two-body problems, namely, those used
for studying the dynamics of comets, minor planets, theMoon, and other natural and artificial satellites.
As applications of the K–S regularization, there are several methods to practically regularize triple collision. The
first one is due to Aarseth and Zare (1974). Their idea is
to introduce two K–S regularizations to two of the threeedges of the triangle formed with three particles. The second one is the so-called “chain method” devised by Mikkola
and Aarseth (1990, 1993). In this method, particles forma chain. To make a chain, we start from a closest pairof particles and form a chain made of a vector that connects the pair. Then, a particle closest to either end of thefirst chain is added to form a new chain. Next, the fourthparticle is selected as the closet particle to the third one.
This procedure is continued until all the particles are exhausted. After the completion of the chain, the K–S regularization is applied to each pair of the chain. The third
example is the global three-body regularization that wasdevised by Heggie (1974). Heggie adopted a time transfor4 In all regularizations described in this paragraph, time t is transformed in the form dt = rdτ near collisions where τ is the newtime variable and r is the mutual distance between colliding bodies.
Usually, the solution passes through singularity within a finite time.
For comparison, we should notice that the McGehee transformation(McGehee, 1974) adopts the form dt = r3/2 dτ so that the solution
does not pass through singularity: it takes an infinite time to arrive
at singularity. See also Section 3.1 for the McGehee transformation(or the McGehee variables).

89

mation dt = r1 r2 r3 dτ where ri (i = 1, 2, 3) are the mutual distances of particles. The new Hamiltonian Γ becomesΓ = r1 r2 r3 (H − E0 ) where H is the original Hamiltonianand E0 is its numerical value. Then the equations of motion are regular for collisions between any pair of particles.
This method is accurate in the case of triple close approach
though it is rather time-consuming. Other characteristicsof this method are that the switching of the closest pairs
is not necessary, and that it can in principle be applied to
the N -body problem for any N . For subsequent developments, readers may consult papers by S. J. Aarseth and
S. Mikkola (e.g. Aarseth, 1985, 1988; Mikkola and Aarseth,
1998, 2002), a review paper (Aarseth, 1999), and a textbook (Aarseth, 2003).
Quite recently, Fukushima (2007) announced that hefound a new scheme to regularize a three-dimensional
two-body problem under perturbations. His method is acombination of Sundman’s time transformation and LeviCivita’s spatial coordinate transformation applied to thetwo-dimensional components of the position and velocity
vectors in the osculating orbital plane. The equations of
motion using his new variables have no singularity evenwhen the mutual distance is extremely small. Hence thenew variables are suitable to deal with close encounters.
The number of dependent variables in the new scheme becomes eight, which is significantly smaller than the existing
schemes to avoid close encounters: for example, the number
is smaller by two compared with the K–S regularization.
5.2. Symplectic integrators
Since many of the systems that celestial mechanics isconcerned with are Hamiltonian systems or their proximity, numerical integration schemes designed specificallyto maintain the Hamiltonian structure are evidently desirable and promising. Symplectic integrators are exactly
the methods that satisfy such requirements. In general,
symplectic integrators conserve the total energy of a system quite well, preventing artificial damping or excitation
due to the accumulation of local truncation errors (e.g.
Yoshida, 1990b; Gladman et al., 1991; Kinoshita et al.,
1991; Yoshida, 1993; Sanz-Serna and Calvo, 1994; Morbidelli, 2002). The mathematical background of symplectic integrators is based on Lie algebra. Its application tophysics is described, for example, in Varadarajan (1974),
Dragt and Finn (1976), Ruth (1983), or Yoshida (2001).
Yoshida (1990b) first noticed that a higher-order symplectic scheme can be composed of several second-orderschemes, and derived some practical pairs of coefficientsfor higher-order integrators. Yoshida also proved that aclass of explicit symplectic integrators rigorously preserve
the total angular momentum within the range of roundoff errors when used for gravitational N -body problems(Yoshida, 1990a).
In practice, general-purpose symplectic integrators

90

Takashi Ito and Kiyotaka Tanikawa

tend to be of low order, mainly because of their complexity in higher orders. Therefore they are not suitable forlong-term numerical simulations. However, the so-called
Wisdom-Holman type symplectic map (“WH map”) devised by Wisdom and Holman (1991), which is specialized
for slightly perturbed motion like solar system planetary
motion, is more accurate by orders of magnitude than thegeneral-purpose symplectic method even when it is of loworder. The principle of the WH map is to split the Hamiltonian into an unperturbed part (for example, Kepler motion when we trace planetary motion) and a perturbation
part (for example, mutual gravitational perturbation). In
each step of the integration, the system is first drifted
forward in time according to unperturbed motion, andthen a kick in momentum is applied that is derived fromthe perturbation part of the Hamiltonian. The first stepis analytic because the Hamiltonian is unperturbed andintegrable. When this part deals with Keplerian drift, a setof canonical coordinates, such as Jacobi coordinates (e.g.
Plummer, 1960), is used to keep the integration scheme
symplectic. Using this scheme, the WH map enhances the
speed of long-term numerical integrations of planetary motion by orders of magnitude. Wisdom and Holman (1992)
investigated the numerical stability of the WH map, anddiscussed the existence of numerical instability within themap. Note that when using the WH map for planetary
dynamics, we have to solve the Kepler equation every timein the drift phase. As for the fast solving methods of the
Kepler problem, we may find it helpful to consult a seriesof papers by T. Fukushima (1996a, 1997a,b,d, 1998).
After the monumental work by Wisdom and Holman(1991), a lot of enhancements have been added to the WH
map. Based on the idea of the WH map, Saha and Tremaine(1992) discovered a special start-up procedure to reduce
the truncation error of angle variables, called the “warmstart”, exploiting one of the characteristics of Hamiltoniansystems—the existence of adiabatic invariants. They also
discovered a symplectic scheme with individual stepsizes,
dividing the total Hamiltonian into each planet’s Kepler
and perturbation parts (Saha and Tremaine, 1994). A technique related to the warm start was devised by Wisdom
et al. (1996) under the name of the “symplectic corrector”.
The symplectic corrector exploits the fact that the valuesof symplectically integrated systems (i.e. coordinates andmomenta) are obtained by canonical transformations of thevariables in real systems. See also Wisdom (2006) for application of the symplectic corrector to the N -body systemdescribed by canonical heliocentric coordinates. Michel andValsecchi (1997) discussed in detail the efficiency of the WH
map compared with traditional integration schemes suchas the Bulirsh-Stoer extrapolation method. For higher order symplectic methods for nearly integrable Hamiltoniansystems, consult recent publications such as Chambers andMurison (2000) or Laskar and Robutel (2001).
Now one of the biggest issues concerning symplectic in-tegrators is how to handle close encounters between particles where the variation of stepsize is inevitable, because
symplectic integrators are inherently not well matched withvariable stepsize schemes (Yoshida, 1993). Based on somepreceding studies (Skeel and Biesiadecki, 1994; Lee et al.,1997b), Duncan et al. (1998) has introduced a method tosolve the problem with the aid of “democratic heliocentric
coordinates”, also known as canonical heliocentric coordinates that use heliocentric positions and barycentric momenta (cf. Charlier, 1902). Their method enables us to handle close encounters among particles by splitting gravitational potentials recursively. The method maintains computational speed equivalent to that of the generic WH map
when no close encounter occurs. When a close encounter
happens, the scheme automatically divides the stepsize andalleviate numerical errors. Duncan et al. (1998) named thescheme SyMBA, or “symplectic massive body algorithm”.
Their program package is named SWIFT, and is available
on their ftp site. One of the major disadvantages of using SyMBA is that it cannot take care of highly eccentric
orbits, i.e. close encounters with central mass. The otherproblem is the complexity of its implementation.
Chambers (1999) originally devised a method similar to
SyMBA but with a slightly different idea. He noticed thatthe drift part of the symplectic scheme does not need to be
calculated analytically, but can be integrated numericallyusing a very accurate method such as the extrapolation
scheme. When a close encounter happens, his scheme moves
the relevant terms in the perturbation Hamiltonian Hpert
to the Keplerian Hamiltonian Hkep , and keeps the contrastof the order of magnitude of Hamiltonians, Hkep À Hpert .
This hybrid scheme enables us to calculate the orbital motion in the system where close encounters frequently occur, with less programming cost than SyMBA. Chambers
et al. (2002) later extended this hybrid method to a scheme
specialized for orbital motions in and around binary stars,
choosing suitably canonical variables. This method is available for integrating planetesimals’ and protoplanets’ orbital
evolution near binary star systems (e.g. Quintana et al.,
2002; Quintana and Lissauer, 2006; Thébault et al., 2004).
Incorporating the above idea of Chambers et al. (2002)
into SyMBA, Levison and Duncan (2000) has developed amodified version of SyMBA that can take account of close
encounters with the central mass, i.e. highly eccentric orbits. Their method is available for systems where close encounters between particles as well as between a particle andthe central mass frequently occur. When neither a close encounter happens nor the orbits of particles are so eccentric,
the modified SyMBA performs as fast as the generic WH
map does. Except for the complexity of its implementation,
the modified SyMBA is one of the ultimate schemes in symplectic integrators that are suitable for planetary dynamics
in planet formation stage.
In a practical use of symplectic integrators for highly
eccentric orbital motions, a series of studies on regularizedTrends in 20th Century Celestial Mechanics
symplectic integrators by Mikkola and collaborators seems
quite promising (Mikkola, 1997, 1999; Mikkola and Innanen, 1999; Mikkola and Tanikawa, 1999a,b; Mikkola andPalmer, 2000). They have brought the K–S regularization
into the WH map, keeping its symplectic structure intact,
though regularization had long been considered incompatible with symplectic integration. Their basic idea is to takethe logarithm of the regularized Hamiltonian. In order toregularize the equations of motion, one multiplies distancer to the original Hamiltonian, which is equivalent to division of the Hamiltonian by the dominant potential. Thismeans that the new Hamiltonian is not separable, i.e. coordinates and momenta are mixed, even when the original Hamiltonian is separable. Taking the logarithm of the
new Hamiltonian, we can separate variables and exploit theleapfrog procedure for symplectic integration.
Mikkola’s (1997) method is not only accurate when
adapted for highly eccentric orbits, but also quite fast, since
there is no need to solve Kepler’s equation which needs
certain iterative procedures involving transcendental functions. Compared with similar previous work, such as Levison and Duncan’s (1994) “RMVS” (regularized mixed variable symplectic integrator) by changing stepsizes automatically and by force-center switching, Mikkola’s method ismore efficient, more accurate, and easier to implement. Breiter (1999) also found a regularized method using K–S variables, which he insists is a similar but complementary approach to the work of Mikkola. See also Preto and Tremaine(1999) for a class of symplectic integrators with adaptive
stepsize for separable Hamiltonian systems.
From this point of view, Rauch and Holman’s (1999)
intensive review is quite worth reading. They have tested
many variants of the WH map and measured their efficiency. They have concluded that an enhanced version of
a potential splitting method such as SyMBA, incorporating time regularization, force-center switching, and an improved kernel function, may be quite efficient when particles are subject to both high eccentricities and close encounters. To integrate eccentric and nearly Keplerian orbits
without close encounters, the time-transformed WH map
by Mikkola is clearly the most efficient and stable.
Symplectic integrators can be extended to include dissipative mechanisms. Malhotra (1994) devised such a method
to integrate the planetary orbital motion dragged by gas,
adding a simple modification to the generic WH map so thatit can include velocity-dependent forces. Mikkola (1998)
discussed the inclusion of non-canonical perturbations in
symplectic integration schemes using the δ-function formalism. J. Touma’s series of work (Touma and Wisdom,
1994, 1998, 2001) demonstrates the feasibility of exploiting
symplectic integrators to calculate systems with long-term
tidal dampings, such as the Earth–Moon system and thecore-mantle coupling of the Earth’s interior.
The research devoted to the enhancement and improvement of symplectic integrators, especially those of91

the WH map, is still going on. For example, Liao (1997)
gave schemes for a low-order mixed symplectic integrator for an inseparable, but nearly integrable Hamiltoniansystem. Although Liao’s schemes are implicit, they havea faster convergence rate of iterative solutions than ordinary implicit integrators do. Fukushima (2001) proposedand examined three kinds of approach to reduce the accumulation of round-off errors in symplectic integrators: thereduction of the number of summations in the implementation of symplectic algorithms, the use of a double-length
routine library in the main summation procedure, and thefull use of the double-length routine library in the entire
procedure of symplectic integration (note that the last
two methods will not be applicable to the WH map unlessdouble-length routines to evaluate trigonometric functions
are provided). Beust (2003) presented a new symplectic
integrator for massive bodies that permits the numericalintegrations of the dynamics of hierarchical stellar systems of any size and shape, provided that the hierarchicalstructure of the systems is preserved along the integration. Chambers (2003) noticed that the previously known
symplectic methods that are higher than second order
include some substeps that travel backward, comparedwith the main integration. To compensate for this, some
substeps in these methods have large coefficients, which
could produce large error terms and reduce the efficiencyof high-order symplectic algorithms. Chambers introduced
complex coefficients in higher order symplectic integrators from the solution of the constraint equations for thesubsteps of symplectic algorithms. He showed that complex integrators with leading error terms that have strictly
imaginary coefficients effectively behave as if they are oneorder higher than expected.
Symplectic integrators can be improved in terms of force
calculation. As is well known, the amount of direct computation of the distances between each pair of gravitationallyinteracting N particles grows as O(N 2 ). To overcome thisdifficulty, a variety of fast (with effort growing more slowly
than O(N 2 )) but approximate force calculation methods
such as those used in the tree code or fast multipole methods have been developed. Wiegert et al. (2004) examinedseveral of these algorithms together with symplectic integrators, and compared their speed in very detail. Also, combination of the special-purpose computer GRAPE (Section 5.4) with symplectic integrators, especially those of
SyMBA-type, might have a possibility to greatly enhance
the computing speed of gravitational N -body interactions
in solar system dynamics.
5.3. Other promising numerical techniques
There will be two major streams of numerical integration schemes used in celestial mechanics. One of them is the
symplectic integrators that we have already mentioned in
the preceding section, and the other is symmetric integra-

92

Takashi Ito and Kiyotaka Tanikawa

tors (cf. Hairer et al., 1993, 2002). Symplectic integrators
will be developed in various ways specialized for each dynamical problem, such as versions for binary systems, threebody systems with regularization, systems with frequent
close encounters, and so on. Symmetric methods have a superiority over the symplectic integrator in that it is easy toconstruct variable stepsize schemes. Among many kinds ofsymmetric integrators, fourth-order time-symmetric methods based on the Hermite integrator seem to be full ofpromise, especially for problems of planetary accretion orstellar dynamics where frequent close encounters require
variability in stepsizes (e.g. Makino, 1991b; Makino andAarseth, 1992; Kokubo et al., 1998; Kokubo and Makino,2004). The Hermite scheme is also a good match with thespecial-purpose GRAPE computer for the gravitational N body problem (Makino, 1991a).
While the Hermite-type integrators are one-step methods, higher-order multistep symmetric integrators are also
widely used in planetary dynamics where a very high accuracy is required. Following a pioneering work by Lambert and Watson (1976), Quinlan and Tremaine (1990) proposed an efficient and useful 12th-order scheme, which is
now a standard of the higher-order linear multistep symmetric method. Later, Evans and Tremaine (1999) as well
as Fukushima (1999b) investigated the numerical instability of the multistep symmetric method independently, andproposed several improved schemes. Though it is not easy
to implement a variable stepsize scheme in multistep symmetric methods, they are quite suitable for planetary dynamics with regular and smooth motion without close encounters between planets. One of the current major problems of multistep symmetric methods is that, they are subject to intrinsic stepsize resonances and instabilities whichwere first recognized by Alar Toomre (cf. Quinlan, 1999;
Fukushima, 1999b). Hence, conventional Störmer–Cowellmethods (cf. Hairer et al., 1993) may still be more widelyused than multistep symmetric methods.
Dynamical models that we deal with in the solar system celestial mechanics is often nearly integrable, such as
the quasi-Keplerian motion. This characteristic has been
employed not only by the Wisdom–Holman symplectic
map, but using traditional methods such as the Encke
method (Encke, 1854; Fukushima, 1996b). Recent works by
T. Fukushima, collectively called “the method of manifold
correction”, also exploits this quasi-Keplerian character ofsolar system dynamics. Fukushima (2003c) proposed a new
approach to numerically integrate quasi-Keplerian orbits.
His method integrates the time evolution of the Kepler
energy and the usual equations of motion simultaneously.
It directly adjusts the integrated position and velocity bya spatial scale transformation in order to satisfy the Kepler energy relation rigorously at every integration step.
By adding the Laplace integral as the second auxiliary
quantity to this method, Fukushima (2003a) extended his
scaling method to integrate quasi-Keplerian orbits in or-

der to suppress the growth of integration errors, not onlyin the semimajor axis, but also in the other orbital elements, especially in the eccentricity and in the longitudeof pericenter. Then again, by adding the orbital angularmomentum vector as another auxiliary quantity to be integrated, Fukushima (2003b) extended his scaling methods
for quasi-Keplerian orbits in order to suppress the growthof integration errors in the inclination and the longitudeof the ascending node. His new approaches provide a fastand high-precision device to simulate the orbital motions
in solar system dynamics at a negligible increase in computational cost. Fukushima’s manifold correction method,combined with the Kustaanheimo–Stiefel regularization,
is now creating a new, substantially large flow in this lineof research (Fukushima, 2004a,b,c,d,e,f, 2005a,b,c,d).
Another important issue concerning numerical integration techniques is parallel computing. Since many problems that are dealt with in celestial mechanics do not involve very large number of bodies such as the solar systemplanetary motion, efficient parallelization of computing is
not easy. However, recent investigation suggests a possibility for us to utilize a large parallel computer for a relatively small number of celestial bodies. For example, Itoand Fukushima (1997) developed a parallelized version of
the extrapolation method to integrate general ordinary differential equations. A simple technique made its load balance among computer processors almost equal, nearly independent of how large amount of data (or how many particles) we use in our numerical model. Also, a method developed by Fukushima (Fukushima, 1997c,e, 1999a) expands
an approximate solution of the equations of motion of celestial particles into Chebyshev polynomials, and approaches
the true solution through the Picard iteration. Though the
Picard–Chebyshev method is not yet being widely applied
to real orbital problems (cf. Arakida, 2006), this method is
expected to be quite suitable for systems with small perturbation, such as orbital and rotational motion in our solar
system calculated with large parallel computers. Similar to
the Picard–Chebyshev method, Saha and Tremaine (1997)
have devised a new parallel computing algorithm for planetary orbital motion, especially suitable for scalar massive
parallel computers. In addition, the near-future version of
the special-purpose computer GRAPE, GRAPE–DR, will
enable us to run a large number of small-scale gravitational
N -body integrations all at once, which will practically bequite a massively parallelized calculation (see also Section
5.4).
Yet another important issue along this line of research
can be the multiple precision arithmetic. Long-term numerical integration potentially and inevitably involves alarge amount of round-off error. As mentioned before, several methods to alleviate the round-off error have been developed by Quinn and Tremaine (1990) and Fukushima(2001). Although their methods are algorithmically realized (i.e. by software), the hardwired implementation ofTrends in 20th Century Celestial Mechanics
multiple precision might now be possible for astronomical application. The implementation of hardwired multiple
precision would also be quite influential on numerical experiments in theoretical studies, such as in helping to gain
information on the detailed structure of various KAM tori,
aiding in inspection of dependency of the Lyapunov exponents on arithmetic precision, and furthering many otherproblems in celestial mechanics and dynamical theory.
5.4. GRAPE: A special-purpose computer for N -bodydynamics
As for numerical research in the fields of celestial mechanics and dynamical astronomy, we cannot avoid mentioning the enormous success of GRAPE (GRAvity piPE),a special-purpose supercomputer for gravitational N -bodyproblems developed by a dedicated research group in Japan(e.g. Sugimoto et al., 1990; Ebisuzaki et al., 1993; Makino
et al., 1997; Makino and Taiji, 1998). GRAPE has created a revolutionary breakthrough in the area of gravitational N -body problems. The peak speed of its latest version, GRAPE–6, has a capacity of more than several tensof TFlops. The near-future goal of the GRAPE project iscalled GRAPE–DR, which will have a peak capacity of 2
PFlops or higher (Makino, 2005).
The essential feature of the GRAPE concept is thepipeline. On GRAPE, computation of the N -body problem
is divided into two parts. The first part is performed on ahost computer (such as a PC), and the other is on a specialpurpose back-end processor (GRAPE). GRAPE calculates
the gravitational force among particles through pipeline(s).
During one timestep, the host computer performs O(N )operations, while GRAPE can do O(N 2 ) because of the
devoted parallel pipeline. Applegate et al. (1986)’s “TheDigital Orrery” also has a parallel architecture in whichsmall computers are connected in a ring. Each processor
takes care of one planet (particle). The Digital Orrery did
not adopt a pipeline architecture because this architecture
lacks the efficiency necessary for small N . Pipelines work
more efficiently when N is sufficiently large.
The GRAPE project team started the development of
GRAPE systems in 1989. The first hardware, GRAPE–1
was a single-pipeline system with low-accuracy force calculation (Ito et al., 1990). It was a machine for demonstratingthat theoretical astrophysicists could construct a practical
hardware. It had a peak speed of 240 MFlops, and turned
out to be useful for a wide range of astrophysical problems
where high accuracy is not very important. GRAPE–2 wasagain a single-pipeline system, but with high accuracy, witha peak speed of 40 MFlops (Ito et al., 1991). GRAPE–3
was the first multiple-pipeline system with a low-accuracy
force calculation pipeline (Okumura et al., 1993). The peak
speed of GRAPE–3 was 380 MFlops/chip.
As a successor of the high-precision machine GRAPE–2,
GRAPE–4 was developed (Makino et al., 1997). The peak93

speed of the GRAPE–4 system is 1.08 TFlops, achieved byrunning 1692 pipeline LSIs in parallel, each providing 640
MFlops. The successor of the low-precision GRAPE–3 machine is GRAPE–5 (Kawai et al., 2000). A GRAPE–5 board
has a peak performance of 38.4 GFlops. The newest editionof the GRAPE system is GRAPE–6 (Makino et al., 2003),
the successor of the high-precision machine (GRAPE–4),
and its single PCI card version, GRAPE–6A (Fukushige
et al., 2005). Similar to GRAPE–4, the primary applicationof GRAPE–6 is simulations of collisional systems, though
it can also be used for collisionless systems. Various improvements have been done for GRAPE–6 from GRAPE–
4, which ended up with a peak speed of 64 TFlops.
We must also mention that the GRAPE system
has the capability of hardwired implementation of particular algorithms for astronomical applications. Thiskind of GRAPE is called PROGRAPE (Hamada et al.,2000). Now PROGRAPE–1 (PROgrammable GRAPE–1), a programmable multi-purpose computer for manybody simulations, is available. The main difference between PROGRAPE–1 and traditional GRAPEs is that
PROGRAPE–1 uses FPGA (Field Programmable Gate
Array) chips as its processing elements, while the traditional GRAPEs rely on a hardwired pipeline processor
specialized for gravitational interactions. GRAPE–7, thesuccessor of the lower-precision machine (GRAPE–5) but
now with FPGA, will also be available shortly (Kawai and
Fukushige, 2006).
The GRAPE project is tremendously successful inthe field of gravitational N -body dynamics, especiallyfor the planet formation and collisional stellar dynamics. Hundreds of academic papers have been publishedusing GRAPE computation results, and the number isever increasing. See the GRAPE project’s webpage for amore detail (http://www.astrogrape.org/ for general information, and http://grape-dr.adm.s.u-tokyo.ac.jp/ for the
GRAPE–DR project).
6. Concluding Remarks
As we have seen, since the beginning of the twentiethcentury, celestial mechanics has pretty much diverged into
many specialized fields that do not communicate much witheach other. Despite the apparent resulting prosperity of celestial mechanics, we suspect that this diversity could produce its potential stagnation through a lack of communication between its various fields. Perhaps it might be uselessto imagine how many people have read and fully understood
the five texts that all have the phrase “celestial mechanics”
in their titles 5 : Smart (1953) Celestial Mechanics, Siegel
5

The fact that popular textbooks such as listed here are mainly
written in English allows them to be widely distributed and read
all over the world. However, we should not forget that there arealso many good textbooks written in other languages. The speed ofunderstanding is clearly fastest when readers obtain information in

94

Takashi Ito and Kiyotaka Tanikawa

and Moser (1971) Lectures on Celestial Mechanics, Taff(1985) Celestial Mechanics: A Computational Guide for
the Practitioner, Brumberg (1995) Analytical Techniquesof Celestial Mechanics, and Poincaré (1892) Les MéthodeNouvelles de la Méchanique Céleste.
At this point, we would like to ask ourselves a very
primitive question to which we have not found an answer:
What is celestial mechanics? What does celestial mechanics
aim at? It seems that there is no general agreement on theanswer to these questions. However, this situation may not
be so strange: Compared with other fields in astronomy,
the disagreement among celestial mechanists on the subject(or target) of celestial mechanics looks remarkably large,
partly because this field has a much longer history thanother fields do. This is what has made celestial mechanics
branch out into so many and such detailed paths.
Astronomy has expanded by an inconceivable degree inthe twentieth century. Researchers extended the frontier ofastronomy to various electromagnetic wavelengths in thefirst half of the twentieth century, beginning from visible
light, reaching radio, infrared, ultraviolet, X-ray, and γ-ray.
Astronomy in the twentieth century was full of discoveries
and astonishment, most of which were brought about using the methods of modern astrophysics. Celestial mechanics has rolled up and down in the midst of the killer waves
of astrophysics. When we remind ourselves of the ages before recent discoveries in astrophysics, it is not wrong to
say that celestial mechanics was the unique field of quantitative science that covered the entire world known to humankind: In the eighteenth century, the solar system wasalmost the universe itself. In the nineteenth century, there
was no notion of extra galaxies, and our galaxy was theuniverse. Spectroscopy of stars started only at the end of
the nineteenth century. No one knew of nuclear reactions,and astrophysical studies of extrasolar stars were just beginning. These facts imply that celestial mechanics was at
the center of the astronomy of those ages. With those centuries behind us, now the world is complicated enough tomake it virtually impossibly to grasp the appropriate position of celestial mechanics in the vast sea of astronomy.
As it is not difficult to imagine, we are not the firstto suffer from this kind of confusion. Already in the early
1920s when Brown et al. (1922) reported the activities ofAmerican celestial mechanists, they were already not suretheir mother language. In this sense, we much prefer the existence ofgood textbooks written in various languages. However, as expected,
non-English textbooks have only a limited number of readers. We
do not think this is desirable. Good textbooks should be read by
as many readers as possible, regardless of the languages in which
they are written. One way to overcome this dilemma is to write
your textbook in your own language as well as in English if you
are the author. It cannot be very difficult to translate your own
book into English. Fortunately, textbooks for celestial mechanics are
generally full of mathematical formulas, which means you likely havelittle to translate into English. This seems only possible through thevoluntary effort of ungrudging scientists.

what the main purpose of celestial mechanics was. Their reports consisted of three parts: the first one was devoted tothe solar system, the second one was to celestial mechanics
as applied to the stars, and the third one was to the theory of the problem of three and more bodies. At the very
beginning of their report, Brown wrote:“Celestial mechanics, broadly interpreted, is involved inpractically all the astronomy of the present time. The
limited meaning of the term now usually adopted refers
only to those problems in which the law of gravitation
plays the chief or only part, and more particularly tothose which deal with motions of bodies about one another and with their rotations.”
which obviously indicates that Brown felt the same perplexity as what we feel now when he tried to define celestialmechanics at that stage of the early twentieth century.
When we are involved with celestial mechanics, an important attitude we should take is to always keep an eye on
the concrete, actual, real phenomena in the universe, evenwhen the methodology we use is highly mathematical or
abstract. We should keep in mind that celestial mechanics
is aimed at better understanding the real structure of thereal world, not at fiddling with complicated equations ormanipulating high-speed computers as a hobby. With thisin mind, we will not lose the loose ends of celestial mechanics in our effort to understand the behavior of this world interms of dynamics, no matter how varied the phenomena
we encounter are.
Before the end of this manuscript, we would like topropose three potential near-future projects that are tightly
related to the issues that we have so far discussed for thefurther development of celestial mechanics in the twentyfirst century:(i) To build a computer chip at least with 10,000-digit
arithmetic, in order to observe complicated chaotic
behaviors in dynamical systems such as the Arnold
diffusion. As an example, take discrete dynamical systems, i.e. maps. At the worst case, numerical errors
can be doubled per one iteration. This means that numerical errors can be multiplied thousand times overten iterations. In other words, we can iterate the map
thirty thousand times before we totally lose accuracy
if we have ten thousand digits. This will enable us todirectly calculate topological entropy of general dynamical systems. We expect that million-digit arithmetic may be possible by the end of the twenty-first
century.
(ii) Occultation observations by the Moon-orbiting satellite on polar orbit, in order to survey multiple stars
in our galaxy (cf. Tanikawa and Mikami, 2000). This
will provide us unbiased data of visual binaries withina relatively short time, which will be the fundamentals not only for star formation study but for the celestial mechanics, particularly the dynamical study

Trends in 20th Century Celestial Mechanicsof binary formation in primordial star clusters.
(iii) To attempt space trigonometry with baseline ofO(100) AU using two pairs of spacecrafts: the firstpair orbits near around the ecliptic, and the other
orbits in a plane that is perpendicular to the ecliptic. The purpose of this project is to make a threedimensional map of our galaxy with the spatial scaleof 104 pc and with the accuracy of 000 .01 or higher,
which will be the basis of future stellar dynamics.
We will be able to get proper motion data over 1–2
years from this mission. Three-dimensional position
data together with velocity data (i.e. proper motion
and radial velocity) will provide us basic information
about the dynamics of our galaxy. The advantage ofthis project is that trigonometric parallaxes can beobtained at once: we do not need to wait half a year.
A disadvantage is that since the spacecrafts move soslowly that the configuration of stars and spacecrafts
does not change so quickly. This disadvantage canbe avoided, however, by increasing the number of
spacecrafts in a pair (a group) from two to four.
In 1974, M. W. Hirsch and S. Smale commented on Newton’s gravitational equations of motion in their celebrated
textbook (Hirsch and Smale, 1974, p. 289) thus:“. . . The flow obtained from this differential equation
then determines how a state moves in time, or the life
history of the n bodies once their positions and velocities are given. Although there is a vast literature of several centuries on these equations, no clear picture has
emerged. In fact it is still not even clear what the basic
questions are for this “problem”.”
We think this sarcastic and rather cynical comment reflects the huge potential of this field that has not yet been
fully examined. The need for celestial mechanics will be
felt to an even greater degree as human-being launch more
artificial satellites and spacecrafts, recalling that the ageof modern celestial mechanics began with Sputnik in the
1950s. Human beings might need to protect the Earth fromthe potential impacts of near-Earth objects or space debris, and again the knowledge and experience of celestialmechanics will be required with urgency. The future of celestial mechanics is closely connected with engineering andspace technology, employed both inside and outside our solar system. We are grateful to mother nature for affording
us the opportunity to witness the dawn of celestial mechanics in the new century.
Acknowledgments
The authors wish to thank Jun Makino and Hidekazu
Ito for giving us numerous comments which suggested directions that made far better the quality of this paper. We also
have greatly benefited from stimulating discussions withand encouragement from Jack J. Lissauer, Arika Higuchi,

95

Yoshihiro Yamaguchi, and Hiroaki Umehara. A review by
Yolande McLean has considerably improved the English
presentation of this paper. The LATEX 2ε class file that theauthors modified and used for this manuscript is based on
elsart5p.cls which is provided and copyrighted by Elsevier.

References
Aarseth, S.J. (1985) Direct method for N -body simulations, in Blackbill, J.U. and Cohen, B.I. eds., Multiple
Time Scales, Academic Press, New York, 377–418.
Aarseth, S.J. (1988) Integration methods for small N -bodysystems, in Valtonen, M.J. ed., The Few Body Problem,
Kluwer Academic Publisher, Dordrecht, 287–307.
Aarseth, S.J. (1999) From NBODY1 to NBODY6: The
growth of an industry, Publ. Astron. Soc. Pac., 111,
1333–1346.
Aarseth, S.J. (2003) Gravitational N -Body Simulations:Tools and Algorithms, Cambridge University Press, Cambridge, England.
Aarseth, S.J. and Lecar, M. (1975) Computer simulationsof stellar systems, Ann. Rev. Astron. Astrophys., 13, 1–
21.
Aarseth, S. and Zare, K. (1974) A regularization of thethree-body problem, Celes. Mech., 10, 185–205.
Abe, Y., Numaguti, A., Komatsu, G., and Kobayashi, Y.
(2005) Four climate regimes on a land planet with wet
surface: Effects of obliquity change and implications forancient Mars, Icarus, 178, 27–39.
Agnor, C.B., Canup, R.M., and Levison, H.F. (1999) On
the character and consequences of large impacts in thelate stage of terrestrial planet formation, Icarus, 142,
219–237.
Alekseev, V.M. (1968a) Quasirandom dynamical systems,
I., JMath. USSR–Sb, 5, 73–128.
Alekseev, V.M. (1968b) Quasirandom dynamical systems,
II., JMath. USSR–Sb, 6, 505–560.
Alekseev, V.M. (1969) Quasirandom dynamical systems,
III., JMath. USSR–Sb, 7, 1–43.
Alekseev, V.M. (1981) Quasirandom oscillations andqualitative questions in celestial mechanics, in Kusnirenko, A.G., Katok, A.B., and Alekseev, V.M. eds.,
Three Papers on Dynamical Systems, Vol. 116 of AMS
Translations, American Mathematical Society, Providence, Rhode Island, 97–169, originally in Russian.
Allen, R.L., Bernstein, G.M., and Malhotra, R. (2002) Observational limits on a distant cold Kuiper belt, Astron.
J., 124, 2949–2954.
Anosova, J.P. (1996) Dynamical evolution of triple systems
with big differences in the masses of bodies: a criterion
for stability, Astrophys. Space Sci., 238, 223–238.
Applegate, J.H., Douglas, M.R., Gürsel, Y., Sussman, G.J.,
and Wisdom, J. (1986) The outer solar system for 200
million years, Astron. J., 92, 176–94.

96

Takashi Ito and Kiyotaka Tanikawa

Arakida, H. (2006) Application of Picard–Chebyshev
method to orbital dynamics, in Proceedings of N -bodyDynamics Symposium, Kashima Space Research Center,
National Institute of Information and Communications
Technology, Itako Fujiya Hotel, Itako, Ibaraki, Japan,
March 1–4, 2006, 114–136.
Arakida, H. and Fukushima, T. (2000) Long-term integration error of Kustaanheimo–Stiefel regularized orbital
motion, Astron. J., 120, 3333–3339.
Arakida, H. and Fukushima, T. (2001) Long-term integration error of Kustaanheimo–Stiefel regularized orbital
motion. II. Method of variation of parameters, Astron.
J., 121, 1764–1767.
Arnold, V.I. (1963) Small denominators and the problemsof stability of motion in classical and celestial mechanics, Russ. Math. Surveys, 18, 85–191, (English translation from Report to the IVth All-Union Mathematical
Congress, Leningrad, 85–191, 1961).
Arnold, V.I. (1964) Instability of dynamical systems withseveral degrees of freedom, Soviet Math. Dokl., 5, 581–
585, (English translation).
Arnold, V.I. (1978) Mathematical Methods in Classical Mechanics, Springer-Verlag, Berlin.
Arnold, V.I., Kozlov, V.V., and Neishtadt, A.I. (1993)
Mathematical aspects of classical and celestial mechanics(Dynamical Systems III), Springer-Verlag, Berlin, second
edition.
Atobe, K. and Ida, S. (2007) Obliquity evolution of extrasolar terrestrial planets, Icarus, 188, 1–17.
Atobe, K., Ida, S., and Ito, T. (2004) Obliquity variations ofterrestrial planets in habitable zones, Icarus, 168, 223–
236.
Aubry, S. and Abramovici, G. (1990) Chaotic trajectories
in the standard map: the concept of anti-integrability,
Physica D, 43, 199–219.
Aubry, S. and LeDaeron, P.Y. (1983) The discrete FrenkelKontorova model and its extensions, Physica D, 8, 381–
422.
Backman, D.E. and Paresce, F. (1993) Main-sequence stars
with circumstellar solid material — The VEGA phenomenon, in Levy, E.H. and Lunine, J.I. eds., Protostars& Planets III, The University of Arizona Press, Tucson,
Arizona, 521–542.
Bailey, M.E. and Emelýanenko, V.V. (1996) Dynamical
evolution of Halley-type comets, Mon. Not. R. Astron.
Soc., 278, 1087–1110.
Bailey, M.E., Chambers, J.E., and Hahn, G. (1992) Origin of sungrazers: a frequent cometary end-state, Astron.
Astrophys., 257, 315–322.
Bangert, V. (1988) Mather sets for twist maps andgeodesics on tori, Dynamics Reported, 1, 1–56.
Battin, R.H. (1987) An Introduction to the Mathematics
and Methods of Astrodynamics, American Institute of
Aeronautics and Astronautics, Inc., New York.
Benedicks, M. and Carleson, L. (1991) Dynamics of theHénon map, Ann. Math., 133, 73–169.
Benettin, G., Galgani, L., and Strelcyn, J.-M. (1976) Kolmogorov entropy and numerical experiments, Phys. Rev.
A, 14, 2338–2345.
Berger, A.L. (1976) Obliquity and precession for the last
5000000 years, Astron. Astrophys., 51, 127–135.
Berger, A.L. (1988) Milankovitch theory and climate, Rev.
Geophys., 26, 624–657.
Berger, A.L. ed. (1989) Climate and Geo-Science — A Challenge for Science and Society in the 21st Century, Kluwer
Academic Publishers, Dordrecht.
Berger, A.L., Imbrie, J., Hays, J., Kukla, G., and Saltzman, B. eds. (1984) Milankovitch and Climate — Understanding the Response to Astronomical Forcing, D. Reidel, Norwell, Mass.
Beust, H. (2003) Symplectic integration of hierarchical stellar systems, Astron. Astrophys., 400, 1129–1144.
Birkhoff, G.D. (1913) Proof of Poincaré’s geometric theorem, Trans. Amer. Math. Soc., 14, 14–22.
Birkhoff, G.D. (1915) The restricted problem of three bodies, Rend. Circ. Mat. Palermo, 39, 265–334.
Birkhoff, G.D. (1920) Surface transformations and theirdynamical applications, Acta Math., 43, 1–119.
Birkhoff, G.D. (1925) An extension of Poincaré’s last geometric theorem, Acta Math., 47, 297–311.
Birkhoff, G.D. (1927) Dynamical Systems, American Mathematical Society, Providence, Rhode Island.
Black, D.C. (1982) A simple criterion for determining the
dynamical stability of three-body systems, Astron. J.,
87, 1333–1337.
Blaes, O., Lee, M.H., and Socrates, A. (2002) The Kozai
mechanism and the evolution of binary supermassive
black holes, Astrophys. J., 578, 775–786.
Boccaletti, D. and Pucacco, G. (1996) Theory of Orbits. 1. Integrable systems and non-perturbative methods,
Springer-Verlag, Berlin.
Boccaletti, D. and Pucacco, G. (1998) Theory of Orbits. 2.
Perturbative and geometrical methods, Springer-Verlag,
Berlin.
Bogard, D.D. (1995) Impact ages of meteorites: A synthesis,
Meteoritics, 30, 244–268.
Bottke, W.F., Vokrouhlický, D., Brož, M., Nesvorny, D.,
and Morbidelli, A. (2001) Dynamical spreading of asteroid families by the Yarkovsky effect, Science, 294, 1693–
1696.
Bottke, W.F., Vokroulický, D., Rubincam, D.P., andBrož, M. (2002) The effect of Yarkovsky thermal forces on
the dynamical evolution of asteroids and meteoroids, in
Bottke, W.F., Cellino, A., Paolicchi, P., and Binzel, R.P.
eds., Asteroids III, The University of Arizona Press, Tucson, Arizona, 395–408.
Brasser, R., Innanen, K.A., Connors, M., Veillet, C.,
Wiegert, P., Mikkola, S., and Chodas, P.W. (2004) Transient coorbital asteroids, Icarus, 171, 102–109.
Breiter, S. (1999) Explicit symplectic integrator for highlyTrends in 20th Century Celestial Mechanics
eccentric orbits, Celes. Mech. Dyn. Astron., 71, 229–241.
Bretagnon, P. (1974) Termes à longues périodes dans lesystème solaire, Astron. Astrophys., 30, 141–154.
Bretagnon, P. and Francou, G. (1992) General theory forthe outer planets, in Chaos, resonance and collective dynamical phenomena in the solar system, Kluwer Academic publishers, Dordrecht, 37–42.
Brouwer, D. (1959) Solution of the problem of artificial
satellite theory without drag, Astron. J., 64, 378–397.
Brouwer, D. and Clemence, G.M. (1961) Methods of Celestial Mechanics, Academic Press, New York.
Brouwer, D. and Hori, G. (1962) The motion of the Moon in
space, in Physics and Astronomy of the Moon, Academic
Press, New York, 1–26.
Brouwer, D. and van Woerkom, A.J.J. (1950) The secularvariations of the orbital elements of the principal planets,
Astron. Pap. Amer. Ephemeris. Naut. Alm., 13, 2, 81–
107.
Brown, E.W., Birkhoff, G.D., Leuschner, A.O., and Russell, H.N. (1922) Celestial Mechanics, Bulletin of the National Research Council, 4, 1–22.
Brown, M.E., Trujillo, C., and Rabinowitz, D. (2004) Discovery of a candidate inner Oort cloud planetoid,, Astrophys. J., 617, 645–649.
Brown, M.E., Trujillo, C.A., and Rabinowitz, D.L. (2005)
Discovery of a planetary-sized object in the scattered
Kuiper belt, Astrophys. J., 635, L97–L100.
Brumberg, V. ed. (1991) Essential Relativistic Celestial Mechanics, Adam Hilger, Bristol.
Brumberg, V.A. (1995) Analytical Techniques of Celestial
Mechanics, Springer-Verlag, New York.
Brumberg, E. and Brumberg, V.A. (1995) Derivative of
Kaula’s inclination function, Celes. Mech. Dyn. Astron.,62, 183–189.
Brumberg, V.A. and Brumberg, E. (2001) Elliptic anomaly
in constructing long-term and short-term dynamical theories, Celes. Mech. Dyn. Astron., 80, 159–166.
Brumberg, E., Brumberg, V.A., Konrad, T., and Soffel, M.
(1995) Analytical linear perturbation theory for highly
eccentric satellite orbits, Celes. Mech. Dyn. Astron., 61,
369–387.
Brunini, A. (2006) Origin of the obliquities of the giant
planets in mutual interactions in the early Solar System,
Nature, 440, 1163–1165.
Bruns, H. (1887) Über die Integrale des VielkörperProblems, Acta Math., XI, 25–96.
Burns, J.A. (2002) Two bodies are better than one, Science,
297, 942–943.
Burns, J.A., Lamy, P.L., and Soter, S. (1979) Radiation
forces on small particles in the solar system, Icarus, 40,
1–48.
Burrau, C. (1906) Uber einige in Aussicht genommene
Berechnun, betreffend einen spezialfall des Dreikorperproblems, Vierteljahrsschrift Astron. Ges., 41, 261.
Butler, R.P., Wright, J.T., Marcy, G.W., Fischer, D.A.,

97

Vogt, S.S., Tinney, C.G., Jones, H.R.A., Carter, B.D.,
Johnson, J.A., McCarthy, C., and Penny, A.J. (2006)
Catalog of nearby exoplanets, Astrophys. J., 646, 505–
522.
Canup, R.M. and Righter, K. eds. (2000) Origin of the
Earth and Moon, The University of Arizona Press, Tucson, Arizona.
Chambers, J.E. (1999) A hybrid symplectic integrator that
permits close encounters between massive bodies, Mon.
Not. R. Astron. Soc., 304, 793–799.
Chambers, J.E. (2001) Making more terrestrial planets,
Icarus, 152, 205–224.
Chambers, J.E. (2003) Symplectic integrators with complex time steps, Astron. J., 126, 1119–1126.
Chambers, J.E. and Murison, M.A. (2000) Pseudo-highorder symplectic integrators, Astron. J., 119, 425–433.
Chambers, J.E. and Wetherill, G.W. (1998) Making theterrestrial planets: N -body integrations of planetary embryos in three dimensions, Icarus, 136, 304–327.
Chambers, J.E., Wetherill, G.W., and Boss, A.P. (1996)
The stability of multi-planet systems, Icarus, 119, 261–
268.
Chambers, J.E., Quintana, E.V., Duncan, M.J., and Lissauer, J.J. (2002) Symplectic integrator algorithms formodeling planetary accretion in binary star systems, Astron. J., 123, 2884–2894.
Charlier, C.L. (1902) Die Mechanik des Himmels, Verlag
Von Veit & Comp., Leipzig.
Chazy, M.J. (1922) Sur l’allure de mouvement dans leproblème de trois corp quand le temps croit indéfiniment,
Ann. Sci. École Norm. Sup., 39, 29–130.
Chenciner, A. and Montgomery, R. (2000) A remarkable
periodic solution of the three-body problem in the caseof equal masses, Annals of Mathematics, 152, 881–901.
Chesley, S.R., Ostro, S.J., Vokrouhlický, D., Čapek, D.,
Girgini, J.D., Nolan, M.C., Margot, J.-L., Hine, A.A.,
Benner, L.A.M., and Chamberlin, A.B. (2003) Direct detection of the Yarkovsky effect by radar ranging to asteroid 6489 Golevka, Science, 302, 1739–1742.
Chirikov, B.V. (1979) A universal instability of manydimensional oscillators system, Physics Reports, 52, 263–
379.
Cohen, C.J. and Hubbard, E.C. (1965) Libration of theclose approaches of Pluto to Neptune, Astron. J., 70,
10–13.
Cohen, C.J., Hubbard, E.C., and Oesterwinter, C. (1973)
Planetary elements for 10000000 years, Celes. Mech., 7,
438–448.
Cohen, B.A., Swindle, T.D., and Kring, D.A. (2000) Support for the lunar cataclysm hypothesis from lunar meteorite impact melt ages, Science, 290, 1754–1756.
Connors, M., Chodas, P., Mikkola, S., Wiegert, P., Veillet, C., and Innanen, K. (2002) Discovery of an asteroid
and quasi-satellite in an Earth-like horseshoe orbit, Meteor. Planet Sci., 37, 1435–1441.

98

Takashi Ito and Kiyotaka Tanikawa

Correia, A.C.M. and Laskar, J. (2001) The four final rotation states of Venus, Nature, 411, 767–770.
Correia, A.C.M. and Laskar, J. (2003) Long-term evolutionof the spin of Venus II. Numerical simulations, Icarus,
163, 24–45.
Correia, A.C.M., Laskar, J., and de Surgy, O.N. (2003)
Long-term evolution of the spin of Venus I. Theory,
Icarus, 163, 1–23.
Daisaka, J.K., Tanaka, H., and Ida, S. (2006) Orbital evolution and accretion of protoplanets tidally interactingwith a gas disk II. Solid surface density evolution with
type-I migration, Icarus, 185, 492–507.
Denjoy, A. (1932) Sur les courbes définies par les équations
différetielles á la surface de tore, J. Math., XI, 333–375.
Deprit, A. (1969) Canonical transformations depending ona small parameter, Celes. Mech., 1, 12–30.
Dermott, S.F., Malhotra, R., and Murray, C.D. (1988) Dynamics of the Uranian and Saturnian systems: a chaotic
route to melting Miranda?, Icarus, 76, 295–334.
Devaney, R.L. (1980) Triple collision in the planar isosceles
three body problem, Inventiones Mathematicae, 60, 249–
267.
Dones, L., Levison, H., Duncan, M., and Weissman, P.
(1998) Formation of the Oort cloud revisited, Bulletin of
the American Astronomical Society, 30, 1113, American
Astronomical Society, DPS Meeting #30, #51.07.
Dones, L., Levison, H., Duncan, M., and Weissman, P.
(2000) Simulations of the formation of the Oort cloud ofcomets, Bulletin of the American Astronomical Society,32, 1060, American Astronomical Society, DPS Meeting#32, #36.02.
Dones, L., Weissman, P.R., Levison, H.F., and Duncan, M.J. (2004) Oort cloud formation and dynamics,
in Festou, M.C., Keller, H.U., and Weaver, H.A. eds.,
Comets II, The University of Arizona Press, Tucson, Arizona, 153–174.
Donnison, J.R. and Mikulskis, D.F. (1992) Three-body orbital stability criteria for circular orbits, Mon. Not. R.
Astron. Soc., 254, 21–26.
Dragt, A.J. and Finn, J.M. (1976) Lie series and invariant
functions for analytic symplectic maps, J. Math. Phys.,
17, 2215–2227.
Drovak, R., Pilat-Lohinger, E., Funk, B., and Freistetter, F.
(2003) Planets in habitable zones:. A study of the binary
Gamma Cephei, Astron. Astrophys., 398, L1–L4.
Duncan, M.J. and Levison, H.F. (1997) A scattered cometdisk and the origin of Jupiter family comets, Icarus, 276,
1670–1672.
Duncan, M.J. and Lissauer, J.J. (1997) Orbital stability of
the Uranian satellite system, Icarus, 125, 1–12.
Duncan, M.J. and Lissauer, J.J. (1998) The effects of postmain-sequence solar mass loss on the stability of our planetary system, Icarus, 134, 303–310.
Duncan, M.J., Quinn, T., and Tremaine, S. (1987) The
formation and extent of the solar system comet cloud,Astron. J., 94, 1330–1338.
Duncan, M., Quinn, T., and Tremaine, S. (1988) The originof short-periodic comets, Astrophys. J., 328, L69–L73.
Duncan, M.J., Levison, H.F., and Budd, S.M. (1995) The
dynamical structure of the Kuiper belt, Astron. J., 110,
3073–3081.
Duncan, M.J., Levison, H.F., and Lee, M.H. (1998) A multiple time step symplectic algorithm for integrating close
encounters, Astron. J., 116, 2067–2077.
Duriez, L. (1977) Théorie générale planétaire en variables
elliptiques 1. Développement des equations, Astron. Astrophys., 54, 93–112.
Duriez, L. (1982) Théorie générale planétaire etendu au casde la resonance et application au system des satellites
Galileens de Jupiter (General planetary theory extended
to the case of resonance and application to the Galilean
satellite system of Jupiter), Celes. Mech., 26, 231–255.
Ebisuzaki, T., Makino, J., Fukushige, T., Taiji, M., Sugimoto, D., Ito, T., and Okumura, S.K. (1993) GRAPE
project: An overview, Publ. Astron. Soc. Japan, 45, 269–
278.
Eckert, W.J., Brouwer, D., and Clemence, G.M. (1951) Coordinates of the five outer planets, Astron. Pap. Amer.
Ephemeris. Naut. Alm., 12.
Edgeworth, K.E. (1949) The origin and evolution of thesolar system, Mon. Not. R. Astron. Soc., 109, 600–609.
Encke, J.F. (1854) Berliner Astronomisches Jahrbuch füar,Königlich Akademie der Wissenschaften, Berlin.
Evans, N.W. and Tremaine, S. (1999) Linear multistep
methods for integrating reversible differential equations,Astron. J., 118, 1888–1899.
Feigenbaum, M. (1978) Quantitative universality for a classof nonlinear transformations, J. Stat. Phys., 19, 25–52.
Fernández, J.A. (1997) The formation of the Oort cloud andthe primitive galactic environment, Icarus, 129, 106–119.
Fernández, J.A. and Ip, W.-H. (1981) Dynamical evolution of a cometary swarm in the outer planetary region,
Icarus, 47, 470–479.
Fernández, J.A. and Ip, W.-H. (1984) Some dynamical aspects of the accretion of Uranus and Neptune: the exchange of orbital angular momentum with planetesimals,
Icarus, 58, 109–120.
Ferraz-Mello, S. (1994) Dynamics of the asteroidal 2/1 resonance, Astron. J., 108, 2330–2337.
Flaschka, H. (1974) Toda lattice II: existence of integrals,
Phys. Rev. B, 9, 1924–1925.
Fleming, H.J. and Hamilton, D.P. (2000) On the origin of
the Trojan asteroids: Effects of Jupiter’s mass accretion
and radial migration, Icarus, 148, 479–493.
Ford, E,B., Kozinsky, B., and Rasio, F.A. (2000) Secular
evolution of hierarchical triple star systems, Astrophys.
J., 535, 385–401.
Froeschlé, Ch. and Morbidelli, A. (1994) The secular resonances in the solar system, in Milani, A., di Martino, M.,
and Cellino, A. eds., Proceeding of the 160th IAU Sym-

Trends in 20th Century Celestial Mechanics
posium, Asteroids, Comets, Meteors 1993,, Kluwer Academic Publishers, Dordrecht, 189–204.
Froeschlé, Ch. and Scholl, H. (1989) The three principal
secular resonances ν5 , ν6 , and ν16 in the asteroidal belt,
Celes. Mech. Dyn. Astron., 46, 231–251.
Froeschlé, Cl., Gonczi, R., and Lega, E. (1997) The fast
Lyapunov indicator: A simple tool to detect weak chaos.
application to the structure of the main asteroidal belt,
Planet. Space Sci., 45, 881–886.
Fujiwara, T. and Hirai, M. (2006) Long-term variabilityof late-type stars: Analysis of photometric records, inSôma, M. and Tanikawa, K. eds., Proceedings of the
Symposium on Historical records and Modern Science,
National Astronomical Observatory of Japan, Tokyo,
Japan, 7–19, (held in Mitaka, Tokyo, June 2–3, 2006).
Fujiwara, T. and Montgomery, R. (2003) Convexity inthe figure eight solution to the three-body problem,
http://arxiv.org/abs/math.DS/0308252/.
Fujiwara, T., Fukuda, H., and Ozaki, H. (2003a) Choreographic three bodies on the lemniscate, J. Phys. A: Math.
Gen., 36, 1–10.
Fujiwara, T., Fukuda, H., and Ozaki, H. (2003b) Evolution of the moment of inertia of three-body figure-eightchoreography, J. Phys. A: Math. Gen., 36, 10537–10549.
Fujiwara, T., Fukuda, H., Ozaki, H., and Yamada, M.
(2004a) Synchronized similar triangles for three-body orbits with zero angular momentum, J. Phys. A: Math.
Gen., 37, 10571–10584.
Fujiwara, T., Yamaoka, H., and Miyoshi, S. (2004b) Survey of long-term variability of stars. I. Reliability of magnitudes in old star catalogues, Astron. Astrophys., 416,
641–646.
Fukushige, T., Makino, J., and Kawai, A. (2005) GRAPE–
6A: A single-card GRAPE–6 for parallel PC-GRAPE
cluster systems, Publ. Astron. Soc. Japan, 57, 1009–1021.
Fukushima, T. (1996a) A fast procedure solving Kepler’s
equation for elliptic case, Astron. J., 112, 2858–2861.
Fukushima, T. (1996b) Generalization of Encke’s method
and its application to the orbital and rotational motionsof celestial bodies, Astron. J., 112, 1263–1277.
Fukushima, T. (1997a) A method solving Kepler’s equation
for hyperbolic case, Celes. Mech. Dyn. Astron., 68, 121–
137.
Fukushima, T. (1997b) A method solving Kepler’s equation without transcendental function evaluations, Celes.
Mech. Dyn. Astron., 66, 309–319.
Fukushima, T. (1997c) Picard iteration method, Chebyshev polynomial approximation, and global numerical integration of dynamical motions, Astron. J., 113, 1909–
1914.
Fukushima, T. (1997d) A procedure solving the extended
Kepler’s equation for the hyperbolic case, Astron. J.,
113, 1920–1924.
Fukushima, T. (1997e) Vector integration of dynamical motions by the Picard–Chebyshev method, Astron. J., 113,99

2325–2328.
Fukushima, T. (1998) A fast procedure solving Gauss’ formof Kepler’s equation, Celes. Mech. Dyn. Astron., 70, 115–
130.
Fukushima, T. (1999a) Parallel/vector integration methods
for dynamical astronomy, Celes. Mech. Dyn. Astron., 73,
231–241.
Fukushima, T. (1999b) Symmetric multistep methods revisited: II. Numerical experiments, in Svoren, J. and Pittich, E.M. eds., Proc. 173rd Colloq. IAU., Evolution andsource regions of asteroids and comets, Astronomical Institute of the Slovak Academy of Sciences, 309, Tatranska Lomnica, Slovak Republic, August 24–28, 1998.
Fukushima, T. (2001) Reduction of round-off error in symplectic integrators, Astron. J., 121, 1768–1775.
Fukushima, T. (2003a) Efficient orbit integration by dual
scaling for consistency of Kepler energy and Laplace integral, Astron. J., 126, 2567–2573.
Fukushima, T. (2003b) Efficient orbit integration by scaling
and rotation for consistency of Kepler energy, Laplace
integral, and angular momentum direction, Astron. J.,
126, 3138–3142.
Fukushima, T. (2003c) Efficient orbit integration by scalingfor Kepler energy consistency, Astron. J., 126, 1097–
1111.
Fukushima, T. (2004a) Efficient integration of highly eccentric orbits by quadruple scaling for Kustaanheimo–
Stiefel regularization, Astron. J., 128, 3108–3113.
Fukushima, T. (2004b) Efficient integration of highly eccentric orbits by scaling methods applied to Kustaanheimo–
Stiefel regularization, Astron. J., 128, 3114–3122.
Fukushima, T. (2004c) Efficient orbit integration by linear transformation for consistency of Kepler energy, full
Laplace integral, and angular momentum vector, Astron.
J., 127, 3638–3641.
Fukushima, T. (2004d) Efficient orbit integration by the orbital longitude method using antifocal anomaly, Astron.
J., 128, 1455–1464.
Fukushima, T. (2004e) Further simplification of the manifold correction method for orbit integration, Astron. J.,
128, 1446–1454.
Fukushima, T. (2004f) Simplification of a manifold correction method for orbit integration, Astron. J., 128, 920–
922.
Fukushima, T. (2005a) Efficient orbit integration by
Kustaanheimo–Stiefel regularization using time element,
Astron. J., 129, 1746–1754.
Fukushima, T. (2005b) Efficient orbit integration by linear
transformation for Kustaanheimo–Stiefel regularization,
Astron. J., 129, 2496–2503.
Fukushima, T. (2005c) Efficient orbit integration by orbital
longitude methods with Sundman transformation on thetime variable, Astron. J., 129, 1171–1177.
Fukushima, T. (2005d) Efficient orbit integration by scaling
and rotation for consistency of Kepler energy, Laplace

100

Takashi Ito and Kiyotaka Tanikawa

integral, and angular momentum direction, Celes. Mech.
Dyn. Astron., 2496–2503.
Fukushima, T. (2007) New two-body regularization, Astron. J., 133, 1–10.
Funato, Y., Makino, J., Hut, P., Kokubo, E., and Kinoshita, D. (2004) The formation of Kuiper-belt binaries
through exchange reactions, Nature, 427, 518–520.
Garcı́a-Sánchez, J., Preston, R.A., Jones, D.L., Weissman, P.R., Lestrade, J-.F., Latham, D.W., and Stefanik, R.P. (1999) Stellar encounters with the Oort cloud
based on Hipparcos data, Astron. J., 1042–1055.
Garcı́a-Sánchez, J., Weissman, P.R., Preston, R.A.,
Jones, D.L., Lestrade, J-.F., Latham, D.W., Stefanik, R.P., and Paredes, J.M. (2001) Stellar encounterswith the solar system, Astron. Astrophys., 634–659.
Gardner, C.S., Greene, J.M., Kruskal, M.D., and Miura, R.
(1967) Method for solving the Korteweg–deVries equation, Phys. Rev. Lett., 19, 1095–1097.
Garfinkel, B. (1959) The orbit of a satellite of an oblate
planet, Astron. J., 64, 353–367.
Getino, J., Ferrándiz, J.M., and Escapa, A. (2001) Hamiltonian theory for the non-rigid Earth: Semidiurnal terms,
Astron. Astrophys., 370, 330–341.
Giacaglia, G.E., Murphy, J.P., and Felsentreger, T.L.
(1970) A semi-analytic theory for the motion of a lunar
satellite, Celes. Mech., 3, 3–66.
Giorgilli, A. and Skokos, C. (1997) On the stability of the
Trojan asteroids, Astron. Astrophys., 317, 254–261.
Gladman, B. (1993) Dynamics of systems of two close planets, Icarus, 106, 247–263.
Gladman, B. (2005) The Kuiper Belt and the solar system’s
comet disk, Science, 307, 71–75.
Gladman, B., Duncan, M., and Candy, J. (1991) Symplectic integrators for long-term integrations in celestial mechanics, Celes. Mech. Dyn. Astron., 52, 221–240.
Goldreich, P. (1966) History of the lunar orbit, Rev. Geophys., 4, 411–439.
Goldreich, P. and Peale, S.J. (1966) Spin-orbit coupling inthe solar system, Astron. J., 71, 425–438.
Goldreich, P., Lithwick, Y., and Sari, R. (2002) Formationof Kuiper-belt binaries by dynamical friction and threebody encounters, Nature, 420, 643–646.
Gomes, R., Levison, H.F., Tsiganis, K., and Morbidelli, A.
(2005) Origin of the cataclysmic late heavy bombardment period of the terrestrial planets, Nature, 435, 466–
469.
Goriely, A. (2000) A brief history of Kovalevskaya exponents and modern developments, Regular and Chaotic
Dynamics, 5, 3–15.
Greenberg, R. (1977) Orbit-orbit resonances in the solar
system: varieties and similarities, Vistas Astron., 21,
209–239.
Greene, J.M. (1979) A method of determining a stochastictransition, J. Math. Phys., 20, 1183–1201.
Gronchi, G.F. (2002) Generalized averaging principle andthe secular evolution of planet crossing orbits, Celes.
Mech. Dyn. Astron., 83, 97–120.
Gronchi, G.F. and Michel, P. (2001) Secular orbital evolution, proper elements, and proper frequencies for nearEarth asteroids: A comparison between semianalytic theory and numerical integrations, Icarus, 152, 48–57.
Gronchi, G.F. and Milani, A. (1999a) Averaging on Earthcrossing orbits, Celes. Mech. Dyn. Astron., 71, 109–136.
Gronchi, G.F. and Milani, A. (1999b) The stable Kozai
state for asteroids and comets. With arbitrary semimajor
axis and inclination, Astron. Astrophys., 341, 928–935.
Gronchi, G.F. and Milani, A. (2001) Proper elements for
Earth-crossing asteroids, Icarus, 152, 58–69.
Hadamard, J. (1898) Les surfaces á courbure oposée etleurlignes géodésiques, Journal de Mathématiques pures
et appliquées, 4, 27–73.
Hagihara, Y. (1970) Celestial Mechanics I. Dynamical Principles and Transformation Theory, MIT Press, Cambridge, Massachusetts.
Hagihara, Y. (1971) Recent advances of celestial mechanics
in the Soviet Union, Vistas in Astronomy, 13, 15–47.
Hagihara, Y. (1972a) Celestial Mechanics II. Perturbation
Theory, part 1, MIT Press, Cambridge, Massachusetts.
Hagihara, Y. (1972b) Celestial Mechanics II. Perturbation
Theory, part 2, MIT Press, Cambridge, Massachusetts.
Hairer, E., Nørsett, S.P., and Wanner, G. (1993) Solving Ordinary Differential Equation I. Nonstiff Problems,
Springer-Verlag, New York, second revised edition.
Hairer, E., Lubich, C., and Wanner, G. (2002) Geometric
Numerical Integration: Structure-Preserving Algorithms
for Ordinary Differential Equations, No. 31 in Springer
Series in Computational Mathematics, Springer-Verlag,
Berlin, Heidelberg.
Hall, G.R. (1984) A topological version of a theorem on
twist maps, Ergod. Th. & Dynam. Sys., 4, 585–603.
Hamada, T., Fukushige, T., Kawai, A., and Makino, J.
(2000) PROGRAPE–1: A programmable, multi-purpose
computer for many-body simulations, Publ. Astron. Soc.
Japan, 52, 943–954.
Han, Y.-B. (1997) The secular variations of the Earth rotation derived from ancient astronomical observations,Astrophys Rep., Publ. Beijing Astron. Obs., 31, 66–70.
Harrington, R.S. (1975) Production of triple stars by the
dynamical decay of stellar systems, Astron. J., 80, 1081–
1086.
Harrington, R.S. (1977) Planetary orbits in binary stars,
Astron. J., 82, 753–756.
Hartmann, W.K., Phillips, R.J., and Taylor, G.J. eds.
(1986) Origin of the Moon, Lunar and Planetary Institute, Houston.
Hartmann, W.K., Ryder, G., Dones, L., and Grinspoon, D.
(2000) The time-dependent intense bombardment of
the primordial Earth/Moon system, in Canup, R. andRighter, K. eds., Origin of the Earth and Moon, The University of Arizona Press, Tucson, Arizona, 493–512.

Trends in 20th Century Celestial Mechanics
Hayashi, C., Nakazawa, K., and Nakagawa, Y. (1985) Formation of the solar system, in Protostars & Planets II,
The University of Arizona Press, Tucson, Arizona, 1100–
1153.
Hayli, A. (1970) More comparisons of numerical integrations in the gravitational N-body problem, Astronomy
and Astrophysics, 7, 249–253, text in French.
Head, J.W., Mustard, J.F., Kreslavsky, M.A., Milliken, R.E., and Marchant, D.R. (2003) Recent ice ages
on Mars, Nature, 426, 797–802.
Heggie, D.C. (1974) A global regularisation of the gravitational N–body problem, Cel. Mech., 10, 217–241.
Heggie, D.C. (1975) Binary evolution in stellar dynamics,
Mon. Not. R. Astron. Soc., 173, 729–787.
Heggie, D.C. and Hut, P. (1993) Binary-single star scattering IV. Analytic approximation and fitting formulae forcross sections and reaction rates, Astrophys. J. Suppl.,85, 347–409.
Heisler, J. and Tremaine, S. (1986) The influence of thegalactic tidal field on the Oort comet cloud, Icarus, 65,
13–26.
Heisler, J., Tremaine, S., and Alcock, C. (1987) The frequency and intensity of comet showers from the Oort
cloud, Icarus, 70, 269–288.
Hénon, M. (1965) Exploration numerique du probleme restreint I., Annals Astrophysique, 28, 499–511.
Hénon, M. (1974) Integrals of the Toda lattice, Phys. Rev.
B, 9, 1921–1923.
Hénon, M. (1976) A two-dimensional mapping with astrange attractor, Commun. Math. Phys., 50, 69–77.
Hénon, M. and Heiles, C. (1964) The applicability of thethird integral of motion: Some numerical experiments,
Astron. J., 69, 73.
Hénon, M. and Petit, J.-M. (1986) Series expansion forencounter-type solutions of Hill’s problem, Celes. Mech.,38, 67–100.
Herman, M.R. (1983) Sur les courbes invariantes par lesdifféomorphisms de l’anneau vol. 1, Astérisque, 103–
104.
Herman, M.R. (1986) Sur les courbes invariantes par lesdifféomorphisms de l’anneau vol. 2, Astérisque, 144.
Higuchi, A., Kokubo, E., and Mukai, T. (2006) Scatteringof planetesimals by a planet: Formation of comet cloudcandidates, Astron. J., 131, 1119–1129.
Hill, G.W. (1897) On the values of the eccentricities andlongitudes of the perihelia of Jupiter and Saturn for distant epochs, Astron. J., 17, 81–87.
Hinkle, K.H., Fekel, F.C., Johnson, D.S., and Scharlach, W.W.G. (1993) The triple symbiotic system CH
Cygni, Astron. J., 105, 1074–1086.
Hirayama, S. (1918) Groups of asteroids probably of common origin, Astron. J., 31, 185–188.
Hirota, R. and Satsuma, J. (1976) A variety of nonlinear
network equations generated from the Bäcklund transformation for the Toda lattice, Prog. Theor. Phys. suppl.,

101

59, 64–100.
Hirsch, M.W. and Smale, S. (1974) Differential Equations, Dynamical Systems, and Linear Algebra, Academic
Press, New York.
Hogg, D.W., Quinlan, G.D., and Tremaine, S. (1991) Dynamical limits on dark mass in the outer solar system,
Astron. J., 101, 2274–2286.
Holman, M.J. and Murray, N.W. (1996) Chaos in highorder mean motion resonances in the outer asteroid belt,
Astron. J., 112, 1278–1293.
Hori, G. (1966) Theory of general perturbations with unspecified canonical variables, Publ. Astron. Soc. Japan,
18, 287–296.
Hori, G. (1967) Non-linear coupling of two harmonic oscillations, Publ. Astron. Soc. Japan, 19, 229–241.
Hori, G. (1970) Comparison of two perturbation theories
based on the canonical transformation, Publ. Astron.
Soc. Japan, 22, 191–198.
Hori, G. (1971) Theory of general perturbations for noncanonical systems, Publ. Astron. Soc. Japan, 23, 567–
587.
Ida, S., Canup, R.M., and Stewart, G.R. (1997) Lunar accretion from an impact-generated disk, Nature, 389, 353–
357.
Ida, S., Larwood, J.D., and Burkert, A. (2000) Evidence
for early stellar encounters in the orbital distributionof Edgeworth–Kuiper Belt objects, Astrophys. J., 528,
351–356.
Innanen, K.A., Zheng, J.Q., Mikkola, S., and Valtonen, M.J. (1997) The Kozai mechanism and the stability of planetary orbits in binary star systems, Astron. J.,
113, 1915–1919.
Ito, T. and Fukushima, T. (1997) Parallelized extrapolation method and its application to the orbital dynamics,
Astron. J., 114, 1260–1267.
Ito, T. and Kojima, S. (2005) Inaccuracies of trigonometricfunctions in computer math libraries, Publ. Natl. Astron.
Obs. Japan, 8, 17–31.
Ito, T. and Malhotra, R. (2006) Dynamical transport ofasteroid fragments from the ν6 resonance, Adv. SpaceRes., 38, 817–825, (astro-ph/0611548).
Ito, T. and Tanikawa, K. (1999) Stability and instabilityof the terrestrial protoplanet system and their possibleroles in the final stage of planet formation, Icarus, 139,
336–349.
Ito, T. and Tanikawa, K. (2001) Stability of terrestrial protoplanet systems and alignment of orbital elements, Publ.
Astron. Soc. Japan, 53, 143–151.
Ito, T. and Tanikawa, K. (2002) Long-term integrations andstability of planetary orbits in our solar system, Mon.
Not. R. Astron. Soc., 336, 483–500.
Ito, T., Makino, J., Ebisuzaki, T., and Sugimoto, D. (1990)
A special-purpose N -body machine GRAPE–1, Computer Physics Communications, 60, 187–194.
Ito, T., Ebisuzaki, T., Makino, J., and Sugimoto, D. (1991)

102

Takashi Ito and Kiyotaka Tanikawa

A special-purpose computer for gravitational many-bodysystems: GRAPE–2, Publ. Astron. Soc. Japan, 43, 547–
555.
Ito, T., Kumazawa, M., Hamano, Y., Matsui, T., and Masuda, K. (1993) Long term evolution of the solar insolation variation over 4 Ga, Proc. Japan Acad., Ser. B, 69,
233–237.
Iwasawa, M., Funato, Y., and Makino, J. (2006) Evolutionof massive blackhole triples. I. Equal-mass binary-single
systems, Astrophys. J., 651, 1059–1067.
Jeffreys, H. and Vicente, R.O. (1957) The theory of nutation and the variation of latitude, Mon. Not. R. Astron.
Soc., 117, 162–173.
Jewitt, D. (1999) Kuiper belt objects, Annu. Rev. Earth
Planet. Sci., 27, 287–312.
Jewitt, D.C. and Luu, J.X. (1993) Discovery of the candidate Kuiper belt object 1992 QB1, Nature, 362, 730–732.
Kaplan, M.L., van Cleve, J.E., and Alcock, C. (2003) Minor body surveyor: A multi-object, high speed, spectrophotometer space mission system employing widearea intelligent change detection, American GeophysicalUnion, Fall Meeting 2003, abstract# P41B–0403.
Katok, A. (1982) Some remarks on Birkhoff and Mather
twist map theorems, Ergod. Th, & Dynam. Sys., 2, 185–
194.
Katok, A. and Hasselblatt, B. (1995) Introduction to theModern Theory of Dynamical Systems, No. 54 in Encyclopedia of Mathematics and Its Applications, Cambridge
University Press, Cambridge.
Kawabata, K., Tanikawa, K., and Sôma, M. (2002) Reliability of astronomical records in the Nihongi, Rep. Natl.
Astron. Obs. Japan, 5, 145–159.
Kawabata, K., Sôma, M., and Tanikawa, K. (2004) A solar
eclipse verified in Shijing, Publ. Astron. Soc. Japan, 869–
877.
Kawai, A. and Fukushige, T. (2006) $158/GFlop astrophysical N -body simulation with a reconfigurable addin card and a hierarchical tree algorithm, International Conference on High Performance Computing, Networking and Storage, November 11–17, 2006, Tampa,
Florida. http://sc06.supercomputing.org/schedule/event
detail.php?evid=5054.
Kawai, A., Fukushige, T., Makino, J., and Taiji, M. (2000)
GRAPE–5: A special-purpose computer for N -body simulations, Publ. Astron. Soc. Japan, 52, 659–676.
Khilmi, G.F. (1961) Qualitative Methods in the Many-Body
Problem, Gordon Breach, New York.
Kimura, H. (1902) A new annual term in the variation oflatitude, independent of the components of the pole’s
motion, Astron. J., XXII, 517, 107–108.
Kinoshita, H. (1977a) Theory of the rotation of the rigid
earth, Celes. Mech., 15, 277–326.
Kinoshita, H. (1977b) Third-Order Solution of an Artificial
Satellite Theory, No. 379 in SAO Special Report, Smithsonian Astronomical Observatory, pp. 104.

Kinoshita, H. and Nakai, H. (1984) Motions of the perihelions of Neptune and Pluto, Celes. Mech., 34, 203–217.
Kinoshita, H. and Nakai, H. (1995) The motion of Pluto
over the age of the solar system, in Ferraz-Mello, S. et al .
ed., Dynamics, ephemerides and astrometry in the solar
system, Kluwer Academic publishers, Dordrecht, 61–70.
Kinoshita, H. and Nakai, H. (1996) Long-term behavior of
the motion of Pluto over 5.5 billion years, Earth, Moon,
and Planets, 72, 165–173.
Kinoshita, H. and Souchay, J. (1990) The theory of thenutation of the rigid earth: Model at the second order,
Celes. Mech. Dyn. Astron., 48, 187–265.
Kinoshita, H., Yoshida, H., and Nakai, H. (1991) Symplectic integrators and their application to dynamical astronomy, Celes. Mech. Dyn. Astron., 50, 59–71.
Kirkwood, D. (1867) Meteoric Astronomy, Lippincott,Philadelphia.
Kiseleva, L.G., Eggleton, P.P., and Anosova, J.P. (1994a)
A note on the stability of hierarchical triple stars withinitially circular orbits, Mon. Not. R. Astron. Soc., 267,
161–166.
Kiseleva, L.G., Eggleton, P.P., and Orlov, V.V. (1994b)
Instability of close triple systems with coplanar initial
doubly circular motion, Mon. Not. R. Astron. Soc., 270,
936–946.
Knez̆ević, Z. and Milani, A. (1994) Asteroid proper elements: The big picture, in Milani, A., di Martino, M.,
and Cellino, A. eds., Asteroids, Comets, Meteors 1993,
Kluwer Academic Publishers, Dordrecht, 143–158.
Knez̆ević, Z., Froeschlé, Ch., Lemaı̂tre, A., Milani, A., andMorbidelli, A. (1995) Comparison between two theories
of asteroid proper elements, Astron. Astrophys., 293,
605–612.
Knez̆ević, Z., Lemaı̂tre, A., and Milani, A. (2002) The determination of asteroid proper elements, in Bottke, W.F.,
Cellino, A., Paolicchi, P., and Binzel, R.P. eds., Asteroids
III, The University of Arizona Press, Tucson, Arizona,603–612.
Kobayashi, H. and Ida, S. (2001) The effects of a stellar
encounter on a planetesimal disk, Icarus, 153, 416–429.
Kokubo, E. and Ida, S. (1995) Orbital evolution of protoplanets embedded in a swarm of planetesimals, Icarus,
114, 247–257.
Kokubo, E. and Ida, S. (1996) On runaway growth of planetesimals, Icarus, 123, 180–191.
Kokubo, E. and Ida, S. (1998) Oligarchic growth of protoplanets, Icarus, 131, 171–178.
Kokubo, E. and Ida, S. (2000) Formation of protoplanetsfrom planetesimals in the solar nebula, Icarus, 143, 15–27.
Kokubo, E. and Ida, S. (2002) Formation of protoplanet
systems and diversity of planetary systems, Astrophys.
J., 581, 666–680.
Kokubo, E. and Makino, J. (2004) A modified Hermite integrator for planetary dynamics, Publ. Astron. Soc. Japan,

Trends in 20th Century Celestial Mechanics
56, 861–868.
Kokubo, E., Yoshinaga, K., and Makino, J. (1998) On
a time-symmetric Hermite integrator for planetary N body simulation, Mon. Not. R. Astron. Soc., 297, 1067–
1072.
Kokubo, E., Ida, S., and Makino, J. (2000) Evolution ofa circumterrestrial disk and formation of a single moon,
Icarus, 148, 419–436.
Kolmogorov, A.N. (1954) Preservation of conditionally periodic movements with small change in the Hamiltonian
function, Dokl. Akad. Nauk. USSR, 98, 527–531.
Kominami, J. and Ida, S. (2002) The effect of tidal interaction with a gas disk on Formation of terrestrial planets,
Icarus, 157, 43–56.
Kominami, J. and Ida, S. (2004) Formation of terrestrial
planets in a dissipating gas disk with Jupiter and Saturn,Icarus, 167, 231–243.
Kominami, J., Tanaka, H., and Ida, S. (2005) Orbital evolution and accretion of protoplanets tidally interactingwith a gas disk: I. Effects of interaction with planetesimals and other protoplanets, Icarus, 178, 540–552.
Kortenkamp, S.J., Malhotra, R., and Michtchenko, T.
(2004) Survival of Trojan-type companions of Neptune
during primordial planet migration, Icarus, 167, 347–
359.
Korycansky, D.G., Laughlin, G., and Adams, F.C. (2001)
Astronomical engineering: A strategy for modifying planetary orbits, Astrophys. Space Sci., 275, 349–366.
Kozai, Y. (1959) The motion of a close Earth satellite, Astron. J., 64, 367–377.
Kozai, Y. (1962a) Second-order solution of artificial satellite theory without air drag, Astron. J., 67, 446–461.
Kozai, Y. (1962b) Secular perturbations of asteroids withhigh inclination and eccentricity, Astron. J., 67, 591–598.
Kozlov, V.V. (1983) Integrability and non-integrability in
Hamiltonian mechanics, Russ. Math. Survey, 38, 1–76.
Kring, D.A. and Cohen, B.A. (2002) Cataclysmic bombardment throughout the inner solar system 3.9–4.0 Ga, J.
Geophys. Res., 107, E2, 4–1–4–6.
Kuchner, M,J., Brown, M.E., and Holman, M. (2002) Longterm dynamics and the orbital inclinations of the classicalKuiper belt objects, Astron. J., 124, 1221–1230.
Kuiper, G.P. (1951) On the origin of the solar system, inHynek, J.A. ed., Astrophysics, McGraw Hill, New York,357–424.
Kustaanheimo, P. and Stiefel, E. (1965) Perturbation theory of Keplerian motion based on Spinor regularization,
J. Reine Angew. Math., 218, 204–219.
Kuwabara, K. and Tanikawa, K. (2006) An extension of
the theorems of three-tangents and three-normals, Phys.
Lett. A, 354, 445–448.
Lagrange, J.L. (1788) Mécanique Analitique, la Veuve Desaint, Paris.
Lambert, J.D. and Watson, I.A. (1976) Symmetric multistep methods for periodic initial value problems, J. Inst.

103

Maths. Applics., 18, 189–202.
Landau, L.D. and Lifshitz, E.M. eds. (1959) Fluid Mechanics, Pergamon Press, Oxford, translated from the Russian by J. B. Sykes and W. H. Reid.
Laplace, P.-S. (1814) Essai Philosophique sur les Probabilités, Courcier Imprimeur, Paris, reprinted in Oeuvres
complètes de Laplace (vol. 8), Gauthier Villars, Paris.
Laskar, J. (1985) Accurate methods in general planetary
theory, Astron. Astrophys., 144, 133–146.
Laskar, J. (1986) Secular terms of classical planetary theories using the results of general theory, Astron. Astrophys., 157, 59–70.
Laskar, J. (1988) Secular evolution of the solar system over
10 million years, Astron. Astrophys., 198, 341–362.
Laskar, J. (1989) A numerical experiment on the chaotic
behaviour of the solar system, Nature, 338, 237–238.
Laskar, J. (1990) The chaotic motion of the solar system:
A numerical estimate of the size of the chaotic zones,
Icarus, 88, 266–291.
Laskar, J. (1993) Frequency analysis of a dynamical system,
Celes. Mech. Dyn. Astron., 56, 191–196.
Laskar, J. (1994) Large scale chaos in the solar system,
Astron. Astrophys., 287, L9–L12.
Laskar, J. (1996) Large scale chaos and marginal stability
in the solar system, Celes. Mech. Dyn. Astron., 64, 115–
162.
Laskar, J. (1997) Large scale chaos and the spacing of theinner planets, Astron. Astrophys., 317, L75–L78.
Laskar, J. and Robutel, P. (1993) The chaotic obliquity ofthe planets, Nature, 361, 608–612.
Laskar, J. and Robutel, P. (2001) High order symplectic
integrators for perturbed Hamiltonian systems, Celes.
Mech. Dyn. Astron., 80, 39–62.
Laskar, J., Froeschlé, C., and Celletti, A. (1992) The measure of chaos by the numerical analysis of the fundamental frequencies. Application to the standard mapping,
Physica D, 56, 253–269.
Laskar, J., Joutel, F., and Boudin, F. (1993a) Orbital, precessional, and insolation quantities for the earth from
−20Myr to +10Myr, Astron. Astrophys., 270, 522–533.
Laskar, J., Joutel, F., and Robutel, P. (1993b) Stabilizationof the Earth’s obliquity by the Moon, Nature, 361, 615–
617.
Laskar, J., Levrard, B., and Mustard, J. (2002) Orbital
forcing of the martian polar layered deposits, Nature,
419, 375–377.
Laughlin, G. and Adams, F.C. (1999) Stability and chaosin the υ Andromedae planetary system, Astrophys. J.,
526, 881–889.
Lazutkin, V.F., Schachmannsky, I.G., and Tabanov, M.B.
(1989) Splitting of separatrices for standard and semistandard mappings, Physica D, 40, 235–248.
Lazzaro, D., Ferraz-Mello, S., and Fernández, J.A. (2006)
Asteroids, Comets, and Meteors, Cambridge University
Press, Cambridge, England, Proceedings of the Interna-

104

Takashi Ito and Kiyotaka Tanikawa

tional Astronomical Union Symposium 229, 7–12 August
2005, Búzios, Rio de Janeiro, Brazil.
Le Calvez, P. (1987) Propriétés dynamiques des régionsd’instabilité, Ann. Scient. Éc. Norm. Sup., 20, 443–464.
Le Calvez, P. (2000) Dynamical properties of diffeomorphisms of the annulus and of the torus, Vol. 4 of
SMF/AMS texts and Monographs, Amer. Math. Soc.,
Province, Rhode Island.
Lecar, M. (1968) A comparison of eleven numerical integrations of the same gravitational 25-body problem, Bulletin Astronomique, 3, 91–104.
Lecar, M. and Franklin, F. (1992) On the original distribution of the asteroids IV. Numerical experiments in theouter asteroid belt, Icarus, 96, 234–250.
Lecar, M., Franklin, F., and Murison, M. (1992) On predicting long-term orbital instability: A relation between
the Lyapunov time and sudden orbital transitions, Astron. J., 104, 1230–1236.
Lecar, M., Franklin, F.A., and Holman, M. (2001) Chaos
in the solar system, Annu. Rev. Astron. Astrophys., 39,
581–631.
Lee, D.-C., Halliday, A.N., Snyder, G.A., and Taylor, L.A.
(1997a) Age and origin of the moon, Science, 278, 1098–
1103.
Lee, M.H., Duncan, M.J., and Levison, H.F. (1997b) Variable time step integrators for long-term orbital integrations, Astron. Soc. Pac. Conf. Ser., 123, 32–37.
Lemaı̂tre, G. (1952) Cordonnées symétriques dans leprobléme des trois corps, Bull. Classe. Sci. Acad. R.
Belg., 38, 582–592.
Levi-Civita, T. (1906) Sur la résolution qualitative duprobléme des trois corps, Acta Math., 30, 305–327.
Levison, H.F. and Agnor, C. (2003) The role of giant planetsin terrestrial planet formation, Astron. J., 125, 2692–
2713.
Levison, H.F. and Duncan, M.J. (1994) The long-term dynamical behavior of short-period comets, Icarus, 108,
18–36.
Levison, H.F. and Duncan, M.J. (1997) From the Kuiper
belt to Jupiter-family comets: The spatial distribution of
ecliptic comets, Icarus, 127, 13–32.
Levison, H.F. and Duncan, M.J. (2000) Symplectically integrating close encounters with the Sun, Astron. J., 120,
2117–2123.
Levison, H.F. and Morbidelli, A. (2003) The formation of
the Kuiper belt by the outward transport of bodies during Neptune’s migration, Nature, 426, 419–421.
Levison, H.F., Lissauer, J.J., and Duncan, M.J. (1998)
Modeling the diversity of outer planetary systems, Astron. J., 116, 1998–2014.
Levison, H., Dones, L., Duncan, M., and Weissman, P.
(1999) The formation of the Oort cloud, Bulletin of theAmerican Astronomical Society, 31, 1079, American Astronomical Society, DPS Meeting #31, #06.05.
Levison, H.F., Done, L., Chapman, C.R., Stern, S.A., Dun-can, M.J., and Zanhle, K. (2001a) Could the lunar “lateheavy bombardment” have been triggered by the formation of Uranus and Neptune?, Icarus, 151, 286–306.
Levison, H.F., Dones, L., and Duncan, M.J. (2001b) The
origin of Halley-type comets: probing the inner Oort
cloud, Astron. J., 121, 2253–2267.
Levison, H,F., Morbidelli, A., Dones, L., Jedicke, R.,
Wiegert, P.A., and Bottke, W.F. (2002) The mass disruption of Oort cloud comets, Science, 296, 2212–2215.
Levison, H.F., Duncan, M.J., Dones, J., and Gladman, B.J.
(2006) The scattered disk as a source of Halley-type
comets, Icarus, 184, 619–633.
Li, T.Y. and Yorke, J.A. (1975) Period three implies chaos,Am. Math. Monthly, 82, 985–992.
Liao, X. (1997) Symplectic integrator for general nearintegrable Hamiltonian system, Celes. Mech. Dyn. Astron., 66, 243–253.
Lie, S. (1888) Theorie der Transformationgruppen I, Teubner, Lipzig.
Liou, J.C., Dermott, S.F., and Xu, Y.L. (1995) The contribution of cometary dust to the zodiacal cloud, Planet.
Space Sci., 43, 717–722.
Lissauer, J.J. (1999) Chaotic motion in the solar system,
Rev. Modern Phys., 71, 835–845.
Lissauer, J.J., Quintana, E.V., Rivera, E.J., and Duncan, M.J. (2001) The effect of a planet in the asteroid belton the orbital stability of the terrestrial planets, Icarus,
154, 449–458.
Lorenz, E.N. (1963) Deterministic nonperiodic flow, J.
Atoms. Sci., 20, 130–141.
Lowry, S.C., Fitzsimmons, A., Pravec, P., Vokrouhlický, D.,
Boehnhardt, H., Taylor, P.A., Margot, J.-L., Galád, A.,
Irwin, M., Irwin, J., and Kusnirák, P. (2007) Direct detection of the asteroidal YORP effect, Science, 316, 272–
274.
Luu, J.X. and Jewitt, D. (2002) Kuiper belt objects: Relics
from the accretion disk of the Sun, Annu. Rev. Astron.
Astrophys., 40, 63–101.
Lyapunov, A.M. (1892) The General Problem of Stabilityof Motion, Kharikov, (in Russian: English translation is
published by Taylor & Francis, London, 1992).
Mackay, R.S. (1992) Greene’s residue criterion, Nonlinearity, 5, 161–187.
Mackay, R.S. (1993) Renormalisation in Area-Preserving
Maps, World Scientific, Singapore.
Mackay, R.S. and Percival, I.C. (1985) Converse KAM:
Theory and practice, Commu. Math. Phys., 98, 469–512.
Makino, J. (1991a) A modified Aarseth code for GRAPE
and vector processors, Publ. Astron. Soc. Japan, 43, 859–
876.
Makino, J. (1991b) Optimal order and time-step criterion
for Aarseth-type N -body integrators, Astron. J., 369,
200–212.
Makino, J. (2005) GRAPE and Project Milkyway, J. Korean Astron. Soc., 38, 165–168.

Trends in 20th Century Celestial Mechanics
Makino, J. and Aarseth, S.J. (1992) On a Hermite integrator with Ahmad-Cohen scheme for gravitational manybody problems, Publ. Astron. Soc. Japan, 44, 141–151.
Makino, J. and Taiji, M. eds. (1998) Scientific Simulations
with Special-Purpose Computers — the GRAPE systems,
John Wiley & Sons, New York.
Makino, J., Taiji, M., Ebisuzaki, T., and Sugimoto, D.
(1997) GRAPE–4: A massively parallel special-purpose
computer for collisional N -body simulations, Astrophys.
J., 480, 432–446.
Makino, J., Fukushige, T., Koga, M., and Namura, K.
(2003) GRAPE–6: Massively-parallel special-purpose
computer for astrophysical particle simulations, Publ.
Astron. Soc. Japan, 55, 1163–1187.
Malhotra, R. (1990) Capture probabilities for secondary
resonances, Icarus, 87, 249–264.
Malhotra, R. (1993) The origin of Pluto’s peculiar orbit,
Nature, 365, 819–821.
Malhotra, R. (1994) A mapping method for the gravitational few-body problem with dissipation, Celes. Mech.
Dyn. Astron., 60, 373–385.
Malhotra, R. (1995) The origin of Pluto’s orbits: implications for the solar system beyond Neptune, Astron. J.,
110, 420–429.
Malhotra, R. and Dermott, S.F. (1990) The role of secondary resonances in the orbital history of Miranda,
Icarus, 85, 444–480.
Mandelbrot, B. (1977) Fractals: Form, Chance, and Dimension, W. H. Freeman & Co., San Francisco.
Mardling, R. and Aarseth, S.J. (1999) Dynamics and stability of three-body systems, in Steves, B.A. and Roy, A.
eds., The Dynamics of Small Bodies in the Solar System,
Kluwer Academic Publishers, Dordrecht, 385–392.
Mardling, R,A. and Aarseth, S.J. (2001) Tidal interactions
in star cluster simulations, Mon. Not. R. Astron. Soc.,
321, 398–420.
Marshal, C. (1990) The Three-Body Problem, Elsevier, Amsterdam.
Marzari, F. and Scholl, H. (1998) Capture of Trojans by a
growing proto-Jupiter, Icarus, 131, 41–51.
Marzari, F., Scholl, H., Murray, C., and Lagerkvist, C.
(2002) Origin and evolution of trojan asteroids, in Bottke, W.F., Cellino, A., Paolicchi, P., and Binzel, R.P.
eds., Asteroids III, The University of Arizona Press, Tucson, Arizona, 725–738.
Mather, J.N. (1982) Existence of quasi-periodic orbits for
twist homeomorphisms of the annulus, Topology, 21,457–467.
Mather, J.N. (1984) Non-existence of invariant circles, Ergod. Th. & Dynam. Sys., 4, 301–309.
Mather, J.N. (1985) More Denjoy invariant sets for area
preserving diffeomorphisms, Comment. Math. Helv., 60,
508–557.
Mather, J.N. (1991) Variational construction of orbits for
twist homeomorphisms, J. Amer. Math. Soc., 4, 203–267.

105

Matsuoka, T. (1983) The number and linking of periodic
solutions of periodic systems, Invent. Math., 70, 319–340.
Matsuoka, T. (1993) The Burau representation of the braid
group and the Nielsen–Thurston classification, Contemporary Mathematics, 152, 229–248.
Matukuma, T. (1933) On the periodic orbits in Hill’s case,Proc. imp. Acad. Japan, 9, 364–366.
May, R.M. (1976) Simple mathematical models with very
complicated dynamics, Nature, 261, 459–467.
McGehee, R. (1974) Triple collision in the collinear threebody problem, Invent. Math., 27, 191–227.
Merman, G.A. (1958) Qualitative investigations of thethree-body problem, Bull. Inst. Teoret. Astronom., 6,
687–712, (in Russian).
Michel, P. and Froeschlé, C. (1997) The location of linear
secular resonances for semimajor axes smaller than 2 AU,
Icarus, 128, 230–240.
Michel, P. and Thomas, F. (1996) The Kozai resonance fornear-Earth asteroids with semimajor axis smaller than 2
AU, Astron. Astrophys., 307, 310–318.
Michel, P. and Valsecchi, G.B. (1997) Numerical experiments on the efficiency of second-order mixed-variable
symplectic integrators for N -body problems, Celes.
Mech. Dyn. Astron., 65, 355–371.
Mignard, F. (1975) Highly eccentric satellite – Application
to Nereid, Astron. Astrophys., 43, 359–379.
Mignard, F. (1982) Long time integration of the Moon’s orbit, in Brosche, P. and Sündermann, J. eds., Tidal Friction and the Earth’s Rotation II, Springer, Berlin, 67–91.
Mikkola, S. (1997) Practical symplectic methods with timetransformation for the few-body problem, Celes. Mech.
Dyn. Astron., 67, 145–165.
Mikkola, S. (1998) Non-canonical perturbations in symplectic integration, Celes. Mech. Dyn. Astron., 68, 249–
255.
Mikkola, S. (1999) Efficient symplectic integration of satellite orbits, Celes. Mech. Dyn. Astron., 74, 275–285.
Mikkola, S. and Aarseth, S.J. (1990) A chain regularization
method for the few-body problem, Celes. Mech. Dyn.
Astron., 47, 375–390.
Mikkola, S. and Aarseth, S.J. (1993) An implementation of
N –body chain regularization, Celes. Mech. Dyn. Astron.,57, 439–459.
Mikkola, S. and Aarseth, S.J. (1998) An efficient integration
method for binaries in N-body simulations, New Astron.,
3, 309–320.
Mikkola, S. and Aarseth, S.J. (2002) A time-transformed
leapfrog scheme, Celes. Mech. Dyn. Astron., 84, 343–354.
Mikkola, S. and Innanen, K. (1999) Symplectic tangent
map for planetary motions, Celes. Mech. Dyn. Astron.,74, 59–67.
Mikkola, S. and Palmer, P. (2000) Simple derivation of
symplectic integrators with first order correctors, Celes.
Mech. Dyn. Astron., 77, 305–317.
Mikkola, S. and Tanikawa, K. (1998a) Does Kozai reso-

106

Takashi Ito and Kiyotaka Tanikawa

nance drive CH Cygni?, Astron. J., 116, 444–450.
Mikkola, S. and Tanikawa, K. (1998b) Does Kozai resonance drive CH Cygni?, Astron. J., 116, 444–450.
Mikkola, S. and Tanikawa, K. (1999a) Algorithmic regularization of the few-body problem, Mon. Not. R. Astron.
Soc., 310, 745–749.
Mikkola, S. and Tanikawa, K. (1999b) Explicit symplectic algorithms for time-transformed Hamiltonians, Celes.
Mech. Dyn. Astron., 74, 287–295.
Mikkola, S., Brasser, R., Wiegert, P., and Innanen, K.
(2004) Asteroid 2002 VE68: A quasi-satellite of Venus,
Mon. Not. R. Astron. Soc., 351, L63–L65.
Mikolajewski, M., Mikolajewska, J., and Khudyakova, T.
(1990) The long-period symbiotic binary CH Cygni I. A
hundred years’ history of variability, Astron. Astrophys.,
235, 219–233.
Milani, A. and Farinella, P. (1994) The age of the Veritas
asteroid family deduced by chaotic chronology, Icarus,
370, 40–41.
Milani, A. and Knez̆ević, Z. (1990) Secular perturbation
theory and computation of asteroid proper elements, Celes. Mech., 49, 347–411.
Milani, A. and Nobili, A.M. (1988) Integration error oververy long time spans, Celes. Mech., 43, 1–34.
Milani, A. and Nobili, A.M. (1997) Stable chaos in the asteroid belt, Icarus, 125, 13–31.
Milani, A., Nobili, A.M., and Carpino, M. (1987) Secular variations of the semimajor axes: Theory and experiments, Astron. Astrophys., 172, 265–279.
Milani, A., Nobili, A.M., and Carpino, M. (1989) Dynamics
of Pluto, Icarus, 82, 200–217.
Milankovitch, M. (1920) Théorie mathématique des phénomènes thermique produis par la radiation solaire,
Gautier–Villars, Paris.
Milankovitch, M. (1930) Mathematische Klimalehre und astronomische Theorie der Klimaschwankungen, SpringerVerlag, in “Handbuch der Klimatologie” (Köppen andGeiger eds.), Band 1. Teil A.
Milankovitch, M. (1941) Kanon der Erdbestrahlung undseine Anwendung auf das Eiszeitproblem, Vol. 133 ofKöniglich Serbische Academie Publication, Königlich
Serbische Academie.
Miller, R.H. (1964) Irreversibility in small stellar dynamicalsystems, Astrophys. J., 140, 250–585.
Miller, R.H. (1971a) Experimental studies of the numericalstability of the gravitational N -body problem, J. Comp.
Phys., 8, 449–463.
Miller, R.H. (1971b) Partial iterative refinements, J. Comp.
Phys., 8, 464–471.
Montenbruck, O. and Gill, E. (2000) Satellite Orbits,
Springer-Verlag, Berlin.
Moons, M. and Morbidelli, A. (1993) The main mean motion commensurabilities in the planar circular and elliptic problem, Celes. Mech. Dyn. Astron., 57, 99–108.
Moons, M. and Morbidelli, A. (1995) Secular resonance inmean motion commensurabilities: the 4/1, 3/1, 5/2, and7/3 cases, Icarus, 114, 33–50.
Moore, G.E. (1965) Cramming more components onto integrated circuits, Electronics, 38, 8, 114–117.
Moore, C. (1993) Braids in classical dynamics, Phys. Rev.
Lett., 70, 3675–3679.
Morais, M.H.M. and Morbidelli, A. (2006) The populationof Near Earth Asteroids in coorbital motion with Venus,
Icarus, 185, 29–38.
Morales-Ruiz, J.J. and Ramis, J.P. (2001a) Galoisian obstruction to integrability of Hamiltonian systems. I,
Methods Appl. Anal., 8, 33–96.
Morales-Ruiz, J.J. and Ramis, J.P. (2001b) Galoisian obstruction to integrability of Hamiltonian systems. II,
Methods Appl. Anal., 8, 97–112.
Morbidelli, A. (1997) Chaotic diffusion and the origin ofcomets from the 2/3 resonance in the Kuiper belt, Icarus,
127, 1–12.
Morbidelli, A. (2002) Modern integrations of solar system
dynamics, Annu. Rev. Earth Planet. Sci., 30, 89–112.
Morbidelli, A. and Froeschlé, C. (1996) On the relationship between Lyapunov times and macroscopic instability times, Celes. Mech. Dyn. Astron., 63, 227–239.
Morbidelli, A. and Giorgilli, A. (1995) Superexponential
stability of KAM tori, Journ. Stat. Phys., 78, 1607–1617.
Morbidelli, A. and Guzzo, M. (1997) The Nekhoroshev
theorem and the asteroid belt dynamical system, Celes.
Mech. Dyn. Astron., 65, 107–136.
Morbidelli, A. and Henrard, J. (1991a) The main secular resonances ν6 , ν5 and ν16 in the asteroid belt, Celes.
Mech. Dyn. Astron., 51, 169–197.
Morbidelli, A. and Henrard, J. (1991b) Secular resonancesin the asteroid belts: theoretical perturbation approachand the problem of their location, Celes. Mech. Dyn.
Astron., 51, 131–167.
Morbidelli, A. and Moons, M. (1993) Secular resonance inmean motion commensurabilities: the 2/1 and 3/2 cases,
Icarus, 102, 316–332.
Morbidelli, A. and Moons, M. (1995) Numerical evidence
on the chaotic nature of the 3/1 mean motion commensurability, Icarus, 115, 60–65.
Morbidelli, A. and Nesvorný, D. (1999) Numerous weak
resonances drive asteroids toward terrestrial planets orbits, Icarus, 139, 295–308.
Morbidelli, A., Petit, J.-M., Gladman, B., and Chambers, J.
(2001) A plausible cause of the late heavy bombardment,
Meteor. Planet Sci., 36, 371–380.
Morbidelli, A., Levison, H.F., Tsiganis, K., and Gomes, R.
(2005) Chaotic capture of Jupiter’s Trojan asteroids inthe early solar system, Nature, 435, 462–465.
Moro-Martı́n, A. and Malhotra, R. (2002) A Study of the
dynamics of dust from the Kuiper belt: Spatial distribution, Astron. J., 124, 2305–2321.
Moro-Martı́n, A. and Malhotra, R. (2003) Dynamical models of Kuiper belt dust in the inner and outer solar sys-

Trends in 20th Century Celestial Mechanics
tem, Astron. J., 125, 2255–2265.
Morse, M. (1946) George David Birkhoff and his mathematical work, Bull. Amer. Math. Soc., 52, 357–391.
Morse, M. and Hedlund, G. (1938) Symbolic dynamics,
Amer. J. Math., 60, 815–866.
Moser, J.K. (1962) On invariant curves of area-preserving
mappings of an annulus, Nachr. Akad. Wiss. Gottingen,
II Math. Phys., 1, 1–20.
Moser, J.K. (1974) Stable and Random Motions in Dynamical Systems, Princeton University Press.
Murray, C.D. and Dermott, S.F. (1999) Solar System Dynamics, Cambridge University Press, Cambridge.
Murray, N. and Holman, M. (1997) Diffusive chaos in theouter asteroid belt, Astron. J., 114, 1246–1257.
Murray, N. and Holman, M. (1999) The origin of chaos in
the outer solar system, Science, 283, 1877–1881.
Murray, N., Holman, M., and Potter, M. (1998) On theorigin of chaos in the asteroid belt, Astron. J., 116, 2583–
2589.
Nagasawa, M. and Ida, S. (2000) Sweeping secular resonances in the Kuiper belt caused by depletion of the solar
nebula, Astron. J., 120, 3311–3322.
Nakagawa, K. and Yoshida, H. (2001a) A list of all integrable two-dimensional homogeneous polynomial potentials with a polynomial integral of order at most four inthe momenta, J. Phys. A: Math. Gen., 34, 8611–8630.
Nakagawa, K. and Yoshida, H. (2001b) A necessary condition for the integrability of homogeneous Hamiltoniansystems with two degrees of freedom, J. Phys. A: Math.
Gen., 34, 2137–2148.
Nakazawa, K., Ida, S., and Nakagawa, Y. (1989a) Collisional probability of planetesimals revolving in the solar
gravitational field I. Basic formulation, Astron. Astrophys., 220, 293–300.
Nakazawa, K., Ida, S., and Nakagawa, Y. (1989b) Collisional probability of planetesimals revolving in the solar
gravitational field II. The validity of the two-body approximation, Astron. Astrophys., 221, 342–347.
Nekhoroshev, N.N. (1977) An exponential estimate of thetime of stability of nearly integrable Hamiltonian systems, Russ. Math. Survey, 32, 1–65.
Nesvorný, D. and Ferraz-Mello, S. (1997) Chaotic diffusionin the 2/1 asteroidal resonance. An application of thefrequency map analysis, Astron. Astrophys., 320, 672–
680.
Nesvorný, D., Bottke, W.F., Dones, L., and Levison, H.F.
(2002) The recent breakup of an asteroid in the mainbelt region, Nature, 417, 720–722.
Nesvorný, D., Bottke, W.F., Levison, H., and Dones, L.
(2003) Recent origin of the solar system dust bands, Astrophys. J., 591, 486–497.
Nobili, A.M., Milani, A., and Carpino, M. (1989) Fundamental frequencies and small divisors in the orbits of theouter planets, Astron. Astrophys., 210, 313–336.
Noll, K.S., Stephens, D.C., Grundy, W.M., Millis, R.L.,

107

Spencer, J., Buie, M.W., Tegler, S.C., Romanishin, W.,
and Cruikshank, D.P. (2002) Detection of two binarytrans-Neptunian objects, 1997 CQ29 and 2000 CF105,
with the Hubble Space Telescope, Astron. J., 124, 3424–
3429.
Novikov, S.P. ed. (1981) Integrable Systems, No. 60 in London Mathematical Society lecture note series, Cambridge
University Press, Cambridge.
Ohtsuki, K. and Ida, S. (1998) Planetary rotation by accretion of planetesimals with nonuniform spatial distribution formed by the planet’s gravitational perturbation,
Icarus, 131, 393–420.
Okumura, S.K., Makino, J., Ebisuzaki, T., Fukushige, T.,
Ito, T., Sugimoto, D., Hashimoto, E., Tomida, K., andMiyakawa, N. (1993) Highly parallelized special-purpose
computer, GRAPE–3, Publ. Astron. Soc. Japan, 45, 329–
338.
Oort, J.H. (1950) The structure of a cloud of comets surrounding the solar system and a hypothesis concerning
its structure, Bull. Astron. Inst. Neth., 11, 91–110.
Paillard, D. (2001) Glacial cycles: Toward a new paradigm,
Rev. Geophys., 39, 3, 325–346.
Painlevé, P. (1897) Léçon sur la Theorie analytique desequation differentielle, A. Hermann, Paris.
Pendleton, Y.J. and Black, D.C. (1983) Further studies on
criteria for the onset of dynamical instability in generalthree-body systems, Astron. J., 88, 1415–1419.
Petit, J.-M. and Hénon, M. (1987) A numerical simulationof planetary rings. I—Binary encounters, Astron. Astrophys., 173, 389–404.
Petit, J.-M., Morbidelli, A., and Valsecchi, G.B. (1999)
Large scattered planetesimals and the excitation of thesmall body belts, Icarus, 141, 367–387.
Plummer, H.C. (1960) An Introductory Treatise on Dynamical Astronomy, Dover, New York.
Poincaré, H. (1890) Sur le problème des trois corps et leséquation de la dynamique, Acta Math., XIII, 1–270.
Poincaré, H. (1892) Les Méthode nouvelles de laMéchanique céleste, Gauthier–Villars, Tome 1 (1892), 2(1892), 3 (1899). English translation is edited and introduced by D. L. Goroff as “New Methods in Celestial Mechanics”, 1993, American Institute of Physics. Subtitlesof the English version are: 1. Periodic and Asymptotic
Solutions, 2. Approximations by Series, and 3. Integral
Invariants and Asymptotic Properties of Certain Solutions.
Poincaré, H. (1912) Sur un théorème de géométrie, Rend.
Circ. Mat. Palermo, 33, 375–407.
Preto, M. and Tremaine, S. (1999) A class of symplectic
integrators with adaptive time step for separable hamiltonian systems, Astron. J., 118, 2532–2541.
Quinlan, G.D. (1999) Resonances and instabilities in symmetric multistep methods, astro-ph/9901136.
Quinlan, G.D. and Tremaine, S. (1990) Symmetric multistep methods for the numerical integration of planetary

108

Takashi Ito and Kiyotaka Tanikawa

orbits, Astron. J., 100, 1694–1700.
Quinn, T. and Tremaine, S. (1990) Roundoff error in longterm planetary orbit integration, Astron. J., 99, 1016–
1023.
Quinn, T., Tremaine, S., and Duncan, M. (1990) Planetaryperturbations and the origin of short-period comets, Astrophys. J., 355, 667–679.
Quinn, T., Tremaine, S., and Duncan, M. (1991) A three
million year integration of the earth’s orbit, Astron. J.,
101, 2287–2305.
Quintana, E.V. and Lissauer, J.J. (2006) Terrestrial planet
formation surrounding close binary stars, Icarus, 185,
1–20.
Quintana, E.V., Lissauer, J.J., Chambers, J.E., and Duncan, M.J. (2002) Terrestrial planet formation in the α
Cenrauri system, Astrophys. J., 576, 982–996.
Rauch, K.P. and Holman, M. (1999) Dynamical chaos in
the Wisdom-Holman integrator: Origins and solutions,Astron. J., 117, 1087–1102.
Richardson, D.L. and Walker, C.F. (1989) Numerical simulation of the nine-body planetary system spanning two
million years, J. Astronaut. Sci., 37, 159–182.
Roy, A.E., Walker, I.W., Macdonald, A.J., Williams, K.,
I.P. Fox, Murray, C.D., Milani, A., Nobili, A.M., Message, P.J., Sinclair, A.T., and Carpino, M. (1988) Project
LONGSTOP, Vistas Astron., 32, 95–116.
Rubincam, D.P. (1995) Asteroid orbit evolution due tothermal drag, J. Geophys. Res., 100, 1585–1594.
Rubincam, D.P. (1998) Yarkovsky thermal drag on small
asteroid and Mars–Earth delivery, J. Geophys. Res., 103,
1725–1732.
Rubincam, D.P. (2000) Radiative spin-up and spin-downof small asteroids, Icarus, 148, 2–11.
Ruelle, D. and Takens, F. (1971) On the nature of turbulence, Comm. Math. Phys., 23, 343–344.
Ruth, R.D. (1983) A canonical integration technique, IEEE
Trans. Nucl. Sci., NS-30, 2669–2671.
Ryder, G. (1990) Lunar samples, lunar accretion and theearly bombardment of the Moon, Eos, 71, 10, 313, 322–
395.
Saari, D.G. (1971) Implications of the inverse-squaremodel, Astrophys. J., 165, 399–1971.
Saha, P. and Tremaine, S. (1992) Symplectic integratorsfor solar system dynamics, Astron. J., 104, 1633–1640.
Saha, P. and Tremaine, S. (1994) Long-term planetary integrations with individual time steps, Astron. J., 108,
1962–1969.
Saha, P. and Tremaine, S. (1997) A parallel integration
method for solar system dynamics, Astron. J., 114, 409–
415.
Saito, M. (2005) Rectilinear three-body problem using symbolic dynamics, PhD thesis, Graduate University for Advanced Studies, Hayama, Kanagawa, Japan.
Sano, M.M. (2004) The classical Coulomb three-body problem in the collinear eZe configuration, J. Phys. A: Math.

Gen., 37, 803–822.
Sanz-Serna, J.M. and Calvo, M.P. (1994) Numerical Hamiltonian Problems, Chapman & Hall, London.
Sasao, T., Okubo, S., and Saito, M. (1980) A simple theory on the dynamical effects of a stratified fluid core
upon nutational motion of the Earth, in Fedrov, E.P.,
Smith, M.L., and Bender, P.L. eds., Nutation andthe Earth’s Rotation, Reidel, Dordrecht, Netherlands,
165–183, International Astronomical Union Symposium
No.78.
Scholl, H., Froeschlé, Ch., Kinoshita, H., Yoshikawa, M.,
and Williams, J.G. (1989) Secular resonances, inGehrels, T. ed., Asteroids II, The University of Arizona
Press, Tucson, 845–861.
Sekiguchi, M. and Tanikawa, K. (2004) On the symmetriccollinear four-body problem, Publ. Astron. Soc. Japan,
56, 235–251.
Sharaf, S.G. and Budnikova, N.A. (1969a) On secular perturbations in the elements of the earth’s orbit and theirinfluence on the climates in the geological past, Bull.
Inst. Theor. Astron., 11, 231–261, in Russian. English
translation is in NASA TT F–12, 467 , 1–37.
Sharaf, S.G. and Budnikova, N.A. (1969b) Secular perturbations in the elements of the Earth’s orbit and the astronomical theory of climate variations, Trudy Inst. Theor.
Astron., 14, 48–84, in Russian.
Shirai, T. and Fukushima, T. (2000) Numerical convolution
in the time domain and its application to the nonrigidearth nutation theory, Astron. J., 119, 2475–2480.
Shirai, T. and Fukushima, T. (2001) Construction of a newforced nutation theory of the nonrigid earth, Astron. J.,
121, 3270–3283.
Shoemaker, E.M., Shoemake, C.S., and Wolfe, R.F. (1989)
Trojan asteroids: Populations, dynamical structure andorigin of the L4 and L5 swarms, in Gehrels, T. ed., Asteroids II, The University of Arizona Press, Tucson, 487–
523.
S̆idlichovský, M. (1999) On stable chaos in the asteroid
belt, Celes. Mech. Dyn. Astron., 73, 77–86.
Siegel, C.L. (1941) Der Dreierstoss, Acta Math., 42, 127–
168.
Siegel, C.L. and Moser, J.K. (1971) Lectures on Celestial Mechanics, Springer-Verlag, Berlin, Original German
edition is C. L. Siegel (1956) “Vorlesungen über Himmelsmechanik”, Die Grundlehren der mathematischen
Wissenschaften in Einzeldarstellungen, Band 85.
Simó, C. (2002) Dynamical properties of the figure eightsolution of the three-body problem, in Chenciner, A.,
Cushman, R., Robinson, C., and Xia, Z. eds., Celestial
Mechanics: Dedicated to Donald Saari for his 60th Birthday, American Mathematical Society, Providence, RI,
209–228.
Sitnikov, K. (1960) The existence of oscillatory motions in
the three-body problem, Dokl. Akad. Nauk USSR, 133,
303–306.

Trends in 20th Century Celestial Mechanics
Skeel, R.D. and Biesiadecki, J.J. (1994) Symplectic integration with variable stepsize, Ann. Numer. Math., 1,
191–198.
Smale, S. (1967) Differential dynamical systems, Bull.
Amer. Math. Soc., 73, 747–817.
Smart, W. (1953) Celestial Mechanics, Longmans, London.
Soffel, M.H. ed. (1989) Relativity in astronomy, celestialmechanics and geodesy, Springer-Verlag, Berlin.
Sôma, M., Kawabata, K., and Tanikawa, K. (2004) Unitsof time in ancient China and Japan, Publ. Astron. Soc.
Japan, 887–904.
Souchay, J. and Kinoshita, H. (1991) Comparison of new
nutation series with numerical integration, Celes. Mech.
Dyn. Astron., 52, 45–55.
Stephens, D.C. and Noll, K.S. (2006) Detection of six transNeptunian binaries with NICMOS: A high fraction ofbinaries in the cold classical disk, Astron. J., 131, 1142–
1148.
Stephenson, F.R. (1997) Historical Eclipses and Earth’s
Rotation, Cambridge University press.
Stephenson, F.R. and Morrison, L.V. (1995) Long-term
fluctuations in the Earth’s rotation: 700BC to AD1990,
Phil. Trans. R. astr. Soc. London, A351, 165–202.
Sterling, D., Dullin, H.R., and Meiss, J.D. (1999) Homoclinic bifurcations for the Hénon map, Physica D, 134,
153–184.
Stern, S.A. (1996) Signatures of collisions in the Kuiper
disk, Astron. Astrophys., 310, 999–1010.
Strom, R.G., Malhotra, R., Ito, T., Yoshida, F., andKring, D.A. (2005) The origin of impacting objects in
the inner solar system, Science, 309, 1847–1850.
Strömgren, E. (1935) Fortgesetzte untersuchungen uber
asymptotische bahnen im probleme restreint, Bulletine
Astronomiques, 9, 87, (also in Publ. Copenhague Obs.
no. 100).
Sugimoto, D., Chikada, Y., Makino, J., Ito, T.,
Ebisuzaki, T., and Umemura, M. (1990) A specialpurpose computer for gravitational many-body problem,
Nature, 345, 33–35.
Sundman, K. (1912) Mémoire sur le probléme des trois
corps, Acta Math., 36, 105–179.
Sussman, G.J. and Wisdom, J. (1988) Numerical evidence
that the motion of Pluto is chaotic, Science, 241, 433–
437.
Sussman, G.J. and Wisdom, J. (1992) Chaotic evolution of
the solar system, Science, 257, 56–62.
Szebehely, V. (1997) Open problems on the eve of the next
millennium, Celes. Mech. Dyn. Astron., 65, 205–211.
Szebehely, V. and Peters, C.F. (1967) Complete solution ofa general problem of three bodies, Astron. J., 72, 876–
883.
Szebehely, V. and Zare, K. (1977) Stability of classicaltriplets and their hierarchy, Astron. Astrophys., 58, 145–
152.
Taff, L.G. ed. (1985) Celestial Mechanics: A Computational

109

Guide for the Practitioner, Wiley, New York.
Tanikawa, K. (1983) Impossibility of the capture of retrograde satellites in the restricted three-body problem,
Celes. Mech., 29, 367–402.
Tanikawa, K. (2000) A Search for collision orbits in the freefall three-body problem II., Celes. Mech. Dyn. Astron.,76, 157–185.
Tanikawa, K. and Mikami, Y. (2000) Space Occultation,
in Proceedings of the 32nd Symposium on Celestial Mechanics, National Astronomical Observatory of Japan,
Mitaka, Tokyo, 181-8555, Japan, 228–233, March 15–17,
Hayam, Kanagawa, Japan.
Tanikawa, K. and Mikkola, S. (2000a) One-dimensional
three-body problem via symbolic dynamics, Chaos, 10,
649–657.
Tanikawa, K. and Mikkola, S. (2000b) Triple collision in the
one-dimensional three-body problem, Celes. Mech. Dyn.
Astron., 76, 23–34.
Tanikawa, K. and Sôma, M. (2001) Reliability of the totality of the eclipse in AD 628 in Nihongi, in Proc. 33th
Symp. Celes. Mech., National Astronomical Observatory,
Mitaka, Tokyo, Japan, 432–446.
Tanikawa, K. and Sôma, M. (2004a) ∆T and the tidal acceleration of the lunar motion from eclipses observed atplural sites, Publ. Astron. Soc. Japan, 879–885.
Tanikawa, K. and Sôma, M. (2004b) On the totality ofthe eclipse in AD 628 in the Nihongi, Publ. Astron. Soc.
Japan, 56, 215–224.
Tanikawa, K. and Yamaguchi, Y. (1987) Coexistence of periodic points in reversible dynamical systems on a surface, J. Math. Phys., 28, 921–928.
Tanikawa, K. and Yamaguchi, Y. (1989) Coexistence ofsymmetric periodic points in the standard map, J. Math.
Phys., 30, 608–616.
Tanikawa, K. and Yamaguchi, Y. (2001) Running homoclinic and periodic points in standard-like mappings,
Progr. Theor. Phys., 105, 399–407.
Tanikawa, K. and Yamaguchi, Y. (2002a) Dynamical ordering of symmetric non-Birkhoff periodic points in reversible monotone twist mappings, Chaos, 12, 33–41.
Tanikawa, K. and Yamaguchi, Y. (2002b) Non-symmetricperiodic points in reversible maps: Examples from thestandard map, Prog. Theor. Phys., 108, 987–997.
Tanikawa, K. and Yamaguchi, Y. (2005) Generalized horseshoes in the standard mapping, Prog. Theor. Phys., 113,
261–281.
Tanikawa, K., Manabe, S., and Broucke, R. (1989) On theorigin of the planetary spin angular momentum by accretion of planetesimals: Property of collision orbits, Icarus,
79, 208–222.
Tanikawa, K., Kikuchi, N., and Sato, I. (1991) On the originof the planetary spin by accretion of planetesimals II.
Collisional orbits at the Hill surface, Icarus, 94, 112–125.
Tanikawa, K., Umehara, H., and Abe, H. (1995) A searchfor collision orbits in the free-fall three-body problem

110

Takashi Ito and Kiyotaka Tanikawa

I. Numerical procedure, Celes. Mech. Dyn. Astron., 62,
335–362.
Taylor, P.A., Margot, J.-L., Vokrouhlický, D.,
Scheeres, D.J., Pravec, P., Lowry, S.C., Fitzsimmons, A.,
Nolan, M.C., Ostro, S.J., Benner, L.A.M., Giorgini, J.D.,
and Magri, C. (2007) Spin rate of asteroid (54509) 2000
PH5 increasing due to the YORP effect, Science, 316,
274–277.
Tera, F., Papanastassiou, D.A., and Wasserburg, G.J.
(1974) Isotopic evidence for a terminal lunar cataclysm,
Earth Planet. Sci. Lett., 22, 1.
Thébault, P., Marzari, F., Scholl, H., Turrini, D., and Barbieri1, M. (2004) Planetary formation in the γ Cephei
system, Astron. Astrophys., 427, 1097–1104.
Thiele, T.N. (1896) Recherches numériques concernant dessolutions périodiques d’un cas spécial du probléme destrois corps, Astron. Nachr., 138, 1.
Thommes, E.W., Duncan, M.J., and Levison, H.F. (2003)
Oligarchic growth of giant planets, Icarus, 161, 431–455.
Tomasella, L., Marzari, F., and Vanzani, V. (1997) Evolution of the Earth obliquity after the tidal expansion of
the Moon orbit, Planet. Space Sci., 44, 427–430.
Touma, J. and Wisdom, J. (1994) Evolution of the Earth–
Moon system, Astron. J., 108, 1943–1951.
Touma, J. and Wisdom, J. (1998) Resonances in the early
evolution of the Earth–Moon system, Astron. J., 115,
1653–1663.
Touma, J. and Wisdom, J. (2001) Nonlinear core-mantle
coupling, Astron. J., 122, 1030–1050.
Tremaine, S. (1995) Is the solar system stable?, CITA–95–
3, in Proceedings of Rosseland Centenary Symposium of
Astrophysics, Oslo, June 16–17, 1994.
Trujillo, C.A. (2003) Discovering the edge of the solar system, American Scientist, 91, 424–431.
Tsiganis, K., Varvoglis, H., and Hadjidemetriou, J.D.
(2002) Stable chaos versus Kirkwood gaps in the asteroid
belt: A comparative study of mean motion resonances,
Icarus, 159, 284–299.
Tsiganis, K., Gomes, R., Morbidelli, A., and Levison, H.F.
(2005) Origin of the orbital architecture of the giant planets of the solar system, Nature, 435, 459–461.
Ueda, Y. (1979) Randomly transitional phenomena in thesystem governed by Duffing’s equation, J. Stat. Phys.,
20, 181–196.
Ueda, Y., Hayashi, C., and Akamatsu, N. (1973) Computer
simulation of nonlinear ordinary differential equations
and nonperiodic oscillations, Electron. Comm. Japan,
56A, 27.
Umehara, H. (1997) The free-fall three-body problem: Escape and collisions, PhD thesis, The Graduate University
for Advanced Studies, Hayama, Kanagawa, Japan.
Umehara, H. and Tanikawa, K. (2000) Binary and triple
collisions causing instability in the free-fall three-bodyproblem, Celes. Mech. Dyn. Astron., 76, 187–214.
Valsecchi, G,B., Jopek, T.J., and Froeschlé, Cl. (1999) Me-teoroid stream identification: a new approach - I. Theory,
Mon. Not. R. Astron. Soc., 304, 743–750.
Varadarajan, V.S. (1974) Lie groups, Lie algebras and theirrepresentation, Prentice-Hall, New Jersey.
Veillet, C., Parker, J.W., Griffin, I., Marsden, B., Doressoundiram, A., Buie, M., Tholen, D.J., Connelley, M.,
and Holman, M.J. (2002) The binary Kuiper-belt object
1998 WW31, Nature, 416, 711–713.
Vokroulický, D., Milani, A., and Chesley, S.R. (2000)
Yarkovsky effect on small near-Earth asteroids: Mathematical formulation and examples, Icarus, 148, 118–138.
von Zeipel, H. (1916) Recherches sur le mouvement des petites planétes, Ark. Astron. Mat. Fys., 11, no. 1, See also
ibid., 11, no. 7, 1–62 (1916); 12, no. 9, 1–89 (1916); 13
no. 3, 1–93 (1917).
Wako, Y. (1970) Interpretation of Kimura’s annual Z-term,
Publ. Astron. Soc. Japan, 22, 525–544.
Wan, X.-S., Huang, T.-Y., and Innanen, K.A. (2001) The
1:1 superresonance in Pluto’s motion, Astron. J., 121,
1155–1162.
Ward, W.R. (1974) Climatic variations on Mars 1. Astronomical theory of insolation, J. Geophys. Res., 79, 3375–
3386.
Ward, W.R. and Canup, R.M. (2006) The obliquity of
Jupiter, Astrophys. J., 640, L91–L94.
Ward, W.R., Murray, B.C., and Malin, M.C. (1974) Climatic variations on Mars 2. Evolution of carbon dioxide
atmosphere and polar caps, J. Geophys. Res., 79, 3387–
3395.
Weidenschilling, S.J. (2002) On the origin of binarytransneptunian objects, Icarus, 160, 212–215.
Weissman, P.R. (1995) The Kuiper belt, Annu. Rev. Astron. Astrophys., 33, 327–357.
Wetherill, G.W. (1967) Collisions in the asteroid belt, J.
Geophys. Res., 72, 2429–2444.
Wetherill, G.W. (1974) Solar system sources of meteorites
and large meteoroids, Annu. Rev. Earth Planet. Sci., 2,
303–331.
Wetherill, G.W. and Stewart, G.R. (1989) Accumulationof a swarm of small planetesimals, Icarus, 77, 330–357.
Whittaker, E.T. (1904) A Treatise on the Analytical Dynamics of Particles and Rigid Bodies, Cambridge University press, (revised in 1952).
Wiegert, P.A., Innanen, K.A., and Mikkola, S. (1997) An
asteroidal companions to the Earth, Nature, 387, 685–
686.
Wiegert, P.A., Innanen, K.A., and Mikkola, S. (1998) The
orbital evolution of near-earth asteroid 3753, Astron. J.,
115, 2604–2613.
Wiegert, P.A., Innanen, K., and Mikkola, S. (2000) The
stability of quasi satellites in the outer solar system, Astron. J., 119, 1978–1984.
Wiegert, P.A., McNeil, D., Duncan, M., and Levison, H.
(2004) Fast force algorithms and solar system integrations, Astron. J., 127, 1218–1226.

Trends in 20th Century Celestial Mechanics
Williams, J.G. (1969) Secular perturbations in the solar system, PhD thesis, University of California, Los Angeles.
Williams, J.G. (1979) Proper elements and family memberships of the asteroids, in Gehrels, T. ed., Asteroids, The
University of Arizona Press, Tucson, 1040–1063.
Williams, G.E. (1990) Tidal rhythmites: Key to the historyof the earth’s rotation and the lunar orbit, J. Phys. Earth,38, 475–491.
Williams, J.G. and Benson, G.S. (1971) Resonances in the
Neptune–Pluto system, Astron. J., 76, 167–177.
Williams, J.G. and Faulkner, J. (1981) The positions ofsecular resonance surfaces, Icarus, 46, 390–399.
Williams, D.M. and Kasting, J.F. (1997) Habitable planetswith high obliquities, Icarus, 129, 254–267.
Williams, D.M., Kasting, J.F., and Wade, R.A. (1997) Habitable moons around extrasolar giant planets, Nature,
385, 234–236.
Wisdom, J. (1980) The resonance overlap criterion and theonset of stochastic behavior in the restricted three-bodyproblem, Astron. J., 85, 1122–1133.
Wisdom, J. (1982) The origin of the Kirkwood gaps: A
mapping for asteroidal motion near the 3/1 commensurability, Astron. J., 87, 577–593.
Wisdom, J. (1983) Chaotic behavior and the origin of the
3/1 Kirkwood gap, Icarus, 56, 51–74.
Wisdom, J. (1987a) Chaotic behavior in the solar system,
Proc. R. Soc. London Ser., 413, 109–129.
Wisdom, J. (1987b) Urey prize lecture: Chaotic dynamics
in the solar system, Icarus, 72, 241–275.
Wisdom, J. (2006) Symplectic correctors for canonical heliocentric N -body maps, Astron. J., 131, 2294–2298.
Wisdom, J. and Holman, M. (1991) Symplectic maps for
the N -body problem, Astron. J., 102, 1528–1538.
Wisdom, J. and Holman, M. (1992) Symplectic maps for
the N -body problem: stability analysis, Astron. J., 104,2022–2029.
Wisdom, J., Peale, S.J., and Mignard, F. (1984) The chaotic
rotation of Hyperion, Icarus, 58, 137–154.
Wisdom, J., Holman, M., and Touma, J. (1996) Symplectic
correctors, in Marsden, J.E., Patrick, G.W., and Shadwick, W.F. eds., Integration Algorithms and ClassicalMechanics, Vol. 10 of Fields Institute Communications,
American Mathematical Society, Providence, Rhode Island, 217–244.
Woolard, E.M. (1953) Theory of the rotation of the Eartharound its center of mass, Astron. Pap. Amer. Ephem.
Naut. Alm., 15, 3–165.
Xia, Z. (1992) The existence of noncollision singularities innewtonian systems, Annals of Mathematics, 135, 411–
468.
Xia, Z. (1994) Arnold diffusion and oscillatory solutions in
the planar three-body problem, J. Differential Equations,
110, 289–321.
Yamaguchi, Y. and Tanikawa, K. (2000) Symmetrical nonBirkhoff period-3 orbits in the standard-like mappings,

111

Prog. Theor. Phys., 104, 943–954.
Yamaguchi, Y. and Tanikawa, K. (2001a) Dynamical ordering of Non-Birkhoff periodic orbits in a forced pendulum,
Prog. Theor. Phys., 106, 1097–1114.
Yamaguchi, Y. and Tanikawa, K. (2001b) Geometrical approach to the splitting of separatrix for standard-likemappings, Phys. Lett. A, 280, 33–36.
Yamaguchi, Y. and Tanikawa, K. (2001c) Non-symmetrical
non-Birkhoff period 2 orbits in the Standard mapping,
Prog. Theor. Phys., 106, 691–696.
Yamaguchi, Y. and Tanikawa, K. (2002a) Dynamical ordering of non-Birkhoff orbits and topological entropy in thestandard mapping, Prog. Theor. Phys., 107, 1117–1145.
Yamaguchi, Y. and Tanikawa, K. (2002b) Dynamical ordering of periodic accelerator modes in the standard mapping, Prog. Theor. Phys., 108, 669–678.
Yamaguchi, Y. and Tanikawa, K. (2002c) Non-Birkhoff periodic orbits in a circle mapping, Prog. Theor. Phys., 108,
231–240.
Yamaguchi, Y. and Tanikawa, K. (2003) Topological entropy in a forced Duffing system, Prog. Theor. Phys.,
110, 861–874.
Yamaguchi, Y. and Tanikawa, K. (2004a) Oscillatory orbits
in the standard mapping, Prog. Theor. Phys., 111, 689–
705.
Yamaguchi, Y. and Tanikawa, K. (2004b) Periodic orbits
in a quadruply reversible non-twist Harper map, Prog.
Theor. Phys. Lett., 111, 757–762.
Yamaguchi, Y. and Tanikawa, K. (2005a) Dynamical ordering of symmetric periodic orbits in the area preservingHénon map, Prog. Theor. Phys., 113, 935–951.
Yamaguchi, Y. and Tanikawa, K. (2005b) Evolution toward
the three-fold horseshoe in the two dimensional cubic
map, Prog. Theor. Phys., 114, 1163–1177.
Yamaguchi, Y. and Tanikawa, K. (2005c) Generalized dynamical ordering and topological entropy in the Hénonmap, Prog. Theor. Phys., 114, 763–791.
Yamamoto, S. and Mukai, T. (1998) Dust production by
impacts of interstellar dust on Edgeworth–Kuiper Belt
objects, Astron. Astrophys., 329, 785–791.
Yamashita, Y. and Maehara, H. (1979) A binary model for
CH Cyg, Publ. Astron. Soc. Japan, 307–316.
Yoshida, J. (1972) Improved criteria for hyperbolic–elliptic
motion in the general three-body problem, Publ. Astron.
Soc. Japan, 24, 391–408.
Yoshida, J. (1974) Improved criteria for hyperbolic–elliptic
motion in the general three-body problem II., Publ. Astron. Soc. Japan, 26, 367–377.
Yoshida, H. (1983a) Necessary conditions for the existenceof algebraic first integrals I., Celes. Mech., 31, 363–379.
Yoshida, H. (1983b) Necessary conditions for the existenceof algebraic first integrals II., Celes. Mech., 31, 383–399.
Yoshida, H. (1990a) Conserved quantities of symplectic integrators for Hamiltonian systems, preprint.
Yoshida, H. (1990b) Construction of higher order symplec-

112

Takashi Ito and Kiyotaka Tanikawa

tic integrators, Phys. Lett. A, 150, 262–268.
Yoshida, H. (1993) Recent progress in the theory and application of symplectic integrators, Celes. Mech. Dyn. Astron., 56, 27–43.
Yoshida, H. (1999) A new necessary condition for the integrability of Hamiltonian systems with a two-dimensional
homogeneous potential, Physica D, 128, 53–69.
Yoshida, H. (2001) Non-existence of the modified first integral by symplectic integration method, Phys. Lett. A,
282, 276–283.
Yoshikawa, M. (1987) A simple analytical model for thesecular resonance ν6 in the asteroidal belt, Celes. Mech.,40, 233–272.
Yoshikawa, M. (1989) A survey of the motions of asteroids
in the commensurabilities with Jupiter, Astron. Astrophys., 213, 436–458.

Yoshinaga, K., Kokubo, E., and Makino, J. (1999) The
stability of protoplanet systems, Icarus, 139, 328–335.
Zappalà, V., Cellino, P., A. Farinella, and Milani, A. (1994)
Asteroid families. II. Extension to unnumbered multiopposition asteroids, Astron. J., 107, 772–801.
Zare, K. and Chesley, S. (1998) Order and chaos in theplanar isosceles three-body problem, Chaos, 8, 475–494.
Zare, K. and Tanikawa, K. (2002) Reversible maps in two
degrees of freedom Hamiltonian systems, Chaos, 12, 699–
705.
Zare, K. and Tanikawa, K. (2005) Reversing symmetries ina two-parameter family of area-preserving maps,, Prog.
Theor. Phys., 113, 251–259.
Ziglin, S.L. (1983) Branching of solutions and non-existenceof first integrals in Hamiltonian mechanics. I, Funct.
Anal. Appl., 16, 181–189.