THE LICK-CARNEGIE EXOPLANET SURVEY: A URANUS-MASS FOURTH PLANET FOR GJ 876 IN AN EXTRASOLAR LAPLACE CONFIGURATION^*

The planetary system orbiting the nearby M4V star GJ 876 (HIP 113020) has proven to be perhaps the most remarkable entry in the emerging galactic planetary census. This otherwise unassuming red dwarf has produced the first example of a giant planet orbiting a low-mass star, the first instance of a mean-motion resonance (MMR) among planets, the first clear-cut astrometric detection of an extrasolar planet, and one of the first examples of a planet in the hitherto unknown mass regime between Earth and Uranus.

The radial velocity (RV) variations of GJ 876 have been monitored by us, the Lick-Carnegie Exoplanet Survey (LCES) team, at the Keck I telescope using the High Resolution Echelle Spectrograph (HIRES) for 12.6 years. Bean & Seifahrt (2009) and Rivera et al. (2005, R05 henceforth) give detailed discussions of the discoveries and some of the studies made for this system. Here, we summarize some of that history.

Marcy et al. (1998) and Delfosse et al. (1998) announced the first companion, "b." They found that it has an orbital period (P_b) of ∼61 days and a minimum mass (m_bsin i_b) of ∼2.1 M_Jup and that it produced a reflex barycentric velocity variation of its host star of amplitude K_b ∼ 240 m s⁻¹. After 2.5 more years of additional Doppler monitoring, Marcy et al. (2001) announced the discovery of a second giant planet in the system. This second companion, "c," has an orbital period, P_c ∼ 30 days, m_csin i_c ∼ 0.56 M_Jup, and K_c ∼ 81 m s⁻¹, and upon discovery of the second planet, the parameters of the first planet were tangibly revised. Marcy et al. (2001) modeled the two-planet system with non-interacting Keplerian orbits and noted both that there was room for improvement in the quality of the fit and that the system's stability depended on the initial positions of the planets.

To a degree that has not yet been observed for any other planetary system, the mutual perturbations between the planets are dynamically significant over observable timescales. This fortuitous situation arises because the minimum masses are relatively large compared to the star's mass, because the planetary periods are near the 2:1 commensurability and because their orbits are confined to a small region around the star.

Laughlin & Chambers (2001) and Rivera & Lissauer (2001) independently developed self-consistent "Newtonian" fitting schemes which incorporate the mutual perturbations among the planets in fitting the RV data. The inclusion of planet–planet interactions resulted in a substantially improved fit to the RV data and suggested that the coplanar inclination of the b–c pair lies near i = 50°. These papers were followed by an article by Nauenberg (2002) who described a similar dynamical fitting method but, in contrast, found a system inclination near i = 90°.

Soon after the announcement of planet "c," Benedict et al. (2002) used the Fine Guidance Sensor on the Hubble Space Telescope to detect the astrometric wobble induced by planet "b," which constituted the first unambiguous astrometric detection of an extrasolar planet. Their analysis suggested that the orbital inclination of planet "b" is close to edge-on (i_b = 84° ± 6°), a result that agreed with the model by Nauenberg (2002), but which was in conflict with the results of Laughlin & Chambers (2001) and Rivera & Lissauer (2001).

R05 analyzed an updated RV data set that included new Doppler velocities obtained at Keck between 2001 and 2005. By adopting self-consistent fits, they were able to announce the discovery of a third companion, "d," with period P_d ∼ 1.94 days, m_dsin i_d ∼ 5.9 M_⊕, and K_d ∼ 6.5 m s⁻¹. Using only the RV data, they were able to constrain the inclination of the system (assuming that all three planets are coplanar) relative to the plane of the sky to 50° ± 3°. For this inclination, the true mass of the third companion is ∼7.5 M_⊕.

Bean & Seifahrt (2009) performed self-consistent fitting of both the Keck RV data from R05 and the astrometry from Benedict et al. (2002). They found that the joint data set supports a system inclination i = 48.9⁺¹⁸₋₁₆, which is in agreement with the i = 50° ± 3° value published by R05. Bean & Seifahrt's analysis (see, e.g., their Figure 1) indicates that the best-fit χ²_ν is not significantly affected by the inclusion of the astrometric data and that the astrometric detection of "b" by Benedict et al. (2002) is not in conflict with the significantly inclined system configuration. Bean & Seifahrt's work also resulted in the first determination of a mutual inclination between the orbits of two planets in an extrasolar planetary system. By allowing the inclinations, i_b and i_c and the nodes, Ω_b and Ω_c to float as free parameters in their three-planet dynamical fit, they derived a mutual inclination, cos Φ_bc = cos i_bcos i_c + sin i_bsin i_ccos(Ω_c − Ω_b), of Φ_bc = 5.0⁺³⁹₋₂₃.

Correia et al. (2009) also performed self-consistent, mutually inclined fits using the RVs from R05 plus 52 additional high precision RV measurements taken with the HARPS spectrograph. They confirmed the presence of companion "d" and further showed that it has a significant orbital eccentricity, e_d = 0.139 ± 0.032. Correia et al. (2009) determined inclinations i_b = 4893 ± 097 and i_b = 4807 ± 206. Their orbital fit yields a mutual inclination Φ_bc = 100, consistent to within measurement uncertainty with a coplanar system configuration. We became aware of Correia et al. (2009) while preparing this article for publication. A. Correia (2009, private communication) graciously shared the HARPS RV data on which their paper is partly based. We have carried out a preliminary analysis that indicates that our new results are not in conflict with these HARPS data. The HARPS RVs do not add much to constrain the system's dynamics because of poor phase coverage at the period of the fourth planet to be discussed below. Throughout this article, however, our analysis and results are based on the Keck RVs alone.

In addition to R05, several other studies, such as West (1999), West (2001), Laughlin et al. (2005), and Shankland et al. (2006), have examined the possibility of transits occurring in this system. The star has also been included in surveys in which it was probed for nearby companions or a circumstellar dust disk. These include Leinert et al. (1997), Patience et al. (2002), Hinz et al. (2002), Trilling et al. (2000), Lestrade et al. (2006), Shankland et al. (2008), and Bryden et al. (2009). Luhman & Jayawardhana (2002) were able to use Keck AO to place limits on potential brown dwarfs in the system at separations of 1–10 AU. Seager & Deming (2009) used the Spitzer spacecraft to search for thermal emission from planet "d." The star was also included in the Gemini Deep Planet Survey (Lafreniere et al. 2007). In addition to some of the work above, a variety of studies have examined the dynamics and/or long-term stability of the system. Examples are Kinoshita & Nakai (2001), Gozdziewski & Maciejewski (2001), Gozdziewski et al. (2002), Snellgrove et al. (2001), Ji et al. (2002), Jones et al. (2001), Zhou & Sun (2003), and Haghighipour et al. (2003). Beauge & Michtchenko (2003) present an analytic model which does a very good job at reproducing the dynamical evolution of the system. Murray et al. (2002), Lee & Peale (2002), Chiang et al. (2002), Thommes & Lissauer (2003), Beauge et al. (2003), Kley et al. (2004), Lee (2004), Kley et al. (2005), and Lee & Thommes (2009) examined scenarios of how the system may have come to be in its current dynamical configuration. Even interior structure models of GJ 876 d have been introduced (Valencia et al. 2007).

In this work, we present a detailed examination of our significantly augmented and improved HIRES RV data set from the Keck observatory. Our analysis indicates the presence of a fourth, Uranus-mass planet in the GJ 876 system. The new planet is in a Laplace resonance with the giant planets "b" and "c," and the system marks the first example of a three-body resonance among extrasolar planets. We expect that the newly resolved presence of the fourth planet will provide additional constraints on the range of formation scenarios that could have given rise to the system and will thus have a broad impact on current theories for planetary formation. The plan for this paper is as follows. In Section 2, we describe the new velocities and the methods used to obtain them. In Section 3, we present our latest two- & three-planet coplanar fits. In Section 4, we thoroughly analyze the residuals of our best three-planet coplanar fits. In Section 5, we discuss our four-planet fits. In Section 6, we explore the possibility of constraining the libration amplitude of the critical angle for the Laplace resonance. In Section 7, we explore the stability of additional planets in the system. Finally, in Section 8, we finish with a summary and a discussion of our findings.

The stellar characteristics of GJ 876 have been described previously in Marcy et al. (1998) and Laughlin et al. (2005). It has a Hipparcos distance of 4.69 pc (Perryman et al. 1997). Its luminosity is 0.0124 L_☉. As in previous studies, we adopt a stellar mass of 0.32 M_☉ and a radius of 0.3 R_☉ based on the mass–luminosity relationship of Henry & McCarthy (1993). We do not incorporate uncertainties in the star's mass (0.32 ±  0.03 M_☉) into the uncertainties in planetary masses and semimajor axes quoted herein. The age of the star exceeds 1 Gyr (Marcy et al. 1998).

We searched for Doppler variability using repeated, high resolution spectra with resolving power R ≈ 70000, obtained with Keck/HIRES (Vogt et al. 1994). The Keck spectra span the wavelength range from 3700 to 8000 Å. An iodine absorption cell provides wavelength calibration and the instrumental profile from 5000 to 6200 Å (Marcy & Butler 1992; Butler et al. 1996). Typical signal-to-noise ratios are 100 per pixel for GJ 876. At Keck, we routinely obtain Doppler precision of better than 3–5 m s⁻¹ for V = 10 M dwarfs. Exposure times for GJ 876 and other V = 10 M dwarfs were typically 8 minutes. Over the past year, these typical exposure times have been raised to 10 minutes.

The internal uncertainties in the velocities are judged from the velocity agreement among the approximately 700 2 Å  chunks of the echelle spectrum, each chunk yielding an independent Doppler shift. The internal velocity uncertainty of a given measurement is the uncertainty in the mean of the ∼700 velocities from one echelle spectrum.

We present results of N-body fits to the RV data taken at the W. M. Keck telescope from 1997 June to 2010 January. The 162 measured RVs are listed in Table 1. The median of the internal uncertainties is 2.0 m s⁻¹. Comparison of these velocities with those presented in Laughlin et al. (2005) and in R05 shows significant changes (typically 3–10 m s⁻¹) in the velocities at several observing epochs. R05 discuss the primary reasons for the improvements—a CCD upgrade and an improved data reduction pipeline. Additionally, many of the velocities presented here are the result of averaging multiple exposures over two-hour time bins. More importantly, the present RV data set is based on a new "super" template. For the new template, 10 exposures, most of which were <900 s, were obtained under very good observing conditions. In contrast, the template used to determine the RVs and uncertainties in R05 was based on only three 900-s exposures taken under relatively poor observing conditions. As a result, the new template has brought the median internal uncertainty down to the current 2.0 m s⁻¹ from the value of 4.1 m s⁻¹ in R05. Additionally, the new template has resulted in significant changes in the old published RVs. Note that since the last observation listed in R05, 36 more high quality observations have been obtained for this system since 2004 December. Our LCES survey assigns high observing priority to GJ 876, given the potential for further characterizing the intricate dynamical configuration between the 30- and 61-day planets and because of the possibility of detecting additional planets in the system.

We used similar procedures to those described in R05 to produce updated two- and three-planet dynamical fits to the GJ 876 Doppler velocities. The epoch for all the fits in this work is JD 2450602.093. As in R05, we initially held the eccentricity of companion "d" (e_d) fixed. We find, however, that fitting for e_d results in a significant improvement in all the fits which include companion "d." We therefore allow e_d to float in our final fits. Uncertainties in the tables listing our best-fit Jacobi parameters are based on 1000 bootstrap realizations of the RV data set. We take the uncertainties to be the standard deviations of the fitted parameters to the 1000 trials. The fitting algorithm is a union of Levenberg–Marquardt (LM) minimization (Section 15.5 of Press et al. 1992) and the Mercury integration package (Chambers 1999).

Our dynamical fits can be integrated to determine the libration amplitudes of the critical arguments of the 2:1 MMR between "b" and "c" and to monitor the linear secular coupling between the pericenter longitudes. For each fit, we integrated the system for 300,000 days (∼821 years) with a time step of 0.05 day. For all simulations in this work, we used the Hybrid symplectic algorithm in the Mercury integration package (Chambers 1999) modified as in Lissauer & Rivera (2001) to include the partial first order post-Newtonian correction to the central star's gravitational potential. The three relevant angles are φ_cb,c = λ_c − 2λ_b + ϖ_c, φ_cb,b = λ_c − 2λ_b + ϖ_b, and φ_cb = φ_cb,b − φ_cb,c = ϖ_b − ϖ_c. For each of these angles, we measure the amplitude as half the difference between the maximum and minimum values obtained during each ∼821-year simulation. As expected, we find that for both two- and three-planet fits to the real data set and for every bootstrapped data set, all three angles librate about 0°. We use the subscript "cb" to indicate the angles relating the longitudes of planets "c" and "b." In Section 5, we will use the subscripts "be" and "ce" to indicate the corresponding angles relating the longitudes of planets "b" and "e," and "c" and "e," respectively.

As in R05, we inclined the coplanar system to the line of sight. We set the inclination relative to the plane of the sky of all the companions to values from 90° to 40° in decrements of 1°. We hold all the longitudes of the ascending node fixed at 0°. We fit for all the remaining parameters. Finally, we examine χ²_ν as a function of the system's inclination. Figure 1 shows the results.

Zoom In Zoom Out Reset image size

The minimum in χ²_ν as a function of the system's inclination is in the range 57°–61°. We find that the location of the minimum has a small dependence on the number of companions in the system. For two planets, we find i = 57° ±  58. For three planets, we find i = 59° ± 32. For four planets, we find i = 61° ± 24. (The results for four-planet fits are presented in Section 5). The error estimates on the coplanar inclination are derived from 100 bootstrap trials of the RV data set. For each bootstrap RV set, we use the best-fit parameters to the real RVs at each inclination value as the initial guess. For each trial, we find the minimum in χ²_ν as a function of the system's inclination. The uncertainties in the fitted inclination of the system are the standard deviations of the locations of the minimum in χ²_ν. Our best three- and four-planet dynamical fits indicate a system with an invariable plane with an inclination relative to the plane of the sky of ∼59°, and we adopt this inclination as a starting point for studying the fitness of models with mutually inclined orbits.

Starting from the coplanar, i = 59°, two-planet fit, we examine the effect of a mutual inclination between companions "b" and "c." We do this by fitting for the inclinations of "b" and "c" (i_b and i_c) and the longitude of the ascending node of "c" (Ω_c). Since the system is rotationally invariant about the line of sight, we leave the longitude of the ascending node of "b" (Ω_b) fixed at 0°. We find χ²_ν = 9.76 (rms = 5.704 m s⁻¹) for the coplanar fit and χ²_ν = 9.91 (rms = 5.706 m s⁻¹) for the mutually inclined fit. This "improvement" is not significant. As we show below, fitting for the mutual inclination of "b" and "c" in the case of four-planet configurations results in a similar "improvement" in χ²_ν. For the two-planet mutually inclined model, the fitted masses and inclinations are m_c = 0.7211 M_Jup, m_b = 2.3689 M_Jup, i_c = 5794, i_b = 5529, and Ω_c = 082. With the approximation that Ω_c ∼ 0°, the invariable plane of this two-planet system has an inclination relative to the plane of the sky of 56°. The mutual inclination for this system is 27.

For the three-planet model with i = 59°, our best coplanar fit has χ²_ν = 3.93 (rms = 3.625 m s⁻¹) and the mutually inclined fit has χ²_ν = 3.99 (rms = 3.620 m s⁻¹). For the mutually inclined fit, the invariable plane of the system has an inclination relative to the plane of the sky of 592. For the mutually inclined fit, we force planet "d" to be in the invariable plane determined by "b" and "c." The mutual inclination between "b" and "c" in this configuration is 07. Since a mutual inclination does not significantly improve the three-planet fit, in Table 2 we show the best-fit parameters for our best coplanar three-planet fit.

Table 2. Three-planet Coplanar Fit for GJ 876 with i = 59°

Parameter	Planet d	Planet c	Planet b
P (days)	1.937781 ± 0.000024	30.0880 ± 0.0091	61.1139 ± 0.0084
ma	6.78 ± 0.51 M⊕	0.7175 ± 0.0042 MJup	2.2743 ± 0.0059 MJup
aa (AU)	0.02080665 ± 0.00000018	0.129590 ± 0.000026	0.208301 ± 0.000019
K (m s−1)	6.60 ± 0.47	88.72 ± 0.52	213.86 ± 0.55
e	0.257 ± 0.070	0.25493 ± 0.00080	0.0292 ± 0.0015
ω (°)	229 ± 28	48.67 ± 0.82	50.7 ± 4.5
MA (°)	358 ± 29	295.0 ± 1.0	325.4 ± 4.6
Offset (m s−1)	50.7 ± 0.4
Fit epoch (JD)	2450602.093
χ2ν	3.9280
rms (m s−1)	3.6254
φcb,c	4.5 ± 07
φcb,b	13.1 ± 32
φcb	11.5 ± 37

Note. ^aQuoted uncertainties in planetary masses and semimajor axes do not incorporate the uncertainty in the mass of the star.

Download table as:  ASCIITypeset image

At i = 59°, the angles φ_cb,c, φ_cb,b, and φ_cb are all in libration for both the model fit (in Table 2) and for all of the bootstrap trials. For the fit in Table 2, these amplitudes are φ_cb,c = 45 ±  07, φ_cb,b = 131 ±  32, and φ_cb = 115 ±  37. Note that the libration amplitudes for the three-planet fits presented in this section are smaller than for the three-planet fits presented in R05.

The critical angles show a complicated evolution for the coplanar three-planet fit with i = 59°. The first angle, φ_cb,c, has peaks at (listed with decreasing power) 1.58, 45.67, and 8.69 years. The second angle, φ_cb,b, has peaks at 8.69, 1.58, and 45.74 years. Finally, the third angle, φ_cb, has two peaks at 8.69 and 45.67 years. These periodicities are also present in the evolution of the eccentricities of the two outer planets. The eccentricities of all three planets show very regular secular evolution.

We performed an error-weighted Lomb–Scargle periodogram (Gilliland & Baliunas 1987) analysis on the residuals of the coplanar three-planet fits for both i = 90° and i = 59°. Both periodograms show a significant peak at ∼125 days. Figure 2 shows the periodogram for the case i = 59°. The periodogram for the case i = 90° looks very similar. For i = 90°, the peak raw power is 311.4. For i = 59°, the peak raw power is 310.2. Even if we allow for a mutual inclination between "b" and "c" and an overall trend as in Correia et al. (2009) (they actually subtract out a trend, due to the geometrical effect of the star's proper motion, prior to fitting the RVs), we still see a strong signal near 125 days in the periodogram of the residuals. For both values of i, folding the three-planet residuals at the peak period shows a strong coherent sinusoidal signature. Figure 3 shows the folded residuals for the case i = 59°. Examination of all the inclined three-planet fits shows that the peak period in the residuals is at 125 days for all values of i ⩾ 40°. For the 1000 fits used in determining the uncertainties in Table 2, the residuals of 941 of them have a peak period near 125 days, with 12 harboring a marginally higher alias peak with a period near 1.0 day. The 125-day periodicity has been adequately sampled at all phases only relatively recently.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also examined subsets of the RV data where we took just the first N observations, assumed i = 59°, and fit for the parameters of the three known planets. We took all values of N from 98 to 162. We then examined the error-weighted Lomb–Scargle periodogram of the residuals for all values of N. All the periodograms have a peak period near 125 days. In Figure 3, we show the folded residuals for the first 98 observations as squares. Comparison between the folded residuals for the first 98 points and that for the full RV set again suggests that the 125-day periodicity has been sampled at all phases only recently. The observed power at 125 days has grown monotonically during the last decade, but the growth rate has been non-uniform. The power grows more quickly whenever we observe around the time when the fourth planet is near the phases in its orbit when its orbital velocity has a large component along the line of sight. The power grows less quickly whenever we observe the system when the orbital motion of the fourth planet is mostly perpendicular to the line of sight. We are constrained to observe near the time of the full Moon, which impedes the efficiency of obtaining uniform phase coverage for periodicities near 120 days. As a consequence, there were several stretches where we consistently observed near the phases when the fourth planet's orbital motion was primarily perpendicular to the line of sight. This circumstance is one of the reasons why so many observations over 12 years were required prior to reliably detecting the fourth periodicity in this system.

As an additional check on the reality of the 125-day signal, we ran a Monte Carlo analysis. We take the coplanar three-planet model at i = 59° and integrate forward to generate a synthetic RV curve. We then generate a set of 1000 synthetic data sets by repeatedly sampling the synthetic RV curve at the times of the Keck observations and adding Gaussian noise perturbations of amplitude equal to the rms of the fit in Table 2. Then, for each RV set, we fit for the three planets and examine the periodogram of the residuals. In not one instance did the power near 125 days exceed the value observed for the real data. Additionally, the peak period is "near" 125 days in only one of the 1000 cases. Thus, the false alarm probability is <0.001.

An alternative to the fourth planet hypothesis is that the 125-day periodicity is due to rotation of the star itself. The rotation period of GJ 876 is at least ∼40 days, based on its narrow spectral lines and its low chromospheric emission at Ca ii H&K (Delfosse et al. 1998). R05 presented photometric evidence of a rotation period of ∼97 days. Our measurements of the Mt. Wilson S activity index also confirms this rotation period. Thus, rotational modulation of surface features cannot explain the 125-day period in the velocities. The periodogram of our measured Mt. Wilson S activity index after subtracting out a 14-year spot cycle has a peak period near 88 days.

The 125-day periodicity in the residuals to our best coplanar three-planet fit motivates exploration of the possibility that a fourth companion exists in the system. For the past several years, as we accumulated new RV measurements, we monitored the quality of dynamical fits containing a fourth planet with P ∼ 125 days. From late 2005 through 2008, these fits preferred a fourth planet with e_e>0.2 and produced system configurations which, when integrated, were dynamically unstable. Only recently, have the four-planet fits resulted in systems which are stable on timescales exceeding a few hundred million years (at least for coplanar configurations). A fourth planet, with a mass similar to Uranus, and on a stable, low-eccentricity orbit, has now definitively emerged from the data.

A primary reason why the latest fits result in a stable system is that we have only recently gathered enough observations to sample all phases of the 125-day periodicity. Figure 3 shows that most of the (small) phase gaps have now been filled in. Additionally, near the points with larger residuals, we also have additional observations with smaller residuals. As a result, the LM routine no longer prefers a fourth planet with period ∼125 days with a relatively large eccentricity. Another reason why it has been difficult to find stable four-planet fits with the fourth planet at 125 days is that the region around the outer 2:1 MMR with companion "b" contains the boundary between stable and unstable orbits for the two-planet system consisting of only "c" and "b" (Rivera & Haghighipour 2007).

To highlight the generally tenuous stability of planets with periods near 125 days, we performed simulations based on the fit in Table 2 in which we added 1200 (massless) test particles on circular orbits in the plane of the system with periods in the range 120–125 days every 0.5 day plus at the fitted period for "e." At each period, we placed 100 test particles equally spaced in initial longitude. In agreement with Rivera & Haghighipour (2007), nearly all are unstable on timescales <10⁶ years. Only one particle survived more than 10⁶ years; however, it was also lost after 7.4 Myr. It started with a period of 122.5 days (astrocentric) at a mean anomaly of 2196. For this long-term survivor the critical angle, φ_Laplace = λ_c − 3λ_b + 2λ_e, librates around 0° with an amplitude of Δφ_Laplace = 22° for the entire time that the particle is stable. The chaotic nature of the particle's orbit becomes apparent when the libration amplitude of the Laplace angle suddenly becomes unconstrained during the last few thousand years of its evolution. In Figure 4, we show the stability of the test particles as a function of initial semimajor axis and longitude. For clarity, we show the astrocentric semimajor axes raised to the sixth power. Unstable initial locations are indicated with open dark gray symbols. The single "stable" initial location, where a test particle survived 7.4 Myr, is indicated with a filled-in black symbol. With an open star, we also mark the approximate fitted location for the fourth planet from the coplanar four-planet fit with i = 59° (for simplicity, we assume zero eccentricity for the fourth planet—the astrocentric eccentricity is actually 0.0448). We also find that if we assign a mass of a few Earth masses to these "test particles," they tend to survive for longer times. The results strongly suggest that a resonant mechanism is required for an object in this region to be stable.

Zoom In Zoom Out Reset image size

Additionally, we find that the orbital parameters of the innermost planet, "d," can influence the stability and chaotic nature of test particles placed in the region near 125 days. Assuming a circular orbit for "d" results in an island of stability near the fitted location of planet "e." Instead of only a single particle which survives for 7.4 Myr, six particles survive for at least 10 Myr in this region if e_d is assumed to be zero. As for the long-lived particle for the case with e_d>0, all the "stable" particles are in the Laplace resonance. The six stable particles under the assumption e_d = 0 are shown as filled-in light gray symbols in Figure 4. We also used the Lyapunov estimator in the SWIFT integration package (Levison & Duncan 1994) to analyze the chaotic nature of test particle orbits with periods near 125 days. We find that the most stable particles in this region have Lyaponov times in the range 10³–10⁴ years.

Following procedures similar to the ones used for the two- and three-planet fits, we worked up four four-planet fits. The best-fit coplanar inclination (i = 59°) with e_d = 0 has χ²_ν = 2.8738 (rms = 3.0655 m s⁻¹). Working from this fit, holding e_d fixed, but fitting for a mutual inclination between "b" and "c" gives χ²_ν = 2.8322 (rms = 3.0219 m s⁻¹). If instead we fit for e_d and assume a coplanar system, the best-fit inclination is at i = 61° and χ²_ν = 2.5991 (rms = 2.9492 m s⁻¹). Finally, fitting for e_d and a mutual inclination between "b" and "c" results in χ²_ν = 2.6098 (rms = 2.9303 m s⁻¹). Since fitting for e_d and the mutual inclination between "b" and "c" results in a significant improvement in χ²_ν, we take this last model to represent our current best fit to the present Keck RV data set. However, simulating this model as well as other mutually inclined models with similar values of χ²_ν shows that the resulting systems are unstable on timescales of ≲1 Myr. Also, with e_d allowed to float, the "improvement" in χ²_ν between the coplanar and mutually inclined models is not significant. The F-test probability that these two fits are statistically the same is 0.88. Thus, Table 3 shows the best-fit parameters for the coplanar four-planet fit with i = 59°, our preferred inclination for the system. The resulting system is stable for at least 400 Myr, and other coplanar fits with similar values of χ²_ν result in systems which are stable for hundreds of millions of years up to at least 1 Gyr.

For completeness, we note that for our best-fit mutually inclined system, the invariable plane of the system is inclined 595 to the plane of the sky and the mutual inclination between "b" and "c" is 37. Fitting for the inclinations of planets "d" and "e" does not improve χ²_ν significantly. This result for the inclination of "d" is in agreement with the similar result by Correia et al. (2009). The F-test probability for our best four-planet (mutually inclined) fit versus the three-planet coplanar fit with i = 59° and with a fitted e_d is 0.0066. This small probability indicates that the four-planet model is significantly different from the three-planet model. It should be noted that for the coplanar case with i = 59°, the fitted longitude of planet "e" is near the location of the long-lived particle discussed above.

The inclination of planet "e" does have a significant effect on the stability of the system based on our best mutually inclined fit. It also influences the libration amplitudes of the critical angles to be discussed below. The orbit of planet "e" must have an inclination and node such that the precession rates of both the nodes and longitudes of periastron of the outer three planets result in a relatively small libration amplitude for the critical angle of the Laplace resonance, φ_Laplace = λ_c − 3λ_b + 2λ_e. In other words, systems with a small forced inclination for "e" should be more stable than systems with a large forced inclination—through experimentation we find this to be true. The coplanar four-planet fits are generally more stable than the mutually inclined versions since a small libration amplitude for φ_Laplace depends only on the precession rates of the longitudes of periastron of the outer three planets and not on the rates of precession of their nodes. For the mutually inclined four-planet fit, we performed a brief search of the parameter space spanned by i_e and Ω_e as well as i_d and Ω_d. Experimentation shows that all four of these can influence the stability (and evolution) of this chaotic system. Also, the libration amplitudes of all the librating critical angles are smaller if we simply place "e" and "d" in the invariable plane determined by "b" and "c." Since we only explored a small region of the parameter space and since we found several equally good fits to the RV data with various values for the forced inclination of "e," other, more stable systems with smaller forced inclinations are likely to exist. One emerging pattern from the search of the parameter space is that smaller forced inclinations, smaller libration amplitudes, and more stable systems are more likely to occur in models in which the orbital plane of "e" is more closely aligned with that of "b." This pattern is also in rough agreement with the finding that coplanar systems are generally more stable than mutually inclined ones.

As for the two- and three-planet fits, the critical angles involving "c" and "b" librate about 0° for all our four-planet fits as well as for all 1000 bootstrap fits used in determining the uncertainties in Table 3. However, their amplitudes are generally larger because of the perturbing influence of the fourth planet.

Since the ∼125-day period is near the 2:1 MMR with planet "b," the 4:1 MMR with planet "c," and the Laplace resonance with both "b" and "c," it is informative to examine all the relevant critical angles associated with these resonances as well as those associated with the resonances involving only "b" and "c." For the 2:1 MMR and the linear secular coupling between "b" and "e," the relevant angles are φ_be,b = λ_b − 2λ_e + ϖ_b, φ_be,e = λ_b − 2λ_e + ϖ_e, and φ_be = φ_be,e − φ_be,b = ϖ_e − ϖ_b. For the 4:1 MMR between "c" and "e," there are four relevant angles: φ_ce0 = λ_c − 4λ_e + 3ϖ_c, φ_ce1 = λ_c − 4λ_e + 2ϖ_c + ϖ_e, φ_ce2 = λ_c − 4λ_e + ϖ_c + 2ϖ_e, and φ_ce3 = λ_c − 4λ_e + 3ϖ_e. The subscripts above correspond to the multiplier in front of ϖ_e. There is also the critical angle involving just the longitudes of periastron of "c" and "e": φ_ce = ϖ_e − ϖ_c. Finally, the critical angle for the Laplace resonance is φ_Laplace = λ_c − 3λ_b + 2λ_e. Table 4 lists the libration amplitudes of all these critical angles for the fit in Table 3. A "C" in the amplitude column indicates that the angle circulates (at least once) during the ∼821-year simulation. The time evolution of the planetary eccentricities and the Laplace angle for the fit in Table 3 are shown in Figures 5 and 6, respectively. The secular evolution of the eccentricities of planets "b," "c," and "d" is apparently slightly less regular as a result of the addition of planet "e," of which the eccentricity evolution is more chaotic than that of the other three planets. Additionally, the evolution of all the librating angles in stable four-planet fits appears chaotic; however, the resulting systems are stable for at least several hundreds of millions of years.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The addition of the low-mass planet "e" to the system model has only a modest effect on the dynamical characterization of the orbits of "b" and "c." Δφ_cb,c is increased from 45 ± 07 to 574 ± 085, and Δφ_cb increases from 115 ± 37 to 225 ± 48. (Additionally, their evolution appears less regular with the addition of "e.") In addition, the 4:1 argument φ_ce0 is observed to librate with Δφ_ce0 = 83° ± 25°, and the 2:1 argument φ_be,b librates with Δφ_be,b = 36° ± 13°. Most significantly, the system also participates in the three-body resonance, with Δφ_Laplace = 40° ± 13°.

Libration of φ_Laplace is required for stability. Long-term simulations based on bootstrap fits associated with four-planet coplanar fits indicate that systems for which φ_Laplace ≳ 140° are unstable on timescales <2 Myr. Other systems with large libration amplitudes for φ_Laplace, 1414, 1291, and 1194, are stable for only 25.9 Myr, 33.6 Myr, and 196.8 Myr, respectively. This suggests that the stability of the system also hinges on the libration amplitude of the critical angle for the Laplace resonance.

Figure 7 plots the planetary positions over 120 one-day intervals starting from JD 2450608.093 in a frame rotating at planet "b" 's precession rate, $\langle \dot{\varpi _b} \rangle =-0\mbox{.\!\!^\circ $}116165\ d^{-1}$. Four snapshots corresponding to successive half-orbits for planet "b" are shown. The plots trace the system's 1:2:4 commensurability at a time when the osculating eccentricity of "e," at e_e ∼ 0.064, is relatively high and when ϖ_b − ϖ_e ∼ 180°. During this phase, triple conjunctions occur when planet "e" is near apastron.

Zoom In Zoom Out Reset image size

For the dynamical configuration given in Table 3, numerical integration indicates that the periastron longitudes ϖ_c and ϖ_b regress more quickly than does ϖ_e. Over a timescale of several years, the three apses approach alignment, angular momentum is transferred to "e," and the eccentricity of "e"'s orbit declines nearly to zero. During this low-eccentricity phase, the arguments φ_be and φ_ce both circulate through 2π. The complicated dynamical environment for "e" results in dramatic changes in its eccentricity over readily observable timescales. Our fit to the data indicates that during the time span of the RV observations, the eccentricity of "e" has varied from e_e ∼ 0 to e_e ∼ 0.075 (see Figure 5).

In the previous section, we showed evidence for a fourth planet in the GJ 876 system, planet "e," which participates in a three-body resonance with planets "c" and "b." Our preferred coplanar model for the GJ 876 planetary system which is consistent with the RV data shows that the critical angle for the Laplace resonance librates about 0° with amplitude 40°. Given the low RV half-amplitude for planet "e," it is natural to question whether there exists sufficient data to truly confirm the existence of the three-body resonance between "b," "c," and "e." In this section, we address the questions concerning the reality of this dynamical configuration. The analysis in this section is another method to address the question of whether the 125-day period is actually due to a planet.

We performed a Monte Carlo simulation in which we took a three-planet model and added a ∼125-day signal plus Gaussian noise. The three-planet model is just the three inner planets from the fit in Table 3. We added a Keplerian model with P = 120.713 days, m = 14.629 M_⊕, e = 0.0448, ω = 25136, and MA = 21688 plus Gaussian noise with rms = 5.9696 m s⁻¹. These (astrocentric) parameters are the result of subtracting off the effect of the three inner planets from the fit in Table 3 and fitting the residuals with a one-planet model (with initial guess based on the parameters of the fourth planet in our preferred four-planet fit) for which only the mean anomaly is allowed to float. We generated 1000 synthetic RV data sets and performed a four-planet Newtonian fit to each data set. We then integrated each fit for 300,000 days (∼821 years) and examined the libration amplitude of the critical angle for the Laplace resonance. Note (again) that the synthetic RVs are based on a Newtonian three-planet configuration plus a Keplerian signal plus Gaussian noise. Thus, we may expect that not too many fits should result in systems with a libration amplitude comparable to the system based on our preferred four-planet model. We find that only 6 simulations have a libration amplitude <40°, and 29 have a libration amplitude <53°, the one-sigma uncertainty on the Laplace angle libration width. Thus, most of these synthetic RVs result in systems which are not comparable to the system based on our preferred four-planet model. These results suggest that we can constrain the libration amplitude of the critical angle for the Laplace resonance. The results also indicate that planet "e" is strongly interacting with "b" and "c." This is in good accord with the eccentricity evolution shown in Figure 5.

Another Monte Carlo experiment we performed is to assume that the system is deep in the Laplace resonance and generate synthetic RVs to see how often we find systems comparable to the system based on our preferred four-planet fit. Among the fits to the bootstrap RVs above, we found one which results in a system in which the libration amplitude of the critical angle for the Laplace resonance is 164. We take this as our model and add Gaussian noise with rms = 3.9604 m s⁻¹, the rms for the fit in Table 3. We generate 1000 synthetic RV data sets and perform four-planet Newtonian fits to them. The 1000 fits are then integrated forward for ∼821 years. We find that 577 systems have a libration amplitude of the critical angle for the Laplace resonance <40°, 799 systems have this libration amplitude <53°, and none have the critical angle in circulation. The largest amplitude is ∼130°. Again, these results indicate that we are constraining the libration amplitude of the critical angle of the Laplace resonance.

Additionally, we examined the long-term stability of several bootstrapped and synthetic systems for which the critical angle for the Laplace resonance is not strictly librating. We find that such systems are unstable on timescales <200 Myr. Thus, our preferred coplanar fit to the current data set shows a system which is participating in a Laplace resonance, and this configuration appears to be required for long-term stability.

As noted in Figure 6, the libration amplitude for the critical angle for the Laplace resonance grows stochastically for long-term simulations performed up to at least 1 Gyr. If this diffusion occurs in the actual system and if we assume that the system were damped down initially, then the currently observed libration amplitude for the critical angle for the Laplace resonance may give a rough estimate of the age of the system.

Rivera & Haghighipour (2007) used massless test particles to show that, except for the innermost region, the entire habitable zone is unstable for small planets because of the perturbing influence of just planets "c" and "b." Based on their work, it is expected that interior to the orbit of "e" there may be a small region of stability just exterior to the orbit of "d." Exterior to the orbit of "e," we expect the region out to the location of the 2:1 MMR with "e" to be unstable. Beyond this location, we expect small planets on low-eccentricity orbits to be stable except possibly around the location(s) of the 3:1 MMR (and perhaps the 4:1 MMR) with "e."

We have performed a similar analysis based on the parameters in Table 3. In fact, due to the strong interactions among planets "e," "b," and "c," we find that particles placed exterior to the orbit of "d" are unstable out to beyond the location of the outer 5:2 MMR with "e." This result is similar to that found by Rivera & Lissauer (2000) in which the interactions between the outer planets in the ϒ Andromedae planetary system clear out a region much larger than what would be cleared out if only the outermost planet were present.

The analyses above were done with test particles placed on circular orbits starting on the x-axis. They suggest that the region among the planets is unstable. This is also suggested by the inherently chaotic nature of the four-planet system. A more thorough exploration of the initial positions of test particles, in which we allow for eccentric orbits with various periastron longitudes and various initial mean anomalies, can be used to find among the planets small islands of stability which are protected by some resonant mechanism. Correia et al. (2009) postulated that a small planet could exist with a period near 15 days. Such a planet is not detectable with the current RV data set(s). The addition of the fourth planet, "e," may have a significant effect on the stability of such a planet. However, we explored the possibility that such a planet could be participating in a four-body resonance with the three outer planets. Short-term simulations of particles placed in this region do indeed show that a small fraction of them are participating in the putative four-body resonance. A few of these test particles have the relevant critical angle, ϕ = λ_f − 4λ_c + 5λ_b − 2λ_e, executing small amplitude libration, Δϕ ∼ 30°, about ϕ ∼ ±70°. This configuration would keep the potential fifth planet away from the line joining the outer three planets and the star when triple conjunctions occur among "c," "b," and "e." Again, the present RV data set(s) do not support such a planet, but it would be interesting to discover a resonant configuration which has no solar system analog.

Using 162 RV observations taken at the Keck I telescope with HIRES, we have shown strong evidence for a fourth periodicity that is present in the RV of GJ 876. This likely corresponds to a planet "e" with period ∼124 days and a minimum mass of ∼13.4 M_⊕. This planet joins the previously known planets "d," "c," and "b" with periods of 1.9378, 30.1, and 61.1 days, respectively. Assuming that the four-planet system is coplanar, we have shown that the inclination of the system relative to the plane of the sky is ∼59°. Our best self-consistent fit shows that the period and mass of the new companion are 124.2 days and 15.4 M_⊕. The value of χ²_ν above 1.0 for our fit in Table 3 suggests that the stellar jitter may be at the 1–2 m s⁻¹ level and holds out the possibility of additional low-amplitude objects in the system which may be revealed with additional observations.

We have also shown that this new planet likely participates in a complicated set of resonances involving planets "c" and/or "b." Simulations indicate that it is in a 2:1 MMR with planet "b," a 4:1 MMR with planet "c," and a Laplace resonance with both "c" and "b." Comparison with prior four-planet fits using fewer RV observations indicates that as more RVs have been accumulated, the "fitted" libration amplitudes of the critical angles associated with these resonances have generally decreased. Additionally, long-term simulations based on our fits indicate that the complicated dynamical structure is necessary for the system's long-term survival. It is important to note that our current best-fit coplanar model may change as we and other groups gather more RV data on this fascinating system. Thus, the system may be deeper in the resonances mentioned above. Additionally, future RVs may eventually constrain the mutual inclination between planets "b" and "c."

Phenomenal signal-to-noise ratio and a decade-plus observational baseline combine to make the GJ 876 system a touchstone for studies of planetary formation and evolution. The presence of well-characterized planet–planet resonances can potentially allow a detailed connection to be made between the properties of GJ 876's protoplanetary disk and the planetary system that eventually emerged. Clearly, the system was prolific, both in terms of the mass and the total number of planets produced.

A zeroth-order prediction of the core-accretion paradigm for giant planet formation is that the frequency of readily detectable giant planets should increase with both increasing stellar metallicity and with increasing stellar mass (Laughlin et al. 2004; Ida & Lin 2004, 2005). During the past decade, both of these trends have been established observationally (see, e.g., Fischer & Valenti 2005 for a discussion of the metallicity trend and Johnson et al. 2008 and Bowler et al. 2010 for discussions of the mass trend.)

Until recently, however, there appeared to be little evidence for the strong expected planet-metallicity correlation among the handful of M dwarf stars that are known to harbor giant planets. In particular, the results of Bonfils et al. (2005) suggested that GJ 876 has subsolar metallicity. This result naturally induces speculation that GJ 876's giant planets might be the outcome of gravitational instability (e.g., Boss 2000) rather than core accretion.

Johnson & Apps (2009) have recently provided an update to the metallicity calibration developed by Bonfils et al. (2005). The Johnson & Apps (2009) calibration indicates that the planet-bearing M dwarfs do appear to be systematically metal-rich, suggesting that there is no breakdown of the planet-metallicity correlation as one progresses into the red dwarf regime. In particular, the new calibration indicates [Fe/H] = 0.39; GJ 876 is more than twice as metal-rich as the Sun.

A supersolar metallicity for GJ 876 dovetails with a history of vigorously efficient planet formation, but there remain a number of fascinating questions with regards to the formation and migration mechanisms that permitted the final configuration to be assembled. We are preparing an in-depth analysis of these issues, even as we continue to gather more observations of this most eminently productive star.

We sincerely thank Jacob Bean for refereeing and improving this work. S.S.V. gratefully acknowledges support from NSF grants AST-0307493 and AST-0908870 and from the NASA KECK PI program. R.P.B. gratefully acknowledges support from NASA OSS grant NNX07AR40G, the NASA Keck PI program, the Carnegie Institution of Washington, and the NAI, NASA Astrobiology Institute. G.L. acknowledges support from NSF grant AST-0449986. N.H. gratefully acknowledges support from NASA EXOB grant NNX09AN05G, the NASA Astrobiology Institute under cooperative agreement NNA04CC08A at the Institute for Astronomy, University of Hawaii, and a Theodore Dunham J. grant administered by Funds for Astrophysics Research, Inc. We also acknowledge the major contributions of fellow members of our previous California-Carnegie Exoplanet team: Geoff Marcy, Debra Fischer, Jason Wright, Katie Peek, and Andrew Howard, in helping us to obtain many of the pre-2008 RVs presented in this paper. The work herein is based on observations obtained at the W. M. Keck Observatory, which is operated jointly by the University of California and the California Institute of Technology, and we thank the UC-Keck, UH-Keck, and NASA-Keck Time Assignment Committees for their support. We also thank those of Hawaiian ancestry on whose sacred mountain of Mauna Kea we are privileged to be guests. Without their generous hospitality, the Keck observations presented herein would not have been possible. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France.

Note added in proof. Some of the scatter in Figure 3 is due to the evolution of the orbit of planet "e". We expect that because of this evolution, the signature of planet "e" may indeed disappear as future RVs are added.