THE IRREGULAR SATELLITES: THE MOST COLLISIONALLY EVOLVED POPULATIONS IN THE SOLAR SYSTEM

The irregular satellites are defined as objects that orbit far enough away from their primary (i.e., Jupiter, Saturn, Uranus, and Neptune) that the precession of their orbital plane is controlled by solar rather than planetary perturbations. This distance was defined to be a semimajor axis a > (2μJ₂R²a³_p)^1/5, where μ is the ratio of the planet's mass to that of the Sun, J₂ is the planet's second zonal harmonic coefficient, R is the planet's radius, and a_p is the planet's semimajor axis (e.g., Burns 1986; Jewitt & Haghighipour 2007).

The semimajor axes a of the satellites around Jupiter to Neptune are more easily described if we scale them by their primary's Hill sphere, defined as R_H = a_p(m_p/(3 M_☉))^1/3, where m_p is the mass of the planet. Collectively, the irregular satellites have a/R_H values between ∼0.1 and 0.5. The prograde objects, however, have a smaller range of values (0.1–0.3) than the retrograde ones (0.2–0.5). Their eccentricity e values range between 0.1 and 0.7, with the retrograde satellites generally having larger e values than the prograde ones. Their inclination i values, to zeroth order, are similar to one another, with the prograde and retrograde populations avoiding 60° < i < 130° where the Kozai resonance is active (Carruba et al. 2002, 2003; Nesvorný et al. 2003). Concentrations of irregular satellites in (a, e, i) space, which are likely families produced by collision events (e.g., Nesvorný et al. 2004), will be discussed in greater detail below.

The orbital stability zones of the irregular satellites were determined numerically by Carruba et al. (2002) and Nesvorný et al. (2003; see also Beaugé et al. 2006). In general, these zones are modestly larger in (a, e, i) space for retrograde objects than for prograde ones. The known irregular satellites fill a large fraction of these zones; there are no obvious voids that require explanation. This implies that the primordial satellite population was once larger, with objects captured onto nearly isotropic orbits. The objects injected onto unstable orbits, however, were either driven inward toward the primary, where they collided with the primary or the primary's regular satellites, or outward to escape beyond the planet's Hill sphere. Numerical simulations indicate that the majority of these entered the zone of the regular satellites (and were presumably lost by collisions) within <10⁴ yr of capture.

The irregular satellites have colors that are consistent with dark C-, P-, and D-type asteroids that are prominent in the outer asteroid Belt and dominate the Hilda and Trojan asteroid populations (e.g., Grav et al. 2003, 2004; Grav & Bauer 2007; Grav & Holman 2004). Spectrally, these objects are a good match to the observed dormant comets. Levison et al. (2009) claim that this may not be a coincidence and use numerical simulations to suggest and that the aforementioned populations may be captured refugees from the massive primordial trans-planetary disk that once existed beyond the orbit of Neptune (see also Morbidelli et al. 2005). This issue will be discussed in greater detail below.

The cumulative SFDs of the irregular satellites are shown in Figure 1. The data on this plot, taken from Scott S. Sheppard's Jovian planet satellite Web site at http://www.dtm.ciw.edu/sheppard/satellites/, contain what is known of the irregular satellite population as of 2009 May 1. Jupiter has 55 irregular satellites (7 prograde, 48 retrograde), Saturn has 38 (9 prograde, 29 retrograde), Uranus has 9 (1 prograde, 8 retrograde), and Neptune has 6 (3 prograde, if one counts Nereid, and 3 retrograde). The decrease in satellite numbers as one moves further from the Sun corresponds to the increasing difficulty observers have in detecting small dark distant objects (e.g., Sheppard et al. 2006). The satellite diameters were computed using Equation (2) of Sheppard et al. (2005, 2006) under the assumption that the objects have albedos of 0.04, values that are typical for dormant comets (e.g., Jewitt 1991) and C-, D-, and P-type asteroids (e.g., Fernández et al. 2003).

It is important to point out here that in terms of detection and discovery, irregular satellites are different from asteroids or TNOs. For the latter, sky coverage, distance from opposition and the observer, recovery, etc. are all important issues that vary from survey to survey. For irregular satellites, however, the Hill spheres for planets like Saturn, Uranus, and Neptune have been completely covered by observations down to the limiting magnitude of the specific survey (e.g., Gladman et al. 2001; Sheppard et al. 2005, 2006). Accordingly, observational incompleteness in each population does not gradually fall off but instead is nearly a step function of magnitude. In other words, the satellites are close to 100% complete down to magnitude M₁ and close to 100% incomplete for magnitude M > M₂ = M₁ + dM, where dM is small. Values of M₁ and M₂ are discussed for individual planets in various irregular satellite detection papers (a list of useful references can be found in Jewitt & Haghighipour 2007 and Nicholson et al. 2008).

Note that Jupiter's Hill sphere has yet to be fully covered, so completeness is not yet at 100%. As described in the introduction, it has also been argued that the prograde and retrograde SFDs, when combined together, are surprisingly similar to one another, particularly when observational incompleteness at Uranus and Neptune is considered (Jewitt & Sheppard 2005; see reviews in Jewitt & Haghighipour 2007 and Nicholson et al. 2008). The combined prograde and retrograde SFDs at each planet have a cumulative power-law index of q ≈ −1 for 20 < D < 200 km (i.e., cumulative number N(>D) ∝ D^q). Using the data from Figure 1, we find that this value also holds, more or less, down to D > 8 km objects (q ≈ −0.9), though Jupiter and Saturn have shallower slopes over this size range (q ≈ −0.8 and −0.7, respectively), and the prograde population for Saturn has a steeper power-law slope for D > 8 km objects than the retrograde ones. The retrograde SFDs for D < 8 km bodies at Jupiter and Saturn, where we have data, then steepen up. The steepest slope, q = −3.3, is found for Saturn's retrograde satellites between 6 < D < 8 km.

Other small body populations, measured over the same size ranges as those described above, cannot rival the shallow power-law slopes seen among the irregular satellites. Instead, most have q values near −2; the outer main Belt, Hilda, and Trojan SFDs for 20 < D < 100 km all have q ≈ −2 (e.g., Figure 2; see also Levison et al. 2009), the main Belt and near-Earth object (NEO) populations for 1 < D < 10 km is from −1.8 to −2.0 (e.g., Bottke et al. 2002, 2005a, 2005b; Stuart & Binzel 2004), and the ecliptic comet SFD for D > 3 km is q = −1.9 ± 0.3 (see the review by Lamy et al. 2004). Even fading or dormant comet populations, which are or have been strongly affected by mass-loss mechanisms and disruption events (e.g., Levison et al. 2002), do not quite reach the shallow slopes of the irregular satellites. Faint Jupiter-family comets (JFCs) that become NEOs and achieve perihelion <1.3 AU have both a steep branch for bright comets (q = −3.3 ± 0.7 for absolute comet nuclear magnitudes H₁₀ < 9) and a shallow branch for dimmer comets (q = −1.3 ± 0.3 for 9 < H₁₀ < 18; Fernández & Morbidelli 2006). Most observed dormant comets have sizes between 2 km < D < 20 km. Those classified as JFCs have q = −1.5 ± 0.3 (Whitman et al. 2006), while those classified as external returning comets and Halley-type comets have q = −1.4 or perhaps even q = −1.2, depending on the modeling parameters used (Levison et al. 2002). Thus, we conclude the irregular satellite for D > 8 km has some of the shallowest, if not the shallowest, SFDs in the solar system.

While we find the similarities intriguing, the various SFDs also have distinct differences that may provide important clues to their origins. For example: (1) Jupiter's retrograde SFD for D < 8 km is shallower than Saturn's, while the shapes of their prograde SFDs are very different from one another. (2) Uranus's prograde population is limited to a single D > 20 km object, while those for Jupiter and Saturn have 3–4 such objects. (3) There is nearly an order of magnitude difference in diameter between Phoebe and the second largest retrograde irregular satellite of Saturn. No comparable size difference can be found among any sub-population in Figure 1 unless one counts the size difference between Nereid, which may be a regular satellite that was scattered during the capture of Triton (Goldreich et al. 1989; Banfield & Murray 1992; Agnor & Hamilton 2006), and its two prograde brethren. (4) Jupiter's and Saturn's prograde populations are comparable to or significantly larger than their retrograde ones for D > 8 km.

Irregular satellite families are clusters of objects with similar proper (a, e, i) parameters with respect to the primary planet. They are produced by cratering or catastrophic impact events, the latter defined as an impact event where 50% of the mass is ejected at escape velocity from the parent body.

Families appear to be an important component of the inventory of irregular satellite populations and can be identified once several tens of objects have been found. Among Jupiter's irregular satellites, Nesvorný et al. (2003, 2004) identified two robust retrograde families using clustering techniques and numerical integration simulations. They are the Carme family, which includes D = 46 km Carme and 13 members with D = 1–5 km, and the Ananke family, which includes D = 28 km Ananke and 7 members with D = 3–7 km. Himalia, a prograde D = 160 km satellite, may also be in a family of four objects with D = 4–78 km. When combined, these families represent about half of Jupiter's known prograde and retrograde satellites. For additional details, see Beaugé & Nesvorný (2007).

A potential problem with the Himalia family is that the dispersion velocity between the members is significantly larger than those inferred for well-studied main Belt asteroid families (Bottke et al. 2001; Nesvorný et al. 2003). One possible solution for this was identified by Christou (2005), who used numerical integration simulations to show that the large dispersion velocities of three members (with D = 4–38 km) could have been produced by numerous close encounters with Himalia itself. If true, some irregular satellites families may be too dispersed to find by clustering algorithms alone. A second potential solution with implications for all irregular satellite is discussed below.

At Saturn, Gladman et al. (2001) identified two prograde groups orbiting Saturn that contain three and four objects, respectively, with D = 10–40 km. If real, these groups constitute ∼80% (seven of nine) of Saturn's known prograde satellites. A retrograde group of four objects with D = 7–18 km has been linked to Phoebe (D = 240 km), though like the Himalia family the dispersion velocities between Phoebe and its putative members are larger than those found among main Belt families (Nesvorný et al. 2003, 2004). If high dispersion velocities are common among irregular satellite families with large parent bodies, perhaps produced by mutual close encounters (e.g., Christou 2005) or some other mechanism (see below), the Phoebe family could have many more members. Thus, like Jupiter, collisional families are likely an important, perhaps dominant, component of Saturn's irregular satellite population.

The irregular satellite populations around both Uranus and Neptune have less than 10 members, too few to probe for families in a meaningful way.

Ideally, we would like to use irregular satellite families to constrain the collisional evolution of irregular satellite systems. After some consideration, however, we decided to shy away from doing so in this paper. Our rationale is that the study of family-forming events among irregular satellites is far less mature than those among main Belt asteroids. At present, it is fair to say that no one yet understands the ejecta mass–velocity distribution function produced by irregular satellites or comet impact events.

To illustrate this complicated and somewhat bizarre issue, consider that while P/D-types are fairly common beyond the outer main Belt, we have yet to identify a robust P/D-type family anywhere in the solar system. Tests indicate P/D-types do not exist in the outer main Belt (Mothé-Diniz & Nesvorny 2008) nor the Hilda/Trojan populations (Brož & Vokrouhlický 2008). For the latter, Brož & Vokrouhlický (2008) investigated the 1200 known Hildas and 2400 known Trojans for asteroid families. Despite the fact that 90% of the Hildas/Trojans are D/P-types, only C-type families were found; two robust families produced by D > 100 km disruption events in the Hilda population and one in the Trojan population. For reference, C-type asteroids have flat, somewhat nondescript spectra and are common across the main asteroid Belt.

At present, we have only begun to explore the physics of impacts into primitive volatile-rich highly porous materials. Preliminary experimental results indicate that when high-velocity projectiles are shot into highly porous targets, impact-generated shock waves may heat up porous material so much that, although the target is under pressure and is in compression, the density decreases because of extensive heat production and resulting thermal expansion (e.g., Holsapple & Housen 2009). This effect could eject fragments at much higher speeds than those inferred from asteroid family-forming events. In turn, this would prevent standard clustering algorithms from finding all of the important members of a family, particularly when few objects are known. All in all, this could provide the easiest solution to the Himalia family conundrum described above. It is even possible that some impact events among irregular satellites launch fragments onto unstable orbits where they could strike the regular satellites (see Section 5.3).

For these reasons, in this paper we refrain from using irregular satellite families as serious constraints for our collisional evolution model.

We infer that while the irregular satellites share many of the same physical characteristics as other small solar system bodies, they evolved in dramatically different ways. Our key points from the data are as follows.

• 1.  
  The irregular satellites were captured on nearly isotropic orbits and were once more numerous than they are today. The size of the initial population is unknown.
• 2.  
  The D > 8 km objects have shallow slopes unlike any small body population yet observed in the solar system. This implies that they were affected by mechanisms in a manner and/or to a degree unlike the C-, P-, and D-type bodies found among the outer main Belt, Hilda, and Trojan asteroid populations or the active/dormant comets.
• 3.  
  The prograde and retrograde satellite SFDs for each planet are similar in certain ways (e.g., population size, overall shape of SFDs) but are distinctly different in other ways (e.g., size and shape of the prograde versus retrograde populations). This makes it difficult to characterize the SFDs without understanding how they reached this point in their evolution.
• 4.  
  Families are an important component of observed irregular satellite populations with several tens of members (i.e., Jupiter, Saturn). We infer from this that impacts have been a key factor in the evolution of the irregular satellite populations. Our understanding of the impact physics controlling irregular family-forming events, however, is still in its infancy.

Nearly all recent papers discussing the origin of the irregular satellites have two sections in common: a section detailing previous attempts to explain the capture of planetesimals from heliocentric to observed irregular satellite orbits and a second section describing the deficiencies of those models. Accordingly, we will try to keep our discussion of these issues brief. The main capture scenarios described to date are: (1) capture by collisions between planetesimals (Colombo & Franklin 1971; Estrada & Mosqueira 2006); (2) capture due to the sudden growth of the gas giant planets, which is often referred to as the "pull-down" capture method (Heppenheimer & Porco 1977), (3) capture of planetesimals due to the dissipation of their orbital energy via gas drag (Pollack et al. 1979; Astakhov et al. 2003; Ćuk & Burns 2004; Kortenkamp 2005), (4) capture during resonance-crossing events between primary planets at a time when gas drag was still active (Ćuk & Gladman 2006); (5) capture in three-body exchange reactions between a binary planetesimal and the primary planet (Agnor & Hamilton 2006; Vokrouhlický et al. 2008), and (6) capture in three-body interactions during close encounters between the gas giant planets within the framework of the so-called Nice model (Nesvorný et al. 2007, and see below).

The problems with most of these models are discussed in several papers; recent summaries can be found in Jewitt & Haghighipour (2007), Nesvorný et al. (2007), Nicholson et al. (2008), and Vokrouhlický et al. (2008). Essentially, all of these models (except possibly step (6); see below) are unsatisfying at some level because they suffer from one or more of the following problems: they are underdeveloped, they are inconsistent with what we know about planet formation processes and/or planetary physics, their capture efficiency is too low to be viable, they require exquisite and probably unrealistic timing in terms of gas accretion processes or the turning on/off of gas drag, they cannot reproduce the observed orbits of the irregular satellites, and they can produce satellites around some but not all gas giants. Moreover, there is the newly recognized problem that if the outer planets migrate after the capture of the irregular satellites, the satellites themselves will be efficiently removed by the passage of larger planetesimals or planets through the satellite system. This means that while different generations of irregular satellites may have existed at different times, the irregular satellites observed today were probably captured relatively late in a gas-free environment.

We argue here that the model with the fewest problems is No. 6 from Nesvorný et al. (2007), with late planetary migration acting not as a "sink" but rather as the conduit for satellite capture. In their scenario, migration leads to close encounters between pairs of gas giant planets over an interval of several millions of years. This allows planetesimals wandering in the vicinity of the encounter site to become trapped onto permanent orbits around the gas giants via gravitational three-body reactions.

To get gas giant close encounters, Nesvorný et al. (2007) invoked the so-called Nice model framework (Tsiganis et al. 2005; Morbidelli et al. 2005; Gomes et al. 2005) that assumes the Jovian planets experienced a violent reshuffling event in the past (presumably ∼3.9 Ga). The starting assumption for the Nice model is that the Jovian planets formed in a more compact configuration than they have today, with all located between 5 AU and 15 AU. Slow planetary migration was induced in the Jovian planets by gravitational interactions with planetesimals leaking out of a ∼35 M_⊕ planetesimal disk residing between ∼16 AU and 30 AU (i.e., known as the primordial trans-planetary disk). Eventually, after a delay of ∼600 Myr (∼3.9 Ga), Jupiter and Saturn crossed a mutual mean motion resonance. This event triggered a global instability that led to a reorganization of the outer solar system; planets moved and in some cases had close encounters with one another, existing small body reservoirs were depleted or eliminated, and new reservoirs were created in distinct locations.

Despite its radical nature, the Nice model has been successfully used to deal with several long-standing solar system dynamics problems. It can quantitatively explain the orbits of the Jovian planets (Tsiganis et al. 2005), a partial clearing of the main Belt via sweeping resonances (Levison et al. 2001; Minton & Malhotra 2009), and the likely occurrence of a terminal cataclysm on the Moon and other terrestrial planets (Gomes et al. 2005; Strom et al. 2005). The timing of the Nice model would then be linked to the formation time of late forming lunar basins like Serentatis and Imbrium ∼3.9 Ga (e.g., Stöffler & Ryder 2001; Bottke et al. 2007). As we will show below, however, the precise timing of the Nice model does not play a large role in determining our results.

Perhaps the most critical test of the Nice model is determining what happens to objects from the primordial trans-planetary disk. According to simulations, giant planet migration led to resonance migration and resonance-crossing events that allowed a small fraction of scattered disk objects to be captured within the outer main Belt, Hilda, Jupiter and Neptune Trojan regions, irregular satellites, and TNO regions (Morbidelli et al. 2005, 2009b; Tsiganis et al. 2005; Nesvorný et al. 2007; Levison et al. 2009; Nesvorný & Vokrouhlický 2009). Thus, the comet-like bodies in these populations all came from the same reservoir and thus should presumably have similar properties.

This prediction appears to hold true from a taxonomy standpoint. The observed small objects in these populations are mainly dormant comet-like objects (i.e., C-, D-, and P-type bodies). We cannot yet rule out the possibility, however, that majority of small planetesimals beyond 2.8 AU formed this way.

More compellingly, this scenario also seems to work from a SFD perspective. Morbidelli et al. (2009b) showed the SFD of the TNOs is similar to that of Jupiter's Trojans, while Levison et al. (2009) showed that the D/P-type objects in the outer main Belt, Hilda, Trojan populations could all have come from a source population with a SFD shaped like the present-day Trojan asteroids. Sheppard & Trujillo (2008) also find that the cumulative luminosity function of the Neptune Trojans may have a turnover at the same place as the Trojans and Kuiper Belt objects (i.e., if converted into diameter, the turnover occurs near D ∼ 100 km). Thus, it is unavoidable; if Nesvorný et al. (2007) are correct and the irregular satellites came from the primordial trans-planetary disk, they should have the same starting population as the Trojans (Figure 2).

Given all this, we will base our calculations on the Nice model framework and will test whether the irregular satellites could have come from an SFD that originally had the same shape as the Trojan asteroids.

As a caveat, we point out that some irregular satellite capture scenarios cannot yet be ruled out based on the modeling done below. It is possible that some of them, if proven true, could change our initial conditions (e.g., sweeping resonances may have slightly modified the nature of the orbital populations; Ćuk & Gladman 2006). At this time, however, we believe no other scenario is mature enough to test within our model.

Nesvorný et al. (2007) modeled satellite capture in several steps. In step 1, they created synthetic Nice model simulations that examined the time shortly after Jupiter and Saturn had migrated through the 2:1 mean motion resonance and initiated the violent reshuffling event. The initial conditions of the gas giants were taken from the Gomes et al. (2005) simulations, with Jupiter, Saturn, Uranus, and Neptune at 5.4 AU, 8.4 AU, 12.3 AU, and 18 AU. The primordial trans-planetary disk was filled with thousands of planetesimals between 21 AU and 35 AU and were given a variety of initial configurations. These systems were tracked for at least 130 Myr using the symplectic integrator SyMBA (Duncan et al. 1998), with close encounters between the gas giants recorded for later use.

In step 2, they sifted the results for planetary systems that resembled our own and then used recorded gas giant close encounter and planetesimal/planet orbital data to create a high-fidelity satellite capture model. First, the orbits of the planetesimals and planets from step 1 at the time of each gas giant close encounter were integrated slightly backward in time to a time just before the encounter. Next, they used the planetesimals to construct orbital distribution maps suitable for creating millions of new test bodies on planetesimal-like orbits. Then, the test bodies and planets were integrated forward in time all the way through the close encounter. Objects deflected by the close encounter into planet-bound orbits were integrated for stability. Those satellites that survived were then followed into the next gas giant encounter, where their orbits could change or they could even be stripped altogether from their primary. This process was repeated until all close encounters within each run were completed.

Overall, Nesvorný et al. (2007) found that planetary encounters can create irregular satellites around Saturn, Uranus, and Neptune with (a, e, i) distributions that are largely similar to the observed ones. A drawback, however, is that Jupiter does not generally participate in close encounters in the Nice model (Tsiganis et al. 2005; Gomes et al. 2005), such that only certain runs were capable of producing Jupiter's irregular satellites. Rather than representing a fundamental flaw in the work, we instead consider this to be an indication that we are still missing some important aspects in our understanding of the Nice model itself. For example, Brasser et al. (2009) find that close encounters between Jupiter and the ice giants are needed to prevent the terrestrial planets from obtaining high eccentricities via secular resonance sweeping. It is also possible that the gas giants were started in a different and even more compact configuration, with the instability triggered by a different resonance (see Morbidelli & Crida 2007; Morbidelli et al. 2007). Regardless, the success of the Nesvorný et al. model for Saturn, Uranus, and Neptune and the orbital similarity of Jupiter's irregular satellite system to the other systems implies to us that Jupiter's irregular satellites were produced by gas giant encounters. On this basis, we will assume that the irregular satellites of all of the planets can be modeled using the Nesvorný et al. (2007) template.

Our collisional modeling simulations will employ Boulder, a new code capable of simulating the collisional fragmentation of multiple planetesimal populations using a statistical particle-in-the-box approach. A full description of the code, how it was tested, and its application to both accretion and collisional evolution of the early asteroid Belt, can be found in Morbidelli et al. (2009a). Examples of its use for the outer main Belt, Hildas, Trojans, and primordial trans-planetary disk can be found in Levison et al. (2009). Boulder was constructed along the lines of comparable codes (e.g., Weidenschilling et al. 1997; Kenyon & Bromley 2001) and can be considered an updated and more flexible version of the well-tested collisional evolution and dynamical depletion model code CoDDEM used by Bottke et al. (2005a, 2005b) to simulate the history of the main Belt.

A major difference between Boulder and other codes like CoDDEM is that Boulder uses the results of smoothed particle hydrodynamics (SPH)/N-body impact experiments to model the fragment SFD produced in asteroid disruption events (Durda et al. 2004, 2007). This is a key factor because, as we will show, accurate descriptions of different kinds of impact events are needed to model the irregular satellite SFDs. The code's procedure for modeling an impact is as follows.

For a given impact between a projectile and a target object, the code computes the impact energy Q, the kinetic energy of the projectile per unit mass of the target, and the critical impact specific energy Q*_D, defined as the energy per unit target mass needed to disrupt the target and send 50% of its mass away at escape velocity (e.g., Davis et al. 2002). For reference, Q < Q*_D events correspond to cratering events, Q ≈ Q*_D events correspond to barely catastrophic disruption events, and Q > Q*_D events correspond to super-catastrophic disruption events.

Numerical hydrocode experiments show that the mass of the largest remnant after a collision follows a linear function of Q/Q*_D (Benz & Asphaug 1999), with the mass of the largest remnant (M_LR) produced by a given impact:

$$\begin{eqnarray} M_{\rm LR}(Q<Q^*_D) & = & \left[-{{1}\over {2}} \left({{Q}\over {Q^*_D}}-1\right)+{{1}\over {2}}\right] M_T, \nonumber\\ M_{\rm LR}(Q>Q^*_D) & = & \left[-0.35 \left({{Q}\over {Q^*_D}}-1\right)+{{1}\over {2}}\right] M_T, \end{eqnarray} \tag{ 1 }$$

where M_T is the target mass.

To choose the fragment SFD ejected from the collision, we take advantage of data derived from the numerical hydrocode impact experiments of Durda et al. (2004, 2007). As part of a project to study asteroid satellite formation, Durda et al. (2004) used smooth particle hydrodynamic (SPH) codes coupled with N-body codes to perform 160 numerical impact experiments where they tracked D = 10–46 km projectiles slamming into D = 100 km basaltic spheres at a wide range of impact speeds (2.5–7 km s⁻¹) and impact angles (15°–75°). Durda et al. (2007) reported that the SFDs produced by these simulations, when scaled to the appropriate parent body size, were a good match to the SFD of the largest members of many observed asteroid families (i.e., the largest bodies had not been seriously affected by comminution since their formation; Bottke et al. 2005a, 2005b).

To include these results in Boulder, the mass of the largest fragment and the slope of the power-law SFD for each of the Durda et al. experiments were tabulated as a function of the ratio Q/Q*_D. Empirical fits to the experimental data indicated that the mass of the largest fragment (M_LF) and slope of the cumulative power-law size distribution of the fragments (q) could be written as

$$\begin{eqnarray} M_{\rm LF} &=& 8\times 10^{-3}\left[{{Q}\over {Q^*_D}}\exp ^{-\big({{Q}\over {4 Q^*_D}}\big)^2}\right]M_{\rm tot}, \nonumber\\ q & = & -\!10+7\left({{Q}\over {Q^*_D}}\right)^{0.4}e^{-{{Q}\over {7 Q^*_D}}} \end{eqnarray} \tag{ 2 }$$

with M_tot being the combined mass of the projectile and target bodies.

Using these equations, Boulder selects a largest remnant, a largest fragment, and the exponent of the power-law fragment SFD. This allows it to accurately treat both cratering and super-catastrophic disruption events in a realistic manner.

Note that in some extreme cases, such as like-sized bodies smashing into one another, M_LF > M_LR. This describes a highly energetic super-catastrophic event capable of pulverizing both the projectile and target bodies. For the runs described here, we assume such impact events create fragments smaller than our resolution limit (D > 0.1 km). Accordingly, we place the total mass of the projectile and target into the code's "trash" bin.

Here, we describe the parameters needed to run our Boulder simulations. Using the Nesvorný et al. (2007) model as our foundation, we tracked model prograde and retrograde irregular satellite populations captured at Jupiter, Saturn, and Uranus, the planets with sufficient numbers of known satellites to constrain our models, over 3.9 Gyr of collisional evolution (i.e., Neptune is left out). We also assumed that the remnants of the primordial trans-planetary disk passing through these irregular satellite systems were capable of striking and disrupting irregular satellites, though we did not model the collisional evolution of the primordial trans-planetary disk for reasons that will be described below.

Several parameters are needed as input into the code: (1) the starting SFDs of the prograde and retrograde irregular satellites, (2) the SFDs and dynamical decay function for the primordial trans-planetary disk, (3) the collision probabilities and impact velocities between the populations, and (4) the disruption scaling laws and bulk densities for the objects. The code follows the evolution of the irregular satellite SFDs while holding most of these parameters constant. An exception to this is that we forced the remnants of the primordial trans-planetary disk to dynamically decay over time according to the numerical results described in Nesvorný et al. (2007) and Nesvorný & Vokrouhlický (2009).

The values of these parameters for each of our populations are discussed in the subsections below.

3.4.1. Initial Size Frequency Distributions of the Irregular Satellites

The shape of the irregular satellites' SFD immediately after capture is assumed to have the same shape as the Trojan asteroid's SFD. The cumulative number N(>D) of known Trojans (from the "astorb.dat" database described above) can be fit by a broken power law with q = −5.5 for D > 100 km and q = −1.8 for D < 100 km, provided the Trojans have a mean albedo of 0.04 (Figure 2).

As a baseline to normalize the Trojan SFD for the irregular satellite populations, we look to the numerical simulations provided by Nesvorný et al. (2007). They predicted that the capture efficiency of irregular satellites from the primordial trans-planetary disk to Uranus and Neptune was (2.7–5.4) × 10⁻⁷. Thus, if the primordial trans-planetary disk had a mass at 3.9 Ga of 35 M_⊕, the mass captured at Uranus and Neptune should have been ∼0.001 lunar masses. This is roughly equivalent to the size of the present-day Trojan population, provided the bulk density of the objects was 1 g cm⁻³. This is consistent with the fact that the largest Trojan, (624) Hektor (D = 270 km; Storrs et al. 2005), is only modestly larger than Phoebe (D = 240 km), Himalia (D = 170 km), and Sycorax (D = 150 km) (Figure 1). Nereid at Neptune, which is D = 340 km, is plausible as well.

The source SFDs used to create our model prograde and retrograde satellite SFDs, defined as populations A and B, are shown in Figure 3. The shape of population A is the Trojan SFD described above, namely, a broken power law with q = −5.5 for D > 100 km and q = −1.8 for D < 100 km. We assumed our objects had bulk densities of 1 g cm⁻³. The starting mass of the SFD was fixed to values near 0.001–0.01 lunar masses by assuming the cumulative number N of D > 250 km objects was 0.1, 0.3, 1.0, and 3.0 (nominally 0.0008, 0.002, 0.008, and 0.02 lunar masses for D > 0.1 km). These values allow us to check whether larger starting populations can be ruled in or out.

Zoom In Zoom Out Reset image size

Population B mimics population A but has a "foot" attached, namely, a shallow branch with q = −2 that starts at D > 250 km. From a probabilistic standpoint, it allows larger objects to be captured into the prograde or retrograde populations. The normalization of population B is the same as above but the amount of mass captured is slightly higher (0.0009, 0.003, 0.009, and 0.03 lunar masses for D > 0.1 km). According to Morbidelli et al. (2009b), the foot might represent objects from the inner component of the primordial trans-planetary disk, which may have been more shallower-sloped than the outer component. Observational evidence for a foot may exist in the Kuiper Belt and other populations with large capture/delivery efficiencies (e.g., perhaps the outer main Belt). This issue will be discussed further in the discussion section.

To make our starting conditions as realistic as possible, we also added stochastic elements to mimic the capture process for each trial case of Boulder.

• 1.  
  The idealized populations A and B were divided into logarithmic intervals dlog D = 0.1. Random deviates were then used to select objects for the source SFD. In other words, if a given size bin in population A contained 0.1 objects, there would be a 10% chance that an object of that size would get into the source SFD. This accounts for the fact that the capture of the largest satellites should be probabilistic. Our testbed runs indicate a stochastic element in the capture process could potentially explain the diameter difference between the largest irregular satellites in each system (i.e., Himalia, Phoebe, Sycorax, and Nereid).
• 2.  
  We defined the parameter f_split, the fraction of satellites going into prograde (and retrograde) populations from the source SFD. Using f_split and random deviates, we assigned objects in the source SFD to the prograde or retrograde populations.

For the irregular satellites orbiting Jupiter and Saturn, we tested f_split = 0.4–0.6, while for those at Uranus, we tested 0.3–0.7 (see below). In the Nesvorný et al. (2007) simulations, captured bodies often had f_split ≈ 0.3–0.5, but values of 0.6–0.7 were also found. Our computational work took advantage of the fact that the collision probabilities and impact velocities for prograde–prograde and retrograde–retrograde collisions were similar to one another (see the following section). This means that a trial case for f_split = 0.3, if the populations are switched, can also be used to examine f_split = 0.7.

In our production runs, we used the following procedure. (1) We chose a value of f_split ⩽ 0.5; we call this α. (2) We tested our prograde and retrograde SFDs against the known prograde and retrograde satellites, respectively. (3) We "switched" prograde and retrograde SFDs, which are the equivalent of testing f_split = 1 − α. (4) We ran the same tests as (2).

3.4.2. Collision Probabilities and Impact Velocities

3.4.3. The Primordial Trans-planetary Disk Population

The other population that can collide with the irregular satellite SFDs is the surviving remnant of the primordial trans-planetary disk. This population readily decays as the gas giant scatters away its members, but not so quickly that it can be ignored. For added realism, we included it in our model.

Figure 4 shows two estimates for this population just prior to Jupiter and Saturn crossing into a mutual mean motion resonance, one containing ≈20 M_⊕ and one with ≈40 M_⊕ (e.g., Tsiganis et al. 2005). As before, we assume that at the moment the instability started, the disk SFD had the same shape as the current Trojan SFD. These populations were normalized by assuming there were 7.5 × 10⁵ and 1.5 × 10⁶ objects with D > 200 km, respectively. Note that most of the mass in each SFD is near the inflection point at D = 100 km, such that the large body end of the SFD does not play a meaningful factor in the evolution of the irregular satellites.

Zoom In Zoom Out Reset image size

To include this population in Boulder, we must account for the fact that the primordial trans-planetary disk is scattered and cleared during the Nice model simulations. This means that the number capable of hitting the irregular satellites decreases dramatically with time. At the same time, we also need to couple this information to the P_i and V_imp values computed between the remnant disk population and the irregular satellite populations.

We characterized the scattering of the disk using new runs based on the Nice model template described in Nesvorný et al. (2007), Nesvorný & Vokrouhlický (2009), and in Section 3.1. Here we started 27,000 disk particles (rather than 7000) and integrated them out to 1 Gyr (as opposed to 130 Myr). Particles were removed from the system if they struck a planet, the Sun, or were thrown out of the outer solar system. For computational expediency, we assumed the same decay curve for objects crossing the orbits of Jupiter and Saturn.

Next, using these data and the techniques described in Charnoz et al. (2009), we computed P_i and V_imp values between disk particles passing near the gas giants and the irregular satellite populations. A running mean of these values was used to eliminate jitter. In some cases, we also interpolated through the values when it was clear that it could be fit to a line segment (e.g., see Bottke et al. 2005b). Beyond 600 Myr, we assumed that the remnant populations shown in the figure were a reasonable "stand-in" for the refugees from the Kuiper Belt and scattered disk that occasionally hit the irregular satellites. Accordingly, we assumed that the impacting disk population did not change in a meaningful way beyond this time.

Our results, which are shown in Figure 5, show that the primordial trans-planetary disk population decays rapidly, with a slower drop off for objects near Uranus than those near Jupiter and Saturn. This occurs because Uranus (and Neptune) are embedded in the disk while Jupiter and Saturn are effectively decoupled from the population. Test results also indicate that the decay is fast enough that the irregular satellites are only seriously affected by disk–satellite collisions during two intervals: the first few tens of Myr of the simulation, where the population is still large, and last several Gyr of the simulation, where the irregular satellite populations have become so decimated that the sporadic impacts of disk bodies actually have some effect.

Zoom In Zoom Out Reset image size

Finally, we find that the impact velocities V_imp between the remnant disk population and the irregular satellite populations do not change very much throughout the simulation. The values used for impacts on the irregular satellites at Jupiter, Saturn, and Uranus are 7.0 km s⁻¹, 4.7 km s⁻¹, and 3.0 km s⁻¹, respectively (Table 2). The higher impact velocities found for impacts at Jupiter and Saturn compensate somewhat for their more rapid decay curves (Figure 5).

3.4.4. Disruption Scaling Laws and Bulk Densities for Irregular Satellites

For a Q*_D function applicable to irregular satellites, we turn to Levison et al. (2009), who recently used Boulder to model the collisional evolution of icy planetesimals in the primordial trans-planetary disk that were scattered into the outer main Belt, Hilda, and Trojan asteroid populations. If the irregular satellites came from that same source population, it is reasonable to assume that the irregular satellites have the same disruption properties as the C-, D-, and P-type objects modeled in Levison et al. (2009).

In Levison et al. (2009), the Q*_D function was assumed to split the difference between SPH impact experiments of Benz & Asphaug (1999), who used a strong formulation for ice, and those of Leinhardt & Stewart (2009), who used the finite volume shock physics code CTH to perform simulations into what they describe as weak ice. To do this, Levison et al. (2009) examined what happened when the Benz & Asphaug (1999) strong ice Q*_D function was divided by a factor, f_Q. They found that the best match to the observed populations came from using f_Q = 3, 5, and 8 (see Figure 6). Note that because we sampled a broad section of parameter space, we chose not to include still more complicating factors (e.g., Q*_D varies with impact velocity, etc.).

Zoom In Zoom Out Reset image size

The bulk densities (ρ) of typical irregular satellites are largely unknown; direct estimates are only available for the largest Jovian and Saturnian satellites. Phoebe, which the Cassini spacecraft encountered in 2004, has ρ = 1.634 g cm⁻³ (Jacobson et al. 2006). Himalia, which has produced measurable gravitational perturbations on nearby satellite Elara, appears to have a comparable or perhaps much higher ρ value, though much depends on its exact shape and size (Emelyanov 2005).

By assuming the irregular satellites came from the primordial trans-planetary Belt, however, we can assume that these objects have comparable ρ values to other objects that came from the disk, namely comets, Kuiper Belt objects, and Trojan asteroids. A survey of the literature produces ρ values of between 0.46 (Centaur 2002 CR₄₆; Noll et al. 2006) and 2.5 g cm⁻³ (Binary Trojan (624) Hektor; Lacerda & Jewitt 2007; see also a recent review by Weissman et al. 2004). Taken together, we choose to make the same assumption made by Levison et al. (2009), namely we will give our objects comet-like bulk densities of 1 g cm⁻³. This ρ value was selected to roughly split the difference between the extremes of the above values.

3.4.5. Caveats

One of the most complex issues in this project is determining when a good match exists between our parameter-dependent model SFDs and the observed irregular satellite SFDs. While there is a strong temptation to compare both "by eye" and simply report the results, this becomes impractical when we are dealing with hundreds of runs that contain millions of output SFDs. Moreover, the results would be highly subjective and would differ between observers.

For this reason, we decided that it would be useful to introduce a scoring (target) function that helps us define the comparison of the modeled and observed SFDs quantitatively. In our first attempts, we tried to apply statistical chi-square-like tests (e.g., Press et al. 1992) to the data. Unfortunately, we found them difficult to apply in our situation because (1) the observed SFDs had poorly defined uncertainties and arbitrary bin sizes, (2) the model SFDs output from Boulder were assigned ever-changing bin center locations as collisional evolution took place, and (3) the uncertainties in the sizes of the observed objects were, in most cases, estimated from an albedo assumption rather than by a rigorous method.

In this situation, we decided to use a much simpler tool where the quantitative threshold of "quality" in the model-to-data comparison was obtained using a predefined metric rather than an exact statistical method. The advantage to this method was that it was fast, easy to evaluate, and did not directly use the unreliable uncertainties in the satellite sizes. The disadvantage is that we cannot assign a statistical level of confidence to our results. With that cautionary note, we describe our method below.

We assume that D^(O)_i, i = 1, ..., N, defines a vector of sizes for the observed population of irregular satellites, while D^(C)_i is the same for the modeled population. No bin has more than one member; if the modeled population has a size bin occupied by several objects, we spread out the respective number of copies in D^(C)_i. The quality of the match between the observed and modeled populations is then defined as the scoring function:

$$\begin{eqnarray} &&S = \frac{1}{N} \sum _{i=1}^{N} \left|{\rm log}\left(\frac{D^{(C)}_i}{D^{(O)}_i} \right)\right|\;. \end{eqnarray} \tag{ 3 }$$

The S-score is evaluated separately for both the prograde and retrograde irregular satellite populations. The modeled population typically contains many more objects than the observed one, N in our case, and we consider the first N largest bodies in the modeled population to evaluate the score function (3). In some cases, though, one of the irregular populations might become decimated by the other, such that the modeled population could contain fewer objects than the observed one. In that case, N in Equation (3) becomes the number of modeled objects. Note that if N ever equals zero, S is assigned an arbitrarily large value indicating a poor match.

As an example of how to use S, consider a modeled SFD whose largest objects get within a factor of 2 in diameter to the members of the observed SFD. In that circumstance, S drops to a value of 0.3 or smaller. Given all other uncertainties and biases, we consider this level of agreement to be reasonable enough that we will assume S ≃ 0.3 is the acceptable threshold of a "good" match. It also produces results that are consistent with our "by eye" evaluations.

We start our description of our model results by showing off a successful trial case. Figure 7 shows eight snapshots from a trial case for Jupiter's irregular satellites with input parameters f_Q = 3, f_split = 0.4, N(D > 250 km) = 0.1, q = −5.5, and M_disk = 40 M_⊕. The initial conditions are shown in the time t = 0 Myr panel (top left). Here f_split specifically corresponds to 40% of the population going into the prograde population. The model retrograde and prograde SFDs are represented by green and magenta dots, respectively, while those of the observed retrograde and prograde SFDs are given by the blue and red dots, respectively.

Zoom In Zoom Out Reset image size

The top end of the model SFDs shows off the stochastic nature of the capture process. For example, for f_split = 0.4, the model prograde objects ended up with 1 D = 180 km object, 4 D = 110 km objects, and 11 D = 90 km objects, while the retrograde object obtained 2 D = 140 km objects, 6 D = 110 km objects, and 18 D = 90 km objects. Thus, while the fraction of objects on prograde versus retrograde orbits is what we would expect (0.38), the mass distribution can be a variable.

The first two time steps, t = 1 Myr and 10 Myr, show that the two populations grind away fastest when they are largest. For example, at 10 Myr, each SFD has decreased by factors of 3–10 over nearly their entire span. We find the shape of the SFDs for D < 10 km largely depends on the nature of the breakup events occurring at larger sizes. Particularly large super-catastrophic disruption events often dominate the population at smaller sizes for extended periods. This allows the prograde and retrograde SFDs for D < 5 km to undergo sudden and sometimes radical changes. In fact, it is common to see the SFDs jump to steeper slopes in the aftermath of a large collision event and then slowly retreat to shallow slopes as widespread grinding beats the new fragment population back down. Both prograde and retrograde scores S for 1–10 Myr are >0.5, indicating that we are still far from a satisfying match to the observational data. At t = 40 Myr, however, the SFDs begin to take on shapes that are similar, in many ways, to the observed populations.

The endgame of the satellite evolution simulation begins near t = 500 Myr. Here the matches between model and observations approach S ≈ 0.3, the minimum threshold value needed for a good fit. At this point, both populations are beaten up and battered. The large body populations have become so depleted that, at any given time, most small fragments are produced by cratering events. Sporadic large-scale catastrophic disruption events, however, can and do dominate the small body populations from time to time. This means the shapes of the SFDs for D < 5 km wiggle and wag up and down for the next several billion years. Occasionally, the match with the observed SFD becomes quite good, as seen for t = 3883 Myr. This cannot last, though, and collisional evolution over several tens of Myr is enough to degrade the fit until the next big disruption event restarts the process. In other words, all of this has happened before and will happen again.

Table 3 shows the parameter sets for the top five runs that produced good fits (S < 0.3) to the observed prograde and retrograde over the last 500 Myr. These parameters produce good fits with the data approximately 50% of the time over the tested interval. We consider this remarkable when one considers all of the possible ways our simulations could have run into trouble. The common theme in the top five runs appears to be f_Q = 3 and an input SFD without a foot (q = −5.5) that was normalized for N(D > 250 km) = 0.1. To some degree, this would argue against input SFDs that contain more large bodies.

We found that using the top five runs alone to analyze trends across 128 runs was generally an unsatisfying way to sift our data. For that reason, we decided to augment our analysis by combining our results across all runs. Here we computed the average number of good fits found over the last 500 Myr in each run and combined these data across our five input parameters. Our net results are shown in Figures 8 and 9.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 8, our Jupiter runs show preferences for f_Q ⩾ 3, with 3 being the favorite, and for smaller values of N(D > 250 km), with 0.1 being the favorite. In Figure 9, the trends are less distinct, but f_split values of 0.5 and q values of −5.5 are generally favored. There does not appear to be any trend for M_disk. These values are in generally good agreement with the trends seen for the top five runs from Table 3, giving us confidence that both approaches yield useful results.

Saturn's irregular satellites were tested against the same parameter suite as that used for Jupiter's runs. Thus, in terms of collisional evolution, the only model differences between the satellite systems should be Saturn's P_i values, which vary slightly from Jupiter's (Table 1), and Saturn's significantly lower V_imp values (Table 2). Figure 10 shows eight snapshots from a trial case for Saturn's irregular satellites. Our initial conditions are shown in the time t = 0 Myr panel. The color code for the different populations is the same as in Figure 7.

Zoom In Zoom Out Reset image size

To demonstrate that there are many ways to get reasonable fits, we show a trial case here with somewhat different parameters than Figure 7: f_Q = 8, f_split = 0.6, N(D > 250 km) = 0.3, q = −5.5, and M_disk = 40 M_⊕. Here an important difference from Figure 7 is that 60% of the starting population goes into the model prograde population rather than the model retrograde one. A larger starting prograde population can be useful because the observed prograde satellites have more D = 10–30 km bodies than the observed retrograde bodies.

While the shapes of the SFDs through the first 500 Myr in Figure 10 are reminiscent of those seen in Figure 7, the factor of 2 decrease in V_imp values at Saturn when compared to Jupiter mean typical collisions are less energetic. The consequence is that the SFDs grind away more slowly, which helps explain the larger abundance of D = 10–30 km bodies present in the prograde population at t = 500 Myr. Low impact energies also mean cratering events are more common than before. As seen in the t = 500 Myr, 1000 Myr, and 1500 Myr time steps, the fragment tail of this ejecta shows up for D < 1–7 km and will actually wither away on timescales of tens of millions of years until a stochastic disruption or cratering event on one of the larger remaining satellites freshens it up with new ejecta. The pattern then begins again. By animating the timesteps, we find the D < 7 km SFDs for both populations wave up and down again and again over billions of years, much like a ringmaster cracking a whip.

Our best score for this case comes at t = 3651 Myr. The S values for the prograde and retrograde SFDs are 0.15 and 0.07, quite good, but with some limitations that we will discuss below. Here, ejecta evolving from a cratering event on the largest remaining retrograde satellite reproduces the observed steep retrograde population near D ≈ 7 km. This excellent match does not last more than a few tens of Myr, but comparable fits are found at several later times in the simulation via the "cracking whip" process described in the previous section.

The main deficiencies in this trial case are twofold. First, our model was unable to reproduce the substantial difference in diameter between Phoebe and the second largest retrograde irregular satellite (Figure 10). Second, our model prograde satellites for D > 30 km are modestly larger than the observed ones. Similar problems show up in other S < 0.3 trial cases.

The culprit is likely to be our initial conditions. Using our selected capture efficiencies, source populations A and B (Figure 3) nearly always produce a continuum of retrograde objects, which makes it difficult to produce significant differences in size between Phoebe and the second-largest retrograde object. Similarly, if all captured objects follow the same Q*_D disruption function, it is difficult to eliminate large numbers of middle-sized objects without affecting Phoebe as well.

The solution to the Phoebe problem is unclear, but we hypothesize that it involves factors that we have yet to consider in our modeling work. For example, Phoebe, classified as a C-type object with water ice and CO₂ on its surface (Johnson & Lunine 2005) as well as a bulk density of 1.634 g cm⁻³ (Jacobson et al. 2006), may be harder to disrupt than typical D- and P-type irregular satellites. In this scenario, the smaller objects would disrupt/erode while Phoebe would be left more or less intact. For a second possibility, it could be that the capture efficiency at Saturn was lower than at the other gas giants, and Phoebe's capture was a fluke at the ∼10% level. This would allow fewer middle-sized objects to be captured. A third possibility is that Phoebe was captured onto its current orbit in a different manner than the other irregular satellites (e.g., perhaps the binary capture mechanism that allowed Triton to be captured at Neptune; Agnor & Hamilton 2006; Vokrouhlický et al. 2008). The problem is finding a mechanism that works during or after the Nice model and allows Phoebe to reach its current orbit (Vokrouhlický et al. 2008). Note that these scenarios do not have to be mutually exclusive; perhaps several of them played a role at some level.

The top five irregular satellite runs for Saturn are somewhat more varied than those for Jupiter (Table 3). Good fits between model and data occurred 25%–40% of the time over the last 500 Myr of the simulations. We believe this fractional value, while lower than Jupiter's by a factor of 1.4, is still highly encouraging, particularly when one considers the Phoebe issue above as well as how difficult it is to fit the always-evolving small end of the SFDs for long intervals. The combined results shown in Figures 8 and 9 indicate that Saturn's runs follow the same overall trends as Jupiter's runs. Moderately favored parameters are f_Q values of 3–8, N(D > 250 km) values of 0.1–1.0, variants of population A (q values of −5.5), f_split values of 0.5, and M_disk values of 40 M_⊕. These results are consistent with the trends seen in the top five runs, where normalization values of N(D > 250 km) = 0.1, population A input SFDs (q = −5.5), and larger disks (M_disk = 40 M_⊕) are favored.

The irregular satellites of Uranus have SFDs that are somewhat different than those examined thus far (Figure 1). The known prograde population is comprised of a single D = 20 km object, while the retrograde SFD is a fairly smooth continuum made up 10 D > 10 km objects and 1 D > 100 km object. This order of magnitude mismatch implies something interesting went on in this system. A partial cause may be the low V_imp values for satellite collisions, which range between 1 km s⁻¹ and 2 km s⁻¹ (Table 2). This is unlikely to be the whole story, though, because Saturn's values are only modestly higher.

Figure 11 shows eight snapshots from a trial case with a very good fit to Uranus's irregular satellites. The input parameters used were f_Q = 8, f_split = 0.3, N(D > 250 km) = 1.0, q = −2.0, and M_disk = 20 M_⊕. The f_split = 0.3 value used here is lower than those used for the Jupiter/Saturn cases but justifiable based on the dynamical capture runs of Nesvorný et al. (2007). We were driven to this value by multiple test runs using f_split values of 0.4 or 0.5; neither choice produced more than a marginal number of successful fits among our many trial cases.

Zoom In Zoom Out Reset image size

We found that in this trial case, the "head start" given to the retrograde SFD over the prograde one allowed it to immediately dominate the collisional evolution of the system. This leads to what we refer to as a runaway collisional process. The prograde SFD quickly becomes decimated, so much so that its fate is controlled by the evolution of the retrograde SFD. This leads to the endgame situation seen at t = 500 Myr, where only a handful of large prograde objects are left standing.

After 500 Myr, the remaining retrograde SFD (and external Kuiper Belt/scattered disk SFDs) are depleted enough that collisional erosion among the prograde objects is dominated by cratering rather than catastrophic disruption events. Bit by bit, these relatively small impacts, over billions of years, grind down the surviving prograde objects. Eventually, at t = 3888 Myr, we end up with an excellent fit between the model and observed SFDs (i.e., S = 0.0005 for prograde objects and S = 0.09 for retrograde objects). The high quality of our fit, and the fact that the prograde population only includes one D ≈ 20 km object, allows us to make the following predictions: (1) by the end of the capture phase, the retrograde population was substantially larger than the prograde one, and (2) the initial difference between the SFDs was greatly exacerbated over billions of years by a runaway collisional process that produced the near-elimination of the prograde population.

The odds of getting a runaway collisional outcome are reasonable but not overwhelming. From our top five irregular satellite runs for Uranus (Table 3), we find success rates of 10%–30% when we checked the last 500 Myr of our simulations. These results also favor larger initial populations (i.e., N(D > 250 km) values of 1.0), population A-type SFDs (i.e., q values of −5.5), and satellites that are easy to disrupt (i.e., f_Q values of 8). Using the combined results shown in Figures 8 and 9, we obtain similar results. There are relatively strong preferences for f_Q values near 8 and N(D > 250 km) values near 1.0. The preferences are more moderate for M_disk values near 20 M_⊕ and q values near −5.5.

An interesting and potentially exciting implication of our model results is that the same collisional cascades that demolish C-, D- and P-type irregular satellites should also create enormous amounts of dark dust. For example, using the starting SFDs shown in Figures 7, 10, and 11, we predict that ∼99% of the mass of the irregular satellite populations at Jupiter, Saturn, and Uranus were lost to comminution over 4 Gyr. Thus, in each case, we have ∼0.001 lunar masses or more of dust to deal with. What happened to this material?

Small particles in the solar system are strongly affected by solar radiation pressure forces and Poynting–Robertson (P–R) drag (for a review of these processes, see Burns et al. 1979). The importance of these mechanisms depends on the size, shape, and composition of the particles. Burns et al. (1979) show that micron-sized dust particles on planetocentric orbits have their eccentricities rapidly pumped up by radiation pressure forces, while P–R drag sends objects up to hundreds of microns in size inward toward the planet. If started from present-day irregular satellite orbits, these particles will eventually achieve crossing orbits with the outermost regular satellites, where they have the potential to strike them.

The collision probabilities between irregular satellite dust and the regular satellites have only been calculated for a limited number of examples, probably because few irregular satellites were known until recently. The most prominent cases in the literature involve dust particles from Phoebe evolving inward to strike Iapetus (D = 1470 km), the outermost regular moon of Saturn with a dark leading side and a bright trailing side. Burns et al. (1996) estimated that D = 20 μm sized particles from Phoebe have a 70% chance of striking Iapetus. Those that get past Iapetus have a 60% chance of striking Hyperion (D = 270 km). The rest are lost to Titan (D = 5150 km), which protects the inner moons from significant dust contamination. Burns et al. report that smaller and faster-moving particles have lower impact probabilities with Iapetus and Hyperion, but few still get past Titan. The equations provided by Burns et al. (1979) indicate comparable impact probabilities should exist between dust particles from irregular satellites and the outermost regular satellites of Jupiter and Uranus. Part of the reason has to do with target size; at Jupiter, Callisto and Ganymede (D = 4820 and 5260 km, respectively) are similar in size to Titan and are much larger targets than Iapetus, while at Uranus, Oberon and Titania (D = 1520 and 1580 km, respectively) are only slightly larger than Iapetus.

The inescapable conclusion from all of this is that the outermost regular satellites of the gas giants should have been blanketed by large amounts of dust produced by irregular satellite comminution. Most of this dark material would have landed within the first few hundreds of Myr after the capture of the irregular satellites. Thus, if the timing of the Nice model is linked with the ages of late-forming lunar basins like Serentatis and Imbrium, which have ages of ∼3.9 Ga (Stöffler & Ryder 2001), most of this dark dust landed on the satellites ∼3.5–3.9 Gyr ago.

Many authors have either suggested or explored variants of this intriguing scenario over the last several decades. While this is by no means a complete list, the interested readers are encouraged to investigate the following papers and references therein (e.g., dynamical links between Phoebe dust and the leading face of Iapetus; Soter 1974; Burns et al. 1996; dark non-icy components with spectroscopic properties similar to C-, D-, and P-type asteroids have been identified on many outermost regular satellites; Cruikshank et al. 1983; Buratti & Mosher 1991, 1995; Buratti et al. 2002; meteoroid contamination from the irregular satellites may be the source of the albedo differences on Ganymede and Callisto; Pollack et al. 1978; Johnson et al. 1983; Bell et al. 1985; see also McKinnon & Parmentier 1986). Presenting all of this material is beyond the scope of this paper, so below we focus on several interesting regular satellites.

6.1.1. Callisto

6.1.2. Iapetus

6.1.3. Titan

6.1.4. Uranian Satellites

6.1.5. Summary and Discussion

Another intriguing implication of this work concerns the crater records of the outermost regular satellites. Most groups have tried to interpret and/or model these crater populations using heliocentric impactors (e.g., refugees from the primordial trans-planetary disk, ecliptic comet, and nearly isotropic comet populations) or planetocentric ones (e.g., planetesimals leftover from satellite accretion; fragments from the catastrophic disruption of a small moon; ejecta from cratering events on the regular satellites; refugees from the asteroid Belt that somehow enter into planetocentric orbits prior to satellite impact; e.g., Strom 1987; Chapman & McKinnon 1986; McKinnon et al. 1991; Levison et al. 2000; Zahnle et al. 1998, 2001, 2003, 2008; Schmedemann et al. 2009; Charnoz et al. 2009). No group has yet considered, however, whether irregular satellites constitute an important source of impactors for the regular satellites.

Immediately after capture, the irregular satellites, whose population should rival the size of the present-day Trojan asteroids (Figure 2), comprise a nearly isotropic distribution of impactors spanning much of the space within each giant planet's Hill sphere (Nesvorný et al. 2007). Several sub-populations within this group are capable of striking the regular satellites. For example, objects embedded in the Kozai resonance will move into the regular satellite region on 10⁴ yr timescales (Carruba et al. 2002; Nesvorný et al. 2003). Another fraction, captured outside the regular satellite-crossing zone but inside the stable irregular satellite zone, will move more slowly onto collision trajectories with the regular satellites (e.g., Carruba et al. 2002; Nesvorný et al. 2003; Haghighipour & Jewitt 2008). Finally, impacts on irregular satellites from irregular or heliocentric impactors may inject new fragments onto unstable orbits where they can then go on to strike the regular satellites. The latter source may create some craters on the regular satellites today.

Irregular satellites capable of striking the regular satellites should have SFDs that change significantly over time via collisional and dynamical evolution (e.g., Figures 7, 10, and 11). This means that if the irregular satellites are important sources of craters, the chronology of events on the outer planet satellites can only be accurately interpreted using a coupled collisional and dynamical evolution model. While a full test of this scenario is beyond the scope of this paper, we decided that, as a proof of concept test, it would be interesting to see whether the shapes of the crater SFDs found on some outer planet satellites are similar to "snapshots in time" of the SFDs found in our best-fit irregular satellite trial cases (Figures 7, 10, and 11).

Our results for Callisto (at Jupiter), Iapetus (at Saturn), and Titania (at Uranus) are shown in Figure 12. The crater counts and error bars were taken from the literature: Callisto (Strom et al. 1981; McKinnon et al. 1991), Iapetus (Jaumann & Neukum 2009), and Titania (Strom 1987). The crater data were converted into cumulative projectile SFDs using the scaling laws described by Zahnle et al. (2001). The impact velocities between regular and irregular satellites were computed using the code of Bottke et al. (1994), with our accepted values essentially the same as those reported in Zahnle et al. (2001) for heliocentric impactors. Note that because the absolute scaling is unknown, the transformed crater SFDs were slid up and down and fit "by eye" to the dominant irregular satellite SFD (usually retrograde) every 1 Myr over 4 Gyr of evolution. The best overall fit between model and data for each case is shown in Figure 12.

Zoom In Zoom Out Reset image size

Overall, the matches are surprisingly good, especially when one considers our unknowns and that we are attempting to compare crater densities, which are integrated over the age of a surface, to snapshots in time from our Boulder trial cases. Thus, if the fits are not flukes, they may be telling us about the evolution of an irregular satellite population captured on unstable orbits (Callisto, Iapetus) or how stochastic breakup events at particular epochs dominate the crater records of some satellites (Titania). They also present us with a tantalizing possibility that the shape of a crater SFD alone may be diagnostic enough to estimate surface ages. On this basis, we analyze our model results below.

For Callisto, the best fit to the crater SFD comes at t = 25 Myr, where t is defined as the time after the initial capture of the irregular satellites (Figure 12). Interestingly, this age is broadly consistent with Levison et al. (2002), who claimed that the clearing of the primordial trans-planetary disk should have produced enough basin-forming impacts to melt Callisto's surface to considerable depth (see also Barr & Canup 2010). If true, our predicted surface age may represent the time Callisto's surface had stabilized enough to support impact structures again.

Iapetus may have the oldest non-saturated surface among the outer planet satellites. Using a projectile SFD for the primordial trans-planetary disk based on Bernstein et al. (2004), Charnoz et al. (2009) showed that heliocentric impactors during the Nice model may be a primary source of basins and D = 10 km craters on Iapetus's surface. The crater SFD for D < 10 km craters, however, is weakly "S"-shaped (Figure 12); our results indicate that such shapes are diagnostic of irregular satellite populations that have experienced a short interval (t ∼ 1 Myr) of intense collision evolution. At this time, we cannot say whether Iapetus's largest basins were produced by heliocentric or planetocentric impactors, but we find it probable that many smaller craters come from irregular satellites.

Titania is the largest moon of Uranus (D = 1578 km) but is comparable in size to Oberon (D = 1523 km) on the outside and is slightly larger than Umbriel (D = 1169 km) on the inside. Despite this, it is far fewer D > 100 km craters than either one (Strom 1987), probably because Titania was extensively resurfaced at an unknown time (Croft & Soderblom 1991). Our results indicate that time would be t ∼ 850 Myr or roughly 3 Ga (Figure 12). The source of geological activity on Titania at such relatively recent times is unknown.

We conclude by pointing out that the flux of projectiles and dark dust from the irregular satellites give us a new and potentially powerful tool for determining the chronology of events on the regular satellites. They also allow us to link these events to those happening elsewhere via the Nice model framework. By knitting all such events together, we hope to eventually create a tapestry that can accurately describe the history of the solar system.