© The Authors 2025

1 Introduction

The Jovian satellite Europa and the Saturnian satellite Enceladus are attractive astrobiological targets due to their activity, subsurface oceans, and mostly young surface ages derived from a lack of impact craters. Yet any meaningful surface age for these bodies has to rely on an accurate derivation of the current impact rate on these bodies.

Zahnle et al. (1998) and Zahnle et al. (2003) used the numerical simulations of the outer Solar System of Levison & Duncan (1997), and a Monte Carlo method to compute the current cratering rates of heliocentric planetesimals on the regular satellites of the giant planets. Zahnle et al. (2003) presented two cases, one based on the abundance of small objects near Jupiter (Case A), and another from the impact crater size-frequency distribution on Triton (Case B). These works assumed that the total number of planetesimals crossing the orbits of the giant planets declined with time as t⁻¹ (Holman & Wisdom 1993). With this assumption the heavily cratered regions of most of the regular icy satellites were computed to be older than 4 Ga, and often as old as the Solar System itself when Case A was considered, but much younger under the scenario Case B. However, the decline in impact flux does not evolve as t⁻¹ but as a stretched exponential (Wong et al. 2021), which can substantially modify the ages, particularly for old surfaces. As such, from these two studies alone there is a large variety in the impact flux onto these satellites and in the resulting surface ages.

More recently Nesvorný et al. (2023) attempts to calculate the impact rate on the regular satellites of the giant planets. They rely on the extensive numerical simulations of the outer Solar System in Nesvorný et al. (2017). They find similar impact probabilities and impact velocities of heliocentric bodies with the satellites compared with Zahnle et al. (2003). Both Nesvorný et al. (2023) and Zahnle et al. (2003) scale the impact rates and impact probabilities of bodies with the satellites to Jupiter.

This study is our attempt to compute the current impact rate based on our own simulations. It builds on the previous studies by Wong et al. (2019, 2021, 2023) and ongoing attempts to compute the surface ages of Enceladus. Since our previous works focused on the satellites of Jupiter, Saturn and Uranus, and because the history of the Neptune system is uncertain, we do not compute the current impact rate onto Neptune’s regular satellites system. The location of Neptune as the outermost planet also makes the impact rate onto its satellites much more sensitive to the underlying trans-Neptunian small body population. Unlike most previous studies we follow Levison et al. (2000) and base our current impact rate onto the giant planets and their satellites from the leakage rate of scattered disc objects beyond Neptune into the Centaur region. Due to the diverse outcomes that exist in the literature of estimates of the current impact rates, we explicitly describe the previous work that we build upon here, our methods and choices, followed by a detailed comparison with other works.

2 Methodology and existing data

There are several different methodologies in the established literature to compute the current impact rate on the icy satellites. This work describes our approach to the problem. We also give a few comparisons to previous works.

2.1 Database of numerical N-body simulations of the outer Solar System

The impact rate of the satellites of the giant planets relies on the dynamical evolution of the outer Solar System. For this work we rely on a database of simulations that are in Wong et al. (2019, 2021, 2023). We summarise these simulations here.

We ran two sets of simulations of giant planet migration. Both sets of migration simulations began with the giant planets on a more compact configuration surrounded by a disc of planetesimals. The first set employed the Graphical Processing Unit (GPU) integrator Gravitational ENcounters with GPU Acceleration (GENGA) (Grimm & Stadel 2014; Grimm et al. 2022) wherein the giant planets feel the gravitational force from the planetesimals, and the migration of the giant planets was the result of energy and angular momentum conservation as they scattered the planetesimals (Fernandez & Ip 1984). These simulations are denoted as GPU M. The second set used the Regularised Mixed Variable Symplectic (RMVS) integrator (Levison & Duncan 1994), version 4, wherein the migration of the planets was orderly due to the application of fictitious forces (Levison et al. 2008), and the planetesimals were modelled as massless test particles; this set is hereafter denoted as RMVS M. From the GENGA simulations we obtained the impact probability of a planetesimal with a giant planet and we recorded the vectors of planetesimals one time step after their closest approach with a giant planet within 40 planetary radii (R_P). From these vectors we enacted a high number of simulations wherein we subjected the regular satellites to the planetesimals that flew past them. We computed the impact probability of a planetesimal with a satellite from these flyby simulations and the encounter probability with the planets Wong et al. (2019).

We discovered that the RMVS migration simulations did not have enough resolution, or remaining particles to accurately compute the impact chronology onto the giant planets over the last 3.5 billion years (Gyr). We ran additional simulations with RMVS and GENGA: we cloned all planetesimals after 1 Gyr of dynamical evolution in the RMVS M simulations. The semi-major axis, perihelion distance, and inclination of the planetesimals after 1 Gyr are shown in Figure 1. We ran these cloned planetesimals for another 3.5 Gyr to obtain a detailed impact chronology with high resolution data and robust impact statistics; these simulations are denoted as RMVS C and GPU C respectively. A summary of the different simulations is presented in Table 1. One additional GENGA simulation, not listed in Table 1, integrated all test particles that remained after 1 Gyr of simulation time for another 3.5 Gyr without cloning (GPU Mixed). From the simulations we computed the impact probabilities (P_P) with the giant planets, and from the GENGA simulations we obtained the encounter probabilities within 40 planetary radii (P_enc) as well.

In what follows we consider the migration phase – the first 1 Gyr – to be separate from the last 3.5 Gyr and we do not discuss the migration phase because it is unimportant for our purpose: Nesvorný et al. (2023) state that the present properties of the ecliptic comet population are not sensitive to the details of Neptune’s early migration. The reason for treating the last 3.5 Gyr differently is that in order to compute the current impact rate onto the satellites we need to calculate the current injection rate of comets from beyond Neptune towards the giant planets. For this we need to simulate a sample of particles whose orbital distribution is similar to the unbiased observed population, or, in our case, the population of particles remaining after 4.5 Gyr of dynamical evolution. This sample has a very different orbital distribution from the initial conditions of the migration simulations, and as such the impact and encounter probabilities with the giant planets of this population are different from during the migration phase.

Fig. 1 Orbital distribution of semi-major axis and perihelion distance of planetesimals after 1 Gyr of evolution of giant planet migration. The colour coding indicates inclination. Data from Wong et al. (2023).

2.2 Computing the satellite impact probabilities

There are two ways to calculate the impact probabilities of the planetesimals with the satellites (P_s). The first method makes use of the relative impact probabilities of the planetesimals on the satellites versus on the host planet and relies on an analytical approximation (Zahnle et al. 2003). We do not use that here.

The second method relies on numerical flyby simulations as detailed in Wong et al. (2019). We can express P_s as the fraction of planetesimals that collided on the satellites and multiply it by the probability of the planetesimals coming within 40 R_P of a giant planet; we call the latter the encounter probability (P_enc).

We have $$P_{\text{s}}=\frac{N_{\text{i}}}{N_{\text{tot}}}{P}_{\text{enc}},$$(1)

where N_i is the recorded number of impacts on the satellites, and N_tot is the total number of planetesimals that come with a distance of 40 R_P. From the flyby simulations of the three planet-satellites systems, we could count the number of impacts on each regular satellite and calculate the absolute impact probabilities from which we can compute the impact rates.

2.3 Re-enacting the close encounters

For this work we decide to redo the flyby simulations using the CPU version of GENGA. The reason is the following.

Following Chambers (1999), the GENGA integrator computes a critical radius around each massive body in the simulation, whose magnitude is given by R_c =max(n₁R_H, n₂υ_Kh). Here R_H is the Hill radius of the body, h is the time step, and n₁ and n₂ are constants set by the user; their default values are 3.0 and 0.4. For example, for Mercury using a time step of h = 0.01 yr we have R_c = 0.04 astronomical units (au) which is dominated by the second term, but for Jupiter with h = 0.4 yr we have R_c = 1.07 au which is dominated by the first term. Within this critical distance two bodies are considered to be within a close encounter, and the integration routine smoothly changes from the Mixed Variable Symplectic (MVS) map (Wisdom & Holman 1991) beyond R_c to the Bulirsch-Stoer routine (Grimm & Stadel 2014; Grimm et al. 2022) when the two bodies are closer than 0.2R_c. The changeover between the MVS and Bulirsch-Stoer regimes needs to be done smoothly so that during a close encounter the integrator should compute at least a few steps inside the changeover regime between 0.2 R_c and 1 R_c.

The default values of n₁ = 3.0 and n₂ = 0.4 in GENGA suffice for many problems: the value of n₁ is identical to the default value in RMVS (Levison & Duncan 1994). It is typically sufficient for low-velocity encounters, or when the mass ratio between the central body and the planets or satellites is greater than O(10⁻⁵). However, when the masses are lower, the key distance factor is not the number of Hill radii, but rather the number of time steps used to sample the changeover region, in other words, the second term. A pair of objects with lunar mass can easily travel 3 mutual Hill radii in a single time-step. In this case the switch from the symplectic regime to the Bulirsch-Stoer regime is instantaneous, leading to large errors during the close encounter. As such, for high-velocity encounters, the n2 parameter becomes dominant (Chambers 1999).

The average encounter speed between the regular satellites of the giant planets and the hyperbolic heliocentric planetesimals can be computed with Pythagoras’ theorem, and is approximately $\sqrt{3}{v}_K$ (Zahnle et al. 2003), so that within a time step the relative distance the satellite and a planetesimal travel is about $\sqrt{3}{v}_{\text{K}}h$. This distance is typically much greater than a few Hill radii of the satellites. As such, the n₂ term dominates. It is likely that n₂ = 0.4 is too low and the simulations may ‘step over’ close encounters as well as collisions between the satellites and the planetesimals. As such, we redid the simulations by setting n₂ = 1. We reckon that this value is high enough because GENGA begins to check for potential close encounters at a distance of $\sqrt{3}{R}_c$ (in other words, at $\sqrt{3}{v}_{\text{K}}h$).

Due to the low impact probabilities with the satellites it is necessary to re-enact many more encounters than are obtained during the giant planet migration simulations. We sampled the velocity distribution of the planetesimals as they penetrated the 40 planetary radii (R_P) sphere, and recorded their values between 38 and 42 planetary radii. We created 25 000 massless planetesimals in this manner per simulation and ran 20 000 simulations per planet for a total of 500 million planetesimals per planet. Our maximum number of planetesimals per simulation was limited by computing time: higher numbers of planetesimals would increase the computing time more sharply than a linear increase due to greater RAM and CPU cache occupation. From testing we found that using 25 000 planetesimals using one CPU core gave excellent results, with each simulation lasting 12–17 seconds. In these flyby simulations all planetesimals passed by the satellites at once.

The isotropic initial position vectors of the planetesimals are computed according to $$\begin{array}{l}\cos \theta =2{\xi }_1-1,\\ \sin \theta =\sqrt{1-{\cos }^2\theta }\\ \varphi =2\pi {\xi }_2,\\ x=S\sin \theta \cos \varphi ,\\ y=S\sin \theta \sin \varphi ,\\ ɀ=S\cos \theta ,\end{array}$$(2)

where ξ₁ and ξ₂ are uniform random numbers in the interval [0, 1], and S is the distance to the planet corresponding to 40 R_P. The angles θ and ϕ are the spherical latitude and longitude, respectively. The assumption of isotopic flybys has been validated by Zahnle et al. (1998) and the output from our GPU simulations confirm this. The velocities of the planetesimals are computed from angular momentum conservation and the requirement that they are initially all moving towards the planet. Their components are computed as $$\begin{array}{l}\psi =2\pi {\xi }_3,\\ v_3=-u\sqrt{{\xi }_4},\\ v_1=\sqrt{u^2-v_3^2}\cos \psi ,\\ v_2=\sqrt{u^2-v_3^2}\sin \psi ,\\ v_x=-v_1\sin \varphi +v_2\cos \theta \cos \varphi +v_3\sin \theta \cos \varphi ,\\ v_y=v_1\cos \varphi +v_2\cos \theta \sin \varphi +v_3\sin \theta \sin \varphi ,\\ v_{ɀ}=-v_2\sin \theta +v_3\cos \theta ,\end{array}$$(3)

where ξ₃ and ξ₄ are two additional uniform random numbers, u is the speed of the particle at 40 R_P, whose magnitude and distribution was computed from the GENGA simulations, and ψ modulates the magnitude of the velocity vector components orthogonal to the direction towards the planet. With this prescription, the impact parameter of the planetesimals, b = L/u, where L is the specific angular momentum, is uniformly distributed in b². Even though this impact parameter distribution is only strictly valid beyond 52 R_U for Uranus, 75 R_S for Saturn and beyond 110 R_J for Jupiter, this prescription does reproduce the trend from the GENGA planet migration simulations that P_enc/P_P ~ 50 (Wong et al. 2019). We computed the hyperbolic velocities at 40 R_P from the migration and post-migration GENGA simulations in Wong et al. (2023), and from there calculated the encounter velocities ‘at infinity’ as $v_{\infty }=\sqrt{u^2-v_{\text{esc}}{\left(r=40R_{\text{P}}\right)}^2}.$. For Jupiter 〈u〉 = 10.1 km s⁻¹ leading to 〈υ_∞〉 = 3.17 km s⁻¹, for Saturn 〈u〉 = 6.85 km s⁻¹ leading to 〈υ_∞〉 = 3.65 km s⁻¹, and for Uranus 〈u〉 = 4.25 km s⁻¹ leading to 〈υ_∞〉 = 2.39 km s⁻¹.

For the purpose of our calculations, we assumed that the satellites reside on their present-day orbits. The initial positions and velocities of the Jovian, Saturnian, and Uranian satellites are from the JPL Horizons website. The time step was 0.01 days or 864 seconds and the maximum simulation time was just 20 days, by which time most of the planetesimals had left. With this time step R_c = 12000 km assuming υ_K = 14 km s⁻¹ for Mimas, and GENGA will scan for encounters at a distance of at least 21 000 km, comparable to the relative distance travelled in a single time step between Mimas and a planetesimal, and will reduce the time step within this distance to the satellites. We therefore expect all collisions to be detected. For those few planetesimals that impacted a satellite their impact velocities were recorded.

3 Results and outcomes

3.1 Fundamental simulation parameters

Many small bodies in the trans-Neptunian region of the Solar System are on meta-stable orbits. These orbits could appear to be stable for billions of years because they are temporarily trapped in mean-motion resonances with Neptune (Duncan et al. 1995; Duncan & Levison 1997). Non-resonant objects with long semi-major axes whose perihelion lies within ~37 au evolve on billion-year timescales, undergoing a random walk in semi-major axis (Fernández et al. 2004), while resonant objects, once dislodged, can become Neptune crossing (Duncan et al. 1995; Duncan & Levison 1997). At short semi-major axis (often ≲100 au) secular effects are important, which alter the eccentricity, and thus the perihelion distance. These secular effects can flip a body from non-crossing to a crossing orbit, changing them from a scattered disc object to a Centaur or ecliptic comet. These combined dynamical effects cause a slow leakage of scattered disc and resonant objects into the Centaur region.

In order to compute the current impact rate onto the regular satellites we need to rely on the distribution of ecliptic comets, and more specifically on the injection rate of scattered disc objects into the realm of the giant planets. Ecliptic comets are scattered disc objects that have left their source region and become Neptune crossers; their motion is thus dominated by encounters with the planets, and by secular oscillations in inclination and eccentricity. A subset of these comets are the short-period Jupiter-family comets, which are concentrated along the ecliptic plane, and whose Tisserand parameter with respect to Jupiter T > 2 (Levison & Duncan 1997). The Tis-serand parameter with respect to a planet is akin to the Jacobi integral in the circular restricted three-body problem, and is roughly conserved during encounters with said planet. We qualified the particles in our simulations that met all the following criteria as an ecliptic comet (EC). Following Wong et al. (2023) our definition focuses on the EC injection rate, and not so much their actual fate. Our criteria are:

• the planetesimal survived the first 1 Gyr during the migration episode;

• its perihelion reached q < 30 au in the last 3.5 Gyr of the simulation, and

• it is lost eventually due to collision with the planets or the Sun, ejection from the Solar System, or survived till the end of simulation at 4.5 Gyr.

As such, only those scattered disc objects that began to cross Neptune are counted as ecliptic comets.

Following Duncan et al. (1995) we computed the injection rate of ecliptic comets from the scattered disc as $$F_{\text{EC}}=N_{\text{SD}}\left|r_{\text{SD}}\right|f_{\text{EC}},$$(4)

where N_SD is the current number of objects in the scattered disc, r_SD is the rate at which the scattered disc population declines, and f_ec is the fraction of scattered disc objects that become ecliptic comets. The rate of decline of scattered disc objects is computed as the number of planetesimals that left the system (ΔN) divided by the final number of the particles at the end (N_fin), and the time interval (Δt) (Duncan et al. 1995). In other words $$r_{\text{SD}}=\frac{1}{N_{\text{fin}}}\frac{\text{Δ}N}{\text{Δ}t}.$$(5)

It is worth mentioning that r_SD depends on the time range over which it is calculated because the rate of decline is not constant with time; instead, the rate of decline is a function of time itself: it decreases with increasing time. To approximate the current injection rate we calculated ΔN and N_fin over the last 0.5 Gyr in the simulations. Doing the same calculation over the last 1 Gyr increases the rate of decline by about 15%, and this increases to 50% if we consider the decline over the last 1.5 Gyr. These differences indicate the non-linearity of the rate of decline. We show this in Fig. 2, where we plotted the remaining percentage of the initial population over the last 3.5 Gyr from one of the GPU simulations. It is clear that the decline is much more rapid between 1 and 2 Gyr than between 3.5 and 4.5 Gyr. The fraction of bodies remaining after 4.5 Gyr is about 20% of the population at 1 Gyr, that is, the population declines by a further 80% between 1 Gyr and 4.5 Gyr of evolution. The percentage remaining after 1 Gyr of migration is typically 0.68% and therefore the total population remaining after 4.5 Gyr is 0.15%. Nesvorný et al. (2023) reports a remaining fraction of 0.25%.

Using Equation (5) to compute r_SD is not without problems, however. The quantity ΔN changes discretely, and Poisson statistics dictate that the standard deviation ${\sigma }_{\text{Δ}N}=\sqrt{\text{Δ}N}$. If ΔN is a small quantity but N_fin is large, then the uncertainty in the ratio between different simulations, that is σ_rSD, will likely be low. But if N_fin is also a small quantity, with standard deviation ${\sigma }_{\text{Nfin}}=\sqrt{N_{\text{fin}}}$, then both the ratio itself, and the uncertainty in the ratio ΔN/N_fin, can take on an (arbitrarily) large range in values due to stochastic variations in ΔN between different simulations. The expectation is that ΔN ≪ N_fin so that if ΔN is a low number σ_rSD ~ |r_SD|(ΔN)^−1/2 ~ |r_SD| . As such, the expectation is that the uncertainties in r_SD will become large when the number of remaining particles in the simulations and their rate of decline is low. These uncertainties cannot be reduced by running more simulations with the same initial number of particles.

To estimate the total number of ecliptic comets we need their mean dynamical lifetime, τ_EC. Levison et al. (2000) argue that $${\tau }_{\text{EC}}\approx {\displaystyle {\int }_0^t\frac{N_{(s)}}{N_{\text{init}}}} \text{d}s\approx \frac{1}{N_{\text{init}}}{\displaystyle \sum _{i=1}^{N_{\text{init}}}{t}_i},$$(6)

where N_init is the initial number of test particles, and N_(s) is the number of test particles remaining at the time s < t. The last step is approximately valid because of the slow decline of planetesimals. The orbital elements or vectors of each body during the simulation are written to disk at regular time intervals: every 0.1 Myr for RMVS and every 0.3 Myr for GENGA. In this manner each simulation produced thousands of data frames. To compute τ_EC we, (a) searched each data frame for a particular particle, (b) determined whether its perihelion q < 30 AU, and then (c) counted the number of data frames the particle was in with q < 30 au. If its perihelion increased again beyond 30 au it was not counted. We computed τ_EC for each particle by multiplying the total number of data frames by the data output interval (either 0.1 or 0.3 Myr). We then averaged these individual times over all simulations per set to compute τ_EC.

The nominal average values of all of these quantities are listed in Table 2 for the different sets of numerical simulations that we performed. All the quantities derived from the simulations, as well as their uncertainties, are in the Supplementary Materials on Zenodo. In what follows we listed all the relevant quantities from all simulation sets, but for the analysis of the impact rate we used only the clone (C) simulations. The reason for this is that in the RMVS M simulations 〈∆N〉( = 1.34 so that σ_∆N = 1.12 and the uncertainty in the rate of decline is σ_rSD ~ |r_SD|, which is too great to derive a meaningful impact rate onto the giant planets.

We computed $r_{\text{SD}}=-{0.179}_{-0.040}^{+0.030}$ Gyr⁻¹ (2σ) from the RMVS C simulations, and r_SD = −0.149 ± 0.004 Gyr⁻¹ (1σ) from the combined GPU simulations. Most of our values fall in between those of previous literature estimates and their uncertainties from Brasser & Morbidelli (2013) (−0.163 ± 0.066), Fernández et al. (2004) (−0.15), Levison et al. (2006) (−0.27), Volk & Malhotra (2008) (−0.15 ± 0.05), but are much lower than that of Di Sisto & Brunini (2007) (−0.52). Some of these studies (Fernández et al. 2004; Di Sisto & Brunini 2007; Volk & Malhotra 2008) compute the injection rate from integrating the observed sample of then-known trans-Neptunian objects, while others (Levison et al. 2006; Brasser & Morbidelli 2013) rely on a distribution of planetesimals that results from long-term integrations of particles subjected to giant planet migration. The agreement between all these results is encouraging and suggests that the leakage rate may only be weakly dependent on the underlying orbital distribution of scattered disc objects.

From the RMVS C simulations we obtained $f_{\text{EC}}={68.1}_{-5.2}^{+20.0}\%$ (2σ), and f_EC = 78.7% ± 1.2% (1σ) from the GPU runs. The mean lifetime of the ecliptic comets in the RMVS clone simulations is ${\tau }_{\text{EC}}={167}_{-164}^{+423}$ Myr (2σ) and in the GENGA simulations the result is ${\tau }_{\text{EC}}={154}_{-151}^{+379}$ Myr (2σ). These large ranges are due to some simulation particles spending a long time on orbits only crossing one or both ice giants, for which the evolution timescale is of the order of 100 Myr (Duncan et al. 1987). We show examples of these long-lived particles in Fig. 3. We plotted the semi-major axis and perihelion distance of two particles that both survive the migration stage. The case displayed in orange spends a shorter time crossing Neptune as an EC than the black case. All of the average EC lifetimes listed above are comparable with that found by Levison et al. (2000), who report τ_EC = 190 Myr, but longer than found by Di Sisto & Brunini (2007) (72 Myr) and Tiscareno & Malhotra (2003) (9 Myr), but the last study only simulated ecliptic comets whose perihelia were close to Jupiter and Saturn, drastically reducing the dynamical lifetime.

For completeness we report that the total fraction remaining for the RMVS C simulations is ${0.19}_{-0.11}^{+0.17}\%$ (2σ) and 0.15% ± 0.05% (1σ) for the GPU C simulations. Last, we report the fraction of the total population that are ECs with q < 30 au, which is computed as N_EC = N_SD |r_SD|f_ECτ_EC. The mean and uncertainties were computed using a Monte Carlo sampling procedure from the underlying distributions in f_EC, r_SD, and τ_EC across all simulations in each set. For the RMVS C simulations we obtained $N_{\text{EC}}={0.024}_{-0.023}^{+0.061}$ N_sd (2σ), and $N_{\text{EC}}={0.018}_{-0.017}^{+0.027}$ N_sd (1σ) from the GPU simulations. The large uncertainties are almost entirely due to the distribution in τ_EC. Before we can calculate the current impact rate onto the giant planets and their satellites we need to know the number of scattered disc objects greater than some threshold diameter, N_SD.

Fig. 2 Decline of the remaining fraction of planetesimals during the last 3.5 Gyr of dynamical evolution. The rate of decline slows down with increasing time. Data are from one of the GPU simulations of Wong et al. (2023).

Fig. 3 Semi-major axis and perihelion distance evolution of two particles that survive for the first 1 Gyr. They spend a long (black) and short (orange) time on Neptune-crossing EC orbits. The blue and cyan lines indicate the semi-major axes of Neptune and Uranus, respectively.

3.2 The current number of scattered disc objects and their size-frequency distribution

Based on Outer Solar System Origins Survey (OSSOS) observations the current estimated number of scattered disc objects with diameter D_i > 10 km is (2.0 ± 0.8) × 10⁷ (Nesvorný et al. 2019), assuming their albedo is 6%. We scaled the number of impactors from the observed number of scattered disc objects with D_i. We have, $$N_{\text{SD},>{\text{D}}_{\text{i}}}=N_{\text{scale }}{\left(\frac{D_{\text{i}}}{10 \text{km}}\right)}^{\alpha }.$$(7)

Observations of distant scattered disc objects in the hot or excited population, as well as Jupiter’s Trojan asteroids, show that the slope α seems to have two values, one for objects with diameter D_i ≳ 100 km, and another one for smaller objects (Shankman et al. 2013; Fraser et al. 2014; Lawler et al. 2018; Wong & Brown 2015; Yoshida & Terai 2017). These studies suggest that for D_i ≳ 100 km |α₁| = 4.5 ± 0.5 and for 1 ≲ D_i ≲ 100 km objects |α₂| = 2.0 ± 0.2. Wong et al. (2023) adopts the nominal values α₂ = −2 and α₁ = −4.5. For impactors with diameter D_i < 1 km, we shall adopt the slope of Singer et al. (2019), in other words α₃ = −0.7.

With the nominal adopted value of α2 by Wong et al. (2023) the expected current number of scattered disc objects with diameter D_i > 1 km is approximately 2.0 × 10⁹, which is consistent within uncertainties with the values obtained by previous numerical and observational works (e.g. Brasser & Morbidelli 2013; Volk & Malhotra 2008; Shankman et al. 2013; Lawler et al. 2018). This estimate varies by at least a factor of three: Shankman et al. (2013) compute 2 × 10⁹ for Di > 1 .5 km (assuming an albedo of 4%), followed by Lawler et al. (2018), who increases it to 3 × 10⁹ for the same diameter. Volk & Malhotra (2008) report a number of 10⁹, probably for objects with D_i > 2 km (the exact diameter is not specified), and Brasser & Morbidelli (2013) estimate 2 × 10⁹ for D_i > 2.3 km. It is possible that these studies adopt different size-frequency distributions at that size range, making scaling all of these values to D_i > 1 km or D_i > 10 km challenging. Can we verify whether any of the above estimates for N_SD make sense?

One source of bodies that we can use are the Jupiter-family comets (JFCs). Following Levison & Duncan (1997) we can define the production rate of JFCs as r_SDf_JFC. Literature estimates show that f_JFC = 16–30% (Levison & Duncan 1997; Di Sisto et al. 2009; Brasser & Morbidelli 2013). The usual reported number of active JFCs with diameter D_JFC > 2 km and perihelion distance q < 2.5 au is 117 (Di Sisto et al. 2009). The dormant to active JFC ratio is about 5–6.5 (Levison & Duncan 1997; Di Sisto et al. 2009; Brasser & Morbidelli 2013), so to first order there are of the order of 700 dormant JFCs. These comets spend 7–10% of their dynamical life with q < 2.5 au (Brasser & Morbidelli 2013), so that there are about 10⁴ JFCs with D_JFC> 2 km in total. The median dynamical lifetime of these comets is 165 kyr (Levison & Duncan 1997; Brasser & Morbidelli 2013), so the loss rate is then approximately 0.06 yr⁻¹. This must match the production rate N_SDf_JFC|r_SD| and results in a nominal value of N_SD = 2 × 10⁹ if we take f_JFC = 20% and |r_SD| = 0.15 Gyr⁻¹, with a range of (1.1 - 2.5) × 10⁹ accounting for the uncertainties of f_JFC and |r_SD|.

From their own analysis of numerical simulations of JFC production, Nesvorný et al. (2017) estimate that at present the number of active JFCs to scattered disc objects is 10⁻⁷ – comparable to 6.8 × 10⁻⁸ of Brasser & Morbidelli (2013) – and they compute that N_SD = (1.5 – 4) × 10⁷ with D_i > 10 km, implying N_SD = (1.5 – 4) × 10⁹ with D_i >1 km if α₂ = −2, consistent on the lower end with earlier estimates. However, Nesvorný et al. (2017) assumes a physical lifetime of 500 orbits. In an independent study on near-earth objects including active JFCs Bottke et al. (2002) obtains N_SD = (2.8 ± 2.0) × 10⁹ objects beyond Neptune with D_i > 1 km. However, most of the estimates for N_SD derived from the JFC population pertain to objects with diameter D_i > 2 km rather than 1 km, which agrees better with the estimates of Brasser & Morbidelli (2013); Shankman et al. (2013); Lawler et al. (2018) than with that calculated by us from the observational constraints by Nesvorný et al. (2019) and the numerical work by Nesvorný et al. (2017) for objects with D_i > 10 km.

A potential problem presents itself when we extrapolate the estimate of the current population back in time before the onset of giant planet migration. The average depletion of the population in the clone simulations is 99.8%, comparable to the 99.7% of Nesvorný et al. (2017). If the current population with D_i > 1 km is 2 billion, the original number was about 10¹², which is higher than that computed by Nesvorný et al. (2023) and Wong et al. (2023) by a factor of 1.5, which is acceptable. But if the current population is 3 or 4 billion then it is difficult to reconcile with the original population unless its size-frequency distribution slopes are on the very steep end of what observational constraints allow.

One way out of this dilemma is to accept that the current population of scattered disc objects with diameter D_i > 10 km is (2 ± 0.8) × 10⁷ (Nesvorný et al. 2019) so that with our nominal depletion the nominal primordial population of such objects is 10¹⁰ – Nesvorný et al. (2017) reports (8 ± 3) × 10⁹. With the primordial disc mass of 18 M_⊕ from Wong et al. (2019, 2023), and adopting the knee size frequency distribution of trans-Neptunian objects from Fraser et al. (2014); Lawler et al. (2018) a primordial number of objects of 10¹⁰ with D_i > 10 km requires that the size-frequency distribution slopes are either (α_1; α₂) = (5.0, 2.1) or (α₁, α₂) = (4.5, 2.2). In either case one of the slopes is at the higher end of literature estimates.

With these faint-end slopes and our depletion amount the current number of scattered disc objects with D_i > 1 km becomes N_SD = (2.5 – 3.2) × 10⁹. In what follows we adopt a nominal value of N_SD = 3 × 10⁹ for D_i > 1 km. The resulting nominal number of ecliptic comets with D_i > 1 km is then N_EC = 5.4 × 10⁷ if we use the outcome of the GENGA simulations.

We now compute the ecliptic comet injection rate. From the RMVS clone simulations we obtain $F_{\text{EC}}={0.366}_{-0.082}^{+0.118}$ yr⁻¹ (2σ), while the combined GENGA simulations yield F_EC = 0.353 ± 0.002 yr⁻¹ (1σ) these estimates do not include the uncertainty in N_SD. We resort to using the outcome of the GPU clone simulations, and we adopt a nominal value of F_EC = 0.36 yr⁻¹. With this approach the largest uncertainty in the impact rate is that from the number of scattered disc objects. In Table 3 we have listed the nominal values of all the relevant parameters that we adopted.

3.3 Impact probabilities with the giant planets and their satellites

The various impact probabilities of planetesimals with the giant planets and the satellites, as well as encounter probabilities with the planets, and impact speeds with the satellites are listed in Tables 4 and 5. The quantities derived from the numerical simulations, as well as their uncertainties, are given in the Supplementary Materials on Zenodo.

The first table lists the impact and encounter probabilities of planetesimals with the giant planets. The first two entries are obtained from the clone simulations run with GENGA for the last 3.5 Gyr (Wong et al. 2023). The third column is P_P for the RMVS clone simulations. The next two columns are the same values as the first two, but are obtained from Wong et al. (2021) when the giant planets were migrating. The last column is the impact probability obtained from the RMVS M simulations. Even though the output from these last simulations is not used in our analysis, we list the impact probabilities here because our values are very similar to those of Levison et al. (2000) and these impact probabilities were used by Zahnle et al. (1998), and Zahnle et al. (2003) in their analysis of the current cratering rates.

Two things are of note. First, the impact probabilities with the planets after migration during the last 3.5 Gyr are lower than that during the migration phase, typically by a factor of 3–4. The reason is the different initial conditions of the migration and clone simulations. During the migration phase the planetesimals are located in a ring surrounding the outermost ice giant. They are excited by the migrating planets and are quickly scattered from one planet to another. During the last 3.5 Gyr most plan-etesimals from this initial ring are either ejected or collide with one of the giant planets. Planetesimals remaining after 1 Gyr of migration reside on meta-stable orbits (Duncan & Levison 1997) whose distribution is shown in Fig. 1, and their time-averaged orbital distribution is very different from the first 1 Gyr. The different initial conditions are also the reason why for Uranus P_enc/P_P < 40 in the M simulations. The reason is that during their migration a significant fraction of planetesimals encounter the ice giants on elliptic orbits, some of which are still recorded in the code as encounters if they venture close enough. For comparison, Zahnle et al. (1998) reports that 20% of planetesimals encounter Jupiter on elliptical orbits.

Second, the impact probabilities we computed with RMVS are systematically higher than with GENGA. The reason for this is likely to be the time step: we chose its value as a compromise between speed and accuracy, with the understanding that most planetesimals were ejected due to the cumulative effect of mostly distant encounters with the giant planets. However, with a large enough time step RMVS does not adequately resolve very close repeated encounters with the giant planets (Levison & Duncan 1994; Michel & Valsecchi 1996), while GENGA, using a Bulirsch-Stoer integrator close to the planet, likely does a better job of resolving such close encounters. The impact probability difference between RMVS and GENGA is a factor of 1–1.5, which is propagated into the impact rates. In what follows we adopt the impact probabilities from the GENGA simulations.

Table 5 contains the impact probabilities of planetesimals with the satellites. The first entry, P_s/P_P, is the impact probability scaled to that of the planet obtained from our flyby simulations, where we counted the number of impacts on a satellite and divided by the number that hits the planet. The next entry, N_i/N_tot, is the impact probability of a planetesimal that ventures within 40 R_P as computed from our flyby simulations This is the number of planetesimals that collide with a satellite divided by the total number of planetesimals (500 million) in the flyby simulations. The next entry, P_s is the absolute impact probability of a heliocentric planetesimal with the satellites, and is computed as N_i/N_tot × P_enc, where we used the P_enc values on the planets from the GENGA simulations during the last 3.5 Gyr. As an example, for Enceladus N_i/N_tot = 0.10 × 10⁻⁶, suggesting that typically we recorded 50 impacts out of 500 million particles in the flyby simulations, and P_enc = 5.95% for Saturn so that the absolute impact probability becomes P_s = 0.1 × 10⁻⁶ × 0.060 = 0.60 × 10⁻⁸, which is a little lower than using the relation P_s/P_P × P_P. The discrepancy between these two values of the absolute impact probability is largest for Mimas and Enceladus. This could be attributed to the flyby simulations missing some encounters and impacts due to the high encounter speeds with these tiny bodies, but several iterations of the flyby simulations with higher values of n₂ > 1 did not remedy that problem. We subsequently list the average impact speed with the satellite obtained from the flyby simulations, followed by the impact rate of objects with D_i > 1 km, and the crater diameter that such an impact produces. The last entry is cratering rate for craters with diameter D_cr > 10 km assuming the crater scaling relations discussed below.

We note that the differences between P_s/P_P and with direct hits during the flyby simulations are typically 30% and are smallest for the Jovian satellites Io, Europa, and Ganymede (about 10%). We do not know what accounts for this large difference. Both the analytical estimate of Zahnle et al. (2003) and our numerical flyby simulations have their own assumptions and limitations, but Zahnle et al. (2003) did say that their Monte Carlo procedures gave results that were 30% different from their analytical solutions. Additional differences could result from numerical factors in the flyby simulations and the fact that the gravitational focusing of the giant planets becomes important for Saturn and Uranus at closer distances than for Jupiter, resulting in a different initial distribution of the impact parameter. The uncertainty in the value of P_P is at least a factor of 1.3.

3.4 Current impact rate on the satellites

We computed the current impact rate on a satellite (or a planet) following the procedure from Duncan et al. (1995) and Levison et al. (2000) based on the injection rate (F_EC) of ecliptic comets with D_i ≥ 1 km originating from beyond Neptune, and on the impact probability with the satellites with the giant planets. The formula for calculating the impact rate is $${\dot{C}}_{\text{s}}=F_{\text{EC}}{P}_{\text{s}}$$(8)

where P_s is the post-migration collision probability onto a satellite as described in Equation (1).

The impact rates listed in the fifth column of Table 5 need to be converted into a cratering rate using crater scaling laws to convert impactor diameter to craters, and vice versa. Wong et al. (2023) employs a crater scaling law based on the π-scaling law, using the simple-to-complex crater diameter relation D_SC = 4 (1 ms⁻²/ɡ) km (Schenk 1991), which applies to all icy satellites apart from Mimas, Enceladus, Tethys, Dione, and Miranda, for which D_SC ~ 15 km (Schenk 1991). For Io, D_SC = 24 (1 m s⁻²/ɡ) km (Pike 1980). The crater scaling law in Wong et al. (2023) for icy bodies is $$\begin{array}{l}{D}_{\text{cr}}=4.66{\left(\frac{D_{\text{i}}}{1 \text{km}}\right)}^{0.917}{\left(\frac{{\rho }_i}{{\rho }_s}\right)}^{0.392}{\left(\frac{v_{\text{imp}}}{1 \text{km} {\text{s}}^{-1}}\right)}^{0.517}{\left(\frac{𝑔}{1 \text{m} {\text{s}}^{-2}}\right)}^{-0.082}\\ D=3.67{\left(\frac{D_{\text{i}}}{1 \text{km}}\right)}^{0.917}{\left(\frac{{\rho }_i}{{\rho }_s}\right)}^{0.392}{\left(\frac{v_{\text{imp}}}{1 \text{km} {\text{s}}^{-1}}\right)}^{0.517}{\left(\frac{𝑔}{1 \text{m} {\text{s}}^{-2}}\right)}^{-0.259}\end{array}$$(9)

where the bottom equation applies to Mimas, Enceladus, Tethys, Dione and Miranda, and the top equation to the other satellites (for Io the coefficient in the top equation is 3.39). For both crater scaling laws we assume that the impact angle is 45°. These scaling laws are somewhat different from those in Zahnle et al. (2003). Wong et al. (2023) assumes an impactor density of ρ_i = 400 kg m⁻³, which is the mean density of comet Tempel 1 (Richardson et al. 2007) and which is a little lower than that of Churyumov-Gerasimenko (Jorda et al. 2016), and a satellite density ρ_s = 1000 kg m⁻³ apart from Io, for which ρ_s = 3520 kg m⁻³. We employ values for the traditional scaling parameters α = 0.85, which describes the simple-to-complex crater diameter ratio (Croft 1985), and β = 0.22, which relates the efficiency of crater excavation in the gravity-dominated regime. The inverse relations for icy bodies are $$\begin{array}{l}{D}_{\text{i}}=0.187{\left(\frac{D_{\text{cr}}}{1 \text{km}}\right)}^{1.09}{\left(\frac{{\rho }_i}{{\rho }_s}\right)}^{-0.427}{\left(\frac{v_{\text{imp}}}{1 \text{km} {\text{s}}^{-1}}\right)}^{-0.564}{\left(\frac{𝑔}{1 \text{m} {\text{s}}^{-2}}\right)}^{0.09}\\ D_{\text{i}}=0.241{\left(\frac{D_{\text{cr}}}{1 \text{km}}\right)}^{1.09}{\left(\frac{{\rho }_i}{{\rho }_s}\right)}^{-0.427}{\left(\frac{v_{\text{imp}}}{1 \text{km} {\text{s}}^{-1}}\right)}^{-0.564}{\left(\frac{𝑔}{1 \text{m} {\text{s}}^{-2}}\right)}^{0.282},\end{array}$$(10)

from which we compute the cratering rate for D_cr > 10 km using the inferred impactor size-frequency distribution slope of α₃ = −0.7 for impactors with D_i < 1 km (Singer et al. 2019), and following Nesvorný et al. (2019) by adopting α₂ = −2.1 for D_imp > 1 km. The uncertainties in the cratering rate are a factors of a few, primarily due to uncertainties in F_EC (including in N_SD). For Io the coefficient of the top equation becomes 0.264.

The nominal value of F_EC and the current impact probability of Jupiter result in a rough estimation of the impact rate on Jupiter, which we compute as 1.0 × 10⁻³ per year for objects with D_i > 1 km. Our estimated Jovian impact rate falls towards the low end, but remains consistent with values from previous studies. The ranges and uncertainty of impact rates in Zahnle et al. (2003) are discussed in Section 4, but suffice to say that the range of these impact rates for D_i > 1 km spans two orders of magnitude, and the uncertainties in these estimates also vary. We also calculated the impact and cratering rates on Enceladus (and the other satellites of all three planets). The impact rate of objects with D_i > 1 km with Enceladus is computed to be 2.12 Gyr⁻¹. For the cratering rate we used the inverse crater scaling law. With our adopted values for the impactor density, satellite density, and impact velocity from Table 5, a crater with diameter D_cr > 10 km on Enceladus requires a projectile diameter of D_i = 0.41 km and adopting the size-frequency distribution slope α₃ = −0.7 from Singer et al. (2019) the current cratering rate becomes 5.0 × 10⁻⁶ km⁻² Gyr⁻¹.

4 Comparison with earlier work

We present three sets of studies that estimate current impact rates and surface age for the icy satellites. These sets include:

• Classic studies (Zahnle et al. 1998, 2003; Kirchoff & Schenk 2009),

• Recent studies by others (Nesvorný et al. 2019, 2023; Bottke et al. 2023, 2024), and

• Recent studies by us (Wong et al. 2019, 2021, 2023), and this work.

Each group of studies represents a series of works that begin with planetary dynamics and impact rate estimations, progressing towards determining crater ages for the icy satellites. These works employ various methodologies and sources, ranging from numerical simulations of planetary migration to observational crater counts. It is important to note that no single paper in any of these sets comprehensively addresses the entire spectrum from planetary dynamics to geological history.

We shall discuss the differences in results from the three studies by listing the impact rates on Jupiter and Enceladus, and the cratering rate on the latter, in Table 6. In the following, we address the large uncertainties and differences between our impact rate estimates, and those from previous studies.

4.1 Classic studies

Zahnle et al. (1998, 2003) constrain the current impact rate on Jupiter through multiple measurements. Zahnle et al. (2003) reports an impact rate of approximately 4 × 10⁻³ yr⁻¹ based on six encounters (smallest ~1 km in diameter) within four Jovian radii over ~350 years. Lamy et al. (2004) observe nine comets crossing Callisto’s orbit over approximately 50 years, suggesting an impact rate of about 1.75 × 10⁻³ per year for D_i > 1 km. Other estimations, using recorded collisions and close encounters, measurements of excess carbon monoxide content, crater counting on Ganymede, and inference from Near Earth Objects, range from 0.4 to 44 × 10⁻³ impacts per year. The ranges Zahnle et al. (2003) discusses are: historical records (4.3 × 10⁻³ to 2.3 × 10⁻² impacts per year), observed impacts (1.4 × 10⁻² to 4.4 × 10⁻² yr⁻¹), excess carbon monoxide content (7.7 × 10⁻³ yr⁻¹), crater counting on heavily cratered terrain on Ganymede (2.0 × 10⁻³ to 4.8 × 10⁻² yr⁻¹), and inferred impacts from Near Earth Objects (1.3 × 10⁻³ yr⁻¹). All impact rates above are scaled to impactors with diameters D_i > 1 km, assuming a cumulative power-law slope of α₂ = −2 (Singer et al. 2019). In any case, the impact rate on Jupiter of Zahnle et al. (1998, 2003) have uncertainties that vary significantly, ranging from a factor of three to an order of magnitude. Zahnle et al. (2003) calibrate the impact rate on Jupiter at $5_{-3}^{+6}\times {10}^{-3}$ yr⁻¹ for impactors with D_i > 1 .5 km.

This commonly cited Jovian impact rate for D_i ≥ 1 .5 km, when translated using their Case A size-frequency distribution, corresponds nominally to 10⁻⁴ yr⁻¹ for impactors with diameter D_i ≥ 10 km, and 7.5 × 10⁻³ yr⁻¹ for D_i ≥ 1 km. Using α₂ = −2.1 we find similar results: 9.3 × 10⁻⁵ yr⁻¹ for D_i ≥ 10 km, and 1.0 × 10⁻² yr⁻¹ for D_i ≥ 1 km. We note that Case A of Zahnle et al. (2003) represents a shallower size-frequency distribution of the impactors, with α₃ = −1 for D_i < 1.5 km, whereas their Case B reflects a steeper distribution closer to a collisionally evolved population deduced from craters on Triton (α₃ = −2.5 for D_i < 1.5 km). Most recent studies adopt Case A as a more accurate representation of heliocentric comet populations in the outer Solar System, so that we shall do the same. The shallow slope is also more closely aligned with the results of Singer et al. (2019).

Zahnle et al. (2003) theoretically assess the impact rate on Saturn by considering several factors: the proportion of comets delivered from beyond Neptune’s orbit to Jupiter’s, the impact probability based on the planets’ relative sizes to their Hill spheres, and the gravitationally enhanced cross-sections, which are proportional to each planet’s mass and escape velocity. From this, they estimate the impact ratio of Jupiter to Saturn to be 1:0.42. This is close to the ratio of 1:0.47 from the average annual Öpik impact probabilities for the three then-known Centaurs with Saturn, for which the impact rate is computed as 2 × 10⁻⁸ yr⁻¹ for a 150 km diameter object (Fernández et al. 2002). Together, this corresponds to a Saturnian impact rate of 3.2 × 10⁻³ yr⁻¹ for D_i ≥ 1 km, using the Case A size-frequency distribution.

With this information we compute that the current impact rate on Enceladus, scaled from Jupiter and using Case A, is approximately 16.7 Gyr⁻¹ for objects with D_i > 1 km. We used P_s/P_P = 2.2 × 10⁻⁶ for the Enceladus to Jupiter impact probability ratio, and that the impact rate on Jupiter is 5 × 10⁻³ yr⁻¹ for D_i ≥ 1 .5 km. This impact rate on Enceladus is a factor of eight higher than our estimate, and is primarily due to Zahnle et al. (2003) taking a factor of six higher impact rate on Jupiter, as well as our different size-frequency distribution between impactor diameters 1 km and 1.5 km. To compute the cratering rate on Enceladus we use the inverse of the crater scaling laws adopted by Zahnle et al. (2003) (see the Appendix). Using the impact velocity of heliocentric planetesimals with Enceladus listed by Zahnle et al. (2003) of υ_i = 24 km s⁻¹ the impactor diameter required to excavate a crater with diameter D_cr > 10 km on Enceladus is D_i = 0.34 km. Adopting the Case A size-frequency distribution implies the current impact rate of such planetesimals on Enceladus is 49 Gyr⁻¹. The resulting cratering rate is then 6.1 × 10⁻⁵ km⁻² Gyr⁻¹, which is a factor of 12 greater than our estimate, with the additional factor of two over the impact rate being caused by our different adopted crater scaling law and slope.

4.2 Recent studies

In the last five years, Nesvorný et al. (2023) have built upon their earlier works of Nesvorný (2015a,b); Nesvorný et al. (2017) to further constrain the initial conditions of the giant planet migration model, specifically the initial semi-major axis of Neptune and its migration pathway. The most recent version of the model assumes that the dynamical instability among the giant planets occurred 10 Myr after the dispersal of the protoplanetary gas disk. The migration rate of Neptune follows an exponential profile (Nesvorný et al. 2017): for the first 10 Myr, the migration e-folding time (τ) is 10 Myr, representing rapid migration. Then, after 10 Myr of dynamical evolution, τ increases to 30 Myr, reflecting a slower migration rate, which lasts up to 500 Myr.

Nesvorný et al. (2023) uses the initial conditions of the simulations of Nesvorný et al. (2017) and reran the last 1 Gyr of evolution. They begin with one million test particles from Nesvorný et al. (2017) and apply ‘on-the-fly’ cloning when any particle ventures closer than 23 AU for the first time, cloning it 50 times to maintain sufficient particle interactions with the planets.

From the simulation’s impact counts, Nesvorný et al. (2023) calculates the impact rate on Jupiter as 3.5 × 10⁻⁵ per year for objects with D_i ≥ 10 km, without considering cometary disruption. They compute the impact rate as Ċ_p = (Δn/t) (n_imp/n_tp), where T is 1 Gyr, and ΔN is essentially equal to the original number of planetesimals in the disc because 99.7% are lost (Nesvorný et al. 2017). In Nesvorný et al. (2023) they adopt ΔN = (8 ± 3) × 10⁹ with D_i > 10 km and they record N_imp = 217 impacts on Jupiter out of a total of N_tp = 5 × 10⁷ test particles; in this manner they obtain their nominal value Ċ_p = 3.5 × 10⁻⁵ per year. With cometary disruption included, the rate decreases to 2.6 × 10⁻⁵ yr⁻¹. Using their adopted impactor power-law slope of α₂ = −2.1, we computed that their impact rate on Jupiter for D_i > 1.5 km becomes 1.9 × 10⁻³ yr⁻¹ without disruption, and 1.4 × 10⁻³ yr⁻¹ with disruption. These values are still higher than ours by factors of a few, but they also fall outside the range of $5_{-3}^{+6}\times {10}^{-3}$ yr⁻¹ of Zahnle et al. (2003). The impact rates for objects with D_i > 1 km are listed in Table 6.

To calculate the collision probability with Enceladus, Nesvorný et al. (2023) record the close encounters within a Hill radius of Saturn and numerically evaluate the collision probability for each Saturnian encounter with Enceladus. Following Zahnle et al. (2003) the impact rate was then scaled to that of Jupiter. Using their nominal impact rate on Jupiter of 3.5 × 10⁻⁵ yr⁻¹ for D_i ≥ 10 km and P_s/P_p = 1.5 × 10⁻⁶ as the Enceladus to Jupiter impact probability ratio – see Table 2 of Nesvorný et al. (2023) – the impact rate on Enceladus becomes 0.053 Gyr⁻¹ for D_i > 10 km, and 6.6 Gyr⁻¹ for D_i > 1 km using α₂ = −2.1; this impact rate is a factor of 2.5 higher than our value.

Nesvorný et al. (2023) do not explicitly compute the crater-ing rate onto the satellites, only the relative impact probabilities. Since some of their methodology follows that of Zahnle et al. (2003) we can naively adopt the same cratering scaling law as Zahnle et al. (2003) and the size-frequency distribution of their Case A for objects with diameter D_i < 1 km. We then computed the cratering rate on Enceladus as 2.4 × 10⁻⁵ km⁻² Gyr⁻¹ for the non-disruption case, and 2.2 × 10⁻⁵ km⁻² Gyr⁻¹ with disruption. These cratering rates are a factor of three to four higher than ours and a factor of three to four lower than that of Zahnle et al. (2003).

5 Discussion

In this section we discuss our results in comparison to previous studies in more detail. To begin, we find that the impact rates onto Jupiter and Enceladus, and by extension the other giant planets and their satellites, is lower than that found by two other studies. We made our methodology transparent, we clarified every choice of every parameter and how they were used. It is strange that despite many of these parameters being similar to some of the previous studies, we still find a very different impact rate. Why is this so? There are several reasons that we can think of.

The first is the outcome of the numerical simulations. The migration simulations in Wong et al. (2023) move Neptune and Uranus out slower than Nesvorný et al. (2017), but Nesvorný et al. (2023) concluded that the details of the migration should not matter for the current impact rates. We concur with this statement based on logical deduction alone. However, Wong et al. (2023) find that the location of the planetesimal disc’s outer edge appears to play a role because it influences the depletion, fraction that becomes ecliptic comets, and the rate of decay, but only when the outer edge is beyond 35 au. Furthermore, the migration occurs very early and the outer Solar System is relatively stable for the past 3.5 Gyr. Therefore the current impact rate – which reflects only the most recent few hundred Myr – should not strongly depend on the specific details of Neptune’s migration, such as its speed or timing if it occurred more than 4 Gyr ago. The adopted disc mass and its width by Wong et al. (2023) are also comparable to Nesvorný et al. (2017), and we obtained the same amount of depletion by 4.5 Gyr of evolution, apart from cases where the disc extended beyond 35 au. The average rate of decline of our scattered disc population is consistent with other works in the literature. If this parameter is fundamentally wrong in our simulations it has to be equally wrong in others; this is highly unlikely. Unfortunately we do not have the numerical resolution to determine f_JFC , and other studies do not report f_EC so that we cannot compare our values to the works of others; the only parameter with which we can do this is r_SD, which agrees across the board.

One difference between the outcome of our numerical simulations, that of Levison & Duncan (1997), and Levison et al. (2000) is that our impact probability with the planets is lower than in these studies, even though the impact probability for the simulations with migration are comparable to those in Levison et al. (2000). Levison & Duncan (1997) do not explicitly report the impact probabilities, yet the simulations of Levison et al. (2000) also find a comparable lifetime for the ecliptic comets than we do, once again lending credibility to our approach. We argued why our clone simulations are a good proxy for the current population and why the impact probabilities from this population should be considered rather than that including the migration; simply put: the migration phase is not a good proxy for the current solar system. In our studies, we derive the impact probability probabilities of the satellites by multiplying the impact count during the flyby simulations with the encounter probability with the planets. Both counts were recorded in dynamic simulations. Nesvorný et al. (2023), however, numerically evaluates the collision probability of each planetesimal encounter with the giant planets and its satellites, yet we find similar impact probability ratios P_s/P_P to Zahnle et al. (2003) and Nesvorný et al. (2023), so that our different approach still yields consistent results.

A second potential reason for disagreement is the current number of scattered disc objects that we use. Neither Zahnle et al. (2003) nor Nesvorný et al. (2023) explicitly mention this number because they have no need for it. Unfortunately this quantity is poorly known and it will probably remain so for the time being. We made clear why we have chosen the number that we did, based on constraints from the total depletion of the population of planetesimals in our simulations as well as constraints on the primordial disc mass and the observed slopes of the size-frequency distribution of these objects, and the observational constraints. We do not foresee that the current population can be much more or much less numerous than what we have used here; factors of two perhaps. As such, this quantity is perhaps also not to blame.

A third possibility of disagreement is the choice of impactor size-frequency distribution. Our choices are certainly different from those of Zahnle et al. (2003) and partially accounts for our lower impact rates, but we deliberately choose the same slope of α₂ = −2.1 as Nesvorný et al. (2023) for objects with diameter 1 ≲ D_i ≲ 100 km to minimise differences in outcome. We extended the population to sub-km objects by adopting α₃ = −0.7 from Singer et al. (2019), but if we restrict ourselves to computing the impact rate of km-sized objects then this choice does not matter. Thus, our difference of a factor of a few with Nesvorný et al. (2023) cannot be explained by the choice of impactor size-frequency distribution slope.

Beyond these relatively obvious potential differences that all yield comparable outcomes, there is one additional divergence in impact rate calculation methods that could contribute to the observed discrepancies in our results: the fundamental difference in approach. Both Zahnle et al. (2003) and Nesvorný et al. (2023) anchor their satellites’ impact rate to that of Jupiter. Whether this is derived from simulations or from observational constraints is less relevant. The satellite-to-Jupiter impact probability ratio and Jupiter’s absolute impact rate yield the resulting satellites’ impact rate. While the impact probability with satellites is implicit in this method, it is indirectly informed by the impact probability with Jupiter, even if this is not explicitly stated in Zahnle et al. (2003).

In contrast, our approach uses the flux of new comets entering the Centaur region from beyond Neptune as the fundamental parameter. This is then multiplied by the absolute impact probability with the satellites – derived from flyby simulations and the encounter probabilities with the giant planets – to yield the satellites’ impact rate. This is where the major differences show up: our impact rate on Jupiter is a factor of eight lower than that of Zahnle et al. (2003), and at least three lower than that of Nesvorný et al. (2023). The discrepancy with Zahnle et al. (2003) is primarily due to their choice of fixing the impact rate on Jupiter from historical observations rather than numerical simulations, while the discrepancy with Nesvorný et al. (2023) is more difficult to explain, but it appears to be due to how they compute the impact rate on Jupiter over the last billion years for objects with D_i > 10 km, and their potentially anomalously low impact probability with this planet.

As a thought experiment we may use the impact rate on Jupiter adopted by Zahnle et al. (2003) at face value and compute what the injection rate of comets with D_i > 1 km from beyond Neptune needs to be to be consistent with this impact rate. We obtain F_EC = 2.6 yr⁻¹. For this high value there is no literature precedent. Even if we were to adopt an impact probability with Jupiter of about 1% (Levison et al. 2000) rather than our value of 0.29%, the injection rate still needs to be 0.75 yr⁻¹, which is still very high. Repeating this check for the impact rate of Nesvorný et al. (2023) with our impact probability and without disruption results in 1.5 yr⁻¹, for which there is also no literature precedent either; using the higher impact probability on Jupiter of 1% lowers this to 0.44 yr⁻¹, which is acceptable, but this cannot be verified without knowing their impact probability with Jupiter.

One potential effect that we did not consider is tidal disruption of comets. From crater chains on Ganymede and Callisto, and assuming their surfaces ages are about 4 Gyr old, Schenk et al. (1996) argue that there are 3.7 × 10⁻³ tidal disruptions per year near Jupiter. Since these disruptions generally happen within 3 R_J for comets with D_i ~ 1 km this rate of disruption is fairly consistent with our estimated impact rate on Jupiter, yet with typically 10 fragments generated these disruptions are insufficient to elevate the impact rate at Jupiter by a factor of a few.

6 Conclusions

In this work, we uses numerical simulations of the evolution of the outer Solar System to compute the impact and cratering rates onto the regular satellites of Jupiter, Saturn, and Uranus, highlighting the complexities involved in determining these rates. We demonstrated how the current impact rate depends strongly on the specific parameters of the dynamical and disc models employed.

We outlined our methodology and parameter choices in detail, allowing future studies to verify or adjust these assumptions and recalculate impact rates. After careful comparison with previous works, we conclude that our impact rates are lower than that the best estimates of previous studies by up to a factor of six, likely due to our reliance on the injection rate of scattered disc objects into the Centaur region rather than anchoring to the impact rate on Jupiter. This difference in both approach and outcome highlights the importance of clearly stating all parameters and assumptions in future studies to enable meaningful comparisons.

Our findings underscore that the uncertainties in impact rates substantially influence surface age estimates, particularly for young, sparsely cratered terrains such as Europa, Enceladus, and young patches on Tethys and Dione. Better constraints on orbital dynamics and impactor size-frequency distributions are needed to refine both current and historical impact rates.

Future advances, such as more precise measurements of small body sizes using stellar occultation techniques, can provide crucial data to constrain the current impact rates and improve our understanding of the historical evolution of small bodies in the outer Solar System. These insights will help bridge gaps between impact dynamics and geological interpretations of icy moon surfaces.

Data availability

Supplementary materials are available on https://doi.org/10.5281/zenodo.14974372

Acknowledgements

This study is supported by the Research Council of Norway through its Centres of Excellence funding scheme, project No. 332523 PHAB. GENGA can be obtained from Bitbucket.

Appendix A Crater scaling

Zahnle et al. (2003) uses the following crater scaling law $$\begin{array}{l}{D}_{\text{cr}}=4.38{\left(\frac{D_{\text{i}}}{1 \text{km}}\right)}^{0.783}{\left(\frac{{\rho }_i}{{\rho }_s}\right)}^{0.333}{\left(\frac{v_{\text{imp}}}{1 \text{km} {\text{s}}^{-1}}\right)}^{0.434}{\left(\frac{𝑔}{1 \text{m} {\text{s}}^{-2}}\right)}^{-0.217}\\ D_{\text{cr}}=5.31{\left(\frac{D_{\text{i}}}{1 \text{km}}\right)}^{0.885}{\left(\frac{{\rho }_i}{{\rho }_s}\right)}^{0.376}{\left(\frac{v_{\text{imp}}}{1 \text{km} {\text{s}}^{-1}}\right)}^{0.490}{\left(\frac{𝑔}{1 \text{m} {\text{s}}^{-2}}\right)}^{-0.245}\\ \times {\left(\frac{D_{\text{SC}}}{1 \text{km}}\right)}^{-0.13}\end{array}$$

where ɡ is the acceleration due to gravity, ρ_i and ρ_t are the densities of the impactor and the satellite, and D_SC is the simple to complex transition crater diameter. Zahnle et al. (2003) uses D_SC = 2.5 km for Europa, Ganymede, Callisto, and Titan (Schenk et al. 2004) and D_SC = 15 km for all other satellites as well as an impactor density of ρ_i = 600 kg m⁻³ and a satellite density of ρ_s = 900 kg m⁻³ for all satellites apart from Io, for which ρ_s = 2700 kg m⁻³. The top equation should be applied when D_cr < D_SC. The inverse equations are then $$\begin{array}{l}{D}_{\text{i}}=0.151{\left(\frac{D_{\text{cr}}}{1 \text{km}}\right)}^{1.28}{\left(\frac{{\rho }_i}{{\rho }_s}\right)}^{-0.425}{\left(\frac{v_{\text{imp}}}{1 \text{km} {\text{s}}^{-1}}\right)}^{-0.554}{\left(\frac{𝑔}{1 \text{m} {\text{s}}^{-2}}\right)}^{0.277}\\ D_{\text{i}}=0.151{\left(\frac{D_{\text{cr}}}{1 \text{km}}\right)}^{1.13}{\left(\frac{{\rho }_i}{{\rho }_s}\right)}^{-0.425}{\left(\frac{v_{\text{imp}}}{1 \text{km} {\text{s}}^{-1}}\right)}^{-0.554}{\left(\frac{𝑔}{1 \text{m} {\text{s}}^{-2}}\right)}^{0.277}\\ \times {\left(\frac{D_{\text{SC}}}{1 \text{km}}\right)}^{0.147}\end{array}$$

References

All Tables

All Figures

	Fig. 1 Orbital distribution of semi-major axis and perihelion distance of planetesimals after 1 Gyr of evolution of giant planet migration. The colour coding indicates inclination. Data from Wong et al. (2023).
In the text

	Fig. 2 Decline of the remaining fraction of planetesimals during the last 3.5 Gyr of dynamical evolution. The rate of decline slows down with increasing time. Data are from one of the GPU simulations of Wong et al. (2023).
In the text

	Fig. 3 Semi-major axis and perihelion distance evolution of two particles that survive for the first 1 Gyr. They spend a long (black) and short (orange) time on Neptune-crossing EC orbits. The blue and cyan lines indicate the semi-major axes of Neptune and Uranus, respectively.
In the text