Stability Constrained Characterization of Multiplanet Systems
MMNRAS , 1–14 (2020) Preprint 3 December 2020 Compiled using MNRAS L A TEX style file v3.0
Stability Constrained Characterization of Multiplanet Systems
Daniel Tamayo ★ † , Christian Gilbertson , , ‡ † , Daniel Foreman-Mackey Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA Department of Astronomy and Astrophysics, 525 Davey Laboratory, The Pennsylvania State University, University Park, PA 16802, USA Center for Exoplanets & Habitable Worlds, University Park, PA 16802, USA Institute for Computational and Data Sciences, The Pennsylvania State University, University Park, PA, 16802, USA Center for Computational Astrophysics, Flatiron Institute, New York, NY 10010, USA
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Many discovered multiplanet systems are tightly packed. This implies that wide parameterranges in masses and orbital elements can be dynamically unstable and ruled out. We present acase study of Kepler-23, a compact three-planet system where constraints from stability, transittiming variations (TTVs), and transit durations can be directly compared. We find that in thistightly packed system, stability can place upper limits on the masses and orbital eccentricitiesof the bodies that are comparable to or tighter than current state of the art methods. Specifically,stability places 68% upper limits on the orbital eccentricities of 0.09, 0.04, and 0.05 for planets 𝑏 , 𝑐 and 𝑑 , respectively. These constraints correspond to radial velocity signals (cid:46)
20 cm/s,are significantly tighter to those from transit durations, and comparable to those from TTVs.Stability also yields 68% upper limits on the masses of planets 𝑏 , 𝑐 and 𝑑 of 2.2, 16.1, and 5.8 𝑀 ⊕ , respectively, which were competitive with TTV constraints for the inner and outer planets.Performing this stability constrained characterization is computationally expensive with N-body integrations. We show that SPOCK, the Stability of Planetary Orbital ConfigurationsKlassifier, is able to faithfully approximate the N-body results over 4000 times faster. We arguethat such stability constrained characterization of compact systems is a challenging “needle-in-a-haystack" problem (requiring removal of 2500 unstable configurations for every stableone for our adopted priors) and we offer several practical recommendations for such stabilityanalyses. Key words: planets and satellites: dynamical evolution and stability
The continually growing sample of discovered exoplanet systemsprovides a rich laboratory in which to test theories of planet forma-tion. Combining mass measurements from radial velocities withradii inferred from transit data yields mass-radius relationships(Weiss et al. 2013; Weiss & Marcy 2014; Wolfgang et al. 2016;Chen & Kipping 2016; Neil & Rogers 2018; Ning et al. 2018; Neil& Rogers 2020). This has provided first constraints on the interiorstructure of planets beyond our solar system (Valencia et al. 2006;Fortney et al. 2007; Seager et al. 2007; Rogers et al. 2011; Lopez &Fortney 2013; Howe & Burrows 2015; Zeng et al. 2016). Orbital ec-centricities provide complementary information on the dynamicalhistories of planetary systems. Interactions with their natal proto-planetary disks tend to circularize exoplanet orbits (Ward 1988), soeccentricities measured at the present day constrain the degree ofdynamical excitation across a given system’s lifetime (Rasio & Ford ★ NHFP Sagan Fellow: [email protected] † Both authors contributed equally to this manuscript ‡ [email protected] Kepler mission (Borucki et al. 2011). Forthe minority of exoplanets near mean-motion resonances (Fabryckyet al. 2014), it is possible to extract extremely precise masses andeccentricities even for low-mass planets (Agol et al. 2005; Holman& Murray 2005; Nesvorn`y & Morbidelli 2008; Lithwick et al. 2012;Wu & Lithwick 2013; Hadden & Lithwick 2014; Jontof-Hutter et al.2015, 2016; Hadden & Lithwick 2017). Additionally, for eccentricgiant planets (Dawson & Johnson 2012), or when strong densityconstraints on the host star are available (e.g., Van Eylen & Al-brecht 2015), orbital eccentricities can be estimated from transitdurations in individual systems. However, only a small fractionof the sub-Neptunes that dominate the current multiplanet sample © a r X i v : . [ a s t r o - ph . E P ] D ec Tamayo, Gilbertson, Foreman-Mackey have observational constraints on their masses ( ≈ ≈ .In this paper we explore the independent constraints providedfrom orbital stability considerations. Many of the known multi-planet exoplanet systems are in dynamically packed configurations(e.g., Lissauer et al. 2011; Fang & Margot 2012). This implies thatmany combinations of masses and orbital eccentricities would leadto rapid dynamical instabilities in those systems. Given the lowlikelihood of discovery just prior to such a violent orbital rearrange-ment, one can thus reject wide areas of parameter space to constrainphysical and orbital parameters that may otherwise be inaccessi-ble observationally. We refer to this process as stability constrainedcharacterization.Several authors have incorporated such stability constraintsthrough direct N-body integrations for important exoplanet discov-eries (e.g., Steffen et al. 2013; Tamayo et al. 2015, 2017; Quarleset al. 2017; Wang et al. 2018; Rosenthal et al. 2019). However,this brute-force approach is orders of magnitude too computation-ally expensive to effectively sample the high-dimensional parameterspaces for timescales comparable to the Gyr ages of most knownmultiplanet systems. To alleviate the computational burden, authorstherefore typically restrict themselves to checking for “immediate"instabilities within 10 − orbits. However, this is not ideal giventhat dynamical instabilities in typical super-Earth compact plane-tary systems tend to be approximately logarithmically spaced intime, with similar numbers of instabilities occurring in each decadeof time out to at least billions of orbital timescales (e.g., Volk &Gladman 2015).There have been extensive efforts to predict such instabilitiesanalytically, yielding many powerful results for two-planet systems(Wisdom 1980; Marchal & Bozis 1982; Gladman 1993; Barnes &Greenberg 2006; Deck et al. 2013; Petit et al. 2017; Petit et al.2018; Hadden & Lithwick 2018; Hadden 2019). In the 3+ planetcase, several authors have run N-body integrations with initial con-ditions drawn from low-dimensional cuts through the full parameterspace, and fitted empirical functional forms to the resulting insta-bility times (Chambers et al. 1996; Yoshinaga et al. 1999; Marzari& Weidenschilling 2002; Zhou et al. 2007; Smith & Lissauer 2009;Funk et al. 2010; Pu & Wu 2015; Obertas et al. 2017; Gratia &Lissauer 2019). Such empirical fits provide insight into the depen-dencies on physical parameters, and complementary analytic inves-tigations have clarified several aspects of the underlying dynamics(Zhou et al. 2007; Quillen 2011; Laskar & Petit 2017; Yalinewich& Petrovich 2019; Petit et al. 2020). However, none of these modelsare yet able to make reliable stability predictions in general, compact3+ planet systems, particularly ones near mean-motion resonances(Tamayo et al. 2020).Recently, Tamayo et al. (2020) presented the Stability of Plane-tary Orbital Configurations Klassifier (SPOCK), a machine learningmodel capable of making reliable stability predictions over 10 or-bits across a wide variety of compact orbital configurations similarto those discovered by the Kepler and TESS missions, up to 10 times faster than direct N-body integration. They argue that for such (cid:46) Neptune-mass planets, short-timescale ( < orbit) instabilitiesare still dominantly driven specifically by the overlap of two-bodymean motion resonances, as in the two-planet case (Wisdom 1980;Obertas et al. 2017). For a given set of initial conditions, they run a 𝑅 ⊕ ), or a mass (or 𝑀 sin 𝑖 ) < 𝑀 ⊕ ). short (10 orbit) N-body integration (see also Tamayo et al. 2016)to generate a set of ten dynamically motivated summary features,two of which involve the MEGNO chaos indicator (Cincotta et al.2003), and the remaining eight derived from analytical two-planetmean-motion resonance models (Hadden 2019). These features arethen passed to a gradient-boosted decision tree machine learningclassifier (Chen & Guestrin 2016), which returns an estimated prob-ability of stability over 10 orbits. Tamayo et al. (2020) show thattheir model not only performs well on a holdout set of resonanttraining examples (not used during the training process), but also tonon-resonant and higher-multiplicity systems.By enabling long-term stability classification of such compactorbital configurations in a fraction of a second (compared to severalhours for a 10 orbit N-body integration), SPOCK computationallyopens up the stability constrained characterization of multi-planetsystems. In this paper we demonstrate how one can use SPOCKto sharpen poorly constrained physical and orbital parameters ofexoplanets and explore the particular combinations of physical andorbital parameters that stability specifically informs.The paper is organized as follows. We begin in Sec. 2.1 withvarious definitions and a discussion of the implicit assumptions andpossible biases introduced by stability constrained characterization.In Sec. 3 we present stability constraints through direct N-bodyintegrations, and contrast them against those from transit timingvariations and transit durations in the compact 3-planet Kepler-23system, where all methods can be directly compared. In Sec. 4 weconsider faster stability constraints using SPOCK, provide practicalrecommendations for its use, and compare these estimates to theN-body results. We conclude in Sec. 5. Before applying stability constraints to planetary systems, we beginwith some definitions and by laying out the underlying assumptions.The orbital evolution of typical multi-planet systems (in par-ticular of any configurations that would go unstable after a fewdynamical cycles) is chaotic. This leads to effectively stochasticevolution, most notably in the orbital eccentricities, which in simplemodels diffuse until orbits begin crossing and cause close encoun-ters between the planets (e.g., Murray & Holman 1997; Zhou et al.2007).Once orbits cross, it can take a long time for close-in low-mass planets to find one another and physically collide (Rice et al.2018). However, we assume that in such crossing configurations,the large eccentricities and strong scatterings would be detectableeither through transit photometry (e.g., Ford et al. 2008; Kippinget al. 2012; Dawson & Johnson 2012; Van Eylen & Albrecht 2015;Price et al. 2015) or by sharp “chopping” variations in their transittimes due to close approaches (Nesvorn`y & Vokrouhlick`y 2014;Deck & Agol 2015). Following previous authors (e.g., Gladman1993; Zhou et al. 2007; Faber & Quillen 2007; Smith & Lissauer2009; Obertas et al. 2017), we therefore operationally define the“instability time” as the time it takes for a pair of planets to comewithin one Hill sphere of one another. Once this occurs, orbits startcrossing almost immediately (on orbital timescales, Gladman 1993),so this criterion is both accurate and simple to compute numericallywith fast N-body algorithms (e.g., Wisdom & Holman 1991).Instability times are most usefully expressed as a number of or-bits. Because point-source Newtonian gravity is scale invariant, we
MNRAS , 1–14 (2020) tability Constrained Characterization can put different systems on equal footing by expressing all massesrelative to that of the central star, and all times and distances in unitsof the innermost planet’s orbital period and semimajor axis, respec-tively. This facilitates comparisons between systems with differentorbital periods and absolute ages. For example, the ∼
40 Myr age ofthe HR 8799 system (Marois et al. 2008) is only one hundred timesyounger than the typical Gyr ages of observed systems. However,the much longer orbital periods of these directly imaged planets ∼
100 yrs, as compared to ∼ . − . ∼ ∼ fewer orbits than the innermost TRAPPIST-1b, whose orbital periodis only 1.5 days (Gillon et al. 2017).The fact that the required timestep for N-body integrationsscales linearly with the innermost orbital period puts suites of directintegrations over the ∼ -orbit dynamical age of the HR 8799 sys-tem within computation reach (Wang et al. 2018). However, statisti-cal exploration of the long-term stability of configurations of typicalmulti-transiting systems with dynamical ages of ∼ − or-bits becomes computationally prohibitive with N-body methods (afew CPU hours per billion orbits per configuration with the fastestavailable algorithm of Wisdom & Holman 1991). The chaotic dynamics leading to such instabilities have at leasttwo important implications. First, chaos renders Newton’s time-reversible equations of motion effectively irreversible. Planets oncrossing orbits will not scatter back onto circular orbits. In particu-lar, integrating a system’s current orbital configuration backward intime does not reveal its past history beyond a few chaotic (Lyapunov)timescales. Rounding errors due to finite floating point precision(though see Rein & Tamayo 2018) renders integrations forward orbackward in time statistically identical, with eccentricities diffusingupward in both time directions (Gaspard 2005; Morbidelli et al.2020). Given this loss of past information, we focus on integrationsforward in time.Second, this chaos implies that small changes to the initialconditions will yield a range of equally valid instability times. Riceet al. (2018) and Hussain & Tamayo (2019) studied the distribu-tions of such instability times in compact planetary configurations,finding them approximately lognormally distributed. While differ-ent orbital configurations in the dataset analyzed by Hussain &Tamayo (2019) had mean instability times spanning four orders ofmagnitude ( ≈ − orbits), the widths of the lognormal dis-tributions imprinted by the chaotic dynamics was much narrower,approximately 0.4 dex. This can be understood as a consequence ofchaotic random walks in action space (Petit et al. 2020). A givenorbital configuration thus has a well-defined mean instability time,and an N-body integration provides a single draw from a relativelynarrow distribution around that value. These widths quantify theerrors on instability times reported from N-body integrations, setthe fundamental limit on instability time predictions (Tamayo et al.2020), and provide important constraints on the dynamics leadingto instability (Hussain & Tamayo 2019; Petit et al. 2020). Eliminating unstable orbital configurations implicitly makes as-sumptions about the formation and evolution of planetary systems.On one extreme, suppose that the orbital architectures of planetarysystems are effectively set during or shortly after the protoplanetarydisk phase, and from then on remain stable for the age of the uni-verse. Under this assumption, one could straightforwardly rule outorbital configurations with orbital lifetimes shorter than a Hubbletime.By contrast, planetary systems may instead be continuallydestabilizing and rearranging themselves into progressively longer-lived configurations (e.g., Laskar 1990; Volk & Gladman 2015; Pu& Wu 2015; Izidoro et al. 2017, 2019). In the traditional paradigmof core accretion, the final stage of mass growth for sub-Neptunesoccurs in a phase of giant impacts (Goldreich et al. 2004; Hansen &Murray 2012, 2013; Dawson et al. 2016; MacDonald et al. 2020),which could continue occurring over timescales comparable withthe age of the system (e.g., Volk & Gladman 2015; Pu & Wu 2015).More detailed models incorporating the effects of pebble accretionand migration (e.g., Bitsch et al. 2019; Lambrechts et al. 2019) tendto capture planets into chains of mean motion resonances. In order toreconcile this with the observed paucity of systems in such MMRs,such resonant chains must destabilize over time (Izidoro et al. 2017,2019), and again the tail of such collisions could continue to thepresent day. Thus, both of these hypotheses would predict that sys-tems should exist with short remaining lifetimes (i.e., compared tothe system’s age). To our knowledge no such systems have beenidentified to date, though quantitative estimates of how many suchunstable systems one should find are not yet clear. Nevertheless, ifsuch unstable systems might exist, one should clearly not rule outunstable configurations out of hand.We are pursuing separately the quantitative question of howfrequently one should expect to discover systems with a given re-maining lifetime, but we can make a useful simplification. One canalways reliably rule out configurations in the limit of instabilitytimes much shorter than the age of the system. For example, onecould confidently rule out configurations with lifetimes shorter than10 orbits in a system with an age of 10 orbits. One would eitherhave to be extremely lucky to catch a system immediately prior tosuch a cataclysm, or such rearrangements would have to be con-stantly reoccurring. Hadden & Lithwick (2017) uniformly analyzed55 multi-planet systems where TTVs allow characterization of thefull orbital architecture, and find that the majority of their solutionsfor each system are stable over at least 10 orbits . The fact that nosystems with short instability times are currently known stronglydisfavors a scenario where planetary systems are continually rear-ranging on short timescales at the present day, and this qualitativepicture agrees with N-body simulations of giant impact accretion(e.g., Dawson et al. 2016).In this work, we choose to uniformly discard orbital configu-rations with lifetimes shorter than 10 orbits, which corresponds tofractions of a few times 10 − to a few times 10 − in typical sys-tems with lifetimes of a few Gyr and innermost orbital periods of0 . − . The exception being resonant chains (see also Mills et al. 2016; Gillonet al. 2017), where a full exploration of the phase space with Markov ChainMonte Carlo becomes difficult (Tamayo et al. 2017)MNRAS , 1–14 (2020)
Tamayo, Gilbertson, Foreman-Mackey not be applicable to the young HR 8799 planets, with a dynamicallifetime of only 10 orbits (Wang et al. 2018).For simplicity, for the remainder of the paper we therefore referto systems surviving for 10 orbits as (long-term) stable , and onesthat suffer close encounters in that timeframe as unstable . Most current efforts to estimate orbital parameters and masses fromobservational data employ Markov Chain Monte Carlo (MCMC)methods that are efficient in the higher dimensional parameterspaces spanned by multiplanet systems, and make it simple to eval-uate confidence intervals and correlations (e.g., Foreman-Mackeyet al. 2013). One approach to stability constrained characterizationwould be to compute the stability probability at each step in anMCMC analysis. Here, we instead apply the stability constraint as apost-processing step performed after generating a set of samples thatdo not incorporate stability. This has the benefit that this procedureis not tightly coupled to a specific analysis pipeline and can be usedwith existing tool chains and workflows that exist for characterizingexoplanet systems.The goal of any Bayesian characterization is to compute pos-terior weighted expectation integrals such as 𝐸 𝑝 ( 𝜃 | data ) [ 𝑓 ( 𝜃 )] = ∫ 𝑓 ( 𝜃 ) 𝑝 ( 𝜃 | data ) d 𝜃 (1)where 𝜃 represents the set of planet parameters (mass, orbital pe-riod, eccentricity, etc.) and 𝑓 ( 𝜃 ) is the target of our inference (seeHogg & Foreman-Mackey 2018, for example, for a more completediscussion). Typically, the “data” in Equation 1 refers to an obser-vational dataset such as a light curve or radial velocity curve. But,here we aim to include an extra piece of data: we have observedthis planetary system, suggesting that it is in a stable configuration.We refer to the collected observational data as 𝑋 and the “obser-vation” of stability as 𝑞 . With this notation, the relevant posteriorprobability is 𝑝 ( 𝜃 | 𝑋, 𝑞 ) = 𝑝 ( 𝜃 ) 𝑝 ( 𝑋 | 𝜃 ) 𝑝 ( 𝑞 | 𝜃 ) 𝑝 ( 𝑋 ) 𝑝 ( 𝑞 ) (2)where, on the right hand side, we have made the reasonable assump-tion that the observed data and stability are independent conditionedon the physical parameters of the system. In other words, the cal-culation of a light curve model (for example) depends only on theparameters of the system and not on whether or not those parametersare stable.Substituting Equation 2 into Equation 1, we find 𝐸 𝑝 ( 𝜃 | 𝑋, 𝑞 ) [ 𝑓 ( 𝜃 )] = ∫ (cid:20) 𝑓 ( 𝜃 ) 𝑝 ( 𝑞 | 𝜃 ) 𝑝 ( 𝑞 ) (cid:21) 𝑝 ( 𝜃 ) 𝑝 ( 𝑋 | 𝜃 ) 𝑝 ( 𝑋 ) d 𝜃 = ∫ (cid:20) 𝑓 ( 𝜃 ) 𝑝 ( 𝑞 | 𝜃 ) 𝑝 ( 𝑞 ) (cid:21) 𝑝 ( 𝜃 | 𝑋 ) d 𝜃 . (3)Now, if we have somehow (using MCMC or otherwise) generatedsamples 𝜃 ( 𝑛 ) ∼ 𝑝 ( 𝜃 | 𝑋 ) from the posterior probability densityonly conditioned on the observational data, we can approximate theintegral in Equation 3 using the usual sampling approximation 𝐸 𝑝 ( 𝜃 | 𝑋, 𝑞 ) [ 𝑓 ( 𝜃 )] ≈ (cid:205) 𝜃 ( 𝑛 ) 𝑤 ( 𝑛 ) 𝑓 ( 𝜃 ( 𝑛 ) ) (cid:205) 𝜃 ( 𝑛 ) 𝑤 ( 𝑛 ) (4)where 𝑤 ( 𝑛 ) = 𝑝 ( 𝑞 | 𝜃 ( 𝑛 ) ) , i.e., the probability the given sample isstable, and the factors of 𝑝 ( 𝑞 ) have canceled. In other words, thestability constraint can be incorporated by re-weighting the samples from an existing MCMC analysis using the stability probabilitycalculated for each sample in the chain.When performing stability constrained characterizationthrough N-body integrations, 𝑝 ( 𝑞 | 𝜃 ) becomes a simple binaryprobability: either stable or unstable with probability one for a givenset of parameters , reducing to a rejection of unstable configurationssince 𝑝 ( 𝑞 | 𝜃 ) = (cid:26) 𝜃 is a stable configuration0 otherwise . (5)By contrast, we advocate for using the continuous probability ofstability 𝑝 ( 𝑞 | 𝜃 ) estimated by SPOCK (Sec. 4.1).Equation 4 becomes exact in the limit of infinite samples. How-ever, rejecting unstable configurations has the potential to substan-tially decrease the effective number of samples, and hence increasethe error introduced by the sampling approximation. This couldhappen, for example, if the stability constraint is significantly moreconstraining than the observational data. Therefore, when using thismethod, care should be taken to ensure that the effective number ofsamples is sufficient to support the results. In this paper, the initialsamples 𝜃 ( 𝑛 ) are independent by design so it is sufficient to justmake sure that the rejection step leaves a sufficiently large numberof samples, but a more thorough analysis would be required whenaccounting for autocorrelation in MCMC analyses. We now consider the constraints imposed by requiring orbital stabil-ity, and compare the improvements to complementary observationalmethods. We focus on constraints from transit photometry that areimmediately available in large transit surveys like Kepler and TESS,and later comment on radial velocity observations. In particular,TTVs can provide information on the planetary masses (Agol et al.2005; Holman & Murray 2005; Nesvorn`y & Morbidelli 2008; Lith-wick et al. 2012; Wu & Lithwick 2013; Hadden & Lithwick 2014;Jontof-Hutter et al. 2015, 2016; Hadden & Lithwick 2017), andboth TTVs and transit duration modeling can constrain orbital ec-centricities (Dawson & Johnson 2012; Kipping 2014; Van Eylen &Albrecht 2015; Xie et al. 2016).We choose the compact, near-resonant, and asteroseismica-clly characterized three-planet Kepler-23 system as a particularlyillustrative example where stability, transit duration (Van Eylen &Albrecht 2015), and TTV (Hadden & Lithwick 2017) analyses canall be applied and compared.Kepler-23 (KOI 168) is a 1 . ± . 𝑀 (cid:12) , 13.4 Kepler-magnitude star (Huber et al. 2013), with an estimated age of 4-8Gyr (Ford et al. 2012). Its three known planets are tightly packed,with orbital periods of 7.107, 10.742 and 15.274 days (with periodsratios between adjacent planets of 1.511 and 1.422, respectively),and radii of 1 . ± . , . ± .
10 and 2 . ± .
09 Earth-radii( 𝑅 ⊕ , respectively (Van Eylen & Albrecht 2015). Its dynamical agethus corresponds to 2 − × inner-planet orbits.Precise stellar parameters from asteroseismology allow foreccentricity constraints from transit durations. A shorter-than-expected transit duration can be explained by its occurrence nearthe pericenter of an eccentric orbit. However, there is a degeneracy In reality, the uncertainty arising from chaotic dynamics (Sec. 2.2) impliesthat 𝑝 ( 𝑞 | 𝜃 ) is not a step function but has finite width. Nevertheless, becausethis width is small, as measured numerically by Hussain & Tamayo (2019),a step function is a good approximation. MNRAS000
09 Earth-radii( 𝑅 ⊕ , respectively (Van Eylen & Albrecht 2015). Its dynamical agethus corresponds to 2 − × inner-planet orbits.Precise stellar parameters from asteroseismology allow foreccentricity constraints from transit durations. A shorter-than-expected transit duration can be explained by its occurrence nearthe pericenter of an eccentric orbit. However, there is a degeneracy In reality, the uncertainty arising from chaotic dynamics (Sec. 2.2) impliesthat 𝑝 ( 𝑞 | 𝜃 ) is not a step function but has finite width. Nevertheless, becausethis width is small, as measured numerically by Hussain & Tamayo (2019),a step function is a good approximation. MNRAS000 , 1–14 (2020) tability Constrained Characterization with the orientation of the pericenter relative to the line of sight(one can observe the same transit duration with a more eccentricorbit by adjusting the orbit orientation such that the transit occursfurther from pericenter, see e.g., Fig. 1 of Van Eylen & Albrecht2015). Additionally, there is a degeneracy with the transit impactparameter, which will also change the transit duration, though thiscorrelation can sometimes be alleviated by modeling the distortionof the transit shape with increasing impact parameter. These effectsgenerically result in eccentricity constraints from transit durationshaving a tail toward large values (Van Eylen & Albrecht 2015).The period ratios between adjacent planets (1.511 and 1.422,respectively) fall near strong MMRs (3:2 and 7:5, respectively). Thisinduces TTVs that have been measured and modeled to put con-straints on the masses and orbital eccentricities (Ford et al. 2012;Hadden & Lithwick 2017). For these TTV and the above tran-sit duration comparisons, we take the publicly available posteriordistributions from Hadden & Lithwick (2017) and Van Eylen &Albrecht (2015), respectively.Finally, we consider the constraints imposed by long-term sta-bility. As argued in Sec. 2.4, stability constraints are complementaryand easily combined with observational constraints. However, fora straightforward comparison, we choose to present stability con-straints independently. Thus, rather than rejecting unstable posteriorsamples 𝜃 ( 𝑛 ) ∼ 𝑝 ( 𝜃 | 𝑋 ) as constrained by TTVs or transit durations(Sec 2.4), we reject unstable configurations starting from the sameprior distribution of orbital configurations as considered by Hadden& Lithwick (2017) and Van Eylen & Albrecht (2015) in their TTVand transit duration analyses.The most reliable (but most computationally expensive) way ofestimating stability is through direct N-body integration. We beginby comparing these “ground truth” results to the other methods, andin Sec. 4 compare these N-body stability constraints to faster resultsobtained with SPOCK.It is important to point out that constraints from each of themethods we compare depend on a variety of physical and observa-tional factors. TTV constraints, in particular, can in some cases bemuch more precise than stability. For reference, out of the 90 planetpairs uniformly analyzed by Hadden & Lithwick (2017), Kepler-23is roughly a typical case, with about 60% of planet pairs yieldingbetter eccentricity constraints (taken as the ratio of the reportederror bars to the peak of the posterior probabilities).Transit durations are less informative for Kepler-23 than av-erage among the sample selected by Van Eylen & Albrecht (2015)with asteroseismology. Approximately 70% of their sample of ≈ To facilitate comparisons, we follow Hadden & Lithwick (2017)and Van Eylen & Albrecht (2015) by drawing all planet parametersindependently, sampling eccentricities uniformly and masses log-uniformly. Eccentricities were drawn from [0,0.9], and masses bysampling bulk densities from 0.3-30 g/cc and using fixed planetary radii of 1.8, 3.2 and 2.3 𝑅 ⊕ for planets, 𝑏 , 𝑐 and 𝑑 , respectively.As typically the case for transiting planets, the orbital periods aremeasured precisely enough that for simplicity we fix them to thevalues listed above. Particularly close to MMRs, these period ratioscan matter for stability, so sampling them from a wider range wouldbe prudent when periods are less certain, for example if there aretransit timing variations over a small set of transits.We sampled the longitudes of ascending node, arguments ofpericenter and mean anomalies uniformly from [ , 𝜋 ] , and incli-nations uniformly from zero to values corresponding to an impactparameter of 0.9 stellar radii. More carefully accounting for cor-relations between orbital eccentricities and impact parameters isimportant for some applications, e.g., modeling transit durations(Van Eylen & Albrecht 2015). By contrast, stability constraintsare very weakly sensitive to small mutual inclinations. Experi-mentation by Tamayo et al. (2020) with inclination-dependent fea-tures provided negligible improvements to their machine learningmodels’ performance on compact systems with mutual inclinations (cid:46) ◦ , and their final SPOCK model has no features that dependon inclinations explicitly. Choices related to sampling inclinationsshould therefore be guided by the original orbital parameter infer-ence problem, and should be negligible for stability constraints innear-coplanar systems over timescales of 10 orbits.Finally, we drew normally distributed stellar masses of 𝑀 ★ = ± . 𝑀 (cid:12) . We note that as mentioned above, gravity’s scale invari-ance implies that the dynamics only depend on the planet-star massratios, rather than their individual masses. An alternative, therefore,that decouples uncertainties in stellar parameters would be to con-strain planet-star mass ratios rather than absolute planetary masses.The star’s absolute mass does change the orbital periods ∝ 𝑀★ / ,so one can convert from 10 orbits to absolute time and from planet-star mass ratios to planetary masses by folding in the uncertaintiesin stellar mass.We ran each configuration for 10 orbital periods of the inner-most planet using a fixed timestep of 0.25 days ≈ WHFast (Rein & Tamayo2015), which is part of the open-source package
REBOUND (Rein &Liu 2012), and is based on Wisdom & Holman (1991). We reject anyconfigurations where any two planets come within the sum of theirindividual Hill radii (see Sec. 2.1) within the span of integration,and accept samples that survive the full 10 orbits. In Fig. 1 we compare the constraints on the planet masses, with68% highest posterior density intervals listed in the legend (seefigure caption). Transit durations hold no information on planetmasses, and we can see that TTVs provide stronger constraints thanstability. Additionally we point out that while TTVs can providelower bounds on planet masses (as for planet 𝑏 in Fig. 1), requiringstability typically only provides upper limits.We see that in this case, stability provides comparable uppermass limits to TTVs for the innermost and outermost planets (seelegend of Fig. 1), but not for the middle planet. We can understandthis qualitatively.At comparable planet masses, there will always be a preferencetoward lower-mass planets, which allow for stability over a broaderrange of orbital eccentricities. However, in the limit where a singleplanet dominates the mass and the other bodies can be treated astest particles, the problem reduces to sets of analytically understoodtwo-planet stability problems, which are much less restrictive. As anillustration, a pair of planets on coplanar and initially circular orbits MNRAS , 1–14 (2020)
Tamayo, Gilbertson, Foreman-Mackey Planet b Mass ( M ) P D F PriorStability[0, 2.23] M Transit Time Variations[0.75, 2.78] M Planet c Mass ( M ) P D F [0, 16.10] M [1.55, 5.34] M Planet d Mass ( M ) P D F [0.01, 5.79] M [0.37, 4.21] M Figure 1.
Comparison of constraints on the masses of the three planets in the Kepler-23 systems from transit time variations and stability. Priors are plotted ingray. Following Hadden & Lithwick (2017), we also list in the legend the 68.3% highest posterior density intervals, i.e., the smallest parameter range containing68% of the distribution. will never undergo close encounters as long as their separationis greater than ≈ . 𝑚 / ,this corresponds to mass constraints that differ by ∼
30 betweentwo and three-planet cases. In the Kepler-23 system, the middleplanet has the largest radius, causing a preference for masses > 𝑀 ⊕ , which are excluded by the density priors on the otherplanets. Stability constraints are therefore strongest in systems withcomparable mass planets, which appears to be a common outcomein observed systems from investigations of correlations betweenplanetary radii and masses in multiplanet systems (Millholland et al.2017; Weiss et al. 2018; Gilbert & Fabrycky 2020; He et al. 2020).Taking this multi-transiting system discovery, we can also askwhat radial velocity precision would be required to reach the upper-mass limits imposed by stability. The 68% upper mass limits of 2.2,16.1, and 5.8 𝑀 ⊕ correspond to radial velocity semi-amplitudes of0.8, 4.8 and 1.5 m/s for planets 𝑏 , 𝑐 and 𝑑 , respectively. We now compare the constraints on the orbital eccentricities forthe Kepler-23 planets from transit durations (Van Eylen & Albrecht2015), TTVs (Hadden & Lithwick 2017) and stability in the top rowof Fig. 2, again listing the 68% highest posterior density intervals.In this case, stability constraints approach those from TTVs, andare much stronger than those from transit durations.Orbital eccentricities introduce radial velocity signals that aresmaller than the observed semi-amplitude by a factor of ∼ 𝑒 (see,e.g., Shen & Turner 2008; Zakamska et al. 2011). The 68% eccen-tricity upper limits from stability of 0 . − . (cid:46)
20 cm/s for the massesfound in Sec. 3.2.
One can show that for a single pair of planets near a j:j-k MMR ,there is a particular combination of the eccentricity vectors e 𝑖 (withmagnitude given by 𝑒 𝑖 and direction given by the longitude ofpericenter, 𝜛 𝑖 ) that is approximately conserved, while a secondcombination Z drives the resonant dynamics, Z 𝑖,𝑖 + ≈ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) e 𝑖 + − e 𝑖 + √ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (6)where the exact expression carries additional coefficients for theeccentricity vectors that depend on the 𝑗 and 𝑘 indices of the reso-nance. Except for the 2:1 MMR, these coefficients are within ≈ Z can further bedecomposed into a component that is forced by the MMR, and afree component set by initial conditions. However, because most ob-served TTV systems are typically far from exact resonance ( ≈ − Z as free eccentricities.When sinusoidal variations in transit times measure eccen-tricities, they dominantly constrain these free eccentricities (e.g.,Lithwick et al. 2012). TTVs thus provide narrower bounds on Z than on the individual eccentricities .At the same time, MMRs are also understood to drive fastinstabilities in compact multiplanet systems (Wisdom 1980; Quillen2011; Deck et al. 2013; Petit et al. 2017; Hadden & Lithwick 2018;Petit et al. 2020). Indeed, Tamayo et al. (2020) showed that their Each j:j-k resonance actually includes k+1 sub-resonances, but these canbe combined into a single resonance via canonical transformation. This wasoriginally shown (exactly) for first-order (k=1) resonances (Sessin & Ferraz-Mello 1984; Wisdom 1986; Henrard et al. 1986), but Hadden (2019) showsthat this also approximately holds for higher order resonances. In principle this is not a one-to-one comparison. If one orbit were eccentricand the other circular, Z would be a factor of √ 𝑒 , then onaverage Z = 𝑒 . MNRAS000
One can show that for a single pair of planets near a j:j-k MMR ,there is a particular combination of the eccentricity vectors e 𝑖 (withmagnitude given by 𝑒 𝑖 and direction given by the longitude ofpericenter, 𝜛 𝑖 ) that is approximately conserved, while a secondcombination Z drives the resonant dynamics, Z 𝑖,𝑖 + ≈ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) e 𝑖 + − e 𝑖 + √ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) , (6)where the exact expression carries additional coefficients for theeccentricity vectors that depend on the 𝑗 and 𝑘 indices of the reso-nance. Except for the 2:1 MMR, these coefficients are within ≈ Z can further bedecomposed into a component that is forced by the MMR, and afree component set by initial conditions. However, because most ob-served TTV systems are typically far from exact resonance ( ≈ − Z as free eccentricities.When sinusoidal variations in transit times measure eccen-tricities, they dominantly constrain these free eccentricities (e.g.,Lithwick et al. 2012). TTVs thus provide narrower bounds on Z than on the individual eccentricities .At the same time, MMRs are also understood to drive fastinstabilities in compact multiplanet systems (Wisdom 1980; Quillen2011; Deck et al. 2013; Petit et al. 2017; Hadden & Lithwick 2018;Petit et al. 2020). Indeed, Tamayo et al. (2020) showed that their Each j:j-k resonance actually includes k+1 sub-resonances, but these canbe combined into a single resonance via canonical transformation. This wasoriginally shown (exactly) for first-order (k=1) resonances (Sessin & Ferraz-Mello 1984; Wisdom 1986; Henrard et al. 1986), but Hadden (2019) showsthat this also approximately holds for higher order resonances. In principle this is not a one-to-one comparison. If one orbit were eccentricand the other circular, Z would be a factor of √ 𝑒 , then onaverage Z = 𝑒 . MNRAS000 , 1–14 (2020) tability Constrained Characterization Planet b Eccentricity ( e b ) P D F PriorTransit Durations[0, 0.31]Stability[0, 0.09]Transit Time Variations[0, 0.08]
Planet c Eccentricity ( e c ) [0, 0.42][0, 0.04][0, 0.05] Planet d Eccentricity ( e d ) [0, 0.31][0, 0.05][0, 0.05] Inner-Pair Free Ecc. ( Z bc ) P D F [0.04, 0.40][0.01, 0.07][0.01, 0.04] Outer-Pair Free Ecc. ( Z cd ) [0.05, 0.38][0.01, 0.04][0.01, 0.04] Center-of-mass Ecc. ( e com ) [0.03, 0.29][0, 0.04][0, 0.05] Figure 2.
Top row: Comparison of constraints on the orbital eccentricities of the three planets in the Kepler-23 systems from transit durations, TTVs, andstability. Priors are plotted in gray. Following Van Eylen & Albrecht (2015) and Hadden & Lithwick (2017), we also list in the legend the 68.3% highest posteriordensity intervals, i.e., the smallest parameter range containing 68% of the distribution. Bottom row: Constraints on the particular eccentricity combinations Z (Eq. 6) that drive the resonant dynamics between the inner (left panel) and outer (middle panel) planet pair. The right panel plots a final combination ofeccentricities that generalizes a conserved quantity in a two-planet, single MMR model (Eq. 7). machine learning model SPOCK, trained only in and near thesediscrete MMRs, generalizes to uniformly spaced systems. Similarto TTVs, one would therefore expect stability to better constrain freeeccentricities between pairs of adjacent planets than their respectiveindividual eccentricities. This has been demonstrated analyticallyfor two-planet systems. For higher multiplicity systems it is difficultto make precise statements without an analytic understanding yetin hand. However, one would qualitatively expect that the closera pair of planets are to a strong MMR, the more stability shouldspecifically constrain the combination Z , rather than the individualeccentricities.Transforming from 𝑁 individual eccentricities to 𝑁 − 𝑒 com = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:205) 𝑁𝑖 = 𝑚 𝑖 e 𝑖 (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:205) 𝑁𝑖 = 𝑚 𝑖 , (7)where the sum runs over all the planets. While this quantity is notstrictly conserved, it reduces to an approximately conserved quantity in both the limit where a pair of near-resonant planets dominate themass, and in the limit of a single dominant planet.In the bottom row of Fig. 2, we show constraints on thesetransformed quantities. We see that stability, and especially TTV,constraints on free eccentricities are tighter. In particular, TTVs areable to exclude zero from the 68% interval. We note that while thestability derived distributions now also peak at non-zero values, thisis an artifact imposed by the adopted prior. Choosing eccentricityvectors randomly and independently results in few samples withcomparable and aligned eccentricity vectors that would yield smallvalues of Z (prior distributions in gray in Fig. 2). This might ap-proximate the random outcomes of giant impacts or planet-planetscattering, but would be a poor prior for, e.g., smooth migrationin a disk with eccentricity-damping, which would act to damp freeeccentricities to zero. In this case a uniform prior in the free eccen-tricities would be more appropriate.Finally, we point out that while TTV modeling can rule out Z = , stability will typically only yield upper limits, since things aregenerally more stable at lower eccentricities. The exception wouldbe near resonances that force non-zero equilibrium eccentricities. In Both the amplitude and phase of the TTV signal depend on Z (e.g.,Lithwick et al. 2012).MNRAS , 1–14 (2020) Tamayo, Gilbertson, Foreman-Mackey
Inner-Pair Free Ecc. ( Z bc ) P D F N-body Orbits[0.03, 0.12]10 N-body Orbits[0.01, 0.07]SPOCK[0, 0.08]
Outer-Pair Free Ecc. ( Z cd ) [0.01, 0.08][0.01, 0.04][0.01, 0.05] Center-of-mass Ecc. ( e com ) [0, 0.10][0, 0.04][0, 0.05] Figure 3.
Comparison of the eccentricity distributions of configurations that remain stable over 10 orbits in N-body integrations (orange) to faster estimateswith SPOCK (purple). Gray histograms show configurations that remain stable over 10 orbits. SPOCK (Tamayo et al. 2020) effectively removes 96% of theunstable configurations in the gray histogram to yield the purple histogram (see end of Sec. 4.5). such cases, circular orbits far from the equilibrium could be stronglyunstable and ruled out. For example, Obertas et al. (2017) foundreductions in instability times by several orders of magnitude nearMMRs in their numerical integrations of initially circular systems. The N-body characterization above required approximately 8500CPU hours in order to identify 833 stable configurations. Giventhat typical MCMC analyses typically aim to have a factor of ∼ − -orbit N-body integrations to rule out only the fastest instabilities, and usingSPOCK, the Stability of Planetary Orbital Configurations Klassifier(Tamayo et al. 2020). SPOCK is a machine learning model thatalso runs a direct 10 -orbit N-body integrations, but is trained tomake stability predictions over 10 orbits from this short time series(Tamayo et al. 2020).As justified below (Sec. 4.3), we begin by rejecting orbit-crossing configurations, the vast majority of which should be un-stable over long timescales. In Fig. 3, we compare the results ofthese faster approaches to the 10 -orbit N-body distributions la-beled ‘Stability’ in Figs 1 and 2. For the SPOCK predictions, weplot histograms of configurations weighted by their estimated prob-abilities of stability from the model following Eq. 4 (see Sec. 4.1).Of the 9414 configurations that survive 10 orbits in gray, over90% go unstable when integrated to 10 orbits. We see in Fig. 3 thatthis leads to eccentricity constraints that are too wide by roughlya factor of two. SPOCK, by contrast, is able to match long-termN-body constraints to within approximately 20% at comparablecomputational cost to the much shorter 10 -orbit integrations. Thisillustrates SPOCK’s potential for the fast characterization of com-pact multiplanet architectures.We now further analyze SPOCK’s performance to understandthe source of deviations with the N-body distributions, and to makerecommendations for its application to new systems. For every input orbital configuration, SPOCK returns an estimatedprobability of stability. For simplicity, Tamayo et al. (2020) consid-ered binary classification into stable and unstable systems, whichone does by choosing a threshold probability to separate the twoclasses (Eq. 5).We advocate instead for approximating 𝑝 ( 𝑞 | 𝜃 ) usingSPOCK’s continuous stability probability estimates directly. Whileone can try to find a threshold that balances the rejection of mostunstable systems without throwing out a significant fraction of thestable configurations one wishes to characterize, the right thresholdwill vary by system, and binary classification always throws out in-formation. By instead weighting all configurations by their stabilityprobabilities estimated by SPOCK (Eq. 4), we allow for more confi-dent classifications from the model to be counted more heavily thanones close to a hand-tuned threshold. This probabilistic approachalso helps avoid potential pitfalls when using SPOCK on sets ofconfigurations distributed very differently from its set of trainingexamples. We illustrate this through an instructive example.Tamayo et al. (2020) present a case study of the three-planetKepler-307 system (third planet is a candidate), where they claimthat SPOCK fails (their Fig. 7). In particular, they drew 1500 sam-ples from the posterior of a TTV analysis of the system and generatedSPOCK probabilities for each configuration. In their earlier analy-sis, Tamayo et al. (2020) held back 20% of their training dataset.Analyzing SPOCK’s performance on this 20% holdout set, theyfound that labeling systems with stability probabilities estimated bySPOCK > .
34 lead to a false positive rate (FPR) of 10%, i.e., only10% of unstable systems were misclassified as stable.However, when they applied this same threshold to Kepler-307, this lead to an FPR of 87%. Tamayo et al. (2020) arguedthat because the TTV analysis had already strongly constrained theresonant dynamics in the system, this caused the 1500 posteriorsamples to strongly cluster in SPOCK’s feature space, making itdifficult for the model to separate stable from unstable systems.While this is certainly a factor, we now instead argue that thisapparent failure is principally a result of their binary classification,and that a probabilistic approach resolves most of the disagreement.For the sake of argument, let us assume that the probabilities
MNRAS000
MNRAS000 , 1–14 (2020) tability Constrained Characterization Inner-Pair Free Ecc. ( Z bc ) P D F [0.006, 0.007][0.006, 0.007] Outer-Pair Free Ecc. ( Z cd ) [0.06, 0.14][0.08, 0.15] Center-of-mass Ecc. ( e com ) N-body Orbits[0.01, 0.06]SPOCK[0.01, 0.07]
Figure 4.
Separate test of SPOCK vs N-body on a different system, Kepler-307, presented by Tamayo et al. (2020) as a case where SPOCK fails. By continuouslyweighting configurations by their probabilities of stability as estimated by SPOCK (purple), we obtain significantly better agreement with stability constrainedeccentricity posteriors performed with N-body integrations (orange) than by performing binary classification into stable and unstable configurations using athreshold optimized on SPOCK’s training dataset as was done by Tamayo et al. (2020) (cf. their Fig. 7). of stability estimated by SPOCK are true probabilities. In SPOCK’straining set, there are many unstable configurations, which SPOCKassigns low stability probabilities. The fact that SPOCK’s FPR onthis distribution of configurations is 10% for a stability threshold of0.34 means that 90% of the unstable systems in its training set areassigned stability probabilities < 0.34.But now imagine giving SPOCK a different set of much morestable orbital configurations, where the SPOCK probabilities areall > .
4. We should still expect many of these configurations tobe unstable (e.g., only half of those with probability 0.5 shouldbe stable). But the fact that none of the systems in this dataset have stability probability below the 0.34 threshold, chosen usingthe training dataset , means that all of the unstable systems aremisclassified as stable, i.e., the FPR is 100%.This is essentially what happened with Kepler-307. The priorTTV analysis removed the vast majority of clearly unstable config-urations, leaving over 95% of samples with stability probabilities >0.34 as estimated by SPOCK. The fact that so few configurations fallbelow the threshold again means that the vast majority of unstableconfigurations are labeled as stable, yielding a misleadingly highFPR of 87%. In summary, a given probability threshold for binaryclassification is only sensible when the data distribution being testedis similar to the training data distribution.The alternative is to weight configurations by their respective,continuous stability probability 𝑝 ( 𝑞 | 𝜃 ) as estimated by SPOCK(Eq. 4). We define the number of effective stable samples returnedby SPOCK as the sum of SPOCK probabilities across the 1500 inputsystems. This yields 1121 effective stable samples from SPOCK,as compared to 967 stable samples as determined through N-bodyintegrations. We therefore see that the probabilities estimated bySPOCK are overall too high, though the discrepancy is at the 10%level.In Fig. 4, we show the equivalent plot to Fig. 3 for Kepler-307,which also shows agreement with N-body to ≈ 𝑐 and 𝑑 than SPOCK. As argued by Tamayo et al. (2020),we suspect that this is due to longer-term secular dynamics (e.g.,Hadden 2019) not effectively captured by SPOCK’s features, which focus on the short-term resonant dynamics. However, the agreementwith N-body is much closer than suggested by the 87% FPR reportedby Tamayo et al. (2020) using a binary classification threshold. Wetherefore recommend weighting configurations by their stabilityprobabilities as estimated by SPOCK.Using continuous probabilities additionally makes the resultdifferentiable with respect to orbital parameters, which allows formore efficient sampling in high-dimensional spaces through Hamil-tonian Monte Carlo techniques (Foreman-Mackey et al. 2020).SPOCK uses non-differentiable gradient-boosted decision trees, butwe are working on differentiable models using neural networks(Cranmer et al., in prep ). The previous section considered an application of SPOCK whereTTV constraints had already ruled out the most unstable configura-tions, leaving behind a mostly stable distribution of samples. Moreinteresting are cases where stability can rule out the majority of con-figurations to significantly constrain orbital parameters and masses,like in Kepler-23. It is instructive to first consider why such a taskmight be challenging for any stability estimator, including SPOCK.In our analysis of Kepler-23, we sampled roughly 2 millionorbital configurations, of which only 837, or ∼ ∼ × ) would lead one to include ∼ ,
000 false positives in the predicted distribution of stable sys-tems. This would swamp the 837 stable systems, yielding unreli-able results. This textbook problem, often framed in terms of animperfect medical test for a rare disease, arises across “needle-in-a-haystack” applications with strong class imbalances. We also notethat the most tightly packed systems, which provide the strongeststability constraints, by definition have the smallest fraction of pos-sible stable configurations. The most interesting systems for stability
MNRAS , 1–14 (2020) Tamayo, Gilbertson, Foreman-Mackey constrained characterization are therefore also the most challengingto accurately characterize.
We found that one helpful way to alleviate this problem is to simplythrow out configurations where any of the orbits cross. Whetheradjacent orbits cross is a non-linear problem involving the orbitalorientations, but for closely spaced planets, the condition for theorbits of planet 𝑖 and 𝑖 + | e i − e i + | > 𝑎 𝑖 + − 𝑎 𝑖 𝑎 𝑖 + , (8)where the 𝑎 denote semimajor axes, and (cid:174) 𝑒 ≡ ( 𝑒 cos 𝜛, 𝑒 sin 𝜛 ) isa vector pointing in the direction of the longitude of pericenter 𝜛 ,with magnitude given by the orbital eccentricity.This is a reasonable cut given that planets on crossing orbitsshould eventually “find" one another and scatter or collide; however,it is not a rigorous step. In particular, orbit-crossing configurationscan be long-lived in MMRs (like those of Neptune and Pluto). Onewould therefore not want to take this step if there are reasons tosuspect such a configuration, e.g., from short transit durations. Wenote that in such cases, sampling eccentricities and pericenter ori-entations independently as done above will be very inefficient insampling the resonant island (see Tamayo et al. 2020, for moreeffective sampling methods).In the typical case, however, where one can reasonably rejectcrossing configurations, this can significantly improve the accuracyof stability-constrained posteriors with SPOCK for two reasons.First, the training dataset for SPOCK was generated drawing or-bital eccentricities log-uniformly up to orbit-crossing values. WhileSPOCK does generally assign low stability probabilities to crossingconfigurations, the lack of examples near and beyond orbit-crossingcan lead to poor extrapolation.More importantly, rejecting crossing orbits significantly cutsdown the number of unstable configurations, alleviating the classimbalance challenge discussed above. In this case, rejecting crossingconfigurations reduces the ∼ , was stableover 10 orbits. A more relevant comparison of these 12 outliers is to thetotal number of 837 stable configurations found. However, the mis-rejection of these 1% of stable configurations is more than compen-sated by the removal of ≈
99% of unstable systems, which we findwould otherwise make up ≈
20% of the final estimated distribution(effective stable samples) in the form of false positives.
SPOCK makes its stability predictions from summary features mea-sured over a short N-body integration of 10 orbits. This allowsSPOCK to catch systems that destabilize quickly, and rigorously In fact, manual examination of the 12 rejections that were long term stablereveal that while all were close, none were actually orbit-crossing—due tothe leading order approximation in Eq. 8. Given the negligible 1% error onthe estimated distributions, we do not pursue this correction. assign them zero probability of stability over 10 orbits . The factthat these highly unstable systems are classified perfectly, and do notcontribute any residual probability to the estimated stable distribu-tion as false positives further alleviates the class imbalance problemdiscussed above.In total, rejecting crossing configurations removed 99% of ourunstable configurations, and the short 10 orbit integrations fromSPOCK further removed approximately 60% as short-lived. Thisleft 9414 non-crossing configurations, surviving at least 10 orbits,of which 825 ≈
10% are stable. While the remaining class imbalanceratio of 10:1 is still a challenge, the improvement from the originalratio of 2500:1 is significant.We note that while in principle the short integrations can alsocatch most of the unstable orbit-crossing configurations, the over-whelming number of crossing configurations causes the small frac-tion of surviving false positives to significantly skew the estimateddistributions. We therefore found that separately rejecting crossingconfigurations significantly improved predictions in this case withwide priors.
We now analyze the ≈
20% discrepancy between SPOCK and N-body shown in Fig. 3. In particular, are errors dominated by residualfalse positives, or are SPOCK probabilities significantly distortingthe distribution of stable systems? The former problem is signifi-cantly preferable to the latter. After all, this is the baseline situationbefore applying any stability constraints, where all unstable con-figurations are included as false positives. Thus, any approximatestability criterion is helpful (and conservative) as long as it removesunstable systems while preserving the distribution of stable config-urations.In Fig. 5, we plot histograms of configurations weighted bytheir SPOCK probabilities as in Fig. 3, only taking the stable config-urations as determined by N-body. This demonstrates that SPOCKis not significantly distorting the distribution of stable systems weare trying to characterize.We note that SPOCK does makes significant absolute errorsin assigning these stable configurations stability probabilities belowunity. Starting from 833 stable configurations (orange), SPOCK re-turns 402 effective samples (purple). However, it is not the absoluteprobabilities that matter, but rather the probabilities assigned to dif-ferent samples relative to one another, and whether they distort theoverall distribution. We see in Fig. 5 that the agreement is excellent,except for deviations in the free eccentricities between the inner pairof planets, where SPOCK skews the distribution by ≈
15% towardlower values.The final question is then how much SPOCK suppresses un-stable systems. As argued above, by assigning zero probability tosystems that do not survive its short 10 orbit integrations (and byrejecting orbit-crossing configurations), SPOCK can perfectly re-ject the vast majority of unstable configurations. The remainder arethe 9414 samples plotted in the gray histogram in Fig. 3, of whichunstable systems outnumber stable ones 10:1.SPOCK reduces these 8589 unstable configurations to 340effective samples. This final suppression of unstable systems by a In reality the orbital evolution is chaotic, leading to a range of equallyvalid instability times, but the probability of one realization going unstablein 10 orbits and another lasting 10 orbits is negligible (Hussain & Tamayo2019). MNRAS000
15% towardlower values.The final question is then how much SPOCK suppresses un-stable systems. As argued above, by assigning zero probability tosystems that do not survive its short 10 orbit integrations (and byrejecting orbit-crossing configurations), SPOCK can perfectly re-ject the vast majority of unstable configurations. The remainder arethe 9414 samples plotted in the gray histogram in Fig. 3, of whichunstable systems outnumber stable ones 10:1.SPOCK reduces these 8589 unstable configurations to 340effective samples. This final suppression of unstable systems by a In reality the orbital evolution is chaotic, leading to a range of equallyvalid instability times, but the probability of one realization going unstablein 10 orbits and another lasting 10 orbits is negligible (Hussain & Tamayo2019). MNRAS000 , 1–14 (2020) tability Constrained Characterization Inner-Pair Free Ecc. ( Z bc ) P D F Stable (N-body)[0.01, 0.07]Stable (SPOCK)[0, 0.06]
Outer-Pair Free Ecc. ( Z cd ) [0.01, 0.04][0.01, 0.04] Center-of-mass Ecc. ( e com ) [0, 0.04][0, 0.04] Figure 5.
Comparison of the eccentricity distributions of only the configurations that remain stable over 10 orbits in N-body integrations (orange) to the samehistograms re-weighted by their estimated probabilities of stability estimated by SPOCK (purple). This shows that SPOCK is not significantly skewing thedistribution of stable systems. factor of 25, while preserving the distribution of stable ones, allowsSPOCK to faithfully approximate the underlying distribution ofstable configurations in Fig. 3. This is a significant achievement in aproblem where unstable configurations initially outnumbered stableones 2500:1. Because the stability probabilities from SPOCK approximate trueprobabilities, we have used the terms interchangeably above forsimplicity. In detail, however, this is not the case. In particular,SPOCK probabilities are fit to capture the model’s uncertainties onthe distribution of example configurations it was trained on. Con-sider creating a large sample of orbital configurations distinct from(but generated in the same way, with different random seeds, as)SPOCK’s training examples (Tamayo et al. 2020). If one collectedall samples where SPOCK estimates probabilities near 60%, to ex-cellent approximation, 60% of those configurations would be stableover 10 orbits when run with N-body, since this was what themodel was trained to achieve.However, while SPOCK’s training set tries to cover the wholeparameter range for observed near-coplanar compact multi-planetsystems, any single system spans a much narrower range in phasespace, where SPOCK might do better or worse than average.In Fig. 6, we evaluate SPOCK’s performance on Kepler-23 bybinning configurations by their probability of stability as estimatedby SPOCK, and comparing it to the fraction of systems in the bin thatwere actually stable in N-body integrations (top panel). A perfectestimator would follow the 1:1 line in blue, and indeed, SPOCKachieves this when tested on holdout sets never seen during trainingof configurations drawn from its training dataset. For our sampleof Kepler-23 configurations, we see that it deviates from the trueprobabilities by (cid:46)
20% (bottom panel).Low probability configurations do not contribute much to thefinal distribution individually. However, the preponderance of un-stable systems means that the vast majority of configurations areassigned low stability probabilities. Over 80% of configurationshave SPOCK probabilities < 0.1, as reflected by the small Poissoncounting errors for the leftmost bins in Fig. 6.The bins in Fig. 6 were therefore chosen so that each bin con- F r a c t i o n s t a b l e Predicted Stability Probability (SPOCK) E rr o r Figure 6.
Predictions, binned by SPOCK’s estimated probabilities of sta-bility, vs. fraction of those configurations that were actually stable over 10 orbits when run with N-body. A perfect model would follow the blue one-to-one line. Error bars denote Poisson counting errors, and bins were chosenso that each bin contributes the same total SPOCK probability to the finalSPOCK distribution (purple histogram) in Fig. 3. tributes the same total probability to the final distribution. Thismakes it a useful diagnostic for analyzing errors in the stability con-strained distributions, since the errors in each bin then contributeequally to the purple histogram in Fig. 3. The 10 orbit integrations of the 9414 non-crossing configurationsthat survive > orbits (Sec. 4.4) required approximately 8500CPU hours. SPOCK predictions took ≈ times faster withSPOCK than N-body (Tamayo et al. 2020). By contrast, a sampleof systems that all go unstable within 10 orbits would take as long MNRAS , 1–14 (2020) Tamayo, Gilbertson, Foreman-Mackey with SPOCK as with direct N-body given that SPOCK is performinga 10 orbit N-body integration to generate its features.Given that one does not know the instability time distributionfor a given set of configurations ahead of time, it is difficult toestimate the time required through direct N-body integrations. Bycontrast, SPOCK requires a flat evaluation time of ≈ . SPOCK was trained on resonant and near-resonant, compact three-planet configurations with masses (cid:46) (cid:38) ◦ ). One particularly interesting application is to resonantchains, where the resonant structure and chaotic separatrices canrender wide swaths of parameter space unstable (e.g. Tamayo et al.2017). SPOCK’s training dataset puts pairs of planets in and nearmean motion resonances, but then initializes the third planet ran-domly. This means that very few training examples are in resonantchains, so SPOCK’s performance on such cases is unclear. At thesame time, such resonant chains, e.g. TRAPPIST-1 (Gillon et al.2017) and Kepler-223 (Mills et al. 2016) are both rare and par-ticularly valuable objects of study, justifying large investments ofCPU time with direct N-body. An interesting direction would be touse suites of N-body integrations for such systems, and use transferlearning (see, e.g., Weiss et al. 2016, for a review) to specializeSPOCK for stability classification for that particular case. Without a full analytic understanding of stability in compact mul-tiplanet systems, it is difficult to provide quantitative criteria forwhen stability constrained characterization will provide strong con-straints. Nevertheless, our partial understanding of the dynamicscan provide some rough guidelines.Quillen (2011) and Petit et al. (2020) argue analytically thatinstabilities in compact, initially circular , 3+ planet systems aredriven by 3-body resonances. Most, if not all, observed multiplanetsystems are at separations where 3-body resonances no longer over-lap (Quillen 2011; Petit et al. 2020) and are stable against this modeof instability for plausible planet masses. This is logical. Systemsunstable through this channel should have already eliminated them-selves, so we should expect to find systems that are stable at zeroeccentricity.Tamayo et al. (2020) suggest that short-lived eccentric config-urations instead destabilize through the overlap of 2-body MMRs,given that their model trained on systems in and near such 2-bodyMMRs generalizes well to uniformly distributed compact config-urations. These resonant widths grow with increasing planet massand orbital eccentricities (Wisdom 1980; Deck et al. 2013; Hadden& Lithwick 2018), so one expects threshold combinations of massesand eccentricities beyond which a system will be short-lived, givena set of orbital periods (which are often the best-known parameters).MMRs of a given order are more closely spaced (and thus eas-ier to overlap) at closer separations (e.g., for first-order resonances, the period ratios of the 7:6 and 6:5 MMR are much closer than thoseof the 3:2 and 2:1). Two-body MMRs always overlap as eccentrici-ties reach orbit-crossing configurations (Hadden & Lithwick 2018;Hadden 2019), so one can expect stable configurations to extend toa fraction of orbit-crossing values, where the exact threshold willbe correlated with the planet masses.Perhaps most importantly, constraints will depend strongly onthe relative spacing between adjacent planets. As argued in Sec. 3.2,numerical experiments suggest that mass constraints can vary byover an order of magnitude between two-planet cases and oneswith 3+ equally spaced planets. Thus, the strongest constraints willnot necessarily occur in the system with the closest pair of adja-cent planets, but in ones with comparably spaced additional bodies.Interestingly, transiting multiplanet systems exhibit a preference to-ward similar period ratios between adjacent planets (Millhollandet al. 2017; Weiss et al. 2018).We have seen that in cases like Kepler-23 where both pairs ofplanets have period ratios at or less than 3:2, one can obtain strongconstraints even for planet masses in the super-Earth regime. Thevalue of SPOCK is that one can quickly check stability constraintson estimated orbital parameters and masses, whereas N-body inte-grations are typically prohibitive.
Kepler-23 is a particularly informative compact three-planet system,where mass and eccentricity constraints from transit timing varia-tions (TTVs), transit durations and stability can be directly com-pared. Through this case study we have shown that in the most com-pact systems, stability can provide comparable constraints to TTVs,and much narrower upper limits than transit durations (Figs. 1 and2). This is particularly relevant for multiplanet systems observedover short time baselines with only a few transits, where TTVsrarely provide strong constraints (Hadden 2019).Such stability constrained characterization is typically pro-hibitive through direct N-body integrations, due to the long systemages and the large number of orbital configurations to be evaluated.We compared these computationally expensive N-body constraintsto much faster ones using the Stability of Planetary Orbital Config-urations Klassifier (SPOCK) from Tamayo et al. (2020).Stability constrained characterization with approximate sta-bility estimators like SPOCK is challenging for compact systemswhere most candidate configurations are unstable (drawing fromour prior resulted in a ratio of 2500 unstable systems for every sta-ble one). Given the preponderance of unstable configurations, inorder for residual false positives not to dominate the final posteriorsrequires strong rejection of unstable systems without skewing thedistribution of stable ones. We show that SPOCK indeed approxi-mately preserves the distribution of stable configurations (Fig. 5),and strongly suppresses unstable configurations to yield good agree-ment with N-body (Fig. 3).We provide several recommendations for stability constrainedcharacterization with SPOCK, most notably first rejecting any cross-ing configurations (Sec. 4.3), and weighting individual configura-tions by their respective probabilities of stability as estimated bySPOCK (Sec. 4.1).While contrasting constraints from transit timing variations,transit durations and stability is informative, we emphasize thatthey are not meant to be mutually exclusive, and are in fact com-plementary (Sec. 2.4). For example, as the GAIA mission increasesthe sample of stars with well measured densities beyond the aster-
MNRAS000
MNRAS000 , 1–14 (2020) tability Constrained Characterization oseismic sample (Van Eylen & Albrecht 2015), joint analyses ofmultiplanet systems with transit durations disfavoring low eccen-tricities and stability cutting off large values could provide improvedconstraints. ACKNOWLEDGEMENTS
We thank the anonymous reviewer for a careful report that improvedthe presentation in this manuscript. We would also like to thank Vin-cent van Eylen and Sam Hadden for sharing their data on Kepler-23for comparison. DT is indebted to the Hews Company for IT sup-port that enabled the completion of this work remotely. Support forDT was provided by NASA through the NASA Hubble Fellowshipgrant HST-HF2-51423.001-A awarded by the Space Telescope Sci-ence Institute, which is operated by the Association of Universitiesfor Research in Astronomy, Inc., for NASA, under contract NAS5-26555. C.G. acknowledges support from the Penn State Eberly Col-lege of Science, Department of Astronomy & Astrophysics, theCenter for Exoplanets and Habitable Worlds, the Center for Astro-statistics, Pennsylvania Space Grant Consortium, and the Institutefor Computational and Data Sciences ( http://ics.psu.edu/ ) atThe Pennsylvania State University, including the CyberLAMP clus-ter supported by NSF grant MRI-1626251, for providing advancedcomputing resources and services that have contributed to the re-search results reported in this paper. The citations in this paper havemade use of NASA’s Astrophysics Data System Bibliographic Ser-vices. This research has made use of the NASA Exoplanet Archive,which is operated by the California Institute of Technology, un-der contract with the National Aeronautics and Space Adminis-tration under the Exoplanet Exploration Program. This researchwas made possible by the open-source projects
REBOUND (Rein &Liu 2012), SPOCK (Tamayo et al. 2020),
Jupyter (Kluyver et al.2016), iPython (Pérez & Granger 2007), (Harris et al. 2020), and matplotlib (Hunter 2007; Droettboom et al. 2016).
DATA AVAILABILITY
To aid in future stability analyses, we provide an accompanyingrepository with Jupyter notebooks to reproduce the figures in thispaper at https://github.com/dtamayo/Stability-Priors . REFERENCES
Adams F. C., Laughlin G., 2003, Icarus, 163, 290Agol E., Steffen J., Sari R., Clarkson W., 2005, Monthly Notices of the RoyalAstronomical Society, 359, 567Barnes R., Greenberg R., 2006, Astrophysical Journall, 647, L163Bitsch B., Raymond S. N., Izidoro A., 2019, Astronomy & Astrophysics,624, A109Borucki W. J., et al., 2011, The Astrophysical Journal, 736, 19Burgasser A. J., Mamajek E. E., 2017, The Astrophysical Journal, 845, 110Chambers J. E., Wetherill G. W., Boss A. P., 1996, Icarus, 119, 261Chatterjee S., Ford E. B., Matsumura S., Rasio F. A., 2008, AstrophysicalJournal, 686, 580Chen T., Guestrin C., 2016, arXiv preprint arXiv:1603.02754Chen J., Kipping D., 2016, The Astrophysical Journal, 834, 17Cincotta P. M., Giordano C. M., Simó C., 2003, Physica D: NonlinearPhenomena, 182, 151Dawson R. I., Johnson J. A., 2012, The Astrophysical Journal, 756, 122Dawson R. I., Lee E. J., Chiang E., 2016, The Astrophysical Journal, 822,54 Deck K. M., Agol E., 2015, The Astrophysical Journal, 802, 116Deck K. M., Payne M., Holman M. J., 2013, Astrophysical Journal, 774,129Droettboom M., et al., 2016, matplotlib: matplotlib v1.5.1,doi:10.5281/zenodo.44579, http://dx.doi.org/10.5281/zenodo.44579
Faber P., Quillen A. C., 2007, Monthly Notices of the Royal AstronomicalSociety, 382, 1823Fabrycky D. C., et al., 2014, The Astrophysical Journal, 790, 146Fang J., Margot J.-L., 2012, The Astrophysical Journal, 761, 92Ford E. B., Quinn S. N., Veras D., 2008, The Astrophysical Journal, 678,1407Ford E. B., et al., 2012, The Astrophysical Journal, 750, 113Foreman-Mackey D., Hogg D. W., Lang D., Goodman J., 2013, Publicationsof the Astronomical Society of the Pacific, 125, 306Foreman-Mackey D., Luger R., Czekala I., Agol E., Price-Whelan A., Barclay T., 2020, exoplanet-dev/exoplanet v0.3.2,doi:10.5281/zenodo.1998447, https://doi.org/10.5281/zenodo.1998447
Fortney J. J., Marley M. S., Barnes J. W., 2007, The Astrophysical Journal,659, 1661Funk B., Wuchterl G., Schwarz R., Pilat-Lohinger E., Eggl S., 2010, As-tronomy & Astrophysics, 516, A82Gaspard P., 2005, in , Chaotic dynamics and transport in classical andquantum systems. Springer, pp 107–157Gilbert G. J., Fabrycky D. C., 2020, The Astronomical Journal, 159, 281Gillon M., et al., 2017, Nature, 542, 456Gladman B., 1993, Icarus, 106, 247Goldreich P., Lithwick Y., Sari R., 2004, arXiv preprint astro-ph/0405215Gratia P., Lissauer J. J., 2019, arXiv preprint arXiv:1908.01117Hadden S., 2019, Astronomical Journal, 158, 238Hadden S., Lithwick Y., 2014, The Astrophysical Journal, 787, 80Hadden S., Lithwick Y., 2017, The Astronomical Journal, 154, 5Hadden S., Lithwick Y., 2018, The Astronomical Journal, 156, 95Hansen B. M., Murray N., 2012, The Astrophysical Journal, 751, 158Hansen B. M., Murray N., 2013, The Astrophysical Journal, 775, 53Harris C. R., et al., 2020, Nature, 585, 357He M. Y., Ford E. B., Ragozzine D., Carrera D., 2020, arXiv preprintarXiv:2007.14473Henrard J., Lemaitre A., Milani A., Murray C., 1986, Celestial mechanics,38, 335Hogg D. W., Foreman-Mackey D., 2018, ApJS, 236, 11Holman M. J., Murray N. W., 2005, Science, 307, 1288Howe A. R., Burrows A., 2015, The Astrophysical Journal, 808, 150Huber D., et al., 2013, The Astrophysical Journal, 767, 127Hunter J. D., 2007, Computing In Science & Engineering, 9, 90Hussain N., Tamayo D., 2019, Monthly Notices of the Royal AstronomicalSociety, 491, 5258Izidoro A., Ogihara M., Raymond S. N., Morbidelli A., Pierens A., Bitsch B.,Cossou C., Hersant F., 2017, Monthly Notices of the Royal AstronomicalSociety, 470, 1750Izidoro A., Bitsch B., Raymond S. N., Johansen A., Morbidelli A., Lam-brechts M., Jacobson S. A., 2019, arXiv preprint arXiv:1902.08772Jontof-Hutter D., Rowe J. F., Lissauer J. J., Fabrycky D. C., Ford E. B.,2015, Nature, 522, 321Jontof-Hutter D., et al., 2016, The Astrophysical Journal, 820, 39Jurić M., Tremaine S., 2008, Astrophysical Journal, 686, 603Kipping D. M., 2014, Monthly Notices of the Royal Astronomical Society,440, 2164Kipping D. M., Dunn W. R., Jasinski J. M., Manthri V. P., 2012, MonthlyNotices of the Royal Astronomical Society, 421, 1166Kluyver T., et al., 2016, Positioning and Power in Academic Publishing:Players, Agents and Agendas, p. 87Lambrechts M., Morbidelli A., Jacobson S. A., Johansen A., Bitsch B.,Izidoro A., Raymond S. N., 2019, Astronomy & Astrophysics, 627,A83Laskar J., 1990, Icarus, 88, 266Laskar J., Petit A., 2017, Astronomy & Astrophysics, 605, A72MNRAS , 1–14 (2020) Tamayo, Gilbertson, Foreman-Mackey
Lin D., Ida S., 1997, The Astrophysical Journal, 477, 781Lissauer J. J., et al., 2011, The Astrophysical Journal Supplement Series,197, 8Lithwick Y., Xie J., Wu Y., 2012, The Astrophysical Journal, 761, 122Lopez E. D., Fortney J. J., 2013, The Astrophysical Journal, 776, 2MacDonald M. G., Dawson R. I., Morrison S. J., Lee E. J., Khandelwal A.,2020, The Astrophysical Journal, 891, 20Marchal C., Bozis G., 1982, Celestial Mechanics, 26, 311Marois C., Macintosh B., Barman T., Zuckerman B., Song I., Patience J.,Lafrenière D., Doyon R., 2008, science, 322, 1348Marzari F., Weidenschilling S. J., 2002, Icarus, 156, 570Millholland S., Wang S., Laughlin G., 2017, The Astrophysical JournalLetters, 849, L33Mills S. M., Fabrycky D. C., Migaszewski C., Ford E. B., Petigura E.,Isaacson H., 2016, Nature, 533, 509Morbidelli A., Batygin K., Brasser R., Raymond S., 2020, Monthly Noticesof the Royal Astronomical Society: LettersMurray N., Holman M., 1997, The Astronomical Journal, 114, 1246Neil A. R., Rogers L. A., 2018, The Astrophysical Journal, 858, 58Neil A. R., Rogers L. A., 2020, The Astrophysical Journal, 891, 12Nesvorn`y D., Morbidelli A., 2008, The Astrophysical Journal, 688, 636Nesvorn`y D., Vokrouhlick`y D., 2014, The Astrophysical Journal, 790, 58Ning B., Wolfgang A., Ghosh S., 2018, The Astrophysical Journal, 869, 5Obertas A., Van Laerhoven C., Tamayo D., 2017, Icarus, 293, 52Pérez F., Granger B. E., 2007, Computing in Science and Engineering, 9, 21Petit A. C., Laskar J., Boué G., 2017, Astronomy & Astrophysics, 607, A35Petit A. C., Laskar J., Boué G., 2018, Astronomy & Astrophysics, 617, A93Petit A. C., Pichierri G., Davies M. B., Johansen A., 2020, arXiv preprintarXiv:2006.14903Price E. M., Rogers L. A., Johnson J. A., Dawson R. I., 2015, The Astro-physical Journal, 799, 17Pu B., Wu Y., 2015, The Astrophysical Journal, 807, 44Quarles B., Quintana E. V., Lopez E. D., Schlieder J. E., Barclay T., 2017,arXiv preprint arXiv:1704.02261Quillen A. C., 2011, Monthly Notices of the Royal Astronomical Society,418, 1043Rasio F. A., Ford E. B., 1996, Science, 274, 954Rein H., Liu S.-F., 2012, Astronomy & Astrophysics, 537, A128Rein H., Tamayo D., 2015, Monthly Notices of the Royal AstronomicalSociety, 452, 376Rein H., Tamayo D., 2018, Monthly Notices of the Royal AstronomicalSociety, 473, 3351Rice D. R., Rasio F. A., Steffen J. H., 2018, Monthly Notices of the RoyalAstronomical Society, 481, 2205Rogers L. A., Bodenheimer P., Lissauer J. J., Seager S., 2011, The Astro-physical Journal, 738, 59Rosenthal M., et al., 2019, The Astronomical Journal, 158, 136Seager S., Kuchner M., Hier-Majumder C., Militzer B., 2007, The Astro-physical Journal, 669, 1279Sessin W., Ferraz-Mello S., 1984, Celestial mechanics, 32, 307Shen Y., Turner E. L., 2008, The Astrophysical Journal, 685, 553Simbulan C., Tamayo D., Petrovich C., Rein H., Murray N., 2017, MonthlyNotices of the Royal Astronomical Society, 469, 3337Smith A. W., Lissauer J. J., 2009, Icarus, 201, 381Steffen J. H., et al., 2013, Monthly Notices of the Royal AstronomicalSociety, 428, 1077Tamayo D., Triaud A. H., Menou K., Rein H., 2015, The AstrophysicalJournal, 805, 100Tamayo D., et al., 2016, The Astrophysical Journal Letters, 832, L22Tamayo D., Rein H., Petrovich C., Murray N., 2017, The AstrophysicalJournal Letters, 840, L19Tamayo D., et al., 2020, Proceedings of the National Academy of Sciences,117, 18194Valencia D., O’Connell R. J., Sasselov D., 2006, Icarus, 181, 545Van Eylen V., Albrecht S., 2015, The Astrophysical Journal, 808, 126Volk K., Gladman B., 2015, The Astrophysical Journal Letters, 806, L26Wang J. J., et al., 2018, The Astronomical Journal, 156, 192Ward W. R., 1988, Icarus, 73, 330 Weidenschilling S. J., Marzari F., 1996, Nature, 384, 619Weiss L. M., Marcy G. W., 2014, The Astrophysical Journal Letters, 783,L6Weiss L. M., et al., 2013, The Astrophysical Journal, 768, 14Weiss K., Khoshgoftaar T. M., Wang D., 2016, Journal of Big data, 3, 9Weiss L. M., et al., 2018, The Astronomical Journal, 155, 48Wisdom J., 1980, The Astronomical Journal, 85, 1122Wisdom J., 1986, Celestial mechanics, 38, 175Wisdom J., Holman M., 1991, Astronomical Journal, 102, 1528Wolfgang A., Rogers L. A., Ford E. B., 2016, The Astrophysical Journal,825, 19Wu Y., Lithwick Y., 2013, The Astrophysical Journal, 772, 74Xie J.-W., et al., 2016, Proceedings of the National Academy of Sciences,113, 11431Yalinewich A., Petrovich C., 2019, arXiv preprint arXiv:1907.06660Yoshinaga K., Kokubo E., Makino J., 1999, Icarus, 139, 328Zakamska N. L., Pan M., Ford E. B., 2011, Monthly Notices of the RoyalAstronomical Society, 410, 1895Zeng L., Sasselov D. D., Jacobsen S. B., 2016, The Astrophysical Journal,819, 127Zhou J.-L., Lin D. N., Sun Y.-S., 2007, The Astrophysical Journal, 666, 423This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000