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

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 quasi-satellite 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.

Appendices

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

$$\begin{aligned} W(\varphi ,x,y,P_h)=\frac{1}{2\pi }\int _0^{2\pi }R(\varphi ,x,y,P_h,\overline{h})\,d\overline{h}\end{aligned}$$

(31)

and of its derivatives

$$\begin{aligned} \frac{\partial W}{\partial \varphi },\quad \frac{\partial W}{\partial x},\quad \frac{\partial W}{\partial y}, \quad \frac{\partial W}{\partial P_{h}} \end{aligned}$$

with no limitations on the asteroid’s eccentricity or inclination.

A1.Change of the variable, over which the integration takes place. To start with, we replace in (31) the integration over the auxiliary variable $\overline{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 $\overline{h}$ and $E$ reads as

$$\begin{aligned} \overline{h}=\varphi - \omega - (E-e\sin E). \end{aligned}$$

(32)

As a result, the expressions for $W$ and its derivative can be written as

$$\begin{aligned} W(\varphi ,x,y,P_h)&= -\frac{1}{2\pi }\int _0^{2\pi }R\frac{\partial \overline{h}}{\partial E}\,dE\nonumber \\ \frac{\partial W}{\partial \varphi }&= -\frac{1}{2\pi }\int _0^{2\pi }\frac{\partial R}{\partial \varphi }\frac{\partial \overline{h}}{\partial E}\,dE,\quad \frac{\partial W}{\partial P_{h}}=-\frac{1}{2\pi }\int _0^{2\pi }\frac{\partial R}{\partial P_{h}}\frac{\partial \overline{h}}{\partial E}\,dE,\nonumber \\ \frac{\partial W}{\partial x}&= -\frac{1}{2\pi }\int _0^{2\pi }\left( \frac{\partial R}{\partial x}\frac{\partial \overline{h}}{\partial E} + R\frac{\partial ^2 \overline{h}}{\partial x \partial E}\right) \,dE,\nonumber \\ \frac{\partial W}{\partial y}&= -\frac{1}{2\pi }\int _0^{2\pi } \left( \frac{\partial R}{\partial y}\frac{\partial \overline{h}}{\partial E} + R\frac{\partial ^2 \overline{h}}{\partial y \partial E}\right) \,dE, \end{aligned}$$

(33)

where

$$\begin{aligned} \frac{\partial \overline{h}}{\partial E}&= -(1 - e\sin E),\\ \frac{\partial ^2 \overline{h}}{\partial E \partial x}&= \frac{\partial ^2 \overline{h}}{\partial E \partial e}\cdot \frac{\partial e}{\partial x}= \cos E \cdot \frac{\partial e}{\partial x},\\ \frac{\partial ^2 \overline{h}}{\partial E \partial y}&= \frac{\partial ^2 \overline{h}}{\partial E \partial e}\cdot \frac{\partial e}{\partial y}= \cos E \cdot \frac{\partial e}{\partial y}. \end{aligned}$$

The “minus” sign before the integrals in formulae (33) appears due to our intention to have the upper limit larger than the lower one.

A2. Derivatives of $R$ with respect to $\varphi ,x,y,P_{h}$. At this point, we need to introduce the uniformly rotating heliocentric reference frame $O\xi \eta \zeta$ with the axis $O\xi$ being directed from the Sun to the planet and the axis $O\zeta$ 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 $\mathbf{{r}}'=\mathbf{{e}}_\xi$, where $\mathbf{{e}}_\xi =(1,0,0)^T$ is the unit vector corresponding to axis $O\xi$. Respectively the expression for disturbing function (2) is reduced to

$$\begin{aligned} R=\frac{1}{|\mathbf{{r}}- \mathbf{{e}}_\xi |}-(\mathbf{{r}},\mathbf{{e}}_\xi ). \end{aligned}$$

Then we get:

$$\begin{aligned} \frac{\partial R}{\partial \varphi }=\frac{\partial R}{\partial \mathbf{{r}}}\frac{\partial \mathbf{{r}}}{\partial \varphi },\quad \frac{\partial R}{\partial x}=\frac{\partial R}{\partial \mathbf{{r}}}\frac{\partial \mathbf{{r}}}{\partial x},\quad \frac{\partial R}{\partial y}=\frac{\partial R}{\partial \mathbf{{r}}}\frac{\partial \mathbf{{r}}}{\partial y}, \quad \frac{\partial R}{\partial P_{h}}=\frac{\partial R}{\partial \mathbf{{r}}}\frac{\partial \mathbf{{r}}}{\partial P_{h}}, \end{aligned}$$

where

$$\begin{aligned} \frac{\partial R}{\partial \mathbf{{r}}}= - \frac{\mathbf{{r}}- \mathbf{{e}}_\xi }{|\mathbf{{r}}- \mathbf{{e}}_\xi |^3}- \mathbf{{e}}_\xi . \end{aligned}$$

A3. Position vector of the asteroid. Let the unit vector $\mathbf{{e}}^*_\xi$ be directed to the pericenter of the osculating orbit, while one more unit vector $\mathbf{{e}}^*_\eta$ 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\xi \eta \zeta$:

$$\begin{aligned} \mathbf{{e}}^*_\xi = \left( \begin{array}{c} \cos \overline{h}\cos \omega -\cos i \sin \overline{h}\sin \omega \\ \sin \overline{h}\cos \omega +\cos i \cos \overline{h}\sin \omega \\ \sin i\sin \omega \end{array}\right) \end{aligned}$$

(34)

and

$$\begin{aligned} \mathbf{{e}}^*_\eta = \left( \begin{array}{c} -\cos \overline{h}\sin \omega -\cos i \sin \overline{h}\cos \omega \\ -\sin \overline{h}\sin \omega +\cos i \cos \overline{h}\cos \omega \\ \sin i\cos \omega \end{array}\right) \end{aligned}$$

(35)

After that the asteroid’s position vector $\mathbf{{r}}$ can be written as the linear combination

$$\begin{aligned} \mathbf{{r}}(i,e,\omega ,\overline{h},E) = \mathbf{{e}}^*_\xi (i,\omega ,\overline{h}) \xi ^*(e,E)+\mathbf{{e}}^*_\eta (i,\omega ,\overline{h}) \eta ^*(e,E) \end{aligned}$$

(36)

with the coefficients (here the orbit with $a=1$ is considered)

$$\begin{aligned} \xi ^*(e,E)=\cos E - e,\qquad \eta ^*(e,E)=\sqrt{1-e^2}\sin E. \end{aligned}$$

A4. Derivatives of the position vector $\mathbf{{r}}$ with respect to $\varphi ,x,y,P_{h}$. Taking into account (32) and (36) we obtain the following relations:

$$\begin{aligned} \frac{\partial \mathbf{{r}}}{\partial x}&= \left( \frac{\partial \mathbf{{r}}}{\partial \overline{h}}\frac{\partial \overline{h}}{\partial e}+\frac{\partial \mathbf{{r}}}{\partial e}\right) \frac{\partial e}{\partial x}+ \left( \frac{\partial \mathbf{{r}}}{\partial \overline{h}}\frac{\partial \overline{h}}{\partial \omega }+\frac{\partial \mathbf{{r}}}{\partial \omega }\right) \frac{\partial \omega }{\partial x}+ \frac{\partial \mathbf{{r}}}{\partial i}\frac{\partial i}{\partial x},\\ \frac{\partial \mathbf{{r}}}{\partial y}&= \left( \frac{\partial \mathbf{{r}}}{\partial \overline{h}}\frac{\partial \overline{h}}{\partial e}+\frac{\partial \mathbf{{r}}}{\partial e}\right) \frac{\partial e}{\partial y}+ \left( \frac{\partial \mathbf{{r}}}{\partial \overline{h}}\frac{\partial \overline{h}}{\partial \omega }+\frac{\partial \mathbf{{r}}}{\partial \omega }\right) \frac{\partial \omega }{\partial y}+ \frac{\partial \mathbf{{r}}}{\partial i}\frac{\partial i}{\partial y},\\ \frac{\partial \mathbf{{r}}}{\partial \varphi }&= \frac{\partial \mathbf{{r}}}{\partial \overline{h}}\frac{\partial \overline{h}}{\partial \varphi }=\frac{\partial \mathbf{{r}}}{\partial \overline{h}},\quad \frac{\partial \mathbf{{r}}}{\partial P_{h}}=\frac{\partial \mathbf{{r}}}{\partial i}\frac{\partial i}{\partial P_{h}}, \end{aligned}$$

where evidently

$$\begin{aligned} \frac{\partial \overline{h}}{\partial \omega }=-1,\qquad \frac{\partial \overline{h}}{\partial e}=\sin E. \end{aligned}$$

A5. Derivatives of the position vector $\mathbf{{r}}$ with respect to $\overline{h}$ and osculating variables $i,e,\omega$. Direct differentiation of (36) results in

$$\begin{aligned} \frac{\partial \mathbf{{r}}}{\partial \overline{h}}&= \frac{\partial \mathbf{{e}}^*_\xi }{\partial \overline{h}}\xi ^* +\frac{\partial \mathbf{{e}}^*_\eta }{\partial \overline{h}}\eta ^*,\qquad \frac{\partial \mathbf{{r}}}{\partial i}=\frac{\partial \mathbf{{e}}^*_\xi }{\partial i}\xi ^* +\frac{\partial \mathbf{{e}}^*_\eta }{\partial i}\eta ^*,\\ \frac{\partial \mathbf{{r}}}{\partial \omega }&= \frac{\partial \mathbf{{e}}^*_\xi }{\partial \omega }\xi ^* +\frac{\partial \mathbf{{e}}^*_\eta }{\partial \omega }\eta ^*, \qquad \frac{\partial \mathbf{{r}}}{\partial e}=\mathbf{{e}}^*_\xi \frac{\partial \xi ^*}{\partial e}+\mathbf{{e}}^*_\eta \frac{\partial \eta ^*}{\partial e},\nonumber \end{aligned}$$

(37)

where

$$\begin{aligned} \frac{\partial \xi ^*}{\partial e}=1,\qquad \frac{\partial \eta ^*}{\partial e}=-\frac{e \sin E}{\sqrt{1-e^2}}. \end{aligned}$$

The derivatives of the unit vectors $\mathbf{{e}}^*_\xi$ and $\mathbf{{e}}^*_\eta$ in (37) can be easily evaluated from (34) and (35) respectively.

A6. Derivatives of the osculating variables $e,i,\omega$ with respect to $x,y,P_{h}$. First we rewrite the expressions for $e,\cos i,\sin i$ in a more compact way [compare with (6)]:

$$\begin{aligned} e&= \frac{\sqrt{s(x^2+y^2)}}{2},\nonumber \\ \cos i&= \frac{2P_h}{s-2},\qquad \sin i = \frac{\sqrt{(s-2)^2 - 4P_h^2}}{s-2}. \end{aligned}$$

(38)

Here

$$\begin{aligned} s=4 - (x^2 + y^2). \end{aligned}$$

After the differentiation of the expressions (38) with respect to $x,y$ we obtain:

$$\begin{aligned} \frac{\partial e}{\partial x}&= \frac{x}{2}\left( \sqrt{\frac{s}{x^2 +y^2}}-\sqrt{\frac{x^2+y^2}{s}}\right) ,\\ \frac{\partial e}{\partial y}&= \frac{y}{2}\left( \sqrt{\frac{s}{x^2 +y^2}}-\sqrt{\frac{x^2+y^2}{s}}\right) , \end{aligned}$$

and

$$\begin{aligned} \frac{\partial i}{\partial x}&= -\frac{4P_h x}{(s-2)\sqrt{(s-2)^2 - 4P_h^2}},\\ \frac{\partial i}{\partial y}&= -\frac{4P_h y}{(s-2)\sqrt{(s-2)^2 - 4P_h^2}},\\ \frac{\partial i}{\partial P_{h}}&= -\frac{2}{\sqrt{(s-2)^2 - 4P_h^2}}. \end{aligned}$$

To finalise we need to find $\partial \omega /\partial x$ and $\partial \omega /\partial y$. Taking into account in what way the variables $x,y$ were introduced [i.e., formulae (5)], the obvious relations can be written:

$$\begin{aligned} \cos \omega = \frac{x}{\sqrt{x^2 + y^2}},\qquad \sin \omega = -\frac{y}{\sqrt{x^2 + y^2}}. \end{aligned}$$

(39)

Differentiating (39) we get:

$$\begin{aligned} \frac{\partial \omega }{\partial x}= \frac{y}{x^2+y^2},\qquad \frac{\partial \omega }{\partial y}=-\frac{x}{x^2+y^2}. \end{aligned}$$

Appendix B: Calculation of the leading term $\overline{W}$ of the averaged disturbing function

As it was mentioned in Sect. 4, the leading term $\overline{W}$ in the expression for disturbing function (19) and its derivatives

$$\begin{aligned} \frac{\partial \overline{W}}{\partial \overline{x}},\quad \frac{\partial \overline{W}}{\partial \overline{y}},\quad \frac{\partial \overline{W}}{\partial \overline{\varphi }} \end{aligned}$$

can be written in terms of complete elliptic integrals $K(k)$ and $E(k)$ of the first and second kind.

B1. Calculation of $\overline{W}$. We begin with a slightly modified version of the formula used above to define $\overline{W}$:

$$\begin{aligned} \overline{W}(\overline{\varphi },\overline{x},\overline{y})=\frac{1}{2\pi }\int ^{\pi }_{-\pi }\frac{du}{\varDelta (u, \overline{\varphi },\overline{x},\overline{y})}, \end{aligned}$$

(40)

where $u=-\overline{h}=\lambda '-h$,

$$\begin{aligned} \varDelta ^2(u,\overline{\varphi },\overline{x},\overline{y})=\overline{x}^2+\overline{y}^2+\overline{\varphi }^2 + 4\overline{\varphi }(\overline{x}\sin u + \overline{y}\cos u)\\ +3(\overline{x}\sin u + \overline{y}\cos u)^2 + [1-(\overline{x}^2+\overline{y}^2)]\sin ^2u. \end{aligned}$$

By means of the standard change of variables $\xi = \mathop {\mathrm{tg}}\nolimits \frac{u}{2}$, we obtain:

$$\begin{aligned} \overline{W}(\overline{\varphi },\overline{x},\overline{y})=\frac{1}{\pi }\int ^{+\infty }_{-\infty } \frac{d\xi }{\sqrt{P_4(\xi )}}. \end{aligned}$$

(41)

Here $P_4(\xi )$ denotes the forth degree polynomial

$$\begin{aligned} P_4(\xi )=d_4\xi ^4 + d_3\xi ^3 + d_2\xi ^2 + d_1\xi + d_0 \end{aligned}$$

with the coefficients rendered by the formulae

$$\begin{aligned} d_0&= (2\overline{y}+\overline{\varphi })^2 + \overline{x}^2,\quad d_1=4\overline{x}(2\overline{\varphi }+ 3\overline{y}),\quad d_2=2(2+5\overline{x}^2-4\overline{y}^2+\overline{\varphi }^2),\\ d_3&= 4\overline{x}(2\overline{\varphi }- 3\overline{y}),\quad d_5=(2\overline{y}-\overline{\varphi })^2 + \overline{x}^2. \end{aligned}$$

The integral (40) has a finite value only in the case $\varDelta (u,\overline{\varphi },\overline{x},\overline{y})\ne 0$ for $\forall u\in [-\pi ,\pi ]$. It takes place when all the roots $\xi _1,\ldots ,\xi _4$ of the equation $P_4(\xi )=0$ are not real:

$$\begin{aligned} \xi _1=a_1+ib_1,\quad \xi _2=a_2+ib_2,\quad \xi _3= \overline{\xi }_1, \quad \xi _4=\overline{\xi }_2\quad (b_1>0,b_2>0). \end{aligned}$$

If $a_{1,2}$ and $b_{1,2}$ are known, the value of $\overline{W}$ is provided by the formula

$$\begin{aligned} \overline{W}(\overline{\varphi },\overline{x},\overline{y})=\frac{4}{\pi \sqrt{d_4\kappa _0}} K\left( \sqrt{\frac{\kappa _1}{\kappa _0}}\right) . \end{aligned}$$

where

$$\begin{aligned} \kappa _0 = (a_1-a_2)^2+(b_1+b_2)^2,\qquad \kappa _1 = (a_1-a_2)^2+(b_1-b_2)^2. \end{aligned}$$

B2. Calculation of the derivatives of $\overline{W}$. By differentiation of (41) we arrive at the evident formulae

$$\begin{aligned} \frac{\partial \overline{W}}{\partial \overline{\varphi }}= \sum _{k=0}^{4}\frac{\partial d_k}{\partial \overline{\varphi }}X_k,\quad \frac{\partial \overline{W}}{\partial \overline{x}}= \sum _{k=0}^{4}\frac{\partial d_k}{\partial \overline{x}}X_k,\quad \frac{\partial \overline{W}}{\partial \overline{y}}= \sum _{k=0}^{4}\frac{\partial d_k}{\partial \overline{y}}X_k, \end{aligned}$$

where

$$\begin{aligned} X_k=\int ^{+\infty }_{-\infty }\frac{\xi ^k d\xi }{P_4^{3/2}(\xi )}. \end{aligned}$$

(42)

Then we apply the change of the variable in (42) as

$$\begin{aligned} \xi = a_1 - b_1\mathop {\mathrm{ctg}}\nolimits \frac{\sigma - \sigma _*}{2} \end{aligned}$$

(43)

The auxiliary quantity $\sigma _*$ in (43) is defined by the relations

$$\begin{aligned} \sin \sigma _*&= \frac{2(a_2-a_1)b_1}{\sqrt{[(a_1-a_2)^2 +(b_1-b_2)^2][(a_1-a_2)^2+(b_1+b_2)^2]}},\\ \cos \sigma _*&= \frac{(a_1-a_2)^2+(b_1^2-b_2^2}{\sqrt{[(a_1-a_2)^2+(b_1-b_2)^2][(a_1-a_2)^2+(b_1+b_2)^2]}}, \end{aligned}$$

The change of variable (43) allows us to express $X_0,\ldots ,X_4$ as linear combinations of the integrals

$$\begin{aligned} I_k=\int _0^{\pi }\frac{\cos ^k \sigma d\sigma }{(\alpha -\beta \cos \sigma )^{3/2}}\quad (k=0,1,2) \end{aligned}$$

with

$$\begin{aligned} \alpha =1+\frac{\kappa _1}{\kappa _0},\qquad \beta = 2\sqrt{\frac{\kappa _1}{\kappa _0}}. \end{aligned}$$

With the aid of the standard handbook Gradshteyn and Ryzhik (2007), one can easily find that

$$\begin{aligned} I_0&= \frac{2}{(\alpha -\beta )\sqrt{\alpha +\beta }} E\left( \sqrt{\frac{2\beta }{\alpha +\beta }}\right) ,\\ I_1&= \frac{\alpha }{\beta }I_0-\frac{2}{\beta \sqrt{\alpha +\beta }} K\left( \sqrt{\frac{2\beta }{\alpha +\beta }}\right) ,\\ I_2&= \frac{2\alpha }{\beta }I_1-\frac{\alpha ^2}{\beta ^2}I_0+ \frac{2\sqrt{\alpha +\beta }}{\beta ^2} E\left( \sqrt{\frac{2\beta }{\alpha +\beta }}\right) . \end{aligned}$$

So the relations connecting the integrals $I_k$ and $X_k$ allow to express $X_k$ in terms of complete elliptic integrals too. Omitting the cumbersome calculations we write down these relations in the recursive form

$$\begin{aligned} X_0&= \frac{2}{b_1^2(d_4\kappa _0)^{3/2}}\left[ (1+\sin ^2 \sigma _*)I_0 - 2\cos \sigma _*\cdot I_1 + \cos 2\sigma _*\cdot I_2\right] ,\\ X_1&= a_1 X_0+\frac{2}{b_1(d_4\kappa _0)^{3/2}}\left[ \frac{1}{2} \sin 2\sigma _*\cdot I_0 + \sin \sigma _*\cdot I_1 - \sin 2\sigma _* \cdot I_2\right] ,\\ X_2&= -a_1^2 X_0 + 2a_1 X_1 + \frac{2}{(d_4\kappa _0)^{3/2}}\left[ \cos ^2\sigma _*\cdot I_0 - \cos 2\sigma _*\cdot I_2\right] ,\\ X_3&= a_1^3 X_0 - 3a_1^2 X_1 + 3a_1 X_2\\&+\frac{2b_1}{(d_4\kappa _0)^{3/2}}\left[ -\frac{1}{2}\sin 2 \sigma _*\cdot I_0 + \sin \sigma _*\cdot I_1 + \sin 2\sigma _*\cdot I_2\right] ,\\ X_4&= -a_1^4 X_0 + 4a_1^3 X_1 - 6a_1^2 X_2 + 4a_1 X_3\\&+\frac{2}{b_1^2(d_4\kappa _0)^{3/2}}\left[ (1+\sin ^2 \sigma _*)I_0 + 2\cos \sigma _*\cdot I_1 + \cos 2\sigma _*\cdot I_2\right] \end{aligned}$$

B3. Computation of the second derivative with respect to $\overline{\varphi }$. To make a conclusion regarding the possibility of a QS-regime at given values of $\overline{x},\overline{y}$, we need to know the sign of $\partial ^2 \overline{W}/\partial \overline{\varphi }^2$ at $\overline{\varphi }=0$ (Christou 2000). The answer can be obtained through the use of the formula

$$\begin{aligned} \left. \frac{\partial ^2 \overline{W}}{\partial \overline{\varphi }^2}\right| _{\overline{\varphi }=0}&= -\frac{1}{\pi }I_0 + \frac{6}{\pi }(x^2+y^2)J_0\\&+\frac{6}{\pi }\frac{\left[ (y^2-x^2)(1+2x^2-4y^2)-12x^2y^2\right] J_1}{\sqrt{20x^2y^2+(1+2x^2)^2+(1-4y^2)^2-1}} \end{aligned}$$

where

$$\begin{aligned} J_0&= \int _0^{\pi }\frac{d\xi }{(\alpha - \beta \cos \xi )^{5/2}} = \frac{2\sqrt{\alpha +\beta }}{3(\alpha ^2-\beta ^2)^2} \left[ 4\alpha E\left( \sqrt{\frac{2\beta }{\alpha +\beta }}\right) -(\alpha -\beta ) K\left( \sqrt{\frac{2\beta }{\alpha +\beta }}\right) \right] ,\\ J_1&= \int _0^{\pi }\frac{cos\xi d\xi }{(\alpha - \beta \cos \xi )^{5/2}} = \frac{\alpha }{\beta }J_0 - \frac{1}{\beta }I_0,\\ \alpha&= \frac{1}{2}+2(\overline{x}^2+\overline{y}^2),\qquad \beta =\frac{1}{2}\sqrt{20x^2y^2+(1+2x^2)^2+(1-4y^2)^2-1}. \end{aligned}$$

We can add also that the value of $\overline{W}$ at $\overline{\varphi }=0$ is provided by the formula

$$\begin{aligned} \overline{W}= \frac{1}{\pi }\frac{2}{\sqrt{\alpha + \beta }} K\left( \sqrt{\frac{2\beta }{\alpha +\beta }}\right) . \end{aligned}$$

(44)

Up to notation, the formula (44) 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 $\varphi$.