Finding Planets with Star Wobbles



Introduction

The last 5-10 years has seen the number of known planets grow to almost two thousand (with almost four thousand waiting to be confirmed). This burst in discovery is due mostly to the Kepler telescope. Kepler stared at roughly 200,000 stars for five years looking for a dip in the light they emit. Planets transiting between us and their host star is one of the easiest and fruitful planet detection methods. However, the first exoplanet was discovered in 1995 by a completely different method; one that looks for the slight wobble of a star due to an orbiting planet. The radial velocity method has since been a consistent method for finding new planets, and confirming candidate planets detected by other methods.

In keeping with our theme of using Python to explore and analyze data, we will be working with real radial velocity data to find exoplanets with just a few built in Python functions.



The Radial Velocity Method

The radial velocity method works on the principle that a star + planet system will orbit about its common center of mass. The figure below depicts this.

alt text alt text
Credit


The speed at which the star moves and its displacement from the center are set by the planet's mass and how close it approaches the star. What radial velocity measures is the velocity of the star in the direction towards us or away from us. To make this measurement, astronomers take a spectrum of the star. If the star is moving towards us, the spectral lines will be blueshifted, whereas if the star is moving away form us, the lines will be redshifted. This shift is directly related to the speed of the star along the line of sight. From the depiction above, it's clear that if the system is orientated face-on (as in the right), there will be no observed motion along the line of sight. On the other hand, if we're viewing the system edge-on, then the radial velocity signal will be the largest. This tilt of the star-planet system on the sky is called inclination, and it is very difficult to measure. If we don't know the inclination of the system, then we can only get a lower bound on the planet's mass.

Radial velocity data might look something like this,

alt text



Here we're plotting the observed radial velocity in meters per second on the vertical axis versus time in days. We can see that the observations for this star occurred over several years. We can see by eye that the data follow roughly a sinusoidal pattern (especially for the last half). The goal, then, is going to be to fit this data to a general sine (or cosine) curve, $$V_r(t) = \text{amplitude} * \cos\left( \text{frequency}*t + \text{phase} \right) + \text{offset}$$ Once we have the best fit period, amplitude, phase, and offset of the data, we can relate them to the physical properties of the planet and its orbit.

From the analytic solution of the gravitational two-body problem we can write the velocity of the star along the line of sight as,

$$ V_r(t) = K \left( \cos( f(t) + \omega ) + e \cos\omega \right)$$

The time parameter here is hidden in the \(f(t)\) angle (also called the true anomaly). If we connect this equation back to the equation for the general cosine wave, we see that the \(K\) parameter is related to the amplitude, \(f(t)\) is related to time and the period, \(\omega\) is related to the phase, and \(e\) (the eccentricity of the orbit) is related to the offset. The semi-amplitude, \(K\), can be approximated as,

$$K = 28.4 \, \text{m} \text{s}^{-1} \, \frac{M_p \sin i}{M_J} \left( \frac{P}{1 \text{yr}} \right)^{-1/3} \left( \frac{M_\star}{M_\odot} \right)^{-2/3} \frac{1}{\sqrt{1-e^2}}$$

This says that maximum speed of the Sun if Jupiter where placed in the Earth's orbit would be about 28 m/s. Astronomers can detect these speeds down to ~ 1 m/s (that's walking speed!). We see that the strength of the radial velocity signal depends on the planet mass, the eccentricity, the orbital period, and the mass of the star. The true anomaly, \(f(t)\) is related to time via the following two equaions,

$$E - e \sin E = n (t-\tau)$$ $$\tan(f/2) = \sqrt{ \frac{1+e}{1-e} } \tan(E/2)$$

Here, \(E\) is the eccentric anomaly, \(n (t-\tau)\) is the mean anomaly, and \( n= 2 \pi / P \) is the mean motion of the planet. Given a time, \(t\), we first have to solve the first equation for the eccentric anomaly, \(E\). This has to be done numerically. In Python we can use the fsolve function to do this easily. Once we have the eccentric anomaly, we can plug it into the second equation to calculate the true anomaly, \(f\). The relationship between \(f\) and \(t\) is shown below for different eccentricities. The straight segments are for a circular orbit, while the most distorted curves are for an eccentricity of 0.9.

alt text



A nice technical summary (but not too technical) on the radial velocity method can be found here. Now that we've laid out our fitting procedure and connected our fitting parameters to our physical parameters, we can get into the coding.



The Python

Given the radial velocity data, Python makes it simple to define a fitting function and then fit the data to that function. We've explored this in one of the previous activities. There we were fitting linear and quadratic functions with one or two parameters. Here, however, we will be fitting five parameters, \(K, \, e, \, \omega, \, \tau, \, n \), to our radial velocity data. We will again be making use of Python's curve_fit function to find the best parameters that describe our data.

Download the two python files rv.py and rv_fit.py. Make sure to download them into the same folder. The first file contains helper functions that solve the orbital equations described above and functions to visualize the data and fit to the data. The second file contains the fitting function and the function to load the data. The second file is missing some code, denoted by FIX MEs.

Download the two data files for the stars HD 10442 and HD 5319 (courtesty of the authors). Make sure to save these to the same folder as the Python files. Once the data files are downloaded, fill in the missing code in rv_fit.py and find the exoplanets!



Results

If you coded everything up correctly, you should see plots that look something like these,

alt text alt text



Are these good fits to the data? The bottom plot shows the residuals of our fit. The closer the residuals are to zero, the better. Take a close look at the residuals for HD 5319. Do you notice anything interesting? It's not too hard to believe that residuals appear to follow another sine curve. We can go through the same procedure as before and fit the rediduals for a planet. In fact, HD 5319 does have two planet with similar periods.

Now that you have measured masses, eccentricities, and periods for these planets, you can compare your results to the accepted values. Note that we can get pretty close given our simple fitting procedure. In reality, astronomers use more sophisticated and custom methods in order to explore all of the uncertainties in the measurements and the models. Nevertheless, this method can prove to be a useful first attempt at discovering planets!



Extras

  • If you have time, dig a little deeper into the helper functions in rv.py. The plotting function demonstrates how to "glue" two plots together.
  • Try modifying the fit_data routine to accept x,y,y_err arguments so that you can fit whatever data you want. Apply this to HD 5319 to fit the residuals for a second planet.
  • Grab some more radial velocity data. One waty to do this you can start on the exoplanet catalogue site, find a planet, and then look up the reference paper. The authors usually provide their data in the text of the article. Here are a few more to try (these can be found here),
  • You can read up on how astronomers use the radial velocity method to find multiple planets here.