Does the gravitational attraction near the surface of dense celestial objects diverge from inverse square? Does the gravitational attraction near the surface of dense celestial objects (neutron star, white dwarf itc) diverge (to infinity) from inverse square?

This question is inspired by the similarity between EM and gravity (inverse square force.). The paper by John Lekner here (doi:10.1098/rspa.2012.0133) shows that there is an electrostatic attraction between charged spheres no matter the polarities of the charges and that it diverges at close separation until electrical short for almost all charge ratios. I am wondering if there is a similar kind of gravitational inverse square divergence for anything other than a black hole. 

Actually make this for a black hole as well, although I know that a black hole is not thought to have a normal surface.
Leckner's paper deals with the effect of induced polarization on the spheres. Electrons are redistributed, making the force different from what one would expect. The gravitational counterpart is tidal distortion: since the gravitational field is non-radial when you have two heavy masses close to each other, matter will move to make the surface an equipotential surface. This means that the gravitational acceleration at the surfaces will not be constant at all locations.

Doing an analytic solution of how two ellipsoids attract each other and deform seems to be tractable (e.g. see this question) but algebraically very tedious and likely involves lots of special functions. See the addendum below for an approximate numeric model.

Black holes produce another complication: since spacetime nearby is curved and expanded the meaning of the distance in the inverse square law becomes problematic. The Paczyński–Wiita potential is an approximation of the potential, and it deviates from the $U=-GM/r$ as $U_{PW}=-GM/(r-R_S)$ (where $R_S$  is the Schwarzschild radius). It makes the force increase faster than the classical potential as we approach $r=R_S$.


Addendum: I did a numeric exploration of the force between two ellipsoidal, self-gravitating masses with centres of mass separated by a given distance. To find the shape I started with spheres and adjusted the semi-major axis (while preserving volume) so the potential along the surface became more equal at the poles. After a few iterations this gives a self-consistent shape. Then I calculated the force (the derivative of the potential) due to this shape on the other mass.



The result is indeed that the force increases faster than $1/r^2$
as the bodies approach each other, since they elongate and eventually merge (a bit before this they will deviate from my ellipsoidal assumption). If one multiplies the force by the squared distance the product should be constant for pure $1/r^2$ forces, but it starts to increase as they approach enough. Note that this is a non-rotating model: with rotation the numbers will change and the ellipsoids will become tri-axial, but I suspect the qualitative behaviour stays the same.