The distance square in the Newton's law of universal gravitation is really a square?

This was suggested by Asaph Hall in 1894, in an attempt to explain the anomalies in the orbit of Mercury. I retrieved the original article in http://adsabs.harvard.edu/full/1894AJ.....14...49H

Interestingly, he mentions in the introduction that Newton himself had already considered in the Principia what happens if the exponent is not exactly 2, and had concluded that the observations available to him strongly supported the exact power 2!

The story is retold, e.g., on p.356 of N.R. Hanson, Isis 53 (1962), 359-378.

See also Section 2 of http://adsabs.harvard.edu/full/2005MNRAS.358.1273V

Let's first see why the inverse square form is special. Betrand's Theorem states that only two types of central potentials will produce stable orbits. The harmonic oscillator potential $V=\frac{1}{2}kr^2$ and the potential $V=-\frac{k}{r}$ that will produce an inverse square force law. Obviously the age of the universe is finite, so the fact that planet's orbits survived until now need not imply it will continue to be so in the future.

Another argument why this type of potential is so common is that, when doing quantum field theory, the propagator (details depend on whether the particle is a (gauge) boson, fermion or scalar, i will stick with scalars for now) has form

$$\frac{1}{q^2+m^2}$$

Thus if this particle where the force carrier of your force with coupling $g$ the potential is basically the fourier transform of the propagator

$$V(r) =-g^2\frac{1}{(2\pi)^3}\int d^3k\frac{4\pi}{q^2+m^2}e^{i\vec{k}\cdot \vec{r}} = -g^2\frac{1}{r}e^{-mr}$$

This is the famous Yukawa potential. For massless force carriers the damping term goes to 1 and the force becomes long range with a inverse square force law. Upto small details this is analogous to the gauge boson case, e.g. the masslessness of the photon makes the EM force long range, where as the massiveness of W,Z bosons make weak forces short-ranged.

Above derivations use the three space dimensions. Theories with extra dimensions have suggested that large extra dimensions will alter the inverse square law at some not-so-short distances (sub-mm range). Published experimental results are to be found e.g. from the Eöt-Wash group ( http://www.npl.washington.edu/eotwash/experiments/shortRange/sr.html ) and are available on the arXiv.

One potential tested here is here $$V(r)=-G\frac{m_1m_2}{r}(1+\alpha\exp(-r/\lambda))$$

The below plot shows the exclusion limits for both parameters $\alpha$ and $\lambda$

But of course, Newton’s theory is not correct; instead Einstein’s theory is correct. If you use general relativity GR, you usually talk about curvature, etc. rather than forces.
Nevertheless, the results can be expressed in terms of an effective force.

This reference http://farside.ph.utexas.edu/teaching/336k/Newton/node116.html gives

$F = -GM/r^2 -3GMh^2/(c^2r^4)$

where h is momentum as the first order correction. Higher orders have been calculated by the PPN and gravitational wave people. This correction is very small except for very fast moving objects. In practise, it applies to bodies orbiting very near to a black hole or neutron star. Famously, it is also responsible for the precession of the perihelion of Mercury.