Tycho: Realistically Simulating Exoplanets within Stellar Clusters. I. Improving the Monte Carlo Approach

Cluster-scale gravitational dynamics is handled by PH4, an MPI-parallel fourth-order Hermite direct N-body integrator with GPU acceleration (McMillan et al. 2011). This AMUSE module is well tested on N-body systems of all sizes, from binaries to clusters containing hundreds of thousands of stars (see Portegies Zwart & McMillan 2018). Like most N-body integrators in AMUSE, PH4 does not include any special treatment of close encounters between stars or binaries. Most AMUSE modules, including PH4, include the possibility of softening the potential, but this is not appropriate for the collisional dynamics of interest here. Instead, Tycho uses the AMUSE Multiples module to manage dynamical close encounters (Pelupessy et al. 2013; Portegies Zwart & McMillan 2018), as we now describe.

Traditional collisional N-body codes devote a substantial fraction of the total code base to modeling the detailed internal motion of multiple systems (where this term henceforth includes stellar binaries) simultaneously with the large-scale motion of the rest of the system (Portegies Zwart et al. 1998; Aarseth 2003). The algorithms involved have varying degrees of efficiency and accuracy, but all gradually relax the accuracy of the perturbed binary integration until "unperturbed" binaries (where the external tidal acceleration is less than ∼10⁻⁶ times any internal acceleration) are treated as isolated unperturbed objects. The alternative approach, widely used in Monte Carlo schemes, treats stable multiples as completely unperturbed until they interact strongly with another object (Fregeau et al. 2004; Tanikawa & Fukushige 2009; Giersz et al. 2013).

As discussed in more detail by Portegies Zwart & McMillan (2018, Section 4.5) the N-body module (here PH4) manages a multiple's center of mass, but is ignorant of its internal dynamics. The Multiples module maintains a database of all stable multiples in the system. Herein, we define that a stable multiple is a binary, a hierarchical triple deemed stable by the Mardling (2007) criterion, or a higher-order multiple system in which the Mardling criterion applied to the outermost orbit(s) indicates approximate stability. Each single or multiple object has an interaction radius R_int assigned to it. For a single star of mass M, R_int ∼ GM/σ², where σ is the cluster velocity dispersion. For a multiple, R_int is twice the outer semimajor axis. For planets, R_int is the Hill radius when bound to a star, and ${R}_{\mathrm{int}}\sim {GM}/{\sigma }^{2}$ when free-floating. When two particles in the N-body code approach within the sum of their interaction radii, the N-body integrator is advanced to and paused at that moment. The particles are then removed from the N-body system, their internal structure is restored, and the interaction is followed as an isolated few-body system using the SmallN integrator, a robust shared-step, time-symmetrized code with analytic extensions (Hut et al. 1995; McMillan & Hut 1996; Portegies Zwart & McMillan 2018). Pure two-body encounters are advanced analytically using the AMUSE Kepler module.

In order to inform the high-fidelity scattering experiments detailed in Section 2.4, we store the encounter history of each star throughout the course of the simulation. At the start of a close encounter, Tycho copies the flattened hierarchical tree of the root particles involved into an AMUSE particle set. This particle set and all required initial conditions are then stored in a chronological database of encounter histories for each star in the simulation. Once the N-body evolution ends, a complete history of all stellar close encounters is available within the database.

The structure of the few-body system is monitored during the SmallN integration and the integration is declared over when the system consists of some number of mutually unbound, receding, stable (as just defined) centers of mass. For stars with planetary systems, the "over" criterion is applied only to the stars. At that point, the internal structure of the newly obtained multiples is stored for future use, corrections are applied to partially account for the tidal field of the rest of the N-body system, the centers of mass are added to the N-body system, and the large-scale integration continues. The internal structure of each center of mass particle is held static until the next interaction occurs. Examples of this hierarchical structure can be found in Figure 2.

Zoom In Zoom Out Reset image size

For our models, all primary stellar masses are drawn from a Kroupa (2001) mass function, as discussed in detail in Section 3.1. In addition, we include primordial binaries and allow dynamical binaries to form naturally, to ensure an accurate depiction of large-scale cluster dynamics. Tycho uses the SeBa stellar evolution module (with solar metallicity), with an additional prescription for mass loss due to winds (Portegies Zwart & Verbunt 1996, 2012; Toonen et al. 2012).

To account for multiple age groups in our stellar population, we implement a procedure for asynchronous evolution within SeBa. This methodology will become increasingly more important for accurate depictions of dynamics as direct simulations of cluster formation from molecular clouds become more widespread (e.g., Wall 2018; Wall et al. 2019). In order to accomplish this, we bin our stars by their age, which allows us to create separate instances of SeBa for each grouping. These separate instances within the AMUSE environment allow for the module to accurately evolve while reducing the computational overhead to a manageable level. See Figure 3 for an example of a Hertzsprung–Russell diagram drawn from a roughly 500 M_⊙ cluster.

Zoom In Zoom Out Reset image size

The evolution of a cluster within the tidal field of its host galaxy is substantially different from one in isolation. The tidal field imposes a limiting Jacobi radius on the cluster, at which stars become unbound (von Hoerner 1957; King 1962). As mass is lost from the cluster, the radius shrinks, leading to a decreased cluster lifetime (Baumgardt & Makino 2003; Gieles & Baumgardt 2008). In order to accurately model the escape of stars from the cluster into the Galactic field, we apply a background gravitational potential which models key aspects of the the Milky Way. As described in detail in Bovy (2015), this static potential features a central bulge, a dusty disk and a dark matter halo modeled after current Milky Way observations. To couple the Galactic potential with the time-dependent evolution of the cluster, we utilize the AMUSE framework's Bridge package, a code-coupling algorithm based on the scheme described by Fujii et al. (2007). A cluster may be placed in any kind of orbit in this potential and the relevant tidal forces will be fully taken into effect.

Of the 4126 confirmed exoplanets to date, 1766 exist within 705 multi-planetary systems.⁶ Accordingly, it is critical to include within Tycho the ability to faithfully model multi-planet systems. In our model, we record the parameters of all close encounters between stellar systems (single or hierarchical) in the cluster and use these parameters to later perform a series of scattering experiments, the procedure for which we will now describe.

Once the full-cluster N-body simulations are completed, we re-simulate our independent encounters in parallel at higher temporal resolution with additional planetary bodies drawn from our mock solar system model described in Section 3.2. We again utilize the specialized few-body integrator SmallN for these scattering experiments. To improve run-time efficiency, we advance the two stellar systems involved in the encounter along their Keplerian orbit to the point, r₁₂, where the gravitational perturbation on the outermost orbit is of order 10⁻³:

$$\begin{eqnarray}&&{r}_{12}\sim 10{a}_{\mathrm{outer}}\displaystyle \frac{{M}_{\mathrm{pert}}}{{M}_{\mathrm{host}}}.\end{eqnarray} \tag{ 1 }$$

Here a_outer is the semimajor axis of the outermost planet, M_pert is the total mass of the perturbing system, and M_host is the mass of the star hosting the planetary system. We then allow SmallN to simulate the close encounter utilizing an adaptive internal time step, storing particle position and velocity vectors every 1 yr until the close encounter is declared over, as described in Section 2.1. We then evolve the system for an additional 100 yr before committing and tabulating the planets' orbital parameters.

To ensure that our dynamical parameters are fully sampled, we simulate each encounter in our database with the host star's planetary system multiple times, with different orientations and mean anomalies. At the start of each scattering experiment instance, we apply a random rotation transformation which samples the unit sphere uniformly to the system using the fast algorithm developed by Arvo (1992). The proof of this method's uniform coverage is outlined by Shoemake (1992). To develop a statistically significant population of encounters from our N-body database, we repeat this procedure 100 times.

As we are interested in the evolution of our systems over their time spent within the parent cluster, we need a quantifiable system for analyzing not only the shifts in a given planet's orbits, but also the overall stability of the system they inhabit. Measuring the angular momentum deficit (AMD) in our systems allows for such a metric. Following Laskar (1997), the ratio of a planet's relative AMD to that of its system as a whole can be interpreted as a secular measure of the perturbation of the orbit. We can then find the regime in which the possibility of strong encounters with other bodies in the planetary system are forbidden, implying long-term stability.

Due to its conservative nature, many different stability criteria can be included into the AMD framework. For our purposes, we specifically include the criteria imposed on the critical AMD (the value at which a given planet becomes unstable) stemming from Hill stability (Petit et al. 2018), orbit-crossing collisions (Laskar & Petit 2017), and first-order mean motion resonance (MMR) overlap (Petit et al. 2017). We present a complete summary of the AMD framework in the Appendix.

In order to allow for an efficient and complete classification of our planetary systems, we adopt the following prescription: if the jth planet's AMD-stability coefficient ${\beta }_{j}={{ \mathcal C }}_{j}/{C}_{c}^{{\rm{H}}}\lt 1$, we do not calculate C_c; if the reverse is true, then we set ${\beta }_{j}={{ \mathcal C }}_{j}/{C}_{c}$. As noted in Section 4 of Petit et al. (2018), this ensures the system's long-term stability or lack thereof is correctly validated. Using β_j as a post-processing indicator, we are able to correctly categorize our systems as stable (where all planets have β_j < 1), unstable (where two or more planets have β_j > 1), or metastable (where the innermost planet is unstable but all other planets are stable).

We draw our clusters' spatial density distributions from King (1962, 1965, 1966) models, which provide rudimentary, if imperfect, fits to observed stellar clusters with masses of a few times 10³M_⊙ (Portegies Zwart et al. 2010). Detailed observations of star-forming regions show stellar distributions which tend to be clumpier in composition than King models (Sánchez & Alfaro 2009; Allison et al. 2010), with stars still embedded in their natal gas for the first few megayears. N-body models of hierarchical or fractal stellar clusters relax to King-like spatial distributions in a short time (Geller et al. 2013; Smith et al. 2013). We synthetically age our stellar population and begin our models at 10 Myr, to avoid this gas-embedded, spatially clumpy phase. In the future, we will couple our models to detailed star formation simulations (such as those presented in Wall et al. 2019) to better capture this initial phase. For convenience, we initialize all of our clusters on circular orbits of Galactocentric radii 9 kpc.⁷

For our stellar masses, we sample a truncated Kroupa initial mass function (IMF; Kroupa 2001). We set our lower and upper stellar mass limits as the lowest mass for an observed red dwarf (M_min = 0.1M_⊙; Dieterich et al. 2014) and the approximate maximum single-component mass following complete gas ejection at ∼5–10 Myr (M_max = 10M_⊙; Lada 2010), respectively. Limiting our lower mass allows for us to more efficiently spend computational resources integrating the cluster's evolution, without significantly altering the overall evolution of the cluster (Kouwenhoven et al. 2014).

We draw the orbital parameters of our primordial stellar binary population from empirical distributions derived from observational surveys. The creation procedure implemented is as follows.

• 1.  
  To determine if a given primary star should have a stellar binary companion, we adopt the empirical binary fractions as a function of primary mass from Figure 12 of Raghavan et al. (2010). If a star is selected to have a secondary component, its position becomes the center of mass of the new binary.
• 2.  
  The secondary mass is selected from the truncated Kroupa IMF and added to the initial mass of the center of mass particle. This allows us to recover the observed IMF when we redistribute the mass between the primary and secondary star.
• 3.  
  Using a uniform mass-ratio distribution (Raghavan et al. 2010; Goodwin 2013), we split the the center of mass particle's new total mass into the primary and secondary star.
• 4.  
  To determine the stellar binary's orbital elements, we use the primary's mass to draw from empirically based distributions. The period is generated from the log-normal distribution of Raghavan et al. (2010) for lower-mass and Sana et al. (2012) for higher-mass systems. Eccentricities are drawn from a uniform distribution well documented in the literature (Raghavan et al. 2010; Duchêne & Kraus 2013; Moe & Di Stefano 2017).
• 5.  
  The orientation of the system with respect to the cluster's coordinate system is chosen randomly and uniformly over the unit sphere as described in Arvo (1992).

Figure 4 provides an overview of the yield of the above-described procedure. Once generated, binaries are added to the N-body particle set to be picked up by Multiples before the cluster is scaled to virial equilibrium. All internal orbital interactions, close encounters, and stellar evolution are tracked and stored by Tycho's bookkeeping procedures for future analysis. There are no initial systems consisting of more than two bodies in our simulated clusters.

Zoom In Zoom Out Reset image size

Within the initial cluster models presented in the previous section, we model the planetary systems after our own solar system—using a stable, well-studied system allows us to differentiate the effects of star–star scattering on multi-planet systems from internal interplanetary dynamics. Thus we adopt a scaled version of the solar system, which contains only analogs for Earth, Jupiter, and Neptune. We refer to these as Terrestrial, Jovian, and Neptunian, respectively. We initialize our systems as coplanar, consistent with leading formation models (Laskar 2000).

In accordance with current formation theory, the placement of each system's Jovian is scaled to the snow line of its host star (Kennedy & Kenyon 2008). Thus, the Jovian planets have semimajor axes ${a}_{\mathrm{jovian}}={a}_{{\rm{J}}}{\left(\tfrac{{M}_{\mathrm{host}}}{{M}_{\odot }}\right)}^{2}$ (Ida & Lin 2005). After scaling the Jovian's semimajor axis, we scale the other planetary orbits so that the system remains long-term stable. By fixing the ratio of periods between each planet and the system's Jovian, we ensure that the solar system's innate stability is preserved regardless of the mass of the host star. A good measure of this is the AMD of the planets (see Section 2.5). Table 3 summarizes our model's parameters and provides the AMD stability coefficient of each planet, averaged over the host star masses in our simulations.

Table 3.  Implemented Planetary Systems

	Mass (MJ)	Eccentricity	Semi-major Axis (au)	Avg. AMD Stability Coefficient, ${\beta }_{\mathrm{AMD}}$
Terrestrial	0.003	0.016	$1.000\times {\left(\tfrac{{M}_{\mathrm{host}}}{{M}_{\odot }}\right)}^{2}$	$(0.9015\pm 4.50)\times {10}^{-4}$
Jovian	1	0.048	$5.454\times {\left(\tfrac{{M}_{\mathrm{host}}}{{M}_{\odot }}\right)}^{2}$	(0.2394 ± 2.35) × 10−4
Neptunian	0.054	0.009	$30.110\times {\left(\tfrac{{M}_{\mathrm{host}}}{{M}_{\odot }}\right)}^{2}$	(0.02361 ± 1.79) × 10−5

Download table as:  ASCIITypeset image

To illustrate our improved approach for simulating planetary systems within clusters, we simulate 40 independently seeded open clusters with initial conditions as described above. These consist of 10 realizations for each of the possible parameter combinations of the number of center-of-mass objects (N ∈ {100, 1000}) and the depth of the King potential (W₀ ∈ {3.0, 6.0}).⁸ Each realization is simulated according to our prescription for approximately 2 Gyr, which is past the point of cluster dissolution in most cases. Planetary systems are generated only around single stars and include only the Jovians in our full-cluster N-body simulations (as discussed above), but we retain the full planetary systems in our scattering experiments. For the purposes of showcasing the dynamics, we limit our close encounter database to those with r_periapsis < 2a_outer.

Our simulations were performed on the Draco supercomputing cluster at Drexel University. Draco consists of 24 compute nodes, each consisting of two Intel Xeon x5650 CPU six-core chips and four to six vintage (Tesla/Titan) NVIDIA GPUs, allowing us to take advantage of PH4's parallel and cuda-enabled design. Our N-body simulations took 553 CPU hours to complete, and our scattering experiments took 1690 CPU hours, for a total of  1.33 weeks distributed across 10 nodes.

A catalog of our simulations has been provided to the publisher, containing a total of 41,313 planetary systems of various configurations. Figure 5 shows the relation between the total change in eccentricity and the total change in semimajor axis for our three planetary populations after cluster dissolution. Most of our stellar population only encountered another stellar system once over the 2 Gyr evolution. Given the much greater scattering cross section of our Neptunian population, our results confirm the expectation that they are far more significantly affected by stellar flybys than the Jovian or Terrestrial planets. Indeed, of the 371 planets that escaped their host system across all encounters, 222 were Neptunians, 77 were Jovians, and 61 were Terrestrials. Further, of these planets that left their systems, the Neptunians were more likely to be captured (26 Neptunians versus fove Jovians and no Terrestrials).

Zoom In Zoom Out Reset image size

We note that most of our systems do not have large changes in eccentricity or semimajor axis even though all of the systems included in our catalog have undergone at least one stellar encounter. Those systems that do experience significant changes in their semimajor axes usually experience changes in their eccentricities as well. We see that more of our systems that retain three planets after the encounter (as plotted in Figures 5–9) experience hardening (i.e., shrinking of the semimajor axis) due to stellar encounters. Apparently, encounters that soften the orbits are either less common, or are more likely to lead to the disruption of the planetary system.

Furthermore, the eccentricity histogram of the left graph of Figure 5 shows that as the system shrinks, our terrestrial population experiences large changes in eccentricity while staying near the original semimajor axis.

Figure 6 compares the change in eccentricity in each of our planetary types drawn from systems that retain all three original planets. The color indicates the long-term stability of the system as determined by the AMD framework. We see that, while all of our systems begin as AMD stable, 37.5% become weakly unstable and 0.81% become fully unstable after 2 Gyr. Of this unstable population, a significant portion of outer planets meet the criteria for eventual collision with an inner body rather than instability generated via the MMR mechanism, with no system becoming Hill unstable. As expected, there is a high concentration of unstable systems in which Jovians and Neptunians both experience large changes in eccentricity. However, very few Jovians reach a point where their AMD stability coefficients exceed unity. In weakly stable systems, most arise from instability in their companion Terrestrial that results in orbital decay and inspiral into the parent star.

Zoom In Zoom Out Reset image size

To provide an example from the catalog, we have included a very unstable system in Figure 6. This system is the result of a binary stellar system (M_b1 = 0.133 M_⊙ and M_b2 = 0.369 M_⊙; period  = 9746 days, a = 7.105 au, and e = 0.248) encountering the host star (${M}_{\star }=0.950\,{M}_{\odot }$). The encounter has a closest approach of 10.001 au over its 35.6 kyr lifetime (a = 1224.27 au and e = 1.008), placing the binary well within the planetary system. A summary of the final orbital elements of the system can be found in Table 4. As a result of this encounter, the planetary system is predicted to result in an eventual collision between the Jovian and Neptunian while the Terrestrial inspirals toward the host star.

Table 4.  Resulting System around Star Number 339

	Eccentricity	Semi-major Axis (au)	Relative Inclination	βAMD
Terrestrial	0.031 (+0.021)	0.898 (+0.005)	−18137	570.68
Jovian	0.752 (+0.704)	4.687 (−0.232)	⋯	3.498
Neptunian	0.99 (+0.981)	363.889 (+336.726)	+23570	1.878

Note. This system can be found in the catalog provided to the publisher under the System Key: Adam_N1000_W6/S1340/Enc-0_Rot-63. The table includes the final value of each parameter and the difference from the initial condition, indicated in parenthesis. We take our inclination measurements from the orbital plane of the most massive planet, which in these proof-of-concept simulations is the system's Jovian.

Download table as:  ASCIITypeset image

Finally, we turn to the relationship between changes in the orbits and the mass of the system's stellar host. As the stellar mass increases so does the change in the semimajor axis, with the Neptunian population strongly following this trend (as seen in Figure 7). Meanwhile, planets experience more dispersion of their eccentricity at lower stellar host masses. As expected, the relative inclinations of the planetary systems resemble that of a Gaussian distribution centered about the co-planar starting condition (σ_T = 47 and σ_N = 59; see Figure 8). Additionally, Figure 9 shows that planets orbiting lower-mass stars have the highest amount of dispersion in their AMD stability coefficient. Note that medians of the three planetary types stay nearly constant across the range of host masses (centered around ${\beta }_{{\rm{T}}}$ = 0.90, β_J = 0.24, and β_N = 0.02).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The results from these simulations are revealing; they predict that a significant portion of planets within gas-ejected open clusters experience long-term stability changes due to stellar close encounters. From planetary ejections to inner-planet inspiral, it is clear that encounters fitting the parameters described in Sections 2.4 and 4.1 can drive dynamical changes in the structure of the planetary system.

As previously noted, we see that systems that become unstable always have a Neptune with a significant eccentricity change. This is expected, as when stellar systems encounter one another, the outer planets will have the largest orbital perturbation. By increasing the orbital eccentricities, these changes lower the outer planets' periapses, resulting in further orbital shifts via planet–planet interactions and, eventually, long-term instability.

We can understand the trends in our planetary systems' orbital changes through the lens of dynamical mass segregation (McMillan et al. 2007). As an open cluster evolves, stellar encounters redistribute energy throughout the cluster, attempting to drive the stars toward energy equipartition. This leads the more massive stars to move more slowly and hence sink toward the center of the cluster. As a result, more massive stars are more likely to encounter other massive stars, leading to larger encounter energies and thus greater perturbations to a planet's semimajor axis. Additionally, lower-mass stars pass through this dense core before being ejected from the cluster via evaporation, often after a single strong encounter. During such strong encounters (which our cuts highlight), even tightly packed planetary systems around these low-mass hosts can experience dramatic orbital perturbations.

Regarding the significant population of unstable Terrestrials, we note that these planets are near the critical point of β = 1 before any encounters. As such, only a small change in the angular momentum of the system is required to push this category of planet toward an instability criterion either via inspiral into its stellar host or a collision in the densely packed system. We find that while the median of our Terrestrial population remains at β = 0.945, they exceed β = 1 at the 61.63 percentile. Conversely, the Jovian and Neptunian populations only exceed β = 1 at the 99.5 percentile despite their tendency toward significant orbital change. Indeed, Dotti et al. (2019) saw similar behaviors in the stability of Terrestrials, where perturbations to the system's Jovian easily caused instability of the inner planets.

In future papers we will improve and expand our initial conditions to better match observed clusters, such as those presented in Table 2, and implement new numerical tools to enhance Tycho's capabilities in order to make more quantitative predictions about the expected observational sample of exoplanets.

As mentioned in Section 2.1, our current methodology does not actively simulate planetary systems between stellar encounters. While this does not affect the validity of our results thanks to the robustness of the AMD stability criterion, integrating a system over long timescales will expose important nuances in its evolution. Recently, a generalized hierarchical approach to the long-term secular integration of planetary systems was described by Hamers & Portegies Zwart (2016). It was further expanded to include a secular treatment of close encounters between such hierarchical systems and is now included in AMUSE (see the SecularMultiples community code; Hamers 2018). To improve both Tycho's physical accuracy and computational efficiency, we plan to integrate this approach into our scattering experiments in the next iteration.

Additionally, Tycho does not follow planetary systems through multiple stellar encounters. Rather, it treats encounters as wholly separate from one another, regardless of a system's history or future. Most stars in these clusters only encounter one other star during the cluster lifetime and thus Tycho's approach is appropriate for the simulations presented here. However, in future larger simulations, the number of stars that undergo several encounters in their lifetime, and particularly those that migrate to the center of the cluster during the core-collapse epoch, will not be insignificant. The implementation of SecularMultiples will allow planetary systems that undergo multiple encounters to be efficiently evolved between encounters, and therefore enable proper handling of such encounter chains.

In this paper, we present a robust methodology, Tycho, that aims to solve the multifaceted problem of evolving planetary systems from their birth in their natal clusters to their eventual dispersion into the Galactic disk. Our proof-of-concept simulations verify that Tycho can simulate such clusters within a reasonable amount of time without compromising the physical integrity required for such work.

Examining the results of our proof-of-concept runs, we find that a significant portion (37.5%) of our remaining solar system analogs become weakly unstable due to close encounters with other stellar systems. A smaller portion become fully unstable, mainly due to the collisional criterion with the inner planet. Additionally, while systems orbiting high-mass stars experience larger changes in their semimajor axes than those around lower-mass stars, the latter experience larger dispersion in β_j. Finally, 0.8% of the planets in our population were ejected from their host system and/or captured by another stellar system due to stellar close encounters.

In summary, we have developed a methodology that allows us to explore how stellar birth environments influence the evolution of planetary systems. Building upon previous numerical solutions, Tycho strikes a balance between computational efficiency of scattering experiments and direct N-body simulations. Furthermore, we plan to continue developing Tycho to enable direct comparisons between simulated and observational planetary populations, the progress of which can be followed on the project's GitHub repository (https://github.com/JPGlaser/Tycho).