Stability chart of the triangular points in the elliptic restricted problem of three bodies
Abstract
The possible observations of Trojanlike extrasolar planets stimulate the deeper understanding of the stability behaviour of the coorbital resonant motion. By using Hill’s equations and the energyrate method an analysis of the stability map of the elliptic restricted threebody 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.
keywords:
celestial mechanics – planets and satellites: dynamical evolution and stability – methods: numerical, analytical1 Introduction
The stability of the Lagrangian triangular points and in the elliptic restricted threebody 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 threebody problem can be described by the wellknown equations of variation (Szebehely, 1967)
(1a)  
(1b) 
where
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 4th order system (1) can be separated into two independent 2nd order systems. From these we can derive an equivalent set of Hill’s equations (Meire, 1980, 1982):
(2a)  
(2b) 
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 Energyrate method
The energyrate method (ERM) (Jazar, 2004) gives the opportunity, beside the wellknown methods (Floquet theory, PoincaréLinstedt series), to investigate the stability boundary of a system
(3) 
where
(4) 
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)
(5) 
where
(6) 
corresponds to excitation in Eq. (3) as one may check in Eqs. (10).
As is wellknown if any nonconservative 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
(7) 
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
(8) 
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)
(9) 
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(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 (logscale) 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 wellknown 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, socalled ”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 subharmonic 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).
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 wellknown 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 Atype and Dtype 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 Vshape 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 Dtype 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 Vshaped 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 twodimensional 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 wellknown 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 subharmonics 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 Atype 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 subharmonics 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 superharmonics 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). Dtype 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.
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 superharmonics are depicted. An even more complete and dense structure could be obtained by overplotting the charts corresponding to different integration times.
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 energyrate 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

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 Atype resonances are situated in the intersection part of the two stable regions (long and short period librations), they are therefore the backbone of the twodimensional stable motion.

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 Btype 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.
Acknowledgments
The author wishes to thank Prof. B. Érdi for helpful discussions. This work was supported by the Hungarian OTKA Grant project number K81373.
References
 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ácsDajka E., Nagy I., Rajnai R. 2009, Celest. Mech. and Dyn. Astron., 104, 145
 Jazar (2004) Jazar G.N., 2004, International Journal of NonLinear 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 threebody problem are
(10) 
where
with
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
(11) 
where
Figure 6 demonstrates and vs. the mass parameter It can be seen that the curves exactly match the wellknown 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
Let us consider, however, a naive analytical approximation to find the resonant curves corresponding to the Btype resonances. The Btype 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
(12) 
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, Btype 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 wellseen 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).
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 Dtype resonances describe a welldefined 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 Dtype 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,Ftypes, (É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