Survivability of planetary systems in young and dense star clusters

TL;DR

This study uses simulations of the Orion Trapezium star cluster to investigate how planetary systems evolve and produce free-floating planets, revealing that planetary loss mechanisms are largely independent of planet mass and initial orbit.

Contribution

It provides the first detailed simulation-based analysis of planetary survivability and free-floating planet production in a young, dense star cluster environment.

Findings

357 free-floating planets formed, with 281 leaving the cluster

Planet loss driven mainly by stellar encounters and planetary scattering

Mass distribution of free-floating planets matches that of bound planets

Abstract

We perform a simulation using the Astrophysical Multipurpose Software Environment of the Orion Trapezium star cluster in which the evolution of the stars and the dynamics of planetary systems are taken into account. The initial conditions from earlier simulations were selected in which the size and mass distributions of the observed circumstellar disks in this cluster are satisfactorily reproduced. Four, five, or size planets per star were introduced in orbit around the 500 solar-like stars with a maximum orbital separation of 400au. Our study focuses on the production of free-floating planets. A total of 357 become unbound from a total of 2522 planets in the initial conditions of the simulation. Of these, 281 leave the cluster within the crossing timescale of the star cluster; the others remain bound to the cluster as free-floating intra-cluster planets. Five of these free-floating intra-cluster planets are captured at a later time by another star. The two main mechanisms by which planets are lost from their host star, ejection upon a strong encounter with another star or internal planetary scattering, drive the evaporation independent of planet mass of orbital separation at birth. The effect of small perturbations due to slow changes in the cluster potential are important for the evolution of planetary systems. In addition, the probability of a star to lose a planet is independent of the planet mass and independent of its initial orbital separation. As a consequence, the mass distribution of free-floating planets indistinguishable from the mass distribution of planets bound to their host star.

Figures

Tables

Table 1. Table 1: Initial cluster model adopted by Portegies Zwart ( 2016 ) ; the final conditions for the disk-size analysis in Portegies Zwart ( 2016 ) , which we adopted as the initial realization for the simulations presented here; and the final conditions. The parameter N total subscript 𝑁 total N_{\rm total} indicates the total number of stars in the simulation; N bnd subscript 𝑁 bnd N_{\rm bnd} , the bound mass of the cluster in solar masses; R vir subscript 𝑅 vir R_{\rm vir} the virial radius of the cluster; Q 𝑄 Q , the virial equilibrium; F 𝐹 F , the fractal dimension of the cluster; N bnd subscript 𝑁 bnd N_{\rm bnd} , the number of bound stars; N bnd subscript 𝑁 bnd N_{\rm bnd} , the number of bound stars with planets; N bnd subscript 𝑁 bnd N_{\rm bnd} , the number of unbound stars with planets; N w / planets subscript 𝑁 w planets N_{\rm w/planets} , the number of stars with planets; N r ≥ 100 ​ a ​ u subscript 𝑁 𝑟 100 𝑎 𝑢 N_{r\geq 100au} the number of stars with disks or planets equal or larger than 100 au; N r ≥ 100 ​ a ​ u subscript 𝑁 𝑟 100 𝑎 𝑢 N_{r\geq 100au} the number of stars with disks or planets equal or larger than 10 au; n bnd subscript 𝑛 bnd n_{\rm bnd} , the number of planets bound to a star; n ff subscript 𝑛 ff n_{\rm ff} , the number of free-floating planets; n unbnd subscript 𝑛 unbnd n_{\rm unbnd} , the number of planets unbound from the cluster; m bnd subscript 𝑚 bnd m_{\rm bnd} , the total planetary mass in Jupiter masses, bound to a star; m ff subscript 𝑚 ff m_{\rm ff} , the total planetary mass in free-floating planets; and m unbnd subscript 𝑚 unbnd m_{\rm unbnd} , the total planetary mass unbound from the cluster.

Parameter	$t=0$ Myr	$t=1$ Myr	$t=11$ Myr
Cluster characteristics
$N_{\mathrm{total}}$	1500	1500	1482
$M_{\mathrm{bnd}}$/$M_{\odot }$	627	618	545
$R_{\mathrm{vir}}$/pc	0.5	0.36	0.32
Q	1.0	0.6	1.0
F	1.6	1.26	0.6
Stellar characteristics
$N_{\mathrm{bnd}}$	1500	977	508
$N_{\mathrm{bnd},w/p}$	0	512	166
$N_{\mathrm{unbnd},w/p}$	0	0	323
$N_{w/\mathrm{planets}}$	0	500	517
$N_{r\ge 100au}$	1500	78
$N_{r\ge 10au}$	1500	578
Planets characteristics
$n_{\mathrm{bnd}}$	–	2522	2165
$n_{\mathrm{ff}}$	–	0	357
$n_{\mathrm{unbnd}}$	–	0	282
$m_{\mathrm{bnd}}$/$M_{\mathrm{Jup}}$	–	3527	2915
$m_{\mathrm{ff}}$/$M_{\mathrm{Jup}}$	–	0	502
$m_{\mathrm{unbnd}}$/$M_{\mathrm{Jup}}$	–	0	395

Table 2. Table 2: Comparison of the distribution of planets at the beginning and at the end of the simulation. In each cell, the count of the number systems with a certain number of planets is given. This count is given per number of planets in the original system. The original distribution 1  M ​ y ​ r 𝑀 𝑦 𝑟 Myr is given in the top summation row; the final distribution at 11  M ​ y ​ r 𝑀 𝑦 𝑟 Myr is given in the last column. For the final distribution of the 6 systems with 3 planets, 2 of these systems originally had 4 planets, 4 originally had 5 planets, and none originally had 6 planets. A total of 5 new planetary systems have been created during the evolution of the cluster, in these systems originally the star had no planets. Of these 5 new planetary system, 3 systems have 1 planet and 2 have gained 2 planets.

	$N_P$
	0	1	2	3	4	5	6	$\sum$
$\sum$	5	0	0	0	109	332	71
0	0	0	0	0	3	25	0	28
1	3	0	0	0	6	35	10	54
2	2	0	0	0	5	17	6	30
3	0	0	0	0	2	4	0	6
4	0	0	0	0	93	11	0	104
5	0	0	0	0	0	240	2	242
6	0	0	0	0	0	0	53	53

Table 3. Table 3: Parameters for planetary systems in which one planet was ejected to a larger distance ( > 800 absent 800 >800  au) from its host star. The second column identified which planet was ejected, followed by the mass of the host, planet mass, and its eventual orbital parameters.

System	Planet	M $\left(M_{sun}\right)$	m $\left(M_{jup}\right)$	a (au)	e
0	a	0.33	1.27	983.2	0.93
1	a	0.16	2.23	900.1	0.89
2	a	0.37	0.85	1107	0.39
3	b	0.37	6.05	1311	0.71
4	a	0.25	1.39	933.6	0.94
5	a	0.20	2.43	1062	0.94
6	a	0.72	6.85	1693	0.90

Table 4. Table 4: Listing of systems that formed by the abduction of a planet from another star. Each of these stars was initially without any planets, but one or two planets were captured from another system. In two cases (#7 and #10) two planets were captured.

System	Planet	M $\left(M_{sun}\right)$	m $\left(M_{jup}\right)$	a (au)	e
7	a	9.16	0.28	1544	0.68
7	b	9.16	0.50	1161	0.22
8	a	6.70	14.69	3332	0.50
9	a	3.61	0.95	891.5	0.77
10	a	0.47	0.03	1100	0.66
10	b	0.47	0.12	1208	0.71
11	a	0.55	0.38	192.5	0.60

Table 5. Table 5: Orbital elements of the merging planets. For each merger the time of the snapshot saved just before the merger is given. For every planet the index of the planetary system is given with a letter denoting the position of the planet in the system (from the innermost planet ’a’ to the outermost planet ’f’). We define the inclination of a planet with respect to the initial orbital plane of the closest planet to the star.

Time (Myr)	id	M ($M_{jup}$)	a (au)	e	i (∘)	id	M ($M_{jup}$)	a (au)	e	i (∘)
3.10	1e	0.30	213.7	0.68	14.4	1d	0.15	105.5	0.18	125.0
3.14	1e	0.46	78.3	0.70	31.8	1c	0.08	33.4	0.38	100.2
3.28	1e	0.54	37.4	0.39	155.2	1b	0.05	13.5	0.29	87.5
3.85	2e	1.30	150.5	0.31	32.4	2d	0.62	118.2	0.21	37.2
4.24	2e	1.93	129.0	0.18	-31.1	2b	0.19	92.1	0.36	37.3
4.96	3e	1.03	234.0	0.50	4.0	3d	0.51	90.9	0.22	7.3
5.13	4f	0.76	401.5	0.51	-38.3	4c	0.13	130.5	0.35	-6.1
5.31	5d	0.58	130.3	0.66	-4.5	5c	0.25	115.6	0.20	-9.4
5.36	6e	0.49	160.6	0.57	7.3	6d	0.24	87.7	0.29	18.9
5.49	7e	1.12	222.3	0.69	38.7	7d	0.56	161.0	0.54	-5.0
5.83	8d	3.93	295.4	0.54	0.0	8a	0.39	96.8	0.99	7.1
5.85	8d	4.31	255.7	0.56	6.1	8c	1.61	296.1	0.52	5.2
5.91	9e	0.37	223.4	0.38	0.0	9a	0.04	215.0	0.84	-10.2
5.98	10e	0.30	101.7	0.39	18.2	10d	0.17	133.3	0.42	-10.5
6.09	5d	0.82	121.9	0.57	-26.1	5b	0.12	44.8	0.63	5.9
6.20	11e	0.81	111.1	0.62	-12.6	11d	0.39	44.4	0.44	-13.9
7.55	12d	0.13	128.1	0.19	-8.9	12b	0.04	98.4	0.60	5.1
7.71	13f	0.17	214.9	0.61	46.2	13e	0.09	42.6	0.92	24.9
8.25	14f	2.29	301.0	0.15	1.1	14e	1.20	127.8	0.01	-1.9

Table 6. Table 6: Orbital elements of the planets resulting from a merger.

Time (Myr)	id a	id b	M $\left(M_{jup}\right)$	a (au)	e	i (∘)
3.10	1e	1d	0.46	78.8	0.70	100.5
3.14	1e	1c	0.54	37.2	0.40	93.6
3.28	1e	1b	0.59	19.0	0.14	92.7
3.85	2e	2d	1.93	129.9	0.23	34.6
4.24	2e	2b	2.12	102.9	0.08	32.6
4.96	3e	3d	1.54	137.5	0.21	5.7
5.13	4f	4c	0.89	256.2	0.26	-9.5
5.31	5d	5c	0.82	105.1	0.55	-7.3
5.36	6e	6d	0.73	105.7	0.27	13.0
5.49	7e	7d	1.69	172.4	0.64	14.4
5.83	8d	8a	4.31	236.3	0.40	5.4
5.85	8d	8c	5.92	220.9	0.50	5.4
5.91	9e	9a	0.41	166.6	0.24	-10.1
5.98	10e	10d	0.46	98.7	0.31	1.4
6.10	5d	5b	0.94	94.8	0.54	-8.0
6.20	11e	11d	1.21	63.2	0.48	-5.4
7.55	12d	12b	0.16	107.4	0.06	-6.3
7.71	13f	13e	0.26	73.0	0.14	27.0
8.25	14f	14e	3.50	255.7	0.10	-1.3

Full text

*Methods. *The initial conditions from earlier simulations were selected in which the size and mass distributions of the observed circumstellar disks in this cluster are satisfactorily reproduced. Four, five, or size planets per star were introduced in orbit around the 500 solar-like stars with a maximum orbital separation of 400 au.

*Results. * Our study focuses on the production of free-floating planets. A total of 357 become unbound from a total of 2522 planets in the initial conditions of the simulation. Of these, 281 leave the cluster within the crossing timescale of the star cluster; the others remain bound to the cluster as free-floating intra-cluster planets. Five of these free-floating intra-cluster planets are captured at a later time by another star.

*Conclusions. * The two main mechanisms by which planets are lost from their host star, ejection upon a strong encounter with another star or internal planetary scattering, drive the evaporation independent of planet mass of orbital separation at birth. The effect of small perturbations due to slow changes in the cluster potential are important for the evolution of planetary systems. In addition, the probability of a star to lose a planet is independent of the planet mass and independent of its initial orbital separation. As a consequence, the mass distribution of free-floating planets is indistinguishable from the mass distribution of planets bound to their host star.

1 Introduction

In recent years several free-floating planets, i.e., planets not orbiting a star, have been discovered by direct infrared imaging (Pacucci et al., 2013) and bycatch in gravitational microlensing surveys (Sumi et al., 2011; Gaudi, 2012; Gould & Yee, 2013). Following star formation theory planets could in principle form in isolation (Gahm et al., 2007; Liu et al., 2013; Haworth et al., 2015), but it seems more likely that they form according to the canonical coagulation process in a disk orbiting a host star (Kant, 1755). If planets are not formed in isolation, there are three major mechanisms by which planets can be liberated. A planet may become unbound as a result of (i) dynamical interaction with another star (Hurley & Shara, 2002; Vorobyov et al., 2017; Cai et al., 2017, 2018; Zheng et al., 2015), (ii) scattering interactions among the planets in a multi-planet system (Veras & Raymond, 2012; Cai et al., 2017, 2018), (iii) copious mass loss in a post-AGB phase (Veras et al., 2015; Veras, 2016) or supernova explosion of the host star (Blaauw, 1961), and (iv) the ejection of fragments when the protoplanetary disk is perturbed (Vorobyov et al., 2017). The relative importance of each of these and other possible processes are hard to assess, but the four listed here are probably most common.

A total of 20 free-floating planet candidates have been identified (Udalski et al., 2008; Wright et al., 2010; Winn & Fabrycky, 2015; Mroz et al., 2018). Two of these orbit each other in the binary-planet 2MASS J11193254-1137466 (Best et al., 2017), but all others are single. Weak micro-lensing searches indicate that the number of free-floating planets with masses exceeding that of Jupiter is about one-quarter of the number of main-sequence stars in the Milky Way Galaxy, whereas Jupiter-mass planets appear to be twice as common as main-sequence stars (Sumi et al., 2011). Interestingly, Earth-mass free floaters are estimated to be only comparable in number to main-sequence stars (Cassan et al., 2012); there appears to be a peak in the number of free-floating planets around the mass of Jupiter.

If rogue planets are liberated upon a strong encounter with another star in a cluster, this process is likely to take place during its early evolution after circumstellar disks have coagulated into planets and most of the primordial gas has been lost. By this time, the stellar density is still sufficiently high that strong encounters between stars are common (Portegies Zwart & Jílková, 2015). Young star clusters may, therefore, make an important contribution to the production of free-floating planets. However, this is at odds with the low number of free-floating planets seen in star clusters. Only one rogue planet was found in the TW Hydra association (Schneider et al., 2016) and a dozen candidates were found in the sigma Orionis cluster (Zapatero Osorio et al., 2013), but no planets were found in the Pleiades cluster despite active searches (Zapatero Osorio et al., 2014). These estimates are in sharp contrast to the number of asteroids and other sōlī lapidēs111sōlus lapis, means “lonely rock” in Latin. expected from the star formation processes (Portegies Zwart et al., 2018b).

The majority of free-floating planets appear as part of the field population, but this may be a selection effect of the methods used to find them (Winn & Fabrycky, 2015). To some degree, however, their relatively high abundance in the field does not come as a surprise. If every star that turns into a white dwarf liberates its planets (and other debris), the number of isolated free floaters should exceed the number of white dwarfs at least by the average number of planets per star. Many of these stars are then already part of the field population once they turn into white dwarfs, giving a natural reduction of free-floating planets in clusters compared to the field population. However, this would mean that dynamical interactions and internal planetary instabilities have a minor contribution to the formation of free-floating planets.

In order to investigate the consequences of stellar evolution and dynamical interactions on the production of free-floating planets, we perform a series of calculations in which we take the relevant processes into account. The main question we address is to what degree the dynamics of a star cluster contribute to the formation and variety of free-floating planets, and what is the relative importance of the various channels for producing these planets.

Planetary systems in our simulation are born stable in the sense that allowing the systems to evolve in isolation would not result in dynamical interactions among the planets. This enables us to study specifically the relative contribution of dynamical interactions on the production of free-floating planets. The stars in our simulations that receive a planetary system are selected such that they remain on the main sequence for the entire duration of the simulation. Stellar mass loss, therefore, does not specifically affect these planetary systems. As a result, in the absence of dynamical interactions these planetary systems are not expected to be affected by either internal planetary dynamics nor by stellar mass loss.

We include, in our simulations, the gravitational interactions between the stars, the interactions inside the planetary systems, and the mass loss due to stellar evolution. In principle, all the three main processes mentioned above are included, although, as mentioned earlier, the effect of stellar evolution is limited by the duration of our simulations. We take all these effects into account as accurately as our computer resources permit, which is particularly important for the long-term dynamical processes among planets orbiting a single star. The simulations are performed using the Astrophysical Multipurpose Software Environment (AMUSE; Portegies Zwart et al., 2009, 2013; Pelupessy et al., 2013). We perform our calculations using a dedicated script, which we call Nemesis, that enables us to integrate the equations of motion of stars with planetary systems and includes the effects of mass loss due to stellar evolution and collisions between stars and planets. Our calculations ignore the primordial gas in the star cluster, but our initial conditions are selected to mimic the initial stellar and planet distribution functions shortly after the primordial gas was expelled and the disks turned into planetary systems. Several example scripts of how AMUSE operates and a more detailed description of the framework is provided in Portegies Zwart & McMillan (2018).

In this work, we focus on the liberation processes and their consequences in a dense star cluster with characteristics comparable to the Orion Trapezium cluster. The majority of the observed field stars and rogue planets may originate from bound clusters, loosely bound associations, and only a minority from isolated stars. Our adopted initial conditions originate from a previous study (Portegies Zwart, 2016) in which the size distribution of circumstellar disks in the Orion Trapezium cluster were reproduced. We considered these conditions suitable for our follow-up study assuming that some of the surviving disks would produce a planetary system. The cluster in the study of Portegies Zwart (2016) was born in virial equilibrium with a fractal density distribution with dimension $F=1.6$. The cluster initially contained 1500 stars with a virial radius of 0.5 pc. At an age of 1 Myr the size distribution of the disks in this cluster is indistinguishable from the observed size distribution of 95 ionized protoplanetary disks larger than 100 au in the Trapezium cluster (Vicente & Alves, 2005).

We adopt the earlier reconstructed initial parameters for the Trapezium cluster and populate the stars that have a surviving disk with a planetary system. The 500 stars with a disk size of at least 10 au at the end of their simulation received either four, five, or six planets with a mean mass of $\sim 0.3\,M_{\rm Jupiter}$. The planets are assumed to have circular orbits in a randomly oriented plane. The correlation between orbital separation and planet mass was selected from the oligarchic growth model for planetary systems by Hansen & Murray (2013) andKokubo & Ida (2002).

After the initialization, we continue the evolution of the star cluster including its planetary systems for 10 Myr to an age of 11 Myr. At that time about half the cluster stars are unbound.

In the following section (§ 2) we describe the setup of our numerical experiment, followed by a description of the initial conditions in § 3. We report on the results in § 4, discuss the results in § 5 and eventually, in § 6, we summarize our findings. In the appendix (§ A) we validate the adopted Nemesis method for integrating planetary systems in stellar clusters.

2 Methods

Integrating planetary systems in star clusters is complicated by the wide range in timescales, ranging from days to millions of years, and the wide range of masses, ranging from Earth-mass up to about 100 ${\rm M}_{\odot}$. The first complication directly indicates that many planetary systems have to be integrated over many orbits, which have to be realized without a secular growth of the error in the energy. The wide range in masses hinders such integrations by introducing round-off and integration errors (Boekholt & Portegies Zwart, 2015). The effect of stellar mass loss complicates the numerical problem. In this section, we describe the methods developed to address these issues.

We use AMUSE for all the calculations presented in this work. This framework is a component library with methods for coupling multi-scale and multi-physics numerical solvers for stellar evolution, gravitational dynamics, hydrodynamics, and radiative transfer. In this paper we incorporate stellar evolution of all the stars in the simulation via the SeBa parametrized stellar evolution code (Portegies Zwart & Verbunt, 1996; Portegies Zwart & Yungelson, 1998; Toonen et al., 2012; Portegies Zwart & Verbunt, 2012). Gravitational interactions between planets are addressed using Huayno, which is a class of a large variety of $N$-body codes based on various kick-drift-kick algorithms via the Hamiltonian splitting strategy of tunable order (Pelupessy et al., 2012). For this work, we adopted the fourth and eighth order shared time-step solvers (Makino & Aarseth, 1992; Nitadori & Makino, 2008). In our case, we adopted the fourth order method for integrating the equations of motion for the stars and the symplectic higher order method for planetary systems.

The computing time for integrating Newtons’ equations of motion of $N$ stars in a cluster scale $\propto N^{2}$. In a relatively small star cluster such as the Trapezium cluster studied in this paper, the integration time step for the top-level parent particles peaks at a fraction of the mean cluster’s crossing timescale, whereas the planetary time step is typically on the order of a few percent of the orbital period around the host star. A multi-time-step approach consequently saves enormously in terms of computer time (see also Aarseth, 1985).

Adding planets to stars increases the number of particles in the system. A more severe performance bottleneck is introduced by the generally tight orbits in which these planets are introduced; i.e., years for planets compared to millions of years for the free-floating stars in the cluster. If all the new objects were introduced in a regular $N$-body code the computation would come to a grinding halt. To prevent this from happening and to reduce the effect of integration errors and round off, we developed the Nemesis package within the AMUSE framework.

The principles that make Nemesis efficient is based on the wide range of scales, which are used as an advantage by separately solving systems that are well separated in terms of temporal or spatial scales. In addition, we introduce the simplification that a planet orbiting one star has a negligible effect on the orbit of a planet around another star in the cluster. This strict separation subsequently allows us to choose different integrators for stars and planetary systems. The latter flexibility allows us to tailor the integration method to the topology of the system. As a consequence, our calculations are naturally parallelized over the many well-separated systems. This results in an enormous acceleration when running on multiple cores because each of the $N$-body integrators can run in parallel for the global intersystem communication timescale. At the same time, energy is conserved per individual system and separately for the global $N$-body system to machine precision. This combination of excellent performance and energy conservation makes Nemesis an ideal tool for integrating planetary systems in star clusters.

2.1 Nemesis module

In Nemesis, planetary systems and stars are integrated together. The underlying assumption is that the entire cluster can be separated into groups. We call these groups “subsystems” or “children” and they can be composed of stars as well as planets that are relatively close together with respect to the size of the cluster. The dynamics in these subsystems is not resolved in the global integrator, which we call the “parent”, but is integrated separately. In many cases, a planetary system is a subsystem, but children may also be composed of several planetary systems that happen to be spatially in close proximity. In this approach, we integrate subsystems separately from the rest of the cluster, but the components of the subsystems and the other cluster objects feel each other’s forces.

2.1.1 Calculating forces

In this section, we explain how the forces in the Nemesis module are calculated. To ease the discussion, we define the term particle. Particles represent the center of masses of a subsystem or of individual objects, such as single stars or free-floating planets. Particles represent the parents in the $N$-body system and are integrated together in one $N$-body code. In practice, the particles are integrated with a fourth- or sixth-order Hermite predictor-corrector method (Makino & Aarseth, 1992; Nitadori & Makino, 2008).

The internal dynamics of each child (the subsystem) is integrated with a separate N-body code. The latter can be a different code, for example, a simple Kepler solver or some high-order symplectic N-body solver. We call this the local subsystem for a particular particle, or the parent’s child. The entire simulation is then composed of as many N-body codes as there are subsystems and one additional code for all the particles that are not part of a subsystem including the center of masses of all the subsystems. The parent system is then composed of subsystems, single stars, or planets.

The gravitational force exerted on each particle is composed of three parts: the forces from all the other objects in the local subsystem, the forces of all the single particles in the global system, and the force of the stars and planets in the other subsystems. In Nemesis we ignore the forces of the individual objects (planets and stars) in the other subsystems. Instead, we take the force from the center of mass of the subsystem into account. As a consequence the stars and planets in a subsystem feel the total force from other subsystems as exerted from the center of mass of that subsystem, but not the individual forces from all the individual components from within that subsystem. Particles in other subsystems, therefore, do not feel the forces of individual planets orbiting a star in the other subsystem. Local particles feel the forces of the other planets and stars in the same system. This procedure, outlined in figure 1, results in a slight error in magnitude and direction of the force on any particle due to the assumption that all objects in another subsystem exert a force from the center-of-mass of that subsystem. As long as a subsystem is composed of a star with some planets, this error remains small, but the error grows when a subsystem is composed of multiple stars. We reduce this error by assuring that subsystems remain small compared to the interparticle distance and that they are not composed of many stars.

2.1.2 Integrating the system

The force calculation in Nemesis is implemented in multiple bridge operations (Fujii et al., 2007; Portegies Zwart & McMillan, 2018). These bridges integrate the equations of motion of the individual components (particles and the subsystems) via a second-order Verlet kick-drift-kick method (see Hut et al., 1995; Jänes et al., 2014).

In the initial kick phase, we accumulate the forces between the single particles and the particles in each of the subsystems. These forces are used to update the velocities of the particles and the objects in each of the subsystems over half a bridge time step, $dt_{\rm bridge}/2$.

In the drift phase, the particles and subsystems are integrated using the forces between the particles in each individual subsystem. Since this is an uncoupled problem, each individual subsystem is integrated in parallel. In this phase, we ignore the forces between the single particles and those that are in subsystems. In the final kick phase, we again calculate the forces between the single particles and the particles in the subsystems based on the new positions after the drift phase, and again update the velocities.

This procedure allows us to integrate particles and subsystems independently. This strict separation of integrating subsystems enables us to adopt a different N-body code for each subsystem, although this is not a requirement. In addition, it makes the concurrent integration of each subsystem possible, which enormously speeds up the procedure for a sufficiently large number of subsystems.

2.1.3 Subsystem dynamics

Subsystems may change their composition at runtime. This can happen when a star or planet is ejected, planets or stars collide, when two or more subsystems merge, or when a single object enters the subsystem. To simplify this process, we recognize two changes to a subsystem:

Merger

Two subsystems are merged to one as soon as their center of masses approaches each other to within the sum of their radii. In this case the radius of each subsystem is the maximum of two radii: it is (1) 5% larger than the distance from the center-of-mass to the outer-most object and (2) the size that corresponds to a likely encounter. The latter is a function of the bridge time step ($t_{\rm nemesis}$), the number of objects in a subsystem, the mass of the subsystem, and a dimensionless factor $\eta$: $t_{\rm enc}=1.0/\eta t_{\rm nemesis}$. We adopt a value of $\eta\simeq 0.2$. Upon the merger of two subsystems, one of the N-body integrators assimilate the other subsystem and the other integrator is terminated. Since both integrators may be different, we assume that the integrator with the largest number of particle survives.

Dissolution

A subsystem can dissolve into individual objects or multiple body parts can split off to form their own separate subsystem. The procedure to decide on the dissolution of a subsystem follows the inverse criteria as for the merger of two or more subsystems. This procedure may lead to the starting of one or more new integrators to take care of the various newly introduced subsystems. Single objects (stars or planets) are incorporated in the global integrator when they escape from a subsystem.

From an astronomical point of view, this procedure looks somewhat arcane, but numerically it has many advantages because it allows us to optimize for efficiency, performance, and accuracy.

2.1.4 Planet and stellar collisions

Apart from the dissolution and merging of dynamical subsystems, we also allow stars and planets to experience physical collisions. Collision can only occur within a subsystem. If two stars in the parent system were to collide, they would first for a separate subsystem within which the collision is handled. Two stars or planets are considered to collide as soon as their mutual distance is smaller than the sum of their radii. A collision always results in a single object, while conserving the mass, volume, and angular momentum in the collision. In principle, it would be relatively easy to perform a hydrodynamics simulation upon each collision, but that is beyond the scope of our current study. A more extensive discussion on such more rewarding events is provided in Portegies Zwart & McMillan (2018).

Isolated stars have a size according to the stellar evolution code, which runs concurrently with the dynamics. The sizes of planets are calculated by assuming a mean planet density of $3$ g/cm3. For improved efficiency, we adopt a special treatment for collisions between planets and the central star of a planetary system. Planets are assumed to collide with their orbiting star as soon as they approach it to within 1 au. This relatively large distance was adopted in order to reduce the computational cost of integrating tight planetary orbits and to minimize the errors associated with their numerical integration. We can easily relax this assumption, but it would result in a considerable increase in computer time.

The new mass of a merged object is the sum of the two individual masses and the new position and velocity are determined by conserving linear momentum and angular momentum. The radius of the collision product of two planets is calculated by conserving the density. A stellar collision acquires its new radius on the stellar evolution track as described in Portegies Zwart & Verbunt (1996).

2.2 Selecting the $N$-body codes in Nemesis

Each subsystem is integrated with a separate $N$-body code. In principle, each of these codes could be different. In practice, however, we use two different techniques to integrate the equations of motion of the stars and planets. The choice of code is based on the requirements for the physics.

For two-body encounters, we adopt a semi-analytic Kepler solver as implemented by Pelupessy & Portegies Zwart (2013). For a typical planetary system in which one particle is much more massive (at least more than 100 times) than the other particles, we use Rebound (Rein & Liu, 2012) with an implementation of a symplectic Wisdom-Holman integrator (WHFAST Rein & Tamayo, 2015). For all other subsystems, we adopt the eighth-order method available in the symplectic integrator Huayno (Pelupessy et al., 2012). The center-of-masses of the subsystems, the single stars, and the free-floating planets are integrated via the Hermite fourth-order predictor-corrector integrator (Makino & Aarseth, 1992; Nitadori & Makino, 2008).

All calculations are executed on a central processing unit (CPU) because the number of particle in each N-body code is relatively small and a graphics processing unit (GPU) would not provide many benefits in terms of speed (Belleman et al., 2008).

2.3 Validation and verification

The performance and accuracy of the Nemesis integrator module is controlled with two parameters: one controls the distance for which individual objects (planets and stars) and subsystems merge or dissolved, and the another controls the time step of the bridge operator. This so-called bridge time step controls the numerical timescale for the interactions between the subsystems and the particles. Both parameters are tuned independently but we choose to express the bridge time step in terms of the encounter distance and the mass of the objects. This adopted scaling leaves only the Nemesis time step, $dt_{\tt Nemesis}$, as a free parameter for integrating the entire $N$-body system. This timescale depends on the topology of the $N$-body system, and we tune its value by performing scaling and validation tests. A detailed analysis of the dependency of the model on the time step in presented in Appendix A. For our choice of initial conditions and integrators we found that an interaction time step of 100 yr gives the most satisfactory results in terms of reproducibility, consistency, energy conservation, and speed.

3 Initial conditions

After developing and validating the numerical framework we can start generating the initial realization for our star cluster with planetary systems. We start the calculations with a cluster of stars, some of which have a planetary system. The initial realization is motivated by Portegies Zwart (2016), who studied the dynamical evolution of the star cluster with 500 to 2500 stars taken from a broken power-law mass function between 0.1 ${\rm M}_{\odot}$  and 100 ${\rm M}_{\odot}$  (Kroupa, 2001). These calculations were performed with a fourth order Hermite $N$-body method including a heuristic description for the size and mass evolution of circumstellar disks. At the start of these calculations each star received a disk with a mass of 1% of the stellar mass and a size of 400 au. During the $N$-body integration the sizes and masses of these disks were affected by close stellar encounters (Jílková et al., 2016). During these simulations the disk size distributions were compared with the protoplanetary disks observed using Hubble Space Telescope WFPC2 of the Trapezium cluster (Vicente & Alves, 2005). In this way Portegies Zwart (2016) was able to constrain the initial cluster parameters. Clusters for which the stars were initially distributed according to a Plummer (1911) distribution did not satisfactorily reproduce the observed disk-size distribution, irrespective of the other parameters, but when the stars were initially distributed according to a fractal with a dimension $F=1.6$ and in virial equilibrium ($Q=0.5$) the simulations satisfactorily reproduced the observed disk size distribution in the Trapezium cluster (KS probability of $\sim 0.8$) in the age range from 0.3 Myr to 1.0 Myr. For our simulations, we adopted the final stellar masses, positions, and velocities for one of these simulations that matched the observed distribution of disk sizes and disk masses best. As a consequence, our initial conditions had already evolved dynamically for 1 Myr before we started our calculation.

In Table 1 we present the initial parameters as adopted by Portegies Zwart (2016) in the left column (indicated by $t=0$ Myr). The third column gives the global cluster parameters at an age of 1 Myr, which are the final conditions for the study performed by Portegies Zwart (2016). We adopted these parameters and in fact, the precise realization of these calculations as initial conditions for our follow-up calculations. The last (fourth) column presents the global cluster parameters at the end of our simulations, at an age of 11 Myr, which is 10 Myr after the introduction of the planetary systems.

During the first 1 Myr of evolution, starting from a fractal spatial distribution (see the leftmost panel in Fig. 2) most of the structure in the initial cluster is lost. The cluster seems to have expanded considerably, as is evidenced by the zoom-out in Fig. 2, but when considering the virial radius has in fact decreased from the initial 0.5 pc to $R_{\rm vir}\simeq 0.36$ pc at an age of 1 Myr. At this moment, we randomly select 500 stars for which the circumstellar disk has survived with a radius of at least 100 au. We subsequently assign a planetary system to 500 of the stars with a surviving disk. The total mass of the planets is identical to the disk mass. The masses and orbital separation of planets are generated using the oligarchic growth model (Kokubo & Ida, 1998) between a distance of 10 au  to 400 au from the host star.

There is no particular reason why we adopted a minimum separation of 10 au, but adopting a smaller minimum separation would have resulted in many more planets with a low mass in very tight orbits. This would have resulted in an enormous increase in computing time. All planets have initially circular orbits with inclination randomly selected from a Gaussian distribution with a dispersion of 1$\degree$  around a plane. This plane is defined as the orbital plane of the planet closest to the star. After the planetary systems are initialized they are rotated to a random isotropic orientation. Each star acquires between 4 and 6 planets with a mass of 0.01 to 130 Jupiter masses (see Fig. 4 and Fig. 6). The total number of planets in the simulation was 2522.

4 Results

When starting the simulation the stars are already 1 Myr old and the stellar density and velocity distribution are the result of the previous calculations reported in Portegies Zwart (2016). We continue to evolve this cluster including its planets for 10 Myr to an age of 11 Myr.

We performed one simulation in which all interactions between stars and planet are taken into account using Nemesis. Snapshots are produced every 1000 yr, but most of the analysis aims at the final snapshot at an age of 11 Myr. A second simulation was performed in which the planetary systems are evolved in isolation without any interactions from other stars. This second run is used for validation purposes only. Even though not explicitly discussed, no free-floating planets were formed in this second run because the initial planetary configurations are intrinsically stable.

4.1 Global evolution of the star cluster

In Fig. 2 we present a projected view of the stars and planets of our simulated cluster at birth (left), at an age of 1 Myr (middle) and at the end of the simulation, at an age of 11 Myr. During the first 1 Myr in which the stars still have circumstellar disks the cluster loses most of its initial fractal structure. During this early phase, the cluster is most dynamically active and the majority of stars experience one or more close encounters with other stars. These encounters cause the truncation of circumstellar disks. By the time we introduce the planetary systems, at an age of 1 Myr, most dynamical interactions have subsided and the cluster has expanded by about an order of magnitude, although the cluster core remains rather compact (see also Table 1). The reduction in density has profound consequences for the survivability of our planetary systems. During the subsequent 10 Myr of evolution the outer parts of the cluster expand by another order of magnitude, but the cluster core remains rather small and bound.

In the overview presented in Table 1 we demonstrate that the cluster hardly loses any mass during its evolution. Mass loss due to stellar winds is rather moderate, reducing the total cluster mass from 618 ${\rm M}_{\odot}$  to 545 ${\rm M}_{\odot}$  in 10 Myr. The majority of this mass loss is caused by the two most massive stars of 73 ${\rm M}_{\odot}$  and 64 ${\rm M}_{\odot}$. These stars experience copious mass loss in the Wolf-Rayet phase followed by a supernova explosion. Such evolution may enrich most of the disk in the cluster by r-processed elements (Portegies Zwart et al., 2018a). The expansion of the cluster by about an order of magnitude and the global mass loss in bound stars cannot be attributed to the stellar mass loss alone. In total, the cluster loses about two-thirds of its stars, one-third in the first Myr, and another third in the following 10 Myr. The structure of the cluster also changes from an initial fractal dimension of $F=1.6$ to $F=1.26$ at 1 Myr and to $F\simeq 0.6$ at the end of the simulation. The eventual cluster, at an age of 11 Myr, can be well described with a Plummer distribution (Plummer, 1911) with a characteristic radius of 0.32 pc. Although, in fig. 2 the cluster appears to expand by two orders of magnitude, the cluster central portion remains rather confined within a parsec.

4.2 Characteristics of the surviving planetary systems

During our calculations, planetary orbits are affected in a number of ways. We start by describing the characteristics of the surviving planetary systems. Later, in § 4.4 and § 4.5 we discuss the planets that are lost due to collisions or ejection from their host star.

In Fig. 3 we present the distribution in eccentricity and semimajor axis of the planets that remain bound up to an age of 11 Myr. About 10% (213 in total) of the planets have experienced considerable orbital variations ($\Delta e>0.1$ or $\Delta a>10\%$) due to a combination of encounters with other stars and internal planetary scattering. We note that in the absence of stellar encounters the planetary systems are not affected by internal scattering. Any changes in the planetary systems in our simulation is therefore the result of interactions with external perturbators (stellar encounters and cluster topology). These interactions put the planets in orbits where internal scattering causes further changes in the orbital parameters.

Some planets acquire eccentricities close to unity, indicating that they may be subject to tidal interactions or even collisions with the host star. Although we ignore tidal effects in our calculations, collisions are taken into account. A total of 75 ($\sim 3.0$ % of the total) planets collided with their parent star and 14 ($\sim 0.6\%$) planets experienced a collision with another planet. We discuss planetary collisions more extensively in § 4.4.

In Fig. 4 we compare the distributions of the number of planets per star in our simulation and compare the distribution with the simulation in which we ignored any stellar encounters. In the latter simulations, the planetary systems are not affected by dynamics and their conditions remain very close to the initial conditions. This indicates that the initial configuration of our planetary systems is stable against internal dynamical evolution.

All stars with planetary systems have either four, five, or six planets initially. In Fig. 4 we subsequently observe that in particular systems with five planets tend to be reduced, whereas only a few stars with four planets or six planets seem to lose any. In addition, by the end of the calculations, the number of systems with three planets seem to be rather small compared to the number of systems with one or two planets. To further quantify the results we also present Table 2, in which we present the number of planets for a star initially (columns) versus the final number of planets (rows).

From Table 2 we see that the systems with three planets by the end of the simulation tend to originate from systems with initially four or five planets. But we find that most systems that initially have five planets reduce directly to one or no planets at all. Curiously enough though, systems that initially have six planets do not lose as many planets, but when they do, they tend to reduce to a single planet, whereas for systems that initially have four planets tend to be rather agnostic about how many planets they lose. Statistically, these changes are significant but much can be attributed to the initial conditions. According to our initial conditions, large disks with a relatively high mass are prone to receiving more planets than small low-mass disks. The large disks tend to be hosted by relatively low-mass stars, and those stars tend to avoid the cluster center, whereas relatively high-mass stars tend to be more abundant in the cluster core. These differences propagate in the distribution of planets and therefore cause an imprint on their future scattering history.

In Fig. 5 we plot the number of planets in a planetary system at the end of our simulations. The majority of stars keep all their planets throughout the calculations, but if a star loses planets, it tends to lose a larger number like three to five rather than just one or two. The lost planets become free-floating or rogue planets, which we discuss in § 4.5.

The redistribution of planets among the stars may also be affected by the masses of the planets. To quantify this we present in Fig. 6 the mean planet-mass as a function of their semimajor axis. The oligarchic-growth model, used to generate the initial planetary systems, leads to more massive planets at larger orbital separation (visible in Fig. 6). To see if there is a mass preference for ejecting planets we also show, in Fig. 6, the final distribution (at an age of 11 Myr). Although the differences between both distributions appear small, the differences at small separation are statistically significant.

To further quantify these findings we present in Fig. 7 the difference in the cumulative distribution of planet mass for the cluster at an age of 1 Myr with respect to 11 Myr. The difference between the two cumulative distributions are small and the fluctuations rather large, but in the final systems, low-mass planets are more abundant than high-mass planets. The turnover occurs near the mean-planet mass in our simulation which is around $1.4$ MJupiter (indicated with the vertical line in Fig. 7). Based on the lack of a correlation between planet mass and orbital separation we argue that the majority of ejections is driven by external perturbations (mostly with other stars) rather than by internal scattering among the planets.

4.3 Migrating and abducted planets

Two rather extreme processes that affect the orbits of planets are their abduction from another star or when a planet is scattered during a close encounter with other planets. In both cases the resulting planet is expected to be parked in a wide orbit with high eccentricity. However, planets that are scattered close to the host star into a parking orbit are expected to have higher eccentricity, on average, than planets abducted from another star (see Jílková et al., 2016).

In our simulations, only a few planets were abducted, and a comparable number of planets were kicked out to the outskirts of their own planetary system by internal scattering. In Fig. 8 we compare the orbital separation and eccentricity of these systems. Although the distributions are rather broad in semimajor axis and in eccentricity, captured planets have on average lower eccentricity and somewhat larger orbital separation compared to ejected planets.

In table 3 we list the migrated planets, and the abducted planets are presented in Table 4. Apart from slight differences in the orbital parameters, the mass of the host star for captured planets tends to be considerably higher than for the migrated planets. This trend is not unexpected because of the stronger gravitational influence of more massive stars whereas low-mass stars are more prone to lose planets.

The abducted planets in Table 4 appear to have large semimajor axes and a broad range in eccentricities. Such abduction explains the observed orbital parameters of the dwarf planet Sedna in the solar system (see Jílková et al., 2015). As an alternative to abduction, a free-floating planet could in principle be captured by a star or planetary system. Capturing free-floating planets was also studied in Goulinski & Ribak (2018), who argued that these systems may not be uncommon, but that they would have a wide range in eccentricities and typically large semimajor axes (Perets & Kouwenhoven, 2012). In our simulations no free-floating planets were captured, and we do not expect this to be a common process because 80% of the ejected planets escape promptly from the cluster (see § 4.5).

4.4 Characteristics of colliding planets

One important aspect of planets is their finite size, which makes them prone to collisions. A total of 75 planets in our simulations collide with another planet or with the parent star. In our simulations, collision with the parent star is not treated realistically in the sense that we ignored tidal effects. We compensate for this by adopting a size of 1 au for planetary-hosting stars. As a result, we overestimate the number of collisions with the host star and we do not acquire hot Jupiter planets. We, therefore, focus on the collisions that occur between planets.

In Table 5 and Table 6 we list the mergers that occurred in our simulations sorted in the moment of the collision. In Table 5 we show the pre-collision parameters of the two planets that participate in the collision, whereas in Table 6 we list the orbital parameters of the merger product.

The orbital parameters for the pre-merger planets are derived from the last snapshot before the merger occurred, which can be up to 1000 years before the actual event. The mean mass of the primary in a colliding planet pair is 1.14 ${\rm M}_{\rm Jup}$  and a secondary of 0.36 ${\rm M}_{\rm Jup}$. The resulting merger product is 1.5 ${\rm M}_{\rm Jup}$. During the calculation, 34 planets collided in a total of 19 events. Several planets experienced multiple collisions, causing the planet mass to increase very effectively and causing the planet to migrate closer toward the host star. These multiple mergers all tend to occur in relatively short succession.

Most mergers tend to occur between neighboring planets, but there are seven occasions where one or more intermittent planets are skipped. In particular the event at $t=5.91$ Myr is interesting because in this case, the outermost planet collides with the innermost planet. Although not taken into account in our calculations, such close encounters among planets could lead to the capture of one planet by the other, giving rise to a binary planet as was observed in Kepler 1625 (Hamers & Portegies Zwart, 2018).

4.5 Production of free-floating planets

By the time the cluster has reached an age of 11 Myr the total mass in bound planets was reduced from $\sim 3527$ $M_{\rm jup}$ to $\sim 2915$ $M_{\rm jup}$ (see Table 1). Planets have been lost by their parent star via encounters with other stars (see § 4.5) , internal planet-scattering ($\sim 60$), by the mass loss of their host stars, and through collisions with the star (75; see § 4.4) of collided with another planet (14). Once liberated, free floaters may remain bound to the cluster (75 planets) or escape its gravitational potential (282, see Table 1). In total 357 planets (out of 2522) were liberated from the gravitational pull of their parent star. In § 4.3 we discussed the possibility of captured planets, but this was not the fate of any of the free-floating planets, because all captured planets were exchanged during a close encounter and always bound to at least one star.

In Fig. 9 we present the number of free-floating planets as a function of time. The majority of free floaters (67%) leave the cluster within a crossing time ($\sim 1$ Myr) after being liberated from their host star. The other $\sim 33$% remain bound to the cluster for an extended period of time and leave the cluster on a much longer timescale, at a typical escape rate of $\sim 8$ planets per Myr.

In Fig. 10 we present the cumulative distributions of the velocity of bound and unbound stars and planets; for the planets we make a distinction between free-floating planets that remain bound to the cluster and those that escape. The distributions for the stars and planets at an age of 11 Myr that are still bound to the cluster show only slight differences (thin lines). Both velocity distributions are statistically indistinguishable (with a KS-statistic of $0.23$). The population of unbound planets, however, tend to have much higher velocities (of $\sim 3^{+5.6}_{-1.2}$ km/s) than the stars ($\sim 1^{+1.3}_{-0.6}$ km/s). This is not unexpected because planets tend to be launched from the stars with their orbital speed, which gives rise to higher mean escape velocity, whereas most stars escape by dynamical evaporation (Fukushige & Heggie, 1995; Portegies Zwart & Takahashi, 1999). This relatively high space motion of the rogue planets is also reflected in the large percentage of liberated planets that escape the cluster.

In Fig. 11 we present the mass distribution of free-floating planets. Those that remain bound to the cluster have statistically the same mass function as those that escape (KS-statistics of $0.14$) and as the global initial planet mass function (KS $=0.11$; see also Fig. 6). Signifying what we already discussed in relation to Fig. 6 and Fig. 7: the ejection of planets is independent of their mass (see also Malmberg et al., 2011; Veras & Moeckel, 2012).

The mass distributions of free-floating planets in the simulation differ considerably from the observed mass distribution. Observational selection effects probably play an important role here because low-mass free-floating planets tend to be very hard to discover. We, therefore, introduce a lower limit of $2.5\,M_{\rm Jup}$ to the mass distribution the simulated distribution of free floaters becomes statistically indistinguishable from the observed sample (KS-statistic is 0.06).

In relation to Fig. 7, we argued that the lack of a mass-dependency of the production of free-floating planets is mainly caused by the importance of strong encounters with other stars rather than internal scattering among planets. To quantify this hypothesis we present in Fig. 12 the cumulative distributions of the number of strong and weak encounters. In this figure strong indicates an encounter within 1500 au. In this analysis, a planet that becomes free floating within 0.5 Myr of such a strong encounter is considered to be liberated as a result of this, otherwise, we consider the planet to be lost as a result of a weak encounter or the internal reorganization of the planetary system.

To further understand the importance of strong encounters we present in Fig. 13 the delay time distribution of liberated planets. The majority of those escape promptly upon a strong encounter with another planetary system or a single star. A considerable number ($\sim 24$ %) require more time (up to about a million years) before they escape from their host star. In this latter population, planetary escape is initiated by the close encounter, but it requires the planetary system to become dynamically unstable before the planet is actually ejected. The timescale for these planetary systems to become unstable appears to be on the order of a million years.

The number of Jupiter-mass free-floating planets have been estimated to about 0.25 of the number of main-sequence stars (Cassan et al., 2012; Mróz et al., 2017). This number is consistent with our findings, even though we adopted that only about one-third of the stars had planets initially. If each star would have a planetary system, our estimates would rise to about $\sim 0.72$  free-floating planets per main-sequence star, which would be on the high side but not inconsistent with the observed estimate of $1.8^{+1.7}_{-0.8}$ (Sumi et al., 2011). Although not taken into account in this work, the number of free-floating planets produced per star depends on the moment circumstellar disks start forming planetary systems, their distribution in mass and orbital parameters, and on the density and velocity distribution of the young cluster.

5 Discussion

We simulated the evolution of a cluster of 1500 stars of which 500 are orbited by a total of 2522 planets (4, 5, or 6 planets of 0.008 ${\rm M}_{\rm Jup}$  to 130 ${\rm M}_{\rm Jup}$  per star in circular planar orbits between 10 au and 400 au). The calculations were performed via the Nemesis script in the Astrophysical Multipurpose Software Environment (Portegies Zwart, 2011; Portegies Zwart et al., 2018; Portegies Zwart & McMillan, 2018) and include the effects of stellar mass loss and the interactions between all objects. We took the initial conditions from earlier calculations that mimic the mass and size distributions of the Orion trapezium star cluster (Portegies Zwart, 2016). We stopped the calculations at an age of 11 Myr, after which we analyzed the population of planets.

In our calculations, we ignored the effect of tidal energy dissipation between stars and planets. When we started this study we argued that this effect had minor consequences, but it turned out that 75 of the planets (3.0 %) have a strong interaction with their host star and 34 planets collide with other planets. Tidal interactions are clearly important and we will improve this in a future version of Nemesis. Considering these systems as resulting either in a collision with the parent star or the formation of a hot Jupiter, we derive a hot-Jupiter formation efficiency of 75 per 500 planetary systems per 10 Myr, or 15% of the planetary systems produce a hot Jupiter, which is not inconsistent with the rate derived by Heller (2018).

Our study mainly focuses on the production of free-floating planets. The planet-ejection probabilities in our simulations are independent of the mass of the planet, which contradicts earlier results of Malmberg et al. (2011); Davies et al. (2014). Part of this result probably depends sensitively on our initial distributions of planet mass and orbital topology. The choice of oligarchic growth causes the more massive planets to be further away from the host star, where they are more vulnerable to perturbations by passing stars. This makes the inner planets more prone to being ejected in the subsequent unstable planetary system that results from an external perturbation.

Our finding that the probability of escaping the parent star is independent of planet mass and the birth distance from the star is a direct consequence of the way in which planets are freed, i.e., in most cases this is the result of a strong encounter between the planetary system and another star or planetary system. In our simulations, interactions between planets and stars lead to a total of 357 free-floating planets from an initial population of 2522 bound planets. This results in 0.24 to 0.70 free-floating planets per main-sequence star, which is consistent with estimates of the number of free-floating planets in the Galaxy by Cassan et al. (2012) and Mróz et al. (2017).

An important reason for the relatively small number of free floaters is their relatively late formation. Most interactions occur in the first 1 Myr of the evolution of the cluster, and strong dynamical encounters drive the size evolution of the circumstellar disks in this phase. By the time we introduced the planets the stellar density had already dropped considerably and the number of strong interactions had subsided. The absence of planets in the first million years enables them to survive to a later epoch. If these disks were already rich in debris or planets they would have been much more vulnerable to external perturbations. The mutual interactions between stars in the earliest cluster evolution \ {\raise-2.15277pt\hbox{\buildrel<\over{\sim}}}\ 1 Myr would have been sufficient for ionizing most planetary systems, leading to a larger population of free-floating objects. Such a sola lapis has recently been found traversing the solar system (Portegies Zwart et al., 2018b).

The distribution of the masses of free-floating planets in our simulation is indistinguishable from the mass distribution of planets bound to their host star. This may have interesting consequences for observations. This comparison may also be made for observed planets. Our cluster is not old enough to produce free-floating planets by the copious stellar mass loss in the post-asymptotic giant branch phase, and it is not a priori clear what effect this would have on the distribution and ejection of multi-planet systems. But to first order we argue that the distribution of free-floating planets is the same as that of bound planets.

6 Conclusions

We simulated the evolution of the Orion Trapezium star cluster including planets. The calculations start with initial conditions taken from earlier calculations at an age of 1 Myr from Portegies Zwart (2016) by converting circumstellar disks into planetary systems and were continued to an age of 11 Myr. Our calculations, performed with AMUSE, include the effects of stellar mass loss, collisions, and the dynamics of the stars and planets. The orbits of the planets are integrated using a symplectic direct $N$-body code whereas the stellar dynamics is resolved using a direct Hermite $N$-body code.

Realizing that we study a chaotic system based on the result of only two simulations, one without stellar interactions and one that included interactions between the planets and the stars, we nevertheless feel sufficiently bolstered by our results to report a number of conclusions. Each of these conclusions is based on the results obtained from the simulation in which all interactions between stars and planets are taken into account. The results enumerated below are therefore rather empirical, although, as argued in the main text, some of these conclusions may be fundamental. All conclusions, however, are a result of the complicated interplay between initial conditions and simulations, and it is sometimes hard to disentangle the two.

• Conclusions regarding planet stability

• $\bullet$

The majority of planets ($\sim 70$%) experience a change in their orbits (in eccentricity or semimajor axis) of less than 5%.

• $\bullet$

A small number of $\sim 10$ % planets acquire a high (\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 0.8) eccentricity. This is not necessarily caused by stars passing closely, but in the majority of cases repeated small perturbations within the cluster and subsequent secular evolution within the planetary system drives these high eccentricities.

• $\bullet$

High eccentricities are also induced by collisions between planets and in the orbits of captured planets.

• $\bullet$

The innermost planets (at 10 au) experience a comparable relative variation in their final orbital parameters (in particular the eccentricity and inclination) due to encounters, perturbations, and internal secular evolution as wider systems.

• $\bullet$

The probability for a planet to escape is independent of its mass or semimajor axis. Low-mass planets that are born relatively close to the parent star are only marginally more prone to ejection than more massive planets born further out (see also, Malmberg et al., 2011; Veras & Moeckel, 2012). This result, however, probably depends sensitively on the initial orbital distributions and masses of the planets. Comparing observed planet-mass distributions and those that survived in a planetary systems may then provide interesting constraints on the initial planet mass function.

• $\bullet$

Seventy-five planets (3.0%) collide with their host star. This number, however, strongly depends on our adopted stellar collision radius and will change when tidal evolution is properly taken into account, but we still expect that collisions between a planet and its host star are rather frequent. Although our collisions are not taken into account realistically because of the large stellar size we adopted, these systems would be eligible to the formation of hot Jupiter planets at a rate of $\sim 0.015$ per star per Myr.

• $\bullet$

The widest planetary systems in our simulations tend to be formed either by ejecting planets on very wide and highly eccentric orbits or by capturing a planet from another star. Both methods seem to be equally important, but the captured planets tend to have somewhat lower eccentricity.

• Conclusions regarding planetary escapers

• $\bullet$

A total of 357 planets (out of 2522 or $\sim 16.5$ %) become unbound from their parent star.

• $\bullet$

Out of 357, 282 ($\sim 80$%) of the free floating planets promptly escape the cluster upon being unbound from their parent stars.

• $\bullet$

The probability for a planet to escape is independent of its mass. As a consequence, the mass function of free-floating planets and the mass function of bound planets are indistinguishable from the initial distribution of planet masses.

• $\bullet$

At the end of our simulations systems with 3 planets were rare compared to systems with 1 or 2 planets, or systems with 4 or more planets. Once a star loses planets, it tends to lose 3 or more (consistent with Table 9 of Cai et al., 2017).

• Conclusions regarding planet collisions

• $\bullet$

Thirty-four planets (1.3%) experienced a collision with another planet.

• $\bullet$

The collision probability between two planets is independent of planet mass.

• $\bullet$

The orbits of planet-planet collisions have a mean eccentricity of $0.33\pm 0.19$ and a relative inclination of $20^{\circ}\pm 35^{\circ}$.

• $\bullet$

Instead of colliding, some of those events may lead to the tidal capture of one planet by another. This would lead to the formation of a binary planet, or moon, as was observed in Kepler 1625B Teachey et al. (2018).

• $\bullet$

It is generally the outermost planet that collides with a planet closer to the parent star. This inner planet is not necessarily the next nearest planet.

• $\bullet$

Planets regularly engage in a cascade of collisions. These chain-collisions are initiated by a dynamical encounter with another star.

• Conclusions regarding the host star clusters

• $\bullet$

The host star ejects 240 (67% of the ejected planets, 10% of all planets) planets with a delay of $0.1$–$0.5$ Myr after the last strong encounter with another cluster member.

• $\bullet$

Young $\sim 10$ Myr old star clusters harbor a rich population of free-floating planets. About one-third of the free-floating planets remain in the cluster for more than a dynamical timescale, up to the end of the simulation. The number of free floaters in these clusters can be as high as $40$ planets for stars between 0.9 ${\rm M}_{\odot}$  and 1.1 ${\rm M}_{\odot}$, or $25$ % of the main-sequence stars (consistent with estimates by Cassan et al., 2012).

A large number (30%) of planetary systems are affected by the presence of the other stars in the cluster, but only $\sim 10$% of those will leave a recognizable trace that allows us to reconstruct the dynamical history based on the topology of the inner planets. For the majority of planetary systems observed today, current instruments are unable to discern the dynamical history because we only observe the inner most planets, rather than the outer parts where dynamical effects are most pronounced. It would require observation of a exo-Kuiper belt to be able to establish the past dynamical history of the planetary system. Possibly the easiest way to perceive the dynamical history of a planetary system is preserved in collision products between planets. We argue that in more than 3% of the planetary systems collisions between planets are initiated by external dynamical perturbations. From the $\sim 4000$ planetary systems known today we then expect more than one hundred to host a collision product.

About 16% of planets eventually become dissociated from their parent star due to interactions with other cluster members or internal reorganization of the planetary system. These ejected planets become free floaters. The majority of those (\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 80%) leave the cluster within a crossing time scale, the rest lingers around the cluster potential and are subject to a slower evaporation process driven by mass segregation. We therefor expect star clusters to be relatively poor in free floating planets. The Galactic field, on the other hand, is contains about 1/4-th of the number of free floating planets as there are main-sequence stars. The Galaxy is then composed of some $5\times 10^{10}$ free floating planets, of which only a dozen are observed.

Acknowledgment

We are grateful to the anonymous referee for many interesting comments that helped us improve the manuscript. In this work we used the following packages: AMUSE (Portegies Zwart, 2011; Portegies Zwart et al., 2018) (which is also available at http://amusecode.org), Hermite0 (McMillan, 2014), Hop (Eisenstein & Hut, 2011), Huayno (Pelupessy et al., 2012), matplotlib (Hunter, 2007), numpy (Oliphant, 2006), Python (van Rossum, 1995), and SeBa (Portegies Zwart & Verbunt, 2012). The calculations ware performed using the LGM-II (NWO grant # 621.016.701) and the Dutch National Supercomputer at SURFSara (grant

15520).

Appendix A Validation

We analyze the accuracy of the hybrid Nemesis strategy as a function of the interaction time step $dt_{\tt Nemesis}$.

A.0.1 Determining the optimal Nemesis time step

This Nemesis time step ($dt_{\tt Nemesis}$) numerically associates two important factors: how often forces between subsystems are calculated and a measure for the interaction distance between individual particles ($d_{\tt Nemesis}$). If two particles are separated by less than this interaction distance ($d_{\tt Nemesis}$) a new subsystem is created within which the interaction between particles is resolved with a separate $N$-body integrator. In principle we create a new subsystem with its own individual $N$-body solver. In practice, however, many of these individual subsystem $N$-body solvers are the same code.

If one particle is spatially separated from several other particles in a subsystem by a distance larger than the interaction distance, $d_{\tt Nemesis}$, this particle is removed from the subsystem and incorporated in the global cluster integration code. For the physics it makes no difference if a particle is part of the global system or of a subsystem. However, the integrator used for any of the subsystems is symplectic and generally more accurate by adopting higher order and a smaller time step, whereas the global $N$-body code adopts larger time steps and is not symplectic.

There is no specific requirement for any particle to be integrated either by the integrator of a subsystem or by the global integrator. The choice of the domain to which the particle belongs is purely based on geometry and the adopted demands for accuracy and precision. In practice, the entire cluster including all the planets could either be integrated by the single global fourth order Hermite code or by one of the symplectic $N$-body codes of the subsystem. The choice of which particle is integrated by what integrator is then only decided on terms of accuracy, precision, and performance.

As a general note, however, the global $N$-body code tends to be less accurate because of larger time stepping and non-symplectic, whereas the subsystem codes adopt rather small shared time steps with a symplectic integrator. As a consequence, we prefer to keep particles that belong to a single planetary system in the same integrator.

The number of stars and planets that are embedded within a single subsystem depends on $dt_{\tt Nemesis}$ (and therefore on $d_{\tt Nemesis}$). In Fig. 14 we show how this number varies as a function of $dt_{\tt Nemesis}$. For very small values of $dt_{\tt Nemesis}$, all the stars and planets are integrated by the global $N$-body integrator, and the number of subsystems $n$ drops to 1, in the extreme. On the other hand, if dt_{\tt Nemesis}\ {\raise-2.15277pt\hbox{\buildrel>\over{\sim}}}\ 200 yr all initial planetary systems are recognized as individual subsystems and assigned their own integrator. In that case, the number of subsystems grows to the actual number of planetary systems we initialized plus one for the global $N$-body system, and $n$ approaches to a value of $501$. We draw a vertical line at $dt_{\tt Nemesis}=100$ yr, which corresponds to our adopted Nemesis time step. For this value, a total of about 400 $N$-body integrators are being initialized and run concurrently.

A.0.2 Subsystem size criterion in Nemesis

The analysis performed in the previous section is calculated on a static initial realization without evolving the cluster dynamically. In § A.0.1 we demonstrated that at a larger time step individual planetary systems are consistently captured in a subsystem. A larger time step is also preferred because this requires fewer interaction steps between the subsystems and the other particles. The evolution of the cluster, however, is dynamic and as a consequence, the value of $dt_{\tt Nemesis}$ should be dynamics to warrant the accuracy and efficiency of the Nemesis method. We tested this hypothesis by integrating the cluster for $0.1$ Myr with various values of $dt_{\tt Nemesis}$. After this time we measured the radius of the largest resolved subsystem. These largest resolved subsystems tend to slow down the integration because they are likely to be composed of a larger number of particles (stars and planets). Such large subsystems may cause the entire calculation to wait for the integration of the large subsystem. the calculation becomes progressively slower when more particles are incorporated in the subsystem. Eventually, this may continue until all the particles are embedded in a single subsystem, which is beyond the purpose of the Nemesis module.

In Fig. 15 we present the measured size of subsystems as a function of $dt_{\tt Nemesis}$. The optimum is reached for a $dt_{\tt Nemesis}\simeq 100$ yr, which results in a maximum radius for subsystems of $\sim 1738$ au. The choice of a time of $dt_{\tt Nemesis}\simeq 100$ yr results in the most efficient calculation of the entire stellar system while at the same time it results in the lowest energy error. With this time step our calculations conserve energy better than one part in $10^{4}$ per planetary system per million years, which is sufficient to preserve the phase space characteristics of $N$-body systems for the 10 Myr over which we performed the simulation (Portegies Zwart & Boekholt, 2014).

The two criteria, i.e., (1) keep each initial planetary system in a single subsystem and (2) prevent subsystems from boundless growth, suggest opposing optimal values for the Nemesis time step $dt_{\tt Nemesis}$. Both criteria appear to match for $dt_{\tt Nemesis}\simeq 100$ yr, which is the value we adopt for all further calculations.

A.0.3 Validation of Nemesis on individual planetary systems

Apart from tuning the performance and accuracy of the compound Nemesis integrator, we also validated this code in a more practical application. For this we opted to study the evolution of a system of five planets that is orbited by another second star of $1$M*⊙* with a semimajor axis of $1500$ au, an eccentricity of $0.5,$ and an inclination of $90\degree$. The planetary system is generated using the oligarchic growth model for a 1 ${\rm M}_{\odot}$  star with a 400 au disk of 0.1 ${\rm M}_{\odot}$. The simulations were performed via two distinct methods: (1) using Nemesis and (2) integrating all objects in a single $N$-body code. The Nemesis method requires two codes: one for the planetary system and one for the center of mass of the planetary system and the orbiting secondary star.

For both integrators, we selected the eighth order symplectic integrator in Huayno. The two codes communicate using a Nemesis time step of $dt_{\tt Nemesis}=100$ yr. For comparison, we also integrated these planetary systems with the same integrator, but all the objects stars and planets are in the same computational domain. In Fig. 16 we present the eccentricities of the planets as a function of the semimajor axis at an age of 0.5 Myr.

Based on the integration of these planetary systems and the earlier tests regarding the migration of planets across integrators, we decided that a Nemesis time step of $dt_{\tt Nemesis}=100$ yr gives satisfactory results in terms of accuracy, precision, and performance.

A.0.4 Energy errors in the composite model

To determine the reliability of the Nemesis for planetary system evolution, we also investigated the evolution of the energy error. We performed this test for the same model as in the previous section, using an isolated planetary system composed of five planets and one perturbing star in a wide orbit. We simulated this system using our method and a fourth order Hermite integrator using a time step of $dt_{\tt Nemesis}=100$ yr. The resulting evolution of the energy error is presented in In Fig. 17.

In Fig. 17 we show the results of the two calculations, one with a fourth order Hermite integrator (green), which is not symplectic. The other calculation (red curve) is performed using Nemesis in which we combine an eighth-order symplectic integrator for the planetary system with the fourth-order Hermite integrator for the binary system. The energy error in the Hermite (green curve) grows monotonically, which is the typical response for a non-symplectic integrator, such as the adopted Hermite scheme. The evolution of the energy error in the hybrid integrator does not grow on a secular timescale. The evolution of the energy error is rather erratic with sharp peaks to low values as well as high values but stays stable overall. The secular growth of Nemesis is much smaller than the single Hermite integrator. This is mainly caused by the fact that the largest energy errors are generated while integrating the planetary system, which, in the Hermite integration (green curve) drives the energy error. An additional advantage is that the calculation with the hybrid Nemesis method took about ten minutes on a workstation, whereas the Hermite scheme (green curve) took 18 hours.

Based on the results presented in Fig. 17, we conclude that in our method the energy error does not grow with time, but remains constant for the duration of the calculation. The Hermite part of the integration does show a monotonic increase of the energy error, but this error remains below the mean error produced in the subsystem code, which is symplectic. The overall energy error, therefore, appears well behaved, but eventually, in the long run, the non-symplectic part of the energy error may start to dominate.