¿Puede la "prolatitud" caracterizarse por un coeficiente J2J2J_2 como la "oblatividad"?

Esta es una ampliación de mi comentario.

Hablando muy groseramente,$J_2$es una medida de lo siguiente: si un punto está sobre el polo del cuerpo y otro punto está sobre el ecuador a la misma distancia$R$desde el centro del cuerpo (¡no la superficie!), ¿cuánto es el primer punto "cuesta arriba" en comparación con el segundo punto con respecto a la gravedad del cuerpo?

Consideremos primero el caso de dos masas puntuales iguales. En cualquier punto que no esté en la misma línea que ellos y no en su plano de simetría, la gravedad total debida a esas masas no apunta a su punto medio, sino que se desvía un poco hacia la masa cercana (ver figura 1 en la imagen adjunta:$A$y$B$son las masas;$O$es su punto medio;$X$es el punto donde medimos la gravedad; Perdón por la mala calidad). Para la suma de las aceleraciones de la gravedad$\mathbf{g}_A$y$\mathbf{g}_B$apuntar a$O$, la razón de sus magnitudes$\frac{g_A}{g_B}$tendria que ser igual que$\frac{|XA|}{|XB|}$(foto 2), es decir,$g_A$tendría que ser más pequeño que$g_B$. Pero$g_A$es de hecho mayor que$g_B$, entonces la suma$\mathbf{g}_A+\mathbf{g}_B$apunta hacia algún punto del segmento$OA$.

Significa que si te mueves a lo largo de un cuarto de círculo con el centro en$O$desde un punto por encima del "ecuador" hasta un punto por encima del "polo" (imagen 3), el ángulo entre la dirección de su movimiento y la dirección de la gravedad local es siempre (excepto en los extremos) menor que$90^\circ$, es decir, te estás moviendo "cuesta abajo". Esto significa que$J_2$es negativo en este caso.

Ahora, si consideramos un disco plano, podemos dividirlo en pares de pequeñas secciones ubicadas simétricamente al centro (foto 4). Para cualquiera de estas dos secciones, la suma de las aceleraciones de gravedad debidas a ellas no apunta al centro del disco, sino que se desvía hacia la sección cercana. Como resultado, la gravedad total del disco apunta a un punto dentro del medio disco cercano. Entonces, si nuevamente se mueve a lo largo de un cuarto de círculo con el centro en el centro del disco desde un punto sobre el "ecuador" hasta un punto sobre el "polo" (imagen 5), el ángulo entre la dirección de su movimiento y la dirección de la gravedad local es mayor que$90^\circ$, es decir, te estás moviendo "cuesta arriba". Entonces$J_2$es positivo en este caso.

Para simplificar el cálculo, simulemos un cuerpo alargado como un potencial central más dos masas más pequeñas al norte y al sur del centro.

Y simulemos un cuerpo achatado como potencial central más un anillo ecuatorial de masas puntuales.

si establecemos$G=1$por conveniencia y calcule el potencial gravitatorio$GM/r$podemos ver que alcanza su punto máximo en el ecuador para el caso achatado y tiene un mínimo en el ecuador para el caso alargado.

Estos movimientos no son sinusoides puros, pero podemos ver que sus comportamientos son en su mayoría como$\sin^2(\theta)$y$-\sin^2(\theta)$y desde$J_2$es el coeficiente delante de tal término (del modelo Geopotencial entre las Ecuaciones 9 y 10):

podemos ver que para un cuerpo achatado$J_2$será positivo y para un cuerpo alargado será negativo.

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()