Existence and Control of Special Orbits around Asteroid 4 Vesta

Abstract

This paper focuses on the existence and control of particular types of orbits around asteroid 4 Vesta, including Sun-synchronous orbits, orbits at the critical inclination, repeating ground-track orbits, and stationary orbits. J₂, J₃, and J₄ terms are considered in the gravity model of Vesta. First, the inclination perturbation caused by solar gravitation is studied, and preset and multiple inclination bias methods are proposed to dampen the local time drift at the ascending node. Compared with Vesta, the control periods of the Sun-synchronous orbits of 21 Lutetia and 433 Eros are much longer. Second, Vesta’s orbits with a critical inclination depend on the semi-major axis and eccentricity. If the eccentricity is not greater than 0.2, inclination decreases slowly and monotonically concerning the semi-major axis. If the eccentricity is not smaller than 0.4, inclination increases rapidly and monotonically. Third, Sun-synchronous repeating ground-track circular orbits of Vesta, which do not exist for Lutetia and Eros, are investigated. Finally, the perturbations of stationary orbits caused by solar gravitation and solar radiation pressure are analyzed.

1. Introduction

Asteroids are terrestrial bodies mainly in the Main Asteroid Belt of our Solar System, primarily located between Mars and Jupiter. The Asteroid Belt (~2.1–3.3 au) is estimated to have about $(1.4\pm 0.5)\times {10}^6$ bodies with a diameter greater than about 1 km [1,2]. Scientists can learn about the origin and evolution of our solar system by studying the geophysics and geochemistry of asteroids. Vesta is the second-largest asteroid and, supposedly, the parent body of meteorites in the Main Asteroid Belt. Coradini et al. [3] theorized that the thermal evolution of Vesta took place at a time contemporary to the formation of Jupiter. There is an assumption that HED (Howardites, Eucrites, and Diogenites) meteorites found on Earth originated from Vesta [4].

Vesta was the fourth body discovered in the Main Asteroid Belt in 1807; observations from NASA’s Hubble Space Telescopes reveal some hints of its surface composition. In recent decades, an increasing number of missions were designed to explore the asteroids in the solar system, including the sample return missions, the OSIRIS-REx (2016) [5] to asteroid Bennu mission, and the Hayabusa2 (2014) [6,7] to asteroid Ryugu mission. The most well-known mission to Vesta is NASA’s Dawn mission [8,9], which was launched on 27 September 2007 and finally entered an orbit around Vesta in July 2011. The spacecraft spent 14 months in a series of orbits of different Vesta altitudes, studying the surface and probing the interior through gravity measurements [10]. Dawn’s observation of Vesta reveals that the surface is dominated by variously sized impact craters [11,12] at a survey orbit altitude of ~2700 km. The gravity field was investigated based on collected data from the Deep Space Network Doppler tracking and optical landmark images [13].

In the past few years, much research on various special spacecraft orbits around planets has been conducted, especially on spacecraft orbits of Earth, Mars, and Jupiter [14,15,16,17,18,19,20]. Usually, these special orbits include stationary orbits, frozen orbits, Sun-synchronous orbits, repeating ground-track orbits, and orbits at the critical inclination. Li et al. [21] proposed an analytic method for frozen orbit design and discussed the properties of frozen orbits around Vesta. Since an asteroid usually has a body with an irregular shape, its gravity field is very different from Earth. Werner et al. [22] introduced a polyhedral method to evaluate the gravitation of irregularly shaped bodies, but the integration is computationally demanding. Moreover, many researchers have conducted much work on equilibrium points of irregularly shaped asteroids [23,24,25,26,27]. Although the harmonic expansions of the gravitational field are always an approximation to reality, the series can be truncated at a limited order to obtain accuracy in modeling. As far as is known, there are few studies about orbits around Vesta using the mean element theory. This paper aims to investigate the existence and control of special orbits around Vesta and mainly considers the perturbation of the non-spherical gravitational effect. It also applies these methods to Lutetia [28] and Eros [29] and presents some interesting results.

The basic parameters of Vesta are listed in

Table 1. The basic parameters of Vesta.

$\mu$ (km3 s−2)	17.288245
Reference radius R (km)	265.0
J2	$7.1060892\times {10}^{-2}$
J3	$-8.7588999\times {10}^{-3}$
J4	$-9.7967997\times {10}^{-3}$
J5	$3.9871881\times {10}^{-3}$
Obliquity of the orbit (deg)	15.66
Mean motion (deg/day)	0.271587
Rotation rate (deg/day)	1617.333128

[13,30,31]. One can see that the term J₂ of Vesta is dominant, but the harmonic coefficients of J₃, J₄, and J₅ are also strong for Vesta. In Earth’s gravity field, the J₂ term is even 1000 times stronger than J₃. The spacecraft’s motion around Vesta will differ due to these gravity properties.

In the inertial frame of the rotational symmetric Vesta model, the gravitational function is given by [32,33]:

$$V=\frac{\mu }{r}\left[1-{{\displaystyle {\sum }_{n=2}^{\infty }{J}_n\left(\frac{R}{r}\right)}}^n{P}_n(\sin \varphi )\right]$$

(1)

where $\mu$ is the gravitational constant of Vesta, R is the radius of Vesta, r is the distance of the spacecraft from the mass center of Vesta, $\varphi$ is the latitude of the spacecraft, and $P_k(k\ge 2)$ is the Legendre polynomial of degree k. When we consider the analytical aspheric perturbations with zonal harmonics up to J₆, we can formalize the gravitational function as follows:

$$\begin{array}{cc}V=& \frac{\mu }{r}{[1-J_2\left(\frac{R}{r}\right)}^2\left(\frac{3}{2}{\sin }^2\varphi -\frac{1}{2}\right)-J_3{\left(\frac{R}{r}\right)}^3\left(\frac{5}{2}{\sin }^3\varphi -\frac{3}{2}\sin \varphi \right)\hfill \\ & -J_4{\left(\frac{R}{r}\right)}^4\left(\frac{35}{8}{\sin }^4\varphi -\frac{15}{4}{\sin }^2\varphi +\frac{3}{8}\right)\hfill \\ & -J_5{\left(\frac{R}{r}\right)}^5\left(\frac{63}{8}{\sin }^5\varphi -\frac{35}{4}{\sin }^3\varphi +\frac{15}{8}\sin \varphi \right)\hfill \\ & -J_6{\left(\frac{R}{r}\right)}^6\left(\frac{231}{16}{\sin }^6\varphi -\frac{315}{16}{\sin }^4\varphi +\frac{105}{16}{\sin }^2\varphi -\frac{5}{16}\right)]\hfill \end{array}$$

(2)

It is important to determine the region where the gravitational influence of a planet prevails over the influence of other planets, for this need scientists introduced the concept of the sphere of influence. This is defined as the region of space in which the planet represents the primary gravitational source. For a spacecraft, P₂, and two planets, P₁ and P₃, in a three-body system, the radius of the sphere of influence of P₁ to P₃ is [14]:

$$r=\rho {\left(\frac{m_1}{m_3}\right)}^{2/5}$$

(3)

where $m_1$ and $m_3$ are the mass of P₁ and P₃, respectively, and $\rho$ is the vector from P₁ to P₃. When P₂ is in the sphere of influence of P₁, the gravitation of P₁ is the main gravity center for P₂. Otherwise, P₃ is the gravity center for P₂. In the solar system, the radius of the sphere of influence of Earth is about $93\times {10}^4$ km from the Sun. For Vesta, the radius of the influence sphere is roughly about $3.9\times {10}^4$ km from the Sun. That is, the asteroid’s gravity is primary compared with the solar gravity at a high altitude of about 1000 km from the mass center of Vesta.

2. Sun-Synchronous Orbits

Sun-synchronous orbits are orbits with a precession rate of the orbital plane equal to the revolution’s angular velocity around the Sun. The spacecraft passes over a given location on the planet each time at the same local solar time so that the remote sensing and optical observation satellites are launched into these orbits. In other words, the rate of node precession is equal to the Vesta mean motion of orbit around the Sun; that is:

$$\dot{\Omega }=n_s,$$

(4)

where $\dot{\Omega }$ is as formulated in Equation (A.1) in Appendix A. Equation (4) can be further written as the following equation:

$$f(\cos i)={\theta }_1{\cos }^{3}i+{\theta }_2\cos i+{\theta }_3=0,$$

(5)

where

$${\theta }_1=\frac{9n{J}_2^2{R}^4}{4p^4}\left[-\frac{35J_4}{6J_2^2}\left(\frac{3}{4}{e}^2+\frac{1}{2}\right)+\frac{3}{2}\sqrt{1-e^2}-\frac{5}{24}{e}^2+\frac{5}{3}\right],$$

(6)

$${\theta }_2=\frac{3n{J}_2{R}^2}{2p^2}\left\{\frac{3J_2{R}^2}{2p^2}\left[\frac{5J_4}{2J_2^2}\left(\frac{3e^2}{4}+\frac{1}{2}\right)-\frac{1}{2}\sqrt{1-e^2}+\frac{3e^2}{8}-\frac{1}{6}\right]+1\right\},$$

(7)

$${\theta }_3=n_s.$$

(8)

and $n_s$ is the mean motion of Vesta around the Sun. We note that Equation (5) is a cubic function of $\cos i$; it has three real roots if the discriminant satisfies the following condition:

$$\Delta =\frac{{\theta }_3^2}{4{\theta }_1^2}+\frac{{\theta }_2^2}{27{\theta }_1^3}\le 0.$$

(9)

If $\Delta >0$, then Equation (5) has one real root and two complex roots. The semi-major axis, a, and eccentricity, e, should satisfy the following:

$$R<a(1-e).$$

(10)

Figure 1a shows the values of the discriminant for a and e. We see the value of the discriminant is always non-negative, so Equation (5) has one real root when a and e are given. Figure 1a indicates that there exists one Sun-synchronous orbit around Vesta at a certain a and e. Figure 1b shows the variations of inclination, i, for different values of a and e.

The spacecraft in Sun-synchronous orbits of Vesta would be perturbed by the Sun. In the latter part of this section, we continue to investigate circular Sun-synchronous-orbit perturbations. The derivative of inclination is formulated using the Lagrange equations, and the secular term of the derivative of inclination in a period caused by solar gravity is [20]:

$$\frac{\mathrm{d}i}{\mathrm{d}t}=\frac{3n_s^2}{8n}\left(\sin i\sin 2\Omega {\sin }^2{i}_s+\cos i\sin \Omega \sin 2i_s\right)$$

(11)

where $i_s$ is the angle between the motion of Vesta around the Sun and the equator of Vesta. For circular Sun-synchronous orbits, Equation (11) is presented as:

$$\frac{\mathrm{d}i}{\mathrm{d}t}=-\frac{3n_s^2}{16n}\sin i{(1+\cos i_s)}^2\sin (2{\beta }_s-2\Omega )$$

(12)

where ${\beta }_s$ is the ecliptic longitude of motions of the Sun. The phrase, ${\beta }_s-\Omega$, in Equation (12) is related to the local time at the descending node, and the derivation of inclination is determined since the local time at the descending node is fixed.

Since Vesta has no atmosphere, the atmospheric drag effect is ignored. For Vesta, the rotation period is 5.3421 h, and a 1-degree drift of the right ascension of ascending node results in a 53.421 s drift of local time at the descending node. Then, the drift of local time at the descending node caused by solar gravitation is:

$$\begin{array}{cc}\Delta T& =\frac{53.421\times 180}{\pi }\Delta \Omega (t)\hfill \\ & =\frac{53.421\times 180}{\pi }\left[\Delta \Omega (t_0)+\frac{\partial \dot{\Omega }}{\partial a}\Delta a(t_0)(t-t_0)+\frac{\partial \dot{\Omega }}{\partial i}\Delta i(t_0)(t-t_0)+\frac{1}{2}\frac{\partial \dot{\Omega }}{\partial i}\frac{\mathrm{d}i}{\mathrm{d}t}{(t-t_0)}^2\right]\hfill \end{array}$$

(13)

where $\Delta a(t_0)$ and $\Delta i(t_0)$ are initial launch biases of a and i, respectively. Local time drift at the descending node comes from four variables in Equation (13): $\Delta \Omega (t_0)$, $\Delta a(t_0)$, $\Delta i(t_0)$, and $\mathrm{d}i/\mathrm{d}t$. Let $\Delta \Omega (t_0)$ be 0, and consider the remaining three items according to the following analysis. According to Equations (A.4) and (A.5), $\partial \dot{\Omega }/\partial a$ is $-3.2463\times {10}^{-10}$ for Vesta; this is about four orders of magnitude lower than $\partial \dot{\Omega }/\partial i$. Hence, $\partial \dot{\Omega }/\partial a\cdot \Delta a(t_0)(t-t_0)$ can be ignored.

Given an initial orbital element of circular Sun-synchronous orbits: a = 508.27 km, e = 0.0001, i = 90.2990 deg, $\Omega =-45$ deg, $\omega =30$ deg, and M = 0 deg, and the ecliptic longitude of the Sun is fixed at ${\beta }_s=90$ deg. We chose a small initial orbital deviation to investigate the local time drift at the descending node during a 2-Earth-year period. For an initial orbital deviation of $\Delta a_0=0.5$ km, $\Delta i_0=0$ deg, and $\Omega \in [0,360]$ deg, Figure 2a shows that the local time drift, $\Delta T$, is positive for $-({\beta }_s-\Omega )\in [0,90]\cup [180,270]$ deg, which means the local time at the descending node will be delayed. In contrast, the local time at the descending node will be brought forward for $-({\beta }_s-\Omega )\in [-90,0]\cup [90,180]$ deg. Therefore, we set $-({\beta }_s-\Omega )$ at −45 deg, where the local time drift will obtain the extremum value at the end of 2 Earth years.

There are two kinds of methods to fix the local time at the descending node. One is the preset inclination bias method, and the other is the multiple inclination bias method [34]. In the case of a preset bias, an initial inclination will be set in place of any other inclination controls during the spacecraft’s lifetime. This is simple and efficient for short-period missions, and the spacecraft will maintain the character of the Sun-synchronous orbits.

When $t_m=-\Delta i(t_0)/\dot{i}$, the local time drift in Equation (13) obtains the extremum:

$$\Delta {\widehat{T}}_m=\Delta \widehat{T}(t_m)=-\frac{24039}{5\pi }\frac{\partial \dot{\Omega }}{\partial i}\frac{\Delta i{(t_0)}^2}{\mathrm{d}i/\mathrm{d}t}$$

(14)

The local time drift at the descending node should be less than the extremum value during the spacecraft’s mission period: $|\Delta \widehat{T}(t_f)|\le \Delta \widehat{T}(t_m)$, where $t_f$ is the end time of the spacecraft’s mission period. Here, we set:

$$t_f=-k\frac{2\Delta i_0(t_0)}{\mathrm{d}i/\mathrm{d}t}$$

(15)

where $k>0$. Then we obtain $\Delta \widehat{T}(t_f)$ from Equation (13) for the following equation:

$$\Delta {\widehat{T}}_f=\frac{24039}{5\pi }(2k^2-2k)\frac{\partial \dot{\Omega }}{\partial i}\frac{\Delta i_0{(t_0)}^2}{\mathrm{d}i/\mathrm{d}t}$$

(16)

Let $\Delta {\widehat{T}}_f=-\Delta {\widehat{T}}_m$. Then, based on Equations (15) and (16), the preset inclination bias is:

$$\Delta i(t_0)=\frac{1-\sqrt{3}}{2}\frac{\mathrm{d}i}{\mathrm{d}t}{t}_f$$

(17)

Figure 3 shows that the preset inclination bias method limits the local time drift at the descending node in a 2-Earth-year mission period.

The multiple inclination bias method utilizes a periodic inclination bias to reduce the local time drift caused by solar gravitation. In a period of $t_c=-2\Delta i(t_{c_0})/(\mathrm{d}i/\mathrm{d}t)$, the extremum of local time drift is equal to a fixed local time drift-bound $\Delta \tilde{T}$:

$$|\Delta \tilde{T}(\frac{t_c}{2})|=|\Delta \tilde{T}|,$$

(18)

Due to the negligible $\partial \dot{\Omega }/\partial a$, using calculation in Equation (13), Equation (18) can be formulated as:

$$\frac{53.421\times 180}{\pi }\left[\Delta \Omega (t_0)-\frac{1}{2}\frac{\partial \dot{\Omega }}{\partial i}\frac{\Delta i{(t_{c_0})}^2}{\mathrm{d}i/\mathrm{d}t}\right]=|\Delta \tilde{T}|$$

(19)

where $t-t_0=t_c/2$. When $\Delta \Omega (t_0)=0$, the following was obtained:

$$-\frac{24039}{5\pi }\frac{\partial \dot{\Omega }}{\partial i}\frac{\Delta i{(t_{c_0})}^2}{\mathrm{d}i/\mathrm{d}t}=|\Delta \tilde{T}|$$

(20)

Then, the multiple inclination bias is as follows:

$$|\Delta i(t_{c_0})|=\sqrt{|\frac{5\pi }{24039}\frac{\mathrm{d}i/\mathrm{d}t\Delta \tilde{T}}{\partial \dot{\Omega }/\partial i}|}$$

(21)

where $\Delta i(t_{c_0})$ has a sign opposite to the sign of $\mathrm{d}i/\mathrm{d}t$. The multiple inclination bias is periodic during a long-period mission. Figure 4 shows the multiple inclination bias method when the limited local time drift bound is 0.2 min. We can see that this method is sensitive to the initial inclination deviation. From Figure 3 and Figure 4, we can see that it is very important to limiting the initial inclination deviation.

We applied the method in this section to the existence and control of Sun-synchronous orbits of Lutetia and Eros. Given the initial orbital elements of the circular Sun-synchronous orbit of Lutetia, a = 186.49 km, e = 0.00001, i = 172.52 deg, $\Omega =-45$ deg, $\omega =30$ deg, and M = 0 deg; and for the initial orbital elements of the Sun-synchronous orbit of Eros, a = 118.72 km, e = 0.60566, i = 165.100 deg, $\Omega =-45$ deg, $\omega =30$ deg, and M = 0 deg. Using the multiple inclination bias method, Figure 5 shows the local time drift change of the Sun-synchronous orbits of Lutetia and Eros, respectively. The local time drifts of the two orbits are not sensitive to the initial inclination deviation but are sensitive to the initial semi-major axis deviation. Compared with the Sun-synchronous orbits’ control period of Vesta, the control periods of the Sun-synchronous orbits of Lutetia and Eros are much longer.

3. Orbits with a Critical Inclination

The orbit with a critical inclination is described as critical as there is zero apogee drift for spacecraft in elliptical orbits at this inclination. This orbit is a practical option for long-term spacecraft flying around oblate planets such as Earth and Mars [16]. The secular perturbations of the first and second order are considered in the investigation of orbits at the critical inclination of Vesta. The mean variation of eccentricity from the secular part of the first and second order are

$$\dot{e}=0.$$

(22)

To meet the condition for orbit with critical inclination, the average variation of the argument of perigee in the second equation of (A.6) needs to be zero; then, we obtain a quadratic equation of variable ${\sin }^{2}i$

$$h({\sin }^{2}i)=A{\sin }^{4}i+B{\sin }^{2}i+C=0$$

(23)

where

$$A=\frac{9J_2^2{R}^4}{4p^4}n\left(\frac{15}{4}\sqrt{1-e^2}-\frac{15}{32}{e}^2+\frac{215}{48}\right)-\frac{35J_4{R}^4}{8p^4}n\left(\frac{81}{16}{e}^2+\frac{21}{4}\right),$$

(24)

$$B=\frac{9J_2^2{R}^4}{4p^4}n\left(-\frac{11}{2}\sqrt{1-e^2}-\frac{3}{8}{e}^2-\frac{103}{12}\right)+\frac{35J_4{R}^2}{8J_2{p}^2}n\left(\frac{9}{4}{e}^2+\frac{31}{14}\right)-\frac{15J_2{R}^2}{4p^2}n,$$

(25)

$$C=\frac{9J_2^2{R}^4}{4p^4}n\left[-\frac{5J_4}{6J_2^2}\left(\frac{9}{2}{e}^2+4\right)+2\sqrt{1-e^2}+\frac{7}{12}{e}^2+4\right]+\frac{3J_2{R}^2}{p^2}n.$$

(26)

Equation (23) has two different real roots of ${\sin }^{2}i$ if the discriminate of $h({\sin }^{2}i)$ is non-negative.

In Figure 6, the variation of critical inclinations as the semi-major axis a increases from 265 km to 663 km are presented. We note that the critical inclination decreases slowly and monotonically concerning the orbit’s altitude for $e\le 0.2$. For $e\ge 0.4$, the critical inclination increases rapidly and monotonically. For Earth, the critical inclinations are fixed when considering only J₂, while in the case of Vesta, critical inclinations depend on the semi-major axis and eccentricity.

We chose two different orbits of critical inclination with a = 492.171 km, e = 0.2, i = 63.499 deg, and a = 522.580 km, e = 0.1, i = 63.519 deg, to evaluate the evolution of $\omega$ and e over a period of 250T, where T is a Vesta day. The zonal harmonic coefficients up to J₆ are considered, and the results are shown in Figure 7. The drift amplitude of $\omega$ is about 2 deg for 250T, and the e is about 0.0023 for the first case in Figure 7. The drift amplitude of $\omega$ and e are slightly smaller in the second case since the altitude of the second orbit is higher near the stationary orbit of Vesta.

Some special properties were found when we applied the method in this section to Lutetia and Eros. For Lutetia, orbits with critical inclination do not exist with the condition of e < 0.6 and a > 153 km. For Eros, orbits with critical inclination do not exist with e < 0.9 and a > 91 km.

4. Repeating Ground-Track Orbits

Repeating ground-track orbits, which have regular ground-track patterns, are defined as orbits with the trajectory ground track repeating after a whole number of revolutions within some days. Repeating ground-track orbits represent an essential role in planet observation. We can monitor the area of Vesta by using spacecraft on these orbits and detect the surface and body of Vesta during a period. The interval of the adjacent ground track in the equator is

$$\Delta \lambda =T_N({\omega }_V-\dot{\Omega }),$$

(27)

where ${\omega }_V$ is the rotational angular velocity of Vesta, and T_N is the nodal period of the motion of the spacecraft, which is presented in Equation (A.11). The condition for achieving a repeating ground track can be written as

$$P\Delta \lambda =2\pi Q,$$

(28)

where P and Q are positive integers. Equation (28) means spacecraft will have completed P revolutions in Q days of Vesta. With the definition of $T_N$ in Equation (27), Equation (28) can be formulated as

$$Z=\frac{\dot{\overline{M}}+\dot{\overline{\omega }}}{{\omega }_V-\dot{\Omega }},$$

(29)

where Z = P/Q is the ground-track repetition parameter, widely used in practice. $\dot{\overline{M}}$ and $\dot{\overline{\omega }}$ are formulated in Equation (A.6). We can see that once the ground-track repetition parameter is given, Equation (29) represents the relation between a, e, and i. The repeating ground-track orbits are often Sun-synchronous, in which Equation (4) is also satisfied. If P and Q are fixed, then a and i can be calculated. Based on Equation (4), we can write Equation (29) as

$$\dot{\overline{M}}+\dot{\overline{\omega }}=Z({\omega }_V-n_s),$$

(30)

for a given ground-track repetition parameter, with Equations (A.1) and (A.6), Equation (30) can be further formulated as a quadratic equation of variable ${\sin }^{2}i$

$$g(\sin i)=\alpha {\sin }^{4}i+\beta {\sin }^{2}i+\gamma =0.$$

(31)

where

$$\begin{array}{ll}\alpha =& \frac{9J_2^2{R}^4}{4p^4}n[\frac{105J_4}{32J_2^2}{e}^2\sqrt{1-e^2}-\frac{35J_4}{18J_2^2}\left(\frac{81}{16}{e}^2+\frac{21}{4}\right)\hfill \\ & +\sqrt{1-e^2}\left(\frac{315e^4}{32(1-e^2)}+\frac{103e^2}{12}+\frac{413}{48}\right)-\frac{51}{32}{e}^2+\frac{269}{48}],\hfill \end{array}$$

(32)

$$\begin{array}{cc}\beta =& \frac{3J_2{R}^2}{4p^2}n\left\{\frac{3J_2{R}^2}{p^2}\left[\frac{35J_4}{6J_2^2}\left(\frac{9e^2}{4}+\frac{31}{14}\right)-\frac{3}{8}{e}^2-\frac{11}{2}\sqrt{1-e^2}-\frac{103}{12}\right]-5\right\}\hfill \\ & +\frac{9J_2^2{R}^4}{4p^4}n\sqrt{1-e^2}\left(\frac{25J_4}{4J_2^2}{e}^2-\frac{35e^4}{4(1-e^2)}-\frac{26}{3}{e}^2-\frac{3}{2}\sqrt{1-e^2}-\frac{19}{3}\right)\hfill \\ & -\frac{9J_2{R}^2}{4p^2}n\sqrt{1-e^2},\hfill \end{array}$$

(33)

$$\begin{array}{cc}\gamma =& \frac{3J_2{R}^2}{4p^2}n\left\{\frac{3J_2{R}^2}{p^2}\left[-\frac{5J_4}{6J_2^2}\left(\frac{9}{2}{e}^2+4\right)+\frac{7}{12}{e}^2+2\sqrt{1-e^2}+4\right]+4\right\}\hfill \\ & +\frac{9J_2{R}^4}{4p^4}n\sqrt{1-e^2}\left(-\frac{5J_4}{4J_2^2}{e}^2+\frac{35e^4}{12(1-e^2)}+\frac{10}{3}{e}^2+\frac{1}{2}\sqrt{1-e^2}+\frac{5}{2}\right)\hfill \\ & +\frac{3J_2{R}^2}{2p^2}n\sqrt{1-e^2}+n+Z(n_s-{\omega }_V).\hfill \end{array}$$

(34)

Equation (31) has four real roots of variable $\sin i$ only if the discriminant of $g(\sin i)=0$ is above zero.

$$\Delta ={\beta }^2-4\alpha [\gamma +Z(n_s-{\omega }_V)]\ge 0.$$

(35)

The rotational angular velocity of Vesta is greater than the mean motion of Vesta around the Sun. Equation (35) indicates that the lower bound of Z is

$$Z\ge \frac{\gamma -{\beta }^2/4\alpha }{{\omega }_V-n_s}.$$

(36)

Vesta rotates on its axis once every 5 h 20 m 30 s, which is denoted as $T_V$; the orbit period of spacecraft around Vesta $T_s$ is 2 h 20 m 48 s at a low altitude of 50 km. The upper bound of Z can be calculated by $T_V/T_s$, which is equal to 2.276. Concerning the semi-major axis for different Z, the change of inclinations of Sun-synchronous repeating ground-track orbits is presented in Figure 8. We can see that for a fixed Z, the inclination decreases with respect to the semi-major axis monotonously, and it is sensitive to the change of the semi-major axis. For Z = 1.1, the maximum meaningful semi-major axis is 540 km near the stationary orbit of Vesta. Compared with major planets in the solar system, such as Earth, Mars, Jupiter, and Saturn, Vesta’s ground-track repetition parameter is small [16,19,35]. For Mars, the repetition parameter is roughly between 12 and 14 [16]. For Jupiter, it is roughly between 2.5 and 3.5 [19]. Larger repetition parameters can offer a larger scale selection for Sun-synchronous repeating ground-track orbits. Figure 9 shows Vesta’s Sun-synchronous repeating ground-track orbit for a given repetition parameter.

We also calculated the discriminant in Equation (35) for Lutetia and Eros. We found that Equation (31) has no real root for Lutetia and Eros in the case of e = 0, which means that there are no Sun-synchronous repeating ground-track circular orbits for Lutetia and Eros. On the other hand, for the case of e > 0, there are a few Sun-synchronous repeating ground-track orbits for Lutetia and Eros when e > 0.9.

5. Stationary Orbits

A stationary orbit, also called a synchronous equatorial orbit, is a type of synchronous orbit around a planet or moon where the orbit is directly over the equator. We call it a synchronous orbit (SO) for Vesta, known as geosynchronous equatorial orbit (GEO) for Earth. Spacecraft in the stationary orbit remain orbiting over the same spot on the planet’s surface and would appear to be standing still above the surface. The most preferred use for spacecraft in this orbit is for communication, navigation, and space science. The spacecraft in the stationary orbit must satisfy very precise conditions to a fixed position related to the orbiting planet, and its orbit positions are given in the form of a station longitude. We will investigate the existence of the stationary orbits of Vesta in this subsection.

Based on Equations (A.12) and (A.13), spacecraft in the stationary orbits of Vesta must satisfy the following condition [19]

$$\{\begin{array}{c}\dot{r}=\ddot{r}=0,\\ \dot{\lambda }=n,\ddot{\lambda }=0,\\ \dot{\phi }=\ddot{\phi }=0.\end{array}$$

(37)

When we set $\phi =0$ and $\dot{\phi }=0$, the first Equation of (A.13) can be written as

$$\ddot{r}-r{\dot{\lambda }}^2+\frac{\mu }{r^2}+\frac{3\mu J_2{R}^2}{2r^4}-\frac{15\mu J_4{R}^4}{8r^6}=0$$

(38)

Solving Equation (38) with the condition in Equation (37), we obtain the radius of the stationary orbit of Vesta, which is 549.74 km.

We continue to analyze the perturbation of the SO and slightly inclined synchronous orbit (ISO) caused by solar radiation pressure, leading to drift of eccentricity. The solar radiation pressure on spacecraft in SO or slightly ISO is

$$F_s=-Kp(\frac{A}{m})S$$

(39)

where K is the reflection parameter of the spacecraft, p is the intensity of solar radiation pressure, A is the area of the spacecraft vertical to the Sun, m is the spacecraft’s mass, and S is the unit vector from the spacecraft to the Sun. We use S in this section as an approximation of the vector from Vesta to the Sun. S can be expressed as

$$S={\left(\begin{array}{ccc}\cos l_s& \sin l_s\cos i_s& \sin l_s\sin i_s\end{array}\right)}^{\prime }$$

(40)

where $l_s$ is the angle between solar radiation direction and X axis and $i_s$ is the obliquity of the ecliptic plane of Vesta, $l=\omega +M+\Omega$. The average distance from Vesta to the Sun is 353.4 million km, and the mean solar irradiance around Vesta is 247.41 W/m² [36]. Then, we obtained the solar radiation pressure of $F_s=8.2527\times {10}^{-9}S$ for a spacecraft with the following parameters: reflection parameter K = 1, average area A = 10 m², and mass m = 1000 kg.

In the Radial–Tangential–Normal coordinate system, the components of $F_s$ are:

$$\{\begin{array}{c}{F}_r=F_s\cdot {\mathrm{r}}_0,\\ F_t=F_s\cdot {\mathrm{t}}_0,\\ F_n=F_s\cdot {\mathrm{n}}_0.\end{array}$$

(41)

where $\left(\begin{array}{ccc}{\mathrm{r}}_0& {\mathrm{t}}_0& {\mathrm{n}}_0\end{array}\right)$ is the unit vector in the kinetic coordinate system [34]. The average derivatives of the eccentricity vector $e_x=e\cos (\omega +\Omega )$, $e_y=e\sin (\omega +\Omega )$ are as follows:

$$\begin{array}{l}\frac{\mathrm{d}{e}_x}{\mathrm{d}t}=\frac{1}{na}\left(F_r\sin l+2F_t\cos l\right),\\ \frac{\mathrm{d}{e}_y}{\mathrm{d}t}=\frac{1}{na}\left(-F_r\cos l+2F_t\sin l\right).\end{array}$$

(42)

We further formulated Equation (42) as [34]:

$$\begin{array}{l}\frac{\mathrm{d}{e}_x}{\mathrm{d}t}=\frac{F_s}{2na}(3\sin l_s\cos i_s-\cos l_s\sin 2l+\sin l_s\cos i_s\cos 2l),\\ \frac{\mathrm{d}{e}_y}{\mathrm{d}t}=\frac{F_s}{2na}(-3\cos l_s+\sin l_s\cos i_s\sin 2l+\cos l_s\cos 2l),\end{array}$$

(43)

When considering the solar radiation pressure perturbation over a long period of time, we obtained the following expressions of the eccentricity vector:

$$\begin{array}{l}{e}_x(t)-e_x(t_0)=\frac{3F_s}{2na}\frac{1}{n_s}[\cos l_s(t)-\cos l_s(t_0)]\cos i_s,\\ e_y(t)-e_y(t_0)=\frac{3F_s}{2na}\frac{1}{n_s}[\sin l_s(t)-\sin l_s(t_0)].\end{array}$$

(44)

Then, the eccentricity vector can be further formulated as:

$$\begin{array}{l}{\left\{e_x(t)-\left[e_x(t_0)-\frac{3F_s}{2na}\frac{1}{n_s}\cos i_s\cos l_s(t_0)\right]\right\}}^2/{\cos }^2{i}_s+\\ {\left\{e_y(t)-\left[e_y(t_0)-\frac{3F_s}{2na}\frac{1}{n_s}\sin l_s(t_0)\right]\right\}}^2={\left(\frac{3F_s}{2na}\frac{1}{n_s}\right)}^2.\end{array}$$

(45)

We observed from Equation (45) that the curve of the eccentricity vector is an ellipse caused by the solar radiation pressure. The eccentricity vector moves along the curve of the ellipse during Vesta’s orbit around the Sun. The semi-major axis of the ellipse is $1.3\times {10}^{-3}$, and the semi-minor axis of the ellipse is $1.2\times {10}^{-3}$.

We then analyzed the inclination change caused by solar radiation pressure. Based on Equation (41) [34], we can write the normal perturbation force caused by solar radiation pressure as:

$$F_n=F_s\left(\sin l_s\sin i_s+\cos l_s\cdot i_x-\sin l_s\cos i_s\cdot i_y\right)$$

(46)

where $i_x=\sin i\sin \Omega ,$$i_y=\sin i\cos \Omega$. The derivatives of the inclination vector are:

$$\begin{array}{l}\frac{\mathrm{d}{i}_x}{\mathrm{d}t}=\frac{F_s}{na}\left(\sin l_s\sin i_s+\cos l_s\cdot i_x-\sin l_s\cos i_s\cdot i_y\right)\sin l,\\ \frac{\mathrm{d}{i}_y}{\mathrm{d}t}=\frac{F_s}{na}\left(\sin l_s\sin i_s+\cos l_s\cdot i_x-\sin l_s\cos i_s\cdot i_y\right)\cos l.\end{array}$$

(47)

For a spacecraft in SO or slightly ISO, $i_x\ll 1$, $i_y\ll 1$. $l_s$ should be constant during one Vesta rotation period; thus, the time integral of the inclination vector in a Vesta rotation period for time t is:

$$\begin{array}{l}{i}_x(t)-i_x(t_0)=-\frac{F_s}{na}\sin l_s\sin i_s[\cos l(t)-\cos l(t_0)],\\ i_y(t)-i_y(t_0)=\frac{F_s}{na}\sin l_s\sin i_s[\sin l(t)-\sin l(t_0)].\end{array}$$

(48)

Equation (48) indicates:

$$\begin{array}{c}{\left\{i_x(t)-\left[i_x(t_0)+\frac{F_s}{na}\sin l_s\sin i_s\cos l(t_0)\right]\right\}}^2+\\ {\left\{i_y(t)-\left[i_y(t_0)-\frac{F_s}{na}\sin l_s\sin i_s\sin l(t_0)\right]\right\}}^2={\left(\frac{F_s}{na}\sin l_s\sin i_s\right)}^2.\end{array}$$

(49)

We learn from Equation (49) that the curve of the inclination vector is a circle caused by the solar radiation pressure; the inclination vector moves on the circle with the period of one Vesta day, and the circle’s radius is less than $1.386\times {10}^{-9}$ deg.

After analyzing the eccentricity drift and inclination drift caused by solar radiation pressure perturbation, we calculate the eccentricity drift and inclination drift caused by solar gravitation perturbation. The average time derivatives of the eccentricity vector $(e_x,e_y)$ are:

$$\begin{array}{ll}\frac{\mathrm{d}{e}_x}{\mathrm{d}t}& =\frac{3r{n}_s^2}{na}[\frac{1}{2}\left({\cos }^2{\beta }_s+{\sin }^2{\beta }_s{\cos }^2{i}_s\right)\sin l+\frac{1}{4}\left(\sin 3l-\sin l\right)\left({\cos }^2{\beta }_s-{\sin }^2{\beta }_s{\cos }^2{i}_s\right)\hfill \\ & +\frac{1}{4}\sin 2{\beta }_s\cos i_s\left(\cos l-\cos 3l\right)+\frac{1}{2}\left(\sin 3l+\sin l\right)\left({\sin }^2{\beta }_s{\cos }^2{i}_s-{\cos }^2{\beta }_s\right)\hfill \\ & +\sin {\beta }_s\cos {\beta }_s\cos i_s\left(\cos l+\cos 3l\right)]-\frac{r{n}_s^2}{na}\sin l,\hfill \\ \frac{\mathrm{d}{e}_y}{\mathrm{d}t}& =\frac{1}{na}[-3r{n}_s^2{\cos }^2\xi \cos l+r{n}_s^2\cos l\hfill \\ & +6r{n}_s^2\cos \xi \sin l\left(-\cos {\beta }_s\sin l+\sin {\beta }_s\cos i_s\cos l\right)].\hfill \end{array}$$

(50)

where $\xi$ is the angle between two unit vectors; one is the unit vector from the center of Vesta to the Sun, and the other is the unit vector from the center of Vesta to the spacecraft in SO, $l=\omega +M+\Omega$. Due to the existence of the trigonometric function of l in the first equation of (50), the drift of $e_x$ caused by solar gravitation is periodic. In the case of the second equation of (50), there is no secular term in the right hand of the equation, so the drift of $e_y$ caused by solar gravitation is also periodic.

Before analyzing the inclination drift caused by solar gravitation perturbation, we obtained the normal perturbation force caused by solar gravitation, which can be written as [37]:

$$F_n=3r{n}_s^2\left(\cos {\beta }_s\cos l+\sin {\beta }_s\cos i_s\sin l\right)\left(\cos {\beta }_s\cdot i_x-\sin {\beta }_s\cos i_s\cdot i_y+\sin {\beta }_s\sin i_s\right).$$

(51)

The derivatives of the inclination vector perturbed by the solar gravitation are:

$$\begin{array}{cc}\frac{\mathrm{d}{i}_x}{\mathrm{d}t}& =\frac{3r{n}_s^2}{2na}\left[\cos {\beta }_s\sin 2l+\sin {\beta }_s\cos i_s\left(1-\cos 2l\right)\right](\cos {\beta }_s\cdot i_x\hfill \\ & -\sin {\beta }_s\cos i_s\cdot i_y+\sin {\beta }_s\sin i_s),\hfill \\ \frac{\mathrm{d}{i}_y}{\mathrm{d}t}& =\frac{3r{n}_s^2}{2na}\left[\cos {\beta }_s\left(1+\cos 2l\right)+\sin {\beta }_s\cos i_s\sin 2l\right](\cos {\beta }_s\cdot i_x\hfill \\ & -\sin {\beta }_s\cos i_s\cdot i_y+\sin {\beta }_s\sin i_s).\hfill \end{array}$$

(52)

The short periodic terms in Equation (52) were neglected for a long period, and the derivatives of the inclination vector were further formulated as:

$$\begin{array}{l}\frac{\mathrm{d}{i}_x}{\mathrm{d}t}=\frac{3r{n}_s^2}{2na}\sin {\beta }_s\cos i_s\left(\cos {\beta }_s\cdot i_x-\sin {\beta }_s\cos i_s\cdot i_y+\sin {\beta }_s\sin i_s\right),\\ \frac{\mathrm{d}{i}_y}{\mathrm{d}t}=\frac{3r{n}_s^2}{2na}\cos {\beta }_s\left(\cos {\beta }_s\cdot i_x-\sin {\beta }_s\cos i_s\cdot i_y+\sin {\beta }_s\sin i_s\right).\end{array}$$

(53)

For a spacecraft in SO or slightly ISO, $i_x\ll 1$, $i_y\ll 1$. The time integral of the inclination vector in a period of Vesta orbiting around the Sun for time t is:

$$\begin{array}{l}{i}_x(t)-i_x(t_0)=\frac{3}{8}\frac{n_s^2}{n}\sin 2i_s\left\{\left(t-t_0\right)-\frac{1}{2}\left[\sin 2{\beta }_s(t)-\sin 2{\beta }_s(t_0)\right]\right\},\\ i_y(t)-i_y(t_0)=-\frac{3}{8}\frac{n_s^2}{n}\left[\cos 2{\beta }_s(t)-\cos 2{\beta }_s(t_0)\right]\sin i_s.\end{array}$$

(54)

We learned from Equation (54) that the drift of $i_y$ caused by solar gravitation is periodic, with a period of half of a Vesta year. If only the secular term of $\mathrm{d}{i}_x/\mathrm{d}t$ is considered, the drift of inclination throughout T is:

$$\Delta i_x=\frac{3}{8}\frac{n_s^2}{n}\sin 2i_{s}T\frac{180}{\pi }.$$

(55)

We obtained $\Delta i_x=0.0119$ deg over one Vesta year from Equation (55). For Earth, we learned that $\Delta i_x$ is 0.2712 deg over one Earth-year period.

6. Conclusions

We analyzed the existence of some types of orbits around Vesta, including the Sun-synchronous orbits, orbits with critical inclination, repeating ground-track orbits, and stationary orbits. Zonal harmonic coefficients J₂, J₃, and J₄ were considered in the gravitation model of these typical orbits. For the evolutionary evaluation of an orbit with a critical inclination, coefficients up to J₆ were considered.

First, we investigated the existence of Sun-synchronous orbits of Vesta, found the variations of inclination for different semi-major axes and eccentricity values, and calculated the local time drift at the descending node caused by solar gravitation. Then, two control methods were given to dampen the local time drift and applied to Lutetia and Eros. As a result, the control periods of Sun-synchronous orbits of Lutetia and Eros were found to be much longer.

Second, we analyzed the orbits with critical inclinations. In the case of Vesta, critical inclinations depend on the semi-major axis and eccentricity. When the zonal harmonic coefficients up to J₆ were considered, the eccentricity and the argument of perigee changed monotonically. Numerical calculation results showed that orbits with critical inclinations at a high altitude do not exist for Lutetia and Eros.

Third, we calculated different repetition parameters of repeating ground-track orbits of Vesta. Compared with major planets in the solar system, Vesta’s ground-track repetition parameters are small. We found that Sun-synchronous repeating ground-track circular orbits do not exist for Lutetia and Eros.

Finally, we calculated the radius of stationary orbits of Vesta in the spherical coordinates of an inertial frame. Then, we analyze the inclination and eccentricity drift caused by solar radiation pressure and solar gravitation, respectively. We found that, in the case of solar radiation pressure perturbation, the eccentricity vector moves on an ellipsis and the inclination vector moves on a circle. For the case of solar gravitation, the drift of the eccentricity vector is periodic, and the drift of inclination in the direction of the X axis is 0.0119 deg over a Vesta year.