Origin and Dynamics of the Mutually Inclined Orbits of Upsilon Andromedae c and d
Rory Barnes, Richard Greenberg, Thomas R. Quinn, Barbara E. McArthur, G. Fritz Benedict
aa r X i v : . [ a s t r o - ph . E P ] N ov Origin and Dynamics of the Mutually Inclined Orbits of υ Andromedae c and d
Rory Barnes , , Richard Greenberg , Thomas R. Quinn , Barbara E. McArthur , andG. Fritz Benedict ABSTRACT
We evaluate the orbital evolution and several plausible origins scenarios forthe mutually inclined orbits of υ And c and d. These two planets have orbital ele-ments that oscillate with large amplitudes and lie close to the stability boundary.This configuration, and in particular the observed mutual inclination, demandsan explanation. The planetary system may be influenced by a nearby low-massstar, υ And B, which could perturb the planetary orbits, but we find it cannotmodify two coplanar orbits into the observed mutual inclination of 30 ◦ . How-ever, it could incite ejections or collisions between planetary companions thatsubsequently raise the mutual inclination to > ◦ . Our simulated systems withlarge mutual inclinations tend to be further from the stability boundary than υ And, but we are able to produce similar systems. We conclude that scatteringis a plausible mechanism to explain the observed orbits of υ And c and d, butwe cannot determine whether the scattering was caused by instabilities amongthe planets themselves or by perturbations from υ And B. We also develop aprocedure to quantitatively compare numerous properties of the observed systemto our numerical models. Although we only implement this procedure to υ And,it may be applied to any exoplanetary system.
1. Introduction
The υ Andromedae ( υ And) planetary system is the first multiple planetary systemdiscovered beyond our own Solar System around a solar-like star (Butler et al. Department of Astronomy, University of Washington, Seattle, WA, 98195-1580 Virtual Planetary Laboratory Lunar and Planetary Laboratory, University of Arizona, Tucson, AZ 85721 Department of Astronomy, University of Texas at Austin, TX 78712 e.g.
Laughlin & Adams 1999; Rivera & Lissauer2000; Lissauer & Rivera 2001; Laskar 2000; Barnes & Quinn 2001; Go´zdziewski et al. e.g.
Stepinski et al. et al. ◦ mutual inclination, via astrometry(McArthur et al. et al. (2010), which place strict constraints on the system’s origin.Such large mutual inclinations among planets are unknown in our Solar System andhence indicate different processes occurred during or after the planet formation process. Weassume υ And c and d formed in coplanar, low eccentricity orbits in a standard planetaryformation models (see e.g.
Hubickyj 2010; Mayer 2010), but then additional phenomena,occurring late or after the formation process, altered the system’s architecture. We considertwo plausible mechanisms: a) If the distant stellar companion υ And B (Lowrance et al. et al. et al. et al. υ And represents the first exoplanetary system with fullthree-dimensional orbits and true masses directly measured (aside from the pulsar systemPSR 1257+12 [Wolszczan 1994]), we may exploit this information when evaluating formationmodels. To that end, we develop a simple metric that quantifies the success of a model atreproducing numerous observed aspects of the υ And planetary system. Although we applythis approach to υ And, it is generalizable to any planetary system.In this investigation we limit the analysis to just υ And c and d, and to the stablefit presented in McArthur et al. (2010). As described in that paper, the inclination andlongitude of ascending node of b are not detectable with the
Hubble Space Telescope , so 3 –rather than explore the range of architecture permitted by RV data, we focus on the knownproperties of c and d. Furthermore, we do not consider the range of uncertainties in c and das the stable fit in McArthur et al. (2010) is surrounded by unstable fits. As we see below,even limiting our scope in this way, we still must perform a large number of simulations overa wide range of parameter space.We first ( §
2) analyze the current orbital oscillations of υ And c and d with an N-bodysimulation. We then use the dynamical properties to constrain our two inclination-raisingscenarios, which we explore through > § υ And B by itself cannot raise the mutual inclinations to 30 ◦ . In § ◦ mutual inclination and extreme proximity to the stability boundary. In § υ And c and dPlanet m (M J ) a (AU) e i ( ◦ ) ̟ ( ◦ ) Ω ( ◦ ) n ( ◦ )c 14.57 0.861 0.24 16.7 290 295 270.5d 10.19 2.70 0.274 13.5 240.8 115 266.1
2. Orbital Evolution of υ And c and d
In this section we determine the orbital evolution of υ And c and d. As the mass andorbit of planet b are unknown and the four-body interactions between υ And A, b, c andd are extremely complex and depend sensitively on a secular resonance, general relativityand the stellar oblateness (McArthur et al. . For this integration we change the coordinatesystem from that in the discovery paper, which is based on the viewing geometry, andinstead reference our coordinates to the invariable (or fundamental) plane. This plane isperpendicular to the total angular momentum vector of the system (although we ignoreplanetary spins), i.e. we rotated the coordinate system. The planets’ orbital elements in thiscoordinate system are listed in Table 1 at epoch JD 2452274.0. The system shown in Table1 was found to be stable in McArthur et al. , but it was also noted that the system is closeto the stability boundary. Therefore we cannot exclude the possibility that other solutions 4 –may result in significantly different dynamical behavior. Here and below we assume thatsuch a situation is not the case.The orbits of planets c and d undergo mutual perturbations which cause periodic varia-tions in orbital elements over thousands of orbits. The long-term changes can be convenientlydivided into two parts: the apsidal evolution (changes in eccentricity e and longitude of peri-astron ̟ ) and the nodal evolution (changes in inclination i , longitude of ascending node Ω,and hence the mutual inclination Ψ). The variations, starting with the conditions in Table 1,are shown in Fig. 1. In this figure, we chose a Jacobi coordinate system in order to minimizefrequencies due to the reflex motion of the star.In Fig. 1 the left panels show the apsidal behavior, and the right show the nodal behav-ior. The lines of apse oscillate about ∆ ̟ = π , i.e. anti-aligned major axes. This revision onceagain changes the expected apsidal evolution. Initially Chiang & Murray (2002) and Malho-tra (2002) found the major axes oscillated about alignment. Then Ford et al. (2005; see alsoBarnes & Greenberg 2006a,c) found the system was better described as “near-separatrix,”meaning the apsides lie close to the boundary between libration and circulation. Now wefind that the system found by McArthur et al. (2010) librates in an anti-aligned sense! Sub-stantial research has examined the secular behavior of exoplanetary systems, yet the storyof υ And shows that predicting the dynamical evolution of a planetary system based onminimum masses and poorly constrained eccentricities is uncertain at best and foolhardy atworst. Even now, without full three-dimensional information about planets b and e (a trendseen in McArthur et al. [2010]) our analysis should only be considered preliminary.Table 2: Dynamical Properties of υ And c and dProperty Value e minc e mind e maxc e maxd ǫ i minc . ◦ i mind . ◦ i maxc . ◦ i maxd . ◦ Ψ min . ◦ Ψ max . ◦ β/β crit ∼ e by more than 0.01). Note also that the ratio of the orbitalperiod, 5.32, is not very close to the low-order 5:1 and 11:2 resonances, so this high frequencyevolution is not due to a mean motion resonance.The numerical integration also allows a calculation of the proximity to the apsidalseparatrix ǫ (Barnes & Greenberg 2006c). When ǫ is small ( < ∼ . υ And was the firstsystem to be identified as near-separatrix (Ford et al. ǫ = 0 .
17 indicating thatthe system is actually not close to the separatrix, according to the McArthur et al. (2010)model.Also of interest is the system’s proximity to dynamical instability, e.g. the ejection ofa planet. Previous studies found planets c and d are close to this limit (Rivera & Lissauer2000; Barnes & Greenberg 2001, 2004; Go´zdziewski et al. β/β crit (Barnes & Greenberg2006b,2007b; see also Marchal & Bozis 1982; Gladman 1993; Veras & Armitage 2004). If β/β crit > <
1, it is unstable. We find β/β crit = 1 .
075 andtherefore these two planets are very close to the dynamical stability boundary. We cautionthat constraints based on β/β crit may be misleading, as Hill stability is strictly only applicableto three-body systems.In Table 2 we list some statistics of our 10 year integration based on 36.5 day outputintervals in astrocentric coordinates. Superscripts min and max refer to the minimum andmaximum values achieved, respectively. Clearly the actual dynamics of this system dependon the presence and properties of planet b, and, in principle, any additional unconfirmedplanets, but without better data on these objects, Table 2 is the best available characteri-zation of the orbital evolution of these two planets. They may also be used to evaluate theorigins scenarios described in the following two sections.
3. Perturbations from υ And B υ And B is a distant M4.5, 0.2 M ⊙ companion star to υ And A (Lowrance et al. et al. et al. [2002]; Raghavan et al. [2006]). Estimates for their 6 –Fig. 1.—
Secular evolution of υ And c (black) and d (red) without planet b.
Top left:
Evolutionof ∆ ̟ . Bottom left:
Evolution of the eccentricities. Black is planet c, red d.
Top right:
Evolutionof Ψ.
Bottom right:
Evolution of the inclinations. Black is planet c, red d. et al. et al. υ And B havepumped up the mutual inclinations of c and d?Given the uncertainty in B’s orbit, we consider a broad parameter space sweep: 13,200simulations that cover the range 500 ≤ a B ≤ . ≤ e B ≤ .
85 and 30 ◦ ≤ i B ≤ ◦ .Here i B is referenced to the initial orbital plane of the planets, not the invariable plane ofthe four-body system. The angular elements were varied uniformly from 0 to 2 π . Notethat this range was chosen to increase the perturbative effects of B and does not reflect anyexpectation of its actual orbit. Each simulation was run for 5 × years, which correspondsto ∼
250 orbits of B. While not a long time, we find that B may destabilize the planetarysystem, which could lead to large mutual inclinations of planet c and d (see § < ◦ . However, 3 simulations ejected planet d; 192 led to planets with Ψ max > ◦ ; 26 withΨ max > ◦ ; and 1 case out of the 13,200 reached Ψ max = 34 ◦ . The 192 non-planar caseswere spread throughout parameter space, with no significant clustering.To explore effects on longer timescales, we integrated 20 cases to 1 Gyr. Four of theprevious simulations that had led to significant mutual inclinations (including that whichled to Ψ max = 34 ◦ ) were tested, in order to examine stability. The other 16 were chosen fromamong those in which Ψ max stayed < ◦ over 250 orbits of B, in order to determine if mutualinclinations could develop over longer timescales. We find that most of these simulations, infact, ejected a planet. Therefore the orbit of υ And B appears able to destabilize a circular,coplanar system.Next we relax the requirement that the planets began on coplanar orbits and ran sim-ulations with initial mutual inclinations Ψ = 3 ◦ , 10 ◦ , and 30 ◦ , and with a B = 700, e B = 0, i B = 30 ◦ . The two planets began with their current best fit semi-major axes and masses,but on circular orbits, and one inclination was set to 3 ◦ , 10 ◦ or 30 ◦ with the masses heldconstant. These systems were integrated for 1 Gyr. In each of these cases, shown in Fig. 2,the initial value of Ψ is maintained for the duration of the simulation. The widths of thelibration increase with Ψ because the interactions among the planets are driving a secularoscillation. It therefore seems that even if the planets began with a nonzero relative incli-nation, υ And B is unable to pump it to the range shown in Fig. 1. These simulations alsodemonstrate that plausible orbits of υ And B will not destabilize the observed system.These simulations indicate that it is unlikely that υ And B could have twisted theorbits into the mutually inclined system we see today. However, it could have destabilizedthe system, which we will see in the next section is a process which can lift coplanar orbits 8 –to Ψ > ∼ ◦ .
4. Planet-Planet Scattering
Our second hypothesis considers the possibility that the planetary system formed inan unstable configuration independent of B, and that encounters between planets ultimatelyejected an original companion leaving a system with high mutual inclinations. Such impulsiveinteractions could drive inclinations to large values, perhaps as large as 60 ◦ (Marzari &Weidenschilling 2002; Chatterjee et al. et al. e.g. one or two additional planets with massesin the range 1 – 15 MJup. At the end of this section we summarize the results of all thesesimulations, but initially we focus on one subset.We completed 5,000 simulations that began with three 10 – 15 MJup mass objects(uniformly distributed in mass) separated by 4–5 mutual Hill radii (Chambers et al. e < . i < ◦ and 0 . < a < years withMercury’s hybrid integrator, conserving energy to 1 part in 10 . About 1% of cases failedto conserve energy at this level and were thrown out. These ranges are somewhat arbitrarybut follow the recent study by Raymond et al. (2010), which considered smaller mass planetsat larger distances. They found that a system consisting of three 3 MJup planets could,after removal of one planet and settling into a stable configuration, end up with Ψ > ◦ about 15% of the time (down from 30% for three Neptune-mass planets). However, theyalso found that only 5% of systems of three 3 MJup planets settled into a configurationwith β/β crit < .
1. Therefore we expect that these two parameters will be the hardest toreproduce via scattering. As we see below, this expectation is borne out by our modeling.In our study, a successful model conserved energy adequately (1 part in 10 ), removedthe extra planet, and the remaining planets all had orbits with a <
10 AU. 2072 trialsmet these requirements (1416 collisions and 656 ejections). We ran each of these final two-planet configurations for an additional 10 years (again validating the simulation via energyconservation) to assess secular behavior for comparison with the system presented in Table2. In Fig. 3 we show the outcome of one such trial in which a hypothetical planet wasejected. The format of this figure is the same as Fig. 1. The behavior is qualitatively similaras in Fig. 1, including anti-aligned libration of the apses, the magnitudes of the eccentricitiesand the inclinations, and the short period oscillation superposed on the longer oscillation.The mutual inclination for this case is even larger than the observed system. This system’s 9 –Fig. 2.— Variation of Ψ for four hypothetical systems with υ And B included (see text for moredetails). The black points are systems in which initially i c = 3 ◦ , blue i c = 10 ◦ , red i c = 30 ◦ , andgreen i d = 30 ◦ .
10 – β/β crit is 1.06, slightly lower than the observed system. This simulation, which is one of theclosest matches to the observed system, shows that the ejection of a single additional planetcould have produced the υ And system.Figure 3 is but one outcome. We next explore the statistics of this suite of simulationsand consider the other orbital elements and dynamical properties. We divide the outcomesinto two cases: Ejections and Collisions. These two phenomena could produce significantlydifferent outcomes. For example, collisions tend to occur near periastron of one planet andapoastron of the other, and we might expect the merged body to have a lower eccentricitythan either of the progenitors. We show the cumulative distributions of the properties listedin Table 2 in Fig. 4. Comparing the values of orbital elements at a given time is not ideal,but as it has been done many times (see e.g.
Ford et al. et al. et al. e , i , and Ψ at the end of the initial 10 year integration.In panels d–h we show the ranges of i min , i max , Ψ min and Ψ max . We find that 8.9%of successful models produced a system with Ψ max > ◦ , consistent with Raymond etal. (2010). Note that ejections produce Ψ max > ◦ about 20% of the time.In Fig. 4i we show the ǫ distribution, which is bimodal with one peak near 0.1 andanother near 10 − . The observed value of 0.17 is not an unusual value, and we find thatsystems with this ǫ value can have appropriate values of e min and e max . We note that thesignificant fraction of systems near the apsidal separatrix contradicts the results of Barnes &Greenberg (2007a), which found that scattering only produced ǫ < .
01 a few percent of thetime. The most likely explanation for this difference is that Barnes & Greenberg consideredcoplanar orbits and forbade collisions, whereas here we explore non-planar motion. Ourresults also indicate that near-separatrix motion is likely a result of collisions, rather thanejections, and ǫ < − (which is unlikely to be measured any time soon) only result fromcollisions.Figure 4j shows the distribution of β/β crit after scattering. Here the difference betweencollisions and ejections is starkest: Ejections have a much broader distribution than collisions.Barnes et al. (2008) noted that systems are “packed” (no additional planets can lie in betweentwo planets) when β/β crit < ∼
2, which is near the peak of the ejection distribution. Our resultsconsistent with Raymond et al. (2009).In this model the mutual inclinations and proximity to instability are the strongestconstraints on the system’s origins, but we may quantify scattering’s ability to reproduce allthe observed features of the system. Most previous studies of scattering have focused on re-producing eccentricity distributions, but all available information should be used. With this 11 –Fig. 3.—
Secular evolution of a simulated system after the ejection of an exterior planet.
Top left:
Evolution of ∆ ̟ . Bottom left:
Evolution of the eccentricities. Black is planet c, red d.
Top right:
Evolution of Ψ.
Bottom right:
Evolution of the inclinations.
12 –goal in mind, we lay out here a simple method to quantify any model’s ability to reproduceobservations. For each simulated system, we calculate its “parameter space distance” ρ fromthe best fit. We define this quantity as ρ = rX (cid:16) η − η j η (cid:17) , (1)where η represents e minj ...β/β crit and j = c,d (see Table 2). This statistic has several limi-tations: It ignores correlations between parameters; is dependent on the coordinate system;ignores uncertainties in the observations (as discussed above), and possibly overweights someparameters by including combinations of variables that are not independent. Although crude, ρ does provide a quantitative estimate of how close a modeled system is to the observed sys-tem. Smaller values of ρ signal a system that is a closer match to the actual system.In Fig. 4k we plot the distributions of ρ . These distributions resemble the β/β crit distributions (panel j), suggesting it is the most important constraint on the system. InTable 3 we show some statistics of our runs, where “min” is the minimum value of the set,“avg” is the mean, σ is the standard deviation, and “max” is the maximum value.For most of the parameters we consider, ejections and collisions do a reasonable jobof producing the observed values. However, this representation does not show any cross-correlation, i.e. does a system with high Ψ also have low β/β crit ? We explore this relationshipin Fig. 5. We see that post-collision systems (blue points) cluster heavily at low β/β crit andΨ max , but post-ejection systems (red points) have a much broader range. Nonetheless, thetwo outcomes seem equally likely to reproduce the system, represented by the “+” (recallthat there are three times as many blue points as red). However, from our models the actualprobability that instabilities can reproduce the υ And system is less than 1%.From our analysis of these 2072 systems, we see that it is possible for planetary ejectionsand collisions to reproduce the observed υ And system, albeit with low probability. Thissuite of simulations is obviously limited in scope, so we performed 36,000 more simulationsrelaxing constraints on planetary mass (allowing uniform values between 1 and 15 MJup),separation (uniform distribution between 2 and 5 mutual Hill radii), and number of planets (3or 4). These other simulations began with two planets with approximately the same semi-major axes as observed, and placed planets interior and/or exterior to these two planets.These simulations show that the equal-mass case we considered here is the best method toachieve large Ψ, as expected from Raymond et al. (2010). For additional planets with massesless than 5 MJup, Ψ values greater than 30 ◦ are very unlikely, but β/β crit ∼ υ And c and d. 13 –Fig. 4.— Cumulative distributions of properties (out of 626 ejections, 1416 collisions) afterscattering of one planet. Black curves refer to cases in which one planet was ejected, red toa collision between two planets. Thick vertical lines correspond to the value for the best fitsystem. If the parameter is measurable for both planets, solid lines indicate planet c, dashedd. 14 –Fig. 5.—
Relationship between proximity to instability ( β/β crit ) and maximum mutual inclination(Ψ max ). Systems which experienced a collision are color-coded blue, ejections red. The best fit tothe υ And system is shown by the + sign.
15 –Table 3: Distribution of Properties After Collision/EjectionProperty Ejection Collisionmin avg σ max min avg σ max e minc . × − × − e mind . × − × − e maxc e maxd ǫ − . × − i minc ( ◦ ) 0.02 7.2 7.50 34.4 0.013 1.2 2.7 30.4 i mind ( ◦ ) 0.002 4.1 5.9 43.7 0.005 0.77 2.1 30.4 i maxc ( ◦ ) 0.08 11.9 10.8 50.8 0.014 1.7 3.8 43.0 i maxd ( ◦ ) 0.03 8.3 8.7 48.7 0.007 1.26 3.4 49.7Ψ min ( ◦ ) 0.06 11.7 9.7 46.8 0.036 2.0 4.5 39.8Ψ max ( ◦ ) 0.1 19.8 15.5 63.9 0.032 2.9 6.7 62.1 β/β crit ρ
5. Conclusions
We have shown that the orbital behavior of the model for Ups And c and d proposed byMcArthur et al. (2010) is quite different from that of orbital models identified by previousstudies that had no knowledge of the inclination or actual mass of the planets. The majoraxes librate in an anti-aligned configuration, and their mutual inclination is substantial andoscillates with an amplitude of about 10 ◦ .We find that the companion star υ And B by itself cannot pump the mutual inclinationup to large values, even if the planets began with a significant relative inclination. However,it may have sculpted the planetary system by inciting an instability that ultimately led toejections of formerly bound planets. The timescale to develop these instabilities is long.The configurations of B that could have such effects are sparsely distributed over parameterspace, and the orbits of previously bound planets cannot be specified. These factors makethe role of υ And B complicated, but suggest an in-depth analysis of its role merits furtherresearch.Even without B, planet-planet scattering may have driven the system to the observedstate. That process can easily reproduce the apsidal motion, but pumping the mutualinclination up to the observed values is difficult and probably requires the removal of a planetwith mass > β/β crit . Collisions may leave a system near that boundary,whereas ejections tend to spread out the planets. Furthermore, we find that collisions tendto produce systems with low β/β crit and low Ψ, while ejections produce a broad range ofΨ, but large values of β/β crit . Nonetheless, Fig. 3 demonstrates that scattering can producesystems similar to υ And.Although scattering is a reasonable process to produce the observed architecture, wecannot determine the triggering mechanism. Did scattering occur because υ And B desta-bilized the planetary system? Or did the planet formation process itself, independent ofB, ultimately lead to instabilities? The presence of υ And B makes distinguishing thesepossibilities very difficult. A larger census of mutual inclinations and stellar companions canresolve this open issue.Alternatively, our decisions about the system at the onset of scattering could be mis-taken. We assumed the planets formed inside the original protoplanetary disk with inclina-tions < ◦ . It may be that larger initial inclinations are possible prior to scattering, in whichcase the planets could be pumped to larger mutual inclinations (Chatterjee et al. et al. i ’s andΩ’s of b and the subsequent orbital evolution of the entire system.Figure 4i shows that scattering tends to produce two types of apsidal behavior: near-separatrix ( ǫ < − ) and motion far from the separatrix ( ǫ > .
01) with a desert in between.Adding a second scatterer to the mix does not erase this bimodality. This result contrastswith the case with no inclinations (Barnes & Greenberg 2007a), in which near-separatrixmotion ( ǫ < .
01) is not a common outcome of scattering. Fig. 4i shows that, in fact, bothoutcomes are likely, at least in systems similar to υ And. Although there are hints of thisstructure in the observed exoplanet population (Barnes & Greenberg 2006c), those resultsare based on radial velocity data. We now know that minimum masses are not necessarily agood indicator of apsidal motion. 17 –For υ And, the large mutual inclination and proximity to instability are strong con-straints on the origin of its planetary system. However, for other systems, this may not bethe case. We outlined in § ρ which provides a statistic for quantitatively comparingmodels. In our analysis we ignored the observational errors, which is regrettable, but neces-sary due to the system’s extreme proximity to dynamical instability (McArthur et al. et al. to find ρ values less that those listed in Table 3, as lower values imply a closer match to the system,assuming that other stable solutions show similar behavior to the one we describe here. Fur-thermore, investigations into exoplanet formation could compare distributions of observedand simulated properties as a quantitative method for model validation.The revisions of McArthur et al. (2010) reveal the importance of the mass-inclinationdegeneracy in dynamical studies of exoplanets. Clearly in some cases masses can be muchlarger than the minimum value measured by radial velocity, which in turn changes secularfrequencies and eccentricity amplitudes. However, large changes in mass due to the mass-inclination degeneracy should be rare, hence, trends using minimum masses may still bevalid. Nonetheless, we urge caution when exploring trends among dynamical properties( e.g. Zhou & Sun 2003; Barnes & Greenberg 2006c), as they may be misleading.Even if 30 ◦ mutual inclinations turn out to be rare, systems with Ψ ∼ ◦ probably arenot (Fig. 4; Chatterjee et al. et al. et al. et al. υ And warns us that habitabilityassessment hinge on the orbital architecture of the entire planetary system.RB acknowledges support from the NASA Astrobiology Institute’s Virtual PlanetaryLaboratory lead team, supported by cooperative agreement No. NNH05ZDA001C. RG ac-knowledges support from NASA’s Planetary Geology and Geophysics program, grant No.NNG05GH65G. BEM and GFB acknowledge support from NASA through grants GO-09971,GO-10103, and GO-11210 from the Space Telescope Science Institute, which is operated bythe Association of Universities for Research in Astronomy (AURA), Inc., under NASA con- 18 –tract NAS5-26555. We also thank Sean Raymond for helpful discussions.
REFERENCES
Armstrong, J.C., Leovy, C.B., & Quinn, T.R. 2004, Icarus, 171, 255Atobe, K., Ida, S. & Ito, T. 2004, Icarus, 168, 223Atobe, K., Ida, S. 2007, Icarus, 188, 1Barnes, R., Go´zdziewski, K., & Raymond, S.N. 2008, ApJ, 680, L57Barnes, R. & Greenberg, R. 2006a, ApJ, 652, L53————. 2006b, ApJ, 638, 478————. 2006c, ApJ, 647, L163————. 2007a, ApJ, 665, L67————. 2007b, ApJ, 659, L53Barnes, R. & Quinn, T.R. 2001, ApJ, 554, 884————. 2004, ApJ, 611, 494Barnes, R. & Raymond, S.N. 2004 ApJ, 617, 569Butler, R.P. et al. et al. et al.
Formation and Evolution of Exoplanets , Rory Barnes (ed), Wiley-VCH, Berlin.Juri´c, M. & Tremaine, S. 2008, ApJ, 686, 603Kozai, Y. 1962, AJ, 67, 591Laskar, J. 2000, PhRvL, 84, 3240Laskar, J., Joutel, F., & Robutel, P. 1993, Nature, 615Laughlin, G. & Adams, F.C. 1999, ApJ, 526, 881Lissauer, J.J. & Rivera, E.J. 2001, ApJ, 554, 1141Lowrance, P.J., Kirkpatrick, J.D., & Beichman, C.A. 2002, ApJ, 572, L79Malhotra, R. 2002, ApJ, 575, L33Marchal, C. & Bozis, G. 1982, CeMDA, 26, 311Marzari, F. & Weidenschilling, S. 2002, Icarus, 156, 570Mayer, L. 2010 in
Formation and Evolution of Exoplanets , Rory Barnes (ed), Wiley-VCH, 19 –Berlin.McArthur, B. et al.
Solar System Dynamics , Cambridge UP, CambridgePatience, J. et al. et al.et al.