<figure><img src="image_1.png" /></figure>

# Innocent Bystanders: Orbital Dynamics of Exomoons During Planet-Planet Scattering

Yu-Cian Hong¹,², Sean N. Raymond³, Philip D. Nicholson¹, and Jonathan I. Lunine¹,²

¹Astronomy Department, Space Sciences Building, Cornell University, Ithaca, NY 14853, USA; ycsylva@gmail.com

²Carl Sagan Institute, Space Sciences Building, Cornell University, Ithaca, NY 14853, USA

³Laboratoire d'astrophysique de Bordeaux, Université de Bordeaux, CNRS, B18N, allée Geoffroy Saint-Hilaire, F-33615 Pessac, France

Received 2017 October 2; revised 2017 December 4; accepted 2017 December 8; published 2018 January 10

## Abstract

Planet-planet scattering is the leading mechanism to explain the broad eccentricity distribution of observed giant exoplanets. Here we study the orbital stability of primordial giant planet moons in this scenario. We use N-body simulations including realistic oblateness and evolving spin evolution for the giant planets. We find that the vast majority (~80%-90% across all our simulations) of orbital parameter space for moons is destabilized. There is a strong radial dependence, as moons past ~0.1 $R_{Hill}$ are systematically removed. Closer-in moons on Galilean-moon-like orbits (<0.04 $R_{Hill}$) have a good (~20%-40%) chance of survival. Destabilized moons may undergo a collision with the star or a planet, be ejected from the system, be captured by another planet, be ejected but still orbiting its free-floating host planet, or survive on heliocentric orbits as “planets.” The survival rate of moons increases with the host planet mass but is independent of the planet’s final (post-scattering) orbits. Based on our simulations, we predict the existence of an abundant galactic population of free-floating (former) moons.

**Key words:** planets and satellites: dynamical evolution and stability

## 1. Introduction

Given that each of the solar system's giant planets hosts at least one large natural satellite, the presence of moons is anticipated around giant exoplanets. The potentially diverse environments on exomoons, and the clues they may provide to planet formation models, make them subjects worthy of research.

The potential habitability of exomoons is affected by a number of processes, including atmospheric dynamics, stellar illumination, tidal heating, planetary magnetic fields, orbital configurations, etc. (Williams et al. 1997; Kaltenegger 2010; Tinney et al. 2011; Heller 2012; Forgan & Kipping 2013; Heller & Barnes 2013, 2014; Heller & Zuluaga 2013; Hinkel & Kane 2013; Dobos et al. 2017). Exomoons have more paths to habitable configurations than planets. Besides being in the liquid water habitable zone of the star (Kasting et al. 1993), effects such as tidal heating provide an opportunity for moons outside the stellar habitable zone to be heated (Reynolds et al. 1987; Scharf 2006; Heller & Zuluaga 2013; Dobos et al. 2017).

Several techniques appear capable of detecting exomoons (see Section 6 of Heller et al. 2014 for a review). The transit technique has the prospect of observing sub-Earth-sized moons (Sartoretti & Schneider 1999; Agol et al. 2005; Holman & Murray 2005; Kipping 2009; Kipping et al. 2009, 2012; Heller 2014; Kipping et al. 2015) and has discovered a Neptune-sized exomoon candidate (Teachey et al. 2017). Microlensing can detect moons down to 0.01 Earth masses (Bennett & Rhie 1996, 2002; Han & Han 2002; Han 2008) and has discovered a sub-Earth mass exomoon candidate around a free-floating planet (Bennett et al. 2014). Its broader sensitivity with heliocentric distances may yield a better chance of exomoon detection. Direct imaging may be able to detect bright, tidally heated exomoons (Peters & Turner 2013).

The frequency of the occurrence of exomoons depends on their formation and dynamical evolution. The Hill stability criterion, Roche limit, stellar tidal stripping, tidal decay, planet migration, direct planet perturbation, and planetary close encounters are among the mechanisms that can destabilize

the orbits of exomoons (Holman & Wiegert 1999; Barnes & O'Brien 2002; Domingos et al. 2006; Donnison 2010; Namouni 2010; Sasaki et al. 2012; Frouard & Yokoyama 2013; Payne et al. 2013; Spalding et al. 2016; Kane 2017). The stability of distant satellites ($a_s > 0.1 R_H$), which would be classified as irregular satellites in the solar system, in planet-planet scattering is investigated by Gong et al. (2013) and Nesvorný et al. (2007).

In this paper, we investigate the orbital dynamical behavior and the stability of primordial satellites as close-in as Io to Jupiter (~0.01 $R_H$), in scenarios where planets gravitationally scatter off each other. Planet-planet scattering has been considered as the most viable candidate mechanism for explaining the prevalence of eccentric orbits of extrasolar planets and reproduce well their observed eccentricity distribution (Rasio & Ford 1996; Weidenschilling & Marzari 1996; Lin & Ida 1997; Marzari & Weidenschilling 2002; Adams & Laughlin 2003; Chatterjee et al. 2008; Ford & Rasio 2008; Jurić & Tremaine 2008; Raymond et al. 2010). Therefore, it is important to test the orbital stability of planetary satellites in this context. In the planet-planet scattering scenario, planets are hypothesized to form closely packed, then they perturb each other and eventually enter an instability phase when they undergo orbit crossing and close encounters. The instability is ended by the removal of some planets in the system by ejection or collision. In the instability phase, the planetary satellites experience strong gravitational perturbations from the perturbing planets. Various sources of perturbations also affect the satellites, such as secular perturbations from the perturbing planets (as opposed to the host planet of the satellites) and stellar perturbations through the change of the host planet's spin and orbits in the instability phase. The orbital parameter space of moons in this work covers close-in regions where planet oblateness plays a major role in moon stability, as an already well-established fact, and the absence of planet oblateness will cause an unrealistic orbital instability effect (Hong et al. 2015). This work also simulates moons within and beyond the critical
semimajor axis (0.04 Hill radii for Jupiter) of the planets, where planet spin can affect moon stability (Tremaine et al. 2009), and moons up to 0.35 Hill radii, where prograde moons can be stable.

Section 2 studies the inner workings of planetary close encounters and discusses factors related to close encounters that affect moon stability, by simulating single planet mutual flyby events. Section 3 discusses the orbital evolution of moons under various sources of perturbations, the moon survival rate and its relevant factors, and the dynamical outcome of moons in planet-planet scattering in full simulations of length 1–100 million years.

## 2. Numerical Method

This work uses the N-body symplectic integrator Mercury (Chambers 1999). All simulations use the Bulirsch-Stoer algorithm, the most accurate for simulating bodies that perturb each other closely, though the slowest in the package, in order to ensure integration accuracy. The conservative version of the Bulirsch-Stoer integrator is used, which is appropriate in cases where the force does not depend on velocity. The spin and oblateness of the planet are adapted into the code in order to correctly simulate the stability of planetary satellites.

### 2.1. Oblateness of the Planet

The planets are simulated as oblate by the addition of the quadrupole moment ($J_2$) 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.

### 2.2. Spin Evolution of the Planet

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

$$F = \frac{-G M m}{r^2} \left\{ \hat{r} - 3 J_2 \left(\frac{R}{r}\right)^2 \left[ (5(\hat{r} \cdot \hat{s})^2 - 1)r - 2(\hat{r} \cdot \hat{s})\hat{s} \right] + \cdots \right\}, \quad (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

$$\frac{\mathbf{r} \times \mathbf{F}}{l} = \frac{d}{dt}(\Omega \hat{s}), \quad (2)$$

where $\Omega$ 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:

$$\frac{d\hat{s}}{dt} = \frac{-3 G M m}{r^3} \frac{J_2}{\lambda \Omega} (\hat{r} \cdot \hat{s})(\hat{r} \times \hat{s}), \quad (3)$$

where $\lambda$ is the normalized moment of inertia. $\lambda$ 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:

$$\frac{d\hat{s}}{dt} = \left( \frac{-3}{2} \frac{n^2}{\Omega} J_2 \lambda \right) (\hat{l} \cdot \hat{s}) (\hat{l} \times \hat{s}) \quad (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.

## 3. Result: Single Close Encounter

We start by considering the effect of a single planet-planet encounter on the stability of one planet's moons.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Planet minimum close encounter distance vs. moon stability boundary. The stability boundary is set where the most distant surviving moon with *e* < 0.5 is located around the planet. Dots with different colors represent simulations with different perturber masses as shown in the legend. Each dot is obtained by taking averages of the simulations with the same *d*<sub>min</sub>. The trend in the plots with perturber mass 0.5–2*M<sub>J</sub>* indicates that *d*<sub>min</sub> plays the most important role in determining the moon stability boundary, and the different slopes of each curve corresponding to different perturber masses show that perturber mass is a secondary factor.</figcaption></figure>

## 3.1. Simulation Setting

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$^{-1}$, 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$^6$ 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.

## 3.2. Results

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$^{-1}$) 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

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Minimum close encounter distance vs. average number of moon captures. Figure legend as in Figure 1.</figcaption></figure>

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.

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_J$ and 2 $M_J$), the stability boundary is 0.6–0.8 times $d_{min}$ for $d_{min} < 0.05$ 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 and 0.35 $R_H$ (~0.003–0.12 au) are stable. For simulations with a 0.1 $M_J$/0.5 $M_J$/1 $M_J$ perturber, when $d_{min}$ is greater than ~0.05/0.1/0.1 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–2 $M_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
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Close encounter geometry vs. moon dynamical outcome. The same planetary close encounter event is split into four subfigures. Each subfigure features a different moon in the system. In each subfigure, the host planet is represented by the cross sign at the origin, the trajectory of the perturbing planet is represented by empty circles, and the trajectory of the moon represented by empty squares. Their trajectories change in color from purple to yellow as time evolves, as is labeled by the vertical color axis (time in days). The solid black dots represent the location of the perturber and the moon at the time of closest approach. See the main text for further simulation details.</figcaption></figure>

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

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.

# 4. Result: Full Integration

## 4.1. Simulation Setting

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 0°.01. 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_H$ to $0.35 R_H$ (0.0016–0.24 au for the entire simulation set) on circular orbits. The inner boundary of $0.01 R_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°–0°.02) 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^{-4}$, 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.

## 4.2. Moon Orbital Evolution Under Different Perturbations

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}$–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.

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_1$ 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
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Evolution of the moon's *a* and *e* early in two early close encounter events. (a) and (b) features the first close encounter, and (c) and (d) feature the second one in the same simulation. The two events are two consecutive effective close encounters. Each colored curve represents a moon with their initial orbital distance in $R_H$ labeled in the legend. All moons orbit the same 0.5 $M_J$ host planet. In (a) and (b), at ~3200 years, a 0.1 $M_J$ perturbing planet (blue diamond) makes a very close approach to the moon-hosting planet with $d_{min} = 0.057$ au. (c) and (b) feature a subsequent close encounter at ~3850 years with $d_{min} = 0.014$ au.</figcaption></figure>

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,

$$a_{\text{crit}} = \left( 2 J_2 R_p^2 a_p^3 (1 - e_p^2)^3 \frac{m_p}{m_*} \right)^{\frac{1}{5}}, \quad (5)$$

in the co-planar case, where $J_2$ 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_{\text{crit}}$ to depend also on their inclinations. Well inside $a_{\text{crit}}$, the orbital angular momentum axis of the moon precesses rapidly around the planet's spin, and well outside $a_{\text{crit}}$, it precesses around the orbital angular momentum axis of the planet. Moons close to $a_{\text{crit}}$ precess about local Laplace planes that lie between the spin and orbital plane of the planet. As the $a_{\text{crit}}$ of a planet coevolves 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_{\text{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_{\text{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

Besides the evolution of planetary orbits, the evolution of the orbits of moons can also destabilize or change themselves by evolving beyond the stability limit. Oftentimes, the orbits of moons evolve when there is an effective close encounter, although mutual interactions among planets when they are outside each other's Hill sphere also affect moons. As their orbits change, they can become exposed outside of the Hill sphere and thus become destabilized or more stable. They can also shift across $a_{\text{crit}}$ as they radially migrate in and out, which causes changes in their inclinations relative to the center of orbital precession (local Laplace plane). Close encounters can

<table><caption>**Table 1**
Probability of Each Dynamical Outcome for All Moons in the 10$^6$–10$^8$ Year Simulations in Set 2</caption><thead><tr><th>Outcome (Surviving)</th><th>Probability(%)</th></tr></thead><tbody><tr><td>Host Planet Bound</td><td>17</td></tr><tr><td>Captured by Perturber</td><td>2</td></tr><tr><td>Free-floater Bound</td><td>2</td></tr><tr><td>Heliocentric</td><td>1</td></tr></tbody><thead><tr><th>Outcome (Removed)</th><th>Probability(%)</th></tr></thead><tbody><tr><td>Ejected</td><td>41</td></tr><tr><td>Collision with Planet</td><td>20</td></tr><tr><td>Collision with Star</td><td>17</td></tr></tbody></table>

throw them onto random inclinations relative to their center of precession and affect their stability.

### 4.2.4. Secular Evolution of Moon Post-instability

After the planetary system has been restored to stability via the removal of one or two planets, there are no more planetary close encounters. Most of the moons simply evolve in a stable manner, with their semimajor axis staying fixed, and eccentricity and inclination undergoing sinusoidal oscillation periodically as a consequence of the precession of their orbital angular momentum around the local Laplace plane. For some moons inside $a_{\text{crit}}$ with inclinations greater than 40° from the equatorial plane, the inclination and eccentricity will be slowly oscillating out of phase with each other by 180°, indicating that the Kozai mechanism is in action. Moons with inclinations under 40° undergo secular evolution, with a fixed semimajor axis and periodically oscillating eccentricities and inclinations. The same argument applies to moons exterior to $a_{\text{crit}}$ inclined relative to the orbital plane. Those that undergo Kozai evolution with high inclinations could become highly eccentric or unstable in the longer term.

### 4.3. Statistics

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^8$ years removed most of these heliocentric “moons,” with the rate reduced to ~1%. The moons removed by $10^8$ years are incorporated into the final statistics in this section.

A simulation of $10^6$ 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^7$ 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^6$ years and $10^6–10^7$ years 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.

## 4.4. Survival Factors

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^6$ 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

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Moon initial semimajor axis ($a_i(R_H)$) vs. moon survival rate. The moon survival rate is predetermined by the moon's initial position in the host planet's Hill sphere. In the semimajor axis range 0.01–0.04 $R_H$, which is roughly where the four Galilean moons are located, the survival rate drops exponentially.</figcaption></figure>

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^6$ 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^6$ years.

$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
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Moon initial semimajor axis $a_i$ ($R_H$) vs. moon dynamical outcome. Where moons will end up in planetary close encounters are also predetermined by their initial semimajor axis, not definitely but probabilistically. In the color legend from top to bottom are moons (1) bound to the original host planet (red), (2) captured by the perturber (green), (3) bound to the free-floating planet (blue), (4) orbiting the star (pink), (5) ejected out of the system as free-floating moons (brown), (6) colliding with the star (purple), (7) colliding with planets (light blue), followed by the overall moon population in black.</figcaption></figure>

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.

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

## 4.5. Final Orbits of Moons

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.

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Planet mass vs. moon stability limit. The axes of the figure represent the cumulative fraction of planets in the semimajor axis of their outermost surviving moon (in $R_H$). Different colors represent planets with different masses in $M_J$.</figcaption></figure>

### 4.5.1. $\Delta a$

Planet-bound moons have low tolerance for radial migration. A large amount of change in the moons' semimajor axis often destabilizes it. Moons that stay bound to stable host planets experience no dramatic changes in semimajor axis. Eighty percent of them migrated outward/inward by less than 46%/32% in semimajor axis (or 90% by less than 73%/45%). Moons stably orbiting ejected free-floating planets (2%) are most tightly clustered at small semimajor axis (in au). The 80th percentile reaches as far as ~0.01 au. Their small orbital distance is very likely a result of the same processes identified in the previous section. In comparison, moons that have ended up on heliocentric orbits have a wider distribution in semimajor axis than the surviving planets.

### 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 \sim 0.6$; as a comparison, the 80th percentile of primordial moons around stable planets reaches $e \sim 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
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Cumulative fraction of final eccentricity of stable planets and moons with different dynamical outcomes. In the legend, different colors represent different populations of moons. For planet-bound moons, those bound to the original host planet are represented in red, those captured by the perturber in green, and those bound to free floaters in blue. Moons on heliocentric orbits by the end of $10^6/10^8$ years are represented in pink/purple. Planets surviving around the central star are represented in black.</figcaption></figure>

end with $e < 0.2$, because their orbits are randomized by the capture process.

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 \sim 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^6$ and $10^8$ years, all of the heliocentric moons in the set form a thermal eccentricity distribution (Figure 8),

$$f(e) = 2e. \qquad (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^8$ 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 $\lesssim 1$ within $10^8$ 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^{-5}$ are ruled

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Cumulative distribution of final inclination of stable planets and moons with different dynamical outcomes. The figure legend is identical to Figure 8. Inclinations are measured relative to the initial zero inclination plane of the system, if not specified.</figcaption></figure>

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^6$ years. In comparison, only ~1% of the planet-bound moons are on retrograde orbits.

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

## 5. Occurrence of Free-floating Exomoons

Here we calculate the abundance of free-floating exomoons and estimate its potential contribution to the vast number of free-floating objects estimated from microlensing observations,
since planet–planet scattering appears to be very efficient at ejecting moons. Sumi et al. (2011) obtained an observed frequency of free-floating planets per main-sequence star of $N_{\text{ff}}/N_{\text{star}} = 1.8^{+1.7}_{-0.8}$, and Mróz et al. (2017) found $N_{\text{ff}}/N_{\text{star}} = 0.25$ for Jupiter-mass free floaters. The equation below adopts a similar method to Veras & Raymond (2012) to estimate the number of free-floating moons:

$$ \frac{N_{\text{moon,ff}}}{N_{\text{star}}} = f_{\text{gp}} \times f_{\text{system,unstable}} \times n_{\text{gp,eject}} \times n_{\text{moon,eject}}, \quad (7) $$

where $N_{\text{star}}$ is the total number of stars, $f_{\text{gp}}$ is the fraction of stars with giant planets, $f_{\text{system,unstable}}$ is the fraction of planetary systems that go unstable, $n_{\text{gp,eject}}$ is the averaged number of ejected giant planets per system, and $n_{\text{moon}}$ is the averaged number of moons per system.

The giant planet occurrence rate $f_{\text{gp}}$ has been measured by radial velocity surveys to be ~0.1–0.2 for Sun-like stars (Cumming et al. 2008; Mayor et al. 2011; Rowan et al. 2016) but with a stellar-mass dependence such that low-mass stars are deficient in gas giants (Johnson et al. 2007; Lovis & Mayor 2007; Winn & Frabrycky 2015; Wittenmyer et al. 2016). Given the predominance of low-mass stars by number, we expect the average $f_{\text{gp}}$ to be in the 1%–10% range. The number of giant planets that participate in a given instability, $n_{\text{gp,unstable}}$, must be at least two and may be larger. We assume an average of 2–4. To match the observed eccentricity distribution of exoplanets, the fraction of unstable planetary systems must be very high. $f_{\text{system,unstable}}$ is at least 75% but is more likely more than 90% (Raymond et al. 2010, 2011). $n_{\text{moon,eject}}$ equals the number of moons per system × ejection rate. The former is unknown, and we assume a value of 1–5; the latter is roughly 40% from the result of simulation set 2. Therefore, the occurrence of moons per star is $O(0.01–1)$, which predicts a galactic population of free-floating (former) moons that may be as abundant as stars. WFIRST can detect such objects down to ~0.1 Earth masses (Spergel et al. 2015).

Another interesting population of moons that sit stably on ejected planets has an occurrence probability of order 1% per system based on the theoretical result in this work. Replacing $n_{\text{moon,eject}}$ with the fraction of free-floating giant planets that carry moons (~0.08) and using the same estimate as above for other parameters yield an occurrence rate of $O(10^{-3}-10^{-2})$ per star. If using the observed frequency of free-floating planets per star (0.25; Mróz et al. 2017) and multiplying it with the fraction of moon-bearing free-floating planets from the simulations (~0.08), the occurrence of moon-bearing free-floating planets per star is $O(10^{-2})$. Both approaches yields probabilities that are not insignificant.

## 6. Conclusion

We have directly simulated the survival of exomoons during giant planet scattering. To summarize, most moons are unstable during planet–planet scattering, and they have rich dynamical outcomes. Planet–planet scattering is a very efficient way to destabilize moons, mainly because of the close approach of the relatively massive perturbing planet compared to moons, and second because planets themselves also experience strong mutual perturbations induced by the instability even outside of close encounters. Only ~10%–20% of the moons within $0.35 R_{\text{H}}$ remain bound to the original host planet; in other

words, only ~10%–20% of the orbital parameter space in the Hill sphere allows for moons to be stable. Close-in moons on Galilean-moon-like orbits ($0.01–0.04 R_{\text{H}}$) have at least twice higher survival rates than the whole population. The majority of moons (~80%–90%) collide with big bodies or become ejected out of the system. Given the high probability of ejection, there is likely a population of free-floating objects in interstellar space that were born as moons of giant planets, with an occurrence rate of $O(0.01–1)$ per star. Future microlensing detection is capable of observing such objects down to ~0.1 Earth masses. A very tiny fraction of moons live on interesting orbits, such as around free-floating planets (occurrence rate $O(10^{-3}-10^{-2})$ per star), around a perturbing planet as a captured object, or on heliocentric orbits. As moons around free-floating planets tend to be very close-in and eccentric in the planet–planet scattering scenario, there may be a good chance of having tidally heated Io-like moons around free-floating planets. The simulated number of moons is not abundant enough to explain the efficiency of producing Trojan-like objects from moon scattering, but intuition and the simulation results for moons on heliocentric orbits suggest that such objects can have a hard time surviving the instability phase. Since planet–planet scattering tends to destroy already-formed terrestrial planets (Veras & Armitage 2005; Matsumura et al. 2013; Carrera et al. 2016) or their building blocks (Raymond et al. 2011, 2012), a possible way to form terrestrial planets could be through the heliocentric moons remaining in the system after planet–planet scattering, although the simulation results do not predict it as a highly efficient mechanism (probability ~0.01).

Properties of the close encounters, planets, and moons all come into play together in determining the stability and dynamical outcomes of moons. Among all factors, the planet mass and the initial orbits of moons stand out as clear predictors of a moon’s final outcome, though only in a probabilistic sense. The source region of moons with different dynamical outcomes can be anticipated to some degree.

As to the final orbits of moons, stable moons experience limited radial migration, so it is not a common phenomenon for moons to swap their orbits. However, moons do move in and out of the planet’s Hill stability limit and the critical semimajor axis during the instability. For the solar system gas giant planets, moons interior to the critical semimajor axis are co-planar with the equator and have tiny eccentricities (Deienno et al. 2011), whereas moons exterior to it are inclined and eccentric (Jewitt & Haghighipour 2007). In contrast, in the bulk simulated extrasolar systems, moons interior to the critical semimajor axis have a tiny probability, $\lesssim 0.1$, of having unperturbed orbits like in the solar system, hinting that solar system gas giant planets might not have experienced planetary close encounters, or that they experienced very few and very mild close encounters. In the region exterior to the critical semimajor axis and interior to the Hill stability limit for prograde moons, moons have a negligible probability, ~0.01, of having unperturbed orbits, just like in the solar system. The solar system irregular moons on prograde orbits could have mixed origins of primordial populations and external capture from circumplanetary material. However, tests in the planet–planet scattering scenario show that moons on retrograde orbits are not produced very efficiently, therefore supporting the origin of retrograde irregular satellites via other mechanisms
such as capture from circumplanetary material (Nesvorný et al. 2007).

We thank John Chambers, Tom Quinn, and John Armstrong for helpful discussions. We are grateful to an anonymous referee for a constructive review. This work was supported by the JWST (NASA grant NNX12AK016) and Cassini (subcontract JPL 64969) projects. S.N.R. thanks the Agence Nationale pour la Recherche for funding via grant ANR-13-BS05-0003-002 (MOJO) and the NASA Astrobiology Institute's Virtual Planetary Laboratory Lead Team, funded via the NASA Astrobiology Institute under solicitation NNH12ZDA002C and cooperative agreement No. NNA13AA93A.

## References

Adams, F. C., & Laughlin, G. 2003, *Icar*, 163, 290

Agol, E., Steffen, J., Sari, R., et al. 2005, *MNRAS*, 359, 567

Barnes, J. W., & O'Brien, D. P. 2002, *ApJ*, 575, 1087

Barnes, R., & Quinn, T. 2004, *ApJ*, 611, 1

Bennett, D. P., Batista, V., Bond, I. A., et al. 2014, *ApJ*, 785, 155

Bennett, D. P., & Rhie, S. H. 1996, *ApJ*, 472, 660

Bennett, D. P., & Rhie, S. H. 2002, *ApJ*, 574, 985

Bills, B. G. 1990, *JGR*, 95, 14137

Bills, B. G. 2005, *Icar*, 175, 233

Canup, R. M., & Ward, W. R. 2009, in Europa, ed. R. T. Pappalardo et al. (Tucson, AZ: Univ. Arizona Press), 59

Carrera, D., Davies, M. B., & Johansen, A. 2016, *MNRAS*, 463, 3226

Chambers, J. E. 1999, *MNRAS*, 304, 793

Chambers, J. E., Wetherill, G. W., & Boss, A. P. 1996, *Icar*, 119, 261

Chatterjee, S., Ford, E. B., Matsumura, S., et al. 2008, *ApJ*, 686, 580

Cumming, A., Butler, R. P., Marcy, G. W., et al. 2008, *PASP*, 120, 531

Deienno, R., Yokoyama, T., Nogueira, E. C., et al. 2011, *A&A*, 536, A57

Dobos, V., Heller, R., & Turner, E. L. 2017, *A&A*, 601, 91

Domingos, R. C., Winter, O. C., & Yokoyama, T. 2006, *MNRAS*, 373, 1227

Donnison, J. R. 2010, *MNRAS*, 406, 1918

Ford, E. B., & Rasio, F. A. 2008, *ApJ*, 686, 621

Forgan, D., & Kipping, D. 2013, *MNRAS*, 432, 2994

Frouard, J., & Yokoyama, T. 2013, *CeMDA*, 115, 59

Gong, Y.-X., Zhou, J.-L., Xie, J.-W., et al. 2013, *ApJL*, 769, L14

Guillochon, J., Ramirez-Ruiz, E., & Lin, D. 2011, *ApJ*, 732, 74

Han, C. 2008, *ApJ*, 684, 684

Han, C., & Han, W. 2002, *ApJ*, 580, 490

Heller, R. 2012, *A&A*, 545, L8

Heller, R. 2014, *ApJ*, 787, 14

Heller, R., & Barnes, R. 2013, *AsBio*, 13, 18

Heller, R., & Barnes, R. 2014, *IJASB*, 14, 335

Heller, R., Williams, D., Kipping, D., et al. 2014, *AsBio*, 14, 798

Heller, R., & Zuluaga, J. I. 2013, *ApJL*, 776, L33

Hilton, J. L. 1991, *AJ*, 102, 1510

Hinkel, N. R., & Kane, S. R. 2013, *ApJ*, 774, 27

Holman, M. J., & Murray, N. W. 2005, *Sci*, 307, 1288

Holman, M. J., & Wiegert, P. A. 1999, *AJ*, 117, 621

Hong, Y.-C., Tiscareno, M. S., Nicholson, P. D., et al. 2015, *MNRAS*, 449, 828

Jeans, J. H. 1919, *MNRAS*, 79, 408

Jewitt, D., & Haghighipour, N. 2007, *ARA&A*, 45, 261

Johnson, J. A., Butler, R. P., Marcy, G. W., et al. 2007, *ApJ*, 670, 833

Jurić, M., & Tremaine, S. 2008, *ApJ*, 686, 603

Kaltenegger, L. 2010, *ApJL*, 712, L125

Kane, S. R. 2017, *ApJL*, 839, L19

Kasting, J. F., Whitmire, D. P., & Reynolds, R. T. 1993, *Icar*, 101, 108

Kinoshita, H., & Nakai, H. 1991, *CeMDA*, 52, 293

Kipping, D. M. 2009, *MNRAS*, 392, 181

Kipping, D. M., Bakos, G. Á, Buchhave, L., et al. 2012, *ApJ*, 750, 115

Kipping, D. M., Fossey, S. J., & Campanella, G. 2009, *MNRAS*, 400, 398

Kipping, D. M., Schmitt, A. R., Huang, X., et al. 2015, *ApJ*, 813, 14

Kouwenhoven, M. B. N., Goodwin, S. P., Parker, R. J., et al. 2010, *MNRAS*, 404, 1835

Lin, D. N. C., & Ida, S. 1997, *ApJ*, 477, 781

Lovis, C., & Mayor, M. 2007, *A&A*, 472, 657

Marzari, F., & Weidenschilling, S. J. 2002, *Icar*, 156, 570

Matsumura, S., Ida, S., & Nagasawa, M. 2013, *ApJ*, 767, 129

Mayor, M., Marmier, M., Lovis, C., et al. 2011, *arXiv:1109.2497*

Morbidelli, A. 2005, *arXiv:astro-ph/0512256*

Mróz, P., Udalski, A., Skowron, J., et al. 2017, *Natur*, 548, 183

Namouni, F. 2010, *ApJL*, 719, 145

Nesvorný, D., Vokrouhlický, D., & Morbidelli, A. 2007, *AJ*, 133, 1962

Nicholson, P. D., Cuk, M., Sheppard, S. S., et al. 2008, in *The Solar System Beyond Neptune*, ed. M. A. Barucci et al. (Tucson, AZ: Univ. Arizona Press), 411

Payne, M. J., Deck, K. M., Holman, M. J., et al. 2013, *ApJL*, 775, L44

Peters, M. A., & Turner, E. L. 2013, *ApJ*, 769, 98

Quinn, T. R., Tremaine, S., & Duncan, M. 1991, *AJ*, 101, 2287

Rasio, F. A., & Ford, E. B. 1996, *Sci*, 274, 954

Raymond, S. N., Armitage, P., & Gorelick, N. 2010, *ApJ*, 711, 772

Raymond, S. N., Armitage, P., Moro-Martín, A., et al. 2011, *A&A*, 530, A62

Raymond, S. N., Armitage, P. J., & Gorelick, N. 2009, *ApJL*, 699, L88

Raymond, S. N., Armitage, P. J., Moro-Martín, A., et al. 2012, *A&A*, 541, A11

Raymond, S. N., & Izidoro, A. 2017, *Icar*, 297, 134

Reynolds, R. T., McKay, C. P., & Kasting, J. F. 1987, *AdSpR*, 7, 125

Rowan, D., Meschiari, S., Laughlin, G., et al. 2016, *ApJ*, 817, 104

Sartoretti, P., & Schneider, J. 1999, *A&AS*, 134, 553

Sasaki, T., Barnes, J. W., & O'Brien, D. P. 2012, *ApJ*, 754, 51

Scharf, C. A. 2006, *ApJ*, 648, 1196

Simon, A., Szatmáry, K., & Szabó, G. M. 2007, *A&A*, 470, 727

Spalding, C., Batygin, K., & Adams, F. C. 2016, *ApJ*, 817, 18

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, *arXiv:1503.03757*

Sumi, T., Kamiya, K., Bennett, D. P., et al. 2011, *Natur*, 473, 349

Teachey, A., Kipping, D. M., & Schmitt, A. R. 2017, *arXiv:1707.08563*

Tinney, C. G., Wittenmyer, R. A., Butler, R. P., et al. 2011, *ApJ*, 732, 31

Tremaine, S. 1993, in *ASP Conf. Ser.* 36, Planets Around Pulsars, ed. J. A. Phillips, S. E. Thorsett, & S. E. Kulkarni (San Francisco, CA: ASP), 335

Tremaine, S., Touma, J., & Namouni, F. 2009, *AJ*, 137, 3706

Veras, D., & Armitage, P. J. 2005, *ApJL*, 620, L111

Veras, D., & Raymond, S. N. 2012, *MNRAS*, 421, L117

Ward, W. 1973, *Sci*, 181, 260

Ward, W. 1974, *JGR*, 79, 3375

Weidenschilling, S. J., & Marzari, F. 1996, *Natur*, 384, 619

Williams, D. M., Kasting, J. F., & Wade, R. A. 1997, *Natur*, 385, 234

Winn, N. J., & Frabrycky, C. D. 2015, *ARA&A*, 53, 409

Wittenmyer, R. A., Butler, R. P., Tinney, C. G., et al. 2016, *ApJ*, 819, 28