COULD JUPITER OR SATURN HAVE EJECTED A FIFTH GIANT PLANET?

We consider a reduced-body subset of the solar system including the Sun, a satellite-hosting gas giant planet, and an IG. The gas giant of mass M_p is initialized on an osculating, circular orbit with semimajor axis a_p. By approximating the gas giant's orbit as circular, we remove any dependence on its orbital phase at the epoch of encounter; t_enc. Because the mass of the ejected IG is not well constrained, we select a fiducial value approximately equal to the mass of Uranus (${M}_{I}=5\times {10}^{-5}{M}_{\odot }$), as was used by NM12. A system of Keplerian satellites is placed in orbit around the gas giant (see Section 3.2 for a detailed description of satellite orbits).

To limit the number of encounter parameters and the computational cost of our survey, we assume a coplanar geometry (we present a more detailed discussion of the effect of inclined encounters in Section 7.1). Thus, at its closest approach the IG's velocity is perpendicular to its separation vector from the gas giant and the encounter is fully specified at this time (t_enc) by the impact parameter b, the relative planet velocity v_rel, and the phase angle θ (see Figure 1).

Zoom In Zoom Out Reset image size

We linearly step the impact parameter outward from ${b}_{\mathrm{min}}=0.02\;{\rm{AU}}$ until encounters no longer lead to ejections. Deienno et al. (2014) found that close encounters with b < b_min excessively excite the Galilean satellites; our estimates in Section 6 for the likelihood of retaining the observed satellite orbits following an ejection therefore represent a conservative overestimate due to our comparatively gentle ejections. However, we note that values smaller than b_min are less likely due to the reduced encounter cross-section at small radii. We estimate the size of this effect by including a reduced simulation sample with b < b_min and find that our final results (Section 6) changed (fractionally) by $\lesssim 12$%.

To determine which encounter parameters are capable of ejecting the IG, we uniformly sample 20,000 parameter combinations within $b\in [{b}_{\mathrm{min}},0.1$ AU], ${v}_{\mathrm{rel}}\in [1,5]{v}_{\mathrm{esc}},$ and $\theta \in [0,2\pi ),$ where ${v}_{\mathrm{esc}}=\sqrt{2{{GM}}_{p}/b}.$ We then remove any unphysical encounters in which the IG is unbound from the Sun prior to the encounter. In the cases of Jupiter and Saturn we find N_sim = 278 and N_sim = 274 as valid parameter combinations, respectively.

To simulate an encounter with the aforementioned parameters from Figure 1, we integrate backward in time until the absolute separation of the planets $| {r}_{I}-{r}_{p}| \geqslant 2\;{\rm{AU}}.$ At this separation a satellite with a Callisto-like or Iapetus-like orbit feels a force from its host that is >10⁴ times greater than that felt from the IG. The positions and velocities of all three massive bodies are then used as initial conditions for the forward simulations that include the satellites. These simulations are integrated forward in time toward the encounter at t = t_enc and are halted at t = 2t_enc. After 2t_enc the influence of the IG on the satellites is again negligible and satellite orbits are no longer perturbed by the IG's influence. In order for close encounters with $b\lt 0.1\;{\rm{AU}}$ to be possible, numerous "soft" encounters between the gas giants and IG, prior to t_enc, were needed in order to build-up the IG's eccentricity. These encounters will supply a small perturbation to regular satellite orbits that, on average, increase satellite eccentricity over time (see Nesvorný et al. 2014, Figure 4). However, the effect of encounters on satellite orbits is strongly dependent on b, which is much larger than 0.1 AU for "soft" encounters. We are therefore only concerned with the strongest (final) encounter that leads to the ejection of the IG as its effect dominates the final satellite orbits.

Finally, we do not include the perturbations from the satellite host's oblateness. This should be a good approximation since the encounter timescale with the IG is much shorter than the precession timescales of the satellites due to the non-spherical shape of its host planet.

In each encounter simulation, we include a ring of Callisto or Iapetus analog satellites around the gas giant planet. The satellite ring consists of N = 100 non-interacting test particles with azimuthal positions randomly sampled from a uniform distribution between 0 and $2\pi$ in order to remove any azimuthal dependence at t_enc. Modeling the satellite system as an ensemble of test particles ensures that the orbital evolution of the satellites is governed solely by gravitational interactions with the host planet and perturbations from the IG and the Sun.

The semimajor axes of the satellites a_s are initialized to ±1% of the current orbital radius of Callisto or Iapetus. This fractional deviation is chosen to be on the order of the observed satellite eccentricities e_s (see Table 1). Because changes in a_s and e_s are related through changes in the satellites' angular momentum, it is unlikely that larger shifts in a_s would be reconcilable with the observed e_s, assuming that each satellite formed out of a circumplanetary disk on a circular orbit (Canup & Ward 2002; Mosqueira & Estrada 2003). Satellites are initialized on nearly coplanar and circular orbits ($i\leqslant 0\buildrel{\circ}\over{.} 1,e\leqslant {10}^{-5}$), as expected from circumplanetary formation scenarios. Satellite particles are not allocated physical sizes as we do not account for particle collisions in our simulations.

We performed our simulations using the REBOUND N-body numerical code (Rein & Liu 2012). We employ the IAS15 integrator (Integrator and Adaptive Step-size control, 15th order; Rein & Spiegel 2015), whose adaptive timestepping ensures optimal resolution of the short satellite orbital periods and rapid encounter timescales.⁵

Figure 2 summarizes the properties of the planetary encounters that result in the ejection of the IG planet. We find N_sim = 278 such encounters. The greatest impact parameter we find capable of ejecting the IG is $\approx 0.05\;{\rm{AU}}.$

Zoom In Zoom Out Reset image size

All successful ejections involve the closest approach in the lower hemisphere of the xy-plane in Jupiter's reference frame ($180^\circ \lt \theta ({t}_{\mathrm{enc}})\lt 360^\circ ;$ see Figure 1). As familiar from spacecraft gravity assists, an IG trailing Jupiter at t_enc will receive a positive velocity kick via their interaction. At these phase angles, the relative velocity vector is rotated by the encounter such that the IG's inertial speed accelerates. Depending on b and v_rel, the encounter can potentially boost the IG to escape velocity from the solar system. From simple vector diagrams, the maximum increase to the IG velocity can be shown to occur when $\theta ({t}_{\mathrm{enc}})=270^\circ .$ Figure 2 highlights this as the distribution of $\theta ({t}_{\mathrm{enc}})$ peaks at ∼270°, where the encounter geometry is most conducive to ejecting the IG.

Low v_rel trajectories more often lead to ejection since Jupiter can more effectively deflect the IG's path. We note that there are many encounters in which the velocity of the IG with respect to Jupiter at t_enc is greater than the escape speed from the solar system at Jupiter's distance from the Sun. This is due to the IG being accelerated upon approach to the encounter, which occurs deep in Jupiter's gravitational well. Such cases consist of the IG becoming unbound from the solar system prior to the realization of the impact parameter at t_enc and being thus accelerated to super-escape speeds even before its closest approach to Jupiter.

The 2D (two-dimensional) histogram in Figure 2 depicts the correlation between the relative velocity of the planets at the epoch of encounter and the phase angle in the Jovian-centric reference frame. Encounters that occur with the IG at a smaller heliocentric distance than Jupiter ($\theta ({t}_{\mathrm{enc}})\lt 270^\circ$) require larger relative velocities in order to eject the IG. As the IG moves outward (increasing θ), the relative velocities necessary to eject the IG decrease because the pull applied by Jupiter following the gravity assist is less efficient at decelerating the IG to sub-escape speeds given the IG's trajectory.

Close planetary encounters between Jupiter and the IG can result in significant perturbations to the orbits of the in situ Jovian satellites. The time evolution of a_s and e_s for five sample Jovian satellites throughout one encounter simulation demonstrates this fact (Figure 3).

Zoom In Zoom Out Reset image size

The five sample satellites that are initialized on nearly circular orbits get kicked to moderately eccentric orbits by a planetary encounter with $b\approx 0.042\;{\rm{AU}}$ and ${v}_{\mathrm{rel}}\approx 13.57$ km s⁻¹. Hence, the characteristic encounter timescale is $b/{v}_{\mathrm{rel}}\sim 5.3$ days, or about one-third of Callisto's orbital period. The orbital perturbation therefore acts roughly as an impulse to the satellites making their azimuthal position at t_enc, a parameter of importance. Indeed, in Figure 3 the satellites whose azimuthal positions are increasingly distant from $\theta ({t}_{\mathrm{enc}})$ are on average less perturbed than satellites that are close to $\theta ({t}_{\mathrm{enc}}).$ In addition, the position of the satellite relative to $\theta ({t}_{\mathrm{enc}})$ determines the direction of the satellite's radial shift. For example, satellites immediately trailing the IG at t_enc will be pulled forward in their orbits (e.g., red curve) increasing a_s, whereas satellites immediately upstream of the IG at t_enc get pulled downward, thus decreasing a_s (e.g., black curve).

The average final orbital elements of all Jovian satellites are depicted in Figure 4 and show close agreement with the results from Deienno et al. (2014) despite the broader investigation of the close encounter parameter space. Changes in a_s and e_s are coupled through the satellite's orbital angular momentum. This effect is highlighted in Figure 4, as satellites that are radially transferred to either a significantly larger or smaller semimajor axis are those which exhibit the largest deviation from a circular orbit.

Zoom In Zoom Out Reset image size

For each satellite we define a reconcilable orbit to be when the satellite's final average eccentricity is $\leqslant 2{e}_{\mathrm{Callisto}}$ (recall e_Callisto ≈ 0.007). The factor of 2 in our definition comes from the subsequent eccentricity evolution due to tidal damping of Callisto (Deienno et al. 2014, see Figure 7) and allows for Callisto to be excited beyond the present e_Callisto at t_enc and consequently settle into its current orbital eccentricity in the 4 Gyr following the solar system's instability phase.

For the remainder of the paper we refer to the event of a simulated satellite's final orbit being reconcilable with Callisto as RS, for "reconcilable satellite." The boundary dividing RS from non-RS satellites is depicted in Figure 4 as a dashed horizontal line.

Final average e_s values are never ≥0.4, implying that no Jovian satellite becomes unbound from Jupiter following the ejection of the IG. That is, the vast majority of planetary encounters that are capable of ejecting an IG from the solar system are not sufficiently violent to strip Jupiter's regular satellites. This favors the possibility of Jupiter being able to eject an IG while retaining a regular satellite whose orbit is Callisto-like. We estimate the likelihood in Section 6.

While the phase angle θ(t_enc) is an important parameter in determining whether the IG is ejected, it has little effect on the fraction of perturbed satellites, due to their uniformly distributed azimuthal positions around Jupiter. Therefore, for each b and v_rel we can marginalize over θ to show in Figure 5 the fraction of satellites that are reconcilable with the orbit of Callisto after each encounter. Interpolating over the fraction of reconcilable Jovian satellites in (b, v_rel) space, the resulting high-resolution contours at 10% and 50% are fitted with cubic functions and overplotted as solid and dashed curves, respectively, to aid in the visualization of the different regions of the parameter space.

Zoom In Zoom Out Reset image size

It is clear that as encounters become closer, the perturbation to e_s increases and the fraction of RSs shrinks. Similarly, as the duration of encounters becomes longer (smaller v_rel), the timescale over which the perturbation is applied grows, thus increasing the deviation of satellite orbits from circular. Hence forming and maintaining a Jovian satellite with a Callisto-like orbit favors encounters that are wide and fast.

We suggest that researchers simulating solar system formation scenarios can use Figure 5 to estimate whether or not a given Jupiter/IG encounter is consistent with Jupiter's current Galilean satellite orbits. We note that our requirement that the IG gets ejected only determines what regions of the plot are populated, hence the contours are useful guides regardless of whether or not the IG survives the encounter (and for any θ). We caution, however, that these will only be approximate due to our assumption of coplanarity. For a given b and v_rel, inclined encounters will increase the fraction of satellites whose eccentricities are reconcilable with current orbits. However, they will also tend to overly excite the satellite inclinations. We discuss this further in Section 7.1. For particular simulations that lie on the boundary of plausibility in Figure 5, one could re-simulate the close approaches, including inclined encounters, following our setup to more accurately quantify the effect.

The properties of planetary encounters between Saturn and an IG are summarized in Figure 6. We find N_sim = 274 simulations that result in the IG being ejected by Saturn. The largest impact parameter that is still capable of ejecting the IG is nearly identical to the Jupiter/IG case; ≈0.05 AU. We perform an additional sampling of encounter parameters with b > 0.05 AU at an increased resolution to confirm that the largest impact parameter capable of ejecting the IG is nearly the same in both the Jupiter and Saturn cases rather than being a statistical anomaly due to the finite sampling of encounter parameters. We find that this limit is real and not an artifact of our sampling procedure. The physical interpretation of this similarity is beyond the scope of this paper.

Zoom In Zoom Out Reset image size

Similar trends to those shown in Figure 2 are observed. This should be expected because, modulo the variations in the physical parameters M_p and a_p, the Saturn/IG encounters investigated are fundamentally equivalent to those described in Section 4.1. With the exception of the decreased magnitude of v_rel by a factor of $\sqrt{{M}_{\mathrm{Sat}}/{M}_{\mathrm{Jup}}}$, on average, the set of successful scattered-by-Saturn simulations are statistically identical to those in the scattered-by-Jupiter simulations; i.e., the histograms in Figures 2 and 6 exhibit the same behavior.

Similar to the case of Jupiter ejecting the IG from the solar system, the Kronian satellite orbits will be perturbed as a result of the encounter between Saturn and the IG. As a result of the wide separation of Iapetus (∼61 R_S) and the similar nature of the encounters between an IG and either Jupiter or Saturn (Figures 2 and 6), it is expected that such orbital perturbations will be more destructive for simulated satellite orbits than for those previously explored in Section 4. This notion is supported in Figures 7–9 when compared to their counterparts in the Jupiter/IG case (Figures 3–5), as the fraction of reconcilable Kronian satellites is systematically lower than in the Jupiter case.

Zoom In Zoom Out Reset image size

The time evolutions of a_s and e_s for five example satellites during a sample planetary encounter are shown in Figure 7. This example is a particularly violent encounter with $b\approx 0.033\;{\rm{AU}}$ and v_rel ≈ 5.87 km s⁻¹, as no satellite in this simulation has a resulting orbit that is reconcilable with Iapetus. One satellite, whose azimuthal position is only 46 from $\theta ({t}_{\mathrm{enc}}),$ is ejected from Saturn by the encounter. The encounter timescale ($b/{v}_{\mathrm{rel}}\sim 9.6$ days) is much less than the orbital period of Iapetus, making the effect of the encounter behave approximately as an impulse.

The average final orbital elements of all Kronian satellites after each unique encounter are depicted in Figure 8. The black contours in Figure 8 are representative of the final orbital elements of the Jovian satellites from Figure 4 to aid in visual comparison between the final orbits of Jovian satellites to Kronian satellites following similar close planetary encounters. This shows that the wide separation of the simulated Kronian satellites makes them more susceptible to excessive orbital perturbations.

Zoom In Zoom Out Reset image size

The event RS in the Saturn/IG case is defined similar to the Jupiter/IG case and is achieved by Kronian satellites whose final average is ${e}_{{\rm{s}}}\leqslant {e}_{\mathrm{Iapetus}}.$ As noted earlier, the eccentricity evolution of Iapetus since the solar system's instability phase is negligible (Castillo-Rogez et al. 2007), thus we assume that resulting orbits from our simulations will go largely unchanged to the present day. The boundary dividing RSs from non-RSs is depicted in Figure 8 as a dashed horizontal line.

As noted in Section 4.2, no Jovian satellite becomes more eccentric than 0.4, whereas a sizable fraction of Kronian satellites become equivalently or excessively excited including a subset of Kronian satellites whose final average e_s is approximately unity. Furthermore, we find that ∼6% of sampled Kronian satellites are sufficiently excited to final e_s > 1. The orbital elements of these ejected satellites are not included in Figure 8. It is clear that the perturbations to the Kronian satellites resulting from close encounters are, on average, much stronger and capable of stripping Iapetus-like satellites from the Kronian system. The large fraction of satellites perturbed beyond the orbit of Iapetus makes it statistically difficult for Saturn to have ejected an IG while retaining an Iapetus-like regular satellite. We formally estimate the likelihood of reconciling Iapetus' orbit in Section 6.

The effect of encounter properties on the resulting Kronian satellite orbits is shown in Figure 9. For each encounter with a given impact parameter and relative planet velocity, the fraction of Kronian satellites whose final orbit is reconcilable with the current orbit of Iapetus is shown. Approximate contours of 10% and 50% fractions are overplotted. The 10% contour is computed identically to the contours in Figure 5 (i.e., cubic interpolation). However, the 50% contour lacks a sufficient number of points to perform a robust cubic interpolation. Therefore we opt for a linear fit in its place. It is clear that the region of the (b, v_rel) parameter space in which ≥50% of the Kronian satellites are reconcilable with Iapetus' orbit is very small compared to the Jovian satellites with only five unique ejections sampled there.

Zoom In Zoom Out Reset image size

The effects of b and v_rel on resulting satellite orbits are nearly identical to those shown in Figure 5. One interesting difference that is unique to Saturn/IG encounters is that regardless of how close the encounter is, if the encounter is sufficiently long (v_rel ≲ 8 km s⁻¹), then <10% of satellites will be reconcilable with Iapetus. This is in contrast to the Jovian satellite case where even the slowest encounters could retain a high fraction of Callisto-like satellites if the encounter's impact parameter is large. Another important feature to note is that no Saturn/IG encounter resulting in the ejection of the latter is guaranteed to preserve an Iapetus-like satellite. That is, even the least violent encounters leading to ejection are only capable of preserving a maximum of ∼70% of the in situ Iapetus-like satellites.

Figure 9 contains information on the likelihood that Iapetus survives an IGE by Saturn with a particular b(t_enc) and ${v}_{\mathrm{rel}}({t}_{\mathrm{enc}})$ when the mutual inclination of the encounter is near zero. Similar to the Jupiter/IG encounter case, we suggest that researchers simulating solar system formation scenarios can use Figure 9 to estimate whether or not a given Saturn/IG encounter is consistent with Iapetus' current orbit in the limit of uninclined planetary encounters.

The resemblance of the final average orbits of the Jovian and Kronian satellites (Figures 4 and 8) to the current orbits of Callisto and Iapetus can in principle be used to compute the likelihood of a fifth giant planet getting ejected by either of the gas giant planets in the early solar system. If such an event were to have occurred, it must be consistent with the satellite orbits presently observed. The likelihood of an IG getting ejected by either gas giant requires the likelihood of the current orbits of Callisto or Iapetus being reconciled by simulated satellites after an IGE. A successful event in which the resulting orbit of a simulated Jovian (Kronian) satellite is reconcilable with Callisto's (Iapetus') current orbit is referred to as RS for "reconcilable satellite."

Given that an IG gets ejected in all simulations (event IGE; "IG ejected"), for the $i\mathrm{th}$ satellite in the $j\mathrm{th}$ simulation, we record whether or not the event RS is achieved using

$$\begin{eqnarray}{p}_{i,j}({\rm{RS}}| {\rm{IGE}})=\left\{\begin{array}{ll}1 & \quad \mathrm{if}\ {\rm{RS}}\\ 0 & \quad \mathrm{if}\ \mathrm{not}\ {\rm{RS}}.\end{array}\right.\end{eqnarray} \tag{ 1 }$$

However, not all planet orbits in our ejection simulations are statistically relevant. Our adopted methodology for determining which encounter parameters result in an ejection is heavily biased toward initially high-eccentricity orbits of the IG (${e}_{I}\lesssim 1$), which are not long-term stable and therefore are uncommon in nature. This bias naturally arises because it is easier to kick the IG to e_I > 1 if e_I is initially very close to unity. Therefore in each simulation j, each satellite's likelihood of an RS must be weighted by the likelihood of the IG having an initially bound orbit with e_I,j, where e_I,j is the IG's initial eccentricity in the $j\mathrm{th}$ simulation.

The corresponding weighting function $W({e}_{I,j})$ is modeled by the distribution of planet eccentricities observed in a statistically significant number of exoplanetary systems; the eccentricity distribution of solar system bodies is insufficient, as it can only be derived in the limit of small number statistics. These data are recovered from the www.exoplanets.org database (Han et al. 2014) and only include RV exoplanet detections with Doppler variation semi-amplitude $K/{\sigma }_{K}\gt 5$ (neglect low signal-to-noise observations). We find that our results do not sensitively depend on whether we focus on RV detections or the full catalog of exoplanets with orbital solutions. The distribution of planet eccentricities is shown in Figure 10. Due to variations in empirically derived eccentricity distributions we consider three proposed analytical forms of the distribution. Namely a Beta probability density function (PDF; Kipping 2013), a Rayleigh PDF plus decaying exponential (Jurić & Tremaine 2008; Steffen et al. 2010), and the model from Shen & Turner (2008). We use a Levenberg–Marquardt least squares algorithm to fit for each PDF's unique parameters. A summary of the adopted distributions including fitted parameters is given in Table 2 and each fitted weighting function is overplotted in Figure 10. Differences in $W({e}_{I,j})$ from adopting three unique PDFs result in ∼6% variance among computed likelihoods.

Zoom In Zoom Out Reset image size

Table 2.  Analytical Fits to the Planet Eccentricity Distribution

PDF Name	PDF	Parameter 1	Parameter 2	Parameter 3
Beta	${W}_{\beta }{(x;a,b)=\frac{1}{B(a,b)}{x}^{a-1}(1-x)}^{b-1}$	0.786 ± 0.055	2.764 ± 0.167	⋯
Rayleigh + exponential	${W}_{{\rm{R}}+\mathrm{exp}}(x;\alpha ,\lambda ,\sigma )=\frac{x(1-\alpha )}{{\sigma }^{2}}\mathrm{exp}(\frac{-{x}^{2}}{2{\sigma }^{2}})+\alpha \lambda \mathrm{exp}(-\lambda x)$	0.782 ± 0.301	5.131 ± 2.342	0.266 ± 0.061
ST08	${W}_{\mathrm{ST}08}(x;k,a)=\frac{1}{k}\left(\frac{1}{{(1+x)}^{a}}-\frac{x}{{2}^{a}}\right)$	0.305 ± 0.012	3.413 ± 0.309	⋯

Note. Parameter columns are written in the order in which they appear in the functional form $W(x;...),$ shown in the PDF column. We write the independent variable of planet eccentricity as x to make a distinction between eccentricity and Euler's number $e\approx 2.7182818.$

Download table as:  ASCIITypeset image

While we take the known exoplanet orbital eccentricities as our nominal distribution, this sample is obviously biased by the requirement that systems be stable. If early on the solar system's IGs moved on significantly more elliptical paths than the known exoplanets, it would be easier to eject planets while keeping Jupiter and Saturn's satellites on nearly circular orbits, raising the likelihood of reconciling satellite orbits. However, this scenario would require a substantial amount of damping to subsequently recircularize the surviving planets orbits, and if eccentricities reached such high values, one might expect to also lose Uranus and Neptune. These considerations may be interesting directions for future work.

One must additionally weigh ${p}_{i,j}$ by the simulation's impact parameter b_j, since wider encounters should occur more frequently. Given our coplanar setup, the differential interaction cross-section for a given impact parameter scales linearly with b_j rather than with b_j², as is true in the full 3D (three-dimensional) case. We find that our results do not sensitively depend on this distinction. We therefore assume $W^{\prime} ({b}_{j})\propto {b}_{j}.$

The distribution of ${e}_{I,j}$ is not explicitly prescribed and is instead determined from the encounter parameters (b, v_rel, θ). The resulting distribution of e_I,j is not sampled uniformly, unlike the distribution in b or θ, implying that the step size ${\rm{\Delta }}{e}_{I}$ when integrating over IG eccentricities is not constant throughout the domain. Hence, simulations sample IG eccentricity bins of various widths, which are taken into account when computing the likelihood function by calculating the width of IG eccentricity bins among the N_sim simulations. In this way, the subspace of IG eccentricities, which is oversampled, gets averaged over. Similar factors ${\rm{\Delta }}{b}_{j}$ and ${\rm{\Delta }}{\theta }_{j}$ are included but are constants because the parameters are sampled uniformly, and therefore do not affect the resulting likelihood.

Combining the aforementioned effects into the likelihood of obtaining an RS orbit given an IGE, we write

$$\begin{eqnarray}{\rm{P}}({\rm{RS}}| {\rm{IGE}}) & = & \displaystyle \frac{1}{\mu }\displaystyle \sum _{i=1}^{N}\displaystyle \sum _{j=1}^{{N}_{\mathrm{sim}}}\;{p}_{i,j}({\rm{RS}}| {\rm{IGE}})\\ & & \times W({e}_{I,j})W^{\prime} ({b}_{j}){\rm{\Delta }}{e}_{I,j}{\rm{\Delta }}{b}_{j}{\rm{\Delta }}{\theta }_{j},\end{eqnarray} \tag{ 2 }$$

where

$$\begin{eqnarray}&&\mu =N\displaystyle \sum _{j=1}^{{N}_{\mathrm{sim}}}W({e}_{I,j})W^{\prime} ({b}_{j}){\rm{\Delta }}{e}_{I,j}{\rm{\Delta }}{b}_{j}{\rm{\Delta }}{\theta }_{j},\end{eqnarray} \tag{ 3 }$$

is the normalization factor of the weighted mean and we are careful to account for the fact that our uniform sampling of (b, v_rel, θ) leads to a non-uniform distribution of ${e}_{I,j}.$

From our ${N}_{\mathrm{sim}}=278$ ejection simulations by Jupiter, each with N satellites, we calculate that the likelihood of Jovian satellite orbits remaining consistent with the observed orbit of Callisto, is ∼42%. This value is the median of the results we obtained from adopting the three distinct eccentricity weighting functions discussed in Section 6.1 (see Table 2). The median absolute deviation among the likelihoods is small; ≲1%.

By contrast, in our N_sim = 274 ejection simulations by Saturn, the likelihood of Kronian satellite orbits remaining consistent with Iapetus' observed orbit is ∼1%, more than an order of magnitude less likely than for Callisto around Jupiter. The main reason for this wide likelihood disparity is that a given encounter will perturb Iapetus more strongly than Callisto because the former is less tightly bound to Saturn than the latter is to Jupiter. In this case, the median absolute deviation among the three eccentricity weighting functions is ≲0.6%, which is comparable to the likelihood itself.

To limit the computational cost of our study, we have restricted our analysis to planetary encounters with no mutual inclination. However, if planetary eccentricities are excited enough to permit orbit-crossing and ejections, one might expect comparably large orbital inclinations. During an inclined encounter, some of the applied torque goes into realigning Callisto's (Iapetus') orbital plane so that, on average, the perturbations to a_s and e_s are reduced as some energy goes into increasing i_s.

A preliminary analysis of inclined encounters with mutual inclination i(t_enc) = 5° effectively revealed no change to the final satellite orbital elements compared to uninclined, but otherwise equivalent encounters. In a more heavily inclined test case with $i({t}_{\mathrm{enc}})=45^\circ ,$ we found that, on average, the final e_s were smaller by ≳0.01 than in the uninclined case.

Therefore, by limiting our investigation to coplanar encounters, we are sampling the largest possible perturbations to e_s without affecting i_s. The introduction of inclined encounters would thus raise the likelihoods quoted in Section 6.2. However, such a 3D case would additionally excite the satellite inclinations. As mentioned by Deienno et al. (2014), these inclinations may provide more rigorous constraints on planetary encounters because even in cases where eccentricity damping is important, the tidal evolution of the inclinations is effectively null. A generalization of this study would therefore consider the combined constraint from the satellites' orbital eccentricities and inclinations. However, the uncertainties in the likelihoods we compute in Section 6.2 are dominated by the large uncertainties in the planetary orbital eccentricities (and inclinations) early in the solar system. We therefore believe that our simplified analysis considering only the orbital eccentricities captures the correct likelihoods at the approximate level allowable by our current state of knowledge.

As a rough check, in a fully 3D case, the differential cross-section for encounters of a given impact parameter would scale as b_j² rather than as b_j. Using the full data set from our coplanar simulations, but adopting this new scaling for W'(b_j) in Equations (2) and (3), we find that our results do not discernibly change. Specifically, in the Jupiter case our results change from ∼42% to 54% and from ∼1.1% to 1.3% in the Saturn case. This is certainly within the errors of our uncertain knowledge of the planetary orbits' initial eccentricities and inclinations and therefore does not appreciably change the interpretation of our results.

We conclude that Jupiter could plausibly have ejected an IG from the solar system. Nearly half of the hypothetical set of ejections that were modeled in Section 6.2 keep Callisto on an orbit that is reconcilable with the one we observe today. Put another way, we conclude that Callisto's orbit cannot meaningfully constrain whether or not Jupiter ejected an additional IG in the early solar system. Nevertheless, this is an important test for the fifth giant planet hypothesis to pass, thus providing more stringent evidence for its plausibility.

Interpretation of our results in the Saturn case is more subtle. We showed in Section 6.2 that an initially circular Iapetus orbit gets overly excited by a single ejection event $\sim 99\%$ of the time. This suggests that Saturn is not capable of ejecting an IG mass planet from the solar system.

But did Iapetus originally move on a circular path, as expected if it formed out of a circumplanetary disk? Starting with a circular orbit, all encounters act to raise the eccentricity. But if the moon's initial path were instead elliptical, some ejection geometries could act to lower the eccentricity, complicating the constraint that we nominally set. Perhaps one reason to doubt that Iapetus formed from a circum-Kronian disk is that it is substantially inclined (≈8°) to the local equilibrium plane that one would expect it to follow.

So what caused this aberrant inclination? Hamilton (2013) recently suggested that Iapetus could indeed have formed on a circular orbit from a disk in its equilibrium plane, and subsequent collisions between Saturn's inner moons could have instead tilted the equilibrium plane for the exterior moons by the requisite amount. In this case, our study implies that Saturn likely did not eject an IG from the solar system because Iapetus' relatively low-eccentricity orbit cannot be reconciled with IGEs by Saturn. On the other hand, Iapetus' inclination could be the signature of an ejection event itself. Nesvorný et al. (2014) recently studied such a scenario, trying to simultaneously match Iapetus' current orbital eccentricity and inclination using simulations of the specific early solar system dynamical instabilities from NM12. They show that some cases are capable of sufficiently exciting Iapetus' orbital inclination while maintaining a low orbital eccentricity even for encounters as close as those considered in this study but not necessarily leading to an ejection of the IG.

In the Kronian case, therefore, the interpretation depends critically on the formation mechanism for Iapetus, which is currently unknown. If Iapetus' inclination is the result of collisions between inner moons (Hamilton 2013), our results show that Saturn is unlikely to have ejected an IG from the early solar system. By contrast, if Iapetus' orbit is fully explainable through close planetary encounters, Nesvorný et al. (2014) showed that this requires many such close approaches that might be capable of ejecting the IG.

7.2.1. Single-encounter Assumption

7.2.2. Close Encounters without Ejection

One sought after quantity relating to this work is the probability of an IG getting ejected by either gas giant given that the orbit of Callisto or Iapetus can be reconciled; ${\rm{P}}({\rm{IGE}}| {\rm{RS}}).$ Using Bayes' theorem, this quantity could in principle be computed with knowledge of the likelihood functions calculated in this paper but also requires the independent probability of the IG getting ejected ($P({\rm{IGE}});$ see NM12) and the normalization factor by the probability of a satellite exhibiting the current orbit of one of the wide-separation satellites (P(RS)). Despite the abundance of work done on the latter (e.g., Canup & Ward 2002; Estrada & Mosqueira 2006; Ward & Canup 2010; Crida & Charnoz 2012; Hamilton 2013; Heller et al. 2015), precise constraints on P(RS) are difficult to compute.

Also, a robust calculation of the likelihood of ejecting an IG from the solar system would require the event to be consistent with a vast number of constraints imposed by solar system bodies, of which the orbital constraint RS imposed by Callisto or Iapetus is just one. Therefore, such a calculation is not practical. However, consistency checks of IGEs by a gas giant, such as those presented in this study, help to substantiate the proposed existence of a fifth giant planet. Based on our current understanding, which includes the main results of our study, there is little evidence demonstrating that an additional IG mass planet could not have existed in the early solar system.