## 2 Stationary Orbits
Therefore, we only need to analyze solutions of Equation (5). From the second equation of (5), we see that only one meaningful latitude of the stationary orbit exists, i.e., \(\phi=0\) (\(\phi=\frac{\pi}{2},\pi\) are also roots of \(\sin(2\phi)=0\), but \(\phi=\frac{\pi}{2}\) makes no sense for stationary orbits, and \(\phi=\pi\) corresponds to a stationary orbit which coincides with that of \(\phi=0\)). Introducing the functions \(f_{0}(r)=\frac{\mu}{r^{3}}+\frac{3\mu_{J}R^{2}}{2r^{5}}-\frac{15\mu_{J}R^{4}} {8r^{\prime}}-\eta_{J}^{2}\) and \(f_{1}(r)=\frac{\mu}{r^{3}}-\eta_{J}^{2}\), the variations of \(f_{0}(r)\), \(f_{1}(r)\), with respect to \(r\), can be seen in the following Figure 1.  
It can be proved that \(f_{0}(r)\) vanishes at two real positive values of \(r\), but only one of them is bigger than \(R\). This solution, which is given by \(r_{0}=2.2414R=1.6024\times 10^{8}\) m, is denoted by E0 in Figure 1. If the perturbation due to the non-spherical gravity field is not considered, one can find that the radius for stationary orbits is about \(2.2381R\), which corresponds to the point E1. Therefore, the existence of the \(J_{2}\) and \(J_{4}\) terms increase the altiti  
Figure 1: Variations of \(f_{0}(r)\), \(f_{1}(r)\) as \(r\) increases from \(2.2R\) to \(2.3R\). To study the stability of the stationary orbits under small perturbations, we mainly used the epicyclic theory ([PERSON] [29]). We can denote by the constant \(h\), the value taken by \(r^{2}\cos^{2}\phi\dot{\lambda}\) for some given initial conditions. Assuming that the deviations from the stationary orbits are small and writing \(r=r_{0}(1+\varepsilon)\), Equation (4) can be linearized as follows ([PERSON] [29], Section 11 of Chapter 6):  
\[\left\{ \begin{array}{l}\ddot{\varepsilon}+(3\alpha_{1}+\alpha_{2}) \varepsilon=0,\\ \ddot{\phi}+(\alpha_{1}+\beta_{2})\phi+\beta_{1}=0,\\ \dot{\lambda}-\frac{h}{r_{0}^{2}(1+\varepsilon)^{2}\cos^{2}\phi}=0\end{array} \right. \tag{7}\]  
where \(\alpha_{1}=-\frac{1}{r_{0}}\frac{\partial U}{\partial r}(r_{0},0)\), \(\alpha_{2}=-\frac{2U}{\partial r^{2}}(r_{0},0)\), \(\beta_{1}=-\frac{1}{r_{0}^{2}}\frac{\partial U}{\partial\phi}(r_{0},0)\), and \(\beta_{2}=-\frac{1}{r_{0}^{2}}\frac{\partial^{2}U}{\partial\phi}(r_{0},0)\). By solving Equation (7), the analytical approximate expressions of \(\varepsilon\), \(\phi\), \(\lambda\) can be formulated as [29]  
\[\left\{ \begin{array}{l}\varepsilon=e\cos k_{1}t,\\ \phi=-\frac{\beta_{1}}{\alpha_{1}+\beta_{2}}+i\cos k_{2}t,\\ \lambda=\sqrt{\alpha_{1}}t-\frac{2\sqrt{\alpha_{1}}}{k_{1}}e\sin k_{1}t,\end{array} \right. \tag{8}\]  
where \(k_{1}^{2}=3\alpha_{1}+\alpha_{2}\), \(k_{2}^{2}=\alpha_{1}+\beta_{2}\), \(i\) is the inclination, and \(e\) is the eccentricity of the orbit. Note that the average of \(\lambda\) is \(\sqrt{\alpha_{1}}\). Let us set \(k_{3}=\sqrt{\alpha_{1}}\). The first two equations in (7) or (8) show that the radial and North-South motions are uncoupled. Using the expression (1), we can get a precise form of the three frequencies \(k_{1},k_{2},k_{3}\)[29]:  
\[\left\{ \begin{array}{l}k_{1}{}^{2}=\frac{\mu}{r_{0}}\left[1-\frac{3}{2}I_{2 }\left(\frac{R}{r_{0}}\right)^{2}+\frac{45}{8}I_{4}\left(\frac{R}{r_{0}} \right)^{4}\right],\\ k_{2}{}^{2}=\frac{\mu}{r_{0}}\left[1+\frac{9}{2}I_{2}\left(\frac{R}{r_{0}} \right)^{2}-\frac{75}{8}I_{4}\left(\frac{R}{r_{0}}\right)^{4}\right],\\ k_{3}{}^{2}=\frac{\mu}{r_{0}}\left[1+\frac{3}{2}I_{2}\left(\frac{R}{r_{0}} \right)^{2}-\frac{15}{8}I_{4}\left(\frac{R}{r_{0}}\right)^{4}\right].\end{array} \right. \tag{9}\]  
Therefore, the satellites will oscillate in the radial and North-South directions with frequencies \(k_{1}\) and \(k_{2}\), respectively. Moreover, the mean motion frequency, \(k_{3}\), in the West-East direction is larger than in the Keplerian case. Namely, for a given semi-major axis, the satellite moves faster than the rate expected at that location in the Keplerian case. In the following, we give some numerical examples to illustrate the above characteristics.