Position and velocity sensitivities at the triangular libration points in the restricted problem of three bodies when the bigger primary is an oblate body

In this paper we have examined the stability of triangular libration points in the restricted problem of three bodies when the bigger primary is an oblate spheroid. Here we followed the time limit and computational process of Tuckness (Celest. Mech. Dyn. Mech. 61, 1–19, 1995) on the stability criteria given by McKenzie and Szebehely (Celest. Mech. 23, 223–229, 1981). In this study it was found that in comparison to other studies the value of the critical mass μc has been reduced due to oblateness of the bigger primary, i.e. the range of stability of the equilateral triangular libration points reduced with the increase of the oblateness parameter I and hence the order of commensurability was increased.

0.038521 . . . , the critical mass. Leontovich (1962) established the non-linear stability of the triangular libration points and proved that the triangular libration points L 4 and L 5 are stable for all values of μ < μ c except for a set of measure zero, where μ c is the Rouths critical mass of the Restricted three-body problem. Breakwell and Pringle (1965) discussed the motion of the third body around L 4 in the Earth-Moon system using Von Zeipel's method with an additional effect of the Sun's gravitation. Also he did a detailed study of the stability of the triangular libration points. Deprit and Palmore (1966), Deprit and Deprit-Bartholome (1967), , Markeev (1969) all discussed stability of triangular libration points with different methods and also established the order of commensurability of the long periodic and short periodic orbits. Alfriend (1970) and Nayfeh (1971) justified the work of Markeev (1969) that L 4 is an unstable equilibrium point where the long periods and short periods of motion around L 4 have commensurability 2 to 1. Henrard (1970) established that the family of long periodic orbits at L 4 does not evolve in a continuous manner with the mass ratio. Discontinuity appears not only at the mass ratios for which the long period of the equilibrium is a multiple of the short periods, but also at mass ratios when the global analysis of the family detects singular bifurcation orbits. Markeev (1973) and Sokol'skii (1975) studied the stability of the Lagrange solutions for the critical mass μ c .
McKenzie and Szebehely (1981) introduced a new idea to test the stability of the equilateral triangular libration point L 4 . They defined the maximum velocity and maximum displacement envelopes within which the third body remains for a long time starting from the suitable initial conditions so that the third body may not cross the x-axis. The conditions of not crossing the x-axis were introduced as the stabil-ity criteria for the third body. By using all the stability criteria of McKenzie and Szebehely (1981), Tuckness (1995) excellently, investigated numerically the sensitivities of the third body around L 4 when it is given positional and velocity deviations away from L 4 with a suitable initial conditions. He used Poincare's surface of sections to compare the periodic, quasi-periodic and stochastic regions to the trajectories with the definitions of stability given by McKenzie and Szebehely (1981). Moreover he also investigated the value of μ (the mass ratio) ranging from zero to linear critical mass 0.038521 . . . . Using the stability criteria he determined some values of μ which are more stable than the others.
Here we extend the study of Tuckness (1995) for analyzing the effect of oblateness of the bigger primary on the stability of L 4 . All the other conditions and criteria will remain the same as those chosen by McKenzie and Szebehely (1981) and Tuckness (1995).

The equations of motion of the third body
The equations of motion of the third body moving in the gravitational field of the two primaries, arë (1) where Ω = n 2 2 (1 − μ)r 2 1 + μr 2 2 + F ; with n is the mean angular velocity of the primaries about their centre of mass; I = A 1 − A 3 = oblateness parameter (0 < I 1); A 1 = a 2 5R 2 = the moment of inertia of the oblate body about the equatorial radius; A 3 = c 2 5R 2 = the moment of inertia of the oblate body about the polar radius; and

Location of libration points
For location of libration points the Jacobi's integral (or Jacobi's manifold) iṡ where C is called Jacobi's constant. Libration points are the singularities of the manifold, Following Szebehely (1967), Bhatnagar and Chawla (1977) the equilateral triangular libration points are given as

Mathematical formulation for critical mass
Let (a, b) be the coordinates of the triangular libration point L 4 . Let ξ, η be the small variations in a & b respectively due to oblateness of the bigger primary, then (a + ξ, b + η) is a point in the vicinity of (a, b). Thus the variational equations of motion of the third body are given bÿ where Let ξ = P e λt &η = Qe λt be two particular solutions of the Eqs. (7), then the characteristic equation for the triangular libration points is where λ 2 = ∧, Let ∧ 1 and ∧ 2 be the roots of the Eq. (8), then i.e.
Equation (13) is a quadratic equation in μ so the roots of this equation will give the critical value of the mass parameter μ.
If μ c is the critical value of μ then If

Numerical integration
The numerical integration in this work were carried out using the classical fourth order Runge-Kutta method with automatic step-size control. The specified accuracy was at least eight or more significant decimal digits and double precision arithmetic was used in all computations. Tests were performed with different levels of accuracy and the specified accuracy was found to be sufficient for this study.
The Jacobian constant was calculated in order to observe the validity of the trajectory generations throughout the numerical integration. This constant of integration was checked to make sure that it remained constant to at least eight decimal places.
For numerical integration, we introduce x = x 1 , y = x 2 ,ẋ = x 3 ,ẏ = x 4 and reduce the two second order differential equations of motion given in (1) to the four first order differential equations as:

Time limit for numerical integration
The study of the time limit was first undertaken by McKenzie and Szebehely (1981). The value of μ used in their study was μ = 0.01214 (Earth-Moon system). The selection of the time interval t f for the numerical integration was determined by examining the effect of t f on the maximum distance that the third body was started from L 4 along two given lines before its motion crossed the x-axis. As seen in Fig. 1, the first line points 120 • away from L 4 and the origin of the coordinate system. The second line is inclined at 300 • with respect to the x-axis, pointing toward the origin of the coordinate system. The dimensionless values of t f they chose to study were: t f = 120, 240, 360, 480, 600, 720, 840, 960, corresponding to 564, 1128, 1692, 2256, 2820, 3384, 3948, and 4512 days for the Earth-Moon system. From Fig. 2, it can be seen that, as the time of integration increases, the maximum distance the third body may be displaced initially from L 4 for librational motion to take place decreases. Also in Fig. 2 the maximum values of the initial velocity as a function of integration time is shown. After the first large change in the distance between 120 and 240 time units, the effect of the final time of integration appears to change the initial velocity by only small amounts. Since accuracy of the numerical integration becomes more questionable with an increase in the time of integration hence a reasonable value of t f = 600 is taken. By using Runge-Kutta method, further investigations were carried out by Tuckness (1995) with two other values of μ = 0.0027 and μ = 0.025, for θ = 10 • , 20 • , 30 • , 40 • , . . . , 360 • . The values of μ = 0.0027 and μ = 0.025 were chosen because of the nature of maximum displacement and maximum velocity envelopes take place, whereas in our case, for the same values of θ and the same time limit (as in Tuckness), the two values of μ were reduced to μ = 0.0026 and μ = 0.024 due to the oblateness parameters I = 0.01. For I = 0.1 the corresponding value of μ is given by μ = 0.024 for maximum velocity and displacement and two values of μ = 0.012 and μ = 0.021 were investigated for nearly non-librational motion around L 4 . In short, the area of the maximum velocity and maximum displacement envelope for μ = 0.0026 is a maximum and a minimum for μ = 0.024. As McKenzie and Szebehely (1981) and Tuck- The value of t f = 100 did not allow enough time for the third body to librate around to the x-axis. This allowed much larger velocity and positional deviations. Therefore, from the studies undertaken by McKenzie and Szebehely (1981), Tuckness (1995) and from our investigations t f = 1000, dimensionless time units would be an acceptable integration time duration for μ = 0.0026 and μ = 0.024.

Comparison of maximum velocity and maximum displacement envelopes using Poincare's surface of sections
In this section the maximum velocity and maximum displacement envelopes boundary values were investigated using time limits for numerical integration given by Tuckness (1995) and Poincare's surface of sections. For this we need Hamilton's equations of motion to be defined as follows: Let at any time t, (x 1 , x 2 ) be the position of the third body with momenta p 1 and p 2 , then Now the time rate of change of phase variables are: Thuṡ The Hamiltonian is given by where F is given in (5).
A solution of the Hamiltonian equations of motion (18) can be represented as a trajectory in a four dimensional phase space. Because of the existence of the integral of motion Eq. (19), the trajectory lies in a three-dimensional subspace (H ) which is equal to a constant of the phase space. The successive intersection of this three-dimensional trajectory with a two dimensional surface is called Poincare's surface of section. The successive intersections of the 3-D trajectory with the surface y = √ 3 2 [1 + 1+2μ 3(1−μ) I ] was taken for Poincare's surface of section in (x,ẋ) or (x 1 , p 1 ) plane, the results can be represented in a much more compact manner. There is no serious loss of information, because the most interesting trajectories of the projection are reflected in the corresponding properties of the set of points.
In case of Tuckness (1995) the phase space can be divided into three regions. One region contains primarily stochastic trajectories. It is connected to a second region containing large embedded islands. A third region consists of primarily regular trajectories and is isolated from the first two regions. From the Fig. 3 the phase space can be divided into three regions. One region contains stochastic trajectories not connected with the second region because of Fig. 3 Surface of Poincare section for μ = 0.001, θ = 108 • rare embedded islands due to oblateness parameter I . The third region contains primarily regular trajectories which are isolated from the first two regions. The maximum velocity and maximum displacement were investigated using the Poincare's surface of sections method for various values of μ. Initial velocities and positions from the maximum velocity and maximum displacement envelopes were numerically integrated using Hamilton's equations of motion from (16) to (19) and their intersection with the surface y = √ 3 2 [1 + 1+2μ 3(1−μ) I ] were plotted. The equations of motion were integrated for 500 to 2000 orbits (or intersection) using the same Runge-Kutta integrator explained above and only the orbits for y > 0 were plotted.
The velocity VL 4 = 0.444 is the value of maximum velocity envelope in the 108 • direction which was found using the stability criteria under consideration. Figure 3 is the Poincare's surface of section investigated for μ = 0.001, θ = 108 • . Velocities greater than VL 4 = 0.444 result in unstable motion according to the boundary envelopes established in the study. Velocities less than 0.444 are contained within the surfaces of boundary envelops depicted in Fig. 3. Also areas of apparent quasi-periodic solutions can be investigated from the Fig. 3 but it is too difficult to show. Figure 4 depicts the Poincare surfaces of section investigated for μ = 0.001, θ = 0 • and maximum initial displacement allowed. In the enlarged view the chaotic regions had been reduced in comparison of Tuckness due to oblatness I .

Determination of stability of the third body using the maximum velocity and displacement envelopes
As seen earlier in this study, the sizes of the maximum velocity and maximum displacement envelopes vary according to the value of μ, I , the final time limit t f and the As seen from the previous section, when with the increase of I (oblateness) μ decreases, the maximum velocity and maximum displacement envelopes become more narrow. The maximum velocity and displacement envelopes consisting of initial velocities and displacements respectively away from L 4, where according to previous definition, stable librational motions result. One could also say that arriving at L 4 with the velocity and direction of the maximum velocity allowed would give the same results. Therefore, the maximum velocity envelope above would be thought of as maximum velocity arrival error envelopes, which could be thought of as the maximum displacement arrival error envelope. One could also think of the value of μ, which allowed the most arrival velocity error at L 4 , as being the most stable value of μ for arriving at velocity errors. Again, the same is true for the maximum displacement envelope.
The value of μ and I which allows the total amount of displacement error (all the displacement vectors of the maximum displacement envelope summed up) will have a maximum displacement envelope with the greatest enclosed area. Tuckness (1995) calculated the area of the displacement envelopes by considering the formula πab, where a and b are semi-major axis and semi-minor axis respectively. Also he considered the maximum displacement envelopes as the symmetric and uniform ellipse but the displacement envelopes are neither symmetric nor uniform. So we calculated the displacement envelope and velocity envelope by considering the formula 1 2 r 2 dθ with the limits of θ , the initial point of envelope to the final point of envelope.
Because the estimated area of the envelopes is a measure of stability, a comparison can be made on how sta- and μ = 0.024 gives the maximum velocity and displacement envelopes for I = 10 −1 . It was also found in Figs. 7 and 8 that μ = 0.024 approximately has nearly zero area in both the maximum velocity and displacement envelopes for I = 10 −2 while μ = 0.012 and μ = 0.023 both represent nearly zero displacement and velocity for I = 10 −1 . This helps to reinforce the ideas in the previous section that mo- approximately are non-librational or unstable for I = 10 −2 , I = 10 −1 respectively. In other words, for μ = 0.024 and 0.023 approximately, very little velocity or displacement away from L 4 will allow the third body to stay around L 4 instead, it will cross the x-axis.
These appear to be distinct minimums which depict that librational motion is difficult to achieve. The minimums marked on Figs. 7 and 8 show a distinct relation to the values of μ and I which are related to the critical values; namely for three to one commensurability μ c = 0.0135160160, 0.0135156540, 0.0135065589, 0.0134260594 and for two to one commensurability the critical masses are μ c = 0.02429389714, 0.0242921132,  . 9 Graph of critical mass versus I 0.0242761497, 0.024125161 corresponding to I = 0, 10 −4 , 10 −3 , 10 −2 respectively. Following Pedersen (1933), Deprit and Deprit-Bartholome (1967) and  the ratio of the long period and short period orbits, s 1 = 1 1 + k 2 &s 2 = k 2 s 1 .
The expression for k 2 as in (11) and value of μ c as in (14) show that the frequency of occurrence of the critical values of the mass parameter increases as μ approaches zero.  Table 1).
The k values of one to seven are depicted in Figs. 7 and 8 which indicate an exponentially decreasing curve as a function of μ and I . From this observation it may be possible to arrive at an analytical expression that is the function of μ and I which governs the stability around the triangular libration point for the specific values of μ and I .
The curves that depict the areas of maximum velocity and displacement envelopes in Figs. 7 and 8 appear to have the same mathematical characteristic differing only in amplitude.

Conclusion
From Figs. 5 & 6 and graphs 7, 8, 9 & 10 we conclude that the range of stability reduces with the increase of oblateness parameter I . In other words we can say that with increase of I , the order of commensurability k increased. Ultimately the long period and short period will be changed due to I . The discussion can also be corroborated by the Table 1.