Numerical Integration Methods in Planetary Science

August 13-16, 2017

University of Toronto at Scarborough


David Hernandez (MIT): Time-symmetric integration in astrophysics

Through study of the modified differential equations, I show that time-symmetric integration in astrophysics has its limitations and can give undesirable errors.  Results from integrations of test problems will be presented.

Seppo Mikkola (Turku): A Comparison of some symplectic methods in very long computations

We compare the performance of four different symplectic integration methods (Mikkola and Palmer 2000, CeMDA 77, 305-317). The integrations were extended to the very long time of 10 billion years. The methods used were the basic Wisdom-Holman with leading  order symplectic corrector (WH), the Gaussian two-point formula, the alternating stepsize WH-method and the modified potential method.

We computed first the outer solar system motions by neglecting the Earth-like planets. The stepsize used was 100 days. The results show the expected energy errors that only fluctuate without any secular trends.

In the other set of computations we included all the planets from Mercury to Pluto. Due to Mercury's short period the stepsize used was 7 days. However, the energy errors behave somewhat differently from the first experiment, The reason is likely the 7 days stepsize witch may have been too long. However, the errors did not grow to large values and the orbital elements behaved in the same way for 30 to 40 million years, after which there were clear phase differences in the variations of the elements. In all but one (with the alternating stepsize method) the planet Mercury escaped after some billions of years. This is not surprising since other authors have arrived at similar conclusions earlier. One interesting result was that the e-folding times of the variational equations solutions seems to be only 3~Myr when Mercury is included."

Jihad Touma (AUB): Lie-Poisson Algorithms: From Planetary Spins to Galactic Nuclei


Hanno Rein (UTSC): JANUS: A bit-wise reversible integrator for N-body dynamics


Sarah Millholland (Yale): New Constraints on the Multi-Resonant Planetary System, Gliese 876

Gliese 876 is one of the most dynamically rich, well-understood exoplanetary systems. The nearby M dwarf hosts four known planets, the outer three of which are in a Laplace mean-motion resonance. A thorough characterization of the complex resonant perturbations exhibited by the orbiting planets, and the chaotic dynamics therein, is key to a complete picture of the system’s formation and evolutionary history. Here we present a reanalysis of the system using six years of new radial velocity data taken with the Levy spectrograph on the Automated Planet Finder (APF) telescope, the Keck HIgh Resolution Echelle Spectrometer (HIRES), and the Carnegie Planet Finder Spectrograph (PFS). These new datasets augment the decades-long collection of existing RV observations. We provide up-to-date estimates of the system parameters by employing a computationally efficient Wisdom-Holman N-body symplectic map, coupled with a Gaussian Process regression model to account for correlated stellar noise and a rapidly converging adaptation of Differential Evolution Markov Chain Monte Carlo. We present our ongoing investigations into the resonant dynamics of the system, including the coupling within the multi-resonant chain, the observability of second-order inclination resonances with current RV precision, the observable consequences of the dynamical chaos, and possible constraints on the system’s formational history.

Raluca Rufu (Weizmann Institute): Triton's Evolution with a Primordial Neptunian Satellite System

The numerical integration scheme in Phantom is leapfrog. This is a good choice for integration when accelerations are dependant upon position alone as you uphold leapfrog's symplectic nature. Is it a good choice when a wider range of physics is involved?

Sam Hadden (Northwestern/CfA): Disturbing Function Integrator

Developing a simplified analytic model is a frequent goal when investigating phenomena in planetary dynamics and usually involves isolating the relevant terms in a disturbing function expansion. I will describe a newly developed Python software package, celmech, that allows the user to easily construct and integrate dynamical equations by selecting individual disturbing function terms.  The package is designed to interface easily with the REBOUND code (Rein & Tamayo, 2015), allowing models to be initialized from, and directly compared with, REBOUND simulations. Celmech also provides transformations to and from orbital elements used in N-body integrations and canonical variables used in Hamiltonian models frequently used to study mean-motion resonances.

Ari Silburt (Toronto/Penn State): HERMES: A hybrid integrator for planetesimal migration

Here I present HERMES, a hybrid integrator designed for planetesimal migration. It builds off two existing rebound integrators, WHFAST and IAS15, using the former to integrate distant particles and the latter to integrate particles undergoing close encounters. In addition, an automatic Hill Switch Boundary has been implemented, ensuring that close encounters are not missed. In this talk I will discuss the strengths and weaknesses of the integrator, and compare its performance to existing integrators.

David Hernandez (MIT): Symplectic integration for the collisional gravitational N-body problem

I discuss efficient exactly symplectic integration for collisional N-body problems which tolerates close encounters.

Erika Nesvold (Carnegie DTM): Adding Collisions to Your N-Body Integrator: Challenges and Benefits

Debris disks, including the asteroid belt and Kuiper belt in our own Solar System, are often modeled as disks of massless test particles via N-body integration. This approach neglects interparticle interactions such as the fragmenting collisions between planetesimals that produce smaller, observable dust grains in debris disks. However, adding fragmenting collisions to an N-body integrator poses a computational challenge, as the number of bodies grows quickly. I will review the various approaches to this problem taken by debris disk researchers, and present results from my own research demonstrating the insights that collisional modeling can provide into the morphology of debris disks.

Martin Duncan (Queens): LIPAD


Simon Grimm (Bern): The GENGA code

GENGA is a GPU N-body code to simulate planetary systems. It uses a hybrid symplectic integrator to handle close encounters between bodies to achieve a good energy conservation. GENGA can simulate up to 30000 fully gravitational interactive bodies to perform planet formation simulations. Or it can use up to a million test particles to simulate the dynamics of small bodies. Alternatively, GENGA can run a large number of small planetary systems in parallel on the same GPU, which is useful to sample a large configuration space of exoplanetary systems or to compute TTVs.

Besides gravity, GENGA includes general relativity corrections, tidal forces, the Yarkovsky effect, Poynting-Robertson drag and a collision and fragmentation model for asteroids. The later can be used to simulate asteroid break-up events and meteorite delivery to terrestrial planets.

During the talk, some of the key problems of the parallelization for a GPU are presented and discussed.  

Junichiro Makino (Kobe): FDPS: Framework for developing large-scale parallel particle-based simulation code.

We have developed FDPS (Framework for developing particle simulators). The basic idea of FDPS is to separate the codes for parallelization and efficient interaction calculation from the codes for pairwise interaction calculation and time integration of particles, and thus to allow researchers try their ideas easily and efficiently on large-scale parallel computers. FDPS provides functions to perform domain decomposition, particle exchange, and interaction calculation, given the data structure of particles and functions to evaluate interactions between particles.

FDPS is implemented as a template library in C++ language, but users can write their programs (and data structures) either in C++ or Fortran. It is also possible to write interaction functions in Cuda or OpenCL, to make efficient use of GPGPUs and other accelerators.  We report the performance and scalability of several particle simulation codes, such as an pure N-body code for galactic dynamics, gravity+SPH, hybrid particle-particle particle-tree code for planetary formation. For many problems, scalability on large machines such as K computer and Sunway TaihuLight is pretty good, up to their full size.

Billy Quarles (U Oklahoma): Probing the early Solar System using GPUs

Studies of the early Solar System are often truncated to the leading order effects and scales due to pragmatic limitations on computing.  Recently, Grimm & Stadel (2014) have developed a GPU solution (GENGA) to this problem that is based upon the widely used code mercury (Chambers 1999).  This code allows for thousands of massive bodies, which makes it ideal for answering questions concerning the influence of an outer planetesimal disk on a primordial system of resonant giant planets (i.e., Nice Model).  We have used GENGA to efficiently examine the timing of a giant planet instability in the early Solar System and found that prior assumptions on interactions within the planetesimal disk appear to underestimate the amount of self-stirring and spreading of the disk.  We will discuss the numerical challenges associated with this problem and how we have used GENGA to investigate a possible Nice-like scenario in the early Solar System.

Pawel Artymowicz & Fergus Horrobin (UTSC): Massively parallel integration on CPU, GPU and MIC

We compare the capabilities of modern CPUs and coprocessors (GPU and MIC/Xeon Phi) to deliver massively parallel simulations. Up to N ~ several billion particles, as well as fluids on giga-cell grids using PPM method can be integrated in several days of computations for astrophysically interesting periods of time (1e3 - 1e5 orbits) on a small hybrid cluster, such as the one we have recently built at UTSC. We discuss in detail a Fortran90 implementation of the 4th order symplectic algorithm of N x 3B problem (a particle disk interacting with a binary system), with an accurate treatment of particles within Roche lobe of the secondary component. This project is part of our investigation of planet migration in disks.  We also discuss a massively parallel statistical study of long timescale (>1e9 of orbits) orbital stability of 1e4 planetary systems. Our general conclusions are that while all of the modern computational platforms are competitive in some applications (such as astro-CFD), MIC processors are particularly suited for accurate, massively parallel particle integrations.

David Tsang (U Maryland): Variational Integrators: Discretizing the Action

I will discuss the variational integrator formulation of symplectic integrators, showing how such integrators can be constructed by discretization of the action integral itself. Such integrators obey Noether's theorem, conserving constants of motion automatically when the (continuous) symmetries of the original action are preserved. I will demonstrate how this allows us to simply construct Wisdom-Holman type integrators for general (velocity dependent) potentials, as well as demonstrate the Wisdom-Holman-Farr method which allows effective adaptive integration while keeping a fixed formal step-size and preserving the symplectic properties.

Jack Wisdom (MIT): The Dynamical Systems Approach to Numerical Integration

I review the dynamical systems approach to numerical integration, and contrast it with the Lie series approach.

Jack Wisdom (MIT): Improving the Wisdom-Holman Integrator

Using the dynamical systems approach, I derive and test improvements in the Wisdom-Holman method.

Matija Cuk (SETI): Numerical implementation of tides for synchronous rotators

Tidal despinning and resulting orbital evolution are the dominant process driving the long-term evolution of planet-satellite systems. A very likely outcome of tidal despinning is the synchronous rotation, and almost all close-in moons in our solar system are currently rotating synchronously, always keeping one face towards their parent planet. However, synchronous rotation is not the end of tidal evolution, and tides within synchronously rotating moons are still important, usually acting to damp the moon's orbital eccentricity and (less commonly) inclination. Since the traditional theory of eccentricity damping is based on analytical description of satellite librations (e.g., Murray and Dermott, 1999), it is not suitable for implementation into a purely numerical integrator where forces should depend only on instantaneous positions, orientations and velocities. This is particularly important in cases when librations are large or synchronous lock is broken, which might have happened to the Moon during the Cassini State transition (Cuk et al, 21016). I will address ways of going around these problems in numerical codes, and how our choices of frequency-dependence may affect numerical integrations.

Terrence Tricco (CITA): Velocity dependant accelerations with leapfrog

Earlier this year we publicly released our smoothed particle hydrodynamics (SPH) code, Phantom. It has support for a variety of physics, ranging from compressible hydrodynamics and self-gravity, to magnetic fields, Navier-Stokes viscosity, Shakura-Sunyaev disc alpha viscosity, dust (one and two fluid approaches), ISM chemistry, numerous potentials, turbulence driving, and point masses. It has been primarily developed and used for stellar, planetary and galactic astrophysics as well as accretion discs.

David Tsang (U Maryland): Slimplectic Integrators: Variational Integrators for Non-Hamiltonian Systems

I will briefly present a recently developed non-Hamiltonian action principle which can capture nonconservative/dissipative behaviour in a variational formalism. Using this nonconservative action, I will apply the variational integrator method to non-Hamiltonian systems. Due to a generalization of Noether's theorem, the resulting "slimplectic" integrators possess the long term stability properties of symplectic integrators, but can be applied to systems where the "constants" of motion are no longer conserved, allowing these powerful numerical techniques to be applied to a much larger range of physical problems.

Dan Tamayo (UTSC/CITA): Operator Splitting Methods for Additional Effects in N-body Simulations

Symplectic integrators split the N-body problem into individually solvable Hamiltonian operators that are then applied in sequence. Their error behavior can be rigorously understood through the Baker-Campbell-Hausdorff (BCH) formula, which relates the full N-body mapping to the split mapping. However, it is not clear how to proceed when velocity-dependent and/or dissipative forces are added to the N-body problem. Such forces are ubiquitous (tides, 1PN general relativity corrections, radiation forces), and many studies have simply added such forces to a Wisdom-Holman scheme.

We will show that since the BCH formula is a statement about general operators (not necessarily symplectic), it still provides a powerful way to analyze the error behavior of N-body simulations including velocity-dependent and/or dissipative forces. We will illustrate and describe the numerical errors introduced when naively adding velocity-dependent forces into a Wisdom-Holman scheme, as well as simple steps to fix the mapping. We will also discuss how this operator splitting formalism can provide intuition and be used to estimate numerical errors that can guide choices for long-term integrations.

Finally, we will present REBOUNDx, a framework for adding additional effects to REBOUND simulations.

Jo Bovy (Toronto): Wendy: A 1-D Nbody Code

Gravity in 1D, planar geometry is remarkably simple, because the gravitational force is independent of distance. The 1D N-body problem can be solved efficiently to machine precision and can act as a testbed for understanding more complicated, realistic simulations. I will discuss the implementation of a new 1D N-body code, wendy, together with some extensions of the standard procedure to unequal masses and external acceleration fields. wendy is available at and can be run directly in your browser.

Julien Salmon (SwRI): HydroSyMBA: a 1d hydrocode coupled to an N-body integrator

The numerical modeling of circumplanetary disks and satellites is particularly challenging because each part of the system requires a very different approach. We have previously developed a 1-dimensional hydrocode that allows the study of the evolution of a disk with a radial grid, allowing one to resolve its time-evolving structure (Charnoz et al. 2010, Salmon et al. 2010). It includes formation of satellites at the disk's edge, but their orbital evolution is only approximated. The orbital evolution of discrete bodies is best modeled by N-body integrators. Of particular interest is the symplectic integrator SyMBA (Duncan et al. 1998) that can accurately treat close encounters and collisions between particles. Current computational capabilities now make it feasible to merge the 2 codes in an efficient manner. Our resulting integrator is constructed around SyMBA. At the beginning and end of the timestep (dt), we call the integration part of the hydrocode for dt/2. The hydrocode solves the equation of mass and angular momentum conservation for the disk, taking into account the viscous torque and the torques associated with the 1st-order Lindblad resonances from outer satellites (Meyer-Vernet and Sicardy 1987). The torque on the satellites is returned to SyMBA and used to compute an acceleration on each object. At the end of the global timestep, we check whether the disk is extending beyond the Roche limit, in which case we spawn a new moonlet that is added to the N-body code, making sure to conserve the angular momentum of the whole system.

Andrew Hesselbrock (Purdue): Planetary Rings that Crash onto Planets, and Crash Your Computer

RING-MOONS is a numerical code that simulates the evolution of viscous planetary rings, the accretion of satellites from these ring systems, and the exchange of angular momentum between satellites, the ring, and the primary body.  The ring is modeled as a 1-D Eulerian finite element model that generates Lagrangian satellites.  While the physics modeled by RING-MOONS is identical to HYDRORINGS, a similar code in the literature, there are some differences in the solution approach.  The evolution of the surface mass density of a ring is best described by a second-order, non-linear, partial differential equation.  The analytical solution to this equation is not known, and explicit solution methods to this equation readily go unstable.  Our model implements a stability criterion that is insufficient to set the timestep and gridspace of the simulation, but directs us to values that do yield stable results.

Damya Souami (UofT): The use of Expectation-Maximisation (EM) statistical algorithms for missing/incomplete data problems.

We have developed a general and rigorous statistical modelling framework for estimating properly parameters of many physical models that depend on noisy, uncertain or missing data.  I will be introducing here the general motivations and methods of the use of statistical EM (Expectation-Maximisation), SEM (Stochastic-Expectation-Maximisation) algorithms, and will share preliminary results of the application of such algorithms in the problem of « Quantisation of the Solar System ». (The associated paper shall be submitted this fall)

Scott Tremaine (IAS, Princeton): Integration Methods for Large N

In an N-body system of fixed total mass and radius, as N grows the

gravitational potential becomes smoother. However the integration

timestep needed to resolve close encounters gets shorter. Is there a

way to resolve this apparent paradox?

Cristobal Petrovich (CITA): Secular dynamics of few-body systems: vectorial formalism and applications

I will discuss the secular evolution of few-body systems using vectorial elements and some of the advantages these elements offer over the more traditional orbital elements. I will discuss the application to four-body systems and triple systems perturbed by external potentials.

Jihad Touma (AUB): Softened Gauss: Tool and Applications


Seppo Mikkola (Turku): Chaos in the Trappist-1 system

I computed the evolution and Lyapunov exponents of the Tappist-1 planetary system. In ten experiments the masses, inclinations and eccentricities were varied within the published error estimates. The system has a very short e-folding time in the variational equations, but it seems to be stable, provided the eccentricities are set to very small values. The masses and inclinations seem to have little effect in the stability. The solution vectors of the variational equations are, most of the time, either parallel to the velocities or in the opposite direction. This suggests that the chaos is mainly in the phases of the planetary motions. One may at least speculate that this makes the orbits stable despite the closeness of resonances. Little disturbances in the periods makes it easily less chaotic, but such systems are often unstable because of eccentricity increase. To obtain the stable orbits it was necessary to adjust the periods to the observed ones. The formation of this stable system is a most interesting problem for celestial mechanics.

David Minton (Purdue): The challenges of numerical modeling of heavily cratered terrains

I will discuss some of the major challenges to modeling the evolution of heavily cratered terrains using a Monte Carlo crater code. I will show that observed equilibrium crater counts can only be reproduced using an anomalous, rather than classical, diffusion model for diffusive downslope topographic erosion. I will discuss how the challenging dynamic range problem of the cratering problem is approached with both sub-pixel and super-domain techniques. I will also discuss the use machine learning techniques to simulate human crater counting.

Gongjie Li (Harvard): Cross-sections for planetary systems interacting with passing stars and binaries

Most planetary systems are formed within stellar clusters, and these environments can shape their properties. This talk considers scattering encounters between solar systems and passing cluster members, and calculates the corresponding interaction cross-sections. The target solar systems are generally assumed to have four giant planets, with a variety of starting states, including circular orbits with the semimajor axes of our planets, a more compact configuration, an ultracompact state with multiple mean motion resonances, and systems with massive planets. We then consider the effects of varying the cluster velocity dispersion, the relative importance of binaries versus single stars, different stellar host masses, and finite starting eccentricities of the planetary orbits. For each state of the initial system, we perform an ensemble of numerical scattering experiments and determine the cross-sections for eccentricity increase, inclination angle increase, planet ejection, and capture. This paper reports results from over 2 million individual scattering simulations. Using supporting analytic considerations, and fitting functions to the numerical results, we find a universal formula that gives the cross-sections as a function of stellar host mass, cluster velocity dispersion, starting planetary orbital radius, and final eccentricity. The resulting cross-sections can be used in a wide variety of applications. As one example, we revisit constraints on the birth aggregate of our Solar system due to dynamical scattering and find N ≲ 10^4 (consistent with previous estimates).

Sarah Millholland (Yale): Constraints on Planet Nine’s Orbit and Sky Position in a Mean-Motion Resonant Framework

A number of authors have proposed that the statistically significant orbital alignment of the most distant Kuiper Belt Objects (KBOs) is evidence of an as-yet undetected planet in the outer solar system, referred to colloquially as “Planet Nine”. Dynamical simulations by Batygin & Brown (2016) have provided constraints on the range of the planet's possible orbits and sky locations. We have recently extended these investigations by exploring the suggestion of Malhtora et al. (2016) that Planet Nine is in small integer ratio mean-motion resonances (MMRs) with several of the most distant KBOs (Millholland & Laughlin 2017). We show that the observed KBO semi-major axes present a set of commensurabilities with an unseen planet at ~654 AU that has a greater than 98% chance of stemming from a sequence of MMRs rather than from a random distribution. We describe the constraints resulting from a Monte-Carlo optimization scheme that used billion-year dynamical integrations of the outer solar system to pinpoint the orbital properties of a planet that is capable of maintaining the KBOs' apsidal alignment. We also demonstrate a connection between the KBOs’ ability to maintain alignment and orbital stability and their participation in resonance. Finally, we discuss constraints on the planet’s current sky position.

Millholland, S. & Laughlin, G. 2017, AJ, 153, 91