When I was in the university (in the late 90s, circa 1995) I was told there had been research investigating the $2$ (the square of distance) in the Newton's law of universal gravitation.
$$F=G\frac{m_1m_2}{r^2}.$$
Maybe a model like
$$F=G\frac{m_1m_2}{r^a}$$
with $a$ slightly different from $2$, let say $1.999$ or $2.001$, fits some experimental data better?
Is that really true? Or did I misunderstand something?
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$
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
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.
There has indeed been such research about the Pioneer anomaly: Two spacecrafts launched in the 1970s into the outer solar system did not move quite as expected (as calculated due to gravity and solar wind) after ca. 1980. Only in/after the 2000-2010 decade the source of the discrepancy, accidental thrust by thermal radiation, became generally accepted consensus. Previously, it was at least conceivable to interprete the data as containing hints of subtle differences between actual observed gravity and our theoretical understanding of gravity.
There was indeed some talk of the exponent on $r$ during the late 90's and the early years of the 21st century. The problem, as I recall, was dark matter, which can only be observed indirectly by looking at the anomalous rotation of galaxies. It was suggested that perhaps Newton's law broke down under certain conditions. Again, as I recall, that while a number of papers were published, nothing much came of the idea.
Building on Jim Graber's answer:
We can absorb the perturbation term into the correction to the power law. $$ F = -GM/r^2 \frac1{(r/r_0)^{\delta(r)}} \approx -GM/r^2 \left( 1 - \delta(r)(r/r_0 - 1) \right)$$
and we get
$$ \delta(r) = -\frac{3h}{r^2c^2}\frac1{r/r_0 - 1}$$
I'm not sure about the physical meaning of $r_0$ though (renormalization scale?).
If $r/r_0 > 1$ then $\delta(r)$ is negative and we have $\alpha$ like 1.9999...
JPL who calculates the orbits of celestial bodies to a high precision uses an expression for the acceleration of one body of negligible mass due to the gravitational force of one other body that looks like:
$\frac{d\bar{v}}{dt}=-\frac{GM}{r^2}{(1-\frac{4GM}{rc^2}}+\frac{v^2}{c^2})\hat{r}+ \frac{4GM}{r^2}(\bar{v}\cdot{\hat{r}})\frac{v^2}{c^2}\hat{r}$
Ignoring velocity dependent parts we have: $\frac{d\bar{v}}{dt}=-\frac{GM}{r^2}\hat{r}+\frac{4(GM)^2}{c^2r^3}\hat{r}$
So actually keeping $a = 2$ but adding a small "inverse cubic" part is actually done do fit experimental values better, even though the inverse cube term is not invented but comes from trying to approximate general relativistic effects using what is known as "the post-Newtonian expansion".
The reason for this is a bit complicated, but basically adding a small inverse cube part as well as velocity-dependent parts explains what it is known as the "anomalous precession of perihelion".
See for instance expression 4-61 in this paper, titled "Formulation for Observed and Computed Values of Deep Space Network Data Types for Navigation"
