Dynamics of Binary Planets within Star Clusters

Abstract

We develop analytical tools and perform three-body simulations to investigate the orbital evolution and dynamical stability of binary planets within star clusters. Our analytical results show that the orbital stability of a planetary-mass binary against passing stars is mainly related to its orbital period. We find that critical flybys, defined as stellar encounters with energy kicks comparable to the binary binding energy, can efficiently produce a wide range of semimajor axes ($a$) and eccentricities ($e$) from a dominant population of primordially tight JuMBOs. Applying our results to the recently discovered Jupiter-Mass Binary Objects (JuMBOs) by the James Webb Space Telescope (JWST), our simulations suggest that to match the observed $\sim$9% wide binary fraction, an initial semimajor axis of $a_{0}=10$–20 au and a density-weighted residence time of $\chi\gtrsim 10^{4}$ Myr pc^-3 are favored. These results imply that the JWST JuMBOs probably formed as tight binaries near the cluster core.

UTF8gbsn

1 Introduction

With an estimated age of $<$1 Myr and a distance at $\approx$470 pc (Hillenbrand & Hartmann, 1998), the Trapezium cluster is one of the youngest star clusters near the Sun. It is located at the center of the Orion Nebular Cluster and has one of the densest stellar environments (McCaughrean & Stauffer, 1994). The extreme conditions of the Trapezium Cluster have made it an ideal target for studying the formation of stars and planetary-mass objects (PMOs), both theoretically and observationally (e.g., Prosser et al., 1994; Muench et al., 2008; Pearson & McCaughrean, 2023).

The recent discovery by Pearson & McCaughrean (2023) of 540 PMO candidates, including 40 Jupiter-Mass Binary Objects (JuMBOs), in the Trapezium cluster poses significant challenges to our conventional understanding of planet formation. In the core accretion model, giant planets form within circumstellar protoplanetary disks. The growth from planetesimals to a planetary core (several $M_{\oplus}$) can last up to a million years, followed by a gas accretion phase lasting up to $\sim$10 Myr (Lissauer, 1993; Pollack et al., 1996). This time span is already a factor of a few longer than the estimated age of the Trapezium Cluster. Furthermore, the timescale for ejecting planets out of embedded systems through planet–planet scatterings or stellar encounters is estimated to be another millions to hundreds of millions of years (Spurzem et al., 2009; Nesvorný & Morbidelli, 2012; Huang et al., 2022).

In addition, the abundance of these PMOs relative to stars is also surprisingly high. The Trapezium cluster is estimated to contain $\sim$2000 stellar members (Morales-Calderón et al., 2011), so the relative abundance of PMOs is $>$$20\%$. This is above the limit on Jupiter-mass free-floating planets (FFPs) that were derived from microlensing surveys (Mróz et al., 2017). It is also in tension with the substellar mass function derived from Euclid, albeit for a less dense environment (Martín et al., 2024).

What is even more surprising is the discovery of 40 JuMBOs (with binary masses of 2–20 $M_{\rm J}$) at wide separations ($d=25$–400 au), which contribute $\sim$$9\%$ of all PMOs discovered in the Trapezium Cluster (Pearson & McCaughrean, 2023). Planet-mass binaries have previously been identified. For example, Best et al. (2017) found 2MASS J11193254–1137466AB to be a pair of $\sim$4 $M_{J}$ objects with $d\approx 4$ au, and Beichman et al. (2013) found WISE 1828+2650 to be a pair of 3–6 $M_{J}$ objects with $d<0.5$ au. Oph 162225–240515 is a wide ($d\approx 242$ au) binary that was originally thought to be in the planetary-mass regime (Jayawardhana & Ivanov, 2006), but later studies have substantially revised the masses upward (Luhman et al., 2007).

Despite these previous discoveries, the existence of such abundant ($\sim$$9\%$) low-mass binaries with wide separations has not been anticipated.

Explaining the formation of binary planets with wide separations is a challenging task. Although two planets (on mutually close astrocentric orbits) can in principle be ejected as a binary pair through stellar encounters (Wang et al., 2024), recent numerical studies suggest that the production rate through this direct ejection is at least two orders of magnitude lower than what is needed to explain the observed number and fraction in the Trapezium Cluster (Yu & Lai, 2024; Portegies Zwart & Hochart, 2023). Another proposed mechanism is through tidal capture during repeated planet-planet scattering in planetary systems (Ochiai et al., 2014; Lazzoni et al., 2023). However, tides generally produce tight binaries rather than wide binaries.

Given the difficulties to form wide binary planets out of the protoplanetary disk (i.e., planet-like formation), star-like formation, in which binary planets form directly within the cluster, may be preferred (see also Portegies Zwart & Hochart, 2023), although the exact formation mechanism remains unclear and requires further investigation.

Once formed, the binary planets in the cluster may evolve as a result of dynamical interactions with passing stars. This is especially the case for such wide binaries as the JWST JuMBOs in the Trapezium Cluster. Due to the extremely high star number density ($n_{\star}\sim 1\times 10^{4}\ \text{pc}^{-3}$), the dynamical lifetime of a typical JWST JuMBO is comparable to the age of the cluster (see Section 2; see also Portegies Zwart & Hochart 2023). Therefore, it is probable that the wide binary planets seen by JWST were born with closer separations and then got softened by stellar flybys. These binary planets with initially close separations may then be the lower mass counterparts of binaries with L-, T-, and Y-type primaries, for which the binary fraction is 10–30% and the binary separation peaks at 3–10 au (Burgasser, 2007; Joergens, 2008).

Regardless of the true nature of the recently discovered PMOs and JuMBOs by JWST, it is necessary to study the dynamical evolution of the planetary binaries in a dense cluster environment to explore their origins.

2 Analytical Method

The dynamical evolution of binary stars in star clusters has been studied extensively in the literature (e.g., Heggie 1975, Heggie & Hut 2003, and Binney & Tremaine 2008). The problem of the binary planet–star interaction is a simplified version, as the binary mass and binding energy are too small to influence the star’s motion in a meaningful way. Therefore, we aim to build a simple analytical model that focuses on the $m\ll M_{\star}$ interaction, where $m$ is the total mass of the binary and $M_{\star}$ is the much more massive third body. For simplicity, only equal-mass binaries are considered in this work.

2.1 Effect of a Single Flyby

Figure 1: Hyperbolic trajectory of the binary planet relative to the passing star. Two planets (blue circles) with a total mass $m$ orbit their mutual barycenter (black cross) with the relative semimajor axis $a$. The periastron distance (closest approach) of the barycenter is denoted as $q_{h}$, which is related to the hyperbolic semimajor axis $a_{h}$ and eccentricity $e_{h}$. The geometry of the flyby trajectory is determined by the impact parameter $b$, the relative velocity between the binary barycenter and the star at infinity $v_{\infty}$, and the stellar mass $M_{\star}$.

As illustrated in Figure 1, the hyperbolic trajectory of the binary barycenter relative to the passing star can be described with the following orbital elements:

$$a_{h} =\frac{\mathcal{G}M_{\star}}{v_{\infty}^{2}},$$ (1)

$$q_{h} =a_{h}(e_{h}-1)=a_{h}\left(\sqrt{1+j_{h}^{2}}-1\right),$$

$$j_{h} =\frac{b}{a_{h}},$$

where $\mathcal{G}$ is the gravitational constant. The hyperbolic semimajor axis $a_{h}$ and the periastron distance $q_{h}$ are determined by the impact parameter $b$, the relative velocity between the star and the binary barycenter at infinity $v_{\infty}$, and the mass of the passing star $M_{\star}$. Assuming $M_{\star}=1M_{\odot}$, the hyperbolic semimajor axis can be quickly estimated with $a_{h}=(30~{}\text{km}~{}\text{s}^{-1}/v_{\infty})^{2}~{}\text{au}$. For young massive clusters with a typical velocity dispersion of $\sigma\sim 1~{}\text{km}~{}\text{s}^{-1}$, $a_{h}$ is of the order of $\sim$1000 au for solar-mass stars. The ratio between $b$ and $a_{h}$ is defined as $j_{h}$, representing the dimensionless angular momentum on a given semimajor axis (see Tremaine, 2023).

The pericenter velocity $v_{q}$ can be obtained via the conservation of angular momentum

$$v_{q}=\frac{bv_{\infty}}{q_{h}}.$$ (2)

The typical timescale of the encounter is

$$t_{\text{enc}}\equiv\frac{2q_{h}}{v_{q}},$$ (3)

and a close encounter is defined as $t_{\text{enc}}\lesssim P$, with the orbital period of the binary $P$ given by

$$P=2\pi\sqrt{a^{3}/(\mathcal{G}m)}.$$ (4)

Here $a$ is the binary semimajor axis and $m$ is the binary total mass. During such encounters, the relative location of the two planets does not change over the span when the gravitational pull from the star is most intense. As a result, the change in the specific orbital energy of the binary, $E=-\mathcal{G}m/(2a)$, can be estimated with an impulse approximation (Farinella & Chauvineau, 1993):

$$\Delta E\simeq\frac{\mathcal{G}^{2}M_{\star}mP}{4abv_{\infty}q_{h}}\int_{-f_{l}}^{f_{l}}\Psi\mathrm{d}f.$$ (5)

Here, the integral term evaluates the total energy change along the hyperbolic path, which is related to the spatial orientation of the hyperbola relative to the binary orbital plane (ascending node $\Omega$ and inclination $i$) as well as the exact locations of the two components. Monte-Carlo simulations of random flyby angles and initial locations show that the orientation integral has an average of zero and a root mean square of the order of unity (Farinella & Chauvineau, 1993). This implies that random close encounters have equal probabilities of softening and hardening the binary. Therefore, the effect of multiple flybys can be modeled as a random walk in $E$ space, with the typical energy change $|\Delta E|$ (i.e., the variance of the $\Delta E$ distribution) being

$$\left|\Delta E\right| \simeq\frac{\mathcal{G}^{2}M_{\star}mP}{4abv_{\infty}q_{h}}.$$ (6)

We define critical flybys as stellar encounters with a typical energy change larger than the binary binding energy $\left|\Delta E\right|>\left|E\right|$. Such encounters would likely induce significant changes to the semimajor axis $a$ and orbital eccentricity $e$ of the binary, or even ionize it. Rearranging the expression, one obtains an inequality for two timescales:

$$t_{\text{cr}}\equiv\frac{2q_{h}}{v_{\infty}}j_{h}<P,$$ (7)

where the left-hand side is defined as the critical flyby timescale $t_{\text{cr}}$.

The encounter timescale $t_{\text{enc}}$ (Equation 3) can be rewritten as $t_{\text{enc}}=(2q_{h}/v_{\infty})(q_{h}/b)$ using Equation (2). The relation between the two timescales, $t_{\text{enc}}$ and $t_{\text{cr}}$, is thus given by

$$t_{\text{cr}}=t_{\text{enc}}\left(\frac{j_{h}^{2}}{\sqrt{1+j_{h}^{2}}-1}\right),$$ (8)

and it is trivial to verify that $t_{\text{cr}}>2t_{\text{enc}}$ holds for all $j_{h}$. As a result, $t_{\text{cr}}<P$ is always a stronger constraint than $t_{\text{enc}}<P$; in other words, all critical flybys are close encounters where the impulse approximation is applicable.

Figure 2: Critical timescale $t_{\text{cr}}$ (Equation 7) computed on the $b$-$v_{\infty}$ space assuming a solar-mass star. Dashed lines denote the orbital periods of 10 $M_{J}$ JuMBOs with a wide range of semimajor axes $a$, of which critical flybys lie on the left-hand side. The solid line represents $j_{h}=1$ ($e_{h}=\sqrt{2}$), which divides the parameter space into the hyperbolic regime ($j_{h}>1$, above the line) and the near-parabolic regime ($j_{h}<1$, below the line).

We compute the critical timescale for a series of $b$ and $v_{\infty}$ using Equation (7) and plot the results in Figure 2. The critical flyby timescale spans several orders of magnitude, from a few days to more than 100 Myr. $t_{\text{cr}}$ also scales differently given the geometry of the flyby trajectory: Near-parabolic flybys with $e_{h}\sim 1$ ($j_{h}\sim 0$) affect the binary more effectively than hyperbolic flybys ($j_{h}>1$) do, as indicated by a slope change separated by the solid line $j_{h}=1$.

To derive the scaling law, one first expands the pericenter distance

$$q_{h} \approx\begin{cases}\cfrac{b^{2}}{2a_{h}}\quad&\text{ if }j_{h}<1,\\ b\quad&\text{ if }j_{h}>1,\\ \end{cases}$$ (9)

and then explicitly writes the fractional orbital energy change for different $j_{h}$

$$|\Delta\mathcal{E}|\equiv\frac{|\Delta E|}{|E|}=\begin{cases}\left(\mathcal{G}M_{\star}\right)^{2}Pb^{-3}v_{\infty}^{-3}\quad&\text{ if }j_{h}<1,\\ \left(\mathcal{G}M_{\star}P/2\right)b^{-2}v_{\infty}^{-1}\quad&\text{ if }j_{h}>1.\end{cases}$$ (10)

Equating $|\Delta\mathcal{E}|$ with 1 yields

$$b_{\text{cr}}=\begin{cases}\left[\left(\mathcal{G}M_{\star}\right)^{2}P\right]^{\frac{1}{3}}v_{\infty}^{-1}\quad&\text{ if }j_{h}<1,\\ \left[\mathcal{G}M_{\star}P/2\right]^{\frac{1}{2}}v_{\infty}^{-\frac{1}{2}}\quad&\text{ if }j_{h}>1,\end{cases}$$ (11)

where $b_{\text{cr}}$ is the impact parameter for critical flybys. Notice $b_{\text{cr}}\propto v_{\infty}^{-1}$ in the near-parabolic regime and $b_{\text{cr}}\propto v_{\infty}^{-\frac{1}{2}}$ in the hyperbolic regime, which corresponds to the slope change in Figure 2.

One can easily estimate the orbital stability of a binary planet against a particular flyby by comparing $t_{\text{cr}}$ with $P$. In Figure 2, we plot orbital periods of 10 $M_{J}$ JuMBOs with $a$ spanning from 0.1 au to 1000 au in dashed lines, and the areas to the left of these lines indicate flyby parameters that would likely induce significant changes in binding energy. For JWST JuMBOs with $a=25$–400 au, critical flybys occur at $b_{\text{cr}}\lesssim 600$–6,000 au for $v_{\infty}=1$ km $\text{s}^{-1}$, and $b_{\text{cr}}\lesssim 300$–3,000 au for $v_{\infty}=2.5$ km $\text{s}^{-1}$. It is worth pointing out that such close flybys may be happening in some of the observed systems (see JuMBO33 and 34 in Figure 3 of Pearson & McCaughrean 2023)

2.2 Dynamical Lifetime Inside a Cluster

Until now, we have only looked at the average effects caused by a single passing star. However, in a cluster environment where thousands of stars are present, it is essential to consider the ongoing interactions from multiple stars, each with varying trajectories and characteristics.

The frequency of encounters with the parameters $b$ and $v_{\infty}$ for a period of time $t$ is given by

$$N=n_{\star}t(\pi b^{2})v_{\infty},$$ (12)

where $n_{\star}$ is the number density of the star cluster, and $\pi b^{2}$ represents the cross section of the encounter. For critical flybys with $b<b_{\text{cr}}$, the induced fractional energy change $\Delta\mathcal{E}$ is of the order of unity. In contrast, the energy drift of non-critical distant flybys can build up, whose cumulative change after $N$ flybys follows the random walk equation $|\Delta\mathcal{E}|_{N}=\sqrt{N}|\Delta\mathcal{E}|$. Apparently, $|\Delta\mathcal{E}|_{N}\propto b^{-2}$ in the near-parabolic case and $\propto b^{-1}$ in the hyperbolic case. In other words, in both cases, the total energy change is dominated by the single closest stellar encounter, rather than the cumulative effects of more frequent distant encounters. As a result, the ionization of a planetary-mass binary mostly results from the closest encounter, and the dynamical lifetime can be defined as the expected time interval between two critical flybys (Equation 13).

$$T \equiv\frac{1}{n_{\star}(\pi b_{\text{cr}}^{2})v_{\infty}}$$ (13)

$$\simeq\begin{cases}0.9~{}\text{Myr}~{}\left(\cfrac{n_{\star}}{10^{4}~{}\text{pc}^{-3}}\right)^{-1}\left(\cfrac{M_{\star}}{M_{\odot}}\right)^{-\frac{4}{3}}\left(\cfrac{m}{10~{}M_{J}}\right)^{\frac{1}{3}}\left(\cfrac{a}{100~{}\text{au}}\right)^{-1}\left(\cfrac{v_{\infty}}{1~{}\text{km s}^{-1}}\right)\quad&\text{ if }j_{h}<1,\\ 1.4~{}\text{Myr}~{}\left(\cfrac{n_{\star}}{10^{4}~{}\text{pc}^{-3}}\right)^{-1}\left(\cfrac{M_{\star}}{M_{\odot}}\right)^{-1}\left(\cfrac{m}{10~{}M_{J}}\right)^{\frac{1}{2}}\left(\cfrac{a}{100~{}\text{au}}\right)^{-\frac{3}{2}}\quad&\text{ if }j_{h}>1.\end{cases}$$

Notice that the discontinuity between the two equations at $j_{h}=1$ arises from different approximations of $q_{h}$ in Equation (9). Adopting the typical number density of $1\times 10^{4}~{}\text{pc}^{-3}$ and the velocity dispersion of $1~{}\text{km s}^{-1}$ for young star clusters, Equation (13) suggests that the dynamical lifetime of a JuMBO with $a=100$ au is only of the order of $\sim$1 Myr. It should be noted that the structure of star clusters dynamically evolves, so the JuMBO stability may vary depending on the cluster density over time.

We also calculate the dynamical lifetime of potential JuMBOs within other types of cluster, assuming a constant star density and relative velocity. For open clusters such as Hyades ($n_{\star}=1~{}\text{pc}^{-3}$ and $v_{\infty}=0.25~{}\text{km s}^{-1}$, Perryman et al. 1998), binary planets similar to the JWST JuMBOs are stable for several billion years, significantly longer than the $\sim$700 Myr age of the Hyades; Conversely, the same type of binary planets would have long been destabilized if they formed in the globular cluster 47 Tucanae ($n_{\star}=1\times 10^{5}~{}\text{pc}^{-3}$ and $v_{\infty}=10~{}\text{km s}^{-1}$, Meylan & Mayor 1986), given their short dynamical lifetime ($\lesssim$1 Myr) compared to the cluster age of 11 Gyr.

2.3 Numerical Validation

We carry out numerical simulations in which binary planets encounter a solar-mass star from random directions. For each simulation, one thousand equal-mass binaries with a binary mass of $10M_{J}$ are initialized on a circular orbit with $a_{0}=15$ au. Their barycenter is then set on a hyperbolic trajectory relative to the passing star with $v_{\infty}=1~{}\text{km}~{}\text{s}^{-1}$. Both the initial binary mean anomalies and the longitude of the ascending nodes are distributed uniformly in $(0,2\pi)$, while the relative inclinations between the hyperbolic trajectory and the binary orbital plane are distributed uniformly in $\cos(i)$, corresponding to the isotropic angles of entry. We then obtain the corresponding critical impact parameter $b_{\text{cr}}=460$ au using Equation (11), and set $b=b_{\text{cr}}/1.5,b_{\text{cr}},$ and $1.5b_{\text{cr}}$, respectively, for the three experiments. All three cases have $j_{h}<1$ and thus belong to the near-parabolic regime. Simulations are performed using the ias15 integrator of rebound (Rein & Spiegel, 2014).

Figure 3: Upper: Orbital distributions and histograms of 1,000 equal-mass binary planets after encountering a solar-mass star with $v_{\infty}=1~{}\text{km}~{}\text{s}^{-1}$. Binaries are initialized with $a_{0}=15$ au, $e_{0}=0$, and $m=10~{}M_{J}$, while the chosen impact parameters correspond to a deep critical flyby ($b_{\text{cr}}/1.5$, left panel), a critical flyby ($b_{\text{cr}}$, middle panel), and a distant flyby ($1.5b_{\text{cr}}$, right panel), respectively. The $x$ axis indicates the fractional energy change $\Delta\mathcal{E}$. The left $y$ axis represents the binary eccentricity post flyby, whereas the right $y$ axis denotes the number of binaries in each bin. The orange dashed line marks the energy change of 0, with a negative $\Delta\mathcal{E}$ corresponding to a decrease in $a$, and vice versa (red arrows). The green shaded region denotes the observable JuMBOs with $a>25$ au. Lower: Histograms of semimajor axes of captured planets ($a_{\text{cap}}$) from ionized JuMBO pairs.

The post-flyby orbital distributions of the three simulations are shown in Figure 3, with $\Delta\mathcal{E}$ denoting the fractional energy change and $e$ denoting the eccentricity. We decide to plot $\Delta\mathcal{E}$ instead of $\Delta a$ because it is a dimensionless scale-free quantity. In other words, even if we ran the same simulations (with corresponding $b_{\text{cr}}$) for binaries with different $m$ and $a_{0}$, their orbital distributions in the $\Delta\mathcal{E}$ space would still look similar to those in Figure 3.

We define the following branching fractions:

1. $f_{<\text{25}}$ — the number of binaries with $a<25$ au divided by the total number of primordial binaries,

2. $f_{>\text{25}}$ — the number of binaries with $a>25$ au divided by the total number, and

3. $f_{\text{ion}}$ — the number of ionized binaries (signified by a positive binary energy post flyby) divided by the total number (i.e., ionization fraction).

With the above definitions, we have $f_{<\text{25}}+f_{>\text{25}}+f_{\text{ion}}=1$. The values of these three fractions are shown in Figure 3.

The left and middle panels of Figure 3 demonstrate the aftermath of critical flybys with different impact parameters. Although both flybys show roughly equal chances of softening and hardening the orbit (that is, a near symmetric distribution centered at $0$), the middle panel shows a strong peak at $\Delta\mathcal{E}=0$, implying that a significant portion of binaries are not strongly affected by the flyby. In comparison, 64% of the binaries are ionized by deeper critical flybys.

When the equal-mass binaries are ionized by the passing star, a common outcome is the capture of one planet by the passing star and the ejection of the other planet. The captured planets are on highly eccentric orbits around the star, with $a_{\text{cap}}$ peaked at 300–400 au and the minimal eccentricity of $e_{\text{min}}=0.5$–0.6 (lower panels in Figure 3). These distant and elliptic planetary orbits are mostly not stable inside a dense cluster environment, and one can expect that subsequent stellar perturbations would easily disrupt them and produce FFPs (Spurzem et al., 2009). It is worth noting that three-body interactions occurred frequently in the early Solar System, both among planetesimals (Funato et al., 2004) and between planetesimals and a planet. For example, Triton – the largest moon of Neptune – is thought to be captured during a close encounter between Neptune and a transneptunian binary (Agnor & Hamilton, 2006).

There is an apparent dichotomy between critical flybys and noncritical flybys. As illustrated in the right panel of Figure 3, the noncritical ($b=1.5b_{\text{cr}}$) flybys induce little variation in both the binary energy and eccentricity. In a cluster environment, such flybys would be roughly twice as likely as critical ($b_{\text{cr}}$) flybys (Equation 12). However, the cumulative effect of two $1.5b_{\text{cr}}$ flybys is by no means comparable to a $b_{\text{cr}}$ flyby. This again validates our previous analysis that the orbital change of a binary planet is dominated by the single closest passing star.

3 Application to JWST JuMBOs

We have demonstrated that starting from a single population, gravitational interactions with passing stars are able to produce a variety of binary systems from a few to several hundred au. The same dynamical process also produces planetary binaries on tighter orbits (close binaries), FFPs, and temporarily captured planets, which are expected to be liberated as FFPs shortly after due to their extremely large semimajor axes.

Limited by the observational capability of the telescope, however, only wide binaries with separations above a certain threshold ($25$ au, for JuMBOs in the Trapezium) can be resolved by JWST, and it is unclear whether the ‘non-binary’ PMOs are isolated FFPs, or unresolved close JuMBOs, or more likely, a combination of both. The applicable constraint from Pearson & McCaughrean (2023) is the $\sim$9% JuMBO-PMO fraction (the number of $a>25$ au JuMBOs divided by the total number of all planetary-mass objects), i.e., the wide binary fraction $f_{\text{wide}}$, which is also given by

$$f_{\text{wide}} =\frac{f_{>25}}{2\times f_{\text{ion}}+f_{<25}+f_{>25}+R_{S:B}}$$ (14)

$$=\frac{f_{>25}}{1+f_{\text{ion}}+R_{S:B}},$$

where $R_{S:B}=N_{\text{single}}/N_{\text{binary}}$ is the primordial ratio between single planets and binary planets governed by the formation process. Since the formation mechanism is still unknown, we simply assume that it only produces binary objects, that is, $R_{S:B}=0$. As a result, Equation 14 gives the upper bound of $f_{\text{wide}}$.

The ionization fractions for the two critical flybys presented in Figure 3 are 64% and 27%, respectively. Deeper critical flybys are also more efficient in making $a>25$ au binaries (6.3% compared to 2.0%), because of stronger energy kicks. The resultant $f_{\text{wide}}$ therefore ranges from 2% to 4%. It is worth pointing out that $f_{\text{wide}}$ is very sensitive to the initial $a_{0}$ and the definition of wide JuMBOs. Therefore, we next explore numerically what initial $a_{0}$ best reproduces the $\sim$9% wide binary fraction, assuming a “primordial JuMBO population” formed at the beginning of the star cluster.

Additional three-body simulations similar to those in Figure 3 were performed. Instead of fixing $b$, we initialize hyperbolic trajectories based on a realistic distribution ($\mathrm{d}N/\mathrm{d}b\propto b$), with the maximum $b$ being what gives $N=1$ in Equation (12). To explore different cluster environments, we define a single parameter $\chi=\int_{0}^{t}n_{\star}\mathrm{d}t$ as the stellar number density-weighted residence time (Batygin et al., 2020), which indicates the time-accumulated effect taking into account the evolution of the cluster. In the lowest order, $\chi\approx n_{\star}t$ assuming a fixed star density after cluster formation.

Figure 4: Ionization fraction $f_{\text{ion}}$ ($+$) and wide binary fraction $f_{\text{wide}}$ ($\times$) from stellar flyby simulations as a function of the binary period $P$. Each symbol represents the result of 5,000 JuMBOs with the same mass (from right to left, $m=1,2,5$, and 10 $M_{J}$, respectively) and $a_{0}$ (grouped by the same color) interacting with a solar-mass star. The impact parameter $b$ is generated from a realistic distribution assuming a density-weighted residence time of $\chi=2\times 10^{4}$ Myr pc^-3 and a fixed relative velocity $v_{\infty}$ of $1~{}\text{km}~{}\text{s}^{-1}$. $N_{\text{cr}}$ on the top of the $x$ axis marks orbital periods where the frequency of critical flybys is $>$$0.1$ and $>$$1$, respectively. The range of $f_{\text{wide}}$ potentially consistent with the JWST discovery (1%–10%) is highlighted in yellow.

In Figure 4, pluses and crosses represent $f_{\text{ion}}$ and $f_{\text{wide}}$ from a simulation set in which 5,000 binary planets encounter random passing stars. The $x$ axis denotes the binary period $P$, a proxy for its dynamical stability (see Equation 7). The chosen $\chi=2\times 10^{4}$ Myr pc^-3 environment correspond to the lower bound of the Trapezium Cluster core density (Heggie & Aarseth, 1992; Hillenbrand & Hartmann, 1998). As shown in both panels, the ionization fraction $f_{\text{ion}}$ grows monotonically as the binary period, in agreement with analytical theory.

For primordially tight JuMBOs ($a_{0}<25$ au), orbital expansion is the main cause for generating wide binaries ranging from 25 to 400 au. The post-evolution $f_{\text{wide}}$ lies between 1% to 10% (yellow highlight) for $a_{0}=10,15$, and 20 au, consistent with the $\sim$9% observed by JWST in terms of order of magnitude. $f_{\text{wide}}$ is merely $\sim$0.1% for $a_{0}=5$ au JuMBOs, due to stronger orbital stability (signified by their low ionization fractions) and fewer chances of getting wide orbits. These two effects can be quantitatively estimated in the following manner: 1) To increase a binary $a$ from 5 au to $>$25 au, a fractional energy change of $\Delta\mathcal{E}=(0.8,1)$ is required. Compared to $\Delta\mathcal{E}=(0.4,1)$ needed for raising 15 au to $>$25 au, it is an event at least a factor of three (0.6/0.2) less likely to occur with random energy kicks. 2) Critical flybys are roughly $5\times$ more frequent for $a_{0}=15$ au than for $5$ au (obtained by combining Equations 11 and 12). Therefore, we expect a factor of $\sim$15 difference in $f_{\text{wide}}$ for $a_{0}=15$ au and $5$ au, in excellent agreement with the left panel of Figure 4. There is also a factor of $\sim$2 difference in $f_{\text{wide}}$ for JuMBOs with different masses: more massive JuMBOs tend to have a smaller wide binary fraction, in line with JWST observation (Figure 4 in Pearson & McCaughrean, 2023). In contrast, Equation (13) shows that if JuMBOs primordially formed as wide binaries ($a_{0}>25$ au), the more massive JuMBOs (which are more stable against stellar flybys) would tend to have higher $f_{\text{wide}}$, contradicting the JWST observation.

Therefore, we conclude that if a primordial JuMBO population formed within the Trapezium Cluster, the scenario in which they were initially on tight orbits ($a_{0}=10$–20 au) reproduces the observed ratio. This scenario requires a high density-weighted residence time of $\chi\gtrsim 10^{4}$ Myr pc^-3, suggesting primordial JuMBOs must have spent most of the time near the core of the cluster until critical flybys significantly changed their orbital properties. This postulated population may also connect to tighter JuMBOs ($a<5$ au, see Beichman et al. 2013 and Best et al. 2017), and our results indicate that the current PMO population observed by JWST is partially made up of unresolved binaries. Follow-up observations that characterize the properties of these PMOs are the key to unveiling their origins.

4 Discussion

In this work, we have built analytical tools to study the orbital evolution and dynamical stability of binary planets within star clusters. By conducting three-body simulations of random encounters between binary planets and stars, we have shown that critical stellar flybys are efficient at producing binary planets with a wide range of $a$ and $e$ from a dominant population.

Applying our results to recently discovered JWST JuMBOs (McCaughrean & Pearson, 2023), we show that their typical dynamical lifetime is comparable to the age of the Trapezium cluster, suggesting that these wide binaries may have encountered stellar flybys and thus dynamically evolved. We propose that the JWST JuMBOs were formed as tighter binaries near/within the core of the Trapezium cluster. To best reproduce the observed $\sim$9% wide binary fraction, an initial semimajor axis of $a_{0}=10$–20 au and a density-weighted residence time of $\chi\gtrsim 10^{4}$ Myr pc^-3 are favored.

The discovered JWST JuMBOs are currently within a radius of $\sim$0.6 pc, larger than the core of the Trapezium cluster, which has a radius of $\sim$0.1–0.2 pc and a central density of $\sim$2–5$\times 10^{4}$ pc^-3 (McCaughrean & Stauffer, 1994; Hillenbrand & Hartmann, 1998).

Given that the dynamical timescale of the ONC (as well as the inner Trapezium Cluster) is comparable to its age, it is likely that the dynamical equilibrium of the cluster was either recently reached or not yet established (Portegies Zwart et al., 2010). Furthermore, three-body simulations of young massive clusters suggest that the stars were more clustered at an early stage, leading to a higher stellar density (see Figure 1 of Fujii, 2015). Taking into account the dynamical evolution of the cluster itself, it is thus reasonable to postulate that the discovered JWST JuMBOs have experienced $\chi\gtrsim 10^{4}$ Myr pc^-3 in the past $\sim$1 Myr.

The current JWST JuMBO distribution beyond the cluster core could also be produced by binary–star encounters: JuMBOs that had critical stellar encounters in the past are also those that had experienced a significant barycentric velocity kick. Specifically, the 90^∘ deflection impact parameter $b_{\text{90}}$ has exactly the same definition as $a_{h}$ (Binney & Tremaine, 2008), suggesting that any stellar encounters in the near-parabolic regime ($b<a_{h}$) would also deflect the relative velocity vector by $>$90^∘. This would induce a barycentric velocity change comparable to its original speed. Therefore, even in a scenario where JuMBOs are hypothesized to have formed closer to the core, it is not so surprising to discover escaped members within 1 pc around the core at 1 Myr, given the typical barycentric velocity of $v=1~{}\text{km s}^{-1}\approx 1~{}\text{pc Myr}^{-1}$. Future optical surveys on the outskirts of the Trapezium Cluster should help test this validity of this scenario.

During the encounters between binary planets and the flyby stars, one of the planets in the binary may get captured by the star into a wide-separation and eccentric orbit. Our simulations show that the production rate of such captured planets is not negligible, especially in the dense environment. This can potentially provide an alternative explanation for the existence of Jupiter-mass planets on extremely wide and eccentric orbits that have been found by direct imaging surveys (e.g., Kraus et al., 2014; Bryan et al., 2016). It may also be connected to the temporarily present rogue planet (Gladman & Chan, 2006; Huang et al., 2022), or the hypothesized “Planet Nine” (Batygin et al., 2019), in the solar system.