Which mass distributions guarantee two bodies have non-Keplerian orbits? Which non-spherical distributions still allow noncircular Keplerian orbits?

Question

The trajectories of two point masses or spherically symmetric masses with respect to their center of mass are conic sections or Kepler orbits.

Consider that the bodies have finite size with respect to their separation and not necessarily uniform, or even spherically symmetric mass distributions.

In that case what are the constraints on their mass distributions and orientations such that their orbits are still Keplerian? Or does any deviation from spherical symmetry of one or both body immediately result in a non-Keplerian orbit?

This answer to the Astronomy SE question Which point in an orbiting body most closely follows its Keplerian trajectory? explains:

If one or both of the bodies have non-spherically symmetric density distribution, the orbits will no longer be keplerian.

but without citing sources nor using math. I can certainly imagine "... are no longer necessarily Keplerian" but is this always true in absolute?

Just for example couldn't two tidally locked ellipsoids have circular orbits about their center of mass?

If so, perhaps that's the only exception, but perhaps not. If they orbited in non-circular orbits would their apparent libration also mean that their centers of mass no longer followed strict Keplerian orbits? Would they differ from ellipses, or still be elliptical but no longer exactly follow equal area per unit time trajectories? Would there be another point within the bodies besides their centers of mass that still did follow a Keplerian orbit?

Question: Which mass distributions guarantee that two bodies will always have non-Keplerian orbits? Which non-spherical distributions still allow for noncircular Keplerian orbits?

Answer 1

This answer is trivial, but Keplerian orbits are assured as long as the gravitational fields are spherically symmetric and thus decay as $1/r^2$. There are an infinite number of mass distributions that can generate a given field. For example a cubic shell of side $l$ with a mass density $(l/2)/((l/2)^2+x^2+y^2)^{(3/2)}$ on each square surface will produce a spherically symmetric field. In this example, the mass density at the corners is about 1/5 of that at the centers of the faces.
Some insight can be gained by analogy with electric field/flux. 'Tubes' of gravitational flux must terminate on elements of mass. If the incoming flux is spherically symmetric, then the mass density of the shell must go as $1/r^2$ where $r$ is its distance from some interior point. In addition, since the surface will generally be at an angle to the incoming radial tube, the mass density gets reduced by the cosine of that angle. This is exactly how the result for the cubic shell is derived.
So, in theory, any shell that completely encloses a central point can produce a spherical field. This result can obviously be extended to solid bodies.
It is interesting to speculate whether a field can go as $1/r^2$ only in one (say equatorial) plane but not in other planes. I think the answer is 'no' because the far-field from any object is spherically symmetric and this boundary condition seems to preclude such solutions except by approximation over a limited range of radii (e.g. perhaps by sticking a dumbell through the hole in a donut).

Answer 2

There are infinitely many such shapes that work, even ones without anything like spherical symmetry at all. To give a simple example, recall the method of images for spheres in electromagnetism. When you move a point charge near a conducting sphere, it induces a non-symmetric charge distribution on the sphere. But the electrostatic field of that charge distribution is exactly mimicked by a fictitious point charge inside the sphere, called the image charge.

That means, conversely, that the gravitational field of a point mass can be exactly mimicked by a non-symmetric mass distribution encasing it. By superposing multiple such distributions, you can get a mass distribution that looks as weird as you want, but whose field has exactly $1/r^2$ behavior outside it.