Stability of the triangular Lagrangian points

Stability chart of the triangular points in the elliptic restricted problem of three bodies

T. Kovács E-mail:
Accepted – Received –

The possible observations of Trojan-like extrasolar planets stimulate the deeper understanding of the stability behaviour of the co-orbital resonant motion. By using Hill’s equations and the energy-rate method an analysis of the stability map of the elliptic restricted three-body problem is performed. Regions of the parameter plane are described numerically and related to the resonant frequencies of librational motion. Stability and instability can therefore be obtained by analysing the two independent frequency modes depending on system parameters. The key role of the long period libration in determining the structure of the stability is demonstrated and also a stability mechanism is found that is responsible for extended life time of the test particle in the unstable domain of the stability map.

celestial mechanics – planets and satellites: dynamical evolution and stability – methods: numerical, analytical
pagerange: Stability chart of the triangular points in the elliptic restricted problem of three bodiesBpubyear: 2002

1 Introduction

The stability of the Lagrangian triangular points and in the elliptic restricted three-body problem (ERTBP) is an extensively studied topic of celestial mechanics. Recent directions in extrasolar planetary dynamics require to rethink our knowledge about the historical results of pure dynamical astronomy. The stability of motion around and depends on two parameters in the ERTBP, namely on the orbital eccentricity () and the mass parameter () of the primaries.

The main results on the stability of and were published in Danby (1964), Bennett (1965), and Tschauner (1971). Numerical studies about the stability was carried out by Lohinger & Dvorak (1993) and Markellos et al. (1995). It was shown by Érdi et al. (2007, 2009) and Rajnai et al. (2010) that there are secondary resonances responsible for the fine structure of the stability chart. Recently Schwarz et al. (2012) demonstrated the stability phenomena in the spatial problem, too.

In this work we present a fast and accurate method to construct the detailed stability map of the problem using Hill’s equations. The results reveal the complex structure of the plane and give evidence for peaks of longer escape times in the unstable region of the parameter space.

2 Equations of motion and methods

The behaviour of the test particle near the and points in the planar elliptic restricted three-body problem can be described by the well-known equations of variation (Szebehely, 1967)



In Eqs. (1) the prime denotes the derivative with respect to the true anomaly The orbital parameter and the dimensionless mass parameter are involved in the problem, and are the masses of the primaries.

2.1 Hill’s equations

By a suitable transformation (Tschauner, 1971), the 4-th order system (1) can be separated into two independent 2-nd order systems. From these we can derive an equivalent set of Hill’s equations (Meire, 1980, 1982):


where are -periodic functions containing the system parameters. Therefore the system (2) can be thought as of two decoupled harmonic oscillators with variable frequencies. The characteristic units in Eqs. (1) and (2) were chosen so that the primaries’ period is i.e. the mean motion is 1. For a detailed description of Eq. (2) see Appendix A.

2.2 Energy-rate method

The energy-rate method (ERM) (Jazar, 2004) gives the opportunity, beside the well-known methods (Floquet theory, Poincaré-Linstedt series), to investigate the stability boundary of a system




and is a single variable function, and is a nonlinear periodic function of the time The functions and may depend on some parameters of the system. It is trivial from Appendix A that Eqs. (2) can be written in the form of Eq. (3)




corresponds to excitation in Eq. (3) as one may check in Eqs. (10).

As is well-known if any non-conservative force is present in the system, the mechanical energy does not remain constant. In other words, due to the effect of in Eq. (3) for a given set of parameters, energy is either pumping into or withdrawing from the system. It can be verified (Jazar, 2004) that in a system such as (3) the energy integral has the following form


where is the kinetic and is the potential energy for free motion.

It is also known that the time derivative of the mechanical energy must be zero over one period for conservative and autonomous systems in a steady state periodic cycle. In order to find the parameters that belong to stable or unstable solutions we may evaluate the energy over a period


where is the period of the external force If provided that is a restoring force, then energy is inserted into the system, and Eq. (3)) is unstable. But, if energy is subtracted and the system is stable. On the common boundary

3 Stability analysis

In the following the stability chart of the and points is studied using the ERM. A common way in stability analysis is to identify the characteristic roots of the problem that characterize the dynamically different domains of the parameter plane. The ERM allows us, if the equations are in a suitable format, a simple and fast method to distinguish the sets of parameters corresponding to stable or unstable dynamics.

In order to get the classical stability map of the plane one has to compute Eq. (8) numerically for both Eqs. (5)


Note that equations (2) have a repulsive force, (see (5)) responsible for free motion. Therefore the stable solution for some pairs can be realized when the averaged energy Consequently, unstable motion appears in the opposite case,

Since Eqs. (2) are linear ODEs with periodic coefficients, the choice of initial conditions is arbitrary, i.e. they do not affect the stability of the solution. Throughout this work initial conditions are used. The control parameter that measures the accuracy is set to

Figure 1: Stability map in the plane. Top: different features of motion around and governed by Eq. (2a). Bottom: stability domains corresponding to Eq. (2b). The contour bar indicates the magnitude of the average energy for given pairs of the parameters The integration time is one period of the primaries. The number of the parameter pairs is set to i.e. 1000 in each direction with equal step size. (Colour figures are available online.)

Figure 1(a) depicts the classical stability map of the ERTBP obtained from Eq. (2a) by the ERM. The light gray region (orange online), where for one period, corresponds to librational motion around and The dark gray area (blue) represents the negative energy domain related to the set of parameters of unstable motion. The black part of the figure (dark blue), situated at the bottom right quarter, belongs to parameters where Eqs. (2) have no solutions. For better visualization (log-scale) we cut the positive energies at and set the negative values to additionally, where the problem has no solution to Therefore, the contour bar is organized as follows: the light gray colour () represents libration, the dark gray () denotes unstable motion, and the black domain () is for complex frequencies. The border of the stability domain is at The charecteristic roots of Eqs. (1) (Bennett, 1965) devide the stability chart exactly into the same parts.

Figure 1(b) shows the stability chart provided by Eq. (2b). It can be seen that the librational domain (light gray (orange online)) is more extended than in panel (a), however, the bottom right part (black/dark blue) is exactly the same. If we want to identify the stability region around the and points in the ERTBP determined by Hill’s equations, we need the intersection of all the individual stable regions in Fig. 1(a) and 1(b). One can immediately verify that the upper panel encompasses the domain in question. In other words, this shows that merely the first Hill’s equation in Eqs. (2) is responsible for the well-known structure of the plane. The possible reason will be discussed in Section 5. Accordingly, we can state that the stability chart obtained from Eq. (2) is consistent with earlier analysis.

One of many advantages of the ERM is that if we integrate the equation(s) of motion until where is an integer number and is the period of the external force, new, so-called ”splitting” curves appear on the stability map. (The name comes from the fact that these resonant periodic lines split the stable regions into several parts.) Actually, these curves are super- and sub-harmonic oscillations and describe higher order periodic responses of the parametric equations.

Figure 2 shows the average energy levels at i.e. for 4 revolutions of the primaries. The similarity with Fig. (1) is obvious but we have several new splitting curves inside the stable regions, both in panel (a) and (b). In Section 4 it will be shown that these resonances correspond to the type and secondary resonances defined in Érdi et al. (2009).

Figure 2: The plane for integration time Top and bottom: stability zones derived from Eqs. (2a) and (2b), respectively. Increasing the computation time by the multiplication of the excitation period, resonant curves turn up inside the stable regions.

4 Resonant frequencies

It was shown by Érdi et al. (2009) that the fine structure of the plane is organized by secondary resonances between the four frequencies in the ERTBP. Two of them, the short and long frequencies (corresponding to the short and long periodic librations around and ), are present also in the circular problem, the additional two (in normalized units: ) appear when the primaries revolve on elliptic orbits.

In general case, Hill’s equations (2) can be considered as of two independent harmonic oscillators with periodic coefficients, i.e. with variable frequencies. Eqs (2) describe the two fundamental frequencies in ERTBP, namely the frequencies corresponding to the short and long periodic components of the libration around and When the problem is circular, Eqs. (2) become simple harmonic oscillators and serve the well-known frequency diagram of the problem. (See Appendix A, Fig. 6.)

The ERM applied to the ERTBP can be used to identify the secondary resonances between the librational frequencies of the test particle and the primaries’ motion. These are the resonant curves on the stability chart defined as A-type and D-type in Érdi et al. (2009). Investigating the problem through Eqs. (2) one do not need to involve the two ”additional” normalized frequencies coming from ellipticity. Let us consider a quick example. The A11 resonance in Érdi et al. (2009) corresponds to the nearly V-shape boundary line of the stability zone. The notation A11 comes from However, after a simple rearrangement we get i.e. the frequency is equal to the half of the normalized frequency of the primaries. In other words, the long period of the librating infinitesimal particle along the A11 resonance curve is twice as the primaries’ period (Danby, 1964). Moreover, this rule holds for any A and D-type resonance. If we integrate Eqs. (2) for and plot the sets of parameters corresponding to condition, the stability map with secondary resonances, of types A and D, can be drawn from.

In Figure 3 the plane is depicted for various integration times. Panel (a) contains the ”original” stability chart of librational motion in the ERTBP with agreement of the energy plot of Fig. 1. Dark blue curves (color online, hereafter dark curves) correspond to the long period motion () while the light blue (color online, hereafter light curves) represent the short period (). We note again, these resonant lines indicate the two components of motion around the points and Accordingly, the nearly V-shaped dark curves describe the 2:1 commensurability of the test particle’s motion, in direction with the primaries’ period. The light line represents the resonant component in the other direction Along these curves the average energy during one period is zero. This picture itself does not say anything about the stability of the two-dimensional motion but only the resonant frequencies. We also need the characteristics of Figs. 1 and 2 to decide whether the motion is stable or not. Fig. 3(a) tells us only the border of the stability region. It has already been shown in Section 3 that libration can occur under the dark curve, where both components of Eqs. (2) are stable.

Other secondary resonances show up on the stability map when the integration time is doubled, as in Fig. 3(b). In this panel two extra lines appear, beside the well-known stability border, corresponding to the 4:1 mean motion resonances . The plot was obtained as follows: the integration time was two periods of the primaries, the pairs were stored when the energy rate was zero during this time. Now, the harmonic resonant curves are the new lines and the old borders are the sub-harmonics of them (2:1 res.). This behaviour is defined as splitting in the ERM, where resonant curves split the stability regions into two parts (Jazar, 2004). The 4:1 resonant curve originating from corresponds to A31 resonance in Érdi et al. (2009). To demonstrate the feature of the resonances Fig. 4 shows two individual solutions of A-type resonances, 2:1 and 4:1.

If one continues the procedure, more and more splitting curves turn up on the stability map. In panel (c) (, 3 revolutions of the primaries) the fundamental resonances are 6:1 – (A51 in Érdi et al. (2009)) and for type D resonance. The sub-harmonics describe the 3:1 (A21) and 2:1 (A11) resonances.

In addition two interesting curves appear in panel (c). Curve a originates from curve d comes from the stability region of and terminates where curve a ends, about at Integrating Eqs. (2) along these curves, one can find that these librations correspond to the super-harmonics of the resonances (6:1). These, in question, are the 2:3 commensurabilities. The role of the meeting point is discussed in Appendix B.

Figure 3(d) demonstrates the resonant frequencies for of the primaries As one may expect various new resonances appeared on the plane. Table 1 collects the resonances corresponding to dark curves on panel (d). D-type resonances are light curves of similar types.

Resonance starting point panel(s)
8:1 0.00228
4:1 0.00876 b
8:3 0.01823
2:1 0.02859 a,b,c
8:5 0.03660

Second column contains the mass parameters along the axis where the resonant lines originate. Third column gives the other panels where the same resonance appears.

Table 1: A-type resonances in Figure 3(d).

To show how complex the stability chart can be, the plane for is plotted in Figure 5. The 10:1 resonances and all their sub- and super-harmonics are depicted. An even more complete and dense structure could be obtained by overplotting the charts corresponding to different integration times.

Figure 3: Resonant frequencies at different integration times. Computations were performed from top to bottom for times: In panel (b)-(d) splitting curves appear. For detailed explanation see the text.

Figure 4: Two individual solutions of Eq. (2a). The dashed line corresponds to the 2:1 resonant frequency for the parameter pair see the circle on the dark curve in Fig. 3(b). While the slolid curve denotes the 4:1 resonance for , this latter point is marked with a cross in the stability chart, see also Fig. 3(b). The dashed and dotted line represents t=25.12 (8), i.e. 4 revolutions of the primaries.

Figure 5: Stability chart for This plot contains the 10:1 resonant curves and all sub- and super-harmonics of them.

5 Summary and conclusions

In the present work the stability analysis of the ERTBP was performed. Using Hill’s equations we showed that the problem can be described as two independent harmonic oscillators with time dependent frequencies. It was also demonstrated that using this approach, combined with the ERM, one can find the same stability chart as many authors published previously.

Beside many advantages, the energy-rate method allows us to obtain the resonant curves on the plane with arbitrary accuracy (obviously depending on numerical truncation errors). Moreover, it is not restricted to moderate values of the parameters as common perturbation technics.

Investigating the librational motion around the and points in the ERTBP by two independent Hill equations, two unresolved problems in stability analysis can be interpreted

  1. Figs. 1 and 2 indicate that the long period component (frequencies from Eq. (2a)) of libration is responsible for the global feature of stability. Indeed, the A-type resonances are situated in the intersection part of the two stable regions (long and short period librations), they are therefore the backbone of the two-dimensional stable motion.

  2. Rajnai et al. (2010) showed that there are significant peaks in the unstable domain responsible for longer escape times. In other words, there may exist some processes that stretches the lifetime of the test particle around and Comparing their result (Rajnai et al. (2010), Fig. 1) with Figure 3, one can see that the longer escape times appear exactly where the resonant frequency curves, corresponding to the shorter frequencies, are situated. An explanation of this behaviour, based on the present work, can be that the resonant curves in the stable domain of that are actually invisible on the classical stability chart, may partially stabilize the whole librational motion around the triangular Lagrangian points.

A direct consequence of this study can be, in future, to extend the investigation to the spatial problem, considering the third equation of motion. Interestingly, the third variational equation beside the and components of (1) is already in form of Hill’s equation, therefore, only a simple variable change is needed to study the extended ERTBP.

The reader may ask why are the B-type resonances (Érdi et al., 2009) omitted from the present study. Since the ERM cannot be used to identify the exact solution of Hill’s equations (2) and therefore one cannot compare the individual frequencies, (perhaps the most important) resonances were not examined in this work thoroughly. However, an interesting analytical approach of B11 resonance is presented in Appendix B.


The author wishes to thank Prof. B. Érdi for helpful discussions. This work was supported by the Hungarian OTKA Grant project number K-81373.


  • Arnold (1978) Arnold V.I., 1978, Mathematical methods of classical mechanics. Springer, New York
  • Bennett (1965) Bennett A., 1965, Icarus, 4, 177
  • Danby (1964) Danby J.M.A., 1964, AJ, 69, 165
  • Érdi et al. (2007) Érdi B., Nagy, I., Sándor, Zs., Süli, Á., Fröhlich, G., 2007, MNRAS, 381, 33
  • Érdi et al. (2009) Érdi B., Forgács-Dajka E., Nagy I., Rajnai R. 2009, Celest. Mech. and Dyn. Astron., 104, 145
  • Jazar (2004) Jazar G.N., 2004, International Journal of Non-Linear Mechanics, 39, 1319
  • Lohinger & Dvorak (1993) Lohinger E., Dvorak R., 1993, A&A, 280, 683
  • Markellos et al. (1995) Markellos V.V., Papadakis K.E., Perdios E.A., 1995. in Roy A.E., Steves B., ed., From Newton to Chaos, Plenum Press, New York, p.371
  • Meire (1980) Meire M., 1980, Bull. Astron. Inst. Czechosl., 31, 312
  • Meire (1982) Meire M., 1982, A&A, 110, 152
  • Rajnai et al. (2010) Rajnai R., Nagy I., Érdi B., 2010, Journal of Physics: Conference Series, 218, 012018
  • Schwarz et al. (2012) Schwarz R., Bazsó Á., Funk B., Érdi B., 2012, MNRAS, 427, 397
  • Szebehely (1967) Szebehely V., 1967, Theory of orbits. Academic Press, New York
  • Tschauner (1971) Tschauner J., 1971, Celest. Mech. and Dyn. Astron., 3, 189
  • Yoshida (1984) Yoshida H., 1984, Celest. Mech. and Dyn. Astron., 32, 73

Appendix A Hill equations in the RTBP

Following Meire’s (Meire, 1982) derivation, the final form of Hill’s equations in the restricted three-body problem are




In Eqs. (10) denotes the eccentricity, is the true anomaly, and are the masses of the primaries, and is the mutual distance between them.

One can obtain the motion around and as a composition of two harmonic oscillations governed by Eqs. (10). Obviously, the form of equations shows that the frequencies of the motions are However, the solution is trivial only in the case when the frequencies do not depend on time, i.e. the orbital parameter is zero. After some computations, one can easily verify that in this case () the frequencies are



Figure 6: Librational frequencies in the circular RTBP. The dark blue (online) curve corresponds to the long period solution or The light blue (online) curve depicts the short period component or Eqs. (11) provide the exact frequencies depending only on the mass parameter as .

Figure 6 demonstrates and vs. the mass parameter It can be seen that the curves exactly match the well-known classical results of the circular RTBP. The lower frequency represents the long period solution while correspondes to short time oscillations. The latter tends to the frequency of the primaries (here normalized to 1) in the limit case when

Appendix B The B11 resonance

In the elliptic case the frequencies of libration Eqs. (10) are not so obvious. In this case we have second order linear differential equations with time dependent coefficients. For an arbitrary given periodic coefficient, there are no general procedures to express the characteristic exponents, unless the fundamental matrix of the problem is obtained by quadratures (Arnold, 1978; Yoshida, 1984). Therefore, it seems to be a hard task to find any analytical expression for the period, and therefore the frequency of the solution

Figure 7: The frequency curve (B11) as the border of stability. A simple averaging by using Eqs. (12) gives the exact curve of the B11 resonance. This curve indicates the parameter pairs where The resonant curves A11, D11, and B11 meet at one point P. Additionally, several resonant points can be identified for different mass parameters at . Crosses along the -axis denote the following B type resonances, from left to right: 5:1, 4:1, 3:1, 2:1, 3:2, 4:3, 1:1, respectively.

Let us consider, however, a naive analytical approximation to find the resonant curves corresponding to the B-type resonances. The B-type resonances are defined as: where and denote the short and long period frequencies, respectively (Érdi et al., 2009). These frequencies, that are actually varying with time, can be easily obtained by evaluating the coefficients in Eqs. (10). Suppose one can describe the solution with frequency changing around the average frequency


where is the period of is a stepsize, and is the number of terms during one period. Consider these averaged frequencies the skeleton of individual librations around and In this sense, B-type resonances can be evaluated simply analytically by using the formulae given by Eqs. (12).

Just for the interesting coincidence with Érdi’s work, who used Fourier analysis to find the exact period, the B11, i.e. resonance is plotted in Figure 7. As is well-seen the curve originates at and follows perfectly the border of stability, touching the critical point and draw out the border line, seen also in Fig. 1, that separates two different kinds of solutions (Eqs. (10)). No analytical solution in linear theory was found so far that can describe this behaviour (Szebehely, 1967).

Figure 8: Various A and D-type resonances near to B11 resonance. The panels (a)-(c) refer to different integration times, (Magnified parts of Figs. 3(c), 3(d), and 5, respectively.) On each panel the touching point is marked with a circle. At these points the two librational frequencies and are equal. Panel (d) verifies that the B11 resonance obtained by using Eqs. (12) contains the circles.

At this point we should explain the role of the pair of the A11 and D11 (P in Fig. 7) resonant curves that touch each other in a common point. These A and D-type resonances describe a well-defined secondary resonance between the librating particle and the primaries’ motion. However, at that particular point where they meet, the two individual libration frequencies become equal. Therefore, the motion at this point also satisfies the condition of B11 resonance, i.e. For different integration times new A and D-type resonant curves appear on the plane that meet in one common point, see Fig. 8. These common points are necessarily situated along the B11 resonance curve due to the condition

Finally, we also note that this method provides the right places of the secondary resonances for the circular RTBP. That is, all the resonant frequencies (with A,B,C,D,E,F-types, (Érdi et al., 2009)) can be identified by the analytical computation Eqs. (11) along the -axis, see Fig. 7. The reason of this is that the coefficients in Eqs. (10) no longer depend on time in the case consequently describe the exactly constant frequencies of the librational components, since in the circular case

Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description