On the survival of resonant and non-resonant planetary systems in star clusters
Katja Stock, Maxwell X. Cai, Rainer Spurzem, M.B.N. Kouwenhoven, Simon Portegies Zwart
MMNRAS , 1–17 (2020) Preprint 24 July 2020 Compiled using MNRAS L A TEX style file v3.0
On the survival of resonant and non-resonant planetarysystems in star clusters
Katja Stock (cid:63) † , Maxwell X. Cai , Rainer Spurzem , , ‡ , M.B.N. Kouwenhoven and Simon Portegies Zwart Astronomisches Rechen-Institut, Zentrum f¨ur Astronomie der Universit¨at Heidelberg, M¨onchhofstraße 12-14, D-69120 Heidelberg,Germany Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands National Astronomical Observatories and Key Laboratory of Computational Astrophysics, Chinese Academy of Sciences,20A Datun Road, Chaoyang District, Beijing 100101, P.R. China Kavli Institute for Astronomy and Astrophysics, Peking University, 5 Yi He Yuan Road, Haidian District, Beijing 100871,P.R. China Department of Physics, School of Science, Xi’an Jiaotong-Liverpool University, 111 Ren’ai Road, Suzhou Dushu Lake Science andEducation Innovation District, Suzhou Industrial Park, Suzhou 215123, P.R. China
Accepted 2020 July 7. Received 2020 June 5; in original form 2020 February 19
ABSTRACT
Despite the discovery of thousands of exoplanets in recent years, the number of knownexoplanets in star clusters remains tiny. This may be a consequence of close stellarencounters perturbing the dynamical evolution of planetary systems in these clusters.Here, we present the results from direct N -body simulations of multiplanetary systemsembedded in star clusters containing N = k, 16k, 32k, and 64k stars. The planetarysystems, which consist of the four Solar system giant planets Jupiter, Saturn, Uranus,and Neptune, are initialized in different orbital configurations, to study the effect of thesystem architecture on the dynamical evolution of the entire planetary system, and onthe escape rate of the individual planets. We find that the current orbital parametersof the Solar system giants (with initially circular orbits, as well as with present-dayeccentricities) and a slightly more compact configuration, have a high resilience againststellar perturbations. A configuration with initial mean-motion resonances of 3:2, 3:2,and 5:4 between the planets, which is inspired by the Nice model, and for which the twooutermost planets are usually ejected within the first yr, is in many cases stabilizeddue to the removal of the resonances by external stellar perturbation and by the rapidejection of at least one planet. Assigning all planets the same mass of 1 M Jup almostequalizes the survival fractions. Our simulations reproduce the broad diversity amongstobserved exoplanet systems. We find not only many very wide and/or eccentric orbits,but also a significant number of (stable) retrograde orbits.
Key words: planets and satellites: dynamical evolution and stability – galaxies:clusters: general – methods: numerical
Studies of nearby giant molecular clouds by Lada, Strom &Myers (1993) suggest that stars generally do not form in iso-lation but also in groups or stellar associations. If clusteredstar formation is the rule rather than the exception, there (cid:63)
E-mail: [email protected] † Fellow of the International Max Planck Research School forAstronomy and Cosmic Physics at the University of Heidelberg(IMPRS-HD) ‡ Research Fellow of Frankfurt Institute for Advanced Studies. is reason to believe that star clusters are promising targetsfor the detection of newborn planetary systems because starand planet formation are closely connected to each other.Despite the large number of 4158 extrasolar planets which were detected in the last 25 yr, only around 30 plan-ets ( < ) have been detected in star clusters so far, andonly one of them has been detected in a globular cluster (seetable 1 in Cai et al. 2019, for a complete list of planet de-tections in star clusters and their corresponding references). As of May 2020, according to the NASA Exoplanet Archivec (cid:13) a r X i v : . [ a s t r o - ph . E P ] J u l K. Stock et al.
Among those planets detected in star clusters are, for ex-ample, 13 planets around 11 stars in the Praesepe (M44)cluster (Quinn et al. 2012; Malavolta et al. 2016; Obermeieret al. 2016; Gaidos et al. 2017; Mann et al. 2017; Rizzuto etal. 2018; Livingston et al. 2019), six planets in four systemsin the Hyades cluster (Sato et al. 2007; Quinn et al. 2014;Mann et al. 2016) with one three-planet system (Mann etal. 2018), and five single-planet systems in the M67 cluster(Brucalassi et al. 2014, 2016, 2017). The origin for the peri-odic RV variations in the giant stars IC 4651 No. 9122, NGC2423 No. 3, and NGC 4349 No. 127, which are all located inan open cluster, is still under debate (Delgado Mena et al.2018). Brucalassi et al. (2017) find a comparable fraction ofgiant planets around stars in the cluster M67 than aroundfield stars but a significantly higher fraction of Hot Jupitersin the cluster compared to the field (see also Brucalassi et al.2016). Although the sample size in these studies is very smalland statistics should therefore be interpreted with caution,the “excess” of Hot Jupiters found in M67 is an indicationfor significant dynamical perturbations from neighbouringstars on the planets in the cluster.Clustered environments pose a threat already for theearly phases of planet formation. Protoplanetary discs maybe photoevaporated by the radiation of nearby massive stars(e.g. St¨orzer & Hollenbach 1999; Armitage 2000; Anderson,Adams & Calvet 2013; Facchini, Clarke & Bisbas 2016) ortruncated due to close encounters (e.g. Clarke & Pringle1993; Olczak, Pfalzner & Spurzem 2006; Portegies Zwart2016; Concha-Ram´ırez, Vaher & Portegies Zwart 2019). Buteven when a planetary system has successfully formed with-out major perturbations, its dynamical fate will still be de-termined by the host star’s position and motion inside thecluster and the properties of the cluster itself like its den-sity (denser clusters, and especially globular clusters, tendto have a more destructive effect on planetary systems thanloosely bound open clusters). Numerous studies have anal-ysed the effect of cluster environments on planetary systemsbeyond the protoplanetary disc phase (e.g. Malmberg et al.2007; Spurzem et al. 2009; Malmberg, Davies & Heggie 2011;Parker & Quanz 2012; Hao, Kouwenhoven & Spurzem 2013;Cai et al. 2017; Cai, Portegies Zwart & van Elteren 2018;Cai et al. 2019; Flammini Dotti et al. 2019; Fujii & Hori2019; van Elteren et al. 2019; Glaser et al. 2020).Spurzem et al. (2009) presented a set of dynamical starcluster models with a large number of planetary systems(consisting of one planet) fully included into the model;they showed that there is a constant rate of planets lib-erated as a result of stellar encounters; they also showedthat stellar encounters act like a diffusive process on plane-tary systems, where changes of semimajor axis and angularmomentum may be directed in both ways. Depending onthe details of the encounter, there is a net flux outward,giving the rate at which free-floating planets are created.There could also be a net flow to the inner boundary, i.e.planets accreted onto the central star, which was not dis-cussed in their paper. Li & Adams (2015) followed anotherapproach – a Monte Carlo model, in which many thousandsof encounters of single objects (single and binary stars) withplanetary systems were modelled. They were able to cover aparameter space substantially larger than that of Spurzemet al. (2009). However, Monte Carlo models suffer from theinaccuracy in the stochastic selection of encounter parame-
Table 1.
Initial conditions for the star cluster simulations.Star cluster 8k 16k 32k 64kNumber of stars 8000 16 000 32 000 64 000Total mass (M (cid:12) ) 4073 7939 16 302 32 619Half-mass radius (pc) 0.78 0.78 0.78 0.78Central density (M (cid:12) pc − ) 3906 6813 13 852 25153Initial tidal radius (pc) 22.58 28.20 35.84 45.16 t age [Myr]10 c [ M p c ] Praesepe Pleiades Hyades IC4651 Ruprecht147 M67 NGC6811 8k 16k 32k 64k100 < M tot /M < 300300 < M tot /M < 600600 < M tot /M < 1000 Figure 1.
Central density ( ρ c ) as a function of the cluster age τ age ,for our four simulated clusters (blue hexagons) at a simulationtime of 100 Myr and the observed clusters from fig. 1 in Fujii &Hori (2019). ters (the impact parameter and the velocity at infinity). Infigs 1 and 2 in Spurzem et al. (2009) one can see that thereal distribution of these parameters in a star cluster differsfrom a random selection, covering the available phase spaceequally.In the work of van Elteren et al. (2019), they adopteda different approach, in which planetary systems were inte-grated together with the stars in the cluster. To reduce thecomputational burden, planets in one system were not af-fecting the orbits of planets in another system. This still ledto an enormous computational burden, which resulted in arather limited parameter study.Note that also very low-mass particles, such as planetes-imals (asteroids and Kuiper belt objects) or comets (Oortcloud objects) are subject to these encounters (see e.g. Veraset al. 2020). A characterization of the importance of suchclose encounters on planetary systems with debris discs ispresented in Portegies Zwart & J´ılkov´a (2015). This pro-cess can lead to flybys of interstellar objects (Torres et al.2019) or to the capture of cometary objects into youngplanetary systems (e.g., Kouwenhoven et al. 2010; Perets& Kouwenhoven 2012; Wang et al. 2015a; Shu et al. 2020).This work is also closely related to a – somewhat less com-prehensive – study by Hands et al. (2019). It was suggestedthat the extraordinary asteroids 90377 Sedna was abductedfrom the debris disc of another star in such a close encounter(J´ılkov´a et al. 2016). The identification of ‘Oumuamua and2I/Borisov as possible interstellar objects in our Solar sys-tem has received much attention recently, and is connectedto the idea that young planetary systems are sources of free- MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs floating comets or planetesimals (see, e.g., Zheng, Kouwen-hoven & Wang 2015; Portegies Zwart et al. 2018; Handset al. 2019; ’Oumuamua ISSI Team et al. 2019; Pfalzner &Bannister 2019, and references therein).In this work, we study the effect of close stellar en-counters on the dynamical architectures of planets that areborn around stars in star clusters. We pay particular at-tention to the dependence on the initial orbital configura-tion of a planetary system before the first encounters withneighbouring stars take place. Our work differs from earlierworks in several aspects. (i) We do not only focus on theeffect of one single encounter on planetary systems but in-stead investigate the cumulative effect of several encounterson planetary systems by following their dynamical evolutionduring a significant fraction of time which they spend in thecluster. This again allows us to compare the distribution oforbital parameters at the end of our simulations with actualobserved properties of planetary systems that are in con-flict with current planet formation theories (e.g. eccentricor retrograde orbits). (ii) Our N -body approach enables arealistic representation of encounters between cluster mem-bers, while many previous works use a Monte Carlo approachwhich typically suffer from inaccuracy by randomly selectingencounter parameters equally from the available parameterspace. (iii) Using a hybrid N -body code allows us to put ev-ery planetary system in different initial configurations whilethe host star’s trajectory through the cluster and thus alsoall external perturbations on the planetary system are thesame for the different system architectures.This paper is organized as follows. Section 2 describesthe computational approach of the simulation of planetarysystems embedded in star clusters and specifies the initialconditions for the star cluster simulation and the simulationof the planetary systems. In Sec. 3 we present the results ofour simulations which are then discussed and summarizedin Sec. 4. Planetary systems evolve through secular evolution, the or-bits of planets being relatively stable for millions, and some-times tens of billions of orbits. Secular evolution is providedby mutual gravitational interaction between the planets, aswell as by external perturbation through passing stars, in astar cluster or in the Galaxy. However, stellar clusters evolvedifferently, namely through two-body relaxation and few-body encounters. Orbits of stars in the system are changedby these processes in less than a single orbital time-scale.The dynamical evolution of star clusters also exhibits de-terministic chaos, so that slightly different initial conditionscan lead to exponentially diverging outcomes in phase spacewithin less than one orbital time (see e.g. Miller 1964; Quin-lan & Tremaine 1992; Boekholt, Portegies Zwart & Valtonen2020).Therefore, a combined simulation of planetary systemsin star clusters is a challenge. The challenge lies not so muchin the different time-scales or hierarchical nature of some ob-jects (in this sense, close stellar binaries and planetary sys-tem are quite similar); rather the problem is to accurately follow resonant and secular effects in the internal evolutionof planetary systems. This is why we simulate star clusterand planetary systems using different simulation codes. Thisis feasible because we assume that the neighbouring stars inthe cluster affect the planets, but the planets have a negli-gible influence on the stellar kinematics.Although currently the decoupled, combined simula-tions of planetary systems and star clusters as describedearlier are state of the art, and fully coupled dynamical sim-ulations of planetary systems in star clusters have only beencarried out for single planetary systems (e.g., Spurzem et al.2009), in the future more development on that side wouldbe important. Using the current
LPS algorithm (see below)neglects the potential effect of more distant perturbers andalso tidal forces of the entire star cluster on the planetarysystem. Also very massive bodies being further away (e.g.stellar or intermediate mass black holes) could have an im-pact on planetary systems which are not taken into accounthere.We first simulate the stellar population in the star clus-ter using
NBODY6++GPU (Wang et al. 2015b, 2016) and in-tegrate the motion of its members inside the cluster us-ing the Hermite scheme.
NBODY6++GPU is a follow-up versionof
NBODY6 (Aarseth 1999) and
NBODY6++ , and has a signifi-cant speedup due to the usage of graphical processing units(GPUs) and parallelization of tasks through a message pass-ing interface (MPI). All required information such as mass,position, velocities, acceleration, and the first time derivativeof the acceleration of all cluster members in our simulationare stored at a high time resolution using the “block time-step” (BTS) storage scheme (Faber et al. 2010; Farr et al.2012; Cai et al. 2015). This scheme allows the reconstructionof stellar encounters in details when planetary systems areassigned to single stars in the cluster at a subsequent step(see Sec. 2.3). The data are stored in
HDF5 format to enablehigh-performance parallel access to the data.The dynamical evolution of the planetary systems issimulated using the LonelyPlanets Scheme ( LPS ). It isbased on the
AMUSE framework (Portegies Zwart 2011; Porte-gies Zwart & McMillan 2018) and uses rebound (Rein & Liu2012) to integrate the planets. Before integrating the plan-ets using the
IAS15 integrator (Rein & Spiegel 2015), allencounters with the next five neighbouring stars are derivedby interpolating the data of the corresponding stars fromthe BTS data (see Cai et al. 2017, 2019, for further expla-nations).
The simulated star clusters in this work contain 8000, 16 000,32 000, and 64 000 stars. We adopt the Kroupa (2001) ini-tial mass function in the mass range of 0.08-100 M (cid:12) . Thestars have an expected average mass of 0.509 M (cid:12) . We drawthe initial positions and velocities for the stars in our clus-ters from the Plummer (1911) model. The initial half-massradius for all clusters is r hm = . . We do not includeprimordial mass segregation and we do not include primor-dial binary systems. All initial parameters for the star clus-ter simulations are listed in Tab. 1. It should be noted that MNRAS000
The simulated star clusters in this work contain 8000, 16 000,32 000, and 64 000 stars. We adopt the Kroupa (2001) ini-tial mass function in the mass range of 0.08-100 M (cid:12) . Thestars have an expected average mass of 0.509 M (cid:12) . We drawthe initial positions and velocities for the stars in our clus-ters from the Plummer (1911) model. The initial half-massradius for all clusters is r hm = . . We do not includeprimordial mass segregation and we do not include primor-dial binary systems. All initial parameters for the star clus-ter simulations are listed in Tab. 1. It should be noted that MNRAS000 , 1–17 (2020)
K. Stock et al. these values are initial cluster properties. After a short phaseof core collapse the clusters rapidly expand and the centraldensities decrease significantly. Simulating the clusters for − Myr leads to central densities that correspond tothose of observed star clusters. Figure 1 shows the centraldensity of our simulated clusters after a simulation time of100 Myr in comparison to the actual observed clusters fromfig. 1 in Fujii & Hori (2019). The central density is not com-parable to the typical density our planetary systems experi-ence during their life in the cluster, and due to the onset ofmass segregation it is unlikely that our 1 M (cid:12) host stars re-main in the small but dense core of a Plummer model clusterfor a long time. The vast majority of all systems experiencemoderate stellar densities of up to a few hundred M (cid:12) pc − .The encounter time-scale (see equation 3 in Malmberget al. 2007) determined for our clusters is comparable to τ enc ≈ . Myr given in Malmberg et al. (2007). Instead ofusing m t = (cid:12) for the total mass of the stars involved inthe encounter as Malmberg et al. (2007) do, we use a valueof m t = . (cid:12) based on the assumption that our 1 M (cid:12) hoststars encounter stars with the average mass of ∼ . M (cid:12) .We set r min = au as the encounter distance to ensurecomparability with Malmberg et al. (2007). For the smallestcluster, we obtain a value of τ enc ≈ . Myr, and τ enc ≈ . Myrfor the largest cluster. This corresponds to encounter ratesof 0.3 (smallest cluster) and 0.9 (largest cluster) encountersper star per Myr.The Lagrangian radii containing different fractions ofthe total cluster mass as a function of time are shown inFig. 2 for the 8k cluster and give an overview of the evolutionof the entire star cluster. For comparison, the initial tidalradius r tid is plotted as well. The half-life of the cluster isdefined as the time at which the 50% Lagrangian radius(half-mass radius) and the tidal radius are equal.We use the standard definition of the tidal radius as r tid = R G (cid:18) M cl M G (cid:19) , (1)where R G and M G are the distance to the Galactic centreand the mass of the Galaxy contained inside R G ; M cl is thestar cluster mass. Our star clusters gradually lose mass overtime due to stellar evolution (see below) which results in ashrinking tidal radius over time.Near the tidal radius stars are typically only marginallybound to the cluster, and may escape from the cluster intothe tidal tails. In reality the situation is much more com-plex, since stars escape through Lagrangian points, and notall stars with positive energy (or outside r tid ) escape imme-diately, some of them may be retained by the cluster. Thisprocess is neatly described in the study of Ernst et al. (2008).Our star clusters are assumed to orbit the Galactic cen-tre in the solar neighbourhood, wherefore the tidal forcesof the galaxy on the cluster are the same as for the solarneighbourhood (Heisler & Tremaine 1986). The formation Note that this definition of the tidal radius is an operationalone, used for example in our N -body code; other definitions usethe truncation of the density profile (King 1962) or the distancebetween the Lagrange point and the cluster centre (see e.g. Justet al. 2009). These definitions differ from ours by a numericalfactor of order unity. r L a g r [ p c ] N = 8000
Figure 2.
Lagrangian radii r Lagr of the 8k star cluster, containingthe indicated fraction of mass, as a function of time. The blackdashed curve shows the initial tidal radius r tide . of tidal tails is observed in our simulations. We do not re-move stars from our simulations even when their position is r (cid:29) r tid . Therefore, we still keep track of the motion of starsin our simulation that have physically already left the clus-ter. Hence, the cluster dissolves faster than Fig. 2 suggest.We assume that the clusters will have reduced their centraldensity significantly after roughly 100 Myr. For example, inthe 8k cluster, almost 20% of the stars are already beyondthe tidal radius after 100 Myr and can be considered to haveleft the cluster. Therefore, we simulate the cluster environ-ment of our planetary systems only for this time span asmost strong encounters will have occurred by that time.The star cluster simulations with NBODY6++GPU includestellar evolution of single and binary stars –– they followthe evolution of masses and radii of all objects accordingto the recipes described in Hurley et al. (2005, and earliercitations of Hurley therein). Since we start without primor-dial binary stars, binary systems are rare — only a few dy-namically formed binaries are found. The stellar evolutionis implemented in the form of parametrized lookup tables;any mass-loss of stars or from binaries is assumed to leavethe cluster instantaneously; mass transfer in a binary is ap-proximately followed. The reader interested in more detailscould have a look into the DRAGON (million body) simu-lations (Wang et al. 2016). In recent years the stellar evo-lution prescriptions for N -body simulations are undergoingconsiderable changes, see for example Khalaj & Baumgardt(2015), Spera, Mapelli & Bressan (2015), and Banerjee etal. (2020) for an overview. The updates according to thatpaper are now also available in NBODY6++GPU , but have notyet been used for the simulations of this paper. Note thatwe select in the
LPS scheme only host stars for planets whichare close to one solar mass – therefore these systems are notsubject to any changes due to stellar evolution, given the rel-atively short time of simulation used here. In future models,we could also initialize planets around more massive stars,which would undergo changes due to stellar evolution (mass-loss due to expansion of the host star on the AGB leads toa loss of planets or wider orbits of those remaining).
MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs Table 2.
Initial orbital parameters of our planetary systems in the different configurations from Li & Adams (2015).Configuration Common parameters Jupiter Saturn Uranus NeptuneStandard e = i = ◦ a = . au a = . au a = . au a = . auCompact e = i = ◦ a = . au a = . au a = . au a = . auResonant e = i = ◦ a = . au a = . au a = . au a = . auEccentric i = ◦ a = . au a = . au a = . au a = . au e = . e = . e = . e = . Eccentric e = . i = ◦ a = . au a = . au a = . au a = . auMassive e = i = ◦ m pl = M Jup a = . au a = . au a = . au a = . au We aim to investigate how the initial configuration of theplanetary systems affects the dynamical evolution of theplanets that are born around stars in clustered environmentsand how it affects the likelihood of the individual planets tosurvive the first tens of millions of years in such a destruc-tive environment. For this purpose, we adopt the six differ-ent initial configurations of Li & Adams (2015) as startingpositions for the planets in our simulations (see Tab. 2).Li & Adams (2015) study scattering encounters betweenSolar system analogues and passing stars (single stars andbinary systems) and determine cross-sections for the dis-ruption of these planetary systems. Their planetary systemscontain the four Solar system giants Jupiter, Saturn, Uranus,and Neptune with their present-day masses. In the “stan-dard configuration” of Li & Adams (2015), they use the cur-rent semimajor axes of the planets but they assign circularand co-planar orbits. Inspired by the Nice model (Gomes etal. 2005), Li & Adams (2015) use two more compact con-figurations in which the three outer planets are closer toJupiter. The first one is referred to as “compact configura-tion”. Although these planetary systems are tightly packed,this configuration is fully stable over 100 Myr. In the sec-ond one, the four planets are in mutual mean-motion res-onance (MMR), wherefore this configuration is called “res-onant configuration”. In this configuration, Jupiter/Saturnand Saturn/Uranus are each in a 3:2 MMR while Uranusand Neptune are in a 5:4 MMR. See Li & Adams (2015) andthe references therein for a further discussion of this initialstate. The initial orbital angles in this configuration playa key role in the question whether not only this system isstable for a certain period of time but also the resonance an-gles librate for similar period of time. In our simulations, wecan fulfil the stability and resonance criterion usually longenough until the first encounters of neighbouring stars startto disturb the planetary systems and break the resonancesbetween the planets. However, it should be mentioned thatthis resonant configuration is generally highly unstable dueto its compactness, and usually at least one of the outerplanets is ejected rapidly when the initial orbital parameterare not chosen properly.Furthermore, Li & Adams (2015) use two eccentric con-figurations (referred to as “Eccentric e = . . While the first eccentricconfiguration is fully stable over 100 Myr, Neptune is ejectedin the second eccentric configuration after 5 Myr if we placethe system in isolation. Therefore, the second eccentric aswell as the resonant configuration both contain an internalinstability leading to a higher vulnerability against externalperturbations. The sixth investigated configuration in Li &Adams (2015) is referred to as “massive configuration” inwhich all planets have Jovian masses instead of their ac-tual masses. Despite the large masses of all four planets, theconfiguration is stable for at least 100 Myr.For all these six configurations, we distribute 200 identi-cal planetary systems around those stars in the cluster whosemasses are closest to 1 M (cid:12) . The host stars within one clus-ter simulation are therefore the same for each configuration.This allows us to work out the differences in vulnerabilityin the clustered environment between those initial configu-rations due to the different positions of the host stars in thecluster. The number of 200 planetary systems per clusterand per configuration is a compromise between computa-tional costs and the possibility to do proper statistics aboutour sample. Since we simulate 200 planetary systems in sixdifferent configurations in all four star clusters, we have atotal number of 4800 different planetary system simulations.On grounds of efficiency, our simulations are therefore car-ried out using the simulation monitor SiMon (Qian et al.2017).Each planetary system is integrated for 100 Myr (as dis-cussed in Sec. 2.2). Planets that are excited to an eccentricity e > . are considered as having been ejected from the sys-tem and are removed from the simulation. The mass-loss ofthe ∼ (cid:12) host stars is negligible during the main-sequencephase and especially during the first 100 Myr which is whyit is not taken into account for the dynamical evolution ofthe planetary systems. For each configuration, we simulate 200 identical planetarysystems and distribute them around ∼ (cid:12) host stars. Aswe do not include primordial mass segregation in our clus-ters the positions of the host stars (and therefore the stellardensities the planetary systems experience) in our clusters MNRAS000
Initial orbital parameters of our planetary systems in the different configurations from Li & Adams (2015).Configuration Common parameters Jupiter Saturn Uranus NeptuneStandard e = i = ◦ a = . au a = . au a = . au a = . auCompact e = i = ◦ a = . au a = . au a = . au a = . auResonant e = i = ◦ a = . au a = . au a = . au a = . auEccentric i = ◦ a = . au a = . au a = . au a = . au e = . e = . e = . e = . Eccentric e = . i = ◦ a = . au a = . au a = . au a = . auMassive e = i = ◦ m pl = M Jup a = . au a = . au a = . au a = . au We aim to investigate how the initial configuration of theplanetary systems affects the dynamical evolution of theplanets that are born around stars in clustered environmentsand how it affects the likelihood of the individual planets tosurvive the first tens of millions of years in such a destruc-tive environment. For this purpose, we adopt the six differ-ent initial configurations of Li & Adams (2015) as startingpositions for the planets in our simulations (see Tab. 2).Li & Adams (2015) study scattering encounters betweenSolar system analogues and passing stars (single stars andbinary systems) and determine cross-sections for the dis-ruption of these planetary systems. Their planetary systemscontain the four Solar system giants Jupiter, Saturn, Uranus,and Neptune with their present-day masses. In the “stan-dard configuration” of Li & Adams (2015), they use the cur-rent semimajor axes of the planets but they assign circularand co-planar orbits. Inspired by the Nice model (Gomes etal. 2005), Li & Adams (2015) use two more compact con-figurations in which the three outer planets are closer toJupiter. The first one is referred to as “compact configura-tion”. Although these planetary systems are tightly packed,this configuration is fully stable over 100 Myr. In the sec-ond one, the four planets are in mutual mean-motion res-onance (MMR), wherefore this configuration is called “res-onant configuration”. In this configuration, Jupiter/Saturnand Saturn/Uranus are each in a 3:2 MMR while Uranusand Neptune are in a 5:4 MMR. See Li & Adams (2015) andthe references therein for a further discussion of this initialstate. The initial orbital angles in this configuration playa key role in the question whether not only this system isstable for a certain period of time but also the resonance an-gles librate for similar period of time. In our simulations, wecan fulfil the stability and resonance criterion usually longenough until the first encounters of neighbouring stars startto disturb the planetary systems and break the resonancesbetween the planets. However, it should be mentioned thatthis resonant configuration is generally highly unstable dueto its compactness, and usually at least one of the outerplanets is ejected rapidly when the initial orbital parameterare not chosen properly.Furthermore, Li & Adams (2015) use two eccentric con-figurations (referred to as “Eccentric e = . . While the first eccentricconfiguration is fully stable over 100 Myr, Neptune is ejectedin the second eccentric configuration after 5 Myr if we placethe system in isolation. Therefore, the second eccentric aswell as the resonant configuration both contain an internalinstability leading to a higher vulnerability against externalperturbations. The sixth investigated configuration in Li &Adams (2015) is referred to as “massive configuration” inwhich all planets have Jovian masses instead of their ac-tual masses. Despite the large masses of all four planets, theconfiguration is stable for at least 100 Myr.For all these six configurations, we distribute 200 identi-cal planetary systems around those stars in the cluster whosemasses are closest to 1 M (cid:12) . The host stars within one clus-ter simulation are therefore the same for each configuration.This allows us to work out the differences in vulnerabilityin the clustered environment between those initial configu-rations due to the different positions of the host stars in thecluster. The number of 200 planetary systems per clusterand per configuration is a compromise between computa-tional costs and the possibility to do proper statistics aboutour sample. Since we simulate 200 planetary systems in sixdifferent configurations in all four star clusters, we have atotal number of 4800 different planetary system simulations.On grounds of efficiency, our simulations are therefore car-ried out using the simulation monitor SiMon (Qian et al.2017).Each planetary system is integrated for 100 Myr (as dis-cussed in Sec. 2.2). Planets that are excited to an eccentricity e > . are considered as having been ejected from the sys-tem and are removed from the simulation. The mass-loss ofthe ∼ (cid:12) host stars is negligible during the main-sequencephase and especially during the first 100 Myr which is whyit is not taken into account for the dynamical evolution ofthe planetary systems. For each configuration, we simulate 200 identical planetarysystems and distribute them around ∼ (cid:12) host stars. Aswe do not include primordial mass segregation in our clus-ters the positions of the host stars (and therefore the stellardensities the planetary systems experience) in our clusters MNRAS000 , 1–17 (2020)
K. Stock et al. T [Myr] F r a c t i o n o f s u r v i v i n g p l a n e t s JupiterSaturnUranusNeptuneAverage
Figure 3.
The survival fractions for the Solar system giant planets as a function of time for the six different initial configurations in a16k Plummer model star cluster. The black dotted curves represent the overall survival fraction averaged over the four planets. are random. However, to ensure comparability between thedifferent initial configurations we use the same 200 host starsfor all planetary systems of the same cluster.In this work we define a planet having “survived” whenit has not been ejected from the planetary system duringthe course of the simulation. This means that the planet’seccentricity has been e ≤ . for at least 100 Myr. For thedetermination of the survival fraction of a certain kind ofplanet we average over all 200 planets of the same type inthe same star cluster.An inspection of the survival fraction as a function oftime for the four different planets in our systems revealslarge differences between the initial configurations. Figure 3shows the survival fraction for all six configurations as afunction of time for the 16k cluster. In Fig. 4, the fraction ofsurviving planets after 100 Myr is plotted against the num-ber of stars, N , in the host star cluster, for the six initialorbital configurations.In all of our simulations, Jupiter is the planet with thehighest survival probability. The reason for this is two-fold:Jupiter is not only the most massive planet (in five of our sixconfigurations) but it is also the innermost planet, so thatits binding energy is by far the largest. The same reason-ing explains why Saturn is usually the second most resistantplanet. Although Neptune is the outermost planet, its bind-ing energy is somewhat larger than that of Uranus due to itslarger mass. This is why in most of our simulations Neptuneis slightly more likely to survive than Uranus.The survival fractions after 100 Myr for the planets instandard configuration in the 16k cluster are 88.0%, 74.5%,62.5%, and 66.0% for Jupiter, Saturn, Uranus, and Nep-tune, respectively. Starting the planets in the 16k cluster ina more compact configuration does not change these valuessignificantly as can be seen in Fig. 3. Also the first eccentricconfiguration in which the planets were assigned their trueeccentricities does not differ significantly from the standard and compact case. The overall survival fraction (averagedover all four planets) is around 72% for these three configu-rations in the 16k cluster.While the differences in the survival fractions of thestandard and compact configurations are negligible, the sur-vival fractions in the resonant case differ significantly fromthose in the compact configuration. Only the fraction of sur-viving Jupiters is comparable to the other configurationsand is 90.0% in the 16k cluster. For Saturn, the survivalfraction decreases from 74.5% in the standard case to 49.0%in the resonant configuration. However, the percentage ofsurviving Uranus- and Neptune-like planets is much lowerand is only 4.0% and 6.5%, respectively. The overall sur-vival fraction in the resonant case is only 37.4% which isthe lowest value for all six configurations in the 16k cluster.Although all planets have initially circular orbits, the effectof planet–planet interaction in this configuration is very de-structive. Due to the compactness of the planetary system,the system is only long-term stable on time-scales of sev-eral ten thousand years. The first encounters have usuallyalready occurred at that time, removing the system fromresonances and exciting the orbital parameters of some orall planets. Uranus and Neptune are the most vulnerableplanets in this configuration. In none of our simulations, allfour planets survived. Usually either Uranus or Neptune (orboth) is ejected latest after several million years. Only in2 out of 800 simulations of the resonant case, Uranus andNeptune survived together. In both cases, Saturn is ejectedwithin the first two million years.In the second eccentric configuration, the survival frac-tions of Jupiter and Saturn are relatively unaffected bythe larger initial eccentricities of all four planets. Only forUranus and Neptune, the differences are significant com-pared to the standard, compact, and first eccentric case.In the 16k cluster, the survival fractions after 100 Myr dropdown to 35% and 48% for Uranus and Neptune, respectively. MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs Standard Compact Resonant
Eccentric1
Eccentric2
Massive N F r a c t i o n o f s u r v i v i n g p l a n e t s JupiterSaturnUranusNeptuneAverage
Figure 4.
The survival fractions for the Solar system giant planets as a function of the number of stars, N , in the host star cluster at t = Myr.
Although Neptune is the outermost planet, it has a signif-icant higher chance to survive in this eccentric planetarysystem than Uranus. In this configuration, Uranus’ fate ismainly determined by secular evolution. Since the planetsalready have an initial eccentricity of e = . , it requires lessangular momentum transfer to another star to trigger de-structive interactions between the planets. Due to its posi-tion between Saturn and Neptune and the fact that it hasthe smallest mass, Uranus is easily excited to highly eccen-tric orbits which often leads to the ejection of the planet.The fractions of surviving planets in the massive config-uration reveal not only the importance of the planetary massduring stellar encounters but also its role during secular evo-lution. The overall survival fraction in the 16k cluster dropsfrom 72.8% in the standard configuration to 63.9% in themassive configuration. While Jupiter had by far the largestlikelihood for survival in the other configurations, the differ-ences between the planets in the massive configuration aresignificantly smaller. The survival fraction for Jupiter in the16k cluster is 88.0% in the standard and compact case, butonly 71.9% in the massive configuration. For Saturn, Uranusand Neptune the survival rates in the massive configurationare all around 60% in the 16k cluster.Our simulated clusters all have the same initial half-mass radius but differ in central density. Therefore, the sur-vival fractions for the different configurations also dependon the number of stars in the host cluster. In general, thesurvival fractions for the different planets decrease with in-creasing stellar density due to an increasing number of closeencounters between cluster members. However, the effect ofan increasing stellar density is larger on the outer planetsof the system since they are more easily liberated by an-other star due to their smaller gravitational binding energy.While the survival fractions for Jupiter, Saturn, Uranus, andNeptune in the standard configuration in the 8k cluster are87.5%, 75.5%, 71.5%, and 70.5%, respectively, these values decrease to 83.0%, 60.0%, 45.0%, and 44.0%, respectively,in the 64k cluster (see Fig. 4). In table 2 in Li & Adams (2015), the authors provide theirejection cross-sections in units of au . In order to normalizethis to obtain an escape or survival fraction, one needs toknow the maximum impact parameter chosen in their mod-els. However, their maximum impact parameter is variable,depending on parameters (e.g. 10 times the semimajor axisof a stellar binary which encounters a planetary system, butnot more than 1000 au). To be able to compare our resultswith those of Li & Adams (2015), we adopt p max = aufor the normalization of their cross-sections. Table 3 lists oursurvival fractions in percentage after integrating the plane-tary systems in the 8k cluster in four of the six different ini-tial configurations for 100 Myr. We assume that our smallestcluster is most similar to the cluster environment simulatedin Li & Adams (2015). However, our models are different inthree aspects — (1) the distribution of impact parametersand relative velocities of encounters is very different to theone assumed in Monte Carlo simulations of encounters asthey did; (2) Li & Adams (2015) stop the planetary systemmodel after the encounter, while we continue all planetarysystems for the entire simulation time of 100 Myr and findmany delayed unstable systems, which reduce the survivalfraction; (3) our simulations take into account the cumula-tive effect of several encounters. Due to these differences, wefind much more ejections of planets in our simulations andhave significantly smaller survival fractions for each planettype than Li & Adams (2015).The ejection of one or several planets can occur eitherduring or directly after the encounter (prompt ejection), orat a later time due to secular evolution (delayed ejection). MNRAS000
Although Neptune is the outermost planet, it has a signif-icant higher chance to survive in this eccentric planetarysystem than Uranus. In this configuration, Uranus’ fate ismainly determined by secular evolution. Since the planetsalready have an initial eccentricity of e = . , it requires lessangular momentum transfer to another star to trigger de-structive interactions between the planets. Due to its posi-tion between Saturn and Neptune and the fact that it hasthe smallest mass, Uranus is easily excited to highly eccen-tric orbits which often leads to the ejection of the planet.The fractions of surviving planets in the massive config-uration reveal not only the importance of the planetary massduring stellar encounters but also its role during secular evo-lution. The overall survival fraction in the 16k cluster dropsfrom 72.8% in the standard configuration to 63.9% in themassive configuration. While Jupiter had by far the largestlikelihood for survival in the other configurations, the differ-ences between the planets in the massive configuration aresignificantly smaller. The survival fraction for Jupiter in the16k cluster is 88.0% in the standard and compact case, butonly 71.9% in the massive configuration. For Saturn, Uranusand Neptune the survival rates in the massive configurationare all around 60% in the 16k cluster.Our simulated clusters all have the same initial half-mass radius but differ in central density. Therefore, the sur-vival fractions for the different configurations also dependon the number of stars in the host cluster. In general, thesurvival fractions for the different planets decrease with in-creasing stellar density due to an increasing number of closeencounters between cluster members. However, the effect ofan increasing stellar density is larger on the outer planetsof the system since they are more easily liberated by an-other star due to their smaller gravitational binding energy.While the survival fractions for Jupiter, Saturn, Uranus, andNeptune in the standard configuration in the 8k cluster are87.5%, 75.5%, 71.5%, and 70.5%, respectively, these values decrease to 83.0%, 60.0%, 45.0%, and 44.0%, respectively,in the 64k cluster (see Fig. 4). In table 2 in Li & Adams (2015), the authors provide theirejection cross-sections in units of au . In order to normalizethis to obtain an escape or survival fraction, one needs toknow the maximum impact parameter chosen in their mod-els. However, their maximum impact parameter is variable,depending on parameters (e.g. 10 times the semimajor axisof a stellar binary which encounters a planetary system, butnot more than 1000 au). To be able to compare our resultswith those of Li & Adams (2015), we adopt p max = aufor the normalization of their cross-sections. Table 3 lists oursurvival fractions in percentage after integrating the plane-tary systems in the 8k cluster in four of the six different ini-tial configurations for 100 Myr. We assume that our smallestcluster is most similar to the cluster environment simulatedin Li & Adams (2015). However, our models are different inthree aspects — (1) the distribution of impact parametersand relative velocities of encounters is very different to theone assumed in Monte Carlo simulations of encounters asthey did; (2) Li & Adams (2015) stop the planetary systemmodel after the encounter, while we continue all planetarysystems for the entire simulation time of 100 Myr and findmany delayed unstable systems, which reduce the survivalfraction; (3) our simulations take into account the cumula-tive effect of several encounters. Due to these differences, wefind much more ejections of planets in our simulations andhave significantly smaller survival fractions for each planettype than Li & Adams (2015).The ejection of one or several planets can occur eitherduring or directly after the encounter (prompt ejection), orat a later time due to secular evolution (delayed ejection). MNRAS000 , 1–17 (2020)
K. Stock et al.
Table 3.
Survival fractions (in percent) at t = Myr in the 8k cluster, in comparison to the results of Li & Adams (2015).Jupiter Saturn Uranus Neptune Jupiter Saturn Uranus Neptune8k Li & Adams (2015)Standard 87.5 75.5 71.5 70.5 98.5 96.6 92.8 88.7Compact 85.0 75.5 68.5 72.0 98.2 96.7 94.3 90.6Resonant 88.5 56.5 6.50 5.0 97.6 96.0 93.9 94.0Massive 74.0 72.0 67.5 70.5 97.6 96.2 92.2 89.5
Table 4.
Fractions of prompt and delayed ejections in the stan-dard configuration of the 16k cluster for the different planet types.Planet Prompt ejection Delayed ejectionJupiter 88% 12%Saturn 61% 39%Uranus 52% 48%Neptune 59% 41%
For the standard case in the 16k cluster, we find a frac-tion of ∼ prompt ejections and ∼ delayed ejections(see Table 4 for a distinction between the different planettypes). However, both events (but especially the latter case)are not well defined in our simulations since the planetarysystems are continuously perturbed by other stars. In manycases, where planetary systems are already moderately orhighly excited, the true source for a planet’s ejection — sec-ular evolution or the next external perturbation — cannotbe clearly identified. The mentioned fraction for the delayedejection should therefore be treated with caution. We of-ten see a strong planet–planet interaction subsequent to anencounter which leaves the system in a highly vulnerablestate. It then only requires a very weak perturbation by an-other star to eject some of the planets which would not havebeen strong enough to disrupt the planetary system withoutthe previous excitation. Those events are counted as delayedejection even though they result from the combined effect ofsecular evolution and (another) prompt ejection due to thenext encounter.Fujii & Hori (2019) perform N -body simulations of dif-ferent cluster types and use a semi-analytical approach forthe calculation of the fraction of ejected planets. The clustermodel which is closest to one of our clusters is a high-densityKing-model cluster (King 1966) with N = and W = .Using the power-law function from equation 10 in Fujii &Hori (2019) and the corresponding best-fitting parameterfor G-type stars, one obtains survival fractions [ − f ejc ( a ) ]for the standard configuration of 93% (Jupiter), 90% (Sat-urn), 82% (Uranus), and 76% (Neptune). These values arehigher than the results for our smallest cluster. The impor-tant difference between our work and Fujii & Hori (2019) isnot so much the different cluster models but the fact that wealso take into account delayed ejections due to planet–planetscattering (for which the multiplicity of our planetary sys-tems plays a crucial role) and the possibility of several strongencounters.In the standard configuration of our 16k cluster the frac-tion of systems in which at least one planet is immediatelyejected after an encounter (regardless of the intruder’s mass) Table 5.
Fractions of planetary systems in which at least oneplanet is ejected during the simulation.Configuration 8k 16k 32k 64kStandard 34% 42% 42% 62%Compact 34% 38% 38% 63%Resonant 100% 100% 100% 100%Eccentric is 26 % . This value is comparable to the study of Malmberg,Davies & Heggie (2011) where they find fractions between15 and 31 % for a mass range of . − . M (cid:12) for the intruderstar. Malmberg, Davies & Heggie (2011) also determine thefractions of systems from which at least one planet has beenejected within 100 Myr after the encounter and find fractionsof 47–69 % for flybys of stars with masses of . − . M (cid:12) . Wefind a corresponding value of 42 % in the standard configu-ration of the 16k cluster (see Tab. 5 for other configurationsand other cluster sizes). The lower value likely stems fromthe shorter integration time. Although we simulate the plan-etary systems for 100 Myr, the remaining simulation timeafter the first strong encounter is shorter wherefore our val-ues cannot be directly compared with those from Malmberg,Davies & Heggie (2011).Although Malmberg et al. (2007) use N -body simula-tions, a direct comparison is also difficult as they only studyencounters between stars but do not explicitly analyze plan-etary systems. Furthermore, their studied clusters are rathersmall in terms of stellar members. They define a star thathas never been part of a binary system or has never under-gone any close encounters with other stars as “singleton”.We calculate the fraction of singletons in our simulationsand obtain values of 50% (8k cluster), 32% (16k cluster),20% (32k cluster), and 6% (64k cluster). Malmberg et al.(2007) provide the fraction of singletons for different half-mass radii and different numbers of cluster members. Tak-ing their largest cluster ( N = ) as reference, our resultsare most similar to the range of initial half-mass radii of0.38–1.69 pc.The most consistent simulations of star clusters withplanetary systems so far have been performed by van Elterenet al. (2019). They adopted the initial conditions from earliersimulations that tried to match the mass and size distribu-tions of circumstellar discs in the Orion Trapezium cluster(Portegies Zwart 2016). In the follow-up calculations by vanElteren et al. (2019), the discs were replaced by planetarysystems, selected according to the Oligarchic growth model MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs (Kokubo & Ida 1998). The parameter search was limited toa cluster of 1500 stars. The initial conditions were generatedfrom the simulation of a star cluster with circumstellar discsof 400 au each for 1 Myr during which the discs were trun-cated and harassed by passing stars. In that time frame,the cluster evolved and the discs were affected by passingstars but not by internal processes. After 1 Myr, a total of977 stars remained bound in the cluster, 512 of which re-ceived a planetary system. The calculation was performedwith 2522 planets with a total mass of 3527 Jovian masses.At an age of 11 Myr, 10 Myr after the birth of the cluster,2165 planets were still bound to their host star: 16.5% ofthe planets became unbound. The majority ( ∼ %) of theejected planets promptly escaped the cluster, the rest lingersaround for at least half a million years before escaping thecluster potential. a - e space We plot the eccentricity as a function of the semimajor axisof the planets for the different cluster sizes in Fig. 5 for the 8kcluster, and in Figs A1, A2, and A3 in the Appendix for the16k, 32k, and 64k clusters, respectively. These figures clearlyshow a trend with increasing cluster size. In the 8k clustermost planets are only excited in eccentricity and just a few ofthem migrate into wider (or sometimes tighter) orbits. Thefraction of highly eccentric orbits and planets that undergosignificant orbital migration increases with increasing clusterdensity. The planets’ distribution in the a - e space is thereforewider for our larger clusters.Having a look on the different initial orbital configu-rations reveals large differences in the a - e space at the endof our simulations. In the 8k cluster (see Fig. 5), the stan-dard and compact configuration look similar. Most planetsroughly retain their initial semimajor axis for 100 Myr. Onlya few have migrated to larger semimajor axes and even lessto tighter orbits (mainly Jovian-like planets, but there is alsoone Uranus-like planet in the standard case with a < au).Especially the outermost planets tend to migrate to verywide orbits of more than 100 au. In the standard configu-ration of the 8k cluster, even one Saturn-analogue can befound beyond a > au. Due to the initially smaller semi-major axis of the outer three planets in the compact config-uration, we observe a smaller number of planets on orbitswith a > au. In the compact configuration of the 16k clus-ter (see Fig. A1 in the Appendix) we can find in general morewide-orbit planets than in the 8k cluster but also wide-orbitplanets with eccentricities e < . which are missing in thestandard configuration. The fraction of planets which re-mains unaffected in their orbital parameters is lower in the32k and 64k clusters.The distribution of planets in the a - e space looks differ-ent for the resonant configuration. In the 8k cluster, mostJovian-like planets were at least excited in eccentricity andsome also migrated within the system (mainly to wider or-bits). None of the Saturn-like planets can retain its initialsemimajor axis and eccentricity, and a clear trend towardswide, eccentric orbits is observable. The few Uranus-likeplanets which survived for 100 Myr all have wide and/oreccentric orbits. Four of these planets have a > au andone even has a > au. All of these four planets have ec-centricities of e > ∼ . . While all Uranus-like planets failed to keep their initial semimajor axis, three of the Neptune-likeplanets succeeded in doing so. However, all of them were atleast slightly excited in eccentricity (as well as the Uranus-like planets). In all of these three systems Uranus was ejectedduring the first few tens of thousands of years after a rela-tively short interaction with Neptune before the first strongencounter happened. Due to the encounters with Neptune,all three Uranus’s migrated to an orbit with a semimajoraxis smaller than that of Jupiter which led to the ejectionof Uranus within the subsequent tens of thousands of years.This fortunate circumstance made the planetary system ro-bust enough to withstand the gravitational perturbations byother stars during the remaining 99 Myr.The three orbital parameters a , e , and i of one of thesethree planetary systems as a function of time are shown inFig. 6. The time is plotted in logarithmic scale to highlightthe planet–planet scattering during the first 100,000 yr. Weadditionally plot the distance of the host star to the clus-ter centre and the distance to the next stellar perturber ingray in the top and middle panel to illustrate the interactionwith the cluster. The cumulative gravitational effects of sev-eral neighbours in distances between 13 000 and 23 000 auremove the resonance of Uranus and Neptune within thefirst 10 000 yr which causes them to slightly interact witheach other and to change their orbital position for a fewhundred years. After migrating back to the second outer-most position Uranus is already excited in eccentricity. Thesubsequent interaction with Jupiter and Saturn and the si-multaneous close approach of a neighbouring star leads tothe prompt ejection of Uranus and the removal of the re-maining resonances in the system. Due to this circumstance,the remaining three planets form a stable system and stayrelatively unperturbed for the rest of the simulation eventhough neighbouring stars closely approach the system sev-eral times.The a - e space for the first eccentric configuration in the8k cluster looks similar to the standard and compact config-uration but there are slight differences. On the one hand, thenumber of Jupiters that have high eccentricities is reduced.On the other hand, the number of Uranus- and Neptune-likeplanets with high eccentricities is larger in the first eccentricconfiguration. While Jupiter, Saturn, and Uranus all startat eccentricities of e ≈ . there are Jupiters and Saturnsthat end up at nearly-circular orbits. However, this is not thecase for Uranus. While some of them keep their initial eccen-tricity, there are no Uranus-like planets that have reducedit after 100 Myr. Most of them are significantly excited ineccentricity.Increasing the initial eccentricity of all planets to e = . makes a large difference in the outcome of our simulations.Especially Uranus and Neptune cover a much wider range inthe a - e space after 100 Myr compared to the first eccentricconfiguration. On the other hand, most of our Jupiters andSaturns “fall back” to circular or almost circular orbits dur-ing the simulation which can be explained by the exchangeof angular momentum between the planets during close en-counters. This effect can be seen in Fig. 7 where Neptune isejected at t = Myr. A 6.9 M (cid:12) star approaches the plane-tary system down to a distance of 310 au and ejects Neptuneout of the system by “kicking” it to the inner regions of theplanetary systems where it transfers angular momentum to
MNRAS , 1–17 (2020) K. Stock et al. Eccentric2 10 Massive a [AU] e JupiterSaturnUranusNeptune
Figure 5.
The a - e space for all planets in the 8k star cluster which are not ejected from their host planetary system at t = Myr, forthe six different initial configurations. A video showing the a - e space for the course of the simulation is available on our Silkroad projectteam webpage: http://silkroad.bao.ac.cn/silkroad-save/a_e_space_N8k.mp4 . T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]051015 i [ d e g ] Jupiter Saturn Uranus Neptune
Figure 6.
The orbital parameters a (top), e (middle), and i (bot-tom) of a resonant planetary system from the 8k cluster as afunction of time. The time is plotted in logarithmic scale. Thescales on the right side correspond to the gray lines in the plotswhich represent the distance of the host star to the cluster centre(top) and the distance to the closest perturber (middle). Jupiter and Saturn. The eccentricity of both planets subse-quently decreases.The massive configuration is more difficult to be excitedin orbital parameters which can be seen the right bottompanel of Fig. 5. There is a clear distinction between thoseplanets that are only slightly perturbed, which are thosewith eccentricities below 0.2, and those which have been suf-ficiently perturbed to trigger fatal planet-planet scattering.Due to the equal mass of all planets in this configuration,the number of highly eccentric planets that have undergoneorbital migration is almost comparable for all kind of plan- T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]051015 i [ d e g ] Jupiter Saturn Uranus Neptune
Figure 7.
Same as in Fig. 6 but for a planetary system withinitial eccentricities of e = . from the 8k cluster. ets. The numbers of almost unexcited planets is only slightlylarger for the inner planets due to their smaller semimajoraxis. Figure 8 shows the orbital elements of a planetary sys-tem in which the planets are mostly unaffected for the firstfew million years despite several encounters. At t = . Myra red dwarf with a mass of 0.3 M (cid:12) approaches the systemcloser than 240 au causing a transfer of energy and angularmomentum from Neptune to the perturber. Due to the in-wards migration of Neptune on an orbit with an eccentricityof around 0.5, a fatal chain reaction with strong planet–planet scattering is triggered in which the eccentricity ofthe three other planets is excited as well. After a very shortchange of position with Saturn, Jupiter migrates inwardsand reaches twice an eccentricity of more than 0.9 beforeit is finally ejected at t = . Myr. During that time and in
MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]0255075 i [ d e g ] Jupiter Saturn Uranus Neptune
Figure 8.
Same as in Fig. 6 but for a massive planetary systemfrom the 8k cluster. the following 7 Myr the remaining planets Saturn, Uranus,and Neptune change their order several times. Due to thatstrong interaction, Uranus migrates outwards to a very wideand eccentric orbit. Saturn follows at t = . Myr which fi-nally leads to the ejection of Uranus at t = . Myr. For theremaining 85 Myr, Saturn retains its wide orbit of a ≈ auwhile Neptune remains at a = . au. a - i space A change in eccentricity is often directly related to a changein inclination since both result from the transfer of angularmomentum. By looking on the a - i space of the planets af-ter 100 Myr for the different cluster sizes in Figs 9, A4, A5,and A6, we can again see a wide distribution in that parame-ter space, although all planets started on coplanar, progradeorbits (Figs A7, A8, A9, and A10 show the e - i space after100 Myr for comparison). Those planets that get excited topolar orbits ( i ≈ ◦ ) or retrograde orbits ( i > ◦ ) are ofspecial interest.In 2006, Remijan & Hollis (2006) found first evidencethat parts of the protoplanetary disc around the binary sys-tem IRAS 16293–2422 are counterrotating which means thatplanets that form in that region would have a retrograde or-bit. The first two detected planets for which a retrograde orpolar orbit is assumed are WASP-17b (Anderson et al. 2010)and HAT-P-7b (Winn et al. 2009).We find planets with inclined orbits of more than ◦ in all of our simulations, independent of the initial configu-ration and stellar density of the host cluster. Some planetsswitch to a retrograde orbit for only a few million years butsome also keep their highly inclined orbit for the rest of thesimulation. An example of the latter case is shown in Fig. 10.At 71 Myr, the encounter of a 1.9 M (cid:12) star of less than 600 aucauses Uranus (the outermost planet at that time) to switchfrom a prograde orbit (inclined by ◦ ) to a retrograde orbitwith i = ◦ . Despite several additional encounters duringthe remaining 29 Myr with periastron distances of less than1000 au, Uranus keeps its retrograde orbit until the end ofthe simulation.There is no clear trend visible in which configuration or with which cluster density we can expect the highest fractionof retrograde orbits. However, in all four cluster simulations,the massive configuration results in the largest number ofhighly inclined orbits with i > ◦ after 100 Myr. We cantherefore conclude that retrograde orbits mainly occur dueto strong external perturbation while planet–planet scatter-ing especially seems to be an additional source for the exci-tation of planetary orbits to the range of i ≈ ◦ − ◦ . In the previous sections, we have shown the differences be-tween the initial configurations averaged over identical 200planetary systems. However, from this we can only have arough estimate of how the dynamical evolution of one plan-etary system looks like if we put it in different initial con-figurations. Therefore, we show the dynamical evolution ofplanetary system t = Myr,Uranus migrates inwards by ∼ . au and Neptune outwardsby ∼ au. Neptune’s increase in eccentricity and inclinationis the largest of all four planets.There is almost no difference in the dynamical evolu-tion of Jupiter, Saturn, and Uranus between the standardand compact configuration. However, instead of migratinginwards, Uranus keeps its semimajor axis in the compactconfiguration unlike Neptune which now migrates outwardsby 3.5 au due to the same encounter as in the standard con-figuration. Neptune is also more excited in eccentricity butless in inclination.In the resonant configuration, Saturn, Uranus, and Nep-tune do not even survive until the first very close encounterat t = Myr which is the only encounter in the standardand compact configuration which affects the system signif-icantly. Uranus and Neptune are both ejected within thefirst 1.3 Myr, whereas Saturn migrates to a = . au andstrongly oscillates in eccentricity within the following mil-lions of years. Jupiter, which migrates outwards to a semima-jor axis of 7.8 au, and Saturn cross their orbits after an addi-tional stellar encounter at t = . Myr which causes the ejec-tion of Saturn at t = . Myr. Due to that interaction, Jupitermigrates back to a smaller orbit of a = . au and stays com-pletely unaffected during additional encounters within therest of the simulation.The dynamical evolution of the first eccentric configu-ration is characterized by the interaction between Uranusand Neptune. The same encounter that affected the stan-dard and compact configuration causes a first close orbitalapproach of Uranus and Neptune after roughly 10 Myr. Dueto that, their initially small eccentricities increase as well astheir inclinations. Additional encounters, especially duringthe last 50 Myr of our simulation cause steady interactionand switch of orbital positions between the two outermostplanets. Jupiter and Saturn stay relatively unaffected in thissimulation.The higher initial eccentricity in the second eccentricconfiguration leads to a quicker interaction between Uranusand Neptune than in the previous case. After a switch of MNRAS000
Same as in Fig. 6 but for a massive planetary systemfrom the 8k cluster. the following 7 Myr the remaining planets Saturn, Uranus,and Neptune change their order several times. Due to thatstrong interaction, Uranus migrates outwards to a very wideand eccentric orbit. Saturn follows at t = . Myr which fi-nally leads to the ejection of Uranus at t = . Myr. For theremaining 85 Myr, Saturn retains its wide orbit of a ≈ auwhile Neptune remains at a = . au. a - i space A change in eccentricity is often directly related to a changein inclination since both result from the transfer of angularmomentum. By looking on the a - i space of the planets af-ter 100 Myr for the different cluster sizes in Figs 9, A4, A5,and A6, we can again see a wide distribution in that parame-ter space, although all planets started on coplanar, progradeorbits (Figs A7, A8, A9, and A10 show the e - i space after100 Myr for comparison). Those planets that get excited topolar orbits ( i ≈ ◦ ) or retrograde orbits ( i > ◦ ) are ofspecial interest.In 2006, Remijan & Hollis (2006) found first evidencethat parts of the protoplanetary disc around the binary sys-tem IRAS 16293–2422 are counterrotating which means thatplanets that form in that region would have a retrograde or-bit. The first two detected planets for which a retrograde orpolar orbit is assumed are WASP-17b (Anderson et al. 2010)and HAT-P-7b (Winn et al. 2009).We find planets with inclined orbits of more than ◦ in all of our simulations, independent of the initial configu-ration and stellar density of the host cluster. Some planetsswitch to a retrograde orbit for only a few million years butsome also keep their highly inclined orbit for the rest of thesimulation. An example of the latter case is shown in Fig. 10.At 71 Myr, the encounter of a 1.9 M (cid:12) star of less than 600 aucauses Uranus (the outermost planet at that time) to switchfrom a prograde orbit (inclined by ◦ ) to a retrograde orbitwith i = ◦ . Despite several additional encounters duringthe remaining 29 Myr with periastron distances of less than1000 au, Uranus keeps its retrograde orbit until the end ofthe simulation.There is no clear trend visible in which configuration or with which cluster density we can expect the highest fractionof retrograde orbits. However, in all four cluster simulations,the massive configuration results in the largest number ofhighly inclined orbits with i > ◦ after 100 Myr. We cantherefore conclude that retrograde orbits mainly occur dueto strong external perturbation while planet–planet scatter-ing especially seems to be an additional source for the exci-tation of planetary orbits to the range of i ≈ ◦ − ◦ . In the previous sections, we have shown the differences be-tween the initial configurations averaged over identical 200planetary systems. However, from this we can only have arough estimate of how the dynamical evolution of one plan-etary system looks like if we put it in different initial con-figurations. Therefore, we show the dynamical evolution ofplanetary system t = Myr,Uranus migrates inwards by ∼ . au and Neptune outwardsby ∼ au. Neptune’s increase in eccentricity and inclinationis the largest of all four planets.There is almost no difference in the dynamical evolu-tion of Jupiter, Saturn, and Uranus between the standardand compact configuration. However, instead of migratinginwards, Uranus keeps its semimajor axis in the compactconfiguration unlike Neptune which now migrates outwardsby 3.5 au due to the same encounter as in the standard con-figuration. Neptune is also more excited in eccentricity butless in inclination.In the resonant configuration, Saturn, Uranus, and Nep-tune do not even survive until the first very close encounterat t = Myr which is the only encounter in the standardand compact configuration which affects the system signif-icantly. Uranus and Neptune are both ejected within thefirst 1.3 Myr, whereas Saturn migrates to a = . au andstrongly oscillates in eccentricity within the following mil-lions of years. Jupiter, which migrates outwards to a semima-jor axis of 7.8 au, and Saturn cross their orbits after an addi-tional stellar encounter at t = . Myr which causes the ejec-tion of Saturn at t = . Myr. Due to that interaction, Jupitermigrates back to a smaller orbit of a = . au and stays com-pletely unaffected during additional encounters within therest of the simulation.The dynamical evolution of the first eccentric configu-ration is characterized by the interaction between Uranusand Neptune. The same encounter that affected the stan-dard and compact configuration causes a first close orbitalapproach of Uranus and Neptune after roughly 10 Myr. Dueto that, their initially small eccentricities increase as well astheir inclinations. Additional encounters, especially duringthe last 50 Myr of our simulation cause steady interactionand switch of orbital positions between the two outermostplanets. Jupiter and Saturn stay relatively unaffected in thissimulation.The higher initial eccentricity in the second eccentricconfiguration leads to a quicker interaction between Uranusand Neptune than in the previous case. After a switch of MNRAS000 , 1–17 (2020) K. Stock et al. Eccentric2 10 Massive a [AU] i [ d e g ] JupiterSaturnUranusNeptune
Figure 9.
The a - i space for all planets in the 8k star cluster which are not ejected from their host planetary system after a simulationtime of 100 Myr for the six different initial configurations. The dotted black line shows the threshold of i = ◦ . Planets near that valuehave polar orbits while those above it have retrograde orbits. T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]050100150 i [ d e g ] Jupiter Saturn Uranus Neptune
Figure 10.
Same as in Fig. 6 but for a planetary system withinitial eccentricities of e = . from the 32k cluster. The time isplotted in linear scale. positions the whole planetary system stays stable for therest of the simulation.The massive configuration reveals the increasing riskfor the innermost planets if all planets in the system havethe same mass. In all of the previous configurations, Jupitersurvived without facing any serious dangers for its orbitalstability while Saturn survived in four configurations. In themassive configuration, these two planets are the only planetsthat get ejected while Uranus and Neptune survive. The ex-ternal perturbation that leads to the ejection of Jupiter andSaturn is the same stellar encounter which is the formativeencounter in the dynamical evolution of the system in thestandard, compact and first eccentric configuration. We have explored the stability and vulnerability of 4800planetary systems, which are exposed to repeated stellar en-counters in the star cluster in which they formed. In each ofour four star clusters, we distribute 200 identical planetarysystems in six different initial configurations, which were in-spired by the Monte Carlo simulations of Li & Adams (2015).All planetary systems were Solar system analogues (withhost star masses of ∼ M (cid:12) ) consisting of the solar system’sgas giant planets Jupiter, Saturn, Uranus, and Neptune. Inthe standard configuration, the planets have their currentsemimajor axes but circular orbits. Two other configurationsare more compact versions of that case which are called thecompact and resonant configurations (due to mutual MMRsbetween the planets). In two additional configurations, theeccentricities of the standard case are increased to the plan-ets’ present-day values and to larger values of e = . (thefirst and second eccentric configurations). The sixth config-uration differed from the standard configuration only in theequal planetary masses of one Jovian mass.Our results for the cluster simulation can be summa-rized as follows: after 100 Myr the star clusters have under-gone the phases of mass segregation, stellar mass-loss andcore collapse, and re-expand again after a time of maxi-mum central density. The maximum central density reachesabout 10 times the initial central density; after mass segre-gation and core collapse the cluster generally re-expands andreaches a quasi-stationary state, where the central densityis about equal to its initial value, and the average densityinside the half-mass radius has dropped by a factor of ap-proximately 10. We observe that at that stage most of thedynamical interactions between planetary systems and pass-ing stars are over, so for the current pilot study we stop ourmodels at 100 Myr.Generally, the most stable planetary systems are the MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]0.02.55.07.5 i [ d e g ] Jupiter Saturn Uranus Neptune 0 20 40 60 80 100 T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]012 i [ d e g ] Jupiter Saturn Uranus Neptune0 20 40 60 80 100 T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]050100 i [ d e g ] Jupiter Saturn Uranus Neptune 0 20 40 60 80 100 T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]0510 i [ d e g ] Jupiter Saturn Uranus Neptune0 20 40 60 80 100 T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]051015 i [ d e g ] Jupiter Saturn Uranus Neptune 0 20 40 60 80 100 T [Myr]10 a [ A U ] R S C [ p c ] e l o g ( r p ) [ A U ] T [Myr]0204060 i [ d e g ] Jupiter Saturn Uranus Neptune
Figure 11.
Comparison of the dynamical evolution of one planetary system in the 32k cluster around the same host star in all sixdifferent initial configurations. Top left: Standard configuration.
Top right:
Compact configuration.
Middle left:
Resonant configuration.
Middle right:
Eccentric
Bottom left:
Eccentric
Bottom right:
Massive configuration.MNRAS000
Massive configuration.MNRAS000 , 1–17 (2020) K. Stock et al. standard and compact ones, and the configuration withsmall (current) eccentricities. The results for the standard,compact, and first eccentric configuration are comparablein fractions of surviving planets and final distribution in a − e and a − i space. However, a trend is observable thatthe compact system becomes slightly more resistant and theeccentric one slightly more vulnerable with increasing stellardensity relative to the standard case.We note that the compact system relative to the stan-dard system shows very little differences — one would expectthat it experiences less strong interactions under the effectof the same encounters as the standard system; our resultof very similar survival fractions can only be explained bystronger internal interactions, which destabilize the systemeven after relative weak perturbations. Furthermore, smallinitial eccentricities seem to not significantly change the vul-nerability of a planetary system.Due to its innermost position and highest mass, Jupiteris generally the planet with the highest chance to survive aperturbation by a stellar encounter of another cluster mem-ber, followed by Saturn. The exact order of the survival frac-tion of the two outermost planets Uranus and Neptune de-pends on the initial configuration and cluster density. How-ever, usually Uranus is slightly more likely to be ejectedfrom the planetary system due to an encounter or secularevolution. Even though Uranus is not the outermost planet,its lower mass makes the planet slightly more vulnerable togravitational perturbations from the host cluster due to itslower gravitational binding energy compared to Neptune.This difference can especially be seen in the survival frac-tions for the second eccentric configurations. In all four clus-ters, Uranus has by far the lowest chance to survive in thesystem if all planets are started with their true semimajoraxes but with an eccentricity of 0.1. If all planets have equalmasses, the differences in survival fractions shrink signifi-cantly. Due to its smallest semimajor axis, Jupiter still has aslightly higher chance for survival while the rates for Saturn,Uranus, and Neptune are almost equal. From this, we candeduce that a planet’s mass (compared to the other planetsin the system) plays a more crucial role for the estimationof its vulnerability than its semi-major axis.The fourth most stable system is the massive configu-ration in the 8k and 16k cluster but the system with initialeccentricities of e = . is instead more resistant in the 32kand 64k cluster. In all clusters, the resonant system is theone with the highest vulnerability. However, the system is aspecial case and very interesting for a certain reason. Our in-tegrations show that without perturbations by passing starsit is generally very short lived, getting unstable after about yr, around that time Uranus and Neptune are inevitablyejected from the planetary system. However, embedded in astar cluster, the system tends to be more stable. We believethat this is due to a process where stellar encounters detuneor break the resonances and thus render the systems morestable. In many of the simulated systems, this is achieved byonly ejecting one of the outer planets (Uranus or Neptune),and then the remaining three-planet system survives muchlonger than in the isolated case.In van Elteren et al. (2019), the authors find that theprobability of a star to lose a planet is independent of theplanet mass and independent of its initial orbital separation.As a consequence, the mass distribution of free-floating plan- ets would be indistinguishable from the mass distribution ofplanets bound to their host star. Our results do not confirmthis. The discrepancy may result from the larger numberof stars in the clusters in our simulations, the longer evo-lutionary time-scales (we integrated for 100 Myr whereas invan Elteren et al. (2019) they integrated up to 10 Myr), andfinally they adopted the Oligarchic growth model for plan-etary systems. In the latter model, planet mass increasesfurther away from the host star. This has interesting con-sequences for the stability of the planetary systems fromperturbations from inside as well as for external perturba-tions. A small perturbation from another star may render anentire planetary system catastrophically unstable, whereasif the outer most planets have low mass, such a system sur-vives more easily in a dense stellar environment.The survival fractions for the different planet types inour simulations are generally smaller than those of Li &Adams (2015). This is due to the different approaches. First,Li & Adams (2015) randomly select their encounter param-eter equally from the available phase space that is not real-istic (see figs 1 and 2 in Spurzem et al. 2009). Secondly, Li &Adams (2015) only focus on the prompt ejections of planetswhile we continue the integration of the planetary systemslong enough to account for secular evolution. Thirdly, ourplanetary systems are exposed to the cumulative effect ofseveral encounters over a significant fraction of the host starcluster’s lifetime. From the reduced survivability of the plan-ets, which we see in our results compared to Li & Adams(2015), we can conclude that the effects of secular evolutionand cumulative encounters are not negligible.We find that passing stars excite mutual inclinationsbetween planets in our planetary systems; quite some caseslead to high values of relative inclination and even tocounter-rotating planets. It is quite impossible to excite sig-nificant inclinations by internal evolution of planetary sys-tems, they are a tell-tale sign of the important role of stellarencounters in shaping the planetary system. While this effecthas been mentioned in previous studies (such as in Spurzemet al. 2009), there is not yet a more quantitative study ofthis process.Our simulations could be and will be refined in futurework in many ways. Planetary systems around more mas-sive stars are subject to orbital changes when the host starbecomes a red giant and finally loses significant mass. Thepresence of many initial binaries, which is expected fromstar and cluster formation, will be an interesting issue —including S- and P-type planetary systems.The Monte Carlo models of Li & Adams (2015) givesome information about encounters between planetary sys-tems and binary stars. Finally, in this work we have onlypresented a limited set of star clusters. A wider parameterstudy may be required to predict the impact of stellar en-counters on the final planetary population in the Galacticfield. Other processes shaping planetary systems in the for-mation process inside a star cluster have also not been takeninto account here.We have, however, clearly shown that encounters ofpassing stars in star clusters have a considerable effect andcontribute to the diversity of planetary systems in all re-spects. MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs ACKNOWLEDGEMENTS
We want to thank Rosemary Mardling, Willy Kley, andGongjie Li for useful discussions. KS and RS acknowledgethe support of the DFG priority program SPP 1992 “Ex-ploring the Diversity of Extrasolar Planets” (SP 345/20-1).RS has been supported by National Astronomical Observa-tories of Chinese Academy of Sciences, Silk Road Project,and by National Natural Science Foundation of China un-der grant no. 11673032. All simulations have been doneon the GPU accelerated cluster “kepler”, funded by Volk-swagen Foundation grants 84678/84680. MBNK was sup-ported by the National Natural Science Foundation of China(grant 11573004) and by the Research Development Fund(grant RDF-16-01-16) of Xi’an Jiaotong-Liverpool Univer-sity (XJTLU).
DATA AVAILABILITY
The data underlying this article will be shared on reasonablerequest to the corresponding author.
REFERENCES
Aarseth S. J., 1999, CeMDA, 73, 127Anderson D. R., et al., 2010, ApJ, 709, 159Anderson K. R., Adams F. C., Calvet N., 2013, ApJ, 774, 9Armitage P. J., 2000, A&A, 362, 968Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R.,Spurzem R., Wang L., 2020, A&A, 639, A41Boekholt T. C. N., Portegies Zwart S. F., Valtonen M., 2020,MNRAS, 493, 3932Brucalassi A., et al., 2014, A&A, 561, L9Brucalassi A., et al., 2016, A&A, 592, L1Brucalassi A., et al., 2017, A&A, 603, A85Cai M. X., Meiron Y., Kouwenhoven M. B. N., Assmann P.,Spurzem R., 2015, ApJS, 219, 31Cai M. X., Kouwenhoven M. B. N., Portegies Zwart S. F.,Spurzem R., 2017, MNRAS, 470, 4337Cai M. X., Portegies Zwart S., van Elteren A., 2018, MNRAS,474, 5114Cai M. X., Portegies Zwart S., Kouwenhoven M. B. N., SpurzemR., 2019, MNRAS, 489, 4311Clarke C. J., Pringle J. E., 1993, MNRAS, 261, 190Concha-Ram´ırez F., Vaher E., Portegies Zwart S., 2019, MNRAS,482, 732Delgado Mena E., et al., 2018, A&A, 619, A2Ernst A., Just A., Spurzem R., Porth O., 2008, MNRAS, 383, 897van Elteren A., Portegies Zwart S., Pelupessy I., Cai M. X.,McMillan S. L. W., 2019, A&A, 624, A120Faber N. T., Stibbe D., Portegies Zwart S., McMillan S. L. W.,Boily C. M., 2010, MNRAS, 401, 1898Facchini S., Clarke C. J., Bisbas T. G., 2016, MNRAS, 457, 3593Farr W. M., et al., 2012, NewA, 17, 520Flammini Dotti F., Kouwenhoven M. B. N., Cai M. X., SpurzemR., 2019, MNRAS, 489, 2280Fujii M. S., Hori Y., 2019, A&A, 624, A110Gaidos E., et al., 2017, MNRAS, 464, 850Glaser J. P., McMillan S. L. W., Geller A. M., Thornton J. D.,Giovinazzi M. R., 2020, arXiv, arXiv:2002.12375Gomes R., Levison H. F., Tsiganis K., Morbidelli A., 2005, Natur,435, 466Hands T. O., Dehnen W., Gration A., Stadel J., Moore B., 2019,MNRAS, 490, 21 Hao W., Kouwenhoven M. B. N., Spurzem R., 2013, MNRAS,433, 867Heisler J., Tremaine S., 1986, Icar, 65, 13Hurley J. R., Pols O. R., Aarseth S. J., Tout C. A., 2005, MNRAS,363, 293J´ılkov´a L., Hamers A. S., Hammer M., Portegies Zwart S., 2016,MNRAS, 457, 4218Just A., Berczik P., Petrov M. I., Ernst A., 2009, MNRAS, 392,969Khalaj P., Baumgardt H., 2015, MNRAS, 452, 924King I., 1962, AJ, 67, 471King I. R., 1966, AJ, 71, 64Kokubo E., Ida S., 1998, Icar, 131, 171Kouwenhoven M. B. N., Goodwin S. P., Parker R. J., DaviesM. B., Malmberg D., Kroupa P., 2010, MNRAS, 404, 1835Kroupa P., 2001, MNRAS, 322, 231Lada E. A., Strom K. M., Myers P. C., 1993, Protostars andPlanets III, Univ. Arizona Press, Tucson, AZ, p. 245Li G., Adams F. C., 2015, MNRAS, 448, 344Livingston J. H., et al., 2019, MNRAS, 484, 8Malavolta L., et al., 2016, A&A, 588, A118Malmberg D., de Angeli F., Davies M. B., Church R. P., MackeyD., Wilkinson M. I., 2007, MNRAS, 378, 1207Malmberg D., Davies M. B., Heggie D. C., 2011, MNRAS, 411,859Mann A. W., et al., 2016, ApJ, 818, 46Mann A. W., et al., 2017, AJ, 153, 64Mann A. W., et al., 2018, AJ, 155, 4Miller R. H., 1964, ApJ, 140, 250Obermeier C., et al., 2016, AJ, 152, 223Olczak C., Pfalzner S., Spurzem R., 2006, ApJ, 642, 1140’Oumuamua ISSI Team, et al., 2019, Nat. Astron., 3, 594Parker R. J., Quanz S. P., 2012, MNRAS, 419, 2448Perets H. B., Kouwenhoven M. B. N., 2012, ApJ, 750, 83Pfalzner S., Bannister M. T., 2019, ApJL, 874, L34Plummer H. C., 1911, MNRAS, 71, 460Portegies Zwart S., 2011, Astrophysics Source Code Library,record ascl:1107.007Portegies Zwart S. F., J´ılkov´a L., 2015, MNRAS, 451, 144Portegies Zwart S. F., 2016, MNRAS, 457, 313Portegies Zwart S., Torres S., Pelupessy I., B´edorf J., Cai M. X.,2018, MNRAS, 479, L17Portegies Zwart S., McMillan S., 2018, Astrophysical Recipes;The art of AMUSE, IOP ebooks, Bristol, UKQian P. X., Cai M. X., Portegies Zwart S., Zhu M., 2017, PASP,129, 094503Quinlan G. D., Tremaine S., 1992, MNRAS, 259, 505Quinn S. N., et al., 2012, ApJL, 756, L33Quinn S. N., et al., 2014, ApJ, 787, 27Rein H., Liu S.-F., 2012, A&A, 537, A128Rein H., Spiegel D. S., 2015, MNRAS, 446, 1424Remijan A. J., Hollis J. M., 2006, ApJ, 640, 842Rizzuto A. C., et al., 2018, AJ, 156, 195Sato B., et al., 2007, ApJ, 661, 527Shu Q., Kouwenhoven M.B.N., Flammini Dotti F., Hong J.,Spurzem R., submitted to MNRASSpera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086Spurzem R., Giersz M., Heggie D. C., Lin D. N. C., 2009, ApJ,697, 458St¨orzer H., Hollenbach D., 1999, ApJ, 515, 669Torres S., Cai M. X., Brown A. G. A., Portegies Zwart S., 2019,A&A, 629, A139Veras D., Reichert K., Flammini Dotti F., Cai M. X., PortegiesZwart S., Kouwenhoven M.B.N., McDonald C. H., ShannonA. B., Mustill A. J., 2020, MNRAS, 493, 5062Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P.,Kouwenhoven M. B. N., Naab T., 2015, MNRAS, 450, 4070MNRAS000
Aarseth S. J., 1999, CeMDA, 73, 127Anderson D. R., et al., 2010, ApJ, 709, 159Anderson K. R., Adams F. C., Calvet N., 2013, ApJ, 774, 9Armitage P. J., 2000, A&A, 362, 968Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R.,Spurzem R., Wang L., 2020, A&A, 639, A41Boekholt T. C. N., Portegies Zwart S. F., Valtonen M., 2020,MNRAS, 493, 3932Brucalassi A., et al., 2014, A&A, 561, L9Brucalassi A., et al., 2016, A&A, 592, L1Brucalassi A., et al., 2017, A&A, 603, A85Cai M. X., Meiron Y., Kouwenhoven M. B. N., Assmann P.,Spurzem R., 2015, ApJS, 219, 31Cai M. X., Kouwenhoven M. B. N., Portegies Zwart S. F.,Spurzem R., 2017, MNRAS, 470, 4337Cai M. X., Portegies Zwart S., van Elteren A., 2018, MNRAS,474, 5114Cai M. X., Portegies Zwart S., Kouwenhoven M. B. N., SpurzemR., 2019, MNRAS, 489, 4311Clarke C. J., Pringle J. E., 1993, MNRAS, 261, 190Concha-Ram´ırez F., Vaher E., Portegies Zwart S., 2019, MNRAS,482, 732Delgado Mena E., et al., 2018, A&A, 619, A2Ernst A., Just A., Spurzem R., Porth O., 2008, MNRAS, 383, 897van Elteren A., Portegies Zwart S., Pelupessy I., Cai M. X.,McMillan S. L. W., 2019, A&A, 624, A120Faber N. T., Stibbe D., Portegies Zwart S., McMillan S. L. W.,Boily C. M., 2010, MNRAS, 401, 1898Facchini S., Clarke C. J., Bisbas T. G., 2016, MNRAS, 457, 3593Farr W. M., et al., 2012, NewA, 17, 520Flammini Dotti F., Kouwenhoven M. B. N., Cai M. X., SpurzemR., 2019, MNRAS, 489, 2280Fujii M. S., Hori Y., 2019, A&A, 624, A110Gaidos E., et al., 2017, MNRAS, 464, 850Glaser J. P., McMillan S. L. W., Geller A. M., Thornton J. D.,Giovinazzi M. R., 2020, arXiv, arXiv:2002.12375Gomes R., Levison H. F., Tsiganis K., Morbidelli A., 2005, Natur,435, 466Hands T. O., Dehnen W., Gration A., Stadel J., Moore B., 2019,MNRAS, 490, 21 Hao W., Kouwenhoven M. B. N., Spurzem R., 2013, MNRAS,433, 867Heisler J., Tremaine S., 1986, Icar, 65, 13Hurley J. R., Pols O. R., Aarseth S. J., Tout C. A., 2005, MNRAS,363, 293J´ılkov´a L., Hamers A. S., Hammer M., Portegies Zwart S., 2016,MNRAS, 457, 4218Just A., Berczik P., Petrov M. I., Ernst A., 2009, MNRAS, 392,969Khalaj P., Baumgardt H., 2015, MNRAS, 452, 924King I., 1962, AJ, 67, 471King I. R., 1966, AJ, 71, 64Kokubo E., Ida S., 1998, Icar, 131, 171Kouwenhoven M. B. N., Goodwin S. P., Parker R. J., DaviesM. B., Malmberg D., Kroupa P., 2010, MNRAS, 404, 1835Kroupa P., 2001, MNRAS, 322, 231Lada E. A., Strom K. M., Myers P. C., 1993, Protostars andPlanets III, Univ. Arizona Press, Tucson, AZ, p. 245Li G., Adams F. C., 2015, MNRAS, 448, 344Livingston J. H., et al., 2019, MNRAS, 484, 8Malavolta L., et al., 2016, A&A, 588, A118Malmberg D., de Angeli F., Davies M. B., Church R. P., MackeyD., Wilkinson M. I., 2007, MNRAS, 378, 1207Malmberg D., Davies M. B., Heggie D. C., 2011, MNRAS, 411,859Mann A. W., et al., 2016, ApJ, 818, 46Mann A. W., et al., 2017, AJ, 153, 64Mann A. W., et al., 2018, AJ, 155, 4Miller R. H., 1964, ApJ, 140, 250Obermeier C., et al., 2016, AJ, 152, 223Olczak C., Pfalzner S., Spurzem R., 2006, ApJ, 642, 1140’Oumuamua ISSI Team, et al., 2019, Nat. Astron., 3, 594Parker R. J., Quanz S. P., 2012, MNRAS, 419, 2448Perets H. B., Kouwenhoven M. B. N., 2012, ApJ, 750, 83Pfalzner S., Bannister M. T., 2019, ApJL, 874, L34Plummer H. C., 1911, MNRAS, 71, 460Portegies Zwart S., 2011, Astrophysics Source Code Library,record ascl:1107.007Portegies Zwart S. F., J´ılkov´a L., 2015, MNRAS, 451, 144Portegies Zwart S. F., 2016, MNRAS, 457, 313Portegies Zwart S., Torres S., Pelupessy I., B´edorf J., Cai M. X.,2018, MNRAS, 479, L17Portegies Zwart S., McMillan S., 2018, Astrophysical Recipes;The art of AMUSE, IOP ebooks, Bristol, UKQian P. X., Cai M. X., Portegies Zwart S., Zhu M., 2017, PASP,129, 094503Quinlan G. D., Tremaine S., 1992, MNRAS, 259, 505Quinn S. N., et al., 2012, ApJL, 756, L33Quinn S. N., et al., 2014, ApJ, 787, 27Rein H., Liu S.-F., 2012, A&A, 537, A128Rein H., Spiegel D. S., 2015, MNRAS, 446, 1424Remijan A. J., Hollis J. M., 2006, ApJ, 640, 842Rizzuto A. C., et al., 2018, AJ, 156, 195Sato B., et al., 2007, ApJ, 661, 527Shu Q., Kouwenhoven M.B.N., Flammini Dotti F., Hong J.,Spurzem R., submitted to MNRASSpera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086Spurzem R., Giersz M., Heggie D. C., Lin D. N. C., 2009, ApJ,697, 458St¨orzer H., Hollenbach D., 1999, ApJ, 515, 669Torres S., Cai M. X., Brown A. G. A., Portegies Zwart S., 2019,A&A, 629, A139Veras D., Reichert K., Flammini Dotti F., Cai M. X., PortegiesZwart S., Kouwenhoven M.B.N., McDonald C. H., ShannonA. B., Mustill A. J., 2020, MNRAS, 493, 5062Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P.,Kouwenhoven M. B. N., Naab T., 2015, MNRAS, 450, 4070MNRAS000 , 1–17 (2020) K. Stock et al.
Wang L., Kouwenhoven M. B. N., Zheng X., Church R. P., DaviesM. B., 2015, MNRAS, 449, 3543Wang L., et al., 2016, MNRAS, 458, 1450Winn J. N., Johnson J. A., Albrecht S., Howard A. W., MarcyG. W., Crossfield I. J., Holman M. J., 2009, ApJL, 703, L99Zheng X., Kouwenhoven M. B. N., Wang L., 2015, MNRAS, 453,2759 MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs Eccentric2 10 Massive a [AU] e JupiterSaturnUranusNeptune
Figure A1.
Same as Fig. 5 but for the 16k cluster. Eccentric2 10 Massive a [AU] e JupiterSaturnUranusNeptune
Figure A2.
Same as Fig. 5 but for the 32k cluster.
APPENDIX A: ADDITIONAL MATERIAL
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000
This paper has been typeset from a TEX/L A TEX file prepared by the author.MNRAS000 , 1–17 (2020) K. Stock et al. Eccentric2 10 Massive a [AU] e JupiterSaturnUranusNeptune
Figure A3.
Same as Fig. 5 but for the 64k cluster. Eccentric2 10 Massive a [AU] i [ d e g ] JupiterSaturnUranusNeptune
Figure A4.
Same as Fig. 9 but for the 16k cluster. MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs Eccentric2 10 Massive a [AU] i [ d e g ] JupiterSaturnUranusNeptune
Figure A5.
Same as Fig. 5 but for the 32k cluster. Eccentric2 10 Massive a [AU] i [ d e g ] JupiterSaturnUranusNeptune
Figure A6.
Same as Fig. 5 but for the 64k cluster.MNRAS000
Same as Fig. 5 but for the 64k cluster.MNRAS000 , 1–17 (2020) K. Stock et al. e i [ d e g ] JupiterSaturnUranusNeptune
Figure A7.
The e - i space for all planets in the 8k star cluster which are not ejected from their host planetary system after a simulationtime of 100 Myr for the six different initial configurations. The dotted black shows the threshold of i = ◦ . Planets near that value havepolar orbits while those above it have retrograde orbits. e i [ d e g ] JupiterSaturnUranusNeptune
Figure A8.
Same as in Fig. A7 but for the 16k cluster. MNRAS , 1–17 (2020) esonant and non-resonant planetary systems in SCs e i [ d e g ] JupiterSaturnUranusNeptune
Figure A9.
Same as in Fig. A7 but for the 32k cluster. e i [ d e g ] JupiterSaturnUranusNeptune
Figure A10.
Same as in Fig. A7 but for the 64k cluster.MNRAS000