Mercury, Vulcan, and General Relativity

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

Year	Anomaly	Prediction	Outcome
1846	Uranus 31.5 arcsec RMS (131.4 peak-to-peak)	New outer planet at ~36 AU	Neptune found ✓
1859	Mercury ~38 arcsec/century reported by Le Verrier (later ~43)	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 — 6-body (Sun + Mercury + Venus + Earth + Jupiter + Saturn) 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,
    dual_render,
    fig_size,
    palette,
)

apply_style()

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

from discoverneptune.data import (
    MERCURY_EPOCH_START,
    MERCURY_EPOCH_STEP,
    MERCURY_EPOCH_STOP,
    PLANET_IDS,
    fetch_planet_vectors,
    find_bundled_data_dir,
)
from discoverneptune.mercury import (
    analytical_gr_precession,
    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

The committed run loads deterministic JPL Horizons-derived J2000 state vectors from data/inner_planet_state_vectors_j2000.csv; set USE_LIVE_HORIZONS = True only when intentionally refreshing them. Venus and Earth dominate Mercury’s classical precession; Jupiter and Saturn supply additional secular perturbations.

Code

# ── Load 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"]),
    )

USE_LIVE_HORIZONS = False
if USE_LIVE_HORIZONS:
    kw = dict(
        start=MERCURY_EPOCH_START,
        stop=MERCURY_EPOCH_STOP,
        step=MERCURY_EPOCH_STEP,
        use_cache=False,
    )
    vectors = {
        planet: fetch_planet_vectors(PLANET_IDS[planet], **kw)
        for planet in ("mercury", "venus", "earth", "jupiter", "saturn")
    }
    vector_source = "Using live JPL Horizons J2000 state vectors"
else:
    bundled_path = find_bundled_data_dir() / "inner_planet_state_vectors_j2000.csv"
    bundled = pd.read_csv(bundled_path)
    vectors = {
        planet: bundled.loc[bundled["planet"] == planet]
        for planet in ("mercury", "venus", "earth", "jupiter", "saturn")
    }
    vector_source = "Using bundled J2000 state vectors"

mercury_sv = _sv_from_df(vectors["mercury"])
venus_sv = _sv_from_df(vectors["venus"])
earth_sv = _sv_from_df(vectors["earth"])
jupiter_sv = _sv_from_df(vectors["jupiter"])
saturn_sv = _sv_from_df(vectors["saturn"])

print(vector_source)
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")

Using bundled J2000 state vectors
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. Le Verrier reported about 38 arcsec/century in 1859; later analyses refined the discrepancy before Einstein’s general relativity explained it 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("(Modern residual estimate: ~43.0 arcsec/century; Le Verrier reported ~38)")

Textbook GR precession:  42.98 arcsec/century
(Modern residual estimate: ~43.0 arcsec/century; Le Verrier reported ~38)

Newtonian integration (500 years)

We integrate Mercury’s orbit for 500 years under purely Newtonian gravity with Venus, Earth, Jupiter, and Saturn as perturbers. The longitude of perihelion (ω̃) precesses at ~525 arcsec/century in this six-body run. Adding the relativistic correction supplies the remaining ~43 arcsec/century.

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",         "31.5 RMS; 131.4 peak-to-peak", "~38 reported; later ~43"),
    ("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 (20.2 % error)",   "N/A — no planet exists"),
    ("Our simulation",  "See current JPL-mode run", 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               31.5 RMS; 131.4 peak-to-peak       ~38 reported; later ~43
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 (20.2 % error)        N/A — no planet exists
Our simulation            See current JPL-mode run   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 — far below this comparison's residual scale.")
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 — far below this comparison's residual scale.
  There, a new planet was the correct (and only) explanation.
============================================================