  Due to the chaotic nature of the Solar System, the question of its long-term
stability can only be answered in a statistical sense, for instance, based on
numerical ensemble integrations of nearby orbits. Destabilization of the inner
planets, leading to close encounters and/or collisions can be initiated through
a large increase in Mercury's eccentricity, with a currently assumed likelihood
of ~1%. However, little is known at present about the robustness of this
number. Here I report ensemble integrations of the full equations of motion of
the eight planets and Pluto over 5 Gyr, including contributions from general
relativity. The results show that different numerical algorithms lead to
statistically different results for the evolution of Mercury's eccentricity
(eM). For instance, starting at present initial conditions (eM ~= 0.21),
Mercury's maximum eccentricity achieved over 5 Gyr is on average significantly
higher in symplectic ensemble integrations using heliocentricthan Jacobi
coordinates and stricter error control. In contrast, starting at a possible
future configuration (eM ~= 0.53), Mercury's maximum eccentricity achieved over
the subsequent 500 Myr is on average significantly lower using heliocentric
than Jacobi coordinates. For example, the probability for eM to increase beyond
0.53 over 500 Myr is >90% (Jacobi) vs. only 40-55% (heliocentric). This poses a
dilemma as the physical evolution of the real system - and its probabilistic
behavior - cannot depend on the coordinate system or numerical algorithm chosen
to describe it. Some tests of the numerical algorithms suggest that symplectic
integrators using heliocentric coordinates underestimate the odds for
destabilization of Mercury's orbit at high initial eM.
