The (In)Stability of Planetary Systems

1 MB·English

by  Rory Barnes

DRAFTVERSIONFEBRUARY5,2008
PreprinttypesetusingLATEXstyleemulateapjv.04/03/99
THE(IN)STABILITYOFPLANETARYSYSTEMS
RORYBARNESANDTHOMASQUINN
Dept.ofAstronomy,UniversityofWashington,Box351580,Seattle,WA,98195-1580
E-mail:[email protected],[email protected]
DraftversionFebruary5,2008
ABSTRACT
4 Wepresentresultsofnumericalsimulationswhichexaminethedynamicalstabilityofknownplanetarysystems,
0 astarwithtwoormoreplanets. Firstwevarytheinitialconditionsofeachsystembasedonobservationaldata.
0 We then determine regions of phase space which produce stable planetary configurations. For each system we
2 perform1000∼106yearintegrations.WeexamineυAnd,HD83443,GJ876,HD82943,47UMa,HD168443,and

n thesolarsystem. Wefindthattheresonantsystems,2planetsinafirstordermeanmotionresonance,(HD82943
a and GJ876) have very narrow zones of stability. The interacting systems, not in first order resonance, but able
J toperturbeachother(υ And,47UMa,andSS)havebroadregionsofstability. Theseparatedsystems,2planets

9 beyond10:1resonance,(weonlyexamineHD83443andHD168443)arefullystable.Furthermorewefindthatthe
bestfitstotheinteractingandresonantsystemsplacethemveryclosetounstableregions.Theboundaryinphase

1 spacebetweenstabilityandinstabilitydependsstronglyontheeccentricities,and(ifapplicable)theproximityof
v the system to perfectresonance. In additionto 106 yearintegrations, we also examinedstability on ∼108 year
1
timescales. Foreachsystemweran∼10longtermsimulations,andfindthattheKeplerianfitstothesesystems
7
allcontainconfigurationswhichmayberegularonthistimescale.
1
1
0 1. INTRODUCTION systemsin whichthe (detected)planets’orbitsare beyondthe
4 10:1resonance. Theseplanetsaremostlikelydynamicallyde-
BySeptember2003,13planetarysystemshadbeendiscov-
0 coupled,howeversomeofthesesystemswarrantinvestigation.
ered(includingourownSolarSystem). Theextra-solarplane-
/ We examine planetary systems on two differenttimescales.
h tary systems(ESPS)are, possiblydueto observationalbiases,
p markedlydifferentfromourowninseveralways. Inoursolar First,weexploreparameterspacein106 yearintegrations. For
- thesesimulationswevaryinitialconditionstodeterminestable
o system(SS)theJovianmassplanetsallorbitatdistanceslarger
regionswithintheobservederrors.Second,wecontinueseveral
str tohtahner5hAaUn,da,ncdonotnanineagrilaynctirpcluanlaertsoribnitas(wei∼<de0r.a0n5g).eEoSfPdSis,toanncthees stablesimulationsforanadditional108-109 years. Fromthese
runswethenlearnhowrobustthepredictedstableregionsare,
a and eccentricities; some are 10 times closer to their primary
: than Mercury, and others orbit with eccentricities larger than and we also determine the mechanisms that lead to stability.
v Specificallyweevaluatethehypothesisthatsomestablesystem
i 0.5. In this paper we attempt to categorize these systems dy-
X requiresecularresonancelocking.
namically, constrainthe errorsof the orbitalparameters,com-
Inmanywaysthispaperperformsthedirectanalysisthatis
r pare our SS to ESPS, explore the long term stability of each
a approximated by MEGNO (Cincotta & Simó 2000, see also
planetary system, and determine the mechanism(s) that main-
Robutel & Laskar 2001, Michtchenko & Ferraz-Mello 2001,
tainstability.
and Goz´dziewski 2002). MEGNO searches parameter space
We examine these systems through numerical simulations.
for chaotic and periodic regions. Our simulations show that
Theintegrationsbeginwithslightlydifferentinitialconditions
to first order MEGNO’s results do uncover unstable regions.
inordertoprobeobservationallyallowedconfigurations. This
In general, however,chaotic systems canbe stable forat least
explorationofparameterspace permitsa quantitativemeasure
109 years, as shown below. Our SS also shows chaos on all
ofthestabilityofeachsystem,and,hence,predictswhichdis-
timescales(foracompletereview,seeLecar,etal.2001),there-
tribution of orbital elementswill most likely result in a stable
foreonlydirectN-bodyintegrationcandeterminestability.
system. Inaddition,acomparisonofstabilitybetweensystems
Thisworkrepresentsthelargestcoherentstudyofplanetary
may revealwhichelementsare mostcriticalto the stability of
system dynamics to date. Our simulations show that the true
planetarysystemsingeneral.
configurations of most planetary systems are constrained by
A quick inspection of the known systems reveals three ob-
justafeworbitalelements(orratiosofelements),andthatsta-
viousmorphologicalclassifications: resonant, interacting,and
ble regions can be identified with integrationson the order of
separated. Resonant systems contain two planets that occupy
106 years. We alsofindthatstability,aswellasconstraintson
orbits very close to a 2:1 mean motion resonance. A 3:1 sys-
stability,arecorrelatedwithmorphology.Resonantsystemsta-
tem, 55Cnc (Marcy et al.2002), has been announced, but its
bility dependsstrongly on the ratio of the periods, interacting
dynamicswillnotbeexaminedhere. Interactingsystemscon-
systemsoneccentricities,andseparatedsystemsarestable.
tain planets which are not in mean motion resonance, but are
This paper is structured as follows. In §2 we describe the
separatedbylessthana10:1ratioinorbitalperiod. Thesesys-
generation of initial conditions, integration technique, and in-
tems are not dynamicallylocked, but the planetsperturb each
troduce the concept of a stability map. In §3 we analyze the
other. TheSSfallsintothiscategory.Thefinalclassificationis
results of the resonant systems HD82943 and GJ876. In §4
1
2 TheStabilityofPlanetarySystems
weexaminetheinteractingsystemsof47UMa,υAnd,andthe thatthisschemerequirestheazimuthalanglesintheplanetary
SS. In §5 the separated systems, specifically HD168443 and systemstobemeasuredfromthelineofsight.
HD83443, are discussed briefly. In §6 we summarize the re-
2.2. Integration
sultsof§3-5,In§7wediscusspossibleformationscenariosand
inconsistenciesbetweenthisworkandotherworkonplanetary With initial conditions determined the systems were then
systems. Finally in §8 we draw generalconclusions,and sug- integrated with the RMVS3 code from the SWIFT suite of
gestdirectionsforfuturework. programs1(Levison & Duncan 1994). This code uses a sym-
plectic integration scheme to minimize long term errors, as
2. NUMERICALMETHODS well as regularization to handle close approaches. The initial
Inthissectionweoutlinethetechniquesusedtoperformand timestep,∆t isapproximately1/30thoftheorbitaltimeofthe
analyze these simulations. This paper follows the example of innermostplanet.Inordertoverifytheaccuracyoftheintegra-
Barnes andQuinn (2001,hereafterPaper I). First we describe tions, the maximum change in energy, ǫ, permitted was 10- 4.
howtheinitialconditionsofeachoftheshorttermsimulations Wedefineǫas
max|E - E |
are determined. In §2.2 we describe our integration method. ǫ= i 0 , (2)
Finally in §2.3 we describe a convenientway to visualize the E0
resultsofthesesimulations,thestabilitymap. where E represents an individual measurement of the energy
i
during the simulation, and E is the energy at time 0. There
0
2.1. InitialConditions are two reasons for using this threshold in ǫ. First, as the in-
tegration scheme is symplectic, no long term secular changes
Inordertoexploreallofparameterspacewemustvary5or-
willoccur,sohighprecisionisnotrequired. Second,thesim-
bital elements, (the period P, the eccentricity e, the longitude
ulationsneededto becompletedin atimelymanner. Ifa sim-
of ascending node Ω, the longitude of periastron, v , and the
ulation did not meet this energy conservationcriterion, it was
inclinationi)andeverymassinthesystem. Theperiodandthe
rerunwiththetimestepreducedbyafactorof10.Theminimum
semi-majoraxisaarerelatedbyKepler’sthirdlaw. Theargu-
timestepweusedwasP /3000. Despitethissmalltimestep,
mentofpericenter,ωisthedifferencev andΩ. inner
afewsimulationsdidnotconserveenergy,andwerediscarded,
Foreachsystemweperform1000simulations,eachwithdif-
exceptthat they were incorporatedinto the errors. Errorsand
ferentinitialconditions. Orbitalelementsthatareeasilydeter-
error bars include information from unconserved simulations.
minedviatheDopplermethod(P,e,v )arevariedaboutaGaus-
Simulations which fail to conserve energy would most likely
sian centered on the nominal value, with a standard deviation
belabeledasunstable,asthefailureofenergyconservationun-
equaltothepublishederror. We donotpermitanyelementto
doubtedlyresults froma close encounterbetweentwo bodies,
bemorethan5σ fromthemean. FortheiandΩelements,flat
which usually results in an ejection. Therefore the estimates
distributions in the ranges 0o <i<5o and 0<Ω<2π were
forstabilityinthesystemspresentedhereshouldbeconsidered
used. NotethatthisdistributionofiandΩpermitsamaximum
upperlimits.
mutualinclinationof 10o. These randomizedorbitalelements
Throughoutthispaperweadoptthenomenclatureofthedis-
arerelativetothefundamentalplane.
coverypapers(planetshavebeenlabeledb,c,d,etc,withorder
TheDopplermethodofdetectiondoesnotproduceanormal
in the alphabetcorrespondingto order of discovery). We will
errordistribution.AsKonacki&Maciejewski(1999)show,the
also introduceanotherschemebasedonmass. Planetswillbe
errorineccentricitycanhavealargetailtowardunity.However
subscripted with a 1, 2, 3, etc, in order of descending mass.
their method, or other statistical methods, such as bootstrap-
Thisnewschemeismoreusefulindiscussingthedynamicsof
ping,requirealltheobservationaldata(includingreflexveloc-
thesystem.
ityerrors)todeterminetheshapeofthiserrorcurve.Asnotall
The short term simulations are integrated until one of the
theobservationsofmultiplanetsystemshavebeenpublishedwe
following criteria is met. 1) The simulation ejects a planet.
areforcedtouseanormaldistributioninorderforcomparison
Ejectionisdefinedasthe osculatingeccentricityofoneplanet
between systems to be meaningful. We therefore encourage
reaching, or exceeding unity. 2) The simulation integrates to
all the observational data to be published, as some of the re-
completionattimeτ. Forthesesimulations,τ isdefinedas
sults presented here (specifically comparingthe percentageof
simulationsthat survived)are less meaningfulbecause we are τ ≡2.8×105P , (3)
1
unable to accurately estimate the error distributions of the or-
bital elements. For completeness we also vary the primaries’ or 280,000times the period of the most massive planet. This
masses,M ,basedonotherobservations(generallydetermined choicecorrespondsto 106 yearsforthe υ Andsystem, aswas
∗
via spectral fitting). However as the stars are all at least 100 simulatedinPaperI.
timesmoremassivethantheirplanetarycompanions,theslight Ifasystemendswithoutejection,thenthestabilityofthesys-
variationsinprimarymassshouldnotaffectstability. temmustbedetermined. Thereareseveralpossibledefinitions
The masses of the planets are determined by the following ofstability. InPaperI,asystemwasstableiftheosculatingec-
relation centricityofeachcompanionremainedbelow1fortheduration
m of the simulation. The most obviousflaw in this definition is
m= obs (1)
|sin(cos- 1(sini)cosΩ)| thataplanetcouldsufferacloseapproachandbethrownoutto
aboundorbitatsomearbitrarilylargedistance. Suchasystem
wheremisthetruemassoftheplanet,andm istheobserved wouldbearnoresemblancetotheobservedsystem,andhence
obs
minimummass. Byvaryingthemassthisway,weaccountfor shouldbelabeledunstable.Hereweadoptamorestringentdef-
allpossibleorientations,andconnecttheinclinationofthesys- inition,namelythatthesemi-majoraxesofallcompanionscan-
teminitsfundamentalplane,toitsinclinationinthesky. Note notchangebymorethanafactorof2. Changesinsemi-major
1SWIFTispubliclyavailableathttp://k2.space.swri.edu/∼hal/swift.html
Barnes&Quinn 3
axisrepresentamajorperturbationtothesystem,thereforethis throughout this paper, instability can arise at any timescale.
secondcutisconservative,andonlyeliminatessystemswhich Therefore,thestabilityzonewillcontinuetoshrinkasthesys-
expelaplanettolargedistanceswithoutfullyejectingit. tem evolves. Second we have chosen a very generous cut in
In addition to these short term simulations, we completed semi-majoraxisspace. Otherstudiespermit∆atobenolarger
simulations to explore longer term stability (∼108 yrs). For than 10% (Chiang, Tabachnik, & Tremaine 2001). Lowering
each system we ran∼10 simulations, chosento covera wide thisvariationwouldundoubtedlyalsoconstrictstabilityzones.
range of stable parameter space. For these runs we started
with the final conditions of stable configurations and contin-
3. RESONANTSYSTEMS
uedthem.Thesesimulationsthereforegiveusahandleonhow
Two systems with orbital periods in 2:1 mean motion res-
the system is likely to evolve on timescales closer to its age
(∼109yrs). Those systems which survived these longer runs onance have been detected: HD829432 and GJ876 (Marcy et
al.2001a). Table1liststheorbitalelementsanderrorsforthe
arethebestcomparisonstothetruesystem. Hencetheyarethe
resonantsystems. For now we do not examinethe 3:1 55Cnc
bestsimulationsfordeterminingthefactorsthatleadtoplane-
system. The currentbest Keplerian fit to the observationsput
tarysystemstability.
HD82943 and GJ876 just beyond perfect resonance. These
Therearetwonotableproblemswiththismethodology.First,
planetsalloccupyhigheccentricityorbitsandhencehavewide
weignoretheeffectsofgeneralrelativity,whichmaybeimpor-
resonancezones. Simulations of these systems show that sta-
tantinsomesystems, specificallyGJ876,andυ And. General
bilityishighlycorrelatedwiththeratiooftheperiods,R,andto
relativitywasincludedinourtreatmentofυAndinPaperI.In
alesserdegreeone . Thesesystemsaretheleaststableasless
PaperI,theinnermostplanethadanegligibleeffectonthesys- 1
than20%ofsimulationssurvivedtoτ.
tem,andwepresumethatgeneralrelativitywillcontinuetobe
unimportantforthesystemsstudiedhere. Second,wetreatall
particlesaspointmasses. Thisisagainespeciallytroublesome
3.1. HD82943
forGJ876andυ And,becauseoftheirproximitytotheir(pre-
sumably)oblateprimaries.Thesphericityofthestarsalsopro- Two planets orbit the 1.05±0.05M⊙ G0 star (Santos et
hibitsanytidalcircularizationofhighlyeccentricplanets(Ra- al.2000) HD82943 at semi-major axis distances of 0.73 and
sio et al.1996). Thismay artificially maintainlarge planetary 1.16AU.Planetbistheinnerandlessmassive,c,theouterand
eccentricities, and increase the frequencyof close encounters. moremassive. Goz´dziewski&Maciejewski(2001)alsoexam-
Howevertheeccentricitiesmustbecomeverylargeforthisef- inedthissystem. Theyfoundthatthesystem ismostlikelyto
fecttobecomeappreciable.Wethereforeassumethatthisphe- bestableinperfectresonance. Thesesimulationsmustrunfor
nomenon will not adversely affect our results. Ignoringthese 340,830 years, τHD82943. The long term simulations were run
twoissuesshouldimpacttheresultsminimallywhilespeeding for100millionyears.
upoursimulationsconsiderably. Forthissystemwefindthat18.8%±4.3%ofthetrialswere
unperturbedtoτ ,and19.3%survivedfor106years,and
HD82943
2.3. StabilityMaps 4.5%failedtoconserveenergy.Fig.1showstheinstabilityrate
to106 years. Mostunstablesimulationsbreakourstabilitycri-
Whenanalyzingthesimulations,itisusefultovisualizethe
terionwithin104years,butotherssurvivedmorethan900,000
results in a stability map. In generala stability map is a 3 di-
yearsbeforebeingejected. Theasymptoticfalloffto106years
mensionalrepresentationofstabilityasafunctionof2parame-
implieswehavefoundmostunstableconfigurations. In91.3%
ters. Inresonantsystems,wefindthatseveralparametersdeter-
of the trials, planet b, the inner and less massive planet, was
mine stability. Themost importantis the ratios ofthe periods
ejected/perturbed.Inordertocheckthesimulations,Fig.2plots
ofthetworesonantplanets,whichwewillcallR:
therateofsurvivalasafunctionofenergyconservation. From
P thisfigureitappearsthatourlimitof10- 4isreasonable,asthere
R≡ outer. (4) appearstobenotrendinstabilityasafunctionofenergycon-
P
inner servation.Thespikeinsurvivabilityat10- 8correspondstoreg-
Incoupledsystems,e ande aretherelevantparameters. The ularorbitsthatwerestableforourinitialchoiceof∆t.Thelack
1 2
advantageto this visualizationis that boundariesbetweensta- ofatrendwithenergyconservation(specificallysurvivalprob-
bilityandinstabilityareeasilyidentified. ability increasing with decreasing ǫ), implies our cutoff value
The disadvantage of this form of visualization is that if the ofǫisstringentenough.
range of parameter space is not uniformly sampled (as it is StabilityinthissystemiscorrelatedwithR,e ,∆Mand∆v ,
c
here), we cannotvisualizeof the errors. It isthereforeimpor- where ∆M is the difference in mean anomaly, and ∆v is the
tanttobearthisdisadvantageinmind. Attheedgesofstability difference in initial longitude of periastron. Slightly beyond
maps, the data are poorly sampled and the information at the perfect resonance is the preferred state for this system. This
edges should largely be ignored. To aid the visualization we system also requiresthe eccentricityof planetc to remainbe-
have smoothed the maps. If a bin contained no data then it low 0.4. These features are shown in Figs.3 and 4. In these
wasgiventheweightedmeanofalladjacentbins,includingdi- greyscaleimages,blackrepresentsunsampledregions,darkest
agonalbins. This methodologycan producesome misleading greymarksregionsinwhichnoconfigurationssurvived,grades
featuresin the stability maps. Mostnotable are tall spikes, or ofstability aredenotedby lightershadesofgrey,andwhite is
deepdepressionsinsparselysampledregions.Wecommenton fullystable.Aswithmoststabilitymapsinthispaper,theouter
thesetypesoferrorswhereappropriate. 2-3gridpointsshouldbeignored. InFig.3,theR- e stability
c
Theprocedureasoutlinedoverestimatesthesizeofstablere- map,thelarge“plateau”atlowe,isthereforepoorlysampled,
gions in two important ways. First, the integration times are asistheislandatR=2.15,e =0.52.Themoststrikingfeatureof
c
generallylessthan0.1%ofthesystems’trueages. Asisshown thisfigureisthatthebestfittothesystem, theasterisk, places
2http://obswww.unige.ch/∼udry/planet/hd82943syst.html
4 TheStabilityofPlanetarySystems
Table1- Initial ConditionsforResonantSystems
System Planet Mass(MJ) Period(d) E

entri
ity Long. of Peri.(Deg) Time of Peri. (JD)
GJ876
0.56 30:12(cid:6)0:02 0:27(cid:6)0:04 330:0(cid:6)12:0 2450031:4(cid:6)1:2
b 1.89 61:02(cid:6)0:03 0:10(cid:6)0:02 333:0(cid:6)12:0 2450045:2(cid:6)1:9
HD82943 b 0.88 221:6(cid:6)2:7 0:54(cid:6)0:05 138(cid:6)13 2451630:9(cid:6)5:9

1.63 444:6(cid:6)8:8 0:41(cid:6)0:08 96(cid:6)7 2451620:3(cid:6)12:0
1
Barnes&Quinn 5
itadjacenttostability. IfRischangedbylessthat5%thesys-
temhasnochanceofsurvivingeven1millionyears. Thismap
showsthatthecurrentfittothesystemisnotcorrect. However
the elements do not need to change by much (specifically e
c
needs to be slightly lower) for the system to have a chance at
stability.
0.1
0.05
0
0 2 4 6
FIG. 3.—TheR- ec stabilitymapforHD82943. Theasteriskrep-
resents the best fit to the system as of 31 July 2002. The data are
FIG. 1.—Thedistributionofinstabilitytimesforunstableconfigu- mostaccurateclosertotheasterisk.Thesystemshowsaclearbound-
rationsofHD82943. Mostunstablesystemssurviveforjust102-104 ary in phase space between unstable (dark grey) regions and stable
yearsbeforeperturbationschangeasemi-majoraxisbymorethana (white)regions. Blackrepresentsunsampleddata. Thestableregion
factorof2. atR=2.15,e =0.52isabininwhich1outof1trialssurvived.
c
The system also shows dependence on mean anomaly and
longitude of periastron. Because of the symmetry of ellipses
wewillintroduceanewvariable,Λ,definedas:
|v -v |, Λ<π
1 Λ≡(cid:26) 3610- |v21-v 2|, Λ>π (5)
0.8 wherethesubscriptsmerelyrepresent2differentplanets,band
cforthissystem. Theorderisunimportant,astheweareonly
concernedwith the magnitude of this angle. In Fig.4, the Λ-
0.6 ∆Mstabilitymapispresented. Theasteriskmarksthebestfit
tothesystem. Stabilityseemstofollowalinerepresentedby
0.4
4 4 π
∆M= Λ+120= (Λ+ ). (6)
3 3 2
0.2
Thisrelationispurelyempirical. Asisshowninthefollowing
0 sections, this interdependencyis unusual for extra-solar plan-
etary systems. Similar plots of ∆M or Λ versus R, show the
-12 -10 -8 -6 -4
sameRdependenceasinFig.3. ThereforeRisclearlythemost
importantparameterinthissystem,buttheseother3alsoplay
an important role in the system. As more observationsof the
FIG. 2.— Survival rate as a function of energy conservation for systemaremade,HD82943shouldfallintotheregiondefined
HD82943. Thelackofatrendimpliestheresultsforthesystemare
by1.95≤R≤2.1,e <0.4,andEq.6.
accurate. c
6 TheStabilityofPlanetarySystems
becomeslockedatΛ≈100o. Althoughnotshown,afterthese
initialperturbationsthesystemsettlesdowntoaconfiguration
inwhicha ≈0.344AU,e ≈0.79anda ≈2400AU,e ≈0.99.
b b c c
1.2 0.8
0.6
1
0.4
0.8
0.2
0 0
Time (Yrs) Time (Yrs)
8 0.08
Outer
6 Inner 0.06
4 0.04
0.02
2
0
0 0 20 40 60
Time (Yrs)
FIG. 5.—OrbitalevolutionofHD82943-348, astable,regularcon-
figuration. Thedata herearesmoothed over a20,000 year interval.
FIG. 4.—TheΛ-∆MstabilitymapforHD82943. Theasteriskrep- TopLeft: Theevolutioninsemi-majoraxis. Theplanetsshowslight
resentsthebestfittothesystemasof31July2002.Thedataaremost variations due to resonant interactions. Top Right: The eccentric-
accurateclosertotheasterisk.Stabilityappearstofollowthelinerep- ities oscillate with a 700 year period, with 0.62 ≤e ≤ 0.85 and
b
resentedbyEq6.Notehoweverthatthesystemisalsoconstrainedto 0.05≤e ≤0.45. This short timescale is not visible in this plot.
c
Λ≤75oand∆M≥30o.TheislandatΛ=80o,∆M=20oisabinin Bottom Left: The inclinations experience oscillations on 2100 year
which1outof1trialssurvived. timescale, again too short to be visible in this plot. The ranges of
inclination are 1o ≤i ≤7o and 0≤i ≤3o. Bottom Right: From
HD82943showsawiderangeofdynamics. Someexamples b c
this distribution function we see the lines of node librate harmoni-
oftheseareshowninFigs.5-8. Theinitialconditionsofthese
cally withan amplitude of 60o. The dashed vertical linerepresents
systems are presented in Table 2. Fig.5 shows an example of Λ ,theinitialvalueofΛ.
theevolutionofaregularsystem,HD82943-348,whichshows 0
noevidenceofchaos. NotethatinsteadofΛ(t),wepresentthe
distribution function of Λ, P(Λ), the probability of Λ, vs. Λ.
This representation of Λ shows that the motion is like that of
aharmonicoscillator;thelongitudesofperiastronarelibrating
withanamplitudeof60o.
Fig.6 (HD82943-382) is a stable case which is clearly
chaotic.Althoughtheeccentricitiesremainclosetotheirinitial 1.2
0.6
values,theinclinationsjumptolargevaluesquickly. Notethat
Λ never exceeds 750, but its motion is slightly nonharmonic,
1
anotherindicationofchaos. 0.4
InFig.7weplottheevolutionofasystemwhichejectsplanet
b. Thissystemexperiencescloseapproachesveryquicklyasis 0.8 0.2
evidencedinthetopleftpanelofFig.7. Theeccentricitiesun-
dergo dramatic fluctuations, and the inclinations nearly reach 0 0
40o. Althoughpoorlysampled,Λsuggestscirculation. Time (Yrs) Time (Yrs)
In Fig.8, the orbital evolution of simulation HD82943-216 30
is shown. This system perturbedplanetc after 280,000years, Outer
Inner
and despite the high eccentricitiesthe system obtained(>0.75 0.02
20
for bothplanets), remainedboundfor 106 years. The inclina-
tions also show large growth. Although initially Λ=41o, the
0.01
system becomes locked at just over 100o. This is misleading 10
as this is the Λ distribution for the full 106 years. It is actu-
allythesuperpositionof2modes.Initially,untilperturbationat 0
280,000years,thesystemcirculateswithaslightpreferenceto- 0 0 20 40 60 80
Time (Yrs)
wardanti-alignment.Howeveraftertheperturbationthesystem
Barnes&Quinn 7
Table 2-Sele
ted Simulationsof HD82943
o a
Trial eb;0 e
;0 R0 (cid:3)0 ( ) (cid:15) Comments
(cid:0)7
216 0.563 0.345 2.0108 42.2 2:6(cid:2)10 C,P(
,299.9)
(cid:0)8
348 0.617 0.419 1.9656 13.6 1:5(cid:2)10 R
(cid:0)9
382 0.516 0.322 2.0303 32.4 9:2(cid:2)10 C
(cid:0)10
698 0.481 0.440 2.0948 68.4 1:3(cid:2)10 C,P(b,7.5),E(b,10.5)
a 3
R=Regular,C=Chaoti
,P=Perturbed(Planet,Time(10 yr)),E=Eje
ted(Planet,
3
Time(10 yr))
1
8 TheStabilityofPlanetarySystems
FIG. 6.—OrbitalevolutionofHD82943-382, astable,chaoticcon-
12
figurationofHD82943. Thedataareaveragedona10,000year in- p 1 p
terval. Top Left: The evolution in Semi-Major Axis. The planets 10
clearlyneverexperienceacloseencounter. TopRight:Theeccentric- 0.8
8
itiesexperiencelowamplitude(≈13%)chaoticoscillations. Bottom
0.6
Left: The inclinations initially jump to large values and experience 6
largeamplitude(≈70%)fluctuations. BottomRight:TheΛdistribu- 0.4
4
tionfunctionsuggeststhatΛisgenerallylibrating,butthatthereare
chaoticfluctuationssuperimposedonthismotion.Thedashedvertical 2 0.2
lineisthevalueofΛ .
0 0 0
Time (Yrs) Time (Yrs)
50
p
Outer
40 Inner
0.1
30
20 0.05
10
0
0 0 50 100 150
25 p 1 p Time (Yrs)
Outer
20 Inner 0.8 FIG. 8.— Orbital evolution of HD82943-216, the perturbation of
HD82943c. The parameters are averaged on 10,000 year intervals.
15
0.6 TopLeft: Theplanets’semi-majoraxesarestableandshownosigns
10 ofcloseencountersuntil260,000years.Atthispointplanetcactually
0.4
crossesb’sorbit. Thisinitialencounter leadstomoreencountersas
5 0.2 acreaches3AUby280,000years,trippingthecriterionforinstability.
TopRight:Theeccentricitiesexperiencesecularchangeuntil210,000
0
0 2000400060008000 0 2000400060008000 years. The system then moves into a lower eccentricity state. The
Time (Yrs) Time (Yrs) eccentricitythengrowstolargevaluesandremainattheirfinalval-
40 p 0.04 uesforanother700,000years. BottomLeft: Aswitheccentricitythe
inclinationsshowslowsecularchangeuntil210,000years.Theincli-
30 0.03 nationsthenleapupto30ointhecaseofplanetb.BottomRight:The
Λdistributionfunctionisthesumof2motions: thepre-perturbation
20 0.02 motion is circulation, the post-perturbation motion is fixed close to
110o.ThedashedlinerepresentΛ .
0
0.01
10
0
1
0 2000400060008000 0 50 100 150
Inner
Time (Yrs)
0.8 Outer
FIG. 7.—Orbital evolutionof HD82943-698, theperturbation and
ejectionofHD82943b.Theorbitalparametersaresampledevery100 0.6
years.TopLeft:Theplanetsexperiencecloseencountersimmediately
aschangesinsemi-majoraxisarevisiblewithin500years.TopRight: 0.4
Theeccentricitiesexperiencelargeamplitudefluctuations,culminat-
0.2
ing inthe ejection of planet b just after 10,000 years. Bottom Left:
The inclinations vary slightly until a close encounter at 1100 years
0
whichsendsbothinclinationsupsubstantially. Theplanetsremainin
thisstateuntiljustpriortoejection,whentheyreturntocoplanarity.
0.8
Bottom Right: The Λ distribution function suggests that Λ remains
below100o,butmayoccasionallycirculate.
0.6
0.4
0.2
We ran 10 simulationsfor108 years. The initial conditions 0
0 20 40 60 80 0 20 40 60 80 100
and results of these simulationsare presentedin Table 3. The
eccentricityevolutionof4simulationsisshowninFig.9.These
simulations, whichranfor100millionyears,showsomecon- FIG. 9.—Long term simulations of the HD82943 system. These
dataareaveragesover50,000yearintervals. TopLeft: Evolutionof
figurations are regular (top left), some are chaotic (top right,
HD82943-035.Anexampleinwhichthesystemisregular.TopRight:
bottomleft),andthatinstabilitycanariseatanytimescale(bot-
EvolutionofHD82943-021,anexampleofchaoticevolution.Bottom
tomright).Theselongtermsimulationsshowthatregionsexist
Right: EvolutionofHD82943-000, another exampleof chaoticmo-
inphasespaceinwhichthissystemcansurviveforbillionsof
tion. Bottom Left: Evolution of HD82943-032. A chaotic system
years. whichejectsplanetbafter2.4millionyears(7τ ).
HD82943
Barnes&Quinn 9
8
Table3- Resultsof LongTerm (10 yr) SimulationsofHD82943
o a
Trial eb;0 e
;0 R0 (cid:3)0 ( ) (cid:15) Result
(cid:0)8
000 0.561 0.395 2.0269 26.5 5:2(cid:2)10 C
(cid:0)8
021 0.571 0.333 2.0615 39.0 1:0(cid:2)10 C
032 0.414 0.556 2.0556 33.3 4.95 (cid:15)(2.4),E(b,2.4)
(cid:0)9
035 0.511 0.236 2.0322 6.7 6:8(cid:2)10 R
(cid:0)9
036 0.534 0.197 2.0080 29.5 6:1(cid:2)10 R
040 0.600 0.393 2.0766 52.8 0.22 (cid:15)(1.005),E(
,1.14)
(cid:0)8
043 0.542 0.386 2.0223 24.8 2:4(cid:2)10 C
(cid:0)8
057 0.549 0.348 2.0382 27.0 1:1(cid:2)10 C
(cid:0)8
073 0.468 0.265 2.0626 56.8 6:3(cid:2)10 R
a 6
R=Regular,C=Chaoti
,(cid:15)=Energy
onservationfailed(Time(10 yr)),E=Eje
tion(Planet,
6
Time (10 yr))
1
10 TheStabilityofPlanetarySystems
Severalpapershavesuggestedthatsecularresonancelocking ofparameterspace isstable. On 106 yeartimescales, 2.4%of
maintainsstabilityinESPSwithlargeeccentricities(Rivera& configurationssurvived,butonly1.7%werestill unperturbed,
Lissauer2000,Rivera&Lissauer2001a,Chiang,Tabachnik,& and 5% failed to conserve energy. In unstable cases planet c,
Tremaine2001). Specifically theysuggestthatthe orientation the inner and less massive, was ejected/perturbednearly 96%
ofthe planets’ellipsesshouldbealigned(Λ≈0). Byexamin- ofthetime. Fig.11showstheinstabilityrate. GJ876showsthe
ingΛinstable,regularlongtermsimulationswecandetermine sametrendasHD82943;mostunstableconfigurationsbreakup
if the longitudesof periastronremain locked. The probability in just a few hundreddynamicaltimes. Again the asymptotic
distribution of Λ for these same four long term simulations is nature of this plot implies most unstable configurations have
presentedin Fig.10. Λ is sampledonceevery100years. The beenidentified. Thissystemshowsthesamesortoftrendinǫ
topleftplotofFig.10issimilartothatofanharmonicoscilla- asHD82943,whichprovesthatourresultsourcorrect.
tor;thisconfigurationislibratingaboutΛ=0withanamplitude
of40o. Theotherplotsaresystemswhicharenotlibrating,but
insteadshowmorerandommotion.Fromthisfigureweseethat
regularmotioniscorrelatedwithlibrationaboutalignment,but
chaoticandunstablegenerallyshowschaoticΛbehavior.
0.4
0.3
0.2
0.1
0
0 2 4 6
FIG.10.—ThedistributionfunctionofΛ(sampledevery100years)
forfourstablecasesofHD82943.Thesefoursystemsarethesameas
inFig.9. ThetopleftplotsasystemlibratingaboutΛ=0. Theother FIG. 11.— The ejection rate for GJ876. In this system unstable
plotsshowthatchaoticmotionisusuallyassociatedwithacirculating configurationsareusuallyejectedwithin100years. Therateasymp-
Λ. toticallyfallstozeroby105years.
3.2. GJ876
Two planets orbit the 0.32±0.05M⊙ (Marcy et al.2001a)
M4starGJ876,alsoknownasGliese876. Thissystemisvery
similar to HD82943, the major differencebeing that the plan-
etslieclosertotheirprimary.Thesemi-majoraxesofthesetwo
planetsare0.13and0.21AU.Notethatinthissystemplanetcis UnlikeHD82943,therearenoobviouszonesofstability.We
theinnerandlessmassiveplanet.Asubstantialamountofwork only see isolated islands in the R- e stability map presented
b
hasalreadybeendoneonthissystem.Notably,inthediscovery inFig.12. Wechoosetheseparametersforourstabilitymapas
paper(Marcyetal.2001)thatstableconfigurationsexistinthe theywerethestrongestindicatorsofstabilityinHD82943,and
system. Rivera&Lissauer(2001b)showthatKeplerianfitting because of the suggestion that the system actually lie in per-
for this system is not precise enough to accurately determine fectresonance. ClosetoR=2.00wesampledtwosimulations
the orbits. They suggest, through N-body fitting, that GJ876 near R=2.02 and e =0.7. Both these were stable. However
b
must actually lie in perfect resonance, and that the orbital el- withsuchpoorstatistics,andatsuchalarge(relative)distance
ements providedin the discoverypaper (whichare used here) fromperfectresonance,we cannotcommentonthelikelihood
suffer from a systematic error. Because P is so short for this thatthesystemwouldbemorestableinperfectresonance. We
1
system(60days)theevolutionoftheorbitalelementshasbeen can, however, point out that there are isolated configurations
observed, andcorroboratedthis theory. Thissection therefore which may hold stable orbits, and prolonged observations of
servesasachecktothishypothesis. thissystem shoulddemonstratewhetherit isindeedin perfect
We ranGJ876for106 years, butτ correspondstoonly resonance. Howeverthislackofalargestableregionstrength-
GJ876
47,000 years. On 47,000 year timescales, only 5.6%±2.8% ensthehypothesisthatthissystemliesinperfectresonance.