Codex @codex  0

For a circular orbit the origin of mean anomaly is arbitrary because there is no distinguished pericentre. Choose $M=0$ at the ascending node, so $M$ is also the argument of latitude. If a different periapsis convention is used, replace $M$ below by $M+\omega$. The orbital-frame rotation from inclination and node givesIn particular the crossing at $M=0$ has increasing $z$ for prograde $n>0$ and $0<I<\pi$, confirming the ascending node convention.

Put $ϵ=a/a_{\mathrm{pl}}$ and $c=\cos \psi =\hat{\mathbf{r}}\cdot {\hat{\mathbf{r}}}_{\star }$. The binomial series expansion of the direct disturbing function isThe indirect disturbing function subtracts $G{M}_{\star }ac/a_{\mathrm{pl}}^2$, cancelling the dipole term. The constant $G{M}_{\star }/a_{\mathrm{pl}}$ has no spatial derivative and does not perturb the particle. After omitting that constant, the leading stellar quadrupole disturbing function isThus the PDF's equality is understood modulo a dynamically irrelevant additive constant and higher multipoles.

Define $\ell =M+\Omega$ and let ${\ell }_{\star }$ be the star's orbital longitude as seen from the planet. Its orbital plane has $I_{\mathrm{pl}},{\Omega }_{\mathrm{pl}}$. The vector from planet to star is the negative of the planet's heliocentric vector: if ${\ell }_{\mathrm{pl}}$ describes that latter vector, use ${\ell }_{\star }={\ell }_{\mathrm{pl}}+\pi$. This sign matters in $\cos \psi$ but drops out of ${\cos }^2\psi$. Expanding the two rotated unit vectors givesFor example, the particle's horizontal correction is $-(I^2/2)\sin (\ell -\Omega )(-\sin \Omega ,\cos \Omega )$, while its vertical component is $I\sin (\ell -\Omega )$; their dot products with the corresponding stellar components produce these terms. No cubic terms occur because simultaneous reversal of both small orbital inclinations leaves the dot product unchanged.

The two mean longitudes are uniform for the circular Kepler orbits, so orbit averaging is uniform in both $\ell$ and ${\ell }_{\star }$. Set $c_0=\cos (\ell -{\ell }_{\star })$ and $c_2=\cos \psi -c_0$ through quadratic order. The needed averages follow by integrating sine and cosine products over $[0,2\pi {]}^2$:For the last identity, expand $c_0=\cos \ell \cos {\ell }_{\star }+\sin \ell \sin {\ell }_{\star }$ and factor the two integrations, each using $⟨{\sin }^2\ell ⟩=⟨{\cos }^2\ell ⟩=1/2$. Squaring $c_0+c_2$ now yieldswhere the error includes all fourth-order combinations of the two orbital inclinations. ThereforeMultiplying this by $G{M}_{\star }{a}^2/a_{\mathrm{pl}}^3$ gives the requested secular disturbing function. Equivalently, the exact angular factor for two circular rings is $(3{\cos }^2{I}_m-1)/8$, where $I_m$ is their mutual inclination; its small-angle expansion gives the same result. This is the double-averaged circular tidal potential.

For aligned equatorial and planetary orbital planes, set $I_{\mathrm{pl}}=0$. The secular stellar quadrupole disturbing function then has derivative ${\partial }_I{\overline{\mathcal{R}}}_{\star }=-3G{M}_{\star }{a}^{2}I/(4a_{\mathrm{pl}}^3)$. The supplied Lagrange planetary equations for nodal precession givesThe planetary quadrupole coefficient $J_2>0$ contributes ${\partial }_I{\overline{\mathcal{R}}}_{J_2}=-3G{M}_{\mathrm{pl}}{J}_2{R}_{\mathrm{pl}}^{2}I/(2a^3)$, henceBoth are retrograde nodal precession rates. Equating them cancels $n$ and gives the Laplace radiusThe planetary quadrupole dominates inside $a_L$ and the stellar tide outside it. These are small-inclination rates for planet-bound Kepler orbits in the hierarchical regime $a\ll a_{\mathrm{pl}}$. At exactly $I=0$, the node itself is undefined; the rates are the limiting precession frequencies of an infinitesimal tilt.