Chaotic Evolution of the Solar System (1992)

What kind of paper is this?

This is a computational/numerical methods paper that combines algorithm development with scientific discovery. The core contribution is both methodological (demonstrating the effectiveness of symplectic mapping for long-term orbital integration) and empirical (providing the first direct numerical confirmation that the entire Solar System exhibits chaotic dynamics). The work validates theoretical predictions through unprecedented computational experiments spanning nearly 100 million years of simulated planetary motion.

What is the motivation?

The authors aimed to address the fundamental question of the long-term stability and predictability of the Solar System. Prior work had limitations:

• Sussman & Wisdom (1988): Found chaos in Pluto’s orbit but did not integrate the full system.
• Laskar (1989): Found evidence for chaos in the whole system (excluding Pluto) but relied on analytically averaged equations, which are perturbative and truncated.

A direct integration of the full system without averaging approximations was required to validate these findings and determine whether the observed chaos was a genuine property of the planetary dynamics or an artifact of the approximation methods.

What is the novelty here?

The study represents the first direct, full-system integration spanning nearly 100 million years. Key innovations included:

• Symplectic Mapping: Application of the Wisdom-Holman mapping method, which allows for much larger time steps (e.g., 7.2 days) compared to multistep methods (which require ~100 steps/orbit) while maintaining long-term energy conservation and numerical stability.
• Custom Hardware: Use of the Supercomputer Toolkit, a multiprocessor system optimized for ODEs, where each processor was 3x faster than the entire previous generation “Digital Orrery”.
• Direct Numerical Validation: This work directly integrates Newton’s equations without analytical approximations, providing independent verification of chaotic behavior, whereas Laskar’s approach relied on secular perturbation theory.

What experiments were performed?

• 100 Myr Integration: Eight separate integrations of the entire planetary system were run for ~100 million years (reversed time) with slightly different initial conditions to measure exponential divergence of trajectories.
• Validation:
  • Compared a 3-million-year segment against the high-precision integration by Quinn, Tremaine, and Duncan (QTD).
  • Compared results with Laskar’s secular resonance angle calculations to verify consistency with the perturbative approach.
• Subsystem Analysis: Additional integrations of just the Jovian planets (outer system) and massless Pluto particles were performed to isolate the source of chaos and determine which subsystems contribute to the overall chaotic behavior.

What outcomes/conclusions?

The key results below are supported by three independent lines of evidence against numerical artifacts (see Is the Chaos Real? below).

• System-wide Chaos: The Solar System is chaotic with an exponential divergence timescale (Lyapunov time) of approximately 4 million years, meaning that initial condition uncertainties grow by a factor of $e$ roughly every 4 million years.
  • Perspective: If the Earth’s position today were shifted by just 15 meters (roughly a bus length), this microscopic error would amplify exponentially. By 100 million years, this uncertainty would span the entire orbit, meaning the Earth could be on the opposite side of the Sun from where it was predicted. The system remains deterministic while being practically unpredictable.
  • Two-segment divergence: The secular phase space divergence is not uniform. The initial segment (roughly the first 40–60 Myr) is dominated by a slower ~12 Myr timescale. The latter segment then transitions to the ~4 Myr component. The inner planets are more sensitive indicators of the faster 4 Myr process; the outer planets show the 12 Myr component for most of the 100 Myr run with the 4 Myr component only appearing in the final 5 Myr. This suggests at least two distinct mechanisms operating simultaneously.
• Jovian Chaos: The Jovian planets (Jupiter, Saturn, Uranus, Neptune) are independently chaotic. This was a surprising result: in the earlier 845 Myr Digital Orrery integration, Neptune’s orbital elements showed discrete spectral lines (unlike Pluto’s broad-band spectrum), which led the authors to dismiss Jovian chaos at that time. A new 1 billion year integration with a slightly perturbed initial position of Neptune showed exponential divergence with a timescale of only 5 million years, overturning the earlier conclusion. Subsequent Störmer predictor checks confirmed a timescale of ~19 Myr. Additional integrations spanning up to 400 Myr at various step sizes found the Lyapunov timescale varies from 3 to 30 Myr and is not a simple function of initial conditions. One specific run (step size ~0.617979 years) yielded quasiperiodic motion, confirming that nearby initial conditions can give qualitatively different long-term behavior.
• Pluto’s Robust Chaos: Pluto’s chaotic motion (10-20 Myr timescale) is robust and persists independently of whether the Jovian planets are themselves chaotic or quasiperiodic. In every outer planet integration, massless Pluto pairs diverged with a timescale between 10 and 20 Myr, regardless of how the Jovians were behaving. This robustness distinguishes Pluto’s chaotic mechanism as independent.
• Secular Resonances: Laskar identified three resonance angles ($\sigma_1$, $\sigma_2$, $\sigma_3$) that alternately librate (oscillate around a fixed value, like a pendulum) and circulate (rotate through all values, like a spinning top), and proposed resonance overlap as the driving mechanism. The authors find $\sigma_1$ and $\sigma_2$ behave consistently with Laskar, but $\sigma_3$ only circulates in their calculation, not librates and circulates as Laskar found. They also identify two additional angles, $\sigma_4 = 3(\varpi_4^\circ - \varpi_3^\circ) - 2(\Omega_4^\circ - \Omega_3^\circ)$ and $\sigma_5 = (\varpi_1^\circ - \varpi_8^\circ) + (\Omega_1^\circ - \Omega_8^\circ)$, both of which alternate between libration and circulation. The polar plot of $\sigma_1$ is more consistent with a chaotic separatrix; the $\sigma_2$ plot resembles a high-dimensional trajectory projected onto a plane, not a chaotic mechanism. Conclusion: resonance overlap is consistent with the data but not unambiguously demonstrated as the sole mechanism.
• Two Independent Mechanisms: The numerical experiments suggest at least two independent sources of chaos: one operating in the Jovian subsystem, and a separate mechanism driving Pluto’s chaos independently in the field of the Jovian planets. The simultaneous presence of two different exponential timescales in the full Solar System integration supports this. The authors speculate that Mercury (high eccentricity and inclination) may be involved in the faster mechanism, analogously to how Pluto’s high eccentricity and inclination couple eccentricity and inclination secular subsystems.
• Predictability Horizon: The chaotic nature fundamentally limits our ability to predict planetary positions beyond roughly 100 million years into the past or future, regardless of improvements in observational precision or computational power.

Is the Chaos Real or a Numerical Artifact?

The authors devote substantial effort to ruling out the possibility that the exponential divergence is a numerical artifact rather than a genuine physical property. Three lines of evidence support the conclusion that it is real:

1. Agreement across radically different methods: The 100 Myr Toolkit integration and Laskar’s secular perturbation theory approach use completely different techniques (direct integration vs. analytically averaged equations, different masses, different initial conditions, different physics). Both find consistent chaotic behavior. This cross-method agreement is strong evidence against a method-specific artifact.
2. Quasiperiodic control integration: The outer planets were integrated for 250 Myr with Uranus removed. Over this period the remaining Jovian planets showed no sign of exponential divergence. This control demonstrates that long-term numerical integrations of planetary systems do not universally produce spurious chaos; the physical configuration matters.
3. Isolated quasiperiodic run: One integration with a particular step size (near 0.617979 years) produced quasiperiodic secular phase space divergence over 300 Myr. This shows the chaos is not a universal artifact of the integration method and that some nearby initial conditions do not lead to chaotic behavior.

The authors conclude that the chaotic character is not sensitively dependent on the precise model or numerical methods used, but acknowledge that definitively ruling out subtle numerical artifacts requires identifying an unambiguous physical mechanism, which remains open.

Reproducibility Details

Data

• Initial Conditions: Derived from JPL Ephemeris DE102, matching the setup used by Quinn, Tremaine, and Duncan (QTD) for direct comparison.
• Masses: Planetary masses consistent with DE102.

Algorithms

Integrator: Symplectic n-body mapping (Wisdom & Holman method)

The Hamiltonian is split into Keplerian motion and planetary interactions:

$$H = H_{\text{Kepler}} + H_{\text{Interaction}}$$

• Structure: Uses a symplectic mapping (Kick-Drift-Kick) where the system evolves via unperturbed Keplerian orbits (Drift) punctuated by instantaneous velocity changes from planet-planet interactions (Kick). This preserves phase space volume, preventing the energy drift common in traditional integrators.
• Step Size: 7.2 days (chosen to align with QTD output timestamps for validation).
• Precision: Pseudo-quadruple precision (128-bit effectively) for position accumulation to minimize roundoff errors, though retrospectively deemed unnecessary for this problem.

Validation Integrator: Traditional linear multistep Störmer predictor (used for 22 Myr validation checks).

Models

• Hamiltonian: The system is modeled using a split Hamiltonian approach separating two-body Keplerian motion from perturbative interactions.
• General Relativity: Modeled using the potential approximation of Nobili and Roxburgh (1986). Why? General Relativity adds velocity-dependent terms that make the Keplerian part of the Hamiltonian non-integrable. To use the fast symplectic mapping (which requires an analytically solvable Kepler step), they needed a potential-only approximation that mimics GR effects (precession) without breaking the Hamiltonian splitting.
• Earth-Moon: Treated via a quadrupole approximation similar to QTD, representing the Earth-Moon system as a single body with appropriate mass distribution.

Evaluation

The primary metric for chaos was the Lyapunov Exponent, estimated via the divergence of nearby trajectories with slightly perturbed initial conditions.

Hardware

• System: Supercomputer Toolkit (MIT/Hewlett-Packard collaboration)
• Configuration: 8-processor configuration used for parallel integrations
• Performance: ~30 years of Solar System evolution per second per processor
• Total Compute: ~1000 hours of Toolkit time for the main experiments

Paper Information

Citation: Sussman, G. J., & Wisdom, J. (1992). Chaotic Evolution of the Solar System. Science, 257(5066), 56-62. https://doi.org/10.1126/science.257.5066.56

Publication: Science 1992

@article{sussmanChaoticEvolutionSolar1992,
  title = {Chaotic {{Evolution}} of the {{Solar System}}},
  author = {Sussman, Gerald J. and Wisdom, Jack},
  journal = {Science},
  volume = {257},
  number = {5066},
  pages = {56--62},
  year = {1992},
  month = {jul},
  doi = {10.1126/science.257.5066.56},
  abstract = {The evolution of the entire planetary system has been numerically integrated for a time span of nearly 100 million years. This calculation confirms that the evolution of the solar system as a whole is chaotic, with a time scale of exponential divergence of about 4 million years. Additional numerical experiments indicate that the Jovian planet subsystem is chaotic, although some small variations in the model can yield quasiperiodic motion. The motion of Pluto is independently and robustly chaotic.}
}

Additional Resources:

• Wisdom-Holman Mapping (Wikipedia)
• Digital Orrery Project