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 𝐽20 and 𝐽40 and the tesseral harmonics 𝐽22 and 𝐽42 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 [114]. 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 𝐽20 and 𝐽40 and the tesseral harmonics 𝐽22 and 𝐽42 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 𝐽22 and 𝐽42, 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×1014 m3/s2, 𝑎, 𝑒, 𝐼, Ω, 𝜔, 𝑀 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𝑒2,𝐻=𝜇𝑎1𝑒2cos(𝐼),=𝑀,𝑔=𝜔,=Ω.(2.3)𝐿, 𝐺, and 𝐻 represent the generalized coordinates, and , 𝑔, and represent the conjugate momenta.

Using the canonical variables, one gets the Hamiltonian 𝐹,𝜇𝐹=22𝐿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𝜇𝑎𝑒𝐿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𝐻=𝜇22𝑋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 𝑛=𝜇2/𝑋3. 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]:𝐻=14𝜋202𝜋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𝐻𝑟=𝜇22𝑋2𝜔𝑒Θ+𝑗=1𝐵2𝑗,0,𝑗,0(+𝑋,𝑌,𝑍)𝑙𝑙=2𝑙𝑚=2𝑝=0𝐵𝑙𝑚𝑝(𝛼𝑚)𝜙(𝑋,𝑌,𝑍)cos𝑙𝑚𝑝(𝛼𝑚).(𝑥,𝑦,𝑧,𝜃)(2.16)

Several authors, [11, 15, 1922], 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 as11𝛼𝑑𝑋+𝑑𝑡𝑑𝑌+𝑑𝑡𝑑𝑍𝑑𝑡=0.(2.24)

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

where 𝐶1 is an integration constant.

Using this first integral, a Mathieu transformation𝑋(𝑋,𝑌,𝑍,Θ,𝑥,𝑦,𝑧,𝜃)1,𝑌1,𝑍1,Θ1,𝑥1,𝑦1,𝑧1,𝜃1(2.26)

can be defined.

This transformation is given by the following equations:𝑋1=𝑋,𝑌1=𝑌,𝑍1=11𝛼𝑋+𝑌+𝑍,Θ1𝑥=Θ,11=𝑥1𝛼𝑧,𝑦1=𝑦𝑧,𝑧1=𝑧,𝜃1=𝜃.(2.27) The subscript 1 denotes the new set of canonical variables. Note that 𝑍1=𝐶1, and the 𝑧1 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, 𝑋1, 𝑌1, 𝑍1, Θ1, 𝑥1, 𝑦1, 𝑧1, 𝜃1, 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=𝜇22𝑋21𝜔𝑒Θ1+𝑗=1𝐵2𝑗,0,𝑗,0𝑋1,𝑌1,𝐶1+𝑙𝑙=2𝑙𝑚=2𝑝=0𝐵𝑙𝑚𝑝,(𝛼𝑚)𝑋1,𝑌1,𝐶1𝜙cos𝑙𝑚𝑝(𝛼𝑚)𝑥1,𝑦1,𝜃1.(2.28)

The Hamiltonian 𝐻rs has all resonant frequencies, relative to the commensurability 𝛼, where the 𝜙𝑙𝑚𝑝(𝛼𝑚) argument is given by𝜙𝑙𝑚𝑝(𝛼𝑚)=𝑚𝛼𝑥1𝜃1+(𝑙2𝑝𝛼𝑚)𝑦1𝜙𝑙𝑚𝑝(𝛼𝑚)0,(2.29)

with𝜙𝑙𝑚𝑝(𝛼𝑚)0=𝑚𝜆𝑙𝑚𝜋(𝑙𝑚)2.(2.30)

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

Each one of the frequencies contained in 𝑑𝑥1/𝑑𝑡, 𝑑𝑦1/𝑑𝑡, 𝑑𝜃1/𝑑𝑡 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𝐻1=𝜇22𝑋21𝜔𝑒Θ1+𝐵1,2010𝑋1,𝑌1,𝐶1+𝐵1,4020𝑋1,𝑌1,𝐶1+𝐵1,2201𝑋1,𝑌1,𝐶1𝑥cos12𝜃1+𝑦12𝜆22+𝐵1,2211𝑋1,𝑌1,𝐶1𝑥cos12𝜃1𝑦12𝜆22+𝐵1,2221𝑋1,𝑌1,𝐶1𝑥cos12𝜃13𝑦12𝜆22+𝐵1,4211𝑋1,𝑌1,𝐶1𝑥cos12𝜃1+𝑦12𝜆42+𝜋+𝐵1,4221𝑋1,𝑌1,𝐶1𝑥cos12𝜃1𝑦12𝜆42+𝜋+𝐵1,4231𝑋1,𝑌1,𝐶1𝑥cos12𝜃13𝑦12𝜆42.+𝜋(2.31)

The Hamiltonian 𝐻1 is defined considering a fixed resonance and three different critical angles associated to the tesseral harmonic 𝐽22; the critical angles associated to the tesseral harmonic 𝐽42 have the same frequency of the critical angles associated to the 𝐽22 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 𝐻1.

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

So, considering (2.31) and (2.32), the Hamiltonian 𝐻4 is found to be𝐻4=𝜇22𝑋24𝜔𝑒Θ42𝑋4+𝐵4,2010𝑋4,𝑌4,𝐶1+𝐵4,4020𝑋4,𝑌4,𝐶1+𝐵4,2201𝑋4,𝑌4,𝐶1𝑥cos4+𝑦42𝜆22+𝐵4,2211𝑋4,𝑌4,𝐶1𝑥cos4𝑦42𝜆22+𝐵4,2221𝑋4,𝑌4,𝐶1𝑥cos43𝑦42𝜆22+𝐵4,4211𝑋4,𝑌4,𝐶1𝑥cos4+𝑦42𝜆42+𝜋+𝐵4,4221𝑋4,𝑌4,𝐶1𝑥cos4𝑦42𝜆42+𝜋+𝐵4,4231𝑋4,𝑌4,𝐶1𝑥cos43𝑦42𝜆42,+𝜋(2.33) with 𝜔𝑒Θ4 constant and 𝐵4,2010=𝜇4𝑋64𝑎𝑒2𝐽2034𝐶1+2𝑋42𝑋4+𝑌42+1431+2𝑌422𝑋4𝑌4𝑋42,(2.34)𝐵4,4020=𝜇6𝑋410𝑎𝑒4𝐽40105𝐶6411+2𝑋42𝑋4+𝑌42232+158𝐶1+2𝑋42𝑋4+𝑌42×1+5𝑌422𝑋4𝑌4𝑋42,(2.35)𝐵4,2201=218𝑋74𝜇4𝑎𝑒2𝐽22𝐶1+1+2𝑋4𝑋4+𝑌42𝑌422𝑋4𝑌4,(2.36)𝐵4,2211=32𝑋74𝜇4𝑎𝑒2𝐽223232𝐶1+2𝑋42𝑋4+𝑌42𝑌422𝑋4𝑌4,(2.37)𝐵4,22213=8𝑋74𝜇4𝑎𝑒2𝐽22𝐶11+2𝑋4𝑋4+𝑌42𝑌422𝑋4𝑌4,(2.38)𝐵4,4211=92𝑋411𝜇6𝑎𝑒4𝐽4235𝐶2711+2𝑋42𝑋4+𝑌42𝐶1+2𝑋4×𝐶1+1+2𝑋4𝑋4+𝑌4𝑋4+𝑌41158𝐶1+1+2𝑋4𝑋4+𝑌42𝑌422𝑋4𝑌4(2.39)𝐵4,4221=52𝑋411𝜇6𝑎𝑒4𝐽42105𝐶1611+2𝑋42𝑋4+𝑌42𝐶131+2𝑋42𝑋4+𝑌42+154154𝐶1+2𝑋42𝑋4+𝑌42𝑌422𝑋4𝑌4,(2.40)𝐵4,4231=𝜇6𝑋410𝑎𝑒4𝐽4235𝐶2711+2𝑋42𝑋4+𝑌42𝐶1+2𝑋4×𝐶11+2𝑋4𝑋4+𝑌4𝑋4+𝑌41158𝐶11+2𝑋4𝑋4+𝑌42×12𝑌422𝑋4𝑌4𝑋4+3316𝑌422𝑋4𝑌4𝑋42.(2.41)

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

The dynamical system described by 𝐻4 is given by𝑑𝑋4,𝑌4=𝜕𝐻𝑑𝑡4𝜕𝑥4,𝑦4,𝑑𝑥4,𝑦4𝜕𝐻𝑑𝑡=4𝜕𝑋4,𝑌4.(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.

The constant of integration 𝐶1 in (2.34) to (2.41) is given, in terms of the initial values of the orbital elements, 𝑎𝑜, 𝑒𝑜, and 𝐼𝑜, by𝐶1=𝜇𝑎𝑜1𝑒2𝑜𝐼cos𝑜2,(2.44) or, in terms of the variables 𝑋4 and 𝑌4,𝐶1=𝑋4𝐼cos𝑜2+𝑌4𝐼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 [2326], 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 𝑑𝑋4𝑑𝑡=𝑃1𝑋4,𝑌4,𝑥4,𝑦4;𝐶1,𝑑𝑌4𝑑𝑡=𝑃2𝑋4,𝑌4,𝑥4,𝑦4;𝐶1,𝑑𝑥4𝑑𝑡=𝑃3𝑋4,𝑌4,𝑥4,𝑦4;𝐶1,𝑑𝑦4𝑑𝑡=𝑃4𝑋4,𝑌4,𝑥4,𝑦4;𝐶1.(3.1)

Introducing𝑋𝑧=4𝑌4𝑥4𝑦4,𝑃𝑍=1𝑃2𝑃3𝑃4.(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 𝐮𝜅(𝑡0)=𝐞𝜅(𝑡0), 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 𝐞𝜅(𝑡0+𝜏) 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 𝐮1,𝐮2,𝐮3,,𝐮𝜅 in constructing the new 𝜅-dimensional subspace spanned by the vectors 𝐞1,𝐞2,𝐞3,,𝐞𝜅.

With new vector 𝐮𝜅(𝑡0+𝜏)=𝐞𝜅(𝑡0+𝜏), the integration is reinitialized and carried forward to 𝑡=𝑡0+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𝑛𝜏𝑛𝑗=1𝑉ln𝜅𝑡0+𝑗𝜏𝑉ln𝜅𝑡0.+(𝑗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, 𝑥4 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: 𝜙2201, 𝜙2211, and 𝜙2221 associated to 𝐽22, and three angles, 𝜙4211, 𝜙4221, and 𝜙4231 associated to 𝐽42, with the same frequency of the resonant angles related to the 𝐽22, but with different phase. The initial conditions corresponding to variables 𝑋4 and 𝑌4 are defined for 𝑒𝑜=0.001, 𝐼𝑜=55, and 𝑎𝑜 given in Table 3. The initial conditions of the variables 𝑥4 and 𝑦4 are 0° and 0°, respectively. Table 3 shows the values of 𝐶1 corresponding to the given initial conditions.

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

Analyzing Figures 58, one can observe a correction in the orbits when the terms related to the tesseral harmonic 𝐽42 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, 𝑥4 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; 𝜙2201, 𝜙2211, and 𝜙2221 associated to 𝐽22 and three angles 𝜙4211, 𝜙4221, and 𝜙4231 associated to 𝐽42. The initial conditions corresponding to variables 𝑋4 and 𝑌4 are defined for 𝑒𝑜=0.01, 𝐼𝑜=55, and 𝑎𝑜 given in Table 5. The initial conditions of the variables 𝑥4 and 𝑦4 are 0° and 60, respectively. Table 5 shows the values of 𝐶1 corresponding to the given initial conditions.

Analyzing Figures 112, one can observe possible irregular motions in Figures 14, specifically considering values for 𝐶1=1.467778013×1011 m2/s and 𝐶1=1.467819454×1011 m2/s, and, in Figures 912, for 𝐶1=1.467765786×1011 m2/s and 𝐶1=1.467821043×1011 m2/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 14 and 912. The dynamical system involves the zonal harmonics 𝐽20 and 𝐽40 and the tesseral harmonics 𝐽22 and 𝐽42. 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 𝐶1, 𝑥4, and 𝑦4 are 𝐶1=1.467778013×1011 m2/s and 𝐶1=1.467819454×1011 m2/s, 𝑥4=0 and 𝑦4=0, respectively. In Figure 14, the initial values for 𝐶1, 𝑥4, and 𝑦4 are 𝐶1=1.467765786×1011 m2/s and 𝐶1=1.467821043×1011 m2/s, 𝑥4=0 and 𝑦4=60, respectively. In each case are used two different values for semimajor axis corresponding to neighboring orbits shown previously in Figures 14 and 912.

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.

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: 𝜙2201, 𝜙2211, and 𝜙2221 associated to 𝐽22 and 𝜙4211, 𝜙4221, and 𝜙4231 associated to the 𝐽42.

In the simulations considered in the work, four cases show possible irregular motions for 𝐶1=1.467778013×1011 m2/s, 𝐶1=1.467819454×1011 m2/s, 𝐶1=1.467765786×1011 m2/s, and 𝐶1=1.467821043×1011 m2/s. Studying the Lyapunov exponents, two cases show chaotic motions for 𝐶1=1.467819454×1011 m2/s and 𝐶1=1.467821043×1011 m2/s.

Analyzing the contribution of the terms related to the 𝐽42, it is observed that, for the value of 𝐶1=1.045724331×1011 m2/s, the amplitudes of the terms 𝐵4,4211, 𝐵4,4221, and 𝐵4,4231 are greater than the other values of 𝐶1. In other words, for bigger values of semimajor axis, it is observed a smaller contribution of the terms related to the tesseral harmonic 𝐽42.

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).