The Orbital Dynamics of Synchronous Satellites: Irregular Motions in the 2 : 1 Resonance

Abstract

The orbital dynamics of synchronous satellites is studied. The 2 : 1 resonance is considered; in other words, the satellite completes two revolutions while the Earth completes one. In the development of the geopotential, the zonal harmonics 𝐽₂₀ and 𝐽₄₀ and the tesseral harmonics 𝐽₂₂ and 𝐽₄₂ are considered. The order of the dynamical system is reduced through successive Mathieu transformations, and the final system is solved by numerical integration. The Lyapunov exponents are used as tool to analyze the chaotic orbits.

1. Introduction

Synchronous satellites in circular or elliptical orbits have been extensively used for navigation, communication, and military missions. This fact justifies the great attention that has been given in literature to the study of resonant orbits characterizing the dynamics of these satellites since the 60s [1–14]. For example, Molniya series satellites used by the old Soviet Union for communication form a constellation of satellites, launched since 1965, which have highly eccentric orbits with periods of 12 hours. Another example of missions that use eccentric, inclined, and synchronous orbits includes satellites to investigate the solar magnetosphere, launched in the 90s [15].

The dynamics of synchronous satellites are very complex. The tesseral harmonics of the geopotential produce multiple resonances which interact resulting significantly in nonlinear motions, when compared to nonresonant orbits. It has been found that the orbital elements show relatively large oscillation amplitudes differing from neighboring trajectories [11].

Due to the perturbations of Earth gravitational potential, the frequencies of the longitude of ascending node Ω and of the argument of pericentre 𝜔 can make the presence of small divisors, arising in the integration of equation of motion, more pronounced. This phenomenon depends also on the eccentricity and inclination of the orbit plane. The importance of the node and the pericentre frequencies is smaller when compared to the mean anomaly and Greenwich sidereal time. However, they also have their contribution in the resonance effect. The coefficients 𝑙, 𝑚, 𝑝 which define the argument 𝜙_{𝑙𝑚𝑝𝑞} in the development of the geopotential can vary, producing different frequencies within the resonant cosines for the same resonance. These frequencies are slightly different, with small variations around the considered commensurability.

In this paper, the 2 : 1 resonance is considered; in other words, the satellite completes two revolutions while the Earth carries one. In the development of the geopotential, the zonal harmonics 𝐽₂₀ and 𝐽₄₀ and the tesseral harmonics 𝐽₂₂ and 𝐽₄₂ are considered. The order of the dynamical system is reduced through successive Mathieu transformations, and the final system is solved by numerical integration. In the reduced dynamical model, three critical angles, associated to the tesseral harmonics 𝐽₂₂ and 𝐽₄₂, are studied together. Numerical results show the time behavior of the semimajor axis, argument of pericentre and of the eccentricity. The Lyapunov exponents are used as tool to analyze the chaotic orbits.

2. Resonant Hamiltonian and Equations of Motion

In this section, a Hamiltonian describing the resonant problem is derived through successive Mathieu transformations.

Consider (2.1) to the Earth gravitational potential written in classical orbital elements [16, 17]𝑉=𝜇2𝑎+^∞_𝑙𝑙=2_𝑙𝑚=0_{𝑝=0−∞𝑞=+∞}𝜇𝑎𝑎𝑎^𝑒^𝑙𝐽_𝑙𝑚𝐹_𝑙𝑚(𝐼)G_𝑙𝑝𝑞(𝑒)cos𝜙_{𝑙𝑚𝑝𝑞}(𝑀,𝜔,Ω,𝜃),(2.1) where 𝜇 is the Earth gravitational parameter, 𝜇=3.986009×10¹⁴ m³/s², 𝑎, 𝑒, 𝐼, Ω, 𝜔, 𝑀 are the classical keplerian elements: 𝑎 is the semimajor axis, 𝑒 is the eccentricity, 𝐼 is the inclination of the orbit plane with the equator, Ω is the longitude of the ascending node, 𝜔 is the argument of pericentre, and 𝑀 is the mean anomaly, respectively; 𝑎_𝑒 is the Earth mean equatorial radius, 𝑎_𝑒=6378.140 km, 𝐽_𝑙𝑚 is the spherical harmonic coefficient of degree 𝑙 and order 𝑚, 𝐹_𝑙𝑚𝑝(𝐼) and G_𝑙𝑝𝑞(𝑒) are Kaula’s inclination and eccentricity functions, respectively. The argument 𝜙_{𝑙𝑚𝑝𝑞}(𝑀,𝜔,Ω,𝜃) is defined by𝜙_{𝑙𝑚𝑝𝑞}(𝑀,𝜔,Ω,𝜃)=𝑞𝑀+(𝑙−2𝑝)𝜔+𝑚Ω−𝜃−𝜆_𝑙𝑚+(𝑙−𝑚)𝜋2,(2.2) where 𝜃 is the Greenwich sidereal time, 𝜃=𝜔_𝑒𝑡 (𝜔_𝑒 is the Earth’s angular velocity, and 𝑡 is the time), and 𝜆_𝑙𝑚 is the corresponding reference longitude along the equator.

In order to describe the problem in Hamiltonian form, Delaunay canonical variables are introduced,𝐿=√𝜇𝑎,𝐺=ℓ=𝜇𝑀𝑎,1−𝑔𝑒²=,𝜔,𝐻ℎ==Ω𝜇.𝑎1−𝑒²cos(𝐼),(2.3)𝐿, 𝐺, and 𝐻 represent the generalized coordinates, and ℓ, 𝑔, and ℎ represent the conjugate momenta.

Using the canonical variables, one gets the Hamiltonian 𝐹,𝐹=𝜇²2𝐿²+^∞_𝑙𝑙=2_𝑚=0𝑅_𝑙𝑚,(2.4)

with the disturbing potential 𝑅_𝑙𝑚 given by𝑅_𝑙𝑚=^𝑙_{𝑝=0+∞𝑞=−∞}𝐵_{𝑙𝑚𝑝𝑞}(𝐿,𝐺,𝐻)cos𝜙_{𝑙𝑚𝑝𝑞}(ℓ,𝑔,ℎ,𝜃).(2.5)

The argument 𝜙_{𝑙𝑚𝑝𝑞} is defined by𝜙_{𝑙𝑚𝑝𝑞}(ℓ,𝑔,ℎ,𝜃)=𝑞ℓ+(𝑙−2𝑝)𝑔+𝑚ℎ−𝜃−𝜆_𝑙𝑚+(𝑙−𝑚)𝜋2,(2.6) and the coefficient 𝐵_{𝑙𝑚𝑝𝑞}(𝐿,𝐺,𝐻) is defined by𝐵_{𝑙𝑚𝑝𝑞}=^∞_𝑙𝑙=2_𝑙𝑚=0_{𝑝=0−∞𝑞=+∞}𝜇²𝐿²𝜇𝐿𝑎^𝑒2^𝑙𝐽_𝑙𝑚𝐹_𝑙𝑚𝑝(𝐿,𝐺,𝐻)𝐺_𝑙𝑝𝑞(𝐿,𝐺).(2.7)

The Hamiltonian 𝐹 depends explicitly on the time through the Greenwich sidereal time 𝜃. A new term 𝜔_𝑒Θ is introduced in order to extend the phase space. In the extended phase space, the extended Hamiltonian 𝐻 is given by𝐻=𝐹−𝜔_𝑒Θ.(2.8)

For resonant orbits, it is convenient to use a new set of canonical variables. Consider the canonical transformation of variables defined by the following relations:𝑋𝑥==ℓ𝐿+,𝑔+𝑌ℎ=,𝐺−𝑦𝐿=,𝑔+𝑍ℎ,=𝐻𝑧−=𝐺ℎ,,Θ𝜃==Θ𝜃,,(2.9)

where 𝑋,𝑌,𝑍,Θ,𝑥,𝑦,𝑧,𝜃 are the modified Delaunay variables.

The new Hamiltonian 𝐻^, resulting from the canonical transformation defined by (2.9), is given by𝐻^=𝜇²2𝑋²−𝜔_𝑒Θ+^∞_𝑙𝑙=2_𝑚=0𝑅^𝑙𝑚,(2.10)

where the disturbing potential 𝑅^𝑙𝑚 is given by𝑅^𝑙𝑚=^𝑙_{𝑝=0+∞𝑞=−∞}𝐵^{𝑙𝑚𝑝𝑞}(𝑋,𝑌,𝑍)cos𝜙_{𝑙𝑚𝑝𝑞}(𝑥,𝑦,𝑧,𝜃).(2.11)

Now, consider the commensurability between the Earth rotation angular velocity 𝜔_𝑒 and the mean motion 𝑛=𝜇²/𝑋³. This commensurability can be expressed as𝑞𝑛−𝑚𝜔_𝑒≅0,(2.12) considering 𝑞 and 𝑚 as integers. The ratio 𝑞/𝑚 defining the commensurability will be denoted by 𝛼. When the commensurability occurs, small divisors arise in the integration of the equations of motion [9]. These periodic terms in the Hamiltonian 𝐻^ with frequencies 𝑞𝑛−𝑚𝜔_𝑒 are called resonant terms. The other periodic terms are called short- and long-period terms.

The short- and long-period terms can be eliminated from the Hamiltonian 𝐻^ by applying an averaging procedure [12, 18]:𝐻^=41𝜋²^02𝜋^02𝜋𝐻^𝑑𝜉_sp𝑑𝜉_lp.(2.13) The variables 𝜉_sp and 𝜉_lp represent the short- and long-period terms, respectively, to be eliminated of the Hamiltonian 𝐻^.

The long-period terms have a combination in the argument 𝜙_{𝑙𝑚𝑝𝑞} which involves only the argument of the pericentre 𝜔 and the longitude of the ascending node Ω. From (2.10) and (2.11), these terms are represented by the new variables in the following equation:𝐻^_lp=^∞_𝑙𝑙=2_𝑙𝑚=0_{𝑝=0+∞𝑞=−∞}𝐵^{𝑙𝑚𝑝𝑞}(𝑋,𝑌,𝑍)cos((𝑙−2𝑝)(𝑦−𝑧)+𝑚𝑧).(2.14)

The short-period terms are identified by the presence of the sidereal time 𝜃 and mean anomaly 𝑀 in the argument 𝜙_{𝑙𝑚𝑝𝑞}; in this way, from (2.10) and (2.11), the term 𝐻^_sp in the new variables is given by the following equations:𝐻^_sp=^∞_𝑙𝑙=2_𝑙𝑚=0_{𝑝=0+∞𝑞=−∞}𝐵^{𝑙𝑚𝑝𝑞}(𝑋,𝑌,𝑍)cos𝑞(𝑥−𝑦)−𝑚𝜃+𝜁_𝑝.(2.15)

The term 𝜁_𝑝 represents the other variables in the argument 𝜙_{𝑙𝑚𝑝𝑞}, including the argument of the pericentre 𝜔 and the longitude of the ascending node Ω, or, in terms of the new variables, 𝑦−𝑧 and 𝑧, respectively.

A reduced Hamiltonian 𝐻_𝑟 is obtained from the Hamiltonian 𝐻′ when only secular and resonant terms are considered. The reduced Hamiltonian 𝐻_𝑟 is given by𝐻_𝑟=𝜇²2𝑋²−𝜔_𝑒Θ+^∞_𝑗=1+𝐵^{2𝑗,0,𝑗,0}(𝑋,𝑌,𝑍)_∞_𝑙𝑙=2_𝑙𝑚=2_𝑝=0𝐵^{𝑙𝑚𝑝(𝛼𝑚)}(𝑋,𝑌,𝑍)cos𝜙_{𝑙𝑚𝑝(𝛼𝑚)}(𝑥,𝑦,𝑧,𝜃).(2.16)

Several authors, [11, 15, 19–22], also use this simplified Hamiltonian to study the resonance.

The dynamical system generated from the reduced Hamiltonian, (2.16), is given by𝑑(𝑋,𝑌,𝑍,Θ)𝑑𝑡=𝜕(𝑥𝜕,𝐻𝑦,𝑧^𝑟,𝜃),𝑑(𝑥,𝑑𝑦𝑡,𝑧,𝜃)=−𝜕(𝑋,𝜕𝑌𝐻,𝑍^𝑟,Θ).(2.17)

The equations of motion 𝑑𝑋/𝑑𝑡, 𝑑𝑌/𝑑𝑡, and 𝑑𝑍/𝑑𝑡 defined by (2.17) are𝑑𝑋𝑑𝑡=−𝛼^∞_𝑙𝑙=2_𝑙𝑚=2_𝑝=0𝑑𝑌𝑚𝐵^{𝑙𝑚𝑝(𝛼𝑚)}(𝑋,𝑌,𝑍)sin𝜙_{𝑙𝑚𝑝(𝛼𝑚)}(𝑥,𝑦,𝑧,𝜃),(2.18)𝑑𝑡=−^∞_𝑙𝑙=2_𝑙𝑚=2_𝑝=0𝑑𝑍(𝑙−2𝑝−𝑚𝛼)𝐵^{𝑙𝑚𝑝(𝛼𝑚)}(𝑋,𝑌,𝑍)sin𝜙_{𝑙𝑚𝑝(𝛼𝑚)}(𝑥,𝑦,𝑧,𝜃),(2.19)𝑑𝑡=^∞_𝑙𝑙=2_𝑙𝑚=2_𝑝=0(𝑙−2𝑝−𝑚)𝐵^{𝑙𝑚𝑝(𝛼𝑚)}(𝑋,𝑌,𝑍)sin𝜙_{𝑙𝑚𝑝(𝛼𝑚)}(𝑥,𝑦,𝑧,𝜃).(2.20)

From (2.18) to (2.20), one can determine the first integral of the system determined by the Hamiltonian 𝐻_𝑟.

Equation (2.18) can be rewritten as1𝛼𝑑𝑑𝑋𝑡=−^∞_𝑙𝑙=2_𝑙𝑚=2_𝑝=0𝑚𝐵^{𝑙𝑚𝑝(𝛼𝑚)}(𝑋,𝑌,𝑍)sin𝜙_{𝑙𝑚𝑝(𝛼𝑚)}(𝑥,𝑦,𝑧,𝜃).(2.21)

Adding (2.19) and (2.20),𝑑𝑌𝑑𝑡+𝑑𝑑𝑍𝑡=(𝛼−1)^∞_𝑙𝑙=2_𝑙𝑚=2_𝑝=0𝑚𝐵^{𝑙𝑚𝑝(𝛼𝑚)}(𝑋,𝑌,𝑍)sin𝜙_{𝑙𝑚𝑝(𝛼𝑚)}(𝑥,𝑦,𝑧,𝜃),(2.22)

and substituting (2.21) and (2.22), one obtains𝑑𝑌𝑑𝑡+𝑑𝑑𝑍𝑡=−(𝛼−1)𝛼1𝑑𝑑𝑋𝑡.(2.23)

Now, (2.23) is rewritten as1−𝛼1𝑑𝑑𝑋𝑡+𝑑𝑑𝑌𝑡+𝑑𝑑𝑍𝑡=0.(2.24)

In this way, the canonical system of differential equations governed by 𝐻_𝑟 has the first integral generated from (2.24):1−𝛼1𝑋+𝑌+𝑍=𝐶¹,(2.25)

where 𝐶₁ is an integration constant.

Using this first integral, a Mathieu transformation(𝑋,𝑌,𝑍,Θ,𝑥,𝑦,𝑧,𝜃)⟶𝑋₁,𝑌₁,𝑍₁,Θ₁,𝑥₁,𝑦₁,𝑧₁,𝜃₁(2.26)

can be defined.

This transformation is given by the following equations:𝑋₁=𝑋,𝑌₁=𝑌,𝑍₁=1−1𝛼𝑋+𝑌+𝑍,Θ¹𝑥=Θ,₁=𝑥−1−1𝛼𝑧,𝑦¹=𝑦−𝑧,𝑧₁=𝑧,𝜃₁=𝜃.(2.27) The subscript 1 denotes the new set of canonical variables. Note that 𝑍₁=𝐶₁, and the 𝑧₁ is an ignorable variable. So the order of the dynamical system is reduced in one degree of freedom.

Substituting the new set of canonical variables, 𝑋₁, 𝑌₁, 𝑍₁, Θ₁, 𝑥₁, 𝑦₁, 𝑧₁, 𝜃₁, in the reduced Hamiltonian given by (2.16), one gets the resonant Hamiltonian. The word “resonant” is used to denote the Hamiltonian 𝐻_rs which is valid for any resonance. The periodic terms in this Hamiltonian are resonant terms. The Hamiltonian 𝐻_rs is given by𝐻_rs=𝜇²2𝑋²¹−𝜔_𝑒Θ₁+^∞_𝑗=1+_∞𝐵_{2𝑗,0,𝑗,0}𝑋₁,𝑌₁,𝐶₁_𝑙𝑙=2_𝑙𝑚=2_𝑝=0𝐵_{𝑙𝑚𝑝,(𝛼𝑚)}𝑋₁,𝑌₁,𝐶₁cos𝜙_{𝑙𝑚𝑝(𝛼𝑚)}𝑥₁,𝑦₁,𝜃₁.(2.28)

The Hamiltonian 𝐻_rs has all resonant frequencies, relative to the commensurability 𝛼, where the 𝜙_{𝑙𝑚𝑝(𝛼𝑚)} argument is given by𝜙_{𝑙𝑚𝑝(𝛼𝑚)}=𝑚𝛼𝑥₁−𝜃₁+(𝑙−2𝑝−𝛼𝑚)𝑦₁−𝜙_{𝑙𝑚𝑝(𝛼𝑚)0},(2.29)

with𝜙_{𝑙𝑚𝑝(𝛼𝑚)0}=𝑚𝜆_𝑙𝑚−(𝑙−𝑚)𝜋2.(2.30)

The secular and resonant terms are given, respectively, by 𝐵_{2𝑗,0,𝑗,0}(𝑋₁,𝑌₁,𝐶₁) and 𝐵_{𝑙𝑚𝑝(𝛼𝑚)}(𝑋₁,𝑌₁,𝐶₁).

Each one of the frequencies contained in 𝑑𝑥₁/𝑑𝑡, 𝑑𝑦₁/𝑑𝑡, 𝑑𝜃₁/𝑑𝑡 is related, through the coefficients 𝑙, 𝑚, to a tesseral harmonic 𝐽_𝑙𝑚. By varying the coefficients 𝑙, 𝑚, 𝑝 and keeping 𝑞/𝑚 fixed, one finds all frequencies 𝑑𝜙_{1,𝑙𝑚𝑝(𝛼𝑚)}/𝑑𝑡 concerning a specific resonance.

From 𝐻_rs, taking, 𝑗=1,2, 𝑙=2,4, 𝑚=2, 𝛼=1/2, and 𝑝=0,1,2,3, one gets𝐻₁=𝜇²2𝑋²¹+𝐵−_1,2201𝜔_𝑒Θ₁+𝐵_1,2010𝑋₁,𝑌₁,𝐶₁+𝐵_1,4020𝑋₁,𝑌₁,𝐶₁+𝐵_1,2211𝑋₁,𝑌₁,𝐶₁cos𝑥₁−2𝜃₁+𝑦₁−2𝜆₂₂+𝐵_1,2221𝑋₁,𝑌₁,𝐶₁cos𝑥₁−2𝜃₁−𝑦₁−2𝜆₂₂+𝐵_1,4211𝑋₁,𝑌₁,𝐶₁cos𝑥₁−2𝜃₁−3𝑦₁−2𝜆₂₂+𝐵_1,4221𝑋₁,𝑌₁,𝐶₁cos𝑥₁−2𝜃₁+𝑦₁−2𝜆₄₂+𝜋+𝐵_1,4231𝑋₁,𝑌₁,𝐶₁cos𝑥₁−2𝜃₁−𝑦₁−2𝜆₄₂+𝜋𝑋₁,𝑌₁,𝐶₁cos𝑥₁−2𝜃₁−3𝑦₁−2𝜆₄₂+𝜋.(2.31)

The Hamiltonian 𝐻₁ is defined considering a fixed resonance and three different critical angles associated to the tesseral harmonic 𝐽₂₂; the critical angles associated to the tesseral harmonic 𝐽₄₂ have the same frequency of the critical angles associated to the 𝐽₂₂ with a difference in the phase. The other terms in 𝐻_rs are considered as short-period terms.

Table 1 shows the resonant coefficients used in the Hamiltonian 𝐻₁.

Table 1 

Resonant coefficients.

Finally, a last transformation of variables is done, with the purpose of writing the resonant angle explicitly. This transformation is defined by𝑋𝑥₄₄=𝑋₁,𝑌₄=𝑌₁,Θ₄=Θ₁+2𝑋₁,=𝑥₁−2𝜃₁,𝑦₄=𝑦₁,𝜃₄=𝜃₁.(2.32)

So, considering (2.31) and (2.32), the Hamiltonian 𝐻₄ is found to be𝐻₄=𝜇²2𝑋²⁴+𝐵−_4,2201𝜔_𝑒Θ₄−2𝑋₄+𝐵_4,2010𝑋₄,𝑌₄,𝐶₁+𝐵_4,4020𝑋₄,𝑌₄,𝐶₁+𝐵_4,2211𝑋₄,𝑌₄,𝐶₁cos𝑥₄+𝑦₄−2𝜆₂₂+𝐵_4,2221𝑋₄,𝑌₄,𝐶₁cos𝑥₄−𝑦₄−2𝜆₂₂+𝐵_4,4211𝑋₄,𝑌₄,𝐶₁cos𝑥₄−3𝑦₄−2𝜆₂₂+𝐵_4,4221𝑋₄,𝑌₄,𝐶₁cos𝑥₄+𝑦₄−2𝜆₄₂+𝜋+𝐵_4,4231𝑋₄,𝑌₄,𝐶₁cos𝑥₄−𝑦₄−2𝜆₄₂+𝜋𝑋₄,𝑌₄,𝐶₁cos𝑥₄−3𝑦₄−2𝜆₄₂+𝜋,(2.33) with 𝜔_𝑒Θ₄ constant and 𝐵_4,2010=𝜇⁴𝑋⁶⁴𝑎_𝑒2𝐽₂₀−43𝐶¹+2𝑋₄²𝑋₄+𝑌₄²+141+23−𝑌⁴²−2𝑋₄𝑌₄𝑋₄₂,(2.34)𝐵_4,4020=𝜇⁶𝑋⁴¹⁰𝑎_𝑒4𝐽₄₀⎛⎜⎜⎝105641−𝐶¹+2𝑋₄²𝑋₄+𝑌₄²²−32+185𝐶¹+2𝑋₄²𝑋₄+𝑌₄²×1+5−𝑌₄₂⎞⎟⎟⎠−2𝑋₄𝑌₄𝑋₄₂,(2.35)𝐵_4,2201=218𝑋⁷⁴𝜇⁴𝑎_𝑒2𝐽₂₂1+𝐶¹+2𝑋₄𝑋₄+𝑌₄²−𝑌₄₂−2𝑋₄𝑌₄,(2.36)𝐵_4,2211=32𝑋⁷⁴𝜇⁴𝑎_𝑒2𝐽₂₂32−23𝐶¹+2𝑋₄²𝑋₄+𝑌₄²−𝑌₄₂−2𝑋₄𝑌₄,(2.37)𝐵_4,2221=−38𝑋⁷⁴𝜇⁴𝑎_𝑒2𝐽₂₂1−𝐶¹+2𝑋₄𝑋₄+𝑌₄²−𝑌₄₂−2𝑋₄𝑌₄,(2.38)𝐵_4,4211=92𝑋⁴¹¹𝜇⁶𝑎_𝑒4𝐽₄₂32571−𝐶¹+2𝑋₄²𝑋₄+𝑌₄²𝐶¹×1+𝐶₁+2𝑋₄+2𝑋₄𝑋₄+𝑌₄𝑋⁴+𝑌₄⁻¹−1581+𝐶¹+2𝑋₄𝑋₄+𝑌₄²−𝑌₄₂−2𝑋₄𝑌₄(2.39)𝐵_4,4221=52𝑋⁴¹¹𝜇⁶𝑎_𝑒4𝐽₄₂110651−𝐶¹+2𝑋₄²𝑋₄+𝑌₄²1−3𝐶¹+2𝑋₄²+15𝑋₄+𝑌₄²4−145𝐶¹+2𝑋₄²𝑋₄+𝑌₄²−𝑌₄₂−2𝑋₄𝑌₄,(2.40)𝐵_4,4231=𝜇⁶𝑋⁴¹⁰𝑎_𝑒4𝐽₄₂−23571−𝐶¹+2𝑋₄²𝑋₄+𝑌₄²𝐶¹×1−𝐶₁+2𝑋₄+2𝑋₄𝑋₄+𝑌₄𝑋⁴+𝑌₄⁻¹−1581−𝐶¹+2𝑋₄×⎛⎜⎜⎜⎝1𝑋₄+𝑌₄²2−𝑌₄₂−2𝑋₄𝑌₄𝑋₄+3316−𝑌⁴²−2𝑋₄𝑌₄𝑋₄₂⎞⎟⎟⎟⎠.(2.41)

Since the term 𝜔_𝑒Θ₄ is constant, it plays no role in the equations of motion, and a new Hamiltonian can be introduced,𝐻₄=𝐻₄+𝜔_𝑒Θ₄.(2.42)

The dynamical system described by 𝐻₄ is given by𝑑𝑋₄,𝑌₄𝑑𝑡=𝜕𝜕𝑥𝐻₄,⁴𝑦₄,𝑑𝑥⁴,𝑦₄𝑑𝑡=−𝜕𝜕𝑋𝐻₄,⁴𝑌₄.(2.43)

The zonal harmonics used in (2.34) and (2.35) and the tesseral harmonics used in (2.36) to (2.41) are shown in Table 2.

Table 2 

The zonal and tesseral harmonics.

The constant of integration 𝐶₁ in (2.34) to (2.41) is given, in terms of the initial values of the orbital elements, 𝑎_𝑜, 𝑒_𝑜, and 𝐼_𝑜, by𝐶₁=√𝜇𝑎_𝑜1−𝑒^2𝑜cos𝐼_𝑜−2,(2.44) or, in terms of the variables 𝑋₄ and 𝑌₄,𝐶₁=𝑋₄cos𝐼_𝑜−2+𝑌₄cos𝐼_𝑜.(2.45)

In Section 4, some results of the numerical integration of (2.43) are shown.

3. Lyapunov Exponents

The estimation of the chaoticity of orbits is very important in the studies of dynamical systems, and possible irregular motions can be analyzed by Lyapunov exponents [23].

In this work, “Gram-Schmidt’s method,” described in [23–26], will be applied to compute the Lyapunov exponents. A brief description of this method is presented in what follows.

The dynamical system described by (2.43) can be rewritten as 𝑑𝑋₄𝑑𝑡=𝑃¹𝑋₄𝑑𝑌₄,𝑌₄,𝑥₄,𝑦₄;𝐶₁,𝑑𝑡=𝑃²𝑋₄𝑑𝑥₄,𝑌₄,𝑥₄,𝑦₄;𝐶₁,𝑑𝑡=𝑃³𝑋₄𝑑𝑦₄,𝑌₄,𝑥₄,𝑦₄;𝐶₁,𝑑𝑡=𝑃⁴𝑋₄,𝑌₄,𝑥₄,𝑦₄;𝐶₁.(3.1)

Introducing𝑧=⎛⎜⎜⎜⎜⎜⎜⎝𝑋⁴𝑌₄𝑥₄𝑦₄𝑍=⎛⎜⎜⎜⎜⎜⎜⎝𝑃₁⎞⎟⎟⎟⎟⎟⎟⎠,𝑃₂𝑃₃𝑃₄⎞⎟⎟⎟⎟⎟⎟⎠.(3.2)

Equations (3.2) can be put in the form𝑑𝑧𝑑𝑡=𝑍(𝑧).(3.3) The variational equations, associated to the system of differential equations (3.3), are given by𝑑𝜁𝑑𝑡=𝐽𝜁,(3.4) where 𝐽=(𝜕𝑍/𝜕𝑧) is the Jacobian.

The total number of differential equations used in this method is 𝑛(𝑛+1), 𝑛 represents the number of the motion equations describing the problem, in this case four. In this way, there are twenty differential equations, four are motion equations of the problem and sixteen are variational equations described by  (3.4)

The dynamical system represented by (3.3) and (3.4) is numerically integrated and the neighboring trajectories are studied using the Gram-Schmidt orthonormalization to calculate the Lyapunov exponents.

The method of the Gram-Schmidt orthonormalization can be seen in [25, 26] with more details. A simplified denomination of the method is described as follows.

Considering the solutions to (3.4) as 𝐮_𝜅(𝑡), the integration in the time 𝜏 begins from initial conditions 𝐮_𝜅(𝑡₀)=𝐞_𝜅(𝑡₀), an orthonormal basis.

At the end of the time interval, the volumes of the 𝜅-dimensional (𝜅=1,2,…,𝑁) produced by the vectors 𝐮_𝜅 are calculated by𝑉_𝜅=‖‖‖‖^𝜅_𝑗=1𝐮_𝑗(𝑡)‖‖‖‖,(3.5) where ⋀ is the outer product and ‖⋅‖ is a norm.

In this way, the vectors 𝐮_𝜅 are orthonormalized by Gram-Schmidt method. In other words, new orthonormal vectors 𝐞_𝜅(𝑡₀+𝜏) are calculated, in general, according to𝐞_𝜅=𝐮^𝜅−∑^{𝜅−1𝑗=1}𝐮_𝜅⋅𝐞_𝑗𝐞_𝑗‖‖𝐮_𝜅−∑^{𝜅−1𝑗=1}𝐮_𝜅⋅𝐞_𝑗𝐞_𝑗‖‖.(3.6)

The Gram-Schmidt method makes invariant the 𝜅-dimensional subspace produced by the vectors 𝐮₁,𝐮₂,𝐮₃,…,𝐮_𝜅 in constructing the new 𝜅-dimensional subspace spanned by the vectors 𝐞₁,𝐞₂,𝐞₃,…,𝐞_𝜅.

With new vector 𝐮_𝜅(𝑡₀+𝜏)=𝐞_𝜅(𝑡₀+𝜏), the integration is reinitialized and carried forward to 𝑡=𝑡₀+2𝜏. The whole cycle is repeated over a long-time interval. The theorems guarantee that the 𝜅-dimensional Lyapunov exponents are calculated by [25, 26]:𝜆(𝜅)=lim_𝑛→∞1𝑛𝜏^𝑛_𝑗=1ln𝑉_𝜅𝑡₀+𝑗𝜏ln𝑉_𝜅𝑡₀+(𝑗−1)𝜏.(3.7)

The theory states that if the Lyapunov exponent tends to a positive value, the orbit is chaotic.

In the next section are shown some results about the Lyapunov exponents.

4. Results

Figures 1, 2, 3, and 4 show the time behavior of the semimajor axis, 𝑥₄ angle, argument of perigee and of the eccentricity, according to the numerical integration of the motion equations, (2.43), considering three different resonant angles together: 𝜙₂₂₀₁, 𝜙₂₂₁₁, and 𝜙₂₂₂₁ associated to 𝐽₂₂, and three angles, 𝜙₄₂₁₁, 𝜙₄₂₂₁, and 𝜙₄₂₃₁ associated to 𝐽₄₂, with the same frequency of the resonant angles related to the 𝐽₂₂, but with different phase. The initial conditions corresponding to variables 𝑋₄ and 𝑌₄ are defined for 𝑒_𝑜=0.001, 𝐼_𝑜=55^∘, and 𝑎_𝑜 given in Table 3. The initial conditions of the variables 𝑥₄ and 𝑦₄ are 0° and 0°, respectively. Table 3 shows the values of 𝐶₁ corresponding to the given initial conditions.

Table 3 

Values of the constant of integration 𝐶₁ for 𝑒=0.001, 𝐼=55^∘ and different values for semimajor axis.

Figure 1 

Time behavior of the semimajor axis for different values of 𝐶1 given in Table 3.

Figure 2 

Time behavior of 𝑥4 angle for different values of 𝐶1 given in Table 3.

Figure 3 

Time behavior of the argument of pericentre for different values of 𝐶1 given in Table 3.

Figure 4 

Time behavior of the eccentricity for different values of 𝐶1 given in Table 3.

Figures 5, 6, 7, and 8 show the time behavior of the semimajor axis, 𝑥₄ angle, argument of perigee and of the eccentricity for two different cases. The first case considers the critical angles 𝜙₂₂₀₁, 𝜙₂₂₁₁, and 𝜙₂₂₂₁, associated to the tesseral harmonic 𝐽₂₂, and the second case considers the critical angles associated to the tesseral harmonics 𝐽₂₂ and 𝐽₄₂. The angles associated to the 𝐽₄₂, 𝜙₄₂₁₁, 𝜙₄₂₂₁, and 𝜙₄₂₃₁, have the same frequency of the critical angles associated to the 𝐽₂₂ with a different phase. The initial conditions corresponding to variables 𝑋₄ and 𝑌₄ are defined for 𝑒_𝑜=0.05, 𝐼_𝑜=10^∘, and 𝑎_𝑜 given in Table 4. The initial conditions of the variables 𝑥₄ and 𝑦₄ are 0° and 60^∘, respectively. Table 4 shows the values of 𝐶₁ corresponding to the given initial conditions.

Table 4 

Values of the constant of integration 𝐶₁ for 𝑒=0.05, 𝐼=10^∘, and different values for semimajor axis.

Figure 5 

Time behavior of the semimajor axis for different values of 𝐶1 given in Table 4.

Figure 6 

Time behavior of 𝑥4 angle for different values of 𝐶1 given in Table 4.

Figure 7 

Time behavior of the argument of pericentre for different values of 𝐶1 given in Table 4.

Figure 8 

Time behavior of the eccentricity for different values of 𝐶1 given in Table 4.

Analyzing Figures 5–8, one can observe a correction in the orbits when the terms related to the tesseral harmonic 𝐽₄₂ are added to the model. Observing, by the percentage, the contribution of the amplitudes of the terms 𝐵_4,4211, 𝐵_4,4221, and 𝐵_4,4231, in each critical angle studied, is about 1,66% up to 4,94%. In fact, in the studies of the perturbations in the artificial satellites motion, the accuracy is important, since adding different tesseral and zonal harmonics to the model, one can have a better description about the orbital motion.

Figures 9, 10, 11, and 12 show the time behavior of the semimajor axis, 𝑥₄ angle, argument of perigee and of the eccentricity, according to the numerical integration of the motion equations, (2.43), considering three different resonant angles together; 𝜙₂₂₀₁, 𝜙₂₂₁₁, and 𝜙₂₂₂₁ associated to 𝐽₂₂ and three angles 𝜙₄₂₁₁, 𝜙₄₂₂₁, and 𝜙₄₂₃₁ associated to 𝐽₄₂. The initial conditions corresponding to variables 𝑋₄ and 𝑌₄ are defined for 𝑒_𝑜=0.01, 𝐼_𝑜=55^∘, and 𝑎_𝑜 given in Table 5. The initial conditions of the variables 𝑥₄ and 𝑦₄ are 0° and 60^∘, respectively. Table 5 shows the values of 𝐶₁ corresponding to the given initial conditions.

Table 5 

Values of the constant of integration 𝐶₁ for 𝑒=0.01, 𝐼=55^∘, and different values for semimajor axis.

Figure 9 

Time behavior of the semimajor axis for different values of 𝐶1 given in Table 5.

Figure 10 

Time behavior of 𝑥4 angle for different values of 𝐶1 given in Table 5.

Figure 11 

Time behavior of the argument of pericentre for different values of 𝐶1 given in Table 5.

Figure 12 

Time behavior of the eccentricity for different values of 𝐶1 given in Table 5.

Analyzing Figures 1–12, one can observe possible irregular motions in Figures 1–4, specifically considering values for 𝐶₁=−1.467778013×10¹¹ m²/s and 𝐶₁=−1.467819454×10¹¹ m²/s, and, in Figures 9–12, for 𝐶₁=−1.467765786×10¹¹ m²/s and 𝐶₁=−1.467821043×10¹¹ m²/s. These curves will be analyzed by the Lyapunov exponents in a specified time verifying the possible regular or chaotic motions.

Figures 13 and 14 show the time behavior of the Lyapunov exponents for two different cases, according to the initial values of Figures 1–4 and 9–12. The dynamical system involves the zonal harmonics 𝐽₂₀ and 𝐽₄₀ and the tesseral harmonics 𝐽₂₂ and 𝐽₄₂. The method used in this work for the study of the Lyapunov exponents is described in Section 3. In Figure 13, the initial values for 𝐶₁, 𝑥₄, and 𝑦₄ are 𝐶₁=−1.467778013×10¹¹ m²/s and 𝐶₁=−1.467819454×10¹¹ m²/s, 𝑥₄=0^∘ and 𝑦₄=0^∘, respectively. In Figure 14, the initial values for 𝐶₁, 𝑥₄, and 𝑦₄ are 𝐶₁=−1.467765786×10¹¹ m²/s and 𝐶₁=−1.467821043×10¹¹ m²/s, 𝑥₄=0^∘ and 𝑦₄=60^∘, respectively. In each case are used two different values for semimajor axis corresponding to neighboring orbits shown previously in Figures 1–4 and 9–12.

Figure 13 

Lyapunov exponents 𝜆(1) and 𝜆(2), corresponding to the variables 𝑋4 and 𝑌4, respectively, for 𝐶1=−1.467778013×1011 m2/s and 𝐶1=−1.467819454×1011 m2/s, 𝑥4=0∘ and 𝑦4=0∘.

Figure 14 

Lyapunov exponents 𝜆(1) and 𝜆(2), corresponding to the variables 𝑋4 and 𝑌4, respectively, for 𝐶1=−1.467765786×1011 m2/s and 𝐶1=−1.467821043×1011 m2/s, 𝑥4=0∘ and 𝑦4=60∘.

Figures 13 and 14 show Lyapunov exponents for neighboring orbits. The time used in the calculations of the Lyapunov exponents is about 150.000 days. For this time, it can be observed in Figure 13 that 𝜆(1), corresponding to the initial value 𝑎(0)=26565.0 km, tends to a positive value, evidencing a chaotic region. On the other hand, analyzing the same Figure 13, 𝜆(1), corresponding to the initial value 𝑎(0)=26563.5 km, does not show a stabilization around the some positive value, in this specified time. Probably, the time is not sufficient for a stabilization in some positive value, or 𝜆(1), initial value 𝑎(0)=26563.5 km, tends to a negative value, evidencing a regular orbit. Analyzing now Figure 14, it can be verified that 𝜆(1), corresponding to the initial value 𝑎(0)=26564.0 km, tends to a positive value, it contrasts with 𝜆(1), initial value 𝑎(0)=26562.0 km. Comparing Figure 13 with Figure 14, it is observed that the Lyapunov exponents in Figure 14 has an amplitude of oscillation greater than the Lyapunov exponents in Figure 13. Analyzing this fact, it is probable that the necessary time for the Lyapunov exponent 𝜆(2), in Figure 14, to stabilize in some positive value is greater than the necessary time for the 𝜆(2) in Figure 13.

Rescheduling the axes of Figures 13 and 14, as described in Figures 15 and 16, respectively, the Lyapunov exponents tending to a positive value can be better visualized.

Figure 15 

Lyapunov exponents 𝜆(1) and 𝜆(2), corresponding to the variables 𝑋4 and 𝑌4, respectively, for 𝐶1=−1.467778013×1011 m2/s and 𝐶1=−1.467819454×1011 m2/s, 𝑥4=0∘ and 𝑦4=0∘.

Figure 16 

Lyapunov exponents 𝜆(1) and 𝜆(2), corresponding to the variables 𝑋4 and 𝑌4, respectively, for 𝐶1=−1.467765786×1011 m2/s and 𝐶1=−1.467821043×1011 m2/s, 𝑥4=0∘ and 𝑦4=60∘.

5. Conclusions

In this work, the dynamical behavior of three critical angles associated to the 2 : 1 resonance problem in the artificial satellites motion has been investigated.

The results show the time behavior of the semimajor axis, argument of perigee and e-ccentricity. In the numerical integration, different cases are studied, using three critical angles together: 𝜙₂₂₀₁, 𝜙₂₂₁₁, and 𝜙₂₂₂₁ associated to 𝐽₂₂ and 𝜙₄₂₁₁, 𝜙₄₂₂₁, and 𝜙₄₂₃₁ associated to the 𝐽₄₂.

In the simulations considered in the work, four cases show possible irregular motions for 𝐶₁=−1.467778013×10¹¹ m²/s, 𝐶₁=−1.467819454×10¹¹ m²/s, 𝐶₁=−1.467765786×10¹¹ m²/s, and 𝐶₁=−1.467821043×10¹¹ m²/s. Studying the Lyapunov exponents, two cases show chaotic motions for 𝐶₁=−1.467819454×10¹¹ m²/s and 𝐶₁=−1.467821043×10¹¹ m²/s.

Analyzing the contribution of the terms related to the 𝐽₄₂, it is observed that, for the value of 𝐶₁=−1.045724331×10¹¹ m²/s, the amplitudes of the terms 𝐵_4,4211, 𝐵_4,4221, and 𝐵_4,4231 are greater than the other values of 𝐶₁. In other words, for bigger values of semimajor axis, it is observed a smaller contribution of the terms related to the tesseral harmonic 𝐽₄₂.

The theory used in this paper for the 2 : 1 resonance can be applied for any resonance involving some artificial Earth satellite.

Acknowledgments

This work was accomplished with support of the FAPESP under the Contract no. 2009/00735-5 and 2006/04997-6, SP Brazil, and CNPQ (Contracts 300952/2008-2 and 302949/2009-7).