Multiplanetary Systems in Star Clusters

Stability of Multiplanetary Systems in Star Clusters

Maxwell Xu Cai ({CJK*}UTF8gbsn蔡栩), M.B.N. Kouwenhoven, Simon F. Portegies Zwart, and Rainer Spurzem
Leiden Observatory, Leiden University, PO Box 9513, 2300 RA, Leiden, The Netherlands
Department of Mathematical Sciences, Xi’an Jiaotong-Liverpool University (XJTLU), 111 Ren’ai Road, Dushu Lake Higher Education
Town, Suzhou Industrial Park, Suzhou 215123, P.R. China
National Astronomical Observatories and Key Laboratory of Computational Astrophysics, Chinese Academy of Sciences, 20A Datun Road,
Chaoyang District, Beijing 100012, P.R. China
Kavli Institute for Astronomy and Astrophysics, Peking University, 5 Yi He Yuan Road, Haidian District, Beijing 100871, P.R. China
Astronomisches Rechen-Institut, Zentrum für Astronomie, University of Heidelberg, Mönchhofstrasse 12-14, 69120 Heidelberg, Germany
E-mail: (MXC)
Accepted 2017 June 9. Received 2017 June 7; in original form 2016 September 27

Most stars form in star clusters and stellar associated. However, only about of the presently known exoplanets are found in these environments. To understand the roles of star cluster environments in shaping the dynamical evolution of planetary systems, we carry out direct -body simulations of four planetary systems models in three different star cluster environments with respectively and stars. In each cluster, an ensemble of initially identical planetary systems are assigned to solar-type stars with and evolved for 50 Myr. We found that following the depletion of protoplanetary disks, external perturbations and planet-planet interactions are two driving mechanisms responsible for the destabilization of planetary systems. The planet survival rate varies from in the cluster to in the cluster, which suggests that most planetary systems can indeed survive in low-mass clusters, except in the central regions. We also find that planet ejections through stellar encounters are cumulative processes, as only of encounters are strong enough to excite the eccentricity by . Short-period planets can be perturbed through orbit crossings with long-period planets. When taking into account planet-planet interactions, the planet ejection rate nearly doubles, and therefore multiplicity contributes to the vulnerability of planetary systems. In each ensemble, of planetary orbits become retrograde due to random directions of stellar encounters. Our results predict that young low-mass star clusters are promising sites for next-generation planet surveys, yet low planet detection rates are expected in dense globular clusters such as 47 Tuc. Nevertheless, planets in denser stellar environments are likely to have shorter orbital periods, which enhances their detectability.

planets and satellites: dynamical evolution and stability – galaxies: star clusters: general – methods: numerical
pubyear: 2017pagerange: Stability of Multiplanetary Systems in Star ClustersStability of Multiplanetary Systems in Star Clusters

1 Introduction

Research of planetary systems can be dated back to at least a few thousand years ago, when astronomers began to observe the night sky and named the planets in the Solar System “wanderers”. Observations of star clusters started several centuries ago, after telescopes became powerful enough to resolve these “clouds of stars”. Only in the recent decades the possible relationship between planetary systems and star clusters was gradually recognised. Thanks to the availability of dedicated observational facilities such as Kepler, more than 3 600 extrasolar planetary systems have now been identified, among which about 610 are multiplanetary systems111 Notably, approximately 20 exoplanets have been detected in star clusters (see Table 1). The discoveries of exoplanets in star clusters show that star clusters contain a variety of celestial bodies, and therefore it is important to further study these star clusters and investigate the planetary systems that they host.

Designation () (days) () DM Cluster Reference
YBP401 b 0.46 4.09 1.14 RV M67 Brucalassi et al. (2016)
Pr0211c 7.9 5 300 0.935 RV M44 Malavolta et al. (2016)
EPIC-210490365 b 3.49 0.29 TS Hyades Mann et al. (2016)
SAND 364 b 1.5 121.7 1.35 RV M67 Brucalassi et al. (2014)
YBP1194 b 0.34 6.96 1.01 RV M67 Brucalassi et al. (2014)
YBP1514 b 0.4 5.11 0.96 RV M67 Brucalassi et al. (2014)
HD 285507 b 0.92 6.08 0.73 RV Hyades Quinn et al. (2014)
Kepler-66 b Nep 17.82 1.04 TS NGC6811 Meibom et al. (2013)
Kepler-67 b Nep 15.73 0.87 TS NGC6811 Meibom et al. (2013)
Pr0201b 0.54 4.33 1.234 RV M44 Quinn et al. (2012)
Pr0211b 1.88 2.15 0.935 RV M44 Quinn et al. (2012)
eps Tau b 7.6 595 2.70 RV Hyades Sato et al. (2007)
NGC 2423 3 b 10.6 714 2.40 RV NGC2423 Lovis & Mayor (2007)
NGC 4349 No 127 b 19.8 678 3.90 RV NGC4349 Lovis & Mayor (2007)
PSR B1620-26 b 2.50 36 525 1.35 TM M4 Backer et al. (1993)
Table 1: List of exoplanet detections in star clusters. DM: detection method; TS: transit; RV: radial velocity; TM: timing; Nep: Neptune-sized; : Stellar mass in solar units ; : planet mass in Jupiter units ; : orbital period in days.

Modern planet formation theories, such as the core accretion scenario (e.g., Mizuno, 1980; Bodenheimer & Pollack, 1986) and the disk-fragmentation scenario (e.g., Boss, 1997), were inspired by the Nebular Hypothesis proposed independently by Immanuel Kant and Pierre-Simon Laplace. According to these theories, planets form in the protoplanetary disks, are essentially byproducts of the star formation process. Consequently, planets are may be common around main sequence stars. The recent discoveries of many exoplanets indeed have confirmed the high frequencies of exoplanets.

Moreover, observational studies suggest that most stars form in groups and/or clusters following the collapse of giant molecular clouds (e.g., Clarke et al., 2000; Lada & Lada, 2003; Porras et al., 2003). Naturally, one would expect that planets also form in star clusters. Indeed, isotope analysis of meteorites suggests that even our own Solar system may have formed in a star cluster (see, e.g., Adams & Laughlin, 2001; Portegies Zwart, 2009; Pfalzner, 2013; Parker et al., 2014; Pfalzner et al., 2015, and references therein). Direct imaging surveys have revealed protoplanetary disks in the Trapezium cluster, which is embedded in the Orion Nebula (McCaughrean & O’dell, 1996; Bally et al., 2000). These disks are likely the precursors of planetary systems. The orbits of outer solar system bodies such as Sedna may also be explained as a consequence of close encounters, (Adams, 2010; Jílková et al., 2015). If planetary systems are formed in star clusters, their dynamical evolution will then be influenced not only by internal processes (e.g. Voyatzis et al., 2013; Zhang et al., 2013; Li et al., 2015, 2016) but also by the dynamical evolution of their host star clusters. The orbital architectures of the planetary systems discovered in the Galactic field can be thus used to constrain the properties of the star clusters in which they were formed.

As of May 2017, only 1% of the known exoplanets are in star clusters (as listed in Table 1). This apparent low frequency of exoplanets in star clusters can partially be attributed to observational selection effects. Nevertheless, theoretical studies indicate that external perturbations due to stellar flybys play a role in the evolution of planetary systems (e.g., Heggie & Rasio, 1996; Laughlin & Adams, 1998; Davies & Sigurdsson, 2001; Bonnell et al., 2001). Naturally, an immediate question arises on how the dynamical evolution of planetary systems depends on the hosting stellar environments. Galaxies and star clusters are distinctively different stellar environments, in the sense that the they have very different relaxation timescales (Binney & Tremaine, 2008). The relaxation time in a typical galaxy of stars is well over a Hubble time, and therefore close encounters are unimportant (except for the galactic center); exoplanets around stars in the galactic field rarely experience external perturbation due to stellar encounters. In contrast, the corresponding relaxation time in open clusters is of the order of  Myr. Should there be any planetary systems in star clusters, they are likely to experience external perturbations. This motivates Portegies Zwart & Jílková (2015) to define the “parking zone” of planetary systems: a range of orbital separations around any star within which the planets native to that star may have been perturbed by close encounters of external perturbers, but are unlikely to be affected by external perturbations once the star becomes a member of galactic field population. A number of recent numerical studies have shown that frequent close encounters that can be destructive for planetary systems (e.g., Spurzem et al., 2009; Malmberg et al., 2011; Hao et al., 2013; Liu et al., 2013; Li & Adams, 2015; Zheng et al., 2015; Wang et al., 2015b; Shara et al., 2016). Therefore, it would be insightful to carry out an exploratory study to understand the role of dynamical stellar environment on the evolution of planetary systems.

Due to the chaotic nature of -body systems with , the research on analytical formulations of generic multiplanetary systems remains extremely challenging. We therefore carry out the investigation with gravitational direct -body simulations. However, given the enormous spatial and temporal differences between the dynamics of star clusters and planetary systems, numerical investigation of planetary systems in star clusters is non-trivial, as it leads to a highly hierarchical and chaotic many-body systems. Substantial progress was made over the years with studies of single-planet systems in star clusters. Using direct -body simulations, Hurley & Shara (2002) investigate the fate of free-floating planets, and find that planets in star clusters, after being liberated by stellar encounters, can remain bound in open clusters for a half-mass relaxation time scale of the star cluster. In another direct -body study, Spurzem et al. (2009) investigate the dynamics of single-planetary systems in star clusters using both direct -body simulations and hybrid Monte-Carlo simulations. Their results indicate that compact and nearly-circular orbits are generally not affected by distant stellar encounters. On the other hand, while short-period planets are difficult to perturb, close encounters can excite these to moderate eccentricities, which may in turn result to orbital decay due to tidal dissipation. Zheng et al. (2015) discuss the effect of the initial virial state and the presence of initial substructure, and find that single-planet systems wider than approximately 200 AU are mostly vulnerable by the time the cluster reaches an age of 50 Myr.

Studies of multiplanetary systems in star clusters have thus far been limited to Monte Carlo scattering experiments. Hao et al. (2013) conclude that multiplanetary systems are affected by both direct interaction with the encountering star and planet-planet scattering. The combined effects can account for the apparent low frequency of exoplanets in star clusters, not only for those on wide orbits that are directly affected by stellar encounters, but also planets close to the star which can disappear long after a stellar encounter has perturbed the outskirts of the planetary system. A more recent study by Li & Adams (2015) extends the simulations to a much wider parameter space with over two million scattering simulations. By varying the compactness of the target solar systems, the velocity dispersion of the host star cluster, stellar host masses, starting eccentricities of planet orbits, and single versus binary perturbers, they characterised the encounter cross sections as a function of stellar host mass, cluster velocity dispersion, semi-major axis, and final eccentricity, and predicted that the Solar System was formed within an open cluster with stars.

In the case of stars with two planetary companions, it is possible to employ a variety of three-body approaches. Mardling (2008) proposes an analytical criterion for determining the stability of arbitrary three-body hierarchies which makes use of the chaos theory concept of resonance overlap. In a follow-up paper, the three-body stability algorithm given in Mardling (2008) is used to determine the stability of an ensemble of mini solar systems with two Jupiter-mass planets in open cluster environments (Shara et al., 2016).

Star cluster environments not only affect the post-formation dynamical evolution of planetary systems, but also affect the planet formation process through the circumstellar disks (e.g., Thies et al., 2010, 2011; Anderson et al., 2013). Adams et al. (2004) points out that the intensive radiation from nearby OB stars can modify the mass-loss rate and evaporation timescales of exposed circumstellar disks, and eventually affect the planet-formation processes and and planet migration through disk-planet interaction. Olczak et al. (2012) carries out simulations to study the mass-loss process driven by intensive radiation at the Arches cluster environment. They find that stellar encounters destroy one-third of the circumstellar disks in the cluster core within the first 2.5 Myr of evolution, and after 6 Myr half of the core population becomes diskless. However they also point out that the disk destruction process ceases after roughly 1 Myr in sparser clusters due to significant cluster expansion (Olczak et al., 2010). A recent study by Portegies Zwart (2016) shows that close encounters result in the truncation of circumstellar disks. This mechanism can be used to reproduce the observed distribution of disk sizes in the Orion Trapezium cluster. Furthermore, he shows that a subvirial () and fractal () initial environment is indicated according to the observed disk size distribution.

In most previous numerical experiments, the star cluster dynamics and planetary system dynamics are integrated by a single code. For example, Spurzem et al. (2009) carried out direct -body simulation with the code NBODY6++ (Spurzem, 1999) by initializing single planets as KS-binary with their host stars. This is highly accurate thanks to the KS regularization (Kustaanheimo & Stiefel., 1965) technique, yet this is prohibitively expensive, even for moderately large systems, and inefficient in handling multiplanetary systems. The Monto-Carlo approach is computationally quite affordable, but the results depend on the assumed distribution of the perturber’s velocity and impact parameter. Finally, free-floating planets (FFP) are expected in star clusters following the ejections from the host stars (see, e.g., Kouwenhoven et al., 2016, and references therein).

The host star clusters themselves are also evolving, driven by internal mechanisms such as two-body relaxation (e.g., Chandrasekhar, 1942; Spitzer & Hart, 1971; Hénon, 1971; Takahashi & Portegies Zwart, 2000) and stellar evolution (e.g., Applegate, 1986; Portegies Zwart et al., 1998; Kouwenhoven et al., 2014), as well as external mechanisms such as the interaction with the Galactic tidal field (e.g., Spitzer, 1987; Cai et al., 2016) and the spiral arms (e.g., Gieles et al. 2007). In order to take all these effects into account, we carry out direct -body simulations with NBODY6++GPU (Wang et al., 2015a, 2016; Aarseth, 2003; Spurzem, 1999), and obtain the properties of all stellar encounters in these simulations. Subsequently, the perturbation data are loaded into the AMUSE framework (e.g., Portegies Zwart et al., 2009, 2013; McMillan et al., 2012; Pelupessy et al., 2013) and are sent to the planetary systems integrator rebound (Rein & Liu, 2012), where they are included in the modelling of the planetary systems.

This paper is organized as follows. The modeling of the perturbations and the implementation of the simulations are presented in Section 2. The initial conditions of the star clusters and planetary systems are described in Section 3. The results are presented and discussed in Section 4. Finally, the conclusions are presented in Section 5.

2 Modeling and Implementation

2.1 Modeling Perturbations

It is convenient to work in a Cartesian coordinate system centered at the host star of the planetary system under consideration, as illustrated in Figure 1. The tidal forces experienced by the planetary systems due to stellar encounters as perturbations. In this frame of reference, the acceleration experienced by a planet is


is the acceleration experienced by planet due to the presence of its host star and the other planets in the system:


where is the mass of the host star, are the masses of the other planets, is the position vector of the -th planet, , and is the gravitational constant.

In a star cluster with stars, suppose that is the position vector of the -th star with respect to the cluster center. The acceleration experienced by the host star due to the other stars in the star cluster is then


where is the relative position vector from the host star to the perturber , and is the mass of the perturber. Finally, the acceleration experienced by planet due to the other stars in the star cluster is calculated as


Consider a simplified scenario where the planetary system is perturbed by only one perturber of mass . The tidal force in the vicinity of the host star is


where is the unit vector of . Therefore, the tidal force is proportional to (for simplicity, hereafter we call it the dependence), which is a strong function of the distance between the targeted planetary system and the perturber. If , to the first order of , the acceleration experienced by planet due to the perturber is


In reality, planetary systems in a star cluster are perturbed simultaneously by member stars. Direct summation of the tidal forces due to these stars is possible but expensive. To explore whether it is possible to simplify the calculation by taking only the closest perturber into account, let us assume that stars are distributed in a virialized () Plummer (1911) sphere, whose gravitational potential exhibits the form:


where is the Plummer scale length, is the total mass of the Plummer sphere (i.e., the total mass of the cluster in our case), and is the distance to the cluster center. As such, around the cluster center where , the potential is roughly constant. The tidal force of the cluster as a function of is therefore


which reaches the maximum at the cluster center (i.e. ). At this point, the relative magnitude of the cluster tide normalized to the magnitude of the perturber tide is


For a typical cluster, the virial radius  pc. Since (Heggie & Hut, 2003), therefore  pc. For a close encounter, we adopt  AU (Adams et al., 2006), therefore . For this reason, we conclude that the perturbations experienced by a targeted planetary system in a star cluster is dominated by the contributions from its neighbor stars222In principle, the Galactic tide acting on the host star cluster is also acting on its planetary systems. In practice, however, we consider the Galactic tide in the solar neighborhood negligible for the evolution of planetary systems..

Figure 1: Schematic illustration of a planetary system consisting of a host star (S) and five planets (P1 to P5), being perturbed by three external perturbers (X1, X2, X3). The dotted circle indicates the boundary of the neighbor sphere, beyond which we ignore the influence of any perturbations (distance not to scale). The arrows indicate the directions and magnitudes of the perturbers’ velocities. In this example, the perturbations due to X1 and X2 are taken into account, whereas the perturbation due to X3 is neglected. Slow perturbers are more likely to remain in the neighbor sphere for longer time, causing stronger perturbations. Note that when a planet’s orbit is sufficiently excited, its orbit may intersect with the orbits of inner planets.

When solving the equation of motion of each planet (Eq. 1), we use rebound (Rein & Liu, 2012) to integrate the term (Eq. 2). The term (Eq. 3) is obtained by summing the pairwise gravitational accelerations from the host star to other stars in the star cluster, which are integrated by NBODY6++GPU (Wang et al., 2015a, 2016; Aarseth, 2003; Spurzem, 1999, and the references therein). The term (Eq. 4) is rather cumbersome to compute, as it requires both the positions of each star in the cluster, and the position of each planet in the planetary systems. Ideally, one would obtain the positions of all stars from NBODY6++GPU and sum up the pairwise accelerations from these stars to a given planet located at a specific position in its orbit. However, the dynamical timescales of star clusters ( Myr) is many orders of magnitudes longer than that of planetary systems (down to a few days for hot Jupiters). Each planetary orbit requires 100 integration steps (depending on the actual integrator) to maintain reasonable accuracy. In each step, the combined tidal force due to other stars needs to be evaluated with respect to the current position of each planet. The coupled system of planetary systems in star clusters is therefore highly hierarchical. It is not practical to integrate such systems with the small timesteps required.

In an analogous scenario where a star cluster evolves in a galactic tidal field, tidal tensors are often used to describe the tidal effects (e.g., Renaud et al., 2011; Rieder et al., 2013), which turns out to be reliable and accurate. Despite this analogy, tidal tensors may not be directly applicable to planetary systems that are perturbed by stellar encounters, due to strong and rapid fluctuations in the tidal field.

As mentioned above, neighbor stars dominate the perturbations of the planetary systems. This property allows us to dramatically reduce the computational costs substantially by considering only the perturbations due to the nearest neighbors, where . In the case where , only the nearest star is taken into account as the perturber, which is analogous to the Monte-Carlo scattering simulations, such as those of Hao et al. (2013), with the important difference that the perturbers are obtained from direct integration of the host star clusters where the multiple mechanisms such as dynamical evolution, stellar evolution, and galactic tidal fields can be taken into account. The case corresponds to an isolated planetary system, which we will use for validation and verification of our methods.

Throughout our study we adopt . The mass of the nearest perturber and its position with respect to the planetary system under consideration are communicated to the planetary system integrator. The planetary system integrator subsequently calculates the tidal force experienced by all planets at each integration step.

2.2 Methods and Optimization

Planetary systems and star clusters are two very different dynamical systems. Planetary systems are AU-scale few-body systems, whereas star clusters are parsec-scale many-body systems. The orbital periods of planets vary from several hours to a few hundred years, whereas the crossing timescale of star clusters are typically  Myr. A simulation of planetary systems are generally considered satisfactory when the relative energy error is , whereas in star cluster simulations this criterion is usually relaxed to . If a star has only one planet, it is possible to integrate the planet as a low-mass binary companion using the regularization technique (e.g. Kustaanheimo & Stiefel., 1965). This technique, however, becomes inadequate in handling multiplanetary systems embedded in star clusters.

It is therefore necessary to integrate these two types of systems using dedicated integrators that are optimized for for each of these dynamical problems, respectively. We employ rebound (Rein & Liu, 2012) for the integration of planetary systems using the IAS15 algorithm (Rein & Spiegel, 2015), which is optimized for handling close encounters. We use the NBODY6++GPU package (with a 4th-order Hermite algorithm) to integrate the star cluster. The perturbations to the targeted planetary system is computed according to the positions of the stars, obtained through the NBODY6++GPU code. Since the tidal force of the perturber depends sensitively on its distance to the targeted planetary system (especially during a close encounter), it is crucial to evaluate the exact positions of the perturbers accurately, with time steps comparable to the planet integration time step ( 1/100th of the shortest orbital period, or even shorter). In practice, it is difficult to determine the perturber’s exact position with such small time steps, as the star cluster integration time steps are orders of magnitude longer than the planet orbit integration time step. Different stars are usually assigned different integration time steps according to their dynamical “activities” – an integration scheme called individual time step scheme, designed to reduce the computational complexity to 333The actual computational complexity depends on the density profile of the star cluster.

In order to bridge the apparent gap of time steps between planetary systems and star clusters, we store the star cluster simulation data, and subsequently make interpolation of the positions of all stars to accommodate the time steps required by the planetary system integration. Farr et al. (2012) propose the Particle Stream Data Format (PSDF) scheme optimized for recording the output of general -body simulations that exploit individual time steps. The PSDF scheme records the data incrementally only when a particle updates its scale, and thereby reduce the redundancy. The data is presented with the YAML444 format, which is highly human readable. Inspired by this idea, we have developed an adaptive storage scheme called Block Time Step (BTS) storage scheme to incrementally store star cluster data at arbitrarily high spatial and temporal resolutions (Cai et al., 2015). The BTS scheme makes it possible to reconstruct the star cluster evolution with full details of stellar encounters with controllable snapshot sizes. To facilitate high-performance parallel access to the data files, we instead store the data with the HDF5555 data format. Accordingly, we carry out star cluster simulations separately and store the integration data with a temporal resolution comparable to 1000 yr per snapshot. Dynamical data such as positions (, velocities (), accelerations () and the first derivative of the accelerations () are stored for the purpose of interpolation.

The coupling between star cluster dynamics and planetary system dynamics is implemented within the AMUSE (Astrophysical Multipurpose Software Environment)666 framework. We construct a GPU-accelerated pseudo-gravitational dynamics interface H5nb6xx, which loads the BTS time series data stored in HDF5 files. At when a simulation starts, H5nb6xx reads the initial state of the star cluster, assigns an ensemble of initially identical multiplanetary systems to solar-type stars. Each planetary system is integrated by one rebound instance implemented in the AMUSE framework. The rebound instances advance their own planetary systems, inquire the positions and masses of the nearby perturbers at time . H5nb6xx responds to these inquiries by loading two adjacent snapshots at and , where . Accordingly, the particle states at are interpolated in parallel on the GPU using a set of seventh-order septic splines. Eventually, accurate positions of the host stars and their neighbors are obtained and transmitted into each rebound instance, which carries out the integration of planetary systems with the additional forces from the perturbers. As such, H5nb6xx and rebound communicate iteratively until the simulation completes. During the simulation, the snapshot of the coupled system is stored at a fixed time intervals for the purpose of optional restarting. The simulations are carried out using SiMon777Available at: (Qian et al. 2017, submitted), an open source Simulation Monitor for computational astrophysics.

A planet is identified as having escaped from its host star when its orbital eccentricity during at least five consecutive integration time steps888Physically speaking, a planet is unbound only if its eccentricity 1.0. In practice, however, this causes numerical difficulties for the integration, and therefore we relax this criterion to in the last five consecutive time steps. . In such a case, the planet is removed from the planetary system and becomes a free-floating planet (FFP). Depending on the escape velocity, the FFP may remain in the host star cluster and get recaptured (Wang et al., 2015b) or escape from the host star cluster (Zheng et al., 2015) .

3 Initial Conditions

3.1 Initial Conditions for Star Clusters

While the total number of stars in star clusters ranges between and , the likelihood for dense globular clusters to bear planets is low (see, e.g., Gilliland et al., 2000; Gonzalez et al., 2001; Masuda & Winn, 2017). We therefore study intermediate-size open clusters consisting of and stars, which are comparable with the total masses of M67 (e.g., Hurley et al., 2005), NGC 6811 (e.g., Meibom et al., 2013) and the Westerlund 1 (e.g., Portegies Zwart et al., 2010) open clusters, respectively. The virial radii of all models are initialized at 1 pc, resulting in a range of central stellar densities that vary by a factor of four. We adopt a Kroupa (2001) stellar initial mass function with a mass range of , which corresponds to a mean stellar mass of . The stellar positions and velocities are sampled from the (Plummer, 1911) model in virial equilibrium (i.e., virial ratio ). We do not include primordial binaries or primordial mass segregation. We evolve our models for 50 Myr, a time during which stellar evolution can be ignored for the mass range under consideration.

3.2 Initial Conditions for Planetary Systems

Current exoplanet data999see e.g., show that planetary systems are immensely diverse: eccentricities are widely spread and the distributions of semi-major axis and mass seem to exhibit complex patterns. In this study we restrict ourselves to systems in which all planets have equal mass, are initially on coplanar and have circular orbits. We further assume that their semi-major axes are equally spaced in terms of their mutual Hill radii (dubbed EMS: Equal Mass and equal Separation in terms of mutual Hill radii, see Zhou et al. 2007; Hao et al. 2013). The scaled separation of the planetary orbits is


where is the semi-major axis of the -th planet (), and is the mutual Hill radius:


The quantity is the mass ratio between a planet and its host star.

Model () [AU] 2k 8k 32k Remarks
Model I 10 1 5 [1, 2.5, 6.3, 15.9, 39.7] 50 50 200 compact multiple-Jupiters
Model II 100 5 [1, 2.5, 6.3, 15.9, 39.7] 50 50 200 compact multiple-Earths
Model III 10 1 5 [5.2, 13.04, 32.7, 82.2, 206.2] 50 50 200 wide multiple-Jupiters
Model IV 100 5 [5.2, 13.04, 32.7, 82.2, 206.2] 50 50 200 wide multiple-Earths
Table 2: Initial conditions of the planetary system models. The host stars are chosen from solar-type stars in the cluster. Model I and Model II are “compact” models, with the inner edge comparable to the Earth’s orbit, and the outer edge comparable to Pluto’s orbit. Model III and Model IV are the “wide” models, with the inner edge comparable to Jupiter’s orbit, and the outer edge of  AU comparable to of Sedna’s orbital semi-major axis. Model I and Model III are multiple-Jupiter models where planet-planet interactions are important; Model II and Model IV are multiple-Earth model where where the gravitational interactions among planets can be safely ignored most of the time. The number of planetary systems in each ensemble are listed in the 2k, 8k, 32k columns, respectively. In total, there are 1200 individual planetary system simulations, and each simulation is carried out for 50 Myr.

We study four EMS models as detailed in Table 2. All models adopt planets orbiting around solar-type stars (). Kokubo & Ida (1998) suggest that separations of are typical, and therefore is adopted in the multiple-Jupiter models (Model I and Model III). The inner edge of Model I is comparable to the Earth’s orbit, while that of Model III is comparable to Jupiter’s orbit. Model II and Model IV are obtained by reducing the planet mass of Model I and Model III by a factor of 1000 while keeping all other parameters unchanged (i.e., and , thus comparable to Earth mass). According to Eq. 10 and 11, the corresponding separation is in this configuration. The configurations serve as comparisons that can be used to disentangle the effects of stellar encounters and planet-planet interactions on the dynamical evolution of planetary systems in star clusters. The wide orbits in Model III and Model IV are useful to provide insights into how external perturbations shape the orbits of objects with large semi-major axes, such as trans-Neptunian objects (TNOs).

Since planetary systems are assigned to solar-type stars, their spatial distribution is random across the cluster. Given that the density profile of Plummer model follows (Plummer, 1911), planetary systems in the central region of the star cluster are immersed in much higher stellar densities than their siblings in the outskirts, therefore they experience much stronger perturbations and more frequent encounters. Our simulations focus on the post-formation regime where the protoplanetary disks have already dissipated. Gas drag is therefore not taken into account. Tidal circularization is important when a planetary orbit is excited to high eccentricity and the periapsis with respect to the host star is small (Chatterjee et al., 2008). Tidal circularization is particularly important when the planetary system is subject to Lidov-Kozai cycles. This is a common mechanism for producing hot Jupiters (e.g., Shara et al., 2016; Hamers et al., 2017). Tidal circularization damps the eccentricity of a perturbed planet, which protects the planet and its planetary system. However, in our -body simulations we do not take this effect into account.

4 Results

4.1 Statistics of Stellar Encounters

Over the timespan of 50 Myr in our simulations, most cluster stars have experienced dozens of crossing times and roughly a half of relaxation time. Therefore, planetary systems within the star cluster will have experienced substantial number of encounters, depending on the neighbor density in the vicinity101010In this paper, we define an encounter as the time span since a perturber becomes the closest star to the targeted planetary system, until it is replaced by another star even closer to the planetary system. In other words, each change of neighbor star is considered as an encounter, regardless of its strength or duration.. Inspired by the approach of modeling the encounter between the host star and a external perturber with a two-body problem (e.g., Spurzem et al., 2009; Heggie, 2006), we carry out our analysis with such a model, and quantify the strength of each encounter with dimensionless parameters and . The quantity is the relative speed of the perturber with respect to the host star, scaled to the average orbital speed of the outermost planet. The parameter is the ratio of the perturbation timescale to the planet orbital timescale, defined as


where is the total mass of the perturbed planetary system, is the mass of the perturber, is the pericenter distance of the perturber with respect to the planetary system center of mass, and is the semi-major axis of the perturbed planetary orbit.

Figure 2 shows the distances (scaled to the semi-major axes of the outermost planet) and the velocities (scaled to the velocity at infinity of the perturber) of three ensembles of stellar encounters in the and clusters during the 50 Myr timespan of simulations whenever a perturber reaches the periapsis.

The frequencies of encounters (including nearly parabolic, hyperbolic, impulsive, tidal encounters; separated by the grey curves in the figure) increase as increases, indicating that more frequent encounters are expected in denser cluster environments. Indeed, as all our modeled clusters have identical initially, a larger results in higher stellar density, and consequently smaller since (see Eq. 12). The strengths of encounters are indicated with the five dashed vertical lines in each panel, in which smaller values are associated with stronger encounters. As compared with hyperbolic encounters, the near parabolic regime is more destructive to the targeted planetary systems – binary systems are formed between the perturber and the host star, with the possibility of triggering Kozai-Lidov oscillations (Naoz et al., 2011; Hamers & Portegies Zwart, 2016). The cumulative frequency distribution of is shown in Figure 3. The cumulative frequency spectra shift leftward with the increment of , showing that the average strengths of encounters increases in dense stellar environments.

We will see below, that the hyperbolic regime, albeit weak, affects the planetary systems cumulatively by series of subsequent relatively weak to moderate encounters (). Stars in their clusters experience Coulomb-like scattering – the contributions of rare and strong encounters are comparable to the contributions of the cumulative effect of a series of frequent and weak encounters.

Figure 2: Ensembles of stellar encounter parameters at the periapsis for the and star clusters (from left to right). The vertical gray line coincides with , where is the semi-major axes of the outermost planets ( AU). If , the encounter is tidal. The diagonal grey line coincides with the hyperbolic eccentricity . Above this line encounters are hyperbolic, below the line they are nearly parabolic. The grey curve separates the regime of adiabatic encounters (, where is the circular velocity at ) and impulsive encounters (). The five blue vertical dashed lines from left to right correspond to and , respectively (assuming that the perturber mass is , and the mass of the host star is according to the initial conditions). The colors of dots in each panel correspond to the counts of encounters. The encounter parameters are collected when the perturber reaches .
Figure 3: Cumulative frequency distribution of the encounter strength parameter for the and hosting clusters. The cumulative curves shift leftward as the total number of particles grows higher, which shows that strong encounters are more common in denser star clusters.

4.2 Perturbed Planetary Systems

Multiplanetary systems are fragile when experiencing stellar encounters in the star cluster environments (e.g., Portegies Zwart & Jílková, 2015; Hao et al., 2013). The planet survival rates depend on the frequency of close stellar encounters, which in turn depends on the stellar density. Planetary systems close to the dense cluster center are more frequently perturbed than those in the cluster outskirts.

Figures (4 - 7) depict the “microscopic” behaviors of four perturbed planetary systems. The first planetary system, shown in Figure 4, is a Model I multiple-Jupiter system in the host cluster, and its evolution is significantly affected by a perturbation at  Myr. The outermost planet P5 is immediately ejected, and the second outermost planet undergoes substantial eccentricity excitation and semi-major axis expansion. As such, the perturber exchanges energy and angular momentum with this planetary system, leading P4 into a retrograde orbit, and P1-P3 tightly coupled as a whole in the inner region of the system (as shown at the bottom panel in particular). For comparison, Figure 5 presents the behavior of a Model II multiple-Earth planetary system orbiting exactly the same host star in the same cluster of Figure 4. The same close encounter at  Myr liberates P4 and P5 simultaneously. However, the system does not exhibit apparent pattern of planet-planet interaction as seen in Figure 4. Rather, planet are evolving mostly independently before and after the close encounter. In the similar way, the evolution of a Model III and a Model IV planetary systems are shown in Figure 6 and Figure 7, respectively. The wide orbits of P5 in these two models are so sensitive to external perturbations that even the mild excitation of eccentricities can be used to trace the history of weak stellar encounters (cf. Portegies Zwart & Jílková, 2015). In both models, eccentricity excitations occur almost exclusively when the planetary system plunges into the dense cluster center, a region where the frequency and strength of encounters significantly increase. Both Figure 4 and Figure 6 exhibit interesting behavior of inclination evolution: slow anti-phase variation of inclinations are seen between the outermost planets and the inner planets combined; the inner planets undergo more rapid anti-phase oscillation of orbital inclinations among themselves. This behavior suggests that the secular evolution of planetary system is largely unaffected by weak encounters, but they can be dramatically changed by strong encounters through the injection of orbital energy and angular momentum.

Figure 4: Orbital element evolution of a Model I planetary system in the cluster. The system has 5 planets initially, which are initially on circular and coplanar orbits. In the top panel, the semi-major axis of each planet is plotted in colored curved as a function of time (in log-scale), while the thick gray curve indicates the distance from the host star to the cluster center (axis on the right, in parsec); in the middle panel, the evolution of the eccentricity of each planet is plotted as a function of time, and the thick gray curve is the distance of the perturber to the host star (in AU, log-scale); at the bottom panel the inclination of each planet with respect to the initial orbital plane is plotted, and the thick gray curve is the acceleration due to the perturber (log-scale, normalized to the acceleration due to the host star). In this particular planetary system, a close encounter encounter occurs at about  Myr, which ejects the outermost planet (P5) immediately, and excites P4 to . This close encounter results in a prograde orbit () of P4, and also strengthens the planet-planet scattering among the remaining planets.
Figure 5: The same host star in the same cluster as in Figure 4, but for a Model II planetary system. P4 and P5 are ejected immediately.
Figure 6: Same as Fig. 5, but for a Model III planetary system in the cluster. The planetary system orbits around the cluster center in approximately eccentric orbits. As it dives into the cluster center, the frequency of perturbations increases significantly. Likewise, the planetary system remains roughly unperturbed when in the outskirts of the cluster. Consequently, excitation are more likely to occur during around the dense cluster center.
Figure 7: Same as Fig. 6, but for a Model IV planetary system in an cluster.

These are only a few examples of the intricate evolution of planetary systems due to both internal and external perturbations. In general we observe in our simulations that these effects depend on the orbit of the planetary system in the star cluster. We can however, deduce several general behaviors. Shortly after the dissipation of protoplanetary disks, planets have obtained nearly circular and coplanar orbits. External perturbations initialized by stellar flybys pump up the eccentricities of planetary orbits. The growth in eccentricity is proportional to the semi-major axes as (Spurzem et al., 2009, Eq. 13), so outer planets with wide orbits are more vulnerable. Indeed, outer planets are excited more quickly, as can be seen in Figures (4 - 7). As the eccentricity of the outer planets grows due to external perturbations, its orbit intersects with the orbits of inner planets. Orbit crossings lead to planet-planet scattering, which in turn result in an inward propagation of eccentricity. This suggests that multiplicity does contribute to the vulnerability of planetary systems. As such, we predict that the stability of a planetary system can be compromised if there exist massive planets with large semi-major axes. Interestingly, while the outermost planets are generally most vulnerable to external perturbations, they may not necessarily be the first ones to get ejected, because the excitation process is complicated by planet-planet interactions and orbital phases.

Moderate encounters (), albeit weak, can still leave marks on the targeted planetary systems (see especially Figure 6). Each encounter results in a small and , causing the orbit to gradually depart from its initial circular and coplanar state, a state corresponding to the maximum possible orbital angular momentum. The orbital angular momentum is partially taken away during the onset of an encounter, and therefore results in the growth of angular momentum deficit (AMD) (see, Laskar, 1997). The AMD of outer planets is partially absorbed by inner planets, consequently causing the inner planets to develop eccentric orbits (Davies et al., 2014). As such, eccentricities are propagated from outer planets to inner planets. Due to the Coulomb-like random scattering among stars in the cluster, the frequency of moderate encounters (as seen in Figure 2) much higher than strong encounters (), especially in denser clusters. These encounters are able to gradually pump up the AMD of a targeted planetary system by repeatedly perturbing the outermost planet. The value of AMD limits the maximum eccentricity and inclination an individual planet can attain. When the AMD is sufficiently large, the planetary system tends to reach an equipartition of AMD (Wu & Lithwick, 2011). Inner planets with small semi-major axis have relatively small orbital angular momentum, and therefore it only contribute to a small fraction of the total AMD. As a consequence of the AMD equipartition, the inner planet is “forced” to contribute as much AMD as possible, which in turn be driven to extreme orbits. It is worthy to point out that absorbing AMD takes time, and therefore if an inner planet is indirectly ejected because it is perturbed by an outer orbits, the ejection may happen well after the perturber departs from the periapsis. Additionally, an inner planet may be ejected earlier than the outer planets, because it has lower AMD capacity. Eventually, all planets in the system obtain high eccentricities, which results in the loss of stability of the entire system.

We conclude that planets can be liberated immediately through very strong encounters (). An alternative channel to eject planets is through the cumulative effect of multiple moderate encounters (). Distance encounters ( > 50) have no direct implication to stability. Furthermore, planet-planet interactions are the catalysts of the destruction of a secularly interacting planetary system, such as in our multiple-Jupiter models.

4.3 Statistical Behavior of Perturbed Planetary System Ensembles

Planetary systems are chaotic few-body systems sensitive to the initial conditions. In order to determine the contributions of external effects (stellar encounters) and internal effect (planet-planet interactions), we plot the planet survival rates111111Survival rates of planets in a given host star cluster: defined as , where is the total number of bound planets at time . as a function of time for different ensembles of simulations. As shown in Figure 8, six ensembles of simulations are plotted in their corresponding panels, with each of the planets distinguished with different colors. The final survival rates of each ensemble of simulations are listed in Table 3. Apparently, for both multiple-Jupiter models (Model I and Model III) and multiple-Earth models (Model II and Model IV), planetary systems in denser stellar environments suffer from higher ejection rates. When keeping the initial arrangement of semi-major axis fixed, the survival rates of multiple-Earth models are substantially higher. When keeping the initial orbital separation fixed, the compact orbit models (Model I and Model II) have significantly higher survival rates compared to the wide orbit models (Model II and Model IV).

Therefore, we can conclude that both internal effect and external effect play important roles in the evolution of planetary systems in star clusters. Furthermore, outer planets tend to be ejected more rapidly; tight orbit inner planets are better protected against ejections.

Models Model I Model II Model III Model IV
2k 0.984 0.988 0.920 0.968
8k 0.940 0.968 0.792 0.900
32k 0.834 0.916 0.550 0.771
Table 3: The overall survival for each ensemble of simulations at  Myr.
Figure 8: Survival rates as a function of time for planetary systems in (top row), (middle row), and clusters. The three panels on the left column are Model III planetary systems, and the three panels on the right column are Model IV planetary systems. Different planets are distinguished with different colors, and the thick black curve is the overall survival rates, defined as the ratio between the total number of ejected planets at the current time and the total number of of planets at .

Table 4 presents the fractions of planetary systems with each of them having exactly surviving planets () by the end of the simulation ( Myr). A comparison between these fractions for the multiple-Jupiter and multiple-Earth models in the and clusters is presented in Figure 9.

Model up to  Myr
0 1 2 3 4 5
I 0.02 0.02 0.0 0.0 0.0 0.96
II 0.00 0.00 0.02 0.0 0.0 0.98
III 0.00 0.06 0.02 0.02 0.06 0.84
IV 0.00 0.00 0.02 0.02 0.06 0.90
I 0.0 0.04 0.02 0.0 0.08 0.86
II 0.0 0.0 0.02 0.04 0.02 0.92
III 0.02 0.06 0.16 0.04 0.14 0.58
IV 0.02 0.02 0.02 0.04 0.18 0.72
I 0.09 0.065 0.05 0.04 0.04 0.715
II 0.015 0.035 0.03 0.025 0.09 0.805
III 0.1 0.185 0.245 0.1 0.11 0.265
IV 0.035 0.04 0.085 0.145 0.265 0.43
Table 4: Fraction of planetary systems with the given number of planets survived in each planetary system , at the end of each simulation ( Myr).
Figure 9: Fractions of planetary systems as a function of the number of surviving planets in each system at  Myr. The data is from Table 4. Circles: multiple-Jupiter models; triangles: multiple-Earth models; green: compact models; blue: wide models.

The excitation process of the ensembles of planetary systems can be inspected with Figure 10, where a grid of snapshots at the space is presented at four time checkpoints ( 1, 5, 10 and 50 Myr) for different models. Planet migrations (changes in ) are seen especially among the planets in the multiple-Jupiter models with moderate to high eccentricities. For comparison, the corresponding snapshots of the multiple-Earth models are shown in the three bottom panels, where planet migrations are less pronounced. The population of highly-eccentric planets increases as increases, until the targeted planetary systems have experienced sufficient perturbations to even induce the eccentricity of the innermost planet, or until the AMD has reached an equipartition across the entire planetary system. The distribution of eccentricities at  Myr is presented in Figure 11.

Figure 10: Time frames of the space at  1 Myr (col.1), 5 Myr (col.2), 10 Myr (col.3) and 50 Myr (col.4). The 1-3 rows are for the Model III wide multiple-Jupiter planetary systems with hosting cluster of , and , respectively. The 4-6 rows are for the Model IV wide multiple-Earth models in , and clusters, respectively.
Figure 11: Same with Figure 10, but for the distribution of eccentricities.

We follow the changes in the orbital eccentricities as a consequence of each stellar encounters. As shown in Figure 12, very strong encounters are rare. Among all models, only of encounters are sufficiently strong to cause excitations of . This suggests that planet ejections is a cumulative process – planets are gradually excited by a number of subsequent encounters with relative small , except for a few very strong encounters that ionize the planets immediately. These findings are consistent with the results in Spurzem et al. 2009 and Hao et al. 2013.

Figure 12: The cumulative frequency of encounters that lead to a given change of eccentricity . Only are plotted. The three curves correspond to three Model IV ensembles in and clusters, respectively. Weak encounters with small dominant the frequency spectra for all models. Very strong encounter causing are rare. Due to the small number of close encounter events, the result is different from the the and results.

Hot Jupiters (HJs) can be produced with the combined effects of stellar encounters and planet-planet interactions. In the cluster, of Model I planets have developed orbital features of  AU and  AU. Tidal circularization can be efficient when these planets are around their orbital periapsis, which in turn provides a mechanism to produce HJs. The HJs rate predicted in (Shara et al., 2016) is , higher than our results. However, we note that our simulations (50 Myr) are much shorter than the simulations in (Shara et al., 2016) (1 Gyr). We therefore suspect that if we were not restricted by the computational costs, our results will be more consistent with (Shara et al., 2016) if we carry out our simulations for longer time.

All planets are initially coplanar, and external torques outside the orbital plane are exerted by perturbers from arbitrary directions. Mutual inclinations form as a natural byproduct of stellar encounters. Therefore, the excitations of orbital eccentricities are usually accompanied by the excitations of inclinations, which is consistent with the results from Parker & Quanz (2012). Figure 13 shows snapshots a grid of space, where the inclinations are measured with respect to the initial orbital planes of planets. A small number of planets have been induced to retrograde orbits. The fraction of retrograde orbits seems to be slightly higher than the multiple-Earth models (Model II and Model IV): there are 4 retrograde planets () in the 32k Model IV ensemble, but only 2 () in the 32k Model III ensemble. Also, in the 8k cluster, the Model IV ensemble has two retrograde orbits, comparing to no retrograde orbits in the corresponding Model III ensemble. Given the low numbers of retrograde orbits, we are not sure yet whether this comparison is statistically significant. We believe that this result is consistent with theoretical understanding: planets in multiple-Earth systems are very weakly coupled (), and carry much less orbital angular momentum than planets in the multiple-Jupiter systems. Consequently, it is easier to flip over the orbits in multiple-Jupiter systems.

Figure 13: Same with Figure 10, but for the space. Each panel is divided into two regimes by the blue horizontal dashed line: prograde orbits () and retrograde orbits ().

5 Conclusions

The collapse of giant molecular clouds triggers star formation in clustered environments. Protoplanetary disks, which are the progenitors of planets, form around newborn stars as byproducts of this process. Both theoretical predictions and observation suggest that planets are common, but only very few exoplanets have been discovered in star clusters. To better understand this apparent dichotomy, we carry out this exploratory study with a grid of simulations to test the dynamical stability of multiplanetary systems in intermediate-mass open clusters. We simulate three host star cluster environments with and Plummer models. Each of these models has the same virial radius  pc. The mass spectra of these clusters are sampled with a Kroupa (2001) initial mass function. For each star cluster model, we distribute an ensemble of an ensemble of equal-mass and equally separated (in terms of mutual Hill radii) planetary systems around solar-type stars. Each planetary system has either five (multiple-Jupiter model) or five (multiple-Earth model) planets distributed in initially circular and coplanar orbits. The star clusters are integrated using the direct -body code NBODY6++GPU (Wang et al., 2015a, 2016; Spurzem, 1999), while the planetary systems are evolved with rebound (Rein & Liu, 2012) using the IAS15 integrator (Rein & Spiegel, 2015). The star cluster data is stored in the Block Time Step storage scheme (Cai et al., 2015). After performing interpolation on the GPU, the perturbation information is passed to rebound within the AMUSE (Portegies Zwart et al., 2009, 2013; McMillan et al., 2012; Pelupessy et al., 2013) framework. Our conclusions can be summarized as follows.

  1. We quantify the strength of each stellar encounter with the dimensionless parameters and , where is essentially the ratio between the perturbation timescale and the orbital timescale, and is the speed of the perturber at infinity. The peak frequency distribution of shifts to a lower values in denser clusters, indicating that the encounters are on average stronger in denser clusters. Moreover, stellar encounters are more frequent in denser and more massive clusters.

  2. The dynamical evolution of planetary systems is sensitive to external perturbations. Consequently, the planet survival rate declines in denser clusters: for clusters with and , the planet survival rates of the compact multiple-Jupiter systems (Model I) are and , respectively, and the survival rates for wide multiple-Jupiter systems (Model III) are and , respectively. Similarly, when evolving the compact multiple-Earth systems (Model II) in the and clusters, the corresponding survival rates are and , respectively, and for wide multiple-Earth systems (Model IV) and , respectively. In terms of the number of planets survived in each planetary system, of wide multiple-Jupiter systems and of wide multiple-Earth systems are able to keep all their planets cluster. This fraction drops to and in the denser cluster. Therefore, we believe that young low-mass star clusters will be prominent sites for next generation planet detection surveys, but the likelihood of detecting planets in dense globular clusters such as 47 Tuc would be low.

  3. External perturbations constrain the maximum sizes of planetary systems. The wider a planet’s orbit is, the more vulnerable it is to external perturbations. In all star clusters environments used in our simulations, the survival rates of the wide models (Model III and Model IV) are much lower than the corresponding compact models (Model I and Model II), even though they evolve in exactly the same host star in the same cluster. As such, we predict that planets in denser star clusters will have smaller orbits, which actually allow them to be detected relatively easily.

  4. Planet-planet interactions are the catalysts of planetary system disruptions. We compare the dynamical evolution of the multiple-Jupiter systems and multiple-Earth systems in identical stellar environments, and found that the multiple-Earth systems in the all clusters have substantially higher survival rates than the multiple-Jupiter systems. In multiple-Earth systems, inner planets absorb the angular momentum deficit of the outer planets through mutual interactions, and consequently leads to transfer of eccentricities from outer planets to inner ones. In multiple-Earth systems, planet mutual interactions are negligible, hence the eccentricity excitations of each planet is solely induced by external perturbations.

  5. The excitation process is cumulative and gradually results in planet ejections. While very strong encounters can cause instantaneous ejection of planets, they are rare: only about of the encounters are strong enough to cause an eccentricity excitation of . In most cases planets are excited gradually by a series of moderate or weak encounters.

  6. One direct consequence of stellar encounters is that they change both the magnitudes and the directions of orbital angular momenta. The changes in angular momenta are primarily caused by perturbers with offsets of the orbital planets, which is usually the case. In the EMS model for planetary systems we distribute planets initially on coplanar circular orbits. As such, the orbital angular momenta are proportional to , where is the semi-major axis. However, outer planets lose their angular momenta quickly though the rapid eccentricity excitations, and therefore they tend to have a higher inclinations than the inner planets. We observe on average planets () with retrograde orbits in each multiple-Jupiter ensemble, and in the multiple-Earth ensembles. The frequency of planets with retrograde orbits increases in denser stellar environments. When planets with highly eccentric orbits approach the host stars, their eccentricities could be damped near the periapsis through the tidal damping effect, resulting to very small orbits (e.g. hot Jupiters). The tidal damping effect does not affect their induced inclinations. We speculate that this may produce an alternative channel of producing large spin-orbit misalignments (as opposed to the Kozai-Lidov mechanism), which is observed among many extrasolar planetary systems through e.g., the Rossiter-McLaughlin effect.

Not all star clusters are as compact as our models. As such, our results provides an upper limit of planet ejection rates in such environments.

It is worthwhile to point out that the background cluster potential may have implications to planetary system stability. For an isolated planetary system, the perturber and the host stars are interacting in the Keplerian potential, and the total energy is conserved. In the star cluster environments (especially in the cluster center), however, stars are being scattered randomly, exhibiting Brownian motion in the background cluster potential. Therefore, our simulations of perturbing planetary systems in star clusters cannot be simplified as isolated planetary systems being perturbed by a series of Keplerian stellar encounters. The background cluster potential affects the planet stability by affecting the trajectories of the perturber. Nevertheless, a quantitative analysis of this effect is beyond the scope of this paper.

This study has been limited to pure dynamical interaction between member stars in star clusters and multiplanetary systems. Planets are assumed to be of equal mass, arranged in initially coplanar circular orbits with equal separation in terms of their mutual Hill radii. Eccentricity damping due to the protoplanetary disks and/or tidal circularization process may contribute to the robustness of planetary systems. In addition, our host star clusters are sampled with the idealized Plummer model in virial equilibrium (). In reality, planetary systems are immensely diverse in terms of orbital architectures and mass spectra. Host star clusters may depart from and exhibit substructures and (Zheng et al., 2015; Portegies Zwart, 2016). It is also possible that the chaotic behavior observed in some mildly perturbed planetary systems only manifests itself when simulating for a more extended time.

Nevertheless, this study aims to highlight the dynamical consequences of stellar encounters for planetary systems in star clusters. The results are also potentially insightful for understanding the frequency of free-floating planets, which are currently subjected to strong observational bias due to the difficulty of observations.


We thank the anonymous referee for her/his comments that helped to improve the manuscript considerably. We thank Adrian Hamers for commenting on the manuscript. We thank Douglas N.C. Lin, Fred C. Adams, Sverre Aarseth, Roberto Capuzzo-Dolcetta, Rosemary Mardling, and Beibei Liu for useful discussions. We thank Long Wang for graciously providing support for the direct -body code NBODY6++GPU. We are grateful for the support from the AMUSE team, in particular Inti Pelupessy, Arjen van Elteren, Nathan de Vries, Michiko Fujii and Lucie Jílková for very helpful discussions. MBNK acknowledges the AMUSE team for supporting his visit to Leiden University. We are grateful for support by Sonderforschungsbereich SFB 881 “The Milky Way System” of the German Research Foundation (DFG), through subproject Z2 and the GPU cluster Milky Way at FZ Jülich, and for the support of the visit of MXC in Heidelberg. We acknowledge support by NAOC CAS through the Silk Road Project and (RS) through the Chinese Academy of Sciences Visiting Professorship for Senior International Scientists, Grant Number . The special GPU accelerated supercomputer laohu at the Center of Information and Computing at National Astronomical Observatories, Chinese Academy of Sciences, funded by Ministry of Finance of People’s Republic of China under the grant , has been used for some of the largest simulations. M.B.N.K. was supported by the National Natural Science Foundation of China (grants 11010237, 11050110414, 11173004, and 11573004). This research was supported by the Research Development Fund (grant RDF-16-01-16) of Xi’an Jiaotong-Liverpool University (XJTLU).

This work was supported by the Netherlands Research Council NWO (grants #643.200.503, #639.073.803 and #614.061.608) by the Netherlands Research School for Astronomy (NOVA). This research was supported by the Interuniversity Attraction Poles Programme (initiated by the Belgian Science Policy Office, IAP P7/08 CHARM) and by the European Union’s Horizon 2020 research and innovation programme under grant agreement No 671564 (COMPAT project).


  • Aarseth (2003) Aarseth, S. J. 2003, Gravitational N-Body Simulations, Cambridge, UK: Cambridge University Press
  • Adams & Laughlin (2001) Adams, F. C., & Laughlin, G. 2001, Icarus, 150, 151
  • Adams et al. (2004) Adams, F. C., Hollenbach, D., Laughlin, G., & Gorti, U. 2004, ApJ, 611, 360
  • Adams et al. (2006) Adams, F. C., Proszkow, E. M., Fatuzzo, M., & Myers, P. C. 2006, ApJ, 641, 504
  • Adams (2010) Adams, F. C. 2010, ARA&A, 48, 47
  • Anderson et al. (2013) Anderson, K. R., Adams, F. C., & Calvet, N. 2013, ApJ, 774, 9
  • Applegate (1986) Applegate, J. H. 1986, ApJ, 301, 132
  • Backer et al. (1993) Backer, D. C., Foster, R. S., & Sallmen, S. 1993, Nature, 365, 817
  • Bally et al. (2000) Bally, J., O’Dell, C. R., & McCaughrean, M. J. 2000, AJ, 119, 2919
  • Binney & Tremaine (2008) Binney, J., & Tremaine, S. 2008, Galactic Dynamics: Second Edition, Princeton University Press, Princeton, NJ USA, 2008
  • Bodenheimer & Pollack (1986) Bodenheimer, P., & Pollack, J. B. 1986, Icarus, 67, 391
  • Bonnell et al. (2001) Bonnell, I. A., Smith, K. W., Davies, M. B., & Horne, K. 2001, MNRAS, 322, 859
  • Boss (1997) Boss, A. P. 1997, Science, 276, 1836
  • Brucalassi et al. (2014) Brucalassi, A., Pasquini, L., Saglia, R., et al. 2014, A&A, 561, LL9
  • Brucalassi et al. (2016) Brucalassi, A., Pasquini, L., Saglia, R., et al. 2016, A&A, 592, L1
  • Cai et al. (2015) Cai, M. X., Meiron, Y., Kouwenhoven, M. B. N., Assmann, P., & Spurzem, R. 2015, ApJS, 219, 31
  • Cai et al. (2016) Cai, M. X., Gieles, M., Heggie, D. C., & Varri, A. L. 2016, MNRAS, 455, 596
  • Chandrasekhar (1942) Chandrasekhar, S. 1942, Chicago, Ill., The University of Chicago press
  • Chatterjee et al. (2008) Chatterjee, S., Ford, E. B., Matsumura, S., & Rasio, F. A. 2008, ApJ, 686, 580-602
  • Clarke et al. (2000) Clarke, C. J., Bonnell, I. A., & Hillenbrand, L. A. 2000, Protostars and Planets IV, 151
  • Davies & Sigurdsson (2001) Davies, M. B., & Sigurdsson, S. 2001, MNRAS, 324, 612
  • Davies et al. (2014) Davies, M. B., Adams, F. C., Armitage, P., et al. 2014, Protostars and Planets VI, 787
  • Faber et al. (2010) Faber, N. T., Stibbe, D., Portegies Zwart, S., McMillan, S. L. W., & Boily, C. M. 2010, MNRAS, 401, 1898
  • Farr et al. (2012) Farr, W. M., Ames, J., Hut, P., et al. 2012, New Astron., 17, 520
  • Fujii et al. (2007) Fujii, M., Iwasawa, M., Funato, Y., & Makino, J. 2007, PASJ, 59, 1095
  • Gieles et al. (2007) Gieles, M., Athanassoula, E., & Portegies Zwart, S. F. 2007, MNRAS, 376, 809
  • Gilliland et al. (2000) Gilliland, R. L., Brown, T. M., Guhathakurta, P., et al. 2000, ApJ, 545, L47
  • Gonzalez et al. (2001) Gonzalez, G., Brownlee, D., & Ward, P. 2001, Icarus, 152, 185
  • Hamers & Portegies Zwart (2016) Hamers, A. S., & Portegies Zwart, S. F. 2016, MNRAS, 459, 2827
  • Hamers et al. (2017) Hamers, A. S., Antonini, F., Lithwick, Y., Perets, H. B., & Portegies Zwart, S. F. 2017, MNRAS, 464, 688
  • Hénon (1971) Hénon, M. H. 1971, Ap&SS, 14, 151
  • Hao et al. (2013) Hao, W., Kouwenhoven, M. B. N., & Spurzem, R. 2013, MNRAS, 433, 867
  • Heggie & Rasio (1996) Heggie, D. C., & Rasio, F. A. 1996, MNRAS, 282, 1064
  • Heggie et al. (1998) Heggie, D. C., Giersz, M., Spurzem, R., & Takahashi, K. 1998, Highlights of Astronomy, 11, 591
  • Heggie & Hut (2003) Heggie, D., & Hut, P. 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics, Cambridge University Press
  • Heggie (2006) Heggie, D. C. 2006, Few-Body Problem: Theory and Computer Simulations
  • Hurley & Shara (2002) Hurley, J. R., & Shara, M. M. 2002, ApJ, 565, 1251
  • Hurley et al. (2005) Hurley, J. R., Pols, O. R., Aarseth, S. J., & Tout, C. A. 2005, MNRAS, 363, 293
  • Jílková et al. (2015) Jílková, L., Portegies Zwart, S., Pijloo, T., & Hammer, M. 2015, MNRAS, 453, 3157
  • Kokubo & Ida (1998) Kokubo, E., & Ida, S. 1998, Icarus, 131, 171
  • Kouwenhoven et al. (2014) Kouwenhoven, M. B. N., Goodwin, S. P., de Grijs, R., Rose, M., & Kim, S. S. 2014, MNRAS, 445, 2256
  • Kouwenhoven et al. (2016) Kouwenhoven, M. B. N., Shu, Q., Cai, M. X., & Spurzem, R. 2016, Mem. Soc. Astron. Italiana, 87, 630
  • Kroupa (2001) Kroupa, P. 2001, MNRAS, 322, 231
  • Kustaanheimo & Stiefel. (1965) Kustaanheimo, P. & Stiefel, E. 1965, J. Reine Angew. Math., 218, 204
  • Lada & Lada (2003) Lada, C. J., & Lada, E. A. 2003, ARA&A, 41, 57
  • Laskar (1997) Laskar, J. 1997, A&A, 317, L75
  • Laughlin & Adams (1998) Laughlin, G., & Adams, F. C. 1998, ApJ, 508, L171
  • Li & Adams (2015) Li, G., & Adams, F. C. 2015, MNRAS, 448, 344
  • Li et al. (2015) Li, Y., Kouwenhoven, M. B. N., Stamatellos, D., & Goodwin, S. P. 2015, ApJ, 805, 116
  • Li et al. (2016) Li, Y., Kouwenhoven, M. B. N., Stamatellos, D., & Goodwin, S. P. 2016, ApJ, 831, 166
  • Liu et al. (2013) Liu, H.-G., Zhang, H., & Zhou, J.-L. 2013, ApJ, 772, 142
  • Lovis & Mayor (2007) Lovis, C., & Mayor, M. 2007, A&A, 472, 657
  • Lynden-Bell & Wood (1968) Lynden-Bell, D., & Wood, R. 1968, MNRAS, 138, 495
  • Masuda & Winn (2017) Masuda, K., & Winn, J. N. 2017, AJ, 153, 187
  • Malavolta et al. (2016) Malavolta, L., Nascimbeni, V., Piotto, G., et al. 2016, A&A, 588, A118
  • Malmberg et al. (2011) Malmberg, D., Davies, M. B., & Heggie, D. C. 2011, MNRAS, 411, 859
  • Mann et al. (2016) Mann, A. W., Gaidos, E., Mace, G. N., et al. 2016, ApJ, 818, 46
  • Mardling (2008) Mardling, R. A. 2008, The Cambridge N-Body Lectures, 760, 59
  • McCaughrean & O’dell (1996) McCaughrean, M. J., & O’dell, C. R. 1996, AJ, 111, 1977
  • McMillan et al. (2012) McMillan, S., Portegies Zwart, S., van Elteren, A., & Whitehead, A. 2012, Advances in Computational Astrophysics: Methods, Tools, and Outcome, 453, 129
  • Meibom et al. (2013) Meibom, S., Torres, G., Fressin, F., et al. 2013, Nature, 499, 55
  • Mizuno (1980) Mizuno, H. 1980, Progress of Theoretical Physics, 64, 544
  • Naoz et al. (2011) Naoz, S., Farr, W. M., Lithwick, Y., Rasio, F. A., & Teyssandier, J. 2011, Nature, 473, 187
  • Olczak et al. (2010) Olczak, C., Pfalzner, S., & Eckart, A. 2010, A&A, 509, A63
  • Olczak et al. (2012) Olczak, C., Kaczmarek, T., Harfst, S., Pfalzner, S., & Portegies Zwart, S. 2012, ApJ, 756, 123
  • Parker & Quanz (2012) Parker, R. J., & Quanz, S. P. 2012, MNRAS, 419, 2448
  • Parker et al. (2014) Parker, R. J., Church, R. P., Davies, M. B., & Meyer, M. R. 2014, MNRAS, 437, 946
  • Pelupessy et al. (2013) Pelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, A&A, 557, A84
  • Pfalzner (2013) Pfalzner, S. 2013, A&A, 549, A82
  • Pfalzner et al. (2015) Pfalzner, S., Davies, M. B., Gounelle, M., et al. 2015, Phys. Scr., 90, 068001
  • Plummer (1911) Plummer, H. C. 1911, MNRAS, 71, 460
  • Porras et al. (2003) Porras, A., Christopher, M., Allen, L., et al. 2003, AJ, 126, 1916
  • Portegies Zwart et al. (1998) Portegies Zwart, S. F., Hut, P., Makino, J., & McMillan, S. L. W. 1998, A&A, 337, 363
  • Portegies Zwart (2009) Portegies Zwart, S. F. 2009, ApJ, 696, L13
  • Portegies Zwart et al. (2009) Portegies Zwart, S., McMillan, S., Harfst, S., et al. 2009, New Astronomy, 14, 369
  • Portegies Zwart et al. (2010) Portegies Zwart, S. F., McMillan, S. L. W., & Gieles, M. 2010, ARA&A, 48, 431
  • Portegies Zwart et al. (2013) Portegies Zwart, S. et al. 2013, Computer Physics Communications, 183, 456
  • Portegies Zwart & Jílková (2015) Portegies Zwart, S. F., & Jílková, L. 2015, MNRAS, 451, 144
  • Portegies Zwart (2016) Portegies Zwart, S. F. 2016, MNRAS, 457, 313
  • Press et al. (2002) Press, W. H., Teukolsky, S. A., Vetterling, W. T., & Flannery, B. P. 2002, Numerical recipes in C++ : the art of scientific computing
  • Proszkow & Adams (2009) Proszkow, E.-M., & Adams, F. C. 2009, ApJS, 185, 486
  • Quinn et al. (2012) Quinn, S. N., White, R. J., Latham, D. W., et al. 2012, ApJ, 756, L33
  • Quinn et al. (2014) Quinn, S. N., White, R. J., Latham, D. W., et al. 2014, ApJ, 787, 27
  • Rein & Liu (2012) Rein, H., & Liu, S.-F. 2012, A&A, 537, A128
  • Rein & Spiegel (2015) Rein, H., & Spiegel, D. S. 2015, MNRAS, 446, 1424
  • Renaud et al. (2011) Renaud, F., Gieles, M., & Boily, C. M. 2011, MNRAS, 418, 759
  • Rieder et al. (2013) Rieder, S., Ishiyama, T., Langelaan, P., et al. 2013, MNRAS, 436, 3695
  • Sato et al. (2007) Sato, B., Izumiura, H., Toyota, E., et al. 2007, ApJ, 661, 527
  • Shara et al. (2016) Shara, M. M., Hurley, J. R., & Mardling, R. A. 2016, ApJ, 816, 59
  • Spitzer & Hart (1971) Spitzer, L., Jr., & Hart, M. H. 1971, ApJ, 166, 483
  • Spitzer (1987) Spitzer, L. 1987, Princeton, NJ, Princeton University Press
  • Spurzem (1999) Spurzem, R. 1999, Journal of Computational and Applied Mathematics, 109, 407
  • Spurzem et al. (2009) Spurzem, R., Giersz, M., Heggie, D. C., & Lin, D. N. C. 2009, ApJ, 697, 458
  • Stoer & Bulirsch (1980) Stoer, J., & Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer)
  • Takahashi & Portegies Zwart (2000) Takahashi, K., & Portegies Zwart, S. F. 2000, ApJ, 535, 759
  • Thies et al. (2010) Thies, I., Kroupa, P., Goodwin, S. P., Stamatellos, D., & Whitworth, A. P. 2010, ApJ, 717, 577
  • Thies et al. (2011) Thies, I., Kroupa, P., Goodwin, S. P., Stamatellos, D., & Whitworth, A. P. 2011, MNRAS, 417, 1817
  • Voyatzis et al. (2013) Voyatzis, G., Hadjidemetriou, J. D., Veras, D., & Varvoglis, H. 2013, MNRAS, 430, 3383
  • Wang et al. (2015a) Wang, L., Spurzem, R., Aarseth, S., et al. 2015a, MNRAS, 450, 4070
  • Wang et al. (2015b) Wang, L., Kouwenhoven, M. B. N., Zheng, X., Church, R. P., & Davies, M. B. 2015b, MNRAS, 449, 3543
  • Wang et al. (2016) Wang, L., Spurzem, R., Aarseth, S., et al. 2016, MNRAS, 458, 1450
  • Wu & Lithwick (2011) Wu, Y., & Lithwick, Y. 2011, ApJ, 735, 109
  • Zhang et al. (2013) Zhang, K., Hamilton, D. P., & Matsumura, S. 2013, ApJ, 778, 6
  • Zheng et al. (2015) Zheng, X., Kouwenhoven, M. B. N., & Wang, L. 2015, MNRAS, 453, 2759
  • Zhou et al. (2007) Zhou, J.-L., Lin, D. N. C., & Sun, Y.-S. 2007, ApJ, 666, 423
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description