Innocent Bystanders: Orbital Dynamics of Exomoons During Planet–Planet Scattering

The planets are simulated as oblate by the addition of the quadrupole moment (J₂) term to the gravitational force between point mass bodies. This treatment ensures that the orbits of moons within the critical orbital distance from a planet are simulated correctly and that unexpected instability will not occur while the planet keeps the close-in moons' orbital angular momentum precessing rapidly (Hong et al. 2015). Numerically, this is done by adding a customized force term that accounts for the quadrupole moment contribution in the gravitational attraction of a planet to other planets and satellites in Mercury. Since the spin of the planet is allowed to change, the direction of the force from the planet's bulge follows the spin evolution of the planet. However, the code assumes that the spin of gas giant planets can react to external torques instantaneously during a planetary close encounter, whereas in reality the timescale for a fluid body to react to torques might be longer than the close encounter time.

The spin of an oblate planet evolves under the influence of the central body and other planetary bodies, especially during a close encounter between planets. In the planet–planet scattering scenario, the spin of a planet can often evolve dramatically, in which the stability of its moons needs to be tested (Tremaine et al. 2009). In this work, the spin evolution model for the planet uses an instantaneous equation of motion, in order to correctly simulate the effect of planetary close encounters on the spin of a planet. The equation of motion for the spin evolution is derived in the following paragraphs.

An extended mass receives different amounts of gravitational attraction from a distant point mass at different radial distances within its interior, which results in a torque on the extended mass, causing its spin to evolve. The force from a point mass M on an extended mass m is

$$\begin{eqnarray}{\boldsymbol{F}} & = & \displaystyle \frac{-\,G\,M\,m}{{r}^{2}}\left\{\hat{r}-3\,{J}_{2}{\left(\displaystyle \frac{R}{r}\right)}^{2}\,[(5{(\hat{r}\cdot \hat{s})}^{2}-1)r\right.\\ & & -2(\hat{r}\cdot \hat{s})\hat{s}]+\cdots \},\end{eqnarray} \tag{ 1 }$$

where G is the gravitational constant, r is the distance between the centers of m and M and $\hat{r}$ points from m to M, R is the radius of m, and $\hat{s}$ the normalized spin vector of m (Hilton 1991). The torque per moment of inertia is

$$\begin{eqnarray}&&\displaystyle \frac{{\boldsymbol{r}}\times {\boldsymbol{F}}}{I}=\displaystyle \frac{d}{{dt}}({\rm{\Omega }}\,\hat{s}),\end{eqnarray} \tag{ 2 }$$

where Ω is the rotation period of the planet, which for simplicity is always assumed to be constant and equal to 10 hr, the same as for Jupiter. The time evolution of the spin components can thus be described in the following equation:

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{s}}{{dt}}=\displaystyle \frac{-3\,G\,M\,m}{{r}^{3}}\displaystyle \frac{{J}_{2}}{\lambda \,{\rm{\Omega }}}(\hat{r}\cdot \hat{s})(\hat{r}\times \hat{s}),\end{eqnarray} \tag{ 3 }$$

where λ is the normalized moment of inertia. λ is taken to be 0.25, close to that of Jupiter, in all simulations. The above equation is derived assuming an instantaneous force and torque within which no variables are averaged over time. It is therefore suitable for the planet–planet scattering case, where planets sometimes interact extremely closely (i.e., ∼0.001 au) over a brief period of time (i.e., ∼1 year). As a side note, equations of spin evolution for secular problems, such as the obliquity evolution of Mars, require the orbit to stay unchanged at least over the timescale of the body's orbital period:

$$\begin{eqnarray}&&\displaystyle \frac{d\hat{s}}{{dt}}=\left(\displaystyle \frac{-3}{2}\displaystyle \frac{{n}^{2}}{{\rm{\Omega }}}{J}_{2}\lambda \right)(\hat{l}\cdot \hat{s})(\hat{l}\times \hat{s})\end{eqnarray} \tag{ 4 }$$

(Ward 1973, 1974; Bills 1990, 2005). The above equation uses averaged quantities such as $\hat{l}$, which is the angular momentum averaged over an orbital period. During a close encounter, in particular, such an assumption lacks accuracy because the small timescale of the variation of perturbation is sometimes comparable to the orbital period of the binary planet pair.

Equation (3), with its dependence on the body's position relative to other bodies, is plugged into the Bulirsch–Stoer integrator alongside the integration for positions and velocities and is updated once per time step. The validity of the integrator is tested with the past obliquity evolution of Mars. The initial conditions are taken from Quinn et al. (1991) at epoch JD 2433280.5, with the obliquity of Mars = 25°. The simulation is integrated backwards and produces results with a close resemblance to those of Ward (1973).

Unless otherwise specified, the initial settings of the simulations are configured as follows: (1) the integration error limit in the orbital energy and angular momentum is ${10}^{-12}$ and (2) the integration time step is 0.1 days in order to accurately integrate satellites with orbital periods as short as a few days like the Galilean satellites.

In each simulation of simulation set 1, the system contains a Sun-like star, a Jupiter-mass planet that hosts moons, and a perturbing planet. The host planet is at 5 au on a circular orbit. Moons are massless test particles placed between 0.01 and 0.35 Hill radii (thereafter R_H) from the host planet on circular and co-planar orbits. Four moons are spaced out evenly in angular position at each of the 10 planetocentric semimajor axes. The planets are set up to undergo close encounters (within each other's Hill sphere) by the following. The perturbing planet is initially placed at 1.2 Hill radii from the host planet. The parameters of the perturbing planet are the only variants across different simulations in set 1. The perturbing planet's mass is equal to 0.1, 0.5, 1.0, or 2.0 M_J (Jupiter mass). Its starting velocity relative to the host planet (V_rel) ranges from 0.001 to 0.009 au day⁻¹, and the direction of its velocity vector is restricted to lie within 60° from the host planet's velocity vector. The range of V_rel is determined from all close encounters in a subset of planet–planet scattering simulations for 10⁶ years. The perturber's starting impact parameter (b) is 0.005–0.1 au (host planet radius = 0.00047 au). The impact parameter is set as a three-dimensional sphere of vectors surrounding the center of the host planet. All simulations last for 10 years. Simulations in which planets collide with each other or in which the perturber starts with e > 1 are excluded. In most of the close encounters during the instability period, the latter configurations are rare.

The minimum close encounter distance d_min plays a major role in determining moon stability. From the initial conditions, it is mainly b that determines d_min and thus moon stability. V_rel plays a minor role because the escape velocity of the host planet (0.035 au day⁻¹) is 1–2 orders of magnitude (∼35–350 times) greater than V_rel. Figure 1 shows the stability boundary of moons around host planets that experience close encounters with different closest approach distances. The stability boundary is set where at least two moons on the same planetocentric distance survive with e ≤ 0.5, because they are anticipated to have high probabilities to survive a full simulation. As will be discussed in Section 4, 80% of the surviving moons have eccentricities under 0.5. The stability boundary for each d_min is an average taken from all of the host planets with the same d_min and with a perturber having the same mass. The larger d_min is, the larger the stability boundary. For simulations with perturber mass 0.5–2 M_J, the stability limits in Figure 1 are positively correlated with d_min in the curves, which demonstrate that d_min plays the most important role. The slope of the curves is determined by the perturber mass, demonstrating that the perturber mass plays the second most important role in determining the stability boundary of a planet. The slope is roughly linear when the perturber mass is less than twice the mass of the host planet, but when the perturber mass reaches twice the mass of the host planet, the linear relation starts to fail.

Zoom In Zoom Out Reset image size

Host planets with lower-mass perturbers have larger stability boundaries. Simulations with a 0.5 M_J perturber have stability boundaries close to or greater than d_min, and simulations with a 0.1 M_J perturber have much larger stability boundaries than d_min. For simulations with larger perturbers (${M}_{{\rm{J}}}$ and 2 ${M}_{{\rm{J}}}$), the stability boundary is 0.6–0.8 times d_min for ${d}_{\min }\lt 0.05\,\mathrm{au}$. Except for simulations with a 2 M_J perturber, there exists a boundary in d_min beyond which almost all moons from $0.01\,\mathrm{and}\,0.35\,{R}_{{\rm{H}}}(\sim 0.003\mbox{--}0.12\,\mathrm{au})$ are stable. For simulations with a 0.1M_J/0.5 M_J/1 M_J perturber, when d_min is greater than $\sim 0.05/0.1/0.1\,\mathrm{au}$, the stability boundary lies close to 0.35 R_H.

Regarding the moon capture rate, like the moon survival rate, the most important determining factor is d_min, and the perturber mass is the second. Figure 2 shows the average number of moons captured by the perturber from the host planet in simulations with different d_min. A perturber that has come closer to the host planet is able to pass by more moons, and because the capture of moons by the perturber requires the perturber to be close to the moons (Figure 3(b)), a smaller d_min yields a higher capture rate. For planet masses ($0.5\mbox{--}2\,{M}_{{\rm{J}}}$) comparable to the host planet, the capture rate is roughly inversely proportional to d_min. A larger perturber mass slightly enhances the capture rate. As seen from Figure 2, the moon capture rate is only significant (>0.1) for d_min under ∼0.03 au for simulations with 0.5–2 M_J perturbers. 0.1 M_J perturbers have nearly no capability to capture moons from a 1 M_J host planet during their encounter, as seen from the average number of captures, but in a few cases the capture rate can reach 1%–5%.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 3 shows how the close encounter geometry can affect the dynamical outcome of moons. The scenario contains an M_J host planet and an M_J perturber. d_min for this encounter event is 0.01 au, shorter than the initial semimajor axis of the four moons (0.014 au/0.042 R_H) with their angular positions evenly spaced out. Despite having the same initial orbital radius, the moons had a range of dynamical outcomes. In (a), the moon is dragged in the forward direction and changes its orbit but stays stable. In (b), the moon heads toward the planet and is able to reach it with the right timing and configuration from the back of the planet for a smooth transfer onto a stable orbit around the perturber to happen. In (c) and (d), the moons are dragged in the backward direction by the planet and thus open up their orbits to become unbound. In summary, moons on the same orbits can end up with different dynamical outcomes based on their relative distance from the perturber and also on the alignment of acceleration from the perturber with the velocity vector of the moon at the time of close encounter.

In set 2, we simulate the dynamical evolution of exomoons in a more general context and over longer timescales. Each simulation consists of a central star, three giant planets, and moons orbiting all of the planets. All bodies in the system are assumed to be fully formed by the time the simulation starts, and the protoplanetary disk is absent.

The star has solar mass and solar radius. Each simulation includes three giant planets, and the combinations of planet masses (chosen from 0.1, 0.3, 0.5, or 1 M_J) for the three planets in each simulation are sampled thoroughly and with equal frequency in the final statistics (two simulations for systems with equal-mass planets and one simulation for each of the other combinations of planet masses, totaling a set of 68 simulations). Simulations where planet collisions occur are rerun, because planetary collisions could introduce complexities to moon survival. This procedure may introduce bias into the final statistics because roughly a quarter of three-planet systems have planet collisions (Raymond et al. 2010), and many of the planetary collisions happen early on. This leaves the third planet—which participates neither in the collision nor very much in the close encounters—with a very high moon survival rate. The innermost giant planet is placed 5 au from the central star, and all planets are separated by 3.5–4.5 mutual Hill radii in order for instability to set in within a reasonable timescale (Chambers et al. 1996; Marzari & Weidenschilling 2002). In addition, they are given initial eccentricities randomly sampled from 0.02 to 0.1, in order to shift the first onset of instability to an earlier time. To avoid immediate collisions between the giant planets upon orbit crossings, they are also given a small mutual inclination of 001. The giant planets are treated realistically as oblate spheroids, with their quadrupole moment ${J}_{2}=0.0147$ equal to that of Jupiter, and the spin of the planet is allowed to evolve in this simulation set.

Each of the three giant planets in a simulation initially has 10 satellites. Each satellite is placed at a different planetocentric orbital distance from $0.01\,{R}_{{\rm{H}}}$ to $0.35\,{R}_{{\rm{H}}}$ (0.0016–0.24 au for the entire simulation set) on circular orbits. The inner boundary of $0.01\,{R}_{{\rm{H}}}$ is close to Io's orbital distance from Jupiter, a distance far enough away from the Roche limit to avoid tidal mass loss and disintegration due to the tidal field of the planet (Guillochon et al. 2011). The outer boundary lies within 0.5 R_H, the Hill stability limit for prograde satellites (Domingos et al. 2006). The moons are initially on circular and nearly co-planar orbits (0°–002) with their host planet's orbital plane and equatorial plane. All simulations in set 2 are integrated to 1 million years. The final results discussed below are drawn from this full set of simulations unless specified. A subset of the above simulations are further integrated up to 10 million years. The energy error $dE/E$ accumulated in the simulation due to the integrator is always under 10⁻⁴, which is an adequate threshold for multiplanet systems (Barnes & Quinn 2004; Raymond et al. 2010). A couple of the simulations with $dE/E\sim {10}^{-4}$ in the above subsets are rerun or ruled out. All simulations in set 2 are integrated to 100 million years to test the stability of "moons" that end up orbiting the star. This set has a higher rate of exceeding the threshold $dE/E$ (∼20%), but those cases are ruled out. All moons are massless test particles that do not perturb each other or the planets.

During the instability phase of the system, when planets sometimes encounter each other closely and have strong secular interactions, the moons' orbits evolve under the influence of various perturbations. The most influential and in most cases the strongest perturbation comes from the close approach of the perturbing planet to the moons. In addition, as the planets' and moons' orbits evolve through the instability phase, the stability of the moons are affected. Changes in the host planet's semimajor axis, eccentricity, and inclination affect the size of the planet's Hill sphere; changes in the moon's own orbits also cause them to become Hill unstable. A rise in the planets' and moons' inclination can put moons under Kozai-like perturbations. The spin evolution of planets could also affect the stability of moons that are close to the critical semimajor axis, as will be further explained below. Mutual planet secular perturbations also perturb moons. Different types of orbital instability may also have been involved in destabilizing the moons.

4.2.1. Direct Effect of Planetary Close Approach

During a close encounter, planets and moons have strong gravitational interactions over a brief period of time, typically on the order of 10⁻¹–1 year. Therefore, the time evolution of the orbital elements of planets and moons experience an instantaneous change at close encounters, as shown in Figure 4. The planet's spin axis also experiences a small sudden change at a close encounter, causing a slight change in the satellite's plane of nodal precession. But the inclination of the moon from the equatorial plane often changes more dramatically during a close encounter. Depending on the geometry of the encounter, their semimajor axis, eccentricity, and inclination may increase or decrease. Although the close encounter geometry is random, the overall eccentricity of the moon tends to increase rather than decrease. Shown in Figure 4 is a typical case of how the moons' eccentricity evolves during the instability phase. At a very close encounter that is strong enough to perturb most of the moons, their eccentricities increase instantaneously by a large amount. Subsequent milder encounters are able to perturb stable moons that are already eccentric, so then the eccentricity increases or decreases depending on the encounter geometry. Sometimes more encounters push the moons onto highly eccentric orbits or destabilize them.

Zoom In Zoom Out Reset image size

Figure 4 also demonstrates how a very close encounter determines the moon's orbit and dynamical outcome. If a perturbing planet has a mass comparable to the host planet, and it approaches the host planet more closely than some of the moons, moons exterior to its closest approach distance all become destabilized. If the perturbing planet is much less massive than the host planet (i.e., at least three times smaller, as in Figure 4), not all moons become destabilized. This demonstrates again how d_min determines the moon stability limit as in Section 3 and also why close-in moons tend to have higher survival rates. In Figures 4(a) and (b), moons interior to the closest approach distance exhibit a trend of eccentricity growth correlated with the semimajor axes. In Figures 4(c) and (d), the close encounter further randomized the outcome of the previous one. The 0.1 M_J perturbing planet pumps up the eccentricity of the majority of moons exterior to and close to its closest approach distance and destabilizes two of the most exterior moons. Close encounter geometry and moon semimajor axes have similar effects to those in the previous close encounter. In between the two featured close encounter events, there are two close encounters with much larger d_min, far away from the moons, and together with the small mass of the perturbing planet, they cause little perturbation on moons. The more distant moons, due to their closer distance to the perturbing planet, tend to become more eccentric. In the case when not all moons are destroyed, moon survival is not as dependent on the semimajor axis; the close encounter geometry (or the relative distance and velocity between the moon and the perturber) plays a more important role than in the previous case, although the strong gravitational binding of closer-in moons could still give them a somewhat better chance of survival. If the closest approach distance is exterior to the moons, chances are higher that the moons could survive, but, as a general trend, those within the sphere of the perturber's gravitational influence will become more eccentric the larger their semimajor axes. However, the encounter geometry adds some uncertainties. In summary, all factors of survival come into play together, but on longer timescales than shown in Figure 4, the seemingly important role of minimum close encounter distance is smeared out.

Destabilized moons can collide with the host planet due to a small pericentric distance. They can also collide with or become captured by the perturber. Moons that are thrown onto hyperbolic orbits can enter heliocentric orbits or get ejected out of the system as free-floating moons. Ejection is by far the most common outcome. Most of the moons that are scattered onto heliocentric orbits experience further close encounters with the giant planets, because the giant planets themselves are on eccentric orbits. In our 100 million year simulations, 94% of the moons on heliocentric orbits become ejected from the system to beyond 1000 au. Depending on the formation efficiency, free-floating "moons" can outnumber free-floating "planets" (Veras & Raymond 2012).

4.2.2. Perturbation from Planet Orbital Evolution

Planetary orbits evolve chaotically as they scatter off each other and secularly interact in the planet–planet scattering phase. Since planetary orbits set the moon stability criteria, moons are hosted in an environment constantly evolving with respect to the stability limit. Different mechanisms leading to the destabilization of moons can be switched on besides the direct perturbation in planetary close encounters. The following subsections discuss relevant stability criteria.

Hill stability—The size of the Hill sphere depends on the planet's semimajor axis, eccentricity, and inclination (Domingos et al. 2006; Donnison 2010). Planets migrate in/out during the instability phase, or their distance from the star varies within an orbit due to a rise in eccentricity. When their Hill sphere shrinks through this path, moons that become exposed outside of the Hill sphere become destabilized without the direct perturbation of a close encounter. The planet's inclination, coupled with eccentricity, also determines the Hill stability.

Critical semimajor axis—The critical semimajor axis of a planet separates the region where the perturbation from planet oblateness and the central star dominates respectively,

$$\begin{eqnarray}&&{a}_{\mathrm{crit}}={\left(2{J}_{2}{R}_{{\rm{p}}}^{2}{a}_{{\rm{p}}}^{3}{(1-{e}_{{\rm{p}}}^{2})}^{3}\displaystyle \frac{{m}_{{\rm{p}}}}{{m}_{* }}\right)}^{\tfrac{1}{5}},\end{eqnarray} \tag{ 5 }$$

in the co-planar case, where J₂ is the quadrupole moment of the planet, R_p the planet's radius, a_p the planet's distance from the star, m_p the planet's mass, and m_* the star's mass (Kinoshita & Nakai 1991; Nicholson et al. 2008; Deienno et al. 2011). In this work, moons often have inclinations high enough for a_crit to depend also on their inclinations. Well inside a_crit, the orbital angular momentum axis of the moon precesses rapidly around the planet's spin, and well outside a_crit, it precesses around the orbital angular momentum axis of the planet. Moons close to a_crit precess about local Laplace planes that lie between the spin and orbital plane of the planet. As the a_crit of a planet co-evolves with its orbit, the center of precession of moons switches between different planes, thus changing the moons' inclinations relative to their local Laplace planes. Some moons wander close to a_crit in systems with high obliquities and could possibly be destabilized through this path (Tremaine et al. 2009), although planets evolving to high obliquities through planet–planet scattering also often acquire high inclinations and eccentricities, making it hard to pin down the exact cause of the moon's orbital instability. Also, in high obliquity systems, the moons' inclination relative to the center of precession could change significantly as a_crit changes, giving them stability or destabilizing them. When a planet evolves to high inclinations, a path may be open to destabilizing its moons. As planet inclination reaches ∼40°, the Kozai mechanism starts to act on the moon that then stays close to and precesses around the orbital plane, causing the eccentricity and inclination to oscillate in a coupled manner.

The scenarios above do not involve any direct perturbation from planetary close encounters but are caused by evolved planetary orbits induced by close encounters plus mutual secular perturbation. As a side note, however, only a small fraction (<5%) of surviving planets or planets within 100 au from the central star in planet–planet scattering have inclinations above 40° (Chatterjee et al. 2008; Raymond et al. 2010), so this is not an important mechanism to destabilize the orbits of moons in the planet–planet scattering scenario.

4.2.3. Evolved Moon Orbit

4.2.4. Secular Evolution of Moon Post-instability

Table 1 summarizes possible dynamical outcomes for moons and the respective probability. The moons that remain planet bound may stay around their original host planet, around a new host planet that captures it away from the original, or around a planet that itself is ejected as a free floater. Other moons may stay on heliocentric orbits, or collide with the host planet, the perturbing planet, or the star. Multiple close encounters throughout the planet–planet scattering phase are very efficient for removing moons. Seventeen percent of the moons in the orbital range of 0.01–0.35 R_H remain bound to the host planet after 1 million years of integration. This survival rate with respect to the entire moon population signifies how much of its parameter space allows for moons to be stable. Inferring from the above numbers and collapsing the parameter space to the innermost region, the moon stability limit is ∼<0.1 R_H. For a similar argument from a different perspective, of the 17%, 88% are initially within 0.1 R_H. In other words, planet–planet scattering truncates the primordial moon disk at least as close as 0.1 R_H. If considering only moons interior to 0.04 R_H, inside which the regular moons of the solar system gas giant planets are located, the moon survival rate amounts to ∼42%, forecasting a decent survival chance for moons on Galilean-moon-like orbits. The rest of the destabilized moons (83%) are removed via various paths, as shown in Table 1. Ejection out of the system is the most common outcome (∼41%), given the small mass of moons relative to the planets (here moons are massless). The collision of moons with the other big bodies is also equally dominant (∼37%); among them, collision with the planets than with the central star are more frequent due to the proximity of the moons to the planets, especially during planetary close encounters. The above two paths remove the majority of the moons. The remaining unstable moons stay in different places in the system; ∼2% of them are captured by the perturbing planet, and ∼2% of them orbit planets that get ejected out of the system as free floaters. The capture of moons via exchange when two planets approach each other closely is a much less efficient way of producing irregular satellites than capture from circumstellar material, probably because of the large number of the latter and the fact that it requires a particular configuration—small d_min. After 1 million years, ∼22% of all moons initially stable around their host planets have destabilized and ended up orbiting the star; however, further integration to 10⁸ years removed most of these heliocentric "moons," with the rate reduced to ∼1%. The moons removed by 10⁸ years are incorporated into the final statistics in this section.

A simulation of 10⁶ years is 1–2 orders of magnitude shorter than required for a subset of simulations to finish the instability phase. A sample subset of simulations taken from those that have not ended the instability phase is further integrated up to 10⁷ years. The survival rate of moons in the subset decreases to half for moons across all orbital distances. The number of close encounters and the magnitude of their perturbations between 0 and 10⁶ years and ${10}^{6}{\textstyle \mbox{--}}{10}^{7}\,{\rm{y}}{\rm{e}}{\rm{a}}{\rm{r}}{\rm{s}}$ are comparable. Although the temporal length of the latter makes the perturbation from close encounters more diffused, the change in survival rate and how perturbed the surviving moons' orbits are not insignificant.

The key factors for moon survival include various system parameters. Parameters not observable in actual systems include the initial orbital distance of moons and properties of planetary close encounters, such as closest approach distance, number of close encounters, and close encounter geometry. Observables include planet mass and planetary orbits (semimajor axis, eccentricity, inclination, and obliquity). All of the above parameters work together in determining the dynamical outcome and stability of moons, so how the moons' orbits would change can only be predicted with a single parameter in a probabilistic sense. Moreover, the long length of the simulations and the large number of close encounters (the average number of close encounters in a 10⁶ year simulation is 224) or the large removal rate of moons probably smears out any clear correlation between some of the parameters and the moon survival rate within the simulated set. Figures 5–7 compare various factors with the moon survival rate and dynamical outcome. The stability limit of moons around the planet is defined at the location of the outermost surviving moon in the Hill sphere. Moon survival rates are nearly equivalent to the stability boundary because moons are all set at different orbital distances around each planet, and due to strong perturbations by close encounters, the survival of an outer moon often ensures that moons interior to it also survive. However, the encounter geometry sometimes adds uncertainties to the relation between moon survival and the stability boundary. Note that parameters related to close encounters do not correlate with moon survival rates in the full-length simulations, unlike what has been shown in single flyby events in Section 3.

4.4.1. Moon Initial Semimajor Axis

a_i versus survival rate—The planetocentric distance of moons determines their survival rate. Quite intuitively, inner moons are more likely to survive because they have several advantages. As a typical case in Figure 4, after the perturber passed by, the inner moons usually experience less change in their orbits than do the outer moons. Closer-in moons also have a lower probability of being perturbed by an encounter since moons are not sensitive to encounters far from them. In Figure 4, in both of the close encounter events, most of the moons interior to the perturbing planet experience little perturbation while those exterior to and near it become destabilized or highly eccentric. Different encounter geometries do add some randomness to the trend.

As shown in Figure 5, the survival rate drops rapidly as a moon's semimajor axis increases. Moons at Io-like distance (0.01 R_H) have a survival probability of ∼0.5–0.6. Beyond 0.1 Hill radii, the chances of survival are very low (<0.1). The survival rate for the Galilean satellites at the end of 10⁶ years could be higher than ∼0.3–0.5 as predicted by the simulations (the outermost satellite, Callisto, is ∼0.04 R_H from Jupiter), because in the simulations, they orbit around planets of all masses but in the dynamical history of the solar system, Jupiter was the most massive planet by the time the satellites formed (Canup & Ward 2009). Their survival rate could also be lower if the solar system has not ended its instability phase by 10⁶ years.

Zoom In Zoom Out Reset image size

a_i versus dynamical outcome—The dynamical outcomes of moons are predetermined in a probabilistic sense by their initial location within the Hill sphere. Figure 6 shows the cumulative distribution in the initial semimajor axis of moons that have different dynamical outcomes. The overall population (black line) divides the planet-bound moons and the unbound moons. Not surprisingly, the planet-bound moons tend to occupy the inner part of the Hill sphere, with their curves consistently lying on the left division of the overall population. Among them, those that stay stable around an ejected free-floating host planet lie in the most interior region, likely due to the fact that they need a sufficiently strong gravitational binding and enough separation with the perturber to survive the planet-ejecting close encounters. By contrast, moons that become heliocentric, including those that later collide with or orbit the star and those that are ejected out of the system, lie consistently on the right-hand side of the overall population, and the three curves are very close to each other, since right after the close encounter that destabilize them, they share the same outcome.

Zoom In Zoom Out Reset image size

4.4.2. Planet Mass

The mass of the host planet determines the extent of its moons' gravitational binding, which if large compared to perturbations, the moons are more likely to survive. Therefore, the more massive the planet, the more stable its moons usually are. Figure 7 shows the cumulative distribution of the outermost surviving moon in semimajor axis on planets with different masses. Moons around more massive planets have wider distributions across the Hill sphere; in other words, they are more likely to have a larger stability boundary. To have moons surviving beyond 0.1 R_H, the probability for planets with M_J, 0.5 M_J, 0.3 M_J, and 0.1 M_J are ∼35%, ∼20%, ∼10%, and ∼5% respectively. Around 0.5 and 1 M_J planets, the 80th percentile of the outermost stable moons reaches 0.1 R_H, whereas around 0.1 and 0.3 M_J planets, they are under a much tighter limit of 0.04 R_H. All planet masses from 0.1 to 1 M_J permit moons on Galilean satellite orbits to have good survival rates.

Zoom In Zoom Out Reset image size

Moons with different dynamical outcomes have experienced different amount of perturbations, so they occupy different parts of the orbital parameter space; in other words, they are dynamically distinguishable. As a general trend, the more they deviate from their initial orbits, the more they have been perturbed.

4.5.1. Δa

4.5.2. Eccentricity

Quite intuitively, moons that stay stable around their original host planet have less excited orbits than those that are destabilized. Moons that remain stable around bound or free-floating planets have lower average eccentricities, and those captured by another planet or go into heliocentric orbits are much more eccentric.

Figure 8 shows the cumulative fraction in eccentricity of moons of different dynamical outcomes. Primordial moons bound to the host planets are, to no surprise, the least eccentric population, since their orbits receive the least perturbation and change the least (Figure 8). Of these, 44% have eccentricities under 0.01, and the average eccentricity of the perturbed ones is 0.31. Moons become eccentric once an encounter strong enough to perturb them occurs, and the probability to recover via subsequent encounters is extremely low; more violent encounters usually only increase eccentricity by a large amount, and milder encounters have a chance to damp down the eccentricity, but only by a moderate amount. The average eccentricity of moons bound to free-floating planets is ∼0.32, higher than moons around stable planets, and the 80th percentile reaches e ∼ 0.6; as a comparison, the 80th percentile of primordial moons around stable planets reaches e ∼ 0.47. Moons that are captured away by a different planet have a higher average eccentricity of 0.57, and they are more evenly distributed across the eccentricity spectrum except at the low end with e < 0.2, because their orbits are randomized by the capture process.

Zoom In Zoom Out Reset image size

Moons stripped away from the planet to orbit the central star have the highest average eccentricity of 0.68; the 80th percentile reaches as high as e = ∼0.91, much higher than all of the planet-bound populations, due to the higher tolerance for eccentricity in the heliocentric system than in the planetocentric system. For both simulations lasting 10⁶ and 10⁸ years, all of the heliocentric moons in the set form a thermal eccentricity distribution (Figure 8),

$$\begin{eqnarray}&&f(e)=2e.\end{eqnarray} \tag{ 6 }$$

This distribution resembles that of wide binary stars and suggests that they have been sufficiently pulled and kicked around for the energy to approximate an equilibrium state (Jeans 1919; Kouwenhoven et al. 2010). However, as a consequence of extremely high eccentricities and the fact that the planets are also on eccentric orbits, interactions are often strong and frequent, making heliocentric moons extremely vulnerable to subsequent removal. By the end of 10⁸ years, only ∼6% of them remain in the system (a < 1000 au). However, if a planetesimal belt exists at the edge of the planetary system, it could help trap moons on the way to ejection, turning them into Kuiper-belt-like (KBO-like) objects (Raymond et al. 2009). Another possibility for producing KBO-like objects out of ejected moons is that, if some instabilities happen before the dissipation of the gas disk, gas drag could act to trap a heliocentric moon by lifting its pericenter from the grasp of the perturbing planets (Raymond & Izidoro 2017). Of the majority of moons venturing beyond 1000 au (classified as "ejected" in this work), ∼35% are able to reach at least 5000 au with eccentricities ≲1 within 10⁸ years. Those moons might have a chance to damp down their eccentricities through the galactic tide and become Oort-cloud-like objects (Tremaine 1993; Morbidelli 2005). Twenty percent of the test simulations for Oort-cloud-like candidates with fractional angular momentum change larger than 10⁻⁵ are ruled out because most of the moons that reach the orbital distance where the galactic tides becomes effective have eccentricities $\lesssim {10}^{-5}$ below 1.

4.5.3. Inclination

Like eccentricity, inclination is an indicator of the dynamical history of moons. Figure 9 shows the cumulative fraction of moons in inclination of different dynamical outcomes. The stable populations, the host-planet-bound and free-floater-bound moons, and the planets are again the least evolved in i, as in a and e. The more unstable moon populations such as the heliocentric moons and the captured ones have larger inclinations; 45% of them are above 40° as seen from the planet's orbital plane. A much larger fraction of them are retrograde from other populations: 23% for moons captured by the perturbing planet and 16% for heliocentric moons by 10⁶ years. In comparison, only ∼1% of the planet–bound moons are on retrograde orbits.

Zoom In Zoom Out Reset image size

Surviving moons have relatively quiet orbits, although still significantly more perturbed than the solar system's regular moons. Of the stable primordial moons, 86% within a_crit are under ∼20° as seen from the host planet's equatorial plane and 73% are above 1°. As for those exterior to a_crit, 81% are under 20° from the host planet's orbital plane. The inner population is significantly less perturbed than the outer one. As for moons stable around ejected free-floating planets, 80% are under 20° from the initial reference plane, but the number of such systems is small, making statistical interpretation difficult. Both of the primordial populations have under 2% on retrograde orbits.

Moons captured by the perturber are more evenly distributed across the inclination spectrum, including a significant portion on the retrograde side, and in general they have higher inclinations than moons bound to their original host planet; this tendency is similar to the simulated capture of satellites from planetesimals (Nesvorný et al. 2007).