Eccentric Extrasolar Planets: The Jumping Jupiter Model

Eccentric Extrasolar Planets: The Jumping Jupiter Model

Icarus 156, 570–579 (2002) doi:10.1006/icar.2001.6786, available online at http://www.idealibrary.com on Eccentric Extrasolar Planets: The Jumping Ju...

864KB Sizes 0 Downloads 56 Views

Icarus 156, 570–579 (2002) doi:10.1006/icar.2001.6786, available online at http://www.idealibrary.com on

Eccentric Extrasolar Planets: The Jumping Jupiter Model F. Marzari Dipartimento di Fisica, Universit`a di Padova, Via Marzolo 8, I-35131 Padova, Italy E-mail: [email protected]

and S. J. Weidenschilling Planetary Science Institute, 620 North 6th Avenue, Tucson, Arizona 85705-8331 Received July 13, 2001; revised November 1, 2001

Most extrasolar planets discovered to date are more massive than Jupiter, in surprisingly small orbits (semimajor axes less than 3 AU). Many of these have significant orbital eccentricities. Such orbits may be the product of dynamical interactions in multiplanet systems. We examine outcomes of such evolution in systems of three Jupitermass planets around a solar-mass star by integration of their orbits in three dimensions. Such systems are unstable for a broad range of initial conditions, with mutual perturbations leading to crossing orbits and close encounters. The time scale for instability to develop depends on the initial orbital spacing; some configurations become chaotic after delays exceeding 108 y. The most common outcome of gravitational scattering by close encounters is hyperbolic ejection of one planet. Of the two survivors, one is moved closer to the star and the other is left in a distant orbit; for systems with equal-mass planets, there is no correlation between initial and final orbital positions. Both survivors may have significant eccentricities, and the mutual inclination of their orbits can be large. The inner survivor’s semimajor axis is usually about half that of the innermost starting orbit. Gravitational scattering alone cannot produce the observed excess of “hot Jupiters” in close circular orbits. However, those scattered planets with large eccentricities and small periastron distances may become circularized if tidal dissipation is effective. Most stars with a massive planet in an eccentric orbit should have at least one additional planet of comparable mass in a more distant orbit. c 2002 Elsevier Science (USA)

1. INTRODUCTION

More than 70 extrasolar planets orbiting solar-type stars have been discovered to date by various surveys (Marcy et al. 2000), by precise Doppler measurements of stellar radial velocities. The masses of these new objects (actually minimum values, as the measured quantity is m sin(i)) range from about 0.4 to 10 Jupiter masses, suggesting that they are jovian-type planets, probably with an icy/rocky core and a massive gaseous envelope. An intriguing aspect of the newly discovered extrasolar planets is that

their orbits are closer—sometimes much closer—to their stars than Jupiter’s distance from the Sun. None is yet known with semimajor axis greater than 4 AU. Some, such as 51 Peg, the first to be found, have extremely small orbits, with semimajor axes less than 0.1 AU and periods of only a few days. This peculiarity facilitates their detection with short-term surveys but probably overestimates their abundance relative to planets with longer periods. Only a few percent of surveyed stars have such close planets. Extrasolar planets at Jupiter’s distance from their stars may well exist and be even more frequent, but they are harder to detect and require longer observation times, comparable to their orbital periods. Another surprising finding is the highly eccentric shape that characterizes the orbits of most of the extrasolar planets. Values of eccentricities larger than 0.9 are observed, and only planets with small semimajor axes (less than about 0.1 AU) have nearly circular orbits, comparable to the maximum value among giant planets in our solar system (e = 0.054 for Saturn). The observed orbits are summarized in Fig. 1. These peculiarities complicate the task of interpreting the origins of the observed extrasolar systems with the conventional core accretion model for the formation of giant planets. Saturn and Jupiter are believed to have formed by gradual accumulation of planetesimals into a core with mass at least a few times that of Earth, followed by gravitational capture of a gaseous envelope from the surrounding nebula (Pollack et al. 1996). The accumulation of such a core would have occurred beyond the nebular “snow line” (believed to be at about 4 AU) where ices condensed in the primordial solar nebula, increasing the amount of solids available for accretion. Jupiter and Saturn did form on almost circular orbits, as expected for protoplanets evolving from an accretion disk. Recently, Bodenheimer et al. (2000) have shown that extrasolar giant planets may form in orbits close to their stars, according to the core accretion model, but this process requires extreme assumptions as to the surface density of the disk and efficiency of accretion. The authors do not exclude that orbital migration was

570 0019-1035/02 $35.00 c 2002 Elsevier Science (USA)  All rights reserved.

THE JUMPING JUPITER MODEL 1

0.8

0.6

0.4

0.2

0 0

1

2

3

a (AU)

FIG. 1. Distribution of eccentricities versus semimajor axes for the extrasolar planets discovered through April 2001. Data from the online Extrasolar Planets Encyclopedia, http://www.obspm.fr/planets.

also an important component; it may possibly explain the low mean density (∼0.4 g cm−3 ) of the planetary companion of star HD 209458 detected by its transit across the star (Charbonneau et al. 1999). Various theories have been proposed based on the assumption that the extrasolar planets originated at large distances from their stars, as did Jupiter and Saturn in our solar system, and subsequently migrated inward to their present positions and acquired large orbital eccentricities. One proposed mechanism for migration is tidal coupling to an evolving protoplanetary nebula (Lin et al. 1996) However, a planet as massive as those observed would open a gap in the nebula (Ward and Hahn 2000); the planet would be constrained to evolve with the disk. That mechanism requires giant planets to form rapidly, on the time scale of nebular viscous evolution, probably less than 106 years. Also, the source of the observed large eccentricities is not explained; while some disk–planet tidal interactions can increase eccentricity (Lin et al. 2000), the expressions break down at large e, and Ward and Hahn (2000) conclude that the net effect of disk–planet interactions may be damping of eccentricity. In this paper we analyze in detail another mechanism proposed to explain the orbital peculiarities of the extrasolar planets: the “Jumping Jupiters” model. It assumes that two or more giant planets did form in the protoplanetary disk and that, subsequently, their orbits became unstable due to their mutual gravitational perturbations (Farinella 1980, Rasio and Ford 1996, Weidenschilling and Marzari 1996, Lin and Ida 1997, Ford et al. 2001). Migration occurs by gravitational scattering during close encounters. The time scale for the onset of instability depends on the initial separation among the bodies. Systems with three or more planets are less stable than systems with only two initial planets; they allow larger initial separations and, consequently, wider feeding zones for the growing protoplanets. Due to mutual

571

perturbations, the eccentricities of the bodies grow and produce crossing orbits. Repeated close encounters among the planets generate a period of chaotic evolution that usually terminates with the ejection of one planet on a hyperbolic trajectory. The two remaining planets are usually on stable eccentric orbits with the inner one close to the star and the other in a much larger orbit, such that close approaches are no longer possible. The “Jumping Jupiter” mechanism has two appealing features. It accounts for inward orbital migration of one of the planets that supplies the orbital energy to the ejected one, and it excites orbital eccentricities consistent with those observed. Since gravitational encounters can produce eccentric orbits with small periastron, tidal dissipation could circularize the orbit near periastron distance and explain the close and circular orbits of 51 Peg type systems. The model can also be invoked together with other evolutionary models to explain planets with eccentric orbits. In this context, the gravitational scattering of planets may account for the excitation of the orbital eccentricity following the formation phase. In this paper we describe in detail the implications of the “Jumping Jupiter” model for different initial conditions. We also explore the predictions of the model in terms of final observable configurations of possible extrasolar planetary systems. 2. THE “JUMPING JUPITER” MODEL IN DETAIL

The basic scenario for the “Jumping Jupiter” model is the evolution of an accretion disk into a multiple planetary system with three or more planets. In the conventional model for runaway growth of protoplanetary embryos, their accretion is a stochastic process, which leads to a series of embryos with comparable masses and more or less equal orbital spacings. (Kokubo and Ida [1998] described this process, which they called “oligarchic growth.”) Multiple cores grow in orbits that are stable during the planetesimal accumulation phase. However, if they reach the critical mass for capture of a massive gaseous envelope from the surrounding nebula, their masses may increase on a short time scale by an order of magnitude or more. This increase enlarges the Hill radius and can render adjacent orbits unstable when the protoplanets reach their final mass at the end of the gas accretion. This evolution was suggested by Weidenschilling and Marzari (1996) for extrasolar planet systems and by Thommes et al. (1999) for the planets of the solar system. The final configuration has to be unstable so that orbits can cross and start the chaotic evolution of the system. The issue of stability in a multiple planetary system is relevant for the model since it gives the order of magnitude of the initial separation that the embryos may have during the formation of the planets and the time scale before the start of the chaotic behavior. Lin and Ida (1997) considered systems with up to nine jovianmass planets, but the formation of so many giant planets in a single system is problematic and the results are harder to interpret. Ford et al. (2001) considered interactions of only two planets. However, three-planet systems allow a much more interesting

572

MARZARI AND WEIDENSCHILLING

range of outcomes of dynamical evolution, as we will demonstrate. 2.1. Instability of a Three-Planet System In the case of two planets of masses m 1 , m 2 orbiting a star of mass M∗ , a general analytical criterion can be invoked to state the stability of the system with initially circular orbits. If we measure the initial orbital separation , defined as the semimajor axis difference a2 − a1 , then stability is guaranteed when √  ≥ 2 3R H (1,2) ,

(1)

where the mutual Hill radius is defined as 

R H (1,2)

m1 + m2 = 3M∗

1/3

(a1 + a2 ) . 2

(2)

If condition (1) is met, the orbits of the two planets will never cross as the result of their mutual perturbations. This criterion was formulated by Gladman (1993) and in a more complete form by Hasegawa and Nakazawa (1990). The situation is more complex for a system containing more than two planets. Chambers et al. (1996) performed numerical integrations of such systems and found that mutual perturbations among three (or more) bodies of comparable mass produce crossing orbits even though each pair of bodies considered separately meets the Hill stability criterion. Instability occurs even √ for initial distances significantly larger than 2 3R H . They also showed that the logarithm of the instability time scale (defined as the time until the first close encounter within one Hill radius) depends linearly on . Their integrations were performed with masses from 10−7 up to 10−5 M , i.e., about 3.3 M⊕ . An increase of the time scale for instability for a given value of  was observed for the cases where larger planetary masses were assumed. Chambers et al. assumed initially circular and coplanar orbits for their integrations. The question of instability was examined further by Ito and Tanikawa (1999), who showed that the time scale for orbit crossing to begin was reduced by up to two orders of magnitude if the planetary embryos were assumed to have initial random velocities (eccentricities and inclinations). The two results are consistent, as the Chambers et al. cases required long times for eccentricities to build up from zero to match the values assumed by Ito and Tanikawa as their initial state. The latter also considered protoplanets with masses 10−7 M . In our model we consider systems of three Jupiter-size planets with masses of the order of 10−3 M , two orders of magnitude larger than the maximum value considered by Chambers et al. and Ito and Tanikawa in their work. We show that for such massive planets, resonances affect the instability time scales such that they are not simply monotonic with separation. To estimate the instability time scale for these systems, we performed numerical integrations of a variety of three planet systems with different initial separations. We followed a similar strategy as in Chambers

FIG. 2. Logarithm of the time to the first close encounter between any two planets versus the initial separation, expressed as K times the mutual Hill radius, for three equal mass planets (m = 10−3 M ). For each value of K at intervals of 0.1, we integrated six cases, varying only the initial position angles of the planets. The continuous line shows the trend we would expect without the resonant perturbations following the linear relation found by Chambers et al. (1996). Our fitting line is approximate and is derived from the data points that are far from the resonance gaps.

et al. starting the planets on circular orbits around a star of one solar mass at semimajor axes a1 , a2 = a1 + KR H (1,2) , and a3 = a2 + KR H (2,3) . Unlike Chambers et al. we did not assume coplanar orbits, but like Ito and Tanikawa we integrated fully tridimensional cases. We assumed initial inclinations of 0.5◦ for the first planet, 1.0◦ for the second, and 1.5◦ for the third. The orbits of the planets have been computed with the mixed variable symplectic integrator (Levison and Duncan, 1994: SWIFT package) adopting a fixed stepsize of 20 days. We first considered a system of three equal mass jovian planets (m = 10−3 M ). We integrated orbits for a range of initial spacings with K between 3 and 5.3, keeping a1 constant at 5 AU. The smallest value of K gives a2 = 6.51 and a3 = 8.47 AU. The largest value gives a2 = 8.01 and a3 = 12.84 AU. In Fig. 2, the time scale for instability, measured as the interval of time from the beginning of the integration until the first close encounter between any pair of planets within a Hill radius, is shown as a function of the parameter K . For each value of K we ran six integrations with different initial values of relative position angles for the planets. Our instability times are significantly longer, for the same K , compared to those found by Chambers et al. (1996) for less massive planets (their largest mass value was m = 10−5 M ). This was an expected result since Chambers et al. also compared planetary systems with different masses and observed an increasing trend in the instability times for more massive planets at large values of K . A surprising feature is instead the strong effect of mean motion resonances among the planets. Chambers et al. observed only weak structures possibly related to resonances in their results. The orbital spacings between their neighboring planetary embryos permitted only

THE JUMPING JUPITER MODEL

573

delayed to 108 –109 years after planetary formation. Thus, the three-planet model is consistent with the core accretion model for the formation of giant planets. However, it does not exclude the possibility of rapid formation by gravitational instability of the disk (Boss 1997), although the formation of more than two planets by that mechanism has not yet been demonstrated. 2.2. The Chaotic Phase

FIG. 3. Logarithm of the instability time versus K for three planets of unequal mass. The planet in the middle has twice the mass of the other two (m 2 = 2 × 10−3 M ). The continuous line has the same meaning as in Fig. 2.

high-order resonant configurations while in our simulations the massive planets encounter strong low-order resonances. This difference is due to the fact that commensurability resonances are at fixed ratios of semimajor axes, while the Hill radius increases with mass; for sufficiently large planetary masses the distances to low-order resonances are a small number of Hill radii. At K  5 we notice a strong reduction in the time scale for instability due to a 2 : 1 resonance between the orbital periods of the first and second body, and between the second and the third. The choice of an equal spacing proportional to R H inevitably leads to this resonance. However, the large range of K values for which the time scale is reduced is surprising, and it proves that even systems near but not precisely in the resonance are strongly affected by it. The signature of a 3 : 1 resonance between the third and first body is also visible around K  3.8. At the maximum initial separation, the time scale for instability can be as large as 109 years. When we consider systems with planets of unequal mass, the pattern is more difficult to interpret. In Fig. 3 we increase the mass of the planet in the intermediate orbit by a factor of 2, and we again see the reduction of the instability time scales caused by the 2 : 1 resonance between the three planets. In this case, K = 3 corresponds to a2 = 6.77, a3 = 9.15 AU. The widest system we considered had K = 5.8, a2 = 9.09, a3 = 16.5 AU. The mutual 5 : 2 resonance is responsible for the plateau in the instability time scale at values of K  5.5. One significant result is the demonstration of the large range of time scales for instability of the orbits with the variation of initial spacing. Ford et al. (2001) argued against systems with three or more planets, on the grounds that they would be unlikely to accrete at precisely the same time. However, we find a range of spacings that allows ample time for core formation and capture of gas during the lifetime of a circumstellar disk (106 –107 years) before the onset of orbital crossings; indeed, instabilities can be

After the onset of instability, the eccentricity of the planets keeps growing until the orbits cross. Close encounters among the planets induce large changes in the orbital elements and cause a prolonged period of chaotic evolution that may last a few million years. The final outcome in the vast majority of cases is that one planet is ejected on an hyperbolic trajectory while two surviving planets are on widely separated stable orbits. If a slight inclination is allowed for the initial orbits, as is reasonable to assume in a realistic system, only in 2 to 5% of the cases we observed a collision, and scattering is by far the preferred outcome. We assumed no enhancement of collisional cross section beyond the usual gravitational focusing and planetary radii equal to Jupiter’s present value; in young systems the radii could be somewhat enhanced before cooling of the planets. Ford et al. (2001) investigated the effects of varying the assumed planetary radius in two-planet systems. For normal-sized “Jupiters,” their collision frequency was somewhat higher (10%), presumably due to lower mean encounter velocities; the collision frequency increases sharply with planetary radius. The numerical computation of the orbits in all the simulations that include the chaotic phase have been performed with a mixed algorithm: the RA15 version of the RADAU integrator (Everhart 1985) during close encounters, and the mixed-variable symplectic integrator from SWIFT elsewhere with a time step of two days. This short time step was required by the highly eccentric orbits obtained by the planets during the scattering phase. This approach is similar to that of Chambers and Wetherill (1998) who instead used a Bulirsch–Stoer method when the massive bodies are close. We believe RADAU is slightly more efficient in handling close encounters than Bulirsch–Stoer. In Fig. 4 we illustrate an example of the “Jumping Jupiter” orbital evolution. Three Jupiter-mass planets initially set on circular orbits with small inclinations (0.5, 1.0, and 1.5 degrees) evolve chaotically after a short (4 × 104 year) phase of noncrossing orbits. The initial semimajor axes are selected in order to shorten the time scale before the onset of crossing orbits, thereby reducing the computational effort. After 2 × 105 years, the planet originally in the innermost orbit is ejected while the planet initially in the middle is left in a distant orbit with a = 32 AU and e = 0.66. The outermost original planet is left in a close orbit with a = 2.41 AU and e = 0.36. These orbits are stable against further close encounters. Assuming in situ formation of giant planets close to the star from massive disks (Bodenheimer et al. 2000), the “Jumping Jupiter” model can still be invoked to excite large orbital eccentricities. As shown in Fig. 5, three planets assumed to form on

574

MARZARI AND WEIDENSCHILLING 1

0.8 0.6 0.4 0.2 0 0 Time (yr) 40

30

20

10

0 0

a3 = 9.50 AU, zero eccentricity and inclinations of 0.5, 1.0, and 1.5 degrees, respectively. The initial orbital angles were generated randomly for each system. The choice of the initial semimajor axes is a good compromise between reasonable initial separation and fast onset of instability. Larger initial separations, while more plausible from the standpoint of planetary formation, would have led to significantly longer integration times, even with the fastest available computers. Each numerical simulation was stopped after one of the planets was ejected from the system. In Fig. 6 we present histograms of the distributions of the orbital elements of the planet injected in the inner orbit. The semimajor axes are concentrated around 2.3 AU, while the eccentricities are peaked at 0.6. The concentration of semimajor axes at about 2.3 AU is related to the conservation of energy and, then, on the initial choice of the planets’ semimajor axes. The initial configuration of three planets has a total energy (neglecting terms of order m/M∗ ) of

Time (yr)

FIG. 4. Typical evolution of the semimajor axis and eccentricity of a system of three jovian-mass planets, with initial semimajor axes of a1 = 5.0, a2 = 7.25, and a3 = 9.5 AU. The orbits become crossing after about 4 × 104 years and, after a series of close encounters, the originally outermost planet is injected into an orbit close to the star while the inner planet is ejected hyperbolically.

initially circular orbits within 2 AU from the star begin to have close encounters until one of the planets is ejected. At the end of the chaotic phase the inner survivor has a semimajor axis of 0.45 AU and an eccentricity of 0.76. The outer one is placed on an orbit with semimajor axis of 6.70 AU and an e = 0.35. Scattering is by far the more common event also at small stellar distances. Lin and Ida (1997) showed that with an enhanced collisional radius by a factor of 2 compared with Jupiter’s size, impacts were rare for initial orbits inside 1 AU. The existence of other planetary companions is a good test for the “Jumping Jupiter” model. The gravitational scattering process requires that the star with one close planet should have in nearly all cases at least one other planet of comparable mass in a more distant orbit. This is in good agreement with the fact that roughly 50% of the stars for which a single planet has been detected show evidence of additional companions, as seen in coherent variations of the residuals to the Keplerian fit (Fischer et al. 2001). 2.3. Final Distribution of the Orbital Elements The dynamical evolution of a three-planet system is highly chaotic and the outcome is very sensitive to the choice of the initial orbital parameters. The chaotic nature is mostly due to the large number of mutual encounters among the planets before the final ejection of one or two planets. To explore the range of possible final configurations, we resorted to a statistical approach. We integrated 64 different systems characterized by fixed initial semimajor axes a1 = 5.00, a2 = 7.25,

  G M∗ m 1 m2 m3 E =− . + + 2 a1 a2 a3

(3)

If there are no collisions to dissipate energy, the final configuration, with one planet in a close orbit, one in a distant orbit, and one escaping on a barely hyperbolic orbit, must have the same energy. If the outer survivor’s orbit is so distant that its energy is

1 0.8 0.6 0.4 0.2 0 0 Time (yr) 8

6

4

2

0 0 Time (yr)

FIG. 5. Orbital evolution of three equal mass jovian planets assumed to have formed within 2 AU from the star. The initial semimajor axes are a1 = 1.00, a2 = 1.50, and a3 = 2.00 AU, respectively. The initial orbital eccentricities are all zero, while the inclinations are 0.5◦ for the first planet, 1.0◦ for the second, and 1.5◦ for the third. Here the planet originally in the middle orbit is ejected, and the originally outermost ends up in the innermost position. The planet left in the inner orbit after the chaotic phase has a = 0.45 AU and e = 0.76. The outer surviving planet has a  6.5 AU.

THE JUMPING JUPITER MODEL

FIG. 6. Histograms showing the statistical distribution of the orbital elements for the planet left in the inner orbit at the end of the chaotic phase of a jovian planet trio. The inclination is relative to the initial reference plane. In each system contributing to the statistical sample the elements of the inner planet are computed just after the ejection of one of the other two bodies. All the planets have the same mass (10−3 M ).

negligible compared with that of the inner planet, then one can show that its final semimajor axis, a f , is af 

−G M∗ m f , 2E

575

0 to 1. The semimajor axes have a broad distribution from about 5 to 90 AU, with peaks around 10, 30, and 60 AU, but it is not clear whether these peaks are significant or if they are related to the choice of initial semimajor axes of the three planets. The breadth of the distribution is due to the small (negative) energy of the outer orbit, which makes only a minor contribution to Eq. (5). In most cases the orbital period of the outer survivor is longer than 30 years. For both inner and outer survivors, the inclinations have a broad range of values with respect to the initial mean orbital plane. Ford et al. (2001) found that for two planets, most final orbits had inclinations <10◦ . We find higher inclinations for the three-planet case. Ford et al., with only two planets, had to start them on orbits with closer spacings in order to initiate crossing. Their encounters were therefore at lower velocities and resulted in a higher fraction of collisions. With three planets, evolution times are longer and there are more scattering events, allowing more evolution of inclinations. Figure 6 shows that the inner planet has a modal inclination of about 15–20◦ , but extending to as much as 120◦ . The outer planet tends to have a lower inclination, ranging up to 45◦ , since it usually carries most of the angular momentum. The escaping planet might however alter the balance of momentum between the two surviving planets. This has immediate consequences on the methods that derive the inclination of extrasolar planet orbits from the detection of a disk around the star (Trilling et al. 2000). According to our model, the planets may have moved out of the formation plane and have a significant inclination respect to the observed dust disk. In the absence of a dust disk, all information about the

(4)

where m f is the mass of the inner planet. From this relation, one can show that af 

mf m1 a1

+

m2 a2

+

m3 a3

,

(5)

or for equal planetary masses a f  (1/a1 + 1/a2 + 1/a3 )−1 . For our choices of a1 , a2 , a3 , this yields a f  0.45 AU and a f = 2.25 AU, in excellent agreement with the results of our integrations. We note that Eq. (5) does not depend on the planet/star mass ratio, which is very small since planets are significantly less massive than stars. Also, our integrations show no correlation between initial and final positions. That is, the evolution is so chaotic that (for equal masses) any of the three initial planets is equally likely to end up in the inner or outer orbit or to be ejected. Due to the slower variation of orbital energy at larger distances from the star, there is much less constraint on the distance of the outer survivor’s orbit; it can be many times the distance of the outermost initial orbit. In Fig. 7 we show histograms of the orbital elements of the outer survivor for the same set of systems shown in Fig. 6. The eccentricities cover the whole range from

FIG. 7. Histograms of orbital parameters for the outer surviving planet of the same cases shown in Fig. 6. Here the inclination is relative to the orbital plane of the inner planet.

576

MARZARI AND WEIDENSCHILLING

original orbital plane is lost and the (potentially) observable quantity is the relative inclination of the planetary orbits. This is shown in Fig. 7. The mutual inclination can be larger than either planet’s inclination to the reference plane, as their (vector) angular momenta tend to be anticorrelated. The difference in orbital planes also has consequences for determining planetary masses by radial velocities, for which the measured quantity is m sin(i). In systems with more than one planet, their relative masses cannot be simply determined by assuming that sin i is the same for all planets. There are ways to generate orbits closer to the star. One is to increase the mass of one planet. From Eq. (5) we see that the inner planet (m 1 , a1 ) has the most effect. The available energy is increased, and if the more massive planet is the one ejected after the chaotic evolution or pushed into the outer orbit, the less massive surviving inner planet can reach a smaller semimajor axis. If one planet is twice as massive as the others and initially in the innermost orbit, then Eq. (5) shows that one of the smaller planets can reach a f  1.6 AU if the larger one is ejected. However, on a statistical basis it is far more probable that one of the smaller planets is ejected. In Fig. 8 we show the distribution of orbital elements of the final inner planet for such a case. Most of the outcomes leave the massive planet in the inner orbit near 3 AU. In only 6 cases out of 64 is the innermost planet closer than 2 AU; however, in 5 of those 6, one of the lighter planets has moved to the inner position. Comparing Fig. 8 to Fig. 6 we notice the reduction of the eccentricity and inclination of the

planet injected in the inner orbit. In Fig. 6 the peak of the eccentricity distribution is at about 0.6 while in the case with a more massive planet (Fig. 8), the eccentricity distribution is peaked at 0.3. Also, the inclinations are lower in this second case. Another possibility to produce a planet closer to the star is a “double jump”; i.e., in some cases the two planets remaining after the first is ejected still can make close approaches and the encounters continue until a second planet is ejected. This process sometimes gives the escaping planet a positive energy, i.e., a hyperbolic excess velocity. In such cases, the semimajor axis of the remaining planet can be as small as about 1/4 of the initial distance of the inner planet. However, this outcome is very rare; two planets are ejected only in about 5% of the cases starting with three equal-mass planets. The probability of multiple ejections increases with the number of starting planets. For example, if a system initially comprising five planets between 5 and 10 AU evolved so that four were ejected or pushed into distant orbits, the innermost remaining planet could reach a semi-major axis near 1 AU. If we assume larger numbers of planets, we have also to assume that the initial protoplanetary disk was more massive than the solar nebula. The average value of disk masses are around ∼0.02M similar to that assumed for the solar nebula; however, there are disks with masses up to 0.5M (Beckwith et al. 1990) that could harbor a large number of planets. Of course, if the in situ formation of giant planets is not restricted to beyond the ice condensation boundary of the circumstellar disk, or if for some disks the “snow line” is significantly closer to the star, then the orbits of the planets were already close to the star, and gravitational scattering allows the final orbit to be closer (Fig. 5). 3. STATISTICS OF THE PLANETARY SYSTEMS

FIG. 8. Histograms of orbital parameters for the inner survivor, when the initially inner planet has a mass twice that of the others. In most of the cases, the more massive planet remains in the inner orbit; the shaded areas indicate cases where one of the less massive planets moved inward. The inclinations are relative to the initial reference plane.

The statistical distribution of the orbital elements shown in Fig. 6 or 8 cannot be directly compared with the presently known distribution of extrasolar planets. We currently observe systems with a random distribution of evolutionary lifetimes while the samples shown correspond to “snapshots” of systems just after one planet is ejected at the end of their chaotic evolution. The two remaining planets continue to interact by long-range perturbations, and in some cases may have large oscillations in eccentricity and inclination on time scales of the order of some 107 years. In Fig. 9 we show the evolution in time of the eccentricities and mutual inclination of two planets of a system after the chaotic phase. This case has long-term oscillations superimposed on short-period oscillations; both are due to the perturbations of the other survivor. The periapsis of the inner orbit is inside 0.05 AU for a small fraction of the time, during which the tidal interaction with the star can circularize the orbit. Figure 10 shows a more extreme example of long-term orbital evolution of two planets after ejection of the third. Here the amplitude of oscillations of mutual inclination gradually increases, as does the eccentricity of the inner planet. After 1.2 × 108 years, the inner planet collides with the star. We repeated the integration

THE JUMPING JUPITER MODEL

577

0.6

0.4 Outer planet

120

90 1

0.8 Inner planet 0.6 0 Time (yr)

FIG. 9. Time evolution of the eccentricities and mutual inclinations of the two surviving planets after the chaotic phase in one system. The semimajor axes, 76 (upper plot) and 2.3 (bottom plot) AU, do not vary significantly. The gravitational interactions between the two surviving bodies cause long-term oscillations of the eccentricity that can lead the periastron of the inner body close to the star so that tidal interactions may circularize the orbit. The dashed line corresponds to that value of eccentricity for which the planet has the periastron inside 0.05 AU.

with a much smaller time step (0.1 day) and the impact on the star is confirmed. General relativistic and stellar quadrupole moment contributions are likely to be negligible in this context and should not affect the overall evolution of the eccentricity of the planet itself. They can cause a slow perihelion shift that is only a fraction of that caused by the secular perturbations of the outer planet, and they do not affect the remaining orbital parameters. Such impact events are of particular interest, in light of the re0.8

0.6 Outer planet 90 60 30 1 0.8 0.6 Inner planet 0.4 0 Time (yr)

FIG. 10. Another example of long-term evolution of two survivors. The outer planet has a  36.9 AU, the inner has a  2.40 AU. Oscillations of the latter’s eccentricity cause impact onto the star (q < 0.005 AU after 1.2 × 108 years.

FIG. 11. As in Fig. 6, the histograms show the statistical distribution of the orbital elements of the inner planet with the difference that in this case the elements are sampled at random evolutionary times between 0 and 2 × 108 years after one planet has been ejected. In this way we account for the postejection orbital evolution, and we can compare our model results with observations. With respect to Fig. 6, two additional cases ended up with a collision and one with a “double jump.”

cent detection of 6 Li in HD82943, indicating that this star has ingested a planet (Israelian et al. 2001). To have a realistic estimate of the distribution of the orbital elements observable for real systems produced by the “Jumping Jupiter” mechanism, we have to further integrate the two-planet systems and extract at random times a set of orbital elements representative of the system. We continued the numerical integration of all the 64 systems for additional 2 × 108 years, and we selected a set of orbital elements for each system at a random time between 0 and 2 × 108 years. The more realistic histograms of the orbital element distributions obtained in this way are shown in Fig. 11. The distribution in eccentricity appears to be bimodal with a first peak near 0.3 and a second one at eccentricity values of about 0.6. This is presumably due to the amplitude of quasi-sinusoidal oscillations of eccentricity about the mean value, with planets generally spending more time near the extreme. Compared with the “snapshot” in Fig. 6, there is a greater proportion of small periastron distances. Gravitational scattering is not in itself able to produce orbits of the type of 51 Peg, with very small semimajor axes and eccentricities. However, it can produce high eccentricities with small periastron distances in a significant fraction of cases, augmented by the continuing perturbations of an outer planet. Our results in this context are consistent with those of Levison et al. (1998). Those planets with periastra close to the star will be affected by the tidal interaction with the outer envelope of the star

578

MARZARI AND WEIDENSCHILLING

planets are known with e ∼ 0.5 and q < 0.1 AU (HD108147 and HD74156). If tidal circularization is effective, then we must conclude that such systems are young and have not had time to be tidally altered. Alternatively, perturbations by a more distant massive planet may mean that such large e and small q are attained for only a small fraction of the time, prolonging their tidal evolution. It would be of great interest to determine the ages of their stars, as well as the existence and orbital parameters of any distant companions. 4. CONCLUSIONS

FIG. 12. Distribution in the a − e plane of the orbital elements of the inner planet at random times. The filled circles are systems in which the inner planet spends a significant amount of time within 0.05 AU from the star so that its orbit may eventually be circularized near the periastron distance.

and by tidal dissipation within the planets themselves. Tides raised on the planet by the star will decrease eccentricity; if the periastron distance is closer than the synchronous point, tides raised on the star by the planet will also decrease e, as well as a. These processes may circularize orbits near the periastron distance, producing “Hot Jupiters” in close orbits. Figure 12 shows values of eccentricity versus semimajor axis for the inner planet of these same 64 cases at the randomly selected times. The clustering of semimajor axes in the range 2.2–2.6 AU is an artifact of our starting conditions, with the inner planet initially at 5 AU. In real systems, there should be some variation of formation distances; also, planetary masses would not be equal as we have assumed. Both of these effects would spread the final distribution of semimajor axes. The distribution of eccentricities is fairly uniform up through e  0.8, in good agreement with the observed distribution. The black dots represent those orbits with minimum periastron distance less than 0.05 AU during oscillations of eccentricity after ejection of one planet and stabilization of the two that remain. The orbits of these planets might be circularized near the periastron distance by tidal dissipation. A quantitative model of this process would be difficult, as the effectiveness of tidal dissipation in stars and gas giant planets is poorly constrained (Lin et al. 2000) and the common expressions for the dissipation rate are not valid for large eccentricities. Also, tidal evolution of the inner orbit would itself alter the interaction with the outer perturbing planet, affecting the periastron distance. More work is needed to develop a quantitative model for circularization of their orbits. The observed distribution of e versus a (Fig. 1) shows an excess of planets in close orbits with small eccentricities (“hot Jupiters”), but there are still several with large eccentricities and small periastron distances, including HD80606 with the largest known value, e = 0.93, with q of only 0.03 AU! Two other

The “Jumping Jupiters” model has some degree of cosmogonical plausibility. Planetary formation has a significant stochastic element. There seems to be no requirement that planets should form only on orbits that are stable on time scales much longer than their formation times; that would imply some a priori knowledge of the final configuration. Oligarchic growth of protoplanets in the core-accretion scenario should lead to orbital instabilities after their masses are increased by capture of gas from the protostellar disk (Weidenschilling and Marzari 1996, Thommes et al. 1999). We have shown that there is a range of initial orbital spacings that allow sufficient time for planetary formation (>106 years), yet lead to the eventual onset of close encounters and orbital chaos, even after delays of 107 –108 years. Our results lead to the prediction that in most cases with a planet in a close orbit with significant eccentricity, there should be another planet of comparable mass in a distant orbit. Its contribution to variations of stellar radial velocity would be much smaller due to its distance, making it harder to detect. However, a greater obstacle to detection is the fact that its orbital period would be many times longer than that of the inner planet, requiring a correspondingly longer interval of observation. If the “Jumping Jupiters” model is correct, distant planets should eventually be detected around most stars with currently known planets. We note that such planets may have orbits with high inclinations relative to the inner ones; in estimating their masses (or mass ratios), one should not assume that all planets in such systems have the same value of sin i. Thus far we have, of necessity, examined only a very limited parameter space of the possible range of planetary masses, initial orbits, and planet/star mass ratios. Our conclusions must therefore be tentative. Dynamical evolution of orbits by gravitational scattering can produce eccentricities of the magnitudes observed for extrasolar planets. Our results imply that the distribution of eccentricities should be fairly uniform, which appears to be consistent with observations. However, it is necessary to invoke tidal dissipation in order to explain the excess of planets in orbits of low eccentricity very close to their stars. Also, the number of observed orbits with very high eccentricities may be affected both by long-term variations due to perturbations by more distant planets and by tidal evolution of planets with small periastron distances. The observed distribution of semimajor axes presents more of a problem. The present results imply that

THE JUMPING JUPITER MODEL

it is hard for scattering to produce a semimajor axis less than about half the initial distance of the innermost planet that formed in a given system. It may be necessary to invoke formation of massive planets significantly closer to their stars than Jupiter’s distance from the Sun, and/or some additional migration mechanism, such as tidal coupling to the protostellar disk, that allows them to move inward from the nebular “snow line” after their formation, before the onset of orbital instability. It would be premature to conclude that there must be only one mechanism for planetary migration. Most of the extrasolar planets detected to date are more massive than Jupiter. Gravitational scattering implies that multiple planets of comparable mass formed about a given star. This would imply circumstellar disks that were considerably more massive than the solar nebula and/or had a higher efficiency of planet formation (possibly due to high abundance of heavy elements). However, we do have the example of the υ Andromeda system with three massive planets, so such systems undoubtedly are possible. It is of interest that the υ And system is close to instability (Rivera and Lissauer 2000, Stepinski et al. 2000, Barnes and Quinn 2001), and the possibility exists that chaotic evolution lies in its future. It should be remembered that only a small fraction of stars examined in radial velocity surveys have shown evidence of such massive planets. As planetary detection becomes more sensitive and is extended to lower masses, it will be interesting to see whether less massive systems tend to be more stable or also show evidence of chaotic orbital evolution. ACKNOWLEDGMENTS We thank E. B. Ford and an anonymous referee for their helpful comments. SJW’s research was supported by NASA’s Planetary Geology and Geophysics Program. This is PSI Contribution 358.

REFERENCES Barnes, R., and T. Quinn 2001. A statistical examination of the short-term stability of the υ Andromedae planetary system. Astrophys. J. 550, 884– 889. Beckwith, S., A. Sargent, R. Chini, and R. Guesten 1990. A survey for circumstellar disks around young stellar objects. Astron. J. 99, 924–945. Bodenheimer, P., O. Hubickyj, and J. J. Lissauer 2000. Models of the in situ formation of detected extrasolar giant planets. Icarus 143, 2–14. Boss, A. P. 1997. Giant planet formation by gravitational instability. Science 276, 1836–1839. Chambers, J., and G. W. Wetherill 1998. Making the terrestrial planets: N-body integrations of planetary embryos in three dimensions. Icarus 136, 304– 327. Chambers, J., G. W. Wetherill, and A. P. Boss 1996. The stability of multi-planet systems. Icarus 119, 261–268. Charbonneau, D., T. M. Brown, D. W. Latham, and M. Mayor 1999. Detection of planetary transits across a Sun-like star. Astrophys. J. 529 L45–L48.

579

Everhart, E. 1985. An efficient integrator that uses Gauss–Radau spacings. In Dynamics of Comets: Their Origin and Evolution (A. Carusi and G. B. Valsecchi, Eds.), pp. 185–202. Farinella, P. 1980. Irregular extrasolar systems and dynamical instability. Moon Planets 22, 25–29. Fischer, D., G. Marcy, P. Butler, S. Vogt, S. Frink, and K. Apps 2001. Planetary companions to HD 12661, HD 92788, HD 38529 and variations in Keplerian residuals of extrasolar planets. Astrophys. J. 551, 1107–1118. Ford, E. B., M. Havlickova, and F. A. Rasio 2001. Dynamical instabilities in extrasolar planetary systems containing two giant planets. Icarus 150, 303–313. Gladman, B. 1993. Dynamics of systems of two close planets. Icarus 106, 247–263. Hasegawa, M., and K. Nakazawa 1990. Distant encounters between Keplerian particles. Astron. Astrophys. 227, 619–627. Israelian, G., N. C., Santos, M. Mayor, and R. Rebolo 2001. Evidence for planet engulfment by the star HD82943. Nature 4311, 163–166. Ito, T., and S. Tanikawa 1999. Stability and instability of the terrestrial protoplanet system and their possible roles in the final stage of planet formation. Icarus 139, 336–349. Kokubo, E., and S. Ida 1998. Oligarchic growth of protoplanets. Icarus 131, 171–178. Levison, H., and Duncan, M. J. 1994. The long-term dynamical behavior of short-period comets. Icarus 108, 18–36. Levison, H., J. J. Lissauer, and M. J. Duncan 1998. Modeling the diversity of outer planetary systems. Astron. J. 116, 1998–2014. Lin, D. N. C., and S. Ida 1997. On the origin of massive eccentric planets. Astrophys. J. 477, 781–791. Lin, D. N. C., P. Bodenheimer, and D. Richardson 1996. Orbital migration of the planetary companion of 51 Pegasi to its present location. Nature 380, 38–40. Lin, D. N. C., J. Papaloizou, C. Terquem, G. Bryden, and S. Ida 2000. Orbital evolution and planet-star tidal interaction. In Protostars and Planets IV (V. Mannings, A. P. Boss, and S. S. Russell, Eds.), pp. 1111–1134. Univ. of Arizona Press, Tucson. Marcy, G. W., W. D. Cochran, and M. Mayor 2000. Extrasolar planets around main-sequence stars. In Protostars and Planets IV (V. Mannings, A. P. Boss, and S. S. Russell, Eds.), pp. 1285–1312. Univ. of Arizona Press, Tucson. Pollack, J., O. Hubickyj, P. Bodenheimer, J. Lissauer, M. Podolak, and Y. Greenzweig 1996. Formation of the giant planets by concurrent accretion of solids and gas. Icarus 124, 62–85. Rasio, F., and E. Ford 1996. Dynamical instabilities and the formation of extrasolar planetary systems. Science 274, 954–956. Rivera, E., and J. Lissauer 2000. Stability analysis of the planetary system orbiting υ Andromedae. Astrophys. J. 530, 454–463. Stepinski, T. F., R. Malhotra, and D. C. Black 2000. The υ Andromedae system: Models and stability. Astrophys. J. 545, 1044–1057. Thommes, E. W., M. J. Duncan, and H. F. Levison 1999. The formation of Uranus and Neptune in the Jupiter-Saturn region of the Solar System. Nature 402, 635–638. Trilling, D., R. H. Brown, and A. Rivkin 2000. Circumstellar dust disks around stars with known planetary companions. Astrophys. J. 529, 499–505. Ward, W. R., and J. Hahn 2000. Disk-planet interactions and the formation of planetary systems. In Protostars and Planets IV (V. Mannings, A. P. Boss and S. S. Russell, Eds.), pp. 1135–1155. Univ. of Arizona Press, Tucson. Weidenschilling, S. J., and F. Marzari 1996. Gravitational scattering as a possible origin for giant planets at small stellar distances. Nature 384, 619–621.