Stability of planets in triple star systems
AAstronomy & Astrophysics manuscript no. A_A_33097_Busetti c (cid:13)
ESO 2018November 21, 2018
Stability of planets in triple star systems
F. Busetti , H. Beust , and C. Harley Univ. of the Witwatersrand, CSAM, Private Bag 3, 2050-Johannesburg, South Africae-mail: [email protected] Univ. Grenoble Alpes, CNRS, IPAG, F-38000 Grenoble, Francee-mail: [email protected]
Received 27 March 2018 / Accepted 27 August 2018
ABSTRACT
Context.
Numerous theoretical studies of the stellar dynamics of triple systems have been carried out, but fewer purely empiricalstudies that have addressed planetary orbits within these systems. Most of these empirical studies have been for coplanar orbits andwith a limited number of orbital parameters.
Aims.
Our objective is to provide a more generalized empirical mapping of the regions of planetary stability in triples by: consideringboth prograde and retrograde motion of planets and the outer star; investigating highly-inclined orbits of the outer star; extending theparameters used to all relevant orbital elements of the triple’s stars and expanding these elements and mass ratios to wider ranges thatwill accommodate recent and possibly future observational discoveries.
Methods.
Using N -body simulations we integrated numerically the various four-body configurations over the parameter space, usinga symplectic integrator designed specifically for the integration of hierarchical multiple stellar systems. The triples were then reducedto binaries and the integrations repeated, to highlight the di ff erences between these two types of system. Results.
This established the regions of secular stability and resulted in 24 semi-empirical models describing the stability boundsfor planets in each type of triple orbital configuration. The results were then compared with the observational extremes discovered todate, to identify regions that may contain undiscovered planets.
Key words. methods: numerical – methods: N -body simulations – planet-star interactions – celestial mechanics – stars: hierarchicaltriples – planetary systems: dynamical evolution and stability
1. Introduction
A large portion of stellar systems are composed of multiple stars.Single stars account for approximately half, binaries a third, andtriples a twelfth (Tokovinin 2014a,b). Over more than two hun-dred years there have been numerous theoretical studies of thethree-body problem and its restricted form. A few recent the-oretical approaches to the stellar dynamics of triple systems,some incorporating numerical checks of their results, includeFord et al. (2000); Mardling & Aarseth (2001); Hamers et al.(2015) and Correia et al. (2016). Theoretical approaches can,however, have limitations in accurately representing some actualsystems.Numerical studies are a powerful tool with which to map theoften discontinuous stability boundaries of systems without anyrestriction on their physical characteristics and configurations.Also, a number of studies have shown that, for systems that con-tain more than one planetary body, the orbits proposed initiallywere simply not dynamically feasible, which has promoted dy-namical analyses as an integral part of the exoplanet discoveryprocess (Horner et al. 2012).A hierarchical triple stellar system consists of three bodiesorganized in two nested orbits — one central binary and a thirdstar orbiting the center of mass of this binary at a larger distance.We shall call the two bodies constituting the binary stars 1 and2, and the outer body star 3. We note that we have made no as-sumptions about the relative masses of the three stars. In manysystems, star 3 is less massive than the central binary, but thereare examples of systems with inverted mass ratios, such as HD 181068, for which star 3 is heavier than the binary (Derekas et al.2011; Borkovits et al. 2013).There are four possible types of planetary orbits in such asystem (see Fig. 1):- S1, S2, and S3 orbits comprise planets that orbit one of theindividual stars (stars 1, 2 and 3 respectively) while being per-turbed by the other stars. S1 and S2 orbits are not distinguishableand will be treated as a single S1 type;- P1 orbits comprise circumbinary planets that orbit the cen-tral binary while being perturbed by the outer third star; and- P2 orbits comprise circumtriple planets that orbit the wholesystem beyond star 3.Planets have been found in 32 triple systems to date. Ofthese, S3 orbits have been found in 29 systems, S1 orbits in twoand a P1 orbit in one. No planet in a circumtriple P2 orbit hasyet been discovered.Numerical studies of planets in multiple systems have fo-cused on binaries. For example, in their key paper, Holman &Wiegert (1999) examined S-type and P-type planetary orbits ineccentric, coplanar binary systems and derived regression equa-tions describing the stability bounds in terms of mass ratiosand eccentricities. Eccentric and inclined (up to 50 ◦ ) P-type or-bits, for equal-mass binaries, are studied by Pilat-Lohinger et al.(2003). Work on S-type and P-type orbits is done by Musielaket al. (2005) for circular, coplanar orbits, defining the regions ofstability as a function of the semi-major axis ratio between thestar and planet and the mass ratio. Mudryk & Wu (2006) alsoanalyze the numerical results of Holman & Wiegert (1999) forS-type orbits in coplanar systems and investigated the instability Article number, page 1 of 12 a r X i v : . [ a s t r o - ph . E P ] N ov & A proofs: manuscript no. A_A_33097_Busetti boundary at a higher resolution. A binary study that addressed afull range of planetary inclinations, for eccentric P-type orbits,was by Doolin & Blundell (2011), which described stability as afunction of mass ratios and eccentricities.One study that addressed retrograde planets in binary, copla-nar, circular, S-type orbits, was by Morais & Giuppone (2012),while Giuppone & Correia (2017) looked at S-type orbits, in-cluding retrograde ones, in a compact binary system, for vari-ous inclinations, eccentricities and orientation angles. Althoughnone have yet been confirmed, planets in retrograde orbits withinand around a triple system should be possible, similar to HAT-P-7b in its binary system.One study of the four-body problem of planets within a triplesystem was by Verrier & Evans (2007). This study was confinedto prograde, coplanar P-type orbits and derived regression rela-tionships analogous to those of Holman & Wiegert (1999).In a more general analysis, in addition to retrograde plane-tary orbits it is also necessary to include highly-inclined orbits ofthe outer star, since the resulting stellar Kozai resonance sculptsthe geometry of the stable planetary region (Kozai 1962; Fordet al. 2000). It is also useful to consider retrograde motion of theouter star. For the 22 triple systems where mutual inclinationshave been established, the mean is a high 63 ◦ and four of thesesystems are retrograde (Tokovinin 2017). One analytical studythat addressed retrograde stellar orbits is by Farago & Laskar(2010). They address two limiting cases, the inner restrictedproblem, where one inner body (i.e., our star 2) has no mass,and the outer restricted problem, where the outer body (star 3)has negligible mass, describing the possible motions of particlesin each regime and providing an expression for the boundary be-tween the regimes.We extended the number of parameters used in previouswork to all the orbital elements of a triple’s stars and widenedthe ranges of both these and the relevant mass ratios, to encom-pass recent observations as well as potential future discoveries.For example, systems with highly inverted mass ratios, such asHD 181068, have been discovered but not yet modeled.The objective of this work was to provide a broader map-ping of the regions of planetary stability in generalized triples.This resulted in 24 semi-empirical models describing the stabil-ity bounds for planets in each possible combination of orbitalmotions in triples and for a range of stellar orbital elements andmass ratios.
2. Methods
Our aim was to determine the orbital stability of both the hoststars and their planets, where collisions between the stars orthe ejection of one of the three stars (typically the least mas-sive body) does not occur over secular timescales that are verylong compared to the orbital periods. Our treatment of secularperturbations was based on classical Newtonian dynamics andassumed that all three bodies are point masses that do not inter-act other than through gravitation and do not evolve in any way.Relativistic and tidal e ff ects, as well as spin, were not addressed.The triple system’s nomenclature is shown in the schematic il-lustration in Fig. 1.There are two S-type orbits, as the S1 and S2 orbits are equiv-alent, and two P-type orbits, one circumbinary and one circum-triple. The stable region for circumbinary (P1) orbits will havean inner and outer bound, while that for circumtriple (P2) orbitswill have an inner bound only. Circumstellar (S1 and S3) regions m m m a a a io a oi P1 orbits
Star 1 Star 3
S2 orbits S3 orbits
Star 2
P2 orbits a io a io S1 orbits
Fig. 1.
Triple system with S-type and P-type orbits. of stability will have their outer edges bounded. The objective isto determine the outer bound a io of the stable regions for S1,S3, and P1 orbits and the inner bound a oi of the P2 orbits. Thefollowing subsections list the selected parameter ranges for thestellar configurations.2.1.1. Semi-major axis ratio a / a The lower limit was selected to be just inside the stability cri-terion for triples as given by Mardling & Aarseth (2001). Thesmallest ratio used was approximately three, lower than thesmallest ratio of 4.1 found in the Eggleton & Tokovinin (2008)catalog of 285 triples and well above the a / a (cid:29) µ = m / ( m + m )The range for this ratio is between zero and one. However, onlythe range 0 < µ < . µ and (1 − µ ) are equivalent, other than for a 180 ◦ change in thelongitude of the ascending node. An upper limit of one and lowerlimits of 0.001 and 0.1 were used for P and S orbits respectively.2.1.3. Outer mass ratio µ = m / ( m + m )Outer mass ratio definitions vary widely. In the Tokovinin(2014a) survey, the largest mass ratio (defined there as m / m )found was 8.9 (for HIP 29860). However, 94% of mass ratiosare below two and 99% are below five. Taking an upper limit offive and assuming that the masses of the binary pair are broadlycomparable, the equivalent value for our ratio definition is 2.5.A mass ratio greater than one implies an inverted system, wherethe outermost star is more massive than the aggregate inner bi-nary. An example is the triple system HD 181068, which has amass ratio of ∼ e and e The only data on inner and outer stellar eccentricities in triplesis from Sterzik & Tokovinin (2002). The mean eccentricity of allorbits is high at 0.39 and in most (70%) of the systems the in-ner orbit is more eccentric than the outer orbit. The di ff erence ineccentricities within these systems ranges from e ff ectively zeroto 0.57. A study of 222 Kepler triples finds outer eccentricities Article number, page 2 of 12. Busetti et al.: Stability of planets in triple star systems spanning the full range, with a broad peak in the middle of therange and an unexplained narrow peak near e (cid:39) .
28 (Borkovitset al. 2015). The integrations therefore covered eccentricities inboth inner and outer orbits ranging from zero to 0.7–0.9 to coverfully this parameter range, with the higher limit required whenKozai resonance occurs.2.1.5. Inner and outer inclinations i and i The inner inclination i was always set to zero, so the outer in-clination i was the mutual inclination i − i of the outer or-bit relative to the inner orbit. For the low-inclination simula-tions, small ranges of mutual inclination were selected, of 0 ◦ –60 ◦ for prograde orbits and 120 ◦ –180 ◦ for retrograde orbits. Inthe high-inclination integrations the full range of 0 ◦ –180 ◦ wasused. Borkovits et al. (2015) found that the distribution of mutualinclinations for 62 Kepler triples had a large peak at 0 ◦ –5 ◦ , indi-cating close-to-coplanar configurations, with a significant 38%portion of the systems in a secondary peak centred at 40 ◦ , sug-gesting Kozai e ff ects.2.1.6. Other orbital elementsThe outer star’s longitude of ascending node Ω and argumentof periapsis ω were both varied from 0 ◦ to 360 ◦ while those forthe inner star were set to zero.2.1.7. Test particlesPlanets were represented by massless test particles. Their ini-tial orbital elements ranged between selected lower and upperlimits — for eccentricity these were 0–0.9, for inclination 0 ◦ –90 ◦ or 90 ◦ –180 ◦ for prograde or retrograde orbits respectivelyand for ω , Ω and mean anomaly M they were 0 ◦ –360 ◦ . Mostinitial orbital elements were randomly generated from uniformdistributions. However, for inclination i the direction of angularmomentum should be uniform over the celestial sphere, whichrequires a uniform distribution in cos i between [-1,1].A log scale was used for semi-major axes. The particles’semi-major axes ranged between 0.02 AU and 0.9 a for S1 orbitsand 0.02 AU and f a for S3 orbits, where f varied between 0.07and 0.60 depending on the stellar configuration. For P1 and P2orbits the lower limit was 0.01 AU while the outer limit was setconservatively at twice the distance at which the 5:1 mean mo-tion resonance occurred. The planet orbiting closest to its centralbody is PSR 1719-14 b, at 0.0044 AU, while Kepler-42c has theclosest orbit to a “normal” star, at 0.006 AU. The distance atwhich a test particle was stopped as being too close to the cen-tral body was set at 0.002 AU. The minimum semi-major axisused must be larger than this, and so 0.01–0.02 AU was selectedfor all orbit types. The furthest planet from its host star yet dis-covered is HIP 77900 b, at 3 200 AU. The distance at which atest particle was assumed to have escaped from the central bodywas set at 10 000 AU. /
20 of the orbital period of the in-ner binary was normally used in the integrations, and varied toverify that results were not a ff ected by numerical errors. Sym-plectic integration schemes usually ensure energy conservationwith 10 − –10 − relative accuracy using this time step size (Lev-ison & Duncan 1994; Beust 2003; Verrier & Evans 2007). Wefound an overall fractional change in the system energy ∆ E / E of ∼ − over a 10 yr integration, although this varied widelydepending on the triple configuration being integrated. Orbitalinstabilities tended to occur quickly – if an orbit did not becomeunstable in as little as 10–100 yr, it was usually stable up to 10 yr. The fact that the stellar systems we investigated were largelycompact was helpful, as 10 yr was often equivalent to manyhundreds of thousands of orbits of the outer star.The secular evolution of orbital parameters was sampled foreach of the configurations used, and an integration time of 10 yr was found to be su ffi ciently long in most cases. Each batchof integrations included at least one run on the least compactconfiguration with an integration time of 10 yr or 10 yr, toconfirm this.The number of planetary test particles used varied from 1 000to 10 000, with higher numbers being used when the rate of insta-bility and hence removal of test particles was high. For coplanarconfigurations, the initial test particle cloud took the form of adisk aligned with the plane of the inner binary, while for con-figurations where the outer star had a high inclination the testparticle disk was aligned with the invariable plane. A procedurefor automatically identifying the edges of the test particle cloudat the end of an integration and measuring its semi-major axiswas incorporated into the HJS algorithm.
3. Results and discussion
For inner (P1) and outer (P2) orbits, the respective outer and in-ner edges of the cleared area of the test particle cloud were stan-dardized by taking them as ratios of the semi-major axis of theouter star, a . These standardized dependent variables a io / a and a oi / a delineate the bounds of the stable regions and are termedthe inner and outer critical semi-major axis ratios respectively.The mean critical semi-major axis ratio for these inner andouter orbits, for the eight possible combinations of orbital mo-tion, are shown in Table 1. These are the average ratios foundover all the combinations used in the parameter space. For P1planetary orbits the mean critical ratio is materially di ff erent(36%) for prograde and retrograde planetary orbits. Both theseratios are slightly ( ∼ ff erentfor prograde and retrograde planetary motions. For retrogradestellar orbits these ratios are slightly ( ∼ ff er-ence made by the orbital direction of the outer star is small forboth prograde and retrograde planetary orbits, but is statisticallymeaningful for inner retrograde planetary orbits and outer pro-grade planetary orbits.The critical semi-major axis ratios were then regressedagainst the parameters discussed in the previous section usinglinear regression to extract semi-empirical relationships of theform a io a , a oi a = C + b µ + b µ + b e + b e + b i + b Ω + b ω (1) Article number, page 3 of 12 & A proofs: manuscript no. A_A_33097_Busetti
Table 1.
Mean critical semi-major axis ratios for all combinations oforbital motions in P1 and P2 orbits.
Orbit Critical Motions Mean critical ratiotype ratio Star 2 Planet Min Mean σ MaxP1 a io / a P P 0.131 0.383 0.147 0.892R 0.113 0.519 0.108 0.828R P 0.110 0.385 0.077 0.897R 0.112 0.537 0.132 0.857P2 a oi / a P P 1.253 2.936 1.449 5.197R 1.184 1.976 0.507 4.916R P 0.924 2.773 0.891 5.384R 0.591 1.960 0.571 4.957
1. P–prograde, R–retrograde for the outer bounds of S1, S3 and P1 orbits, and the innerbounds of P2 orbits, respectively, where C and b i are the regres-sion constant and coe ffi cients. In all the regressions of P-typeand S-type orbits, the univariate relationships between the crit-ical ratios and the independent variables were linear, with oneexception – the relationship between a io / a and µ in S3 orbits.The results of the regressions are tabulated in Table 2. Inthis and subsequent tables of regression results, the various pa-rameters are defined as follows: σ -standard deviation of criticalratios; R -coe ffi cient of determination, F - F -statistic for overallsignificance, S E -standard error and
MAPE -mean average per-centage error, all being for the overall regression fit; t - t -statisticfor individual coe ffi cients and N -number of data points.The various regressions resulted in 24 regression constantsand 192 coe ffi cients. The semi-major axis ratio a is included asan error check as the critical ratio scales with it and its coe ffi cientshould be zero.Examining the constants, for circular, coplanar outer stellarorbits, the inner and outer stability boundaries for prograde plan-etary orbits are found at around 0.4 times and 2.4 times the dis-tance of the outer star respectively, and for retrograde planetaryorbits, at around 0.6 times and 1.5 times respectively. The pa-rameters i , Ω and ω have a negligible influence on the criticalratio in these low-inclination cases.For P1 orbits, the dominant influence on the critical semi-major axis ratio is only the outer star’s eccentricity, particularlyfor retrograde planetary orbits. For P2 orbits this e ff ect is evenstronger, for both planetary motions. In the case of a retrogradeouter star these influences are largely unchanged for P1 orbits,but for P2 orbits the inner mass ratio becomes important for aprograde planet.The outer star dominates the regions of stability in a triplesystem, with its eccentricity having by far the largest influence.The configuration of the inner binary has little e ff ect on eitherinner or outer orbits. This results from the fact that, as a con-sequence of the Mardling stability limit, the outer star is suf-ficiently far away that the inner binary e ff ectively resembles asingle point mass. These conclusions apply to both prograde andretrograde planetary orbits. However, the greater stability of ret-rograde planetary orbits results in inner and outer bounds thatare closer to the outer star compared with the prograde case.The regression coe ffi cient of the outer eccentricity e is large– for highly eccentric orbits of the outer star, the critical semi-major axis ratio can expand by over 80% for prograde plane-tary orbits and more than double for retrograde orbits. The innerbound can shrink by a quarter for prograde planetary orbits andby over 80% for retrograde orbits. The e ff ect of the direction of motion of the outer star on theplanetary stability bounds is shown in Table 3. The significancelevel is shown for the mean critical ratio, and di ff erences in theregression constant are shown for comparison. Although the ab-solute di ff erences in mean critical ratio are small, two are statis-tically significant.The e ff ect of the direction of motion of planets on their sta-bility bounds is shown in Table 4. The absolute di ff erences are alllarge, averaging 37% for P1 orbits and 31% for P2 orbits, withthe signs being consistent with the greater stability of retrogradeorbits. Since the eccentricity of the outer star is by far the most in-fluential variable on both the inner and outer planetary stabilitybounds, a series of integrations was run to examine the relation-ship between these two variables. Typical relationships for theouter bounds are illustrated in Fig.2, for one stellar configura-tion. a oi /a e Fig. 2.
Outer stability bound for P2 orbits as a function of outer stareccentricity, for a =
100 AU, µ = µ = . , e = i = Ω = ω = i = ◦ and 180 ◦ . For the outer bound the critical ratio should be an increasingfunction of outer eccentricity e . This is true for all four casesshown, although for the two retrograde planet cases the criti-cal ratio flattens out for e (cid:38) .
6. There are discontinuities inthe critical ratio, resulting from resonances, in each case. This ismost visible in the prograde star and prograde planet case, withgaps occurring at eccentricities of 0.10, where the critical ratiojumps from 2.37 to 2.70; and at 0.39, with the ratio undergoinga step change from 2.82 to 3.05. These instabilities are far lesspronounced in the two retrograde planet cases.A retrograde outer body allows the outer bound for a pro-grade planet to move substantially closer to it, with a criticalratio at e = ∼ ∼ ff erence diminishes with increasing outereccentricity, with both critical ratios converging toward 3.2 asthis eccentricity approaches unity. Article number, page 4 of 12. Busetti et al.: Stability of planets in triple star systems T a b l e . S u mm a r yo fr e g r e ss i on s t a ti s ti c s f o r a llt r i p l e o r b it a l c on fi gu r a ti on s , f o rr e g r e ss i on e qu a ti on s a i o / a , a o i / a = C + b µ + b µ + b e + b e + b i + b Ω + b ω . O r b it C r it . H i gh M o ti on s C r it . s e m i - m a j o r a x i s r a ti o R e g r e ss i on c o e ffi c i e n t s I n t e g r a ti on s M od e l fi t t yp e r a ti o i n c l . S t a r P l a n e t M i n M ea n σ M a x C a µ µ e e i Ω ω N o . B ound s R a t e R F S E M A P E f ound ( % )( % ) P a i o / a PP . . . . . . - . - . - . - . . - . - . . . R . . . . . . . - . - . - . . . . . . R P . . . . . . . - . - . - . - . . . . . R . . . . . . - . - . . - . . . . . . NK PP . . . . . . - . - . . . - . . . . . K P . . . . . . - . . . . . . . . . NK R P . . . . . - . . . . . - . . . . . K P . . . . . . - . - . . . - . . . . . P a o i / a PP . . . . . - . . . . . . - . - . . . R . . . . . - . . - . - . . . . . . . R P . . . . . - . . . . . - . . . . . R . . . . . - . . . . . - . . . . . NK PP . . . . . - . . . . . . . . . . K P . . . . . - . . . . . - . . . . . NK R P . . . . . - . . - . . . . . . . . K P . . . . . - . . - . . . . . . . . S a i o / a PP . . . . . . - . . - . - . . . . . . R . . . . . . - . . - . . . . . . . R P . . . . . . - . - . - . - . . . . . . R . . . . . . - . - . - . - . . . . . . S a i o / a PP . . . . . . - . . - . - . . . . . . R . . . . . . . . - . - . . . . . . R P . . . . . . - . . - . - . . . . . . R . . . . . . - . . - . - . . . . . . A v e r a g e / t o t a l . . . . . . . . . . . . . . --- N o t e s . ( ) NK –non - K o za i , K – K o za i ( ) P –p r og r a d e , R – r e t r og r a d e ( ) O f a b s o l u t e v a l u e s f o rr e g r e ss i on c o e ffi c i e n t s ( ) O f i n t e g r a ti on s a ndbound s ( ) T h e d i r ec ti on a l m ov e m e n t s o f t h e s t a b ilit ybound s f o r p l a n e t a r yo r b it s a r e g e n e r a ll y a s e xp ec t e d , e x ce p t f o rr e t r og r a d e i nn e r p l a n e t a r yo r b it s , f o r w h i c h t h e bound i n c r ea s e s i n s t ea do f s h r i nk i ng . T h i s m a yb ea tt r i bu t a b l e t o t h e s m a ll s a m p l e , s i n ce v e r y f e w r e t r og r a d e s t a r / r e t r og r a d e p l a n e t c o m b i n a ti on s w e r e s t a b l e . I t a li c s d e no t ec o e ffi c i e n t s w it h s i gn s t h a t w e r e no t a s e xp ec t e d . I gno r i ng c o e ffi c i e n t ss m a ll e r t h a n0 . , on l y f ou r d i dno t h a v e t h ee xp ec t e d s i gn a nd t h e s e w e r ea ll s m a ll ( < | . | ) . ( ) T h e s e m i - m a j o r a x i s r a ti o a i s i n c l ud e d a s a n e rr o r c h ec k a s t h ec r iti ca l r a ti o s ca l e s w it h it a nd it s c o e ffi c i e n t s hou l db eze r o . I t s a v e r a g ea b s o l u t e v a l u e o f < . t h e r e f o r e p r ov i d e s a n i nd i ca ti ono f t h e i n t r i n s i ce rr o r o f t h e m e t hodo l ogy . Article number, page 5 of 12 & A proofs: manuscript no. A_A_33097_Busetti
Table 3. Di ff erences in mean critical semi-major axis ratios and regres-sion constants for P1 and P2 orbits, for retrograde relative to progradestellar motion, with the same planetary motion. Orbit Ratio Planet Di ff erence in Di ff erence intype motion mean critical regressionratio* constant(%) (%)P1 a io / a P 1 6R 3* 15P2 a oi / a P -6* -12R -1 15 * significant at the 5% level
Table 4. Di ff erences in mean critical semi-major axis ratios and regres-sion constants for P1 and P2 orbits, for retrograde relative to progradeplanetary motion, with the same stellar motion. Orbit Ratio Star 3 Di ff erence in Di ff erence intype motion mean critical regressionratio* constant(%) (%)P1 a io / a P 36* 31R 39* 42P2 a oi / a P -33* -38R -29* -18 * significant at the 5% level
Typical relationships for the inner bounds are illustrated inFig. 3, for the same stellar configuration. For the inner boundthe critical ratio is expected to be a decreasing function of outereccentricity. This is the general trend in each case, albeit withhigh scatter, which reflects that the inner orbit bounds are alwaysmore di ff use than the outer orbit bounds. This is a result of thegreater sparsity of surviving test particles compared with outerorbits, since inner orbits tend to be more chaotic. No stable innerorbits are found for outer eccentricities above 0.7. Prograde star, prograde planetRetrograde star, prograde planetPrograde star, retrograde planetRetrograde star, retrograde planet a io /a e Fig. 3.
Inner stability bound for P1 orbits as a function of outer stareccentricity, for a =
100 AU, µ = µ = . , e = i = Ω = ω = i = ◦ and 180 ◦ . The only P1 orbit found in a triple system is in HW Virginis, witha semi-major axis ratio of 0.37 (Beuermann et al. 2012). No cir-cumtriple P2 orbits have been discovered to date. The planet withthe smallest P1 orbit in a binary, in terms of absolute semi-majoraxis, is Kepler 47b (Orosz et al. 2012), with a semi-major axisratio (here relative to the binary) of a io / a = .
54. The small-est semi-major axis ratio of a planet’s orbit is for Kepler 16b,with a io / a = .
14 (Doyle et al. 2011). Both lie well outsidethe smallest mean critical semi-major axis ratios of ∼ To highlight how the P1 and P2 planetary stability bounds in atriple di ff er from those in a binary, the integrations were repeatedafter reducing the triple to a binary by condensing the inner bi-nary into a single body. The results of this approximation for themean critical semi-major axis ratios for the prograde stellar caseare compared in Table 5. Table 5. Di ff erence between mean critical semi-major axis ratios intriples and binaries, for P1 and P2 orbits. Planetary P1: a io / a P2: a oi / a orbit motion Triple Binary ∆ (%) Triple Binary ∆ (%)Prograde 0.383 0.364 -5* 2.936 2.703 -8*Retrograde 0.519 0.471 -9* 1.976 1.691 -14* ∆ (%) 35 29 - -33 -37 - * significant at the 5% level Both the inner and outer mean critical semi-major axis ratiosmove toward the central star, for both prograde and retrogradeplanetary orbits. While these movements are statistically signifi-cant, they are nevertheless small, confirming that the influence ofthe outer star is dominant. The regression coe ffi cients are com-pared in Table 6.In all cases, aside from the constant the only variable of ma-jor influence is the eccentricity of the outer star, e . Generally,the e ff ect of the outer star’s eccentricity is far larger on the outerstability bounds than on the inner stability bounds. For inner or-bits its influence is larger in binaries than in triples, for both pro-grade and retrograde planetary orbits. For outer orbits, its e ff ectin triples and binaries is essentially the same for retrograde or-bits, but for prograde orbits it has a larger influence in triples thanin binaries. Compared with binaries, the influence of the outerstar’s eccentricity is significantly smaller for inner prograde or-bits and materially larger for outer prograde orbits; there are nosimilar di ff erences for retrograde orbits. Overall, the relativelysmall di ff erences between triples and binaries result from theMardling stability limit for triples, which precludes them frombecoming too compact. Here we examine the stable planetary region for large inclina-tions of the outer star relative to the invariable plane of the triplesystem. Once the orbit of the outer star is no longer coplanarwith the orbit of the inner binary, in other words, the mutual in-clination of the two orbits is no longer zero, the stable planetary
Article number, page 6 of 12. Busetti et al.: Stability of planets in triple star systems
Table 6.
Regression coe ffi cients and model fits for P1 and P2 orbits,triples compared with binaries, for regression equations a io / a , a oi / a = C + b µ + b µ + b e + b e + b i + b Ω + b ω . Parameter Triple Binary ∆ Coe ff .Coe ff . t Coe ff . t (%)Inner orbits, prograde planet C a µ -0.005 -0.4 - - - µ -0.020 -7.1 -0.006 -1.8 -68 e -0.004 -0.5 - - - e -0.114 -8.9 -0.435 -39.4 283 i Ω -0.001 -6.1 - - - ω -0.001 -6.3 - - - R C a -0.003 -23.8 0.000 -25.8 -98 µ µ e e i Ω -0.002 -15.9 - - - ω -0.002 -13.4 - - - R C a µ µ -0.043 -19.2 -0.026 -4.9 -39 e -0.010 -1.4 - - - e -0.468 -31.3 -0.567 -34.9 21 i Ω ω R C a -0.004 -40.4 0.000 -45.7 -99 µ µ -0.008 -2.8 0.023 7.3 -375 e -0.036 -3.5 - - - e i Ω ω R x – y plane is aligned withthe invariable plane. The outer limit of the test particle cloud isof no relevance here, since this is simply where it has been trun-cated. Fig. 4.
Position in space of surviving test particles after 100 kyr. P1 andP2 prograde orbits, with a =
40 AU, e = e = i = ◦ , i = ◦ , µ = . µ = . Examining the paths of individual test particles over timeyields the following observations. Firstly, their orbits librate,with the argument of periapsis ω oscillating. Secondly, thesemi-major axis of their orbits, a , remains e ff ectively constant.Thirdly, their orbits remain circular, with their eccentricity stay-ing very close to zero. Fourthly, their orbits vary in inclination.In an integration they are initially coplanar with the invariableplane and then become more inclined over time, eventually os-cillating around a final non-zero inclination with a period muchlonger than those of the other orbital elements. The typical annu-lar vertical “chimney” shown in Fig 4 therefore consists of testparticles in circular orbits of constant semi-major axis but vary-ing inclination. Fifthly, some test particles develop Kozai reso-nance and some of these eventually undergo orbit flipping intoretrograde motion.There is a limit to the maximum inclination of planetary or-bits. The number of stable orbital bounds falls exponentially withincreasing inclination, as shown in Fig. 5, where the lines endwhen there are no more stable orbital bounds. inclination Both prograde and retrograde stellar orbits were used, and plan-etary motions were prograde. In practical terms, we know thatfor real bodies, as opposed to test particles, the Kozai critical in-clination is larger than the theoretical 39 ◦ , for example Grishinet al. (2017). Also, if the ratio of the initial angular momentumof the outer orbit to that of the inner orbit is (cid:38) ◦ or 270 ◦ . For the integrations the sam-ple was first split into prograde and retrograde stellar cases. Eachcase was then split into subsets where Kozai resonance did or did Article number, page 7 of 12 & A proofs: manuscript no. A_A_33097_Busetti
Number of orbits remainingInitial orbit inclination (deg)
No. μ μ i a Outerstar(-) (-) (deg) (-) orbit1 0.50 0.50 30 40 Prograde2 0.50 0.50 65 40 Prograde3 0.50 0.50 65 20 Prograde4 0.50 0.50 90 40 Prograde5 0.50 2.00 65 40 Prograde6 0.50 2.00 65 40 Retrograde7 0.25 2.00 65 40 ProgradeStellar parameters
Fig. 5.
Number of stable P1 and P2 prograde orbits remaining after1 Myr as a function of initial planetary inclination, for various tripleconfigurations. not occur. Figure 6 shows the critical semi-major axis ratios forthe inner and outer planetary orbits against the “absolute” incli-nation i (i.e., i for prograde orbits and 180 ◦ − i for retrogradeorbits), for both Kozai and non-Kozai regimes. As usual, there ismore scatter for inner orbits than outer orbits.For inner orbits, the critical semi-major axis ratio has almostno dependence on the outer star’s inclination for prograde stel-lar orbits, and only a weak one for retrograde orbits. Retrogradestellar orbits allow the stable planetary region to extend furtherout (with a regression constant of 0.69 versus 0.46), but as in-clination approaches 90 ◦ , this e ff ect decreases. (The regressionlines for inner and outer orbits should ideally coincide at an in-clination of 90 ◦ .)For outer orbits, the critical ratio also has minimal depen-dence on inclination for prograde stellar orbits. However, it canbe seen in the bottom panel of Fig. 6 that this low average depen-dence is made up of two parts. For low inclinations the criticalratio is relatively constant at 1.9–2.0, but as inclinations increasethe critical ratio begins to decline, in other words, the orbital sta-bility bound begins to move inwards. This makes intuitive sensebecause once Kozai resonance begins, as the outer star’s incli-nation increases, the eccentricity of the inner binary decreases,which allows this to occur. This divergence begins at around 45 ◦ and is probably a manifestation of the start of Kozai resonance.This suggests that for real bodies which may also be relativelylarge, Kozai resonance is more likely to begin when i (cid:38) ◦ than for i > ◦ .For retrograde stellar orbits dependence on i is muchstronger. In this case planetary orbits can again approach muchcloser to the outer star (with a regression constant of 1.2 versus1.9), but this increased stability declines as inclination rises to90 ◦ . Also of interest are the regions of stability interspersed withgaps of instability.The mean semi-major axis ratios for the non-Kozai andKozai cases are shown in Table 7. Generally, with a retrogradeouter star the inner orbits move outwards toward this star and theouter orbits move inwards, reflecting the greater stability of ret-rograde orbits. For a prograde outer star the existence of Kozai a io /a = 0.0003 i + 0.456R² = 0.003 a io /a = -0.0021 i + 0.688R² = 0.2070.00.10.20.30.40.50.60.70.80.91.0 0 10 20 30 40 50 60 70 80 90 Prograde star Retrograde star a io /a i a oi /a = -0.0005 i + 1.975R² = 0.006 a oi /a = 0.0081 i + 1.173R² = 0.5711.01.21.41.61.82.02.22.42.62.8 0 10 20 30 40 50 60 70 80 90 Prograde star Retrograde star a oi /a i Fig. 6.
Critical semi-major axis ratios versus outer star inclination: panela) inner stability bound; panel b) outer stability bound.
Table 7.
Non-Kozai and Kozai cases – mean critical semi-major axisratios for prograde and retrograde stellar orbits, for P1 and P2 orbits.
Orbit Regime N Mean ∆ (%) σ Min Maxtype criticalratioPrograde outer starInner NK 256 0.484 0.099 0.049 0.837 a io / a K 47 0.495 2.2 0.100 0.288 0.748Outer NK 341 1.941 0.051 1.755 2.242 a oi / a K 317 1.933 -0.4 0.176 1.759 2.532Retrograde outer starInner NK 177 0.666 0.078 0.423 0.853 a io / a K 101 0.531 -20.4* 0.096 0.105 0.818Outer NK 272 1.335 0.126 1.095 1.820 a oi / a K 204 1.753 31.3* 0.195 1.363 2.510
NK–non-Kozai, K–Kozai* significant at the 5% levelArticle number, page 8 of 12. Busetti et al.: Stability of planets in triple star systems resonance makes no di ff erence to the critical ratios of either theinner or outer bounds, which remain virtually identical. For a ret-rograde outer star, however, the bounds are very di ff erent. WhenKozai resonance occurs the inner bound contracts inwards, whilethe outer bound moves outwards, and by a larger proportionalamount. For both inner and outer orbits these movements are in adirection opposite to that usually induced by a retrograde stellarorbit. This shows, unsurprisingly, that Kozai resonance increasesplanetary instability.At the 5% significance level, the critical semi-major axis ra-tios for planetary orbits under non-Kozai and Kozai regimes donot di ff er for a prograde outer body, but for the retrograde casethey are substantially di ff erent, by an absolute 20%–30%. Thenumber of stable inner orbits found under Kozai resonance witha prograde outer star is small.The regression data for these four cases is shown in Table2. None of the orbital parameters have any significant influenceon the critical semi-major axis ratio, with the possible excep-tion of µ in the prograde star / Kozai case. This stability boundis therefore e ff ectively represented by the constant. The data forthe constant in the four cases is shown in Table 8. Table 8.
Regression constants for non-Kozai and Kozai regimes, for P1and P2 orbits.
Orbit Regime N Constant ∆ Std. t
95% boundstype (%) Err. Lower UpperPrograde outer starInner NK 256 0.545 0.032 17.2 0.483 0.607 a io / a K 47 0.303 -44 0.130 2.3 0.041 0.566Outer NK 341 1.92 0.01 207.6 1.90 1.94 a oi / a K 317 2.35 22 0.05 47.7 2.25 2.44Retrograde outer starInner NK 177 0.685 0.024 28.0 0.637 0.733 a io / a K 101 0.738 8 0.086 8.6 0.567 0.908Outer NK 272 1.29 0.02 67.4 1.25 1.33 a oi / a K 204 1.36 5 0.07 18.3 1.21 1.50
NK–non-Kozai, K–Kozai
For a prograde outer star, Kozai resonance results in the con-stant being 44% smaller for inner orbits and 22% larger for outerorbits. A retrograde outer body does not result in a significantchange in either constant.
The S1 and S2 orbits around stars 1 and 2 respectively are inter-changeable, so only one case is examined, denoted S1. The meancritical semi-major axis ratios for S1 orbits, for the four possiblecombinations of orbital motion, are shown in Table 9.
Table 9.
Mean critical semi-major axis ratios for all combinations oforbital motions in S1 orbits.
Orbit Critical Motions Mean critical ratiotype ratio Star 2 Planet Min Mean σ MaxS1 a io / a P P 0.015 0.180 0.049 0.760R 0.015 0.185 0.122 0.771R P 0.015 0.197 0.059 0.772R 0.028 0.196 0.157 0.772
1. P–prograde, R–retrograde
For S1 planetary orbits in triples, the mean critical semi-major axis ratio is in the range 0.180–0.197. S1 planetary or-bits on average extend ∼
8% further out when the outer star is ina retrograde orbit. For prograde outer stellar orbits, retrogradeplanetary orbits extend slightly (3%) further out while for ret-rograde stellar orbits they are virtually the same. The results ofthe corresponding regressions are shown in Table 2. The key de-terminants of the critical semi-major axis ratio for planetary S1orbits are only the inner binary’s mass ratio and eccentricity, andfor a retrograde outer star the influence of eccentricity is muchweaker.
The two S1 orbits in triple systems found to date, HD 126614A b and HD 2638 b, have semi-major axis ratios of 0.0649 and0.00172 respectively. There is some uncertainty whether a third,Fomalhaut b, is a planet or a dust cloud or disk.Among the closest S1 orbits found in a binary are 0.7 AU forOGLE-2013-BLG-0341L b, with a semi-major axis ratio a io / a of 0.058–0.041 (Gould et al. 2014); 1.09 AU for HD 59686 b,with a io / a ∼ a io / a ∼ The triple was reduced to a binary by setting the outer star’s massto a negligible value and the integrations were repeated. The re-sults of this approximation for the mean critical semi-major axisratios for the prograde stellar case are compared in Table 10.
Table 10. Di ff erence between mean critical semi-major axis ratios intriples and binaries, for S1 orbits. Planetary S1: a io / a orbit motion Triple Binary ∆ (%)Prograde 0.180 0.257 43*Retrograde 0.185 0.289 57* ∆ (%) 2 12 - * significant at the 5% level The outer star of a triple has a significant constraining influ-ence on S1-type orbits compared with the binary case. In triplesthe added influence of the outer star makes S1 orbits move in-wards, with the mean critical semi-major axis ratio reducing byover 30%, from ∼ ∼ ff erence in criticalratio between prograde and retrograde planetary orbits shrinksfrom 12% to 2%. The regressions for the binary case are com-pared with those for the triple case in Table 11.The results for prograde and retrograde planetary orbits intriples and binaries are qualitatively similar, in that the only vari-ables of influence are the inner mass ratio µ and eccentricity e .This e ff ect is not attributable to Kozai, as the mutual inclinationhas no influence. However, the e ff ect of the mass ratio is signifi-cantly larger in triples than in binaries, while that of eccentricity(and the constant) are e ff ectively the same. Article number, page 9 of 12 & A proofs: manuscript no. A_A_33097_Busetti
Table 11.
Regression coe ffi cients and model fits for S1 orbits, triplescompared with binaries, for regression equations a io / a , a oi / a = C + b µ + b µ + b e + b e + b i + b Ω + b ω . Parameter Triple Binary ∆ Coe ff .Coe ff . t Coe ff . t (%)S1 orbits, prograde planet C a µ -0.746 -32.4 -0.497 -50.7 -33 µ e -0.398 -19.7 -0.396 -58.8 0 e -0.006 -0.5 -0.011 -1.6 67 i Ω ω R C a µ -0.667 -29.8 -0.491 -39.6 -26 µ e -0.381 -18.9 -0.411 -47.9 7 e i Ω ω R The mean critical semi-major axis ratios for S3 orbits, for thefour possible combinations of orbital motion, are shown in Table12.
Table 12.
Mean critical semi-major axis ratios for all combinations oforbital motions in S3 orbits.
Orbit Critical Motions Mean critical ratiotype ratio Star 3 Planet Min Mean σ MaxS3 a io / a P P 0.009 0.289 0.098 0.893R 0.010 0.361 0.158 0.920R P 0.007 0.287 0.096 0.853R 0.009 0.360 0.156 0.908
1. P–prograde, R–retrograde
For S3 planetary orbits in triples, the critical semi-major axisratio is in the range 0.287–0.361 (compared with 0.180–0.196for S1 orbits) with wide ranges around these values dependingon the parameters of the triple, and is essentially independentof the direction of motion of the outer star. The critical semi-major axis ratio is around 25% larger for retrograde planetaryorbits than for prograde orbits, significantly larger than the 8%di ff erence found for S1 orbits. The smallest and largest critical semi-major axis ratios a io / a calculated for planets in S3 orbits, based on data in Wagner et al.(2016), are 0.00061 and 0.265 respectively. Recently discovered KELT-4 Ab appears to have a very small a io / a of 0.000132(Eastman et al. 2016), as does Proxima Centauri b’s 3.3 × − (Anglada-Escudé et al. 2016). Although the integrations foundsome stable bounds with semi-major axis ratios as low as 0.009for a few stellar configurations, there were very few critical ra-tios lying below 0.1 and most ranged from 0.1–0.6, with a fewas high as 0.9. This suggests that there may be planets at muchgreater distances from their outer host stars than discovered todate. For S3 orbits the triple was reduced to a binary by merging theinner binary into a single star, as for the P1 / P2 case, after whichthe integrations were repeated. The results of this approximationfor the prograde stellar case are compared in Table 13.
Table 13. Di ff erence between mean critical semi-major axis ratios intriples and binaries for S3 orbits. Planetary S3: a io / a orbit motion Triple Binary ∆ (%)Prograde 0.289 0.309 7*Retrograde 0.361 0.368 2 ∆ (%) 25 19 - * significant at the 5% level The S3 orbits in a triple system have a smaller mean criticalsemi-major axis ratio than in a binary, since the inner binary in atriple has a larger e ff ect on the planet than a single body of equalmass at the same distance. The average di ff erence in S3 orbits forbinaries and triples is around 5%, which is much smaller than theaverage 50% di ff erence found for S1 orbits and more comparablewith the average (absolute) 9% di ff erence for P1 / P2 orbits. Theregressions for the binary case are compared with those for thetriple case in Table 14.Aside from the constant, the only orbital element of impor-tance is e , the eccentricity of the outer star containing the S3orbit. The next largest influence, that of the outer mass ratio µ ,is far smaller. The e ff ect of the mass ratio µ is smaller still, un-like in the S1 case, where it was of equal quantitative importanceto e . These coe ffi cients are virtually identical for prograde andretrograde planetary motions as well as the direction of motionof the outer star. The only empirical work on triples to compare with is that ofVerrier & Evans (2007). Our prograde P1 mean critical ratio of0.383 is 18% smaller than the 0.466 from that study, while ourP2 mean critical ratio of 2.94 is close to their 2.92. The reasonfor the P1 discrepancy is not clear. Although a larger parameterspace and larger clouds of test particles were used in our study,the boundary of P1 orbits is always much more tenuous than thatfor P2 orbits, and the selection of a density cuto ff to define itsedge is necessarily arbitrary. So while results should be consis-tent within a study, a di ff erence of this magnitude across studiesis quite possible for P1 orbits.A larger but still sparse area of overlap exists between our“triple-reduced-to-binary” cases and previous work on binaries.Regression constants from this study are compared with those ofprevious empirical studies, where provided, in Fig. 7.For prograde S1 and P1 orbits our results are very close tothose in previous work. There have been no empirical studies Article number, page 10 of 12. Busetti et al.: Stability of planets in triple star systems
Table 14.
Regression coe ffi cients and model fits for S3 orbits, triplescompared with binaries, for regression equations a io / a , a oi / a = C + b µ + b µ + b e + b e + b i + b Ω + b ω . Parameter Triple Binary ∆ Coe ff .Coe ff . t Coe ff . t (%)S3 orbits, prograde planet C a µ -0.075 -5.2 - - - µ e -0.020 -1.9 - - - e -0.604 -58.6 -0.585 -83.1 -3 i Ω ω R C a µ µ e -0.011 -1.3 - - - e -0.580 -64.6 -0.599 -96.7 2 i Ω ω R Regression constant
P - prograde, R - retrograde
Fig. 7.
Comparison of results for binaries – regression constants. on S1 retrograde orbits. For P1 retrograde orbits there are twoother results; ours coincides with the lower one, from Doolin& Blundell (2011). The critical semi-major axis ratios found inprevious research are shown in Fig. 8, together with our ratiosand regression constants for comparison.For binary S1 prograde planetary orbits, previous results forthe critical semi-major axis ratio varied widely, ranging from0.22–0.46 and averaging 0.38. Our corresponding mean criti-
Critical semi-major axis ratio
P - prograde, R - retrograde
Fig. 8.
Comparison of results for binaries – critical semi-major axisratios and regression constants. cal ratio of 0.26 and regression constant of 0.36 are within thisrange. For P1 prograde orbits the range is 1.60–3.85, the studyby Holman & Wiegert (1999) being the lowest, and with an av-erage of 2.44. Our mean critical ratio and constant of 2.70 and2.71 respectively are 11% higher than this.Two empirical studies have addressed retrograde planetaryorbits in binaries. The one for S1 orbits finds a critical ratio of0.60, with our results for the mean critical ratio and constant, of0.29 and 0.38 respectively, being around 45% lower. The studyfor P1 orbits, by Doolin & Blundell (2011), finds critical ratios inthe range 1.3–2.7, comparable with our results that have a meancritical ratio of 1.69 and a regression constant of 1.25. Despitethe small sample sizes, the above results generally compare well.The one relatively large di ff erence, for retrograde S1 orbits, sug-gests further investigation.
4. Conclusions
We have investigated numerically the long-term stability of plan-ets in a large number of possible orbits around the stars in a triplesystem and provided a generalized mapping of the regions of sta-bility. This was done for prograde and retrograde motion of theplanets and for the outer body of the triple, as well as for highly-inclined orbits of the outer body. The construction of multipleregression equations resulted in 24 semi-empirical models, onefor each type of orbital configuration.The greatest influences on a planet’s critical semi-major axisratio are: for P1 and P2 orbits, the eccentricity of the outer bodyand, to a far lesser extent, the inner mass ratio; for S1 orbits,the inner binary’s eccentricity and inner mass ratio; and for S3orbits, the outer body’s eccentricity and outer mass ratio. Furtherdetail can be found in Busetti (2018).We extended the number of parameters used to all relevantorbital elements of the triple’s stars and their mass ratios. Wealso expanded these parameters to wider ranges that will accom-modate recent and possibly future observational discoveries. Theinvestigation of how the regions of planetary stability in triplesdi ff ered from those in binaries enabled comparison with, and Article number, page 11 of 12 & A proofs: manuscript no. A_A_33097_Busetti confirmation of, some previous studies. It also contributed somenew data points.These generalized results can be useful in the investigation ofobserved systems, providing a fast method of determining theirstability bounds within the large parameter space that resultsfrom observational uncertainties. The relationships expressed inthe regression models can also be used to guide searches forplanets in triple systems and to select candidates for surveysof triple systems. The geometry of the stable zone indicates notonly where to look for planets but also the most appropriate tech-nique to search for them.
References
Anglada-Escudé, G., Amado, P. J., Barnes, J., et al. 2016, Nature, 536, 437Beuermann, K., Dreizler, S., Hessman, F., & Deller, J. 2012, A&A, 543Beust, H. 2003, A&A, 400, 1129Beust, H., Bonfils, X., Montagnier, G., Delfosse, X., & Forveille, T. 2012, A&A,545Borkovits, T., Derekas, A., Kiss, L. L., et al. 2013, MNRAS, 428, 1656Borkovits, T., Rappaport, S., Hajdu, T., & Sztakovics, J. 2015, MNRAS, 448,946Busetti, F. 2018, PhD thesis, University of the WitwatersrandCorreia, A. C. M., Boué, G., & Laskar, J. 2016, CeMDA, 126, 189Derekas, A., Kiss, L. L., Borkovits, T., et al. 2011, Science, 332, 216Doolin, S. & Blundell, K. M. 2011, MNRAS, 418, 2656Doyle, L. R., Carter, J. A., Fabrycky, D. C., et al. 2011, Sci, 333, 1602Eastman, J. D., Beatty, T. G., Siverd, R. J., et al. 2016, AJ, 151Eggleton, P. P. & Tokovinin, A. A. 2008, MNRAS, 389, 869Farago, F. & Laskar, J. 2010, MNRAS, 401, 1189Ford, E. B., Kozinsky, B., & Rasio, F. A. 2000, ApJ, 535, 385Giuppone, C. A. & Correia, A. C. M. 2017, A&A, 605, A124Gould, A., Udalski, A., Shin, I. G., et al. 2014, Sci, 345, 46Grishin, E., Perets, H. B., Zenati, Y., & Michaely, E. 2017, MNRAS, 466, 276Hamers, A. S., Perets, H. B., & Portegies Zwart, S. F. 2015, MNRAS, 455, 3180Holman, M. J. & Wiegert, P. A. 1999, AJ, 117, 621Horner, J., Hinse, T., Wittenmyer, R., et al. 2012, MNRAS, 427, 2812Kozai, Y. 1962, AJ, 67, 591Levison, H. F. & Duncan, M. J. 1994, Icarus, 108, 18Levison, H. F. & Duncan, M. J. 2013, Astrophysics Source Code Library, 1,03001Mardling, R. A. & Aarseth, S. J. 2001, MNRAS, 321, 398Morais, M. H. M. & Giuppone, C. A. 2012, MNRAS, 424, 52Mudryk, L. R. & Wu, Y. 2006, ApJ, 639, 423Musielak, Z. E., Cuntz, M., Marshall, E. A., & Stuit, T. D. 2005, A&A, 434, 355Orosz, J. A., Welsh, W. F., Carter, J. A., et al. 2012, Sci, 337, 1511Pilat-Lohinger, E., Funk, B., & Dvorak, R. 2003, A&A, 400, 1085Santerne, A., Hébrard, G., Deleuil, M., et al. 2014, A&A, 571, A37Sterzik, M. F. & Tokovinin, A. A. 2002, A&A, 384, 1030Tokovinin, A. 2014a, AJ, 147, 86Tokovinin, A. 2014b, AJ, 147, 87Tokovinin, A. 2017, ApJ, 844, 103Trifonov, T., Lee, M. H., Re ff ert, S., & Quirrenbach, A. 2018, AJ, 155, 174Verrier, P. & Evans, N. 2007, MNRAS, 382, 1432Wagner, K., Apai, D., Kasper, M., et al. 2016, Sci, 353, 673ert, S., & Quirrenbach, A. 2018, AJ, 155, 174Verrier, P. & Evans, N. 2007, MNRAS, 382, 1432Wagner, K., Apai, D., Kasper, M., et al. 2016, Sci, 353, 673