Chaotic dynamics of the planet in HD 196885 AB

5.8 MB·

by  Suman Satyal

Mon.Not.R.Astron.Soc.000,1–??(2002) Printed19June2014 (MNLaTEXstylefilev2.2)
Chaotic dynamics of the planet in HD 196885 AB
S. Satyal,1(cid:63) T.C. Hinse2,3 B. Quarles,4 J.P. Noyola1
1The Department of Physics, University of Texas at Arlington, Arlington, Texas, 76019, USA.
2Korea Astronomy and Space Science Institute, 304-358 Daejeon, Republic of Korea.
3Armagh Observatory, College Hill, BT61 9DG, Armagh, UK.
4
4Space Science and Astrobiology Division 245-3, NASA Ames Research Center, Moffett Field, CA 94035, USA.
1
0
2
Released2014January6
n
u
J
ABSTRACT
8
Depending on the planetary orbit around the host star(s), a planet could orbit ei-
1
ther one or both stars in a binary system as S-type or P-type, respectively. We have

analysed the dynamics of the S-type planetary system in HD 196885 AB with an
]
P emphasis on a planet with a higher orbital inclination relative to the binary plane.
E The mean exponential growth factor of nearby orbits (MEGNO) maps are used as
. an indicator to determine regions of periodicity and chaos for the various choices of
h
the planet’s semimajor axis, eccentricity and inclination with respect to the previ-
p
ouslydeterminedobservationaluncertainties.Wehavequantitativelymappedoutthe
-
o chaotic and quasi-periodic regions of the system’s phase space which indicate a likely
r regime of the planet’s inclination. In addition, we inspect the resonant angle to deter-
t
s minewhetheralternationbetweenlibrationandcirculationoccursasaconsequenceof
a Kozai oscillations, a probable mechanism that can drive the planetary orbit to a very
[
large inclination. Also, we demonstrate the possible higher mass limit of the planet

2 and improve upon the current dynamical model based on our analysis.
v
Key words: methods: numerical - celestial mechanics - binaries: general - planetary
8
6 systems - stars: individual: HD 196885
2
1
.
1
0 1 INTRODUCTION alistofcandidatehostsofsuchmoonsaroundtransitingex-
4 oplanets. Noyola et al. (2013) recently calculated the mini-
1 The study of exoplanets has made a significant leap since
mumrequiredfluxfromtheexoplanetaryradioemissionsin
: thefirstplanetaroundasolartypestar,51Pegasib(Mayor
v ordertodetecttheexomoonsusingseveralgroundbasedra-
& Queloz 1995), was discovered almost two decades ago.
i diotelescopes.Althoughanexomoonhasyettobedetected,
X The number of confirmed exoplanets has already surpassed
theorists have already determined relations to constrain an
most peoples’ expectations with a count of 1792 as of June
r exomoon’s habitability by considering energy flux (i.e., ra-
a 2, 20141. New exoplanets are being added almost everyday
diative and tidal) and orbital stability (Heller 2012; Cuntz
into the database. The ultimate conquest of the contem-
et al. 2013).
porary exoplanet research is to find a planet within its re-
spective habitable zone with a stable planetary orbit and a The focus of this paper, however, is the study of the
terrestrial mean density. Various ground based surveys and dynamicsofanexoplanetwhoseorbitaleccentricityandin-
spacebasedtelescopes,suchastheKepler2 missionandfu- clinationrelativetothebinaryplanearesignificantlyhigher
ture missions such as the transiting exoplanet survey satel- thanthemedianofthosediscovered.Theorbitalstabilityof
lite(TESS)(Rickeretal.2010)andtheJamesWebb3 space suchasystemdependsonvariousorbitalparametersandthe
telescope are set for the detection and characterization of long-term stability is vital for life to develop under current
exoplanetary systems which will deliver a promising future theories. Thus, understanding the complete dynamics of a
in the search of terrestrial planets which could support life. systembearsvitalrequirement.Variousnumericaltoolshave
Meanwhile,thesearchforexomoonsisalsoanongoingpro- been developed and used in the past to address the orbital
cessandrecentworkbyKippingetal.(2013)hasproduced configuration leading to the stability, instability or chaos
of planets in a binary system. Quarles et al. (2011) used
the maximum Lyapunov exponent (MLE), originally devel-
oped by Lyapunov (1907), to determine the orbital stabil-
(cid:63) e-mail:[email protected]
1 www.exoplanet.eu ityorinstabilityforthecircularrestrictedthree-bodyprob-
2 http://kepler.nasa.gov/ lem (CRTBP) case. The mean exponential growth factor
3 http://www.jwst.nasa.gov/ of nearby orbits (MEGNO) (Cincotta & Simo´ 1999) maps
(cid:13)c 2002RAS
2 Satyal et al.
havebeenusedtostudythedynamicalstabilityofirregular
HD 196885 B Ab
satellites(Hinseetal.2010)andextrasolarplanetdynamics
(Go´zdziewski&Maciejewski2001;Go´zdziewskietal.2001). Mass(m) 0.45M(cid:12) 2.98MJ
Recently, Satyal et al. (2013) used MLE, MEGNO and the SemimajorAxis(a) 21±0.86AU 2.6±0.1AU
Hillstabilitytimeseriesmethods(Szenkovits&Mako´2008) Eccentricity(e) 0.42±0.03 0.48±0.02
Inclination(i) 0◦ [0◦ -90◦]
tostudytheorbitalperturbationofaplanetduetothestel-
lar companion in the binary systems of γ Cephei and HD ArgumentofPeriapsis(ωp) 241.9◦ ±3.1 93.2◦ ±3.0
Ascendingnode(Ω) 79.8◦ ±0.1 0.0◦
196885.
MeanAnomaly(M) 121◦ ±45 349.1◦ ±1.8
Planets that are formed from a planetary disk are ex-
pected to have near circular orbits due to tidal circulariza-
Table1.Best-fitorbitalparametersoftheHD196885systemas
tion,yetroughly19%ofknownexoplanetswithwell-known
orbital parameters have an eccentricity greater than 0.44. o=bt1a.3in3eMd(cid:12)fr.omTh(eCphlaaunveitn’seitncalli.na2t0i1o1n).(iMpla)sissomfeparsiumreadryresltaatri,vemtAo
Thisestimate,however,couldbebiasedandmaynotinclude thebinaryorbit,thenfori =0◦,theplanetaryorbitiscoplanar
pl
the true population. The actual intrinsic phenomenon on withthebinaryorbit.
howsuchhighlyeccentricorbitsareformedwhilemaintain-
ing orbital stability is not quite fully understood, however
some theories do exist that postulate such origins. Kozai totalsimulationtime(lifetime map),weconsideritastable
(1962) and Lidov (1962) proposed a mechanism, now re- system.Iftheplanetgetsejectedfromthesystemorcollides
ferred to as Lidov-Kozai Mechanism, to explain the varia- withthecentralbodywithinthesimulationtime,itwouldbe
tions of a test particle’s eccentricity and inclination. This anunstablesystem.Forsomeinitialconditionstheplanetary
mechanism was first applied to the exoplanets by Holman orbitdisplaysachaoticmotion,whichweshallconfirmfrom
et al. (1997) while studying the chaotic variations in the the evaluation of the time evolution of the resonant angle.
eccentricity of the planet in 16 Cygni AB. Furthermore, its Eventhoughthemotionischaoticitdoesnotmeanthatits
applicationtotheplanetsincludestheworkdonebyInnanen unstableandthesystemcanmaintain(insomecasesithas
et al. (1997), Wu & Murray (2003), Fabrycky & Tremaine maintained) a stable configuration for millions of years.
(2007),Wuetal.(2007),Veras&Ford(2010),Correiaetal. Thispaperisoutlinedasfollows.InSection2wepresent
(2011), Lithwick & Naoz (2011), Naoz et al. (2011), Katz thebasictheoryoftheLidov-Kozai mechanism,ournumer-
et al. (2011), Naoz et al. (2012), Naoz et al. (2013) and Li ical approach, a discussion of MEGNO as a chaos indicator
et al. (2013). Some other mechanisms have also been pro- andtheinitialsetupofthesystem.InSection3wepresent
posed,suchasplanet-planetscatteringtoexplainsuchvaria- ourresultsfromarepresentativesampleofsinglyandmulti
tions(seeforexample,Rasio&Ford(1996);Chatterjeeetal. integratedorbitsandMEGNOmapsfollowedbydiscussion.
(2008); Nagasawa et al. (2008)). In this work, we make use WeconcludeinSection4withabriefoverviewofourresults.
of the Lidov-Kozai mechanism to explain such highly ec-
centric but stable orbits of the S-type planet, HD 196885
Ab. The planet is part of a binary star system (HD 196885
2 THEORY
AB) whose eccentricity is observationally determined to be
0.48andtheplanet’sinclinationisunconstrainedwhichim- 2.1 The Lidov-Kozai Mechanism
pliesavalueanywherebetweenzeroandninetydegrees(for
prograde motion) with the binary plane. Kozai (1962) developed an analytical theory to explain the
The dynamics of the planet in HD 196885 have been secularperturbationsinducedbyJupiterontheasteroidsin
studied in a previous work by Satyal et al. (2013). In that the Solar System. Similar theory was developed by Lidov
study,investigationswereperformedinarestrictedmanner (1962) to study the evolution of orbits of artificial satellites
where the planet’s inclination, (i ) (cid:54) 25◦ with the binary of planets that are directly influenced by the gravitational
pl
plane, and the dynamics of the system was not fully ex- perturbationsoftheSun.Forthis,Kozaiconsideredtheper-
plored. For this paper, we have considered a full range of turber’s (Jupiter) orbit to be circular. Thus the asteroid’s
prograde orbits in i value, from 0◦ to 90◦, and have made vertical angular momentum and the secular energy is con-
pl
useofthechaosindicator,MEGNO,toproducemapswhich servedandthesystemisintegrable.Asaconsequenceofthe
(cid:113)
demonstrateregionsofperiodicityandchaosforavarietyof conservationofthequantity 1−e2 cosi thetimeevolu-
pl pl
initialconditions.Forsimilarinitialconditionsadynamical tion of eccentricity and inclination of the planetary orbit is
lifetime mapisproducedbyusingtheinformationaboutthe anti-correlated:whentheeccentricityissmall,thentheincli-
planet’s ejection time from the system or the collision time nationishighandviceversa.Forrelativeorbitalinclinations
with the stellar host. Also, the dynamics of the system is less than 39.2 degrees the argument of pericentre circulates
studiedintermsoftheplanet’smaximumeccentricity.Then, between 0◦ and 360◦ and the planet’s orbit is precessing.
theresonantangleisanalyzedforevidenceofameanmotion This property remains true for other values of initial ec-
resonanceatthebest-fitlocationoftheplanetandforpossi- centricity.However,forrelativeinclinationslargerthan39.2
blealternationbetweenlibrationandcirculationarisingdue degrees,theplanetbecomesaKozailibrator(Kozairegime)
to chaos induced by Kozai oscillations. withit’sargumentofpericentrelockedandexhibitingeither
If a system demonstrates quasi-periodic orbits (as in librations (oscillations) around 90 or 270 degrees. These li-
the case of MEGNO maps) and/or the planet survives the brationcentersareknowntobestablefixedpoints.AKozai
librator with an initially circular orbit will undergo a large
variationineccentricityandinclinationwithinKozaicycles.
4 www.exoplanet.eu Forrelativeinclinationsclosetooraroundthecriticalvalue
(cid:13)c 2002RAS,MNRAS000,1–??
Monthly Notices: LaTEX2ε guide for authors 3
of39.2degrees,theargumentofpericentreexhibitsintermit- 2.2 Numerical Approach and methods
tent behaviour displaying alternations between circulations
Ourapproachhasconsideredthemotionofaplanetofmass,
and librations. In general, the coupling between eccentric-
m aroundastarofmass,m .Whencalculatingtheinitial
ityandinclinationprovidesaneffectiveremovalmechanism. pl A
conditions in position and velocity, we have used the best-
For a large initial relative inclination the amplitude of the
fit orbital elements, semi major axis (a), eccentricity (e),
planet’s eccentricity variation increases. Following Innanen
inclination (i), argument of periapsis (ω), ascending node
etal.(1997)themaximumextentineccentricityforaKozai
(Ω) and mean anomaly (M), that are obtained from the
cycle can be expressed as,
radialvelocitymeasurements(Chauvinetal.2011).Forany
unconstrained parameters, the values are taken in a range
(cid:114)
5 orconsideredzero.ThelistoforbitalparametersoftheHD
(e ) = 1− cos2(i ) , (1)
max pl 3 o pl 196885systemaregiveninTable1.Thesubscriptsbin and
pl are used to denote the secondary binary component and
where (i ) is the initial relative inclination of the
o pl theplanet,respectively.Themeanlongitude(λ)thatisused
planet.Eventually,foralargeenougheccentricity,theplanet
to calculate the resonant angle (Φ) is calculated from the
iseitherejectedfromthesystemorcollideswiththeprimary
longitudeoftheperiapsis((cid:36)=Ω+ω)andtheorbit’smean
component rendering its orbit to be unstable.
anomaly (M), which is given as λ = (cid:36) + M.
However, the vertical angular momentum and the sec-
Using the orbital integration package MERCURY
ular energy are conserved only when the pertuber’s orbit is
(Chambers&Migliorini1997;Chambers1999),thebuilt-in
circular.Whentheperturberiseccentricoriftheplanethas
Radaualgorithmwasusedtointegratetheorbitsofthesys-
non-negligiblemass,theverticalangularmomentumcompo-
tem in astrocentric coordinates when investigating the evo-
nent of the planet and the perturber is not conserved and
lutionoforbitalelementsforasingleinitialconditionandto
theplanetshowsqualitativelydifferentbehaviour.Forsuch
producealifetime mapformultipleinitialconditions.MER-
cases, Lithwick & Naoz (2011) and Naoz et al. (2013) have
CURY was effective in monitoring the ejection/collision of
furtherdevelopedtheLidov-Kozaimechanismdiscussedear-
aplanetduetocloseencounterwiththesecondarystarand
lier and formulated for the eccentric perturber, which they
provided robust results for the purpose of determining the
havecalledittheeccentricKozaimechanism (EKM).When existence of orbital resonance. A time step of (cid:15) = 10−3
the perturber has an eccentric orbit the quadrupole-order
year/stepwasconsideredtohavehighprecisionbecausethe
approximationisinsufficientandrequirestheoctupole-order
errorinchangeintotalenergyandtotalangularmomentum
term as well (Ford et al. 2000). Following Lithwick & Naoz
wascalculatedateachtimestepwhichfellwithintherange
(2011) and Naoz et al. (2013), the energy function and the of10−16 to10−13 inbothcasesduringthetotalintegration
constant of motion, F, in quadrupole and octupole term is,
period of 50 Myr. The data sampling (DSP) was done per
F ≡ F + (cid:15)F and the constant is given by
quad oct yearforshorterintegrationperiods(upto100kyr)andper
thousand years for longer integration periods. The lifetime
(cid:15)≡ apl eper , (2) map and the maximum eccentricity map are generated for
aper 1−e2per 12,000initialconditionsinapl andipl,andsimulatedfor50
Myr.
where a and a are the semimajor axis of a planet
pl per WeusedtheMECHANIC5software(S(cid:32)loninaetal.2012,
and pertuber respectively, and e is the pertuber’s eccen-
per
2014) optimised to N-body code to calculate the orbits of
tricity.InabinarystarsystemwithanS-typeplanetaround
the given masses and the MEGNO maps on a multi-CPU
the primary star, as in our case, the secondary star would
computing environment. The MEGNO criterion is used to
beaperturber.Whene =0,(cid:15)=0,thusFreducestothe
per
differentiate between the periodic and chaotic dynamics of
case with quadrupole term only. Further mathematical set
a system. Originally developed by Cincotta & Simo´ (1999,
up of equations of motion concerning the quadrupole and
2000), MEGNO has been widely used in the study of dy-
octupole order terms can be found in Naoz et al. (2013).
namical astronomy (Go´zdziewski et al. 2001; Go´zdziewski
During their study of the EKM, some of the remarkable
&Maciejewski2001;Hinseetal.2008,2010,andreferences
results were obtained including the flipping of the plane-
therein). For reason of completeness we outline some de-
tary orbit from prograde to retrograde and its eccentricity
tailsoftheMEGNOformalism.FollowingCincotta&Simo´
reaching to the extreme values (∼ 1). Two main cases have
(2000) MEGNO (Y) is defined as
been discussed regarding the initial conditions in e and
pl
i : high eccentricity, low inclination (HeLi) and low eccen-
pl
tricity,highinclination(LeHi).Primarily,onlytheoctupole 2(cid:90) t (cid:107)δ˙(s)(cid:107)
Y(t)= s ds, (4)
term is in play during HeLi flip (Li et al. 2013) and both t (cid:107)δ(s)(cid:107)
t0
terms (quadrupole and octupole) are in play during LeHi
along with its time-averaged mean value
flip (Naoz et al. 2011, 2013). From Li et al. (2013), the flip
criterion is given as, 1(cid:90) t
(cid:104)Y(cid:105)(t)= Y(s) ds. (5)
t
t0
(cid:15)> 8 1−e2per , (3) The quantity δ is the variational state vector and t is time.
57−eper(4+3e2per)cos(ωper+Ωper) The absolute norm of δ measures the distance between two
nearby points in phase space using a euclidian metric. The
where (cid:15) is the constant parameter (Eq. 2), e , ω
per per
and Ω are the perturber’s initial eccentricity, argument
per
of periapsis and ascending node, respectively. 5 https://github.com/mslonina/Mechanic
(cid:13)c 2002RAS,MNRAS000,1–??
4 Satyal et al.
(a)Initiali =0◦ (b)Initiali =20◦ (c)Initiali =35◦
pl pl pl
(d)Initiali =39.2◦ (e)Initiali =39.7◦ (f)Initiali =45◦,60◦,70◦,80◦
pl pl pl
Figure 1.Evolutionofplanet’seccentricityandargumentofperiapsisforvariouschoicesofplanet’sorbitalinclinationwiththebinary
plane, simulated for 1 ×105 years. The Savitzky-Golay smoothing function with third degree polynomial is used to smooth the high
frequency data as shown in (a),(b) and (c). For the first three choices of i values, the ω is circulating within 0◦ to 360◦. In contrast
pl
to the coplanar case, the amplitude of circulation is varying and is maximum at 90◦ and 270◦ (b and c). When i = 39.2◦ (d), the ω
pl
shows circulation in the region of phase space with eccentricity amplitude maximum at 90◦ and 270◦, and for even higher i (39.7◦)
pl
thephasespacedividesintotwodistinctregionswhereω showscirculationandlibrationaround90◦ and270◦ (e).Thisis,however,an
unstable region where the amplitude of eccentricity oscillation reaches ∼1. When i is set to 45◦, 60◦, 70◦ and 80◦ (f), the ω librates
pl
around90◦ andthelibrationamplitudeineccentricityisfoundtoincreasewithincreasingi untiltheplanetcollideswiththecentral
pl
bodyfori greaterthan80◦.
pl
variationalvectorisfoundfromsolvingthevariationalequa- MEGNOanditsrelationshipwithLyapunovexponents,see
tions of motion in parallel to the newtonian equations of (Hinse et al. 2010, and references therein).
motion (Mikkola & Innanen 1999; Go´zdziewski et al. 2001; In our study, each map is generated with 1.5×105 ini-
Hinse et al. 2010). Equations 4 and 5 can be rewritten tial conditions, for a resolution of (300×500) in various
(Go´zdziewskietal.2001)astwofirst-orderdifferentialequa- orbital elements and simulated up to 200 Kyr. The pur-
tions ple/blue color in the map denotes the quasi periodic region
andtheyellowistheregionofchaos.TheMEGNOquantity
dx δ˙ ·δ dw x (cid:104)Y(cid:105) is colour-coded in the colour bar.
= t and =2 , (6)
dt (cid:107)δ(cid:107)2 dt t
and are solved in tandem alongside the equations of mo- 3 RESULTS AND DISCUSSION
tionandvariationalequationsofmotion.Attheendofeach
3.1 The Phase Space Evolution
integration step the quantities Y(t) and (cid:104)Y(cid:105)(t) can be de-
termined from Y(t) = 2x(t)/t and (cid:104)Y(cid:105)(t) = w(t)/t. When To observe the system that is evolving in the phase space,
generating MEGNO maps we chose to plot the mean value the planet’s eccentricity (e) is plotted versus the argument
(cid:104)Y(cid:105). of periapsis (ω), for various choices of the planet’s rela-
In general, MEGNO has the parameterisation (cid:104)Y(cid:105) = tive orbital inclination with the binary plane (Fig. 1). The
α×t+β (Cincotta & Simo´ 2000; Go´zdziewski et al. 2001). Savitzky-Golaysmoothingfunctionisusedtofilterthedata
For a quasi-periodic initial condition, we have α (cid:39) 0.0 and as we are interested in the secular dynamics of the sys-
β (cid:39) 2.0 (or (cid:104)Y(cid:105) → 2.0) for t → ∞ asymptotically. If the temandthesecularperturbationtheorywhichinitiatesthe
orbit is chaotic, then (cid:104)Y(cid:105) → λt/2 for t → ∞. Here λ is Kozaimechanisminthecaseofinclinedorbits.Wefindthat,
the maximum Lyapunov exponent (MLE) of the orbit. In below 39◦ and for our choice of i values, the phase space
pl
practice,whengeneratingourMEGNOmaps,weterminate starts to evolve where ω starts to circulate between 0◦ and
agivennumericalintegrationofachaoticorbitwhen(cid:104)Y(cid:105)> 360◦. The eccentricity, however, oscillates with large and
12.0. In contrast, quasi-periodic orbits have |(cid:104)Y(cid:105)−2.0| (cid:54) increasing amplitude for increasing i (for example 20◦ to
pl
0.001. For more details on the mathematical properties of 35◦)whilecomparedtothecoplanarcase(seeFigs.1a,b,c).
(cid:13)c 2002RAS,MNRAS000,1–??
Monthly Notices: LaTEX2ε guide for authors 5
(a)Initiali =10◦ (b)Initiali =30◦ (c)Initiali =39.2◦
pl pl pl
Figure 2. Planet’s eccentricity and inclination time series when planet’s initial orbital inclination with the binary plane is set at 10◦,
30◦ and 39.2◦. In these double-axis plots, the left y-axis has eccentricity and right y-axis has inclination plotted versus time along the
commonx-axis.Theorbitalintegrationwascarriedoutfor50Myrbut(a)and(b)aretruncatedat2×104 yearstoclearlydemonstrate
theoscillationsofbothelements(e andi).(c)isplotteduptotheinstabilitypoint(∼4×104 years)afterwhichtheplanetcollidedwith
thecentralbody.
(a)Initiali =50◦ (b)Initiali =60◦ (c)Initiali =83◦
pl pl pl
Figure 3. Planet’s eccentricity and inclination time series when planet’s initial orbital inclination with the binary plane is set at 50◦,
60◦ and83◦.Thetimeseriesat(a)and(b)aretruncatedat2×104 yearsandthetimeseriesat(c)isplotteduptotheinstabilitypoint
(∼1.5×104 years)afterwhichtheplanetcollidedwiththecentralbody.
Whentheinitiali issetat39.2◦,theω showscircula- 3.2 The Kozai Resonance
pl
tion only with eccentricity amplitude maximum at 90◦ and
270◦ (Fig. 1d), but the system becomes orbitally unstable
The evolution of a planet’s eccentricity and inclination was
as the e reaches its maximum value causing the planet to
pl furtherexploredforawiderangeofi values.Theconserva-
collide with the central body within 50 Kyr. Further incre- pl
tionoftheKozaiintegralterm(seeSec.3.2)suggeststhatas
ment inthei valueto 39.7◦ makes thephase space divide
pl thesystemevolvesintime,thedecrease/increaseinplanet’s
into two distinct regions where the ω is found to circulate
eccentricity is compensated by the increase/decrease in it’s
as well as librate (Fig. 1e). This is also a transition regime
inclination relative to the binary plane, also shown by Hol-
of the phase space beyond which the ω displays libration
man et al. (1997) for circular perturber and by Lithwick &
only. Holman et al. (1997) have shown a similar behaviour
Naoz (2011) for an eccentric perturber. For small initial i
in the momemtum variable (1-e2) plotted versus its conju- pl
values (i.e., 10◦) we found a constant-amplitude oscillation
gate angle (ω) for a highly inclined (i = 60◦) planetary
pl of e and i at the average values of 0.55 and 12◦, respec-
orbit. pl pl
tively(Fig.2a)fortheintegrationperiodof50Myr.Asthe
initialvalueofi isincreasedtheamplitudeofoscillationsof
pl
Furthermore,whenthei issetinincreasingorderfrom both parameters increase as well. In Fig. 2b we have shown
pl
45◦ to80◦,theamplitudesofthee oscillationsisfoundto a case for i = 30◦. In this double-axis plot, the y-axis on
pl pl
increaseaswell,andinthisinclinationregimetheargument theleftsideindicatesthetimeevolutionofeccentricityfrom
of periapsis only shows libration around 90◦ (Fig. 1f). The its best-fit value of 0.48. As the system evolves, the eccen-
librationamplitudeincreaseswiththeincreaseintheinitial tricity varies between 0.25 and 0.55. The y-axis on right
i values,eventuallyenteringaregimewheretheplanetcan side illustrates the variation in inclination from its initial
pl
escapefromthesystemorcollidewithitshoststarwheni assigned value of 30◦. The i value oscillates between 30◦
pl pl
reaches a value greater than 82◦ and e approaches 1.0 and40◦ andthevariationinitsamplitudefromtheinitially
max
(Eq. 1). assigned value behaves in an anti-correlated fashion to that
(cid:13)c 2002RAS,MNRAS000,1–??
6 Satyal et al.
of e . Also, the variation in the e agrees with the case
pl pl
when the perturber has an eccentric orbit, the inner orbit’s
eccentricity can reach extremely large values (Ford et al.
2000; Naoz et al. 2013; Teyssandier et al. 2013). The x-axis
inFigs.2aand2baretruncatedat20Kyrtoclearlydemon-
strate the oscillations of the two parameters, (e and i ).
pl pl
Theamplitudeoftheoscillationsmerelychangedduringits
evolution for given total integration time.
Whentheinitiali valuewasfurtherincreasedto39.2◦,
pl
theplanetfollowsapathtoinstabilitywithin40kyrthrough
a growing value of eccentricity. It is also found that there
is a small instability window (most of the time the planet
collided with the central body) when i is set at ∼(39◦ - Figure 4. A maximum eccentricity (emax) map of the planet,
pl HD 196885 Ab, for varying i and a , simulated for 50 Myr.
40◦). The time series plot of e and i (Fig. 2c) shows in- pl pl
pl pl The colour bar indicates the emax reached by the planet during
creasing values of the planet’s eccentricity and inclination
thetotalsimulationtime,whichalsoincludesthecaseswhenthe
until it reaches an extreme value ∼1 and i oscillates up
pl planetsuffersanejectionorcollision(especiallywhenepl reaches
to 80◦, hence minimizing the periastron distance and even- avaluegreaterthan0.9).Thedarkercolourrepresentsthebest-
tually forcing the planet to collide with its host star. From fit e parameter (0.48) and the lightest colour represents the
pl
a theoretical aspect it can be inferred that when the rel- emaxvaluetheplanetattainedfortherespectivechoicesofinitial
ative inclination hits the critical angle mark (ic = 39.2◦), conditions in ipl and apl. The red line at 2.6 au is the best-fit
thelong-periodoscillationsbetweeneccentricityandinclina- semimajoraxisoftheplanet.Thismapcomplimentsthelifetime
tionensue.Theinitialeccentricitybecomesinsensitivelead- map(Fig.6)inawaythattheverticalinstabilitybarsarelocated
ing to forced eccentricity, which is the basis for the Kozai atthesameregionswheretheplanet’semax exceeds0.7.
resonance. The other factor causing the planet to exhibit
Kozai cycles is the libration of the argument of periapsis
(ω) around 90◦, whose significant evolution is observed in maximum eccentricity reaches the extreme values, from 0.7
Fig. 1 for i greater than 39◦. We have not considered the to ∼ 1, (see the maximum eccentricity map, Fig. 4) for the
pl
tidaldissipationwithintheregimesofcloseapproachwhich orbital inclination greater than 55◦; nonetheless, the planet
mayaswellproducemisalignment(Winnetal.2010).Also, survives the total integration time until the i is raised to
pl
generalrelativityisnotincludedinourcalculationsthepre- 82◦ . The e stays relatively low, between 0.48 and 0.68,
max
cession of which can lead to large eccentricity oscillations for the i values less than 55◦ excluding the case when i
pl pl
(Naoz et al. 2013; Antonini & Merritt 2013). is∼39◦.Beforereachingtheflippingpoint(90◦)theplanet
Forgreaterorbitalinclinationvalues,suchas,whenini- gets ejected from the system when ipl (cid:62) 82◦ and the emax
tial i is set to 50◦ (or, 60◦) (Fig. 3a,b) the amplitudes of reachescloseto1.0withintheintegrationtime.Wenotethat
pl
eccentricity and inclination oscillations reach as high as ∼ emax appearstodecreaseinFig.4forveryhighinclinations
0.6and53◦(or,∼0.8and64◦),respectively,andthesystem (85◦−90◦), but this is due to the approach to instability
continuestomaintaintheperiodicorbitsthroughoutthein- occurring faster than our data output frequency.
tegration period of 50 Myr. Also, further increasing the i
pl
valuepushese toitsextremelimitresultingintheplanet-
pl
star collision when the i value is increased beyond 82◦.
pl 3.3 The MEGNO and the MMR Analysis
Above the critical value of the planet’s orbital inclination
withthebinaryplane(ic =39.2◦),theprecessionofthear- We have generated the MEGNO maps considering a wide
gument of pericentre is replaced by libration around 90◦ as parameterspaceoftheplanet’sinitialinclinationandsemi-
discussed earlier. Satyal et al. (2013) encountered a similar majoraxiswhichillustratesvariationsinthequasi-periodic
dynamicalbehaviourintheirstudyoftheplanetintheHD and chaotic phase space for varying i and a values (Fig.
pl pl
196885 system for the planet’s orbital inclination less than 5). We found that, for i less than 39◦, the best-fit value
pl
25 degrees. For higher initial ipl, as shown here, the secular of apl (the vertical line at 2.6 au) lies within the chaotic
perturbation causes the epl and ipl values to oscillate with region. However, it’s important to note that a chaotic orbit
higheramplitudes,whichcaneventuallycauseinstabilityin does not necessarily imply an unstable one. This is shown
the system due to the collision of the planet with the cen- above in Sec. 3 through the orbital integration of select ini-
tral body or escape due to a parabolic or hyperbolic orbit. tial conditions. Also, the global lifetime map of i vs. a
pl pl
Forexample,whentheinitialipl =83◦ andtheepl isforced (Fig.6)showswhichinitialconditionssurviveforthefullin-
beyond1withinashortintegrationperiod(15kyr,Fig.3c). tegrationtime.Thecolourcodeinthelifetime mapisbased
For specific initial conditions, Lithwick & Naoz (2011) ontheplanet’ssurvival,ejectionorcollisiontimewherethe
haveshownthattheplanetaryorbitcanflipwhenthecrite- lightestcolourrepresentsthefullsurvivaloftheplanet.And
rion given in Eqn. 3 is satisfied. But the initial parameters the dark coloured vertical bands signifies instability, which
from our system does not satisfy this criterion despite the means the planet at these initial conditions did not survive
presenceofthehighlyeccentricperturber.Forthegivenini- the total simulation time of 50 Myr.
tialparametersofHD196885AB,(cid:15)=0.063((cid:15)isaconstant The mean motion resonance (MMR) associated with
term given by Eqn. 2) while the value obtained from Eqn. the chaotic regions in the MEGNO maps can be estimated
3 is 0.174. Thus, the planet is not expected to flip its orbit bycalculatingitspositionfromtheperturbationtheory.For
which is in agreement with Fig. 3c and Fig. 4. The planet’s thebest-fitvaluesoftheplanet’sandbinarysemimajoraxes,
(cid:13)c 2002RAS,MNRAS000,1–??
Monthly Notices: LaTEX2ε guide for authors 7
Figure 5.MEGNOmapsindicatingthequasi-periodicandchaoticregionsofthegiantplanetintheHD196885systemforinclination
varyingfrom(left)0◦to90◦andsemimajoraxisrangingfrom[2.3-2.9]au,simulatedfor2×105years.Verticallineindicatesthelocation
ofthebest-fitsemimajoraxisoftheplanetata =2.6au.Asmallareaalongtheverticallineandforloweri regionismagnifiedto
pl pl
clearlydifferentiatebetweenthequasi-periodicandchaoticregions(rightpanel)andsimulatedfor8×105 years.Theplanet’snear 39:2
MMRislocatedslightlyleftoftheplanet’sbest-fitlocation.TheMMRstopswheni reaches∼30◦,withasmallsemiquasi-periodic
pl
regionextendingfrom30◦ to34◦ whichweexpecttovanishforlongersimulationtime.Yellowregionsignifiesthechaoticorbitswithin
the total simulation time and dark purple region signifies the periodic orbits. The colour bar indicates the strength in the value of
MEGNO,<y>.
andmassesofthestars,theratioofthe(p+q)/pcommensu-
rabilitycanbecalculatedas(seeMurray&Dermott(1999)),
(cid:18) (cid:19)3 (cid:18) (cid:19)1
p+q = abin 2 mA 2, (7)
p a m +m
pl A B
Wherem andm representmassesoftheprimaryand
A B
secondary star respectively; q is the order of the resonance
and p is its degree. To search for the chaos associated with
(p+q)/pMMR,weconsideredtheresonantangle(Φ)ofthe
form
Φ=k λ +k λ +k (cid:36) +k (cid:36) +k Ω +k Ω ,
1 bin 2 pl 3 bin 4 pl 5 bin 6 pl Figure 6. A global dynamical lifetime map of the planet, HD
(8)
196885 Ab, for varying i and a , simulated for 50 Myr. The
pl pl
Wherethecoefficientsfollowtherelationsk =p+q,k colourbarindicatesthesurvivaltime,wheredarkercolourrepre-
1 2
(cid:80) sentstheinstability(ejectionorcollision)andthelightestcolour
=-p,and k =0.Themeanlongitudeλisafunctionof
i i representthestability(survival)uptotheintegrationperiod.The
themeananomaly,M,andlongitudeofpericenter,(cid:36).Also,
solidredlineat2.6auisthebest-fitsemimajoraxisoftheplanet.
thesumofk andk isevenasrequiredbysymmetryinthe
5 6 Two dotted lines indicate the observational uncertainty of ±0.1
d’Alembert rules.
au (see Table 1). The vertical bars evident in Figs. 4 and 5 at
Fortheplanetat2.6aufromtheprimarycomponentan 2.45,2.55,2.65,2.75aualsoappear.
interactionwiththemeanmotionresonance(MMR)isfound
tobenear a39:2commensurabilitywiththesecondarycom-
ponent(someoftheanglesinFig.7showatemporarylock). Initial conditions for different i values, but along the
pl
The time evolution of the angles (cid:36) −(cid:36) (secular reso- best-fitsemimajoraxisline,werepickedwherethemeasure
pl bin
nance,Fig.7,f)showscirculationbetween0◦and360◦when of MEGNO is chaotic (Fig. 5, right panel) and the orbits
thei issetatlowervalues(forexample0◦)whichsuggests were integrated for a short period with a high data sam-
pl
thattheobservedchaosiscausedbytheMMR.Secularres- pling frequency. We would expect the resonance angle (Φ)
onance at the semi quasi-periodic region (i = 35◦) shows time series to alternate between modes of circulation and
pl
circulationaswell,butwithadifferentorientation(Murray librationforinitialconditionstakenfromthechaoticregion
& Dermott 1999). Then, the MMR value is used to search and Φ to circulate only when the initial condition is taken
foralibratingresonantangle(Eq.8),fork =39,k =-2, from the quasi-periodic region. At 0◦, 10◦, 31◦ and 35◦, Φ
1 2
and various other combinations of k , k , k and k due to is found to librate and circulate as expected for chaotic be-
3 4 5 6
theinfluenceofsecularcomponents.Inourcase,thei was haviour (Fig. 7, a-d). This also confirms the earlier indi-
pl
variedwithrespecttothebinaryorbitalplaneandtheorbits cation of the planet’s interaction with a near 39:2 MMR
were integrated using astrocentric coordinates, considering with the secondary star. One of the factors inducing chaos
theprimaryattheorigin.Thus,k inEq.8wouldbezeroas in the region of low inclination orbit is possibly due to this
5
Ω becomesundefined(i =0◦)andconventiondictates near39:2MMR.TheMMRinthecaseofhigherinclination
bin bin
that it be set to zero. orbit changes its dynamical character which results in the
(cid:13)c 2002RAS,MNRAS000,1–??
8 Satyal et al.
(a)Planet’sInclination(i )=0◦ (b)Planet’sInclination(i )=10◦ (c)Planet’sInclination(i )=31◦
pl pl pl
(d)Planet’sInclination(i )=35◦ (e)Planet’sInclination(i )=50◦ (f)Planet’sInclination(i )=0◦,35◦
pl pl pl
Figure 7. The time evolution of resonant angles (Φ) calculated for different i values but along the best-fit location of semimajor
pl
axis (a = 2.6 au). The location for Φ angle test is chosen from MEGNO map (see Fig. 5) where measure of MEGNO (< y >) is
pl
quasi-periodic(e)andchaotic(a,b,candd).(f)showsthetimeevolutionoftheangles((cid:36) −(cid:36) )fori valuesat0◦ and35◦.
pl bin pl
periodic orbits. The location of the planet (vertical line at choiceofa ,from2.3auto2.9au.Theinstabilityatthisre-
pl
2.6au)isalmostatthecenteroftheMMR,thusminimizing gionisalsoconfirmedbythelifetime mapwheretheplanet
the amplitude of the libration. The 39:2 MMR has a min- isfoundtocollidewithitshoststar.TheMEGNOhasbeen
imal effect for i greater than 40◦ and the planet enters a successfulindeterminingthechaoticregionsurroundingthis
pl
quasi-periodic region for i up to 55◦, a region where the instability point. The minimum critical angle (i (cid:62) 39.2◦)
pl pl
e deviatesleastfromitsbest-fitvalue.Theresonantan- required to initiate the Kozai mechanism lies within this
max
gleisfoundtocirculateonlyinthesequasi-periodicregions. inclination regime. Once the Kozai resonance phenomenon
Φplottedforinitialconditionsat50◦ hasshowncirculation is triggered, it drives the planet into higher inclination or-
in Fig. 7e. This gives a strong evidence that one of the fac- bits(hencehighereccentricity)causingtheplanettocollide
tors responsible for the chaos below 39◦ is produced by the withtheprimary.However,beyond40◦andlessthan55◦the
near 39:2 MMR interaction. Also, the secular perturbation planetmaintainsquasi-periodicorbitseventhoughitissub-
inducedbytheprecessionofω andΩisfoundtocontribute ject to the Kozai resonance. The observed horizontal stripe
totheobservedchaoticdynamicsofthesystemintheKozai intheMEGNOmapati ∼39◦canbeconsideredasasep-
pl
regime.Thelibrationandcirculationofresonantangleseen aratrixlinewhichseparatesthephasespacefromcirculation
in Fig. 7 is obtained only when the coefficients of (cid:36) and when i < 39.2◦, and libration when i (cid:62) 39.2◦.
pl pl pl
Ω are included in the Φ term.
pl TheMEGNOmapshavedisplayedquasi-periodicorbits
A small region of possible quasi-periodicity appears along the best-fit a value ( black line, Fig. 5) and for i
pl pl
when the initial i lies between 30◦ and 34◦ (Fig. 5 right values greater than 40◦ and until it rises as high as ∼55◦
pl
panel).Thiswouldseemlikeaperiodicregionbutthereso- abovethebinaryplane.Beyondthisinclinationregime,the
nance angle test (Fig. 7c) shows that the Φ is found to cir- maps have displayed the chaotic orbits. The MEGNO have
culate and librate, which indicates this region as a chaotic shown an overall effectiveness at indicating the MMR loca-
zone. The small black and red dots in this region are in- tions, and the chaotic and periodic regions of the system’s
dicative that many of these initial conditions have not had phasespace.Inaddition,thedynamicalmaximumeccentric-
enough time to reveal their true nature (quasi-periodic or ity and lifetime maps have addressed the caveats of chaos
chaotic).Whatwecalledbeforeaquasi-periodicregionand seen in the MEGNO maps to discriminate which regions
expected the circulation of Φ has revealed itself as a region are truly unstable or merely chaotic within the integration
ofcomplexdynamicswithmanyoverlappingandinteracting period.
resonances.
The planet’s interaction near the 39:2 MMR with the
Astripofunstablechaosisobservedwhentheinitiali secondaryisobservedwhenthemeanbest-fitvalueofsemi-
pl
rangesbetween∼(39◦-40◦)(Fig.5,rightpanel)andforany major axis is 2.6 au. When the planet’s location is moved
(cid:13)c 2002RAS,MNRAS000,1–??
Monthly Notices: LaTEX2ε guide for authors 9
to2.58auand2.68au(butwithintheobservationaluncer-
tainty limit), the planet is found to exhibit oscillations due
to the 20:1 and 19:1 MMRs, respectively.
Theplanet’sorbitalmotionisalsoanalysedinthevary-
ingi ande plane(Fig.8)todisplaypossibleperiodicand
pl pl
chaoticregimes.Theverticallineate =0.48indicatesthe
pl
planet’sbest-fitlocation,and,likeinthepreviouscase(Fig.
5), the planet clearly lies in chaotic zone for e > 0.4 and
pl
i < 30◦. Driven by the Kozai resonance, a similar chaotic
pl
stripe appears at i ∼ 40◦. Periodic regions continue be-
pl
yondthisangle,howeversomechaoticislandsdoappearfor
i greater than 55◦.
pl
The discovery of the planet in the HD 196885 system
(Chauvin et al. 2011) using the radial velocity (RV) tech-
nique has only allowed the determination of the planet’s
minimum mass of 2.98M . We have performed a test to in- Figure 8. MEGNO maps indicating the periodic and chaotic
J
orbits of the giant planet in HD 196885 for inclination varying
dicatetheperiodicand/orchaoticregionsofthesystemfor
from0◦ to90◦,plottedversustheeccentricityandsimulatedfor
variations in the planet’s mass (m ) ranging from 2M to
pl J 2×105 years.Verticallineindicatestheplanet’sbest-fitlocation
20M and i varying from 0◦ to 90◦ (Fig. 9). The map
J pl ate =0.48.
pl
indicates that quasi-periodic regions exist for various com-
binationofm andi .Thesystemdisplaysperiodicorbits
pl pl
form mostlybetween3M and6M withasmallchaotic
pl J J
regionaround4.5M plustheexclusionofchaoticregionat
J
∼ 40◦. The chaotic region continues below 30◦ at the best-
fitsemimajoraxislocationoftheplanet.Themapindicates
otherpossibleperiodicislandsforhighermassvaluesofthe
planet specifically for i between (40◦ - 50◦) and m be-
pl pl
tween ∼ (8.5 - 11.5)M , (12.5 - 15)M and (16 - 19)M .
J J J
4 CONCLUSIONS
Basedonouranalysis,thebestconfigurationfortheplanet
to display quasi-periodic orbits is when the planet’s orbital
inclinationrelativetothebinaryplaneliesbetween41◦ and
55◦. Chauvin et al. (2011) made a similar suggestion re- Figure 9. MEGNO maps indicating the periodic and chaotic
orbits of the giant planet in HD 196885 for inclination varying
gardingtheplanet’sinclinationthata highrelativeinclina-
from0◦to90◦,plottedversustheplanet’smassandsimulatedfor
tion favours for the periodic orbits. Also, Giuppone et al.
2×105 years.Verticallineindicatestheplanet’sbest-fitlocation
(2012)havefoundthattheplanethaseithercoplanarorbits
atmpl =2.98MJ.
(progradeandretrograde)orhighlyinclinedorbitsnearthe
Lido-Kozai equilibriumpoints,i-90◦=±47◦.Thisrangeof
anglesgetssmallerwhentheplanet’smassisincreasedmore have forced the planetary embryo or protoplanet into pre-
than the current value of 2.98M . The planet can have an cession during its early formation stage, eventually driving
J
interaction with the near 39:2 MMR from perturbations of the planet into high inclination orbits.
thesecondarycomponentofthebinarywhichappearstobe The planets higher mass is possibly constrained be-
responsible for the chaotic region as seen in the MEGNO tween (3 - 6)M ; however, these choices also depend on
J
maps below 39◦. The plots of the resonance angle time se- the choice of planet’s orbital inclination and it is mostly
riesforthenear39:2MMRclearlydemonstratethecircula- favorable when m is less than 6M and i less than
pl J pl
tionandlibrationofΦbasedonthechoicesofinitialcondi- ∼ 55◦, with some exceptions of chaotic islands (Fig. 9).
tion. The planet’s orbital configuration below 39◦ is proven The planetary mass higher than 9M and up to 19M is
J J
chaotic, and hence lesser possibility for the planet to be in likely when the planet’s orbital inclination lies in a smaller
thatregimeunderthecurrentassumptionsofplanetaryfor- regime, somewhere between 40◦ and 50◦. A possibility of
mation. havinganyotherterrestrialplanetswithinthe2.6auregion
Thebest-fitvalueoftheplanet’seccentricitywassetat around the primary is less likely. The circumprimary disc
0.48.Thenthesystemwasevolvedfor50Myrandthemaxi- can be strongly hostile to planetesimal accretion in this
mume valuewascalculatedina vs.i parameterspace. region (Thebault 2011).
pl pl pl
Thee isfoundtodeviateleastfromitsbest-fitvaluewhen
pl
i is at ∼35◦ and ∼(40◦-50◦), thus, these regions are ex- Acknowledgement:SSandJNwouldliketothankThe
pl
pected to demonstrate the most quasi-periodic regimes for DepartmentofPhysicsatUTArlington,ZdzislawMusielak
the planet. The high amplitude oscillation in e below 35◦ and Manfred Cuntz for their continuous support and guid-
pl
mainly arises due to the planet’s secular interactions with ance. BQ gratefully acknowledges support from the NASA
the secondary. The perturbation from the secondary may postdoctoralprogram.TCHgratefullyacknowledgesfinan-
(cid:13)c 2002RAS,MNRAS000,1–??
10 Satyal et al.
cial support from the Korea Research Council for Funda- MikkolaS.,InnanenK.,1999,CelestialMechanicsandDy-
mentalScienceandTechnology(KRCF)throughtheYoung namical Astronomy, 74, 59
Research Scientist Fellowship Program and financial sup- MurrayC.D.,DermottS.F.,1999,Solarsystemdynamics
portfromKASI(KoreaAstronomyandSpaceScienceInsti- Nagasawa M., Ida S., Bessho T., 2008, ApJ , 678, 498
tute)grantnumber2014-1-400-06.Numericalcomputations NaozS.,FarrW.M.,LithwickY.,RasioF.A.,Teyssandier
werepartlycarriedoutusingtheSFI/HEAIrishCentrefor J., 2011, Nature , 473, 187
High-EndComputing(ICHEC)andthePLUTOcomputing NaozS.,FarrW.M.,LithwickY.,RasioF.A.,Teyssandier
clusterattheKoreaAstronomyandSpaceScienceInstitute. J., 2013, MNRAS , 431, 2155
AstronomicalresearchattheArmaghObservatoryisfunded Naoz S., Farr W. M., Rasio F. A., 2012, ApJL , 754, L36
by the Northern Ireland Department of Culture, Arts and Naoz S., Kocsis B., Loeb A., Yunes N., 2013, ApJ , 773,
Leisure(DCAL).Wewouldalsoliketothanktheanonymous 187
referee for his/her comments and suggestions on eccentric NoyolaJ.P.,SatyalS.,MusielakZ.E.,2013,ArXive-prints
Lidov-Kozai mechanism which significantly improved this Quarles B., Eberle J., Musielak Z. E., Cuntz M., 2011,
work. A&A, 533, A2
Rasio F. A., Ford E. B., 1996, Science, 274, 954
Ricker G. R., Latham D. W., Vanderspek R. K., Ennico
K. A., Bakos G., Brown T. M., Burgasser A. J., Char-
REFERENCES
bonneauD.,ClampinM.,DemingL.D.,DotyJ.P.,Dun-
Antonini F., Merritt D., 2013, ApJL , 763, L10 ham E. W., Elliot J. L., Holman M. J., Ida S., Jenkins
Chambers J. E., 1999, MNRAS , 304, 793 J. M., 2010, in American Astronomical Society Meeting
Chambers J. E., Migliorini F., 1997, in AAS/Division for Abstracts #215 Vol. 42 of Bulletin of the American As-
PlanetarySciencesMeetingAbstracts#29Vol.29ofBul- tronomicalSociety,TransitingExoplanetSurveySatellite
letin of the American Astronomical Society, Mercury - A (TESS). p. 450.06
New Software Package for Orbital Integrations. p. 1024 Satyal S., Quarles B., Hinse T. C., 2013, MNRAS , 433,
ChatterjeeS.,FordE.B.,MatsumuraS.,RasioF.A.,2008, 2215
ApJ , 686, 580 Sl(cid:32)oninaM.,Go´zdziewskiK.,MigaszewskiC.,2012,inAre-
Chauvin G., Beust H., Lagrange A.-M., Eggenberger A., nou F., Hestroffer D., eds, Proceedings of the workshop
2011, A&A , 528, A8 ”Orbital Couples: Pas de Deux in the Solar System and
Cincotta P., Simo´ C., 1999, Celestial Mechanics and Dy- theMilkyWay”.HeldattheObservatoiredeParis,10-12
namical Astronomy, 73, 195 October2011.Editors:F.Arenou,D.Hestroffer.ISBN2-
Cincotta P. M., Simo´ C., 2000, aaps, 147, 205 910015-64-5, p. 125-129 Mechanic: a new numerical MPI
Correia A. C. M., Laskar J., Farago F., Bou´e G., 2011, framework for the dynamical astronomy. pp 125–129
Celestial Mechanics and Dynamical Astronomy, 111, 105 Sl(cid:32)onina M., Go´zdziewski K., Migaszewski C., 2014, ArXiv
CuntzM.,QuarlesB.,EberleJ.,ShukayrA.,2013,PASA e-prints
, 30, 33 SzenkovitsF.,Mako´Z.,2008,CelestialMechanicsandDy-
Fabrycky D., Tremaine S., 2007, ApJ , 669, 1298 namical Astronomy, 101, 273
FordE.B.,KozinskyB.,RasioF.A.,2000,ApJ ,535,385 Teyssandier J., Naoz S., Lizarraga I., Rasio F. A., 2013,
Giuppone C. A., Morais M. H. M., Bou´e G., Correia ApJ , 779, 166
A. C. M., 2012, A&A , 541, A151 ThebaultP.,2011,CelestialMechanicsandDynamicalAs-
Go´zdziewski K., Bois E., Maciejewski A. J., Kiseleva- tronomy, 111, 29
Eggleton L., 2001, A&A , 378, 569 Veras D., Ford E. B., 2010, ApJ , 715, 803
Go´zdziewskiK.,MaciejewskiA.J.,2001,ApJL,563,L81 WinnJ.N.,FabryckyD.,AlbrechtS.,JohnsonJ.A.,2010,
Heller R., 2012, A&A , 545, L8 ApJL , 718, L145
Hinse T. C., Christou A. A., Alvarellos J. L. A., Wu Y., Murray N., 2003, ApJ , 589, 605
Go´zdziewski K., 2010, MNRAS , 404, 837 Wu Y., Murray N. W., Ramsahai J. M., 2007, ApJ , 670,
Hinse T. C., Michelsen R., Jørgensen U. G., Go´zdziewski 820
K., Mikkola S., 2008, A&A , 488, 1133
Holman M., Touma J., Tremaine S., 1997, Nature , 386,
254
Innanen K. A., Zheng J. Q., Mikkola S., Valtonen M. J.,
1997, AJ , 113, 1915
KatzB.,DongS.,MalhotraR.,2011,PhysicalReviewLet-
ters, 107, 181101
KippingD.M.,HartmanJ.,BuchhaveL.A.,SchmittA.R.,
Bakos G. A´., Nesvorny´ D., 2013, ApJ , 770, 101
Kozai Y., 1962, AJ , 67, 591
Li G., Naoz S., Kocsis B., Loeb A., 2013, ArXiv e-prints
Lidov M. L., 1962, PlanetaryandSpaceScience , 9, 719
Lithwick Y., Naoz S., 2011, ApJ , 742, 94
Lyapunov A., 1907, Annales de la facult`e des sciences de
Toulouse, 2:9, 203
Mayor M., Queloz D., 1995, Nature , 378, 355
(cid:13)c 2002RAS,MNRAS000,1–??