Forcing Planets to Evolve: Interactions Between Uranus and Neptune at Late Stages of Dynamical Evolution

TL;DR

The work tackles the challenge of exploring chaotic late-stage outer Solar System evolution by large N-body simulations. It introduces a method to artificially evolve planetary orbital elements independently through user-defined velocities and accelerations, implemented in REBOUND and Mercury6.2, to emulate the dynamical influence of a massive planetesimal disk without simulating it. The authors validate the approach with one- and two-planet tests, showing independent control is robust for a single planet but couples through secular interactions in multi-planet systems, and they reproduce the Uranus-Neptune migration and eccentricity-damping episodes described in Tsiganis et al. (2005) by using dynamical-friction-based timescales and exponential damping. They discuss limitations (notably nonconservation of angular momentum and resonant effects) and highlight the method's potential to systematically explore orbital histories and their impact on trans-Neptunian populations.

Abstract

In early Solar System numerical simulations, where chaos is a primary driver, it is difficult to explore parameter space in a systematic way. In such simulations, stable configurations are hard to come by, and often require special fine-tuning. In addition, it is infeasible to run suites of well-resolved, realistic simulations with a disk of massive particles to drive planetary evolution where enough particles remain to represent the transneptunian populations to robustly statistically compare with observations. To complement state of the art full N-body simulations, we develop a method to artificially control each planet's orbital elements independently from each other, which when carefully applied, can be used to test a wider suite of models. We modify two widely used publicly available N-body integrators: (1) the C code, \texttt{REBOUND} and (2) the FORTRAN code, \texttt{Mercury6.2}. We show how the application of specific fictitious forces within numerical integrators can be used to tightly control planetary evolution to more easily explore migration and orbital excitation and damping. This tool allows us to replicate the impact a massive planetesimal disk would have on the planets, without actually including the massive planetesimals, thus decreasing the chaos and simulation runtime. We demonstrate the utility of this tool by applying it to the coupled orbital evolution of Uranus and Neptune, and show that Neptune's eccentricity damping and radial outward migration have the appropriate affect on Uranus' eccentricity.

Figures (10)

• Figure 1: Fifty million year evolution of the semimajor axis ($a$ (au); top left), eccentricity ($e$; top middle), inclination ($i$ (deg); top right), argument of pericenter ($\omega$ (deg); bottom left), longitude of ascending node ($\Omega$ (deg); bottom middle), and true anomaly ($f$ (deg); bottom right) of a Jupiter-mass planet orbiting a solar-mass star. The planet is initialized with the orbital elements shown in Table \ref{['all-sim-params']}. The orbital elements $a$, $e$, $i$, $\omega$, and $\Omega$ are artificially time evolved by the arbitrary functions (dashed line) in Equation (\ref{['aeifuncs']}), denoted by the "Eq" column in Table \ref{['all-sim-params']}. Each orbital element in fact evolves independently from each other and follows the functional form prescribed in the code.
• Figure 2: One hundred fifty million year evolution of the semimajor axis ($a$), eccentricity ($e$), inclination ($i$) (orange) of a Jupiter-mass (top row) and Neptune-mass (bottom row) planet orbiting a solar-mass star. The orbital elements are initialized and artificially evolved by an exponential (dark maroon dashed line) with parameters described in Table \ref{['all-sim-params']}. Similar to the one planet example in Figure \ref{['fig:1planetex']}, the additional user velocity and acceleration prescription evolves the orbital elements analogous to a desired function. However, the evolution of the eccentricity and inclination of each planet have a short timescale oscillation on top due to the secular effects between the planets. In contrast to the one planet example, where we can evolve the element to a final value, in this case, the end result is near, but not exactly the final value. This is due to the two planets causing additional oscillations of each other's orbital elements.
• Figure 3: Fifty million year evolution of Neptune's inclination for three different simulations (solid lines) initialized as described by in Table \ref{['all-sim-params']}. We artificially evolve Neptune's inclination (dashed lines) as a sinusoid, varying $\Delta i$ and $\tau_i$ relative to the secular amplitude, $\mathcal{A}_{\rm sec}$ and oscillation timescale, $\tau_{\rm sec}$ with three examples: (1) $\Delta i > \mathcal{A}_{\rm sec}$ and $\tau_i > \tau_{\rm sec}$ (top panel), (2) $\Delta i > \mathcal{A}_{\rm sec}$ and $\tau_i < \tau_{\rm sec}$ (middle panel), and (3) $\Delta i < \mathcal{A}_{\rm sec}$ and $\tau_i > \tau_{\rm sec}$ (bottom panel).
• Figure 4: Two, fifty million year evolutions of Jupiter's $\omega$ under the gravitational influence of a solar-mass star and Neptune-mass planet. Both simulations are initialized with the parameters described in Table \ref{['all-sim-params']}. The first simulation does not include any additional velocities or accelerations (orange). In the second simulation, we add a linear perturbation (maroon) over time to Jupiter's $\omega$, making Jupiter precess faster (blue).
• Figure 5: Eccentricity for Uranus (top panel) and Neptune (bottom panel) calculated by using the semimajor axis, pericenter, and apocenter curves from Figure 1 of Tsiganis:2005. The best-fit exponential damping timescale at the early stage (before 20 Myr) is 1.5 Myr for both Uranus and Neptune. At late times (after 30 Myr), the best-fit exponential timescales are 20 Myr and 6 Myr for Uranus and Neptune, respectively (see Section \ref{['sec:bestfit_eccentricity']} for the best fit final and initial eccentricities).