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 to simultaneously 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.
We investigate planetesimal accretion via a direct N-body simulation of an annulus at 1 au orbiting a 1 M? star. The planetesimal ring, which initially contains N = 106 bodies is evolved into the oligarchic growth phase. Unlike previous lower resolution studies, we find that the mass distribution of planetesimals develops a bump at intermediate mass after the oligarchs form. This feature marks a boundary between growth modes. The smallest planetesimals are packed tightly enough together to populate mean motion resonances with the oligarchs, which heats the small bodies, enhancing their growth. If we depopulate most of the resonances by decreasing the width of the annulus, this effect becomes weaker. To clearly demonstrate the dynamics driving these growth modes, we also examine the evolution of a planetary embryo embedded in an annulus of collisionless planetesimals. In this case, we find that the resonances push planetesimals away from the embryo, decreasing the surface density of the bodies adjacent to the embryo. This effect only occurs when the annulus is wide enough and the mass resolution of the planetesimals is fine enough to populate the resonances. The bump we observe in the mass distribution resembles the 100 km power-law break seen in the size distribution of asteroid belt objects. Although the bump produced in our simulations occurs at a size larger than 100 km, we show that the bump location is sensitive to the initial planetesimal mass, which implies that this feature is potentially useful for constraining planetesimal formation models.
Proxima Centauri b provides an unprecedented opportunity to understand the evolution and nature of terrestrial planets orbiting M dwarfs. Although Proxima Cen b orbits within its star’s habitable zone, multiple plausible evolutionary paths could have generated different environments that may or may not be habitable. Here, we use 1-D coupled climate-photochemical models to generate self-consistent atmospheres for several evolutionary scenarios, including high-O2, high-CO2, and more Earth-like atmospheres, with both oxic and anoxic compositions. We show that these modeled environments can be habitable or uninhabitable at Proxima Cen b’s position in the habitable zone. We use radiative transfer models to generate synthetic spectra and thermal phase curves for these simulated environments, and use instrument models to explore our ability to discriminate between possible planetary states.
The obliquity of the Earth, which controls our seasons, varies by only ~2fdg5 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 semianalytic 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.
The υ Andromedae system is the first exoplanetary system to have the relative inclination of two planets’ orbital planes directly measured, and therefore offers our first window into the three-dimensional configurations of planetary systems. We present, for the first time, full three-dimensional, dynamically stable configurations for the three planets of the system consistent with all observational constraints. While the outer two planets, c and d, are inclined by ~30°, the inner planet’s orbital plane has not been detected.
We explore the impact of obliquity variations on planetary habitability in hypothetical systems with high mutual inclination. We show that large-amplitude, high-frequency obliquity oscillations on Earth-like exoplanets can suppress the ice-albedo feedback, increasing the outer edge of the habitable zone. We restricted our exploration to hypothetical systems consisting of a solar-mass star, an Earth-mass planet at 1 AU, and 1 or 2 larger planets. We verified that these systems are stable for 108 years with N-body simulations and calculated the obliquity variations induced by the orbital evolution of the Earth-mass planet and a torque from the host star. We ran a simplified energy balance model on the terrestrial planet to assess surface temperature and ice coverage on the planet’s surface, and we calculated differences in the outer edge of the habitable zone for planets with rapid obliquity variations. For each hypothetical system, we calculated the outer edge of habitability for two conditions: (1) the full evolution of the planetary spin and orbit and (2) the eccentricity and obliquity fixed at their average values. We recovered previous results that higher values of fixed obliquity and eccentricity expand the habitable zone, but we also found that obliquity oscillations further expand habitable orbits in all cases. Terrestrial planets near the outer edge of the habitable zone may be more likely to support life in systems that induce rapid obliquity oscillations as opposed to fixed-spin planets. Such planets may be the easiest to directly characterize with space-borne telescopes. Key Words: Exoplanets—Habitable zone—Energy balance models. Astrobiology 14, 277–291.
We study the influence of outer solar system architecture on the structural evolution of the Oort Cloud (OC) and the flux of Earth-crossing comets. In particular, we seek to quantify the role of the giant planets as “planetary protectors.” To do so, we have run simulations in each of four different planetary mass configurations to understand the significance of each of the giant planets. Because the outer planets modify the structure of the OC throughout its formation, we integrate each simulation over the full age of the solar system. Over this time, we follow the evolution of cometary orbits from their starting point in the protoplanetary disk to their injection into the OC to their possible re-entry into the inner planetary region. We find that the overall structure of the OC, including the location of boundaries and the relative number of comets in the inner and outer parts, does not change significantly between configurations; however, as planetary mass decreases, the trapping efficiency (TE) of comets into the OC and the flux of comets into the observable region increases. We determine that those comets that evolve onto Earth-crossing orbits come primarily from the inner OC but show no preference for initial protoplanetary disk location. We also find that systems that have at least a Saturn-mass object are effective at deflecting possible Earth-crossing comets but the difference in flux between systems with and without such a planet is less than an order of magnitude. We conclude by discussing the individual roles of the planets and the implications of incorporating more realistic planetary accretion and migration scenarios into simulations, particularly on existing discrepancies between low TE and the mass of the protoplanetary disk and on determining the structural boundaries of the OC.
Recent numerical simulations have demonstrated that the Sun’s dynamical history within the Milky Way may be much more complex than that suggested by its current low peculiar velocity (Sellwood, J.A., Binney, J.J. [2002]. Mon. Not. R. Astron. Soc. 336, 785–796; Roškar, R., Debattista, V.P., Quinn, T.R., Stinson, G.S., Wadsley, J. [2008]. Astrophys. J. 684, L79–L82). In particular, the Sun may have radially migrated through the galactic disk by up to 5–6 kpc (Roškar, R., Debattista, V.P., Quinn, T.R., Stinson, G.S., Wadsley, J. [2008]. Astrophys. J. 684, L79–L82). This has important ramifications for the structure of the Oort Cloud, as it means that the Solar System may have experienced tidal and stellar perturbations that were significantly different from its current local galactic environment. To characterize the effects of solar migration within the Milky Way, we use direct numerical simulations to model the formation of an Oort Cloud around stars that end up on solar-type orbits in a galactic-scale simulation of a Milky Way-like disk formation.
We investigate the evolution of grains composed of an ice shell surrounding an olivine core as they pass through a spiral shock in a protoplanetary disk. We use published three-dimensional radiation-hydrodynamics simulations of massive self-gravitating protoplanetary disks to extract the thermodynamics of spiral shocks in the region between 10 and 20 AU from the central star. As the density wave passes, it heats the grains, causing them to lose their ice shell and resulting in a lowering of the grain opacity. In addition, since grains of different sizes will have slightly different temperatures, there will be a migration of ice from hotter grains to cooler ones. The rate of migration depends on the temperature of the background gas, so a grain distribution that is effectively stable for low temperatures can undergo an irreversible change in opacity if the gas is temporarily heated to above ~150 K. We find that the opacity can drop more and at a significantly faster rate throughout the spiral shocks relative to the prediction of the standard dust grain model adopted in hydrodynamical calculations of protoplanetary disks. This would lead to faster gas cooling within spiral arms. We discuss the implications of our results on the susceptibility of disks to fragment into sub-stellar objects at distances of a few tens of astronomical units.
We evaluate the orbital evolution and several plausible origin scenarios for the mutually inclined orbits of υ And c and d. These two planets have orbital elements that oscillate with large amplitudes and lie close to the stability boundary. This configuration, and in particular the observed mutual inclination, demands an explanation. The planetary system may be influenced by a nearby low-mass star, υ And B, which could perturb the planetary orbits, but we find it cannot modify two coplanar orbits into the observed mutual inclination of 30°. However, it could incite ejections or collisions between planetary companions that subsequently raise the mutual inclination to >30°. Our simulated systems with large mutual inclinations tend to be further from the stability boundary than υ And, but we are able to produce similar systems. We conclude that scattering is a plausible mechanism to explain the observed orbits of υ And c and d, but we cannot determine whether the scattering was caused by instabilities among the planets themselves or by perturbations from υ And B. We also develop a procedure to quantitatively compare numerous properties of the observed system to our numerical models. Although we only implement this procedure to υ And, it may be applied to any exoplanetary system.
We present N-body simulations of planetary accretion beginning with 1 km radius planetesimals in orbit about a 1 M⊙ star at 0.4 AU. The initial disk of planetesimals contains too many bodies for any current N-body code to integrate; therefore, we model a sample patch of the disk. Although this greatly reduces the number of bodies, we still track in excess of 105 particles. We consider three initial velocity distributions and monitor the growth of the planetesimals. The masses of some particles increase by more than a factor of 100. Additionally, the escape speed of the largest particle grows considerably faster than the velocity dispersion of the particles, suggesting impending runaway growth, although no particle grows large enough to detach itself from the power law size-frequency distribution. These results are in general agreement with previous statistical and analytical results.
The water content and habitability of terrestrial planets are determined during their final assembly, from perhaps 100 1,000-km “planetary embryos” and a swarm of billions of 1–10-km “planetesimals.” During this process, we assume that water-rich material is accreted by terrestrial planets via impacts of water-rich bodies that originate in the outer asteroid region. We present analysis of water delivery and planetary habitability in five high-resolution simulations containing about 10 times more particles than in previous simulations. These simulations formed 15 terrestrial planets from 0.4 to 2.6 Earth masses, including five planets in the habitable zone. Every planet from each simulation accreted at least the Earth’s current water budget; most accreted several times that amount (assuming no impact depletion). Each planet accreted at least five water-rich embryos and planetesimals from the past 2.5 astronomical units; most accreted 10–20 water-rich bodies. We present a new model for water delivery to terrestrial planets in dynamically calm systems, with low-eccentricity or low-mass giant planets—such systems may be very common in the Galaxy. We suggest that water is accreted in comparable amounts from a few planetary embryos in a “hit or miss” way and from millions of planetesimals in a statistically robust process. Variations in water content are likely to be caused by fluctuations in the number of water-rich embryos accreted, as well as from systematic effects, such as planetary mass and location, and giant planet properties. Key Words: Planetary formation—Water delivery—Extrasolar planets—Cosmochemistry. Astrobiology 7(1), 66–84.
We present results of numerical simulations that examine the dynamical stability of known planetary systems, a star with two or more planets. First we vary the initial conditions of each system on the basis of observational data. We then determine regions of phase space that produce stable planetary configurations. For each system we perform 1000 ~ 106 yr integrations. We examine υ And, HD 83443, GJ 876, HD 82943, 47 UMa, HD 168443, and the solar system. We find that the resonant systems, two planets in a first-order mean motion resonance (HD 82943 and GJ 876) have very narrow zones of stability. The interacting systems, not in first-order resonance, but able to perturb each other (υ And, 47 UMa, and the solar system), have broad stable regions. The separated systems, two planets beyond 10 : 1 resonance (we examine only HD 83443 and HD 168443) are fully stable. We find that the best fits to the interacting and resonant systems place them very close to unstable regions. The boundary in phase space between stability and instability depends strongly on the eccentricities and (if applicable) the proximity of the system to perfect resonance. Furthermore, we also find that the longitudes of periastron circulate in chaotic systems but librate in regular systems. In addition to 106 yr integrations, we also examined stability on ~108 yr timescales. For each system we ran ~10 long-term simulations, and find that the Keplerian fits to these systems all contain configurations that are regular on this timescale.
We present results from 44 simulations of late stage planetary accretion, focusing on the delivery of volatiles (primarily water) to the terrestrial planets. Our simulations include both planetary “embryos” (defined as Moon to Mars sized protoplanets) and planetesimals, assuming that the embryos formed via oligarchic growth. We investigate volatile delivery as a function of Jupiter’s mass, position and eccentricity, the position of the snow line, and the density (in solids) of the solar nebula. In all simulations, we form 1–4 terrestrial planets inside 2 AU, which vary in mass and volatile content. In 44 simulations we have formed 43 planets between 0.8 and 1.5 AU, including 11 “habitable” planets between 0.9 and 1.1 AU. These planets range from dry worlds to “water worlds” with 100+oceans of water (1 ocean=1.5×1024 g), and vary in mass between 0.23M⊕ and 3.85M⊕. There is a good deal of stochastic noise in these simulations, but the most important parameter is the planetesimal mass we choose, which reflects the surface density in solids past the snow line. A high density in this region results in the formation of a smaller number of terrestrial planets with larger masses and higher water content, as compared with planets which form in systems with lower densities. We find that an eccentric Jupiter produces drier terrestrial planets with higher eccentricities than a circular one. In cases with Jupiter at 7 AU, we form what we call “super embryos,” 1–2M⊕ protoplanets which can serve as the accretion seeds for 2+M⊕ planets with large water contents.
Because of the high eccentricities (~0.3) of two of the possible planets about the star υ Andromeda, the stability of the system requires careful study. We present results of 1000 numerical simulations which explore the orbital parameter space as constrained by the observations. The orbital parameters of each planet are chosen from a Gaussian error distribution, and the resulting configuration is integrated for 1 Myr. We find that 84% of these integrations are stable. Configurations in which the eccentricity of the third planet is lesssim0.3 are always stable, but when the eccentricity is gsim0.45, the system is always unstable, typically producing a close encounter between the second and third planets. A similar exercise with the gas giants in our solar system sampled with the same error distribution was performed. Approximately 81% of these simulations were stable for 106 yr.