Jupiter Ejected: Numerical Dynamics in Simulations of the Outer Solar System

By Junyu Chen

The outer solar system planets—Jupiter, Saturn, Uranus, Neptune, and Pluto (downgraded to dwarf planet in 2006)—travel in stable elliptical orbits around the Sun. However, numerical simulations of these trajectories—which use computational methods to solve the ordinary differential equations (ODEs) for the system—can cause error and instability. In the late 1980s, Edward Lorenz demonstrated that ODE solvers can introduce computational chaos [2]; since then, many studies have explored these effects in the context of the n-body problem [1, 3]. Here we examine the numerical dynamics of a six-body system that consists of the Sun, Jupiter, Saturn, Uranus, Neptune, and Pluto. We focus specifically on the ejection of Jupiter due to numerical effects.

One can derive the ODEs for the n-body problem by treating the bodies as point masses \(M_i\) that are governed by Newton’s universal law of gravitation in three-dimensional (3D) space: 

\[\vec{F}_{M_2M_1} = G \frac{{-M}_1M_2\ \vec{r}}{|r^3|},\]

where \(\vec{F}_{M_2M_1}\) is the force that \(M_2\) exerts on \(M_1\) and \(\vec{r}\) is the vector from \(M_1\) to \(M_2\). In our problem, this calculation produces a 36-dimensional ODE system; there are six second-order ODEs (one for every body), each of which involves a 3D \(\vec{r}\) vector.

Dozens of well-known algorithms can solve these types of systems, beginning with Euler’s method:

\[x_{t+1}=x_t+h\dot{x}_t.\]

Euler’s method is part of the Runge-Kutta family of methods, which use an estimate that is constructed via a power series to extrapolate the system state. The \(n^{\textrm{th}}\)-order Runge-Kutta truncates the series after the \(n^{\textrm{th}}\) term, rather than after the \(h\dot{x}_t\) term as in Euler’s method. This truncation is one source of numerical error — and numerical dynamics. Of course, a larger number of terms in the series begets a smaller truncation error, meaning that Euler’s method is the least accurate member of the Runge-Kutta family.

Figure 1. Outer solar system simulation by (1a) Euler’s method and (1b) fourth-order Runge-Kutta (RK4) ordinary differential equation solvers. The step size is set at 100 days for both solvers. Figure courtesy of Junyu Chen.

In our study, we utilize Euler’s method, fourth-order Runge-Kutta (RK4), and two symplectic solvers — common methods for the simulation of astronomical objects under gravitational interactions. The first symplectic solver is the leapfrog method, which preserves angular momentum; the second is WHFast (a Wisdom-Holman integrator), which preserves energy [3]. All of our solver implementations fix the time step at a constant value rather than adapting it on the fly, as time-step effects are of primary interest in our study.

Figure 2. Three outer solar system simulations by fourth-order Runge-Kutta (RK4) that demonstrate the ejection of Jupiter. 2a. Simulation with a solver step size of 120 days and a timespan of 34 kiloyears (kyr). 2b. Simulation with a solver step size of 150 days and a timespan of 15 kyr. 2c. Simulation with a solver step size of 170 days and a timespan of eight kyr. Figure courtesy of Junyu Chen.

The final matter involves choosing the value of time step \(h\). In general, a smaller time step increases the computational effort but yields a more accurate solution; indeed, Euler’s method limits to the first theorem of calculus as \(h\) goes to zero. However, the effects of floating-point arithmetic will come into play long before one reaches that point, which is another source of numerical dynamics. When simulating planetary orbits, it might seem reasonable to choose a step size that correlates with the orbital period. Unfortunately, this approach is problematic here because of the different periods of the planets in the study, which range from around 4,300 days for Jupiter to more than 90,000 days for Pluto. We thus expect the numerical error to affect Jupiter the most at any given step size (since all 36 ODEs are solved in parallel, we cannot use different time steps for different planets).

We begin by comparing 270-year simulations of the outer solar system that are produced by Euler’s method and RK4, both with a step size of 100 days (see Figure 1). The first part of Animation 1 gives the Euler’s method simulation.

Unsurprisingly, the Euler trajectories in Figure 1a deviate quite rapidly from the true elliptical orbit and spiral outwards because the solver uses the tangent to extrapolate the planets’ future positions. As we conjectured, this effect is strongest for Jupiter, which passes Pluto’s orbit after 170 years. In contrast, Neptune completes less than two revolutions during the experimental period and only deviates from its true orbit by a small amount.

Figure 3. Relationship between step size and time to eject Jupiter (in kiloyears). Ejection is defined as surpassing 50 astronomical units from the Sun (Pluto’s aphelion). Figure courtesy of Junyu Chen.

As we expected, the Runge-Kutta trajectories in Figure 1b are far better; they closely follow the true elliptical paths of the planets throughout the experiment. However, longer simulations reveal a different story; Jupiter first spirals in towards the Sun and then is ejected from the solar system (see Figure 2). Part 2 of Animation 1 portrays a simulation of this ejection (Jupiter also leaves the solar system in the Euler trajectory, but with a slow spiraling motion rather than an abrupt ejection).

The mechanism behind this effect is interesting. After spiraling within one astronomical unit of the Sun, Jupiter has only three to four time steps to complete an elliptical orbit with \(h = 100\) days, which is far too coarse for RK4 to extrapolate the trajectory with any reasonable accuracy. When the planet is in the vicinity of the Sun, it therefore takes a relatively gigantic step and picks up enough speed to eject itself from the solar system (the additional energy essentially comes from the large step size when compared to the distance between Jupiter and the Sun). 

The progression of images in Figure 2 shows that the time step of the simulation affects this behavior, though Jupiter is always the first planet to be ejected. Two interrelated reasons are responsible for this outcome. First, Jupiter has a tighter orbit and a shorter orbital period than the other planets, so the same \(h\) will have stronger effects. Second, the ejection mechanism depends on a planet’s proximity to the Sun. Indeed, missions like Voyager and New Horizons leveraged these dynamics to essentially use a close approach to a large body as a gravitational “slingshot.” Since Jupiter is the closest planet to the Sun in this six-body simulation, it feels these effects more strongly than the other four planets.

Such a combination of numerical and dynamical effects explains Jupiter’s proclivity for ejection. Following this line of reasoning, one would expect longer time steps to produce faster ejections (see Figure 3). Part 3 of Animation 1 similarly demonstrates this effect via three RK4 simulations with different time steps.

Figure 4. Three-dimensional trajectories of outer solar system simulation by (4a) WHFast and (4b) leapfrog. Figure courtesy of Junyu Chen.

The two symplectic solvers produce much more accurate trajectories (see Figure 4). This result makes logical sense, as these solvers are expressly designed to abide by the laws of classical mechanics and preserve energy (WHFast) or momentum (leapfrog). Part 4 of Animation 1 compares the animated trajectories of the four solvers—Euler’s method, RK4, WHFast, and leapfrog—using the same time step and total time duration.

In contrast, the RK4 trajectories do not preserve energy or angular momentum (see Figure 5). The resulting losses in energy and momentum are correlated to the inward spirals in Figures 1b and 2. Meanwhile, Figure 5 also proves that the conservation of these quantities in the symplectic solvers is imperfect.

WHFast conserves energy quite effectively. While the momentum of the trajectories does oscillate, that oscillation is fairly small; both the energy and momentum track the actual observations from NASA’s Horizons System quite well (see the red curves in Figure 5). As expected, the energy of the leapfrog trajectories is not conserved. Surprisingly though, neither is the momentum — though the oscillations are also small. These physical quantities do not match the observations as well as WHFast. Yet Jupiter is not ejected in either case regardless of the simulation length, since these algorithms are designed for conservative systems.

Figure 5. Energy and angular momentum of the simulated system versus time. In all simulations, \(h = 200\) days. We obtained the data labeled “Horizons” from NASA’s Horizons System with a time span of 1900 to 2022. Figure courtesy of Junyu Chen.

Numerical dynamics can have striking effects, as evident from the spiral trajectories of Euler’s method, the ejection of Jupiter by RK4, and the energy and angular momentum differences between the two symplectic solvers. Even though lowering the time step helps to mitigate the spiral effects, Euler’s method is certainly unacceptable for n-body simulations. RK4 is not suited for anything beyond a short-term simulation, as it eventually ejects Jupiter—an interesting but obviously nonphysical effect. Symplectic solvers like leapfrog and WHFast—which are designed for conservative systems—are the obvious choice and produce accurate simulations of the outer solar system.

References
[1] Hernandez, D.M., Hadden, S., & Makino, J. (2020). Are long-term N-body simulations reliable? Mon. Not. Roy. Astron. Soc., 493(2), 1913-1925.
[2] Lorenz, E.N. (1989). Computational chaos – a prelude to computational instability. Phys. D: Nonlin. Phenom., 35(3), 299-317.
[3] Wisdom, J., & Holman, M. (1991). Symplectic maps for the N-body problem. Astron. J., 102, 1528-1538.

Junyu Chen is an undergraduate student at the University of Colorado Boulder majoring in computer science, astrophysics, and mathematics. This work is from the final project in the Chaotic Dynamics course taught by Elizabeth Bradley.