A new perturbative solution to the motion around triangular Lagrangian points in the elliptic restricted three-body problem

The equations of motion of the planar elliptic restricted three-body problem are transformed to four decoupled Hill’s equations. By using the Floquet theorem, a perturbative solution to the oscillator equations with time-dependent periodic coefficients are presented. We clarify the transformation details that provide the applicability of the method. The form of newly derived equations inherently comprises the stability boundaries around the triangular Lagrangian points. The analytic approach is valid for system parameters 0<e≤0.05\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 < e \le 0.05$$\end{document} and 0<μ≤0.01\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 < \mu \le 0.01$$\end{document} where e denotes the eccentricity of the primaries, while μ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} is the mass parameter. Possible application to known extrasolar planetary systems is also demonstrated.


Introduction
In the era of exoplanets and specifically designed space missions, the co-orbital motion in the vicinity of the equilateral points L 4 and L 5 has become again the focus of attention (Kumar and Ishwar 2015;Elshaboury et al. 2016;Singh and Tyokyaa 2016;Barbosu et al. 2017;Zahra et al. 2017;Qian et al. 2018;Singh and Amuda 2019;Suraj et al. 2020). Since the seminal work of Szebehely (1967), the orbits near the libration points have been discussed extensively by the community. The analytic description of the Trojan-like resonant dynamics in the elliptic case of the restricted three-body problem (ERTBP) is based mainly on the averaged motion (Duffy 2012;Liang et al. 2018). Erdi (1977Erdi ( , 1978 showed the perturbation effects up to second order in Jupiter's eccentricity, perihelion and ascending node precession by using a three-parameter expansion. Morais (2001) considered an averaged disturbing potential to describe the secular variation of the Trojans' orbital elements in case of an oblate primary. Recently, Robutel et al. (2016) and Páez et al. (2016) have investigated the co-orbital resonance based on Hamiltonian formalism whereby the fast angles have been averaged out. These latter analytic studies are also capable to locate higher-order resonances as well as very slow secular frequencies.
It has been demonstrated (Tschauner 1971;Erdi 1974;Meire 1981;Matas 1982) that the coupled equations of the ERTBP can be written in the form of independent ordinary differential equations with variable coefficients. The primary goal of these studies is to explore the stability map of eccentricity-mass parameter dated back to Danby (1964). Interestingly, the analysis given by Erdi (1977) and Robutel et al. (2016) also leads to a pendulum-like equation; however, they do not attempt to solve it by classical techniques such as Floquet theorem (Lichtenberg and Lieberman 1983). Here we propose a detailed derivation of Hill's equations of the ERTBP and make a comprehensive analysis of their applicability which is still out of literature. Furthermore, analytic expressions for the solution of Hill's equations are given in the regime of moderate eccentricities and mass parameter with good agreement with numerical calculations.

Basic context
In this paper, we mainly follow the notations used in, e.g., Tschauner (1971), Meire (1980), Meire (1982). We use a dimensionless non-uniformly rotating coordinate system (x 1 , x 2 ) fixed to the primaries. The primaries having the mass m 1 and m 2 are located at (−μ, 0) and (1 − μ, 0), respectively, where μ is the mass ratio, μ = m 2 /(m 1 + m 2 ) (the origin coincides with the center of mass). The dimensionless (x 1 , x 2 ) coordinates are obtained by dividing the original dimensional coordinates of the problem by the variable distance between the primaries. The linear motion of the third body around the Lagrangian points L 4 and L 5 is determined by the coupled differential equations (Szebehely 1967) where x 1 and x 2 are the synodic Cartesian rectangular dimensionless coordinates of the third body; furthermore, Here primes denote the derivation with respect to the true anomaly v. A new form of Eqs. (1) and (2) can be introduced as By introducing the notationx = (x 1 , x 2 ) T and the vector x = (x,x ) T , the system of Eqs. (1) and (2) can be written in a more compact hypermatrix form as where 1 2 is the two-dimensional identity matrix and 0 is the zero matrix.

Hill's equation
Hill's equation has the form of where J (v) is a periodic function-in our case-with the minimal period of 2π. It is important to note that the first derivative term, ξ , is missing in Eq. (7). In order to obtain Eq. (7) from a general second-order ODE the term of y can be eliminated by the substitution where A and v 0 are arbitrary constants. We will show that Eqs. (1) and (2) can be rewritten as four coupled second-order differential equations. Let us introduce four new functions y 2 ) T , i = 1, 2 and y = (ỹ 1 ,ỹ 2 ) T with the following two properties. The first property makes possible that the new system of differential equations of y-s splits into two independent parts, namelỹ The separation is more obvious if we write Eq. (10) in the form of The second property links the new and original variables, or simply x i = y i . We will see that by this choice, the equations of the system can be rewritten into the form of Hill's equation. From the first (10) and second (12) properties, we havex =ỹ 1 +ỹ 2 = P 1 ·ỹ 1 + P 2 ·ỹ 2 = P 1 P 2 · ỹ 1 y 2 , which means that from Eq. (12) and Eq. (13) the relationship between x and y can be written as We have to calculate the elements of the matrices P i . On the one hand from Eq. (6) and Eq. (14), on the other hand from the derivatives of Eq. (11) and Eq. (14), x = 1 2 1 2 P 1 P 2 y + 1 2 1 2 P 1 P 2 y = 0 0 P 1 P 2 y + 1 2 1 2 P 1 P 2 P 1 0 0 P 2 y = P 1 P 2 P 1 + P 2 1 P 2 + P 2 2 y.
In the last two equations, the multiplication factors of y must be equal; thus, from the equality of the elements in the second rows we can write which are matrix differential equations of Riccatti type. Based on Tschauner's argument (Tschauner 1971), the following matrix elements satisfy Eq. (17) where Using these results, the inverse of the matrix T (T −1 ) can be calculated. This matrix is necessary to get the vector y from the vector x (y = T −1 · x), see Eq. (14). It is easy to verify that and Since all elements of P i contain a multiplicative factor of r , we can define the matrices Q i as Q i = P i /r with the following elements Let det Q i be the determinant of matrix Q i with the above elements. It can be shown that According to Eq. (10) If we use the relations and derived from Eq. (11), we get We have to emphasize that these two equations are not independent from each other. Equations (32) and (33) 2 are possible without any integration.
According to Eq. (9), the elimination of the first-order terms in Eqs. (34) and (35) 1 . Furthermore, as we know that 2q

, one can write
2 ). Equation (36) describes the original transformation used by Tschauner (Tschauner 1971). The derivatives of y By using the above equations, and similar derivatives for y (i) 2 , finally, Eqs. (34) and (35) can be written in the form of Hill's equation as where the coefficients J (i) j (i, j = 1, 2) are periodic functions with period of 2π. At this point, it is necessary to consider the numerical solvability of these equations. In this study, we focus on the parameter range only where the value of c(μ, e) in Eq. (22) is a strictly positive real number (c > 0) and μ < 1/3. This parameter range is the shaded region in Fig. 1.
The numerical solvability of Eqs. (38) and (39)    In contrast to the analysis of q (i) 21 , the shaded region can be divided into three subregions, which are denoted by I, II and III in Fig. 1. For brevity, Table 1 shows the signs of the minimum and maximum values in the different subregions. We can see that q (1) 21 changes its sign (and therefore goes through zero) in the subregions II and III during the integration, which means that it is not recommended to solve the equation for ξ (1) 2 in these subregions. It is interesting to note that in the subregion I the sign is always negative, which means that q (1) 21 is purely imaginary. Because of the former equality y (1) 2 must also be purely imaginary. Fortunately, this does not mean that the equation for ξ (1) 2 is unsolvable but we have to consider that ξ (1) 2 remains always purely imaginary; thus, multiplying the equation by the imaginary unit solves this problem. The situation is similar in the case of q (2) 21 , but the sign change happens only in the subregion III (see Table 1). Following the above analysis, the numerical integration has been performed in order to validate the results. Figure 2 depicts the trajectory for e = 0.048, μ = 0.000954 (the case of Jupiter). The solution of Eqs. (1) and (2) and Eqs. (38) originating from the appropriate initial conditions perfectly overlap. This means that Hill's equations can be applied to solve the equations of motion around the L 4 and L 5 points.
ficients, can be solved by Floquet theorem (Hagel 1992). We seek the solution in the form of where w(v) is the so-called Floquet function, which has the same period as J (v). Constants a and b are determined by the initial conditions. Since the derivation for both ξ (1) 1 and ξ (2) 1 are the same, we omit the indices in the rest part of the paper. Let us rewrite Eqs. (38) for w(v) and ψ(v) The above differential equations split into two parts with the coefficients of sin and cos From Eq. (44), we obtain where C is a constant of integration. At the end of the calculation, we will see that a 2 = e C that is With this equation, Eq. (43) becomes Now we are looking for the solution of w(v) in a third-order Taylor series in the eccentricity e where w (i) are 2π periodic functions of v.

Taylor series of J (1) 1 (v) and J (2) 1 (v)
The periodic coefficients to be solved have complicated forms; therefore, the solution can be obtained by a third-order Taylor expansion in the eccentricity. Let us utilize J by using the earlier introduced notations. Useful expressions will be B i ≡(2c 2 + 1 + (−1) i λ) −1 and λ ≡ √ 1 − 9g. The third-order Taylor expansion of J (i) 1 is then where Let us write back the results of the Taylor expansions into Eq. (46), and use the fact that Then we can collect the terms for e 0 , e 1 , e 2 and e 3 ; thus, 4 new differential equations can be obtained (also for i = 1, 2, so for simplicity we omit index i) for the terms of w(v): It can be seen that a particular solution corresponding to 2π periodicity for Eq. (51) is The differential equations (52) ih (v). The homogeneous part of Eq. (52) is which is a harmonic oscillator with frequency (4α) 1/2 ; thus, the solution of the equation is Generally, in order to fulfill the 2π periodicity of w(v), the constants must be chosen as follows K 1 = K 2 ≡ 0. For the inhomogeneous solution, we use the following trial function where w 1,1 are w 1,0 constants. By calculating the derivatives from the coefficients, we can simply obtain the values of w 1,1 and w 1,0 , namely We use the same steps for the solution of Eq. (53). By using trigonometric identities, it can be seen that the differential equation has the following form Like in the previous case the solution of the homogeneous part is where again K 1 and K 2 must disappear for the 2π periodicity, K 1 = K 2 ≡ 0. The trial function of the particular solution of the inhomogeneous equation is: Again by calculating the appropriate derivatives, the equality of the coefficients implies: Only the solution of Eq. (54) is left The homogeneous solution reads where again the constants are K 1 = K 2 ≡ 0 due to the periodicity of w(v). The trial function for the particular solution of the inhomogeneous equation is where the forms for w 3,1 and w 3,3 coefficients are w 3,1 = − w 0,0 ε + w 1,1 γ + 1 2 w 1,1 δ + w 2,0 β + 1 2 w 2,2 β − 12αw 1,1 w 2,0 w 0,0 − 6αw 1,1 w 2,2 w 0,0 + 15αw 3 1,1 2w 2 Then by using the fact that ψ (v) = a 2 w −2 (v) (see Eq. (46)), ψ(v) can be calculated if we again expand ψ (v) into Taylor series in e up to third order Now we have the expressions for w(v) and ψ(v); thus, ξ(v) = aw(v) cos ψ(v) + b can be calculated. It is left to determine the constants a and b, which are controlled by the initial conditions ξ(0) ≡ ξ 0 and ξ (0) ≡ ξ 0 . As the differential equations are second-order linear differential equations with periodic coefficients, the initial conditions can be arbitrary; therefore, we use the simple conditions of x 0 = 1, y 0 = 1, v x0 = 0, v y0 = 0, from which ξ 1,0 , ξ 1,0 , ξ 2,0 and ξ 2,0 can be easily achieved. By using the values ξ 0 and ξ 0 , At the end, the only task is to use the transformations detailed in Eqs. (37), calculate y

Illustrations and discussion
The prominent example of co-orbital dynamics is the Sun-Jupiter-Trojan configuration in our own Solar system. We apply the perturbative solution described in Sect. 3.1 to this structure first. Figure 3 depicts the trajectory around the Sun-Jupiter triangular Lagrangian point. The integration time is 20 periods of Jupiter (ca. 240 years). The analytic and numerical solutions match perfectly, although after some time (∼ 38 − 40 periods) they start to deviate.
Recently, Lillo-Box et al. (2018) have studied the physical parameters and dynamical properties of possible exo-Trojans in systems with close-in (orbital period < 5 days) planets. We selected two of them, HAT-P-20b (e = 0.015, μ = 0.0091) and WASP-36b (e = 0.0, μ = 0.0021), to provide the analytic solution in these regimes 1 . The orbits are plotted in Fig. 4a and b, respectively. The panels show the paths for T =20 periods again. Due to the zero eccentricity of the planet, the analytic solution for WASP-36b remains very close to the numerical outcome for much longer times. Considering the Earth-Moon system with e = 0.054 and μ = 0.012, it falls close to the limit of third-order solution. The analytic solution diverges after 5-6 revolutions (∼ 130−150 days) of the Moon. We have seen that for the Sun-Jupiter system the analytic curve traces the numerical method reasonably well, while the eccentricity falls into the same range. In addition, we have found that the rather large mass parameter-compared to planetary systems-does not affect the precision of the analytic solution provided the eccentricity is small enough, practically zero. This is, however, not the case for the Moon. Consequently, systems with moderate nonzero eccentricity and mass parameter of the same size require an improved analytical solution, e.g., higher-order expansion in mass.
In this work, we described the motion around the triangular Lagrangian points using Hill's equations. As a perturbative solution, a third-order expansion of Floquet function w(v) in eccentricity was presented. This method is capable to follow analytically the orbit of a massless particle around the equilibrium points L 4 and L 5 in the ERTBP. A precise trajectory forecast for moderate eccentricity (e ≤ 0.05) and mass parameter (μ ≤ 0.005) is achievable for tens of secondary's orbital periods.