It's well known that in low orbits around astronomical bodies with uneven mass distribution, orbits will shift around, and the orbiting object may crash into the surface.
Why is this the case? Newtonian gravity seems to dictate that we can treat all masses as point masses. In this case, a satellite orbiting a planet would simply be two point masses orbiting around their common center of mass. Why would the mass distribution matter at all? The planet, I thought, could even be some really really long cylinder and a satellite would still orbit in perfect ellipses around its center of mass.
Best Answer
It is incorrect to think that
Indeed, you can't do that, and one of the consequences of that fact is the orbital decay that you refer to.
The theorem you're thinking of holds only for bodies with perfect spherical symmetry, and states that for such a body, the gravitational field outside it is identical to that of a point mass at the body's centre.
However, you run into trouble immediately after you leave that setting. For example, the interactions between two spherical bodies cannot, in general, be reduced to the interactions between point masses. You can indeed substitute one body for a point mass to find its field, but then you need to find the force on the second body caused by this spatially inhomogeneous field.
It is, sometimes but not always, a good approximation to treat this as a total force acting on its centre of mass:
If the second body is very small compared to the bodies' separation, then the relative variation of the first body's gravitational field over the second body's spatial extent will probably be small, and it can be neglected. As a better approximation, you can include the first-order variation of the first body's gravitational field, which is called a tidal force.
If the second body is perfectly rigid, it will respond to any inhomogeneous force field by means of displacements of the COM caused by the total force, which is again the point-particle model. As a better approximation, you can include linear deformations caused by first-order variations in the force field. If you couple this with the tidal forces, you'll have a nice basic description of tidal drag.
Notice, for example, that the Moon obeys both conditions pretty well.
For a body with an arbitrary mass distribution, you can probably now guess that it gets a lot hairier. If the distribution is really wild, then there is nothing for it but to calculate the field numerically, and orbits will be hard to come by. If the body is mostly spherical, then there is some hope.
The Earth and the Moon are sort of on the latter case, but their distributions are far from uniform. The Earth's geoid, which describes its shape as far as gravity goes, is not spherical:
Image source. Note the exaggerated scale: 50 mGal corresponds to a change in $g$ of 0.005%.
For the Moon, mass concentrations underneath mare make the gravity even more uneven.
The tool to deal with these departures from sphericality is borrowed from electrostatics: the multipole expansion. That is, if the body is mostly spherical, then you can treat its gravitational field as a superposition of a point mass plus a point dipole, quadrupole, octupole, and so on.
$$\Phi(\mathbf r)=\sum_{l,m}Q_{lm} \frac{Y_l^m(\theta, \phi)}{r^{l+1}}.$$
The good news is that the contribution of each succeeding term gets weaker and weaker, and most importantly they get weaker the further away the test mass is. As a bonus, you can usually make the dipolar terms to zero by putting the origin at the centre of mass.
For the Earth, the first non-point-mass term that you need to account for is the zonal quadrupole term, with $l=2$ and $m=0$, which describes the Earth's bulge at the equator and the corresponding flatness at the poles. Including that as in this article, you get a gravitational potential that looks like $$ \Phi=-\frac{GM}{r}\left[1-J_2 \left(\frac{a_1}{r}\right)^2P_2(\cos\theta)\right]. $$ Here $a_1$ is the Earth's equatorial radius, and $J_2\approx0.001$ is called the quadrupole moment coefficient. For low Earth orbit, this effect is on the order of 0.1%, which is small but it's a perturbation that needs to be dealt with. For orbits which are further than 3 or 4 $a_1$ from Earth, the effect becomes correspondingly smaller, in a quadratic fashion. This is on the order of 0.01% of the Keplerian force, but again it's a perturbation that needs to be accounted for and which can add up over many orbits to give sizeable effects.