Three-Body System: Stable and Chaotic Orbits

Three bodies moving in the closed orbits, exposed to the mutual gravitational interactions only, enter gradually, as a rule, into the gravitational resonances: periods of their rotations become related as the rational fractions.

In long term, this phenomenon may lead to either and for the most part, stabilization of the orbits or, for some period ratios, be the cause of the chaotic motions, also.

Investigation of that phenomenon in this work was based on the, so cold, Newton’s three-body system consisting of one massive and 2 considerably smaller bodies (m_0 \gg m_1, m_2) turning around it. The same procedure may be used for examination of the correspondent restricted three-body problem (m_0 \gg m_1, m_2 \approx 0).

Numerical analysis of the problem was carried out in the following way. After discretizing the time sub-space, the differential equations of motions were replaced by the corresponding difference equations. Their approximate solutions were procured by use of the perturbation technique.

It was revealed that the chaotic motions zones are situated around the T_2:T_1 ~ 1:3 and 3:1 resonances. This outcome is in accordance with the Wisdom’s result for the Sun-Jupiter-asteroid triad /9/, as well as with the observational data. Extensions of the chaotic motions zone depend on the mass ratios \frac {m_1} {m_0} and \frac {m_2} {m_1}.

Neglecting all the other gravitational influences, the adopted model may be applied to the Sun, together with any pair of the heavenly bodies, provided that they do not form the binary system. Also, it is not necessary that the small bodies form the mutually closest pair.

Thus, it was possible to point out, roughly, some of the Solar System regions in which the orbits may become chaotic. Orbital eccentricities of numerous meteoroids, comets and asteroids in these zones steadily increase and their paths began to cross the solar system planets’ orbits.

Regarding the planets, it seems that Venus, Mars, Saturn and Uranus orbits lie in the unstable orbits zones of the Sun-Venus-Mars and Sun-Saturn-Uranus triads, respectively. In the long term, this fact may produce some changes in the orbit geometry of these heavenly bodies, so that they may enter, eventually, into the stable orbit regions.


We consider one large and two small bodies (m_0 \gg m_1, m_2) exposed to no influence other than that of their mutual gravitation. They are moving in the bounded orbits.

Aiming simplification of the calculus, we apply non-dimensionalisation and remove fundamental units by suitable scaling of these quantities. We shall assume


Bar over the symbol denotes the real quantity and G is the gravitational constant. \bar{R}_1(0) and \bar{R}_2(0) denote initial positions of the small bodies with respect to the massive one.

From now on, all kinematical and dynamical quantities (without bar over them) will be dimensionless. In particular


In the perturbation method we use in this work, the number \varepsilon \ll 1 is to be taken for the small parameter.

Since the small bodies m_1 and m_2 perform relative rotations around m_0, we introduce three kinds of reference frames lying in the Laplace’s plane. The first one Cxy is the inertial frame of reference, the second represents transport r.f. 0x’y’ (0x’, 0y’ parallel with 0x and 0y), while the polar coordinates R_1,\varphi_1 and R_2,\varphi_2 defining positions of the small bodies with respect to the large one, are the relative reference frames (Fig. 1).

Reference Frames

Fig.1 – Reference Frames

Gravitational Forces and Harmonization of Motions

Resulting gravitational forces acting on the bodies are represented in the Fig.2.

Gravitational Interactions

Fig.2 – Gravitational Interactions

As known, in the three-body case, these forces meet at one point called center of the gravitational attraction /2, 6/. Here follow their intensities:


Evidently, all dynamical and, therefore, corresponding kinematical characteristics of the system will have extremes when three bodies become aligned,


that is, when the small bodies are in conjunctions (n = 0,2,4,…), or in oppositions (n = 1,3,5,…) with respect to the large one. Again, we underline the fact that the small body’s binary systems moving in the 1:1 resonance that cannot take both aligned positions with respect to the massive body had to be excluded from the model and, consequentially, our solution does not contain the capture into binary system case.

If we have on our disposition, the small bodies conjunctional and oppositional angular velocities \omega_k^c and \omega_k^0, k = 1,2, then, assuming linear change of these variables between two alignments in the time interval h, we may put


where t04-for07represent the mean angular velocities. Now, expression (1) will have the form


Taking into account that


where T_k is the orbital period of the body k, one now obtains


Since the sum, difference, product and quotient of two non-zero real numbers from which one is rational and the other irrational must be irrational and the right hand side of this equation is the rational number, it comes out that both ratios on the left hand side have to be eider rational, eider irrational.

Is it possible to prove mathematically anyone of these assertions? For now, it is hard to say, but the observational evidence of the Newtonian three body systems in our surrounding space indicates that corresponding periods are related as rational numbers. As said, the term resonance is used in celestial mechanics to denote this phenomenon.

It is well-known fact that relation between Earth – Venus periods is 13:8, Saturn – Jupiter is 5:2, Pluto – Neptune 3:2 and so on, not to mention the numerous solar system moons orbital resonance.

Combining the perturbation technique and method of the successive approximations, motion of the same three-body model was investigated in the work /6/. As a rule, the outcome, for the series of arbitrary chosen initial conditions, was resonant motion of the bodies.

Our assumption is that the cause of resonance is alternation of the small bodies’ aligned positions with respect to the massive one. Linear constellations of the bodies produce extremes in the gravitational field and consequentially, extremes in the related kinematical quantities. These field peaks and bottoms require adjustment of the mean orbital velocities in such a way that orbital periods of the bodies in the time interval h between two alignments become proportional.

Capture Position and Initial Conditions

If, given the masses of the bodies, one adopts initial positions and velocities arbitrary, the outcome may be either motions with one, or three open orbits, collision of two, or three bodies, either motions in closed orbits.

It makes no sense to study the three – body system allowing possibility of its decomposition, so we are going to consider only these initial conditions producing bounded orbits, without presumption that these orbits are going to be stable.

Initial Constellation, Conjunctions and Oppositions

Fig.3 – Initial Constellation, Conjunctions and Oppositions

In order to simplify the initial condition problem we shall take that the bodies are aligned in the moment t = 0. We choose conjunction of the small bodies with respect to the large one for the initial („capture“) constellation (Fig.3).

It was taken that Cx axis of the inertial reference frame xCy (C is the common center of masses) lies along the aligned bodies.

Distances from the massive body R_k(0), k = 1,2 can be chosen arbitrary. On the other hand, in order that the orbits become bounded, the velocities \vec{v}_k (0) have to be orthogonal to the Cx axis, while their intensities have to be smaller than the escape (parabolic /8/) velocities.

Aiming to obtain the orbits, which are close to the regular paths of m_1 and m_2, we shall take that in the assumed capture constellation motions of the small bodies are dynamically balanced with the massive body attraction force, only. Therefore, initial angular velocities of these bodies correspond to their circular velocities. Consequentially,


Variables and their Expansion into the Power Series of epsilon

It’s convenient to adopt


to be variables in the exposed problem. They can be expanded into the power series of small number \varepsilon, representing the body 1 and body 0 mass relation:


We shall take that the leading terms R_{k0} and \omega_{k0} are constants while the other terms are homogeneous functions of time


Equations of Motions

Let’s represent differential equations of the small bodies motions in the inertial reference frame, as follows


From now on, we consider nothing but the linear constellations of the bodies, that is, their positions for t = nh (n = 0, 1, 2, 3,…). Multiplying equations (4) by the unit vectors \vec{e}_{Rk}, one obtains the equations of motions


The minus sign corresponds to the body 2, in opposition (Figures 4 and 5). If we multiply the same equation by \vec{e}_{\varphi k} we come to the alignment conditions


and to the conclusion that, in fact, these conditions represent differential form of the angular momentum conservation law


Now, we can introduce perturbation series (3) into the equations (6)


and obtain the alignment conditions arranged by order of \varepsilon: \varepsilon^0, \varepsilon^1, \varepsilon^2, …


We can take that the leading terms of the series (2) are equal to the initial positions and initial angular velocities of the bodies 1 and 2. In fact, these terms represent the kinematic quantities corresponding to the two-bodies systems 0, 1 and 0, 2. Consequentially the sums of the following terms will contain the small bodies’ mutual „perturbations“.

In accordeance with the initial condition (2), we may now write down


Absolute Accelerations

Now, we shall write down absolute accelerations of the bodies and, using perturbation series (3), arrange their components by order 0, 1 and 2 of the small parameter \varepsilon.




Fig.4 – Conjunction





Fig. 5 – Opposition


Differential Component of the Radial Acceleration

Approximately, the second derivative of the variable R with respect to time may be written in the finite-difference form




and h, as said, represents the time interval between two alignments.

Taking into account the fact that


we can also write this derivative in the form


The conjunction – opposition time interval h is equal


By use of the perturbation series (3), \ddot{R}_k obtains the following form






Here follow the terms of the given expression classified by order 0, 1 and 2 of the parameter \varepsilon:


Centripetal Acceleration

Let’s now introduce the perturbation series (3) into the centripetal component of the radial acceleration and arrange its terms by order of the small parameter.


Using (6a), (6b) and (6c) one obtains



We have assumed conjunction as the capture moment and the dynamical balance between the massive and small bodies, neglecting the small bodies’ interactions. Thus, the initial conditions for both bodies correspond to the perturbations of the order (\varepsilon^0):

t = 0 (conjuction):


t = nh (n = 1,3…) (oppositions):


t = nh (n =0, 2, 4….) (conjunctions):




Here follows procedure of the calculus.

Introducing the initial conditions (8) one finds the solutions R_{kj}(h), j > 0 of the equations (9). Then, by use of the perturbation series (3a), he comes to the positions of the bodies in opposition R_1(h) and R_2(h). The obtained values R_{k1} and R_{k2} are input in (10) and solutions of these equations R_1(2h) and R_2(2h) go again in (9) and so on (we emphasize the fact that, theoretically, number of iterations – that is, length of the motion simulation – is unlimited) until the stabilized or chaotic solutions appear.

Stable and Chaotic Orbits

The coefficients A_k, C_{kj}, O_{kj} depend on the \frac {R_{20}} {R_{10}} ratio. In other words, since non-dimensional R_{10} = 1, on R_{20}, only.

Equations (9) and (10) will always have the solutions, excepting the case when A_k(R_{20}) = 0. Simply, the existing solutions correspond to the stable, while the solutions corresponding to the coefficients A_k \to 0 go together with the chaotic orbits.

Thus, according to our calculus, the body 1 will have chaotic orbits around the resonance \frac {T_2} {T_1} = \frac {1} {1+\pi/\sqrt{2}} = 0,31 \approx \frac {1} {3}, while the orbits of the body 2 are going to be chaotic around \frac {T_2} {T_1} = \frac {1+\pi/\sqrt{2}} {1} = 3,22 \approx \frac {3} {1}. The corresponding initial positions R_{20} of the body 2 are 0,46 (~ 0,48) and 2,18 (~ 2,08).

It is important to notice that the obtained results do not depend on the power series approximations of the difference equations solutions but only on the preciseness of the finite -difference approximation (7) of the derivative \ddot{R}. The more accurate approximation would approach the result to the exact resonances, producing orthogonality of the conjunction – opposition directions.

The first of these results coincides with the the one derived by Wisdom in the case of the Sun-Jupiter-asteroid system. He was studying the most prominent Main Belt (Kirkwood) gap and obtained the maximal Lyapunov exponent (LE) for the corresponding trajectories near \frac {T_A} {T_J} = \frac {1} {3} commensurability /8/.

Now, the question of convergence, or divergence emerges.

Of course, the calculus possibilities implies assumption of the finite power series order, as well as the limited number of iterations.

Regarding the power series (2a), there is no problem with their convergence toward r_K, provided that \varepsilon (= \frac {m_1} {m_0}) is sufficiently small and that perturbation order is sufficently high.

As long as the bodies move in the stable orbits zones (LE < 0), the obtained results are satisfactorily correct. Since every irrational number may be represented, with the chosen accuracy, by the rational one, in these domains, given R_{10}(=1), R_{20} , the outcome will be the stable and resonant mean motions of two bodies in the orbits neighbouring the initial ones.

However, it was necessary to establish some criterion for the limits of these regions. Rational pressumption is that, as long as

R_{kN} = R_{kN-2} (k=1,2; Nthe number of iterations)

the small bodies move in the stable orbit zones.

On the other hand, our results are, obviously, worthless in the bifurcation (LE ~ 0) and chaotic (LE > 0) domains.

The problem with the bifurcation orbit zones lies in the fact that the orbit’s nature manifestation in these domains requests tens of millions of years and consequently, convergence (or divergence) of the iterational procedure is accordingly slow.

Still, in order to obtain, at least, a rough sketch of this bifurcation zone, we had to establish a criterion for its limits. Therefore, taking into account the adopted calculus precision, we assumed that the body moves in this transitional to the chaos zone as long as

\mu_1R_{kn+2}R_{kn} | \le \mu_2   (k=1,2; n = 0,1,2,…, N),

where \mu_1 and \mu_2 are the small numbers.

In this work, it was chosen that \mu_1 = 0,001 and \mu_2 = 0,01.

The chaotic orbits zone is clearly distinctive because following the limit \mu_2 the results begin to behave markedly erratically. Observational evidence of the chaotic motion in the “Keplerian” orbit is the sudden increase of its eccentricity appearance /8/.

As a conclusion, we may say that these solutions:

  • are valid in the stable orbits zones,
  • they point up the chaotic motions producing resonance,
  • roughly register the limits of stable, bifurcation and chaotic orbit zones, as well as the capture into the binary systems domains.

Extensions of the chaotic orbit and binary systems zones depend mainly on the mass ratios \varepsilon and \lambda. On the other hand, it seems that these ratios do not affect appreciably the orbits in the stable orbit zones.

We shall illustrate the obtained result with one example. Let us take that the parameter
\varepsilon = \frac {m_1} {m_0} is equal 0,001, while the parameter \lambda = \frac {m_2} {m_1} = 0,5.

Functions representing the mean orbits


and the mean resonance


are shown in the Fig. 6. Out of the limited bifurcation and chaotic motion zones and out of the capture into the binary systems domains, orbits are stable.

Stable, Bifurcation and Chaotic Orbit Zones & Resonance Fig. 6 – Stable, Bifurcation and Chaotic Orbit Zones & Resonance

Some of the Solar System Chaotic Domains

Neglecting all the other gravitational influences and assuming that the Solar System bodies move in the same, ecliptic plane, the Sun, together with any couple (not necessarily mutualy closest) heavenly bodies, provided that they do not form the binary system, correspond to our model.

In such a way, combining the Sun, planet and asteroid into the Newton’s restricted three body system, it was possible to point out, roughly, some of the Solar System regions in which the orbits may become chaotic. Orbital eccentricities of the numerous meteoroids, comets and asteroids in these zones steadily increase and their paths began to cross the solar system planets’ orbits.

There is a lot of the small Solar System bodies orbiting the Sun. Most of these bodies inhabit two asteroid belts. The first one is between Mars and Jupiter, the Main Belt, situated on 2 to 3,3 AU distance from the Sun, while the second one, the Kuiper Belt, lies in the region beyond the planets extending from the orbit of Neptune, at 30 AU, to, approximately, 55 AU from the Sun.

The inner and the outer planets bifurcation and chaotic orbit zones are represented in the Figures 7 and 8. Under the line are the bifurcation and chaotic orbit domains of the asteroids moving near 1/3, while over the line are domains of these bodies moving near 3/1 resonance with the corresponding planet.

As said, the most prominent the Main Belt gap lies around the orbits corresponding to 1/3 resonance with Jupiter. On the other hand there is, practically, no small bodies in the Kuiper belt moving around the 3/1 resonance with Uranus.

Regarding the Solar system planets motions stability, the Venus, Mars and the Saturn, Uranus mean periods ratios are neighboring the chaotic orbital resonance 1/3 : T_v / T_{Ma} = 0,327 ; T_{Sat} / T_U = 0,351. Extensions of the Venus and the Saturn unstable motion zones were determined by use of the SM_aV and SUS_a (=0,1,2) triads, while extensions of the Mars and the Uranus unstable motion domains, by use of the SVM_a and SS_aU (= 0,1,2) triads, respectively.

In addition, the planets Jupiter and Saturn move in the bifurcation zones of the corresponding triads. Extensions of the planets unstable orbits zones are represented in the Figs 7 and 8, as well. Concerning the zones where the influences of two planets interleave (the most distinct are these in the Main Belt and the one in which Saturn moves), it is evident that the three body model is not adequate there.

Inner Planets: Bifurcation and Chaotic Orbit Domains

Fig. 7 – Inner Planets: Bifurcation and Chaotic Orbit Domains

Outer Planets: Bifurcation and Chaotic Orbit Domains

Fig. 8 – Outer Planets: Bifurcation and Chaotic Orbit Domains

Does that mean that orbits of these planets will relay become chaotic in the nearby future? Of course not. The time measure in the long term dynamics of the heavenly bodies chaotic motions is, at least, tens of millions of years /8,10/. Our assumption is that in the future, these four planets will, eventually, migrate inward, or outward, into the stable orbits regions and assume harmonized motions with the neighboring heavenly bodies.


Phenomenon of the resonance in the Newton’s three-body system was investigated in this work. As known, resonant motions of the bodies produce, for the most part, stable orbits, but some period ratios are causes of the chaotic motions, also.

It was demonstrated that the chaotic motions zones are situated around the T_1 : T_2 ~ 1:3 and 3:1 resonances and that their extensions mostly depend on the small bodies and the large body mass ratios.

The obtained results were used to examine some of the Solar System regions in which the orbits may become chaotic.

It was revealed that Venus, Mars, Saturn and Uranus orbits lie in the unstable orbits zones of the Sun-Venus-Mars and Sun-Saturn-Uranus systems, respectively. In the long term, this fact may produce some changes in the orbital geometry of these heavenly bodies, so that they may, eventually, enter into the stable orbit regions.

Concerning the small Solar System bodies moving through the interplanetary space along the paths crossing the planets’ orbits, most of them come from the Main and the Kuiper asteroid belts. Namely, parts of the Main Belt lie in the chaotic domains of the Sun-Earth-asteroid, Sun-Mars-asteroid and Sun-Jupiter-asteroid systems, while the central part of the Kuiper Belt is occupied by the chaotic zone of the Sun-Uranus-asteroid system.


1. M. Milanković: Osnovi nebeske mehanike, Naučna knjiga, Beograd, 1988.

2. T. Anđelić: Uvod u astrodinamiku, Matematički institut, Beograd, 1983.

3. Hawking S., Israel W. 1989. 300 Years of Gravitation, Cambridge University Press.

4. Pars L. A. 1956. A Treatise of Analytical Dynamics, Heinemann, London.

5. R. Bellman: Perturbation Techniques in Mathematics, Physics, and Engineering,

6. Marjanov M.: One Solution of the Newton’s Three-Body Problem, 2nd International
Congress of SSM, Palic, 2009.

7. Marjanov M. 2004. Homogeneous and Inhomogeneous Gravitational Fields,
Zbornik radova Gradjevinskog fakulteta u Subotici.

8. Wisdom J.: Chaotic Behavior and the Origin of the 3/1 Kirkwood Gap, Department of
Physic, University of California, Santa Barbara, California 93106

9. Pivko S.: Mehanika III – Dinamika, Naučna Knjiga, Beograd 1974.

10. John D. Hadjidemetriou: On periodic orbits and resonance in extrasolar planetary
systems, Department of Physics, University of Thessaloniki, Thessaloniki, Greece.

11. ORBITAL RESONANCES AND CHAOS IN THE SOLAR SYSTEM, Renu Malhotra Lunar and Planetary Institute, 3600 Bay Ar Blvd, Houston.



Odnosi perioda rotacija tri tela koja se kreću u zatvorenim orbitama isključivo pod dejstvom uzajamnih gravitacionih interakcija postaju vremenom, po pravilu, racionalni brojevi. Kaže se da njihova kretanja ulaze u rezonancu.

Dugoročno, posledice tog fenomena mogu biti usklađena kretanja i stabilne orbite, ali postoje i neki odnosi perioda pri kojim kretanja postaju haotična.

Ispitivanje tih pojava se u ovome radu osniva na tzv. Njutnovom sistemu tri tela: jednog masivnog i dva tela znatno manjih masa (m_0 \gg m_1, m_2), koja kruže oko njega. Postupak se može primeniti i na odgovarajući ograničeni problem tri tela (m_0 \gg m_1, m_2 \approx 0).

Numerička analiza problema u radu je obavljena na sledeći način. Prvo je vremenski podprostor diskretizovan, pa su diferencijalne jednačine kretanja aproksimirane odgovarajućim diferencnim. Potom je za rešavanje tih jednačina primenjena perturbaciona tehnika.

Pokazuje se da su zone haotičnih kretanja malih tela u okolinama rezonanci T_2:T_1 ~ 1:3 i 3:1, što je u skladu s rezultatom Džeka Vizdoma (J. Wisdom) za sistem Sunce-Jupiter-asteroid /9/, kao i sa opservacionim podacima. Širine tih oblasti zavise od odnosa masa \frac {m_1} {m_0} i \frac {m_2} {m_1}.

Zanemare li se svi drugi gravitacioni uticaji, usvojeni model može se primeniti na Sunce i bilo koji par nebeskih tela, uz uslov da ona (u okviru sistema tri tela) ne formiraju binarni sistem. Pri tome, nije neophodno da par malih nebeskih tela bude međusobno najbliži.

To dopušta da se ukaže, doduše grubo, na neke od zona Solarnog sistema u kojima orbite mogu postati haotične. Ekscentriciteti orbita brojnih meteoroida, kometa i asteroida u tim zonama neprekidno rastu, pa njihove putanje počinju da presecaju orbite planeta sunčevog sistema.

Što se tiče planeta, izgleda da se Venera, Mars, Saturn i Uran kreću u zonama nestabilnih orbita trijada Sunce-Venera-Mars, odnosno, Sunce-Saturn-Uran i to bi vremenom moglo dovesti do izvesnih promena geometrija orbita tih nebeskih tela, kako bi se one našle u regionima stabilnih kretanja.

This entry was posted in Celestial Mechanics. Bookmark the permalink.