Astronomical calculations reveal the solar system's dynamical evolution, including its chaoticity, and represent the backbone of cyclostratigraphy and astrochronology. An absolute, fully calibrated astronomical time scale has hitherto been hampered beyond $\sim$50 Ma, because orbital calculations disagree before that age. Here we present geologic data and a new astronomical solution (ZB18a), showing exceptional agreement from $\sim$58 to 53 Ma. We provide a new absolute astrochronology up to 58 Ma and a new Paleocene-Eocene boundary age (56.01 $\pm$ 0.05 Ma). We show that the Paleocene-Eocene Thermal Maximum (PETM) onset occurred near a 405-kyr eccentricity maximum, suggesting an orbital trigger. We also provide an independent PETM duration (170 $\pm$ 30 kyr) from onset to recovery inflection. Our astronomical solution requires a chaotic resonance transition at $\sim$50 Ma in the solar system's fundamental frequencies.

I report results from accurate numerical integrations of Solar System orbits over the past 100Myr with the integrator package HNBody. The simulations used different integrator algorithms, step sizes, initial conditions, and included effects from general relativity, different models of the Moon, the Sun's quadrupole moment, and up to sixteen asteroids. I also probed the potential effect of a hypothetical Planet 9, using one set of possible orbital elements. The most expensive integration (Bulirsch-Stoer) required 4~months wall-clock time with a maximum relative energy error <~3e{-13}. The difference in Earth's eccentricity (DeE) was used to track the difference between two solutions, considered to diverge at time tau when max|DeE| irreversibly crossed ~10\% of mean eE (~0.028x0.1). The results indicate that finding a unique orbital solution is limited by initial conditions from current ephemerides and asteroid perturbations to ~54Myr. Bizarrely, the 4-month Bulirsch-Stoer integration and a symplectic integration that required only 5~hours wall-clock time (12-day time step, Moon as a simple quadrupole perturbation), agree to ~63Myr. Internally, such symplectic integrations are remarkably consistent even for large time steps, suggesting that the relationship between time step and tau is not a robust indicator for the absolute accuracy of symplectic integrations. The effect of a hypothetical Planet~9 on DeE becomes discernible at ~65Myr. Using tau as a criterion, the current state-of-the-art solutions all differ from previously published results beyond ~50Myr. I also conducted an eigenmode analysis, which provides some insight into the chaotic nature of the inner Solar System. The current study provides new orbital solutions for applications in geological studies.

Due to the chaotic nature of the Solar System, the question of its dynamic long-term stability can only be answered in a statistical sense, e.g. based on numerical ensemble integrations of nearby orbits. Destabilization, including catastrophic encounters and/or collisions involving the Earth, has been suggested to be initiated through a large increase in Mercury's eccentricity (eM), with an estimated probability of ~1%. However, it has recently been shown that the statistics of numerical Solar System integrations are sensitive to the accuracy and type of numerical algorithm. Here I report results from computationally demanding ensemble integrations (N=1,600 with slightly different initial conditions) at unprecedented accuracy based on the full equations of motion of the eight planets and Pluto over 5Gyr, including contributions from general relativity. The standard symplectic algorithm produced spurious results for highly eccentric orbits and during close encounters, which were hence integrated with a suitable Bulirsch-Stoer algorithm, specifically designed for these situations. The present study yields odds for a large increase in Mercury's eccentricity that are less than previous estimates. Strikingly, in two solutions Mercury continued on highly eccentric orbits (after reaching eM values >0.93) for 80-100Myr before colliding with Venus or the Sun. Most importantly, none of the 1,600 solutions led to a close encounter involving the Earth or a destabilization of Earth's orbit in the future. I conclude that Earth's orbit is dynamically highly stable for billions of years, despite the chaotic behavior of the Solar System.

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.