Research

Snapshot of REBOUND simulation of Saturn's Rings.

REBOUND code

Many N-body codes have been developed over the years. But REBOUND is different. REBOUND is built around a few core principles:

• Open Source. REBOUND is open source and will always be open source. Furthermore, all development is done in the open too. There are no secret development branches only available to colloaborators. I don't think secrets and science work well together.
• Easy to use. The python interface to REBOUND allows you to access all the features of REBOUND with minimal effort. Setting up a simulation often takes only a few lines of code. Built-in interactive visualization modules let you monitor the simulations in real-time.
• Accuracy. The integrators that come with RBEOUND, WHFast and IAS15 are among the world's most accurate integrators for long-term N-body simulations.
• Speed. Many of the core routines in REBOUND have been heavily optimized.

The REBOUND project started in summer 2011, when Shang-Fei Liu and I started to write REBOUND during the ISIMA summer school. It is a highly modular code that can handle a variety of problems such as long term orbit integrations, planetary rings and even debris disks. We use a set of highly efficient numerical algorithms and integrators. For example, the Symplectic Epicyclic Integrator (SEI) that I developed together with Scott Tremaine, is ideally suited for studying ring dynamics. The WHFast integrator, on the other hand is well suited for studies of the Solar System. The self-gravity between particles can be calculated using a distributed tree, an FFT method or even a GRAPE. The collision detection makes use of a tree code or an efficient plane-sweep algorithm.

In recent years REBOUND has become a community project with many collaboratos. Most notably Dan Tamayo and Dave Spiegel. Many of my students have also contributed to REBOUND. The code is still being actively developed. Check out its github page to see what is new. Although REBOUND is a very young code, it has already lead to numerous publications.

Planet 9 and Kuiper Belt objects (Caltech/R. Hurt/WorldWide Telescope)

Numerical integrators

Numerical integrators are algorithms that can be used to solve equations of motions. These are differential equations which can be hard to solve accurately and efficiently. I've developed several numerical integrators that can be used for systems which interact gravitationally such as planetary systems, moons or particles in Saturn's rings. All of these integrators are part of the REBOUND package.

IAS15

IAS15 is a 15th-order integrator specifically tuned for high accuracy in gravitational dynamics. The integrator is based on the well known Gauss-Radau quadrature first used for Solar System integrations by Everhardt. Dave Spiegel and I provide an implementation that has been significantly optimized and can handle conservative as well as non-conservative forces. It comes with an automatic step-size control to choose an optimal timestep. The algorithm can handle close encounters and high-eccentricity orbits. The systematic errors are kept well below machine precision and long-term orbit integrations over a billion orbits show that IAS15 is optimal in the sense that it follows Brouwer's law, i.e. the energy error behaves like a random walk. Our tests show that IAS15 is often superior to a mixed-variable symplectic integrator (MVS) and other popular integrators, including high-order ones, if high accuracy is the goal. In fact, IAS15 often preserves the symplecticity of Hamiltonian systems better than the commonly-used nominally symplectic integrators.

WHFast

WHFast is a fast and accurate implementation of a Wisdom-Holman symplectic integrator for long-term orbit integrations of planetary systems. WHFast is significantly faster and conserves energy better than all other Wisdom-Holman integrators that we tested. Dan Tamayo and I have achieve this by significantly improving the Kepler-solver and ensuring numerical stability of coordinate transformations to and from Jacobi coordinates. These refinements allow us to remove the linear secular trend in the energy error that is present in other implementations. For small enough timesteps we achieve Brouwer's law, i.e. the energy error is dominated by an unbiased random walk due to floating-point round-off errors. We implement symplectic correctors up to order eleven that significantly reduce the energy error. We also implement a symplectic tangent map for the variational equations. This allows us to efficiently calculate two widely used chaos indicators the Lyapunov characteristic number (LCN) and the Mean Exponential Growth factor of Nearby Orbits (MEGNO).

WHFast is machine independent. Surprisingly most other numerical integrators are not machine independent because they rely on mathematical function such as sin and cos which vary slighlty depending on the compiler, operating system or library version. This means that rerunning the exact same simulation on a different computer might result in completely different results. We achieve machine independence by not using any mathematical functions other than additions, multiplications and square roots. We implement all other functions ourselves as serious expansions.

The original version of WHFast uses Jacobi coordinates. We recently implemented WHFastHelio, using the same efficient Kepler solver as WHFast but using democratic heliocentric coordinates that are better suited for systems in which close encounters might occur.

SEI

The shearing sheet is a model dynamical system that is used to study the small-scale dynamics of astrophysical disks. Numerical simulations of particle trajectories in the shearing sheet usually employ the leapfrog integrator, but this integrator performs poorly because of velocity-dependent (Coriolis) forces. SEI is a new integrators for this purpose; it is symplectic, time-reversible and second-order accurate, and can easily be generalized to higher orders. Moreover, the integrator is exact when there are no small-scale forces such as mutual gravitational forces between disk particles. In numerical experiments Scott Tremaine and I have shown that this integrator has errors that are often several orders of magnitude smaller than competing methods.

Artist's conception of Kepler-11 (NASA/Tim Pyle)

Formation of planetary systems

The Kepler spacecraft has discovered thousands of planet candidates. Within that sample, there are many multi-planetary systems. These are especially interesting because they can constrain planet formation scenarios. For example, one can see statistically significant peaks in the period ratio distirbution near integer ratios such as 2:1. These locations are mean motion resonances. For example, if a system is in an exact 2:1 resonance, then one planet orbits the stars exactly twice as often as the other. The same features have been found in systems discovered by radial velocity.

The origin of these peaks is linked to the planets' formation history and most liekly planetary migration. As planets form, they interact with the proto-planetary disk and exchange angular momentum. The planets then start to move around and can capture into resonances with other planets. However, this is only part of the story. Because most of the planets in the Kepler sample are very small (sub-Jupiter mass), they are especially sensitive to random fluctuations in the proto-planetary disk. These fluctuations are, among other effects, created by the magneto rotational instability (MRI). They make planets undergo a random walk. Thus, there is a competition between convergent migration pushing planets into resonance and the random walk moving them away from the exact resonance.

I simulated all discovered Kepler systems, including migration forces. By taking canonical values for both the smooth and stochastic component, the final period ratio distribution (dashed curve in the figure) looks surprisingly similar to the observed distribution. The results show that planetary migration and not tides are responsible for shaping the architecture of most planetary systems. In the future, we will be able to constrain migration even more, by comparing planets of different location and size.

Mean motion resonances (MMR) are common among extra-solar planets. The best example is the Laplace 4:2:1 resonance in the GJ876 system. These systems are generally taken as strong evidence that planet migration exists.

The system HD200964 is also in a mean motion resonance, more precisely a 4:3 MMR. And this leads to a problem. Because HD200964 consists of planets that are very massive, the systems should have captured into the 2:1 resonance during migration. In a study with Matthew Payne, Dimitri Veras and Eric Ford, we found that it is practically impossible to form this system by any means!

This is a somewhat unsatisfying result because this planetary system should not exist! And it's not the only one. There are several systems that look exactly the same.

A (very big) moonlet inside Saturn's rings (NASA/JPL/SSI)

Stochastic migration of moonlets in Saturn's rings

As the name suggests, moonlets are small moons. In fact they are so small that they are embeped in Saturn's rings, not being massive enough to open a complete gap. Many of these objects have been observed by the Cassini Spacecraft. The following picture shows several of them.

Interestingly, these moonlets don't have fixed Keplerian orbits. They move around through interactions with the rings. Proto-planets experience a very similar effect soon after they form. They also interact with a disk (the proto-planetary disk) and begin to move. That is called planetary migration.

During my PhD I studied both planetary migration and the migration of moonlets. In both cases, stochastic migration will be important, especially for small mass bodies. John Papaloizou and I have shown that the stochasticty in Saturn's are generated by the finite size of ring particles and the interaction with gravitational wakes. In a collaboration with Margaret Pan, we were able to show that the observed motion is indeed compatible with such a random walk.