We describe a software package called VPLanet that simulates fundamental aspects of planetary system evolution over Gyr timescales, with a focus on investigating habitable worlds. In this initial release, eleven physics modules are included that model internal, atmospheric, rotational, orbital, stellar, and galactic processes. Many of these modules can be coupled simultaneously to simulate the evolution of terrestrial planets, gaseous planets, and stars. The code is validated by reproducing a selection of observations and past results. VPLanet is written in C and designed so that the user can choose the physics modules to apply to an individual object at runtime without recompiling, i.e., a single executable can simulate the diverse phenomena that are relevant to a wide range of planetary and stellar systems. This feature is enabled by matrices and vectors of function pointers that are dynamically allocated and populated based on user input. The speed and modularity of VPLanet enables large parameter sweeps and the versatility to add/remove physical phenomena to assess their importance. VPLanet is publicly available from a repository that contains extensive documentation, numerous examples, Python scripts for plotting and data management, and infrastructure for community input and future development.

Using an energy balance model with ice sheets, we examine the climate response of an Earth-like planet orbiting a G dwarf star and experiencing large orbital and obliquity variations. We find that ice caps couple strongly to the orbital forcing, leading to extreme ice ages. In contrast with previous studies, we find that such exo-Milankovitch cycles tend to impair habitability by inducing snowball states within the habitable zone. The large amplitude changes in obliquity and eccentricity cause the ice edge, the lowest latitude extent of the ice caps, to become unstable and grow to the equator. We apply an analytical theory of the ice edge latitude to show that obliquity is the primary driver of the instability. The thermal inertia of the ice sheets and the spectral energy distribution of the G dwarf star increase the sensitivity of the model to triggering runaway glaciation. Finally, we apply a machine learning algorithm to demonstrate how this technique can be used to extend the power of climate models. This work illustrates the importance of orbital evolution for habitability in dynamically rich planetary systems. We emphasize that as potentially habitable planets are discovered around G dwarfs, we need to consider orbital dynamics.

The obliquity of the Earth, which controls our seasons, varies by only ~2.5 degrees over ~40,000 years, and its eccentricity varies by only ~0.05 over 100,000 years. Nonetheless, these small variations influence Earth's ice ages. For exoplanets, however, variations can be significantly larger. Previous studies of the habitability of moonless Earth-like exoplanets have found that high obliquities, high eccentricities, and dynamical variations can extend the outer edge of the habitable zone by preventing runaway glaciation (snowball states). We expand upon these studies by exploring the orbital dynamics with a semi-analytic model that allows us to map broad regions of parameter space. We find that in general, the largest drivers of obliquity variations are secular spin-orbit resonances. We show how the obliquity varies in several test cases, including Kepler-62 f, across a wide range of orbital and spin parameters. These obliquity variations, alongside orbital variations, will have a dramatic impact on the climates of such planets.