{'fraction_non_alphanumeric': 0.0362085838109884, 'fraction_numerical': 0.025630271680879817, 'mean_word_length': 4.250203881911597, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamics and Origin of the 2:1 Orbital Resonances of the GJ 876 Planets

(Abridged) A dynamical fit has placed the two planets about the star GJ 876 in coplanar orbits deep in 3 resonances at the 2:1 mean-motion commensurability with small libration amplitudes. The libration of both lowest order mean-motion resonance variables, theta_1 and theta_2, and the secular resonance variable, theta_3, about 0 deg. differs from the familiar geometry of the Io-Europa pair, where theta_2 and theta_3 librate about 180 deg. By considering a condition for stable simultaneous librations of theta_1 and theta_2, we show that the GJ 876 geometry results because of the large orbital eccentricities e_i, whereas the very small e_i in the Io-Europa system lead to the latter's geometry. Surprisingly, the GJ 876 resonance configuration remains stable for e_1 up to 0.86 and for amplitude of libration of theta_1 approaching 45 deg. with the current e_i. We find that inward migration of the outer planet of the GJ 876 system results in certain capture into the observed resonances if initially e_1<0.06 and e_2<0.03 and the migration rate |(da_2/dt)/a_2|<0.03(a_2/AU)^{-3/2} yr^{-1}. The bound on the migration rate is easily satisfied by migration due to planet-nebula interaction. If there is no eccentricity damping, eccentricity growth is rapid with continued migration within the resonance, with e_i exceeding the observed values after a further reduction in the semi-major axes a_i of only 7%. With eccentricity damping (de_i/dt)/e_i = -K|(da_i/dt)/a_i|, the e_i reach equilibrium values that remain constant for arbitrarily long migration within the resonances. The equilibrium e_i are close to the observed e_i for K=100 (K=10) if there is migration and damping of the outer planet only (of both planets). It is as yet unclear that planet-nebula interaction can produce the large value of K required to obtain the observed eccentricities.


INTRODUCTION
have discovered two planets about the nearby M dwarf star GJ 876. A preliminary fit of the stellar radial velocity (RV) variations due to two unperturbed Kepler orbits implies that the orbital periods of the two planets are nearly in the ratio 2:1. This resonance is an analog to the orbital resonances among the satellites of Jupiter and Saturn (e.g., Peale 1999), but it is the first to be discovered among extrasolar planetary systems. Table 1 shows the system parameters from the Marcy et al. analysis based on data taken at both the Keck and Lick Observatories and an adopted stellar mass of 0.32M ⊙ . The orbital periods are approximately 30 and 60 days. The reduced chi-square statistic of χ 2 ν = 1.88 indicates that the two-Kepler system is an adequate fit to the data. However, the rather large minimum masses of the planets of 0.56 and 1.89 M J (Jupiter masses) mean that the two-Kepler fit may not be a good determination of the system characteristics because the large mutual perturbations will ensure that the orbits deviate from Kepler orbits. In fact, substitution of the Marcy et al. parameters from Table 1 (with sin i = 1 for both planets, where i is the inclination of the orbital plane to the plane of the sky) into a calculation of the perturbed orbital motions beginning at the specified initial epoch leads to large variations in the orbital elements that were assumed constant in the fit. On the other hand, projecting the orbital parameters in Table 1 to those at a different epoch within the time span of the data set leads to much less variations in the orbital elements and apparent long term stability of the system (Marcy et al. 2001). For some choices of epoch, the system is not in the 2:1 resonance, and it eventually becomes unstable (Laughlin & Chambers 2001). Clearly, a dynamical fit of system characteristics to the RV observations is necessary to constrain the parameters that define an RV curve for GJ 876.
The publication of the data set by Marcy et al. (2001) allowed Laughlin & Chambers (2001) to perform such a fit. They assumed that the two planets are on coplanar orbits, and used two methods to minimize χ 2 ν as a function of the initial system parameters. Starting with the 2-Kepler fit, a Levenberg-Marquardt minimization scheme driving an N -body integrator was used to find a local minimum in χ 2 ν . A second method uses a genetic algorithm combined with a simple model for the variations of the orbital elements due to the 2:1 orbital resonance to search for the global minimum in χ 2 ν . Both methods converged to similar solutions. The dynamical fits are sensitive to sin i since the determined masses (and perturbations) grow as i is decreased. The best-fit solutions obtained using the Levenberg-Marquardt method are included in Table 1 for both the Keck+Lick data and Keck data alone. The Keck data is much more accurate than the Lick data and the solution for this data alone has a minimum χ 2 ν = 1.59 with an rms scatter of 6.86 m s −1 for sin i = 0.55. Although the Lick data is less accurate than the Keck data, the longer time span of observation that includes the Lick data may give the better solution. The solution for the Keck+Lick data has a minimum χ 2 ν = 1.46 with an rms scatter of 13.95 m s −1 for sin i = 0.78, but the χ 2 ν minimum is broad, with several nearby solutions giving nearly as good a fit. The conclusion of Laughlin & Chambers is that the broad minimum in χ 2 ν allows probable values of 0.5 < sin i < 0.8.  examined mutually inclined orbits and got similarly small χ 2 ν values for several different configurations, but they confirmed the Laughlin & Chambers results for coplanar orbits.) Figure 1 shows the trajectories of the motions for 3100 days (the average periapse precession period for both planets) for the two Laughlin-Chambers best-fit solutions in plots of e j sin θ j versus e j cos θ j , where e j is the eccentricity of the jth planet (with j = 1 and 2 for the inner and outer planets, respectively), are the two lowest order, eccentricity-type meanmotion resonance variables at the 2:1 commensurability, and λ j and ̟ j are the mean longitude and longitude of periapse of the jth planet. The two Laughlin-Chambers solutions place the system deep in both mean-motion resonances, with θ 1 and θ 2 librating about 0 • with remarkably small amplitudes. The simultaneous librations of θ 1 and θ 2 about 0 • mean that the secular resonance variable also librates about 0 • . Although the parameters Fig. 1.-Small amplitude librations of the two 2:1 mean-motion resonance variables, θ 1 = λ 1 − 2λ 2 + ̟ 1 and θ 2 = λ 1 −2λ 2 +̟ 2 , about 0 • for the GJ 876 planets. Trajectories for 3100 days (the average periapse precession period for both planets) are shown in plots of e j sin θ j versus e j cos θ j for the Laughlin & Chambers (2001) best-fit solutions to the Keck data alone and the combined Keck and Lick data.
will be more tightly defined as more data is acquired, the almost singular nature of the low amplitude librations indicates that the real system is most likely indeed locked in multi-resonance librations at the 2:1 mean-motion commensurability. This means that the line of apsides of the two orbits are nearly aligned with conjunctions of the two planets always occurring very close to the periapse longitudes. The resonance configuration ensures that the planets can never make close approaches in spite of their large masses, and barring external perturbations or an unreasonably high dissipation of tidal energy in the planets, the system should be stable for the lifetime of the star. In fact, we find that the librations are stable even for amplitudes of the inner planet resonance variable θ 1 approaching 45 • if we induce larger amplitudes of libration of the resonance variables by changing the initial value of the mean anomaly of the inner planet. The existence of the mean-motion resonances means that the assumption that the Fig. 2.-Variations in the semi-major axes, a 1 and a 2 , and eccentricities, e 1 and e 2 , of the GJ 876 planets for the Laughlin-Chambers Keck+Lick solution. The small amplitude librations of the resonance variables ensure that the semi-major axes and eccentricities have little variation.
orbits are nearly coplanar should be correct, and such coplanarity is consistent with the accretion of planetary bodies in a nebular disk surrounding the forming star. We shall assume that the orbits are coplanar hereinafter. Furthermore, there is not much to distinguish the two Laughlin-Chambers solutions in Figure 1, and we shall assume the parameters appropriate to the solution based on both the Keck and Lick data. Figure 2 shows the variations in the semi-major axes and eccentricities of both planets for 10 4 days for the Laughlin-Chambers Keck+Lick solution. The small amplitude librations of the resonance variables about 0 • ensure that both the eccentricities and semi-major axes have little variation in spite of the large mutual perturbations. This is why the 2-Kepler fit of Marcy et al. (2001) could produce a fit that was not too bad.
The fact that both 2:1 mean-motion resonance variables librate about 0 • in Figure 1 contrasts with the geometry of the Io-Europa 2:1 orbital resonances at Jupiter, where the resonance variable involving Io's longitude of periapse, θ 1 , li-brates about 0 • but the resonance variable involving Europa's longitude of periapse, θ 2 , librates about 180 • (e.g., Peale 1999). This means that in the Io-Europa case the lines of apsides are antialigned, with the periapses 180 • apart (i.e., θ 3 librates about 180 • ). Conjunctions occur when Io is near periapse and Europa is near apoapse. In §2 we explain why the GJ 876 and Io-Europa systems have different resonance configurations by considering a condition for stable simultaneous librations of the two mean-motion resonance variables. In particular, the longitudes of periapse should precess in the same direction at the same average rate, so that the relative alignment of the lines of apsides is maintained. In the Io-Europa case, the orbital eccentricities are sufficiently small that the precession rate d̟ i /dt is dominated by a single term lowest order in eccentricities and containing cos θ i . The overall sign of the coefficient for this term is < 0 for i = 1 (inner satellite) but > 0 for i = 2 (outer satellite), and retrograde precessions of both satellites require θ 1 to librate about 0 • and θ 2 to librate about 180 • . In the case of GJ 876, the eccentricities are sufficiently large that there are large contributions to the precession rates from higher order terms whose cosine arguments are linear combinations of the resonance variables. Summing over all contributing terms, coincident retrograde precessions of the GJ 876 planets require both θ 1 and θ 2 to librate about 0 • .
If the orbits of the two planets about GJ 876 were originally much further apart, with the ratio of their mean motions considerably greater than 2:1, any process that drives the orbits toward each other could lead to capture into the 2:1 orbital resonances that we observe. In §3 we describe a particular migration process due to the gravitational interaction between the planets and the nebular disk surrounding the young GJ 876. We note the simulations by Bryden et al. (2000) and Kley (2000), which show that two planets that are massive enough to open gaps individually in the gas disk can rather quickly clear out the disk material between them, if they are not separated too far. Disk material outside the outer planet exerts torques on the planet that are not opposed by disk material on the inside, and the outer planet migrates toward the star on the disk viscous time scale. Any disk material left on the inside of the inner planet exerts torques on the inner planet that push it away from the star. Thus the condition of approaching orbits necessary to form the resonances is established. We discuss the effect of planet-nebula interaction on orbital eccentricities and note that the GJ 876 system is in the interesting regime where it is uncertain whether eccentricity damping or growth is expected because of the rather large planet-star mass ratio of the outer planet. (Although the mass of the outer planet, M 2 , for the Laughlin-Chambers Keck+Lick solution is only 2. In §4 we present the results of a series of numerical orbit integrations in which the orbits of the GJ 876 planets, initially far from the 2:1 meanmotion commensurability, are forced to approach each other. (The numerical methods are described in the Appendix.) We show that capture into the two 2:1 mean-motion resonances and the secular resonance is certain if the orbital eccentricities start reasonably small and the rate of migration is not too fast. If there is no eccentricity damping, the eccentricities of both planets increase rapidly after resonance capture and exceed the observed values after a very short migration of the resonantly locked planets. Eccentricity damping of the formė i /e i ∝ȧ i /a i , where a dot over a symbol denotes d/dt andȧ i is the forced migration rate, leads to a termination in the eccentricity growth, and the eccentricities reach equilibrium values that remain constant for arbitrarily long migration in the resonances. We find that significant eccentricity damping with |ė i /e i | ≫ |ȧ i /a i | is required to produce the observed eccentricities of the GJ 876 system. In §5.1 we show that alternative damping of the eccentricities by tidal dissipation within the star or planets is insignificant, and in §5.2 we discuss other related studies. Our conclusions are summarized in §6.

COMPARISON WITH THE IO-EUROPA SYSTEM
As shown in Figure 1, the GJ 876 system has both lowest order, eccentricity-type mean-motion resonance variables at the 2:1 commensurability, θ 1 = λ 1 − 2λ 2 + ̟ 1 and θ 2 = λ 1 − 2λ 2 + ̟ 2 , and hence the secular resonance variable, θ 3 = ̟ 1 − ̟ 2 , librating about 0 • . This resonance configuration is different from that of Io and Europa, where θ 2 and θ 3 librate about 180 • . It should be noted that the differences are not due to the additional resonances involving Ganymede in the Io-Europa case. In the scenario where the resonances among the inner three Galilean satellites of Jupiter are assembled by differential tidal expansion of the orbits (Yoder 1979;Yoder & Peale 1981), Io is driven out most rapidly and the resonance variables involving Io and Europa only are captured into libration first. These resonance variables have the same centers of libration before and after the 2:1 commensurability between Europa and Ganymede is encountered. The differences in the resonance configurations of GJ 876 and Io-Europa are instead due to the magnitudes of the eccentricities involved and can be understood from a condition for stable simultaneous librations of the two mean-motion resonance variables.
This condition for stable simultaneous librations of θ 1 and θ 2 is that the longitudes of periapse, ̟ 1 and ̟ 2 , on average precess at the same rate. For coplanar orbits, the equations for the variation of ̟ i , a i , and e i in Jacobi coordinates are (e.g., Peale 1986) where is the mass of the star (or Jupiter), M i is the mass of the ith planet (or satellite), the Hamiltonian and the disturbing potential If we neglect terms of order (M 1 /M 0 ) 2 and higher and assume that a 1 < a 2 , Φ can be expanded to the form where and the summation is over the range 0 ≤ j ≤ ∞, −∞ ≤ k, ℓ ≤ ∞, and 0 ≤ m, n ≤ ∞. The coefficient C j kℓmn (β), where β = a 1 /a 2 , can be written in terms of the Laplace coefficient b j 1/2 (β) and its derivatives (Brouwer & Clemence 1961). For a given cosine argument φ jkℓ , the term lowest order in eccentricities is of order e |k| 1 e |ℓ| 2 . Within the 2:1 mean-motion resonances and the secular resonance, the perturbations are dominated by the terms whose arguments are nearly fixed because of the resonances, i.e., those terms involving the resonance variables θ 1 , θ 2 , θ 3 , and their linear combinations. Since the terms in equation (4) for da i /dt and equation (5) for de i /dt are proportional to ∂H/∂λ i and ∂H/∂̟ i and λ i and ̟ i appear in the cosine arguments only, the cosines are changed into sines and there is no secular change in a i and the forced e i if the resonance variables θ 1 , θ 2 , and θ 3 librate about either 0 • or 180 • . (Note, however, that if the eccentricities are not small and da i /dt and de i /dt are not dominated by a single term, there is also the possibility that the sums of all contributing terms are zero for resonance variables having values other than 0 • and 180 • .) We have calculated the contributions to the precession rate d̟ i /dt from secular terms in equation (8) with argument of the form kθ 3 (including the k = 0 non-resonant secular term), up to fourth order in eccentricities, and from meanmotion resonance terms with argument of the form kθ 1 − ℓθ 2 (k = ℓ), up to third order in eccentricities, and they are listed in Table 2. We set the resonance variables to their average values, i.e., the libration center values. In addition, we ignore the small deviation of β from 2 −2/3 (since the commensurability is not exact) and evaluate C j kℓmn (β) at β = 2 −2/3 . For Io and Europa, we adopt e 1 = 0.0026 and e 2 = 0.0013, which are the equilibrium eccentricities before the resonances with Ganymede are encountered in the tidal scenario (Yoder & Peale 1981). The numerical values in Table 2 are for the current semi-major axis of Io, but they can be scaled to other values of the semi-major axis since the precession rates are proportional to a −3/2 1 . As we can see in Table 2, because the orbital eccentricities of Io and Europa are small, the precession rate d̟ i /dt is dominated by a single term lowest order in the eccentricities. This is the term with argument θ i and proportional to e i in the disturbing potential (eq. [8]), and whose contribution to d̟ i /dt is proportional to 1/e i (eq. [3]). Since the coefficients of the e 1 cos θ 1 and e 2 cos θ 2 terms are C 2 1000 = −1.19 and C 1 0−100 = +0.43, respectively, retrograde precessions of both Io and Europa require θ 1 = 0 • and θ 2 = 180 • . Furthermore, the requirement that the regression rates are identical implies a simple relationship between the eccentricities: where n i is the mean motion and the secular motion̟ sec,i includes contributions from secular terms in the disturbing potential and the additional secular motion induced by the oblateness of Jupiter. As we can see in Table 2, the contributions from the secular terms in the disturbing potential and even the secular motion induced by the oblateness of Jupiter (which is +0.12 • day −1 for Io's orbit) are small compared with the ∼ −1.5 • day −1 precession induced by the first order mean-motion resonance terms with arguments θ 1 and θ 2 . Thus e 1 /e 2 ≈ −β −1/2 (C 2 1000 /C 1 0−100 )(M 2 /M 1 ) = 1.9, where the last equality is for the masses of Io and Europa.
For GJ 876, we know from numerical orbit integration of the Laughlin-Chambers Keck+Lick solution that on average d̟ 1 /dt = d̟ 2 /dt = −0.116 • day −1 . The numerical values in Table 2 are obtained using e 1 = 0.255, e 2 = 0.035, and a 1 = 0.130 AU, which are the average values for the Laughlin-Chambers Keck+Lick solution (see Fig. 2). As we would expect from the analysis above for the Io-Europa case, with θ 1 = θ 2 = 0 • , the contributions to d̟ 1 /dt and d̟ 2 /dt from the e 1 cos θ 1 and e 2 cos θ 2 terms, respectively, have opposite signs. However, because the eccentricities are large, these terms no longer dominate the precession rates. In particular, the largest contri-bution to d̟ 2 /dt comes from terms second order in eccentricities -the largest of which being the e 1 e 2 cos(θ 1 + θ 2 ) term, which is of order e 1 (= 0.255) lower than the e 2 cos θ 2 term, but whose coefficient is large and negative (C 3 1−100 = −4.97). We can also see from Table 2 that both d̟ 1 /dt and d̟ 2 /dt converge only slowly with the order of the kθ 1 − ℓθ 2 terms. With the inclusion of terms up to e 3 , the precession rates have the correct, negative signs when θ 1 = θ 2 = θ 3 = 0 • , but the magnitude of the precession rate for the inner (outer) planet is significantly larger (smaller) than that for both planets from numerical orbit integration. Based on the magnitudes and the alternating signs of the contributions from the lower order terms, the expected contributions from the e 4 (and higher order) terms are consistent with their bringing the results into agreement with the actual precession rates. There is no simple analytic expression relating e 1 to e 2 when the eccentricities are large, but as in the small eccentricity limit, the eccentricities are related by the requirement that the periapses precess at the same rate.
Although the analysis in this section is able to explain why the GJ 876 and Io-Europa systems have different resonance configurations, the slow convergence of the series for moderate to large eccentricities means that it has limited usefulness if one is interested in the more general question of what stable configurations are possible for different periapse precession rates and masses. A practical way to investigate this latter question is through numerical migration calculations like those in §4.1, which drive a system through a sequence of configurations with different precession rates. As we shall see, for systems with masses like those in GJ 876, there are stable configurations with θ 1 , θ 2 , and θ 3 librating about 0 • for 0.15 e 1 0.86 (see Fig. 3). We have also found that configurations with resonance variables librating about angles other than 0 • and 180 • are possible when the masses are different from those in GJ 876; these configurations will be discussed in a subsequent paper.

MIGRATION SCENARIO FOR ORI-GIN OF RESONANCES
We now turn to the question of the origin of the two 2:1 mean-motion resonances and the secular resonance in the GJ 876 system. A scenario for the origin of this resonance configuration is that it was assembled by differential migration of orbits that were initially much further apart. We shall see in §4 that libration of all three resonance variables is easily established by any process that drives the orbits toward each other. Although the results presented in §4 are quite general, they will be interpreted in the context of a particular differential migration process, namely via the gravitational interaction between the planets and the nebula from which they form.
It is generally accepted that planets form in a disk of gas and dust that surrounds a young star. Since the two planets around GJ 876 are of order Jupiter mass and are presumably gas giants, they must have formed within the lifetime of the gas disk. If a planet is sufficiently massive, the torques exerted by the planet on the gas disk can open an annular gap in the disk about the planet's orbit. The conditions for gap formation are that the Roche radius of the planet, r R = (M/3M 0 ) 1/3 a, exceeds the scale height H of the disk, or equivalently, and that the viscous condition is satisfied (e.g., Lin & Papaloizou 1993). In the above equations, M and M 0 are the masses of the planet and the star, a is the semi-major axis of the planet's orbit, Ω is the angular Kepler speed at a, and the kinematic viscosity ν is expressed using the Shakura-Sunyaev α prescription: ν = αH 2 Ω. Since the planet-star mass ratios of the planets around GJ 876 are M 1 /M 0 = 2.28 × 10 −3 and M 2 /M 0 = 7.17 × 10 −3 , these planets are expected to open gaps individually during their growth to their final masses if their orbits are sufficiently far apart, unless H/a 0.09 or α 1/2 (H/a) 7.5 × 10 −3 .
The numerical values α = 4 × 10 −3 and H/a = 0.05 are typical of models of protoplanetary disks (e.g., Bryden et al. 2000;Kley 2000;Papaloizou et al. 2001). Hartmann et al. (1998) have also inferred from observed properties of T Tauri disks that α ∼ 10 −2 . The most promising source of an effective viscosity in accretion disks is magnetohydrodynamic (MHD) turbulence initiated and sustained by the magnetorotational instability, which is capable of producing an effective α as large as 0.1 (e.g., Stone et al. 2000). However, protoplanetary disks are probably too weakly ionized for MHD turbulence to develop fully, except at small radii ( 0.1 AU) and possibly in a layer near the surface of the disk at larger radii (Gammie 1996). Another possible source of effective viscosity is damping of density waves excited by many small (of order Earth mass) planets, which is capable of producing an effective α 10 −3 (Goodman & Rafikov 2001). Bryden et al. (2000) and Kley (2000) have performed hydrodynamic simulations of a system consisting of a central star, a gas disk and two planets. The planets are massive enough (M/M 0 ≈ 10 −3 ) to open gaps individually, and the ratio of their initial semi-major axes is a 1 /a 2 ≈ 1/2. They found that the planets clear out nearly all of the nebular material between them in a few hundred orbital periods. Then the torques exerted by the nebular material outside the outer planet drives that planet toward the star, whereas any nebular material left on the inside of the inner planet drives that planet away from the star. The depletion of the inner disk means that the inner planet may not move out very far, but the net effect is always to drive the orbits of two massive planets toward each other. The time scale on which the planets migrate is the disk viscous time scale, whose inverse is (Ward 1997 where a dot over a symbol denotes d/dt. In addition to its effect on the semi-major axis, planet-nebula interaction can also affect the orbital eccentricity of a planet (e.g., Goldreich & Tremaine 1980;Artymowicz 1992Artymowicz , 1993Pa-paloizou et al. 2001). A planet interacts with a disk in the vicinity of Lindblad and corotation resonances. The leading contribution toȧ is due to Lindblad resonances with the l = m Fourier components of the planet's perturbation potential. (The indices l and m for the Fourier series in time and azimuthal angle in this context should not be confused with ℓ and m in eq. [8].) For a planet orbiting in a disk gap, the so-called co-orbital Lindblad and corotation resonances are not important, and the leading contributions toė are damping due to l = m ± 1 corotation resonances and excitation due to l = m − 1 outer and l = m + 1 inner Lindblad resonances. The net effect of the corotation resonance damping and the Lindblad resonance excitation depends on the distribution of the nebular material, which is itself determined by the interaction with the planet. If the gap is not too wide and many resonances of both types are present, previous calculations indicate that there is net eccentricity damping. For example, if we consider an outer disk of constant surface mass density Σ that extends radially from a + ∆ to ∞, with ∆ ≪ a, then integration of equations (30) and (31)  In the case of GJ 876, since the planets clear out nearly all of the nebular material between them and the inner disk is likely to be depleted, the dominant resonant interactions are expected to be those between the outer disk and the outer planet. However, the estimate (14) forė is not adequate because, in addition to the fact that the mass distribution near the inner edge of the outer disk is not modeled properly, the outer planet is sufficiently massive that only a few low-m resonances are likely to be present in the disk. In fact, if the outer planet is able to open a gap out to the 2:1 commensurability, there would be no corotation resonances and only one Lindblad resonance with l = m − 1 = 1 at the 3:1 commensurability, and there would be eccentricity growth instead of damping (Artymowicz 1992). (In addition, the disk can become eccentric and the growth of the planet's orbital eccentricity can be enhanced by the interaction with the eccentric disk, if the planetary mass is comparable to a characteristic mass of the disk; Papaloizou et al. 2001.) From a comparison of the resonant and viscous torques, an approximate condition for opening a gap out to the 2:1 commensurability is (Artymowicz 1992;Lin & Papaloizou 1993) Since the planet-star mass ratio of the outer planet around GJ 876 is M 2 /M 0 = 7.17 × 10 −3 and close to the critical value in equation (16), without knowing the exact values of the disk parameters, it is unclear whether one should expect eccentricity damping or growth. However, as we shall see in §4, significant eccentricity damping with |ė/e| ≫ |ȧ/a| is required to produce the observed eccentricities of the GJ 876 system, unless the migration after resonance capture is severely limited. Therefore, at least the condition (16) must not be satisfied, implying that for the outer disk of the young GJ 876 system. A more detailed analysis using hydrodynamic simulations is necessary to determine whether sufficient eccentricity damping can be produced by such an outer disk.

NUMERICAL RESULTS FOR MI-GRATION SCENARIO
In this section we present the results of a series of numerical orbit integrations designed to determine the conditions under which the dynamical properties of the current GJ 876 system could be produced by any process (such as the planetnebular interaction discussed in §3) that drives the orbits of the two planets toward each other. We consider a system consisting of a central star and two planets, where the stellar and planetary masses are those for the Laughlin-Chambers Keck+Lick solution (Table 1). Unless stated otherwise, the planets are initially on coplanar, circular orbits, with the mean longitudes differing by 180 • and the ratio of the semi-major axes a 1 /a 2 = 1/2, far from the 2:1 mean-motion commensurability.
In addition to the mutual gravitational interactions of the star and the planets, we force the osculating semi-major axis of planet i to migrate at a rateȧ i . In most cases, we assume that only the outer planet is forced to migrate inward, with a migration rate of the formȧ 2 /a 2 = constant. The effects of adopting a more general form of the migration rate (e.g.,ȧ 2 /a 2 being a function of a 2 ) or forcing the inner planet to migrate outward are discussed in §4.2. For the calculations in §4.2, we also damp the osculating eccentricity at a rateė i . We adopt a damping rate of the forṁ e i /e i = −K|ȧ i /a i |, where K is a positive constant. As we shall see, this form of eccentricity damping has the convenient property that the eccentricities reach equilibrium values after capture into the 2:1 resonances. We note that the relationė/e ∝ȧ/a is valid for the simple estimates (14) and (15) for the interaction of a massive planet with an outer disk if ∆/a ≈ constant and that the same relation between the damping and migration rates also holds for a planet that is too small to open a gap in a gas disk with constant H/a and undergoes the socalled type I migration (Artymowicz 1993;Ward 1997).
The numerical orbit integrations were performed using the symplectic integrator SyMBA (Duncan, Levison, & Lee 1998) modified to include the orbital migration and eccentricity damping terms. SyMBA is based on a variant of the Wisdom-Holman (1991) method and employs a multiple time step technique to handle close encounters. (The latter feature is not essential for the integrations presented here.) We confirmed the results by repeating some of the integrations using a Bulirsch-Stoer integrator that has also been modified to include migration and damping. We also checked that both modified integrators give the correct exponential decays in a and e when they are used to integrate the orbit of a single planet withȧ/a = constant andė/e = constant. We describe how the integrators were modified in the Appendix.

Migration without Eccentricity Damping
We consider first inward migration of the outer planet without eccentricity damping. In Figure 3 we show the results of a calculation with a 1 = 0.5 AU and a 2 = 1.0 AU initially andȧ 2 /a 2 = −5 × 10 −5 yr −1 . The migration rate is consistent with that due to planet-nebular interaction at ∼ 1 AU (eq. [13]). Figure 3a shows the time evolution of the semi-major axes and eccentricities, and Figure 3b shows the evolution of the two 2:1 mean-motion resonance variables, θ 1 and θ 2 , in plots of e i sin θ i versus e i cos θ i . Initially, only the outer planet migrates inward at the prescribed rate. When the 2:1 mean-motion commensurability is encountered, both θ 1 and θ 2 (and hence the secular resonance variable θ 3 ) are captured into libration about 0 • . We do not see libration of θ 2 and θ 3 about 180 • , which is expected for small eccentricities (see §2), because the planetary masses are so large that fairly large eccentricities (with e 1 reaching about 0.15) are generated before the 2:1 commensurability is encountered. Because of the forced migration, the centers of libration are actually slightly offset from 0 • . The resulting resonant interaction slows down the migration of the outer planet and forces the inner planet to migrate inward while keeping the ratio of the semimajor axes nearly constant (a 1 /a 2 ≈ 2 −2/3 ); it also causes the eccentricities to increase rapidly. The centers of libration remain near 0 • and the amplitudes of libration remain small as the eccentricities increase.
The sequence of configurations with increasing eccentricities that the system is driven through after resonance capture is in fact a sequence with increasingly less negative periapse precession rates. As we discussed in §2, the forced eccentricities (and the libration centers) for a system with stable simultaneous librations of θ 1 and θ 2 are determined by the requirement that the longitudes of periapse on average precess at the same rate (which in turn is determined byλ 1 − 2λ 2 +̟ i = 0). Although a longer integration indicates that the system does eventually become unstable when e 1 ≈ 0.86, the existence of stably librating configurations with e 1 up to 0.86 is remarkable. In particular, the configurations with e 1 0.71 have prograde periapse precessions. We have integrated the configurations at t = 2×10 4 yr (with e 1 ≈ 0.59 and retrograde periapse precessions) and 5×10 4 yr (with e 1 ≈ 0.81 and prograde periapse precessions) forward with migration turned off and confirmed that the system remains stable as seemingly indefinitely repeating configurations.
As we can see in Figure 3a, without eccentricity damping, the eccentricities exceed the observed values for the GJ 876 planets (dashed lines) when the semi-major axes of the resonantly locked planets have decreased by only ≈ 7% after capture into the resonances. This result is insensitive to the adopted parameters. If the migration rate is not too fast (see below), the evolution of a system with different initial a 2 (but for convenience, the same a 1 /a 2 ) orȧ 2 /a 2 is essentially identical to that shown in Figure 3 if we plot the semi-major axes in units of initial a 2 and time in units of (ȧ 2 /a 2 ) −1 . Therefore, unless by coincidence resonance capture occurs just before migration stops because of, e.g., nebula dispersal, eccentricity damping is necessary to produce the observed eccentricities of the GJ 876 system. From a set of calculations similar to that shown in Figure 3, with initial a 1 /a 2 = 1/2, a 2 = 1, 2, and 4 AU, and differentȧ 2 /a 2 , we find that certain capture of both 2:1 mean-motion resonance Fig. 3.-(a) Time evolution of the semi-major axes and eccentricities and (b) evolution of the mean-motion resonance variables θ 1 = λ 1 − 2λ 2 + ̟ 1 and θ 2 = λ 1 − 2λ 2 + ̟ 2 in plots of e i sin θ i versus e i cos θ i for a calculation where the outer planet is forced to migrate inward withȧ 2 /a 2 = −5 × 10 −5 yr −1 and there is no eccentricity damping. Both θ 1 and θ 2 are captured into small amplitude libration about 0 • , but the eccentricities exceed the observed values for the GJ 876 planets (dashed lines) shortly after resonance capture, when the semi-major axes of the resonantly locked planets decrease by only ≈ 7%.
variables requires where a 2 is the semi-major axis of the outer planet when the 2:1 commensurability is encountered. For migration rate within a factor of a few of, but below, the above limit, both resonance variables are captured into libration, but the centers of libration can be significantly different from 0 • and the amplitudes of libration can be large. The condition (18) is satisfied by a migration rate due to planet-nebular interaction with the nominal parameter values in equation (13) by almost 3 orders of magnitude.
To examine the effects of orbital eccentricities on capture into the 2:1 resonances, we performed a series of calculations with non-zero initial eccentricities. Four types of initial conditions were considered: the outer planet was initially at either periapse or apoapse, and either the initial e 1 or e 2 was fixed at 0.01 and the other initial eccentricity was varied. The inner planet was always started at periapse, with its mean longitude differing from that of the outer planet by 180 • , and the remaining parameters were identical to those used in the calculation shown in Figure 3. Since gravitational interaction between the planets causes the eccentricities to fluctuate even when a 1 and a 2 are close to their initial values, the eccentricities quoted below are the maximum eccentricities when a 1 and a 2 are close to their initial values and not the initial eccentricities. We find that certain capture of both 2:1 mean-motion resonance variables requires e 1 0.06 and e 2 0.03.
For eccentricities above these limits, there is nonzero probability for the planets to be captured into higher order resonances (e.g., 5:2) encountered before the 2:1 commensurability. Figure 4 shows the results of a calculation similar to that shown in Figure 3, but with eccentricity damping of the formė 2 /e 2 = −K|ȧ 2 /a 2 |, where K = 100. The capture of the resonance variables θ 1 and θ 2 into libration about 0 • and the  Fig. 3, but for a calculation with eccentricity damping of the formė 2 /e 2 = −K|ȧ 2 /a 2 |, where K = 100. After resonance capture, in addition to librations of θ 1 and θ 2 about 0 • , the eccentricities reach equilibrium values close to the observed values for the GJ 876 planets (dashed lines) and remain constant for arbitrarily long migration in the resonances. (The jagged nature of the plot (b) prior to equilibrium is due to sparse sampling.) initial evolution after resonance capture are similar to the case without damping. However, the eccentricity growth eventually terminates when the damping balances the excitation due to resonant interaction between the planets. The eccentricities reach equilibrium values that remain constant for arbitrarily long migration in the resonances. With K = 100, the equilibrium eccentricities are close to the observed eccentricities of the GJ 876 system. At the end of the calculation shown in Figure 4, when t = 4.6×10 4 yr, the semi-major axes are also similar to those of the GJ 876 planets.

Migration with Eccentricity Damping
We have integrated the configuration at the end of the calculation shown in Figure 4 forward with migration and damping turned off. The system remains stable, with almost no change in the amplitudes of libration (compare Figs. 4b and 5), which are somewhat smaller than those of the GJ 876 system (compare Figs. 1 and 5). Although the Fig. 5.-Continued small amplitude librations of the mean-motion resonance variables θ 1 and θ 2 after termination of planet migration and eccentricity damping. The configuration at the end of the calculation shown in Fig. 4 is integrated forward for 500 yr with migration and damping turned off.
libration amplitudes are fairly similar in the two Laughlin-Chambers best-fit solutions, we may find in the future, as more data is acquired, that improved best-fit solutions will have smaller libration amplitudes. Alternatively, larger amplitudes could be generated after termination of migration and damping by encounters with remaining planetesimals.
In Figure 6 we show the time evolutions of the eccentricities for a set of calculations withȧ 2 /a 2 = −5 × 10 −5 yr −1 ,ė 2 /e 2 = −K|ȧ 2 /a 2 |, and different K. The equilibrium eccentricities decrease with increasing K and are significantly different from the observed eccentricities of the GJ 876 system if K is more than a factor 2-3 larger or smaller than K = 100.
As long asė 2 /e 2 = −K|ȧ 2 /a 2 |, where K is a positive constant, the equilibrium eccentricities are determined by the value of K only and are insensitive to either the magnitude or functional form ofȧ 2 /a 2 . This is demonstrated in Figure 7, where we plot the results of a calculation with the Fig. 6.-Decrease of the equilibrium eccentricities with increasing eccentricity damping rate. Time evolutions of the eccentricities for a set of calculations withȧ 2 /a 2 = −5 × 10 −5 yr −1 ,ė 2 /e 2 = −K|ȧ 2 /a 2 |, and different K are shown.
same K(= 100) as that shown in Figure 4 but with a different form of the migration rate:ȧ 2 = constant orȧ 2 /a 2 ∝ 1/a 2 . The equilibrium eccentricities are identical in Figures 4 and 7. Ifė 2 /e 2 is not exactly proportional to |ȧ 2 /a 2 |, the eccentricities would decrease slowly (after reaching maxima) or increase slowly as the resonantly locked planets migrate, but the ratio of damping to migration rates, |(ė 2 /e 2 )/(ȧ 2 /a 2 )|, just before migration and damping stop should be close to 100 for the final eccentricities to be close to the observed values of the GJ 876 system. Thus far we have considered forced inward migration of the outer planet only. This is the most likely situation in the planet-nebular interaction scenario discussed in §3, since the inner disk is likely to be depleted and the dominant interactions are expected to be those between the outer disk and the outer planet. However, if the inner disk is not depleted, at least initially, or migration and damping are due to another process, the inner planet may also be forced to migrate outward. In Figure 8 we show the results of a calculation where, for simplicity, the inner planet is forced Fig. 7.-An example demonstrating that the equilibrium eccentricities are determined by the ratio of damping to migration rates, K, and do not depend on the functional form of the migration rate. The results of a calculation with the same ratio (K = 100) as that shown in Fig. 4 but with a different form of the migration rate (ȧ 2 = constant) are shown.
to migrate outward at the same rate that the outer planet is forced to migrate inward (ȧ 1 /a 1 = −ȧ 2 /a 2 = 5 × 10 −5 yr −1 ) andė i /e i = −K|ȧ i /a i |. Initially, the inner (outer) planet migrates outward (inward) at the prescribed rate. After resonance capture, resonant interaction overcomes the forced outward migration of the inner planet, and both planets migrate inward slowly. The equilibrium eccentricities are close to the observed eccentricities of the GJ 876 system when K = 10. Therefore, even with eccentricity damping of both planets, significant eccentricity damping with |ė i /e i | ≫ |ȧ i /a i | is required to produce the observed orbital eccentricities of the GJ 876 planets.

Eccentricity Damping by Tidal Dissipation in the Star and Planets
We have found in §4 that significant eccentricity damping with |ė i /e i | ≫ |ȧ i /a i | is required to pro- Fig. 8.-Time evolution of the semi-major axes and eccentricities for a calculation where the inner planet is also forced to migrate outward, withȧ 1 /a 1 = −ȧ 2 /a 2 = 5 × 10 −5 yr −1 ,ė i /e i = −K|ȧ i /a i |, and K = 10. Even with eccentricity damping of both planets, K ≈ 10 is required to produce the observed eccentricities of the GJ 876 planets (dashed lines).
duce the observed eccentricities of the GJ 876 system if the migration has been at all extensive after resonance capture. As we discussed in §3, it is as yet unclear whether sufficient damping could be produced by planet-nebula interaction, even if the condition (17) is satisfied and eccentricity damping is expected (see also §5.2). In this subsection we show that alternative eccentricity damping by tidal dissipation within the star and planets during planet migration is completely negligible.
The star GJ 876 would most likely still be in its pre-main-sequence contracting phase during disk evolution and planet migration. To estimate the rate of circularization of an orbit due to tidal dissipation in the star, we use the stellar radius when the time t since initial contraction is approximately the migration time scale of |a/ȧ| = 2×10 4 yr (eq. [13]) and keep the planets at their current distances from the star. The planets would have migrated for longer than the migration time scale (see, e.g., Fig. 4) and it would take time for the planets to form. But by adopting this minimum time to yield the maximum probable stellar radius and by keeping the planets at their current distances, we maximize our estimate of the tidal rate of circularization of an orbit.
Since the mass of GJ 876 is only 0.32M ⊙ , we assume that GJ 876 follows the nearly vertical Hayashi track in the HR diagram and remains fully convective during its contracting phase. Hence the star remains a polytrope of index n = 1.5, with effective temperature T e ≈ 3500 K (e.g., D'Antona & Mazzitelli 1994). From the virial theorem, the stellar luminosity where σ is the Stefan-Boltzmann constant, G is the gravitational constant, M 0 and R 0 are the stellar mass and radius, and the form of the gravitational energy for a polytrope of index n is given by Chandrasekhar (1939). Integration of equation (20), with 3/2(5 − n) = 3/7 for a polytrope of index n = 1.5, gives R 0 ≈ 3.85 × 10 11 cm after t = 2 × 10 4 yr, leading to a ratio of the stellar radius to the orbital semi-major axis R 0 /a ≈ 0.197 for the inner planet and R 0 /a ≈ 0.125 for the outer planet.
Since the present orbital period of the inner planet of the GJ 876 system is about 30 days, which is long compared to the periods of free oscillation of the star, we can use the rate for circularizing an orbit corresponding to dissipation dominated by the equilibrium tide for the star, which is given by (Zahn 1989) where q = M/M 0 is the planet-star mass ratio and t f = (M 0 R 2 0 /L) 1/3 . We neglect the reduction of the turbulent viscosity when the tidal period is short compared to the turnover time of the largest convective eddies and adopt the maximum λ circ ≈ 0.048 for a fully convective star. Substitution of R 0 /a from the previous paragraph, along with q for the Laughlin-Chambers Keck+Lick solution (Table 1), t f = (M 0 /4πσT 4 e ) 1/3 ≈ 0.574 yr, and λ circ ≈ 0.048, into equation (21) gives 1 t circ ≈ 9.1 × 10 −9 yr −1 for inner planet, ≈ 7.6 × 10 −10 yr −1 for outer planet.
Since we have assumed that the planets are at their closest proximity to the star for the whole time the tidal effects are damping the eccentricities and that the star is inflated to a maximum size, these values are extreme upper bounds on the rate of eccentricity damping by tides raised on the star. Even in its inflated state, the dissipation in the star can have essentially no effect on the orbital eccentricities. The theoretical circularization rate in equation (21) has been tested by comparing with the observed circularization rates of binaries. It is in good agreement with the observed rates for binaries containing giant stars (Verbunt & Phinney 1995), but even with the maximum λ circ , it is about 50-100 times slower than the observed rates for binaries containing main-sequence solar-type stars, which have radiative cores (Claret & Cunha 1997;Goodman & Oh 1997). It is not yet clear why there is a discrepancy in the latter case or that this discrepancy is relevant for a fully convective pre-main-sequence star. Nevertheless, even if we increase the circularization rates in equation (22) by a factor 100, they are still much smaller than the migration rate due to planet-nebula interaction.
Tidal dissipation within a gaseous planet damps its orbital eccentricity at a rate given by (e.g., Peale, Cassen, & Reynolds 1980) where R, k 2 , and Q are the radius, the potential Love number, and the dissipation function of the planet. If we adopt values for R, k 2 , and Q similar to those of Jupiter, with R = 7 × 10 4 km, k 2 = 0.38 (Gavrilov & Zharkov 1977), and Q = 5 × 10 4 (which is approximately the lower bound on Q for Jupiter ;Yoder & Peale 1981), |ė/e| ≈ 1.5 × 10 −12 yr −1 for the inner planet of the GJ 876 system. The rate for damping the outer planet's orbital eccentricity from dissipation within itself is of course even smaller. The planets would most likely be larger and contracting during their migration. However, unless the planetary radii are unrealistically large and comparable to the Roche radii, the eccentricity damping rate due to tidal dissipation within the planets is smaller than the migration rate due to planet-nebula interaction by several orders of magnitude.

Other Studies
After completing our calculations (Lee & Peale 2001), two papers with complementary calculations came to our attention (Snellgrove, Papaloizou, & Nelson 2001;Murray, Paskowitz, & Holman 2001). Snellgrove et al. (2001) also find from numerical orbit integrations with forced migration and eccentricity damping of the outer planet that the orbital eccentricities of the GJ 876 planets require a short time scale for eccentricity damping compared to the migration time scale. Note, however, that they adopt the minimum planetary masses from the two-Kepler fit by Marcy et al. (2001) and, not surprisingly, have difficulties matching both e 1 and e 2 from their calculations to those from the same fit, since the eccentricities from the two-Kepler fit are not consistent with smallamplitude simultaneous librations of θ 1 and θ 2 . Snellgrove et al. also develop an analytic resonance theory that is first order in the eccentricities. As we showed in §2, a first order theory is inadequate for understanding the current GJ 876 resonance configuration, which has θ 1 , θ 2 , and θ 3 all librating about 0 • . Furthermore, as we found in §4 (see, e.g., Fig. 3), because the eccentricities are excited to large values in the GJ 876 evolution before resonance capture, there is no time during the evolution when a first order theory can be useful for this system. Thus all predictions coming from this theory should be viewed with caution. Snellgrove et al. also present a hydrodynamic simulation of the planet-nebula interaction, where both planets are inside a cavity with almost no disk material. At the end of this simulation, e 1 ≈ 0.34 and the semi-major axes of the planets have decreased by about 13% since resonance capture. We have performed a numerical orbit integration without eccentricity damping similar to that shown in Figure 3 but with the planetary masses adopted by Snellgrove et al. and find that e 1 ≈ 0.38 for the same reduction in semi-major axes. Thus it is not clear that the eccentricities have reached equilibrium values at the end of this simulation; even if they have, these equilibrium values are too large for the GJ 876 system. It appears that the disk model used in this hydrodynamic simulation is not very effective in damping the eccentricities, which emphasizes the uncertainty in such damping. Murray et al. (2001) consider mainly an alternative scenario in which the migration and eccentricity damping of the outer planet are due to scattering of planetesimals in the disk population. Their numerical simulations confirm our results (and those of Snellgrove et al. 2001) of easy capture of an inner planet into resonance, the dual migration of both planets in the resonance, and the growth of the eccentricities, although they limit their numerical and analytical studies to an outer planet of Jupiter mass and inner planets of several Earth masses and devote much of the discussion to resonances of higher order than those at the 2:1 mean-motion commensurability. Thus their results are not directly applicable to the GJ 876 system.

CONCLUSIONS
Radial velocity measurements by Marcy et al. (2001) have revealed two planets in resonant orbits about the star GJ 876. The remarkable orbital fit obtained by Laughlin & Chambers (2001), which finds both lowest order, eccentricity-type mean-motion resonance variables at the 2:1 commensurability librating with small amplitudes, means that the resonances are almost certainly real and indefinitely stable. The existence of the eccentricity-type resonances implies that the assumed coplanarity of the orbits is probably close to reality.
The GJ 876 planetary system has revealed properties of the 2:1 orbital resonances that have not been observed nor analyzed before. The libration of both lowest order mean-motion resonance variables, θ 1 = λ 1 − 2λ 2 + ̟ 1 and θ 2 = λ 1 − 2λ 2 + ̟ 2 , and the secular resonance variable, θ 3 = ̟ 1 − ̟ 2 , about 0 • was not anticipated, as the familiar Io-Europa 2:1 resonance has θ 1 librating about 0 • , but θ 2 and θ 3 librating about 180 • -a configuration that would persist in the absence of Ganymede. Thus conjunctions for the Jovian satellites occur when Io is near periapse and Europa is near apoapse, whereas conjunctions of the two planets about GJ 876 occur when both planets are near periapse. We understood this to be mainly a function of the eccentricities of the orbits, where the resonance configuration with θ 1 ≈ 0 • and θ 2 ≈ θ 3 ≈ 180 • must obtain when the eccentricities are small, but the resonance configuration with θ 1 ≈ θ 2 ≈ θ 3 ≈ 0 • prevails for a system with masses like those in GJ 876 when the eccentricities are large. A necessary condition for stable simultaneous librations of both mean-motion resonance variables is that̟ 1 =̟ 2 on average, so that the relative alignment of the lines of apsides of the two orbits is maintained. The periapse precessions are dominated by resonant terms in the disturbing potential whose arguments are θ 1 , θ 2 , θ 3 , and their linear combinations. The dominance of the lowest order terms when the eccentricities are small allows equal precession rates only if the lines of apsides are anti-aligned (θ 1 ≈ 0 • , θ 2 ≈ θ 3 ≈ 180 • ), whereas dominance of higher order terms when the eccentricities are large results in equal precession rates for a system with masses like those in GJ 876 when θ 1 ≈ θ 2 ≈ θ 3 ≈ 0 • . The equality of the precession rates also determines a relationship between the eccentricities of the two orbits, although there is no simple analytic expression for this relationship when the eccentricities are large, since the first order theory is not a good representation. The stability of the GJ 876 resonance configuration for values of e 1 up to 0.86 was also a surprise.
Any process that drives the two originally widely separated orbits toward each other can result in capture of the planets into orbital resonances at the 2:1 commensurability. The naturally occurring situation where nebular disk material is cleared between two planets sufficiently massive to individually open gaps in the disk (Bryden et al. 2000;Kley 2000) leads to inward migration of the outer planet and possibly outward migration of the inner planet. Thus the likely origin of the resonances in the GJ 876 system is this differential planet migration due to torques induced by the planet-nebula interaction. We have shown that forced inward migration of the outer planet of the GJ 876 system results in certain capture of θ 1 , θ 2 and hence θ 3 into libration if initially e 1 0.06 and e 2 0.03 and |ȧ 2 /a 2 | 3 × 10 −2 (a 2 / AU) −3/2 yr −1 . The lat-ter rate is three orders of magnitude higher than the likely rate of ∼ 5 × 10 −5 (a 2 / AU) −3/2 yr −1 due to planet-nebular interaction. The bounds on the eccentricities result not so much from the transition from certain to probabilistic capture at the 2:1 resonances but from likely capture into higher order resonances such as 5:2 before the 2:1 commensurability is encountered. Continued migration of the planets while locked in the 2:1 resonances leads to rapid growth in the orbital eccentricities that exceed the observed eccentricities of the GJ 876 system after only a further decrease in the semi-major axes of about 7% if there is no eccentricity damping. So unless resonance capture occurred near the end of migration, the observed values of the eccentricities require eccentricity damping. With damping of the formė i /e i = −K|ȧ i /a i |, where K is a positive constant, eccentricity growth is terminated at values of the eccentricities that increase with decreasing K, and the eccentricities remain constant for indefinite duration of the migration. The observed eccentricities result for K ≈ 100 if there is forced migration and eccentricity damping of the outer planet only, but for K ≈ 10 if there is also forced migration and eccentricity damping of the inner planet. This result is independent of the magnitude or functional form ofȧ i /a i as long aṡ e i /e i = −K|ȧ i /a i | is preserved. Relaxing the last condition leads to a slow drift in the eccentricities during migration and would require the migration to terminate as the eccentricities pass through the observed values.
Existing analytic estimates of the effects of planet-nebular interaction are consistent with eccentricity damping of the formė i /e i = −K|ȧ i /a i |, if the planet-star mass ratio is not too large (e.g., Goldreich & Tremaine 1980;Artymowicz 1992Artymowicz , 1993Lin & Papaloizou 1993;Ward 1997). However, the planet-star mass ratio of the outer planet of the GJ 876 system is sufficiently close to the critical value separating eccentricity growth from damping for nominal values of the disk parameters that it is uncertain whether such damping would occur. Even if the disk parameters are such that eccentricity damping would occur, it is not clear that the magnitude of K would be sufficiently large to constrain the eccentricities in the GJ 876 system to the observed values. We have shown that the alternative eccentricity damp-ing mechanism involving the dissipation of tidal energy within the star and the planets is completely negligible. Further long-term hydrodynamic simulations with different physical assumptions and parameters are required to determine whether planet-nebular interaction could produce sufficient eccentricity damping to allow arbitrary migration of the planets within the resonances in the young GJ 876 system while preserving eccentricities comparable to those observed. If not, the migration must have been finely tuned to stop when the system had progressed to its observed state, although this latter constraint is too ad hoc to be believable.
It is a pleasure to thank Greg Laughlin for furnishing the details of the dynamical fit by Laughlin & Chambers before publication and Lars Bildsten for pointing out the enhancement of tidal dissipation in a pre-main-sequence star. We also thank D.N.C. Lin, J.J. Lissauer, and W.R. Ward for informative discussions. This research was supported in part by NASA grants PGG NAG5-3646 and OSS NAG5-7177.

A. NUMERICAL METHODS
In this appendix we describe how the symplectic integrator SyMBA and the Bulirsch-Stoer integrator used for the numerical orbit integrations presented in §4 were modified to include the forced orbital migration and eccentricity damping terms.
A second-order symplectic integrator for a Hamiltonian of the form H = H 0 + H 1 , where H 0 and H 1 are separately integrable, can be represented as The three operators in equation (A1) represent a single step of an algorithm that starts with evolving the system under the influence of H 0 only for half a time step τ /2, then evolving it for a full time step τ under the influence of H 1 , and then evolving it for another half a time step τ /2 under H 0 . For example, in the Wisdom-Holman (1991) method for the gravitational N -body problem, the Hamiltonian in Jacobi coordinates is divided into H 0 that describes the Keplerian motion of the planets around a central star and H 1 that describes the perturbation of the planets on one another and is a function of the positions q i only. A step of the Wisdom-Holman method is thus: (1) Each planet evolves along a Kepler orbit for time τ /2; (2) each planet receives a kick to its momentum of the amount −τ ∂H 1 /∂q i while q i is unchanged (since H 1 does not involve the canonical momenta); (3) each planet evolves along a Kepler orbit for time τ /2, starting with the new momentum after the kick. Recursive application of the basic algorithm (A1) allows one to construct symplectic integrators for Hamiltonians that consist of more than two integrable parts. For example, a single step of an algorithm for a Hamiltonian of the form The symplectic integrator SyMBA (Duncan et al. 1998) is based on a variant of the Wisdom-Holman method, with the gravitational N -body Hamiltonian written in terms of positions relative to a central star and barycentric momenta, and employs a multiple time step technique to handle close encounters. In the SyMBA algorithm, the Hamiltonian is divided into more than two parts and the recursive application of the algorithm (A1) discussed in the previous paragraph is utilized. Although the additional forced orbital migration and eccentricity damping terms are not Hamiltonian, they can be incorporated in an analogous manner. Our modified algorithm is where E grav (τ ) denotes a complete time step for the conservative gravitational N -body problem using the SyMBA algorithm, and E a (τ /2) and E e (τ /2) denote changing the canonical variables according to the imposedȧ i andė i terms, respectively, for time τ /2. During the application of theȧ i term, all of the other orbital elements are constant, and a i,1 = a i,0 exp(τȧ i /2a i ), where a i,0 and a i,1 are a i at the beginning and end of the step, ifȧ i /a i = constant (this can be easily generalized for, e.g.,ȧ i /a i ∝ a γ i ). Note that we do not use truncated approximation such as a i,1 = a i,0 (1 + τȧ i /2a i ). Similarly, the E e (τ /2) step changes the eccentricities according to e i,1 = e i,0 exp(−τ K|ȧ i /a i |/2) ifė i /e i = −K|ȧ i /a i |. By modifying the algorithm in a symmetric manner and using exact solutions in the E a (τ /2) and E e (τ /2) parts, there should be little (if any) secular growth in the energy error (e.g., Mikkola 1998).
For the Bulirsch-Stoer integrator, additional terms must be included in the equations of motion in Cartesian coordinates to account for the forced orbital migration and eccentricity damping. In the following, we simplify the notation by considering a specific planet and dropping the subscript. Let (x, y, z) be the Cartesian coordinates of the planet with respect to the star and r be the distance of the planet from the star. The osculating orbital elements a, e, i, f , ω, and Ω are the semi-major axis, eccentricity, inclination, true anomaly, argument of periapse, and longitude of the ascending node on the xy plane, respectively. The additional terms in the equations of motion due to the forced migrationȧ and eccentricity dampingė terms with similar expressions for the other coordinates. To evaluate the partial derivatives in equations (A3) and (A4) and similar expressions for the other coordinates, we need to express the position and velocity in terms of the osculating orbital elements: x = r cos Ω cos (ω + f ) − r cos i sin Ω sin (ω + f ), y = r sin Ω cos (ω + f ) + r cos i cos Ω sin (ω + f ), where r,ṙ, and rḟ are in terms of a, e, and f (e.g., Murray & Dermott 1999 where equation (A7) should be used to get the functional forms. Unfortunately, the orbital elements must be calculated at each call to the differential equations when there is eccentricity damping. a Two-Kepler fit by Marcy et al. 2001 and dynamical fit using a Levenberg-Marquardt N -body integration scheme by Laughlin & Chambers 2001. b The parameters are the planetary mass M in terms of Jupiter mass M J , the period P , the semi-major axis a, the orbital eccentricity e, the longitude of the periapse ̟, the epoch T , and the mean anomaly M. The stellar mass is 0.32M ⊙ . , with |k| + |ℓ| + 2m + 2n = N , in the disturbing potential (eq. [8]) are grouped together under order e N .