The origin and evolution of wide Jupiter Mass Binary Objects in young stellar clusters

The recently observed population of 540 free-floating Jupiter-mass objects, including 40 dynamically soft pairs in the Trapezium cluster have raised interesting questions on their formation and evolution. We test various scenarios for the origin and survivability of these free floating Jupiter-mass objects and Jupiter-mass Binary Objects (JuMBOs) in the Trapezium cluster. The numerical calculations are performed by direct N-body integration of the stars and planets in the Trapezium cluster starting with a wide variety of planets in various configurations. We discuss four models: SPP, in which selected stars have two outer orbiting Jupiter-mass planets; SPM, where selected stars are orbited by Jupiter-mass planet-moon pairs; ISF in which JuMBOs form in situ with the stars, and FFC, where we introduce a population of free-floating single Jupiter-mass objects, but no initialized binaries. Models FFC and SPP fail to produce enough JuMBOs. Models SPM can produce sufficient JuMBOs, but requires unusually wide orbits for the planet-moon system around the star. The observed JuMBOs and free-floating Jupiter-mass objects in the Trapezium cluster are best reproduced if they formed in pairs and as free-floaters together with the other stars in a smooth (Plummer) density profile with a virial radius of 0.5pc. A fractal stellar distribution also works, but requires relatively recent formations (>0.2Myr after the other stars formed) or a high (50%) initial binary fraction. This would make the primordial binary fraction of JuMBOs even higher than the already large observation fraction of 8%. The fraction of JuMBOs will continue to drop with time, and the lack of JuMBOs in Upper Scorpius could then result in its higher age, causing more JuMBOs to be ionized. We then also predict that the interstellar density of Jupiter-mass objects (mostly singles with 2% lucky surviving binaries) is 0.05/pc$^{3}$.


Introduction
Recently [1] reported the discovery of 42 Jupiter-Mass Binary Objects (JuMBOs) in the direction of the Trapezium cluster.Their component masses range between 0.6 M Jup and 14 M Jup and they have projected separations between 25 au and 380 au.Two of these objects have a nearby tertiary Jupiter-mass companion.They also observed a population of 540 single objects in the same mass range.This discovery initiates discussions on the origin and survivability of weakly bound Jupitermass pairs in a clustered environment.
Free-floating Jupiter-mass objects have already been detected, with the first sightings in the direction of the Trapezium cluster more than twenty years ago [2][3][4].Since then, many more have been found, for example, in the young clustered environment of Upper Scorpius [5], and through gravitational microlensing surveys in the direction of the Galactic bulge [6].Their abundance may be as high as 1.9 +1. 3  −0.8 per star [6], although a considerable fraction of these could be in wide orbits around a parent star, or have masses ≪ M Jup .
The origin of these free-floating planets has been debated in [7].Generally, star formation via the collapse of a molecular cloud through gravitational instability leads to objects considerably more massive than Jupiter [8,9], and in disks, planets tend to form with lower masses.The formation of extremely low-mass stars, however, remains an active field of research.Theory developed by [10] introduces the possibility for in situ formation of objects with Jupiter masses.Even so, the current consensus remains that the large population of Jupiter-mass free-floaters arises from fully packed planetary systems [11], or as ejectees from their birth planetary system following encounters with other stars in the cluster [12].Single Jupiter-mass free-floating objects then originally formed in a disk around a star to become single later in time [12][13][14][15][16][17].The number of super-Jupiter mass free-floating objects formed in this way are predicted to be on the order of one (∼ 0.71) per star [12], but lower-mass free-floaters orphaned this way may be much more abundant [15]; The origin of relatively massive free-floaters through dynamical phenomena is further complicated by the tendency for lower-mass planets to be more prone to ejections [17][18][19][20].
Explaining the observed abundance and mass-function of single free-floating Jupiter-mass objects is difficult.In particular, the large population of objects in Upper Scorpius challenges the formation channels.The recent discovery of a large population of paired free-floaters complicates matters even further and puts strong constraints on their origin.So far, binary free-floating planetary-mass objects have been rare, and were only discovered in tight (few au) orbits [21], including: • 2MASS J11193254-1137466 AB: a 5 to 10 M Jup primary in a a = 3.6 ± 0.9 au orbit [22].
Such tight pairs could have formed as binary planets (or planet-moon pairs) orbiting stars, before being dislodged from their parents [27].If only a few of these objects were discovered in tight orbits, such an exotic scenario would explain their existence reasonably well, but the discovery of a rich population of 42 wide JuMBOs [1] requires a more thorough study of their origin.
Assuming a dynamical history, we perform direct N-body simulations of a Trapezium-like star cluster with primordial Jupiter-mass objects (JMO) and JuMBOs.Our simulations focus on four models that could explain the abundance, and properties (mass and separation distributions) of the observed JMOs in the cluster.Alternative to forming in situ (scenario ISF ), one can naively imagine three mechanisms to form JuMBOs. [28] argued that hierarchical planetary systems could explain these binaries.Here, the outer two planets get stripped by a passing star during a close encounter.The two ejected planets would lead to a population of free-floating planets, but could also explain the observed population of JuMBOs.We call this scenario SPP (for star planetplanet).
Another possibility is that JuMBOs result from the ejection of planet-moon pairs (or binary planets) originally orbiting some star.We call this scenario SPM, for star planet-moon.Finally, we explore the hypothesis whether a sufficiently large population of free-floating JMOs could lead to a population of JuMBOs by dynamical capture of one JMO by another.We call this scenario F F C (free-floating capture).A similar scenario was proposed by [29] for explaining very wide stellar pairs, but the model also works for wide planetary orbits [30,31] We start by discussing some fundamental properties of the environmental dynamics in section 2, followed by a description of models in section 3, the numerical simulations to characterize the parameters of the acquired JuMBOs in section 4, and the resulting occurrence rates in section 5. We conclude in section 6.We initialise our cluster using parameters found by [32], who numerically modelled disk-size distributions and concluded that the Trapezium cluster was best reproduced for a cluster containing ∼ 2500 stars with a total mass of ∼ 900 M ⊙ and a half-mass radius of ∼ 0.5 pc.The results were inconsistent with a Plummer [32] distribution, but match the observations if the initial cluster density distribution represented a fractal dimension of 1.6 which we adopt here (see [33]).For consistency with earlier studies, we also perform our analysis for Plummer models.
Adopting a Plummer distribution of the Trapezium cluster (with virial radius r vir = 0.5 pc), the cluster core radius becomes r c ≃ 0.64r vir ∼ 0.32 pc with a core mass of 250 M ⊙ .This results in a velocity dispersion of v disp ≡ GM/(r 2 c + r 2 vir ) 1/2 ≃ 0.97 km/s.Assuming a mean stellar mass in the cluster core of 1 M ⊙ the unit of energy expressed in the kinematic temperature kT becomes ∼ 8 • 10 42 erg.
Observed JuMBOs are found in the mass range of about 0.6 M Jup to 14 M Jup and have a projected separation of 25 au to 380 au.The corresponding averages are r i j = 200 ± 109 au, ⟨M prim ⟩ = 4.73 ± 3.48 M Jup , and ⟨M sec ⟩ = 2.81 ± 2.29 M Jup .The median and 25 % to 75 % percentiles are r i j = 193.8+78.2  −114.1 au M prim = 3.67 +1.31 −1.57M Jup , and M sec = 2.10 +1.05 −1.05 M Jup .Assuming a thermal distribution in eccentricities, random projection and arbitrary mean anomalies, we expect the two objects to be bound in orbits with a typical semi-major axis of a expected ∼ 220 au.
To simplify our analysis, let us assume that the observed variation in projected distances between the two JMOs corresponds to an orbital separation, and express distances in terms of semimajor axis (see Appendix A for motivation).In practice, the differences between the projected separation and the actual semi-major axis of the orbit is small.Adopting a statistical approach, a thermal distribution in eccentricities and a random projection on the sky, the semi-major axis is statistically ∼ 1.2 times the projected separation.Keeping in mind that we do not definitively know whether the observed JuMBOs are truly bound, and even if they were, their underlying eccentricity distribution remains unknown.Nevertheless, in practice, this difference between projected separation and actual semi-major axis of a bound population is negligible compared to the 25% to 75% uncertainty intervals derived from the simulations.
To first order, the binding energy of JuMBOs ranges between ∼ 5•10 37 erg and 1.4•10 41 erg (or at most ∼ 0.02 kT ).This makes them soft upon an encounter with a cluster star.On average, soft encounters tend to soften these binaries even further [34], although an occasional soft encounter with another planet may actually slightly harden the JuMBO.
The hardest JuMBO, composed of two 14 M Jup planets in a 25 au orbit would be hard for another encountering object of less than 17 M Jup .For an encountering 1 M Jup object, a 25 au orbit would be hard only if the two planets are about three times as massive as Jupiter.This implies that JuMBOs are also generally soft for any encountering free-floating giant planet unless they are in tight enough orbits or the perturber is of low enough mass.Overall, independent of how tight the orbit, JuMBOs are expected to be relatively short-lived since they easily dissociate upon a close encounter with any other cluster member.The JuMBO ionization rate is then determined by the encounter probability, rather than the encounter parameters.
Once ionized they contribute to the population of free-floating single objects.Note that in the Trapezium cluster, even the orbits of 2MASS J11193254-1137466AB, and WISE 1828+2650 would be soft ( > ∼ 0.025 kT); they could be the hardest survivors of an underlying population.To further understand the dynamics of JuMBOs in a clustered environment, and to study the efficiency of the various formation scenarios we perform direct N-body calculations of the Trapezium star cluster with a population of JMOs in various initial configurations.
For each of our proposed models, ISF (in situ formation of JuMBOs ), SPP (JuMBOs formed via ejections of a host stars' outer planets), SPM (as planet-moon pairs orbiting a star), and F F C (mutual capture of free-floaters) we perform a series of N-body simulations with properties consistent with the Trapezium cluster.
Each cluster starts with 2500 single stars taken from a broken power-law mass-function [35] with masses between 0.08 M ⊙ and 30 M ⊙ distributed either in a Plummer sphere (model Pl) or a fractal distribution with a fractal dimension of 1.6 (model Fr).All models start in virial equilibrium.We run three models for each set of initial conditions, with a virial radius of 0.25 pc, 0.5 pc and 1.0 pc, called model R025, R050 and R100, respectively.We further assume stellar radii to follow the zero-age main sequence, and the radius of JMOs based on a density consistent with Jupiter (∼ 1.3 g/cc).
For our proposed models, we initialize a population of single JMOs and/or JuMBOs (binary JMOs).The single (and the primaries in planet pairs) are selected from a power-law mass function between 0.8 M Jup and 14 M Jup , which is consistent with the observed mass function [1].We fitted a power-law to the primary-planet mass function, which has a slope of α JuMBO = −1.2(considerably flatter than Salpeter's α Salpeter = −2.35).This choice follows from the observed mass distribution and is further motivated by the fact that the first dozen discovered free-floaters having a similarly flat mass function [3].Indeed, large statistics provided by gravitational microlensing surveys allowed a reliable measure of the slope, with α = −1.3+0.3 −0.4 [6].This mass function is slightly steeper than the slope derived for lower-mass ( < ∼ 1 M Jup ) free-floaters (α = −0.96+0.47 −0.27 [36]).
For each model, we have a special set of initial configurations.The clusters all have the same statistical representation.But the distribution of JMOs and JuMBOs varies per model.In figure 1 we sketch the various models.

Model F F C: JMOs as free-floating among the stars
For the models with free-floating JMOs, model F F C, we sprinkle the single objects, with mass taken from a power-law distribution of slope α = −1.2,into the cluster potential as single objects using the same initial distribution function as we used for the single stars (either Plummer or fractal).These models were run with ∼ 600 objects with a mass > 0.8 M Jup .
We performed additional runs with 10 4 free-floaters.Some of these runs have a different lower limit to the mass function, to keep the number of objects with a mass > 0.8 M Jup at ∼ 600 (assuming that lower-mass objects are unobservable).Each simulation is evolved for 1 Myr, after which we study the population of free floating JMOs and the population of JuMBOs.A few simulations were extended to 10 Myr, to study the long-term survivability of JuMBOs.

Model SPP: Star hierarchically orbited by two planets
For scenario SPP, 150 stars below and above 0.6 M ⊙ are selected to host planetary systems.The mid-mass point (of 0.6 M ⊙ ) represents twice the mean stellar mass in the mass function.
As the case for the model described in section 3.1, the mass of the primary planet, M prim , was chosen from a power-law distribution with slope α = −1.2.The mass of the secondary planet, M sec , was selected randomly from a thermal distribution between 0.2M prim and M prim .The more massive planet can therefore either be the inner planet or the outer one.A consequence of our mass-ratio distribution is that we have a slight preference for planets of comparable mass (as observed), and that we have a population of < ∼ 0.8 M Jup objects.This low-mass population contributes to ∼ 7.3 % of the total.

Submission
The distance from the first planet a 1 and the second planet a 2 (such that a 2 > a 1 ) are selected according to various criteria.The inner orbit a 1 was selected randomly between 25 au and 400 au from a flat distribution in a.The outer orbit, a 2 , was typically chosen to be five times larger than the inner planet's Hill radius.This guarantees the stability of the planetary systems if isolated.
Both planetary orbits are approximately circular, with a random eccentricity from the thermal distribution between circular and 0.02.The two planets orbit the star in a plane with a relative inclination randomly between −1 • and 1 • .The other orbital elements are randomly taken from their isotropic distributions.The system's orientation in space then gets randomized.We perform an additional series of simulations with pre-specified orbital separations for the two planets a 1 and a 2 , to follow the model proposed in [28].The results of these runs are presented in figure 11.

Model SPM: Star orbited by a pair of planet-mass objects
In the SPM models we initialize planet pairs (or planet-moon pairs) in orbit around a star.The masses of the stars, planets and moons are selected as in the SPP model.The planet-moon system's orbit was selected from a flat distribution in a between 25 au and 200 au, and with an eccentricity from the thermal distribution with a maximum of 0.02.
To warrant the stability of the star-planet-moon system, we choose an orbital separation such that the planet-moon pair stays within 1/3rd of its Hill radius in orbit around the star.The planetmoon system is randomly oriented.These systems tend to be dynamically stable, but some fraction may be subject to von Zeipel-Lidov-Kozai cycles [37][38][39].
With the adopted range of masses and orbital parameters, the time-scale for a cycle is of the order of a few Myr.Two JMOs in a circular 100 au orbit around a 1 M ⊙ star would lead to a circumstellar orbit of ∼ 3490 au.The von Zeipel-Lidov-Kozai cycle period of such a system is < ∼ 1.6 Myr, depending on the eccentricity of the planet-moon system around the star.

Model ISF : JMOs in weakly bound orbits
Primordial JuMBOs (model ISF ) are initialized with semi-major axis following a flat distribution between 25 au and 1000 au, an eccentricity from the thermal distribution between 0 and 1.The masses are selected as in model SPP although we investigate both the case where the mass-ratio follows a uniform and thermal distribution.Each system is subsequently rotated to a random orientation.The binaries are scattered in the cluster potential as single objects using the same initial distribution function as used for the stars.With the arrival of results, we expand our analysis on this scenario by exploring a more contrived parameter space (see section 4.5).
Figure 1 illustrates the four models with a schematic diagram.

The simulations
All calculations are performed using the 4-th order prediction-correction direct N-body integrator ph4 [40] through the Astrophysical MUltipurpose Software Environment, or AMUSE for short [41][42][43].The data files are stored in AMUSE formatted particle sets, and available via zenodo 10.5281/zenodo.10149241;the source code is available at github https://github.com/spzwart/JuMBOs.The script to run the simulations is fairly simple.It is essentially the same script from chapter 2 in [43], including the collision-detection stopping condition.Runs are performed with the default time-step parameter η = 0.03, which typically leads to a relative energy error < 10 −8 per step, and < ∼ 10 −5 at the end of the run.The fractal initial conditions with small virial radius (0.25 pc) are, not surprisingly, the most taxing.The relative energy error in these runs can be somewhat higher at times, but never exceed 10 −3 , which according to [44], suffices for a statistically reliable result.Snapshots were stored every 0.1 Myr and we analyse the data at an age of 1 Myr.Although incorporating stellar evolution, general relativity and the Galactic tidal field would be straightforward in the AMUSE framework, we decided to ignore those processes.We do not expect any stars to effectively lose much mass during the short duration of the simulations.Meanwhile, incorporating general relativity would have made the calculations expensive without much astrophysical gain, and the tidal field hardly has any influence on the close encounter dynamics in the cluster central portion.Moreover, in the observations, JuMBOs seem to be confined towards the cluster's central region.However, this is in part an observational selection effect.Namely, it is difficult to identify JMOs away from the cluster center due to intercluster extinction and the lower coverage at the periphery (S.Pearson, private communication).

Finding JuMBOs
Considering the average local kinetic energy of surrounding objects (stars and planets), JuMBOs are soft, complicating their identification in the numerical models.Generally, one considers hard binary pairs or multiples in direct N-body simulations, and finding soft pairs requires some extra effort.We search for JuMBOs by first finding every individual objects' nearest neighbors, and determining their binding energy.If an individual's nearest neighbour has another object identified as its own nearest neighbor, we adopt that as the close pair, and the initially selected object as a tertiary.Afterwards, we order the particles in terms of distance and binding energy, on which the eventual designation is based.
Instead of identifying JuMBOs as bound pairs, we also analyse the data only considering nearest neighbors using connected components.With this method we do not establish JuMBOs as bound objects, but as close pairs.This second method mimics observations, in which boundness cannot yet be established.
We denote single stars s, and planets p. Pairs of objects are then placed in parenthesis, for binary stars we write (s, s), a planetary system with one planet can be (s, p).A system with two planets then either becomes ((s, p), p), for a hierarchy of planets, or (s, (p, p)) for a planet-pair orbiting a star as in the SPM model.A JuMBO in this nomenclature becomes (p, p).

Results
The main results are presented in table 1 and table 2, but also in the more fine-tuned simulations presented in table 3, and table 4.
Given the large parameter space available regarding how to distribute JMO's among stars, we start by exploring part of this space with a selection of simulations in which we vary the way JMOs and JuMBOs are distributed among the stars.In addition, we cover a small portion of the cluster parameter space, including the virial radius and the density profile.All the other parameters we keep constant.
The results of these simulations are reported in section 4.1, and in the tables 1 and 2 in section 4.4.Section 4.5 further explores our favorite model, one in which JuMBOs form as isolated pairs together with JMOs and stars.

Distinguishing between the various models
Table 1 summarises the outcomes of our simulations for each scenario.Rows are named after their model designation followed by either the letter "Pl" for the Plummer model, or "Fr" for the Fractal model.The model name ends with the virial radius "R" in parsec, here R025 indicates 0.25 pc, R050 for 0.5 pc and R100 for 1 pc virial radius.
The SPP and F F C models systematically fail to reproduce the observed population of JuMBOs by a factor of 50 to 400.Changing the initial distribution in the semi-major axis of the inner orbit from a uniform distribution to a logarithmic distribution reduces the formation rate of JuMBOs even further.There are several systematic trends in terms of cluster density that depend on one's choice of a Plummer or fractal distribution, but it is not clear how these models can lead to JuMBOs.
Both models produce a considerable population of binary stars and single planetary systems, in particular the fractal distributions, where the typical fraction of dynamically formed binary stars is around 4 %, and the fraction of JMOs captured by a star is 0.7 % per star.Interestingly, models that already start with some paired configuration with a star (SPP and SPM) tends to produce more binaries and planetary systems than the models where planet-mass objects do not orbit stars (F F C).The high abundance of hierarchical multiple planets, ((s, p), p), in the SPP Pl and SPM Pl models reflects some of the initial conditions.These models also tend to produce a relatively rich population of single planet systems (s, p).
The only models that produce a considerable population of JuMBOs are the SPM and ISF models.The former however requires the binaries to orbit a ≳ 900 au, something we deem infeasible given observations and the size distribution of circum-tellar disks in the Trapezium cluster.For the latter, the Plummer distributions tend to produce a sufficient number of JuMBOs whereas the fractal model produces too few.
Figure 2 shows the cumulative distribution of the two-point correlation function, ζ(r i , r j ), between stars and JMOs for models ISF Pl R050 and ISF Fr R050.Note that for calculating the nearest mutual distance, JuMBOs are treated as two separate JMO.This gives rise to the left shoulder in the blue curve in figure 2. This shoulder is also visible in the star-star curve (black), whose height exceeds that of JMO-JMO above ∼ 100 au, and falls below within a few tens of au, where JMO-JMO pairs become more abundant since these correspond to the tightest initialised JuMBOs who are able to survive.The broader and larger JMO-JMO shoulder observed for the Plummer model results from its ability to better preserve wide JuMBOs.The distribution for the mutual distance between stars and JMOs does not show such a pronounced shoulder.Even so, Table 1: The average number of systems per simulation categorized into groups.The possible outcomes are the number of single stars (n s ), binaries (n (s,s) ), star orbited by a single planet (n (s,p) ), star orbited by two planets (n ((s,p),p)) ), single isolated planets (n p ) and JuMBOs (n (p,p) ).Note that we only list those objects with a mass > 0.8 M Jup .As a consequence, the total number of planetary mass objects do not always add up to 600.The number of stars also do not always add up to 2500 because of collisions and hierarchies not listed in the  at distances r i j > ∼ 10 4 au both distributions converge.The more violent nature of fractal distributions is also shown here in two ways.Namely, the fractal model exhibits a larger proportion of ≲ 10 3 au detections, implying numerous high-energy encounters.Moreover, the smaller proportion of JMO-JMO detections within ≲ 10 3 au highlights the tendency for the fractal model to more efficiently ionise JuMBOs.
Having established the existence of an overpopulation of nearest neigbors among stars and JMOs, we further explore the orbital characteristics of the surviving JuMBOs.Table 2 lists the median, the 25% and 75% quartiles for the JuMBO, (p, p), distribution in primary mass, secondary mass, semi-major axis and eccentricity.The table omits models that produce too few JuMBOs to calculate the median and quartiles.We compare observations to the medial values, citing the quartiles as they better portray the skewness of distributions.This is noticeable in the large difference between the range of the lower and the upper quartiles.
The primary masses produced in the models SPM and ISF , tend to be on the high side, but the secondary masses are in the observed range.SPM models tend to lead to orbits that are too tight.Omitting the orbits with a < 25 au JuMBOs does not improve the median orbital separation.In terms of the orbital separation (or projected distance) the best model seems to be ISF Plummer with a 0.5 pc virial radius, but the distributions are wide, and although the 1 pc ISF fractal model exhibits a low formation rate, it does reproduce the observed separation distribution.Meanwhile, models ISF Fr at an age of ∼ 50 kyr to ∼ 0.2 Myr compare quite favourably to the observations (see section 5.6) and could hint at the idea of JuMBOs forming later in the cluster evolution when the system has started to relax.
The eccentricities in model SPM Pl are generally smaller than in the SPM Fr models.Recall that here, planet-moon pairs started in nearly circular orbits.In the fractal models, the eccentricities are more effectively perturbed and thermalized, whereas in the Plummer models, this does not happen.In model ISF , JuMBOs start with higher average eccentricities, in which case the difference in eccentricity between the Plummer and fractal models is less pronounced.Later, in section 4.5, we also explore ISF models with initially circular orbits, and models including a

Submission
Table 2: Simulation results for the models that produce a sufficiently large population of JuMBOs to be considered feasible (mainly models SPM and ISF ).We present the median values and the quartile intervals for 25 % and 75 %.The observed median inter-JMO distance in the Trapezium cluster is r i j = 193.8population of free-floating JMOs alongside the JuMBOs .

Stellar and planetary collisions
We encountered several collisions in the simulation.The majority occur between two stars (83%), with the rest between a star and a JMO.Most collisions happen in the fractal models.In figure 3, we present the cumulative distribution of collisions in the models ISF Fr with a virial radius of 0.25 pc (solid blue), 0.5 pc (orange) and 1.0 pc (green), and for the equivalent models F F C Pl with the thin dash-dotted curves.The higher density naturally leads to more collisions, which tend to occur at earlier times.Plummer models typically yield fewer collisions.The only Plummer models in which collisions among stars were common is in model F F C with 64, 20 and 18 collisions on average per cluster, for those models with a viral radius of 0.25 pc, 0.5 pc and 1.0 pc, respectively.Interestingly, the fractal models from the same series and virial radii only experience 30, 17 and 14 collisions.It came as a bit of a surprise that models F F C Pl lead to so many collisions, throughout the investigation, no JMO-JMO collisions occured.

Runaway JuMBOs
In addition to mergers, ejection events also occurred.We adopt the classical definition of identifying runaway objects if their velocity with respect to the center of mass of the entire stellar system (dominated by the bound cluster) exceeds 30 km/s [45].
Not surprisingly, no JuMBO escaped the cluster with a high velocity, but there are some slow runaways, which were born in the cluster periphery and never experienced an encounter with a nearby star.Adding a tidal field to our calculations may cause this population to increase.
Single runaway stars and JMOs are also rare in our simulations.We attribute this to the lack of hard binaries in the initial conditions.We encounter on average one runaway JMO for each of the fractal models, and typically twice as many runaway stars.The JMOs, however, have on average a velocity of 110 km/s, whereas the stars escape with ∼ 63.3 km/s.

Fine-tuning for the binary fraction and separation distribution
With results suggesting an ISF origin, we further explore its consequences, and try to derive some of the earlier JuMBO properties to see if those are reconcilable with our understanding of planet and star formation in section 4.5.
Figure 4 shows how quickly the JuMBOs population decreases in time.The binary fraction among JMOs initially drops quickly (even exponentially in the fractal models), before slowing down to a survival fraction of 2% and 10% after 0.2 Myr in the 0.5pc and 1.0pc virial radius configurations respectively.For the latter model, the fraction of JuMBOs drops eventually to about 4 %; lower than the observed 8 %.
Figure 5 shows how the orbital separation of JuMBOs evolves in time and its dependency on the initial conditions.Note how quickly the distribution for model ISF Fr R050 drops with time, before converging in < ∼ 0.1 Myr to a median separation < 100 au.In all fractal runs, the distribution is much narrower and typically has a smaller mean than the observed properties.
Changing the initial eccentricity and/or semi-major axis distribution have a negligible effect on the survival rate of JuMBOs, and is unable to salvage the fractal models for explaining the observed population of JuMBOs.
Starting with a tighter population of JuMBOs will help delay their ionization, but hardly affects the eventual distribution in orbital separation.On the other hand, it will help in making the fraction of JuMBOs consistent with the observations (see 5).Such consistency is reached ∼ 0.05 Myr to 0.2 Myr after the start of the simulations.At this time, the fraction of JuMBOs to JMOs, as well as the separation distribution of the surviving JuMBOs are consistent with the observed population.We consider this a strong argument in favor of JuMBOs as late formed objects (see 5).
Although 0.05 Myr is a blink on cosmic scales, and it would be curious why the observed population all formed so near in time, one can envisage that JuMBOs take longer to form.In doing, JuMBOs will emerge after the initial violent phase of a young cluster and exhibit a less  rapid depletion in their population to that observed in figure 4. It follows that this would provide a wider range in time frame with which JuMBOs may emerge from.
In figure 6 we present the orbital distribution for JuMBOs in model ISF Fr050 (with and without free-floating JMOs) at an age of 50 kyr and likewise for ISF Fr100 at an age of 0.2 Myr.Shortly after the JuMBOs are introduced in the simulation, the widest orbits are being ionized, and by the time the cluster reaches an age of 1.0 Myr, the orbital distribution is skewed to values much lower than observations (see figure 7).At this age, however, the number of JuMBOs as well as their orbital separations matches the observed populations in the Trapezium cluster.

Introducing a population of free-floaters in model ISF
With the current results, we can further constrain the simulation parameters by repeating several calculations for better statistics and exploring other parts of the parameter space.We perform this analysis for model ISF since this model seems most promising in producing a sufficiently large population of JuMBOs in the range of observed parameters.
Once more, we restrict ourselves to a virial radius of 0.5 pc and 1.0 pc for both the Plummer and fractal models.Some models were run with an additional population of single JMO's (denoted with an 'ff' in the models name).The masses of these free-floaters are identical to the JuMBO primaries, and they are distributed in the same density profile.These models differ from F F C since they are initialised with a population of JuMBOs .Each model calculation was repeated 10 times to build up a more reliable statistical sample.Table 3 summarises our results.Models in the top segment of table 3 have more flexible initial conditions and are initialised with 500 JuMBOs , and when applicable, 500 free-floating JMOs.The JuMBOs themselves are initialised with semi-major axis taken from a uniform distribution and ranging between 0 ≤ a [au]≤ 1000.Primary masses are also taken from a flat distribution with the secondary having masses drawn from a mass ratio with uniform distribution.In all cases, the mass ranges between 0.8 ≤ M [M Jup ] ≤ 14.
Unlike the top segment, these 'O' models fix the number of JMOs to 600, (n JuMBO + n FF = 600) reflecting the observations [1].In cases where the number of free-floaters, n FF , is not zero, we take into account the typical survival rate of JuMBOs after 1 Myr (first column of the first segment of table 3) to get the correct proportion.
Finally, we run an additional series of simulations in which the primordial JuMBOs have an eccentricity between 0.0 and 0.2 sampled from a uniform distribution, rather than the usual thermal eccentricity distribution (model Fr 050ffOC where the 'C' denotes circular).Models ending by the letter 'L', such as model FR R050ffL, are those for which we extended the simulation time to 10 Myr.
Overall, we note that although the orbital distribution of the surviving JuMBOs in the Plummer models better coincide with observations than those in the fractal models (see figure 7), the tail to wider orbits are more pronounced, potentially highlighting problems with initial conditions.Reducing the maximum orbital separation from 1000 au to 400 au helps resolve this discrepancy, but it seems a bit unfair to tune the initial conditions to mimic the observed parameters if the majority of primordial JuMBOs survives.We apply this small but important change in models in the bottom segment (and identified with an 'O' as their final letter); these models have their JuMBOs initialised with 25 ≤ a [au]≤ 400 and masses drawn from a power-law with a slope of α JuMBO = −1.2(as was the case for PPM and SPM) and with a thermalised mass ratio.It would be relatively straightforward to acquire a satisfactory comparison between simulations and observations if we start the simulations with a population of JuMBOs that reflects the observations.Since JuMBOs in Plummer models only marginally evolve, here we focus mostly on fractal runs which exhibits a natural way of trimming the wide JuMBO population.
Adding a population of free-floaters to the Plummer models, makes considerable difference in the evolution of the JuMBO fraction, but eventually, after about 1 Myr their fraction converges to roughly the same value.They also influence the orbital distribution of JuMBOs since the interaction between a relatively tight JuMBO and a relatively low-mass JMO could be hard.This is reflected in the survival rate of JuMBOs.The consequence of this hardening process is also visible in figure 7 where the JuMBOs in the models that included free-floaters are, on average, Table 3: Statistics on the surviving JuMBOs.⟨...⟩ gives the median, while the ± denote the lower and upper quartiles.Col. 1: The fraction of JuMBOs present at the end of the simulation relative to the number initialised.Col. 3: The mass ratio of JuMBO systems.Col. 4: The primary mass of the JuMBO system.Col. 5: The semi-major axis of the JuMBO system.Col. 6: The eccentricity of the system.
Model tighter.Although one can reference to the fact that ionisation will also lead to a tighter orbital distribution, the enhanced survival rate between Pl 050 and Pl 050ff (and Fr 050 vs. Fr 050ff) suggests the ionization plays a secondary role to the observed trend.Although less obvious in fractal runs, we find that JuMBOs are more likely to survive if simulations include a population of free-floating JMOs.This tightening of the orbits establishes itself already at a very early age, as shown in figure 4 and figure 6.The JuMBOs in these runs, however, remain soft for encounters with any of the stellar-mass objects in the cluster.The hardening then reduces the interaction cross-section of the JuMBO , making it less vulnerable to any interaction, including ionization.
The fractal distribution efficiently prunes off any wide orbits since its violent nature provokes many encounters resulting in ionization.The tendency for JuMBOs to ionize at any encounter (JMO or stellar) is reflected by the little variation between runs of the same virial radius and the similar survival rates between fractal runs evolved till 1 Myr, and Fr 050ffL which, we recall, evolves until 10 Myr.The ease of ionisation is also reflected by the models with free-floating JMOs having smaller orbital separations on average.This trend is even most apparent in the Plummer models (see figure 7), where the orbital separations are ⟨a⟩ ∼ 296 au for Pl 050 compared to ⟨a⟩ ∼ 187 au for model Pl 050ff.
Overall, we find that free-floating JMOs lead to more JuMBOs.Their presence tends to yield tighter surviving systems through hardening, although they also effectively ionize the wider JuMBO systems, this latter case being no different to any stellar-JuMBO interaction.
Increasing the number of JMOs also enhances the chances of two free-floaters settling into a newly formed binary, on average, only 0.40 new JuMBO systems emerge in models Fr 050 compared to the 0.65 in Fr 050ff.For the Plummer models the increased JuMBO formation rate is even more striking by increasing from 1.75 for model Pl 050 to 3.15 for model Pl 050ff.This result sharply contrasts the F F C models, where, irrespective of how many JMOs were added to the cluster (up to 10 4 ), no JuMBOs formed in any of our simulations (see section 5.2 for a discussion).The primordial presence of JuMBOs mediates their further formation through interactions with JMOs.
In figure 8 we present the probability distribution of primary JuMBO mass and mass ratio, q, for Plummer and fractal models.The observed points (red) seem to cover a similar parameter space, but the overall distribution does not seem to be consistent with the observations.For the Plummer models this is, in principle, relatively simple to resolve, by using the observed distribution as input.For the fractal models such fine-tuning will be considerably harder because the fraction of survivors is only ∼ 4 %.
The overabundance of equal-mass JuMBOs at < ∼ 4 M Jup primary masses in the observations is striking, and hard to explain.If not just statistics or observational bias, this could indicate some unexplored formation mechanism.The high-mass-ratio population JuMBOs is further illustrated in figure 9, where we plot the observed JuMBO population (data from [1]).
Over plotted in figure 9 is a dotted curve which represents a separation between the high-mass ratio (red) and the low-mass ratio (blue) population.We empirically draw this curve which follows Here we adopted a distribution in q ranging from 0.2 to 1.0 following a thermal distribution such that the lowest mass primary has an equal-mass secondary, and a ratio of 0.2 for the most massive primary.The dotted curve seems to separate high-mass ratio JuMBOs from the low massratio cases.One interpretation for this curve is related to the binding energy, which is proportional to the total JuMBO mass and inversely proportional to orbital separation.The dotted curve then represents a constant binding energy for a ∼ 4.7 M Jup mass objects in a a = 1000 au orbit around a 1 M ⊙ -star.Or equivalently, two 10 M Jup -objects in a 25 au orbit.
We could continue calibrating the initial conditions for a more consistent comparison with the observations, but the global parameter tuning seems to indicate that either the Plummer or the fractal models with a 0.5 pc virial radius compare most favorably to the observations.At this point, we do not see a natural mechanism to remove the tail of very wide orbits among the JuMBOs, and it is a surprising that the observed orbital distribution seems to cut off rather sharply at about 400 au.
In figure 10 we show the primary mass function of the surviving JuMBOs for two models ISF Fr R050.The black curves are the initial distributions.Coloured lines of the same linestyle represent the corresponding final distribution.The difference between the initial and final mass primary function is minimal; one tends to lose some objects in the low-mass end, causing the curve to become flatter, and the mean mass to increase.Although the global trends are similar, the steeper mass function leads to a better comparison with the observations.We do not present mass functions from the other simulations, because the trends are essentially equivalent, once more highlighting that JuMBOs are prone to ionisation no matter with whom they interact.

Discussion
We explore the possible origin of the rich population of Jupiter-mass binary objects (JuMBOs) in the direction of the Trapezium cluster.The main problems in explaining the observations hides in the large number of Jupiter-mass objects (JMOs), their large binary fraction, and the wide separations.Assuming that they are bound, their orbits would be soft for any encounter in the cluster, theory dictates that they should not survive for more than a few hundred kyr (∼ 0.4 Myr for Plummer models, and < ∼ 0.1 Myr for the 0.5 pc ISF fractal models).The ease at which JuMBOs are ionized is illustrated in figure 4. It may be clear that the fractal models, due to their high frequency of strong encounters in the earliest phase of the cluster lifetime, have difficulty preserving wide planet-mass pairs.Plummer models are less dynamically interactive, and the fraction of JuMBOs remains much higher, where even relatively wide pairs can survive.
An alternative explanation for the large population of pairs among the free-floating JMOs might be that they form late.If JuMBOs only formed after ∼ 0.2 Myr, they have a better chance of surviving in the harsh cluster dynamical environment.Such a late formation would hardly affect the estimated mass of the objects, because the cooling curves used to estimate their mass from the observed temperature and luminosity are roughly flat at such a young age [3].
The binary fraction continues to drop well after 0.2 Myr for all models, and by the time the  cluster is 1 Myr old only ∼ 4% of the binaries in the fractal models survive.The survival fraction in the Plummer models is considerably higher.The fraction of JuMBO continues to drop, and by the time the cluster is ∼ 10 Myr the fraction of JuMBOs is < ∼ 2 %.Interestingly enough, [5], reported the detection of 70 to 170 single JMOs in Upper Scorpius, which has an age of about 10 Myr.None of the objects in Upper Scorpius is paired, although this observation could be biased in terms of missing close pairs due to relatively low resolution of the observations.Our calculations did not include primordial stellar binaries (or higher order systems), nor did we take the effect of stellar evolution and supernovae into account.Those processes may have a profound effect on the fraction of JuMBOs, tending to reduce, rather than increase their number.
Starting with a large population of ( > ∼ 600) single free-floating planetary-mass objects among the stars (but without JuMBOs ) would grossly overproduce the expected number of free-floaters, and consequently fail to reproduce the observed number of JuMBOs .This model, however, naturally leads to a mass-ratio distribution skewed to unity, as is observed.We consider this model undesirable by the lack of a large population of free-floating planets in the Trapezium cluster and no JuMBOs.This could indicate the existence of a large population of unobservable low-mass objects, but we consider this a rather exotic possibility.

Failure of model SPP: star with a hierarchical planetary system
The SPP model systematically fails to reproduce the observed population of JuMBOs by a factor of 50 to 400, leading us to rule out these scenarios as their possible origins.Changing the initial distribution in semi-major axis of the inner orbit from a uniform distribution to a logarithmic distribution reduces the formation rate of JuMBOs even further.
To further explore the failure of model SPP, we perform an additional series of simulations in the Plummer distribution with virial radii of 0.25 pc, and 0.50 pc.According to [28], the eventual orbital separation of the JuMBO would be consistent with the difference in orbital separation between the two planets when orbiting the star.We performed additional simulations with a mutual  The bullet points along each line correspond with the adopted orbital separation of the two planets (a 1 and a 2 ).The red symbols indicate an average orbital separations for the JuMBOs between 25 au and 380 au.The symbol sizes give the number of JuMBOs , in these simulation linearly scales with a maximum of 9. separation a 2 − a 1 = 100 au and a 2 − a 1 = 200 au, expecting those to lead to consistent results relative to the observed range in orbital separation for the JuMBOs, as was argued in [28].The other orbital parameters, for the planet masses, their eccentricities and relative inclinations, in these models were identical to the other SPP models.
The results of these simulations are presented in figure 11.The JuMBO-formation efficiency for these models peaks for an orbital separation a 1 > ∼ 1000 au, but steeply drops for smaller values of a 1 .From a total of 45 calculations with various ranges of a 1 and a 2 , 39 produced a total of 910 JuMBOs.Although results are not inconsistent with [28], we find ∼ 50 % wider distribution in separations than [28], who argued that the initial orbital distance a 2 − a 1 would be preserved.
The JuMBO formation rate is found to be orders of magnitude smaller than what [28] expected.They calculate the rate by means of 4-body scattering experiments, in which a star which hosts two equal-mass planets that orbit with semi-major axes a 1 and a 2 (a 1 < a 2 ), encounters a single star.Their largest cross-section of roughly a 2 1 is obtained if the encounter velocity 0.8v ⋆ /v 1 , for an encounter with a star with velocity v ⋆ .For an encounter at the simulated cluster's velocity dispersion the inner planet would then have an orbital separation of ∼ 900 au around a 1 M ⊙ star.
Note that an inner orbital separation of a 1 = 900 au for a 10 M Jup planet leads to a Hill radius of about 160 au.An orbit with a 2 = 1100 au, then is unstable.Still, even in the runs where we use these parameters, the total number of JuMBOs remains small compared to the number of free-floaters.Even if each star in the Trapezium cluster was born with two such planets at most one-third of the 42 observed JuMBOs could conceivably be explained, and the number of free-floating JMOs would run in the thousands.A stable hierarchical system of two Jupiter-mass planets in a circular orbit around a 1 M ⊙ star, would be dynamically stable if a 1 ≃ 120 au and a 2 ≃ 210 au (which is hard for an encounter with a JMO).
The results of the cross-section calculations performed by [28] are consistent with our direct 5.2 Failure of model F F C: Free-Floating Jupiter-Mass Objects Submission N-body simulation, but requires initial orbital separation too wide in comparison with a realistic population of inner-planetary orbits for JMOs.Indeed, observational constraints on the existence of > ∼ 900 au JMOs are quite severe, and we consider it unrealistic to have 300 out of 2500 stars hosting such wide planetary systems.This is further motivated when one considers the small sizes of the observed disks are smaller than 400 au in the Trapezium cluster [46].

Failure of model F F C: Free-Floating Jupiter-Mass Objects
In scenario, F F C, we initialize 600 to 10 4 single JMOs in a cluster of stars without JuMBOs (see section 3.1), expecting that some soft pairs form naturally through interactions with the stars.Soft binary formation, through three-body interactions is not expected to be very effective [47], but with a sufficiently large population, one might expect a few JuMBOs to form.
A JuMBO can form in models F F C, when two JMO and a single star occupy the same phase space volume.In such a scenario, the star can escape with the excess angular momentum and energy, leaving the two planetary-mass objects bound.The distance at which a JMO with mass m and a star with mass M can be considered bound can be estimated from the 90 • turn-around distance, which is r 90 = G(M + m)/v 2 .For our adopted clusters r 90 ≃ 900 au.Following [47], we estimate the probability of this to happen at ∼ 10 −2 per JMO per relaxation time.The relaxation time of our Trapezium model cluster is approximately t rlx ∝ N/(6 ln(N))t cross ∼ 64t cross .With a crossing time of about 1 Myr (roughly the crossing time for our 1 pc models) we expect ∼ 0.5 JuMBOs to form.The estimates by [47] adopted equal mass objects, but the more detailed numerical study, carried out by [48] arrives at a similar number of soft binaries.The latter study, however, focused on post-core collapsed clusters, which is not appropriate for our Plummer models, but more in line with our fractal models.Stellar sub-clumps collapse in the fractal models within ∼ 0.2 Myr, mimicking the post-collapse evolution as addressed in [48].Therefore, in principle, their model is appropriate for our F CC models.
Interestingly enough, the F F C models produce quite a rich population of stellar pairs (82.8 ± 6.6) and several cases where a JMO is captured by a star (6.5 ± 3.3) for the Plummer as well as for the fractal models.But no JuMBOs formed.A higher abundance of stellar pairs compared to planetary captures is somewhat unexpected.In figure 12 we present the cumulative distribution of the eccentricities found in several of our model calculations.Each of them is consistent with the thermal distribution (indicated as the black dotted curve).
The secondary masses in the stellar pairs that formed are statistically indistinguishable from the primary masses of the stars that captured a planet, and their orbits are wider; 182±67 au for the binaries and 288 ± 164 au for the captured JMOs.With the higher masses the binaries are roughly 100 times harder than the JMO captured systems; straddling the hard-soft boundary.

The long-term survivability of JuMBOs
To study the long-term survivability of JuMBOs we execute 10 runs for ISF Fr 050 until an age of 10 Myr.Our aim is to look at the survival of JuMBOs in older clusters, such as Upper Scoprius.Overall, the JuMBO survival rate decreases rapidly with a half-life < 1 Myr, the survivors have tighter orbits.The population of JuMBOs eventually settles at a population of dynamically hard pairs, in which case the mean orbital separation ⟨a⟩ < ∼ 20 au for two 3 M Jup objects.The hardness of these pairs is mostly the result of the decrease in the cluster density with time, rather than in the shrinking of the surviving JuMBOs.
In figure 13 we present the distribution in primary mass and mass ratio for simulation ISF Fr R050ffL at an age of 10 Myr (adopting the uniform primary mass function).An equivalent diagram for our 1 Myr runs is shown in the right panel of figure 8.The most likely survivors have high primary and secondary masses, as could be expected based on the the hardness of these systems.Note that here we adopted a uniform mass function for the primaries and secondaries scaling their masses through a uniformly sampled mass-ratio, q. 2), and to equal mass systems (0.2 ≤ q ≤ 1.0, sampled from a thermal distribution).

Observational Constraints
No matter the initial conditions, f surv and e barely changes, while a shows only marginal differences.In turn, the fractal models exhibit a natural tendency for trimming out wide binaries.Similarly, JuMBOs are found not to rely on their mass parameters.Indeed, no matter the distribution masses are drawn from, the evolution of the primary mass tends to flatten with a very weak trend favouring larger M prim (see figure 10).This is reflected in table 3 since models initialised with low primary masses (bottom segment) exhibit upper quartile ranges spanning larger values and the tendency for JuMBOs to skew rightwards in the M prim vs. q diagram shown in figure 14.
Although we aimed to reproduce the primary mass-function and mass ratio distribution, we clearly under-represent high mass primaries with a low mass ratio.For the Plummer models it will be relatively straightforward to reproduce both populations (the high mass ratio, as well as those with a low mass ratio) since JuMBOs in these runs are rarely ionized, but in the fractal models the survival rate is too low to still reflect the initial conditions.
In figure 15, we present the distribution of the dynamically formed JMO-star systems in semimajor axis and eccentricity.At any given instance in time, a typical fractal run has ⟨N⟩ ∼ 4 of these systems present in the cluster.The parameter space is widely covered, with signs of low eccentricity but very wide (a ≥ 700 au) binaries.The vast majority exhibit large eccentricities and semi-major axis, reflecting their dynamical origin.The non-negligible amount of these systems emerging provides an interesting prospect of detecting ultra-cold Jupiters orbiting stars that have recently fostered them.
Our slight preference for the fractal models stems from their natural consequence of producing higher-order hierarchical systems, which on rare occasions produce triple JMO systems (Jupiter Mass Triple Objects, JuMTOs) as was observed in [1], and its capability of removing wide JuMBOs.Contrariwise no triples were detected in Plummer models and after 1 Myr.Nevertheless, if formed in situ, the Plummer may have a sizable fraction of primordial triples survive: we are then back at weighing the relative importance of initial conditions versus dynamical evolution and emergence.

Assessment on the origin of JuMBOs
In the ISF Plummer models, depending on the clusters configuration, ∼ 50% − 70% of the JuMBOs are ionized within 1 Myr, compared to ∼ 96 % for the fractal models.The observed population of free-floaters and JuMBOs can then be reproduced if the cluster was born with half of its population of JMOs as free-floaters and the other half as pairs.The current observed primary and secondary masses of JuMBOs would reflect the conditions at birth, while the semi-major axis and eccentricity distributions would have been affected considerably by encounters with other cluster members.These processes tend to drive the eccentricity distribution to resemble the thermal distribution (probably with an excess of > ∼ 0.7 eccentricities [49]).The semi-major axes of the JuMBOs would have widened, on average by approximately 5% due to encounters with free-floating planet-mass objects.
Alternative to a Plummer initial stellar distribution we consider fractal distributions, which are also able to satisfactorily reproduce the observed populations.In the fractal models, > ∼ 90 % of the primordial JuMBOs become ionized, and in principle the entire observed populations of free-floating JMOs and JuMBOs can be explained by a 100% initial binarity among the JuMBOs.We then conclude that JMOs are preferentially born in pars with a rather flat distribution in orbital separations with a maximum of ∼ 400 au.Higher order multiplicity (JuMTOs and JuMQOs) form naturally from interactions between two or more JuMBOs in the fractal models.
This model (ISF Fr R050) satisfactorily explains the observed orbital separation distribution, with a ∼ 15% excess of systems with a separation > ∼ 400 au.We do not expect a rich population of orbits with separation smaller than the observed 25 au.Indeed, when processing results we found that fractal models whose JuMBOs are initialised with 10 ≤ a [au] ≤ 1000, ∼ 65% of JuMBOs have orbital separations larger than 25 au (∼ 90% when we restrict the initial We have a slight preference for the ISF fractal models with 0.5 pc virial radius because hierarchical triple JMOs form naturally in roughly the observed proportion (on average ∼ 4 triples among ∼ 40 pairs and ∼ 500 single JMOs).The singles then originate from broken-up pairs, and the triples form in interactions between two JMO pairs.The dynamical formation of soft triple JMOs is quite remarkable, and observational follow-up would be of considerable interest.
The mass function of single JMOs should resemble the combined mass functions of the primary and secondary masses of the JuMBOs since the ionization probability for a JuMBO does not depend on its mass, but instead on the chance encounter with a star.In this sense, the majority of JMOs would originate from ionized JuMBOs.

Fine tuning of model ISF R050
The short-periods and the small binary fraction of the fractal models could be salvaged if the JuMBOs form late.If the majority of the observed population formed > ∼ 0.2 Myr later than the stars, the cluster's density profile would already have been smoothed out, leading to fewer strong dynamical encounters.The cluster would somewhat resemble a JuMBO-friendly Plummer-like structure.Not only would that mediate the survival of JuMBOs but it also would allow them to preserve their orbital characteristics.
In table 4 we present a small experiment to support the argument of JuMBOs forming late.The table shows the fraction of surviving JuMBOs, and their orbital parameters.We omitted the eccentricities, as they are consistent with the thermal distributon.
The experiment starts with simulation model ISF Fr R050.Instead of running with JuMBOs we run with 900 test masses distributed in the same fractal structure as the stars.At various moments in time, we stop the run, and restart it with JuMBOs and JMOs.Each restart is initialized with 300 JuMBOs and 600 JMOs.The JuMBOs are taken from the same distribution functions for primary mass, secondary mass and orbital parameters from the models ISF ; Their maximum orbital separation is 400 au.
It should not come as a suprise that the fraction of surviving JuMBOs increases when they form later.We illustate this further in figure 16, where we plot the surviving fraction of JuMBOs as a function of the time they were introduced in the simulation.We started introducing them at the same time as the stars (left, t birth = 0 Myr), and as last point they were introduced together with the stars (right most point at t birth = 1 Myr).The left-most point essentially replicates our earlier ISF Fr configurations, and the right-most point reflects the initial conditions.
In figure 17, we present the orbital distribution for t birth = 0.2 of model ISF Fr R050 (red), and for model ISF Pl R050 for which the JuMBOs were initialized at birth (t birth = 0).The black curve shows the observed projected separation distribution of the JuMBOs from [1].
The consistency of the orbital separation in the fractal model in figure 17 is in part by construction, as this was also the input initial condition.But similar input parameters were adopted for the Plummer model and the fractal model, but in both these cases the JuMBOs were initialized together with the stars.

Conclusions
The discovery of 40 relatively wide pairs, two triples and 540 single JMOs in the Trapezium cluster emphasizes our limited understanding of low-mass star and high-mass planet formation.To derive characteristics for their origin we performed simulations of Trapezium-equivalent stellar clusters (2500 stars in a virialized 0.25 pc to 1.0 pc radius) with various compositions of JMOs and stars.
Models in which planets form in wide hierarchical circumstellar orbits (model SPP), as proposed by [28], produce many single free-floating planets, but insufficient numbers of pairs.The ratio of single to pairs of planet-mass objects in these models are too low by a factor of 50 to 400, irrespective of the initial stellar distribution function.
The models in which pairs of planetary-mass objects orbit stars in the form of a planet-moon system (or binary planets, model SPM), produce a sufficent number of free-floating planetary pairs, and cover the proper range of orbits.In particular the models that start with fractal initial conditions tend to produce a sufficent fraction of JuMBOs among free-floating objects O(0.1), which is close to the observed value of 0.078 ± 0.012.In the Plummer distribution, the number of stars that survive with at least one planet-mass objects is considerable.These cold Jupiters have a typical orbital separation of ⟨a⟩ = 382 ± 75 au, and rather high ⟨e⟩ = 0.74 ± 0.08 eccentricity.For the JuMBOs these models generally predict low-eccentricities (e < ∼ 0.4), whereas others lead to thermalized distributions (⟨e⟩ ∼ 0.6). Figure 17: Distribution of orbital separation for the observed JuMBOs (black) and those from the two models ISF Fr R050 with t birth = 0.2 Myr (red), and adopting t birth = 0 Myr we plot both ISF Pl R050 (blue) and ISF Fr R050.The fractal model compares well with the observed distribution (KS p-value = 0.939) compared to the Plummer and fractals distributions (KS p-value< 0.003).The blue and orange curves are also presented in figure 7, as the red dashed curves (but now in log-scale).
For model SPM to produce a sufficient number of JuMBOs it requires planet-moon pairs to form in > ∼ 900 au orbits around their parent star.Such wide orbits are exotic since the circumstellar disks observed in the Trapezium cluster tend to be smaller than 400 au.We therefore do not see how such wide planet-moon pairs can form around stars.If, however, a population of cold a > ∼ 100 au JMOs are found in the Trapezium cluster, we do consider this model a serious candidate for producing JuMBOs.Investigating some of the observed JuMBOs in the images published in [1], we get the impression that some JuMBOs may have nearby stars, but a thorough statistical study to confirm this correlation is necessary.Model SPM can easily be confirmed or ruled out by establishing the existence of those left over (or dynamically formed) planetary systems.
Ruling out models F F C, SPP, and possibly SPM, we are left with the simplest solution; JuMBOs form together with the single stars in the cluster.This model reproduces the observed rates and orbital characteristics (a bit by construction, though); it can also be used to further constrain the initial conditions of the cluster as well as the JuMBOs.
We tend to prefer model ISF with a 0.5 pc Plummer sphere because it can be tuned rather easily to reproduce the observed population of JuMBOs and JMOs.On the other hand, we also consider the equivalent fractal model (ISF Fr R050) a good candidate so long as JuMBOs were formed somewhat later in the clusters' evolution, once the rate of violent encounters reduced.If JuMBOs form within ∼ 0.1 Myr of the stars, our results diverge with observations, forming too few JuMBOs and with too tight orbits.Even so, one strong point of the fractal models is their natural ability to prune wide JuMBOs and their natural ability to form Jupiter-Mass Triple Objects (JuMTOs).
Single free-floating planetary objects were discovered in abundance before in the Upper Scorpius association (between 70 and 170 candidates) [5], but these were considered to be single free-floaters.With an age of about 11 Myr [5], Upper Scoprius is expected to be rich in single Jupiter-mass free-floating planets, but binaries will be rare as the majority will be ionized.We still could imagine that a few JuMBOs have survived until today.
Finally, we would like to comment briefly on the nature of the objects observed.We wonder that, if these objects formed in situ, and therefore not around a star, they would be deprived of a rocky core.JuMBOs and JMOs would then more resemble a star in terms of the structure, rather than a planet.This may have interesting consequences on their dynamics, their evolution, and when they encounter another star (collisions in our simulations are relatively frequent).In those terms, we also wonder to what degree the term "planet" is rectified at all, and maybe it is time to revive the IAU discussion on the definition of a planet.

Energy Consumption
For 8400 hours of CPU time and an energy consumption rate of 12 Watt/hr [54] and an emission intensity of 0.283 kWh/kg [55] for performing all calculations in this manuscript, we emitted ∼ 357 kg of CO2.This amounts to roughly 32.48 tree-years of Carbon sequestration.Enough to charge 43 426 smartphones.A Similarity between r i j and a The comparison between simulations and observations is somewhat hindered by the different perspectives.Whereas dynamicists prefer to use Kepler orbital elements, from an observational perspective such data is not always available.In our current study, we try to compare populations of binaries with observed objects.The latter are projected separations, which do not directly translate in orbital elements without full knowledge of the 6-dimension phase space of the orbit.We therefore have to compare projected separation with what we prefer to use, the semi-major axis of a bound two body orbit.Both panels in figure 18 motivate our choice of analysing results in terms of the semi-major axis given the similarity between the curves.In all cases, r i j exhibits longer tails at shorter separations/orbits.However, these differences are so small, especially in the fractal case, and considering the width of the distribution (illustrated with the large values for the quartile intervals) we can safely interchange between one and the other.In doing so, we assume that the observed projected separation of JuMBOs are equivalent to their semi-major axis, easing our discussion.

Figure 1 :
Figure 1: Illustration of the four configurations for JMOs in the stellar cluster.Stars are represented with yellow bullets, and JMOs in red.From top left to bottom right we have (as indicated): model F F C for the free-floating single planets; SPP for outer orbiting planets; SPM as bound planetmoon pair orbiting a star, and model ISF in situ formation of jumbos.

Figure 2 :
Figure 2: Mutual nearest neighbor distribution function for model ISF Fr R050 (solid lines) and ISF Pl R050 (dash-dotted lines).Our analysis only considers single objects such that JuMBOs are deconstructed into two JMOs.

Figure 3 :
Figure 3: Number of collisions as a function of time for model ISF Fr (solid lines) and for models F F C Pl (dash-dotted lines).

Figure 4 :
Figure4: Fraction of JuMBOs as a function of time for models ISF Pl R050 and ISF Pl R100 (left panel), ISF Fr R050 and ISF Fr R100 (right panel) In both cases, we also perform a calculation initialised with a population of single JMOs.The number of free-floating objects is the same as the number of primordial JuMBOs and their masses are taken from the primary mass function.Scatter points denote where in time the model is nearest to the 9% observed benchmark for fraction of JuMBOs relative to JMO free-floaters.

Figure 5 :Figure 6 :
Figure 5: Evolution of the median orbital separation and the 20% and 75% quartiles of the distribution for different fractal runs.Model ISF Fr 050 looks at the case where JuMBOs are initialised with 0 ≤ a [pc]≤ 1000 while ISF Fr 050O and ISF Fr 050OC reduce the range to 0 ≤ a [pc]≤ 400.In addition, Fr 050OC considers initially circular JuMBOs, with eccentricities ranging between 0 ≤ e ≤ 0.2.

Figure 8 :
Figure 8: Probability distribution of primary mass versus mass ratio for the simulated JuMBOs in model ISF Pl R050ff (left) and ISF Fr R050ff (right) The red crosses show the observed JuMBOs.The primary masses for these runs ware generated from a uniform distribution, which is still evident in the Plummer model (left), but in the fractal model the initial conditions are lost.

Figure 9 :
Figure 9: Distribution function of the observed JuMBOs, for primary mass, observed projected separation and mass-ratio (colors).The red bullet points indicate the high mass-ratio population, whereas the blue bullets show the low-mass ratio population.The dotted curve separating the two populations equals 10 3 au/(M prim (1 + q)) as indicated.

Figure 10 :
Figure 10: Distribution of the primary mass of the surviving JuMBO for models Fr 050 with an without an adapted mass function.The black curves give the initial mass distribution and the colored ones the mass distribution after 1 Myr.

Figure 11 :
Figure 11: The number of JuMBOs produced in model SPP, as fraction of the number of freefloating planets for various simulations starting with a Plummer sphere and virial radius of 0.5 pc.The bullet points along each line correspond with the adopted orbital separation of the two planets (a 1 and a 2 ).The red symbols indicate an average orbital separations for the JuMBOs between 25 au and 380 au.The symbol sizes give the number of JuMBOs , in these simulation linearly scales with a maximum of 9.

Figure 12 :Figure 13 :
Figure 12: Eccentricity distribution for several models (indicated in the lower right inset) for stellar binaries and captured planet-mass objects.Overplotted, for comparison, is the thermal distribution.

Figure 14 :
Figure14: Contours of constant density in the plane of primary mass versus mass-ratio for model ISF 050ff.Red crosses denote observed JuMBOs.In this model, the initial mass function and mass-ratio distributions for the JuMBOs were skewed to lower mass primaries (α = −1.2),and to equal mass systems (0.2 ≤ q ≤ 1.0, sampled from a thermal distribution).

5. 5 Figure 15 :
Figure15: Contours of constant density for semi-major axis and eccentricities for JMO-stellar systems for model ISF Fr 050ff.These systems are surprisingly common in our simulations, and they even tend to be relatively tight, though the actual orbits remain dynamically soft.

Figure 16 :
Figure16: Surviving fraction of JuMBOs for model ISF Fr R050 as a function of the moment they were introduced in the simulation (t birth ).The orange bullets present the actual measurements, and the curve is a least squares fit, represented by f JuMBOs ≃ 0.33 + 0.18 log 10 (t birth /Myr).

Table 4 :
Final distribution of fraction of JuMBOs and orbital parameters for models ISF Fr R050 but with the time of birth for the JuMBOs at t birth .The best combination of values of survival fraction, and orbital separation are achieved if JuMBOs are initialized somewhat later than the stars by between 0.2 Myr or 0.4 Myr.The initial fraction of JuMBOs over JMOs is typically 50 %.t birth /Myr f (p,p) ⟨M prim ⟩/M Jup ⟨M sec ⟩/M Jup axis to 25 ≤ a [au]≤ 400).This increases to ∼ 90% for Plummer models.The fraction of pairs among the wide systems is already high, and there are not enough single JMOs observed to accommodate this tight binary population, unless a considerable fraction of the single observed JMOs are in fact such tight binaries.