What is the sign of Earth's J₂?

Question

From Wikipedia's Geopotential model; The deviations of Earth's gravitational field from that of a homogeneous sphere and the question For the mathematical relationship between J2 (km^5/s^2) and dimensionless J2 - which one is derived from the other? I have

$$\tilde{J_n}=-\frac{J_n}{\mu R^n}$$

which has a prominent minus sign. If the dimensionless $\tilde{J_n}$ (roughly 0.001) is positive, then the dimensioned $J_n$ (roughly 1.8E+25 m⁵/s²) should be negative.

A quick test in a numerical integrator also shows me that if I don't make it negative then the nodal precession is the wrong direction.

So does the linked Wikipedia article give the sign incorrectly twice? Instead of +1.7555E+10 km⁵/s² should it be -1.7555E+10?

---

Quotes from the Wikipedia article

According to JGM-3 one therefore has that $J_{2} = 0.1082635854 \cdot 10^{-02} \cdot 6378.1363^{2} \ \cdot \ 398600.4415 \ \text{ km5/s2 } = 1.75553 \ \cdot \ 10^{10}$ and $J_{3} = -0.2532435346 \cdot 10^{-05} \cdot \ 6378.1363^{3} \cdot \ 398600.4415 \text{ km6/s2 } = -2.61913 \cdot \ 10^{11} \text{ km6/s2 }$

and

For JGM-3 the values are:

μ = 398600.440 km3⋅s−2

$J_{2}$ = 1.75553 × 1010 km5⋅s−2

$J_{3}$ = −2.61913 × 1011 km6⋅s−2

---

Quick numerical check using $J_{2}$ = +1.7555E+25 m⁵/s²

$$\omega_p = -\frac{3}{2}\frac{R_E^2}{(a(1-e^2))^2} J_2 \omega \cos i$$

gives a precession of -360 degrees in 89.4 days which matches the numerical check.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
from scipy.interpolate import interp1d
from scipy.optimize import brentq

def deriv(X, t):
    xyz, v = X.reshape(2, -1)
    acc_0 = -GMe * xyz * ((xyz**2).sum())**-1.5
    acc_2 = acc_J2(xyz)
    return np.hstack((v, acc_0+acc_2))

def acc_J2(xyz):
    x, y, z = xyz
    x2, y2, z2 = xyz**2
    r2 = (xyz**2).sum(axis=0)
    r = np.sqrt(r2)
    r7 = r**7
    u = J2 * r**-5 * 0.5 * (3*z2 - r**2)
    ax = J2 * (x/r7) * (6*z2 - 1.5*(x2+y2)) 
    ay = J2 * (y/r7) * (6*z2 - 1.5*(x2+y2)) 
    az = J2 * (z/r7) * (3*z2 - 4.5*(x2+y2))
    return np.stack([ax, ay, az], axis=0)

halfpi, pi, twopi = [f*np.pi for f in (0.5, 1, 2)]
to_degs, to_rads = 180/pi, pi/180

Re = 6378137. # meters
J2 = +1.7555E+25 # m^5/s^2
GMe = 3.986E+14 # m^3/s^2

a = 6378137. + 400000.
v0 = np.sqrt(GMe/a)

times = np.arange(0, 100*24*3600, 60)
inc = 60 * to_rads
sinc, cinc = [f(inc) for f in (np.sin, np.cos)]

X0 = np.array([a, 0, 0, 0, v0*cinc, v0*sinc])

answer, info = ODEint(deriv, X0, times, full_output=True)

x, y, z = answer.T[:3]
crossings = np.where((z[1:] > 0) * (z[:-1] <= 0))[0]
Cux = interp1d(times, x, kind='cubic', assume_sorted=True)
Cuy = interp1d(times, y, kind='cubic', assume_sorted=True)
Cuz = interp1d(times, z, kind='cubic', assume_sorted=True)
times_zc = np.array([brentq(Cuz, times[i-1], times[i+1]) for i in crossings[1:-1]])
zcs = np.array([F(times_zc) for F in (Cux, Cuy, Cuz)])
longinodes = np.arctan2(zcs[1], zcs[0])

if True:
    plt.figure()
    plt.plot(times_zc/3600/24, to_degs * longinodes)
    plt.title('400 km orbit at 60 degrees inclination \n assuming J2 = +1.7555E+25 m^5/s^2', fontsize=14)
    plt.ylabel('longitude of ascending node (deg)', fontsize=14)
    plt.xlabel('time (days)', fontsize=14)
    plt.show()