Console Tutorial #2
greg posted in Uncategorized on November 21st, 2005
Updated: Jan. 12th, 2007.
The star Upsilon Andromedae lies roughly 40 light years away in the constellation Andromeda, and it is visible to the naked eye if you know where to look. The star is roughly 2.5 billion years old (about half the age of our Solar System), and is about 30% more massive (and roughly twice as luminous) as the Sun. Because of Upsilon Andromedae’s relative brightness and proximity, and because it is fairly similar to the Sun, it was among the original group of 50-odd stars surveyed for planets by Marcy and Butler with the radial velocity technique. It was also one of the first stars found to be accompanied by an extrasolar planet (see Butler et al. 1997, “Three New 51-Peg Type Planets”, Astrophysical Journal Letters, 474, 115-118).
One of the recently published Upsilon Andromedae data sets (Fischer et al. 2003) contains 253 individual radial velocity measurements that span more than 15 years of observation, and is therefore capable of giving the systemic console a very informative workout. In this tutorial, we step throught the procedure for “discovering” the planets that exist in the Fischer et al (2003) radial velocity data. We assume that you have worked through the first tutorial, which analyzed the much less extensive data set for HD 4208, and which introduced the basic concepts of radial velocity fitting.
To begin, bring up the console, and select “upsand” from the system menu. You’ll need to scroll down a little ways to find it:
A click displays the complete Upsilon Andromedae radial velocity data set (obtained with the 3-meter telescope on Mt. Hamilton, see photo below) in the console window. Notice how the data-taking rate has increased with time. After Upsilon Andromedae was discovered to harbor an interesting planetary system (in early 1999), it has been observed more frequently. The radial velocity surveys, with their limited budgets of observing time, continually need to balace the confliciting demands of obtaining as many observations as possible of particularly interesting systems, while simultaneously devoting as uniform a coverage as possible to as many different stars as possible.
In the case of the star HD 4208, we could use the console to find a decent one-planet fit to the radial velocity data set simply by manipulating the sliders, and then applying a final “polish” to get a best fit. This was possible because the planet accompanying HD 4208 has only been observed for two orbits, and because it is apparently the only large planet orbiting that star. For Upsilon Andromedae, however, the complexity of the radial velocity data set hints at the presence of several significant planets with different periods. To get a foothold, we can compute a “Lomb-Scargle” periodogram of the data. This tells us which frequencies (that is, which periods) are most prevelant in the data. Compute the periodogram by clicking the “show” button. The Lomb-Scargle algorithm involves a fair amount of work for the computer. On a Mac 2.5Ghz G5, it takes a few seconds for the completed, densely sampled periodogram to appear in its own window:
The periodogram indicates that a number of different periodicities are present in the data. The most important peaks are listed to the right of the plot. The higher the power, the more important the peak. In this case, the strongest peak is at 4.617 days. We make the hypothesis that the highest peak is due to a planet, and we therefore activate the first row of orbit sliders (using the button to the far left of the row).
By typing in the box, enter 4.617 into the period controller for planet #1. Note that if you use the slider, you will have difficulty getting this exact value to appear.
With a 4.617 day period, the full radial velocity plot on the console window becomes totally useless. The prospective planet has executed more than a thousand orbits over the length of time (~15 years) that the star has been observed. Use the zoom and scroll sliders beneath the plot to focus in on the well-sampled region near 11,400 days (JD-2440000):
Arghh! You’ve uncovered a deficiency in the initial “beta” release of the program. The stellar radial velocity curve arising from a single planet on a circular orbit should be a smooth sine-wave, and yet the plot displays an unattractive sawtooth-like pattern. This results from undersampling of the plotted radial velocity waveform. You can circumvent this problem by adjusting the number of points that are plotted in order to display a smooth radial velocity curve under high magnification. Entering a value of 25000 into the “Set crvSamplePoints” box and hitting return will generate a smooth curve at the expense of slightly degraded console performance.
Use the zoom button and back out to a slightly lower level of magnification:
In the above screenshot, the curve now looks smooth, but its envelope is still affected by undersampling. If you see this artifact, don’t worry. It’s associated only with the “Set crvSamplePoints” value that you give to the plotting algorithm. The value of crvSamplePoints doesn’t affect the correct operation of the fitting procedure. Under the hood, the program always correctly computes the radial velocity curve wherever it needs to be evaluated.
A quick inspection of the radial velocities suggests that more than one planet may be contributing to the overall radial velocity signal. There are trends in the distribution of points which clearly unfold over periods of hundreds of days and more. Ignore these longer-term trends for the moment, and focus on pinning down the orbit of the 4.617 day planet.
Click on the “Mean Anomaly” line minimization button for planet one:
This automatically varies the mean anomaly (defined at the epoch of the first radial velocity data point) through the full 360 degree range, and selects the value that gives the best fit to the data. The chi-square of the fit drops to 77.49, and the RMS scatter drops to 60.77 meters per second.
Next, by clicking on the line minimization buttons for Period, Mass, Mean Anomaly, and Stellar Offset, you can try to improve the one-planet fit. A few cycles of clicking brings the chi-square down to 65.8 (the exact values that you calculate might deviate slightly from this, depending on history of operations that you’ve performed). A final round of improvement can be obtained by sending the Period, Mass, Mean Anomaly, and Stellar Offset parameters to the Levenberg-Marquardt multi-parameter minimization algorithm, (i.e. by checking the square left hand buttons and then clicking the “polish” button):
The resulting model is by no means a satisfactory fit to the data! Upsilon Andromedae is clearly accompanied by a ~0.7 Jupiter mass planet on a 4.167 day orbit, but there is more going on. Indeed, when Marcy and Butler announced their discovery of the 4.167 d planet in 1997, on the basis of the first 36 radial velocities, they noted that,
“The 36 velocities show evidence for variability in the gamma velocity of the sine wave on a timescale of about 2 yr with amplitude of ∼50 m s-1.”
Back in 1997, that was all they could say, but with the large number of additional velocities that have been published since then, we are in much better shape to figure out what is going on. Indeed, there are several ways that one can proceed to improve the fit. One strategy might be to explore non-zero orbital eccentricities for the planet. We choose not to do this for the following reason: If a 4.167 day planet were given a significant (say e>0.05) eccentricity, then it would experience intense tidal heating via the same process that forces Jupiter’s moon Io to be incredibly volcanically active, and the orbital eccentricity would be driven toward zero on a timescale that is very short in comparison to the lifetime of the star.
An alternate way to proceed is to search for additional planets in the system. To do this, we want to subtract off the radial velocity signal arising from the short-period 4.617 day planet, and look for the remaining periodicities in the system. This is accomplished by clicking on the “periodogram of the residuals” button, which operates just as advertised, yielding:
In this periodogram, the strongest peak in the data occurs at 242.394 days, with the second-highest peak falling at 1290 days. There are also peaks near exactly one day, which are aliases arising from the fact that the star was often observed with a 24-hour night-to-night cadence. Let’s assume that there is a second planet in the system with a 242.394 day period. Activate the second planet, enter the period, and then use the line minimizers to improve the fit (leave the eccentricities of both planets at zero for the moment). After a few rounds of clicking, the two planet verison of the system rolls into a local minimum of parameter space.
The fit can be improved a bit more by polishing with the Levenberg-Marquardt algorithm, but clearly, however, there may be additional bodies contributing to the radial velocity variation of the star. To find what might still be lurking, compute a periodogram of the residuals to the 2-planet fit (by clicking “update” in the residuals periodogram window).
The 242 day peak has been completely eliminated, leaving a single peak at 1290.0 days that contains a lot of power. Activiate a third planet with this period, and polish the fit by line minimizations (start with the Mean Anomaly) followed by multi-parameter minimization (polish) using Levenberg-Marquardt.
With three planets on the books, the fit has improved significantly, but is not yet adequate. At this point, one could try to add a fourth planet based on the result of the residuals periodogram. Alternately, one could open the system up to eccentric orbits. To explore eccentricity, pull the eccentricity sliders of the outer two planets up to modest non-zero values (for example e=0.1 for both). Then improve the radial velocity fit with line minimization followed by another application of the Levenberg-Marquardt algorithm.
A successful model of the system looks like this (full-res .pdf here):
By clicking on “Orbital View”, you can see the actual top-down view of the co-planar model system:
