ABSTRACT

For certain values of semi-major axis and eccentricity, orbit plane precession caused by Earth oblate is synchronous with the mean orbital motion of the apparent Sun (a sun-synchronism). Many forces cause slow changes in the inclination and ascending node of sun-synchronous orbits. In this work, we investigate the analytical perturbations due to the direct solar radiation pressure and gravitational waves effects. A full analytical solution is obtained using technique of canonical Lie-transformation up to the order three in $J_2$ (the oblateness of the Earth). The solar radiation pressure and gravitational waves perturbations cause second order effects on all the elements of the elliptic orbit (the eccentricity, inclination, ascending node, argument of perigee, and semi-major axis) consequently these perturbations will cause disturbance in the sun-synchronism. Also we found that the perturbation or the behavior of gravitational waves almost the same as the perturbation or the behavior of solar radiation pressure and their coupling will incorporate the sun-synchronism through the secular rate of the ascending node precession.

Keywords:

Orbital Mechanics, Celestial Mechanics, Astrodynamics, Gravitational Waves, Perturbations

1. Introduction

It is well known that the sun-synchronous orbits are useful for remote sensing mission, because the illumination of the sub satellite point (the point where a straight line drawn from a satellite to the center of the Earth intersects the Earth’s surface) is nearly constant on each successive pass. Also sun-synchronous orbits must be retrograde (i.e. inclination greater than 90˚). If the sun-synchronous orbit is high enough and the sub satellite point is initiated outside the Earth’s shadow, then the illumination of the sub satellite point will be remain continuously throughout the year. Since a sun-synchronous orbits follows the Sun, so the solar pressure perturbations cannot be neglected, if the orbit is un-shadowed, long period occur in all the orbital elements except semi-major axis [1] [2] . It is well known that Earth oblateness causes the orbital angular momentum vector of earth satellite to process about the Earth’s polar axis. For certain value of semi-major axis and eccentricity orbit plane precession caused by oblateness is synchronous with the mean orbital motion of the apparent Sun (a sun-synchro- nism). Gravitational perturbations arising from a spherical Earth, Moon, and Sun cause long period fluctuations in the mean argument of perigee, eccentricity, inclination, and ascending node [3] . Relativistic effects with SRP incorporate the sun-synchronism [4] . The short and long period’s perturbations due to the effect of GW on the orbital elements of Jupiter are studies by [5] . Semi-analytical theory of the mean orbital motion due to the effect of gravitational waves is investigated by [6] . Many authors discussed the effect of solar radiation pressure and the effects of gravitational waves separately but the purpose of the current contribution is to investigate the coupling effects of solar radiation pressure and gravitational waves on sun-synchronism of elliptic orbits, using the perturbation technique based on the canonical Lie-transformation. Two canonical transformations obtained to eliminate the short-period and long-period terms. The solution is obtained up to order three and two in secular and periodic terms respectively assuming the effect of direct solar radiation pressure and the gravitational waves of the same order of magnitude and equal to the order of $J_2^2$ (the second order zonal harmonic of the geo potential i.e. 10⁻⁶).

2. The Formulation

The Hamiltonian which describes the moving Earth’s satellite under the influence of SRP and the GW can be expressed as

$H=H_g+H_s+H_w$ (1)

where $H_g$ represents the Hamiltonian of a satellite in the Earth’s gravitational field

$H_g=vP-L$ (2)

v is the velocity of satellite, P is the canonical momentum, and L is the Lagrangian of a satellite moving about a central body in Einstein form of the post Newtonian gravitational field [7]

$L=\frac{1}{2}{v}^2-U_N+\frac{1}{c_l^2}\left\{\frac{v^4}{8}-\frac{1}{2}{U}_N^2-\frac{3}{2}{v}^2{U}_N\right\}$ (3)

where $c_l$ is the speed of light and $U_N$ the Earth’s gravitational potential is developed in series up the fourth harmonic as [8]

$\begin{array}{c}{U}_N=-\frac{\mu }{r}\{1-\frac{J_2}{2r^2}{R}^2\left(3{\sin }^2\delta -1\right)-\frac{J_3}{2r^3}{R}^3\left(5{\sin }^3\delta -3\sin \delta \right)\\ -\frac{J_4}{8r^4}{R}^4\left(35{\sin }^4\delta -30{\sin }^2\delta +3\right)\}\end{array}$

Using $\sin \delta =\sin i\sin \left(f+\omega \right)$ ,$\frac{p^2}{2}=\frac{\mu }{r}-\frac{\mu }{2a}$ , and $\frac{p^4}{8}=\frac{1}{2}\left(\frac{{\mu }^2}{r^2}-\frac{{\mu }^2}{ra}+\frac{{\mu }^2}{4a^2}\right)$

where R is the equatorial radius of the Earth, δ is the latitude referred to the equatorial plane, i is the inclination of the orbit plane, f the true anomaly, ω the argument of perigee, and is the semi-major axis as in Figure 1. Then the Hamiltonian $H_g$ will be in the form

$\begin{array}{c}{H}_g=-\frac{\mu }{2a}+\frac{J_2\mu }{4r^3}{R}^2\left\{{\theta }_0+{\theta }_1\cos 2\left(f+\omega \right)\right\}\\ +\frac{J_3\mu }{8r^4}{R}^3\left\{{\theta }_2\sin \left(f+\omega \right)+{\theta }_3\sin 3\left(f+\omega \right)\right\}\\ +\frac{J_4\mu }{64r^5}{R}^4\left\{{\theta }_4+{\theta }_5\cos 2\left(f+\omega \right)+{\theta }_6\cos 4\left(f+\omega \right)\right\}\end{array}$

$\begin{array}{l}+\frac{1}{c_l^2}\left(-3\frac{{\mu }^2}{r^2}+2\frac{{\mu }^2}{ra}-\frac{{\mu }^2}{8a^2}\right)+\frac{J_2}{c_l^2}\left(\frac{{\mu }^2{R}^2}{2r^4}-\frac{3{\mu }^2{R}^2}{8a{r}^3}\right)\left\{{\theta }_0+{\theta }_1\cos 2\left(f+\omega \right)\right\}\\ +\frac{J_3}{c_l^2}\left(\frac{{\mu }^2{R}^3}{4r^5}-\frac{3{\mu }^2{R}^3}{16a{r}^4}\right)\left\{{\theta }_2\sin \left(f+\omega \right)+{\theta }_3\sin 3\left(f+\omega \right)\right\}\\ +\frac{J_4}{c_l^2}\left(\frac{{\mu }^2{R}^4}{32r^6}-\frac{3{\mu }^2{R}^4}{128a{r}^5}\right)\left\{{\theta }_4+{\theta }_5\cos 2\left(f+\omega \right)+{\theta }_6\cos 4\left(f+\omega \right)\right\}\\ +\frac{J_2^2}{c_l^2}\left(\frac{{\mu }^2{R}^4}{32r^6}\right)\left\{{\theta }_7+{\theta }_8\cos 2\left(f+\omega \right)+{\theta }_9\cos 4\left(f+\omega \right)\right\}\end{array}$

In terms of the Delaunay set of canonical variables $\left(l,g,h,L,G,H\right)$ which are defined in terms of Keplerian elements $\left(a,e,i,\Omega ,M\right)$ by

$\begin{array}{l}l=M=n\left(t–\tau \right);L=\sqrt{\mu a}\\ g=\omega ;G=\sqrt{\mu a}\left(1-e^2\right)\\ h=\Omega ;H=G\cos i\end{array}$ (4)

$l$ is the mean anomaly, $\omega$ argument of perigee, $\Omega$ longitude of the node, and $a,e,i$ are the semi-major axis, the eccentricity, and the inclination respectively. The Hamiltonian $H_g$ will be in the form

$\begin{array}{c}{H}_g=-\frac{\mu }{2L^2}+\frac{J_2{\mu }^4}{4L^6}{R}^2{\left(\frac{a}{r}\right)}^3\left\{{\theta }_0+{\theta }_1\cos 2\left(f+g\right)\right\}\\ +\frac{J_3{\mu }^5}{8L^8}{R}^3{\left(\frac{a}{r}\right)}^4\left\{{\theta }_2\sin \left(f+g\right)+{\theta }_3\sin 3\left(f+g\right)\right\}\\ +\frac{J_4{\mu }^6}{64L^{10}}{R}^4{\left(\frac{a}{r}\right)}^5\left\{{\theta }_4+{\theta }_5\cos 2\left(f+g\right)+{\theta }_6\cos 4\left(f+g\right)\right\}\\ +\frac{{\mu }^4}{c_l^2{L}^4}\left\{-3{\left(\frac{a}{r}\right)}^2+2\left(\frac{a}{r}\right)-\frac{1}{8}\right\}\end{array}$

$\begin{array}{l}+\frac{J_2{\mu }^6}{c_l^2{L}^8}{R}^2\left\{\frac{1}{2}{\left(\frac{a}{r}\right)}^4-\frac{3}{8}{\left(\frac{a}{r}\right)}^3\right\}\left\{{\theta }_0+{\theta }_1\cos 2\left(f+g\right)\right\}\\ +\frac{J_3{\mu }^7}{c_l^2{L}^{10}}{R}^3\left\{\frac{1}{4}{\left(\frac{a}{r}\right)}^5-\frac{3}{16}{\left(\frac{a}{r}\right)}^4\right\}\left\{{\theta }_2\sin \left(f+g\right)+{\theta }_3\sin 3\left(f+g\right)\right\}\\ +\frac{J_4{\mu }^8}{c_l^2{L}^{12}}{R}^4\left\{\frac{1}{32}{\left(\frac{a}{r}\right)}^6-\frac{3}{128}{\left(\frac{a}{r}\right)}^5\right\}\left\{{\theta }_4+{\theta }_5\cos 2\left(f+g\right)+{\theta }_6\cos 4\left(f+g\right)\right\}\\ +\frac{J_2{}^2}{c_l^2}\left(\frac{{\mu }^8{R}^4}{32L^{12}}\right){\left(\frac{a}{r}\right)}^6\left\{{\theta }_7+{\theta }_8\cos 2\left(f+g\right)+{\theta }_9\cos 4\left(f+g\right)\right\}\end{array}$ (5)

where

${\theta }_0=1-3c^2$ , ${\theta }_1=-3s^2$ , ${\theta }_2=15s^3-12s^2$ , ${\theta }_3=-5s^3$ , ${\theta }_4=9-90c^2-105c^4$ , ${\theta }_5=-20+160c^2-140c^4$ , ${\theta }_6=35s^4$ , ${\theta }_7=\frac{27}{2}{c}^4-15c^2+\frac{11}{2}$ , ${\theta }_8=-6+24c^2-18c^4$ , ${\theta }_9=\frac{9}{2}{s}^4$ , and $c=\cos i$ and $s=\sin i$ .

The Hamiltonian of the solar radiation pressure $H_s$ in terms of Delaunay elements is given by [9]

$H_s=n_{\odot }{L}_{4s}-\beta \frac{L^2}{\mu a_{\odot }^2}\frac{r}{a}\left\{T_c\cos f+T_s\sin f\right\}$ (6)

where

$\begin{array}{c}{T}_c=\frac{1}{4}\left(1+c\right)\left(1+{\theta }_{\odot }\right)\cos \left(g+h-l_{4s}\right)\\ +\frac{1}{4}\left(1+c\right)\left(1-{\theta }_{\odot }\right)\cos \left(g+h+l_{4s}\right)\\ +\frac{1}{4}\left(1-c\right)\left(1+{\theta }_{\odot }\right)\cos \left(g-h+l_{4s}\right)\\ +\frac{1}{4}\left(1-c\right)\left(1-{\theta }_{\odot }\right)\cos \left(g-h-l_{4s}\right)\\ +\frac{1}{2}s\sin є\left\{\cos \left(g-l_{4s}\right)-\cos \left(g+l_{4s}\right)\right\}\end{array}$

and

$\begin{array}{c}{T}_s=-\frac{1}{4}\left(1+c\right)\left(1+{\theta }_{\odot }\right)\sin \left(g+h-l_{4s}\right)\\ -\frac{1}{4}\left(1+c\right)\left(1-{\theta }_{\odot }\right)\sin \left(g+h+l_{4s}\right)\\ -\frac{1}{4}\left(1-c\right)\left(1+{\theta }_{\odot }\right)\sin \left(g-h+l_{4s}\right)\\ -\frac{1}{4}\left(1-c\right)\left(1-{\theta }_{\odot }\right)\sin \left(g-h-l_{4s}\right)\\ -\frac{1}{2}s\sin є\left\{\sin \left(g-l_{4s}\right)-\sin \left(g+l_{4s}\right)\right\}\end{array}$

$l_{4s}=k_s=n_{\odot }t$ ( $n_{\odot }$ is the solar mean motion of order unity rather than $\text{O}\left(H_s\right)$ ), $L_{4s}=K_s$ is the conjugate of $L_{4s}$ which will be found from

$L_{4s}=-\frac{\partial H_s}{\partial l_{4s}}$ , ${\theta }_{\odot }=\cos є$ , $є$ is the obliquity of the ecliptic, and β is constant

which represents satellite area to mass ratio corresponding to a spherical satellite.

Finally the Hamiltonian of the gravitational waves $H_w$ in terms of Delaunay elements is given by

$\begin{array}{c}{H}_w=n_w{L}_{4w}+ϵL^4\frac{r^2}{a^2}\{\frac{s^2}{8}[\cos \left(2h+l_{4w}+{\alpha }_1\right)-\cos \left(2h-l_{4w}-{\alpha }_1\right)\\ -A_g\left[\sin \left(2h+l_{4w}+{\alpha }_2\right)+\sin \left(2h-l_{4w}-{\alpha }_2\right)\right]]\\ +\left(\frac{c^2}{16}+\frac{c}{8}-\frac{1}{16}\right)[\cos \left(2f+2g+2h+l_{4w}+{\alpha }_1\right)\\ +\cos \left(2f+2g+2h-l_{4w}-{\alpha }_1\right)-A_g[\sin \left(2f+2g+2h+l_{4w}+{\alpha }_2\right)\end{array}$

$\begin{array}{l}+\sin \left(2f+2g+2h-l_{4w}-{\alpha }_2\right)+\sin \left(2f+2g-2h+l_{4w}+{\alpha }_2\right)\\ +\sin \left(2f+2g-2h-l_{4w}-{\alpha }_2\right)]]+\left(\frac{c^2}{16}-\frac{c}{8}+\frac{1}{16}\right)\\ \begin{array}{c}\\ \end{array}\times \left[\cos \left(2f+2g-2h+l_{4w}+{\alpha }_1\right)+\cos \left(2f+2g-2h-l_{4w}-{\alpha }_1\right)\right]\}\end{array}$ (7)

where $ϵ=\frac{1}{2}{n}_w^2{h}_x$ and $A_g=-\frac{h_x}{h_{+}}$ , $n_w$ is the frequency of the wave. ${\alpha }_1$ and

${\alpha }_2$ are the phase difference, $h_{+}$ and $h_x$ are the amplitudes of the wave in the two orthogonal directions in the transverse plane. $l_{4w}=k_w=n_{w}t$ , $L_{4w}=K_w$ is

the conjugate of $l_{4w}$ which will be found from $L_{4w}=-\frac{\partial H_s}{\partial l_{4w}}$ .

According to Equations (5), (6), and (7) we can write the Hamiltonian H in (1) as

$H=H_0+{\displaystyle {\sum }_{n=1}^3\frac{{є}^n}{n!}{H}_n}$

where є is the small parameter of the expansion which represents the second zonal harmonic $J_2$ of the Earth’s potential, and

$H_0=-\frac{{\mu }^2}{2L^2}$ (8)

$H_1=-\frac{A_1}{L^6}{\left(\frac{a}{r}\right)}^3\left\{{\theta }_0+{\theta }_1\cos 2\left(f+g\right)\right\}+n_{\odot }{L}_{4s}+n_w{L}_{4w}$ (9)

$H_2=H_{21}+H_{22}+H_{23}+H_{24}+H_{25}$ (10)

$H_{21}=\frac{A_{21}}{L^8}{\left(\frac{a}{r}\right)}^4\left\{{\theta }_2\sin \left(f+g\right)+{\theta }_3\sin 3\left(f+g\right)\right\}$ (11)

$H_{22}=\frac{A_{22}}{L^{10}}{\left(\frac{a}{r}\right)}^5\left\{{\theta }_4+{\theta }_5\cos 2\left(f+g\right)+{\theta }_6\cos 4\left(f+g\right)\right\}$ (12)

$H_{23}=\frac{A_{23}}{L^4}\left\{-3{\left(\frac{a}{r}\right)}^2+2\left(\frac{a}{r}\right)-\frac{1}{8}\right\}$ (13)

$H_{24}=A_{24}{L}^2\frac{r}{a}\left\{T_c\cos f+T_s\sin f\right\}$ (14)

$\begin{array}{c}{H}_{25}=A_{25}{L}^4\frac{r^2}{a^2}\{\frac{s^2}{8}[\cos \left(2h+l_{4w}+{\alpha }_1\right)-\cos \left(2h-l_{4w}-{\alpha }_1\right)\\ -A_g\left[\sin \left(2h+l_{4w}+{\alpha }_2\right)+\sin \left(2h-l_{4w}-{\alpha }_2\right)\right]]\\ +\left(\frac{c^2}{16}+\frac{c}{8}-\frac{1}{16}\right)[\cos \left(2f+2g+2h+l_{4w}+{\alpha }_1\right)\\ +\cos \left(2f+2g+2h-l_{4w}-{\alpha }_1\right)-A_g[\sin \left(2f+2g+2h+l_{4w}+{\alpha }_2\right)\end{array}$

$\begin{array}{l}+\sin \left(2f+2g+2h-l_{4w}-{\alpha }_2\right)+\sin \left(2f+2g-2h+l_{4w}+{\alpha }_2\right)\\ +\sin \left(2f+2g-2h-l_{4w}-{\alpha }_2\right)]+\left(\frac{c^2}{16}-\frac{c}{8}+\frac{1}{16}\right)\\ \begin{array}{c}\\ \end{array}\times \left[\cos \left(2f+2g-2h+l_{4w}+{\alpha }_1\right)+\cos \left(2f+2g-2h-l_{4w}-{\alpha }_1\right)\right]\}\end{array}$ (15)

$H_3=\frac{A_3}{L^8}\left\{\frac{1}{2}{\left(\frac{a}{r}\right)}^4-\frac{3}{8}{\left(\frac{a}{r}\right)}^3\right\}\left\{{\theta }_0+{\theta }_1\cos 2\left(f+g\right)\right\}$ (16)

For different orders of $J_2$ , the coefficients A’s are

$A_1=\frac{{\mu }^4}{4}{R}^4$ , $A_{21}=\frac{J_3{\mu }^5}{4J_2{}^2}{R}^3$ , $A_{22}=\frac{J_4{\mu }^6}{32J_2{}^2}{R}^4$ , $A_{23}=\frac{2!{\mu }^6}{c_l^2{J}_2{}^2}{R}^2$ ,

$A_{24}=-\frac{\beta }{J_2{}^2\mu a_{\odot }^2}$ , $A_{25}=\frac{1}{2}\frac{n_w^2{h}_x}{J_2{}^2\mu }$ , $A_3=\frac{3!{\mu }^6}{c_l^2{J}_2{}^2}{R}^2$ .

The second-order term in the Hamiltonian H which is given in Equation (10) is consisting of the non-relativistic terms $H_{21}$ , $H_{22}$ , the relativistic term $H_{23}$ , the terms of solar radiation pressure $H_{24}$ and the terms of gravitational waves $H_{25}$ .The third order term $H_3$ is representing the mixed perturbation due to Newtonian quadruple field and Schwarzschild acceleration.

3. The Algorithm of Solution

Let the considered system of differential equations in terms of canonical elements under the effect of solar radiation pressure and GW are

$\frac{\text{d}{L}_i}{\text{d}t}=-\frac{\partial H}{\partial l_i};\frac{\text{d}{l}_i}{\text{d}t}=\frac{\partial H}{\partial L_i}\left(i=1,2,3,4\right)$ (17)

where $L_i$ and $l_i$ are the set of the Delaunay variable in Equation (4). H is the Hamiltonian in equations from (9) to (16).We use Lie transform to solve our problem [10] [11] . This method is one of the canonical perturbation methods proposed to build successive canonical transformation for Hamiltonian systems depending on a small parameter based on the consideration of Lie series and Lie transform. The system of Equation (17) will be solved up to order three in secular terms by performing two successive canonical transformation $\left(l_i,L_i;є\right)\to \left({l^{\prime }}_i,{L^{\prime }}_i\right)$ and $\left({l^{\prime }}_i,{L^{\prime }}_i;є\right)\to \left({l^{″}}_i,{L^{″}}_i\right)$ analytic in $є$ at $є=0$ to eliminate the short and long period terms, respectively from the Hamiltonian, where the primes indicate the transformed variables. The transformed Hamiltonians are in the form

$H^{*}\left(-,l_2,-,L_i;є\right)=H_0^{\ast }\left(L_i\right)+{\displaystyle {\sum }_{n=1}\frac{{є}^n}{n!}{H}_n^{\ast }\left(-,L_i\right)}$ (18)

$H^{**}\left(-,L_i;є\right)=H_0^{\ast \ast }\left(L_i\right)+{\displaystyle {\sum }_{n=1}\frac{{є}^n}{n!}{H}_n^{\ast \ast }\left(-,L_i\right)}$ (19)

With generators $W$ and $W^{*}$ expandable as

$W\left(l_i,L_i;є\right)={\displaystyle {\sum }_{n=0}\frac{{є}^n}{n!}{W}_{n+1}\left(l_i,L_i\right)}$ (20)

$W^{*}\left(-,l_2,l_3,L_i;є\right)={\displaystyle {\sum }_{n=0}\frac{{є}^n}{n!}{W}_{n+1}^{\ast }\left(-,l_2,l_3,L_i\right)}$ (21)

The algorithm for eliminating the short-period terms (i.e. those depending on $l$ and $l_4$ ) determining the new Hamiltonian $H^{*}$ and the generator $W$ is

Zero order: $H_0^{\ast }=H_0$

First order: $H_1^{\ast }=H_{1s}$

$W_1=\frac{{L^{\prime }}^3}{{\mu }^3}{\displaystyle \int H_{1p}\text{d}{l}^{\prime }}$

Second order: $H_2^{\ast }=H_{2s}+{\left\{H_1+H_1^{\ast },W_1\right\}}_s$

$W_2=\frac{{L^{\prime }}^3}{{\mu }^3}{{\displaystyle \int }}^{\text{}}\left(H_{2p}+{\left\{H_1+H_1^{\ast },W_1\right\}}_p\right)\text{d}{l}^{\prime }$

Third order: $\begin{array}{c}{H}_3^{\ast }=H_{3s}+{\left\{H_2,W_1\right\}}_s+2{\left\{H_2^{\ast },W_1\right\}}_s+{\left\{H_1,W_2\right\}}_s\\ +{\left\{H_1^{\ast },W_2\right\}}_s-{\left\{\left\{H_1^{\ast },W_1\right\},W_1\right\}}_s\end{array}$

where the brackets represent Poisson bracket and subscripts s and p indicates the secular (averaged over ${l^{\prime }}_1$ and ${l^{\prime }}_4$ ) and periodic parts in $l^{\prime }$ . The remaining long-period problem is then governed by a Hamiltonian depending only on ( $L^{\prime }$ , $G^{\prime }$ , $H^{\prime }$ and g). The second algorithm of canonical transformation to remove the long-period terms (i.e. those depending on $g^{\prime }$ ) and determining the new Hamiltonian $H^{**}$ and the generating function $W^{*}$ is given by

Zero order: $H_0^{\ast \ast }=H_0$

First order: $H_1^{\ast \ast }=H_1^{\ast }$

Second order: $H_2^{\ast \ast }=H_2^{\ast }$

$W_1^{\ast }=-\frac{1}{2}\frac{\partial H_1^{-1}}{\partial G}{{\displaystyle \int }}^{\text{}}{H}_{2p}^{\ast }\text{d}{g}^{″}$

Third order:

$H_3^{**}=H_3^{*}+{\left\{H_2^{*},W_1^{\ast }\right\}}_s+2{\left\{H_2^{**},W_1^{\ast }\right\}}_s+2{\left\{H_1^{\ast },W_2^{*}\right\}}_s-{\left\{\left\{H_1^{**},W_1^{\ast }\right\},W_1^{\ast }\right\}}_s$

$\begin{array}{c}{W}_2^{*}=-\frac{1}{3}\frac{\partial H_1^{-1}}{\partial G}{\displaystyle \int ({\left\{H_2^{*},W_1^{\ast }\right\}}_p+2{\left\{H_2^{**},W_1^{\ast }\right\}}_p+2{\left\{H_1^{\ast },W_2^{*}\right\}}_p}\\ -{\left\{\left\{H_1^{**},W_1^{\ast }\right\},W_1^{\ast }\right\}}_p)\text{d}{g}^{″}\end{array}$

where the subscripts s and p indicates the secular (averaged over $g^{\prime }$ ) and the periodic parts in $g^{″}$ . Finally, the remaining secular problem, being independent of all angle variables. The elements of the short period are obtained from

$\begin{array}{l}{l}_i={l^{\prime }}_i+є{l^{\prime }}_i^1+\frac{1}{2}{є}^2{l^{\prime }}_i^2\left(i=1,2,3,4\right)\\ L_i={L^{\prime }}_i+є{L^{\prime }}_i^1+\frac{1}{2}{є}^2+{L^{\prime }}_i^2\end{array}$ (22)

where ${l^{\prime }}_i^1=\frac{\partial W_1}{\partial {L^{\prime }}_i}$ and ${L^{\prime }}_i^1=-\frac{\partial W_1}{\partial {l^{\prime }}_i}$

${l^{\prime }}_i^2=\frac{\partial W_2}{\partial {L^{\prime }}_i}+\left\{{l^{\prime }}_i^1,W_1\right\}$ and ${L^{\prime }}_i^1=-\frac{\partial W_2}{\partial {l^{\prime }}_i}+\left\{{L^{\prime }}_i^2,W_1\right\}$ (23)

the prime indicate the transformed terms. For the inverse transformations

$l_i^1=-{l^{\prime }}_i^1$

$l_i^2=-{l^{\prime }}_i^2+2\left\{{l^{\prime }}_i^1,W_1\right\}$

$L_i^1=-{L^{\prime }}_i^1$

$L_i^2=-{L^{\prime }}_i^2+2\left\{{L^{\prime }}_i^1,W_1\right\}$ (24)

Similarly for the elements of the long period. The secular terms will obtained from the substitution of the new transformed Hamiltonian $H^{*}$ and $H^{**}$ in Equations (18) and (19) into the Equation (17). These algorithms are actually very simple, but calculating them by hand is laborious, therefore, all computations were carried out by computer program, using the algebraic manipulation language MATHMATICA V10.

4. Third Order Solution

We employ the previous algorithms to obtain analytical expressions for the short-period and long-period terms of the orbital elements; the solution is expressed in terms of the Delaunay variables. We obtain the following results.

4.1. The Short-Period Terms

Zero order: $H_0^{\ast }=-\frac{{\mu }^2}{2L^2}$ (25)

First order: $H_1^{\ast }=A_1\frac{{\theta }_0}{L^3{G}^3}+n_{\odot }{L}_{4s}+n_w{L}_{4w}$ (26)

$\begin{array}{c}{W}_1=\frac{A_1}{{\mu }^2{G}^3}\{{\theta }_0\left(f-l+e\sin f\right)+\frac{{\theta }_1}{2}[\sin 2\left(f+g\right)\begin{array}{c}\\ \end{array}\\ +\frac{e}{3}\sin \left(3f+2g\right)+e\sin \left(f+2g\right)]\}\end{array}$ (27)

Second order: $H_2^{\ast }=H_{20}^{\ast }+H_{21}^{\ast }+H_{22}^{\ast }+H_{23}^{\ast }+H_{24}^{\ast }+H_{25}^{\ast }$

where $H_{20}^{\ast }$ , $H_{21}^{\ast }$ , and $H_{22}^{\ast }$ are coming from Newtonian effect, $H_{23}^{\ast }$ re- presents the relativistic effect, $H_{24}^{\ast }$ is coming from the solar radiation pressure, and $H_{25}^{\ast }$ is coming from GW

$\begin{array}{c}{H}_{20}^{\ast }=\frac{A_1^2}{{\mu }^2}\{\left(-\frac{105}{4}{c}^4-\frac{15}{2}{c}^2+\frac{15}{4}\right){\eta }_{37}+\left(-\frac{15}{4}{c}^4-\frac{27}{2}{c}^2+\frac{15}{4}\right){\eta }_{55}\\ +\left(-27c^4+18c^2-3\right){\eta }_{46}+e^2\left(-\frac{45}{4}{c}^4+24c^2-\frac{3}{2}\right){\eta }_{37}\cos 2g\}\end{array}$ (28)

$H_{21}^{\ast }=A_{21}\left\{3e{\eta }_{35}s\left(1-5c^2\right)\right\}\sin g$ (29)

$\begin{array}{l}{H}_{22}^{\ast }=A_{22}\{\left(105c^4-90c^2+9\right)\left(\frac{5}{2}{\eta }_{37}-\frac{3}{2}{\eta }_{55}\right)\\ \begin{array}{c}\\ \end{array}+e^2\left(-105c^4+120c^2-15\right){\eta }_{37}\cos 2g\}\end{array}$ (30)

$H_{23}^{\ast }=A_{23}\left(3{\eta }_{31}-\frac{15}{8}{\eta }_{40}\right)$ (31)

$H_{24}^{\ast }=\frac{3}{2}{A}_{24}e{L}^2{T}_c$ (32)

$\begin{array}{c}{H}_{25}^{\ast }=A_{25}{L}^4\{\left(1+\frac{1}{2}{e}^2\right)\frac{s^2}{8}[\cos \left(2h+l_{4w}+{\alpha }_1\right)-\cos \left(2h-l_{4w}-{\alpha }_1\right)\\ -A_g\left[\sin \left(2h+l_{4w}+{\alpha }_2\right)+\sin \left(2h-l_{4w}-{\alpha }_2\right)\right]]+\left(\frac{c^2}{16}+\frac{c}{8}-\frac{1}{16}\right)\\ \times \left(-\frac{1}{2}+2e^2\right)\left[\cos \left(2g+2h+l_{4w}+{\alpha }_1\right)+\cos \left(2g+2h-l_{4w}-{\alpha }_1\right)\right]\\ +\left(\frac{c^2}{16}-\frac{c}{8}+\frac{1}{16}\right)\left(-\frac{1}{2}+2e^2\right)[\cos \left(2g-2h+l_{4w}+{\alpha }_1\right)\\ \begin{array}{c}\\ \end{array}+\cos \left(2g-2h-l_{4w}-{\alpha }_1\right)]\}\end{array}$ (33)

where ${\eta }_{ij}=L^{-i}{G}^{-j}$ (i and j are integers), and the corresponding generators are

$W_{20}=\frac{A_1{}^2}{{\mu }^4}\left\{{\displaystyle \sum }_{n=1}^6{\displaystyle \sum }_{m=-1}^2{b}_{1nm}^s\sin \left(nf+2mg\right)+\left(f-l\right){\displaystyle \sum }_{n=0}^3{\displaystyle \sum }_{m=0}^1{b}_{1nm}^c\cos \left(nf+2mg\right)\right\}$ (34)

$W_{21}=\frac{A_{21}}{{\mu }^2}\left\{\underset{n=1}{\overset{5}{{\displaystyle \sum }}}\underset{m=0}{\overset{2}{{\displaystyle \sum }}}{b}_{2nm}^c\cos \left(nf+\left(2m-1\right)g\right)+e\left(f-l\right)\left(15s^3-12s\right){\eta }_{05}\sin g\right\}$ (35)

$\begin{array}{c}{W}_{22}=\frac{A_{22}}{{\mu }^2}\{\underset{n=1}{\overset{7}{{\displaystyle \sum }}}\underset{m=-1}{\overset{2}{{\displaystyle \sum }}}{b}_{3nm}^s\sin \left(nf+2mg\right)+\left(f-l\right)[\left(105c^4-90c^2+9\right)\begin{array}{c}\\ \end{array}\\ \times \left(\frac{5}{2}{\eta }_{07}-\frac{3}{2}{\eta }_{55}\right)+e^2\left(-105c^4+120c^2-15\right){\eta }_{07}\cos 2g]\}\end{array}$ (36)

$W_{23}=\frac{A_{23}}{{\mu }^2}\left\{3{\eta }_{01}\left(f-l\right)-2{\eta }_{10}\left(E-l\right)\right\}$ (37)

$\begin{array}{c}{W}_{24}=-\frac{A_{24}}{{\mu }^2}{L}^5\{\left[\left(1-\frac{1}{2}{e}^2\right)\sin E-\frac{e}{4}\sin 2E\right]T_c\\ +\left[\left(\frac{1}{2}{e}^2-1\right)\cos E+\left(\frac{e}{4}-\frac{e^3}{8}\right)\cos 2E\right]T_s\}\end{array}$ (38)

$\begin{array}{c}{W}_{25}=A_{25}{L}^4\{\left[\left(1+\frac{1}{2}{e}^2\right)\sin E+\frac{1}{2}\sin 2E-2e\sin E\right]\frac{s^2}{8}[\cos \left(2h+l_{4w}+{\alpha }_1\right)\\ -\cos \left(2h-l_{4w}-{\alpha }_1\right)-A_g\left[\sin \left(2h+l_{4w}+{\alpha }_2\right)+\sin \left(2h-l_{4w}-{\alpha }_2\right)\right]]\\ +\left(\frac{c^2}{16}+\frac{c}{8}-\frac{1}{16}\right)\{\left[\left(\frac{1}{2}+2e^2\right)\sin E+\frac{1}{2}\sin 2E-2e\sin E\right]\\ \times \left[\cos \left(2g+2h+l_{4w}+{\alpha }_1\right)+\frac{c}{8}\cos \left(2g+2h-l_{4w}-{\alpha }_1\right)\right]\end{array}$

$\begin{array}{l}-{\left(1-e^2\right)}^{\frac{1}{2}}\left[2e\cos E-\frac{1}{2}\cos E\right][\sin \left(2g+2h+l_{4w}+{\alpha }_1\right)\\ +\sin \left(2g+2h-l_{4w}-{\alpha }_1\right)]+2A_g{\left(1-e^2\right)}^{\frac{1}{2}}\left[2e\cos E-\frac{1}{2}\cos E\right]\\ \times [\cos \left(2g+2h+l_{4w}+{\alpha }_2\right)+\cos \left(2g+2h-l_{4w}-{\alpha }_2\right)\\ \begin{array}{c}\\ \end{array}+\cos \left(2g-2h+l_{4w}+{\alpha }_2\right)+\cos \left(2g-2h-l_{4w}-{\alpha }_2\right)]\}\end{array}$

$\begin{array}{l}+\left(\frac{c^2}{16}-\frac{c}{8}+\frac{1}{16}\right)\{\left[\left(\frac{1}{2}+2e^2\right)\sin E+\frac{1}{2}\sin 2E-2e\sin E\right]\\ \times \left[\cos \left(2g-2h+l_{4w}+{\alpha }_1\right)+\cos \left(2g-2h-l_{4w}-{\alpha }_1\right)\right]\\ +{\left(1-e^2\right)}^{\frac{1}{2}}\left[2e\cos E-\frac{1}{2}\cos E\right][\sin \left(2g-2h+l_{4w}+{\alpha }_1\right)\\ \begin{array}{c}\\ \end{array}+\sin \left(2g-2h-l_{4w}-{\alpha }_1\right)]\}\}\end{array}$ (39)

where E is the eccentric anomaly and the (b) coefficients are coefficients depend on the eccentricity, c and s (i.e. the inclination) and these coefficients is

coming from the expansion of $\frac{a}{r}$ , $\frac{a}{r}\cos f$ , $\frac{a}{r}\sin f$ in the Hamiltonian H

Third order:

$\begin{array}{c}{H}_3^{\ast }=A_3\left\{\left(1-5c^2\right)\left(\frac{3}{4}{\eta }_{35}-\frac{5}{8}{\eta }_{53}\right)+\frac{3}{8}{e}^2\left(c^2-1\right){\eta }_{35}\cos 2g\right\}\\ +\frac{A_1{A}_{21}}{{\mu }^2}\{\left(1-3c^2\right)\left(-45{\eta }_{35}+36{\eta }_{53}-27{\eta }_{44}\right)+[27e^2\left(1-c^2\right){\eta }_{35}\underset{{}_{}^{}}{\overset{}{{}_{}^{}}}\\ +\frac{6\left(c^2-1\right)}{e^2}\left(3{\eta }_{35}+3{\eta }_{53}-9{\eta }_{44}+3{\eta }_7\right)]\}\cos 2g+\frac{A_{24}}{{\mu }^2}{n}_{\odot }{L}_1^5\left(\frac{3e}{2}-\frac{3e^3}{4}\right){T^{\prime }}_s\left(l_{4\odot }\right)\\ -A_{25}{L}^4{n}_w\{[\left(2+e^2\right)\frac{s^2}{8}\left[-A_g\left[\cos \left(2h+l_{4w}+{\alpha }_2\right)+\cos \left(2h-l_{4w}-{\alpha }_2\right)\right]\right]\end{array}$
$\begin{array}{l}+\left(\frac{c^2}{16}+\frac{c}{8}-\frac{1}{16}\right)\{-4e{\left(1-e^2\right)}^{\frac{1}{2}}{\left(1-e^2\right)}^{\frac{1}{2}}[\cos \left(2g+2h+l_{4w}+{\alpha }_1\right)\\ \begin{array}{c}\\ \end{array}+\cos \left(2g+2h-l_{4w}-{\alpha }_1\right)]\}+\left(\frac{c^2}{16}-\frac{c}{8}+\frac{1}{16}\right)\{-4e{\left(1-e^2\right)}^{\frac{1}{2}}\\ \begin{array}{c}\\ \end{array}\times \left[\cos \left(2g-2h+l_{4w}+{\alpha }_1\right)+\cos \left(2g-2h-l_{4w}-{\alpha }_1\right)\right]\}\}\end{array}$ (40)

This order represents the joint effects of the relativistic effect and the solar radiation pressure and GW effects, ${T^{\prime }}_s\left(l_{4\odot }\right)$ is the differentiation of $T_s$ with respect to $l_{4\odot }$ . From the generators $W_1$ and $W_2$ we can determine the elements of the short-period transformation and its inverse according to the Equations (22) and (23).

4.2. The Long-Period Terms

We now proceed to performing a second canonical transformation to eliminate $g^{\prime },h^{\prime },{l^{\prime }}_{4\odot },{l^{\prime }}_{4w}$ (the argument of perigee, the longitude of node and the augmentation variable of the mean motion of the Sun $n_{\odot }$ and the frequency of the gravitational waves $n_w$ ) from the Hamiltonian. This transformation contains the long-periodic perturbations, whereas the secular terms will be derived from the canonical equations of motion and the final Hamiltonian $H^{**}$ . For the simplicity of writing we dropped the primes.

Zero order: $H_0^{\ast \ast }=-\frac{{\mu }^2}{2L^2}$ (41)

First order: $H_1^{\ast \ast }=A_1\frac{{\theta }_0}{L^3{G}^3}$ (42)

Second order:

$\begin{array}{l}{H}_{20}^{\ast \ast }=\frac{A_1^2}{{\mu }^2}\{\left(-\frac{105}{4}{c}^4-\frac{15}{2}{c}^2+\frac{15}{4}\right){\eta }_{37}+\left(-\frac{15}{4}{c}^4-\frac{27}{2}{c}^2+\frac{15}{4}\right){\eta }_{55}\\ \begin{array}{c}\\ \end{array}+\left(-27c^4+18c^2-3\right){\eta }_{46}\}\end{array}$ (43)