Urbain Le Verrier entered the École Polytechnique in 1831 and initially pursued chemistry. By the late 1830s he had turned to celestial mechanics, publishing memoirs on the secular perturbations of planetary orbits and the stability of the solar system.
His early contributions extended Laplace’s Mécanique Céleste, improving the convergence of perturbation series and producing new tables for the Connaissance des Temps. These seemingly routine calculations built the precise mathematical toolkit that would later make the Neptune prediction possible.
Every planet pulls on every other planet. In a solar system with only one planet, Kepler’s laws give an exact elliptical orbit. Add a second planet and the orbits are no longer exactly elliptical — each planet perturbs the other.
The key insight of perturbation theory: if the perturbing planet is much less massive than the Sun, its effect is small and can be computed as a correction to the Keplerian orbit. The correction is expressed as a series expansion in the mass ratio \(m'/M_\odot\).
The gravitational potential of a perturbing planet can be expanded in a Fourier series involving the orbital elements of both planets. For circular, coplanar orbits this simplifies to:
\[ R = \frac{G m'}{a'} \sum_{j=-\infty}^{\infty} b_{1/2}^{(j)}(\alpha) \cos(j(\lambda - \lambda')) \]
where \(\alpha = a/a'\) is the ratio of semi-major axes, \(\lambda\) and \(\lambda'\) are the mean longitudes, and \(b_{1/2}^{(j)}(\alpha)\) are the Laplace coefficients — integrals that Le Verrier and his contemporaries computed by hand.
The Laplace coefficient is defined as:
\[ b_s^{(j)}(\alpha) = \frac{1}{\pi} \int_0^{2\pi} \frac{\cos(j\theta)}{(1 - 2\alpha\cos\theta + \alpha^2)^s} \, d\theta \]
Le Verrier evaluated these by expanding the integrand in power series and summing term by term — tedious but tractable with patience. We can compute them numerically in a line of code.
import numpy as np
from discoverneptune.perturbation import laplace_coefficient
# Laplace coefficients for α = a_Jupiter/a_Saturn ≈ 0.545
alpha_js = 5.2 / 9.5
for j in range(5):
b = laplace_coefficient(0.5, j, alpha_js)
print(f" b_{{1/2}}^{{({j})}}({alpha_js:.3f}) = {b:.6f}") b_{1/2}^{(0)}(0.547) = 2.181903
b_{1/2}^{(1)}(0.547) = 0.623691
b_{1/2}^{(2)}(0.547) = 0.259915
b_{1/2}^{(3)}(0.547) = 0.119477
b_{1/2}^{(4)}(0.547) = 0.057494The \(j=0\) term gives the secular (long-term average) perturbation. The \(j=1\) term drives the dominant periodic perturbation at the synodic frequency — the beat frequency between two planets’ orbital periods.
We can test the first-order analytical prediction against a full REBOUND N-body integration. To keep the comparison honest, we use circular coplanar orbits, removing every eccentricity- and inclination-driven term — what remains is the synodic perturbation the first-order theory predicts, plus the near-resonant amplification it cannot capture.
Consider a toy system: a “Jupiter” at 5.2 AU and a “Saturn” at 9.5 AU, both on circular orbits. The analytical perturbation theory predicts a periodic longitude perturbation on Jupiter at the synodic period.
from discoverneptune.perturbation import (
first_order_longitude_perturbation,
predict_longitude_perturbations,
)
# Saturn's perturbation on Jupiter (circular coplanar)
m_saturn = 2.858e-4 # solar masses
a_jup, a_sat = 5.2, 9.5
P_jup, P_sat = 11.86, 29.46
amplitude_rad, synodic_period = first_order_longitude_perturbation(
m_perturber_solar=m_saturn,
a_inner_au=a_jup,
a_outer_au=a_sat,
period_inner_yr=P_jup,
period_outer_yr=P_sat,
)
print(f"Predicted perturbation amplitude: {np.degrees(amplitude_rad) * 3600:.1f} arcsec")
print(f"Synodic period: {synodic_period:.2f} years")Predicted perturbation amplitude: 33.7 arcsec
Synodic period: 19.85 yearsNow compare against a REBOUND simulation of the same circular system:
import rebound
import matplotlib.pyplot as plt
# Build a circular coplanar Jupiter-Saturn system in REBOUND
sim = rebound.Simulation()
sim.units = ("yr", "AU", "Msun")
sim.add(m=1.0) # Sun
sim.add(m=9.547e-4, a=5.2, e=0.0, inc=0.0, hash="Jupiter")
sim.add(m=2.858e-4, a=9.5, e=0.0, inc=0.0, hash="Saturn")
sim.move_to_com()
sim.integrator = "whfast"
sim.dt = 0.05
# Integrate and record Jupiter's heliocentric longitude
times = np.linspace(0, 200, 2000)
jup_lon = np.empty(len(times))
for i, t in enumerate(times):
sim.integrate(t)
sun = sim.particles[0]
p = sim.particles["Jupiter"]
jup_lon[i] = np.arctan2(p.y - sun.y, p.x - sun.x)
# Subtract the mean motion to isolate the perturbation
jup_lon_unwrapped = np.unwrap(jup_lon)
mean_motion_fit = np.polyfit(times, jup_lon_unwrapped, 1)
jup_perturbation_arcsec = (
(jup_lon_unwrapped - np.polyval(mean_motion_fit, times)) * np.degrees(1) * 3600
)
# Analytical prediction
analytical_perturbation_arcsec = (
predict_longitude_perturbations(
times, m_saturn, a_jup, a_sat, P_jup, P_sat
) * np.degrees(1) * 3600
)
# The raw detrended signal is dominated by the Jupiter–Saturn "Great
# Inequality" (the slow, large beat from their near-5:2 period ratio), which
# first-order theory deliberately omits. Isolate the synodic-frequency band
# for a like-for-like comparison.
dt = times[1] - times[0]
freqs = np.fft.rfftfreq(len(times), dt)
spectrum = np.fft.rfft(jup_perturbation_arcsec)
f_synodic = 1.0 / synodic_period
band = (freqs > 0.5 * f_synodic) & (freqs < 1.5 * f_synodic)
jup_synodic_arcsec = np.fft.irfft(np.where(band, spectrum, 0.0), n=len(times))def _build(theme="light"):
apply_style(theme)
pal = palette(theme)
cycle = category_cycle(theme)
fig, ax = plt.subplots(figsize=fig_size("wide"))
ax.plot(
times, jup_perturbation_arcsec,
color=cycle[0], lw=0.8, alpha=0.25,
label="REBOUND — full signal (incl. Great Inequality)",
)
ax.plot(
times, jup_synodic_arcsec,
color=cycle[0], lw=1.2,
label="REBOUND — synodic component",
)
ax.plot(
times, analytical_perturbation_arcsec, "--",
color=cycle[1], lw=1.2,
label="First-order theory",
)
ax.set_xlabel("Time (years)")
ax.set_ylabel("Jupiter longitude perturbation (arcsec)")
ax.set_title("Analytical vs Numerical Perturbation Theory")
ax.legend()
return fig
display(dual_render(_build, alt="First-order theory vs the synodic-frequency component of a REBOUND integration, with the full Great-Inequality-dominated signal faint in the background"))At the synodic frequency the comparison becomes meaningful: both curves oscillate at the ~20-year Jupiter–Saturn beat, but the first-order theory still underpredicts the numerical amplitude by roughly a factor of three — Jupiter and Saturn sit near a 5:2 period ratio, and the near-resonance amplifies the perturbation beyond what the lowest-order expansion captures. The full numerical signal is larger still: the resulting “Great Inequality” — a slow beat with a period near 900 years — dominates the raw perturbation. Capturing that required pushing the expansions to higher order, which is exactly the hand computation Laplace and Le Verrier mastered.
By the early 1840s, Le Verrier had published refined orbital tables for several planets. This work established two things: his mastery of perturbation theory at a practical level, and his sensitivity to small residuals between theory and observation. When François Arago pointed him toward the Uranus anomaly in 1845, Le Verrier had exactly the right tools for the job.