Quasi-satellite orbits in the general context of dynamics in the 1:1 mean motion resonance: perturbative treatment

Celestial Mechanics and Dynamical Astronomy manuscript No. (will be inserted by the editor) Quasi-satellite orbits in the general context of dynamics in the 1:1 mean motion resonance. Perturbative treatment ⋆ Vladislav V. Sidorenko · Anatoly I. Neishtadt · Anton V. Artemyev · Lev M. Zelenyi Abstract Our investigation is motivated by the recent discovery of asteroids orbiting the Sun and simultaneously staying near one of the Solar System planets for a long time. This regime of motion is usually called the quasi-satellite regime, since even at the times of the closest approaches the distance between the asteroid and the planet is significantly larger than the region of space (the Hill’s sphere) in which the planet can hold its satellites. We explore the properties of the quasisatellite regimes in the context of the spatial restricted circular three-body problem “Sun-planet-asteroid”. Via double numerical averaging, we construct evolutionary equations which describe the long-term behaviour of the orbital elements of an asteroid. Special attention is paid to possible transitions between the motion in a quasi-satellite orbit and the one in another type of orbits available in the 1:1 resonance. A rough classification of the corresponding evolutionary paths is given for an asteroid’s motion with a sufficiently small eccentricity and inclination. ⋆ Results of this paper were partially presented as Paper DDA 101.04 at the 44th Annual Meeting of the Division on Dynamical Astronomy of the American Astronomical Society, 2013, Paraty, Brazil. V.V.Sidorenko Keldysh Institute of Applied Mathematics Russian Academy of Sciences, Miusskaya Sq., 4, 125047 Moscow, RUSSIA Moscow Institute of Physics and Technology Institutskiy S-Str., 9, 141700 Dolgoprudny, RUSSIA E-mail: vvsidorenko@list.ru A.I.Neishtadt Loughborough University, LE11 3TU Loughborough, UK A.I.Neishtadt · A.V.Artemyev · L.M.Zelenyi Space Research Institute Russian Academy of Sciences, Profsoyuznaya Str., 84/32, 117997 Moscow, RUSSIA 2 V.V.Sidorenko et al. 1 Introduction During recent decades, the properties of the so-called quasi-satellite orbits (QSorbits) have been intensively explored. Within the scope of the three body problem “Sun-planet-asteroid”, the motion of an asteroid in a QS-orbit corresponds to the 1:1 mean motion resonance, with the resonant argument ϕ = λ − λ′ librating around 0 ( λ and λ′ being the mean longitudes of the asteroid and the planet, respectively). The asteroid motion in QS-orbit is bounded to the planet’s neighborhood of a size which can be small enough in comparison with the value of semimajor axis a′ of the planet (Fig. 1). Nevertheless, the trajectory of the asteroid will never cross the Hill sphere of the planet, wherefore the asteroid cannot be considered as a satellite in the usual sense. Fig. 1 The orbital motion of a quasi-satellite and its host planet. Panel a is a Sun-centered reference frame that preserves the orientation in the absolute space. The quasi-satellite and the planet move around the Sun with the same orbital period in elliptic and in circular orbits respectively. Panel b is a Sun-centered frame rotating with the mean orbital motion of the planet. To our knowledge, for the first time the existence of QS-orbits was discussed by Jackson (1913). Long enough the studies of this class of orbits were limited to the consideration of the periodic motions classified as “f-family” by Strömgren. As remarkable from the different points of view we can mention in this context the investigations by Broucke (1968), Henon (1969) and Benest (1977). Carrying out his research on the NASA contract Broucke (1968) actually anticipated the application of QS-orbits in astrodynamics. Henon (1969) proved the (planar) stability of f-family in the Hill approximation and conjectured the existence of the natural retrograde ”satellites” that are much farther from the host planet than the collinear libration points L1 and L2 . Benest (1977) established the conditions for three-dimensional stability of f-family periodic orbits in the frame CR3BP. At the end of 80th the opportunity to insert a spacecraft into QS-orbit around the Martian moon Phobos was thoroughly studied and finally realized in the former USSR (Kogan 1989; Lidov and Vashkovyak 1993, 1994). Quasi-satellite orbits. Perturbative treatment 3 Phobos was the goal of the last Soviet interplanetary mission (Sagdeev and Zakharov 1989). Since the Hill sphere of Phobos is very close to its surface it is impossible to circumnavigate it in a Keplerian-type way. Inspired by the results of Henon (1969) and Benest (1976) A.Yu.Kogan (then a mission specialist at Lavochkin Aerospace company, USSR) proposed a QS-orbit as a solution. Although this mission was only partially successful one of the two launched spacecrafts attained stable QS-orbit (A.Yu.Kogan, private communication). The related activity was reviewed in (Kogan 1990), where in particular the currently popular definition of QS-orbit was formulated probably for the first time: quasi-satellite orbits are the trajectories of restricted three-body problem which are (1) located far beyond the Hill’s sphere surrounding the minor primary body and (2) much less distant from it than from the major primary. Later the application of QS-orbits in astrodynamics was considered by Tuchin (2007), Gil and Schwartz (2010) and many other specialists. Since outside the Hill sphere the gravity field of the planet is weak enough, a QS-orbit can be treated as a slightly perturbed heliocentric Keplerian ellipse. It offers great opportunities for analytical consideration of the motion in QS-orbits. For the first time such a strategy was applied by Mikkola and Innanen (1997). A little bit later the perturbation approach was used to study the properties of OS-orbits in the papers by Namouni (1999) and Namouni et al. (1999). These papers provided the greatest progress in understanding of the key dynamical structures responsible for the long-term evolution at 1:1 mean motion resonance. New types of orbital behavior were described (in particular, the so called compound orbits). The terminology (except for the minor modifications) and the concepts introduced in (Namouni 1999) and (Namouni et al. 1999) became standard for further theoretical research of QS-orbits (Christou 2000; Brasser et al. 2004; Mikkola et al. 2006, etc.)1 . Below we will be dealing also with many of these concepts. In the last decade the interest to QS-orbits has increased due to the discovery of the actual quasi-satellites for Venus (Mikkola et al. 1997), Earth (Connors et al. 2004; Wajer 2009, 2010), Jupiter (Kinoshita and Nakai 2007) and Neptune (Fuente Marcos and Fuente Marcos 2012). The role of the quasi-satellite dynamics in the early Solar nebula was discussed by Kortenkamp (2005). Guipone et al. (2010) investigated the properties of the QS-motions in extrasolar planetary systems. An important phenomenon revealed by Namouni (1999) and Namouni et al. (1999) is a possibility for an asteroid in the 1:1 mean motion resonance to change, from time to time, the qualitative character of its orbital motion. In particular, under a special choice of the initial conditions, transitions between the motion in 1 Being evidently unaware about the studies of QS-orbits by specialists in astrodynamics, Namouni (1999) and Namouni et al. (1999) used for this type of the orbital motion the term ”retrograde satellite orbit” taken from the classical investigations of the periodic solutions in RC3BP. Looking through the literature it is easy to note that the community of celestial mechanicians is not uniform regarding what is more preferable. Since our activity was stimulated by the long standing discussions with specialists who coined the term ”QS-orbit” in 80th (A.Yu.Kogan, M.L.Lidov and M.A.Vashkovyak), it predetermined our choice. 4 V.V.Sidorenko et al. a QS-orbit and that in a horseshoe orbit (HS-orbit) take place. More complicated scenaria are possible too (Brasser et al. 2004; Namouni 1999; Namouni et al. 1999). To study the secular evolution of the resonant motion with qualitative changes in the behavior of the argument ϕ , one can apply an approach developed by Wisdom (1985) in his investigation of the 3:1 mean motion resonance. This approach is, in fact, general enough and contained no restrictions on the type of the resonant orbit to model (e.g., Yokoyama 1996). In essence it is based on the presence of the adiabatic invariant in the asteroid’s dynamics at the resonance. Calculating level curves of the adiabatic invariant one can draw phase portraits characterizing the secular evolution of the motion. The consideration of the adiabatic invariance violation (due to the transitional phenomena) allows to understand the appearance of the chaos in asteroid’s motion. The first step in application of Wisdom’s method to the 1:1 resonance was carried out in (Nesvorny et al. 2002). In the system’s phase space the resonant phenomena are localized: they occur in narrow resonance regions (Arnold et al. 2006). It always interesting to compare the properties of the resonant motions with the properties of non-resonant motions when corresponding phase trajectories are close to the border of the resonance region. In the case of 1:1 mean motion resonance a lot of information for such a comparison can be found in the papers on the non-resonant motions of the asteroid and the planet body in the close orbit (e.g., Lidov and Ziglin 1974; Gladman 1993). The goal of our paper is twofold: (i) to establish the conditions at which the motion in QS-orbit is possible, and (ii) to explore when this regime of orbital motion is perpetual and when it is temporary. Wisdom’s scheme of the meanmotion resonance analysis allows to notice what was not noticed in the previous studies on quasi-satellite dynamics. Section 2 begins with the description of an averaging procedure used to determine the secular evolution in a mean motion resonance. The first averaging is carried out over the orbital motion, whereafter the phase variables are rescaled, and the problem is shaped into a form called a “slow-fast” system (SF-system). This is a two degrees of freedom Hamiltonian system with the variables evolving at different rates: some variables are “slow”, while the other are “fast”. The second averaging is then performed over the “fast” motions of the SF-system. This provides us the evolutionary equations describing the secular effects in the asteroid’s motion. Section 3 is devoted to the analysis of these secular effects, for various regimes of motion. The transitions between different regimes of orbital motion (QS → HS, HS → QS, etc) are discussed. In Section 4, the consideration is restricted to motion in orbits with a small inclination and eccentricity. In this case, the asteroidal dynamics demonstrates some kind of scaling. Section 5 provides an example of the dynamics of an actual asteroid in a QS-orbit. In Section 6, the summary of the main results is presented. Details of the averaging procedures are elucidated in Appendices A and B. Quasi-satellite orbits. Perturbative treatment 5 2 Double averaged motion equations for investigation of dynamics at 1:1 mean-motion resonance 2.1 Averaging over orbital motions Through all stages of our analysis we shall use the motion equations written in the Hamiltonian form. The units are chosen so that the distance between the primaries (i.e., between the Sun and the planet) is unity, the sum of their masses is also unity, while the period of their rotation around the system’s barycenter is 2π . Since the mass of the planet µ is substantially smaller than the mass of the Sun, the quantity µ is a small parameter of the problem. To start with, we introduce the Delaunay canonical variables (Murray and Dermott 1999) L, G, H, l, g, h, where l is the mean anomaly of the asteroid, the other variables being related to the asteroid’s osculating elements (the semimajor axis a , the eccentricity e , the inclination i , the argument of pericenter ω , the longitude of the ascending node Ω ) by the formulae L= p (1 − µ)a, G=L g = ω, p 1 − e2 , H = G cos i, h = Ω. The canonical relations give the equations of motion of the asteroid: ∂K dL =− , dt ∂l dG ∂K =− , dt ∂g ∂K dl = , dt ∂L dg ∂K = , dt ∂G dH ∂K =− , dt ∂h dh ∂K = , dt ∂H with the Hamiltonian K defined as: (1 − µ)2 − µR(L, G, H, l, g, h − λ′ ), (1) 2 L2 R being the disturbing function. For the restricted circular three-body problem the disturbing function admits the form K=− R= 1 |r − r′ | − (r, r′ ) (2) with r = r(L, G, H, l, g, h) and r′ = r′ (λ′ ) being the position vectors of the asteroid and the planet relative to the Sun. The brackets around the second term in (2) denote the scalar product. The Hamiltonian K is a function of the time t through its dependence upon the planet’s mean longitude λ′ (in our units, λ′ = t + λ′0 ). As it follows from the formula (1), it is reasonable to replace the variable h with the variable e h = h−λ′ . e will be time-independent: The new system’s Hamiltonian K e = K(L, G, H, l, g, e K h) − H. Our next step is to perform the canonical transformation 6 V.V.Sidorenko et al. (L, G, H, l, g, e h) 7−→ (Pϕ , Pg , Ph , ϕ, g, h) (3) defined by the generating function S = Pϕ l + [Pg + Pϕ − 1]g + [Ph + Pϕ − 1]e h. The relations between new and old variables are: L = ∂S = Pϕ , ∂l ϕ = ∂S = l + e h + g, ∂Pϕ G = ∂S = Pg + Pϕ − 1, g = ∂S = g, ∂g ∂Pg h. H = ∂S = Ph + Pϕ − 1, h = ∂S = e ∂e h ∂Ph The purpose of transformation (3) is to introduce the resonant phase ϕ into the consideration. It is straightforward to verify that ϕ = l + (h − λ′ ) + g = λ − λ′ , where λ = l + h + g is the asteroid’s mean longitude. Once the transformation (3) will be accomplished, the Hamiltonian assumes the form of 2 e = − (1 − µ) − Ph − Pϕ − µR , K 2Pϕ2 where an insignificant constant term has been dropped. Let R be a region in the system’s phase space, defined by the condition 1 /2 |Pϕ − 1| < ∼µ . (4) We shall call it the resonant region, since the inequality (4) is equivalent to the inequality 1/ 2 |n − n′ | < , ∼ µ where n and n′ are the mean motions of the asteroid and of the planet respectively 2 . In the resonant region, the phase variables evolve at different rates. The variables Pϕ , Pg , Ph , g are the “slow” ones: dPϕ dPg dPh dg , , , ∼µ, dt dt dt dt 2 The given definition of the resonance region is a standard one for the investigations of the resonant phenomena in multifrequency Hamiltonian system obtained from integrable one by a small perturbation of order µ (Arnold et al. 2006). The asteroid entering into the Hill sphere results in the violation of the last condition (if the integrable Hamiltonian corresponds to asteroid’s motion in Keplerian orbit around Sun). For this reason we will consider only those trajectories which are definitely away from the planet’s Hill sphere. The perturbation theory should be developed in quite another way if one would like to take into account the dynamical effects close or inside the Hill’s sphere (e.g. Robutel and Pousse 2013). Quasi-satellite orbits. Perturbative treatment 7 while the variable ϕ is the “semi-fast”: dϕ ∼ µ1 / 2 , dt with h being the only “fast” variable in R : dh ∼1. dt As usual, the investigation of the secular effects in the motion of the asteroid begins with the averaging over the fast variable. The thus-averaged Hamiltonian eavr = 1 K 2π Z 2π 0 e dh K corresponds to the Hamiltonian system with two degrees of freedom, and it depends upon Ph as a parameter. Instead of the variables Pg , g , it is now convenient to introduce the variables x= p p 2(1 − Pg ) cos g , y = − 2(1 − Pg ) sin g . (5) They are defined on the disk D (Ph ) = {x2 + y 2 ≤ 2(1 − Ph )}. With an accuracy of O(µ1/2 ) , the relations between x, y, Ph and the osculating eccentricity and inclination take the form of e2 = 1 2 (x + y 2 )[4 − (x2 + y 2 )] , 4 cos i = (6) 2Ph . 2 − (x 2 + y 2 ) From Eq. (6), it follows that the centre of the disk D(Ph ) corresponds to the motion in a circular orbit. A similar relation for the argument of the pericenter can be written as ω=  2π − arccos p x ,y≥0   x2 + y 2     arccos p x , x2 + y 2 y<0 (x2 + y 2 6= 0). (7) For a given value of Ph , the eccentricity e and the inclination i of the asteroid do not exceed, at any moment of time, the values emax and q = 1 − Ph2 + O(µ1/2 ) imax = arccos Ph + O(µ1/2 ) , respectively. (8) (9) 8 V.V.Sidorenko et al. It would be worth dwelling upon the special case of Ph ≈ 1 . As follows from relations (8)-(9), this condition implies emax ≪ 1, imax ≪ 1 . Therefore this is the situation where the asteroid moves in the orbit of both a small eccentricity and a small inclination (Section 4). To simplify application of the perturbation technique to the analysis of the long-term evolution of motion, it will be instrumental to introduce the auxiliary quantity q σ = 1 − Ph2 . Under the condition Ph ≈ 1 , we clearly have σ≈ p e2 + i2 ≪ 1 . (10) As it follows from (10), the quantity σ characterizes in a certain sense how close the orbit of the asteroid is to the orbit of the planet. At σ → 0 the orbit of the asteroid becomes more and more close to the orbit of the planet. 2.2 The “slow-fast” system Our next step is standard for the analysis of resonant phenomena in the multifrequency systems (Arnold et al. 2006). Following the prescription of this theory, we undertake the scale transformation τ = εt, Φ = (1 − Pϕ )/ε ≈ (1 − a)/2ε , where ε = µ1/2 . Without loss of accuracy, it is possible to rewrite the averaged equations of motion in R as follows: dϕ = 3Φ, dτ dx ∂W =ε , dτ ∂y dΦ ∂W =− , dτ ∂ϕ (11) dy ∂W = −ε . dτ ∂x Here W (ϕ, x, y, Ph ) = 1 = 2π Z 2π 0 R[1, Pg (x, y ), Ph , ϕ − g (x, y ) − h, g (x, y ), h] dh. The expressions for Pg (x, y ) and g (x, y ) can be obtained easily from Eq. (5). Following (Schubart 1964) we apply the numerical integration to compute the averaged disturbing function W (ϕ, x, y, Ph ) .3 The dynamical system (11) will be called below the “slow-fast” system (or SF-system). Thus we would like to emphasise the existence of a timescale separation: the equations describing the behaviour of the variables ϕ, Φ make a “fast” 3 Some authors use for ing (Namouni et al. 1999). W (ϕ, x, y, Ph ) the term “ponderomotive potential” follow- Quasi-satellite orbits. Perturbative treatment 9 dΦ subsystem ( dϕ dτ , dτ ∼ 1 in general), while the “slow” subsystem consists of the dy equations for the variables x, y ( dx dτ , dτ ∼ ε ). Remark. Previously, the variable ϕ and its conjugate momentum were classified as “semi-fast” (Sec. 2.1). Averaging allows us to “forget” about the processes on the orbital motion time scale. So in the subsequent analysis of the evolutionary equations (11) we shall regard ϕ as a “fast” variable without risk of a confusion. The differential form Ψ = ε−1 dy ∧ dx + dΦ ∧ dϕ defines a symplectic structure in the phase space of the SF-system. The Hamiltonian of this system is 3Φ2 + W (ϕ, x, y, Ph ). 2 Finally, it is worth mentioning that differentiation of the averaged disturbing function W (ϕ, x, y, Ph ) with respect to Ph results in the following equation describing the evolution of the longitude of the ascending node: Ξ= dΩ ∂W = −ε . dτ ∂Ph 2.3 Properties of the fast subsystem Some of the results discussed in this section are not new - they were obtained in (Namouni et al. 1999). We rederive them here not only for the purposes of self-containedness: to apply Wisdom’s approach in the following sections we need to know more about the properties of the fast subsystem than we were able to learn from the preceding publications. At ε = 0 , the fast dynamics is governed by the one-degree-of-freedom Hamil- tonian system dϕ = 3Φ, dτ dΦ ∂W =− , dτ ∂ϕ (12) which depends on x, y, Ph as parameters. Let ϕ(τ, x, y, Ph , ξ ) , Φ(τ, x, y, Ph , ξ ) (13) denote a solution of equations (12), laying on the level set Ξ = ξ : Ξ (ϕ(τ, x, y, Ph , ξ ), Φ(τ, x, y, Ph , ξ ), x, y, Ph ) = ξ. In general, the angle ϕ oscillates or rotates: ϕ(τ + T, x, y, Ph , ξ ) = ϕ(τ, x, y, Ph , ξ ) mod 2π , T (x, y, Ph , ξ ) being the fast motion period. The solution (13) is rotational in the case of max W (ϕ, x, y, Ph ) < ξ ϕ (14) 10 V.V.Sidorenko et al. The rotation of the resonant phase ϕ corresponds to the motion of the asteroid in the non-resonant orbit. Following (Namouni et al. 1999) we shall refer to it as the passing orbit or, briefly, the P-orbit. The set of points (x, y ) satisfying inequality (14) for given ξ, Ph will be denoted with DP (ξ, Ph ) . Interpretation of oscillatory solutions requires caution since on the same level set Ξ = ξ different types of such solutions may exist. Emergence of such variety depends upon the number and location of the local maxima and minima of the function W (ϕ, x, y, Ph ) , for given values of the “parameters” x, y, Ph . Fig. 2 Behaviour of the averaged disturbing function W with fixed x, y, Ph . The following values are taken: Ph = 0.96825 ( σ = 0.25 ) in all cases, x = 0.18509, y = −0.10686 ( e = 0.21250, ω = 30◦ , i = 7.74554◦ ) in the upper panel, x = 0.04542, y = −0.20884 ( e = 0.21250, ω = 77.73110◦ , i = 7.74554◦ ) in the middle panel, x = 0.05415, y = −0.03127 ( e = 0.06250, ω = 30◦ , i = 14.03624◦ ) in the lower panel To illustrate the situation, we present, in Fig. 2, several examples demonstrating various types of behaviour of the function W . Similar to (Garfinkel 1977, Fig.2) and (Namouni et al. 1999, Fig.1) the abbreviations near the level lines characterize the secular evolution of the corresponding motion of the asteroid on the intermediate time scale (i.e., on the interval of order 1/µ1/2 in the initial units Quasi-satellite orbits. Perturbative treatment 11 of time): QS - quasi-satellite orbit, HS - horseshoe orbit, P - passing orbit, T tadpole orbit, QS+HS - a certain compound orbit (the existence of a variety of compound orbits at 1:1 mean motion resonance was revealed by Namouni et al. (1999), more details can be found in (Namouni 1999; Christou 2000); first example of QS+HS orbit appeared in (Wiegert et al. 1998)). All of the above mentioned orbits (except for the P-orbit only) are described by oscillatory solutions. The middle panel in Fig. 2 corresponds to the asteroid moving in the orbit which crosses the planetary orbit (the function W thus being unlimited). In the limit µ = 0 the orbits crossing at the exact mean motion resonance ( n = n′ ) takes place when asteroid’s osculating elements satisfy the condition cos ω = e. (15) On the plane (ω, e) the crossing condition (15) defines a curve, separating asteroid’s orbits linked and unlinked with the planet’s orbit (Fig. 3). The motion in the crossing orbits without a collision is possible due to the resonance. Fig. 3 The partition of the set of asteroid’s orbits ( n = n′ ) provided by the crossing condition (15): linked orbits - one node inside the planet’s orbit, the other is outside; unlinked orbits - both nodes inside the planet’s orbit. The blue line is a crossing curve defined by (15). The diagram is drown in two ways. In what follows we find it is more clear to present the phase portraits in the plane (e, ω) . Circular diagrams are convenient when we need to show regions satisfying a certain condition. The lower panel in Fig. 2 demonstrates that at some values of x, y, Ph motion in a QS-orbit is impossible. The relation between regions where such regimes are possible and the regions where such regimes are impossible is illustrated by 12 V.V.Sidorenko et al. Figs 4,5. Specifically, Fig. 4 provides the examples of the set BQS (Ph ) ⊂ D(Ph ) consisting of the elements (x, y ) for which, at a given value of Ph , the QS-regime is possible. Similar diagrams in terms of other variables can be found in (Christou 2000, Fig. 2) and (Mikkola et al. 2006, Fig. 6). Below we suppose to classify the qualitative properties of the asteroid’s motion depending on the value ξ of the Hamiltonian Ξ . We will use the designation DQS (ξ, Ph ) to identify the subset of BQS (Ph ) with elements for which there is a solution (11), corresponding to the motion in QS-orbit. It is easy to show that BQS (Ph ) = ∪ξ≥ξmin (Ph ) DQS (ξ, Ph ), where ξmin (Ph ) denotes the minimal value of ξ , for which the motion in QS-orbit is possible for a given Ph . We found numerically that DQS (ξ, Ph ) → ∂D(Ph ) as ξ ց ξmin (Ph ) . Fig. 5 demonstrates typical structural changes of DQS (ξ, Ph ) , which take place as ξ varies. If ξ is slightly larger than ξmin (Ph ) , the set DQS (ξ, Ph ) has a ringlike shape (Fig. 5,b). At ξ = ξh (Ph ) , a bifurcation takes place: holes emerge in the ring. As ξ increases, the number of holes decreases from four (Fig. 5,c) to two (Fig. 5,d). Simultaneously, the size of the holes increases, and at ξ = ξb (Ph ) their borders (the “internal” borders of DQS (ξ, Ph ) ) approach the outer border of DQS (ξ, Ph ) . The shapes of DQS (ξ, Ph ) for ξ > ξb are shown in Fig. 5,e and Fig. 5,f. The difference of the last Figures stems from the existence of the bifurcation at ξ = ξs (Ph ) , when the triangular regions in the central part of D(Ph ) separate from the peripheral part of DQS (ξ, Ph ) . Fig. 4 The diagram demonstrates at what values of the slow variables the QS-motion is possible for a given value of Ph = 0.96825 ( σ = emax = 0.25 ). The set BQS (Ph ) is painted in dark yellow. Fig. 6 and Fig. 7 provide more examples of the behaviour of the function W . Fig. 6 illustrates the formation of a singularity of the type W ∼ 1/σ at ϕ = 0 , Quasi-satellite orbits. Perturbative treatment 13 Fig. 5 A possible structure of the set DQS (Ph , ξ) for Ph = 0.96825 ( ξmin = 1.7735 , ξh = 2.3849 , ξb = 2.5470 , ξs = 4.0606 ): a - ξ < ξmin ( DQS (ξ, Ph ) = ⊘ ), b - ξ = 2.35 ∈ (ξmin , ξh ) , c - ξ = 2.40 ∈ (ξh , ξb ) , d - ξ = 2.45 ∈ (ξh , ξb ) , e - ξ = 3.50 ∈ (ξb , ξs ) , f - ξ = 4.50 > ξs Fig. 6 Behaviour of the averaged disturbing function W at fixed x, y, Ph . Here Ph = 0.99499 ( σ = 0.1 ) in both cases, x = 0.03898, y = −0.02251 ( e = 0.045, ω = 30◦ , i = 5.12871◦ ) in the left panel, x = 0.07368, y = −0.04254 ( e = 0.085, ω = 30◦ , i = 3.03062◦ ) in right panel as σ → 0 (or Ph → 1 ). In Section 4, this property will be used to analyze the asteroid motion in a QS-orbit with a small eccentricity and inclination. The right panel in Figure 7 demonstrates that at some x, y, Ph we can encounter three different periodic solutions on the same level of the Hamiltonian of the “fast” subsystem. 14 V.V.Sidorenko et al. Fig. 7 Behaviour of the averaged disturbing function W at fixed x, y, Ph . Here Ph = 0.86603 ( σ = 0.5 ) in both cases, x = 0.06500, y = −0.03753 ( e = 0.075, ω = 30◦ , i = 29.71851◦ ) in the left panel, x = 0.37711, y = −0.21772 ( e = 0.425, ω = 30◦ , i = 16.91652◦ ) in right panel. In this paper, the analysis is limited to the case where the behaviour of the function W is described by graphs similar to those presented in Fig. 2 and Fig. 6 (i.e. where a QS-orbit shares the same energy level with a HS-orbit). The examination of the function W led us to conclude that this kind of behaviour is taking place for all (x, y ) ∈ D(Ph ) at Ph > Ph∗ ≈ 0.95924 (or σ < σ ∗ ≈ 0.28258 ). As it follows from relations (8) and (9), this corresponds to the motion of the asteroid in orbits satisfying, for all moments of time, the inequalities e(t) < e∗max ≈ 0.28258, i(t) < i∗max ≈ 16.414◦ . So the rest of the paper is devoted to an exploration of the properties of QS-orbits with small and medium values of inclinations and eccentricities. The general situation can be treated in the same way, but it will be a rather cumbersome investigation. 2.4 Averaging over the resonant phase variations. Evolutionary equations To study the long-term behaviour of slow variables, we use evolutionary equations obtained by averaging of the right hand sides of equations (11) over the solutions of the “fast” subsystem: dx =ε dτ  ∂W ∂y  , D ∂W dy = −ε dτ ∂x where  ∂W ∂ζ  = 1 T (x, y, Ph , ξ ) Z T (x,y,Ph ,ξ ) 0 E , ∂W (ϕ(τ, x, y, Ph , ξ ), x, y, Ph )dτ, ∂ζ ζ = x, y. (16) (17) Quasi-satellite orbits. Perturbative treatment 15 In the region DQS (ξ, Ph ) , the averaging procedure provides us with two vector fields, depending on what periodic solution of the “fast” subsystem was chosen in Eq. (17). One of them describes the evolution of slow variables in the case of asteroid motion in a QS-orbit, the other characterizes the evolution of an HS-orbit. Applying the averaging procedure (17) we should take care of the situation when the solution (13) is non-periodic and corresponds to the separatrix on the “fast” subsystem’s phase portrait. Taking this in mind we define in the disk D(Ph ) the so-called uncertainty curve Γ (ξ, Ph ) (Wisdom 1985; Neishtadt 1987; Neishtadt and Sidorenko 2004; Sidorenko 2006): Γ (ξ, Ph ) = {(x, y ) ∈ D(Ph ), ∃ϕ∗ (x, y, Ph ) which is a point of a local or global maximum of the function W (ϕ, x, y, Ph ) ( x, y, Ph being treated as parameters) satisfying the condition W (ϕ∗ , x, y, Ph ) = ξ} As it is evident, for a given value of Ph the “separatrix” solution exists on the level Ξ = ξ only when (x, y ) ∈ Γ (ξ, Ph ) . For ξ ≥ ξh (Ph ) the uncertainty curve consists typically of one or two components: ΓP (ξ, Ph ) = ∂DP (ξ, Ph ) and ΓQS (ξ, Ph ) ⊂ ∂DQS (ξ, Ph ) (Fig. 8). For ξ ∈ [ξmin (Ph ), ξh (Ph )) it does not exist. Fig. 8 The location of the uncertainty curve components ΓQS and ΓP for Ph = 0.96825, ξ = 3.5 . Green represents the areas of the asteroid motion in the passing orbits at these values of Ph and ξ . When the projection of the system (11) phase point on the plane (x, y ) approaches the uncertainty curve Γ (ξ, Ph ) the qualitative changes in the behavior of the resonant phase ϕ take place (corresponding, for example, to a transition from a QS-orbit to an HS-orbit or back). To prolongate the solutions of the averaged equations (16) across the uncertainty curve, one can follow a rather straightforward strategy based on the matching of solution with the same limit values at Γ (ξ, Ph ) . In some cases, this matching is not unique. For example, in Fig. 9,a we present a situation at the border of DQS (ξ, Ph ) when one option for continuation results in exiting this region, with a further motion in HS+QS-orbit, while the second option can be described as reflection from the border, with a transition from QS-motion 16 V.V.Sidorenko et al. to HS-motion. Since both options are possible in the system (11), the evolution of the motion near the uncertainty curve in this example has a probabilistic nature4 (Arnold et al. 2006). Similar example is provided by Fig. 9,d. Fig. 9,b and Fig. 9,c present situations where the solution leaves the neighbourhood of ΓQS (ξ, Ph ) in a unique way. Fig. 9 Matching of solutions to the evolutionary equations (16) at the uncertainty curve Γ (ξ, Ph ) 2.5 Adiabatic approximation The evolutionary equations (16) provide so called “adiabatic approximation” of the system’s long-term dynamics (Wisdom 1985; Yokoyama 1996). Indeed the first two equations in (11) correspond to 1DOF Hamiltonian system depending on slowly varying parameters x, y . Therefore this system approximately preserves the value of the adiabatic invariant (Henrard 4 If we consider the solution of the non-averaged system with the initial conditions x(0), y(0) exactly on Γ the dynamics is definitely deterministic and depends mainly on to what part of fast subsystem separatrix variables ϕ(0), Φ(0) belong. But far from the uncertainty curve (i.e., deep inside in the yellow zone in our phase portraits) the initial conditions corresponding to the crossing the above mention parts of the separatrix in the plane (ϕ, Φ) are mixed strongly. And a small uncertainty in the initial conditions does not allow us to predict uniquely the qualitative behavior of the system when the phase point leave the vicinity of Γ . Nevertheless, for the set of the possible initial conditions we can compare the measures of subsets resulting in the different qualitative behavior later on and then characterize the dynamics in the probabilistic way. Quasi-satellite orbits. Perturbative treatment 17 1982; Arnold et al. 2006) I (x, y, Ph , ξ ) = I Φ dϕ. (18) For averaged equations (16) the adiabatic invariant (18) is an exact integral. So the phase trajectories x(τ ), y (τ ) are laying on its level curves. The matching of the phase trajectories at Γ satisfies the conditions IHS +QS = IQS + IHS and 2IP = IHS +QS or 2IP = IHS for the components ΓQS and ΓP respectively (indices denote the orbital regimes used to compute the adiabatic invariant (18)). The validity of the adiabatic approximation is limited due the violation of the adiabatic invariance in the vicinity of the uncertainty curve Γ . The related phenomena are discussed in Sec. 3.2. 3 Investigation of secular effects on the basis of evolutionary equations 3.1 Comparative analysis for different level sets of the Hamiltonian Ξ A good way to understand the evolution of slow variables is to consider the phase portraits of system (16). As it follows from Sec. 2.5 the drawing of the phase portraits reduces to the drawing of the level curves of the adiabatic invariants IP,QS,... (x, y, Ph , ξ ) for given values Ph , ξ . Several examples of phase portraits are given in Fig. 10. For better visualization, these portraits present the behaviour of the averaged osculating elements ω and e given by formulae (6) and (7). Dependent on the values of ξ, Ph , these phase portraits differ in the number of equilibrium points and/or the separatrices placement. There exist some other important features which depend on ξ, Ph also. It is interesting to compare these phase portraits with the evolutionary diagrams obtained in (Namouni 1999) by means of direct integration of motion equations. The phase portrait in Fig. 10 resembles Fig. 12 from (Namouni 1999) and represents the situation when the motion in a QS-orbit is impossible ( ξ < ξmin (Ph ) ). All trajectories are related to HS-orbits with a circulating argument of the pericenter, ω . The relatively simple portrait in Fig. 10,b is typical for ξ ∈ (ξmin , ξh ) . It is the case when transitions between the motion in a QS-orbit and that in an HS-orbit are impossible. As a remarkable property of this case, we can mention the opposite directions of the pericenter circulation for these orbits. There is no similar diagram in (Namouni 1999), because the motion equations were integrated only with the initial conditions corresponding to P-orbits and HS-orbits. More complex behaviour of phase trajectories is presented on the phase portrait in Fig. 10,c (compare with Fig. 13 in (Namouni 1999)). One of the most interesting properties of asteroid motion is associated with the closed contours composed by the phase trajectories related to a QS-orbit and an HS-orbit. Such a contour 18 V.V.Sidorenko et al. corresponds to the alternation of the QS- and HS-regimes in orbital motion over a large enough time span. Nevertheless, at a certain moment, this alternation can be violated due to the escape from DQS (ξ, Ph ) . The opportunity of escape is provided by outgoing trajectory matched to the described contour (Fig. 9,a or d). Numerically, this phenomenon was established by Namouni (1999). The phase portrait in Fig. 10,d demonstrates the further complication of asteroid dynamics as ξ increases. The fission of DQS (ξ, Ph ) (the separation of triangular areas from peripheral part) at ξ = ξs is accompanied by the appearance of a stable equilibrium corresponding to motion in passing orbit with ”frozen” pericenter: ω = 90◦ or ω = 270◦ (Lidov-Kozai resonance). In (Namouni 1999) the passing orbits with ω librating around Lidov-Kozai resonance at ω = 90◦ or ω = 270◦ can be seen in Fig. 18. Fig. 10 Examples of the secular evolution at 1:1 mean motion resonance. Thick lines correspond to motions in QS-orbits. Two sets of phase trajectories in the domain DQS (ξ, Ph ) (painted yellow) demonstrate the dependence of secular effects on the character of motion (HS-orbit or QS-orbit) on intermediate time scale. In all cases Ph = 0.96825 ( emax = 0.25 ). The values of Hamiltonian of SF-system (11): a - ξ = 1.75 , b - ξ = 2.35 , c - ξ = 3.5 , d ξ = 4.5 In Fig. 11 we tried to summarize some information about the motion in QSorbits. The choice of parameters ( σ and ξ = ξ · σ ) is justified by our intention to simplify the structure of the presented diagram. As one can see, there are practically straight-line borders between regions with different properties of QSorbits, and the following limits exist: ξ min = lim σ · ξmin (Ph ) = 0.68644..., σ→0 ξ h = lim σ · ξh (Ph ) = 0.8346..., σ→0 ξ b = lim σ · ξb (Ph ) = 0.886... σ→0 Quasi-satellite orbits. Perturbative treatment 19 The regularity of the system properties with respect to proper scaling of parameters stimulates an application of the appropriate perturbation technique for analysis of asteroid dynamics in case σ ≪ 1 (Section 4). Fig. 11 Dependence of QS-orbit properties on parameters σ, ξ . Area A: any motion in QSorbit preserves its character forever. Area B: due to the holes in DQS (ξ, Ph ) (e.g., Fig.3,c and Fig.3,d) for some initial values the alternating escapes from QS-orbits and returns back become possible. Area C: the configuration of DQS (ξ, Ph ) (e.g., Fig.3,e and Fig.3,f) allows permanent motion in QS-orbit only at small inclination 3.2 Chaotization Matching of phase trajectories on the uncertainty curve Γ (ξ, Ph ) provides “zeroorder” theory of slow variables evolution in the case of qualitative transformation of fast variables behaviour. Actually in the neighbourhood of Γ (ξ, Ph ) the adiabatic invariance is violated. As a consequence the projection of a phase point of the system (11) onto the plane x, y jumps in quasi-random way from the incoming trajectory of the averaged system (16) to some outgoing trajectory located at the distance of the order ε| ln ε| from the outgoing trajectory obtained by a formal matching. The accurate estimations of the quasi-random jumps were given in (Timofeev 1978; Cary et al. 1986; Neishtadt 1986, 1987). This mechanism of chaotization (giving rise to the so-called adiabatic chaos in the asteroid’s dynamics) is typical for mean-motion resonances (Wisdom 1985; Neishtadt 1987; Neishtadt and Sidorenko 2004; Sidlichovsky 2005; Sidorenko 2006; Batygin and Morbidelli 2013). It is likely the responsible for the increment of the eccentricity observed by Namouni (1999) in the series of the transitions QS → HS → . . . . 20 V.V.Sidorenko et al. 4 Dynamics of the asteroid moving in an orbit of a small eccentricity and inclination In the case of the asteroid moving in an orbit with a small eccentricity and inclination, the resonance condition 1/2 |n − n′ | < ∼ µ implies Ph ≈ 1 or σ = p 1 − Ph2 ≪ 1 (Section 2.1). Using different simplifying assumption (Namouni 1999) derived analytical expressions characterizing the secular evolution evolution of QS-orbits, HS-orbits and P-orbits. Our goal is to provide a qualitative description of long-term dynamics for any resonant orbit with a small eccentricity and inclination satisfying the only restriction imposed at the beginning of the paper: do not approach the Hill sphere of the planet. In particular it implies µ1/3 ≪ σ ≪ 1 . To achieve our goal we must consider the motions which do not satisfy the simplifying assumptions applied in Namouni (1999). 4.1 The leading term in the expression for the averaged disturbing function W Graphs presented in Fig. 6 demonstrate two important properties of the averaged disturbing function at σ → 0 ( Ph → 1 ): (i) the increase of its values for a resonant phase ϕ close to zero, and (ii) the decrease of the profile “thickness” (i.e., the decrease of the interval of resonant phase values at which W ∼ 1/σ ). Analytically, it implies the following structure for the function W : W (ϕ, x, y, σ ) = 1 σ W   ϕ x y +{Residual term which is in general of the order of 1}, , , σ σ σ where ϕ ∈ [−π, π ] and (19) 1 W (ϕ, x, y ) = 2π ∆(ϕ, x, y, h) = lim σ→0 1 σ Z π −π dh , ∆ |r(ϕ, x, y, h) − r(λ′ )| ϕ=σϕ, x=σx, y =σy = 2 ϕ − 4ϕ(x sin h − y cos h) + 3(x sin h − y cos h)2 + 2 2 2 2 (x + y ) + [1 − (x + y )] sin2 h. Now we shall address the approximate evolutionary equations obtained by approximating W with its leading term only. To that end, it will be useful to discuss the main properties of the function W (ϕ, x, y ) . The function W is defined on the set R1 × D2 \ D0 , where D 2 = {(x, y ) ∈ R2 | x2 + y 2 ≤ 1}, D0 = {ϕ, x, y| ∃h : ∆(ϕ, x, y, h) = 0}. The following symmetries are easy to verify: W (ϕ, x, y ) = W (ϕ, −x, y ) = W (−ϕ, x, −y ) = W (−ϕ, x, y ). Quasi-satellite orbits. Perturbative treatment 21 It is noteworthy that, for a given values x, y , the function W is even with respect to ϕ , while the function W does not, in general, possesses this property. In the studies of QS-orbits, the excessive symmetry of the dynamical model based on truncated expression for the disturbing function were mentioned previously by Mikkola et al. (2006). This symmetry should not be surprising: in a non-explicit way, here we apply the Hill’s approximation of the three body problem; the excessive symmetry of Hill’s approximation (in comparison with the original problem) is well known (Henon 1997). Fig. 12 The behaviour of the function W at some fixed values of x, y : x = 0.7, y = 0.0 (red line); x = 0.85, y = 0.0 (green line); x = 1.0, y = 0.0 (blue line) The function W (ϕ, x, y ) and its derivatives ∂W , ∂y ∂W , ∂x ∂W ∂ϕ can be expressed in terms of the elliptic integrals of the first and second kind (Appendix B). One can use these expressions to accelerate numerical computations at the stage of averaging over the resonant phase oscillations/rotations. In the case of 2 2 ϕ ∼ ( x + y ) 1 /2 ≪ 1 , the following approximate formula can be derived: W (ϕ, x, y ) ≈ 1 ln 2π h 256 2 2 2 ϕ + (x + 4y ) i2 2 2 . (20) − 16ϕ y Expression (20) allows one to understand better the behaviour of the function W (ϕ, x, y ) for small enough ϕ, x, y , i.e., when this function becomes singular. 22 V.V.Sidorenko et al. Fig. 13 Examples of the set D QS (ξ) structure: a - ξ < ξ min ( D QS (ξ) = ⊘ ), b - ξ = ∗ ∗ 0.82 ∈ (ξ min , ξ h ) , c - ξ = 0.84 ∈ (ξ h , ξ b ) , d - ξ = 1.2 > ξ b 4.2 Derivation of the evolutionary equations in the case of the leading term approximation for W If the leading term approximation is used for W , the equations for fast subsystem (12) can be rewritten as dϕ = 3Φ, dτ ∂W dΦ =− , dτ ∂ϕ (21) where τ = σ −3/2 τ, Φ = σ 1/2 Φ, ϕ = σ −1 ϕ. We are interested in oscillatory solutions to Eqs. (21), as this naturally provides the approximation for solutions to Eqs. (12), corresponding to QS-regimes of the orbital motion. Let the pair ϕ(τ , x, y, ξ ) , Φ(τ , x, y, ξ ) , denote a solution to Eqs. (21), satisfying the equality 3 2 Φ (τ , x, y, ξ ) + W (ϕ(τ , x, y, ξ ), x, y ) = ξ 2 (22) Quasi-satellite orbits. Perturbative treatment 23 at any moment τ ∈ R1 ( x, y being considered as fixed parameters). This solution exists only in the case of (x, y ) ∈ DQS (ξ ) , where   ∂2W (0, x, y ) > 0 DQS (ξ ) = (x, y ) ∈ D2 : W (0, x, y ) < ξ, 2 ∂ϕ Fig. 13 illustrates the dependence of DQS (ξ ) on ξ . The changes in the structure of this set resemble the changes of the set DQS (ξ, Ph ) when ξ varies (Fig. 5), although in general the leading term approximation of the disturbing function results in the loss of certain fine details (some of them we discuss below). For ξ ∈ (ξ min , ξ h ) , the set DQS (ξ ) has a ring-like structure ( ξ min , ξ h have been defined at the end of Section 3.1). Using the expression for W in terms of elliptic integrals, one obtains: 1 ξmin = π r ! 3 4 K ∗ ∗ At ξ = ξ h , holes emerge in DQS (ξ ) ; at ξ = ξ b ( ξ b = 0.89542... ) the borders of these holes approach the “external” border of DQS (ξ ) . The difference between ∗ the values of ξ b and ξ b is probably the consequence of the excessive symmetry of Hill’s approximation mentioned above (one more consequence being the existence ∗ of only two holes in DQS (ξ ) for all ξ ∈ (ξ h , ξ b ) ). The solution (22) is used to obtain an approximate expression for the righthand sides of equations (16) describing the evolution of the slow variables on the level set Ξ = ξ/σ for the case of the asteroid moving in a QS-orbit:  ∂W ∂ζ  ≈ 1 T (ξ/σ, Ph (σ )) ≈ Z T (ξ/σ,Ph (σ )) σ 2 ∂ζ 0 1 σ 2 T (ξ ) 1 ∂W Z T (ξ ) 0 ∂W ∂ζ    ϕ σ ϕ (τ , ξ ) , ξ τ, Ph (σ ), σ x y , σ σ   x y , , σ σ  dτ ≈ (23) dτ . Here T (ξ ) denotes the period of oscillatory solutions (22), while ζ = x, y and ζ = x, y . In the case of the asteroid moving in HS- or P-orbits, construction of similar approximate expression is a more delicate procedure since there is no suitable periodic solution to the auxiliary system (21). However one can use non-periodic solutions to this system. Specifically, let us consider the solution to Eqs (12), corresponding to motion in an HS-orbit with Ξ = ξ/σ ≫ 1 . A crude (but sufficient for our purposes) estimate of its period is T ≈π r 8σ . 3ξ (24) Due to the afore-mentioned symmetry of the leading term in the expression for the disturbing function, averaging can be performed over the reduced time interval:  ∂W ∂ζ  ≈ 4 T Z T /4 0 1 ∂W σ 2 ∂ζ   ϕ σ ξ τ, Ph (σ ), σ  x y , , σ σ  dτ (25) 24 V.V.Sidorenko et al. Then we can replace ϕ(τ, ξ/σ, Ph (σ )/σ ) in the expression (25) with the function ϕ(τ , ξ ) , which provides (in combination with Φ(τ , ξ ) = dϕ/dτ ) a non-periodic solution to (21), satisfying the conditions:  W ϕ(0, ξ ), This results in  ∂W ∂ζ  ≈ 4 T x y , σ σ Z  = ξ, ϕ(τ , ξ ) → +∞ at τ → +∞ . σ −3/2 T /4 1 ∂W σ 1/2 ∂ζ 0  ϕ(τ, ξ ), x y , σ σ  dτ . (26) To finish up, we replace the period T in Eq. (26) with its approximate expression (24) and expand the integration over the positive semiaxis (since σ −3/2 T (ξ/σ ) approaches infinity as σ → 0 ):  ∂W ∂ζ  Z 1 p ≈ 23ξ πσ +∞ 0 ∂W ∂ζ  ϕ(τ , ξ ) , x y , σ σ  dτ . (27) The convergence of the integral in Eq. (27) can be proven rigorously: it follows from the rapid decrease of the function W at ϕ → ∞ . Physically, the approximation (27) means that in the case of motion in an HS-orbit the evolution of slow variables ( e, ω, i, Ω ) is substantial only during close approaches of the asteroid to the planet. In a similar way, the “averaging” can be done for the motion in a P-orbit. Formally, the ensuing expression coincides with (27), but the function ϕ(τ , ξ ) is rendered in this case by a solution to Eq. (21), satisfying the conditions: ϕ(0, ξ ) = 0, ϕ(τ , ξ ) → +∞ at τ → +∞ . 4.3 Dynamics of the asteroid. Estimations based on the approximate evolutionary equations Insertion of expressions (23) and (27) into Eqs. (16) furnishes approximate equations characterising the secular evolution of asteroid motion in a near-circular low-inclination orbit in the 1:1 resonance. In Fig. 14, we present examples of the phase portraits corresponding to the approximate evolutionary equations. The behaviour of the phase trajectories has a remarkable similarity: it does not depend on σ , if we use emax as a unit for length along the vertical axis. Fig. 14,a is practically identical to Fig. 10,b – it illustrates that the approximate equations are best of all suited for the description of the motion without qualitative changes in the behavior of the resonant phase ϕ (i.e., asteroid moves permanently in a QS- or HS-orbit). Fig. 14,b resembles Fig. 10,c except for the location of the regions corresponding to P-orbits (these regions are painted in green). The difference emerges due to the absence of compound QS+HS-regimes in the model based on the leadingterm approximation of W (this is another consequence of its excessive symmetry). So the model discussed above should be applied with caution for the analysis of motions with transitions between different types of orbits. Quasi-satellite orbits. Perturbative treatment 25 Fig. 14 Examples of phase portraits obtained in the leading term approximation for W : a ξ/σ = 0.82 < ξ h , b - ξ/σ = 1.2 > ξ b The approximate expressions (23) and (27) for the averaged derivatives of W allow one to estimate the rates of the evolution of the orbital elements: µ dΩ dω µ de di , ∼ 2, , ∼ 3 dt dt σ dt dt σ (28) for the case of motion in a QS-orbit, and de di µ dΩ dω µ , ∼ , , ∼ 2 , dt dt σ dt dt σ (29) for other types of orbital motion. As it follows from Eqs. (28),(29), evolution is faster for motion in a QS-orbit. This is not surprising since for such an orbit the mean distance from the asteroid to the planet is smaller than for other types of orbits and, consequently, the perturbations due to the attraction of the asteroid by the planet are more substantial. For motion with transitions between different types of orbital behaviour, formulae (28) and (29) yield the following estimates TQS ∼ σ3 , µ THS,P,... ∼ σ2 , µ (30) where TQS and THS,P,... denote the duration of motion in a (temporary) QSorbit and in orbits of other types, respectively. An evident consequence of (30) is: TQS ≪ THS,P,... . 26 V.V.Sidorenko et al. As it was mention at the beginning of this Section under some additional assumptions Namouni (1999) derived analytical expressions describing the secular evolution of the orbits under investigation. In particular, his formulas for QS-orbits are limited to the case when the resonant phase ϕ does not oscillate - i.e., the “fast” subsystem is in “quasi-steady” equilibrium close to local minimum of the averaged disturbing function W (ϕ, x, y, Ph ) (in our phase portraits the motions in these orbits correspond to trajectories at the lower boundary of the yellow areas). The HS-orbits in (Namouni 1999, Sec. 4.1) approach the close vicinity of the planet’s Hill sphere what is beyond the scope of our analysis. The consideration of P-orbits in (Namouni 1999, Sec. 4.2) will be commented in the Section 4.4. 4.4 Some remarks about the properties of the motion in the passing orbits with small eccentricity and inclination The relation of the Fig. 3 and Fig. 10,d leads to a conclusion that there are two topologically different types of the passing orbits. The passing orbits with ω librating around 0◦ or 180◦ are linked with the orbit of the planet. The libration of ω around 90◦ or 270◦ takes place in the case of the motion in unlinked passing orbits. To study the motion in linked passing orbits Namouni (1999) applied the averaging over the resonant phase ϕ neglecting the variations in its rate. Formally it is identical to the averaging over motion in the close nonresonant orbits undertaken by Lidov and Ziglin (1974). As a result, the expressions for double-averaged disturbing function in these papers differ by designations only. But then Lidov and Ziglin (1974) only mentioned the √ existence of the “centre” type equilibrium at e = σ/ 2 , ω = 0◦ or ω = 180◦ , while Namouni (1999) derived the analytical solution to the system of the evolutionary equations. For small but non-zero values of σ the regions of the unlinked orbits can be described as the tiny triangles encompassing the lines ω = 90◦ and ω = 270◦ . As a consequence of the results presented in (Lidov and Ziglin 1974) one can obtain for unlinked p passing orbits the existence of the “centre” type equilibrium at e ≈ σ 23 , ω = 90◦ or ω = 270◦ , what is in a good agreement with the numerical investigations. 5 Example: the future escape of asteroid 2004GU9 from the QS-orbit Now we apply the evolutionary equations (16) for the analysis of secular effects in orbital motion of the near-Earth asteroid 2004GU9. The restricted circular threebody problem is insufficient for accurate investigation of the actual asteroid dynamics. So we are interested more in understanding of the time scales involved and of some other quantitative characteristics of the phenomena discussed in previous Sections. Currently, asteroid 2004GU9 is moving along a QS-orbit with osculating elements presented in Table 1 (Mikkola et al. 2006; Wajer 2010). We chose it among the other quasi-satellites of the Earth due to the absence of this object’s close Quasi-satellite orbits. Perturbative treatment 27 encounters with Venus and Mars – this restriction justifies, to some extent, the investigation of the secular effects under the scope of RC3BP. Fig. 15 demonstrates the behaviour of the resonant phase ϕ according to the results of direct numerical integration of the equations of motion corresponding to RC3BP with the initial values provided by the elements in Table 1 and the mass parameter µ = 3.04 · 10−6 (we added the mass of the Moon to the mass of the Earth). The motion in a QS-orbit will last for approximately 500 years, with a subsequent transition to an HS-orbit. In (Wajer 2010), similar conclusion has been achieved via taking into account the gravitational pull of all other Solar system planets, as well as the Moon and Pluto. Table 1 Osculating orbital elements of asteroid 2004GU9. Epoch: March 14, 2012 (JD2456000.5) Element Value a(AU ) 1.001056350821795 e .1362904920360489 i (◦ ) 13.64944749947083 Ω (◦ ) 38.74489028357296 ω (◦ ) 280.6255989836612 M (◦ ) 217.2153150601352 Note. The values of the orbital elements were borrowed from the JPL Small-Body Database Fig. 15 The behaviour of the resonant phase ϕ of 2004GU9 The graphs in Fig. 16 allow us to compare the long term evolution of the osculating elements, according to Eqs. (16), with the results of a direct numerical 28 V.V.Sidorenko et al. integration. In Fig. 17 we present the projection of the asteroid phase trajectory onto the plane ω, e (red line). Black lines correspond to the appropriate phase trajectories of the system (16) when the matching at the uncertainty curve is accomplished according to the procedure described in Sec.2.4. Fig. 16 and Fig. 17 demonstrate that our approach based on double averaging of the equations of motion provides a really accurate description of the secular evolution. Fig. 16 The evolution of the orbital elements of 2004GU9. Blue lines correspond to the result of integration of the equation of motion; red lines characterise the secular behaviour of the orbital elements according to the evolutionary equations (16) Conclusion In this paper we have considered the three body system “Sun-planet-asteroid” in the 1:1 mean-motion resonance. A special attention was given to the motion of the asteroid in a quasi-satellite orbit, when the asteroid is located (permanently or for a long enough time) in the vicinity of the planet, though being out of the planet’s Hill sphere. As usual, in the mean-motion resonance three dynamical time scales can be distinguished. The “fast” process corresponds to the planet and asteroid orbital motions. The “semi-fast” process is a variation of the resonant phase (which, in a certain sense, describes the relative position of the planet and asteroid in their orbital motion). Finally, the “slow” process is the secular evolution of the orbit’s shape (the eccentricity) and orientation (the longitude of the ascending node, the inclination, and the argument of the pericenter). To study the “slow” process, we constructed the evolutionary equations by means of numerical averaging over the “fast” and “semi-fast” motion. As a spe- Quasi-satellite orbits. Perturbative treatment 29 Fig. 17 Phase portrait of the system (16) at ξ = 2.87506, Ph = 0.96269 ( σ = 0.27061 ). The red line characterises the evolution of ω and e in the solution of the equations of motion with the initial conditions corresponding to the current values of the orbital elements of 2004GU9 (marked with the red dot) cific feature of these evolutionary equations, we should point out the ambiguity of their right-hand sides for some values of the “slow” variables. The ambiguity appears because the averaging can be performed over the “semi-fast” processes with qualitatively different properties (in other words, it can be done over a QS-orbit, or an HS-orbit, with the same values of the Hamiltonian Ξ ). Consideration of this ambiguity provided us with an opportunity to predict whether the motions in QSor HS-orbits are permanent or not. For non-permanent motions in QS-orbits, the conditions of capture into this regime and escape from it have been established. For clarity, our analysis was restricted to the case of the asteroid motion in orbits with small and intermediate values of the inclinations and eccentricities, when maximum two regimes of the “semi-fast” evolution are possible. The general case can be studied similarly. The evolutionary equations were used to draw the phase portraits demonstrating in a clear way the secular effects in the asteroid motion (for example, the libration or oscillation of the pericenter). Some kind of scaling has been found in the asteroid dynamics at small inclinations and eccentricities. It has turned out that in this case many properties of the motion (for example, the duration of the QS-regime) can be established using the Hill approximation in the restricted circular three-body problem. The excessive symmetry of this approximation creates a difficulty for correct description of all possible modes of motion, but hopefully it can be overcomed after taking into account the high order terms of perturbation theory. To illustrate the typical rates of the orbital elements’s secular evolution, the dynamics of the near-Earth asteroid 2004GU9 was considered. This asteroid will keep its motion in a QS-orbit for the next several hundred years. For better prediction, our model should be improved by taking into account the influence of the 30 V.V.Sidorenko et al. other planets (Venus, Mars, Jupiter, etc), along with the appropriate modification of the averaging procedure. Acknowledgments The work was supported in part by the Russian Foundation for Basic Research (projects 13-01-00251 and NSh-2519.2012.1) and by the Presidium of the Russian Academy of Sciences under the scope of the Program 22 “Fundamental problems of Solar system investigations”. We are grateful to M.A.Vashkovyak and M.Efroimsky for reading the manuscript and useful discussions. We thank also anonymous referees for all their corrections and suggestions. Quasi-satellite orbits. Perturbative treatment 31 Appendix A. Calculation of the function W and its derivatives Below we describe some technical details of the computational procedure which has been applied to obtain numerically the values of the averaged disturbing function Z 1 W (ϕ, x, y, Ph ) = 2π and of its derivatives ∂W , ∂ϕ ∂W , ∂x 2π R(ϕ, x, y, Ph , h) dh (A.1) 0 ∂W , ∂y ∂W ∂Ph with no limitations on the asteroid’s eccentricity or inclination. 1.Change of the variable, over which the integration takes place. To start with, we replace in (A.1) the integration over the auxiliary variable h with integration over the asteroid’s eccentric anomaly E . This allows us to avoid the necessity to solve the Kepler equation in our calculation of the asteroid’s position. The relation between the variables h and E reads as h = ϕ − ω − (E − e sin E ). (A.2) As a result, the expressions for W and its derivative can be written as 1 W (ϕ, x, y, Ph ) = − 2π ∂W 1 =− ∂ϕ 2π Z 2π 0 ∂R ∂h dE, ∂ϕ ∂E 1 ∂W =− ∂x 2π ∂W 1 =− ∂y 2π Z 2π 0 Z 2π 0   Z 2π R 0 ∂h dE ∂E ∂W 1 =− ∂Ph 2π Z 2π 0 ∂R ∂h ∂2h +R ∂x ∂E ∂x∂E ∂R ∂h ∂2h +R ∂y ∂E ∂y∂E   (A.3) ∂R ∂h dE, ∂Ph ∂E dE, dE, where ∂h = −(1 − e sin E ), ∂E ∂2h ∂ 2 h ∂e ∂e = · = cos E · , ∂E∂x ∂E∂e ∂x ∂x ∂ 2 h ∂e ∂e ∂2h = · = cos E · . ∂E∂y ∂E∂e ∂y ∂y The “minus” sign before the integrals in formulae (A.3) appears due to our intention to have the upper limit larger than the lower one. 2. Derivatives of R with respect to ϕ, x, y, Ph . At this point, we need to introduce the uniformly rotating heliocentric reference frame Oξηζ with the axis Oξ being directed from the Sun to the planet and the axis Oζ being aligned with the normal to the plane of the primary’s orbital motion. In this reference frame the position vector of the planet is r′ = eξ , where eξ = (1, 0, 0)T is the unit vector 32 V.V.Sidorenko et al. corresponding to axis Oξ . Respectively the expression for disturbing function (2) is reduced to 1 − (r, eξ ). R= |r − eξ | Then we get: ∂R ∂R ∂r = , ∂ϕ ∂r ∂ϕ ∂R ∂R ∂r = , ∂x ∂r ∂x where ∂R ∂R ∂r = , ∂y ∂r ∂y ∂R ∂R ∂r = , ∂Ph ∂r ∂Ph r − eξ ∂R =− − eξ . ∂r |r − eξ |3 3. Position vector of the asteroid. Let the unit vector e∗ξ be directed to the pericenter of the osculating orbit, while one more unit vector e∗η is parallel to this orbit’s minor semi-axis (and directed in the same way as asteroid’s velocity vector at the pericenter). It is not difficult to write down the expressions for introduced unit vectors through their projections on the axes of the reference frame Oξηζ :   cos h cos ω − cos i sin h sin ω ∗  eξ = sin h cos ω + cos i cos h sin ω  sin i sin ω and   − cos h sin ω − cos i sin h cos ω ∗  eη = − sin h sin ω + cos i cos h cos ω  sin i cos ω (A.4) (A.5) After that the asteroid’s position vector r can be written as the linear combination r(i, e, ω, h, E ) = e∗ξ (i, ω, h)ξ ∗ (e, E ) + e∗η (i, ω, h)η ∗ (e, E ) (A.6) with the coefficients (here the orbit with a = 1 is considered) ξ ∗ (e, E ) = cos E − e, η ∗ (e, E ) = p 1 − e2 sin E. 4. Derivatives of the position vector r with respect to ϕ, x, y, Ph . Taking into account (A.2) and (A.6) we obtain the following relations: ∂r = ∂x  ∂r ∂h ∂r + ∂e ∂h ∂e  ∂e + ∂x  ∂r ∂h ∂r + ∂ω ∂h ∂ω  ∂ω ∂r ∂i + , ∂x ∂i ∂x ∂r = ∂y  ∂r ∂h ∂r + ∂e ∂h ∂e  ∂e + ∂y  ∂r ∂h ∂r + ∂ω ∂h ∂ω  ∂r ∂i ∂ω + , ∂y ∂i ∂y ∂r ∂h ∂r ∂r , = = ∂ϕ ∂h ∂ϕ ∂h ∂r ∂r ∂i = , ∂Ph ∂i ∂Ph where evidently ∂h = − 1, ∂ω ∂h = sin E. ∂e Quasi-satellite orbits. Perturbative treatment 33 5. Derivatives of the position vector r with respect to h and osculating variables i, e, ω . Direct differentiation of (A.6) results in where ∂e∗ξ ∗ ∂e∗η ∗ ∂r = ξ + η , ∂h ∂h ∂h ∂e∗ξ ∗ ∂e∗η ∗ ∂r = ξ + η , ∂i ∂i ∂i ∂e∗ξ ∗ ∂e∗η ∗ ∂r = ξ + η , ∂ω ∂ω ∂ω ∂ξ ∗ ∂r ∂η ∗ = e∗ξ + e∗η , ∂e ∂e ∂e ∂ξ ∗ = 1, ∂e (A.7) ∂η ∗ e sin E = −√ . ∂e 1 − e2 The derivatives of the unit vectors e∗ξ and e∗η in (A.7) can be easily evaluated from (A.4) and (A.5) respectively. 6. Derivatives of the osculating variables e, i, ω with respect to x, y, Ph . First we rewrite the expressions for e, cos i, sin i in a more compact way (compare with (6)): e= Here p s ( x2 + y 2 ) 2P h cos i = , s−2 2 sin i = p (A.8) , (s − 2)2 − 4Ph2 . s−2 s = 4 − ( x2 + y 2 ) . After the differentiation of the expressions (A.8) with respect to x, y we obtain: and x ∂e = ∂x 2 r ∂e y = ∂y 2 r s − x2 + y 2 s x2 + y 2 − r x2 y2 s ! , r x2 + y 2 s ! , + ∂i 4P x p h , =− ∂x (s − 2) (s − 2)2 − 4P 2 h ∂i 4P y p h , =− ∂y (s − 2) (s − 2)2 − 4P 2 h 2 ∂i . = −p ∂Ph (s − 2)2 − 4P 2 h To finalise we need to find ∂ω/∂x and ∂ω/∂y . Taking into account in what way the variables x, y were introduced (i.e., formulae (5)), the obvious relations can be written: x y , sin ω = − p . (A.9) cos ω = p x2 + y 2 x2 + y 2 Differentiating (A.9) we get: ∂ω y = 2 , ∂x x + y2 ∂ω x =− 2 . ∂y x + y2 34 V.V.Sidorenko et al. Appendix B. Calculation of the leading term W of the averaged disturbing function As it was mentioned in Sec. 4, the leading term W in the expression for disturbing function (19) and its derivatives ∂W , ∂y ∂W , ∂x ∂W ∂ϕ can be written in terms of complete elliptic integrals K (k) and E (k) of the first and second kind. B1. Calculation of W . We begin with a slightly modified version of the formula used above to define W : 1 W (ϕ, x, y ) = 2π Z π −π du , ∆(u, ϕ, x, y ) (B.1) where u = −h = λ′ − h , 2 2 2 ∆2 (u, ϕ, x, y ) = x + y + ϕ + 4ϕ(x sin u + y cos u)+ 2 2 3(x sin u + y cos u)2 + [1 − (x + y )] sin2 u. By means of the standard change of variables ξ = tg W (ϕ, x, y ) = 1 π Z +∞ −∞ p u 2 dξ P4 ( ξ ) , we obtain: (B.2) . Here P4 (ξ ) denotes the forth degree polynomial P4 (ξ ) = d4 ξ 4 + d3 ξ 3 + d2 ξ 2 + d1 ξ + d0 with the coefficients rendered by the formulae 2 d0 = (2y + ϕ)2 + x , d1 = 4x(2ϕ + 3y ), d3 = 4x(2ϕ − 3y ), 2 2 2 d2 = 2(2 + 5x − 4y + ϕ ), 2 d5 = (2y − ϕ)2 + x . The integral (B.1) has a finite value only in the case ∆(u, ϕ, x, y ) 6= 0 for ∀u ∈ [−π, π ] . It takes place when all the roots ξ1 , . . . , ξ4 of the equation P4 (ξ ) = 0 are not real: ξ1 = a1 + ib1 , ξ2 = a2 + ib2 , ξ3 = ξ 1 , ξ4 = ξ 2 (b1 > 0, b2 > 0). If a1,2 and b1,2 are known, the value of W is provided by the formula W (ϕ, x, y ) = 4 √ K π d4 κ 0 r κ1 κ0  . where κ0 = ( a1 − a2 ) 2 + (b 1 + b 2 ) 2 , κ1 = (a1 − a2 )2 + (b1 − b2 )2 . Quasi-satellite orbits. Perturbative treatment 35 B2. Calculation of the derivatives of W . By differentiation of (B.2) we arrive at the evident formulae 4 X ∂dk ∂W = Xk , ∂ϕ ∂ϕ 4 X ∂dk ∂W = Xk , ∂x ∂x k=0 4 X ∂dk ∂W = Xk , ∂y ∂y k=0 where Xk = Z +∞ −∞ ξ k dξ 3 /2 P4 (ξ ) k=0 (B.3) . Then we apply the change of the variable in (B.3) as ξ = a1 − b1 ctg σ − σ∗ (B.4) 2 The auxiliary quantity σ∗ in (B.4) is defined by the relations 2(a2 − a1 )b1 , sin σ∗ = p 2 [(a1 − a2 ) + (b1 − b2 )2 ][(a1 − a2 )2 + (b1 + b2 )2 ] cos σ∗ = p (a1 − a2 )2 + (b21 − b22 [(a1 − a2 )2 + (b1 − b2 )2 ][(a1 − a2 )2 + (b1 + b2 )2 ] , The change of variable (B.4) allows us to express X0 , . . . , X4 as linear combinations of the integrals Ik = Z π 0 cosk σdσ (α − β cos σ )3/2 (k = 0, 1, 2) with κ1 , α=1+ κ0 β=2 r κ1 . κ0 With the aid of the standard handbook Gradshteyn and Ryzhik (2007), one can easily find that r  2β 2 √ I0 = E , α+β (α − β ) α + β r  2β 2 α K , I1 = I0 − √ β α+β β α+β r  √ 2α 2β 2 α+β α2 I2 = E . I1 − 2 I0 + β β β2 α+β So the relations connecting the integrals Ik and Xk allow to express Xk in terms of complete elliptic integrals too. Omitting the cumbersome calculations we write down these relations in the recursive form X0 = 2 b21 (d4 κ0 )3/2 X 1 = a1 X 0 +   (1 + sin2 σ∗ )I0 − 2 cos σ∗ · I1 + cos 2σ∗ · I2 , 2 b1 (d4 κ0 )3/2 h i 1 sin 2σ∗ · I0 + sin σ∗ · I1 − sin 2σ∗ · I2 , 2 36 V.V.Sidorenko et al. X2 = −a21 X0 + 2a1 X1 + 2 ( d4 κ 0 )3/2   cos2 σ∗ · I0 − cos 2σ∗ · I2 , X3 = a31 X0 − 3a21 X1 + 3a1 X2 + i h 1 2b1 − sin 2σ∗ · I0 + sin σ∗ · I1 + sin 2σ∗ · I2 , 2 ( d 4 κ0 ) 3/ 2 X4 = −a41 X0 + 4a31 X1 − 6a21 X2 + 4a1 X3 + 2 b21 (d4 κ0 )3/2  (1 + sin2 σ∗ )I0 + 2 cos σ∗ · I1 + cos 2σ∗ · I2  B3. Computation of the second derivative with respect to ϕ . To make a conclusion regarding the possibility of a QS-regime at given values of x, y , we need to know 2 the sign of ∂ 2 W /∂ϕ at ϕ = 0 (Christou, 2000). The answer can be obtained through the use of the formula ∂2W ∂ϕ 2 ϕ=0 1 6 = − I0 + (x2 + y 2 )J0 π π   2 2 2 2 2 2 6 (y − x )(1 + 2x − 4y ) − 12x y J1 + p π 20x2 y 2 + (1 + 2x2 )2 + (1 − 4y 2 )2 − 1 where J0 = Z π 0 √ dξ (α − β cos ξ )5/2 J1 = α=  2 α+β = 4αE 3(α2 − β 2 )2 1 2 2 + 2(x + y ), 2 Z π 0 r 2β α+β  − (α − β )K r 2β α+β  1 cosξdξ α = J 0 − I0 , 5 / 2 β β (α − β cos ξ ) β= 1p 20x2 y 2 + (1 + 2x2 )2 + (1 − 4y 2 )2 − 1. 2 We can add also that the value of W at ϕ = 0 is provided by the formula 2 W = √ K π α+β 1 r 2β α+β  . (B.5) Up to notations the formula (B.5) coincides with the formula (46) for the secular potential in Namouni (1999) which was used there to study the limiting case of QS-orbits with non-oscillating resonant phase ϕ . , Quasi-satellite orbits. Perturbative treatment 37 References Arnold, V.I., Kozlov, V.V., Neishtadt, A.I.: Mathematical Aspects of Classical and Celestial Mechanics, 3rd Edition, Springer New York (2006) Batygin, R., Morbidelli, A.: Analytical treatment of planetary resonances. Astron. Astrophysics. 556, A28 (2013) Benest, D.: Libration effects for retrograde satellites in the restricted three-body problem. I Circular plane Hill’s case. Celest. Mechanics 13 (1976), 203-215 Benest, D.: Effects of the mass ratio on the existence of retrograde satellites in the circular restricted problem. IV. Three-dimensional stability of plane periodic orbits. Astron. Astrophysics. 54, 563-568 (1977) Brasser, R., Innanen, K.A., Connors, M., Veillet, C., Wiegert, P., Mikkola, S., Chodas, P.W.: Transient coorbital asteroids. Icarus 171, 102-109 (2004) Broucke, R.A.: Periodic orbits in the restricted three-body problem with Earth-Moon masses. NASA Technical Report 32-1168 (1968) Cary, J.R., Escande, D.F., Tennyson, J.L.: Adiabatic invariant change due to separatrix crossing. Phys. Rev. A 34, 4256-4275 (1986) Christou, A.A.: A numerical survey of transient co-orbitals of the terrestial planets. Icarus 144, 1-20 (2000) Connors, M., Veillet, C., Brasser, R., Wiegert, P., Chodas, P., Mikkola, S., Innanen, K.: Discovery of Earth’s quasi-satellite. Meteorit. Planet. Sci. 39, 1251-1255 (2004) de la Fuente Marcos, C., de la Fuente Marcos, R.: (309239) 2007RW10: a large temporary quasi-satellite of Neptune. Astron. Astrophysics 545 L9 (2012) Garfinkel, B.: Theory of the Trojan asteroids.I. Astron. J. 82, 368-379 (1977) Gil, P.J.S., Schwartz, J.: Simulations of quasi-satellite orbits around Phobos. J. of Guidance, Control, and Dynamics 33, 901-914 (2010) Gladman, B.: Dynamics of systems of two close planets. Icarus 106, 247-263 (1993) Gradshteyn, I.S., Ryzhik,I.M.: Table of integrals, series, and products, 7th Edition, Academic Press (2007) Guippone, C.A., Beaugé, C., Michtchenko, T. A., Ferraz-Mello, S.: Dynamics of two planets in co-orbital motion. MNRAS 407, 390-398 (2010) Henon, M.: Numerical exploration of the restricted problem. V. Hill’s case: periodic orbits and their stability. Astron. Astrophysics. 1, 223-238 (1969). Henon, M.: Generating Families in the Restricted Three Body Problem, Springer, Berlin Heidelberg (1997) Henrard, J.: The adiabatic invariant - The use in celestial mechanics. NATO ASI Proc. 82: Applications of Modern Dynamics to Celestial Mechanics and Astrodynamics 153-171 (1982) Jackson, J.: Retrograde satellite orbits. MNRAS 74, 62-82 (1913) Kinoshita, H., Nakai, H.: Quasi-satellites of Jupiter. Celest. Mech. Dyn. Astr. 98, 181-189 (2007) Kogan, A.Yu.: Distant satellite orbits in the restricted circular three-body problem. Cosmic Research 24, 705-710 (1989) Kogan, A.Yu.: Quasi-satellite orbits and their applications. Proceedings of the 41st Congress of the International Astronautical Federation, edited by R.Jehn. Paper 90-307 (1990) Kortenkamp, S.J.: An efficient, low-velocity, resonant mechanism for capture of satellites by a protoplanet. Icarus 175, 409-418 (2005) Lidov, M.L., Vashkovyak, M.A.: Perturbation theory and analysis of the evolution of quasisatellite orbits in the restricted three-body problem. Cosmic Research 31, 187-207 (1993) Lidov, M.L., Vashkovyak, M.A.: On quasi-satellite orbits in the restricted elliptic three-body problem. Astronomy Letters 20, 676-690 (1994) Lidov, M.L., Ziglin, S.L.: The analysis of restricted circular twice-averaged three body problem in the case of close orbits. Celest. Mechanics 9, 151-173 (1974) Mikkola, S., Innanen, K.: Orbital stability of planetary quasi-satellites. In: The Dynamical Behavior of Our Planetary System (R.Dvorak and J.Henrard, Eds). Kluwer, Boston, 345355 (1997) Mikkola, S., Brasser, R., Wiegert, P., Innanen, K.: Asteroid 2002VE68, a quasi-satellite of Venus. MNRAS 351, L63-65 (2004) Mikkola, S., Innanen, K., Wiegert, P., Connors, M., Brasser, R.: Stability limits for the quasisatellite orbit. MNRAS 369, 15-24 (2006) 38 V.V.Sidorenko et al. Murray, C.D., Dermott, S.F.: Solar System Dynamics, Cambridge University Press (1999) Namouni, F.: Secular interactions of coorbiting objects. Icarus 137, 293-314 (1999) Namouni, F., Christou, A.A., Murray, C.D.: Coorbital dynamics at large eccentricity and inclination. Phys. Rev. Lett. 83, 2506-2509 (1999) Neishtadt, A.I.: Change of an adiabatic invariant at a separatrix. Soviet Journal of Plasma Physics 12, 568-573 (1986) Neishtadt, A.I.: Jumps of the adiabatic invariant on crossing the separatrix and the origin of the 3:1 Kirkwood gap. Sov. Phys. Dokl. 32, 571-573 (1987) Neishtadt, A.I.: On the change in the adiabatic invariant on crossing a separatrix in systems with two degrees of freedom. J. Appl. Math. Mech.-USSR 51, 586-592 (1987) Neishtadt, A.I., Sidorenko, V.V.: Wisdom system: dynamics in the adiabatic approximation. Celest. Mech. Dyn. Astron. 90, 307-330 (2004) Nesvorny, D., Thomas, F., Ferraz-Mello, S., Morbidelli,A.: A perturbative treatment of the co-orbital motion. Celest. Mech. Dyn. Astron. 82, 323-361 (2002) Robutel, P., Pousse, A.: On the co-orbital motion of two planets in quasi-circular orbits. Celest. Mech. Dyn. Astron. 117, 17-40 (2013) Sagdeev, R.Z., Zakharov, A.V.: Brief history of the Phobos mission, Nature 341, 581-585 (1989) Sidlichovsky, M.: A non-planar circular model for the 4/7 resonance. Celest. Mech. Dyn. Astron. 93, 167-185 (2005) Sidorenko, V.V.: Evolution of asteroid orbits at resonance 3 : 1 of their mean motions with Jupiter (planar problem). Cosmic Research 44, 440-455 (2006) Schubart, J.: Long-period effects in nearly commensurable cases of the restricted three body problem. Smithson. Astrophys. Obs. Spec. Report 149 (1964) Timofeev, A.V.: On the problem of constancy of the adiabatic invariant during change of the type of motion. Sov. Phys., JETP 75, 1303-1308 (1978) Tuchin, A.G.: Quasi-synchronous orbits and their employment for the approach of a spacecraft to Phobos, Cosmic Research 45, 131-136 (2007) Vashkovyak, M.A.: Numerical averaging method in the problem of the evolution of resonance asteroid orbits. Cosmic Research 27, 6-11 (1989) Wajer, P.: 2002AA29: Earth’s recurrent quasi-satellite? Icarus 200, 147-153 (2009) Wajer, P.: Dynamical evolution of Earth’s quasi-satellites: 2004GU9 and 2006FV35. Icarus 209, 488-493 (2010) Wiegert, P.A., Innanen, K.A., Mikkola, S.: The orbital evolution of near-Earth asteroid 3753. Astron. J. 115, 2604-2613 (1998) Wisdom, J.: A perturbative treatment of motion near the 3/1 commensurability. Icarus 63, 272-286 (1985) Yokoyama, T.: A simple generalization of Wisdom‘s perturbative method. Celest. Mech. Dyn. Astron. 64, 243-260 (1996)