Может ли «вытянутость» характеризоваться коэффициентом J2J2J_2, как «сплющенность»?

Это расширение моего комментария.

Очень грубо говоря,$J_2$является мерой следующего: если одна точка находится над полюсом тела, а другая точка находится над экватором на том же расстоянии$R$от центра тела (не поверхности!), насколько первая точка находится «в гору» по сравнению со второй точкой относительно силы тяжести тела?

Сначала рассмотрим случай двух одинаковых точечных масс. В любой точке, которая не находится с ними на одной линии и не на их плоскости симметрии, общая гравитация этих масс не указывает на их среднюю точку, а несколько отклоняется в сторону ближней массы (см. рисунок 1 на прилагаемом изображении:$A$и$B$это массы;$O$их середина;$X$точка, в которой мы измеряем гравитацию; извините за плохое качество). Для суммы ускорений свободного падения$\mathbf{g}_A$и$\mathbf{g}_B$указать на$O$, отношение их величин$\frac{g_A}{g_B}$должно быть так же, как$\frac{|XA|}{|XB|}$(рис. 2), т. е.$g_A$должен быть меньше, чем$g_B$. Но$g_A$на самом деле больше, чем$g_B$, поэтому сумма$\mathbf{g}_A+\mathbf{g}_B$указывает на некоторую точку на отрезке$OA$.

Это означает, что если вы двигаетесь по четверти окружности с центром в$O$от точки над «экватором» до точки над «полюсом» (рис. 3) угол между направлением вашего движения и направлением местной силы тяжести всегда (кроме конечных точек) меньше$90^\circ$, т. е. вы движетесь «под гору». Это означает, что$J_2$в данном случае является отрицательным.

Теперь, если рассматривать плоский диск, мы можем разделить его на пары небольших участков, расположенных симметрично относительно центра (рис. 4). Для любых двух таких участков сумма ускорений силы тяжести, обусловленных ими, не указывает на центр диска, а отклоняется в сторону ближнего участка. В результате полная гравитация диска указывает на точку внутри ближнего полудиска. Так, если вы снова будете двигаться по четверти окружности с центром в центре диска из точки над «экватором» в точку над «полюсом» (рис. 5), угол между направлением вашего движения и направление местной силы тяжести больше, чем$90^\circ$, т. е. вы движетесь «в гору». Так$J_2$в данном случае положителен.

Для простоты расчета давайте смоделируем вытянутое тело как центральный потенциал плюс две меньшие массы к северу и югу от центра.

А давайте смоделируем сплюснутое тело как центральный потенциал плюс экваториальное кольцо точечных масс.

Если мы установим$G=1$для удобства и рассчитать гравитационный потенциал$GM/r$мы можем видеть, что он достигает пика на экваторе для сплюснутого случая и имеет минимум на экваторе для вытянутого случая.

Эти покачивания не являются чистыми синусоидами, но мы можем видеть, что их поведение в основном похоже на$\sin^2(\theta)$и$-\sin^2(\theta)$и с тех пор$J_2$- коэффициент перед таким членом (из модели геопотенциала между уравнениями 9 и 10):

мы можем видеть, что для сплюснутого тела$J_2$будет положительным, а для вытянутого тела отрицательным.

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