We combine analytical understanding of resonant dynamics in two-planet systems with machine learning techniques to train a model capable of robustly classifying stability in compact multi-planet systems over long timescales of $10^9$ orbits. Our Stability of Planetary Orbital Configurations Klassifier (SPOCK) predicts stability using physically motivated summary statistics measured in integrations of the first $10^4$ orbits, thus achieving speed-ups of up to $10^5$ over full simulations. This computationally opens up the stability constrained characterization of multi-planet systems. Our model, trained on $\approx 100,000$ three-planet systems sampled at discrete resonances, generalizes both to a sample spanning a continuous period-ratio range, as well as to a large five-planet sample with qualitatively different configurations to our training dataset. Our approach significantly outperforms previous methods based on systems' angular momentum deficit, chaos indicators, and parametrized fits to numerical integrations. We use SPOCK to constrain the free eccentricities between the inner and outer pairs of planets in the Kepler-431 system of three approximately Earth-sized planets to both be below 0.05. Our stability analysis provides significantly stronger eccentricity constraints than currently achievable through either radial velocity or transit duration measurements for small planets, and within a factor of a few of systems that exhibit transit timing variations (TTVs). Given that current exoplanet detection strategies now rarely allow for strong TTV constraints (Hadden et al., 2019), SPOCK enables a powerful complementary method for precisely characterizing compact multi-planet systems. We publicly release SPOCK for community use.

Instabilities in compact planetary systems are generically driven by chaotic dynamics. This implies that an instability time measured through direct N-body integration is not exact, but rather represents a single draw from a distribution of equally valid chaotic trajectories. In order to characterize the "errors" on reported instability times from direct N-body integrations, we investigate the shape and parameters of the instability time distributions (ITDs) for ensembles of shadow trajectories that are initially perturbed from one another near machine precision. We find that in the limit where instability times are long compared to the Lyapunov (chaotic) timescale, ITDs approach remarkably similar lognormal distributions with standard deviations ~0.43 $\pm$ 0.16 dex, despite the instability times varying across our sample from $10^4-10^8$ orbits. We find excellent agreement between these predictions, derived from ~450 closely packed configurations of three planets, and a much wider validation set of ~10,000 integrations, as well as on ~20,000 previously published integrations of tightly packed five-planet systems, and a seven-planet resonant chain based on TRAPPIST-1, despite their instability timescales extending beyond our analyzed timescale. We also test the boundary of applicability of our results on dynamically excited versions of our Solar System. These distributions define the fundamental limit imposed by chaos on the predictability of instability times in such planetary systems. It provides a quantitative estimate of the intrinsic error on an N-body instability time imprinted by chaos, approximately a factor of 3 in either direction.