Computational Celestial Mechanics means purely analytical, computer-assisted and numerical methods with the typical feature of the necessity of carrying out a vast amount of calculations, aimed to determine dynamical features of bodies of planetary systems.

Introduction

An accurate prediction of the dynamics of the objects of the solar system often requires very long computations. Even the easiest problem provided by the two-body model deserves computational skill in solving Kepler's equation, which allows to derive the Keplerian elements of the orbit as a function of time. In passing from the two to the three body problem, an increasing complexity due to nonintegrability and chaos is a trademark of Celestial Mechanics; the three body problem motivated the development of perturbation theories aimed to find approximate solutions of the equations of motion. Refined analytical perturbative techniques, such as KAM or Nekhoroshev theory, can be applied to some problems of Celestial Mechanics under suitable assumptions; most likely, effective results often require very lengthy computations which can be implemented through computer-assisted techniques. Models of Celestial Mechanics can be studied also by numerical integrations, eventually using frequency map analysis or synthetic theories; however, due to the numerical errors introduced by the machine particular care must be taken when running over long time scales, as over the age of the solar system. Both analytical and numerical techniques are worldwide used to compute accurate ephemerides.

Keplerian motion and Kepler's equation

The two-body problem is the study of the motion of two material points, e.g. a planet and a satellite, subject to the mutual gravitational attraction. Its solution is provided by Kepler's laws, according to which the motion takes place on an ellipse. Assume that the planet is located at one focus, which is taken as the origin of a reference frame whose abscissa coincides with the perihelion line; in this frame let \(r\) (the orbital radius) and \(f\) (the true anomaly) be the polar coordinates of the satellite. Let \(a\) and \(e\) be the semimajor axis and eccentricity of the ellipse. The solution of the two-body problem (see, e.g., [Roy]) is provided by the set of formulae

\[r=a(1-e\ cos\ E)\] \[f=2 \arctan\Big(\sqrt{{{1+e}\over {1-e}}}\ \tan{E\over 2}\Big) \] \[\ell=E-e\sin E\ ,\]

where \(\ell=\ell(t)\) is the mean anomaly related to the time by \(\ell(t)=n t +\ell(0)\ ,\) \(n\) denoting the frequency of revolution. Last formula, known as Kepler's equation, must be solved to provide the eccentric anomaly \(E\) as a function of \(\ell\ ,\) and therefore of the time. Inserting such solution in the first two equations, one obtains the variation with time of the orbital radius and of the true anomaly. Its solution can be found numerically, up to a given precision, through an iterative algorithm e.g. using a Newton's method, or it can be expressed analytically by means of the Bessel's functions \(J_k(x)\) as \[ E=\ell+\sum_{k=1}^\infty {1\over k}\Big[J_{k-1}(ke)+J_{k+1}(ke)\Big]\sin(k\ell)\ . \]

Perturbation theory : The precession of the perihelion as an example

Canonical perturbation theory ([FM]) is related with nearly integrable Hamiltonian systems and provides efficient tools for Celestial Mechanics. A typical application is the computation of the precession of the perihelion of Mercury.

Consider a nearly-integrable Hamiltonian function of the form \[ H({\underline I},{\underline \varphi})=h({\underline I})+\varepsilon f({\underline I},{\underline \varphi})\ , \] where \(h\) and \(f\) are analytic functions depending on the actions \({\underline I}\) (varying on an open set of \( R^n\)) and on the angles \({\underline\varphi}\) (belonging to the standard \(n\)-dimensional torus); \(\varepsilon>0\) is a small parameter which measures the strength of the perturbation. The aim of classical perturbation theory ([FM]) is to construct a near-to-identity canonical transformation, say \(C:({\underline I},{\underline \varphi})\rightarrow ({\underline I'},{\underline \varphi'})\ ,\) which allows to remove the perturbation to higher orders in the perturbing parameter; the transformed Hamiltonian takes the form \[ H_1({\underline I'},{\underline \varphi'})\ =\ h_1({\underline I'})+\varepsilon^2 f_1({\underline I'},{\underline\varphi'})\ , \] where \(h_1\) and \(f_1\) denote the new unperturbed Hamiltonian and the new perturbing function. The integrable part of the transformed Hamiltonian is simply given by the sum of the old integrable Hamiltonian and the average of the perturbation over the angle variables. The explicit derivation of the expressions for \( h_1\) and \( f_1\) often presents a considerable computational complexity.

An example of the implementation of classical perturbation theory is the computation of the precession of the perihelion of Mercury in the Mercury-Sun-Jupiter system (within the framework of the restricted, planar, circular, three-body model). Delaunay action-angle variables are introduced, where the actions are related to the osculating semimajor axis and eccentricity of the Keplerian orbit, while the angles are the mean anomaly and the difference between the argument of perihelion and the time. Denoting by \(\varepsilon\) the mass-ratio of the primaries, the problem is described by the two degrees-of-freedom Hamiltonian \[ H(L,G,\ell,g)=-{1\over {2L^2}}-G+\varepsilon R(L,G,\ell,g)\ , \] where the perturbing function \(R=R(L,G,\ell,g)\) represents the interaction between Mercury and Jupiter. Long computations finally lead to a development of the perturbing function; the first few terms are given by \[ R(L,G,\ell,g)=R_{00}(L,G)+R_{10}(L,G)\cos\ell+R_{11}(L,G)\cos(\ell+g)+R_{12}(L,G)\cos(\ell+2g))+...\ , \] where, denoting by \(e=\sqrt{1-{G^2\over L^2}}\ ,\) the coefficients \(R_{ij}\) are given by the expressions \[ R_{00}=-{L^4\over 4}(1+{9\over {16}}L^4 +{3\over 2}e^2)+...\ ,\quad R_{10}={{L^4e}\over 2}(1+{9\over {8}}L^4)+...\ , \] \[ R_{11}=-{3\over 8}L^6(1+{5\over {8}}L^4)+...\ ,\quad R_{12}={{L^4e}\over 4}(9+5L^4)+...\ ,\quad \] To compute the precession of the perihelion, a first order perturbation theory is implemented, which provides a new integrable Hamiltonian function of the form \[ h_1(L',G')=-{1\over {2L'^2}}\ -\ G'\ +\ \varepsilon\ R_{00}(L',G')\ . \] From Hamilton's equations one obtains that, to the lowest order, the rate of precession of the perihelion is given by \[ \dot\gamma=\varepsilon\ {{\partial R_{00}(L',G')}\over {\partial G'}} ={3\over 4}\varepsilon L'^2G'\ . \] From astronomical data it turns out that \(\varepsilon=9.49\cdot 10^{-4}\ ;\) setting to one the Jupiter-Sun distance, the semimajor axis of Mercury amounts to 0.07440; taking into account that the eccentricity is 0.2056 and that the orbital period of Jupiter is about 11.86 years, one obtains \[ \dot\gamma=154.65\ \ {{\rm arcsecond}\over {\rm century}}\ , \] which represents a first order contribution due to Jupiter to the precession of perihelion of Mercury. More refined evaluations must be obtained performing the above computations to higher orders. The value found by Leverrier on the basis of the data available in the year 1856 was of 152.59 arcsecond/century.

Computer-assisted stability results

Computer-assisted proofs

The construction of high order perturbation theories requires very long computations, which can be treated through a computer implementation. Nevertheless the computer introduces rounding-off and propagation errors; the development of the interval arithmetic technique allows to make such computations rigorous. Computer-assisted proofs has been successfully applied in the context of the four-color problem, the existence of the Lorenz attractor or the computation of effective stability estimates through KAM theory.

KAM theory

KAM theory ([K], [A], [M]) can be implemented to investigate the stability of some models of Celestial Mechanics, most notably the N-body problem. In the framework of the three-body problem, Arnold ([A]) proved the following result: "If the masses, eccentricities and inclinations of the planets are sufficiently small, then for the majority of initial conditions the true motion is conditionally periodic and differs little from Lagrangian motion with suitable initial conditions throughout an infinite interval of time". Arnold gave a complete proof for the case of three coplanar bodies; the spatial three-body problem was studied in [LR], making use of Poincaré variables, the Jacobi's reduction of the nodes and Birkhoff's normal form. The full proof of Arnold's theorem was given by Féjoz on the basis of Herman's results on the planetary problem. Concrete estimates on the perturbing parameter (i.e., the Sun-Jupiter mass-ratio) ensuring the existence of KAM invariant tori were computed by M. Hénon ([H]); nevertheless, direct applications of KAM theory to the three-body problem give analytical results which are much smaller than the astronomical data. In particular, the application of Arnold's theorem to the restricted three-body problem is valid provided the mass-ratio of the primaries is less than \(10^{-333}\ ,\) while Moser's theorem yields \(10^{-48}\) (see [H]), still far from the present Jupiter-Sun mass-ratio amounting to about \(10^{-3}\ .\)

A real challenge comes from the interaction between KAM theory and computer-assisted techniques: in the restricted, planar, circular three-body problem, the existence of invariant tori has been proved for the actual value of the Jupiter-Sun mass ratio ([CC]). This results stems from a combination of very long effective estimates with a 15,000 line (Fortran) program, which is used to compute an approximate solution using interval arithmetic.

KAM estimates and computer-assisted implementations have been used to study several models of Celestial Mechanics, like the spin-orbit model ([C]) or the planetary problem ([LG]), providing results in agreement with the astronomical observations.

Nekhoroshev theory

The stability of a dynamical system for exponentially long times can be investigated through Nekhoroshev theorem ([N]). In order to obtain optimal estimates, computer-assisted techniques have been implemented in a number of examples. Interesting applications in Celestial Mechanics are found in connection to the stability of the triangular Lagrangian points, which can be exploited through Nekhoroshev theory combined with computer-assisted implementations of suitable Birkhoff normal forms.

Long-time integrations

Stability of the solar system

The first numerical integrations of the solar system dynamics were performed at the end of last century. G.J. Sussman and J. Wisdom (MIT, USA) built the Digital Orrery, a computer devoted to the simulation of a planetary system. Their model ([SW]) was composed by the four giant planets (Jupiter, Saturn, Uranus, Neptune) and by Pluto. The numerical integration of the equations of motion gives that the external planets are essentially stable over a time of 845 million years, while Pluto shows a chaotic behavior. At the same time J. Laskar (Bureau des Longitudes, Paris) applied a numerical procedure based on Leverrier perturbation theory (see [L]). The results, also extended and refined in later works, concerned the dynamics of all planets except Pluto over an interval from -15 to +10 Gyr. The giant planets (from Jupiter to Neptune) are confirmed to be regular, a moderate chaotic behavior is observed for Venus and the Earth, while Mercury and Mars provide a definite chaotic dynamics with large excursions of the eccentricity of Mercury, possibly increasing beyond 0.6 in 5 Gyr. Laskar results show that one cannot make long predictions on the dynamics of the inner solar system: after 100 million years, an error of 15 meters in the initial position of the Earth leads to an error comparable to the Earth-Sun distance. On the other hand, T. Ito and K. Tanikawa, performing direct numerical integrations over \(\pm 4\) Gyr, found an overall stability of the planetary motion with just a moderate chaotic diffusion of the eccentricity of Mercury. The contribution of general relativity was also considered and led to substantial changes as far as the perihelion frequency of Mercury is concerned. In contrast to the Newtonian model where large chaotic zones are observed, general relativity does not present orbits trapped in a \(g_1=g_5\) resonance (i.e., a resonance between the frequencies of perihelion of Mercury and Jupiter).

Chaotic obliquity of the planets

The change of orientation of the planetary axes of rotation was investigated in [L], as a result of the effect of the precessional motion and of the gravitational interaction with the other bodies of the solar system. Numerical investigations provide results on the obliquity, defined as the angle between the equatorial and orbital planes; its value varies as a consequence of the perturbations induced on the orbital plane. A frequency map analysis shows a large chaotic zone of the obliquity of Mars, ranging from \(0^o\) to \(60^o.\) As for the Earth, the presence of the Moon stabilizes the rotation axis, since a chaotic zone is found well beyond the actual obliquity of the Earth (about \(23^o\)), i.e. between \(60^o\) and \(90^o.\) A striking remark is that without the Moon the chaotic region would extend from \(0^o\) to \(85^o,\) thus provoking drastic effects on the Earth's climate.

Synthetic theories

Synthetic theories allow to study the long-time behavior of the solar system dynamics, through the analysis of a reduced number of quantities (frequencies, amplitudes and phases) with respect to the large amount of data provided by a numerical integration of the equations of motion. The term synthetic is intended as opposed to analytic, since it uses a compressed output of the numerical integration. The output is typically represented by a discrete time series, built from one observable sampled at a given interval of time. Using spectral theory one computes the leading frequencies and the corresponding Fourier coefficients. Therefore a model of the dynamics is represented by a finite Fourier series with numerical coefficients whose arguments are given by linear functions of the time. Synthetic theories have been successfully used, for example, to compute ephemerides, to analyze satellite's perturbations or to investigate asteroid's resonances.

Ephemerides computation

Planetary ephemerides are computed by several institutions of the world, eventually using different techniques; a deep computational effort is performed with the goal to obtain the best accuracy on the basis of the available observational data. The ephemerides are then used for several purposes: astronomical studies, dynamical investigations and comparison with analytic theories, spacecraft navigation and planning, observational predictions and reductions. The JPL ephemeris program started in the 1960's; the ephemerides are computed through a numerical integration of the equations of motion of the main bodies of the solar system, taking into account tides, librations and relativistic effects. Over the years, the different versions are known with the acronym DE (Development Ephemeris) followed by a number. The HM Nautical Almanac Office together with the US Nautical Almanac Office regularly publish ephemerides concerning the positions of Sun, Moon, planets, satellites, asteroids as well as the main astronomical phenomena on the basis of ephemerides computations made in collaboration with JPL. Numerical ephemerides are also computed at the Institute of Applied Astronomy of the Russian Academy of Sciences (IAA-RAS) on the basis of a dynamical model similar to the one used at JPL and with computational accuracy comparable to that of JPL. Their ephemerides are known by the acronym EPM (Ephemerides of Planets and the Moon). Analytical theories of the planets and the Moon have been developed at the Bureau des Longitudes and at IMCCE (Institut de mécanique céleste et de calcul des éphémérides ) in Paris. The VSOP (Variations Séculaires de l'Orbite Planétaires) ephemerides include the effect of solar oblateness and relativistic perturbations. However, over long time scales the discrepancy with the numerical ephemerides is quite relevant. With this motivation new planetary ephemerides were developed at IMCCE, known with the akronym INPOP (Intégration Numérique Planétaire de l'Observatoire de Paris). They are based on the numerical integration of the equations of motion of the planets, including the Earth's orientation and the Moon's rotation.

Choreographies of the \(N\)-body problem

Remarkable new solutions of the \(N\)-body problem have been found in [CM], [Mo], in which all masses move on the same curve without experiencing collisions. Such solutions are named choreographies; the first one was discovered by J.-L. Lagrange in 1772 for \(N=3\) and it corresponds to the solution in which the 3 bodies form an equilateral triangle in the rotating frame. Always in the case \(N=3\), a figure eight solution was established in [CM], [Mo]. A very large number of (Lissajous and not especially Lissajous-like) choreography solutions was found in [CGMS] by merging variational techniques, symmetry methods and numerical computations. Massive computations have been performed by C. Simò, who found hundreds of simple Newtonian choreography solutions ([S], see also the animations in [Sl]). The fexistence of choreographies in the \(N\)-body problem has been shown in [KS] through the implementation of a computer-assisted proof. In particular, non-symmetric choreographies for 6 and 7 bodies have been shown to exist. The method allows also to study the linear stability of the figure eight solution, restricted to the plane.

References

[A] V.I. Arnold, Small denominators and problems of stability of motion in classical and celestial mechanics, Uspehi Mat. Nauk, 18, no. 6 (114), 91-192 (1963)

[C] A. Celletti, Analysis of resonances in the spin-orbit problem in Celestial Mechanics: The synchronous resonance (Part I), Journal of Applied Mathematics and Physics (ZAMP), 41, 174-204 (1990)

[CC] A. Celletti, L. Chierchia, KAM Stability and Celestial Mechanics, Memoirs of the American Mathematical Society, 187, 878 (2007)

[CGMS] A. Chenciner, J. Gerver, R. Montgomery, C. Simò, Simple Choreographic Motions of \(N\) bodies: A preliminary study, in Geometry, Mechanics, and Dynamics, volume in honor of the 60th birthday of J.E. Marsden, P. Newton, P. Holmes, A. Weinstein, ed. , Springer-Verlag (2002)

[CM] A. Chenciner, R. Montgomery, A remarkable periodic solution of the three body problem in the case of equal masses, Annals of Math. 152, 881-901 (2000)

[FM] S. Ferraz-Mello, Canonical Perturbation Theories, Springer-Verlag, Berlin, Heidelberg, New York (2007)

[H] M. Hénon, Explorationes numérique du probléme restreint IV: Masses egales, orbites non periodique, Bullettin Astronomique, 3, 1, fasc. 2, 49-66 (1966)

[KS] T. Kapela, C. Simò, Computer assisted proofs for nonsymmetric planar choreographies and for stability of the Eight, Nonlinearity 20, 1241-1255 (2007)

[K] A.N. Kolmogorov, On the conservation of conditionally periodic motions under small perturbation of the Hamiltonian, Dokl. Akad. Nauk. SSR, 98, 527-530 (1954)

[L] J. Laskar, Large scale chaos and marginal stability in the solar system, Cel. Mech. Dyn. Astr., 64, 115-162 (1996)

[LR] J. Laskar, P. Robutel, Stability of the planetary three-body problem. I. Expansion of the planetary Hamiltonian, Celestial Mechanics and Dynamical Astronomy, 62, no. 3, 193-217 (1995)

[LG] U. Locatelli, A. Giorgilli, Invariant tori in the secular motions of the three-body planetary systems, Celestial Mechanics and Dynamical Astronomy, 78, 47-74 (2000)

[M] J. Moser, On invariant curves of area-preserving mappings of an annulus, Nach. Akad. Wiss. Gottingen, Math. Phys. Kl. II, 1, 1-20, (1962)

[Mo] C. Moore, Braids in Classical Gravity, Physical Review Letters, 70, 3675–3679 (1993)

[N] N.N. Nekhoroshev, An exponential estimate of the time of stability of nearly integrable Hamiltonian systems, Russian Math. Surv., 32, 1-65, (1977)

[Roy] A.E. Roy, Orbital motion, Adam Hilger (1991)

[S] C. Simò, New families of Solutions in $N$-Body Problems, in Proceedings of the 3rd European Congress of Mathematics, ed. C. Casacuberta, R. M. Miró-Roig, J. Verdera, S. Xambò, Progress in Mathematics series, Vol 201, 101-115, Birkäuser, Basel (2001)

[Sl] C. Simò, http://www.maia.ub.es/dsg/nbody.html

[SW] G.J. Sussman, J. Wisdom, Chaotic evolution of the solar system, Science, 241, 56-62 (1992)

Internal references

• Alain Chenciner (2007), Three Body Problem, Scholarpedia, 2(10):2111.
• Luigi Chierchia and John N. Mather (2008), Kolmogorov-Arnold-Moser Theory, Scholarpedia, 5(9):2123.
• Sylvio Ferraz-Mello (2008), Celestial Mechanics, Scholarpedia, 4(1):4416.
• Jacques Laskar (2008), Stability of the solar system, Scholarpedia.
• James Meiss (2007), Hamiltonian Systems, Scholarpedia, 2(8):1943.
• Richard Montgomery (2010), N-body choreographies, Scholarpedia, 5(11):10666.
• Michele Trenti and Piet Hut (2008), N-body simulations (gravitational), Scholarpedia, 3(5):3930.
• Ferdinand Verhulst (2007), Hamiltonian Normal Forms, Scholarpedia, 2(8):2101.

See also

Celestial Mechanics, Hamiltonian Normal Forms, Hamiltonian Systems, Kolmogorov-Arnold-Moser Theory, N-body simulations (gravitational), Stability of the solar system, Three Body Problem