IRREGULAR SATELLITE CAPTURE BY EXCHANGE REACTIONS

Gomes et al. (2005) proposed a new planetary-orbit evolution model to explain LHB (PC phase in our context). Because their numerical simulation contained a planetesimal disk of a limited number of objects, Nesvorný et al. (2007) improved resolution by cloning planetesimals at about 20 Myr before Jupiter and Saturn entered mutual 2/1 mean motion resonance. The starting data of Nesvorný et al. (2007) contained 6868 planetesimal particles, each of ∼1.6 Pluto mass. The advantage of this initial data modification was twofold: (i) the more massive particles in the original disk of Gomes et al. could produce an unrealistically granulated planet migration, and (ii) a variety of planetary migration scenarios can be studied by using different disk clones. It is important to verify the dependence of results on (ii). Our starting orbits of outer planets are those of Nesvorný et al. (2007) and correspond to the situation described by Gomes et al. about 20 Myr before PC: namely, they are near-circular and coplanar at 5.4, 8.4, 12.3, and 18.0 AU (terrestrial planets were not included in our work).

These initial data were propagated using a symplectic integrator SyMBA (Duncan et al. 1998). SyMBA is an efficient numerical integrator that can handle close encounters between massive bodies. While all particles are assumed to be massive in the simulation, only the planets gravitationally interact with each other and gravitationally influence planetesimal disk particles. Planetesimals gravitationally interact with planets, making them migrate, but the code neglects mutual gravitational effects between the planetesimals themselves.

Using 50 initial disk variants, Nesvorný et al. (2007) obtained 14 cases of plausible planetary evolution, each resulting in a solar system architecture approximately matching the current state. Here we selected seven of them and integrated the orbital evolution following the planet-crossing phase of the Nice model (Tsiganis et al. 2005; we used those that best matched the final architecture of the solar system). During this phase, planets no longer undergo close approaches with each other, but slowly migrate to attain their final positions. In each of the selected simulations we recorded parameters of close planetesimal flybys of the planets during a control interval of 100 Myr during this post-PC phase. We believe this time interval is representative enough to characterize such encounters and still be amenable to our computer power (in Section 4 we argue how our results could be extrapolated to later periods of the solar system evolution). In particular, we recorded data for every approach when a planetesimal entered the planetary Hill sphere. At the instant of the closest approach to the planet, we used the relative position and velocity of the planetesimal with respect to the planet to compute three parameters of the hyperbolic-flyby approximation: (i) the orbital pericenter q with respect to the planet, (ii) the velocity of the planetary approach at infinity V_∞, and (iii) the inclination of the hyperbolic plane with respect to the ecliptic. This latter parameter was found to be isotropic in space. For purposes of scaling our simulations to the real number of planetesimals in the disk, we also recorded the initial number of disk particles N_TP after planets stopped crossing their orbits and the total number of planetesimal encounters N_enc with each of the planets.

Table 1 provides N_TP and N_enc values for each of our seven runs (note N_TP is different in each of the runs because it gives the number of particles at the moment when planets stopped crossing their orbits). The class designations DE and MA in Table 1 denote the type of planetary evolution defined in Nesvorný et al. (2007): (i) DE stands for "direct emplacement" indicating that both Uranus and Neptune attained their near-final orbits nearly at the end of the mutual-orbit crossing phase, while (ii) MA stands for "Malhotra-class evolution," indicating the post-PC planetary system was still considerably compact and Uranus and Neptune underwent significant migration evolution, similar to that originally proposed by Malhotra (1995) and Hahn & Malhotra (1999). Note, however, that Malhotra (1995) started her planetary migration simulations with zero eccentricity and inclination orbits of both planets and planetesimals in the disk, while in our MA class the orbits have large eccentricity and inclination values after the planet crossing phase. Job no. 4 is somewhat in between these two classes and we designate it MA/DE. Following Tsiganis et al. (2005), we favor DE-type migrations, because MA-type migrations nearly always result in an unrealistic damping of Uranus's and Saturn's orbital eccentricity. For the sake of comparison, though, we keep two MA-type migrations among our selected cases (runs 2 and 7).

Figure 1 shows planetary-orbit evolution for two of our simulations (runs 1 and 2). The 100 Myr test interval over which we recorded planetesimal flybys near the migrating planets is denoted by the black arrows. The bottom panel, corresponding to MA-type run 2, indicates that Uranus's eccentricity became effectively zero, in contrast to the currently observed value.

Zoom In Zoom Out Reset image size

Figure 2 shows the evolution of the planetesimal disk mass in runs 1 and 2 from the previous figure, as well as the time-resolved number of planetesimal flybys near the four planets (note the initial disk mass before PC phase was ∼25 Earth masses). In both cases, the disk mass dropped by ∼65% during the post-PC interval of 100 Myr; this is because the majority of particles have been swept out of the solar system by encounters with Jupiter, or placed in outer solar system structures, such as the Oort cloud or the scattered disk, by encounters with Jupiter and other planets. As expected, the number of planetary encounters decreases in time roughly proportionally to the overall disk-mass decay. Some fine structures, such as a drop of Uranus encounters in run 1 at about 10 Myr, can be explained by a closer analysis of the simulations. In this case, the decrease of Uranus's eccentricity (Figure 1) decoupled this planet from the denser outer parts of the planetesimal disk.

Zoom In Zoom Out Reset image size

Figure 3 shows the distribution of the recorded values of q and V_∞ for all planets in two runs 1 and 2. While the overall results are broadly similar in the two cases, there are small differences such as systematically slightly lower encounter speeds V_∞ in the MA-type run (run 2). For instance, the distribution maximum of V_∞ in the case of Neptune occurs at ∼3 km s⁻¹ (run 1) and ∼2 km s⁻¹ (run 2). The maximum shifts to larger values of V_∞ for interior giant planets, up to ∼6 km s⁻¹ (run 1) and ∼5 km s⁻¹ (run 2) for Jupiter. These rather high velocities are due to the fact that the planetesimal disk is highly excited during the post-PC phase within the Nice model (see, e.g., Figure 3 in Nesvorný et al. 2007). Note that small relative velocities (with respect to the planets) that are generally favorable for irregular capture by exchange reactions (e.g., Agnor & Hamilton 2006 assumed encounter speeds ⩽0.5 km s⁻¹) occur for only a very small fraction of recorded flybys.

Zoom In Zoom Out Reset image size

The distribution of planetocentric pericenter values q is nearly linear up to the Hill radius, with only small variations at the largest values due to details of the migration run (Figure 3). For instance, the drop at ∼0.5 AU for Uranus in run 1 is due to the decrease of the orbital eccentricity at ∼10 Myr (Figure 1) and the lack of interaction with denser regions of the planetesimal disk. Similarly, the decrease of q-distribution for Neptune after ∼0.55 AU in run 2 reflects the variations of the instantaneous Hill radius at different heliocentric distances during the migration of this planet.

The capability to dissociate a given binary planetesimal system is determined mainly by q and V_∞. In particular, it is easy to show that the efficiency of exchange reaction is directly related to the quantity ϕ = GM'/(qV²_∞), where M' is the mass of the planet (e.g., Agnor & Hamilton 2006). In analytic estimations by Agnor & Hamilton the other two fundamental parameters are q itself and ψ = (M/M')(q/a_B); M is the total mass of the binary system and a_B is the pre-encounter semimajor axis of the binary orbit (assumed circular). Note that the exchange reaction only weakly depends on the mass ratio of the binary components. We postpone the analysis of ψ to later sections, because ψ depends on additional parameters in the model such as the distribution of a_B.

The distributions of q and V_∞ characterizing planetary encounters in our migration runs are independent and determined from our migration runs. This makes it easy to construct the expected distribution of ϕ in our model for different planets. Results are shown in Figure 4 for migration runs 1 and 2. The higher the ϕ value is, the higher the planet-intrinsic capability to dissolve the binary.³ Despite systematically larger relative velocities V_∞ (Figure 3), we note that Jupiter has about an order of magnitude larger ϕ values. This is due to the larger mass of the planet and the statistically smaller q for planetesimals encountering Jupiter. In the overall statistics of the expected population of irregular satellites created by this mechanism, the larger ϕ values of Jupiter will be compensated by a much smaller number of encounters N_enc (Table 1 and Section 3). It is worth noting that Agnor & Hamilton (2006) only considered encounters with values of ϕ larger than unity in their model of Triton's capture. Figure 4 indicates such values are rarely attained in the post-PC phase.

Zoom In Zoom Out Reset image size

Having characterized the planetary migration and encounters within the planetesimal disk, we now describe our approach to modeling the dissociation of binary systems. To handle this problem we constructed a numerical code that propagates a binary system during its approach to a planet and monitors the resulting configuration of the three bodies. In this part of our work we use the three-body problem (planet-binary), neglecting the gravitational tidal field of the Sun in the planet's vicinity. We use a high-accuracy Burlish–Stoer integrator with automatically controlled absolute error (e.g., Press et al. 1992). The initial conditions and model parameters are specified as follows.

• 1.  
  We assume that the two components in the binary system have initially a circular orbit. The orbital plane is chosen to be random in space, while the distribution of their separation, a_B, is given by dN ∝ da_B/a²_B. This distribution corresponds to an uniform distribution of binding energies. We note that compact systems, which favor binary dissociation followed by a satellite capture, are preferred in this way. The minimum and maximum values of a_B are set so that (i) the most compact systems correspond to the Roche distance of bodies with a bulk density of 2 g cm⁻³ (e.g., Bertotti et al. 2003, p. 436), and (ii) the most distant systems are at 0.3 of the Hill radius with respect to the solar tidal forces. The nodal positions of the orbital plane of the binary components relative motion in space are chosen randomly.
• 2.  
  The barycenter of the binary system is initially assumed to approach the planet along a hyperbolic orbit with an isotropic distribution of incoming directions, and q and V_∞ parameters are assumed to have distributions identical to those from our migration runs (see, e.g., Figure 3).
• 3.  
  Each simulation uses fixed values of the total mass of the binary system, M, and of the secondary-to-primary mass ratio, μ. In this respect, we consider two different cases for which we conduct two different sets of simulations:

  • (a)  
    Origin of numerous populations of small irregular satellites with sizes ≲100 km (note only a few of them are larger than this limit); we probe the capture efficiency by assuming μ = 0 and M spans a large interval of values, such that the secondary component is always massless and the primary component has various values of mass (see Table 2); and
  • (b)  
    Origin of Triton such that one of the binary components has a mass of Triton, that sets a minimum value of M ∼ 1.64 Pluto masses, while the second component has again a spectrum of masses from zero to large values (see Table 3).

We note that the energy and angular momentum excess of the incoming pair of planetesimals, which needs to be overcome to transfer one of the components onto an orbit bound to the planet, depends on the total mass M of the binary, while it is basically independent of their mass ratio μ. Dependence on this latter parameter may start to play a role in a statistically very infrequent case of q comparable to or smaller than a_B (note, moreover, that for low-mass binaries, a_B is much smaller than the planetocentric distance of regular satellites). As a result, our first set of simulations with a massless secondary (μ = 0) are basically valid for an arbitrary mass ratio. This is important because observations so far indicate that trans-Neptunian binaries have components of comparable mass (or typically within an order of magnitude; e.g., Noll et al. 2007). To verify our conclusions, we also performed a series of simulations for run 1 (Section 3.1) with equal binary components (μ = 0.5) of sizes ranging from 10 to 100 km. None of them provided any captured satellite, even with a large number (10¹⁵) of trials. This is consistent with the results of μ = 0 binaries reported below.

In practice, we start our integration at distance 1.5 AU from the planet and compute the state vector of an hyperbolic orbit with respect the planet. This sets the location and speed of the binary's center-of-mass. The motion of the planet and the two components of the binary system are then numerically propagated through the encounter until the planetocentric distance of any of the components reaches a distance of 1.5 AU again. We then determine the relative state vectors of the primary and secondary components with respect to the planet and check if any of them is on a planet-bound orbit. If so, we record its Keplerian orbital elements with respect to the planet.

We conducted the analysis of satellite capture for each of the four giant planets and for each of the selected migration runs listed in Table 1. In each case we performed a large number of trial simulations, N_trials, changing randomly (but respecting determined distributions) the initial parameters: q, V_∞, a_B, and angular variables. Each simulation assumed fixed values for M and μ.

The efficiency of satellite capture from low-mass binaries is expected to be very low (e.g., Agnor & Hamilton 2006). On the other hand, these systems may be more abundant in the disk than high-mass binaries due to the steep size distribution. Therefore, we must perform a large number of trial runs. We used N_trials ≳ 10¹⁰ in each simulation. To cope with such a huge number of numerical integrations, we adopt an additional filter and integrate only those configurations that satisfy

$$\begin{equation} q \le \kappa a_{\rm B}\left(3 \frac{M^{\prime }}{M}\right)^{1/3}, \end{equation} \tag{ 1 }$$

where the coefficient κ is of the order of unity. The condition in Equation (1) means the binary system is at the Hill stability limit when passing through the pericenter q of its orbit relative to the planet. The coefficient κ>1 is a safety parameter. In our simulations we chose κ = 1.5, except for those of Triton capture, where we set κ = 2.5. All trial cases that violate Equation (1) are not integrated and are assumed to fail in producing satellites. They count, however, toward the overall statistics. Using a limited number of simulations, we have verified that cases violating Equation (1) do not really produce captured satellites. Additionally, in all runs we verified that captured orbits have q well below the critical limit $\kappa \,a_{\rm B} \big(3 \frac{M^{\prime }}{M}\big)^{1/3}$ from the right-hand side of Equation (1) for the κ value used.

The orbital parameters of captured objects recorded by our code are their initial osculating elements. To estimate the long-term stability of these orbits and/or perform comparison with the observed population of irregular satellites, we would ideally need to perform a full-fledged integration of their subsequent evolution. With our large number of trial cases, however, we obtain tens to hundreds of thousands of irregular satellite orbits in each simulation (making in total millions of satellite orbits). Their numerical integration over appropriate time interval would surpass our computer capacity. Therefore, we use an alternative approach to deal with this problem. The main sources of instability for the captured irregular satellites are due to solar perturbations, such as the evection instability and Kozai-cycle instability. Typically, the amplitude of the former is smaller than that of the latter. We thus adopt the simplest approach and compute the minimum and maximum values of the pericenter distance during the Kozai cycle (obtaining thus also minimum and maximum values of the inclination variations). This can be achieved analytically because Kozai dynamics implies two integrals of motion (e.g., Kozai 1962; Kinoshita & Nakai 1999; Carruba et al. 2002; Nesvorný et al. 2003):

$$\begin{equation} \begin{array}{c} K_1 = \cos i\sqrt{1-e^2}, \end{array} \end{equation} \tag{ 2 }$$

$$\begin{equation} \begin{array}{c} K_2 = (2+3e^2)\,(3\cos ^2 i -1) + 15 e^2\sin ^2 i \cos 2\omega , \end{array} \end{equation} \tag{ 3 }$$

where (a, e, i) are the osculating semimajor axis, eccentricity and inclination, and ω is the osculating argument of the pericenter. In addition to K₁ and K₂, the semimajor axis a is also conserved by the Kozai dynamics (unless the dissipation effects become important). With the initial orbital parameters (a, e, i) known, and ω statistically modeled, we determine (K₁, K₂) and compute analytically the minimum/maximum values of e, i, and the orbital pericenter q during the Kozai cycle.⁴ We then apply the following filters to eliminate unstable or inappropriate orbits.

• 1.  
  Small irregular satellites. We eliminate all orbits whose pericenter during the Kozai cycle becomes smaller than the planetocentric distance of the outermost regular satellite, which turns to be 0.013 AU for Jupiter (Callisto), 0.023 AU for Saturn (Iapetus), 0.004 AU for Uranus (Oberon) and 0.003 AU for Neptune (Triton in this case).
• 2.  
  Capture of Triton. We eliminate all orbits whose pericenter during the Kozai cycle becomes smaller than Neptune's radius. We also eliminate those orbits for which the pericenter is always larger than ∼0.004 AU.

In both cases we also eliminate all satellites whose semimajor axis would exceed the long-term stability zone determined by the numerical experiments of Nesvorný et al. (2003). This typically means a ≳ 0.4 R_Hill for prograde satellites and a ≳ (0.5–0.6) R_Hill for retrograde satellites, depending on the time-averaged eccentricity of the orbit (see Figures 9–12 of Nesvorný et al. 2003; R_Hill is the radius of the Hill sphere for the planet).

The first condition described above is based on the assumption that a population of small irregular satellites crossing the orbits of more massive regular satellites would be destroyed by collisions or would be dynamically eliminated from the system. It basically results in the elimination of orbits in some interval of inclination value around 90° (e.g., Carruba et al. 2002; Nesvorný et al. 2003). The more complicated second condition for Triton's capture is based on the recent orbital evolution model of Ćuk & Gladman (2005). These authors postulate that Triton quickly evolved from its initial (capture) orbit by collisions with a pre-existing population of regular satellites similar to those around Uranus (and thus extending to ∼0.004 AU planetocentric distance). Debris from those satellites produced a massive circum-Neptunian disk that made Triton's orbit evolve into its current orbit in only ∼0.1–1 Myr. In particular, the orbital eccentricity of Triton has shrunk from its initial value to ∼0, while its inclination has not likely undergone significant evolution (Ćuk & Gladman 2005). Note, however, that we do not eliminate captured satellites with inclinations inconsistent with Triton's orbit (∼157° with respect to Neptune's equator). This is because we are interested in the intrinsic efficiency of the capture of satellites with Triton-like orbit by the binary splitting process, of which the actual Triton is just one possible realization.

To evaluate the number N_cap of irregular satellites captured in stable orbits with our method, we must include the value of the argument of the pericenter ω. This angle is random at the moment of capture. For this reason, with each captured orbit, characterized by planetocentric (a, e, i) orbital elements, we associate a statistical weight W by randomly choosing 10³ values of ω and evaluating the fraction of (a, e, i; ω) orbits that are stable according to our criteria; in particular, if M such orbits are stable than W = M/10³.

To compare the orbits of the captured satellites with those of existing irregular satellites, we use the mean orbital elements rather than the osculating ones (e.g., Nesvorný et al. 2003). This is because the Kozai cycle may significantly change the osculating elements on a timescale of ∼0.1–1 ky. It is possible to determine time-average eccentricity and inclination values over a Kozai cycle semi-analytically (e.g., Kinoshita & Nakai 1999). Such a computation, however, is not straightforward. We thus use a simpler approximate approach by taking an average of minimum and maximum Kozai-cycle values of e and i instead of their true time averages. They are sufficiently close to the time averages of e and i for the purposes of our paper.

Because of computational limitations, our migration runs in Section 2.1 represented the planetesimal disk by only a few thousands of equal-mass planetesimals (enough, though, to reasonably represent planetary migration and obtain a good statistical description of q and V_∞). Thus, the characteristic particle masses in these simulations were ∼1.6 of Pluto mass, or, assuming a bulk density of 2 g cm⁻³, bodies of ∼2, 920 km in size. In reality, however, the planetesimal disk contains bodies with sizes ranging over the whole spectrum from objects that are very small to objects larger than the ones we used in our migration runs. Bodies in each size range, assuming that some fraction of them are binaries, may be potential sources of irregular satellites. We take this fact into account in the following.

We start with a description of our best guess for the size distribution (SFD) of the planetesimals in the disk. We then assume that the encounter parameters, such as distributions of q or V_∞ determined in Section 2.1, are the same for bodies of different sizes (at least at some plausible interval we shall use) and that the captured population of irregular satellites scales with the estimated number of planetesimals of some particular size. This will allow us to quantify the number N_sat of captured irregular satellites from a disk with a realistic size distribution of bodies (Section 2.4).

Both the observational constraints (e.g., Bernstein et al. 2004; Elliot et al. 2005; Petit et al. 2006, 2007) and the evolutionary models (e.g., Stern & Colwell 1997; Gladman et al. 2001) are, to some degree, uncertain and yield only an approximate insight into the SFD of the planetesimal disk prior to or just after the planet-crossing phase. For simplicity, we shall thus assume a broken power law size distribution with a roll-off at size D₀ and differential distribution exponents q₁ for D ⩾ D₀ and q₂ for D ⩽ D₀. As far as the current population is concerned, Bernstein et al. (2004) found that D₀ ∼ 50 km, q₁∼ 4–5, and q₂∼ 2–3 in the current trans-Neptunian population. While finding nearly the same value of q₁, Gladman et al. (2001) and Petit et al. (2006) argue, however, that the roll-off size D₀ happens for smaller bodies, perhaps between 10 and 25 km. The robust finding that q₁ > 4 and q₂ < 4 implies that bodies with diameters D₀ dominate the mass of the system. We shall assume that the size distributions of the planetesimals in the trans-Neptunian disk and in the post-PC phase were the same as now, only with their total number scaled to characterize its initial mass in our simulations, which is tied to the assumed 35 Earth mass of the primordial planetesimal disk beyond Neptune.

Our assumptions about SFD mean that the number dN of bodies in each size bin dD is dN ∝ D^−q dD. In particular, for D ⩾ D₀

$$\begin{equation} dN = (q_1-1)\,\left(\frac{D_1}{D}\right)^{q_1}\, \frac{dD}{D_1}, \end{equation} \tag{ 4 }$$

and for D ⩽ D₀

$$\begin{equation} dN = (q_1-1)\,\left(\frac{D_1}{D_0}\right)^{q_1}\, \left(\frac{D_0}{D}\right)^{q_2}\,\frac{dD}{D_1}. \end{equation} \tag{ 5 }$$

Here D₁ is an auxiliary size threshold for which N(D > D₁) = 1; we also denote the mass of a body with size D₁ to be M₁. The total mass, M_disk, of the planetesimal disk is given by

$$\begin{equation} M_{\rm disk} = M_0\,\frac{(q_1-1)(q_1-q_2)}{(q_1-4)(4-q_2)}\, \left(\frac{D_1}{D_0}\right)^{q_1-1}. \end{equation} \tag{ 6 }$$

The shape of the assumed size distribution is determined by the (q₁, q₂, D₀) parameters, and the absolute number of objects in the population is defined by the D₁ value. The latter parameter is not known a priori, especially in the early post-PC phase. We, however, circumvent this parameter by using Equation (6) that relates D₁ to M_disk known from the initial data of our migration simulations. In particular, we have that

$$\begin{equation} D_1 = D_0\,\left[\frac{(q_1-4)(4-q_2)}{(q_1-1)(q_1-q_2)}\, \frac{M_{\rm disk}}{M_0}\right]^{1/(q_1-1)}. \end{equation} \tag{ 7 }$$

As an example of the (q₁, q₂, D₀; M_disk) parametrization of the disk we may use M_disk ∼ 12.4 Earth masses (the case of run 1). To get the feeling how the (q₁, q₂, D₀; M_disk) parametrization determines our results, we use for definiteness q₁ = 4.5, q₂ = 2.8 and D₀ = 50 km. With these values, Equation (7) yields M₁ ≃ 44.6 Pluto masses or about one body of Mars mass in the disk population. We would thus expect, as a result of our mathematical exercise, about one body of that mass in the disk population and about 90 bodies of the Pluto mass.

It is interesting to note that (see Equation [7])

$$\begin{equation} M_1\propto M_0^{(q_1-4)/(q_1-1)}\;. \end{equation} \tag{ 8 }$$

Since q₁ > 4, the exponent on the right-hand side of Equation (8) is positive. This means that smaller D₀ values would result in smaller D₁ for a fixed value of M_disk. For example, with the same slope exponents as above, M_disk ∼ 12.4 Earth masses and D₀ = 25 km, the disk would hold one body with ∼33 Pluto masses and some 60 bodies of the Pluto mass. In addition, with D₀ fixed, larger q₁ values result in smaller D₁ for a fixed value of M_disk.

Finally, we note that the size distribution discussed in this section characterizes the total population of disk objects. Our work further considers a fraction of them are binaries. This subpopulation is assumed to obey the same mass distribution as the whole system, or in other words the fraction of binaries is assumed to be roughly constant in the size interval of interest. In further work we actually assume all planetesimals are binaries, and only the final results would eventually be multiplied by the expected binary fraction (e.g., Noll et al. 2007).

In the final step, we use all previously collected information to estimate the expected population of the irregular satellites captured by dissociation of the planetesimal binaries. If the equal-mass population of objects in our simulation were to represent the reality, N_cap satellites become captured in stable orbits from N_trials trial cases (described in Section 2.2), and f_bin is the fraction of binaries, we would expect to capture N_sat ∼ f_bin N_cap(N_enc/N_trials) satellites. With a continuum of planetesimal sizes (Section 2.3) such that at each infinitesimal size bin dD we have dN objects (Equations [4] and [5]), we obtain

$$\begin{equation} N_{\rm sat}\sim \int f_{\rm bin}\,N_{\rm cap}\, \left(\frac{N_{\rm enc}}{N_{\rm trials}}\right) \left(\frac{d N}{N_{\rm TP}}\right), \end{equation} \tag{ 9 }$$

where integration is performed over bodies of all sizes D. Note that N_cap = N_cap(D) depends on the size/mass of the binary system which underwent dissociation in the gravity field of the planet. Similarly, f_bin = f_bin(D) may likely be a size-dependent parameter.

In practice, we choose finite size/mass bins (D, D + ΔD) and replace integration in Equation (9) with a finite summation:

$$\begin{equation} N_{\rm sat} \sim \sum _{\rm bins}\, f_{\rm bin}\,N_{\rm cap}\, \left(\frac{N_{\rm enc}}{N_{\rm trials}}\right) \left(\frac{\Delta N}{N_{\rm TP}}\right). \end{equation} \tag{ 10 }$$

The value ΔN of number of objects in the diameter bin of size ΔD around D is given by Equations (4) and (5).⁵ Our choice for the size bins is given in Table 2 for experiments with the capture of small (<50 km size) irregular satellites, and in Table 3 for the capture of Triton. In both cases we specify the bins by minimum D_l and maximum D_r diameters. Moreover, since M ∝ D³, we define the effective diameter $D_{\rm eff}=\int _{D_{\rm l}}^{D_{\rm r}} \xi ^4\,d\xi \big/\!\int _{D_{\rm l}}^{D_{\rm r}} \xi ^3\,d\xi = (4/5) \big(D_{\rm r}^5-D_{\rm l}^5\big)\big/ \big(D_{\rm r}^4-D_{\rm l}^4\big)$ and the corresponding mass M_eff, assuming a bulk density of 2 g cm⁻³.

Our simulations of the capture of small irregular satellites assume a massive primary and massless secondary (thus μ = 0). The number of trials N_trials in each size bin is given in Table 2. The simulations of Triton's capture must take into account that one component of the binary is a ∼1.64 Pluto mass object. This also sets the minimum total mass of the binary system and fixes the mass ratio μ for each of the bins. For Triton's capture simulations, we used N_trials = 2.5 × 10¹⁰ (Table 3).

In order to verify the robustness of our results on changes of the chosen bin sizes, we also performed a test for run 1 where we took 20 bins for the capture of small irregular satellites and 25 bins for the capture of Triton (the size range between the smallest value of D_l and the largest value D_r was the same). In both cases we also took twice as many trial integrations N_trial as indicated in Tables 2 and 3.

We first analyze the ability to capture small irregular satellites from our runs where the secondary in binary systems was considered massless and the mass of the primary sampled a broad interval of values (Table 2).

Figures 5–8 show the values of the semimajor axis a and the Kozai-cycle mean value of the eccentricity e for captured satellites that meet stability criteria from Section 2.2. We give results for bins 5–8 (i.e., large primary components in the binary systems), because bins 1–4 provide only very limited (or even zero) number of captured satellites, typically with orbits of large a (close to stability limit) and very high e (typically e > 0.8). We nevertheless include them in the statistical analysis below. In all cases, capture from binary systems results in orbits with systematically higher values of the orbital eccentricity than that of the observed satellites (red triangles). The discrepancy is the largest for Jupiter, but even for Neptune the match is not good.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next we used Equation (10) to obtain the total number N_sat of small irregular satellites by adding captured populations in each of the size bins. Table 4 gives results for selected disk parameters (D₀, q₁, q₂). We also explored (D₀, q₁, q₂) more extensively, namely assuming D₀ in the (10, 50) km range, q₁ in the (4.1, 4.9) range, and q₂ in the (2.3, 3) range. Performing 10⁶ random trials, the maximum number of captured irregular satellites around Neptune (in a statistical sense explained in Section 2.4) was ≃0.25. We also found that a smaller population (by about an order of magnitude) is expected around the other giant planets. This result is in contrast with the finding of Sheppard et al. (2006); see also Jewitt & Sheppard (2004) or Nicholson et al. (2007), who show that the D ⩽ 10 km population of irregular satellites is numerous and broadly similar for each of the giant planets. Our predicted ∼0.01–0.1 captured satellites is about three orders of magnitude below the observed population.

In order to check how sensitive our results are on the particular choice of 8 size bins from Table 2, we repeated the simulation with 20 size bins that covered the same interval (50–10, 000 km) for the primary size. Secondaries were again massless and we also performed twice as many trial integrations. For a given set of disk parameters (D₀, q₁, q₂) we obtained ≃30% fewer captured satellites N_sat. This difference is not significant. Because of much lower computer time expense, we shall use the 8-bin set in the following.

We took advantage of our test run to probe which binary systems contribute most to the population of the captured irregular satellites. Figure 9 shows the contribution ΔN_sat = N_cap(N_enc/N_trials)(ΔN/N_TP) of each of the bins in the final value N_sat for all of the planets. Note that much more massive primaries are necessary to provide satellites around Jupiter than around Neptune. The primaries with a peak efficiency for the delivery of small irregular satellites have ∼5, ∼1, ∼0.05, and ∼0.002 Pluto-masses for Jupiter, Saturn, Uranus, and Neptune, respectively. Comparison of ΔN_sat for the simulations with 8 and 20 size bins also allows us to understand why the nominal case with 8 bins overestimates the number of captured satellites (Figure 9).

Zoom In Zoom Out Reset image size

3.1.1. Capture of Triton

Figure 10 shows osculating orbital parameters for satellites meeting the criteria of Triton capture (Section 2.2). Each of the panels corresponds to one individual size bin outlined in Table 3 (the bin identifier is labeled in the right bottom corner). Obviously, a smaller number of Triton-like objects are captured from systems of lower total mass (the initial bins 1, 2, 3...), while much more are captured when the mass of the primary increases (the last bins ..., 9, 10). In the first case the initial eccentricity is systematically very high (≳0.9), while only massive primaries allow lower eccentricity capture orbits. The importance of each of the bins for Triton capture probability depends not only on the number of the captured satellites N_cap, but also on the estimated number ΔN of the binary systems in the planetesimal disk (see Equation [10]). Figure 11 shows the contribution ΔN_sat = N_cap(N_enc/N_trials)(ΔN/N_TP) of each of the bins in the final value N_sat for three different values of (D₀, q₁, q₂) triples. Interestingly, bins 2 and 3, corresponding to about 2 to 3 mass ratio between the primary and secondary (Triton), are preferred for the Triton capture. Systems with a lighter primary are more abundant in the disk but provide comparably less captured satellites. On the other hand, very massive binary systems yield a lot of captured satellites, but are very rare in the disk.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 12 shows the expected number of satellites N_sat meeting the criteria of Triton capture (Section 2.2) when the parameters of the planetary-disk distributions are varied in their admissible intervals. In particular, we considered D₀ values in the interval (10, 50) km, q₁ values in the interval (4.1, 4.9), and q₂ values in the interval (2.3, 3). Overall, the chance to capture Triton during the 100 Myr interval during the post-PC phase is less than 2%. The highest values of N_sat are obtained for q₁ ∼ 4.1–4.3, while observations tend to prefer higher values of this parameter (e.g., Gladman et al. 2001; Bernstein et al. 2003; Elliot et al. 2005; Petit et al. 2006). The right panel in Figure 12 shows a correlation between N_sat and q₁, which is mainly due to the strong influence of q₁ on the number of planetesimals ΔN in the size bins used in our simulation. We found that N_sat is significantly less dependent on D₀, and nearly independent of q₂.

Zoom In Zoom Out Reset image size

We have re-run our simulation with 25, instead of 10, size bins and we obtained results that predict only 7% fewer satellites (see Figure 11). This accuracy is enough for our conclusions and it justifies our usage of coarser bins from Table 3 for further runs.

The orbits of the captured massless irregular satellites in run 2 indicate the same inconsistency with the observed population, as shown in Figures 5–8 for run 1. In particular, orbits with high eccentricities are preferentially populated leaving the space of low-eccentricity orbits, where the observed satellites reside, depleted or entirely void. As expected, the total number of captured satellites is up to one order of magnitude larger than in run 1. This is especially true for Uranus and Neptune, because these two planets penetrate into a heavier residual disk with less scattered orbits of planetesimals, such that the efficiency factor of capture ϕ is a factor larger than in run 1 (see Figure 4). Nevertheless, N_sat for all planets is still significantly smaller than unity and several orders of magnitude smaller than the observed population. The exception is Neptune, for which we would predict about one captured irregular satellite, obviously not enough to explain its debiased population of small irregular satellites (e.g., Jewitt & Sheppard 2004). Moreover, the predicted population of small irregular satellites at Saturn is more than three orders of magnitude smaller than that at Neptune, contradicting the observations (e.g., Jewitt & Sheppard 2004). In the same way, although our expected likelihood of capturing one Triton-like orbit increases by a factor of a few, yet, at its maximum, it reaches only a few percent.

The planet-crossing phase of the Nice model is violently chaotic. The planetary configuration at its end, the degree of the planetesimal disk depletion, and its dynamical excitation at different heliocentric distances is unknown and strongly dependent on the sequence of planet encounters. As a result, the realistic uncertainty of the irregular satellite capture does not derive only from parameter details of simulations described above (such as the number of size bins, the number of trial integrations, the size-distribution parameters D₀, q₁, or q₂ etc), but is likely dominated by variations in the possible planetary evolution scenarios. To face this fact, we used the results of Nesvorný et al. (2007), who determined 14 different possibilities of the planet evolution after Jupiter and Saturn entered the 2/1 mean motion resonance. Here we selected seven of them and applied our method to estimate the number of captured irregular massless satellites and likelihood of capturing Triton during the first 100 Myr of the post-PC phase.

Table 1 indicates that the mass of the residual trans-Neptunian disk and the number of planetesimal encounters with planets varied by nearly a factor of 2 in these different variants of planetary evolution. Overall, the resulting population of irregular satellites does not deviate from that in runs 1 and 2 (described above) by more than a factor of 2 (see Table 4). In none of the cases did the number of captured satellites come close to the observed population. In addition, the likelihood of Triton's capture remained between 0.1 and 3%. For the small satellites, we also noted in all runs that the eccentricity of the captured satellites is significantly smaller than observed.

We note this orbital inconsistency occurred in all our simulations and presents, in our opinion, the strongest argument against the capture mechanism. It might, in principle, only be removed by assuming subsequent weak orbital evolution that could pull-down orbital eccentricities while preserving the semimajor axis after satellite capture (e.g., McGleam et al. 2006). However, in our mind, no such universal process has been found as yet. In the absence of gas drag at the late stage of planetary evolution described in this paper, tidal evolution remains the most important long-term effect. This is too weak and it would also primarily drive the evolution along lines of constant pericenter distance q. Another possibility is mutual collisions between the satellites, which may be an interesting evolutionary process yet to be studied in depth, but below we argue that the initially captured population is not populous enough to explain the discrepancy via this mechanism.

Note that we always assumed f_bin = 1, which likely makes our population of captured satellites an overestimate. We purposely avoid closer analysis of how f_bin(D) affects our results because too little is known about this function. Current observations suggest a 22⁺¹⁰₋₅% abundance of binaries in the classical Kuiper belt, with smaller fractions observed in the hot Kuiper belt and the scattered disk (Noll et al. 2007). This fraction might be, however, larger among the largest trans-Neptunian objects (e.g., Brown et al. 2006), and if these objects possess more than one satellite, one may effectively have also f_bin ⩾ 1. This would slightly enhance the role of the satellite capture from the population of the largest planetesimals. Another process that might slightly increase the role of satellite captures from these largest planetesimals is that their orbits could be dynamically less excited than those of the smaller disk components, such that their relative encounter velocity with planets might be systematically smaller. We have seen, however, that the largest objects do not contribute dominantly to the overall capture probability. We thus do not expect that these omitted effects would significantly increase the efficiency of irregular satellites production from the exchange reactions that was several orders of magnitude low in our simulations.

All the above-mentioned conclusions apply to the interval of 100 Myr following the planet-crossing phase of the solar system evolution according to the Nice model. We now argue how they might be extrapolated further in time, but also to the pre-PC phase of the Nice model. First, we comment on the possible satellite population captured during the phase post-dating the 100 Myr analyzed previously in Section 3. In spite of the fact that this phase lasted another ∼3.7 Gyr, we see from Figure 2 that the number of planetesimal encounters with planets decreases very fast. It is not a straightforward exercise to extrapolate this observed trend because a quasi-exponential match of the 100 Myr integrated time interval would provide an unrealistically low-mass trans-Neptunian population today. Clearly, at some time during the late evolution the rapid decay rate of the disk mass (e.g., Figure 2) and the number of planetary encounters had to change into a slower evolution. However, we argue that rather than extrapolating our previous results using a ratio of timespan it makes more sense to use the ratio of mass losses from the planetesimal disk. Indeed, if the number of planetary encounters in Figure 2 decays proportionally to the disk mass, their ratio is constant. So the accumulated number of the planetary encounters during the past ∼3.7 Gyr might be, in the first approximation, proportional to the remaining mass loss from the disk. In run 1, for instance, the final disk mass of ∼4 Earth masses should be next decreased to today's value by another ∼4 Earth masses (note, for the sake of our argument, we do not actually need to know the small residual mass of the trans-Neptunian population today). This is a factor of ∼2 less than what the planetesimal disk lost during the above analyzed 100 Myr interval. As a result, we would expect about half the number of irregular satellites created in this late phase than during the time interval investigated in Section 3. Overall, this would mean a negligible population.

Next, we comment on the pre-PC phase of the Nice model. It is important to note that this interval of time does not include as yet the planetary formation phase during which the giant planets formed within a comparably cold planetesimal disk and were surrounded by a residual gas disk. What we have now in mind is the phase covered by analysis of Gomes et al. (2005), notably starting from a compact configuration of outer planets that already have cleared most of their orbital space. An edge of a ∼35 Earth mass planetesimal disk is located about ∼1–1.5 AU beyond the outer planet, which was found to provide a ∼600–800 Myr long quasi-stable period before LHB takes place. Note, however, that our analysis remains valid even if the PC phase of the Nice model occurred earlier in the evolution of the solar system and was not related to the LHB. During the pre-PC time interval outer planets experience only little migration evolution, preserving their compact configuration until Jupiter and Saturn reach their mutual 2/1 mean motion resonance. This is reflected by the fact that disk particles only occasionally reach the planet zone.

If we apply the disk mass-depletion scaling for the number of planetary encounters, we may expect a comparable number or fewer satellites created during the pre-PC phase than in the 100 Myr interval after planets stopped crossing their orbits (note the initial disk mass had 30 Earth masses and it decreased to about 24 Earth masses before the PC phase). This assumes that the encounter parameters, such as distribution of V_∞ or q, are the same as described in Section 3. Indeed, we verified this using two simulations. At first, we started again the initial ∼20 Myr of run 1 (Figure 1), during which the gas giants did not yet reach their mutual 2/1 mean motion resonance. During this time interval we recorded parameters of planetary encounters and verified that distributions in Figure 3 still hold. Moreover, the frequency of the encounters was about two orders of magnitude less than during the post-PC phase, as we expected. We also performed a series of test simulations putting the initial planetary system described by Gomes et al. (2005) interior to massive disks with various surface density profiles. We again integrated the evolution for about 20 Myr and verified that even in this initial phase of Gomes et al. evolution, V_∞ and q, with respect to planets, have approximately the same distribution as in Figure 3.

However, we have to point out that a comparable, or even a factor of few larger, captured population of irregular satellites before the PC phase is not enough to explain their observed population today. We recall the results of Nesvorný et al. (2007), who determined that the planet-crossing phase should decimate the pre-existing irregular population by a factor of ∼1000 or more. As a result, the pre-PC captured satellites from binaries should be even much less populous than those created later. This argument does not apply to Triton's capture, because once evolved into its current orbit it would not likely be destabilized during the planetary-encounter phase. Yet, the small factor of increase does not bring the probability of capturing Triton in the pre-PC phase of the Nice model any closer to unity.

We may thus conclude that since the giant planets completed their formation, the surrounding planetesimal disk was cleared, and gas was dissipated from the solar system (the initial moment of the current version of the Nice model), irregular satellite capture from disk binaries is not a viable mechanism to provide the overall populations of irregulars or capture of Triton: the captured orbits are inconsistent with the observed ones and the efficiency is low. It is fair to recall that Jupiter may be a special case, since this planet does not typically participate in close planetary encounters within the current version of the original Nice model (although see a novel route by Morbidelli et al. 2007 that may possibly involve Jupiter's participation in planetary close-encounters and a discussion in Nesvorný et al. 2007). It is possible that its population of irregular satellites was created in part by dissociation of the primordial binary systems in the planetesimal disk very early after (or simultaneously with) Jupiter's formation. A number of Jupiter's irregulars in secular resonances are considered to suggest their early formation in the gaseous environment. Our work, however, suggests that the capture from binaries is likely ineffective in later epochs because of the poor match of orbital characteristics (Figure 5) and a small number of captured satellites.

As far as Triton is concerned, our simulations above make a capture possible onto orbits that could then evolve to the precursor of its currently observed orbit. However, they also indicate that it is unlikely to happen after the planet-crossing phase within the Nice model. In fact we only see a possibility of its capture from a binary system very early after planet formation. In this early phase, orbits of both growing planets and planetesimals in the disk were kept cold by the action of surrounding gas. This means that the relative velocity at the encounters between the planetesimals and the planets was much lower than shown in Figure 3 for the later phase and more consistent with assumptions in Agnor & Hamilton (2006). Whatever formation process of disk-binaries dominated, it is likely that it was operational during this early phase of planet's evolution (in fact, might have been at the peak of its efficiency, e.g., Goldreich et al. 2002; Astakhov et al. 2005). Given the time constraints for this phase, it would be interesting to study constraints on disk parameters (such as its state of excitation, its mass and planetesimal size distribution) to attain an order of unity probability to capture Triton.