Free Access
Issue
A&A
Volume 624, April 2019
Article Number A120
Number of page(s) 16
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201834641
Published online 24 April 2019

© ESO 2019

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 by catch 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; Mróz et al. 2019). 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ēs1 expected from the star formation processes (Portegies Zwart et al. 2018a).

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 offree-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 orderto 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 ~0.3 MJupiter. 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) and Kokubo & 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 Sect. 2 we describe the setup of our numerical experiment, followed by a description of the initial conditions in Sect. 3. We report on the results in Sect. 4, discuss the results in Sect. 5 and eventually, in Sect. 6, we summarize our findings. In 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 M. 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, 2012; Portegies Zwart & Yungelson 1998; Toonen et al. 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 ∝ N2. 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.

thumbnail Fig. 1

Diagram of Nemesis method. (A) A particle is an individual object or a subsystem consisting of multiple individual objects. In this study, the individual objects are either stars or planets. (B) The gravitational force on a particle is the sum of the force from other particles [1] and the forces from the individual objects in the subsystems.(C) The gravitational force on individual objects is the sum of the force from the particles [1] and the forces from the other individual objects in the containing subsystem [2], but not from those individual objects in other subsystems [3]. The forces from A.1 and B.2 are each in a self-contained system and can be calculated in an N-body code; the forces in A.2 and B.2 are connected to the self-contained systems and are evolved with a leapfrog algorithm.

Open with DEXTER

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 Fig. 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 forcesare used to update the velocities of the particles and the objects in each of the subsystems over half a bridge time step, dtbridge∕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 (tnemesis), the number of objects in a subsystem, the mass of the subsystem, and a dimensionless factor η: tenc = 1.0∕ηtnemesis. We adopt a value of η ≃ 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 assoon 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).

Isolatedstars 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 cm−3. 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 adoptedin 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, d tNemesis, 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 M and 100 M (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 constrainthe 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 ~ 0.8) in the age range from 0.3 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 Rvir ≃ 0.36 pc at an age of1 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 to 400 au from the host star.

There isno 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° 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 Figs. 4 and 6). The total number of planets in the simulation was 2522.

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.

thumbnail Fig. 2

Projected view of the simulated star cluster at t = 0 Myr (initial conditions adopted by Portegies Zwart 2016; left panel), at t = 1 Myr (middle panel and the adopted initial conditions), and at t = 11 Myr (right panel, our final conditions). Stars are red bullets, single free floating planets black triangles.

Open with DEXTER

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 M to 545 M in 10 Myr. The majority of this mass loss is caused by the two most massive stars of 73 M and 64 M. 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. 2018b). 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 ≃ 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 Sects. 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 (Δe > 0.1 or Δ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 (~ 3.0% of the total) planets collided with their parent star and 14 (~0.6%) planets experienced a collision with another planet. We discuss planetary collisions more extensively in Sect. 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 toreduce 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 Sect. 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 theinitial 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.

thumbnail Fig. 3

Eccentricity as a function of the semimajor axis the planets that survive up to an age of t = 11 Myr.

Open with DEXTER
thumbnail Fig. 4

Histogram for the number of systems with a certain number of planets. The dotted curve gives the initial distribution with either 4, 5, or 6 planets per star. The final conditions for the simulation without stellar dynamics are identical to this initial distribution. The distribution of the simulation in which we included the stellar encounters at an age of 11 Myr is presented as the solid curve with slanting lines.

Open with DEXTER
Table 2

Comparison of the distribution of planets at the beginning and at the end of the simulation.

thumbnail Fig. 5

Histogram of the number of systems with a certain number of lost planets. Of the original 500 planetary systems the majority (386) do not lose any planets. Only 25 systems lose 1 or 2 planets and 102 systems lose 3 or more planets.

Open with DEXTER
thumbnail Fig. 6

Mean planet mass as a function of semi major axis (in a moving bin of 50 planets). The initial (at 1 Myr, in black) and the final (at 11 Myr in red) mean mass only differ slightly. The mean mass, 1.4 MJupiter, is depicted with a green horizontal line.

Open with DEXTER

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 Sect. 4.5).

thumbnail Fig. 7

Relative difference between the cumulative distributes of the masses of the bound planets initially and at 11 Myr. Positive values indicate an overabundance at the end of the simulation. Initially, the mean planet mass is 1.398 ± 1.05 MJup at 11 Myr the mean mass is only fractionally different at 1.404 ± 4.191 MJup; the latter value is indicated by the vertical green line.

Open with DEXTER
thumbnail Fig. 8

Eccentricity as a function of the semimajor axis for captured planets (black diamonds) and migrated planets with semimajor axis larger than 800 au (red dots) at an age of t = 11 Myr. The mean and standard deviation for both sets are also plotted. The mean orbital elements for the captured planets is ac = 1539 ± 824 au and ec = 0.6 ± 0.2, and ac = 1141 ± 258 au and ec = 0.8 ± 0.2 for the migrated planets.

Open with DEXTER
Table 3

Parameters for planetary systems in which one planet was ejected to a larger distance (> 800 au) from its host star.

Table 4

Listing of systems that formed by the abduction of a planet from another star.

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 Tables 5 and 6 we list the mergers that occurred in our simulations sorted in the momentof 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 yr before the actual event. The mean mass of the primary in a colliding planet pair is 1.14 MJup and a secondary of 0.36 MJup. The resulting merger product is 1.5 MJup. 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 timethe cluster has reached an age of 11 Myr the total mass in bound planets was reduced from ~ 3527 Mjup to ~ 2915 Mjup (see Table 1). Planets have been lost by their parent star via encounters with other stars (see Sect. 4.5), internal planet-scattering (~60), by the mass loss of their host stars, and through collisions with the star (75; see Sect. 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 Sect. 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 (~ 1 Myr) after being liberated from their host star. The other ~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 ~ 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 ) than the stars (). 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 Figs. 6 and 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 MJup 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. Themajority of those escape promptly upon a strong encounter with another planetary system or a single star. A considerable number (~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 ~ 0.72 free-floating planets per main-sequence star, which would be on the high side but not inconsistent with the observed estimate of (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 youngs cluster.

Table 5

Orbital elements of the merging planets.

Table 6

Orbital elements of the planets resulting from a merger.

thumbnail Fig. 9

Number of free-floating planets (Nfp) as a function of time. The solid curve (black) indicates all free planets; the red curve indicates the subset of free floaters that also escape the cluster.

Open with DEXTER
thumbnail Fig. 10

Cumulative distribution (normalized) of the velocity of planets (black curves) and stars (red curves) at an age of 11 Myr. Planets and stars bound to the cluster are plotted with a thin line; the thick curves indicate the unbound objects.

Open with DEXTER
thumbnail Fig. 11

Cumulative distribution of the masses of all planets (red), the free-floating planets that are bound to the cluster (green), and those unbound from the cluster (blue). These three curves are statistically indistinguishable. The dotted curve indicates the mass distribution of 16 observed potential free-floating planets from Luhman et al. (2005); Marsh et al. (2010); Zapatero Osorio et al. (2000); Delorme et al. (2012); Liu et al. (2013); Gagné et al. (2014a,b,c, 2015); Schneider et al. (2014, 2016); Luhman (2014); Liu et al. (2016); Kellogg et al. (2016). For a different comparison, we introduce a lower mass cutoff to the initial sample of planets of 2 MJup and compare this with the observed sample (black).

Open with DEXTER

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 MJup–130 MJup per star in circular planar orbits between 10 and 400 au). The calculations were performed via the Nemesis script in the Astrophysical Multipurpose Software Environment (Portegies Zwart 2011; Portegies Zwart et al. 2018c; 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 independentof 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–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 ≲ 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. 2018a).

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.

thumbnail Fig. 12

Number of planets that became unbound from their host star as a function of time. The number of planets that escaped their host within 0.5 Myr following a strong encounter (within 1500 au, red curve) is about twice as large as the planets thatescape without evidence of having experienced a strong encounter (black curve). The dotted black curves indicate thedependency on the timescale within which a strong encounter is supposed to lead to ejected planets; the lower curves indicate the cumulative distribution for planets that are liberated within 1 Myr of a close encounter, whereas the upper curve is for 0.2525 Myr).

Open with DEXTER
thumbnail Fig. 13

Number of planets that escape from their host star as a function of the time between a close encounter (within 1500 au) and the moment of escape. The majority of the planets escape promptly upon an encounter, but a considerable number require more time, up to about a million years.

Open with DEXTER

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 andwere 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

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

  • A small number of ~10% planets acquire a high (≳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.

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

  • 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.

  • 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.

  • 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 ~0.015 per star per Myr.

  • 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

  • A total of 357 planets (out of 2522 or ~16.5%) become unbound from their parent star.

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

  • 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.

  • 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

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

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

  • The orbits of planet-planet collisions have a mean eccentricity of 0.33 ± 0.19 and a relative inclination of 20° ± 35°.

  • 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).

  • 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.

  • 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

  • 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 thelast strong encounter with another cluster member.

  • Young ~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 M and 1.1 M, 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 ~ 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 ~ 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 (≳ 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 × 1010 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. 2018c; 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 dtNemesis.

A.1 Determining the optimal Nemesis time step

This Nemesis time step (dtNemesis) numerically associates two important factors: how often forces between subsystems are calculated and a measure for the interaction distance between individual particles (dNemesis). If two particles are separated by less than this interaction distance (dNemesis) 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, dNemesis, 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 dtNemesis (and therefore on dNemesis). In Fig. A.1 we show how this number varies as a function of dtNemesis. For very small values of dtNemesis, 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 dtNemesis ≳ 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 dtNemesis = 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.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 Sect. A.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 dtNemesis should be dynamicsto 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 dtNemesis. 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.

thumbnail Fig. A.1

Number of the initial intact planetary systems as a function of the Nemesis time step. The chosen time step of 100 yr is shown as a green vertical line. The time step is not optimal for this criterion, but was chosen as it gives better accuracy and a higher computational speed.

Open with DEXTER

In Fig. A.2 we present the measured size of subsystems as a function of dtNemesis. The optimum is reached for a dtNemesis ≃ 100 yr, which results in a maximum radius for subsystems of ~1738 au. The choice of a time of dtNemesis ≃ 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 104 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 dtNemesis. Both criteria appear to match for dtNemesis ≃ 100 yr, which is the value we adopt for all further calculations.

A.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°. The planetary system is generated using the oligarchic growth model for a 1 M star with a 400 au disk of 0.1 M. 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.

thumbnail Fig. A.2

Sizeof the largest subsystem as a function of dtNemesis. after 0.1 Myr of evolution. The vertical green line indicates the adopted value of 100 yr.

Open with DEXTER
thumbnail Fig. A.3

Eccentricity as a function of semimajor axis for a planetary system orbited by a secondary star of 1 M after 0.5 Myr of integration using Nemesis (big black bullets) and the single 8th order symplectic integrator in Huayno (smaller white bullets). The final eccentricity of the planets in the direct integration and the component method are indistinguishable in the figure, with an absolute mean error < 2 × 10−4 for each of the planets.

Open with DEXTER

For both integrators, we selected the eighth order symplectic integrator in Huayno. The two codes communicate using a Nemesis time step of dtNemesis = 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. A.3 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 dtNemesis = 100 yr gives satisfactory results in terms of accuracy, precision, and performance.

A.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 dtNemesis = 100 yr. The resulting evolution of the energy error is presented in Fig. A.4.

thumbnail Fig. A.4

Total energy error as a function of time for a validation simulation consisting of a single star orbiting a system of 5 planets. The energy error of our method (in red) is compared to the results obtained using a 4th order Hermite code for all particles (green). The time evolution of the energy error is more erratic in the Nemesis method because of the close interactions of the orbiting star. The overall error, however, remains rather constant over a long timescale, whereas for the Hermite method the energy error is smoother but clearly grows with time.

Open with DEXTER

In Fig. A.4 we show the results of the two calculations, one with a fourth order Hermite integrator (green), which is notsymplectic. 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. A.4, 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.

References

  1. Aarseth, S. J. 1985, Multiple Time Scales (Cambridge: Academic Press), 377 [Google Scholar]
  2. Belleman, R. G., Bédorf, J., & Portegies Zwart, S. F. 2008, New Astron., 13, 103 [NASA ADS] [CrossRef] [Google Scholar]
  3. Best, W. M. J., Liu, M. C., Dupuy, T. J., & Magnier, E. A. 2017, ApJ, 843, L4 [NASA ADS] [CrossRef] [Google Scholar]
  4. Blaauw, A. 1961, Bull. astr. Inst. Netherlands, 15, 265 [NASA ADS] [Google Scholar]
  5. Boekholt, T., & Portegies Zwart, S. 2015, Comput. Astrophys. Cosmol., 2, 2 [NASA ADS] [CrossRef] [Google Scholar]
  6. Cai, M. X., Kouwenhoven, M. B. N., Portegies Zwart, S. F., & Spurzem, R. 2017, MNRAS, 470, 4337 [NASA ADS] [CrossRef] [Google Scholar]
  7. Cai, M. X., Portegies Zwart, S., & van Elteren, A. 2018, MNRAS, 474, 5114 [NASA ADS] [CrossRef] [Google Scholar]
  8. Cassan, A., Kubas, D., Beaulieu, J.-P., et al. 2012, Nature, 481, 167 [NASA ADS] [CrossRef] [Google Scholar]
  9. Davies, M. B., Adams, F. C., Armitage, P., et al. 2014, Protostars and Planets VI (Tucson: University of Arizona Press), 787 [Google Scholar]
  10. Delorme, P., Gagné, J., Malo, L., et al. 2012, A&A, 548, A26 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  11. Eisenstein, D., & Hut, P. 2011, ApJ, 498, 137 [NASA ADS] [CrossRef] [Google Scholar]
  12. Fujii, M., Iwasawa, M., Funato, Y., & Makino, J. 2007, PASJ, 59, 1095 [NASA ADS] [Google Scholar]
  13. Fukushige, T., & Heggie, D. C. 1995, MNRAS, 276, 206 [NASA ADS] [Google Scholar]
  14. Gagné, J., Faherty, J. K., Cruz, K., et al. 2014a, ApJ, 785, L14 [NASA ADS] [CrossRef] [Google Scholar]
  15. Gagné, J., Lafrenière, D., Doyon, R., et al. 2014b, ApJ, 792, L17 [NASA ADS] [CrossRef] [Google Scholar]
  16. Gagné, J., Lafrenière, D., Doyon, R., Malo, L., & Artigau, É. 2014c, ApJ, 783, 121 [NASA ADS] [CrossRef] [Google Scholar]
  17. Gagné, J., Burgasser, A. J., Faherty, J. K., et al. 2015, ApJ, 808, L20 [NASA ADS] [CrossRef] [Google Scholar]
  18. Gahm, G. F., Grenman, T., Fredriksson, S., & Kristen, H. 2007, AJ, 133, 1795 [NASA ADS] [CrossRef] [Google Scholar]
  19. Gaudi, B. S. 2012, ARA&A, 50, 411 [NASA ADS] [CrossRef] [Google Scholar]
  20. Gould, A., & Yee, J. C. 2013, ApJ, 764, 107 [NASA ADS] [CrossRef] [Google Scholar]
  21. Goulinski, N., & Ribak, E. N. 2018, MNRAS, 473, 1589 [NASA ADS] [Google Scholar]
  22. Hamers, A. S., & Portegies Zwart, S. F. 2018, ApJL, submitted [arXiv:1810.11060] [Google Scholar]
  23. Hansen, B. M. S., & Murray, N. 2013, ApJ, 775, 53 [NASA ADS] [CrossRef] [Google Scholar]
  24. Haworth, T. J., Facchini, S., & Clarke, C. J. 2015, MNRAS, 446, 1098 [NASA ADS] [CrossRef] [Google Scholar]
  25. Heller, R. 2018, A&A, submitted [arXiv:1806.06601] [Google Scholar]
  26. Hunter, J. D. 2007, Comput. Sci. Eng., 9, 90 [Google Scholar]
  27. Hurley, J. R., & Shara, M. M. 2002, ApJ, 565, 1251 [NASA ADS] [CrossRef] [Google Scholar]
  28. Hut, P., Makino, J., & McMillan, S. 1995, ApJ, 443, L93 [NASA ADS] [CrossRef] [Google Scholar]
  29. Jänes, J., Pelupessy, I., & Portegies Zwart, S. 2014, A&A, 570, A20 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  30. Jílková, L., Portegies Zwart, S., Pijloo, T., & Hammer, M. 2015, MNRAS, 453, 3157 [NASA ADS] [CrossRef] [Google Scholar]
  31. Jílková, L., Hamers, A. S., Hammer, M., & Portegies Zwart, S. 2016, MNRAS, 457, 4218 [NASA ADS] [CrossRef] [Google Scholar]
  32. Kant, I. 1755, Allgemeine Naturgeschichte und Theorie des Himmels (Topeka, KS: Petersen) [Google Scholar]
  33. Kellogg, K., Metchev, S., Gagné, J., & Faherty, J. 2016, ApJ, 821, L15 [NASA ADS] [CrossRef] [Google Scholar]
  34. Kokubo, E., & Ida, S. 1998, Icarus, 131, 171 [NASA ADS] [CrossRef] [Google Scholar]
  35. Kokubo, E., & Ida, S. 2002, ApJ, 581, 666 [NASA ADS] [CrossRef] [Google Scholar]
  36. Kroupa, P. 2001, MNRAS, 322, 231 [NASA ADS] [CrossRef] [Google Scholar]
  37. Liu, M. C., Magnier, E. A., Deacon, N. R., et al. 2013, ApJ, 777, L20 [NASA ADS] [CrossRef] [Google Scholar]
  38. Liu, M. C., Dupuy, T. J., & Allers, K. N. 2016, ApJ, 833, 96 [NASA ADS] [CrossRef] [Google Scholar]
  39. Luhman, K. L. 2014, ApJ, 786, L18 [NASA ADS] [CrossRef] [Google Scholar]
  40. Luhman, K. L., Adame, L., D’Alessio, P., et al. 2005, ApJ, 635, L93 [NASA ADS] [CrossRef] [Google Scholar]
  41. Makino, J., & Aarseth, S. J. 1992, PASJ, 44, 141 [NASA ADS] [Google Scholar]
  42. Malmberg, D., Davies, M. B., & Heggie, D. C. 2011, MNRAS, 411, 859 [NASA ADS] [CrossRef] [Google Scholar]
  43. Marsh, K. A., Kirkpatrick, J. D., & Plavchan, P. 2010, ApJ, 709, L158 [NASA ADS] [CrossRef] [Google Scholar]
  44. McMillan, S. L. W. 2014, AAS/Division of Dynamical Astronomy Meeting, 45, 303.01 [NASA ADS] [Google Scholar]
  45. Mróz, P., Udalski, A., Skowron, J., et al. 2017, Nature, 548, 183 [NASA ADS] [CrossRef] [Google Scholar]
  46. Mróz, P., Udalski, A., Bennett, D. P., et al. 2019, A&A, 622, A201 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  47. Nitadori, K., & Makino, J. 2008, New Astron., 13, 498 [NASA ADS] [CrossRef] [Google Scholar]
  48. Oliphant, T. E. 2006, A Guide to NumPy (USA: Trelgol Publishing), 1 [Google Scholar]
  49. Pacucci, F., Ferrara, A., & D’Onghia, E. 2013, ApJ, 778, L42 [NASA ADS] [CrossRef] [Google Scholar]
  50. Pelupessy, F. I., & Portegies Zwart, S. 2013, MNRAS, 429, 895 [NASA ADS] [CrossRef] [Google Scholar]
  51. Pelupessy, F. I., Jänes, J., & Portegies Zwart, S. 2012, New Astron., 17, 711 [NASA ADS] [CrossRef] [Google Scholar]
  52. Pelupessy, F. I., van Elteren, A., de Vries, N., et al. 2013, A&A, 557, A84 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  53. Perets, H. B., & Kouwenhoven, M. B. N. 2012, ApJ, 750, 83 [NASA ADS] [CrossRef] [Google Scholar]
  54. Plummer, H. C. 1911, MNRAS, 71, 460 [NASA ADS] [CrossRef] [Google Scholar]
  55. Portegies Zwart S. 2011, Astrophysics Source Code Library [record ascl:1107.007] [Google Scholar]
  56. Portegies Zwart S. F. 2016, MNRAS, 457, 313 [NASA ADS] [CrossRef] [Google Scholar]
  57. Portegies Zwart, S., & Boekholt, T. 2014, ApJ, 785, L3 [NASA ADS] [CrossRef] [Google Scholar]
  58. Portegies Zwart, S. F., & Jílková, L. 2015, MNRAS, 451, 144 [NASA ADS] [CrossRef] [Google Scholar]
  59. Portegies Zwart, S., & McMillan, S. 2018, Astrophysical Recipes: the Art of AMUSE (Bristol,UK: Iop Publishing Ltd) [Google Scholar]
  60. Portegies Zwart, S. F., & Takahashi, K. 1999, Celest. Mech. Dyn. Astron., 73, 179 [NASA ADS] [CrossRef] [Google Scholar]
  61. Portegies Zwart, S. F., & Verbunt, F. 1996, A&A, 309, 179 [NASA ADS] [Google Scholar]
  62. Portegies Zwart, S. F., & Yungelson, L. R. 1998, A&A, 332, 173 [NASA ADS] [Google Scholar]
  63. Portegies Zwart, S. F., & Verbunt, F. 2012, Astrophysics Source Code Library [record ascl:1201.003] [Google Scholar]
  64. Portegies Zwart, S., McMillan, S., Harfst, S., et al. 2009, New Astron., 14, 369 [NASA ADS] [CrossRef] [Google Scholar]
  65. Portegies Zwart, S., McMillan, S. L. W., van Elteren, E., Pelupessy, I., & de Vries, N. 2013, Comput. Phys. Commun., 183, 456 [NASA ADS] [CrossRef] [Google Scholar]
  66. Portegies Zwart, S., Torres, S., Pelupessy, I., Bédorf, J., & Cai, M. X. 2018a, MNRAS, 479, L17 [NASA ADS] [CrossRef] [Google Scholar]
  67. Portegies Zwart, S., Pelupessy, I., van Elteren, A., Wijnen, T. P. G., & Lugaro, M. 2018b, A&A, 616, A85 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  68. Portegies Zwart, S., van Elteren, A., Pelupessy, I., et al. 2018c, AMUSE: the Astrophysical Multipurpose Software Environment [CrossRef] [Google Scholar]
  69. Rein, H., & Liu, S.-F. 2012, A&A, 537, A128 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  70. Rein, H., & Tamayo, D. 2015, MNRAS, 452, 376 [NASA ADS] [CrossRef] [Google Scholar]
  71. Schneider, A. C., Cushing, M. C., Kirkpatrick, J. D., et al. 2014, AJ, 147, 34 [NASA ADS] [CrossRef] [Google Scholar]
  72. Schneider, A. C., Windsor, J., Cushing, M. C., Kirkpatrick, J. D., & Wright, E. L. 2016, ApJ, 822, L1 [NASA ADS] [CrossRef] [Google Scholar]
  73. Sumi, T., Kamiya, K., Bennett, D. P., et al. 2011, Nature, 473, 349 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  74. Teachey, A., Kipping, D. M., & Schmitt, A. R. 2018, AJ, 155, 36 [NASA ADS] [CrossRef] [Google Scholar]
  75. Toonen, S., Nelemans, G., & Portegies Zwart, S. 2012, A&A, 546, A70 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  76. Udalski, A., Szymanski, M. K., Soszynski, I., & Poleski, R. 2008, Acta Astron., 58, 69 [NASA ADS] [Google Scholar]
  77. van Rossum, G. 1995, Extending and embedding the Python interpreter, Report CS-R9527 [Google Scholar]
  78. Veras, D. 2016, R. Soc. Open Sci., 3, 150571 [NASA ADS] [CrossRef] [Google Scholar]
  79. Veras, D., & Moeckel, N. 2012, MNRAS, 425, 680 [NASA ADS] [CrossRef] [Google Scholar]
  80. Veras, D., & Raymond, S. N. 2012, MNRAS, 421, L117 [NASA ADS] [CrossRef] [Google Scholar]
  81. Veras, D., Eggl, S., & Gänsicke, B. T. 2015, MNRAS, 451, 2814 [NASA ADS] [CrossRef] [Google Scholar]
  82. Vicente, S. M., & Alves, J. 2005, A&A, 441, 195 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  83. Vorobyov, E. I., Steinrueck, M. E., Elbakyan, V., & Guedel, M. 2017, A&A, 608, A107 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  84. Winn, J. N., & Fabrycky, D. C. 2015, ARA&A, 53, 409 [NASA ADS] [CrossRef] [Google Scholar]
  85. Wright, E. L., Eisenhardt, P. R. M., Mainzer, A. K., et al. 2010, AJ, 140, 1868 [NASA ADS] [CrossRef] [Google Scholar]
  86. Zapatero Osorio, M. R., Béjar, V. J. S., Martín, E. L., et al. 2000, Science, 290, 103 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
  87. Zapatero Osorio, M. R., Béjar, V. J. S., & Peña Ramírez, K. 2013, Mem. Soc. Astron. It., 84, 926 [NASA ADS] [Google Scholar]
  88. Zapatero Osorio, M. R., Gálvez Ortiz, M. C., Bihain, G., et al. 2014, A&A, 568, A77 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
  89. Zheng, X., Kouwenhoven, M. B. N., & Wang, L. 2015, MNRAS, 453, 2759 [NASA ADS] [CrossRef] [Google Scholar]

1

sōlus lapis, means “lonely rock” in Latin.

All Tables

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.

Table 2

Comparison of the distribution of planets at the beginning and at the end of the simulation.

Table 3

Parameters for planetary systems in which one planet was ejected to a larger distance (> 800 au) from its host star.

Table 4

Listing of systems that formed by the abduction of a planet from another star.

Table 5

Orbital elements of the merging planets.

Table 6

Orbital elements of the planets resulting from a merger.

All Figures

thumbnail Fig. 1

Diagram of Nemesis method. (A) A particle is an individual object or a subsystem consisting of multiple individual objects. In this study, the individual objects are either stars or planets. (B) The gravitational force on a particle is the sum of the force from other particles [1] and the forces from the individual objects in the subsystems.(C) The gravitational force on individual objects is the sum of the force from the particles [1] and the forces from the other individual objects in the containing subsystem [2], but not from those individual objects in other subsystems [3]. The forces from A.1 and B.2 are each in a self-contained system and can be calculated in an N-body code; the forces in A.2 and B.2 are connected to the self-contained systems and are evolved with a leapfrog algorithm.

Open with DEXTER
In the text
thumbnail Fig. 2

Projected view of the simulated star cluster at t = 0 Myr (initial conditions adopted by Portegies Zwart 2016; left panel), at t = 1 Myr (middle panel and the adopted initial conditions), and at t = 11 Myr (right panel, our final conditions). Stars are red bullets, single free floating planets black triangles.

Open with DEXTER
In the text
thumbnail Fig. 3

Eccentricity as a function of the semimajor axis the planets that survive up to an age of t = 11 Myr.

Open with DEXTER
In the text
thumbnail Fig. 4

Histogram for the number of systems with a certain number of planets. The dotted curve gives the initial distribution with either 4, 5, or 6 planets per star. The final conditions for the simulation without stellar dynamics are identical to this initial distribution. The distribution of the simulation in which we included the stellar encounters at an age of 11 Myr is presented as the solid curve with slanting lines.

Open with DEXTER
In the text
thumbnail Fig. 5

Histogram of the number of systems with a certain number of lost planets. Of the original 500 planetary systems the majority (386) do not lose any planets. Only 25 systems lose 1 or 2 planets and 102 systems lose 3 or more planets.

Open with DEXTER
In the text
thumbnail Fig. 6

Mean planet mass as a function of semi major axis (in a moving bin of 50 planets). The initial (at 1 Myr, in black) and the final (at 11 Myr in red) mean mass only differ slightly. The mean mass, 1.4 MJupiter, is depicted with a green horizontal line.

Open with DEXTER
In the text
thumbnail Fig. 7

Relative difference between the cumulative distributes of the masses of the bound planets initially and at 11 Myr. Positive values indicate an overabundance at the end of the simulation. Initially, the mean planet mass is 1.398 ± 1.05 MJup at 11 Myr the mean mass is only fractionally different at 1.404 ± 4.191 MJup; the latter value is indicated by the vertical green line.

Open with DEXTER
In the text
thumbnail Fig. 8

Eccentricity as a function of the semimajor axis for captured planets (black diamonds) and migrated planets with semimajor axis larger than 800 au (red dots) at an age of t = 11 Myr. The mean and standard deviation for both sets are also plotted. The mean orbital elements for the captured planets is ac = 1539 ± 824 au and ec = 0.6 ± 0.2, and ac = 1141 ± 258 au and ec = 0.8 ± 0.2 for the migrated planets.

Open with DEXTER
In the text
thumbnail Fig. 9

Number of free-floating planets (Nfp) as a function of time. The solid curve (black) indicates all free planets; the red curve indicates the subset of free floaters that also escape the cluster.

Open with DEXTER
In the text
thumbnail Fig. 10

Cumulative distribution (normalized) of the velocity of planets (black curves) and stars (red curves) at an age of 11 Myr. Planets and stars bound to the cluster are plotted with a thin line; the thick curves indicate the unbound objects.

Open with DEXTER
In the text
thumbnail Fig. 11

Cumulative distribution of the masses of all planets (red), the free-floating planets that are bound to the cluster (green), and those unbound from the cluster (blue). These three curves are statistically indistinguishable. The dotted curve indicates the mass distribution of 16 observed potential free-floating planets from Luhman et al. (2005); Marsh et al. (2010); Zapatero Osorio et al. (2000); Delorme et al. (2012); Liu et al. (2013); Gagné et al. (2014a,b,c, 2015); Schneider et al. (2014, 2016); Luhman (2014); Liu et al. (2016); Kellogg et al. (2016). For a different comparison, we introduce a lower mass cutoff to the initial sample of planets of 2 MJup and compare this with the observed sample (black).

Open with DEXTER
In the text
thumbnail Fig. 12

Number of planets that became unbound from their host star as a function of time. The number of planets that escaped their host within 0.5 Myr following a strong encounter (within 1500 au, red curve) is about twice as large as the planets thatescape without evidence of having experienced a strong encounter (black curve). The dotted black curves indicate thedependency on the timescale within which a strong encounter is supposed to lead to ejected planets; the lower curves indicate the cumulative distribution for planets that are liberated within 1 Myr of a close encounter, whereas the upper curve is for 0.2525 Myr).

Open with DEXTER
In the text
thumbnail Fig. 13

Number of planets that escape from their host star as a function of the time between a close encounter (within 1500 au) and the moment of escape. The majority of the planets escape promptly upon an encounter, but a considerable number require more time, up to about a million years.

Open with DEXTER
In the text
thumbnail Fig. A.1

Number of the initial intact planetary systems as a function of the Nemesis time step. The chosen time step of 100 yr is shown as a green vertical line. The time step is not optimal for this criterion, but was chosen as it gives better accuracy and a higher computational speed.

Open with DEXTER
In the text
thumbnail Fig. A.2

Sizeof the largest subsystem as a function of dtNemesis. after 0.1 Myr of evolution. The vertical green line indicates the adopted value of 100 yr.

Open with DEXTER
In the text
thumbnail Fig. A.3

Eccentricity as a function of semimajor axis for a planetary system orbited by a secondary star of 1 M after 0.5 Myr of integration using Nemesis (big black bullets) and the single 8th order symplectic integrator in Huayno (smaller white bullets). The final eccentricity of the planets in the direct integration and the component method are indistinguishable in the figure, with an absolute mean error < 2 × 10−4 for each of the planets.

Open with DEXTER
In the text
thumbnail Fig. A.4

Total energy error as a function of time for a validation simulation consisting of a single star orbiting a system of 5 planets. The energy error of our method (in red) is compared to the results obtained using a 4th order Hermite code for all particles (green). The time evolution of the energy error is more erratic in the Nemesis method because of the close interactions of the orbiting star. The overall error, however, remains rather constant over a long timescale, whereas for the Hermite method the energy error is smoother but clearly grows with time.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.