## Console Tutorial #3

The star gl 876 (also known as IL Aqr, and sometimes as GJ 876) is an exceedingly ordinary red dwarf. With its mass of roughly 1/3 that of the Sun, it is a typical representative of the most common population of stars in the Galaxy. Like all red dwarfs, gl 876 is tiresomely thrifty with its nuclear fuel. It fuses hydrogen into helium only 1% as quickly as the Sun, and because it is convective throughout its interior, it will shine steadily for nearly a trillion years before finally runnning short of hydrogen and swelling to become a modest red giant. At 15 light years distance, Gl 876 is currently the fortieth nearest known stellar system to the Sun. It lies in the constellation Aquarius, and is approximately 100 times too faint to be visible to the naked eye.

Despite being a very ordinary star, Gl 876 harbors an absolutely extraordinary planetary system. We can use the systemic console to explore what is going on.

The radial velocity data sets for Gl 876 have been published by Marcy et al. 2001, and by Rivera et al. 2005. In the published data sets, there are a total of 171 velocities, 16 taken at Lick Observatory during the 1990s, and 155 obtained at Keck observatory over the past eight years:

The Gl 876 observations come from two different telescopes (Lick and Keck). Because the two telescopes have different iodine cells and spectrographs, they have independent radial velocity zero points. this means that both data sets can be adjusted by adding or subtracting a constant baseline velocity from each. This is accomplished using the Stellar Offset slider (introduced in tutorial 1), as well as the gj876_1.vels slider. Both velocity offsets need to be treated as free parameters in a planetary fit to the data, and as such, are accompanied by their own line-minimization and Levenberg-Marquardt buttons:

The run of radial velocities over the past nine years of observations indicates that the star is accompanied by at least one planet that induces very significant radial velocity variations. A periodogram of the data shows that there is a very large amount of power at a period of 60.993 days:

Based on the suggestion provided by the periodogram, activate a single planet in orbit around the star. Using the line minimization buttons, optimize the planetary mass, period, mean anomaly, and telescope offset settings. Once the values have converged, send all of these parameters to the Levenberg-Marquardt algorithm (that is, check the left hand boxes by each slider display, and then click on the “polish” buton). The end result should be a 1.96 Jupiter mass planet with a period of 60.953 days, and a mean anomaly at the epoch of the first Lick radial velocity of 314 degrees:

Using the zooming and scrolling sliders, it is clear that many sections of the radial velocity data are fit quite well:

Overall, however, the one-planet fit has problems, and the reduced Chi-square is still absolutely terrible. Clearly, a more complicated fit to the data is required.

One possibility would be to improve the fit by allowing the orbit to become eccentric. This offers some immediate help, but a more effective strategy emerges if we compute a periodogram of the residuals:

The periodogram indicates that there is a huge power spike at 30.332 days. Interestingly, this spike was *not* present in the original periodogram of the radial velocities themselves. It only shows up in the residuals periodogram. This is a very strong indication that a second body is present in orbit around the star.

Activate a second planet with a period of 30.332 days, and use the line-minimization buttons (starting with the mean anomaly for the 30.332-day period planet) to work up a 2-planet fit to the data. For the time being, continue to leave the eccentricity of both planets at zero. After the fit stops changing when the line-minimization buttons are pressed, submit the active system parameters (that is, the periods, the mean anomalies, and the masses of both planets, along with the telescope offsets) for a final “polish” by the Levenberg-Marquardt algorithm. You should end up with something that looks like this:

The fit is considerably better. The presence of the two planets creates an undulating envelope in the radial velocity waveform, which does a good job of following the long-timescale trends in the data. The Chi-square value is much lower, but it is still well above unity. Clicking on orbital view shows an interesting configuration:

At the time when the first Lick radial velocity was obtained, the planets were on nearly opposite sides of the star, and because the period of the outer planet is nearly twice that of the inner planet, they will return to this configuration approximately every two months. The gradual drift in the phase of the alignments causes a beat frequency that generates the overall envelope of the radial velocity curve.

So far, all of our exploration with the console has assumed that the planetary motion consists of independent elliptical orbits. With this approximation, the model radial velocity curves are computed very quickly by solving Kepler’s Equation. In reality, however, the situation is more complex. In addition to feeling the gravitational force from the central star, the planets also pull on one another gravitationally, and perturb each other’s orbits. In particular, for a configuration like the one we’ve uncovered for the Gl 876 system, the planet’s gravitational pulls tend to add up over time because the planets repeatedly encounter one another at nearly the same point in their orbits.

Just as a small series of pushes can lead to a large-amplitude swing for a child on a playground, an ordered sequence of planet-planet pulls can cause a profound effect on the planetary orbits in which the planets participate in *resonance*. A classical example of gravitational resonance occurs between Io, Europa, and Ganymede in the Jovian satellite system. Io orbits two times for every one orbit of Europa, and Europa in turn orbits two times for every one orbit by Ganymede. The gravitational tugs between the moons conspire in a remarkable way to enforce the maintenance of a 4:2:1 orbital ratio, in which the orbital periods are driven through small oscillations around perfect commensurabiltiy. The mathematical details of the situation were first worked out by Laplace in the 1700s, shortly before the French Revolution.

When the “Integrate” button is clicked, the console switches from the approximation of perfect, fixed Keplerian orbits to a full integration (using a Runge-Kutta 4th/5th-order scheme with adaptive stepsize control) of Newton’s equations of motion for the planets and the star.

This switch makes the computed radial velocity curve fully self-consistent, with the drawback of being much more time-consuming to compute accurately. After the console has been switched to provide the radial velocity curve via integration, it will also display the radial velocity curve that it computed using the Keplerian approximation. By comparing the two curves, one can evaluate how much of an effect the inclusion of planet-planet interactions is having. For example, if you zoom in on the first five 60-day cycles of the system shown above, the difference between the two curves is noticeable, but still relatively small:

After a few tens of cycles, however, the curves are completely different, and the self-consistent radial velocity curve does not fit the data at all:

Can you recover the best two-planet fit to the system using the console? This requires patience (and an unbothered browser window) since the Levenberg-Marquardt algorithm takes a fair amount of time to run when full integration is enabled. The 2-planet co-planar fit to the gl 876 system is described fully in a recent paper by Laughlin et al. (2005). During the past nine years, the planets have swept through a full precessional cycle, as shown by this figure:

as well as this .mpg animation, which animates the shapes of the orbits over a period of 100 years.

For reference, the fit is shown below (and here with a higher-resolution .pdf file):

What do you see in the periodogram of the residuals? For several years, a very interesting signal was sitting in the published radial velocity data, and anyone (armed with the systemic console) could have taken advantage of it.