It is incorrect to think that
Newtonian gravity seems to dictate that we can treat all masses as point masses.
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.
Best Answer
Addressing only the question of "Will the Earth grow cold?", the fate of the universe turns on cosmological considerations.
Short answer: Yes.
Slightly longer answer: By measuring the thermal gradient of the crust in deep caves and boreholes at many places on the continental and oceanic crust it is possible to find an approximation to total geothermal power.
Aside: A substantial portion of the geothermal flux in this era is from radioactive decay rather than from gravitational potential.
Results published (the link is the the arXiv preprint, but the paper also appeared in Phys. Rev. B) by the Borexino collaboration in 2010 and by the KamLAND collaboration in 2011 (Nature Geoscience) are consistent with roughly half of the geothermal power of the Earth being due to radiological decay.
These measurements also put strict upper limits on the power of a theorized natural nuclear reactor at the core (the data are now consistent with zero reactor power).
Both Borexino and KamLAND are large anti-neutrino detectors and are directly sensitive to the anti-neutrino emissions of radioactive beta decays such as those found in the Uranium and Thorium chains (but not to Potasium-40 on account of using inverse beta decay as the detection mechanism). From this data we can reconstruct the overall radioactive decay activity in the deep Earth, and compute the total power represented.
Disclaimer: I worked on KamLAND for 3 years, but am not named as an author on the paper cited herein.
Further disclaimer: much of the text here is adapted from my earlier answer on Skeptics.SE.