In the last post we looked at modelling the radial velocity of a star as the reflex motion due to an orbiting planet’s gravitational influence. This is the most prominent way used to model extrasolar planets to fit RV data. The Keplerian solution offers a lot in its simplicity, but it has a major shortcoming in that it is not particularly dynamic, and really only approximates any system with more than two bodies.

A planetary system may be thought of as an *n*-body problem, where you have *n* gravitating objects all influencing each other. This is different from the Keplerian models that take into account only the star-planet gravitational interaction. This will be a rather rudimentary look at simulating a gravitational *n*-body model.

First let’s look at vectors. In this case, we’ll consider a velocity vector defined simply with an *x* and a *y* component.

This simply allows us to get some sense of where our velocity is toward in a two-dimensional cartesian plane. We’re used to saying that an object has a velocity of, for example, 5 m s^{-1}, but this alone tells you little useful in the context of modelling a system – 5 m s^{-1} in *what direction*? By defining a vector with an *x* and *y* component, we can specify the direction that the velocity is in. For example, , then the velocity of the object is to the right, whereas implies moving down. If , then the object is moving both to the right and up, with the actual total velocity being determined by Pythagorean Theorem.

This is called the *magnitude* of the velocity vector.

Now each particle will have its own gravitational field. To determine the attraction to a particle from a given coordinate, we can use Newton’s Law of Universal Gravitation.

Where *m* is the mass of the particle we’re considering, let’s call it , and *r* is the distance between it and the other particle, let’s call it . Since we’re considering in this context only the gravitational attraction of on (not both particles on each other), we need not take into account the mass of both particles at this time.

Once you have determined the gravitational attraction of toward , you will need to determine which *direction* the gravitational attraction is toward. The angle between two points is simply given by

The *x* and *y* components of the gravitational attraction may be determined with and , respectively. These *x* and *y* components may be added to the *x* and *y* components of the velocity vector to determine the new .

This will have to be done for every particle. You will have to calculate the gravitational influence of all particles on each other particle. For *i* particles, the *new* velocity of any given particle *j*, , is given by

Where *r* is the separaton between the two particles. For most chosen values for , the system will have a systemic velocity (corresponding to the *γ* term in the Keplerian radial velocity equation). If you want to determine , you can simply sum up all of the initial velocity vectors and divide it by the total system mass *M*.

N-body simulations are calculation-intensive and are typically run on computers, where gravitational attractions and their modifications of velocity vectors may be calculated over and over for each new position. Particles are typically assigned a mass, a starting coordinate and a starting velocity vector, and a computer will calculate new positions as dictated by Newtonian physics for as long as the user wants.

Let’s look at a specific example. For simplicity, I set the gravitational constant to unity and use no units. I simulate three particles, a star and two planets, with positions of (0,0), (-30, 0), (60, 0), starting velocity vectors of (0.1, 0.1), (0, 0.5), (0, -0.22), and masses of 10, 0.1 and 0.1, respectively. Letting the simulation run for 3,000 time steps and plotting the y-axis component of the star’s velocity vector provides us with the stellar radial velocity.

What is going on? The RV behaviour is decidedly non-Keplerian – there are no Keplerian models that would be able to fit this data. If we look at a diagram of the system, it becomes clear what is happening here.

Our radial velocity plot is showing us a nightmare in planetary system architectures: planet-planet scattering. The outer planet (whose positions are shown in green) is clearly interacting very strongly with the inner planet (whose positions are marked in red).

Not all dynamically interacting multi-planet systems are so chaotic. Systems often have planets in close but stable orbits where they gravitationally tug on each other, causing the stellar radial velocity profile to deviate from a simple Keplerian model.

While Keplerian models are not capable of fitting to planetary systems with interacting planets, this Newtonian fit can. It is dynamical in that it simulates the influences of the planets on each other, however it can be particularly computer-intensive. Newtonian fits are often run of systems where Keplerian fits suffice simply to probe the dynamical behaviour of a system, such as its long-term stability or possible resonances between the orbits. Knowing whether to use Keplerian or Newtonian models for planets will depend on the planetary system architecture, but it is clear that Newtonian models are powerful tools for predicting the behaviour of planets in multi-planet systems.