{'fraction_non_alphanumeric': 0.035240503168211756, 'fraction_numerical': 0.029341074382744952, 'mean_word_length': 4.432931999321689, 'pattern_counts': {'":': 0, '<': 8, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The GJ 876 Planetary System -- A Progress Report

We present an updated analysis of the GJ 876 planetary system based on an augmented data set that incorporates 65 new high-precision radial velocities obtained with the Keck telescope from 2001 to 2004. These new radial velocities permit a more accurate characterization of the planet-planet interactions exhibited by the system. Self-consistent three-body orbital fits (which incorporate both the estimated instrumental uncertainties and 6 m/s Gaussian stellar jitter) continue to show that GJ 876 b and GJ 876 c are participating in a stable and symmetric 2:1 resonance condition in which the lowest order, eccentricity type mean-motion resonance variables are all librating. The planets are also locked in a secular resonance which causes them to librate about apsidal alignment. The small libration widths of all three resonances likely point to a dissipative history of differential migration for the two planets in the system. (for full abstract see paper).


Introduction
GJ 876 (HIP 113020) is the lowest mass star currently known to harbor planets, and it is accompanied by perhaps the most remarkable exoplanetary system discovered to date. In 1998, Marcy et al. (1998) and Delfosse et al. (1998) announced the discovery of a P ∼ 60 d companion orbiting the star. This planet, designated GJ 876 "b" is a super-Jovian object, with M sin(i) = 1.9M JUP , and it induces a large (K ∼ 210 m s −1 ) radial velocity variation in its red dwarf companion. After continued Doppler monitoring of GJ 876, Marcy et al. (2001) announced the discovery of a second (M sin(i) = 0.6 M JUP ) planet in the system. This object (designated GJ 876 "c") has a P ∼ 30 d orbit, and was identified to be participating in a 2:1 mean-motion resonance with the outer planet "b".
The red dwarf GJ 876 (RA=22 53, DEC=14 16) is observable from both hemispheres, and is distinguished by being the fortieth-nearest stellar system, with a Hipparcos-determined distance of 4.69 pc (Perryman et al. 1997). Its spectral type is M4 V. Using the bolometric correction of Delfosse et al. (1998), the Hipparcos-estimated parallax indicates a stellar luminosity of 0.0124L ⊙ . The red-dwarf mass-luminosity relation of Henry & McCarthy (1993) therefore implies a mass of 0.32 M ⊙ , and an estimated radius of R ⋆ = 0.3 R ⊙ .
A definitive identification of the resonance conditions obeyed by the planets is made possible by the large dynamic range of the GJ 876 radial velocities. Among the 113 Doppler velocities obtained with the Keck telescope, the individually estimated instrumental errors have an average precision of 4.65 m s −1 , with individual precision estimates ranging as low as 2.3 m s −1 . The two planets induce velocity swings in the star of nearly 0.5 km s −1 , and thus allow us to take full advantage of the fine Doppler precision. Furthermore, the outer planet has been observed for more than forty orbital periods. These fortuitous circumstances allow the planet-planet interactions to be probed with a degree of refinement that is exceeded only for the planets in the solar system (e.g. Laplace 1799Laplace -1802 and by the planets orbiting the 6.2 ms radio pulsar PSR B 1257+12 (Wolszczan & Frail 1992, Konacki & Wolszczan 2003. The gravitational perturbations exerted by the planets on each other induce a significant non-Keplerian component to the orbital motion. In particular, the periastra of the planets precess at a rate̟ ∼ −41 • yr −1 . The non-Keplerian aspects of the motion lead to a relatively high best-fit value (currently χ 2 = 2.81) for models that use dual-Keplerian fitting functions to model the observed radial velocity variation. The strong planet-planet perturbations do, however, enable the construction of dynamical fits to the radial velocity data that both improve the χ 2 statistic of the orbital fit, and place the planets into the secular apsidal alignment resonance and deeply within the two co-planar 2:1 mean-motion resonances (Laughlin & Chambers 2001, Nauenberg 2002. The existence of this multiply resonant configuration can be understood as the consequence of differential migration of the two planets within GJ 876's protoplanetary disk (e.g. Lee & Peale 2001, 2002Lee 2004), and the presence of strong mutual interactions between the two planets leads to a partial removal of the so-called sin(i) degeneracy. System configurations in which the planetary orbits are inclined by less than 30 • to the plane of the sky exhibit significantly worse fits to the radial velocity data.
The importance of the GJ 876 system arises because it can provide interesting constraints for theories of planetary formation. Because the current resonant state is sensitive to details of the system's history, it is worthwhile to evaluate the degree of confidence that can be placed in the present best-fit orbital parameters. The plan of this short paper is thus as follows: In §2, we describe co-planar fits to the radial velocity data set. These fits allow us to construct a detailed model of the system, and equally importantly, allow us to evaluate the confidence for which we can determine the orbital parameters. In §3, we broaden our analysis to include the possibility of system configurations in which the planets do not orbit in the same plane. In this case, the orbital angular momentum vectors of the planets can precess, and the planets can be potentially observed to evolve through transiting configurations. In §4, we briefly discuss how our results bear on current theoretical studies of the nascent GJ 876 planetary system.

Co-planar configurations of the Two Planets
We first assume that the planets GJ 876 "b" and "c" are in a co-planar configuration perpendicular to the plane of the sky (i s = i b = i c = 90 • ), and obtain self-consistent three-body fits to the combined Lick-Keck radial velocity data set This set includes the 16 Lick velocities listed in Marcy et al. (2001), and the 113 Keck velocities listed in Table 1. All of our orbital fits are obtained using a Levenberg-Marquardt multi-parameter minimization algorithm (Press et al. 1992) driving a three-body integrator as described in Laughlin & Chambers (2001). The best edge-on coplanar fit is listed in Tables 2 and 3. This reference fit has twelve free parameters, including the planetary periods, P b and P c , the mean anomalies, M b and M c at epoch JD 2449679.6316, the orbital eccentricities, e b and e c , the longitudes of periapse ̟ b and ̟ c , the planetary masses, m b and m c , and two velocity offsets, o 1 , and o 2 . The quantity o 1 is an offset velocity added to the first GJ 876 radial velocity, v L1 , obtained with the Lick 3-meter telescope (t = JD2449679.6316, v L1 = 58.07 m s −1 ; see Marcy et al. 2001). The parameter o 2 is an offset velocity added to all of the radial velocities taken with the Lick Telescope. It accommodates the different velocity zero-points of the Lick and Keck telescopes. The mass of the star is fixed at 0.32 M ⊙ . The mean longitudes, λ i (used to compute the mean-motion resonant arguments) are related to the mean anomalies and longitudes of periapse The orbital elements listed in Table 2 are osculating orbital elements at epoch JD 2449679.6316 (the epoch of the first radial velocity point taken at Lick Observatory in 1994.9 listed by Marcy et al. 2001) and are expressed in Jacobi coordinates. As explained in Lee &Peale 2003, Jacobi coordinates are the most natural system for expressing multiple-planet fits to radial velocity data. For reference, in Table 3, we express the orbital configuration of the system (again at JD 2449679.6316) in Cartesian coordinates. In the Cartesian system, the line of sight from the Earth to the Star is in the negative y-direction, and the y-component of velocity for the star relative to the system center of mass is measured, by convention, as a negative radial velocity.
The uncertainties in the orbital fit are estimated using Monte Carlo simulation of synthetic data sets (Press et al. 1992). In this procedure, we assume that the true configuration of the system is that given by the orbital parameters listed in Tables 2 and 3 (i.e. the best-fit co-planar, edge-on system). We produce 100 synthetic data sets by integrating this assumed planetary configuration forward in time, sampling the stellar reflex velocity at all of the observed epochs, and adding (in quadrature) noise drawn from Gaussian distributions corresponding to (1) an assumed σ = 6 m s −1 stellar jitter and (2) the individual velocity errors. Chromospherically quiet G and K dwarfs in the ongoing radial velocity surveys typically show RMS scatters σ ∼ 3−5m s −1 arising from stellar jitter (Saar, Butler, & Marcy 1998). Nauenberg (2002) argued that excess scatter in the dynamical fits to the GJ 876 radial velocity data should be attributed to a stellar jitter of 2 − 4m s −1 . Preliminary work by Wright & Marcy (2004) indicates that M3-M4 dwarfs with chromospheric activity similar to GJ 876 display typical jitter values σ = 4 ± 2 m s −1 , motivating our conservative choice of σ = 6 m s −1 .
The assumption of a 6 m s −1 stellar jitter gives an average χ 2 = 1.51 ± 0.09 for the fits to the synthetic data sets, consistent with the value χ 2 = 1.535 obtained from fitting to the actual data set. (All the χ 2 values that we quote are computed using only the instrumental uncertainties, and do not include the scatter expected to arise from stellar jitter). The uncertainty quoted for each orbital parameter is the variance computed for the parameter from the fits to the 100 synthetic data sets. We find that the distributions of parameter estimates are generally consistent with underlying Gaussian distributions. We note, however, that significant co-variation does exist between some of the orbital parameters (e.g. M i and ̟ i ), making it impossible to generate systems that are fully consistent with the radial velocity data by independently sampling orbital parameters from the inferred underlying distributions.
We conclude that the nominal, edge-on, coplanar two-planet model of Table 2 is fully consistent with the current set of radial velocity measurements of the star. If the actual stellar jitter for GJ 876 is smaller than 6 m s −1 , then one can contemplate extracting additional information (related, say, to the inclinations and nodes of planets "b" and "c", or to additional bodies) from the lists of radial velocities in Table 1 and in Marcy et al. (2001).
The three-body fit to the radial velocities indicates that the two major planets in the GJ 876 are locked in a symmetric configuration, with the resonant arguments θ 1 = λ c − 2λ b + ̟ c , and θ 2 = λ c − 2λ b + ̟ b both librating about zero degrees. The orbital configuration of the best-fit edgeon coplanar model of the GJ 876 data set is shown in Figure 1. In this figure, the positions of the planets are plotted as filled circles at 60 successive one-half-day intervals beginning on JD 2449710, when the planets were both near periastron. The positions are plotted in the frame centered on the star. Also plotted (as clouds of dots) are the positions of the planets at every one-half-day interval since the epoch of the first Lick Observatory data point taken on JD 2449679.6316. The figure shows that the orbits of the planets do not close, while examination of the time-dependent osculating orbital elements shows that the periapses of the planets are precessing at a rate oḟ ̟ = −41 • yr −1 . This rapid precession is the primary reason why Keplerian fits to the data show higher χ 2 values than the self-consistent three-body fits. Figure 2 shows the fitted reflex velocity of the star in comparison with the radial velocity data. The most striking feature of this figure (aside from the dominant ∼ 60 d periodicity) is the modulation arising from the 8.7 year precession period for the planets' joint line of apsides. This precession has now been observed for more than one full period, and the planets have completed a full librational cycle for both the secular |̟ c − ̟ b | resonance argument, as well as the 2:1 resonance arguments θ 1 and θ 2 . These librations are manifest in the slightly non-sinusoidal envelope of the overall stellar reflex velocity. The non-Keplerian aspect of the motion is also evident in the wave of small-amplitude velocity reversals running through the radial velocity curve. In the summed Keplerian model, this wave has an asymmetric shape, and is produced (along with the overall modulation) by the inner planet "c" having a fixed period P c = 30.12 d that is slightly less than half the P b = 61.02 d period of the outer planet. In the self-consistent fit, the small velocity reversals display a symmetric waveform, and arise largely from the precession of the inner eccentric orbit and the librations about the three resonances. For additional related discussion of the manifestation of the orbital dynamics in the radial velocity curve, see Nauenberg (2002) The primary assumption underlying the fit given in Tables 2 and 3 is that the planetary orbits are co-planar and are being viewed edge-on. While there is no a-priori observational evidence to indicate that the system is co-planar, it is likely that the planets arose from a relatively flat protoplanetary accretion disk. Numerical integrations of the differential migration of the system which assume this scenario show that the eccentricity must in general be forced to higher values than observed before significant mutual inclination is excited (see also Thommes & Lissauer 2003). Hence it makes dynamical sense to prefer co-planar models. We also note that astrometric evidence obtained by Benedict et al. (2002) suggests that the outer planet in the system is being viewed fairly close to an edge-on configuration.
If we assume co-planar inclinations with i s < 90 • and construct a succession of fits, we obtain the run of best-fit χ 2 values shown by the thick dashed line of Figure 3. This sequence shows a very slight decline in the value of χ 2 as the system is tilted from i s = 90 • ( χ 2 = 1.535) to i s = 59 • ( χ 2 = 1.525). For co-planar inclinations having i s < 38 • , however, χ 2 experiences a rapid rise. Similar behavior in the χ 2 profile was observed by both Laughlin & Chambers (2001) and , although with more radial velocity data, the dip in χ 2 has grown shallower. Laughlin & Chambers (2001), and  both intepreted the configuration with the minimum χ 2 as representing the likely coplanar inclination of the system, whereas Nauenberg (2002) suggested that the improvement found by those authors in going from i s = 90 • to i s ∼ 45 • was not significant. Our primary aim is thus to ascertain what significance can be ascribed to this trend in χ 2 . That is, which co-planar inclinations can be ruled out by the dynamical fits to the data?
In Figure 4, we plot the best-fit osculating eccentricities, e b and e c as a function of co-planar inclination 90 − i s . As i s decreases from 90 • , the fitted planetary masses increase by sin(i) −1 , and the fitted eccentricities also increase. The inner planet eccentricity, for example, increases from e c = 0.22 at i s = 90 • to e c = 0.38 at i s = 20 • . The best fit χ 2 value, however, changes very little in the face of this large eccentricity increase. This behavior occurs because the primary non-Keplerian interaction between the planets is the̟ = −41 • yr −1 precession (see e.g. Ford 2003). As the masses of the planets increase, the precession rate also increases. This increase, however, can be essentially exactly offset by an increase in the orbital eccentricities, which act to decrease̟.
For each of the 100 Monte-Carlo realizations of synthetic data sets which were previously generated for the edge-on coplanar system listed in Tables 2 and 3, we perform the same procedure of incrementing the co-planar inclination and obtaining fits. The results are shown as the cloud of dots in Figure 3, in which thirteen randomly selected sequences are also plotted as dark lines in order to give a representative idea of the trends for particular realizations. These fits show that the shallow minimum observed near i ∼ 60 • in the fits to the actual data cannot be believed.
Figures 5, 6 and 7 show the fitted values of the 2:1 and secular libration widths (|θ 1 | max , |θ 2 | max , and |̟ b − ̟ c | max ) as a function of co-planar inclination for the Monte-Carlo realizations. The fits to the actual data (heavy dashed lines) are fully consistent with the behavior observed in the Monte-Carlo realizations, providing further evidence that the inclination of the co-planar system cannot be confidently extracted from the data (assuming Gaussian stellar jitter with σ = 6 m s −1 ).
The libration width figures indicate that as the masses of the planets are increased (i.e. as 90 − i s increases) the libration widths |θ 1 | max , |θ 2 | max , and |̟ b − ̟ c | max all show a decrease, reaching minimum values near i s ∼ 45 • . This phenomenon occurs because the librations are more readily sensed in an radial velocity data set for planets of larger mass. Hence, a given observed perturbation must arise from a smaller libration if the planet masses are increased. The increase in χ 2 observed for systems with i s < 30 • is associated with the inability to match the observed perturbations with further decreases in the libration widths of the resonances. We note that simulations of resonant capture (Kley et al 2003) and differential migration favor narrow libration widths. These scenarios would therefore favor the prediction that the system will eventually be found to lie in the neighborhood of i s ∼ 45 • .

The Prospects for Observing GJ 876 "c" in Transit
The a priori probability that a planet on a Keplerian orbit transits its parent star as seen from the line of sight to Earth is given by, where a is the semi-major axis of the orbit, R ⋆ and R pl are the radii of the star and planet, respectively, e is the orbital eccentricity, and ̟ is the argument of periastron referenced to the intersection of the plane of the sky with the orbital plane, namely the line of nodes. For the inner planet of the GJ 876 system, this probability is only ∼ 1%, if we assume a stellar radius of R = 0.3R ⊙ . Planet-planet interactions in the GJ 876 system, however, allow the nodal line of the inner, less massive planet to precess into transit for a significantly wider variety of observationally consistent non-coplanar configurations. Transits of GJ 876 by the inner planet "c", if they occur, are therefore likely to be visible for a period of order two years as the node of the planetary orbit sweeps across the face of the star. The scientific opportunities from such transits would be somewhat analogous to the opportunities provided by the series of mutual eclipses observed in the Pluto-Charon system in the 1980s (Binzel 1989). Such configurations require a mutual inclination between planets "b" and "c". Because a transiting configuration is relatively easy to observe (the transit depth is expected to be of order 10%, this system makes an interesting photometric target for small-aperture telescopes (e.g. Seagroves et al. 2003).
Using the planetary evolution models computed by Bodenheimer, Laughlin, & Lin (2003), and assuming i s = 90 • , we estimate that the planetary radius of GJ 876 "c" should be 0.93 R JU P if the planet has a solid core, and 1.03 R JU P if it does not. Insolation-driven atmospheric-interior coupling, which can lead to an increased radius (see e.g. Guillot & Showman 2002), is expected to be negligible for planet "c". For an assumed tidal quality factor Q = 10 6 , the eccentricity damping timescale is of order 250 Gyr (Goldreich & Soter 1966), indicating that the energy generated by interior tidal heating should not affect the planetary radius. We estimate that the effective temperature at the planet's τ = 1 surface is 210K, assuming an albedo a = 0.4. Benedict et al. (2002) used the FGS instrument on HST to obtain a preliminary measurement of the inclination of the outer planet in the GJ 876 system, obtaining a value i b = 84 ± 6 • . In order to illustrate the possibility that the inner planet may periodically experience transit epochs, we assume that the orbital plane of the outer planet is coincident with the line of sight at JD 2449679.6316 (i b = 90 • ). We then choose (1) a specific value for the osculating inclination of the inner planet at the epoch of the first radial velocity point, as well as (2) the osculating value of the difference in nodal longitudes at the first radial velocity epoch. With these parameters fixed, we then obtain a self-consistent fit to the radial velocity data to determine all the other orbital parameters. When an acceptable fit is obtained, we integrate the system forward to check for the occurrence of (inner planet) transits within the next 100 years.
The results are shown in the lower left hand panel of Figure 8, which shows the result of 1296 such separate self-consistent fits. In the figure, the fits are organized by the choice of osculating starting inclination of the inner planet orbit (y-axis of the figure panels), and by the initial angle between the two ascending nodes (x-axis of the figure panels). The nominal edge-on co-planar system therefore corresponds to the bottom row of cells. Scenarios where the inner planet was transiting during the last season of observations (and specifically during the transit epoch near JD 2453000.57) have their cells colored black. Systems that start transiting within 100 years of the last radial velocity observation in Table 1 are indicated by dark gray (transits to start very soon) to light gray (transits starting in the year 2103). Some regions of the diagram contain systems which were transiting during the past ten years, but which have by now moved out of alignment. On average, over the range of configurations plotted in Figure 8, the line of nodes of the inner planet precesses at rates of order −4 • yr −1 .
The lower right panel of Figure 8 maps the distribution of χ 2 values obtained for the 1296 separate self-consistent fits. The lowest values found are χ 2 = 1.52, matching the best-fit coplanar i s = 90 • model. All of the models have χ 2 < 1.65. The Monte Carlo analysis of the previous section thus indicates that they are acceptable fits to the radial velocity data, assuming σ = 6 m s −1 . In the top two panels of Figure 8, we plot the libration widths of the secular apsidal alignment, and the θ 1 resonance argument for each fit. These panels show that, for the range of mutual inclinations sampled, the resonant conditions are always fulfilled.

Discussion
Our analysis continues to show that the Non-Keplerian interaction between the two planets in the GJ 876 system indicates that the planets are participating in both the 2:1 mean motion resonances, as well as in the secular apsidal resonance. Radial velocities accumulated over the last four years show that the libration widths of all three resonances are narrow, which argues for a dissipative history of differential migration for the system.
It is interesting to note, however, that the planet-planet interactions are in a sense quite subtle, and suffer from a degeneracy which prevents simultaneously accurate measurement of the eccentricity of the inner planet and the overall inclination of the system. Extensive Monte-Carlo simulations suggest that the eccentricity of the inner planet lies in the range 0.2 < e c < 0.35, and that the system has i s > 20 • . This situation is based on an assumption for the stellar jitter of 6 m s −1 . If this assumption turns out to be conservative, and the actual jitter is less, then it will be possible to obtain considerably better constraints on the orbital parameters of the system, and as more radial velocities are obtained, perhaps confirm or rule out the presence of additional small bodies in this remarkable exoplanetary system. Plausible and detailed histories for the origin of the resonances in the GJ 876 system were proposed by Lee & Peale (2001, 2002. In their scenario, the planets originally formed in low eccentricity orbits with semi-major axes larger than those currently observed, and with a larger period ratio than the present-day 2:1 commensurability. The planets then grew large enough to open gaps in the protoplanetary disk. Hydrodynamic simulations by Bryden et al. (2000), and Kley (2000) suggest that a residual ring of disk material between two massive planets is rapidly cleared as a consequence of repeated spiral shock passages from the protoplanetary wakes. This clearing process appears to require only several hundred orbits after the planets have been established. After the ring of gas between the planets has vanished, the planets will experience differential migration. The spiral wake driven through the outer disk will exert a negative torque on the outer planet, causing it to spiral inward. The inner planet will either be pushed outward by a remnant inner disk, or more likely, will retain a more or less constant semi-major axis. The inward-migrating outer planet then captures the inner planet into a low-order mean motion resonance (which in the case of GJ 876 was the 2:1) and the planets migrate in together. Lee & Peale (2002) demonstrated this mode of resonant capture for GJ 876 through the use of torqued three-body simulations. Additional N-body simulations of the GJ 876 precursor system were performed by a number of authors including Snellgrove, Papaloizou & Nelson (2001), Murray, Paskowitz & Holman (2002), Nelson & Papaloizou (2002), and Beaugé, Ferraz-Mello, & Michtchenko (2004). More recently, full hydrodynamical simulations by Papaloizou (2003), and Kley, Peitz & Bryden (2004) have also demonstrated capture of GJ 876 "b" and "c" into the observed resonances as a consequence of differential migration driven by disk torques.
Once the planets are migrating in resonance in response to outer disk torques, the orbits lose angular momentum and energy at different rates. In the absence of a dissipative mechanism, this mismatch causes the planetary eccentricities to increase. Lee & Peale (2002) introduced an ad-hoc eccentricity damping term to the migration. In cases where eccentricity damping was not used, they found that the semi-major axes decreased by only 7% before the eccentricities were pumped to their observed nominal values (e c = 0.22, and e b = 0.03). They therefore suggested that either (i) the disk dissipated before the planets were able to migrate very far, or alternately, that (ii) an effective mechanism exists for eccentricity damping during resonant migration.
Option (i) appears to require fine-tuning in order to provide an explanation for the current state of the GJ 876 system. The GJ 876 red dwarf, with M = 0.3M ⊙ , is nearly one hundred times less luminous than the Sun. The inner planet, GJ 876 "c", with its surface temperature T ∼ 210 K, is not far inside the location of the current snow-line of the GJ 876 system. For GJ 876 "b", located at a = 0.2AU , we estimate a temperature at τ = 1 of T b ∼ 160 K, which places it at or beyond the present snowline. The stellar evolution models of Baraffe et al. (2002), however, indicate that during contraction phases between 1 and 10 million years when giant planet formation likely took place, GJ 876 was more than ten times as luminous as it is now. The possibility of nearly in situ formation for the GJ 876 planets is therefore unlikely, but not fully out of the question (see e.g. the accretion models of Bodenheimer, Hubickyj, & Lissauer 2000). Certainly, the comparatively luminous early phases of M star evolution pose interesting tests for theories of planet formation.
Option (ii) may also be problematic. Recent 2D hydrodynamical simulations, such as those of Kley, Peitz, & Bryden (2004) are able to follow the planet-planet-disk evolution over secular timescales t > 5 × 10 4 yr These simulations self-consistently model both the resonance capture and differential migration processes, and show that eccentricity damping arising from the disk gas is much smaller than that required by Lee & Peale (2002) to explain the current state of the GJ 876 system as arising from significant differential migration. Kley, Peitz, & Bryden (2004) remark, however, that it remains to be seen whether 3-D hydrodynamic calculations, which incorporate a more realistic equation of state, and which adequately resolve the gas flow close to the planets, will provide the needed increase in the eccentricity damping rate.
The low expected temperature of GJ 876 "c" leads naturally to speculation that a potentially habitable terrestrial world might exist in the system. The usual definition of the planetary habitability zone, as given in Kasting et al. (1993), combined with the stellar properties of GJ 876, suggests that the habitable zone of GJ 876 is located interior to the orbit of planet "c" (a c = 0.13AU ) at a radius r h ∼ 0.1AU . Menou & Tabachnik (2003) report that terrestrial planets placed in habitable circular orbits with 0.1 AU < a < 0.2 AU are rapidly ejected by the outer two planets. We have verified this conclusion using the updated orbital parameters given in Table 2. We remark, however, that the clear history of resonant capture and inward dynamical migration in this system suggests that a terrestrial-mass object orbiting interior to the two gas giant planets may have been captured into a 2:1 resonant orbit with GJ 876 "c", leading to a high-eccentricity analog of the Laplacian resonant condition observed among Io, Europa and Ganymede. Such an object would have an orbital period of order P ∼ 15 days, and a semi-major axis of a t = 0.08AU . Numerical experiments show that stable systems of this sort are readily found in which the resonant argument θ 2 between the planet "c" and the putative interior terrestrial planet is librating, and where the eccentricity of the terrestrial planet is e t ∼ 0.3. If such a system is not fully co-planar, then one can expect precession of the nodal line, and hence periodically recurring transits. An Earth-size planet transiting GJ 876 would produce a transit depth of 0.3%, which is readily detectable with modest-aperture telescopes from the ground (Henry 1999). A habitable planet in the GJ 876 system would display a maxium separation from the primary star of approximately 20 mas, which places the system within the top 300 candidates among the 1139 nearby stars currently being considered as potential targets for NASA's TPF mission. 1 It is likely that the GJ 876 system will reveal further surprises as it is studied photometrically and spectroscopically from the ground and from space. Furthermore, even the present radial velocity data set may harbor much additional information if the stellar jitter turns out to be smaller than σ = 6 m s −1 that we have assumed in this study.
We thank Drs. Peter Bodenheimer, Eric Ford, Man Hoi Lee, Jack Lissauer, Mike Nauenberg and Eugenio Rivera for useful discussions. This material is based upon work supported by the National Aeronautics and Space Administration under Grant NNG-04G191G issued through the Terrestrial Planet Finder Foundation Science Mission (to GL). We acknowledge support by NSF grant AST-0307493 (to SSV), NSF grant AST-9988087 and travel support from the Carnegie Institution of Washington (to RPB), NASA grant NAG5-8299 and NSF grant AST95-20443 (to GWM), and by Sun Microsystems. This research has made use of the Simbad database, operated at CDS, Strasbourg, France.   Tables 2 and 3. Velocities obtained at Lick Observatory (listed in Marcy et al. 2001) and the velocities taken at Keck Observatory (listed in Table 1) are shown as small solid circles. The plotted Lick velocities include a fitted offset between the telescopes which resulted in o 2 = 44.476m s −1 being added to each of the 16 Lick Observatory measurements. Bottom Panel: Residuals to the orbital fit.