© The Authors 2025

1 Introduction

All major planets in the Solar System except Mercury feature asteroidal companions in the so-called tadpole or Trojan configuration, where a particle is confined to the vicinity of the Lagrangian triangular equilibrium locations that lead or trail the planet by 60^∘ (Murray & Dermott 1999; Pan & Gallardo 2025). Numerical simulations demonstrate the transient nature of many of these Trojans and, by implication, their likely origin as dynamically unstable planet-crossing small bodies (Christou 2000; Morais & Morbidelli 2002; Connors et al. 2004).

Trojans of Jupiter, Neptune, and Mars have gigayear dynamical lifetimes and likely date from events early in Solar System history (Levison et al. 1997; Scholl et al. 2005; Lykawka & Horner 2010) with Mars featuring the only dynamically long-lived Trojans among terrestrial planets. Unlike Main Belt asteroids, the Mars Trojans (MTs hereafter) experience a relatively benign collisional environment while their proximity to the Sun and long dynamical lifetime make them a useful proving ground for models of asteroid orbital-rotational evolution by radiation forces (Rivkin et al. 2003; Ćuk et al. 2015; Christou et al. 2017).

Their small sizes notwithstanding (≤2 km; Trilling et al. 2007), MTs exhibit some unique features: a marked asymmetry favouring L₅ vs. L₄ residents (de la Fuente Marcos & de la Fuente Marcos 2013; Christou 2013); a strong orbital concentration near (5261) Eureka, suggesting that most L₅ Trojans are products of YORP-induced rotational spin-up and breakup (‘YORPlets’; Christou et al. 2017); and the olivine-dominated composition of this Eureka family (Borisov et al. 2017; Polishook et al. 2017), one of only two known to exist (Galinier et al. 2024). Available compositional clues and dynamical models argue equivocally for either an asteroidal or a Martian origin for MTs (Polishook et al. 2017; Wargnier et al. 2025). Their continued study helps determine the nature of the link – if any – between MTs and the Martian system and our understanding of the early evolution of Mars by implication.

The Minor Planet Center website¹ lists 18 MT asteroids, of which all but four have been confirmed as stable residents by numerical integration of the orbit (Scholl et al. 2005; Christou 2013; de la Fuente Marcos & de la Fuente Marcos 2013; Ćuk et al. 2015; Christou et al. 2020; de la Fuente Marcos & de la Fuente Marcos 2021). Of the two exceptions, 2023 FW₁₄ is the second L₄ Trojan of Mars to be discovered (de la Fuente Marcos et al. 2024). Its short, ∼ 10⁷ yr dynamical lifetime points to a likely origin in the Mars-crosser population. The other asteroid, 2018 FM₂₉, was discovered in 2018 and observed again in 2020, 2022, and 2024. The orbital semimajor axis uncertainty of this asteroid is now several orders of magnitude smaller than the radial width of the MT region (see Sect. 2.2) and we confirmed by numerical integration its status as a stable resident at L₅. The recovery of single-opposition asteroids 2015 TL₁₄₄ and 2016 AA₁₆₅ in 2022 and 2024 respectively place them firmly within the Trojan clouds of Mars and we have similarly confirmed their stability with numerical integration of the orbits. Inclusion of the new Trojans raises the tally of MTs to 18 and of L₅ residents in particular to 16.

Using the methods described in the following Section of the paper, we calculated proper elements – libration width L, eccentricity e_P and inclination I_P – (Table 1) from the latest osculating elements of the objects available through the AstDys² online service (Knežević & Milani 2012). Their distribution, shown in Fig. 1, exhibits several interesting features. Ćuk et al. (2015) note a strong correlation between libration width and inclination among Eureka family asteroids, showing it to be a consequence of the solar-thermal Yarkovsky effect (Bottke et al. 2006). We emphasise it is that correlation, rather than the direct detection of Yarkovsky acceleration in the orbital motion achieved recently for the near-Earth asteroids (e.g. Greenberg et al. 2020), that supports a Yarkovsky-dominated orbital evolution for MTs. Ćuk et al. further note that the asymmetric orbital distribution of family Trojans with respect to Eureka is compatible with a non-gravitational acceleration acting against the direction of motion. This might arise if bi-directional diurnal Yarkovsky is suppressed relative to the strictly dissipative seasonal variant through the coupling of spin and shape evolution (Rubincam 2000; Vokrouhlický & Čapek 2002; Čapek & Vokrouhlický 2004; Statler 2009; Cotto-Figueroa et al. 2015; Ćuk et al. 2015). Secondly, the new MTs have similar eccentricity and inclination to the Eureka family but much higher libration width. Their orbits are in fact similar to that of 2009 SE, an L₅ Trojan that might be unrelated to Eureka (de la Fuente Marcos & de la Fuente Marcos 2021).

The availability of a larger MT sample than was available in earlier works motivates a fresh look at their orbital distribution using tools from dynamical systems and from statistical analysis. The present study uncovers evidence for a more complex evolutionary history of this population than recognised previously, including secondary rotational breakup of YORPlet asteroids and a role for collisions as a competing process for MT population evolution alongside thermal-rotational effects. In the next Section we expose the methods used to compile a new catalogue of MT proper elements, extract interesting features from it, and evaluate their statistical significance. Results are described in Sect. 3 while in Sect. 4 we examine origin scenarios for the new clusters. In Sect. 5 we present our conclusions.

2 Methods

2.1 Confirming the stability of MTs 2009 SE, 2015 TL₁₄₄, 2016 AA₁₆₅ and 2018 FM₂₉

To determine the long-term stability of candidate Trojans, we generated 50 dynamical clones of each asteroid from the state covariance at epoch MJD 60 600.0 available through AstDys using the method of Duddy et al. (2012). The clones plus the nominal orbit were then numerically integrated for 10⁹ yr with the MERCURY package (Chambers 1999) under the gravity of the eight major planets and an integration time step of 4d. For all clones of each asteroid, the critical angle l_r=λ − λ_Mars remained in libration around −60^∘ until the end of the integration, indicating dynamical lifetimes comparable to the Solar System age. Examples of individual clone orbital evolution are shown in Fig. 2 where we also included 2009 SE. In the case of 2015 TL₁₄₄ we observe libration width variations of 5^∘−10^∘ over timescales of 10⁸ yr. We observe such variations in most clones of the asteroid and attribute them to the proximity of the orbit to the chaotic region surrounding the tadpole-horseshoe separatrix at a libration width of ∼ 156^∘ (Morais 1999). Otherwise, the clone integrations imply similar dynamical stability to the other Trojans.

2.2 Modelling tadpole libration of MTs

In the context of the circular restricted Sun-planet-asteroid problem, the evolution of the guiding centre a_r(l_r) in the vicinity of the planetary orbit (|a_r| = |a−a_P| ≪ 1) is described by the integral (Namouni et al. 1999) $$a_r^2=C-\frac{8\mu }{3}S(l_r,e,I,\omega )$$(1) $$S(l_r,e,I,\omega )={\frac{1}{2\pi }{\int }_{-\pi }^{\pi }({|\mathbf{r}-{\mathbf{r}}_P|}^{-1}-\mathbf{r}\cdot {\mathbf{r}}_P)\mathrm{d}\lambda |}_{a_r=0}$$(2)

where C is a constant, r and r_P are the heliocentric position vectors of the asteroid and planet respectively, l_r = λ−λ_P and distances are in units of a_P. As S has a minimum S_min ≃ 1/2 at l_r≃± π/3, Eq. (1) can be written as $$a_r^2=a_0^2-(8\mu /3)\mathrm{\Delta }S$$(3)

where ΔS = S − S_min and $a_0^2$ = C − (8μ/3) S_min. The radial width a₀ for tadpole libration was obtained by evaluating Eq. (3) at the turning points for 0<ΔS < 1. For e ∼ 0, I ∼ 0^∘ and Δ S = 1 we obtain the classical result (Yoder et al. 1983; Murray & Dermott 1999; Morais 1999, 2001) that $$a_{0,max}\simeq \sqrt{(8/3)\mu }{a}_P.$$(4)

For the case of Mars, we used μ ∼ 3.2 × 10⁻⁷, a_P ∼ 1.524 au and found $$a_{0,max}\simeq 0.00141\mathrm{a}\mathrm{u}.$$(5)

At large libration amplitude, the motion of the guiding centre around the equilibrium point is strongly asymmetric (Érdi 1988; Milani 1994). For this reason, we chose to parameterise the resonant state of the asteroid orbit by the full width of l_r libration rather than the libration amplitude D (Milani 1993; Ćuk et al. 2015). We refer to this width as L so that, for small libration amplitude (S ≳ S_min), L ≃ 2 D. We also chose d̄ = a₀/a_{0, max} = Δ S as the action-like quantity, obtained by numerical evaluation of Eq. (3) at a_r = 0 with ω = 0^∘, e = 0, and I = 0^∘. Our approach in determining the resonant elements is similar to that of Vokrouhlický et al. (2024) for the Jupiter Trojans. We report proper element estimates for the MTs in Table 1 while pointing out that the moderate inclination of MTs modifies the effective potential with respect to the planar case (Fig. 3). The principal effect of the change is a planetward displacement of the turning point of the tadpole by up to ∼ 4^∘ and an increase in L by a similar amount, manifesting as a fixed offset in the libration widths of the clusters investigated in Sect. 3.2. Including the correction does not alter the orbital dispersion of MT families (Sects. 3.1 and 3.2) and we chose not to include it in our computation of L to maintain consistency with existing proper element catalogues (e.g. Vokrouhlický et al. 2024).

2.3 Proper element estimation

We calculated the proper eccentricity and inclination by carrying out a ∼7.3 × 10⁵ yr integration of the nominal orbit as described in Sect. 2.1 sampled every 128d. The complex quantities z=e exp i ϖ and ζ = I exp iΩ formed from the output were then fitted to Fourier series of the form $$z(t)=A_0{e}^{\mathrm{i}{g}_{0}t}+\sum _{i=1}^n{A}_i{e}^{\mathrm{i}{\mathrm{v}}_{\mathrm{i}}\mathrm{t}},\zeta (t)=B_0{e}^{\mathrm{i}{s}_0\mathrm{t}}+\sum _{i=1}^n{B}_i{e}^{\mathrm{i}{\mu }_{\mathrm{i}}\mathrm{t}}$$(6)

using the Frequency-Modified Fourier Transform algorithm (FMFT; Sidlichovský & Nesvorný 1997). The coefficients A_i and B_i in Eq. (6) are complex quantities. The frequencies v_i, μ_i with i ≠ 0 are linear combinations of the fundamental Solar System secular eigenmodes (Laskar et al. 2004) while ν₀ and μ₀ are the free precession eigenmodes of the asteroid and $\mathrm{i}=\sqrt{-1}$. Proper elements were recovered from the fit as e_P=|A₀| and I_P=|B₀|.

For the libration width L we adopted the average of the difference between the maximum and minimum value of l_r within a 120-sample window applied to the numerical time series. An example of this procedure is shown in Fig. 4. Finally, we obtained d̄ from L by evaluating Eq. (3) for a circular, planar orbit as described in the previous Section.

2.4 Statistical tests of cluster significance

To establish that any particular asteroid grouping represents a real association, we used Monte-Carlo simulations to form the distribution of an appropriate test statistic. That distribution was then compared against the value of the statistic for the candidate cluster to estimate the probability – or p-value – that the claimed cluster is a chance arrangement between unrelated objects. Two different statistical functions were used, one based on the Minimum-volume Enclosing Sphere (MES) and the other on the Cluster Index (CI).

Given a set of m points in a multi-dimensional Euclidean space, the MES of that set is defined by the property that $$V\left(P_{\text{MES\hspace{0.17em}}}\right)<V(P),\mathrm{\forall }P:{\left\{{\mathbf{x}}_i\right\}}_{i=1\cdots m}\in P$$(7)

where P(x₀, r₀) is the sphere with centre x₀ and radius r₀ while V(P) is the sphere volume, used here as a measure of compactness for a given m-subset of M proper orbits.

On each MC trial, we generated a synthetic set of M random points from some reference distribution representing the null hypothesis of uniformity. We then calculated the MES volume V(P_MES) for every possible m-subset of points and compared it with the MES volume for the candidate asteroid cluster V(P_MES,c)=V_c, recording a ‘hit’ when the condition min_j V(P_j) < η V_c is satisfied, where j runs through the set of subsets. The quantity $$p(M,m)=\frac{\text{Number of hits}}{\text{Number of trials}}$$(8)

then estimates the probability that a cluster as compact as our candidate cluster could have arisen by chance. We introduced the η factor as a convenient way to incorporate the different sources of uncertainty in MES volume calculation. In this paper, we used η = 3.0 as a compromise between applying an overly strict criterion and the risk of underestimating V_c due to the statistical uncertainty in the proper orbits.

The Cluster Index (Liu et al. 2008; Huang et al. 2015, see also Markwardt et al. 2023) is defined as $$CI=\frac{\sum _{j=1}^N\sum _{i\in C_j}{\parallel {\mathbf{x}}_i^{(j)}-{\overline{\mathbf{x}}}^{(j)}\parallel }^2}{\sum _{i=1}^M{\parallel {\mathbf{x}}_i-\overline{\mathbf{x}}\parallel }^2}$$(9)

where x̄^(j) is the mean for the j-th out of N clusters within a sample of M points and x̄ is the overall sample mean. A CI value ≪ 1 implies a strongly-clustered set. The synthetic set of orbits in each Monte Carlo trial was split into N=2 clusters (of size M₁ and M₂ respectively) by the kmeans unsupervised clustering search algorithm (MacQueen 1967) implemented via the FindClusters function in Mathematica v12.0 (Wolfram 2003). The procedure initially assigns random centroids to the k clusters with each point subsequently assigned to its nearest centroid. The CI values for the randomly-generated sets form the cumulative distribution function to be compared with the CI value for the candidate asteroid cluster. The p-value was then estimated in a manner analogous to the MES method, that is $$p(M_1,M_2)=\frac{\text{Number of hits}}{\text{Number of trials}}.$$(10)

Kmeans requires the number of clusters as input. Here we used the Elbow method (Thorndike 1953) to heuristically determine the optimum cluster number for our sample by evaluating the numerator in Eq. (9) as a function of k. This is the Within-Cluster Sum of Squares (WCSS; MacQueen 1967) with the property that it vanishes as k → M. The evaluation at k=2 (Fig. 5) corresponds to an ‘elbow’ – the largest change between consecutive values as well as the largest change between consecutive changes – as the optimum cluster number in this case. We present a full listing of our statistical tests and their outcomes in Table 2.

3 Results

3.1 Eureka family structure and the role of Yarkovsky

A notable feature of the Eureka family is a strong anti-correlation between libration width and proper inclination (Fig. 1, top panel). We attribute this to Yarkovsky-driven coupled evolution of L and sin I, while the evolution of the eccentricity is dominated by random diffusive growth (Ćuk et al. 2015; Christou et al. 2020). Numerical models of orbital evolution via either mode suggest an approximate family age of 1 Gyr when compared to the observed orbital dispersion (Ćuk et al. 2015). Here we revisit these conclusions with the benefit of a larger sample of family asteroids and improved orbit determination due to longer observational arcs. The shape and orientation of the family may be computed from the sample covariance matrix (Johnson & Wichern 2007) $$\mathbf{\Sigma }(\mathbf{x})=\frac{1}{N-1}(\mathbf{x}\mathbf{-}\overline{\mathbf{x}}){(\mathbf{x}\mathbf{-}\overline{\mathbf{x}})}^T$$(11) where x = {[d̄_i, e_i, sin I_i]}_{i=1⋯s M} is a 3 × M matrix with rows the proper elements for the N family members, x̄ is the sample mean and the superscript ‘T’ indicates the transpose. The eigenvalues and eigenvectors of Σ, reported in Table A.1 define the shape and orientation of Gaussian elliptical iso-contours as shown in Fig. 6. The coordinate frame defined by the ellipse principal axes amounts to an anticlockwise rotation by an angle θ = 163.7^∘ from the positive d̄ axis. A statistical measure of association between d̄ and sin I is the Pearson correlation coefficient $$n_{13}=c_{13}/\sqrt{c_{11}{c}_{33}}$$(12)

with values between +1 or −1 (Johnson & Wichern 2007). The two extremes correspond to points lying along a line. For the sample of M = 12 Eureka family members we find $$n_{13}=-0.939,$$(13)

that is a strong inverse association between d̄ and sin I. Using only the seven asteroids considered in Ćuk et al. (2015) we obtain the slightly higher negative correlation of −0.971. In contrast, we find no evidence of correlation of either of those elements with e(|n₂₃|, |n₁₂| <0.05).

These results reinforce the conclusions of Ćuk et al. (2015) on the controlling influence of non-gravitational forces in the orbital distribution of MTs. Moreover, all of the five additional Trojans in this study (blue points in Fig. 6) are located either near Eureka or to higher d̄ and lower sin I. Therefore, the asymmetric distribution of family Trojans relative to Eureka attributed to the dominance of the seasonal Yarkovsky effect in the 2015 study is also confirmed for our 1.7 × larger sample.

Agglomerative hierarchical clustering (Ward 1963) was applied in de la Fuente Marcos & de la Fuente Marcos (2021) to a set of 44 Mars co-orbitals in the space of Tisserand constant (Tisserand 1896) vs. the critical angle l_r. Those authors find that all family members known at that time clustered near Eureka. We applied the same procedure on the Eureka family but now in the space of proper elements with values listed in Table 1, using the same s/w implementation of the algorithm as those authors from the SciPy Python library (Virtanen & et al 2020) through the linkage and dendrogram functions and with a distance function of the form $$\begin{array}{rl}{D}^2=& A\ast (8\mu /3)\ast {(\overline{d}-{\overline{d}}_{PB})}^2+B\ast {(e_P-e_{P,PB})}^2\\ & +C\ast {(\sin I_P-\sin I_{P,PB})}^2.\end{array}$$(14)

Here the subscript ‘PB’ represents the chosen parent body and the first term represents the quantity (a−a_PB)/a_PB for our chosen set of variables. We used the following three choices for the vector (A, B, C):(5/4, 2, 2), (1/4, 2, 2) and (1, 1, 1). The last choice corresponds to the standard Euclidean distance used in de la Fuente Marcos & de la Fuente Marcos, while the other two choices correspond to the traditional Hierarchical Clustering Method (HCM) as applied to Main-Belt asteroids (Zappalà et al. 1990, 1994) and to the Jupiter Trojans (Milani 1993) when single linkage is used. In all our runs, the clustering procedure was started from (5261) Eureka.

Figure 7 shows the dendrogram output for both single and complete linkage. In the single linkage runs, asteroid 2011 SC₁₉₁ was always separated from the rest of the group, probably reflecting its status as the family member with the lowest inclination and highest libration width (Fig. 6). The similarity of dendrograms produced by the Zappalà and Milani metrics suggests a relatively minor contribution from the semimajor axis term in Eq. (14). A compact grouping identified in all cases is composed of 2011 SP₁₈₉, 2018 EC₄ and 2018 FM₂₉ with δ [L, e, I]=0.38^∘ ± 0.21, 0.0017 ± 0.0021, 0.18^∘ ± 0.14^∘ (errors reported at 2-σ significance). Hereafter we refer to this grouping as the 2018 EC₄ cluster after its brightest member and investigate it more closely in the next section.

3.2 New MT clusters

To investigate the statistical significance of the 2018 EC₄ cluster we applied the MES compactness criterion (Sect. 2.4) with a Gaussian null distribution defined from the sample of M = 12 Eureka family asteroids. The cumulative distribution of MES volumes for 1000 Monte Carlo trials and for all possible m = 3 subsets is shown in Fig. 8, top panel. We find a ∼4% probability that a grouping more compact than these three asteroids could have occurred by chance within the Eureka family. If we include nearby asteroid (311999) to this group and repeat the procedure, we find that this larger grouping is no longer significant (p-value >0.5). We therefore conclude that the EC₄ cluster is the largest statistically significant sub-grouping of family asteroids under the current observational completeness of the family. We return to the relation between this cluster and (311999) in Sect. 4.

Turning now to Trojans outside of the Eureka family, we observe that the libration widths of two of the new stable Trojans highlighted in Table 1 are much higher than any Eureka family member but similar to that of 2009 SE, with the proper orbit of 2016 AA₁₆₅ in particular being almost identical to that of 2009 SE (δ [L, e, I] = 0.36^∘ ± 0.96^∘, 0.0002 ± 0.0021, 0.23^∘ ± 0.14). Asteroid 2015 TL₁₄₄ has the highest libration width of any stable MT to-date, including (121514) at L₄ (Table 1), but otherwise its eccentricity and inclination are similar to those of the other two asteroids (δ [e, I] = 0.0027 ± 0.0021, 0.34^∘ ± 0.14^∘). The strongly elongated shape of this orbital grouping (σ_{d̄ ≫ σe}, σ_{sin I}) renders the MES approach unsuitable and we resort to the CI-based approach to assess its statistical significance as a cluster of related asteroids.

Application of the kmeans algorithm with k = 2 separated the population of M₁ + M₂ = 15 L₅ asteroids into the candidate cluster and the Eureka family with a CI value >99.9% significant against the null hypothesis (Table 2). To test the robustness of the result, we replaced the Gaussian d̄ variate with a uniform variate drawn between 0 and 0.6 and repeated the procedure. We then substituted the full Eureka family with only the 3 brightest family members to force size parity between the two clusters. We again find that kmeans separates the clusters with >95% significance (Fig. 8, middle panel). Identification of this grouping as a separate cluster is therefore robust and these small asteroids must be related in some way. For the remainder of the paper, we identify this cluster by its brightest member, 2009 SE.

In additional clustering experiments, we found that reapplication of kmeans but now with k=3 again separates out the 2009 SE cluster and additionally a Trojan pair formed by the Eureka family members with highest e_P in Table 1. Applying instead the DBSCAN clustering algorithm (Ester et al. 1996) yielded the same result except that the new cluster was composed of the three highest e_P Trojans. The large eccentricity and inclination spread among members of this cluster (δ e_P ≃ 0.023, δ I_P ≃ 2.6^∘) however indicates that its identification was spurious and we do not consider it further here.

We further investigated the very close pair 2016 AA₁₆₅ − 2009 SE by using the MES approach. As for the full cluster (Table 2), its statistical significance depends on the initial sample size M. For M=15 we obtained a significance of 91%, however this assumption ignores the strong statistical association between Eureka family asteroids. Reducing the sample size to M=6 (Fig. 8, bottom panel) and then M=4 as done previously increased significance to 96 and 99.5% respectively, on a par with our earlier conclusions for either cluster.

As a final note in this Section, we point out that the asteroid sizes in these clusters are comparable to – though at the lower end of – members of identified pairs and clusters in the Main Belt (Vokrouhlický & Nesvorný 2008; Pravec et al. 2010, 2018). MB clusters are recognised by virtue of remarkably similar osculating orbits between members. The new Trojan clusters do not show this degree of similarity (Table 3) and are likely older than their MB counterparts (see next section). We postulate that they remained compact and identifiable by virtue of similar proper orbits owing, on one hand, to the overall sparsity of MT asteroids and, on the other, their spatial confinement 60^∘ ahead or behind Mars which mitigates differential orbital evolution by inverse-distance perturbations from the other planets (Christou 2013).

4 Discussion: cluster origin scenarios

4.1 2018 EC₄ cluster: A potential product of cascade disruption?

Asteroids in both new groups occupy a relatively narrow range in H implying similar-sized objects. We can estimate asteroid size through the expression $$D(\mathrm{k}\mathrm{m})=\frac{1329}{\sqrt{p_V}}{10}^{-H/5}$$(15)

if the geometric albedo p_v is known (Bowell et al. 1989). Here we assumed p_V=0.2 similar to other Eureka family asteroids (Christou 2013). For the EC₄ cluster we obtained nominal diameter estimates of 276 m (2018 EC₄), 191 m (2011 SP₁₈₉) and 175 m (2018 FM₂₉), yielding same-density mass ratios in the range 0.25–0.77, generally high for products of YORP-induced rotational fission (mass ratio ∼0.2; Pravec et al. 2010, 2018).

Alternatively, the parent body of this cluster is another Eureka family asteroid. The nearest family member is (311999) (Fig. 1) with D ∼ 0.7 km. A genetic relationship with this asteroid is more compatible with the rotational fission mechanism. Even if the group separated from (311999) as a single object, the inferred mass ratio of 0.094 fits comfortably within the range predicted by fission theory. A relatively young age for this cluster as implied by its apparent compactness is also consistent with separation from (311999). This asteroid is likely a Eureka YORPlet itself, therefore a link to the cluster implies an ongoing process of cascade disruption of the Eureka progenitor, similar to that invoked for some Main Belt clusters (Fatka et al. 2020).

By comparing the orbital distribution of numerically propagated test particles under the Yarkovsky effect to the observed distribution of Trojans, Ćuk et al. (2015) obtain an approximate upper limit on the family age of 1 Gyr, corresponding to the separation event responsible for 2011 SC₁₉₁ (point at extreme bottom right in Fig. 6). We make use of that result here to estimate the age of the EC₄ cluster. We assumed that the separation distance s in d̄ − sin I space³ of the asteroid from the parent body after time t_age is related to the separation rate ṡ by $$s/t_{age}=(\dot{s}-{\dot{s}}_{PB}).$$(16)

If the same inverse-linear size dependence holds for ṡ as for the semimajor axis drift ȧ in the absence of resonance, Eq. (16) becomes $$s/t_{age}=\dot{s}(1-D/D_{PB}).$$(17)

From the analysis in Sect. 3.1 we have s_SC191 ≃ 0.0536. Adopting p_V = 0.4 for Eureka family asteroids as in Ćuk et al., Eqs. (15) and (17) yield $${\dot{s}}_{SC191}=0.070{\mathrm{G}\mathrm{y}\mathrm{r}}^{-1}$$(18)

and, by applying Eq. (17) to (311999) and each of the EC₄ cluster asteroids we obtain t_age = 44−127 Myr. As the Yarkovsky drift rate depends on a number of unknown bulk and surface characteristics of the asteroid (Bottke et al. 2006), the size-corrected value obtained here may differ from that of 2011 SC₁₉₁. An additional uncertainty factor of ∼2 arises from the statistical uncertainty of the computed proper elements (Table 1). However, even if the true separation ages are several times shorter than our estimated lower bound, this cluster appears to be significantly older than other previously identified clusters or pairs of small main-belt asteroids (≲ 4 Myr; Pravec et al. 2018, 2019).

4.2 2009 SE cluster as the outcome of collisional fragmentation

The similar eccentricity and inclination between this cluster and the Eureka family on one hand but much higher libration width on the other, could arise if this cluster is a remnant of an ejecta cloud from a past collision on an MT. Such ejecta would escape from the surface of a parent body of diameter D and bulk density ρ if v_ej > v_esc where $${\mathrm{v}}_{\mathrm{e}\mathrm{s}\mathrm{c}}=D\sqrt{(2/3)G\pi \rho }$$(19)

and G is the universal gravitational constant. The separation velocity of magnitude ${\mathrm{v}}_{\text{sep}}=\sqrt{{\mathrm{v}}_{\text{ej}}^2-{\mathrm{v}}_{\text{esc}}^2}$ for escaping fragments would cause differences between fragment and parent body orbits given by the Gauss equations (e.g. Nesvorný et al. 2006): $$\delta a/a=\frac{2}{{\mathrm{v}}_{\text{circ\hspace{0.17em}}}{(1-e^2)}^{1/2}}[(1+e\cos f){\mathrm{v}}_t+(e\sin f){\mathrm{v}}_r]$$(20) $$\delta e=\frac{\sqrt{1-e^2}}{{\mathrm{v}}_{\text{circ\hspace{0.17em}}}}[\frac{e+2\cos f+e{\cos }^{2}f}{1+e\cos f}{\mathrm{v}}_t+(\sin f){\mathrm{v}}_r]$$(21) $$\delta I=\frac{\sqrt{1-e^2}}{{\mathrm{v}}_{\text{circ\hspace{0.17em}}}}\frac{\cos (\omega +f)}{1+e\cos f}{\mathrm{v}}_{\mathrm{w}}$$(22)

where the orbit elements are those of the parent asteroid, ${\mathrm{v}}_{\text{circ}}=\sqrt{G{M}_{\odot }/a}$ is the heliocentric speed for a circular orbit of radius a while v_t, v_r and v_w are the respective transverse, radial and out-of-plane components of v_sep. For e∼0 and a = a_P, Eqs. (4) and (20) give $${\mathrm{v}}_{\text{sep},0}\simeq 11.2\mathrm{m}{\mathrm{s}}^{-1}$$(23)

for the minimum velocity ‘kick’ needed to evict a small-amplitude (d̄∼0) MT. A kick of this magnitude would, at the same time, have a negligible effect on the other elements (δe, δI = O(10⁻³) from Eqs. (21) and (22)). Figure 9 shows the distribution of randomly-generated particles that separated from Eureka and (311999) with fixed separation speeds of 5, 20 and 200 m s⁻¹. Their orbits were obtained from Eqs. (20)–(22) by sampling ω and f uniformly over the range 0-2π and for an isotropic separation direction. Values for the remaining elements were set to be those of the respective parent asteroids while L was calculated from δa using the method of Vokrouhlický et al. (2024).

We observe that a separation speed of 5 m s⁻¹ from either (5261) or (311999) is sufficient to increase L to the values in the 2009 SE cluster (top panel). This scenario is compatible with separation from (311999) in particular, except that an additional mechanism is then required to increase e from 0.045 to ∼ 0.06 to match that of (311999) (bottom panel). Delivering such a change in fragment eccentricity from an impact on the asteroid requires a separation speed of 200 m s⁻¹, much higher than typical fragment ejection speeds from collisional disruption of km-sized bodies (≤10 m s⁻¹; Jutzi et al. 2010).

An additional factor to consider is the timing of the impact event and subsequent orbital evolution of the fragments. We have seen (Sect. 3.1 and Ćuk et al. 2015) that small MTs would gradually migrate to orbits with lower inclination and higher libration width due to the Yarkovsky effect. Low-speed ejection of the 2009 SE cluster asteroids from (5261) sometime in the past readily delivers the asteroids to their observed orbits without the need to invoke unusually high ejection speeds. Since the cluster’s proper inclination is significantly lower from Eureka, the time elapsed since the collision event would be comparable to the Eureka family age. Assuming a similar rate of orbital migration and working as in Sect. 4.1, we find t_age ≃ 400 Myr.

The origin scenario we propose for this cluster is compatible with the one available model of the likely collisional history of the Trojans (Christou et al. 2017). The model gives about even odds that Eureka experienced no catastrophic disruption (Q ∼ Q^* where Q is the impactor kinetic energy per unit target mass and Q^* is an impactor- and target-dependent threshold; Housen & Holsapple 1990; Benz & Asphaug 1999) over the age of the Solar System. Impacts with Q < Q^* should occur more frequently and we propose that at least one such event capable of imparting a speed of ∼ 5 m s⁻¹ to fragments 1/10 the size of Eureka occurred within the last ∼1 Gyr of Solar System history.

Interestingly, this timeframe is also the expected collisional lifetime of D∼300m MTs such as 2009 SE and 2018 EC₄ (Christou et al. 2017). The existence of either clusterings might therefore signify their relatively recent creation from progenitor bodies derived from Eureka, (311999) or both. In some cases the member separation epoch in young asteroid families may be constrained by backwards propagation of the orbits (Nesvorný et al. 2002; Carruba et al. 2018). We relegate such an exercise for the MT families to future work.

5 Conclusions

• We identify three new faint (H=20−21) L₅ MTs, raising the total number of trailing Trojan cloud residents to 16. Among the new finds is 2015 TL₁₄₄ with the highest width of tadpole-type libration (∼ 83^∘) among stable MTs.

• The new Trojans have similar eccentricity and inclination to the previously identified Eureka family of MT asteroids (Christou 2013; de la Fuente Marcos & de la Fuente Marcos 2013). However, two of the objects have quite high libration width and may not be rotational fission products (‘YORPlets’) as suggested for Eureka family asteroids (Christou 2013; Ćuk et al. 2015; Christou et al. 2020). We propose instead that these objects, together with 2009 SE identified previously by de la Fuente Marcos & de la Fuente Marcos (2021), originated as collisional impact ejecta from a Eureka family asteroid with the most likely parent being (5261) Eureka itself. We also identify an additional statistically significant cluster of similarly faint Eureka family members that may have formed as YORPlets of asteroid (311999) sometime in the last 10⁸ yr.

• We revisit the significance of the correlation between libration amplitude and inclination attributed to Yarkovsky-driven orbital evolution (Ćuk et al. 2015) by calculating the Pearson correlation coefficient for the Eureka family asteroids. Our estimate of 94% implies a strong association between these two orbital properties. We further confirm the lack of Eureka-derived YORPlets with orbital evolution dominated by a positive Yarkovsky acceleration, consistent with the proposed dominance of the seasonal over the diurnal Yarkovsky effect over the ∼ 1 Gyr age of the family (Ćuk et al. 2015).

While many families of 1 km across or larger Main Belt asteroids are now known (Nesvorný et al. 2015, 2024), the high MB asteroid number density and the overlap between families makes it difficult to identify small asteroid clusters older than ∼ 10⁷ yr (Pravec et al. 2018; Vokrouhlický et al. 2024). At the same time, MT proximity to Earth combined with the relative sparsity of their number and an environment free from frequent collisions and from planetary close approaches highlight this population as a unique natural laboratory to study small asteroid evolution.

The decadal Legacy Survey of Space and Time (LSST; Jones et al. 2016) will discover and catalogue approximately 5 million main belt asteroids in addition to the ≳1 million currently known. Assuming similar discovery efficiency and a size distribution n (< D) ∝ D^−0.45 for MTs (Christou et al. 2020), LSST should find ∼ 70 new MTs as faint as H=23 (D∼75 m for a geometric albedo of 0.2), allowing to test the conclusions of this paper by comparing the orbital distributions of those smaller asteroids and of existing MTs. The outcome of this exercise will have wider implications for the efficiency of diurnal vs. seasonal Yarkovsky and the relative roles of radiation forces vs. collisions for sub-km asteroids on timescales 10⁷−10⁹ yr. This is a regime where few observational constraints exist, yet extremely relevant to the important coupled problems of NEA/meteorite source regions and the migration routes of asteroid fragments to the Earth.

Acknowledgements

A.D.O. was supported by a grant of the Italian National Institute for Astrophysics for fundamental researches projects (INAF, act n. 38/2023). Work by A.H. was supported by the Science and Technology Facilities Council (STFC). We would like to thank the High Performance Computing Resources team at New York University Abu Dhabi and especially Jorge Naranjo for helping us with our numerical simulations. Insightful comments by the anonymous reviewer significantly improved the manuscript. Astronomical research at the Armagh Observatory & Planetarium is grant-aided by the Northern Ireland Department for Communities (DfC).

Appendix A Statistical analysis of Eureka family proper orbits

Table A.1 shows the covariance and correlation of the sample of proper orbits in Table 1 as well as the covariance matrix eigenvalues and eigenvectors used in Sect. 3.1.

Appendix B Uncertainty modelling of MT proper element estimates

Formal uncertainties reported in Table 1 were calculated as $${\sigma }_x=\sqrt{{\sigma }_{x,\mathrm{o}}^2+{\sigma }_{x,\mathrm{m}}^2}$$(B.1)

where x = L, e_P or I_P, the first term under the square root refers to the orbit uncertainty and the second to the fitting model. Estimates of σ_{x, o} and σ_{x, m} for the MTs are shown in Table B.1.

For the libration width, we adopted as model uncertainty σ_{L, m} the 2-σ dispersion of the 120-sample estimates obtained from the numerical time series (Sect. 2.3). To calculate the error introduced by the filtering procedure in estimating σ_{e, m} and σ_{I, m}, we applied FMFT to selected MTs and for different Fourier polynomial order n (2, 5, 10, 15 and 20). Fig. B.1 shows these estimates shifted vertically to match the value at n = 10. Consequently, we adopted the values 0.0015 and 0.1^∘ containing 95% of e_P and I_P estimates as respective model errors for all MTs.

To estimate the uncertainty arising from the finite knowledge of the orbit, the asteroids were first ranked according to the semimajor axis uncertainty σ_a, as the coorbital state is most sensitive to this element. This information was obtained from the state covariance for epoch JD2460600.5 from AstDys. These values ranged between ∼2 × 10⁻⁹ au (5261) to ∼2 × 10⁻⁶ au (2009 SE). Each asteroid was then assigned to one of three bins demarcated by the values σ_a = 5 × 10⁻⁷ au and σ_a = 5 × 10⁻⁸ au. We selected one representative for each bin, calculated proper elements for 20 dynamical clones of each and adopted 95% confidence intervals from the 20-sample statistics as our estimates of σ_{x, 0} for all asteroids in the respective bins. The bin representative asteroids in order of increasing σ_a were (5261), 2009 SE and 2015 TL₁₄₄.

References

---

All Tables

All Figures