Mercury, Vulcan, and General Relativity

Urbain Le Verrier applied the same inverse-problem method twice:

Year	Anomaly	Prediction	Outcome
1846	Uranus ~133 arcsec RMS	New outer planet at ~36 AU	Neptune found ✓
1859	Mercury ~43 arcsec/century	New inner planet (Vulcan) at ~0.14 AU	Never found ✗

The companion notebooks reproduce the Neptune prediction computationally. This notebook demonstrates the Mercury case:

Newtonian simulation — 5-body (Sun + Mercury + Venus + Earth + Jupiter) WHFast integration; measure pomega (longitude of perihelion) drift rate.
GR-corrected simulation — same system with reboundx post-Newtonian corrections.
The GR contribution should match the textbook formula: ≈ 43 arcsec/century.
Vulcan analysis — show what mass and detectability a hypothetical inner planet would require.
Side-by-side comparison — Neptune success vs. Vulcan failure.

Both 500-year integrations take ~30 seconds total (WHFast, dt = 0.001 yr ≈ 9 hours).

Setup

Code

from IPython.display import display
from discoverneptune.plot_style import (
    apply_style, palette, category_cycle, dual_render,
    SEQUENTIAL_CMAP, truth_line, annotate_stat, fig_size,
)
apply_style()

import matplotlib.pyplot as plt
import numpy as np

from discoverneptune.data import (
    MERCURY_EPOCH_START,
    MERCURY_EPOCH_STOP,
    MERCURY_EPOCH_STEP,
    PLANET_IDS,
    fetch_planet_vectors,
)
from discoverneptune.mercury import (
    MERCURY_A,
    MERCURY_PERIOD_YR,
    analytical_gr_precession,
    build_mercury_simulation,
    fit_precession_rate,
    measure_pomega_series,
    vulcan_mass_estimate,
)
from discoverneptune.simulation import StateVector

%matplotlib inline

print("Imports OK")

Imports OK

Fetch J2000 state vectors

We fetch initial conditions for Mercury, Venus, Earth, and Jupiter from JPL Horizons at J2000.0. Venus and Earth are included because their secular perturbations dominate Mercury’s classical precession (~532 arcsec/century from planetary perturbations alone).

Code

# ── Fetch J2000 state vectors for Mercury, Venus, Earth, Jupiter, Saturn ────

def _sv_from_df(df) -> StateVector:
    """Extract a StateVector from row 0 of a JPL Horizons DataFrame."""
    row = df.iloc[0]
    return StateVector(
        x=float(row["x"]),   y=float(row["y"]),   z=float(row["z"]),
        vx=float(row["vx"]), vy=float(row["vy"]), vz=float(row["vz"]),
    )

kw = dict(start=MERCURY_EPOCH_START, stop=MERCURY_EPOCH_STOP, step=MERCURY_EPOCH_STEP)

mercury_sv = _sv_from_df(fetch_planet_vectors(PLANET_IDS["mercury"], **kw))
venus_sv   = _sv_from_df(fetch_planet_vectors(PLANET_IDS["venus"],   **kw))
earth_sv   = _sv_from_df(fetch_planet_vectors(PLANET_IDS["earth"],   **kw))
jupiter_sv = _sv_from_df(fetch_planet_vectors(PLANET_IDS["jupiter"], **kw))
saturn_sv  = _sv_from_df(fetch_planet_vectors(PLANET_IDS["saturn"],  **kw))

print(f"Mercury  x={mercury_sv.x:.4f} AU, y={mercury_sv.y:.4f} AU")
print(f"Venus    x={venus_sv.x:.4f} AU, y={venus_sv.y:.4f} AU")
print(f"Earth    x={earth_sv.x:.4f} AU, y={earth_sv.y:.4f} AU")
print(f"Jupiter  x={jupiter_sv.x:.4f} AU, y={jupiter_sv.y:.4f} AU")
print(f"Saturn   x={saturn_sv.x:.4f} AU, y={saturn_sv.y:.4f} AU")

Mercury  x=-0.1407 AU, y=-0.4439 AU
Venus    x=-0.7186 AU, y=-0.0225 AU
Earth    x=-0.1685 AU, y=0.9688 AU
Jupiter  x=4.0035 AU, y=2.9354 AU
Saturn   x=6.4086 AU, y=6.5680 AU

Textbook GR precession formula

Einstein’s formula gives the excess perihelion advance per orbit: Δφ = 6πGM/(ac²(1−e²)). For Mercury this predicts ~43 arcsec/century — the famous anomaly that Le Verrier identified in 1859 and that Einstein’s general relativity explained in 1915.

Code

# ── Textbook GR formula ───────────────────────────────────────────────────────
#
#  Δφ/orbit = 6π G M_sun / (a_M c² (1−e_M²))
#
#  In REBOUND units (G = 4π²  AU³ yr⁻² M_sun⁻¹):
#  Δφ/orbit = 6π × 4π² / (a_M × c²_AU_per_yr × (1−e_M²))
#
#  Rate = Δφ/orbit × (100 yr/century / T_Mercury) × (180°/π × 3600 arcsec/°)

gr_formula = analytical_gr_precession()
print(f"Textbook GR precession:  {gr_formula:.2f} arcsec/century")
print(f"(Le Verrier's anomaly:  ~43.0 arcsec/century — literature value)")

Textbook GR precession:  42.98 arcsec/century
(Le Verrier's anomaly:  ~43.0 arcsec/century — literature value)

Newtonian integration (500 years)

We integrate Mercury’s orbit for 500 years under purely Newtonian gravity with Venus, Earth, and Jupiter as perturbers. The longitude of perihelion (ω̃) precesses at ~517 arcsec/century from classical perturbations. The 43 arcsec/century anomaly is the gap between this rate and the observed value.

Code

# ── Newtonian integration: 500 years ─────────────────────────────────────────
print("Running Newtonian simulation (500 yr) …")
times_n, pomega_n = measure_pomega_series(
    mercury_sv, venus_sv, earth_sv, jupiter_sv,
    saturn_sv=saturn_sv,
    n_years=500,
    sample_every_yr=1.0,
    enable_gr=False,
)
rate_newt = fit_precession_rate(times_n, pomega_n)
print(f"Newtonian precession rate: {rate_newt:.2f} arcsec/century")

Running Newtonian simulation (500 yr) …
Newtonian precession rate: 524.69 arcsec/century

GR-corrected integration (500 years)

The same integration with reboundx post-Newtonian corrections enabled. The difference between the GR-corrected and Newtonian rates isolates general relativity’s contribution.

Code

# ── GR-corrected integration: 500 years ──────────────────────────────────────
print("Running GR-corrected simulation (500 yr) …")
times_g, pomega_g = measure_pomega_series(
    mercury_sv, venus_sv, earth_sv, jupiter_sv,
    saturn_sv=saturn_sv,
    n_years=500,
    sample_every_yr=1.0,
    enable_gr=True,
)
rate_gr = fit_precession_rate(times_g, pomega_g)
gr_contribution = rate_gr - rate_newt

print(f"GR-corrected rate:    {rate_gr:.2f} arcsec/century")
print(f"Newtonian rate:       {rate_newt:.2f} arcsec/century")
print(f"GR contribution:      {gr_contribution:.2f} arcsec/century")
print(f"Textbook prediction:  {gr_formula:.2f} arcsec/century")

Running GR-corrected simulation (500 yr) …
GR-corrected rate:    567.65 arcsec/century
Newtonian rate:       524.69 arcsec/century
GR contribution:      42.96 arcsec/century
Textbook prediction:  42.98 arcsec/century

Perihelion precession plots

Left panel: the raw longitude of perihelion over 500 years, showing the Newtonian and GR-corrected curves diverging. Right panel: the GR−Newtonian difference grows linearly at ~43 arcsec/century, matching Einstein’s prediction.

Code

# ── Plot: pomega drift — GR minus Newtonian ──────────────────────────────────
delta_arcsec = (np.unwrap(pomega_g) - np.unwrap(pomega_n)) * np.degrees(1.0) * 3600.0

def _build_precession(delta_arcsec, times_n, pomega_n, times_g, pomega_g, gr_contribution, *, theme="light"):
    apply_style(theme)
    pal = palette(theme)
    fig, axes = plt.subplots(1, 2, figsize=fig_size("double"))

    # Left: raw pomega series
    ax = axes[0]
    ax.plot(times_n, np.degrees(np.unwrap(pomega_n)), color=pal.before, label="Newtonian", lw=1.0)
    ax.plot(times_g, np.degrees(np.unwrap(pomega_g)), color=pal.after, label="GR-corrected", lw=1.0, ls="--")
    ax.set_xlabel("Time (yr)")
    ax.set_ylabel("ω̃ — longitude of perihelion (°)")
    ax.set_title("Perihelion precession: Newtonian vs GR")
    ax.legend()

    # Right: GR − Newtonian difference
    ax = axes[1]
    ax.plot(times_g, delta_arcsec, color=pal.after, lw=1.0)
    # Overlay the fitted linear trend
    trend = gr_contribution / 100.0 * times_g   # arcsec, since gr_contribution is per century
    ax.plot(times_g, trend, color=pal.truth, ls="--", lw=0.8, label=f"Fit: {gr_contribution:.1f} arcsec/century")
    ax.set_xlabel("Time (yr)")
    ax.set_ylabel("Δω̃ = GR − Newtonian (arcsec)")
    ax.set_title("GR contribution to Mercury precession")
    ax.legend()

    return fig

display(dual_render(
    _build_precession,
    delta_arcsec, times_n, pomega_n, times_g, pomega_g, gr_contribution,
    alt="Mercury perihelion precession: Newtonian vs GR-corrected",
))
print("Figure rendered.")

Figure rendered.

Vulcan: the planet that wasn’t

If a hypothetical inner planet (Vulcan) were responsible for Mercury’s anomaly instead of GR, how massive would it need to be? First-order secular perturbation theory gives the required mass as a function of orbital distance. The result: any such planet would be far too massive to have escaped detection during solar observations and transits.

Code

# ── Vulcan table ──────────────────────────────────────────────────────────────
#
# What mass would a hypothetical inner planet (Vulcan) need to explain
# Mercury's 43 arcsec/century anomaly?
#
# Formula (first-order secular perturbation, inner perturber):
#   dω/dt ≈ (3/4) × n_M × (m_V / M_sun) × (a_V / a_M)²
#
# Solving for m_V:  m_V = dω/dt × (4/3) / (n_M × (a_V/a_M)²)

M_EARTH_IN_SOLAR = 1.0 / 332_946.0  # solar masses per Earth mass

print(f"{'a_Vulcan (AU)':>14}  {'Mass (M_sun)':>14}  {'Mass (M_Earth)':>14}  {'Detectability':>30}")
print("-" * 80)

for a_v in [0.10, 0.15, 0.20, 0.25, 0.30]:
    m_sun = vulcan_mass_estimate(a_v)
    m_earth = m_sun / M_EARTH_IN_SOLAR

    # Rough detectability: Venus (0.72 AU) has max elongation ~47°
    # Mercury (0.387 AU) has max elongation ~28°
    # Inner planets would be observable during transits / solar observations
    if m_earth > 0.1:
        note = "≫ Moon — easily detected by transits"
    elif m_earth > 0.001:
        note = "Moon-class — detectable by transits"
    else:
        note = "Sub-lunar — still no transit observed"

    print(f"{a_v:>14.2f}  {m_sun:>14.4e}  {m_earth:>14.3f}  {note:>30}")

print()
print("Conclusion: any planet massive enough to produce 43 arcsec/century")
print("would be far too large to have escaped detection. GR is the only answer.")

 a_Vulcan (AU)    Mass (M_sun)  Mass (M_Earth)                   Detectability
--------------------------------------------------------------------------------
          0.10      1.5966e-06           0.532  ≫ Moon — easily detected by transits
          0.15      7.0960e-07           0.236  ≫ Moon — easily detected by transits
          0.20      3.9915e-07           0.133  ≫ Moon — easily detected by transits
          0.25      2.5546e-07           0.085  Moon-class — detectable by transits
          0.30      1.7740e-07           0.059  Moon-class — detectable by transits

Conclusion: any planet massive enough to produce 43 arcsec/century
would be far too large to have escaped detection. GR is the only answer.

Neptune vs Vulcan: two predictions, two outcomes

A side-by-side comparison of Le Verrier’s two uses of the inverse method. For Uranus, the anomaly was caused by a real planet and the method succeeded. For Mercury, the anomaly was caused by spacetime curvature — no planet could explain it.

Code

# ── Side-by-side comparison table ────────────────────────────────────────────

print("\n" + "=" * 70)
print("  Le Verrier's Two Predictions: Uranus / Neptune vs Mercury / Vulcan")
print("=" * 70)
print()

rows = [
    ("Anomaly",         "~133 arcsec RMS",           "~43 arcsec/century"),
    ("Observed body",   "Uranus (19.2 AU)",           "Mercury (0.39 AU)"),
    ("Le Verrier's",    "New outer planet",           "New inner planet (Vulcan)"),
    ("prediction",      "at ~36 AU",                  "at ~0.14 AU"),
    ("True explanation","New planet — Neptune",       "General Relativity (GR)"),
    ("Prediction result","Neptune found 1846 ✓",     "Vulcan — never found ✗"),
    ("True distance",   "30.07 AU (2.4 % error)",    "N/A — no planet exists"),
    ("Our simulation",  "a ~32.7 AU, RMS ~2.5 arcsec", f"GR contrib = {gr_contribution:.1f} arcsec/cy"),
]

col_w = 28
header = f"{'':20}  {'Uranus / Neptune':>{col_w}}  {'Mercury / GR':>{col_w}}"
print(header)
print("-" * len(header))
for label, col1, col2 in rows:
    print(f"{label:20}  {col1:>{col_w}}  {col2:>{col_w}}")
print()
print("The same method — predict a body from orbital anomalies — succeeded for")
print("Uranus and failed for Mercury. The difference: GR is not negligible at 0.39 AU.")


======================================================================
  Le Verrier's Two Predictions: Uranus / Neptune vs Mercury / Vulcan
======================================================================

                                  Uranus / Neptune                  Mercury / GR
--------------------------------------------------------------------------------
Anomaly                            ~133 arcsec RMS            ~43 arcsec/century
Observed body                     Uranus (19.2 AU)             Mercury (0.39 AU)
Le Verrier's                      New outer planet     New inner planet (Vulcan)
prediction                               at ~36 AU                   at ~0.14 AU
True explanation              New planet — Neptune       General Relativity (GR)
Prediction result             Neptune found 1846 ✓        Vulcan — never found ✗
True distance               30.07 AU (2.4 % error)        N/A — no planet exists
Our simulation         a ~32.7 AU, RMS ~2.5 arcsec   GR contrib = 43.0 arcsec/cy

The same method — predict a body from orbital anomalies — succeeded for
Uranus and failed for Mercury. The difference: GR is not negligible at 0.39 AU.

Summary

Code

# ── Summary ───────────────────────────────────────────────────────────────────

print("\n" + "=" * 60)
print("  Summary")
print("=" * 60)
print(f"  Textbook GR formula:        {gr_formula:.2f} arcsec/century")
print(f"  Newtonian simulation rate:  {rate_newt:.1f} arcsec/century")
print(f"  GR-corrected rate:          {rate_gr:.1f} arcsec/century")
print(f"  GR contribution (sim):      {gr_contribution:.1f} arcsec/century")
print(f"  GR contribution (formula):  {gr_formula:.1f} arcsec/century")
print()
print("  → GR closes the Mercury anomaly exactly.")
print("  → No planet (Vulcan) is needed or consistent with observations.")
print()
print("  This contrasts with Neptune: GR corrections at 19+ AU are only")
print("  ~0.01 arcsec/century — five orders of magnitude below the noise floor.")
print("  There, a new planet was the correct (and only) explanation.")
print("=" * 60)


============================================================
  Summary
============================================================
  Textbook GR formula:        42.98 arcsec/century
  Newtonian simulation rate:  524.7 arcsec/century
  GR-corrected rate:          567.6 arcsec/century
  GR contribution (sim):      43.0 arcsec/century
  GR contribution (formula):  43.0 arcsec/century

  → GR closes the Mercury anomaly exactly.
  → No planet (Vulcan) is needed or consistent with observations.

  This contrasts with Neptune: GR corrections at 19+ AU are only
  ~0.01 arcsec/century — five orders of magnitude below the noise floor.
  There, a new planet was the correct (and only) explanation.
============================================================