Can "prolateness" be characterized by a $J_2$ coefficient like "oblateness"?

Question

When an imperfect sphere can be obtained by rotating an ellipse around an axis, it is called a spheroid. There are two types of spheroids, oblate ones and prolate ones.

Most solar bodies can be assimilated to an oblate spheroid, for the purpose of first approximation computations, such as the prediction of satellite movements orbiting that body. The Earth is one of them. On the other hand, some moons in the solar system have a prolate approximation.

When a planet is approximated by an oblate spheroid, its oblateness is characterized by a single coefficient often referred to as $J_2$. When $J_2$=0, the planet's shape is close to a sphere. High positive $J_2$ means that the eccentricity of the ellipse of revolution is high. For example the J2 of Mars is almost twice that of Earth, according to this NASA site.

Can $J_2$, and it alone, be used to characterize the shape of an prolate body? If yes, is there any concrete example of such a characterization?

This question was triggered by Any exact analytical solution for non-Keplerian orbits ..., where it can be read

These refer to Kelperian-like orbital parameterization of orbits around an oblate/prolate spheroid characterized by $J_2$.

At first reading, my interpretation is that the sentence implies that oblate and prolate spheroids differ only by the value of a single characteristic coefficient. In addition, in further comment threads, @uhoh seems to claim that the sign of $J_2$ is what makes the difference.

Answer 1

This is an expansion of my comment.

Very roughly speaking, $J_2$ is a measure of the following: if one point is over the body's pole, and another point is over the equator at the same distance $R$ from the body's center (not surface!), how much is the first point "uphill" compared to the second point w.r.t. the body's gravity?

Let us first consider the case of two equal point masses. At any point which is not on the same line with them and not on their plane of symmetry, the total gravity due to those masses does not point at their midpoint, but deviates somewhat towards the near mass (see picture 1 in the attached image: $A$ and $B$ are the masses; $O$ is their midpoint; $X$ is the point where we measure the gravity; sorry for the poor quality). For the sum of the gravity accelerations $\mathbf{g}_A$ and $\mathbf{g}_B$ to point to $O$, the ratio of their magnitudes $\frac{g_A}{g_B}$ would have to be the same as $\frac{|XA|}{|XB|}$ (pic. 2), i.e., $g_A$ would have to be smaller than $g_B$. But $g_A$ is in fact greater than $g_B$, so the sum $\mathbf{g}_A+\mathbf{g}_B$ points toward some point on the segment $OA$.

It means that if you move along a quarter-circle with the center at $O$ from a point above the "equator" to a point above the "pole" (pic. 3), the angle between the direction of your movement and the direction of local gravity is always (except at the endpoints) smaller than $90^\circ$, i.e., you're moving "downhill". It means that $J_2$ is negative in this case.

Now, if we consider a flat disc, we can divide it into pairs of small sections located symmetrially w.r.t. to the center (pic. 4). For any two such sections, the sum of the gravity accelerations due to them does not point to the center of the disc, but deviates toward the near section. As a result, the total gravity of the disc points to a point inside the near half-disc. So if you again move along a quarter-circle with the center at the center of the disc from a point above the "equator" to a point above the "pole" (pic. 5), the angle between the direction of your movement and the direction of local gravity is greater than $90^\circ$, i.e., you're moving "uphill". So $J_2$ is positive in this case.

Answer 2

For simplicity of calculation let's simulate a prolate body as a central potential plus two smaller masses north and south of center

And let's simulate an oblate body as a central potential plus an equatorial ring of point masses.

If we set $G=1$ for convenience and calculate the gravitation potential $GM/r$ we can see that it peaks at the equator for the oblate case and has a minimum at the equator for the prolate case.

These wiggles are not pure sinusoids, but we can see that their behaviors are mostly like $\sin^2(\theta)$ and $-\sin^2(\theta)$ and since $J_2$ is the coefficient in front of such a term (from Geopotential model between Equations 9 and 10):

$$u_{J2} = J_2 \frac{1}{r^3}\frac{1}{2} \left( 3 \sin^2(\theta)-1 \right)$$

we can see that for an oblate body $J_2$ will be positive and for a prolate body it will be negative.

---

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D 


def phi_central(x, m):
    """central potential"""
    x0 = np.zeros(3)
    r = np.sqrt(((x - x0)**2).sum(axis=-1))
    return m / r, x0

def phi_polar(x, m, h):
    """quadrupole potential (axial)"""
    zhat = np.array([0, 0, 1])
    x1, x2 = h * zhat, -h * zhat
    r1 = np.sqrt(((x - x1)**2).sum(axis=-1)) # top
    r2 = np.sqrt(((x - x2)**2).sum(axis=-1)) # bottom
    return m * (1/r1 + 1/r2), x1, x2
                 
def phi_equatorial(x, m, r, N=100):
    """quadrupole potential (equatorial (xy))"""
    theta = np.linspace(0, 2*np.pi, N+1)[:-1]
    # ring in the xy plane
    ring = np.vstack([r * f(theta) for f in (np.cos, np.sin, np.zeros_like)])
    r = np.sqrt(((x[..., None] - ring)**2).sum(axis=-2)) # equator
    return m / r.mean(axis=-1), r, ring


mcen, mpol, meq = 1, 0.1, 0.3
hpol, req = 0.3, 0.9

theta = np.linspace(0, 2*np.pi, 361)
R = 1.0
N = 1000

xz_plane = np.stack([R * f(theta) for f in (np.sin, np.zeros_like, np.cos)], axis=-1)

phi_c, xcen = phi_central(xz_plane, mcen)

phi_p, xp1, xp2 = phi_polar(xz_plane, mpol, hpol)

phi_e, r, ring = phi_equatorial(xz_plane, meq, req, N=N)
                    
fig, ax = plt.subplots(1, 1)
ax.plot(np.degrees(theta), phi_c)
ax.plot(np.degrees(theta), phi_p)
ax.plot(np.degrees(theta), phi_e)
ax.set_ylim(0, 1.1)
ax.set_xlabel('theta (deg)')
ax.set_ylabel('gravitational potential')
plt.show()

fig = plt.figure()
ax  = fig.add_subplot(1, 1, 1, projection='3d', proj_type = 'ortho')
ax.plot([0], [0], [0], 'ok')
x, y, z = zip(xp1, xp2)
ax.plot(x, y, z, 'or')
x, y, z = ring
ax.plot(x, y, z, '-b')
x, y, z = xz_plane.T
ax.plot(x, y, z, '-g')
plt.show()