This paper investigates the secular motion of a massless asteroid within the framework of the doubleaveraged elliptic restricted threebody problem. By employing Poincar variables we analyze the stability properties of asteroid orbits in the presence of planetary perturbations. Our study reveals that periodic orbits identified in the planar configuration maintain stability in the spatial perturbed problem across a wide range of parameter values. These findings supported by numerical simulations contribute to a deeper understanding of asteroid dynamics and have implications for studying exoplanetary systems with highly eccentric host stars.
['Astrophysics - Earth and Planetary Astrophysics', 'Mathematics - Dynamical Systems']
Stability analysis of spatial perturbed elliptic restricted 3body problem with doubleaveraging
{'Stability analysis of spatial perturbed elliptic restricted 3-body problem with double-averaging': 'Yan Luo 2 , Kaicheng Sheng 1 ∗  \n1  \nSchool of Mathematics, Shandong University, Jinan, China 2 Research Centre for Mathematics and Interdisciplinary Sciences, Shandong University, Qingdao, China', 'Abstract': "This paper investigates the secular motion of a massless asteroid within the framework of the double-averaged elliptic restricted three-body problem. By employing Poincar'e variables, the stability properties of asteroid orbits in the presence of spatial perturbations were analysed. The study reveals that periodic orbits identified in the planar configuration maintain stability in the spatial perturbed problem across a wide range of parameter values. These findings, supported by numerical simulations, contribute to a deeper understanding of asteroid dynamics and have implications for studying exoplanetary systems with highly eccentric host stars.  \nKeywords : 3-body problem, Double-averaging, Periodic orbit, Linear stability.", '1 Introduction': "The study of the stability of celestial bodies within gravitational systems is a foundational aspect of celestial mechanics. The study of the 3-body problem formed the foundation for understanding the complex dynamic behaviours observed in natural and artificial systems. The restricted 3-body problem plays an important model role in the 3-body problem model. This paper focuses on the double-averaged spatial perturbed elliptic restricted 3body problem (ER3BP) involving a star, a planet, and an asteroid. The study is under the assumption that the mass of the planet is much smaller than the mass of the star. The averaging over fast phases corresponds to the motions of the star-planet system and the asteroid. The classical restricted 3-body problem has been instrumental in developing our understanding of orbital mechanics, particularly in cases where the orbits of the primary bodies are assumed to be circular. However, in reality, most celestial bodies, including stars, planets and asteroids, follow elliptical orbits, which introduces significant complexities to the problem. These complexities necessitate more advanced mathematical techniques and have led to the development of the ER3BP, where the primary bodies follow elliptical orbits.  \nDouble-averaging method simplifies the complex gravitational interactions, resulting in a simplified system that retains the essential dynamical features of the original problem while being more tractable for analysis over long timescales. The seminal work of Aksenov [1979] on the double-averaged elliptical restricted 3-body problem laid the foundation for understanding the long-term evolution of orbits in such systems. His work demonstrated how the double-averaging method could be used to reduce the complexity of the problem, making it possible to derive analytical expressions for the evolution of orbital elements over time.  \nThe significance of stability in celestial mechanics, particularly in the context of Hamiltonian systems, has been explored by Arnold [1961], who established fundamental results on the stability of equilibrium positions. This work is further elaborated in the comprehensive treatment of classical mechanics by Arnold et al. [2006]. Moreover, the secular evolution of orbits, especially in hierarchical systems, has been investigated by Lidov [1962] and Kozai [1962], whose work on the Lidov-Kozai mechanism has profoundly influenced our understanding of orbital dynamics in N-body systems. The relevance of the Lidov-Kozai mechanism has been further underscored by recent studies, such as those by Katz et al. [2011] and Lithwick and Naoz [2011], who explored its effects in systems with eccentric perturbers.  \nThe stability of planar orbits in double-averaged circular problems was established by Neishtadt [1975]. His study indicates that planar orbits also remain stable in the linear approximation of the double-averaged elliptic problem with a sufficiently small eccentricity of the perturber's orbit. The stability of periodic solutions in the restricted 3-body problem, particularly in the presence of elliptical orbits, has also been explored by Leontovich [1962] and Moser [1968], whose work on Hamiltonian systems with two degrees of freedom provided crucial insights into the resonance phenomena that often govern the stability of such systems. Harrington [1968] and Ziglin [1975] analyzed the dynamical evolution of the 3-body problem, while more recent studies, such as those by Naoz et al. [2013] and Sidorenko [2018], have extended these ideas to modern exoplanetary systems. The double-averaged ER3BP specifically benefits from these approaches, offering insights into the stability conditions that govern the long-term behaviour of small bodies in elliptical orbits.  \nIn recent years, there has been renewed interest in the study of the double-averaged ER3BP, particularly in the context of exoplanetary systems. The discovery of numerous exoplanets with highly eccentric orbits has prompted researchers to revisit the problem of stability in such systems. Studies by Huang and Lei [2022] and Lei [2022] have focused on the linear stability of the inner case of the double-averaged spatial elliptic restricted 3-body problem, providing new insights into the conditions under which stability can be achieved in these systems. A completed numerical analysis, including all possible situations for the evolution of planar orbits in the double-averaged ER3BP, can be found in the work by Vashkovyak [1982].  \nThe investigation of equilibria within the double-averaged ER3BP, as studied by Palaci'an et al. [2006], and the analysis of apsidal alignment by Neishtadt et al. [2021] analyzed both the linear and nonlinear stability of apsidal alignment in the spatial double-averaged ER3BP, further illustrated the rich dynamical behaviour that emerges in such systems. Linear stability of the equilibria and periodic orbits of the asteroid in spatial double-averaged ER3BP in the inner case with arbitrary inclination is studied in Huang et al. [2024].  \nThis paper aims to provide a comprehensive analysis of the stability of the spatial  \nperturbed double-averaged ER3BP, building on the foundational work of Aksenov [1979], Vashkovyak [1981, 1982], Neishtadt et al. [2021] and others. By combining analytical methods with numerical simulations, the long-term stability of the asteroid's orbits in doubleaveraged ER3BP is investigated in linear approximation. The results have important implications for both theoretical and practical aspects of celestial mechanics, providing new insights into the stability of orbits in a wide range of astrophysical systems. Through this work, we hope to contribute to the ongoing exploration of the dynamic and ever-evolving nature of the cosmos.  \nA notable characteristic of many exoplanets is their orbits, which often exhibit large eccentricities and inclinations, deviating significantly from the near-circular and coplanar trajectories observed in the solar system [Stephen et al., 2012, Team, 2017]. Our research holds the potential to make substantial contributions to the domain of exoplanetary dynamics [Shevchenko, 2016]. Furthermore, the secular evolution of motions within the non-restricted three-body problem has been extensively studied through the application of averaging techniques [Harrington, 1968, Ziglin, 1975, Michtchenko and Malhotra, 2004], one can extend our work into the non-restricted 3-body problem.", '2 Statement of the problem': "Figure 1: Coordinate frames.  \n<!-- image -->  \nThe considered spatial perturbed elliptic restricted three-body problem involves a star S , a planet J , and an asteroid A [Brouwer and Clemence, 1961]. Adopting a coordinate system analogous to that in Neishtadt et al. [2021], the star is positioned at the origin O of a Cartesian frame Oxyz . The Oxy plane of this system is defined by the orbital plane of the star-planet system. Consequently, the coordinates of the planet and the asteroid are ( x J , y J , 0) and ( x, y, z ), respectively. A rotating Cartesian frame Ox ' y ' z ' is introduced, where the Ox ' y ' plane aligns with the asteroid's osculating orbital plane. In this frame, the asteroid's coordinates are ( x ' , y ' , 0). The asteroid's orbital motion is characterized by standard osculating elements: semi-major axis ( a ), mean anomaly ( l ), eccentricity ( e ), argument of periapsis ( ω ), inclination ( i ), and longitude of the ascending node (Ω). Refer to Fig. 1 for a visual representation of the system geometry. According to the rotation of the Cartesian  \nframe, the transformation between frames Oxyz and Ox ' y ' z ' is  \n[ x y z ] T = T Ω ×T i ×T ω × [ x ' y ' 0 ] T , (2.1)  \nwhere  \nT Ω =   cos (Ω) -sin (Ω) 0 sin (Ω) cos (Ω) 0 0 0 1   , T i =   1 0 0 0 cos( i ) -sin( i ) 0 sin( i ) cos( i )   , T ω =   cos ( ω ) -sin ( ω ) 0 sin ( ω ) cos ( ω ) 0 0 0 1   . (2.2)  \nConsequently, the expressions of x , y and z are  \nx = (cos Ω cos ω -cos i sin Ω sin ω ) x ' +( -cos Ω sin ω -cos i sin Ω cos ω ) y ' , y = (sin Ω cos ω +cos i cos Ω sin ω ) x ' +( -sin Ω sin ω +cos i cos Ω cos ω ) y ' , z = (sin i sin ω ) x ' +(sin i cos ω ) y ' . (2.3)  \nIt follows from Shevchenko [2016], the planet moves in a prescribed elliptic orbit:  \nx J = a J (cos E J -e J ) , y J = a J √ 1 -e J 2 sin E J , l J = E J -e J sin E J . (2.4)  \nHere a J , e J , E J , l J are the semi-major axis, the eccentricity, the eccentric anomaly, and the mean anomaly of the planet's orbit. We put a J = 1 for convenience in the following.", '3 Hamiltonian of the system': "We introduce the canonical Delaunay elements ( l, g, h, L, G, H ), where g ( ≡ ω ) and h ( ≡ Ω) are the argument of pericenter and the ascending node of the asteroid. The elements L = √ (1 -µ ) a , G = L √ 1 -e 2 , and H = G cos i correspond to the Keplerian energy, total angular momentum, and z -component of angular momentum, respectively [Brouwer and Clemence, 1961]. To facilitate the dynamical analysis of the asteroid, canonical Poincar'e variables are denoted by  \np 1 = L, q 1 = l + g + h, p 2 = √ 2 ( L -G ) cos ( g + h ) , q 2 = -√ 2 ( L -G ) sin ( g + h ) , p 3 = √ 2 ( G -H ) cos h, q 3 = -√ 2 ( G -H ) sin h. (3.1)  \nAssuming the planet and star have masses µ and 1 -µ , respectively, such that the total system mass equals unity. µ ≪ 1. The Hamiltonian of the asteroid is given by:  \nF = -(1 -µ ) 2 2 L 2 + µU -µ ( x x J + y y J ) , (3.2)  \nwhere  \nU = -V = -1 √ ( x -x J ) 2 +( y -y J ) 2 + z 2 (3.3)  \nis the perturbing gravitational potential which is expressed in terms of the asteroid's coordinates ( x, y, z ). These coordinates can be transformed into Poincar'e variables with equations (2.3) and (3.1). The asteroid's motion in its elliptic orbit is governed by the equations:  \nx ' = a (cos E -e ) , y ' = a √ 1 -e 2 sin E, l = E -e sin E, (3.4)  \nwhere E denotes the eccentric anomaly [Shevchenko, 2016]. The planet's coordinates, x J and y J , are considered prescribed functions of time.  \nThe double-averaged Hamiltonian ¯ F of the asteroid is defined as:  \n¯ F = 1 (2 π ) 2 ∫∫ [0 , 2 π ] 2 Fdldl J = -(1 -µ ) 2 2 L 2 + µ ¯ U = -(1 -µ ) 2 2 L 2 -µ ¯ V , (3.5)  \nwhere  \n¯ U = 1 (2 π ) 2 ∫∫ [0 , 2 π ] 2 Udldl J , (3.6)  \nand µ ¯ V = -µ ¯ U is the double-averaged force function of gravity of the planet.  \nPerforming the double-averaging over the mean anomalies of the planet and the asteroid, the Hamiltonian becomes independent of q 1 . Consequently, the conjugate momentum p 1 (equivalent to L ) constitutes a first integral of the system. The first term in the doubleaveraged Hamiltonian, ¯ F , thus remains constant. The dynamics of the remaining variables, ( p 2 , q 2 ), ( p 3 , q 3 ), are therefore governed by a two-degree-of-freedom (2-DOF) Hamiltonian system with Hamiltonian µ ¯ U .  \nBy introducing the slow time variable τ = µt , the Hamiltonian transforms to ¯ U . Notably, ¯ U depends on the planet's eccentricity, e J , as parameters, and the ratio between the semimajor axis a of the asteroid and a J = 1 of the planet.  \nTaking i = 0, the spatial perturbed ER3BP reduces to a planar problem on the invariant plane p 3 = 0, q 3 = 0. This planar system is described by a one-degree-of-freedom (1-DOF) Hamiltonian, ¯ R Θ ( p 2 , q 2 ), where ¯ R Θ ( p 2 , q 2 ) = ¯ U ( p 2 , q 2 , p 3 = 0 , q 3 = 0) is independent of p 3 and q 3 . Following Aksenov [1979] and Neishtadt et al. [2021], the Hamiltonian of the double-averaged planar ER3BP is given by:  \n¯ R Θ = 1 (2 π ) 2 ∫∫ [0 , 2 π ] 2 Rdldl J , (3.7)  \nwhere  \nR Θ = -1 √ ( x -x J ) 2 +( y -y J ) 2 . (3.8)  \nNumerical investigations by Vashkovyak [1982] have demonstrated the existence of stationary solutions (equilibria)  \nq 2 = 0 , p 2 = p 2 ∗ , (3.9)  \nfor some domains in the plane of parameters a , e J in the double-averaged planar ER3BP. Periodic orbits, separatrix and orbits of circulating are included as well. Fixing the parameter a and e J in ¯ R Θ , one can obtain the figure of Θ and e . See Fig. 2.  \nFigure 2: Trajectories of Θ, e in double-averaged planar ER3BP.  \n<!-- image -->  \np 2 ∗ is the root of the equation  \n∂ ¯ R ∂p 2 ∣ ∣ ∣ ∣ q 2 =0 = 0 . (3.10)  \nThe value of the eccentricity at the equilibrium (3.9) is provided by the equation  \ne ∗ = √ 1 -( 1 -p 2 2 ∗ 2 √ a ) 2 . (3.11)  \nAs detailed in Vashkovyak [1982], the equilibria defined by equation (3.9) correspond to an apsidal alignment scenario where ω +Ω = 0. (Fig. 3a) The study also identifies periodic orbits in terms of Θ = ω +Ω and the asteroid's eccentricity ( e ) surrounding these equilibria. (Fig. 3b)  \n̸  \nTo comprehensively assess the linear stability of asteroid orbits in the spatial perturbed ER3BP, it is essential to consider non-equilibrium cases (i.e., Θ = 0), especially to consider the cases of periodic orbits. The small perturbation case is considered in this paper, which means the inclination i is regarded as a small perturbation. In this case, the formula of x , y and z can be obtained with i = 0 and ω +Ω=Θ:  \nx = cos(Θ) x ' -sin(Θ) y ' , y = sin(Θ) x ' +cos(Θ) y ' , z = 0 . (3.12)  \nFor small i , proper approximations of formulas (3.1) are obtained:  \ncos (Ω) ≈ p 3 √ G · i 2 / 2 , sin (Ω) ≈ -q 3 √ G · i 2 / 2 , cos ( ω ) ≈ p 2 p 3 -q 2 q 3 √ G · i 2 √ p 2 2 + q 2 2 , sin ( ω ) ≈ -p 2 q 3 + p 3 q 2 √ G · i 2 √ p 2 2 + q 2 2 . (3.13)  \nFigure 3: Case of equilibria (apsidal alignment) and periodic orbits.  \n<!-- image -->  \nSubstituting (3.13) into (2.3) and then taking in (3.3), the function U = U Θ with increasing order of p 3 and q 3 can be written as  \nU Θ = R Θ + W Θ + O ( p 4 3 + q 4 3 ) , W Θ = 1 2 ( A Θ p 2 3 +2 B Θ , p 3 q 3 + C Θ q 2 3 ) (3.14)  \nwhere  \nA Θ = yy J G [ ( x -x J ) 2 +( y -y J ) 2 ] 3 / 2 , B Θ = ( xy J + yx J ) 2 G [ ( x -x J ) 2 +( y -y J ) 2 ] 3 / 2 , C Θ = xx J G [ ( x -x J ) 2 +( y -y J ) 2 ] 3 / 2 . (3.15)  \nHere all quantities are calculated at i = 0 and Θ = ω +Ω. Calculations are similar to the procedure in Neishtadt et al. [2021].  \nAveraging W Θ over the mean anomaly of the asteroid l and the mean anomaly of the planet l J , the averaged value of W Θ is  \n¯ W Θ = 1 2 ( ¯ A Θ p 3 2 +2 ¯ B Θ p 3 q 3 + ¯ C Θ q 3 2 ) , (3.16)  \nwhere  \n¯ A Θ = 1 4 π 2 ∫∫ [0 , 2 π ] 2 ( A Θ ) d l d l J , ¯ B Θ = 1 4 π 2 ∫∫ [0 , 2 π ] 2 ( B Θ ) d l d l J , ¯ C Θ = 1 4 π 2 ∫∫ [0 , 2 π ] 2 ( C Θ ) d l d l J (3.17)  \nare the averaged values of A Θ , B Θ and C Θ respectively.  \nThe stability of the periodic orbits of the asteroid is guaranteed in linear approximation if the sequential principal minor ¯ A Θ > 0 and  \nD Θ = det [ ¯ A Θ ¯ B Θ ¯ B Θ ¯ C Θ ] > 0 . (3.18)  \nThe stability of the periodic orbits of the asteroid will be proved analytically with some conditions and then shown numerically for the general case with the help of Matlab.", '4 Limiting cases': "When the distance between the asteroid and the star is much smaller than between the planet and the star, an expansion over the ratio between a and a J can be performed. This is called the inner case of the double-averaged spatial ER3BP. Huang et al. [2024] considered the double-averaged value of the coefficients in the quadratic form of p 3 and q 3 and averaged W Θ with the same procedure in (3.17). The sequential principal minor  \n¯ A Θ , inner ≈ 3 ( 1 + 4 e 2 -5 e 2 cos (Θ) 2 ) 4 G (1 -e 2 J ) 3 2 a 2 , D Θ , inner ≈ 9 (1 + 3 e 2 -4 e 4 ) 16 G 2 (1 -e 2 J ) 3 a 4 , (4.1)  \nare calculated to be positive for a is small and e J ↛ 1. The orbits of the asteroid in the double-averaged model are stable in linear approximation for the inner case.  \nOn the contrary, when the distance between the asteroid and the star is much larger than the distance between the planet and the star (outer case), an expansion of d = 1 /a can be performed analytically. Taking  \nM 1 = x 2 J + y 2 J x 2 + y 2 + z 2 , M 2 = -2 ( xx J + yy J ) x 2 + y 2 + z 2 , (4.2)  \nand  \nM = M 1 + M 2 , (4.3)  \nthe approximate formula of the force function of gravity U = U Θ , outer of the planet up to a 3 is expanded as  \nU Θ , outer = -1 √ x 2 + y 2 + z 2 ( 1 -1 2 M + 3 8 ( M 1 M 2 + M 2 2 ) -5 16 M 3 2 ) . (4.4)  \nThen, the double-averaged value of the coefficients in the quadratic form of p 3 and q 3 is  \n¯ W Θ , outer = 1 2 ( ¯ A Θ , outer p 3 2 +2 ¯ B Θ , outer p 3 q 3 + ¯ C Θ , outer q 3 2 ) . (4.5)  \nAveraging over mean anomalies is difficult. By Kepler's equation  \nl = E -e sin E, l J = E J -e sin E J , (4.6)  \nthe double-averaged value of ¯ W Θ , outer can be derived by  \n¯ W Θ , outer = 1 4 π 2 ∫∫ [0 , 2 π ] 2 W Θ , outer d l d l J = 1 4 π 2 ∫∫ [0 , 2 π ] 2 W Θ , outer d l d E d l J d E J d E d E J = 1 4 π 2 ∫∫ [0 , 2 π ] 2 W Θ , outer (1 -e cos E ) (1 -e J cos E J ) d E d E J . (4.7)  \n¯ A Θ , outer , ¯ B Θ , outer and ¯ C Θ , outer are derived similarly. Substituting x , y , x J , y J into the doubleaveraged value of the coefficients, the coefficients up to order d 4 are  \n¯ A Θ , outer = 3 (1 -e 2 J ) 4 (1 -e 2 ) 3 2 G d 3 -75 e e J cos (Θ) (1 -e 2 J ) 32 (1 -e 2 ) 5 2 G d 4 + O ( d 5 ) , ¯ B Θ , outer = 15 e e J sin (Θ) (17 e 2 J -24) 128(1 -e 2 ) 5 2 G d 4 + O ( d 5 ) , ¯ C Θ , outer = 3 (1 + 4 e 2 J ) 4 (1 -e 2 ) 3 2 G d 3 -15 e e J cos (Θ) (43 e 2 J +34) 64 (1 -e 2 ) 5 2 G d 4 + O ( d 5 ) . (4.8)  \nIn the above formulas, as d is small enough, for e ↛ 1, all the coefficients of d 4 in (4.8) are bounded, thus terms of d 4 , as well as higher order terms, are omitted. It is obvious that ¯ A Θ , outer > 0 and ¯ C Θ , outer > 0. The determinant of the quadratic form (4.5) up to order a 7 is  \nD Θ , outer = 9 (1 -e 2 J ) (1 + 4 e 2 J ) 16 (1 -e 2 ) 3 G 2 d 6 + 45 e (44 e J +39 e 3 J -83 e 5 J ) cos (Θ) 256 (1 -e 2 ) 4 G 2 d 7 (4.9)  \nSimilarly, when d is small, the coefficient of d 7 in (4.9) is bounded as far as e ↛ 1. For small values of d , the sequential principal minor ¯ A Θ , outer > 0 and D Θ , outer > 0, thus (4.5) is a positive definite quadratic form. Consequently, for the outer case, periodic orbits of the double-averaged planar elliptic restricted 3-body problem are stable in the linear approximation with spatial perturbation.", '5 General case': "When the ratio between a and a J = 1 is arbitrary (general case), the orbits of the asteroid and the planet are assumed not in collision, i.e.,  \na (1 + e ) < 1 -e J for 0 < a < 1 , a (1 -e ) > 1 + e J for a > 1 . (5.1)  \nFor periodic orbits around the equilibria in Vashkovyak [1982], the orbits of the asteroid sway up and down periodically. When the sum of the argument of the pericenter and the ascending node of the asteroid is Θ = ω +Ω, denote the distance between the asteroid and the planet by  \nd 3 Θ ( l, l J ) = [ ( x -x J ) 2 +( y -y J ) 2 ] 3 / 2 . (5.2)  \nThen substituting (3.12) into (3.15), we get  \nyy 1  \nΘ Θ  \nA Θ = J Gd 3 Θ ( l, l J ) = Gd 3 Θ ( l, l J ) [sin(Θ) x ' y J +cos(Θ) y ' y J ] , B Θ = ( xy J + yx J ) 2 Gd 3 Θ ( l, l J ) = 1 2 Gd 3 Θ ( l, l J ) [ sin(Θ) x ' x J +cos(Θ) y ' x J +cos(Θ) x ' y J -sin(Θ) y ' y J ] , C Θ = xx J Gd 3 ( l, l J ) = 1 Gd 3 ( l, l J ) [cos(Θ) x ' x J -sin(Θ) y ' x J ] . (5.3)  \nThe double-averaged value of the quadratic form of p 3 and q 3 is  \n¯ W Θ = 1 2 ( ¯ A Θ p 3 2 +2 ¯ B Θ p 3 q 3 + ¯ C Θ q 3 2 ) , (5.4)  \nwhere  \n¯ A Θ = 1 4 π 2 G ∫∫ [0 , 2 π ] 2 1 d 3 Θ ( l, l J ) [sin(Θ) x ' y J +cos(Θ) y ' y J ] d l d l J , ¯ B Θ = 1 8 π 2 G ∫∫ [0 , 2 π ] 2 1 d 3 Θ ( l, l J ) [ sin(Θ) x ' x J +cos(Θ) y ' x J +cos(Θ) x ' y J -sin(Θ) y ' y J ] d l d l J , ¯ C Θ = 1 4 π 2 G ∫∫ [0 , 2 π ] 2 1 d 3 Θ ( l, l J ) [cos(Θ) x ' x J -sin(Θ) y ' x J ] d l d l J (5.5)", '5.1 General case with small Θ': "In the real cosmos, the Lidov-Kozai effects mostly occur with small changes in eccentricities, thus it is important to consider the small periodic orbits around the equilibrium in Vashkovyak [1982]. In this case, it is reasonable to assume that Θ is small, then  \nd 3 Θ ( l, l J ) = [ ( x -x J ) 2 +( y -y J ) 2 ] 3 / 2 = d ' 3 Θ ( l, l J ) + 3 d ' Θ ( l, l J ) ( y ' x J -x ' y J ) Θ + O (Θ 2 ) , (5.6)  \nwhere  \nd ' Θ ( l, l J ) = √ ( x ' -x J ) 2 +( y ' -y J ) 2 . (5.7)  \nThis subsection will prove the positive definiteness of the sequential principal minor of the quadratic form ¯ W Θ with small Θ. The proof is under the assumption that the orbit of the asteroid is inside the planet's orbit. The outside case can be proved similarly. Firstly, let us show  \nand  \n∫∫ [0 , 2 π ] 2 x ' y J d 3 Θ ( l, l J ) dldl J ∼ O (Θ) , (5.8)  \n∫∫ [0 , 2 π ] 2 y ' x J d 3 Θ ( l, l J ) dldl J ∼ O (Θ) . (5.9)  \nFigure 4: New orbit and symmetry of orbit.  \n<!-- image -->  \nThe orbit of the asteroid A concerning ( x ' , y ' ) can be regarded as an orbit of A 0 with x = x ' and y = y ' , in detail, the blue orbit of A 0 can replace the orange orbit of A in Fig. 4a. It can be found that the new orbit (orbit of A 0 ) and the orbit of the planet are symmetrical about the y -axis, See Fig 4b, then  \nd ' Θ ( l, l J ) = d ' -Θ (2 π -l, 2 π -l J ) , d ' Θ (2 π -l, l J ) = d ' -Θ ( l, 2 π -l J ) . (5.10)  \nAlso, distances with expansions of Θ are  \nd 3 Θ ( l, l J ) = d ' 3 Θ ( l, l J ) + 3 d ' Θ ( l, l J ) ( y ' x J -x ' y J ) Θ , d 3 Θ (2 π -l, 2 π -l J ) = d ' 3 Θ ( l, l J ) -3 d ' Θ ( l, l J ) ( y ' x J -x ' y J ) Θ , d 3 Θ ( l, 2 π -l J ) = d ' 3 Θ ( l, 2 π -l J ) + 3 d ' Θ ( l, 2 π -l J ) ( y ' x J + x ' y J ) Θ , d 3 Θ (2 π -l, l J ) = d ' 3 Θ ( l, 2 π -l J ) -3 d ' Θ ( l, 2 π -l J ) ( y ' x J + x ' y J ) Θ . (5.11)  \nThen formula (5.8) is proved by  \n∫∫ [0 , 2 π ] 2 x ' y J d 3 Θ ( l, l J ) dldl J = ∫∫ [0 ,π ] 2 [ x ' y J d 3 Θ ( l, l J ) -x ' y J d 3 Θ (2 π -l, 2 π -l J ) + x ' y J d 3 Θ (2 π -l, l J ) -x ' y J d 3 Θ ( l, 2 π -l J ) ] dldl J = ∫∫ [0 ,π ] 2 x ' y J [ d 3 Θ (2 π -l, 2 π -l J ) -d 3 Θ ( l, l J ) d 3 Θ ( l, l J ) d 3 Θ (2 π -l, 2 π -l J ) + d 3 Θ ( l, 2 π -l J ) -d 3 Θ (2 π -l, l J ) d 3 Θ (2 π -l, l J ) d 3 Θ ( l, 2 π -l J ) ] dldl J =6Θ ∫∫ [0 ,π ] 2 x ' y J [ -d ' Θ ( l, l J ) ( y ' x J -x ' y J ) d 3 Θ ( l, l J ) d 3 Θ (2 π -l, 2 π -l J ) + d ' Θ ( l, 2 π -l J ) ( y ' x J + x ' y J ) d 3 Θ (2 π -l, l J ) d 3 Θ ( l, 2 π -l J ) ] dldl J ∼ O (Θ) .  \nSimilarly,  \n∫∫ [0 , 2 π ] 2 y ' x J d 3 Θ ( l, l J ) dldl J =6Θ ∫∫ [0 ,π ] 2 y ' x J [ -d ' Θ ( l, l J ) ( y ' x J -x ' y J ) d 3 Θ ( l, l J ) d 3 Θ (2 π -l, 2 π -l J ) + d ' Θ ( l, 2 π -l J ) ( y ' x J + x ' y J ) d 3 Θ (2 π -l, l J ) d 3 Θ ( l, 2 π -l J ) ] dldl J ∼ O (Θ) ,  \nwhich proved formula (5.9). Consequently, for small Θ,  \n¯ B Θ = Θ 8 π 2 G ∫∫ [0 , 2 π ] 2 x ' x J -y ' y J d 3 Θ ( l, l J ) d l d l J + cos(Θ) 8 π 2 G ∫∫ [0 , 2 π ] 2 y ' x J + x ' y J d 3 Θ ( l, l J ) d l d l J = Θ 8 π 2 G ∫∫ [0 , 2 π ] 2 x ' x J -y ' y J d 3 Θ ( l, l J ) d l d l J + O (Θ) 8 π 2 G ∼ O (Θ) . (5.12)  \nThus, ¯ B Θ is of order Θ in the case of small Θ. Numerical checks of ¯ B Θ in Matlab for some values of a , e J are performed with variable Θ and e taken from 0 to 0 . 1 and 0 to 1 respectively. In such a way we verified that ¯ B is small for Θ ∈ [0 , 0 . 1].  \nFigure 5: Value of ¯ B Θ for small Θ.  \n<!-- image -->  \nBy formula (5.8), the double-averaged value ¯ A Θ is  \n¯ A Θ = cos(Θ) 4 π 2 G ∫∫ [0 , 2 π ] 2 y ' y J d 3 Θ ( l, l J ) d l d l J + O (Θ 2 ) 4 π 2 G = 1 2 π 2 G ∫∫ [0 ,π ] 2 [ d 3 Θ ( l, 2 π -l J ) -d 3 Θ ( l, l J ) d 3 Θ ( l, l J ) d 3 Θ ( l, 2 π -l J ) ] y ' y J dldl J . (5.13)  \nFor l, l J ∈ (0 , π ), it is obvious that  \nd 3 Θ ( l, 2 π -l J ) -d 3 Θ ( l, l J ) d 3 Θ ( l, l J ) d 3 Θ ( l, 2 π -l J ) > 0 . (5.14)  \nThus ¯ A Θ is positive. Numerical results of A in Matlab with a = 0 . 3, e J = 0 . 2 and a = 4, e J = 0 . 4 coincide with the analytical conclusion.  \nFigure 6: Value of ¯ A Θ for small Θ.  \n<!-- image -->  \nBy formula (5.9), the double-averaged value ¯ C Θ is  \n¯ C Θ = cos(Θ) 4 π 2 G ∫∫ [0 , 2 π ] 2 x ' x J d 3 Θ ( l, l J ) d l d l J + O (Θ 2 ) 4 π 2 G = 1 2 π 2 G ∫∫ [0 ,π ] 2 x ' x J d 3 Θ ( l, l J ) dldl J . (5.15)  \nTherefore, proving formula (5.15) is positive is equivalent to proving  \n∫∫ [0 ,π ] 2 xx J d 3 Θ ( l, l J ) d l d l J > 0 , (5.16)  \n̸  \nwhere x = x ' , i.e. the orbit of asteroid A can be considered as an orbit of A 0 in apsidal alignment case. The case of e J = 0 is the circular problem proved in [Neishtadt, 1975], thus e J = 0 can be assumed.  \nFigure 7: Prove of C .  \n<!-- image -->  \nTaking x > -ae , which means A 0 ( x, y ) is on the right of the orbit of the asteroid, denoted by A 0 ∈ I R , then A ' 0 ( -2 ae -x, y ) is on the left, denoted by A ' 0 ∈ I L . Similarly,  \ntaking x > e J , then J ( x J , y J ) ∈ I R J and J ' ( -2 e J -x J , y J ) ∈ I L J . See Fig. 7. The distance A 0 J = d 1 , A ' 0 J ' = d 2 , A 0 J ' = d 3 and A ' 0 J = d 4 . It is shown in Fig. 7a, when d 2 < d 4 , d 2 < d 3 , d 1 < d 4 and d 1 < d 2 , the considered function  \nˆ C = xx J d 3 1 + ( -2 ae -x )( -2 e J -x J ) d 3 2 + x ( -2 e J -x J ) d 3 3 + ( -2 ae -x ) x J d 3 4 > xx J d 3 2 + ( -2 ae -x )( -2 e J -x J ) d 3 2 + x ( -2 e J -x J ) d 3 2 + ( -2 ae -x ) x J d 3 2 = 4 aee J d 3 2 > 0 . (5.17)  \nIf replacing d 2 < d 4 by d 2 > d 4 , see Fig. 7b, then by rearrangement inequality,  \nˆ C ≤ xx J d 3 1 + ( -2 ae -x )( -2 e J -x J ) d 3 4 + x ( -2 e J -x J ) d 3 3 + ( -2 ae -x ) x J d 3 2 > xx J d 3 4 + ( -2 ae -x )( -2 e J -x J ) d 3 4 + x ( -2 e J -x J ) d 3 4 + ( -2 ae -x ) x J d 3 4 = 4 aee J d 3 4 > 0 . (5.18)  \nFor other cases, one can prove them by rearrangement inequality directly. Thus, the integral  \n∫∫ [0 ,π ] 2 xx J d 3 Θ ( l, l J ) d l d l J = [ ∫∫ I L × I L J + ∫∫ I L × I R J + ∫∫ I R × I L J + ∫∫ I R × I R J ] xx J d 3 Θ ( l, l J ) d l d l J = ∫∫ I R × I R J ˆ C d l d l J > 0 (5.19)  \nAs a result, ¯ C Θ is positive. Meanwhile, this is also an analytical proof of variable C in Neishtadt et al. [2021] instead of the numerical work. Numerical calculations of ¯ C Θ were performed in Matlab for some values of a , e J . The variables Θ and e are taken from 0 to 0 . 1 and 0 to 1 respectively with some grids. In such a way we verified that ¯ C is always positive. Fig. 8 are some numerical results with a = 0 . 3, e J = 0 . 2 and a = 4, e J = 0 . 4.  \nSince ¯ B Θ ∼ O (Θ) is rather smaller than ¯ A Θ and ¯ C Θ , the determinant D Θ = ¯ A Θ ¯ C Θ -O (Θ 2 ) > 0. Numerical checks are shown in Fig. 9. Thus, ¯ A Θ > 0, D Θ > 0, and ¯ W Θ is a positive definite quadratic form. Hence, small periodic orbits around the equilibria of the double-averaged planar elliptic restricted 3-body problem are stable in the linear approximation as small periodic orbits of the double-averaged spatial perturbed elliptic restricted 3-body problem for small Θ.", '5.2 General case with arbitrary Θ': 'When Θ is not small, expansion of Θ can not be performed. Calculations of the coefficients of the quadratic form ¯ W Θ have been done numerically in Matlab. With the averaging  \n<!-- image -->  \nFigure 8: Value of ¯ C Θ for small Θ.  \n<!-- image -->  \n<!-- image -->  \nFigure 9: Value of D Θ for small Θ.  \n<!-- image -->  \nprocedure in (4.7), detailing the expression of double-averaged values of coefficients in the quadratic form:  \n¯ A Θ = 1 4 π 2 G ∫∫ [0 , 2 π ] 2 yy J (1 -e cos E ) (1 -e J cos E J ) G [ ( x -x J ) 2 +( y -y J ) 2 ] 3 / 2 d E d E J , ¯ B Θ = 1 8 π 2 G ∫∫ [0 , 2 π ] 2 ( xy J + yx J ) (1 -e cos E ) (1 -e J cos E J ) G [ ( x -x J ) 2 +( y -y J ) 2 ] 3 / 2 d E d E J , ¯ C Θ = 1 4 π 2 G ∫∫ [0 , 2 π ] 2 xx J (1 -e cos E ) (1 -e J cos E J ) G [ ( x -x J ) 2 +( y -y J ) 2 ] 3 / 2 d E d E J . (5.20)  \nTaking a J = 1, for chosen values of a and e J , each pair of (Θ , e ) belongs to one certain periodic orbit which is determined by  \n¯ R = -1 4 π 2 ∫∫ [0 , 2 π ] 2 (1 -e cos E ) (1 -e J cos E J ) √ ( x -x J ) 2 +( y -y J ) 2 d E d E J . (5.21)  \nConsidering values of variables satisfying conditions (5.1), and then plotting the figures of (Θ , e, ¯ A Θ ) and (Θ , e, D Θ ) numerically for some considerable a and e J . The numerical work takes a = 0 . 3 and e J = 0 . 2 in Fig. 10, a = 0 . 2 and e J = 0 . 4 in Fig. 11, a = 4 and e J = 0 . 4 in Fig. 12, a = 10 and e J = 0 . 1 in Fig. 13. All the numerical results demonstrate that the sequential principal minor ¯ A Θ > 0 and D Θ > 0, thus the quadratic form ¯ W Θ is positive defined. Therefore, the numerical calculation gives stability of the periodic orbits of the asteroid in the linear approximation in the double-averaged spatial perturbed elliptic restricted 3-body problem for all values of parameters.  \nFigure 10: Value of ¯ A Θ and D Θ for a = 0 . 3 and e J = 0 . 2.  \n<!-- image -->', '6 Conclusion': "A complete analysis of secular effects on the motion of a massless asteroid within the framework of the spatial perturbed, double-averaged elliptic restricted three-body problem has been conducted. The stability of the asteroid's orbits was investigated with linearization.  \nFigure 11: Value of ¯ A Θ and D Θ for a = 0 . 2 and e J = 0 . 4.  \n<!-- image -->  \nFigure 12: Value of ¯ A Θ and D Θ for a = 4 and e J = 0 . 4.  \n<!-- image -->  \nFigure 13: Value of ¯ A Θ and D Θ for a = 10 and e J = 0 . 1.  \n<!-- image -->  \nNotably, periodic orbits originating from the planar problems are stable within the spatial perturbation across all parameter regions. Numerical simulations corroborate these findings for a variety of periodic orbits. The model's applicability to systems with highly eccentric planets renders it particularly valuable for exoplanet studies.", 'Acknowledgements': 'The authors express their gratitude to Prof. Anatoly Neishtadt for suggestions on the topic and to Prof. Xijun Hu for discussions. Kaicheng Sheng thanks the National Natural Science Foundation of China (NSFC) for the support of this research (Grant: 12371192 & 12271300).'}