This just sets up the plotting to be "inline" rather than opening separate windows for each plot, and boosts the resolution (I have a "retina" display on my mac, so you need pretty high resolution to get good-looking plots).

You will need the rebound and seaborn packages to run (but not view) this notebook. You can create a fresh environment called trappist and install them on a machine with Anaconda in the following way:

conda create -n trappist
pip install rebound
pip install seaborn

Then issuing the command source activate trappist (mac or linux) or activate trappist (windows) will activate this environment. From here you can run python code that uses these modules, or you can start up jupyter notebooks (like this one!):

source activate trappist
jupyter notebook TRAPPIST1.ipynb

Note that rebound will not install on Windows; if you do not have access to a mac/linux machine (Y2 labs have linux machines, for example, as do some of the other computing labs across campus) you will only be able to look at this code, not run it. Don't worry---there are many ways to get 10 marks on this assignment, so you will still be able to complete the work.

If you need a tutorial on Jupyter notebooks, you can find one here: http://jupyter.org

These are a few of my favourite seaborn plotting settings:

A couple of utility functions for dealing with angles that wrap:

Set up the TRAPPIST-1 parameters. See the paper on Canvas or here. epoch is a time (in days) that sits in the middle of the sequence of transits reported there.

A word about rebound: it stores the planetary system it is currently simulating in a "simulation" object. Between integrations, you can add and remove planets from the sim, change their properties, etc. At the beginning of the integration, you set up your sim with all the planets you want to include in it. You can see examples of this in both the C and Python code at the rebound git repository. There are even iPython notebooks with examples, too! There is a nice feature in the interface as well, which is called a "simulation archive." It stores the states of a simulation at various points through the integration and allows for replay.

Here we write some code to set up the TRAPPIST system. We take the convention that the $x$-axis points toward the Earth, so when planets cross the $x$ axis in the $x$-$y$ plane, they transit. We don't know the exact orientation of the orbits in the system (i.e. the angles $\omega$, $\Omega$ or $i$), so we choose these randomly, and then work out the orbital phase that is necessary to ensure that the object transits at the given time for the orbit we have chosen.

Because we don't know the exact masess, periods, etc, (well, periods are actually really precisely determined), we draw them randomly within the reported uncertainties. Therefore, every time we run the code below we are drawing a Monte-Carlo sample of a possible TRAPPIST1 system that is consistent with the observations.

• For 3 Marks: Comment the code below, explaining what each statement is doing and how we ensure that the planets transit at the appropriate time.

A note on units: we choose $G = 4\pi^2$, which means we can measure masses in $M_\odot$, time in years, and lengths in AU.

• For 1 Mark: Explain why we can do this. Kepler's third law will be helpful here.

The bits at the beginning where we are modifying parameters of the simulation are just set up to optimize the performance of the simulation. It will use a type of integrator called "Wisdom-Holman" (the same Wisdom as explained the Kirkwood gaps we talked about earlier in class). You can read about W-H integrators here; essentially, they exploit the fact that the planets move along Keplerian orbits about the host star with only small perturbations to make the integration much more efficient. Rebound also uses a "corrector" to improve the accuracy of the simulations; those with interest in Hamiltonian dynamics should read about correctors here.

Rebound has the nice feature that it can exit (stop the simulation) if any body gets sufficiently far from the central star; here we set a boundary at $100 \, \mathrm{AU}$ to detect if any bodies have been ejected from the simulatino.

For those without access to rebound or who just prefer C or C++, here is a suite of initial conditions (1k) in all. The format is

• One line with the 8 masses of the star and planets (in orbital order)
• Eight lines giving "x y z vx vy vz" for each object (in the COM frame)

You can use these initial conditions to fulfill parts of the project using your own code if you like.

Note that you can access all the particles in a simulation via sim.particles; each particle has its own orbital elements (p.e, p.omega, p.pomega, p.Omega, p.a, p.l (mean longitude), p.i, p.P (period), etc) in addition to the six-dimensional $\vec{x}$ and $\vec{v}$ that we use below

• For 5 marks Write your own code (i.e. not rebound) that loads one of these initial conditions and integrates the system for 1000 years. Remember: this is just Newtonian gravity: $$ \vec{a}_i = \sum_{j \neq i} \frac{G M_j \left( \vec{r}_j - \vec{r}_i \right)}{\left| \vec{r}_j - \vec{r}_i \right|^3} $$ Plot the trajectories of the system in the $x$-$y$ plane for the first and last 10 years. Is your system unstable over 1000 years?

Let's do the above 5-mark part with rebound. First, we set up a simulation, and initialise an archive stored on disk in "1000yr.bin". We will have the archive store its state every year (or so).

Now we integrate the system to 1000 years. The state of the system will be stored 1 year or so; we don't force the integrator to stop evenly on 1-year boundaries, but rather let it keep its step-size and just record the state at the closest step to 1-year boundaries. (Note that we don't have to make the simulation archive if we don't want to save the intermediate states of the system---we can just call sim.integrate without setting up the archive.)

Now we can load the archive and make some plots:

Let's check the energy error in the simulation (note that energy should be exactly conserved, so any difference between the initial energy and subsequent energies is down to integration error).

• For 2 marks if you have written your own code to do the simulations, make a similar plot of the energy error over your 1000 yr simulation.

Let's plot the state of the system at the beginning and end of our simulation:

Now let's plot the eccentricity of each planet as a function of time over our simulation.

• For 2 marks If you have written your own code to simulate the system, make this plot. Note that you can calculate the eccentricity using the following two formulas for the energy and angular momentum of an orbit: $$ E_\mathrm{orb} = -\frac{G M m}{2 a} $$ and $$ \frac{L}{L_\mathrm{circ}} = \sqrt{1-e^2} $$ with $$ L_\mathrm{circ} = m \sqrt{G M a} $$ you can calculate the energy and angular momentum of the orbit at any instant by computing the separation between a planet and the star and the relative velocity between these two, and then $$ E \simeq \frac{1}{2} m v_\mathrm{rel}^2 - \frac{G M m}{r} $$ and $$ L \simeq m \left| \vec{v}_\mathrm{rel} \times \vec{r} \right| $$ since $m \ll M$.

It looks like, in this particular simulation, there are a number of quasi-periodic oscillations in the eccentricity of some of the planets. This is a standard feature of planets in disks with small eccentricity; essentially we are seeing the normal modes of the system. If you want to know more, look up "Laplace-Lagrange theory," or see Murray & Dermott, Chapter 7.

What about resonances? The discovery paper states that the inner two planets are near a 5:8 commensurability; is there a corresponding resonant angle: $$ \theta_\mathrm{res} = 5 \lambda_1 - 8 \lambda_2 + 3 \begin{cases}\varpi_1 & \mathrm{inner} \\ \varpi_2 & \mathrm{outer}\end{cases} $$ that oscillates?

Apparently not.

From here, there are a number of ways to proceed in the assignment; pick whichever combination of tasks you would prefer. I will mark everything you do, but the maximum score is 10 marks; only do significantly more than 10 marks worth of work if you are excited to explore these topics (good for you!).

• For 2 marks Find an initial condition that goes unstable (a planet escapes) within 10k years. Plot the initial configuration; can you guess that this would be unstable based on the planetary spacings?
• For 2 marks Estimate the fraction of draws that result in unstable systems on 10k year timescales.
• For 2 marks Plot the 7th planet mass and initial period for a ensemble of simulations integratied over 10k yr. Indicate on your plot which values resulted in stable systems and which gave unstable systems; do you see a pattern?
• For 3 marks A recent paper, arXiv:1703.04166 claims a $P = 18.764 \, \mathrm{d}$ orbit for the outer planet in the system based on further observations of the TRAPPIST1 system with during campaign 12 of the K2 mission. Modify the initial-condition-generating code to produce this period for the outer planet but otherwise without changing the distributions of the other parameters, and investigate the stability of systems where the outer planet has this period. Are they generally stable? Are there any you can find that are unstable? Can you estimate or bound the probability of instability?
• For 3 marks Check that the initial conditions we generated have the correct transit times. Look up the times of each planet's transit in the discovery paper, and integrate the system to a transit (via sim.integrate((t_transit - epoch)*days) or using your own code). Verify that the appropriate planet is close to the $x$-axis at the transit time.
• For 4 marks Look for resonances in some integrations of the system (see the discovery paper for a list of possible commensurabilites). You will have the best luck looking for lower order resonances (the inner two planets could be in a third order resonance---$5$:$8$---whose strength will scale as $e^3$, which is pretty weak). You don't have to find any, but show at least ten different searches if you don't find one.
• For 5 marks For an ensemble of simulations over 10k yr (at least ~100), record the planet masses from the stable simulations. Plot a histogram of the masses of the 7th planet; compare to the Gaussian distribution used to draw simulations above. Does the demand for stability constrain the mass of the seventh planet beyond what is already known from the transit observations? Plot such a histogram for all the planets; do you obtain better constraints on the masses of the inner planets, too?
• For 5 marks One of the signs of near-resonant or resonant orbits are transit timing variations (TTVs). Given that transits occur when a planet crosses the $x$-axis (here assume for simplicity 2D orbits, so this is equivalent to $y = 0$), write code that integrates the system for 10s of years and extracts the times of all transits of one (or several, if you like) of the planets. Plot $t_\mathrm{transit} \mod P$ for all the transits; if the orbit were strictly periodic, you would see a consistent orbital phase at transit (i.e. a horizontal line). Do you? If not, how significant are the transit time variations?
• For 5 marks A recent paper, arXiv:1703.04166 claims that the all planets in the TRAPPIST1 system are in three-body resonances based on further observations of the TRAPPIST1 system with during campaign 12 of the K2 mission. See Table 1 from arXiv:1703.04166 for the parameters of the resonances. We have not discussed three-body resonances in class, but the resonant angle is $$ \Theta = p \lambda_1 - (p+q) \lambda_2 + q \lambda_3. $$ Note that this sum of angles satisfies rotational invariance (i.e. the sum of the coefficients vanishes) by including the mean motions of three planets (hence "three-body resonance") rather than a second angle (e.g. $\varpi$) from one of two planets' orbits. Plot this resonant angle for several triples of planets in several simulations. Can you find a simulation in which the angle librates (oscillates) rather than circulating? Does the central value of the libration match the constraints reported in Table 1 of arXiv:1703.04166?
• For 5 marks Download and run some of the simulations reported in Tamayo, Rein, Petrovich & Murray (2017) and verify the resonances they describe. For an additional 10 marks, reproduce their musical composition linked to from Tamayo's website.