Real planets are not points: multipoles, Legendre, and secular perturbation

Idea. The recurrence builds the polynomials; the generating function at a special \(x\) collapses to an elementary function whose Taylor series is known; Legendre's equation degenerates at \(x = 1\) into an equation for the first derivative.

1. Bonnet at \(n = 4\): \(5P_5 = 9xP_4 - 4P_3\), so

\[ 5P_5 = \tfrac98(35x^5 - 30x^3 + 3x) - 2(5x^3 - 3x) = \tfrac18\big(315x^5 - 270x^3 + 27x - 80x^3 + 48x\big) = \tfrac18(315x^5 - 350x^3 + 75x), \] \[ P_5(x) = \tfrac18\big(63x^5 - 70x^3 + 15x\big). \]

\(P_5(1) = (63 - 70 + 15)/8 = 1\) ✓, and only odd powers appear ✓. The leading coefficient \(63/8\) equals \(\binom{10}{5}2^{-5} = 252/32 = 63/8\) ✓.

2. \(g(0,t) = (1+t^2)^{-1/2} = \sum_m c_m(-1)^mt^{2m}\), so

\[ P_{2m}(0) = (-1)^m\binom{2m}{m}4^{-m},\qquad P_{2m+1}(0) = 0 . \]

\(m = 1\): \(-2/4 = -\tfrac12\), and \(P_2(0) = \tfrac12(0-1) = -\tfrac12\) ✓. \(m = 2\): \(6/16 = \tfrac38\), and \(P_4(0) = \tfrac38\) ✓. \(m = 3\): \(P_6(0) = -20/64 = -\tfrac5{16}\).

3. \(P_n\) is a polynomial, so \(P_n''(1)\) is finite and \((1-x^2)P_n''\) vanishes at \(x = 1\). The equation \((1-x^2)P_n'' - 2xP_n' + n(n+1)P_n = 0\) at \(x = 1\) reads \(-2P_n'(1) + n(n+1)P_n(1) = 0\), and \(P_n(1) = 1\), so \(P_n'(1) = \tfrac12n(n+1)\). On \(P_5\): \(P_5' = \tfrac18(315x^4 - 210x^2 + 15)\), so \(P_5'(1) = 120/8 = 15 = \tfrac12\cdot5\cdot6\) ✓.

4. Over the pole \(\cos\theta = 1\) and \(P_2 = 1\); over the equator \(\cos\theta = 0\) and \(P_2 = -\tfrac12\). The ratio is \(-2\). Since \(J_2 \gt 0\), \(\Phi_2\) is positive over the pole and negative over the equator: the equator is the more strongly bound — the potential well is deepened where the bulge is.

Check. The values \(|P_{2m}(0)| = \tfrac12, \tfrac38, \tfrac5{16}\) are all at most \(1\), as the bound \(|P_n| \le 1\) requires, and they decrease, as they must if \(\sum P_{2m}(0)t^{2m}\) is to converge at every \(|t| \lt 1\). Part 3 on \(P_2\): \(P_2' = 3x\), \(P_2'(1) = 3 = \tfrac12\cdot2\cdot3\) ✓.

Idea. Polynomials of degree \(\lt n\) form an \(n\)-dimensional space, and \(P_0,\dots,P_{n-1}\) are \(n\) elements of it with distinct degrees, hence a basis. Orthogonality then does the rest.

1. The matrix expressing \(P_0,\dots,P_{n-1}\) in the monomial basis \(1,\dots,x^{n-1}\) is triangular with nonzero diagonal (the leading coefficients), hence invertible. So \(q = \sum_{k \lt n}b_kP_k\) and \(\int qP_n = \sum_{k\lt n}b_k\int P_kP_n = 0\) by orthogonality. \(\square\)

2. \(x^2 = \tfrac13P_0 + \tfrac23P_2\). From \(P_3 = \tfrac12(5x^3 - 3x)\), \(x^3 = \tfrac25P_3 + \tfrac35x = \tfrac35P_1 + \tfrac25P_3\). From \(P_4 = \tfrac18(35x^4 - 30x^2 + 3)\),

\[ x^4 = \tfrac8{35}P_4 + \tfrac67x^2 - \tfrac3{35} = \tfrac8{35}P_4 + \tfrac47P_2 + \big(\tfrac27 - \tfrac3{35}\big)P_0 = \tfrac15P_0 + \tfrac47P_2 + \tfrac8{35}P_4 . \]

Check the \(P_0\) coefficient by orthogonality: \(a_0 = \tfrac12\int_{-1}^1x^4\,\mathrm dx = \tfrac12\cdot\tfrac25 = \tfrac15\) ✓. Only polynomials of the parity of \(x^k\) appear, as they must.

3. \(\cos^2\theta = \tfrac13 + \tfrac23P_2(\cos\theta)\), so

\[ r = R\big(1 - \tfrac f3\big) - \tfrac{2f}{3}R\,P_2(\cos\theta), \qquad R_0 = R\big(1-\tfrac f3\big) = 6371.009\ \mathrm{km},\quad R_2 = -\tfrac{2f}3R = -14.256\ \mathrm{km}. \]

\(R_0\) is the average of \(r\) over the sphere of directions, because \(\int_{-1}^1P_2 = 0\) kills the second term under the average in \(\cos\theta\) (and \(\mathrm d\Omega = \sin\theta\,\mathrm d\theta\,\mathrm d\phi\) is uniform in \(\cos\theta\), unit 05). It is \((2R + c)/3\) with \(c = R(1-f)\): the geodesists' mean radius, published as \(6371.0088\) km ✓.

Check. At \(x = 1\) every expansion must give \(1\), since \(P_n(1) = 1\): \(\tfrac15 + \tfrac47 + \tfrac8{35} = \tfrac{7+20+8}{35} = 1\) ✓, and \(\tfrac35 + \tfrac25 = 1\) ✓. In part 3, at the pole \(r = R_0 + R_2 = 6356.752\) km \(= R(1-f)\) ✓, and at the equator \(r = R_0 - \tfrac12R_2 = 6378.137\) km ✓.

Idea. Rodrigues' formula puts \(P_n\) in the form "\(n\)-th derivative of a function vanishing to order \(n\) at both ends", and that is precisely what integration by parts \(n\) times needs to transfer all the derivatives to the other factor with no boundary terms.

1. By Leibniz, \(D^k[(x-1)^n(x+1)^n] = \sum_j\binom kj D^j(x-1)^n\,D^{k-j}(x+1)^n\). With \(k \lt n\), each \(D^j(x-1)^n\) is a multiple of \((x-1)^{n-j}\) with \(n - j \ge 1\), so every term vanishes at \(x = 1\). The same with the roles of the factors swapped at \(x = -1\).

2. Integrate by parts \(k\) times, for \(k = 1,\dots,n\):

\[ \int_{-1}^1q\,D^nv_n = (-1)^k\int_{-1}^1(D^kq)\,D^{n-k}v_n , \]

each boundary term \(\big[D^{j}q\;D^{n-1-j}v_n\big]_{-1}^{1}\) vanishing by part 1 since \(n-1-j \lt n\). At \(k = n\), \(D^nq = 0\) because \(\deg q \lt n\). For \(m \lt n\), take \(q = P_m\). \(\square\)

3. \(P_n = \frac{(2n)!}{2^n(n!)^2}x^n + (\text{degree} \lt n)\) — the leading coefficient of \(D^n x^{2n}/(2^nn!)\). By part 2 only the top term of the first factor survives:

\[ \int P_n^2 = \frac{(2n)!}{2^n(n!)^2}\cdot\frac1{2^nn!}\int x^n\,D^nv_n = \frac{(2n)!}{4^n(n!)^3}(-1)^n\int (D^nx^n)\,v_n = \frac{(2n)!}{4^n(n!)^2}\,I_n , \]

using \(D^nx^n = n!\) and \((-1)^nv_n = (1-x^2)^n\). For \(I_n\), integrate by parts with \(u = (1-x^2)^n\), \(\mathrm dv = \mathrm dx\):

\[ I_n = \big[x(1-x^2)^n\big]_{-1}^1 + 2n\int x^2(1-x^2)^{n-1} = 2n\int\big(1 - (1-x^2)\big)(1-x^2)^{n-1} = 2n(I_{n-1} - I_n), \]

so \(I_n = \frac{2n}{2n+1}I_{n-1}\) with \(I_0 = 2\), giving \(I_n = 2\cdot\frac{2\cdot4\cdots(2n)}{3\cdot5\cdots(2n+1)} = \frac{2\cdot4^n(n!)^2}{(2n+1)!}\). Therefore

\[ \int_{-1}^1P_n^2 = \frac{(2n)!}{4^n(n!)^2}\cdot\frac{2\cdot4^n(n!)^2}{(2n+1)!} = \frac{2}{2n+1}. \qquad\square \]

Check. \(n = 1\): \(\int x^2 = \tfrac23\) ✓. \(n = 2\): \(\tfrac14\int(9x^4 - 6x^2 + 1) = \tfrac14\big(\tfrac{18}5 - 4 + 2\big) = \tfrac25\) ✓. And \(I_1 = \int(1-x^2) = \tfrac43 = \tfrac23I_0\) ✓.

Another way

The deck's route: Legendre's equation in self-adjoint form and two integrations by parts give orthogonality, and squaring the generating function gives the norm. The two proofs use disjoint facts — the equation there, the explicit formula here — which is a genuine check that the three definitions of \(P_n\) agree.

Idea. Everything is along one line, so \(Q\) is a multiple of the single trace-free form built from \(\hat{\mathbf e}\); the multiple is a moment of inertia, and it is the reduced mass that appears, exactly as in unit 08's two-body reduction.

1. Dipole: \(m_1\vec r_1 + m_2\vec r_2 = \big(-\frac{m_1m_2}M + \frac{m_2m_1}M\big)d\,\hat{\mathbf e} = 0\). Quadrupole: with \(\vec r_k = s_kd\,\hat{\mathbf e}\), \(s_1 = -m_2/M\), \(s_2 = m_1/M\),

\[ Q(\vec u,\vec v) = \Big(\sum_km_ks_k^2\Big)d^2\big[3(\vec u\cdot\hat{\mathbf e})(\vec v\cdot\hat{\mathbf e}) - \vec u\cdot\vec v\big], \qquad \sum_km_ks_k^2 = \frac{m_1m_2^2 + m_2m_1^2}{M^2} = \mu . \]

2. \(Q(\hat{\mathbf e},\cdot)\): the matrix is \(\mu d^2(3\hat{\mathbf e}\hat{\mathbf e}^\top - \mathbb 1)\). Eigenvector \(\hat{\mathbf e}\), eigenvalue \(2\mu d^2\); any vector perpendicular to \(\hat{\mathbf e}\), eigenvalue \(-\mu d^2\), twice. Trace \(2 - 1 - 1 = 0\) ✓.

3. \(Q(\hat{\mathbf r},\hat{\mathbf r}) = \mu d^2(3\cos^2\psi - 1) = 2\mu d^2P_2(\cos\psi)\), \(\psi\) the angle from \(\hat{\mathbf e}\), so

\[ \Phi_2 = -\frac{G\,Q(\hat{\mathbf r},\hat{\mathbf r})}{2r^3} = -\frac{G\mu d^2}{r^3}P_2(\cos\psi). \]

It is most negative along the line of the pair — a prolate distribution, the opposite sign from an oblate planet's \(J_2\) term.

4. \(\mu = \frac{5.9722\times10^{24}\cdot7.342\times10^{22}}{6.0456\times10^{24}} = 7.253\times10^{22}\) kg, so \(\mu d^2 = 7.253\times10^{22}\times(3.844\times10^8)^2 = 1.072\times10^{40}\ \mathrm{kg\,m^2}\). The Earth's \(C - A = 1.08263\times10^{-3}\times5.9722\times10^{24}\times(6.378137\times10^6)^2 = 2.630\times10^{35}\ \mathrm{kg\,m^2}\). The pair's is larger by a factor \(4.07\times10^4\): from far away the Earth's bulge is invisible beside the Moon's orbit.

Check. Dimensions: \(\mu d^2\) is \(\mathrm{kg\,m^2}\), like \(C - A\) ✓. Limit \(m_2\to0\): \(\mu\to0\) and a lone point mass at the origin has no quadrupole ✓. Equal masses: \(\mu = M/4\), and each mass at distance \(d/2\) contributes \(\frac M2\cdot\frac{d^2}4\), total \(Md^2/4\) ✓.

Idea. Each moment is an integral of a polynomial in \(\vec r\,{}'\); shifting \(\vec r\,{}'\) mixes in the lower moments, and only the lowest one has nothing below it to mix in.

1. \(\int\rho\,(\vec r\,{}' - \vec c) = \vec d - M\vec c\).

2. Expanding,

\[ 3\big(\vec u\cdot\vec r\,{}' - \vec u\cdot\vec c\big)\big(\vec v\cdot\vec r\,{}' - \vec v\cdot\vec c\big) = 3(\vec u\cdot\vec r\,{}')(\vec v\cdot\vec r\,{}') - 3(\vec u\cdot\vec r\,{}')(\vec v\cdot\vec c) - 3(\vec u\cdot\vec c)(\vec v\cdot\vec r\,{}') + 3(\vec u\cdot\vec c)(\vec v\cdot\vec c). \]

Integrating against \(\rho\) and subtracting \(\big(\int\rho r'^2 - 2\vec c\cdot\vec d + Mc^2\big)\vec u\cdot\vec v\) gives the stated formula term by term. Trace: put \(\vec u = \vec v = \hat e_i\) and sum over \(i\). The first term gives \(0\); the bracket gives \(-3\cdot2(\vec c\cdot\vec d)\); the next \(+2\cdot3(\vec c\cdot\vec d)\); the last \(M(3c^2 - 3c^2)\). Total \(0\) ✓.

3. If \(M \ne 0\), \(M\) is lowest and trivially origin-independent; the new dipole \(\vec d - M\vec c\) vanishes if and only if \(\vec c = \vec d/M\), the centre of mass — exactly one origin. If \(M = 0\) and \(\vec d \ne 0\), part 1 gives \(\vec d'' = \vec d\): the dipole is the lowest moment and does not move. If \(M = 0\) and \(\vec d = 0\), every correction term in part 2 carries a factor \(\vec d\) or \(M\), so \(Q'' = Q\). \(\square\)

Check. Take the origin at the centre of mass (\(\vec d = 0\)): the formula reduces to the parallel-axis shift \(Q'' = Q + M(3\vec c\,\vec c^{\,\top} - c^2\mathbb 1)\). For a single point mass \(M\) at \(O\) (so \(Q = 0\)), seen from \(O + \vec c\), it sits at \(-\vec c\) and the definition gives directly \(M[3(\vec u\cdot\vec c)(\vec v\cdot\vec c) - c^2\,\vec u\cdot\vec v]\) ✓. Dimensions: every term is mass times length squared ✓. This is why gravity's expansion always starts monopole-then-quadrupole: mass is positive, so the dipole can always be removed, and it is also why P4's answer did not depend on anything but \(\mu\) and \(d\).

Idea. A linear map carries the ball to the ellipsoid and multiplies volumes by its determinant, so every second moment of the ellipsoid is a second moment of the ball, scaled. \(J_2\) then measures only the bulge, and anything spherical dilutes it.

1. With \(x = Au\), \(A = \operatorname{diag}(a_1,a_2,a_3)\), \(\mathrm dV_x = a_1a_2a_3\,\mathrm dV_u\) (unit 05):

\[ \int_Ex_3^2\,\mathrm dV = a_1a_2a_3\,a_3^2\int_Bu_3^2\,\mathrm dV = a_1a_2a_3\,a_3^2\cdot\tfrac13\int_0^1 r^2\,4\pi r^2\,\mathrm dr = \tfrac{4\pi}{15}a_1a_2a_3\,a_3^2 , \]

using \(\int_Bu_1^2 = \int_Bu_2^2 = \int_Bu_3^2\) by symmetry of the ball. With \(\rho = M/(\tfrac43\pi a_1a_2a_3)\), \(\rho\int_Ex_3^2 = \tfrac M5a_3^2\), and likewise for \(x_2\), so \(I_{11} = \tfrac M5(a_2^2 + a_3^2)\), and cyclically.

2. \(C = I_{33} = \tfrac25MR^2\), \(A = I_{11} = \tfrac15M(R^2 + c^2)\) with \(c = R(1-f)\):

\[ J_2 = \frac{C - A}{MR^2} = \frac{R^2 - c^2}{5R^2} = \frac{1 - (1-f)^2}{5} = \frac{2f - f^2}{5} = 1.33888\times10^{-3} \]

for \(f = 1/298.257223563\). The measured \(1.08263\times10^{-3}\) is smaller by a factor \(0.8086\): a uniform Earth would have \(23.7\%\) too much quadrupole.

3. The spherical part adds \(M_{\mathrm c}\) to the mass and equal amounts to \(A\) and \(C\), so \(C - A = \tfrac15M_{\mathrm s}(R^2 - c^2)\) is unchanged while \(M = M_{\mathrm s} + M_{\mathrm c}\):

\[ J_2 = \frac{M_{\mathrm s}}{M}\cdot\frac{2f - f^2}{5}, \qquad \frac{M_{\mathrm s}}{M} = \frac{1.08263}{1.33888} = 0.8086,\qquad \frac{M_{\mathrm c}}{M} = 0.191 . \]

Check. Dimensions: \(J_2\) is a ratio of squared lengths ✓. Limits: \(f \to 0\) gives \(J_2 \to 0\), a sphere ✓; \(M_{\mathrm c}\to0\) recovers part 2 ✓; \(M_{\mathrm c}\to\infty\) gives \(J_2\to0\), a point mass with a massless envelope, whose bulge weighs nothing ✓. A number to compare. \(M_{\mathrm c}\) is the core's mass in excess of the surrounding density. For the real Earth — core radius \(3485\) km, core mass \(1.94\times10^{24}\) kg, so mean core density \(10.9\) and mean mantle density \(4.45\ \mathrm{g\,cm^{-3}}\) — that excess is \((10.94 - 4.45)\times10^3\times\frac43\pi(3.485\times10^6)^3 / M_\oplus = 0.193\) of the mass. Agreement to a percent is partly luck: the model ignores the density rising continuously through the mantle and the core's own flattening. But it shows how \(J_2\) and \(f\) together weigh an interior nobody can visit.

Idea. Set the secular node rate equal to the Sun's apparent eastward motion and solve for the one unknown.

1. \(\dot\Omega = 2\pi/31556926\ \mathrm s = 1.99106\times10^{-7}\ \mathrm{s^{-1}}\).

• Sentinel-2: \(a = 7164.137\) km, \(n = \sqrt{GM/a^3} = 1.041173\times10^{-3}\ \mathrm{s^{-1}}\), \(\cos i = -\frac{2(1.99106\times10^{-7})(7164.137)^2}{3(1.08263\times10^{-3})(6378.137)^2(1.041173\times10^{-3})} = -0.148570\), so \(i = 98.54^\circ\). Published \(98.62^\circ\): the prediction is \(0.08^\circ\) low.
• Landsat 8: \(a = 7083.137\) km, \(n = 1.059084\times10^{-3}\ \mathrm{s^{-1}}\), \(\cos i = -0.142773\), \(i = 98.21^\circ\). Published \(98.2^\circ\): agreement to the digits given.

2. \(\cos i = -\frac{2\dot\Omega}{3J_2R^2}\cdot\frac{a^2}{n}\) and \(n \propto a^{-3/2}\), so \(\cos i \propto -a^{7/2}\). Higher orbits need a larger \(|\cos i|\), so \(i\) moves further past \(90^\circ\): it grows with altitude, as the two satellites show ✓.

3. Sun-synchronism needs \(\dot\Omega \gt 0\) — the node must move eastward with the Sun — and every factor in \(-\tfrac32J_2(R/a)^2n\cos i/(1-e^2)^2\) other than \(\cos i\) is positive. So \(\cos i \lt 0\): the orbit must be retrograde.

Check. Dimensions: \(\dot\Omega a^2/(R^2n)\) is \(\mathrm{s^{-1}}/\mathrm{s^{-1}}\), dimensionless ✓. Scaling check with the numbers: \((7164.137/7083.137)^{3.5} = 1.0406\), and \(0.148570/0.142773 = 1.0406\) ✓. The residual \(0.08^\circ\) for Sentinel-2 is the size one expects from what \(J_2\) alone omits — the deck's notes find \(0.04^\circ\) from the mean-versus-osculating distinction alone, and "mean altitude" over an oblate Earth is itself defined to a few km.

Idea. On a circular orbit the right-hand sides are trigonometric polynomials in the fast angle \(u\). The constant term integrates to a ramp; the rest integrates to bounded oscillations. That is the whole distinction between secular and periodic, done without the averaging theorem.

1. Put \(k = \tfrac32J_2(R/a)^2\). Dividing by \(L\sin i\) and by \(-L\):

\[ \dot\Omega = -k\,n\cos i\,(1 - \cos2u), \qquad \dot i = -k\,n\sin i\cos i\,\sin2u . \]

At first order \(i\) and \(a\) on the right may be held at their initial values, and \(u = nt + u_0\). Integrating,

\[ \Omega(t) = \Omega_0 - k\cos i\Big[nt - \tfrac12\big(\sin2u - \sin2u_0\big)\Big], \qquad i(t) = i_0 + \tfrac12k\sin i\cos i\,\big(\cos2u - \cos2u_0\big). \]

2. \(\Omega\): secular part \(-k\,n\cos i\,t\), the deck's \(\dot\Omega\); periodic part \(\tfrac12k\cos i\,\sin2u\), amplitude \(\tfrac12k|\cos i| = \tfrac34J_2(R/a)^2|\cos i|\), twice per orbit. \(i\): no secular part at all; periodic amplitude \(\tfrac12k\sin i\,|\cos i|\). Numbers, with \((R/a)^2 = 0.789522\) and \(\cos i = -0.149589\):

• \(\Omega\) wobble: \(\tfrac34(1.08263\times10^{-3})(0.789522)(0.149589) = 9.590\times10^{-5}\) rad \(= 0.00549^\circ\);
• \(i\) wobble: that times \(\sin i = 0.98875\), \(9.482\times10^{-5}\) rad \(= 0.00543^\circ\);
• \(\Omega\) drift per orbit: \(2\pi k|\cos i| = 1.2051\times10^{-3}\) rad \(= 0.0690^\circ\), eastward since \(\cos i \lt 0\).

3. Drift after time \(t\) is \(k|\cos i|\,nt\); the wobble amplitude is \(\tfrac12k|\cos i|\). They are equal when \(nt = \tfrac12\), i.e. after \(\frac{1}{4\pi} = 0.0796\) of an orbit — for every orbit, every inclination, every \(J_2\), because both are the same prefactor \(k\cos i\) and the ratio is pure geometry of \(\sin^2u\). After a hundred orbits the drift is \(6.90^\circ\) and the wobble still \(0.0055^\circ\).

Check. Averaging \(\dot\Omega\) over \(u\) gives \(-kn\cos i\) exactly, the deck's formula at \(e = 0\) ✓. \(i\) has no secular term, in agreement with \(L_z = L\cos i\) being exactly conserved for an axisymmetric potential ✓. At \(i = 90^\circ\) everything vanishes — a polar orbit feels no net twist ✓. Dimensions: \(k\) is dimensionless, so amplitudes are angles and \(kn\) is a rate ✓.

Idea. \(\vec A\) is conserved by the Kepler flow, so the deck's lemma gives its rate of change from the perturbation alone. Its component across itself turns the ellipse; its component along itself changes the eccentricity. Averaging with the correct weight keeps the first and kills the second.

1. The deck's \(\vec a_{\mathrm p} = \frac{3GMJ_2R^2}{2r^5}\big[(5z^2/r^2 - 1)\vec r - 2z\hat{\mathbf z}\big]\) at \(z = 0\) is \(-\frac{3GMJ_2R^2}{2r^5}\vec r = -\frac{3GMJ_2R^2}{2r^4}\hat{\mathbf r}\).

2. \(\vec A = \vec v\times(\vec r\times\vec v) - GM\hat{\mathbf r}\). Replacing \(\vec v\) by \(\vec v + \varepsilon\vec a\) and keeping the first order in \(\varepsilon\): \(\vec a\times(\vec r\times\vec v) + \vec v\times(\vec r\times\vec a)\), which is the claim. The \(GM\hat{\mathbf r}\) term has no \(\vec v\) in it. For \(\vec a = f\hat{\mathbf r}\) the second term vanishes, and the lemma gives \(\dot{\vec A} = f\,\hat{\mathbf r}\times\vec L\).

3. In the orbit plane, with \(\vec L = L\hat{\mathbf z}\) and the right-handed triple \((\hat{\mathbf r},\hat{\boldsymbol\theta},\hat{\mathbf z})\), \(\hat{\mathbf r}\times\hat{\mathbf z} = -\hat{\boldsymbol\theta}\), so \(\dot{\vec A} = -fL\hat{\boldsymbol\theta} = fL(\sin\nu, -\cos\nu)\). With \(\vec A = GMe\,\hat{\mathbf x}\) at this instant,

\[ \dot e = \frac{\dot A_x}{GM} = \frac{fL\sin\nu}{GM}, \qquad \dot\varpi = \frac{\dot A_y}{GMe} = -\frac{fL\cos\nu}{GMe}. \]

Average in time, \(\langle X\rangle = \frac1T\int_0^{2\pi}X\,\frac{r^2}{L}\mathrm d\nu\). Since \(fr^2 = -\tfrac32GMJ_2R^2(1+e\cos\nu)^2/p^2\) is even in \(\nu\), \(\langle\dot e\rangle \propto \int(\text{even})\sin\nu = 0\): no secular change in \(e\). And

\[ \langle\dot\varpi\rangle = \frac{3J_2R^2}{2Tp^2e}\int_0^{2\pi}(1+e\cos\nu)^2\cos\nu\,\mathrm d\nu = \frac{3J_2R^2}{2Tp^2e}\cdot2\pi e = \frac32\,n\,J_2\Big(\frac Rp\Big)^2 , \]

using \(\int\cos\nu = \int\cos^3\nu = 0\), \(\int2e\cos^2\nu = 2\pi e\), and \(2\pi/T = n\). It is positive: the apse line advances with the motion.

4. At \(i = 0\): \(\dot\Omega + \dot\omega = \tfrac{J_2(R/a)^2n}{(1-e^2)^2}\big[-\tfrac32 + \tfrac34(5-1)\big] = \tfrac32\,n\,J_2\,\frac{R^2}{a^2(1-e^2)^2} = \tfrac32nJ_2(R/p)^2\) ✓. This confirms the deck's \((1-e^2)^{-2}\) at every \(e\), in the equatorial case.

5. \(\dot\omega = 0\) when \(\cos^2 i = \tfrac15\): \(i = 63.4349^\circ\) or \(116.5651^\circ\). For the Molniya orbit, \(n = 1.45842\times10^{-4}\ \mathrm{s^{-1}}\), \((R/a)^2 = 0.057660\), \((1-e^2)^{-2} = 4.886\), so \(\tfrac34J_2(R/a)^2n(1-e^2)^{-2} = 3.336\times10^{-8}\ \mathrm{s^{-1}} = 0.16515^\circ\)/day. Then:

• \(i = 60^\circ\): \(5\cos^2i - 1 = 0.25\), \(\dot\omega = 0.0413^\circ\)/day;
• \(i = 63.4^\circ\): \(5\cos^2i - 1 = 0.002441\), \(\dot\omega = 4.03\times10^{-4}\,{}^\circ\)/day, one degree in \(2480\) days, \(6.8\) years.

Above \(63.43^\circ\) the sign flips and perigee regresses: the critical inclination separates advance from regress.

Check. Dimensions: \(J_2(R/p)^2\) is dimensionless, so \(\dot\varpi\) has those of \(n\) ✓. \(J_2\to0\): no precession ✓. \(e\to0\): the rate stays finite, \(\tfrac32nJ_2(R/a)^2\), agreeing with the deck's near-circular effective-potential calculation at \(i = 0\) ✓. The factor \(e\) cancelled between \(|\vec A| = GMe\) and \(\int(1+e\cos\nu)^2\cos\nu\): the apse line of a nearly circular orbit turns at a finite rate even though it is barely defined.

Another way

Unit 08's Binet equation with the extra force: in the equatorial plane \(F/m = -GM/r^2 - \tfrac32GMJ_2R^2/r^4\), so \(u'' + u = \frac1p + \tfrac32J_2R^2u^2/p\) to first order. Linearising about the circular orbit \(u_0 = 1/p\) gives \(u'' + (1 - 3J_2R^2/p^2)u = \text{const}\): radial oscillations at angular frequency \(1 - \tfrac32J_2(R/p)^2\) per radian of \(\theta\), so perigee recurs after \(2\pi(1 + \tfrac32J_2(R/p)^2)\) and advances \(3\pi J_2(R/p)^2\) per orbit — \(\tfrac32nJ_2(R/p)^2\) per unit time. It only covers small \(e\); the \(\vec A\) route covers all.

Idea. The formula binds three independently measured numbers. Where it holds, it confirms the surface is a fluid level surface; where it fails, it says the planet is not in hydrostatic equilibrium.

1. Jupiter: \(q = (1.7585\times10^{-4})^2(7.1492\times10^7)^3/1.26686534\times10^{17} = 0.08919\), so \(f = \tfrac32(0.0146965) + \tfrac12(0.08919) = 0.02204 + 0.04460 = 0.06664\). Measured \(0.06487\): \(2.7\%\) high.

2. Mars: \(\omega_{\mathrm s} = 7.0882\times10^{-5}\ \mathrm{s^{-1}}\), \(q = 0.004595\), so \(f = 0.002941 + 0.002298 = 0.005238 = 1/190.9\). Measured \(1/169.8\): \(11\%\) low.

3. Uniform Jupiter: \((2(0.06487) - 0.06487^2)/5 = 0.02511\). Measured \(0.01470\), a ratio of \(0.585\) against the Earth's \(0.809\). Jupiter's mass is far more concentrated towards its centre than the Earth's — a compressible hydrogen envelope over a dense interior.

4. Mars. The formula assumes the surface is a level surface of the rotating-frame potential, which is true for a fluid. Mars has a rigid lithosphere carrying enormous non-hydrostatic loads — the Tharsis volcanic province above all — which contribute to \(J_2\) and to the shape without being in fluid equilibrium. Jupiter, a fluid planet, misses only by the second-order terms, which are large because \(q\) is not small.

Check. \(q\) is dimensionless: \(\mathrm{s^{-2}\,m^3/(m^3\,s^{-2})}\) ✓. The size of Jupiter's miss is consistent with second order: \(q^2 \approx 0.008\), and the miss is \(0.0018\), a fraction of that ✓. For the Earth the same formula missed by \(0.055\%\), with \(q^2 \approx 1.2\times10^{-5}\) ✓.

Idea. Differentiating a harmonic family with respect to a parameter gives harmonic functions, because the \(\vec r\)-Laplacian and the parameter derivative commute. The generating function supplies the family.

1. Fix \(\vec r \ne 0\). On the set \(\{(\vec r,s) : |s| \lt |\vec r|\}\) the function \(F = 1/|\vec r - s\hat{\mathbf z}|\) is smooth (its denominator is \(\ge r - |s| \gt 0\)) and \(\nabla_{\vec r}^2F = 0\), since \(\vec r \ne s\hat{\mathbf z}\). By the law of cosines and the generating function with \(x = \cos\theta\), \(t = s/r\), \(F = \sum_ns^n\,r^{-n-1}P_n(\cos\theta)\), so \(r^{-n-1}P_n(\cos\theta) = \frac1{n!}\partial_s^nF\big|_{s=0}\). Mixed partial derivatives of a smooth function commute (Clairaut, unit 03), so

\[ \nabla^2\big(r^{-n-1}P_n(\cos\theta)\big) = \frac1{n!}\,\partial_s^n\big(\nabla_{\vec r}^2F\big)\Big|_{s=0} = 0 .\qquad\square \]

2. Polynomial: \(P_n\) has the parity of \(n\), so \(r^nP_n(z/r) = \sum_ja_j\,z^{n-2j}r^{2j}\), and \(r^{2j} = (x^2+y^2+z^2)^j\) is a polynomial. Harmonic: let \(G(\vec r,w) = w/|w\vec r - \hat{\mathbf z}|\) on \(\{|w|\,|\vec r| \lt 1\}\), where the denominator is \(\ge 1 - |w||\vec r| \gt 0\), so \(G\) is smooth. By the chain rule \(\nabla_{\vec r}^2\big[1/|w\vec r - \hat{\mathbf z}|\big] = w^2\,(\nabla^2h)(w\vec r)\) with \(h(\vec y) = 1/|\vec y - \hat{\mathbf z}|\), and \(w\vec r \ne \hat{\mathbf z}\), so \(G\) is harmonic in \(\vec r\). Expanding with \(t = wr\): \(G = w\sum_n(wr)^nP_n(\cos\theta) = \sum_nw^{n+1}r^nP_n(\cos\theta)\). Taking \(\frac1{(n+1)!}\partial_w^{n+1}\) at \(w = 0\) and commuting as in part 1 gives \(\nabla^2(r^nP_n) = 0\). \(\square\)

3. \(r^2P_2 = \tfrac12(3z^2 - r^2) = z^2 - \tfrac12(x^2+y^2)\), with Laplacian \(2 - 1 - 1 = 0\) ✓. For \(\phi = 3z^2/r^5 - 1/r^3\) the deck gives \(\partial_x\phi = 3x/r^5 - 15z^2x/r^7\) and \(\partial_z\phi = 9z/r^5 - 15z^3/r^7\). Differentiating again with \(\partial_jr = x_j/r\):

\[ \partial_x^2\phi = \frac3{r^5} - \frac{15x^2}{r^7} - \frac{15z^2}{r^7} + \frac{105x^2z^2}{r^9},\qquad \partial_z^2\phi = \frac9{r^5} - \frac{90z^2}{r^7} + \frac{105z^4}{r^9}, \]

and \(\partial_y^2\phi\) as \(\partial_x^2\phi\) with \(y\) for \(x\). Summing: \(\frac{15}{r^5} - \frac{15(x^2+y^2)}{r^7} - \frac{120z^2}{r^7} + \frac{105z^2(x^2+y^2+z^2)}{r^9} = \frac{15}{r^5} - \frac{15(x^2+y^2+z^2)}{r^7} = 0\) ✓.

4. The tidal potential inside a body is, at order \(\ell = 2\), \(-\frac{GM}{d^3}r'^2P_2(\cos\psi)\), a multiple of \(r'^2P_2\), which is harmonic by part 2 (with the axis along \(\hat{\mathbf s}\), the direction of the other body, and \(\psi\) the angle from it). The tidal tensor is its Hessian, and the trace of the Hessian is the Laplacian (unit 05): zero. \(\square\)

Check. \(n = 0\): part 1 says \(1/r\) is harmonic, which was the input ✓. \(n = 1\): \(z/r^3\) and \(z\) are harmonic — the field of a dipole, and a linear function ✓. Parts 1 and 2 are related by \(\vec r \mapsto \vec r/r^2\) with a factor \(1/r\): \(r^{-n-1}P_n = r^{-1}\cdot(r^nP_n)(\vec r/r^2)\), a map called the Kelvin transform. That it carries every harmonic function to a harmonic function is a general theorem, not proved here and not needed: parts 1 and 2 were proved separately.

Idea. Tides go as mass over distance cubed, not squared, because the uniform \(\ell = 1\) part of the field falls with the Earth and only its variation across the Earth remains.

1. \(GM_{\mathrm{Moon}}/d^3 = 8.632\times10^{-14}\ \mathrm{s^{-2}}\) and \(GM_\odot/d_\odot^3 = 3.964\times10^{-14}\ \mathrm{s^{-2}}\): the Sun's tide is \(0.459\) of the Moon's, although the Sun's pull on the Earth is \(179\) times the Moon's.

2. Take a point on the equator with both bodies in the equatorial plane, and \(\alpha\) the angle from the point to a body. Then \(P_2(\cos\alpha) = \tfrac12(3\cos^2\alpha - 1) = \tfrac14 + \tfrac34\cos2\alpha\): as the Earth turns, each body contributes a wave in \(2\alpha\), twice a day, of amplitude proportional to its tidal strength. At new or full Moon the two bodies are at \(\alpha\) differing by \(0\) or \(\pi\), so the waves are in phase and add: \(1 + 0.459\). At quarter Moon they differ by \(\pi/2\), so the waves differ in phase by \(\pi\) and subtract: \(1 - 0.459\). The ratio of ranges is \(1.459/0.541 = 2.70\).

3. Consecutive terms carry an extra \(r'/d\) (times a ratio of \(P_\ell\) values of order one): \(R/d_{\mathrm{Moon}} = 0.0166\) and \(R/d_\odot = 4.3\times10^{-5}\). The lunar \(\ell = 3\) tide is a percent-level correction; the solar one is negligible.

4. On the undisturbed surface \(\Phi_\oplus = -GM_\oplus/R\). Raising the surface by \(h\) changes \(\Phi_\oplus\) by \(gh\), \(g = GM_\oplus/R^2\). A level surface of \(\Phi_\oplus + \Phi^{(2)}\) therefore has \(gh + \Phi^{(2)}(R) = \text{const}\), and choosing the constant so that \(h\) has zero average over the sphere,

\[ h = -\frac{\Phi^{(2)}}{g} = \frac{GM_{\mathrm b}R^2/d^3}{GM_\oplus/R^2}P_2(\cos\psi) = \frac{M_{\mathrm b}}{M_\oplus}\Big(\frac Rd\Big)^3R\,P_2(\cos\psi). \]

For the Moon the prefactor is \(\frac{4.9029\times10^{12}}{3.9860\times10^{14}}(0.016592)^3(6.378\times10^6\ \mathrm m) = 0.358\) m; with \(P_2\) ranging from \(1\) to \(-\tfrac12\), the range is \(1.5\times0.358 = 0.54\) m.

Check. Dimensions: \((M_{\mathrm b}/M_\oplus)(R/d)^3\) is dimensionless, so \(h\) is a length ✓. Limit \(d\to\infty\): \(h\to0\) ✓. The number: real open-ocean tides are of order half a metre to a metre, with larger values at coasts where resonance amplifies them ✓. The bulge is at \(\psi = 0\) and \(\psi = \pi\), because \(P_2\) is even — the deck's two bulges again.