© The Authors 2025

1 Introduction

Jupiter trojans (JTs or trojans hereafter) are populations of small asteroids that share the orbit of Jupiter around the Sun. At a mean distance of about 5.2 au, JTs occupy two stable regions around the L4 and L5 Lagrange equilibrium points of the restricted three-body problem in the rotating system of Jupiter around the Sun. One of the most intriguing aspects of JT is the puzzle behind its origin. This population is believed to have formed in the early Solar System, with its current orbits being the result of a complex interplay in its gravitational interaction with Jupiter and other giant planets (Nesvorný et al. 2018; Morbidelli & Nesvorný 2020).

Various mechanisms have been proposed regarding the origin of the trojan population (see Emery et al. (2015) for a review). One of the prevailing theories explaining the origin of JTs involves the capture of planetesimals during the giant planets instability phase (Nesvorný et al. 2018; Bottke et al. 2023a). One particularly successful hypothesis is the ’jumping capture’ model proposed by Nesvorný et al. (2013). Such a model builds upon the ’chaotic capture’ model previously proposed by Morbidelli et al. (2005). According to Nesvorný et al. (2013), Jupiter orbit was scattered multiple times by encounters with an ice giant, causing the planet Lagrangian regions to shift in semi-major axis (Nesvorný & Morbidelli 2012). As a result, these regions eventually covered a populated area of planetesimals, capturing them in the process. This model not only reproduces the orbital distribution of JTs, but also offers a potential explanation for the asymmetric abundances of the leading and trailing populations confirmed by Uehata et al. (2022). An alternative model proposed by Pirani et al. (2019a,b) suggests that Jupiter formed at a distance of about 20 au or further and migrated to its current location while growing in mass, sweeping up planetesimals, and capturing JT in the 1:1 resonance. This model accounts for the asymmetry in the size of the two trojan swarms, but the captured bodies do not match the observed masses or inclinations of the observed trojans. The trojans will be visited by the Lucy space mission (NASA) to help unravel the origin of such an intriguing population (Levison et al. 2021).

Jupiter trojans have experienced significant collisional evolution since their capture, evidenced by the presence of nine collisional families in the L4 swarm and four in the L5 swarm, according to a recent update by Vokrouhlický et al. (2024). Although previous studies investigated the collisional dynamics of JTs, many were conducted before the publication of the capture model (i.e. Marzari et al. 1997; de Elía & Brunini 2007). Consequently, the initial conditions assumed previously, including the initial size distribution of the population and the physical parameters governing collisional responses, should be revisited. To address this, recent studies such as those of Marschall et al. (2022) and Bottke et al. (2023b) have examined the collisional evolution of JTs within the context of the giant planet instability scenario. Using the Boulder code (Morbidelli et al. 2009), both investigations modelled the evolved size-frequency distribution (SFD) of trojans. Marschall et al. (2022) revealed that an unconventional weak collisional strength, with a minimum at 20 m, could ensure the survival of satellites like Queta (3548 Eurybates), over the age of the Solar System. Lately, Bottke et al. (2023b) suggested that an even weaker scaling law is necessary to match the SFD of a trojan population originating from the primordial outer belt, located just behind Neptune when the planet started migrating.

When investigating the collisional evolution of JTs, it is important to consider that its initial population may have been eroded early, before it left the primordial outer belt (in the framework of JT captured by Jupiter during the giant planet instability phase). Benavidez et al. (2025; B25 hereinafter) employed ALICANDEP-22 to model the collisional evolution of the primordial belt and the subsequent evolution of the Kuiper Belt populations, including insights from recent findings by the Outer Solar System Survey (OSOSS) (Petit et al. 2023). The OSOSS analysis showed a striking similarity in the shape of the SFD for the dynamically hot and cold components of the Kuiper Belt within a specific magnitude range H_r ≈ 5.3 to 8.3 (Kavelaars et al. 2021; Petit et al. 2023). As proposed by these authors, this hints at a potential common formation mechanism in the corresponding diameter range of 100 < D < 400 km. Simultaneously, Lyra et al. (2023) highlighted the significance of poly-disperse pebble accretion following streaming instability in planetesimal formation, indicating that such a process led to the accretion of bodies up to several hundred kilometres in diameter at heliocentric distances reaching the current location of the cold Kuiper Belt (CKB). B25 blended such improvements and demonstrated that under suitable conditions, it is feasible to reproduce the present-day size-frequency distribution of both the hot and cold components of the Kuiper Belt.

Building upon the work of B25, we propose to continue their study in this paper. Here, we model the collisional evolution of JT using as input a fraction of the early evolved outer belt at the onset of the giant planets instability. Additionally, we investigate the history of cratering on bodies with sizes equivalent to the targets of the Lucy mission.

This paper is organised as follows. Section 2 introduces the methodology used to model the collisional evolution, giving an overview of our ALICANDEP-JT code. Section 3 compares ALICANDEP-JT results with previous models under the same or equivalent conditions. In Section 4, we present our new results for modelling the collisional evolution of JT under the capture hypothesis. Finally, Section 5 summarises our results and discusses their implications.

2 Methods

In this study, we used ALICANDEP-JT, a simplified version of the code ALICANDEP (Benavidez & Campo Bagatin 2009; Campo Bagatin & Benavidez 2012), adapted to track the evolution of a single zone. The collisional evolution code ALICANDEP-JT is divided into two components: the first component is an algorithm determining the impact outcome of cratering, fragmentation, and re-accumulation, following collisions between two bodies of any given size, within the size range under study. Then, the main code component models the collisional evolution of a given initial population over time.

2.1 Impact outcome

In this section, we outline the essential features of the first component of the model, based on the Petit & Farinella (1993) algorithm to model the cratering and fragmentation processes with the improvements introduced by O’Brien & Greenberg (2005). For further details, we recommend referring to Petit & Farinella (1993).

The main output of the algorithm is the fragmentation matrix f(i,j,k), which is used as input in the following evolution equations. Such a matrix is a three-dimensional (3D) array containing the number of fragments of any size resulting from all possible collisions between two bodies of diameters d(i) and d(j). For any single collision between two bodies of given masses (M_i, M_j), at some relative velocity (V_rel), their relative kinetic energy (E_rel) can be compared with the threshold shattering energy of each colliding body to determine the outcome: cratering or shattering. The relative kinetic energy is: $$E_{\mathrm{r}\mathrm{e}\mathrm{l}}=\frac{M_i{M}_j{V}_{\mathrm{r}\mathrm{e}\mathrm{l}}^2}{2(M_i+M_j)}.$$(1)

The impact kinetic energy is assumed to be equally partitioned between the two colliding bodies. Fragmentation occurs if: $$E_{\mathrm{r}\mathrm{e}\mathrm{l}}>2Q_S^{\ast }{M}_i,$$(2)

where $Q_S^{\ast }$ is the threshold specific energy for shattering, defined as the kinetic impact energy per unit mass of the target such that the impact produces the largest fragment with half the mass of the original target. The occurrence of a cratering event is a direct result of reversing the condition in Eq. (2). Thus, an abrupt transition from cratering to shattering regime is assumed, in agreement with a wealth of literature on shattering experiments, as in Fujiwara et al. (1977).

Generally speaking, when two bodies collide, there are three possible outcomes: (i) both objects are shattered; (ii) one body is shattered, while the other is cratered; or (iii) both objects are cratered. The latter case is limited to low-velocity collisions, which are not expected to take place in the JT system. In each case, the size distribution of fragments from each body needs to be calculated and combined to determine the final size distribution after computing the eventual fragment re-accumulation by self-gravity.

In impacts with a specific energy larger than the critical specific energy threshold for shattering, the mass fraction of the largest fragment with respect to the target mass is given by $$f_{l,i}=M_{max,i}/M_i,$$(3)

where f_l,i is a function of the specific impact energy ($Q_S^{\ast }$), the target mass (M_i), and the relative kinetic energy (E_rel). It can be determined as follows.

Shattering events (E_rel >2$Q_S^{\ast }$ M_i):

• $f_{l,i}=\frac{1}{2}{\left(\frac{Q_{S,i}{M}_i}{E_{\mathrm{r}\mathrm{e}\mathrm{l}}/2}\right)}^{1.24}$,

• The mass distribution of fragments follows a power law with an exponent, b_i, which is related to f_l,i and estimated by a numerical approach (see Eq. (9) in Petit & Farinella 1993);

Cratering events (E_rel < 2$Q_S^{\ast }$ M_i):

• The largest surviving body always has a mass of M_i − M_crat, where M_crat is the mass excavated from the crater and converted into ejected fragments. This is calculated based on the values of E_rel and α, the excavation coefficient of the impact, which depends on the material properties. Following Petit & Farinella (1993), we assume α = 1 × 10⁻⁴ s²/m². Thus M_crat can be calculated as: $$\begin{array}{rl}& M_{\text{crat }}=\alpha E_{\text{rel }}\text{ for }{E}_{\text{rel }}\le \frac{M_1}{100\alpha },\\ & M_{\text{crat }}=\frac{9\alpha }{200Q_S\alpha -1}{E}_{\text{rel }}+\frac{M_i}{10}\frac{1-20Q_S\alpha }{1-200Q_S\alpha }\text{ for }\frac{M_i}{100\alpha }>E_{\text{rel}};\end{array}$$

• The mass distribution of fragments follows a single slope power law with an exponent b_i with a constant value of 0.8 (Melosh 1989);

• The mass fraction of the largest fragment originated from the crater is $$f_{l,i}=1-b_i+b_i{f}_{l,i}^{b_i}{\left(\frac{M_{dust}}{M_i}\right)}^{1-b_i},$$

  where M_dust is the lower cutoff of the fragment mass produced distribution.

2.1.1 Re-accumulation model

A collision is always assumed to produce a remnant body comprising the two largest fragments from the colliding bodies. Apart from the largest remnant of each body, a collision also results in a distribution of smaller fragments which acquire non-negligible velocities. Thus, after any collision, each small fragment may either be re-accumulated back on the largest remnant or escape its gravitational potential.

Thus, as explained in Petit & Farinella (1993), using the energy balance equation allows us to compute the escape velocity, as $$\frac{1}{2}(M_1-M_{max,1}+M_2-M_{max,2})v_{esc}^2+W_{tot}=W_{fr,1}+W_{fr,2}+W_h,$$(4)

where W_tot is the total potential energy of the two bodies just before fragmentation, given by: $$W_{tot}=-\frac{3G{M}_1^{5/3}}{5Q_1}-\frac{3G{M}_2^{5/3}}{5Q_2}-\frac{G{M}_1{M}_2}{Q_1{M}_1^{1/3}+Q_2{M}_2^{1/3}},$$(5)

where Q_i is given by $$Q_i={\left(\frac{4\pi {\rho }_i}{3}\right)}^{-1/3}.$$

W_fr,i is the sum of the self-gravitation potential energy of each fragment of body, i, given by $$\begin{array}{rl}{W}_{fr,i}& =\frac{3G}{5Q_i}{\int }_{m=0}^{m=\mathrm{\infty }}{m}^{5/3}{n}_i(m)dm\\ & =-\frac{3G}{5Q_i}\frac{5M_{max,i}^{5/3}-3b_i{M}_{max,i}^{b_i}{M}_{dust}^{5/3-b_i}}{5-3b_i},\end{array}$$

and W_h is an estimate of the potential energy when the fragments are separated by a distance of the order of the Hill’s radius.

Here, we include the corrections made by O’Brien & Greenberg (2005) to compute W_h according to the impact outcome. Thus, if both bodies are shattered, we have $$W_h=-\frac{3G{(M_1+M_2)}^{5/3}}{5}\frac{{\left(3M_{sun}\right)}^{1/3}}{R_0}.$$(6)

However, if one body is shattered and the other is cratered, we have $$W_h=-\frac{3G{(M_1-M_{\text{crat},1})}^{2/3}(M_2+M_{\text{crat},1})}{2}\frac{{\left(3M_{sun}\right)}^{1/3}}{R_0}.$$(7)

In the case where both bodies are cratered, we have $$\begin{array}{rl}{W}_h=& -3G{(M_1+M_2-M_{\text{crat},1}-M_{\text{crat},2})}^{2/3}\\ & \frac{\times (M_{\text{crat},1}+M_{\text{crat},2})}{2}\frac{{\left(3M_{sun}\right)}^{1/3}}{R_0},\end{array}$$(8)

where M_sun and R₀ are the mass of the Sun and the orbital radius, which is on the order of 5.2 au for JTs.

Then, it is time to compute the number and mass of fragments in bins representing the distribution function. Mass bins are organised into N_bin logarithmic, discretised bins, having their central values scaled by a factor of 2, and their limits are scaled accordingly. The central mass values are given by the following geometric sequence: $$m_{bin,k}=g^{k-1}{m}_{bin,1},$$(9)

where g = 2 and k is the bin number. Each bin spans the range of $$\frac{m_{bin,k}}{\sqrt{g}},m_{bin,k}\sqrt{g}.$$(10)

Therefore, neighbouring mass bins scale by a factor of 2, and bodies in bin k + 3 are twice the size of those in bin k.

Petit & Farinella (1993) defined the total kinetic energy of ejecta from body i as E_fr,i = f_ke E_rel/2, where f_ke is the fraction of the impact energy that is partitioned into kinetic energy of the fragments (see Sect. 2.1.3 for a new approach to this parameter).

Our model assumes a cumulative distribution for the fraction of mass larger than a given value, V: $$f(>V)={\left(\frac{V}{V_{min}}\right)}^{-K},$$(11)

where the exponent takes the value K = 9/4 (Gault et al. 1963) and V_min is a fixed lower cutoff to the ejection velocity. Choosing a V value equal to the escape velocity allows us to find the fraction of escaping mass from the target (Campo Bagatin et al. 2001).

The fraction of re-accumulated mass is determined in the following way. If the number of fragments in a given size bin is less than 100, the velocity of each fragment is drawn at random from a Maxwellian velocity distribution: $$P(v;V_0)=\sqrt{\frac{2}{\pi }}\frac{\sqrt{3}{v}^2}{V_0^3}\exp (-\frac{3v^2}{2V_0^2}),$$(12)

where v is the impact speed and $V_0^2=⟨v^2⟩$ is the rms velocity value. Otherwise, if the number of fragments in the bin is larger than 100, the fraction of non-escaping fragments is calculated as: $$f={\int }_0^{V_{esc}}P(v;V_0)dv,$$(13)

where V_esc is the escape velocity.

2.1.2 Scaling laws for shattering and disruption

Scaling laws are mathematical expressions derived from observations, experiments, and shattering modelling, intended to extend the validity of experimental and modelling outcomes beyond the original studied size range. A scaling law depends on multiple parameters involved in the collision process: impact speed, material properties, density, gravitational forces, and so on. It comprises a set of parameters that is a function of the target diameter.

Over the years, several scaling laws were proposed, either in terms of threshold impact strength (S₀) for shattering (i.e. Davis et al. 1985; Housen & Holsapple 1990) or in terms of disruption specific energy ($Q_D^{\ast }$). The former was mainly derived from experiments, the latter was derived from numerical simulations. Here, $Q_D^{\ast }$ represents the threshold specific energy required to disperse the target into a spectrum of individual fragments but possibly re-accumulated, the largest fragment having exactly half of the mass of the original target. Benz & Asphaug (1999), through smoothed-particle hydrodynamics (SPH), self-consistently modelled this $Q_D^{\ast }$ for a wide range of diameters. The result takes the functional form of $$Q_D^{\ast }=A{\left(\frac{R_T}{1\mathrm{c}\mathrm{m}}\right)}^a+B\rho {\left(\frac{R_T}{1\mathrm{c}\mathrm{m}}\right)}^b,$$(14)

where A, a, B, b are constants determined for different target compositions and impact speeds, while R_T is the target radius. The first term in Eq. (14) corresponds to the strength regime, where fragments cannot re-accumulate (the fragments are too dispersed and moving too fast to be held together by gravity). On the other hand, the second term corresponds to the gravity regime. This regime comes into play when the mass of the shattered target is large enough so that most ejected fragments cannot escape the gravitational potential well. In the gravity regime, fragments can re-accumulate and form gravitational aggregates (rubble piles).

O’Brien & Greenberg (2005) and de Elía & Brunini (2007) carried out collisional evolution simulations using the specific energy required for shattering a target into a spectrum of fragments, the largest one having exactly half the mass of the original target, denoted as $Q_S^{\ast }$. Alternatively, a number of studies (Davis et al. 1985; Petit & Farinella 1993; Campo Bagatin et al. 1994, 2001; Benavidez & Campo Bagatin 2009, among others), have used a comparable approach focusing on impact strength, namely, $Q_S^{\ast }$ = S₀/ρ.

The units of $Q_S^{\ast }$ and $Q_D^{\ast }$ are energy per unit mass, which allows us to make a rough comparison of the two scaling laws, with caveats. In the strength regime, both coincide, as gravity does not contribute. The meaning of the rightmost branch (e.g. Fig. 1) of the two functions is different: In the case of $Q_S^{\ast }$, the increase with size corresponds to the effect of gravitational self-compression. In the case of $Q_D^{\ast }$, the energy required to disperse fragments against the gravitational potential adds to self-compression, producing a noticeable difference between $Q_S^{\ast }$ and $Q_D^{\ast }$.

In an attempt to understand the link between $Q_S^{\ast }$ and $Q_D^{\ast }$, de Elía & Brunini (2007) parameterised Eq. (14) as $$Q_D^{\ast }=C_1{D}^{-{\lambda }_1}(1+{\left(C_{2}D\right)}^{{\lambda }_2}).$$(15)

To reproduce the Benz & Asphaug (1999) scaling law for icy bodies with impact speeds of 3 km/s, the constant parameters C₁, C₂, λ₁, and λ₂ were set to values 24, 2.3, 0.39, and 1.65, respectively (to provide $Q_D^{\ast }$ in SI units). The idea was to have the possibility to examine a family of $Q_S^{\ast }$ identical to $Q_D^{\ast }$ in the strength regime, but with a lower value in the gravity regime. By varying the values of C₂ and λ₂ in Eq. (15), it is possible to create a family of $Q_S^{\ast }$ that could be related to $Q_D^{\ast }$. Therefore, $Q_S^{\ast }$ is expressed as: $$Q_S^{\ast }=C_1{D}^{-{\lambda }_1}[1+{\left(C_2^{\mathrm{\prime }}D\right)}^{{\lambda }_2^{\mathrm{\prime }}}],$$(16)

where $C_2^{\mathrm{\prime }}$ and ${\lambda }_2^{\mathrm{\prime }}$ are slight variations of the parameters C₂ and λ₂. In this work, we updated the Petit & Farinella (1993) algorithm, including Eq. (16) for the scaling law and tested values for $C_2^{\mathrm{\prime }}$ and ${\lambda }_2^{\mathrm{\prime }}$. The goal of this approach is to crumble the $Q_D^{\ast }$ into a combination of factors as $Q_S^{\ast }$, f_ke (fraction of kinetic energy transferred to fragments) and reaccumulation processes. We give more details on this approach in Sect. 3.3.

Fig. 1 Comparison of two selected scaling laws from EB07, $Q_{S,a}^{\ast }$ = (1.2, 1.75) and $Q_{S,b}^{\ast }$ = (2.3, 1.25), constructed using Eq. (16). For reference, the standard scaling law for $Q_D^{\ast }$ from Benz & Asphaug (1999) is also shown.

2.1.3 Fraction of kinetic energy transferred to fragments

An important parameter in the re-accumulation process is f_ke, defined as the fraction of impact energy that is transferred as kinetic energy of the produced fragments. This is directly related to its re-accumulation chances, after comparison with the gravitational potential energy of the rest of the mass.

Nonetheless, our knowledge of this parameter is very limited, as the available data are restricted to laboratory experiments on centimeter to decimeter targets sizes (Nakamura & Fujiwara 1991; Giblin et al. 2004), and to hydrocode experiments on targets with a diameter over a few hundreds of metres, such as that of Melosh & Ryan (1997). O’Brien & Greenberg (2005) suggested a power-law function for f_ke, depending on size (diameter): $$f_{ke}=f_{k{e}_0}{\left(\frac{D}{1000\mathrm{k}\mathrm{m}}\right)}^{\gamma },$$(17)

where γ is always between 0 and 1 and f_ke0 is the value at 1000 km. Suitable values of f_ke0 are in the range of 0.05 to 0.3 (Davis et al. 1989). de Elía & Brunini (2007) attempted to model the collisional evolution of the L₄ swarm, using γ = 0.7 and f_ke0 = 0.35. These values generate an f_ke at smaller sizes that is not aligned with the experimental data (Nakamura & Fujiwara 1991; Giblin et al. 2004) (see Fig. 5). Therefore, in this work, we experimented with variations of Eq. (17), while comparing our model with other collisional models (see Sect. 3.3). Such a formulation allows for a more flexible use of both parameters, $Q_S^{\ast }$ and f_ke.

2.2 Collisional evolution

The collisional evolution of a population of interacting bodies, including their mutual fragmentation, is performed by a set of n non-linear, second-order differential equations. Such equations can be numerically solved. The change per time step in the number of objects belonging to a given bin k is calculated as: $$\mathrm{\Delta }N(k)=\sum _i^n\sum _j^i{P}_{int}\cdot f(i,j,k)\cdot s(i,j)\cdot f_g\cdot \frac{(N(i)-{\delta }_{ij})N(j)}{1+{\delta }_{ij}}\mathrm{\Delta }t,$$(18)

where P_int is the intrinsic probability of collision inside the JT swarm, s(i, j) = 1/4 · (D_i + D_j)² is the impact cross-section and f_g ≈ 1 + (V_esc/V_rel)² accounts for the gravitational focusing increasing the cross-section. Therefore, the number of objects in bin k at time step t + δt is $$N_k(t+\delta t)=N_k(t)+\delta t\cdot d{N}_k(t).$$(19)

The bin population is integrated in time through the Euler method, satisfying the following relation at each time step, $$d{N}_k(t+\delta t)<ϵ\cdot N_k(t),$$(20)

where 0 < ϵ < 1 is the parameter that specifies the fraction of the population in any bin k that should not be increased/decreased at each time step. Then, the code readjusts the time step (increasing or decreasing it) for the next iteration according to the actual difference in the number of objects. This feature avoids divergence in the evolution algorithm.

Campo Bagatin et al. (1994) noted a fictitious wavy pattern that can potentially extend across all size ranges in collisional evolution system. In their study, the authors elucidated that this structure arises due to a sharp cutoff or discontinuity at some small size. To neutralise the effect of such a wavy pattern, an artificial tail at small sizes was introduced in the distribution. The tail extends from the smaller limit assumed for fragments (M_dust) up to the first size bin selected as the lower limit for evolution. The role of the tail is only to provide projectiles to the smallest size range of bodies. Thus, the evolution of the tail itself is not performed. Instead, it is extrapolated through a power law, with a slope calculated from a set of the smallest size bins (e.g. 15, 20) in the evolving populations. This methodology is efficient and removes any artificial wavy structure.

In addition, the model includes the possibility to have either a constant relative impact velocity or a Maxwellian velocity distribution function. This feature can be toggled on or off according to specific requirements.

We considered only mutual collisions among L4 trojans, as previous studies have shown that the contribution of other projectile populations, such as Hilda asteroids and Jupiter-family comets (JFCs), is negligible. The intrinsic collision probabilities between these populations and the L4 trojans are more than an order of magnitude lower than those for trojan-trojan collisions – for instance, a factor of 24 and 26 lower for JFCs and Hildas, respectively (Dell’Oro & Paolicchi 1998; Dell’Oro et al. 2001). In addition, the number of potential impactors is significantly smaller: recent estimates suggest that the number of Hildas larger than 3 km (Zain et al. 2025; Vokrouhlický et al. 2025) is at least one order of magnitude smaller than the corresponding L4 trojan population (Wong & Brown 2015). As a result, their contribution to the collisional evolution is negligible compared to trojan-trojan collisions. While cometary impacts may have played a more significant role during a brief early phase after implantation (~10 Myr, Marchi et al. 2023), our model aims to characterise the long-term collisional evolution, during which any potential evidence of that transient phase may have been smoothed out.

Although our model does not explicitly track the dynamical evolution of fragments after collisions, we note that some ejected fragments may escape the L4 region if their post-impact orbits fall outside the stability zone. This loss mechanism has been investigated in previous studies. Marzari et al. (1995) initially suggested that up to ~20% of collisional fragments might leave the trojan clouds but later refined this estimate to ~10% based on improved modelling (Marzari et al. 1997). More recently, de Elía & Brunini (2007) applied a detailed orbital analysis and concluded that the escape rate is significantly lower, with about 50 bodies larger than 1 km leaving the L4 region per Myr. These results suggest that, although this mechanism exists, it likely represents a minor contribution to the overall mass loss and does not significantly impact the long-term collisional evolution. A broader discussion of this process can be found in the review by Emery et al. (2015).

3 Testing ALICANDEP-JT model

In this section, we describe how we used ALICANDEP-JT to reproduce the modelling results of two previous studies on the collisional evolution of JTs, namely, de Elía & Brunini (2007) (EB07, hereafter) and Marschall et al. (2022) (M22, hereafter). By comparing our results with EB07, we aim to validate our numerical model, while the comparison with M22 allows us to assess different approaches to the implemented scaling laws.

3.1 Comparison with previous JT collisional evolution models

First, we compared ALICANDEP-JT with de Elía & Brunini (2007). Since EB07 adopts a similar approach to the fragmentation and re-accumulation process to ALICANDEP-JT (as described in Sect. 2.1), this comparison serves to validate our code. Once this test was done, the next step was to understand how the combination of the $Q_S^{\ast }$ scaling law and the subsequent re-accumulation (influenced by f_ke) can produce results comparable to those of a model using $Q_D^{\ast }$, such as M22. In all simulation runs discussed in Sect. 3, we used the same boundary conditions and parameters outlined in the referenced studies to perform analogous numerical simulations with our model, granting a consistent comparison.

We carried out Kolmogorov–Smirnov (KS) tests to compare our results with previous studies. Such test compares two populations and determines whether two samples come from the same underlying population or have significantly different distributions. The test takes the two samples as input and returns two outcome values: namely, the test statistics (D) and the p-value. The D value is the maximum difference between the cumulative distribution functions of the two samples. It represents a measure of the discrepancy between the two distributions. The p-value represents the probability of observing the obtained D, or an even more extreme value, given that the null hypothesis is true. The null hypothesis is the assumption that the two samples are drawn from the same underlying population or have the same distribution. A minimum p-value significance level of 0.05 is commonly considered sufficient to retain the null hypothesis. However, a visual analysis is also useful because some patterns and details are invisible to a statistical test.

3.2 Comparison with de Elía & Brunini (2007)

Here, the initial population is defined by a cumulative SFD described by a broken power law with a transition diameter at 60 km. Thus, the starting cumulative SFD has a shallow slope, q = 3, for small bodies, and a steeper slope, q = 4.7, for larger bodies. The result is an initial population with a mass of 8.37 × 10⁻⁵ M_⊕ (assuming a bulk density of 1500 kg m⁻³ for bodies), which is significantly more massive than the value predicted nowadays for the JT (~10⁻⁵M_⊕) (Morbidelli et al. 2005; Vinogradova & Chernetenko 2015; Pitjeva & Pitjev 2020). To model the collisional evolution of the L4 swarm, EB07 considered a constant value for relative impact speed, 4.66 km s⁻¹ and an intrinsic impact probability value: 7.79 × 10⁻¹⁸ km⁻² yr⁻¹ (Dell’Oro & Paolicchi 1998).

EB07 developed an algorithm based on Petit & Farinella (1993) to model the impact outcome, and used as reference the Benz & Asphaug (1999) scaling law corresponding to icy targets at 3 km/s impact velocity. To compare our results to those of EB07, we perform runs using several of their sixteen different scaling laws. Figure 2 illustrates two specific cases where the selected scaling laws, $Q_{S,a}^{\ast }$ and $Q_{S,b}^{\ast }$, were derived from Eq. (16) with parameter values ($C_2^{\mathrm{\prime }}$, ${\lambda }_2^{\mathrm{\prime }}$) set to (1.2, 1.75) and (2.3, 1.25), respectively. These $Q_S^{\ast }$ scaling laws are plotted alongside the standard Benz & Asphaug (1999) scaling law for $Q_D^{\ast }$ (icy composition). It is important to note that $Q_D^{\ast }$ and the $Q_S^{\ast }$ family are not directly comparable as they describe different definitions (threshold specific energy for disruption vs. threshold specific energy for shattering).

Table 1 presents the D-statistic and the p-value for each case. Since the p-values for both cases are above a significance level of 0.05, we accept the null hypothesis, indicating that the distributions obtained using EB07 and ALICANDEP-JT are statistically indistinguishable. Furthermore, Fig. 2 shows a visual comparison of the cumulative SFD of two tested cases. The SFDs of both numerical models show similar features, with dips and peaks occurring at the same size ranges.

The selected sampled cases from EB07 provide highly consistent results. Both scaling laws $Q_{S,a}^{\ast }$ and $Q_{S,b}^{\ast }$ demonstrate excellent agreement with EB07, reinforcing the validation of our model.

3.3 Comparison with Marschall et al. (2022)

The initial distribution used by M22 is characterised by a cumulative SFD given by a power law with a slope of q = 2.1 for small bodies, and a steeper slope (q ~ 4.5) for bodies larger than 100 km. That produces a less massive initial population compared to EB07 (4 × 10⁻⁶ M_⊕, assuming a bulk density of 1000 kg/m³). Note that, at this stage, we are assuming identical initial conditions only for comparison purposes. Thus, the intrinsic impact probability is set to 7 × 10⁻¹⁸ km⁻² yr⁻¹ and the impact velocity is set to a unique value of 4.6 km s⁻¹ (Davis et al. 2002; Nesvorný et al. 2018).

In Marschall et al. (2022) the impact outcome is based on SPH modelling combined with Parallel K-D tree GRAVity code (pkdgrav) following the methodology described in Bottke et al. (2005). The authors examined two distinct characteristic scaling laws based on Eq. (14). The first was inspired by the scaling law for basalt material proposed by Benz & Asphaug (1999) with a minimum at ~250 m, which they also scale by different factors. The second, referred to as the ‘unconventional’ scaling law, has a minimum at 20 m. This unconventional choice for the minimum of Q_D was originally proposed to better match the crater distribution observed on irregular and regular satellites of giant planets (see Bottke et al. 2023b; Bottke et al. 2024).

In this section, we consider two scaling laws for comparison with Marschall et al. (2022). In case 1, the applied scaling law is indicated by Q₀, which is approximately 100 times weaker than the scaling law for basalt proposed by Benz & Asphaug (1999), with a minimum at 250 m. In case 2, the applied scaling law corresponds to the ‘unconventional’ one in M22, with a minimum at 20 m.

M22 employed a $Q_D^{\ast }$ scaling law (Eq. (14)), whereas our code utilises a $Q_S^{\ast }$ scaling law (Eq. (16)); however, as it was stated in Sect. 2.1.2, $Q_D^{\ast }$ sets an upper limit for $Q_S^{\ast }$. Here we test different values of $C_2^{\mathrm{\prime }}$ and ${\lambda }_2^{\mathrm{\prime }}$ in Eq. (16), combined with the approach for f_ke described in Sect. 2.1.3, until we find a good match for the final SFD. A verification that $Q_S^{\ast }$ is consistent with $Q_D^{\ast }$ is always necessary (i.e. $Q_S^{\ast }$ must always be smaller than $Q_D^{\ast }$ in the gravity regime). Figures 3 and 4 show the set of tested scaling laws. To check the influence of different parameters, one of the parameters is kept constant while the other is given different values. Figure 3a the shows scaling law, Q₀, where $C_2^{\mathrm{\prime }}$ is constant and ${\lambda }_2^{\mathrm{\prime }}$ is given three different values, while Fig. 3b shows the outcome corresponding to a constant ${\lambda }_2^{\mathrm{\prime }}$ and $C_2^{\mathrm{\prime }}$ taking three different values. Figure 4 is the same as Fig. 3, but for the ‘unconventional’ scaling law. The influence of the explored scaling laws is explained in more detail later in this section.

In Sect. 2.1.3, we mention that f_ke can be given as a function of diameter (Eq. (17)). Here, we further explore this concept by adding a constant term to the expression: $$f_{ke}=f_{k{e}_0}{\left(\frac{D}{1000\mathrm{k}\mathrm{m}}\right)}^{\gamma }+\beta ,$$(21)

where β is a positive constant added to match the results of the impact experiment (Nakamura & Fujiwara 1991; Giblin et al. 2004). Such an expression allows for a small but not negligible value of f_ke at small sizes, as opposed to the former formulation of that parameter, which vanished quickly within the same size range. In addition, γ can be adjusted to set the diameter at which the function separates the constant regime and the exponential regime.

Four combinations of parameters are explored to assess the influence of f_ke. Such combinations are listed in Table 2 and are represented in Figure 5. The first case corresponds to the parameters used by de Elía & Brunini (2007) (Eq. (17)). The second case uses the same equation with values of f_ke0 and γ in accordance with Melosh & Ryan (1997). Through hydrocode calculations, Melosh & Ryan (1997) studied the impact energy in the context of fragment dispersion of large bodies. They also present constraints on f_ke, namely a value of 0.1 at 1000 km and a D^1/2 decrease for small bodies (O’Brien & Greenberg 2005). In the third case, the β term is implemented to fit the experimental data at small sizes (Eq. (21)) (Nakamura & Fujiwara 1991; Giblin et al. 2004), but still using the results from Melosh & Ryan (1997), i.e., the f_ke value at 1000 km and the γ value are the same. (At 1000 km: f_ke = f_ke0 +β = 0.1). Nakamura & Fujiwara (1991) carried out an experimental study on the velocity distribution of fragments after a collision with basalt and alumina spherical targets of 6 cm in diameter, and they found f_ke ~ 0.01. Giblin et al. (2004), in a similar study on icy, 15 cm diameter targets, found a wider range of f_ke values, still in the same order of magnitude as in Nakamura & Fujiwara (1991). Giblin et al. (2004) manufactured three different types of icy targets: two types of ‘rubble piles’ and one more compact snowball; combinations of these three types of target (each type including targets with slightly different masses in each collision event) combined with three types of projectile generate various initial conditions and, as a result, different f_ke. In the fourth case, an increase in the value of γ is explored to better fit the data from Marschall et al. (2022). Appendix A shows the results for all different combinations of parameters in Table 2, highlighting the need for the non-negligible β parameter used in this section.

When comparing our results with those of Marschall et al. (2022), it is important to acknowledge that multiple sets of runs for the same case may show changes due to random seeds employed in the Boulder code. Consequently, the single SFD provided by R. Marschall and used for comparison in our study represents a potential outcome of the same initial system state. Thus, Figures 10 and 13 in M22 show slight differences (up to a factor of 2) in the number of bodies at sizes smaller than a few kilometres. This dispersion is represented as a shaded region in Figures 6 and 7, emphasising the essential agreement of our results.

We conducted a Kolmogorov-Smirnov test to provide a robust comparison of the results obtained by both models. We notice that, for comparison, we used an SFD falling roughly in the middle of the evolved distributions by M22. The best agreement with the reference simulations is achieved with the parameter values $C_2^{\mathrm{\prime }}$ = 0.3 and ${\lambda }_2^{\mathrm{\prime }}$ = 1.6 for case 1 (Q₀), and $C_2^{\mathrm{\prime }}$ = 8 and λ′2 = 1.8 for case 2 (the unconventional scaling law). In both cases, with a suitable choice of the parameters defining f_ke and $Q_S^{\ast }$, the evolved SFDs closely match those in M22. Specifically, the corresponding $Q_S^{\ast }$ exhibits a minimum shifted toward larger sizes by a factor of 2 and a strength reduced by a factor of 2 relative to the reference $Q_D^{\ast }$ scaling law (see Figures 3b and 4b). The following two sections give a detailed comparison with M22 for case 1 and case 2, guiding in the identification of the combination of scaling law parameters that best reproduces the reference runs from M22 summarised above.

Fig. 2 Cumulative SFDs for cases $Q_{S,a}^{\ast }$ and $Q_{S,b}^{\ast }$. Dashed lines correspond to the initial population, while solid lines indicate evolved populations with ALICANDEP-JT code, and dotted lines correspond to evolved populations in de Elía & Brunini (2007).

Fig. 3 $Q_S^{\ast }$ scaling laws tested to study case 1. (a) parameter $C_2^{\mathrm{\prime }}$ is kept constant at value 0.85, while ${\lambda }_2^{\mathrm{\prime }}$ takes three different values. (b) parameter ${\lambda }_2^{\mathrm{\prime }}$ is kept constant at 1.60, while $C_2^{\mathrm{\prime }}$ takes three different values. The solid black line corresponds to the Q0 case, namely the $Q_D^{\ast }$ used by Marschall et al. (2022), for a reference.

Fig. 4 $Q_S^{\ast }$ scaling laws tested to study case 2. (a) parameter $C_2^{\mathrm{\prime }}$ is kept constant at value 10, while ${\lambda }_2^{\mathrm{\prime }}$ takes three different values. (b) parameter ${\lambda }_2^{\mathrm{\prime }}$ is kept constant at 1.80, while $C_2^{\mathrm{\prime }}$ takes three different values. The solid black line corresponds to the unconventional scaling law, namely, that of $Q_D^{\ast }$, used by Marschall et al. (2022) for a reference.

Fig. 5 Functions of fke tested and listed in Table 2. Black symbols correspond to values from Nakamura & Fujiwara (1991); Melosh & Ryan (1997); and Giblin et al. (2004), as indicated in the legend.

Fig. 6 Cumulative SFDs for runs corresponding to case 1, using the scaling laws showed in Figs. 3a and 3b). The fke function is determined by fke0 = 0.085, γ = 0.5 and β = 0.015. The solid black line is the evolved cumulative SFD using the scaling law, Q0, from Marschall et al. (2022) for visual comparison.

3.3.1 Case 1: Q₀ scaling law (D_min = 250 m)

The scaling law, Q₀, can be expressed by replacing parameters A, a, B and b in Eq. (14) by 1.50 × 10⁶ erg g⁻¹, −0.40, 1.30 erg cm³ g⁻², and 1.30, respectively (Marschall et al. 2022). The same scaling law can be expressed as Eq. (15) using 1.9, 0.4, 2 and 1.7 in place of C₁, λ₁, C₂ and λ₂, respectively. For this scaling law, parameter $C_2^{\mathrm{\prime }}$ is tested for values ranging from 0.3 to 1.2, combined with ${\lambda }_2^{\mathrm{\prime }}$, ranging from 1.4 to 1.85, in Eq. (16). The constant values for the f_ke function are f_ke0 = 0.085, γ = 0.5 and β = 0.015, which correspond to the best combination for this scaling law that match the reference run by M22 (Appendix A shows the results for all different parameter combinations in Table 2).

Table 3 shows the combinations of $C_2^{\mathrm{\prime }}$ and ${\lambda }_2^{\mathrm{\prime }}$ explored in this study, along with the D and p-values for all runs. The statistical test outcomes indicate a substantial agreement between the evolved populations using ALICANDEP-JT and M22. The values of D are low, and the p-values consistently exceed the significance level, often approaching 1. As a result, comparing p-values may not help in determining the best run. Then, a visual inspection of the SFD is convenient at this stage.

Figure 6 shows the cumulative SFDs evolved under the $Q_S^{\ast }$ scaling laws presented in Fig. 3. All modelled SFD obtained for different combinations of $C_2^{\mathrm{\prime }}$ and ${\lambda }_2^{\mathrm{\prime }}$ can reproduce the dip-like feature as the reference SFD modelled by M22 (solid black line). Smooth differences in the evolved population arise from variations in these parameters. In the cases where $C_2^{\mathrm{\prime }}$ are kept constant, there are only negligible variations in the evolved cumulative SFD (see Fig. 6a). In the cases where ${\lambda }_2^{\mathrm{\prime }}$ is kept constant, more pronounced variations can be observed. As $C_2^{\mathrm{\prime }}$ decreases, the dip feature becomes deeper (as shown in Fig. 6b). In summary, the scaling law parameters $C_2^{\mathrm{\prime }}$ = 0.3 and ${\lambda }_2^{\mathrm{\prime }}$ = 1.60 yield an SFD that closely aligns with the reference SFD in M22.

Fig. 7 Cumulative SFD for runs corresponding to case 2, using the scaling law showed in Figs. 4a and 4b. The fke function is determined by fke0 = 0.090, γ = 0.5, and β = 0.010. The solid black line is the evolved cumulative SFD from Marschall et al. (2022) for a visual comparison.

3.3.2 Case 2: Unconventional scaling law (D_min = 20 m)

The unconventional $Q_D^{\ast }$ can be produced by replacing parameters A, a, B and b in Eq. (14) by 7.65 × 10⁵ erg g⁻¹, −0.36, 1.40 erg cm³ g⁻² and 1.36, respectively (Marschall et al. 2022). The same scaling law can be written as Eq. (16), using 1.55, 0.36, 24 and 1.71 in place of C₁, λ₁, C₂ and λ₂, respectively. In this case, parameter $C_2^{\mathrm{\prime }}$ is tested for values ranging from 4 to 12 combined with ${\lambda }_2^{\mathrm{\prime }}$, ranging from 1.40 to 1.80, in Eq. (16). The constant values for the f_ke function are f_ke0 = 0.090, γ = 0.5 and β = 0.010, which correspond to the best combination for this scaling law that matches the reference run from M22 (Appendix B shows the results for all different parameter combinations of Table 2).

The statistical analysis presented in Table 4 includes D and p-values for all explored cases. The results reveal a significant similarity between the cumulative distributions obtained with ALICANDEP-JT and M22 considering their unconventional scaling law. Figure 7 shows the evolved SFDs under the $Q_S^{\ast }$ scaling laws shown in Figure 4.

The analysis of Figure 7 provides insights into the impact of the scaling law $C_2^{\mathrm{\prime }}$ and ${\lambda }_2^{\mathrm{\prime }}$ parameters on the resulting cumulative SFDs. Specifically, Figure 7b shows that increasing $C_2^{\mathrm{\prime }}$ (while keeping ${\lambda }_2^{\mathrm{\prime }}$ constant) enhances the production of fragments in the 0.1–1 km range, aligning our modelled SFD more closely with the reference distribution in M22. However, this adjustment also leads to a slight overabundance of bodies relative to the initial distribution, particularly near 1 km. This feature, a direct consequence of the unconventional scaling law applied, was also observed by M22 and attributed to the limited yet non-negligible fragmentation of bodies larger than 10 km. Notably, the extreme case tested with $C_2^{\mathrm{\prime }}$ = 4 exhibits significant deviations from the reference run provided by R. Marschall. Overall, our model is consistent with the results of M22, with a preference for $C_2^{\mathrm{\prime }}$ = 8 and ${\lambda }_2^{\mathrm{\prime }}$ = 1.80.

4 Results

As outlined in the introduction, the simulations of this section build upon recent work (B25), which successfully modelled the debiased SFD of the hot and cold populations within the Kuiper Belt. That study encompassed several key components, described below.

The SFD of the primordial outer belt is given by: $$\begin{array}{rl}N(<H_r)=& N_{0_1}{10}^{{\alpha }_s(H_r-H_o)}\exp [-{10}^{\beta (H_r-H_B)}]\\ & +N_{0_2}{10}^{{\alpha }_l(H_r-H_o)},\end{array}$$(22)

where N₀₁ and N₀₂ are constant values. For the parameter values, they followed Petit et al. (2023) using H₀ = −2.6, β = 0.25, and H_B = 8.1. They assumed a representative albedo of the hot population within the Kuiper belt of p_r = 0.08 and a shallow power law slope (α_l = 0.18) for large bodies. Furthermore, the slope of small bodies was examined under two distinct values; in case A, a steep slope was considered, α_S = 0.4, whereas in case B, a shallower slope was chosen, α_S = 0.3. They examined three scaling laws: the nominal one determined by the optimal parameters identified for the unconventional scaling law in Section 3.3.2, and two variations of this scaled by factors of 0.1 and 10 (see Fig. 8). A Maxwellian distribution for impact velocities and different timings for the onset of the giant planets instability (t_i) were used.

In the present work, we modelled the collisional evolution of JTs using ALICANDEP-JT with boundary conditions corresponding to each run in B25. Given the negligible exchange of bodies between the two swarms of JTs and the absence of significant collisional interactions, the modelling process was focussed on the most densely populated group, namely, the L4 swarm. The key inputs for the corresponding numerical simulations are as follows.

• Initial SFD: the initial size-frequency distribution was set as a fraction of the evolved population of the outer trans-Neptunian disc. Here, we relied on the assumption that since the capture process is size-independent, the captured JT population would be expected to maintain the same SFD shape as the evolved outer belt population at the onset of instability. Consequently, we normalised the whole SFD based on the number of bodies in the size bin corresponding to ~200 km, approximately the largest size for the JT population. This operation leads to a capture probability of 6.5 × 10⁻⁷, which falls within the uncertainties of the jump capture probability of (5 ± 3) × 10⁻⁷ estimated by Nesvorný et al. (2013). Thus, the initial SFD preserves all the features of the evolved outer belt at the time of dynamical instability, since the probability of capture at any size is the same. Through this procedure, the resulting total mass of the JT population is about 3 × 10⁻⁶ M_⊕, slightly different depending on the specific initial SFD conditions. We note that we performed runs starting with an initial population refereed as case A (α_S = 0.4) and case B (α_S = 0.3). It is worth noting that the timing of the onset of the instability in the primordial outer belt is of little relevance to the evolution of JTs. This is a consequence of the relatively low-intensity collisional evolution in the primordial disc (B25). For the modelling in this study, we considered the onset of the giant planet instability to take place at t_i = 20 Myr.

• Scaling law: the same three scaling laws used for the primordial outer belt were also used here (see Fig. 8).

• Relative speed and impact probability: The most probable relative speed in the Maxwellian distribution for impact velocities is 4.6 km/s, and the intrinsic impact probability of 7 × 10⁻¹⁸ km⁻² yr⁻¹ are adopted, in agreement with previous estimations (Dell’Oro & Paolicchi 1998; Marschall et al. 2022).

Fig. 8 Scaling laws tested for the collisional evolution of JTs. As reference, the solid black line represents the unconventional scaling law $Q_D^{\ast }$ from Marschall et al. (2022). The solid blue line (SL1) is the $Q_S^{\ast }$ that best reproduces the evolved SFD by Marschall et al. (2022) (see Sect. 3.3.2). The dashed and dash-dotted blue lines represent SL1 scaled up and down by a factor of 10 and 0.1, respectively.

4.1 Final size-frequency distribution

Figure 9 shows the initial and final cumulative SFDs for the captured population evolved under cases A and B. Each panel shows the results for each tested scaling law. Solid lines show the final cumulative distribution after ~4500 Myr of evolution. The black line shows the observed cumulative SFD for the L4 swarm constructed using the data from the Minor Planet Center¹ and the grey dashed line is the observed cumulative SFD for the L4 swarm constructed from Wong & Brown (2015). The size distributions were derived from absolute magnitude (H) data using the standard relationship (Pravec & Harris 2007) as $$D=\frac{1329\times {10}^{-H/5}}{\sqrt{p_V}},$$(23)

where D is the diameter in kilometres and p_V is the geometric albedo in the V band. Moving from absolute magnitudes to diameters may introduce significant uncertainties, as real objects have a range of albedo values (Grav et al. 2012). This introduces systematic uncertainties into size estimates based solely on photometric data. For this calculation, we adopted an average geometric albedo of p_V = 0.07 for JTs larger than approximately 10 km, as estimated by Grav et al. (2011).

A noticeable gap emerges between the SFDs evolved from the two different α_S slope values (referred to as cases A and B). Using the same scaling law, case A (α_S = 0.4) keeps its steep SFD that persists over the observed distribution. On the contrary, case B, corresponding to a shallower slope (α_S = 0.3), underestimates the population, particularly at sizes smaller than ~20 km. This finding is consistent with the conclusions by B25, who determined that a slope of α_S = 0.4 in the outer primordial distribution achieves an excellent match with the de-biased populations of the Kuiper belt.

A common feature of cases A and B is an excess of bodies around 1 km that is not usually observed when traditional scaling laws (which have a minimum at hundreds of kilometres) are applied. This feature is a by-product of the unconventional scaling law used and was also noticed by M22. Here, we see that this feature is more pronounced for runs corresponding to the scaling law for strong bodies, and it is smoothed for the weakest bodies. For instance, Fig. 9a shows an excess of bodies – up to a factor of 2 over the corresponding initial population – in a range spanning up to one order of magnitude around 1 km. This is due to inefficient collisional grinding (in the relatively rarefied JT environment) of newly produced fragments from the disruption and cratering of bodies larger than 10 km.

Focusing on the results of case A, which are more consistent with the observational data, and comparing the final SFD for different scaling laws (represented by the solid blue lines in Figure 9), it is evident that each scaling law produces a distinct shape in the final SFD. Thus, the weaker the scaling law, the smaller the gap between the modelled and the observed JT SFDs. Conversely, decreasing the strength (e.g. to the nominal scaling law, shown in the middle panel) heightens the fragmentation of bodies around 10 km, giving rise to a dip in the SFD in the 5 to ~30 km range. Hence, the evolved SFD approximates the observed distribution in that size range. Notably, the SFD modelled using the scaling law that represents the weakest bodies (bottom panel) aligns even more closely with the observed SFD. This is due to the fragility of bodies, which allow for fragmentation up to ~100 km, thereby increasing the depth of the dip, which in this case extends up to ~60 km. Consequently, the scaling law SL1 × 0.1, combined with a suitable initial SFD (α_S = 0.4), matches the observed cumulative size distribution of JT above the completeness limit (~10 km).

Fig. 9 Initial (dashed lines) and final (solid lines) cumulative SFDs of the modelled JT population from each scaling law shown in Fig. 8. Blue and red lines correspond to cases A and B, corresponding to αS = 0.4 and αS = 0.3, respectively. The black solid line represents the cumulative SFD of data from the Minor Planet Center for the L4 swarm, and the dashed gray line corresponds to the SFD estimated by Wong & Brown (2015), including an extrapolation to smaller sizes.

4.2 Collisional families

An additional way to constrain our modelling results is by comparing the number of collisional families whose parent bodies were larger than ~100 km. This provides a valuable observational constraint, as only a few such large families are currently identified among the JTs. Based on a detailed analysis of the trojan population, Vokrouhlický et al. (2024) identified three to five families whose parent bodies likely exceeded 100 km in diameter, using the methodology of Durda et al. (2007) to estimate progenitor sizes. This range offers a useful benchmark to test our collisional evolution model. In our best-case scenario (SL1×0.1, case A), the initial SFD has 17 bodies with diameters larger than D ~ 90 km. After ~4.5 Gyr of evolution, five of them underwent catastrophic disruption: three with D ~ 95 km and two with D ~ 120 km. While the diameters of our disrupted bodies are somewhat smaller than those reported by Vokrouhlický et al. (2024), the agreement is overall satisfactory, especially considering the different approaches used to define parent body sizes. Our results are therefore consistent with the observed number of large trojan families and support the validity of the size distribution and collisional parameters adopted in our model. This supports the ‘weak’ scaling law as the most suitable one for modelling the evolution of the JTs.

4.3 Comparison with the colour sub-populations of JTs

The bimodal colour distribution observed among JTs has been interpreted as evidence of two distinct subpopulations: a red (R) group, considered to be the remnant of the captured primordial bodies, and a less red (LR) group, proposed to consist mainly of fragments produced through collisional evolution after emplacement in the region. This hypothesis, widely discussed in Wong & Brown (2015); Wong & Brown (2017), is supported by the different SFDs observed for the two colour groups. In particular, the LR population exhibits a significantly steeper slope and becomes dominant at smaller sizes, whereas the R population follows a broken power-law SFD consistent with a preserved primordial distribution at large sizes. In this framework, collisions play a crucial role in shaping the present-day colour and size distribution of the trojan population, suggesting that the surface properties of small objects are largely determined by collisional processing. To test the consistency of this scenario with our collisional evolution models, we made use of a tool in our code that allows us to track the number of primordial bodies throughout the evolution. We applied this tool to our reference case (SL1×0.1, case A) and examined the SFD of the total evolved population, the remaining primordial objects, and the collisional fragments produced by the evolution. The results are presented in Fig. 10. The initial SFD, which resembles that of the trans-Neptunian source region at the time of capture, is shown as a black dashed line. The solid black line shows the total population at the end of the simulation. The remnant primordial bodies, corresponding to the R population, are plotted by a dashed red line, while the fragments (interpreted as the LR population) are shown by a solid blue line. As expected, the primordial population remains largely unaltered at large sizes (D > 70 km) but becomes strongly depleted below a few kilometres. In contrast, the fragment population follows a steep slope down to small sizes and begins to dominate the total SFD around D ~ 14 km. This behaviour is in good agreement with the observed trends for the LR population described by Wong & Brown (2015).

Such results provide quantitative support for the hypothesis that the colour bimodality of JTs reflects their collisional history. The R and LR groups, as defined spectrally, can be naturally mapped onto categories of primordial bodies and collisional fragments, respectively, in our model. The fact that the transition between R- and LR-dominated regimes occurs around ~14 km and that the LR population dominates entirely below a few kilometres strengthens the view that small trojans are primarily produced when their parent bodies are located in the L4 trojan region.

Fig. 10 SFDs from our best case run (SL1×0.1, case A). The black dashed line shows the initial SFD of the primordial population at the time of JT emplacement. The solid black line corresponds to the total evolved population after 4 Gyr of collisional evolution. The red line represents the remnant primordial population, while the blue line shows the population of fragments produced during collisions.

4.4 Expected craters SFD of Lucy Mission targets

This section presents the predicted crater size distribution for the Lucy mission targets at the equilibrium point L4 corresponding to the evolution conditions that better match the observed JT SFD, as shown in Sect. 4.1 (case A and SL1 × 0.1). To provide context, we begin by briefly describing the key characteristics of (15094) Polymele, (11351) Leucus, (21900) Orus, and (3548) Eurybates. Afterwards, we introduce the cratering law used in this work and present results for the cratering record on the mentioned targets.

4.4.1 Lucy Mission targets in the L4 swarm

• (15094) Polymele is the smallest Lucy target in L4 and it probably has an ellipsoidal shape (Levison et al. 2023). Its estimated effective diameter is 21 km, with a rotational period of ~11.5 h (Levison et al. 2021). It is a dark, reddish, and P-type asteroid. Stellar occultations revealed that Polymele is accompanied by a small satellite tentatively named Shaun, about 5–6 km in diameter (Buie et al. 2022).

• (11351) Leucus exhibits interesting features, as revealed by stellar occultations (Buie et al. 2021). This trojan slow rotator appears to have an elongated, prolate shape with approximate sizes 63.8 × 36.6 × 14.8 km, corresponding to an equivalent spherical diameter of 41 km. Lucus is a D-type asteroid with a very dark surface and apparently significant topographic features that may suggest a complex collisional history.

• (21900) Orus photometric observations and stellar occultations revealed updated features of this asteroid (Mottola et al. 2023; Buie et al. 2020). Orus D-type asteroid has a retrograde spin with a sidereal rotation period of 13.486 hours, as determined from photometric light-curve observations (Mottola et al. 2023). Stellar occultation data suggest a surface-equivalent spherical diameter of 60.5 ± 0.96 km. Orus is characterised by its very red colour and it appears to have a convex shape, indicative of its surface morphology and potential collisional history.

• (3548) Eurybates photometric observations and stellar occultations provided refined details about this asteroid (Mottola et al. 2023; Buie et al. 2020). This C-type JT has a retrograde spin with a sidereal rotation period of 8.73 h, as determined from photometric light curve observations. Stellar occultation data estimate its surface-equivalent spherical diameter as 69.3 ± 1.4 km. As the largest member of the Eurybates collisional family, it apparently has a convex shape. Additionally, the Hubble Space Telescope identified a 1-km satellite orbiting Eurybates (Noll et al. 2020), likely a remnant of its collisional history.

4.4.2 Cratering scaling law

To estimate any cumulative crater size distribution, we need to calculate the transient crater diameter (D_t) caused by a projectile with size d and the corresponding final crater size (D_c). The Pi-group scaling law formulation provided by Holsapple & Housen (2007) allows us to obtain D_t as a function of the impact conditions and material properties. After the collapse of the transient crater, the final crater is usually estimated as D_c = λD_t (λ ≥ 1 is determined empirically). In fact, in the cratering process, immediately after the crater excavation, a transient crater is formed. Afterwards, the outer structure collapses due to gravity, and the crater is ultimately modified to its final form; the resulting crater is wider than its transient form. Typical values for λ are usually assumed to be in the range 1.2–1.5. Thus, the following relationship allows for the calculation of the diameter of the crater produced by any specific collision. $$D_c=\lambda k_{1}d{[\frac{gd}{2V_P^2}{\left(\frac{\rho }{\delta }\right)}^{\frac{2\nu }{\mu }}+k_2{\left(\frac{Y}{\rho V_P^2}\right)}^{\frac{2+\mu }{2}}{\left(\frac{\rho }{\delta }\right)}^{\frac{\nu (2+\mu )}{\mu }}]}^{-\frac{\mu }{2+\mu }},$$(24)

where V_P is the projectile velocity component perpendicular to the surface, δ is the bulk density, Y is the strength of the target material, k₁, k₂, μ, and ν account for the nature of the target terrain. Finally, g = (4/3) (D/2) πGρ, is the surface gravity, where D is the diameter of the target, G is the gravitational constant, and ρ is the density of the target, while considering that it is spherical and homogeneous.

The knowledge of some parameters in Eq. (24) is limited, especially regarding the composition of JTs, which are likely a combination of rock and ice (Wong et al. 2014). Then, we need to make some assumptions about material properties. Specifically, we set both the density of the material (ρ) and the bulk density (δ) to 1 g cm⁻³, consistent with the assumptions of the collisional evolution model. Holsapple & Housen (2007) pointed out that 1/3 < μ < 2/3 (depending on the energy dissipation of the material) and the experimental results suggest that ν ≃ 0.4 for most materials. Consequently, parameters (k₁, k₂, ν, μ) are assigned suitable values for the target composition of rocks (0.93, 1, 0.4, 0.55) and highly porous materials (0.725, 1, 0.4, 0.4). Finally, the Y parameter remains largely unknown. In this study, we explore a range of strength values for rocky materials: 10 kPa (soft rock) and 1 MPa (hard rock). This choice is justified by the wide range of estimates found in the literature. To provide context for these values, some regions of the lunar regolith have low tensile strengths ranging 1–12 kPa (Slyuta 2014; Ding et al. 2024), while harder rock materials can reach strengths in the 10–18 MPa range (Holsapple 1993). Additionally, for highly porous materials, we assumed 10 kPa.

4.4.3 Cratering record

Here, we present the results for the Lucy targets located in L4, showing the cumulative number of craters per unit area (N_c(> D) [km⁻²]) as a function of crater size. Normalization to the unit area is achieved by calculating the surface area of each body, considered as a sphere of radius equal to its equivalent radius. The effective diameter of each target were taken from Table 3 in Levison et al. (2021). The figures in this section show the outcome corresponding to our preferred parameter choice discussed in Section 4.1 (case A and SL1 × 0.1), for the two tested cratering scaling laws: rock and highly porous material, as specified earlier in this paper.

Figure 11 presents the cumulative crater distribution for the four targets. As expected, all targets follow a similar trend for small craters. However, noticeable differences arise for large craters due to the dependence of crater size on Y and μ in Eq. (24), which governs the transition between strength– and gravity–dominated regimes. For rocky targets, the influence of the strength parameter is particularly evident. Increasing the material strength by a factor of two not only reduces the overall crater size but also extends the range of large craters across different target sizes. In the case of low-strength material (Y = 10 kPa), the surface may become fully saturated with small craters (<0.2 km) and are more than 10% saturated for large craters. Conversely, for high-strength material (Y = 1 MPa), the surface may have significant crater deficiency at all size ranges. Notably, the cratering scaling law for hard rock produces a crater size distribution similar to porous targets with Y = 10 kPa.

Those results highlight the critical role of material properties in shaping the expected crater size distribution on JT surfaces. To explore this further, we focus on Eurybates and apply a cratering scaling law previously used for the icy surfaces of small and mid-sized Saturnian satellites, characterised by (κ₁, κ₂, ν, μ) = (1.67, 0.351, 0.397, 0.38) and Y = 15 kPa (Di Sisto & Zanardi 2013; Rossignoli et al. 2019). As shown in Figure 12, if JTs have icy compositions, their surfaces may reach approximately 10% saturation for craters in the ~1–5 km size range, with even higher saturation levels at both small and large craters.

Fig. 11 Cumulative number of craters per unit surface, on (15094) Polymele, (11351) Leucus, (21900) Orus and (3548) Eurybates. Results gathered from a run under case A (αS = 0.4). The parameter for rocky material are (κ1, κ2, ν, μ)=(0.93,1,0.4,0.55), using Y=10 kPa for soft rock (dashed lines) and Y=1 MPa for hard rock (solid lines). The parameter for porous material are (κ1, κ2, ν, μ)=(0.725,1,0.4,0.4) (dotted lines), grey lines indicate crater saturation levels, with a dashed line indicating 10% and a dotted line 100% of geometric saturation.

Fig. 12 Cumulative number of craters per unit surface area on (3548) Eurybates testing different materials, as indicated in the legend. Results are obtained from a run under case A (αS = 0.4), gray lines indicate crater saturation levels, with the dashed line indicating 10% and the dotted line full (100%) of geometric saturation.

5 Summary and conclusions

In this study, we examine the collisional evolution of the JTs population using a customised version of ALICANDEP called ALICANDEP-JT. Our investigation involves a comprehensive comparison with established models, such as those proposed by de Elía & Brunini (2007) and Marschall et al. (2022), enabling us to both validate our model and explore the relationship between shattering and disruption threshold specific energy. Following the validation and determination of suitable values for the parameters in the shattering scaling laws, we employed our code package to model the collisional evolution of JT under the capture hypothesis. We used an initial SFD derived from recent studies (B25) on the collisional evolution of the primordial trans-Neptunian disc as a baseline. The key findings of this study are as follows.

• Code validation and link between $Q_D^{\ast }$ and $Q_S^{\ast }$: Overall, our comparison of ALICANDEP-JT with previous research yielded positive results. First, we showed that our model can reproduce similar outcomes to de Elía & Brunini (2007), under equivalent fragmentation model and evolutionary conditions. Furthermore, we compared our model with the evolved SFD obtained by Marschall et al. (2022) using the Boulder code. The main difference between the two models is that the Boulder code does not deal with the fragmentation and re-accumulation process; instead, it directly includes the results of previous collision modelling via SPH codes and scales the fragment distribution accordingly. A satisfactory agreement between the two models was achieved for the considered cases, establishing a link between $Q_D^{\ast }$ and $Q_S^{\ast }$ under consistent conditions for the re-accumulation model and the f_ke parameter. The plots in Appendix A highlight the need of including a constant β in the f_ke function to improve the match with M22 in the size range from 100 m to a few kilometres, especially in the case the scaling law has a minimum at 250 m. Our analysis revealed that the location of the minimum in $Q_S^{\ast }$ is shifted by a factor of 2 toward large sizes, compared to the corresponding $Q_D^{\ast }$ minimum.

• Collisional evolution of JT under the capture scenario: Our findings support the hypothesis that JTs originated in the outer disc, before being captured by Jupiter during the instability phase of the giant planets. To model this scenario, we picked a fraction of the evolved outer population at the onset of the instability period, and let it evolve collisionally under the JT boundary conditions. Such fraction corresponds to a capture probability of 6.5 × 10⁻⁷, which is consistent with the capture probability of (5 ± 3) × 10⁻⁷ estimated by Nesvorný et al. (2013). The optimal agreement between the modelled SFDs of JTs after 4500 Myr and the observed population within the L4 swarm allowed us to identify the most suitable scaling law governing the fragmentation process and the preferred slope of the SFD at the small-size end, for the primordial outer population. Our analysis revealed that weak bodies necessarily underwent collisional grinding to match the observed JT SFDs closely. Furthermore, we determined that only a primordial population with a steep slope of the cumulative magnitude distribution (α = 0.4) can accurately reproduce the observed SFD of JTs. These results support the match between the B25 findings and the present-day debiased Kuiper Belt populations. This study shows that the JT underwent a more intense collisional evolution than Kuiper Belt objects. Such a collisional process modified the JT SFD, giving it a more pronounced wave-like structure than the one left in the outer belt by collisional processes during the pre-instability phase.

  Although our findings are in promising agreement with the observed JT population, we acknowledge that other possible combinations of parameters (within limits) may also yield satisfactory results. Multiple factors may cause such a degeneracy. For instance, to build the observed SFD, one needs to make assumptions, for instance, for the mean albedo value. Several intrinsic parameters related to the fragmentation and evolution processes must also be assumed in the numerical modelling, such as the bulk density, impact speed distribution, intrinsic impact probability, and so on (even if for many of these, enough information is already available to reduce the degeneracy).

• Collisional families: among all the simulations performed, only our best-fit scenario (SL1 × 0.1, case A) reproduced the observed number of large collisional families. In this case, five parent bodies with diameters above 90 km underwent catastrophic disruption over 4.5 Gyr. This is in good agreement with the three to five families inferred from the current trojan population (Vokrouhlický et al. 2024). This outcome further supports the plausibility of the SL1 × 0.1 configuration as the most realistic scenario explored in this work.

• Colour sub-populations. Our modelling supports the hypothesis proposed by Wong & Brown (2015); Wong & Brown (2017) that the observed colour bimodality among JTs reflects their collisional evolution. By tracking the primordial and fragment populations separately, we show that the less red (LR) component (interpreted as collisional fragments) follows a steep size-frequency distribution and becomes dominant below ≈14 km in diameter. Conversely, the red (R) component – related to primordial bodies – retains the shape of the initial SFD at large sizes and becomes progressively depleted below a few kilometre sizes due to disruptive collisions. This dual behaviour reproduces the main observational trends and reinforces the interpretation that the R population traces the original captured bodies, while the LR population is largely the result of subsequent collisional processing.

• Cratering records: after our code was validated and the collisional evolution of JT was explored under the capture hypothesis, the ALICANDEP-JT package was used to make predictions about the surface characteristics of the objects that will be visited by the Lucy spacecraft in the L4 swarm. Specifically, we analysed the crater distribution on bodies of sizes comparable to Lucy’s targets. For rocky JTs, low-strength materials (soft rock, Y = 10 kPa) result in a geometric saturation larger than 10%. In contrast, high-strength materials (hard rock, Y = 10 MPa) lead to crater-deficient surfaces. For porous materials, the crater distribution is similar to that of hard rock, while icy compositions imply a distribution that falls between soft and hard rock, reaching a geometric saturation close to 10%.

  Further analyses of Eurybates show that changes in the cratering scaling law parameters, such as μ and k, produce different features in the crater size distribution. Such findings emphasize the need for precise knowledge of the surface composition properties to accurately constrain the population of small JTs from the cratering record, ensuring a proper interpretation of Lucy’s mission data.

  Former studies of the cratering records on JT Lucy targets, such as that of Marchi et al. (2023), have provided valuable information. The parameter sets chosen by Marchi et al. (2023) to represent cohesive soils and porous bodies produce distributions that are relatively similar to each other. Nevertheless, they seem to overlook the substantial variability introduced by significant uncertainties on critical parameter values for the cratering scaling laws. In the present work, we tried to cover this issue by introducing a comprehensive exploration of the parameter space, crucial for capturing the diversity of surface properties of JTs.

Acknowledgements

We thank Raphael Marschall for kindly providing the SFDs from his work, and Romina Di Sisto for her helpful suggestions. PGB and ACB acknowledge support from the Generalitat Valenciana through grant CIAICO/2021/328. PGB also acknowledges funding from the European Space Agency under the Science Faculty program (Funding reference ESA-SCI-SC-LE-191). VSC thanks the University of Coimbra for supporting his stay at the University of Alicante through the Erasmus+ program. TSR acknowledges funding from Ministerio de Ciencia e Innovación (Spanish Government), PGC2021, PID2021-125883NB-C21. This work was (partially) supported by the Spanish MICIN/AEI/10.13039/501100011033 and by ‘ERDF A way of making Europe’ by the ‘European Union’ through grant PID2021-122842OB-C21, and the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia ‘María de Maeztu’) through grant CEX2019-000918-M. We also thank reviewer Vokrouhlický for helpful suggestions on the submitted manuscript.

References

Appendix A Testing f_ke parameters for case 1: Q₀ scaling law (D_min = 250 m)

In this appendix, we show the SFD for various combinations of f_ke parameters listed in Table 2 for each scaling law refereed as case 1. Figures A.1 and A.2 illustrate that different combinations of $C_2^{\mathrm{\prime }}$ and ${\lambda }_2^{\mathrm{\prime }}$, with β = 0 in the f_ke function do not produce a depth as significant as that achieved in the reference run from M22, across the size range of a few meters to several kilometres. However, by assuming a constant value of β = 0.015, our modelled SFDs align closely with the reference SFD from M22 within the same size range (see Figures A.3 and A.4).

Fig. A.1 Evolved SFDs for case 1 (Q0 scaling law with Dmin = 250 m) with parameters indicated in the legend. fke parameters: fke0 = 0.35, γ = 0.7, and β = 0. The black solid line is the reference case of M22.

Fig. A.2 Same as A.1 using fke parameters: fke0 = 0.1, γ = 0.5, and β = 0.

Fig. A.3 Same as A.1 using fke parameters: fke0 = 0.085, γ = 0.5, and β = 0.015.

Fig. A.4 Same as A.1 using fke parameters: fke0 = 0.085, γ = 0.8, and β = 0.015

Appendix B Testing f_ke parameters for case 2: Unconventional scaling law (D_min = 20 m)

Here, we show the SFD for different combinations of f_ke parameters listed in Table 2 for each scaling law refereed as case 2. Figures B.1, B.2, B.3, and B.4 show that the β value is not very critical when the minimum of the scaling law is at ~20 m.

Fig. B.1 Evolved SFDs for case 2 (unconventional scaling law with Dmin = 20 m) with parameters indicated in the legend. fke parameters: fke0 = 0.35, γ = 0.7, and β = 0. The black solid line is the reference case of M22.

Fig. B.2 Same as B.1 using fke parameters: fke0 = 0.1, γ = 0.5, and β = 0.

Fig. B.3 Same as B.1 using fke parameters: fke0 = 0.090, γ = 0.5, and β = 0.010.

Fig. B.4 Same as B.1 using fke parameters: fke0 = 0.090, γ = 0.8, and β = 0.010.

All Tables

All Figures

	Fig. 1 Comparison of two selected scaling laws from EB07, $Q_{S,a}^{\ast }$ = (1.2, 1.75) and $Q_{S,b}^{\ast }$ = (2.3, 1.25), constructed using Eq. (16). For reference, the standard scaling law for $Q_D^{\ast }$ from Benz & Asphaug (1999) is also shown.
In the text

	Fig. 2 Cumulative SFDs for cases $Q_{S,a}^{\ast }$ and $Q_{S,b}^{\ast }$. Dashed lines correspond to the initial population, while solid lines indicate evolved populations with ALICANDEP-JT code, and dotted lines correspond to evolved populations in de Elía & Brunini (2007).
In the text

Fig. 3 $Q_S^{\ast }$ scaling laws tested to study case 1. (a) parameter $C_2^{\mathrm{\prime }}$ is kept constant at value 0.85, while ${\lambda }_2^{\mathrm{\prime }}$ takes three different values. (b) parameter ${\lambda }_2^{\mathrm{\prime }}$ is kept constant at 1.60, while $C_2^{\mathrm{\prime }}$ takes three different values. The solid black line corresponds to the Q0 case, namely the $Q_D^{\ast }$ used by Marschall et al. (2022), for a reference.

In the text

Fig. 4 $Q_S^{\ast }$ scaling laws tested to study case 2. (a) parameter $C_2^{\mathrm{\prime }}$ is kept constant at value 10, while ${\lambda }_2^{\mathrm{\prime }}$ takes three different values. (b) parameter ${\lambda }_2^{\mathrm{\prime }}$ is kept constant at 1.80, while $C_2^{\mathrm{\prime }}$ takes three different values. The solid black line corresponds to the unconventional scaling law, namely, that of $Q_D^{\ast }$, used by Marschall et al. (2022) for a reference.

In the text

	Fig. 5 Functions of fke tested and listed in Table 2. Black symbols correspond to values from Nakamura & Fujiwara (1991); Melosh & Ryan (1997); and Giblin et al. (2004), as indicated in the legend.
In the text

	Fig. 6 Cumulative SFDs for runs corresponding to case 1, using the scaling laws showed in Figs. 3a and 3b). The fke function is determined by fke0 = 0.085, γ = 0.5 and β = 0.015. The solid black line is the evolved cumulative SFD using the scaling law, Q0, from Marschall et al. (2022) for visual comparison.
In the text

	Fig. 7 Cumulative SFD for runs corresponding to case 2, using the scaling law showed in Figs. 4a and 4b. The fke function is determined by fke0 = 0.090, γ = 0.5, and β = 0.010. The solid black line is the evolved cumulative SFD from Marschall et al. (2022) for a visual comparison.
In the text

Fig. 8 Scaling laws tested for the collisional evolution of JTs. As reference, the solid black line represents the unconventional scaling law $Q_D^{\ast }$ from Marschall et al. (2022). The solid blue line (SL1) is the $Q_S^{\ast }$ that best reproduces the evolved SFD by Marschall et al. (2022) (see Sect. 3.3.2). The dashed and dash-dotted blue lines represent SL1 scaled up and down by a factor of 10 and 0.1, respectively.

In the text

Fig. 9 Initial (dashed lines) and final (solid lines) cumulative SFDs of the modelled JT population from each scaling law shown in Fig. 8. Blue and red lines correspond to cases A and B, corresponding to αS = 0.4 and αS = 0.3, respectively. The black solid line represents the cumulative SFD of data from the Minor Planet Center for the L4 swarm, and the dashed gray line corresponds to the SFD estimated by Wong & Brown (2015), including an extrapolation to smaller sizes.

In the text

Fig. 10 SFDs from our best case run (SL1×0.1, case A). The black dashed line shows the initial SFD of the primordial population at the time of JT emplacement. The solid black line corresponds to the total evolved population after 4 Gyr of collisional evolution. The red line represents the remnant primordial population, while the blue line shows the population of fragments produced during collisions.

In the text

Fig. 11 Cumulative number of craters per unit surface, on (15094) Polymele, (11351) Leucus, (21900) Orus and (3548) Eurybates. Results gathered from a run under case A (αS = 0.4). The parameter for rocky material are (κ1, κ2, ν, μ)=(0.93,1,0.4,0.55), using Y=10 kPa for soft rock (dashed lines) and Y=1 MPa for hard rock (solid lines). The parameter for porous material are (κ1, κ2, ν, μ)=(0.725,1,0.4,0.4) (dotted lines), grey lines indicate crater saturation levels, with a dashed line indicating 10% and a dotted line 100% of geometric saturation.

In the text

	Fig. 12 Cumulative number of craters per unit surface area on (3548) Eurybates testing different materials, as indicated in the legend. Results are obtained from a run under case A (αS = 0.4), gray lines indicate crater saturation levels, with the dashed line indicating 10% and the dotted line full (100%) of geometric saturation.
In the text

	Fig. A.1 Evolved SFDs for case 1 (Q0 scaling law with Dmin = 250 m) with parameters indicated in the legend. fke parameters: fke0 = 0.35, γ = 0.7, and β = 0. The black solid line is the reference case of M22.
In the text

	Fig. A.2 Same as A.1 using fke parameters: fke0 = 0.1, γ = 0.5, and β = 0.
In the text

	Fig. A.3 Same as A.1 using fke parameters: fke0 = 0.085, γ = 0.5, and β = 0.015.
In the text

	Fig. A.4 Same as A.1 using fke parameters: fke0 = 0.085, γ = 0.8, and β = 0.015
In the text

	Fig. B.1 Evolved SFDs for case 2 (unconventional scaling law with Dmin = 20 m) with parameters indicated in the legend. fke parameters: fke0 = 0.35, γ = 0.7, and β = 0. The black solid line is the reference case of M22.
In the text

	Fig. B.2 Same as B.1 using fke parameters: fke0 = 0.1, γ = 0.5, and β = 0.
In the text

	Fig. B.3 Same as B.1 using fke parameters: fke0 = 0.090, γ = 0.5, and β = 0.010.
In the text

	Fig. B.4 Same as B.1 using fke parameters: fke0 = 0.090, γ = 0.8, and β = 0.010.
In the text