W. Kley1 - M. H. Lee2 - N. Murray3 - S. J. Peale2
1 - Institut für Astronomie & Astrophysik, Abt. Computational Physics, Universität Tübingen, Auf der Morgenstelle 10, 72076 Tübingen, Germany
2 - Department of Physics, University of California, Santa Barbara, CA 93106, USA
3 - Canada Research Chair in Astrophysics; CITA, University of Toronto, 60 St. George Street Toronto, Ontario M5S 3H8, Canada
Received 7 January 2005 / Accepted 24 March 2005
The two planets about the star GJ 876 appear to have undergone extensive migration from their point of origin in the protoplanetary disk - both because of their close proximity to the star (30 and 60 day orbital periods) and because of their occupying three stable orbital resonances at the 2:1 mean-motion commensurability. The resonances were most likely established by converging differential migration of the planets leading to capture into the resonances. A problem with this scenario is that continued migration of the system while it is trapped in the resonances leads to orbital eccentricities that rapidly exceed the observational upper limits of and . As seen in forced 3-body simulations, these lower eccentricities would persist during migration only for an eccentricity damping rate exceeding . Previous theoretical and numerical analyses have found or even eccentricity growth through disk-planet interactions.
In an attempt to find effects that could relax the excessive eccentricity damping requirement, we explore the evolution of the GJ 876 system using two-dimensional hydrodynamical simulations that include viscous heating and radiative cooling in some cases. Before we evolve the whole system, the disk with just the outer planet embedded is brought into equilibrium. We find that the relaxed disk remains circular in all models for low planet-to-star mass ratios q2, but becomes eccentric for high mass ratios for those models with fixed temperature structure. The disk in models with full radiative thermodynamics remains circular for all q2 considered due to the larger disk temperatures. Given the small stellar mass, the mass ratio for the GJ 876 system is rather high (with minimum ), and so the GJ 876 disk may have been slightly eccentric during the migration.
With a range of parameter values, we find that a hydrodynamic evolution within the resonance, where only the outer planet interacts with the disk, always rapidly leads to large values of eccentricities that exceed those observed - similar to the three-body results. The resonance corresponding to the resonant angle (involving the inner planet's periapse longitude, ) is always captured first. There is no additional delay in capturing into resonance that is attributable to the secular prograde contribution to the precession of from the interaction with the disk, but an eccentric disk can induce a large outer planet eccentricity e2 before capture and thereby further delay capture of for larger planetary masses. The delay in capturing into libration, while delaying the resonance-induced growth of e2, has no effect on the forced eccentricities of both planets, which are uniquely determined by the resonance conditions, once both are librating.
Only if mass is removed from the disk on a time scale of the order of the migration time scale (before there has been extensive migration after capture), as might occur for photoevaporation in the late phases of planet formation, can we end up with eccentricities that are consistent with the observations.
Key words: accretion, accretion disks - planets and satellites: formation - hydrodynamics - methods: N-body simulations
Both theoretical and numerical analyses verify the ubiquity of planet migration due to interaction with the protoplanetary disk of gas and dust (e.g. Nelson et al. 2000; Ward 1997). The clearing of disk material between two planets both capable of opening a gap in the disk leads to differential migration of the two planets as at least the outer planet is forced in by the material remaining outside its orbit. The convergence of the orbits naturally allows capture into stable resonant orbital configurations. Hydrodynamical simulations of two embedded planets capturing each other into resonance have been performed by several groups (Kley et al. 2004; Snellgrove et al. 2001; Kley 2000; Papaloizou 2003; Bryden et al. 2000). The resonance capture has also been analyzed with extensive three-body calculations, where migration is simply imposed with either the semi-major axis migration rate and the eccentricity damping rate specified explicitly or ad hoc forces added to produce migration and eccentricity damping (Kley et al. 2004; Snellgrove et al. 2001; Lee & Peale 2002; Nelson & Papaloizou 2002; Lee 2004).
Depending on the migration rates, masses, and initial orbital separations and eccentricities of the two planets, capture can occur in different resonances. For planets as massive as those in GJ 876, capture into the 2:1 resonances is robust if the initial and the initial eccentricities are small, and it is most likely that the 2:1 resonances for GJ 876 were established by converging differential migration (Lee & Peale 2002). The sequence of 2:1 resonance configurations that a system with initially nearly circular orbits is driven through by continued migration depends mainly on the planetary mass ratio M1/M2. For as in the GJ 876 system, a system is first captured into antisymmetric configurations with librating about and (and hence ) librating about . Continued migration forces e1 to larger values and e2 from increasing to decreasing until when Then the system converts to symmetric configurations like that of GJ 876, with both and librating about (Lee 2004; see also Fig. 16 below). There are other types of 2:1 resonance configurations (with both and librating about or asymmetric librations of and about values far from either or ) for , but they are either unstable for planets as massive as those in GJ 876 or not reachable by convergent migration of planets with nearly constant masses and coplanar orbits (Lee 2004). Asymmetric libration configurations can result from convergent migration for 2:1 resonances with and for 3:1 resonances for a wider range of M1/M2 (Ferraz-Mello et al. 2003; Beaugé et al. 2003; Kley et al. 2004; Lee 2004).
While it is easy to understand the symmetric resonance configuration of GJ 876 from convergent migration, the small eccentricities of this system are a puzzle. The continued migration after capture into resonance drives the pair of planets deeper into the resonances ( 2 n2 - n1 < 0 increasing, with being the mean orbital motions). The identical mean retrograde precessions of the must therefore decrease in magnitude to keep near zero, i.e., to maintain the resonance configuration. This decrease is effected by increasing the orbital eccentricities when the system is in the configuration with both and librating about . The three-body calculations have shown that the eccentricity growth is rapid if there is no eccentricity damping, with the average eccentricities ( , ) of the fit of Laughlin & Chambers (2001) exceeded after only a 7% decrease in the orbital radii after capture (Lee & Peale 2002). The eccentricities can be maintained at the observed low values as migration continues if there is significant eccentricity damping of (Lee & Peale 2002). In contrast the full hydrodynamical calculations typically give similar timescales for both migration and eccentricity damping (Kley et al. 2004). The low eccentricities become even more problematic if disk-planet interactions drive eccentricity growth. Ogilvie & Lubow (2003) and Goldreich & Sari (2003) have suggested that the eccentricity-damping corotation torques can be reduced sufficiently to trigger eccentricity growth if the eccentricity is above a critical value, and the eccentricity of the outer planet in GJ 876 is well above their estimates for this critical value.
Hydrodynamical simulations of two planets interacting with their protoplanetary disk are computationally expensive, and the number of such simulations in the literature is relatively small (Kley et al. 2004; Snellgrove et al. 2001; Kley 2000; Papaloizou 2003; Bryden et al. 2000). In addition, often the masses used in these simulations are not appropriate for the GJ 876 system or the simulations did not continue for a very long time. While there has been more extensive three-body calculations with imposed migration and eccentricity damping (Kley et al. 2004; Snellgrove et al. 2001; Lee & Peale 2002; Nelson & Papaloizou 2002; Lee 2004), they do not model disk-planet interactions self-consistently and, in particular, the effects of the apsidal precession induced by the disk have not been considered.
In this paper we present a more comprehensive set of numerical calculations treating the evolution of two planets interacting with their protoplanetary disk, with special emphasis on the resonant system GJ 876. For this purpose we solve the full hydrodynamical equations, including viscous heating and radiative cooling effects in some cases, and follow the joint planet-disk evolution over several thousand orbits of the planets. We find that we cannot reproduce the observed low eccentricities using these straightforward assumptions. Taking into account mass loss from the disk, as might be the case in the late dissipation stages of protoplanetary disks subject to, e.g., photoevaporation, we are able to reduce the final values of the eccentricities close to the observed values.
In Sect. 2 we give an overview of the observational data for GJ 876. In Sect. 3 we explain our physical and numerical approach, and in Sect. 4 we present the results of our simulations. This is followed in Sect. 5 with three-body calculations and analytic theory to interpret the results of Sect. 4. We state our conclusions in Sect. 6.
Table 1: Orbital parameters of the two planets of the planetary system GJ 876 at epoch JD 2 449 679.6316 assuming co-planarity and as given by Laughlin et al. (2004). The adopted stellar mass is . P denotes the orbital period, M the mass of the planet, a the semi-major axis, e the eccentricity, the angle of periastron at epoch, and q the mass ratio M/M*.
Although a dynamical fit can in principle yield the orbital
inclinations and masses of the planets, the
value of the fit
for GJ 876 is relatively insensitive to the inclination i of the
assumed coplanar orbits for
rapidly only for
(Laughlin et al. 2004).
In Fig. 1 we show the planetary mass ratio M1/M2and the average orbital eccentricities
for the best-fits found by
Laughlin et al. (2004) for
The deviation from Keplerian motion that is most strongly constrained
by the available observations of GJ 876 is the resonance induced
retrograde apsidal precession of the orbits at an average rate of
and indeed this precession
has now been observed for more than one full period.
As i decreases, the individual planetary mass Mj increases
(and M1/M2 is nearly constant and 0.31). The increase in planetary mass acts to increase
This increase can be offset by an increase in
which acts to decrease
observed value without any significant change
To allow for the mass uncertainties, we shall consider systems with
different planetary masses in our numerical models.
|Figure 1: Planetary mass ratio, M1/M2, average orbital eccentricities, and , and the ratio of eccentricity damping rate to migration rate, K, for the GJ 876 best-fit solutions of Laughlin et al. (2004) with coplanar inclination . The value of is for the equilibrium eccentricities to match if there is inward migration and damping of the outer planet only.|
|Open with DEXTER|
The most striking feature of GJ 876 is the exact 2:1 orbital resonance of the planets. The angles , , and describe the state of the system for coplanar orbits, and all three angles librate about in the GJ 876 system. The libration amplitudes are , , and for the fit shown in Table 1, and the smallest libration amplitudes in the range are at and and at (Laughlin et al. 2004). The small libration amplitudes indicate that the system is deep in the resonance, which occurs naturally in the differential migration scheme for forming the resonant structure if the planets approach the resonance slowly with initially small eccentricities (Lee & Peale 2002).
Lee & Peale (2002) have shown using three-body integrations with imposed migration and eccentricity damping that the eccentricities reach equilibrium values that remain constant for arbitrarily long migration within the resonances if . For inward migration and damping of the outer planet only, is required to produce equilibrium eccentricities that match the eccentricities of the (or ) fit of Laughlin & Chambers (2001). We have repeated the calculations in Lee & Peale (2002) for the updated fits in Laughlin et al. (2004), and the value of K as a function of iis shown in Fig. 1. For , , which is close to the value inferred by Lee & Peale (2002). Benedict et al. (2002) have claimed from HST astrometric measurements that the inclination is very close to , which would require a significantly larger amount of eccentricity damping ( ). The increase in with decreasing i leads to a decrease in K, but K must still be greater than about 40 if .
We use cylindrical coordinates ( ) and consider a vertically averaged, infinitesimally thin disk located at z=0, where the origin of the coordinate system is at the position of the star. The basic hydrodynamic equations (mass and momentum conservation) describing the time evolution of such a two-dimensional disk with embedded planets have been stated frequently and are not repeated here (see Kley 1999). Additional information on the treatment of embedded planets is given in Kley et al. (2004).
We also present radiative
models with an improved thermodynamic treatment using
the thermal energy equation
In addition to the simplified treatment of radiation as expressed in Eq. (1)
we also have considered models including radiative diffusion in the
In this case a term
|Figure 2: The global grid structure with gray-shaded initial surface density superimposed. Every second grid points is shown. The dots denote the location of the star and the two planets and the oval line refers to the Roche lobe of the outer planet.|
|Open with DEXTER|
We work in a rotating reference system, rotating approximately with the initial period of the outer planet. As the coordinate system is accelerated and rotating, we take care to include the indirect terms.
The initial hydrodynamic structure of the disk (density, temperature, velocity) is axisymmetric with respect to the location of the star, and the surface density scales as , with a superimposed initial gap (Kley 2000). The initial velocity is pure Keplerian rotation ( ), and the initial temperature stratification is given by which follows from an assumed constant H/r. For the isothermal models the initial temperature profile is left unchanged, while for the radiative cases it is evolved according to Eq. (1).
The kinematic viscosity
is parameterized by an
where the isothermal sound-speed is given by
and H(r) is either held fixed (for the isothermal models)
or, for the radiative models, calculated from the temperature/sound-speed as
At the inner radial boundary outflow conditions are applied; matter may flow out, but none is allowed to enter. This procedure mimics the accretion process onto the star. The density gradient is set to zero at , while the angular velocity there is fixed to be Keplerian. In the azimuthal direction, periodic boundary conditions for all variables are imposed.
These specified boundary conditions allow for a well defined quasi-stationary state if the planets are not allowed to respond to the disk.
After an initial relaxation procedure, where the orbital parameters of the planets are not changed and the disk structure is brought to an equilibrium, the planets are "released'' in all cases and are allowed to migrate (change their semi-major axes) in accordance with the gravitational disk forces/torques exerted on them.
During the evolution, material may enter from the disk into the Roche lobe of a planet. This material is partly removed from the simulation to model accretion onto the planet; for the detailed procedure see Kley (1999). In some models this mass is not added to the dynamical mass of the planets, so that mass is not strictly conserved. In other simulations the mass is added to that of the planet, maintaining conservation of mass. In Table 2 below it is indicated that only one of the presented models (h8a) has a varying dynamical mass. In general it does not change the outcome of a simulation noticably whether mass is accreted onto a planet or not.
For the viscosity a value of is used for all models. This is probably on the large side for protoplanetary disks, but it allows for a rapid evolution of the system and hence a reasonable computational effort; a larger speeds up the evolution and migration of the planet. It has been shown earlier that the migration speed has no influence on final magnitude of the eccentricities (Lee & Peale 2002) which tends to justify this approach. In addition, a larger -value leads to a gap that is not so well cleared (Kley 1999) which will tend to induce a larger periastron advance ( ) and an increased eccentricity damping. Both effects presumably serve to minimize the final eccentricity of the outer planet.
The density in the system is adjusted such that in the relaxed initial state there is approximately of material within the computational domain (see below).
The viscous terms, including all necessary tensor components, are treated explicitly. To ensure stability in the gap region, where there are very strong gradients in the density, an artificial bulk viscosity has been added, with a coefficient . For a detailed discussion of the viscosity related issues and tests, see Kley (1999).
The energy equation Eq. (1) is solved explicitly applying operator-splitting. The heating and cooling term D-Q is treated as one sub-step in this procedure, see Günther et al. (2004). The additional radiative diffusion part in the energy Eq. (5) is solved applying an implicit method to avoid possible time step limitations. To solve the resulting matrix equations we use Successive Over Relaxation (SOR) with an adaptive relaxation parameter (Kley 1989).
A larger mass ratio M/M* induces stronger torques and produces low densities in the gap region. To prevent numerical instabilitites caused by too large gradients we have found it preferable to work with a density floor, where the density cannot fall below a specified minimum value . For our purpose we use a value of in dimensionless values, where the typical (initial) density is of .
The forces of the disk are taken into account in a first order
approximation to reduce the computational effort.
To avoid problems with under-resolved material close to the planet,
a torque cutoff radius of
is applied where
material closer to the planet than
does not contribute to
the force acting on the planet.
The problem of choosing the optimum value for
Using a cutoff radius prevents large unphysical variations of
the forces acting on the planet. Very high resolution
simulations (D'Angelo et al. 2004,2002)
show a nearly symmetric distribution
of material close to the planet indicating that those regions
do not contribute too much to the total torque.
As these regions cannot be resolved in our simulations we have to use
a torque cutoff instead.
For all isothermal models we use
where the Hill radius is given by
In calculating the gravitational potential of the disk and planet the vertical extent of the disk is taken into account, assuming a vertically isothermal structure, i.e. Gaussian density distribution. For the smoothing length of the potential we choose . For a planet this is equivalent to only about .
In Fig. 3 the azimuthally averaged profile is shown for the relaxed equilibrium disk states for three different masses q2 of the outer planet starting from (black solid line), over (gray dotted line) upto (dashed line), with the last value close to the minimum mass of the outer planet of the GJ 876 system. One notices that for the two smaller planetary masses the outer disk is much more "quiescent'' than for the large mass. Larger planetary masses apparently lead to a restructuring of the disk. The wave-like disturbances in the disk (seen as the oscillatory behavior in the curves) are significantly stronger for .
|Figure 3: The azimuthally averaged density profile for the relaxed configurations for three different masses of the outer planet: and . The planet is located at a2 = 0.35 with a fixed semi-major axis, and is allowed to accrete. In these models the mass of the inner planet has been switched off. The density is given in dimensionless units.|
|Open with DEXTER|
|Figure 4: Gray scale plots of the surface density for the relaxed state with no inner planet for two different masses: top: and bottom: . Due to the higher planetary mass much stronger wave-like disturbances are created in the density, and the disk becomes eccentric with very small pattern speed in the inertial frame.|
|Open with DEXTER|
The existence of the two different equilibrium states of the disk is further illustrated in Fig. 4 where we display gray scale plots of the surface density for the relaxed state with no inner planet for two different masses (q2 = 3.5 and ). In both cases the small but non-vanishing eccentricity of the outer planet ( e2 = 0.01) leads to wave-like disturbances in the disk oscillating with the planetary period. But while for the lower mass case ( ) the disk structure remains quite regular, the second high mass case ( ) shows a strongly disturbed disk which has gained a significant eccentricity of about 0.25 near the gap edge. The transition from the non-eccentric state to the eccentric state which is here a function only of the planetary mass may depend also on the viscosity and temperature on the disk which we have held fixed. To check if this effect is induced by the non-zero eccentricity of the planet we have run a model with a vanishing eccentricity (e2 = 0) and found the same behavior. We find that the eccentric disk is nearly stationary in the inertial frame. The mass flow onto the planet and through the gap is higher than the flow for the non-eccentric disk that prevailed for the lower mass planet. While it would be interesting, a further exploration of these features is beyond the scope of this paper, and we leave it to a subsequent study. This disk eccentricity may be closely related to the disk eccentricity induced by resonant interaction of the outer 3:1 eccentric Lindblad resonance of planet 2 with the disk, as discussed by Papaloizou et al. (2001), and to the problem of the inner disk in cataclysmic variables by Lubow (1991). This violent interaction of a massive planet with the disk explains also some of the outlined numerical difficulties and requirements (such as density floor and artificial viscosity).
|Figure 5: The time evolution of the semi-major axis and eccentricity for three models with different masses of the outer planet, and zero mass inner planet.|
|Open with DEXTER|
|Figure 6: The time evolution of the periastron for three models with different masses of the outer planet. The thick dashed line is a fit corresponding to a 1 deg yr-1 shift.|
|Open with DEXTER|
In Figs. 5 and 6 the time evolution of the orbital elements of the outer planet for three models with different planet masses are displayed. All models have been relaxed before releasing the planet. For direct comparison the physical disk mass in all three cases has been scaled to be identical ( ). With respect to the reference model (c3) with q2 = 0.0059, this implies for models (d3a) and (e4a) a density reduction factor of 0.98 and 0.9, respectively (cp. Fig. 3). The migration rate of the planets depends primarily on the amount of mass near the 2:1 Lindblad resonance. We observe that the lowest ( q2 = 0.001) and highest ( q2 = 0.0059) mass planet have very similar initial migration rates, while the intermediate mass planet has a faster migration rate. The intermediate mass ( ) has the largest rate because there is more mass very close to the location of the 2:1 resonance, which lies here at r = 0.56(note the steep density gradient in Fig. 3). From the average density profile one might have suspected a smaller migration rate for the massive planet. However, due to the eccentric disk (Fig. 4), material gets close to the planet every orbit, which increases the migration rate slightly. Thus, the non-monotonic variation of migration speed with planet mass is related to the different equilibrium disk states, circular and eccentric. The right panel of Fig. 5 shows a temporary increase of the eccentricity followed by a decline for the larger mass planet case, while the smaller masses both show declining eccentricities at nearly the same rates indicating the standard eccentricity damping in planet-disk interaction. This different behavior of the eccentricity evolution is definitely caused by the eccentric disk state for the higher planet mass case. It remains to be studied if this effect may have some relation to the observed large eccentricities of the extrasolar planets.
In Fig. 6 the evolution of the periastron is displayed for the three cases. For comparison the superimposed black dashed-line refers to a periastron advance of exactly yr-1. This positive yr-1 (for the small and high mass planet) is solely due to the interaction of the disk with outer planet alone. The amount of planetary precession induced by the presence of the disk is determined primarily by the material being very close to the planet. Here, the behaviour of the low and high mass models are comparable again, while for the intermediate mass, there is less material very close to the planet and the precession rate is smallest.
|Figure 7: The relaxed azimuthally averaged density profile for two isothermal models with and without considering the inner planet (solid, dotted lines), and two radiative models including only heating and cooling (short-dashed), and additionally radiative diffusion (long-dashed). The mass of the inner planet is and that of the outer . The gray short-dashed line is identical to the corresponding line in Fig. 3 (labeled 5.9).|
|Open with DEXTER|
|Figure 8: The relaxed azimuthally averaged temperature profile for the isothermal (solid line) and two radiative models including only heating and cooling (short-dashed), and additionally radiative diffusion (long-dashed). The mass of the inner planet is and that of the outer .|
|Open with DEXTER|
Table 2: Parameters of the models used to study the evolution of a planetary system under the action of an outer disk. Given are a model name, the used thermodynamics (TD), isothermal (iso) or radiative (rad), the type of outer boundary condition (obc) (open or closed), the mass-ratio of the inner planet (q1), the mass ratio or evolution of the outer planet (q2), the figure in which the model is displayed. The displayed radiative models include heating and cooling only as the inclusion of the diffusion does not change results significantly.
In Fig. 8 the change in temperature due to the inclusion of radiative effects is displayed. Due to the relatively large density of disk and the high viscosity, the temperature in the radiative cases rise considerably above the isothermal case, and is given in the disk region by the equilibrium of heating and cooling D=Q. The disk thickness increases from H/r = 0.05 to about H/r = 0.15, for the full radiative models. Including radiative diffusion in the plane of the disk (long-dashed line) leads to higher temperature in the gap region. The full thermodynamic models are no longer eccentric even for the large planetary mass. The included radiation and larger H/r leads to a narrower gap and additional damping of the modes. We do not investigate at this point the detailed dependence of the disk thermodynamics on variations of the physical viscosity.
These fully relaxed models including the inner planet serve as the starting point for the dynamical evolution of the planets. The periastron advance for the outer planet due to the inner planet having q1 = 0.00175 and a1=0.20 alone, with no disk forces is found to be yr-1. Thus, the disk and the planet generate a shift of similar magnitude.
Having constructed several equilibrium models holding the orbital elements of the planets fixed, we now release the planets and follow the evolution of their orbital elements. For an overview Table 2 lists the models and their most important parameters.
In the first instance we keep the damped and reflecting boundary conditions at the outer boundary, i.e. we model the situation where the disk remains in full contact with the planet during the evolution. In Figs. 9 to 11 we display the evolution of the semi-major axes, the eccentricities, and the two resonant angles and for three models, a test case (model h4) with a lighter outer planet ( q1 = 0.00175, q2 = 0.0035), a case (h2) resembling the edge-on ( q1 = 0.00175, q2 = 0.0059) system, and finally the ( q1 = 0.0021, q2 = 0.0071) case (h9), respectively. In all cases the outer planet captures the inner one into a 2:1 resonance. Although the orbits are initially quite far from 2:1 mean-motion commensurability, the resonant angle is already in libration, with amplitude as large as . However, this libration of has little dynamical consequence and the increase in e1 is slow, until the mean motions approach the 2:1 commensurability. There is a delay in the capture of (and ) into resonance, and the libration amplitude of is smaller than that of . Note that in the first test-case for smaller planet mass, when the initial disk state was not eccentric, the capture of occurs faster and the libration amplitudes of both and are much smaller.
|Figure 9: The time evolution of the orbital elements ( ) for an isothermal model (h4) with and (test case with lower outer planet mass).|
|Open with DEXTER|
|Figure 10: The time evolution of the orbital elements ( ) for an isothermal model (h2) with and (as GJ 876 edge on).|
|Open with DEXTER|
|Figure 11: The time evolution of the orbital elements ( ) for an isothermal model (h9) with and (as GJ 876 at inclination).|
|Open with DEXTER|
|Figure 12: The time evolution of the orbital elements ( ) for a radiative model including heating and cooling (m1) with and (as GJ 876 edge on). Note that at around t=380 the outer planet leaves the computational grid and the results are not reliable anymore.|
|Open with DEXTER|
The fact that (and even and ) can be librating when the orbital mean-motions are far from the 2:1 commensurability is due to the 1/ej dependence of the resonance-induced retrograde precession of at small ej. The relative timing of the capture of and into resonance is affected by the masses and initial eccentricities of the planets. We shall discuss these points further in Sect. 5.
Once both and are librating about , the eccentricities rise rapidly. In cases h2 and h9, applicable to GJ 876, the eccentricity of the inner planet rises above 0.4, which clearly exceeds the upper limit of 0.31 for GJ 876 (see Fig. 1).
The evolution with the inclusion of heating and cooling (m1) is displayed in Fig. 12. Again the system is caught in a 2:1 resonance, and the eccentricities rise to high values. Here, due to the smoother initial state, the alignment of the resonant angles is much stronger and the libration amplitudes for both angles are reduced strongly. The evolution beyond t=380 yr is unreliable because both planets are outside the computational domain by then. The models with radiative diffusion do not yield any different results from models with only heating and cooling. The equilibrium configurations as shown in Figs. 7 and 8 above are very similar already.
In summary, there exist two different types of evolutions of the orbital elements depending on the state of disk. For a nearly circular disk (models h4, m1) the delay in the capture of (or equivalently ) allows e2 to be initially damped. After capture in this second resonance e2 rises rapidly and increases far beyond the upper limit inferred from the observations. The libration of the resonant angles after capture is very small. This result confirms earlier findings of hydrodynamical evolutions of resonant planets (Kley et al. 2004). For an eccentric disk state, i.e. for larger outer planet masses, capture of is more delayed and the libration of the two resonant angles remains quite large. Nevertheless, the eccentricities quickly exceed the upper limits for GJ 876.
|Figure 13: The time evolution of the orbital elements ( ) for an isothermal model (k3b) with and , including disk dispersal.|
|Open with DEXTER|
In Figs. 13 and 14 two isothermal cases
are presented; in the first (k3b), the mass of the outer planet remains
constant, while in the second model (h8a) the planet is allowed to accrete
from its surroundings and to grow in mass. In the second case the mass
of the inner planet was chosen higher
q1 = 0.0021 to compensate for
the higher final mass of the outer planet. Again, capture into the
2:1 resonance occurs, but the eccentricities do not rise to
such high values. The loss of mass from the disk turns the migration
off before the eccentricities get large. The alignment of the
resonant angles occurs on
longer timescales than before. Similar behavior can be seen in the
evolution of the radiative model (m2) displayed in Fig. 15. Again
the alignment of the resonant angles occurs faster and
libration of the resonant angles is smaller than in the isothermal
models due to the non-eccentric disk state.
|Figure 14: The time evolution of the orbital elements ( ) for an isothermal model (h8a) with and a growing outer planet . The model includes disk dispersal.|
|Open with DEXTER|
|Figure 15: The time evolution of the orbital elements ( ) for a radiative model (m2) with and , including disk dispersal.|
|Open with DEXTER|
|Figure 16: Time evolution of the orbital elements ( ) for the baseline three-body model with and . The planets are initially on coplanar circular orbits. The outer planet is forced to migrate inward at the rate , and there is no eccentricity damping or additional apsidal precession.|
|Open with DEXTER|
Figure 16 shows the evolution of the semimajor axes, eccentricities, and resonant angles and for the baseline three-body model with and . The outer planet is forced to migrate inward at the rate (which matches the migration rate at for the h2 model in Fig. 10), and there is no eccentricity damping or additional apsidal precession. The planets are initially on coplanar circular orbits, with and . The circular orbits and larger a2 compared to the hydrodynamical simulations (where ) reduce the initial eccentricity variations. As in the hydrodynamical simulations, the resonant angle is already librating about at the start of the evolution, even though the planets are started even further from the 2:1 mean-motion commensurability than in the hydrodynamical simulations. With the small initial eccentricity variations and the use of Jacobi coordinates, it can be seen in Fig. 16 that (and ) are first captured into libration about at and that spends more time around even before the capture. The system passes smoothly over to the configuration with both and librating about when (see also Lee 2004). The small libration amplitudes of the configuration with both and librating about are similar to those found in the radiative model m1 (Fig. 12).
For two planets orbiting a star in coplanar orbits, the equations of
motion for the periapse longitude and eccentricity are
When the eccentricities are very small, the relevant terms in are
and we have
The relative timing of the capture of and into resonance needs some clarification. Due to the 1/e1 dependence of the resonance-induced precession of at small e1 (Eq. (11)), a small initial value of e1 leads to an extremely mobile longitude of periapse , such that the large mass of the outer planet can induce sufficient retrograde motion of to cause to already be librating at the start of the evolution when the orbital mean motions are far from the 2:1 commensurability. This initial libration is evident in the three-body model shown in Fig. 16 and in the hydrodynamical simulations shown in Sect. 4.2. It should be noted, however, that this libration of has little dynamical consequence and the increase in e1 (see next paragraph) is slow, until the mean motions approach the 2:1 commensurability. In spite of the relatively small initial value of e2 in Fig. 16 and the same 1/e dependence of the resonance-induced precession (Eq. (12)), and hence are not librating initially, although they do spend more time near . This is due to the inner planet's lower mass being insufficient to perturb the more massive outer planet's enough to cause to librate that far from the 2:1 commensurability of the mean motions.
The reason the eccentricities increase with continued migration can be understood as follows. Continued inward migration of the outer planet while the system is within the resonances means a2 is consistently slightly smaller than it would have been without the migration. The increased value of n2 therefore means is slightly larger at any instant than it would be without the migration. The argument of the sine in Eq. (13), , is then slightly greater than on average. Since C1 < 0, and the eccentricity must grow. A similar argument applies to if is also librating about (see Eq. (14)). It can be shown from the equations for the evolution of the orbital energy and angular momentum that at least one, but not necessarily both, of the eccentricities must be increasing at any given time for continued migration within the resonance, if e22 < (3 + e12)/4for the 2:1 resonance and there is no eccentricity damping (see Eq. (A10) of Murray et al. 2002).
For , the resonant interaction is no longer dominated by the lowest order resonant terms when , and e2 changes from increasing to decreasing with continued migration until when . This change in the evolution of e2 is not obvious in Fig. 16 because during this phase; see Fig. 5 of Lee (2004). Then the libration center of and changes from to (with a slight phase shift due to the continued migration), and both eccentricities continue to grow because of the slight phase shift in the arguments of the terms involving linear combinations of and . However, because there are now many terms in the disturbing potential that contribute to the resonant interaction, a simple demonstration of the co-precession of the periapses with both and librating about and of the increase in both eccentricities with continued migration is elusive. Continued migration causes the eccentricities to quickly exceed those observed without a large eccentricity damping as discussed above. Although the hydrodynamical simulations self-consistently produce , the magnitude of this damping is far less than that necessary to maintain the small observed eccentricities in the GJ 876 system during continued migration.
|Figure 17: Time evolution of the orbital eccentricities for the three-body models with additional apsidal precession and p = 1 and 10, compared to the case without additional apsidal precession (p = 0) from Fig. 16.|
|Open with DEXTER|
We have performed three-body integrations similar to that in Fig. 16, but with additional apsidal precession and p = 1 and 10. The much smaller disk-induced precession of the inner planet's orbit was neglected. Except for a larger libration amplitude for in the p = 10 case, the evolution of and is similar to that shown in Fig. 16 for p = 0, and there is no additional delay in the capture of . Figure 17 shows the evolution of the eccentricities for p = 0, 1, and 10. At a given time, the additional positive apsidal precession results in an increase in e1 and a decrease in e2, but the effect is small for e1 < 0.31 (the upper limit for GJ 876; Fig. 1) even for p = 10. The additional apsidal precession and forced migration rates used in the three-body integrations are proportional to a2-3/2 and hence the mean motion n2. A decrease in the positive, disk-induced precession rate due to, e.g., disk dispersal would cause the eccentricities to adjust to values appropriate for the precession rate at a given time and approach the p = 0 case without disk-induced precession.
The analytic theory developed by Yoder & Peale (1981) for the 2:1
resonances of Io and Europa and the Laplace resonance takes into
account the apsidal precession induced by the oblateness of Jupiter
(which is largest for the innermost satellite Io), and it can be
adapted to understand the effects of the disk-induced apsidal
precession (which is larger for the outer planet).
As the orbits converge toward the 2:1 mean-motion commensurability
2 n2 - n1 < 0 increases), the resonance
would be encountered before the
apsidal precession is dominated by
However, because the resonance-induced retrograde apsidal precession
is proportional to 1/ej and much larger in magnitude than the
disk-induced prograde precession for small ej,
can be captured into
libration in a sequence that differs little in order or timing from
the case where there is no disk-induced precession.
When we include
Eqs. (11) and (12) for small ej, stable
retrograde precessions of both orbits still require
requirement that the orbits precess at the same rate on average
implies the following relationship between the forced eccentricities:
|Figure 18: Same as Fig. 16, but for the three-body model with the following initial conditions: e1 = 0.01, e2 = 0.05, and the orbits are antialigned, with the inner planet at periapse and the outer planet at apoapse.|
|Open with DEXTER|
The "isothermal'' simulations with fixed disk temperature structure have shown that more massive planets with a planet-to-star mass ratio of are able to perturb the disk sufficiently to make it eccentric (even if the planetary orbit is circular). This effect may be caused by an eccentric instability driven by an interaction of the m=2 mode at the outer eccentric 3:1 Lindblad resonance with a slightly eccentric disk (see Papaloizou et al. 2001). For smaller masses in the isothermal models and for all masses in the radiative models, the planet-disk interaction produces the typical gap and spiral arms, but the disk remains otherwise circular.
The simulations also confirm that for non-eccentric, circular disks the eccentricity growth of the outer planet is suppressed until is locked into libration about (see models h4, m1, m2). However, and are always eventually captured into libration, and the subsequent increase in the eccentricities of both planets past the observational upper limits in GJ 876 ( , ) occurs on a time scale shorter than the expected viscous time scale of the disk. The final eccentricities of both planets are then substantially larger than those seen in GJ 876, a result which was found already in previous hydrodynamic (Kley et al. 2004; Papaloizou 2003) and three-body simulations (e.g. Lee & Peale 2002).
In addition to causing an increased delay in the capture of and , eccentric disks (models h2, h9) lead to significantly larger libration amplitudes for ( ) than for the circular disk case. Interestingly, it appears that in the case of GJ 876 the libration of is in fact larger ( for ) than that of ( for ), as shown recently by Laughlin et al. (2004).
Lee & Peale (2002) have shown that if sufficient eccentricity damping ( ) is applied during migration, the observed configuration of GJ 876 can be maintained for arbitrary migration times. We have updated the amount of eccentricity damping required using the updated dynamical fits by Laughlin et al. (2004). For the best fits with coplanar inclination , , , and must exceed . However, no such rapid eccentricity damping mechanism is presently known. The hydrodynamic simulations typically give comparable time scales for eccentricity damping and semi-major axis decrease. Some current theories of planet-disk interactions indicate that eccentricity driving should occur (Goldreich & Sari 2003; Ogilvie & Lubow 2003). Such eccentricity behavior would not be consistent with the observed state of GJ 876 if there were extensive migration after capture into the resonances.
If the disk is removed soon after the planets are captured into resonance, as in Figs. 13 to 15, the driving disappears and the final eccentricities are smaller, like those observed. However, the disk must vanish before the post-capture orbits shrink by 10%, similar to the result found in 3-body integrations by Lee & Peale (2002). We have used in the present simulations a high viscosity and high surface density to increase the migration rates, and to be able to perform the simulations in a reasonable amount of computer time. In turn this implies a need for rapid disk dispersal to halt migration. If more realistic values were used for the viscosity and the disk mass, the necessary disk dispersal time scale could be much more extended. What appears to be somewhat fine tuned in the present simulations would be perhaps more reasonable. Our scenario only requires that capture occurred in the final formation phase of GJ 876, and that the planets did not migrate over a very long distance while locked in the resonances.
Additional effects that might influence the eccentricity damping in these type of simulations are possibly three-dimensional. Fully 3D-simulations of circular planets on a fixed orbit (Bate et al. 2003; D'Angelo et al. 2003b) have shown that the spiral structures are not so clear and material may enter the Roche lobe from above the midplane, an effect which may alter the eccentricity damping properties. The inclusion of magnetic fields might lead to a magnetic coupling of the planetary field with that of the disk, which will have an influence on the eccentricity evolution. But if the consideration of these and perhaps other processes fail to produce sufficient eccentricity damping, the elimination of the disk shortly after the capture of the GJ 876 planets into the 2:1 resonances is a possible, albeit perhaps unsatisfying means of accounting for the observed low-eccentricity configuration.
We would like to thank Greg Laughlin for stimulating discussions during the course of this project and for furnishing the details of the dynamical fits by Laughlin et al. (2004). We also thank D. N. C. Lin and M. Nagasawa for informative discussions. The work was sponsored by the KITP Program "Planet Formation: Terrestrial and Extra-Solar'' held in Santa Barbara from January 2004 until March 2004. We thank the organizers for providing a pleasant and very stimulating atmosphere. This research was supported in part by the National Science Foundation under Grant No. PHY99-0794 and by NASA grants NAG5-11666 and NAG5-13149.