Impact Rates in the Outer Solar System

Comets are icy objects that reach the inner solar system after leaving distant reservoirs beyond Neptune and are dynamically evolving onto elongated orbits with small perihelion distances (see Dones et al. 2015 and Kaib & Volk 2023 for recent reviews). Their activity, manifesting itself by the presence of a dust/gas coma and tail, is driven largely by solar heating and sublimation of various ices. Once they reach the inner solar system, comets are short-lived, implying that they must be resupplied from external reservoirs. Here we focus on the ecliptic comets (ECs, low- to moderate-inclination planet-crossing bodies in the region of the giant planets; see Levison & Duncan 1997 and Fraser et al. 2023 for further discussion) because (1) the population of ECs is relatively well characterized from observations and allows us to construct a realistic model (Section 2), and (2) ECs and their Centaur precursors are more important for impact cratering in the outer solar system than other types of comets (e.g., long-period comets), escaped asteroids, and Trojans (e.g., Zahnle et al. 1998; Levison et al. 2000; Zahnle et al. 2001; and Zahnle et al. 2003, hereafter Z03).

Levison & Duncan (1997, hereafter LD97) considered the origin and evolution of ECs. They assumed that the classical Kuiper Belt 30–50 au from the Sun was the main source of ECs. They showed that small Kuiper Belt objects (KBOs) reaching Neptune-crossing orbits can be scattered by encounters with the outer planets to small perihelion distances (q < 2.5 au), at which point they are assumed to become active and visible Jupiter-family comets (JFCs). The JFCs reaching q < 2.5 au for the first time have a narrow inclination distribution in the LD97 model because the orbits were assumed to start with low inclinations in the Kuiper Belt, and the inclinations changed little before the comets reached Jupiter-crossing orbits. ⁶ Levison & Duncan (1994) and LD97 pointed out that the inclination distribution of JFCs widens over time due to scattering encounters with Jupiter. LD97 found their best fit to the observed inclination distribution of JFCs (median ≃13°) when they assumed that ECs remain active for ≃12,000 yr (with a range of 3000–30,000 yr) after first reaching q < 2.5 au.

The escape of ECs from the classical Kuiper Belt is driven by slow chaotic processes in various orbital resonances with Neptune. Because these processes affect only part of the belt, with most orbits being stable for billions of years, comet delivery from the classical Kuiper Belt is inefficient (Nesvorný et al. 2017). Duncan & Levison (1997) suggested that the scattered disk, ⁷ whose prototype is (15874) 1996 TL66 (semimajor axis a = 134, q = 35 au; Luu et al. 1997), should be a more prolific source of ECs. This is because scattered disk objects (SDOs) can approach Neptune during their perihelion passages and be scattered by Neptune to orbits with shorter orbital periods, implying a faster loss rate for SDOs than for classical KBOs (Duncan et al. 2004). The population of SDOs is also inferred to be larger than that of classical KBOs (∼0.05 M_⊕ versus ∼0.01 M_⊕; Fraser et al. 2014; Nesvorný 2018), thus representing a large source for ECs.

Z03 used the spatial distribution of ECs from LD97 to determine impact rates in the outer solar system. They estimated that the current rate at which EC nuclei with diameters D > 1.5 km strike Jupiter is ${0.005}_{-0.003}^{+0.006}$ yr⁻¹. ⁸ Their rate for Jupiter was based on considerations such as the four known passages of JFCs within 3 R_J in the past 150 yr and models of the population of near-Earth JFCs from LD97. Z03 then calculated impact rates on Saturn, Uranus, and Neptune by scaling from the rate on Jupiter using LD97's model. Impact rates on the satellites of the giant planets were, in turn, calculated by Z03 with an Öpik-like formalism (Öpik 1951). ⁹ The main uncertainties in the impact rates given in Z03 arise from (1) the assumed source reservoir in LD97 and hence the spatial distribution of ECs (see discussion above), (2) our poor understanding of the sizes of the nuclei of active ECs, (3) the approximate nature of the impact flux calculation in LD97 (see Levison et al. 2000), and (4) the small number of direct impacts recorded on each planet in LD97. As for item 2, a common procedure is to infer the diameter D of a comet nucleus from its absolute magnitude H_T, which quantifies the total brightness of a comet, including coma, at a standard distance of 1 au from the Earth and Sun (e.g., Brasser & Morbidelli 2013). The nature of the relation between H_T and D is, however, uncertain because comets vary greatly in their levels of activity. ¹⁰

An important prerequisite for the interpretation of the cratering record in the outer solar system is to have a reliable model of ECs. ¹¹ We developed such a model in Nesvorný et al. (2017, hereafter N17). In N17, we performed end-to-end simulations in which cometary reservoirs, including the scattered disk and Oort cloud, were produced in the early solar system and evolved for 4.5 Gyr. The simulations included the effects of the four giant planets. The model was calibrated from the observed population of active ECs but includes other types of comets as well. We considered different scenarios for the duration of cometary activity, including a simple model in which comets remain active for N_p(q) perihelion passages with perihelion distance q < 2.5 au. To constrain N_p(2.5), we compared the orbital distribution and number of active comets produced in the model to observations. The observed distribution was well reproduced with N_p(2.5) ≃ 300–800 (see Figure 1). Here we adopt N_p(2.5) ≃ 600 as a reference value for comets with D ∼ 1 km (N17 inferred that N_p(2.5) should increase with the size of the nucleus). As the median orbital period of ECs is ≃8 yr (N17), our reference value corresponds to ≃4800 yr—a factor of 2.5 shorter than the nominal LD97 estimate. Ultimately, this is a consequence of new ECs (i.e., bodies reaching Jupiter-crossing orbit for the first time) having a wider inclination distribution in the N17 model because they start with larger inclinations in the scattered disk.

Zoom In Zoom Out Reset image size

The main uncertainty in the N17 model was related to the absolute calibration. They used Jupiter Trojans, whose size distribution is well characterized from observations down to D ∼ 3 km (Wong & Brown 2015; Yoshida & Terai 2017), for this purpose. This method relies on the assumption that the Trojan implantation efficiency from the original planetesimal disk is relatively well determined (Morbidelli et al. 2005; Nesvorný et al. 2013). We modeled the collisional evolution of Jupiter Trojans and SDOs and found that the size distribution changes after their implantation were insignificant (Nesvorný et al. 2018; Bottke et al. 2023). In addition, Pluto/Charon craters indicate that the size distribution slope of impactors with D > 1 km in the Kuiper Belt is similar to that of Jupiter Trojans (Singer et al. 2019). The use of the Trojan size distribution to model comets may therefore be justified.

A better calibration of planetary impactors in the outer solar system is provided by the Outer Solar System Origins Survey (OSSOS; Bannister et al. 2018) observations of Centaurs (Cabral et al. 2019; Dorsey et al. 2023). In Nesvorný et al. (2019), we used the N17 model to predict the orbital distribution and number of Centaurs. The model distribution was biased by the OSSOS simulator and compared with the OSSOS Centaur detections. We found a good fit to the observed orbital distribution, including the wide range of orbital inclinations, which was the hardest characteristic to fit in previous models (see Figure 2 and Marsset et al. 2019 for how the orbital inclinations of KBOs and Centaurs correlate with photometric color and Nesvorný et al. 2020 for modeling of the inclination–color relationship). The N17 dynamical model, in which the original population of outer disk planetesimals was calibrated from Jupiter Trojans (see above), implies that OSSOS should have detected 11 ± 4 Centaurs with semimajor axes a < 30 au, perihelion distances q > 7.5 au, and diameter D > 10 km (absolute magnitude H_r < 13.7 for a 6% geometric albedo). This is consistent with 15 actual OSSOS Centaur detections with H_r < 13.7.

Zoom In Zoom Out Reset image size

By slightly adjusting the N17 model to accurately match the OSSOS Centaur detections, Nesvorný et al. (2019) inferred that the inner scattered disk at 50 au < a < 200 au, from where most ECs evolve (N17), should contain (2.0 ± 0.8) × 10⁷ D > 10 km objects at the present time. This is consistent with independent estimates by Lawler et al. (2018) and Di Sisto & Rossignoli (2020). We can then rewind the history of the population using the N17 simulations to estimate that the original trans-Neptunian disk contained (8 ± 3) × 10⁹ planetesimals with D > 10 km. Reference population estimates for smaller-diameter cutoffs can be obtained from the size distribution of Jupiter Trojans: cumulative N(D) ∝ D^−2.1 down to D ∼ 3–5 km (Grav et al. 2011; Wong & Brown 2015; Yoshida & Terai 2017; Uehata et al. 2022; see discussion in Section 3.1) and perhaps even down to D ∼ 1–2 km.

With the orbital model and absolute calibration in place, we can estimate the impact rates of ECs and Centaurs on the outer planets and their moons. For that purpose, we repeat the simulations in N17. Specifically, we select the best case from N17, ¹² in which Neptune started at a_N,0 = 22 au; migrated at a rate proportional to $\exp (-t/{\tau }_{1})$, with τ₁ = 10 Myr, for the first 10 Myr; and then migrated at a rate $\propto \exp (-t/{\tau }_{2})$, with τ₂ = 30 Myr, for 10 Myr < t < 500 Myr (case 2 in N17). Time t is measured from the dispersal of the protoplanetary gas disk some 4.56 Gyr ago. The dynamical instability happened at t = 10 Myr in this case. The present properties of the EC population are not sensitive to the details of Neptune's early migration (e.g., N17). The galactic tide and stellar encounters were included in the simulation, while the putative Planet Nine (Trujillo & Sheppard 2014; Batygin & Brown 2016, N17) was not.

Since we are interested in the current impact rate of ECs/Centaurs (see Wong et al. 2019, 2021 and Bottke et al. 2023 for impact rates in the early solar system), we repeated the last segment of the N17 simulation, in which cometary reservoirs were evolved from 1 Gyr ago to the present time. By slicing this wide time interval into smaller segments and comparing the results, we verified that the EC/Centaur population changed little over the past 1 Gyr (e.g., the impact rate on Jupiter changed by <10% in the past 1 Gyr). We thus used the full statistics in the 1 Gyr long interval to infer the current impact rate. The original simulations started with 10⁶ test planetesimals in the disk between 24 and 30 au at t = 0. To further improve the model statistics, we cloned test bodies as they evolved toward the inner solar system. Specifically, in the last 1 Gyr interval, we monitored the heliocentric distance r of each body and cloned it 50 times when r first dropped below r* = 23 au. ¹³ This effectively corresponds to N = 5 × 10⁷ (initial) test planetesimals. The simulations were performed on 2000 Ivy Bridge cores of the NASA Ames Pleiades Supercomputer. We used the swift_rmvs4 integrator (Levison & Duncan 1994) and a 0.2 yr time step. All impacts of bodies on the outer planets were recorded by the N-body integrator.

To compute impact rates on the outer planets' moons, we first recorded all encounters of model bodies within ≃1 Hill radius of each planet. For every encounter, we then used Equation (3) from Nesvorný et al. (2004) to compute the collisional probability with moons. ¹⁴ The computation of collisional probabilities accounts for the moons' real orbits, including their (typically small) orbital eccentricities and inclinations. We verified that, for a moon on a circular orbit, the results are consistent with the Öpik equations (Öpik 1951; Zahnle et al. 1998; Z03; see footnote 8 in Nesvorný et al. 2004). The bodies that were bound to the planet were treated separately; their collisional probability was computed from Equation (13) in Nesvorný et al. (2003). In most cases, the dynamical and physical characteristics of moons (orbital elements, physical radii, surface gravities, etc.) were taken from the NASA JPL Planetary Satellites site. ¹⁵ We assumed spherical shapes for all moons and accounted for the gravitational focusing of comets by moons. The reduction in outbound impactor flux due to collisions with the planet ("shielding") was also accounted for (Lissauer et al. 1988). We accumulated the collisional probability over all recorded encounters and, following Z03, expressed it as a fractional impact probability relative to Jupiter (Tables 1–3).

We considered two ways an EC can become inactive: it becomes dormant, either because a refractory mantle forms or because the comet loses all volatiles, or it breaks into small pieces and disappears. (We ignore the possibility that the comet breaks into several large fragments.) We must distinguish between these possibilities because they have different implications for the impact flux. The impact flux is expected to be higher if the comets become dormant because the nucleus is still intact. If, instead, the ECs disrupt, they must be removed from the pool of impactors. Studies of meteor orbits suggest that disruption is the main physical loss mechanism for JFCs (Ye et al. 2016). Nonetheless, we considered two end-member models. In each, we assumed that comets have a physical lifetime of N_p(2.5) = 600 perihelion passages within 2.5 au, as in N17 (Section 1). ¹⁶ In the first case (Sections 3.1.1 and 3.2.1), we assumed that the comets became dormant and retained them in the impact flux calculation. In the second case (Sections 3.1.2 and 3.2.2), we assumed that the comets were disrupted and removed them from the impact flux calculation. We only model comet fading or disruption at small perihelion distances in this work (as parameterized by N_p (2.5)). Comet fading and/or disruption at larger perihelion distances (e.g., due to activity driven by supervolatiles or tidal encounters with planets) is not accounted for.

3.1. Planetary Impacts

3.1.1. Results without Comet Disruption

In total, our model recorded 217 direct impacts on Jupiter, 70 impacts on Saturn, and 62 impacts on Uranus. We do not consider impacts on Neptune and its moons in this work. Neptune and its moons are bombarded not only by Centaurs but also by SDOs. We do not have complete statistics for SDOs because we cloned the bodies in our model at r* = 23 au. Interestingly, 24%, 11%, and 6% of the bodies were bound to Jupiter, Saturn, and Uranus, respectively, at the time of impact. ¹⁷ For Jupiter, this is close to the 21% reported in Levison et al. (2000) and 15% ± 2% given by Kary & Dones (1996). For Jupiter, three out of 217 impactors (∼1.4%) had Tisserand parameters with respect to Jupiter T_J < 2, indicating that only 1%–2% of the impactors were near-isotropic comets (LD97) from the Oort cloud (Vokrouhlický et al. 2019). The great majority of impactors (98%–99%) were ECs with 2 < T_J < 3.05. The EC precursors—the low-inclination Centaurs evolving from the scattered disk—also dominate impacts on Saturn and Uranus (e.g., only one of 62 Uranus impactors had a retrograde heliocentric orbit).

The impact rate of ECs/Centaurs on Jupiter is computed as follows. Having effectively 5 × 10⁷ test bodies in the original trans-Neptunian disk (10⁶ original planetesimals × 50 clones for bodies that reached r* = 23 in the past billion years; see Section 2), we recorded 217 impacts in 1 Gyr. According to the calibration discussed above, there were (8 ± 3) × 10⁹ planetesimals with D > 10 km when Neptune began to migrate. The current rate of impacts of D > 10 km bodies on Jupiter (i.e., the average for the last billion years; see Section 2) is therefore 217 × 8 × 10⁹/(5 × 10⁷ × 10⁹) = 3.5 × 10⁻⁵ yr⁻¹, implying a timescale for impacts by D > 10 km bodies of ≃2.9 × 10⁴ yr. If we adopt the reference cumulative size distribution N(D) ∝ D^−2.1 from N17 (Section 2) and assume that it extends down to D = 1 km, we estimate a timescale for impacts by D > 1 km bodies of ≃230 yr. We use this timescale as a baseline in the rest of this paper. For comparison, Dones et al. (2009) inferred the impact flux on Jupiter from the historical record of close approaches and impacts (Schenk & Zahnle 2007). They estimated an impact rate of at least 4 × 10⁻³ yr⁻¹ for D > 1 km, implying a timescale of ≲250 yr, which is consistent with our model results. Saturn and Uranus receive 0.32 and 0.29 of the Jupiter impact flux, respectively, implying impact timescales of ≃720 and ≃790 yr, respectively, for D > 1 km impacts.

One major source of uncertainty of our model estimates is the extrapolation of the impact rate from D > 10 to 1 km. Above, we assumed a single power-law slope of 2.1 for simplicity. Subaru telescope observations of Jupiter Trojans indicate that their cumulative power index changes from ≃2.25 for D > 5 km to ≃1.8 for D < 5 km, assuming that albedo is independent of size (e.g., Uehata et al. 2022). If the shallower slope is extended all the way down to 1 km, we would infer a ≃340 yr timescale for impacts of D > 1 km bodies on Jupiter, i.e., an impact rate of ∼2/3 of our previous estimate. ¹⁸ It may also be that the size distribution of Jupiter Trojans is not a good proxy for cometary impactors. If we generously assume that the cometary size distribution is steeper, N(D) ∝ D^−2.4, for 1 km < D < 10 km, based on craters on Iapetus's dark terrain (Kirchoff & Schenk 2010), we obtain a ≃120 yr timescale for impacts of D > 1 km bodies on Jupiter.

3.1.2. Results with Comet Disruption

Figure 3 shows how the number of planetary impacts is reduced when we assume that comets disrupt after a certain number of orbits with q < 2.5 au. As expected, comet disruption has the largest effect for Jupiter, where the impact rate, assuming N_p(2.5) = 600, is reduced by 27% from its value neglecting disruption. For Saturn and Uranus, the corresponding factors are 6% and 2%. ¹⁹ The average time between impacts of D > 1 km bodies on Jupiter, Saturn, and Uranus is then 320, 760, and 810 yr, respectively (here and elsewhere in the paper, we scale from the nominal reference timescale of 230 yr for Jupiter discussed above). For all planets, the reduction factor decreases if we assume long physical lifetimes for comets. For example, for N_p(2.5) = 3000, as appropriate for D ∼ 10 km ECs from N17, the impact flux on Jupiter is reduced by only 8%.

Zoom In Zoom Out Reset image size

The fraction of bound impactors is larger when comet disruption is included in the model. For example, without accounting for comet disruption, we find that 24% of Jupiter impactors were bound to Jupiter when they impacted. With N_p(2.5) = 600, which should be appropriate for D ∼ 1 km comets (N17), we find that the fraction of bound impactors is ≃30%. Bound Jupiter impactors are slightly less likely to have many perihelion passages below 2.5 au than the population of ECs as a whole. Bound objects typically have T_J ≈ 3 (Kary & Dones 1996). Bodies that encounter Jupiter at higher velocities, and hence have smaller Tisserand parameters, are more easily scattered by Jupiter into orbits with small perihelion distances (Levison & Duncan 1997; see Figure 9 of Fernández et al. 2018). This means that events such as Shoemaker–Levy 9, which was tidally disrupted by Jupiter in 1992 and collided with the planet in 1994, would happen slightly more often—relative to impacts from unbound orbits—if comets disrupt. ²⁰

3.2. Impacts on Moons

3.2.1. Results without Comet Disruption

The impact probabilities on the outer planets' moons, both with and without comet disruption, are reported in Tables 1–3. There are some notable differences from Z03. We obtain higher impact speeds for the outermost moons. For example, we find 〈v_imp〉 = 7.1 km s⁻¹ for Phoebe, whereas Z03 reported 〈v_imp〉 = 3.2 km s⁻¹ (here v_imp is the impact speed of an individual body, and 〈v_imp〉 is the mean impact speed computed over all recorded impacts). Ultimately, Z03 based their encounter velocity (v_∞, the velocity at "infinity"; in practice, the planet's Hill sphere) distribution on objects that struck Jupiter in LD97. That distribution is biased toward low-velocity encounters because of gravitational focusing by Jupiter. In reality, there is a wide range of values of v_∞ at each planet (Figure 4). The outermost satellites of the giant planets have orbital velocities v_orb < 〈v_∞〉, where 〈v_∞〉 is the average encounter velocity of bodies at the planet's Hill sphere (v_orb is the orbital velocity of a moon around its parent planet). Therefore, gravitational focusing by the planet is a minor effect for them. This implies that impacts on the outer moons sample the high-velocity portion of the background velocity distribution more than the inner moons and the planets themselves (Zahnle et al. 1998, Section 2.1).

Zoom In Zoom Out Reset image size

By contrast, the inner satellites of the giant planets have orbital velocities v_orb > 〈v_∞〉. As a result, Centaurs/ECs are more strongly focused gravitationally. In the limit v_orb ≫ 〈v_∞〉, the impact rate per unit area on a moon ${{ \mathcal F }}_{\mathrm{imp}}\mathop{\propto }\limits_{\sim }1/({a}_{\mathrm{moon}}{v}_{\infty }^{2}$), where a_moon is the moon's semimajor axis (see, e.g., Equation (4) in Z03). For reference, 〈v_∞〉 = 4.1 km s⁻¹ for bodies that impact Jupiter, while 〈v_∞〉 = 8.8 km s⁻¹ for bodies that cross Jupiter's orbit. These values are, respectively, approximately 30% and 70% of Jupiter's mean orbital speed.

Other differences may arise from the dependence on the structure of the source regions assumed in different works. Z03 adopted the results from Levison & Duncan (1997) and Levison et al. (2000), where bodies started with low inclinations in the classical Kuiper Belt. This leads to a dynamically colder population of Centaur/EC impactors and lower impact speeds for more distant satellites. Here, instead, Centaurs/ECs evolve from the scattered disk in our model (N17) and so have higher inclinations and larger impact speeds.

The mean impact speeds for the inner moons that we find in this work are typically slightly smaller than the mean impact speeds given in Z03. For a monodisperse encounter velocity (v_∞), higher impact speeds for the outer moons would imply higher impact speeds for the inner moons as well. If we approximate the distribution of v_∞ as a Maxwellian, the impact probability weighted mean encounter velocity is $\tfrac{1}{2}{v}_{\infty }$ in the limit v_orb ≫ 〈v_∞〉. The slightly lower impact speeds for the inner moons probably reflect this sampling and the different velocity distributions we and Z03 assume (see Zahnle et al. 2001, Equation (1)). Our low-velocity tail probably extends to slightly lower speeds than that of Z03, and this tail is preferentially sampled by the inner moons.

We find higher impact rates for distant moons and lower impact rates for the innermost moons (Tables 1–3; recall that these rates are normalized to Jupiter) relative to Z03. Therefore, our impact probability profiles with orbital radius are flatter than in Z03. For example, the outer moons of Saturn have impact probabilities that are up to ∼twice as high as in Z03 (e.g., Phoebe has P_imp = 1.8 × 10⁻⁸ here and P_imp = 8.7 × 10⁻⁹ in Z03). The higher (normalized) impact probabilities for the outermost moons are a consequence of higher encounter velocities in our model, which result in lower gravitational focusing on the planets.

The mid-sized moons of Saturn have impact probabilities that are slightly smaller than reported in Z03 (e.g., Mimas has P_imp = 1.7 × 10⁻⁶ in Z03 and P_imp = 1.1 × 10⁻⁶ in this work; the impact probability P_imp is relative to Jupiter). Part of this difference results from the lower impact rates on Saturn relative to Jupiter that we find here (0.32 versus 0.42 in Z03).

For the innermost moons (Metis–Thebe for Jupiter, Pan–Janus for Saturn, and Portia and Puck for Uranus), the differences with Z03 are larger. These result from, e.g., differences in assumed sizes for the satellites and shielding by the planet (Lissauer et al. 1988), which we account for but Z03 did not. As an independent check, we verified that our impact probabilities for the innermost moons follow the expected scaling. For these moons, the value of P_imp should be approximately ${({R}_{\mathrm{moon}}/{R}_{\mathrm{planet}})}^{2}({R}_{\mathrm{planet}}/a)$, where R_moon and R_planet are the physical radii of the moon and planet, and a is the semimajor axis of the moon. The first term is the ratio of cross sections, while the second term accounts for gravitational focusing by the planet.

3.2.2. Results with Comet Disruption

The impact flux is reduced when the effects of comet disruption are included in the model (Tables 1–3). The flux is reduced by ∼25%–45% for the Jovian moons, ∼5%–15% for the Saturnian moons, and <1% for the Uranian moons. This trend, with the reduction factor decreasing with heliocentric distance, is expected because comet disruption mainly reduces the impactor population at smaller heliocentric distances. Compared to planetary impacts (Section 3.1), the reduction factors for moons are larger. For example, the impact flux on Jupiter was reduced by 27% from the original flux for N_p(2.5) = 600 (appropriate for ∼1 km impactors). For Jovian moons, however, the impact flux is reduced by ∼25%–45% for N_p(2.5) = 600, and there is a trend with smaller reduction for inner moons and larger reduction for outer moons.

The impact speeds on the outer moons of Jupiter are slightly lower in the case with disruption. For example, the mean impact speeds on Himalia are 8 and 9 km s⁻¹ with and without disruption, respectively (Table 1). This most likely happens because in the case with disruption, there is not enough time to excite the heliocentric orbits of EC impactors as much. The effects of disruption on the impact speeds for the moons of Saturn and Uranus are small (Tables 2 and 3).

Table 2. The Impact Probabilities P_imp and Average Impact Speeds 〈v_imp〉 for Saturnian Moons

	Zahnle et al. (2003)		This Work
			No Disruption		With Disruption
	Pimp	〈vimp〉	Pimp	〈vimp〉	Pimp	〈vimp〉	τimp
		(km s−1)		(km s−1)		(km s−1)	(Myr)
Saturn	0.42	⋯	0.32	⋯	0.30	⋯	⋯
Pan	⋯	⋯	7.9 × 10−9	30.1	7.4 × 10−9	30.6	⋯
Atlas	⋯	⋯	8.0 × 10−9	28.5	7.0 × 10−9	29.3	⋯
Prometheus	1.7 × 10−7	32	6.5 × 10−8	29.6	5.9 × 10−8	30.4	3800
Pandora	1.0 × 10−7	31	5.6 × 10−8	29.6	5.1 × 10−8	30.3	4400
Epimetheus	1.8 × 10−7	30	1.3 × 10−7	29.1	1.2 × 10−7	30.0	1900
Janus	4.5 × 10−7	30	2.7 × 10−7	28.6	2.5 × 10−7	29.4	920
Mimas	1.7 × 10−6	27	1.1 × 10−6	26.2	9.7 × 10−7	26.6	230
Enceladus	2.2 × 10−6	24	1.5 × 10−6	23.1	1.4 × 10−6	23.3	160
Tethys	7.9 × 10−6	21	6.0 × 10−6	21.0	5.4 × 10−6	20.9	41
Telesto	3.5 × 10−9	21	3.3 × 10−9	21.0	3.0 × 10−9	20.9	⋯
Calypso	3.5 × 10−9	21	2.4 × 10−9	21.0	2.2 × 10−9	20.9	⋯
Dione	7.1 × 10−6	19	5.6 × 10−6	18.7	5.1 × 10−6	18.4	45
Helene	5.8 × 10−9	19	6.2 × 10−6	18.7	5.6 × 10−6	18.6	⋯
Rhea	9.6 × 10−6	16	7.8 × 10−6	15.8	7.0 × 10−6	15.8	32
Titan	5.4 × 10−5	10.5	4.8 × 10−5	11.2	4.5 × 10−5	11.2	5.0
Hyperion	1.0 × 10−7	9.4	1.0 × 10−7	10.4	9.2 × 10−8	10.3	2500
Iapetus	1.4 × 10−6	6.1	1.6 × 10−6	7.9	1.5 × 10−6	7.7	150
Phoebe	8.7 × 10−9	3.2	1.8 × 10−8	7.1	1.6 × 10−8	7.0	⋯

Note. The impact probabilities are given with respect to the impact probability on Jupiter for the case with no comet disruption. To infer the rate of impacts on a moon, P_imp needs to be multiplied by the impact rate of bodies on Jupiter (see Section 3.1.1 and the caption of Table 1). For example, the rate of impact of 1 km bodies on Enceladus is 6.7 × 10⁻⁹ yr⁻¹, or one such impact every ∼150 Myr, for our nominal size distribution and no cometary disruption. For the same size distribution but with cometary disruption, the corresponding rate is 6.4 × 10⁻⁹ yr⁻¹ and a timescale of ∼160 Myr. The last column reports the average time between impacts of D > 1 km bodies (τ_imp) for the case with comet disruption and our nominal size distribution (timescales exceeding the age of the solar system are not shown).

Download table as:  ASCIITypeset image

Table 3. The Impact Probabilities P_imp and Average Impact Speeds 〈v_imp〉 for Uranian Moons

	Zahnle et al. (2003)		This Work
			No Disruption		With Disruption
	Pimp	〈vimp〉	Pimp	〈vimp〉	Pimp	〈vimp〉	τimp
		(km s−1)		(km s−1)		(km s−1)	(Myr)
Uranus	0.25	⋯	0.29	⋯	0.28	⋯	⋯
Portia	4.7 × 10−7	18	4.5 × 10−7	16.9	4.5 × 10−7	16.9	500
Puck	6.6 × 10−7	15	4.9 × 10−7	14.1	4.9 × 10−7	14.1	460
Miranda	5.2 × 10−6	12.5	3.5 × 10−6	12.7	3.5 × 10−6	12.7	65
Ariel	2.1 × 10−5	10.3	1.4 × 10−5	10.6	1.4 × 10−5	10.6	16
Umbriel	1.6 × 10−5	8.7	1.3 × 10−5	9.2	1.3 × 10−5	9.2	18
Titania	1.8 × 10−5	6.8	1.3 × 10−5	7.3	1.3 × 10−5	7.3	17
Oberon	1.3 × 10−5	5.9	1.0 × 10−5	6.5	1.0 × 10−5	6.5	23
Sycorax	⋯	⋯	2.6 × 10−8	3.9	2.5 × 10−8	3.9	⋯

Note. The impact probabilities are given with respect to the impact probability on Jupiter for the case with no comet disruption. To infer the rate of impacts on a moon, P_imp needs to be multiplied by the impact rate of bodies on Jupiter (see Section 3.1.1 and the caption of Table 1). For example, the rate of impact of 1 km bodies on Miranda is 1.5 × 10⁻⁸ yr⁻¹, or one such impact every ∼65 Myr, for our nominal size distribution (the results with and without cometary disruption are the same). The last column reports the average time between impacts of D > 1 km bodies (τ_imp) for the case with comet disruption and our nominal size distribution (timescales exceeding the age of the solar system are not shown).

Download table as:  ASCIITypeset image

The lower cratering rates in the case that includes cometary disruption imply somewhat older surface ages for lightly cratered terrains on outer planet satellites such as Enceladus and Europa. For example, Z03 estimated a surface age of Europa of between 30 and 70 Myr (more recent papers infer slightly older ages: 60–100 Myr in Zahnle et al. 2008 and 40–90 Myr in Bierhaus et al. 2009). Scaling from Z03, our lower impact flux in the model with comet disruption would imply a surface age of between 45 and 105 Myr. The catastrophic disruption timescales of small moons reported in Z03, which assume that a crater with a diameter equal to the diameter of the satellite dooms the moon, would also be longer.

Figure 5 shows how the impact flux depends on the orbital radius of a moon. Here we disregard the moons' sizes by normalizing the impact flux of D > 1 km comets per 10⁶ km² of the moon's surface (per Myr). When comet disruption is not accounted for, the normalized impact flux is ${{ \mathcal F }}_{\mathrm{imp}}\simeq$ 0.0005–0.002/(10⁶ km² Myr) for the outer moons and ${{ \mathcal F }}_{\mathrm{imp}}\simeq 0.01$–0.05/(10⁶ km² Myr) for the inner moons. The higher flux on the inner moons is due to gravitational focusing by the planets. This trend persists when we account for comet disruption (bottom panel of Figure 5). Given that the flux reduction is larger for the outer moons, however, the overall slope of the ${{ \mathcal F }}_{\mathrm{imp}}$ dependence on a_moon becomes slightly steeper with comet disruption. For example, the outer/irregular moons of Jupiter have ${{ \mathcal F }}_{\mathrm{imp}}\,\simeq$ 0.0004/(10⁶ km² Myr) with N_p(2.5) = 600 compared to ${{ \mathcal F }}_{\mathrm{imp}}\,\simeq$ 0.0006/(10⁶ km² Myr) without comet disruption.

Zoom In Zoom Out Reset image size

The impact flux on the satellites depends on whether comets undergo disruption or become dormant. The flux is higher if comets become dormant (Tables 1–3). We expect that both cases apply to some degree. For example, if all comets become dormant, the ratio of the number of dormant to active ECs should be roughly 30–60. This is because the dynamical lifetime of ECs with q < 2.5 au is t_dyn ∼ 3 × 10⁵ yr (LD97), while their physical lifetime is t_phy ∼ 5000–12,000 yr (LD97, N17). The ratio of the number of dormant to active ECs should be ∼t_dyn/t_phys.

Licandro et al. (2016) obtained WISE observations of asteroids on (mostly Jupiter family–like) cometary orbits (ACOs) and inferred their size distribution. They found that the number of ACOs is smaller than the number of JFCs for diameters of <10 km. This suggests that most ECs end their lives by being disrupted instead of becoming dormant (see also Ye et al. 2016). Therefore, the impact fluxes that we calculate for the case with cometary disruption are likely to be more realistic. This may not be true for Centaurs, which are less likely to be disrupted because of their greater distances from the Sun.

The main results of this work can be summarized as follows.

• 1.  
  We determined the current rates at which comets and Centaurs strike the moons of Jupiter, Saturn, and Uranus. Compared to Z03, we find a higher impact flux on the outer moons and a smaller impact flux on the inner moons. This is a consequence of the larger orbital inclinations in our model, in which the scattered disk is the main source of comets/Centaurs. The impact speeds on the outer moons are significantly higher (e.g., 7.1 km s⁻¹ for Phoebe compared to 3.2 km s⁻¹ reported for Phoebe in Z03), implying larger craters for a given impactor size.
• 2.  
  When comet disruption is accounted for, the impact flux on the outer planets and their moons is reduced. The reduction factor depends on the adopted disruption model. For example, for N_p(2.5) = 600, which was the preferred value for D ∼ 1 km comets in N17, the impact flux on Jupiter is reduced to 73% of the original value (the mean interval between impacts of D > 1 km comets on Jupiter is ∼230 yr without comet disruption and ∼315 yr with comet disruption).
• 3.  
  For N_p(2.5) = 600, the impact flux is reduced by ∼25%–45% for Jovian moons, ∼5%–15% for Saturnian moons, and <1% for Uranian moons. The lower cratering rates in the case that includes cometary disruption implies somewhat older surface ages for lightly cratered terrains on outer planet satellites such as Enceladus and Europa. For example, the reduced impact fluxes with comet disruption would imply older surfaces (e.g., an ∼45–105 Myr surface age for Europa).
• 4.  
  The impact rate on Jupiter changes by less than 10% in the past 1 Gyr. Wong et al. (2019) found a much stronger decay during the first billion years of the solar system. Together, these results show that most of the cratering occurred early on (Wong et al. 2021; Bottke et al. 2023).
• 5.  
  The average time between impacts of D > 1 km bodies for the case with comet disruption and our nominal size distribution is 2.7, 5.1, 3.2, and 5.7 Myr for Io, Europa, Ganymede, and Callisto, respectively (Table 1). It is 42, 45, 32, and 5.0 Myr for Tethys, Dione, Rhea, and Titan (Table 2) and 16, 18, 17, and 23 Myr for Ariel, Umbriel, Titania, and Oberon (Table 3).

The simulations were performed on the NASA Pleiades Supercomputer. We thank the NASA NAS computing division for continued support. The work of D.N. was supported by the NASA solar system Workings program. The work of D.N., L.D., and M.D.P. was supported by the NASA Cassini Data Analysis Program. This material is based in part on work done by M.W. while serving at the National Science Foundation. We thank the reviewers, Darryl Seligman and Wes Fraser, for helpful suggestions on the submitted manuscript.