Numerical-Analytical Study of Linked Orbits in the Restricted Elliptic Doubly Averaged Three-Body Problem

The restricted elliptic doubly averaged three-body problem is considered. The terms up to the second order inclusive with respect to the orbital eccentricity of the perturbing body are retained in the expansion of the perturbing function. The elements of the elliptical orbit for the perturbed body are deemed arbitrary. The special, so-called linked orbits of a negligible-mass body are investigated by numerically integrating the averaged equations in Keplerian elements. For such orbits the points of their intersection with the orbital plane of the perturbing body are on different sides of this orbit. The evolution of hypothetical and some real cometary orbits is described in the simplest Sun-Jupiter-comet model; their differences from the corresponding orbits in the circular problem have been revealed.


INTRODUCTION AND FORMULATION
OF THE PROBLEM Studies of the long-term evolution of orbits in the restricted elliptic three-body problem are generally carried out in the doubly averaged formulation. The integrable case of a circular orbit of the perturbing body is widely used. The studies of this case initiated by the renowned scientists H. von Zeipel (von Zeipel, 1910) and N.D. Moiseev (Moiseev, 1945) were elaborated significantly by M.L. Lidov (Lidov, 1961(Lidov, , 1962 and I. Kozai (Kozai, 1962). These studies are described in detail in the monographs by Shevchenko (2017) and Ito and Ohtsuka (2019). Note that the extensive paper by H. von Zeipel has become deservedly famous and is reflected in the above scientifichistorical study by Ito and Ohtsuka (2019) only due to the reference to it in Baily and Emel'yanenko (1996) in connection with the study of the evolution of one of the types of cometary orbits. The works by Moiseev, Lidov, and Kozai performed later were both the results of a qualitative study of the averaged three-body problem and the necessity of investigating the orbital dynamics of artificial satellites of planets and the dynamics of asteroids.
In his extensive paper Zeipel (1910) identified and qualitatively studied three main cases of the orbit location for the perturbed body in the doubly averaged circular problem: the inner one, the outer one, and the case of intersecting, in particular, the so-called linked orbits. The spatial location of linked orbits is characterized by the fact that one of the points of intersection of the orbit of a negligible-mass body with the plane of the orbit of the perturbing body is located inside it, while the other one is outside. Such a classification in the restricted circular doubly averaged three-body problem for uniformly close orbits, along with an analysis of the conditions for their intersection, was proposed by Lidov and Ziglin (1974). The topology of two linked and unlinked Keplerian orbits of all those types was described in detail by Kholshevnikov and Titov (2007).
Note that the phrase "like the rings of a chain" proposed in the monograph by Ito and Ohtsuka (2019) as the English analog of the French term "comme les anneaux d'une shaine" used by von Zeipel (1910) is more sensible and geometrically understandable than the term "linked orbits". The linked orbits in the restricted elliptic doubly averaged three-body problem are the subject of our study. Consider the motion of particle Р of negligible mass under the attraction of a central point S of mass m and a perturbing point J of mass m 1 m moving relative to S in an elliptical orbit with a semimajor axis а 1 and eccentricity е 1 . Let us introduce a rectangular coordinate system Oxyz with the origin at point S whose reference plane xOy coincides with the orbital plane of point J. Let the Ox axis be directed to the pericenter of the orbit of point J, the Oy axis be in the direction of its motion from the pericenter in the reference plane, and the Oz axis complements the coordinate system to a right-handed one. The perturbed orbit of point Р is characterized by osculating Keple-rian elements: the semimajor axis а, the eccentricity е, the inclination i, the argument of pericenter ω, and the longitude of the ascending node Ω. In the chosen coordinate system the perturbing point J has coordinates x 1 , y 1 , and z 1 = 0.
The secular part W of the complete perturbing function is used to investigate the orbital evolution of point Р: (1) Here, is the distance between the perturbed and perturbing points, λ and λ 1 are the mean longitudes of these points, and f is the gravitational constant. The absence of low-order commensurabilities between the mean motions of points J and Р is assumed. In the function W a 1 and e 1 act as parameters of the evolution problem.
The first integrals of the equations of perturbed motion in elements are (2) while one more first integral exists in the case of е 1 = 0 (Moiseev, 1945): (3) THE AVERAGED PERTURBING FUNCTION AND EVOLUTION EQUATIONS Another expression, equivalent to (1), for the function W via well-known formulas is also commonly used in analytical studies: where V is the attractive force function of the elliptical Gaussian ring simulating the averaged influence of the perturbing point and Е is the eccentric anomaly of point Р.
Since for the orbits under consideration the distance r can be both less than and greater than r 1 as they evolve, the commonly used expansions of the inverse distance 1/Δ into series in Legendre polynomials are inapplicable. Therefore, in this paper we will use the analytical expression for the function V, though with a limited accuracy up to inclusive, from Vashkov'yak (1986) for a nearly coplanar system of N of weakly elliptical Gaussian rings, but valid for any relation between r and r 1 . Taking into account the orientation of the introduced coordinate system and assuming that N = 1, i 1 = ω 1 = Ω 1 = k 1 = u 1 = v 1 = 0, and h 1 = e 1 in Eqs. (6)-(8) of this paper, we will obtain (5) and, for coherence, we will permit ourselves to also reproduce the basic simplified computational formulas from this paper by supplementing them with new analytical relations.
The function Ф dependent on the rectangular coordinates x, y, z is expressed via hypergeometric Gaussian functions F:  Here, The constant coefficients B n and H n are defined by the recurrence relations (9) while the empirical value of ζ* is taken to be 0.5.
Remark. The function V depends only on the squares of the coordinates y and z, while the coordinate x enters into this function both quadratically and linearly (into the numerator of ε and via it into ζ). The existence of two planar particular solutions, y = 0 and z = 0, in the singly averaged (only over λ 1 ) evolution problem follows from such a double symmetry of the function V with respect to y and z.
The rectangular coordinates x, y, z are expressed via Е by the well-known formulas for unperturbed Keplerian motion (10) Below, it will be convenient to introduce a new independent variable, a "dimensionless time" τ, according to the formula cos cos sin sin cos , cos sin sin cos cos , sin sin , sin cos cos sin cos , sin sin cos cos cos , cos sin .
is the mean motion of point Р, and a normalized perturbing function (13) To describe the evolution of orbits, we will use the Lagrange equations in elements with the function w that is their first and unique integral: The existence of stationary solutions of these equations is possible if the conditions are fulfilled.

PARTIAL DERIVATIVES OF THE NORMALIZED FUNCTION w WITH RESPECT TO THE ELEMENTS
Generally, for arbitrary orbits of point Р a solution of Eqs. (14) can be found apparently only by a numerical method, while the process of calculations can be controlled by the constancy of the function w along this solution. In Vashkov'yak (1986) the partial derivatives of the function w with respect to the elements were calculated by a difference method. Here, we use a combined method, in which the quadratures (15) are found numerically by the Gaussian method, while the derivatives of the normalized function (16) are found analytically. For the completeness of the set of formulas, we give the necessary expressions to calculate the derivatives of the function with respect to the orbital elements: (20) Relations (7)-(11) and (15)-(22) represent a complete set of formulas to calculate the right-hand sides of the evolution equations (14).
In view of the properties of the function V (see the remark in the previous section), system (14) has two particular solutions or integrable cases.
Case 1. If sini = 0, then the orbital plane of point Р coincides with the orbital plane of the perturbing point J. By symmetry, it turns out that di/dτ = 0. However, in this planar solution at any а, а 1 , and e 1 , in addition to regular orbits, there exist irregular ones intersecting (but not linked) with the orbit of point J (Vashkov'yak, 1982).
Case 2. If cos i = 0 and sinΩ = 0, then in the elliptic problem under consideration the orbital plane of point Р is orthogonal to the orbital plane of the perturbing point J and passes through its line of apsides. By symmetry, it turns out that di/dτ = 0 and dΩ/dτ = 0. The above formulas of the quadratic approximation in е 1 It is easy to verify that at y = 0 the derivatives of the functions ε, ζ, μ, and ν with respect to y become zero, so that and Equations (14) simplified for this case of orthogonal-apsidal orbits take the form (25) A general qualitative study of this case by taking into account the possible intersections of the orbits of points Р and J was carried out by a numerical-analytical method in Vashkov'yak (1984) for arbitrary a, а 1 , and е 1 . In this paper greater attention is given to linked orbits and, in particular, to the stationary solutions of Eqs. (25) existing at ω 0 = 0, π. It can be shown that, in this case, and the stationary eccentricities themselves are determined as the roots of the transcendental equation ON THE EVOLUTION OF SOME HYPOTHETICAL AND REAL COMETARY ORBITS IN THE SOLAR SYSTEM MODEL First we will turn to the linked orbits of point Р highly inclined to the reference plane. In the integra-  , , , , 0, ; 0, 0. w a a e e e ble case of orthogonal-apsidal orbits, a numerical solution of Eq. (26) makes it possible to find the stationary values of the eccentricity е 0 at given ω 0 = 0°, 180°, Ω 0 = 180° and fixed parameters а, а 1 , е 1 . In addition, it can be shown that the values of е 0 depend on ω 0 and Ω 0 only through the combination δ 1 = sign(cosω 0 cosΩ 0 ) = ±1. However, the difference in е 0 at different signs of δ 1 is fairly small and ~.
At an inclination and longitude of the ascending node different from those adopted in case 2, all of the orbital elements will change with time. Our numerical integration of the evolution system (14) by the Runge-Kutta method at а 1 = 5.2 AU, е 1 = 0.048, and the mass ratio m/m 1 in the Sun-Jupiter system makes it possible to estimate such changes for fictitious (or hypothetical) cometary orbits. Tables 1 and 2 give such an estimate in the time interval t* = 500000 years for orbits with a semimajor axis а = 10а 1 = 52 AU, e 0 (δ 1 = 1) = 0.9890, е 0 (δ 1 = -1) = 0.9905. These tables present Table 1 was compiled for three values of i 0 = 90°, 85°, 95°. In the exact solution of Eqs. (14) i 0 = 90°, Ω 0 = 0°, 180° the deviations are zero and, therefore, the corresponding rows are omitted. At Ω 0 = 0 and 180° the results for i 0 = 95° identical to the corresponding data for i 0 = 85° are also omitted.  The initial deviations in i 0 and Ω 0 from their equilibrium values at t = t* lead to insignificant changes in the shape of the orbit (Δ 1 , Δ 3 ), but to a noticeable change in its orientation (Δ 2 reaches 24°, Δ 4 is about 32°). Table 2 was compiled for more significant initial deviations of the orbit from the orthogonal one, i 0 = 90° ± 30°. At Ω 0 = 0° and 180° it also contains the results for i 0 = 120° identical to the corresponding data for i 0 = 60°.
In this case, the inclination also remains close to the initial one, but the deviations of the remaining elements are tens of degrees, reaching approximately 75°.
Below, we will consider the orbits of point Р linked with the orbit of point J, but with an arbitrary spatial orientation. The evolution of such orbits is usually studied using numerical methods even in the integrable doubly averaged circular problem (е 1 = 0) due to the absence of a rigorous analytical expression for the averaged perturbing function. The families of phase trajectories in the (еcosω, esinω) plane constructed as equipotential contours of the function W are presented in Ito and Ohtsuka (2019, Section 5.8, Fig. 24). These families correspond to the hypothetical linked orbits of point Р for three pairs of values of the ratio a/a 1 and the constant of the integral с 1 . In all three cases considered there exist stationary center-type singular points and closed periodic trajectories enclosing them in the phase plane.
The ellipticity of the orbit of the perturbing point J, naturally, leads to qualitative changes of the families of trajectories in the circular problem. Since the equa-

I I I III
tions for е and ω are not decoupled from the remaining ones, as is the case at е 1 = 0, due to the absence of the integral с 1 in the elliptic problem, the evolution of these elements can be traced only in projection of the phase trajectory onto the (еcosω, esinω) or (ω, е) plane. In this paper we numerically integrated system (14) for а 1 = 5.2 AU, е 1 = 0.048, and the mass ratio m/m 1 in the Sun-Jupiter system. The difference of this ratio from the one adopted in the calculations by Ito and Ohtsuka (2019) for the Sun-(Earth + Moon) system does not affect the structure of the phase trajectoriesне, but leads only to a change in the time scale.
For a comparison with the results of the circular problem in cases I, II, and III, out of all families of integral curves in the circular problem we chose the trajectories with ω 0 = 180° and initial е 0 = 0.55, 0.4, and 0.65, respectively. The initial inclinations were calculated from е 0 and с 1 as , while the initial longitude of the ascending node Ω 0 was taken to be zero. Figures 1-3 show the trajectories for cases I, II, and III, respectively, and fragments of the polar diagrams with plotted numerical values of the angles ω and radii е. The circles and triangles mark the initial and final points, respectively. The time intervals are 100000 years for cases I, II and 500000 years for case III. The dashed lines are the so-called "separatrices", which are not integral curves and correspond to the intersections of the orbits of points P and J. In all three cases, the closed periodic trajectories of the circular problem are modified and become nonperiodic, nevertheless, retaining an oscillatory pattern and remaining in the regions of linked orbits. The amplitude of these oscillations can both decrease (Figs. 1, 3) and increase ( Fig. 2) with time.
However, the eccentricity of Jupiter's orbit can also lead to qualitative changes in the behavior of the trajectory. Figure 4 shows a chaotic trajectory that begins in the ω libration region at ω 0 = 180°, then passes into the circulation region, and returns to the libration region, but relative to ω = 0. The initially linked orbit of point Р intersects the orbit of the perturbing point J as it evolves, exiting from one linking region, then traverses the region of unlinked orbits and enters another linking region.
Interestingly, real cometary orbits almost orthogonal to the ecliptic, with the above types of evolution, are also found within the model under consideration (Sun-Jupiter-comet). If we make a selection by perihelion distance q < 2 AU and inclination 85° < i < 95°o n the set of cometary orbits presented in the JPL database (https://ssd.jpl.nasa.gov/sbdb_query.cgi#x), then only ten orbits remain in this sample. The evolution of seven of them is reduced to successive intersections of Jupiter's orbit with passages of the linking regions. Figures 5-7 give an idea of the evolution of the remaining three orbits. In contrast to the previous figures, fragments of the (ω, е) plane are shown here not in the polar coordinate system, but in the rectangular one, which is more convenient for highly elliptical orbits. The up-to-date (initial) values of the elements for our numerical integration were taken from the above JPL database. The initial and final values in the (ω -е) plane are marked by the circles and triangles, respectively.  26) is е 0 = 0.9818. The phase point has no time to make one revolution relative to the libration center, while the trajectory passes through the "separatrix" corresponding to the intersection of the cometary orbit with Jupiter's orbit. After the relatively short time interval of ~200000 years, there occur the second intersection of these orbits and the entry into another linking region with a ω libration center offset by 180° from the initial one. In this case, the initially prograde motion of the comet becomes retrograde already in 100000 years, while the inclination, increasing monotonically, reaches about 115°. Figure 6 presents the variations in orbital elements for comet C/1861 J1 (Great comet) in a time interval of 3 Myr. The motion of the phase point begins and remains in the region of orbital linking and ω libration, while the corresponding stationary value of the eccentricity is е 0 = 0.9820. In this case, the initially prograde motion of the comet becomes retrograde in 2 Myr, while the inclination reaches about 132° at a minimum of 31°. Figure 7 presents the variations in orbital elements for comet 122P/de Vico in a time interval of 3 Myr. The motion of the phase point also begins and remains in the region of orbital linking and ω libration, while the corresponding stationary value of the eccentricity In conclusion, it should be recalled that all the above examples of orbital evolution and, in particular, the examples of libration of the argument of perihelion were constructed within the adopted model (Sun-Jupiter-comet). However, the real motions of comets in the Solar System can differ from their model motions even qualitatively. These differences are due to the neglect of the influence of the remaining plan- ets, except for Jupiter, and the averaged model of the problem. For example, in the rigorous solution of the complete system of non-averaged differential equations the libration of the argument of perihelion can change to its circulation. The orbit of comet 122P/de Vico is an example of chaotic orbital evolution in a real cometary medium, but this property is found only in the non-averaged Solar System model including both several perturbing bodies (Baily et al. 1992) and one body (Ito and Ohtsuka, 2019; Section 5.8, Fig. 25).
Regarding the methodology of the work described in this paper, note that owing to the studies performed relatively recently and reflected in the monographs and the paper by Kondrat'ev (2007Kondrat'ev ( , 2012 and Antonov et al. (2008), the representation of the force function for an elliptical Gaussian ring was obtained in a closed form without expansions in terms of any parameters. However, its practical use involves wellknown difficulties and may in future be the subject of a special study.

OPEN ACCESS
This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.