# Does apsidal precession rate have any dependence on the argument of periapsis? Perhaps higher order?

Writing this comment got me thinking, which invariably leads to confusion in my case.

@DavidHammen's answer to show that for Juno's near 90° polar orbit the major precession is apsidal; periapsis is moving away from the equator towards the north pole. Part of that answer is:

\begin{aligned} \dot\omega &= \phantom{-} \frac 3 4 J_2 \left(\frac R p\right)^2 n\,(5 \cos^2 i -1) \\ \dot\Omega &= - \frac 3 2 J_2 \left(\frac R p\right)^2 n \cos i \end{aligned} where $$R$$ is the equatorial radius of the planet in question, $$J_2$$ is the planet's second dynamic form, $$p=a(1-e^2)$$ is the semi-latus rectum, $$a$$ is the semi-major axis length of the orbit, $$e$$ is the eccentricity of the orbit, $$n$$ is the mean motion, and $$i$$ is the inclination of the orbit.

The JPL News item NASA’s Juno Mission Expands Into the Future is accompanied by an interesting graphic shown below using the steady advance of periapsis towards the pole both as a "calendar" to map out events in mission history and to show when this advance will naturally bring Juno's orbit to intersection with some of Jupiter's moons; flyby(s) of which are an important component of Juno's extended mission.

What surprises me is that the equation suggests that for a purely polar orbit ($$\i=\pi/2$$) around a pure ellipsoid the rate of apsidal precession $$\dot{\omega}$$ seems to be constant even though the orbit spends different times near the poles and equator and at different distances as it precesses around.

Question: Is this true exactly for these conditions? Is the rate of apsidal precession $$\dot{\omega}$$ really constant and independent of the argument of periapsis $$\omega$$, or is this just a first order term and there are smaller higher order corrections to $$\dot{\omega}$$ that depend on $$\omega$$ explicitly?

I am sure that there are $$\dot{e}$$ terms so I'm not asking about long term evolution of a given orbit. Instead suppose one starts a set of otherwise identical orbits that differ only in initial $$\omega_0$$, would $$\dot{\omega}$$ be exactly the same for all of them?

below: From Juno's Mission Goes On

$$\newcommand{\d}{\partial} \newcommand{\r}{\mathbf{r}}$$

Yes, there are many other terms involving higher powers of $$\omega$$, though the powers of $$\omega$$ are more often encountered as $$\cos 2 \omega$$, $$\sin 3 \omega$$, etc., and the size of most of them is usually ignorable, because they are multiplied by things which are small, like $$J_2$$ cubed or $$J_5$$. There are recipes for generating however many terms you want to use, with arbitrarily high powers of all the orbital elements, but they are complicated to explain. For the first part of the derivation, you might read Robert Battin's An Introduction to the Mathematics and Methods of Astrodynamics (1999), where this material is most of Chapter 10. It's also in pretty much every other solid astrodynamics text you might like, though all with slightly different notation. The second part is not covered well in Battin, Vallado, or other such books I've read, with the exception of William Kaula's Theory of Satellite Geodesy, so I recommend you try to find a copy of that.

I tried to boil down the derivation to something I could post here, but I couldn't manage it. The algebra is just too dense. The main idea is to write the full dynamics as the simple Keplerian two-body solution plus something else, where all of the extra accelerations are known and the task is to find the change to the two-body motion which matches those perturbations. This is called variation of parameters, which is a general method for solving certain types of differential equations, invented by Euler and extended by Lagrange and Gauss, specifically to treat this problem in determining the orbits of Jupiter and Saturn. They figured out the equations describing how the orbital elements change with time, based on all of the things other than the simple Keplerian two-body point-mass that you want to consider. That results in Lagrange's planetary equations. Of the six, the one you asked about is

$$\frac{d \omega}{d t}=\frac{-\cos i}{\sqrt{\mu a (1-e^2)} \sin i} \frac{\d R}{\d i} + \frac{\sqrt{1-e^2}}{e \sqrt{\mu a}} \frac{\d R}{\d e}$$

in Lagrange's form, where $$R$$ is the ''disturbing function'', whose derivatives give the accelerations we want to study. In the case where the disturbance is caused by a third gravitating body, such as the sun or moon, $$R$$ equals

$$\frac{Gm}{\rho}\left(\frac{\rho}{d} - \frac{r \cos \alpha}{\rho}\right)$$

where $$m$$ is the third body's mass, $$\rho$$ is its distance from the central body, $$d$$ is its distance from the orbiting body, and $$\alpha$$ is the angle between the vectors $$\r$$ and $$\rho$$. Finding an $$R$$ that gives rise to the acceleration you wish to study can be a difficulty, so Gauss's alternate form uses the accelerations directly, but I'm not going to introduce it here. Instead, let's look at where we get $$R$$ and how we take its derivatives. If we're looking at the effect of the central body not being a point mass, then it is traditional to write $$R$$ as

$$\frac{GM}{r} \sum_{n=1}^\infty\left(\frac{a}{r}\right)^n \sum_{m=0}^n \bar{P}_{nm}(\sin\phi)[\bar{C}_{nm}\cos(m\lambda)+\bar{S}_{nm}\sin(m\lambda)]$$

as described in this answer. Keep in mind that we didn't have to do it this way. We chose to do it this way, because this is actually a lot easier than many other ways we could have chosen! The coefficients $$\bar{C}_{nm}$$ and $$\bar{S}_{nm}$$ are defined as the result of integrating appropriate weighting functions times the gravitational potential, with domain the mass distribution of the central body, but they are actually measured by comparing this formula to the observed motion of satellites. This process is the "satellite geodesy" of Kaula's title.

What I called "the second part" above is the next step, which is actually taking the derivatives of the disturbing function, and then fiddling around until we come up with some approximation we can solve. Much of the work at this point is classical mechanics at the level of Herbert Goldstein's book, making a bunch of "canonical" transformations, which really means changing variables until you transform the Hamiltonian into something with only trivial solutions. This is the origin of mean element theories like Brouwer's work that began with Delaunay variables and went on from there. The other half, which I've only seen worked out to this level of detail in Kaula, is to rewrite various pieces of the disturbing function as things whose derivatives can more easily be taken. It so happens that the sine of the latitude $$\phi$$ equals $$\sin i \sin(\omega + f)$$, where $$f$$ is the true anomaly, and $$\sin \phi$$ is the argument of the associated Legendre polynomials! Furthermore, one can rewrite the terms in sin and $$\cos m \lambda$$ (longitude) as polynomials in sin and cos of ($$m(\Omega-\theta)+c(\omega+f)$$), where c is an integer depending on the power of cos or sin, and $$\theta$$ is Greenwich hour angle. This means every term in the expansion contains $$\omega$$ explicitly, except for the few that are symmetric about the rotational axis of the spherical coordinate system.

================================================================

As an aside, something I find even more interesting happens at the so-called critical inclination: if $$\cos^2 i = 1/5$$, which happens at approximately 63.435 degrees, then by that formula $$\dot{\omega}$$ is exactly zero, no matter what value $$J_2$$ takes, and thus the argument of perigee does not change, no matter what body is being orbited.

This is used practically for the so-called Molniya orbit, a highly-eccentric (typically $$e$$ is 0.6 to 0.75) orbit with argument of perigee at 270 degrees (near the south pole), so that the apogee is near the north pole, and could thus give much better geographic coverage of the Soviet Union than could a (necessarily near-equatorial) geostationary orbit. High eccentricity exploits Kepler's equal area equal time law to arrange it so that the satellite spends most of its 12-hour orbit moving very slowly across the northern sky, and then zooms very quickly around the southern half of the globe before popping back up over the north where it wants to be and slowing down again. In such an orbit, you need $$\dot{\omega}$$ to vanish for the whole concept to work, or else your satellite that starts out dwelling over the northern hemisphere will start rotating around so that apogee becomes first equatorial and then over the southern hemisphere, where you don't need the communication coverage.

Now, of course $$\dot{\omega}$$ is never exactly zero, so these vehicles do have to do their own kind of stationkeeping, but the forces on them that they maneuver to cancel are very different than those found in other kinds of named orbits. The exact nature of this oddity has attracted a great deal of research interest over the years, some with a rather frightening degree of mathematical sophistication. To give a flavor of the sort of thing one may find, here are links to three very different papers that touch on interesting parts of the issue.

"the critical inclination in the main problem of artificial satellite theory is an intrinsic singularity. Its significance stems from two geometric events in the reduced phase space on the manifolds of constant polar angular momentum and constant Delaunay action. In the neighborhood of the critical inclination, along the family of circular orbits, there appear two Hopf bifurcations, to each of which there converge two families of orbits..."

and gets harder from there.

• @uhoh copied in some equations, and added a tearline between my direct answer to you and my additional ramblings. Battin's book has a fair amount of material on how to derive these expressions for yourself (Chapter 10, pages 471-504), but he leaves equations that I wanted to post here as homework problems, as is sadly typical. Vallado Chapter 9 covers the same material with many worked examples, but none of them are quite right to go here. Amusingly, the very next topic is SGP4, because setting up alternate equations so that these things cancel is exactly what those mean elements are for. Feb 4, 2021 at 14:22
• – uhoh
Feb 7, 2021 at 1:14
• @uhoh expanded answer is in work. I've written much of it already, but some of the final pieces are quite complicated. hopefully by tomorrow... Feb 7, 2021 at 1:56
• @uhoh I also didn't realize you had made edits to the references. Now that you mention it, I have no problem linking abstracts instead --- my "rejection" was just an accidental side effect of both of us editing the post at the same time, sorry. :( Feb 7, 2021 at 2:06
• @uhoh that was just a failure to paste in the top couple lines of the $\LaTeX$. It's working now, because I reinserted the missing \newcommand{\d}{\partial} \newcommand{\r}{\mathbf{r}} Feb 11, 2021 at 2:15