Gravitational $n$-body problem with tidal forces This year I'm working on modelling the gravitational $n$-body problem using Newton's law of gravity where I assume that for large enough distances, planetary bodies can be modelled as point masses. But, it's not clear to me at what point this approximation is a good one. 
If we assume that our bodies can be modelled as point masses then we can assume that these bodies have zero rotation. However, in the case of $n=2$ we already have a problem for short distances:


*

*Over millenia, the geometry of the planets will actually deform in ways similar to a fluid would, if we assume zero rotation. 

*Because of $(1)$, we would have to drop the spherical approximation over large time scales and hence the point mass approximation. 

*Suddenly rotation matters. 


For this reason, it seems prudent to use Gauss's Law of Gravitation which takes into account the planetary surface. However, I wonder whether I can rigorously demonstrate that for sufficiently large distances, Newton's approximation works very well. 
This leads me to the question of what distances would be sufficiently large for the approximation to work well (i.e. accurately model orbits) for millions of years?
 A: The standard way to represent the gravitation of a not quite spherical object is to use spherical harmonics. I'll just look at the first non-spherical term of the spherical harmonic expansion, commonly called the $J_2$ term. This term results from an object's oblateness. The gravitational potential of an oblate spheroid as a function of distance $R$ from the object's center of mass and an angle $\theta$ between the position vector and the object's equatorial plane is:
$$\phi(R,\theta) = -\frac{\mu} R \left(1-J_2\left(\frac{r}{R}\right)^2 \frac{3\sin^2\theta - 1}2\right)$$
where $\mu\equiv GM$ is the object's gravitational parameter, $r$ is the object's  equatorial radius, and $J_2$ is the object's dimensionless second dynamic form factor. The two gas giant planets have the largest $J_2$ values, 0.016297 and 0.014733 for Saturn and Jupiter, respectively.
For the Earth, $J_2$ is just 0.001083, over an order of magnitude smaller than that for the gas giants. This is more than large enough to perturb the orbits of satellites in low Earth orbit; the sun synchronous satellites that various space agencies use to observe the Earth depend on this perturbation. The Earth's oblateness has observable effects as far out as the orbit of the Moon.
So does one need to model oblateness (i.e., $J_2$) for the gravitational interactions between planets? The answer is maybe, and only when the planets in question are closest to one another (inferior conjunction or opposition). The reason is the $\mu_s/{R_s}^2$ gravitational acceleration of a planet toward the Sun generally overwhelms the $\mu_p/{R_p}^2$ acceleration toward even the closest other planet by many orders of magnitude. The $J_2$ term is orders of magnitude smaller yet. In particular, to be able to ignore the oblateness of a planet, we need
$$\begin{aligned}
\frac{\mu_p}{{R_p}^2} \left(\frac{r_p}{{R_p}}\right)^2 3 J_2 & \frac{|3\sin^2\theta -1|} 2 \lll \frac {\mu_s} {{R_s}^2} \\
&\text{or}
\end{aligned}$$
$$\frac{\mu_p}{\mu_s} \left(\frac{R_s r_p}{{R_p}^2}\right)^2 3 J_2 \frac{|3\sin^2\theta -1|} 2 \lll 1 \tag{1}$$
Jupiter and Saturn are closest to one another when Jupiter is at aphelion, Saturn is at perihelion, and the two are aligned just right. In this case the unitless "oblateness factor" described by equation (1) of Jupiter's oblateness on Saturn's orbit and of Saturn's oblateness on Jupiter's orbit are about $3.2\times10^{-12}$ and $2.5\times10^{-13}$, respectively, assuming the planets are aligned just right.
These are small values, but not ideally small values ("ideally small" would be smaller than 10-16). The uncertainties in the gravitational parameters and in the planet's position vectors overwhelm these very small factors. Moreover, the moments when Saturn and Jupiter can "feel" the oblateness of one another occurs but briefly, when the Sun, Jupiter, and Saturn are nearly in line with one another. Mars can similarly briefly "feel" the oblateness of Jupiter, and Uranus, that of Saturn.
These effects are very small, much smaller than the uncertainties in the gravitational parameters and in the propagated positions of the planets. They can safely be ignored. On the other hand, if you are modeling the orbits of satellites about a planet, you very much need to account for the non-spherical nature of the planet. Oftentimes, even the $J_2$ term isn't enough.
