A pair of Jovian Trojans at the L4 Lagrange point
Timothy R. Holt, David Vokrouhlický, David Nesvorný, Miroslav Brož, Jonathan Horner
aa r X i v : . [ a s t r o - ph . E P ] S e p MNRAS , 1–22 (2020) Preprint 16 September 2020 Compiled using MNRAS L A TEX style file v3.0
A pair of Jovian Tro jans at the L4 Lagrange point
Timothy R. Holt, , ⋆ David Vokrouhlick´y, David Nesvorn´y, Miroslav Broˇz andJonathan Horner Centre for Astrophysics, University of Southern Queensland, Toowoomba, Queensland 4350, Australia Southwest Research Institute, Department of Space Studies, 1050 Walnut St., Ste 300, Boulder, CO-80302, USA Institue of Astronomy, Charles University, V Holeˇsoviˇck´ach 2, Prague 8, CZ-180 00, Czech Republic
Accepted XXX. Received YYY; in original form ZZZ
ABSTRACT
Asteroid pairs, two objects that are not gravitationally bound to one another, but sharea common origin, have been discovered in the Main belt and Hungaria populations.Such pairs are of major interest, as the study of their evolution under a variety ofdynamical influences can indicate the time since the pair was created. To date, noasteroid pairs have been found in the Jovian Trojans, despite the presence of severalbinaries and collisional families in the population. The search for pairs in the JovianTrojan population is of particular interest, given the importance of the Trojans astracers of planetary migration during the Solar system’s youth. Here we report adiscovery of the first pair, (258656) 2002 ES and 2013 CC , in the Jovian Trojans.The two objects are approximately the same size and are located very close to theL4 Lagrange point. Using numerical integrations, we find that the pair is at least Myr old, though its age could be as high as several Gyrs. The existence of the(258656) 2002 ES –2013 CC pair implies there could be many such pairs scatteredthrough the Trojan population. Our preferred formation mechanism for the newlydiscovered pair is through the dissociation of an ancient binary system, triggered bya sub-catastrophic impact, but we can not rule out rotation fission of a single objectdriven by YORP torques. A by-product of our work is an up-to-date catalog of JovianTrojan proper elements, which we have made available for further studies. Key words: minor planets , asteroids: general
The discovery of asteroid pairs, two objects sharing avery similar heliocentric orbit, recently brought yet an-other piece of evidence into the mosaic of small So-lar system bodies’ evolution on short timescales (e.g.,Vokrouhlick´y & Nesvorn´y 2008). Examples of these coupleshave been found in the Main belt and Hungaria populations(Vokrouhlick´y & Nesvorn´y 2008; Pravec & Vokrouhlick´y2009; Ro˙zek et al. 2011; Pravec et al. 2019). The similaritybetween the heliocentric orbits of the two members of anidentified asteroid pair hints at a common and recent ori-gin for the objects, that most likely involves their gentleseparation from a parent object. Indeed, backward orbitalpropagation of heliocentric state vectors of the componentsin many pairs has allowed researchers to directly investigatethe possibility of their past low-velocity and small-distanceapproach (see Vokrouhlick´y et al. 2017, for the most out-standing example discovered so far). ⋆ E-mail: [email protected] (TRH)
The well-documented cases of pairs among asteroidsidentified to date all feature separation ages of less thana million years. Vokrouhlick´y & Nesvorn´y (2008) speculatedabout three processes that could have led to the formation ofthose paris: (i) collisional break-up of a single parent object,(ii) rotational fission of such an object driven by radiationtorques, and (iii) instability and separation of the compo-nents of a binary system. Whilst each of these possibilitiescan explain the origin of asteroid pairs, with some beingmore likely than others for individual pair cases, evidencehas been found that the majority of currently identified pairswere probably formed through the rotational fission of theirparent object (e.g., Pravec et al. 2010, 2019). It is worthnoting that Main belt binaries in the same size category(i.e., with primary diameters of one to a few kilometers), arealso believed to be primarily formed through the rotationalfission of their parent body (e.g., Pravec & Harris 2007;Margot et al. 2015a). This is an interesting population-scaleresult that informs us about a leading dynamical process forfew-km size asteroids in the Main belt. It would certainly be © T. R. Holt et al. desirable to extend this knowledge to other populations ofsmall Solar system bodies.Attempts to detect orbital pairs in other populationshave, to date, either failed or were not strictly convincing.For instance, the orbital evolution of bodies in the near-Earth population is very fast and chaotic and, at the sametime, the number of known objects is limited (see, e.g.,Moskovitz et al. 2019, and references therein). Searches inpopulations beyond the Main belt were not successful fordifferent reasons. Whilst dynamical chaos could also be rel-evant, a more important factor concerns the smallest sizeof bodies found at larger distance from the Sun. The small-est bodies found in Cybele zone, and amongst the Hildasor Jovian Trojans, are about an order of magnitude largerthan the smallest known asteroids in the inner Main beltor the Hungarias (e.g. Emery et al. 2015). The proposedpair-formation processes have a characteristic timescale thatrapidly increases as a function of parent body size. For thatreason, it is no surprise that, to date, no recently formed( ≤ My) traditional pairs sharing the same heliocentric or-bit have been detected beyond the Main belt. If any pairs doexist in these distant small-body populations, they should berevealed by their tight configuration in proper element spaceand long-term backward orbital propagation, if the stabilityin that particular zone of orbital phase space allows. Withthat guideline in mind, we focus here on the Jovian Trojanpopulation. The leap to the Trojan population might appearto contradict the logical steps of gradually extending ourknowledge of Main belt pairs by searches among the Cybeleor Hilda populations first. However, we argue that the caseof possible Jovian Trojan pairs is actually more interestingbecause of that population’s entirely different origin.The Jovian Trojan population consists of two swarmsof objects, librating on tadpole trajectories about the Jo-vian L4 and L5 Lagrange points. Indeed, 588 Achilles Wolf(1907) was the first discovered object to serve as an exampleof a solution to the restricted three body problem (Lagrange1772). Whilst originally considered to be just an extension ofthe main belt, and particularly the Hilda and Thule popula-tions, towards the orbit of Jupiter, the Jovian Trojans weresoon realised to be a totally distinct group of objects, witha unique history(see Emery et al. 2015, for a review). Mostimportantly, the majority of the Jovian Trojans are thoughtto have formed in a vast trans-Neptunian disk of planetesi-mals, at a heliocentric distance beyond ≃ au, and becamecaptured onto their current orbits during the planetesimal-driven instability of giant planets (see Nesvorn´y 2018, for re-view). The physical properties of the Trojans, such as theirmaterial strength or bulk density, are therefore most likelydifferent from most of the asteroidal populations, resem-bling rather those of comets and Centaurs with which theyshare the birth-zone. Though relatively stable, the JovianTrojans can escape their stable region (e.g. Di Sisto et al.2014; Holt et al. 2020, and references therein), and con-tribute to other populations, most notably the Centaurs (seeDi Sisto et al. 2019, and references therein). An example ofthis, (1173) Anchises, exhibits significant dynamical insta-bility on timescales of hundreds of millions of years, with theresult that it will likely one day escape the Jovian Trojanpopulation and become a Centaur before being ejected fromthe Solar system, disintegrating, or colliding with one of theplanets (Horner et al. 2012). Despite their importance as a source of information onthe Solar system’s past evolution, fact that the Jovian Tro-jans are markedly farther from Earth than the Main Belt hasmade them significantly more challenging targets for study.As a result, our knowledge of the collisional history, binarity,and the presence/absence of pairs in the Trojan populationremains far smaller than our knowledge of the main Asteroidbelt (e.g. Margot et al. 2015b). In fact, to date, no confirmedTrojan pairs have been discovered, and the true level of bi-narity in the population remains to be uncovered. The mostfamous confirmed binary in the Trojan population is (617)Patroclus, accompanied by a nearly equal size satellite Me-noetius (both in the 100 km range; e.g., Marchis et al. 2006;Buie et al. 2015). The Patroclus-Menoetius system is fullyevolved into a doubly synchronous spin-orbit configuration(see Davis & Scheeres 2020, and references therein), and rep-resents an example of the kind of binary systems which areexpected to be common among Trojans. A number of suchbinaries, comprising two components of almost equal size,have been found amongst the large trans-Neptunian objects(e.g., Noll et al. 2020). This comparison is of particular in-terest, given that the Patroclus system was, in all likelihood,implanted to the Trojan region from the trans-Neptunian re-gion source zone (e.g., Nesvorn´y et al. 2018). It seems likelythat the Patroclus system represents the closest example ofan Edgeworth-Kuiper belt binary system. Further informa-tion on the Patroclus system will become available in thein coming decades, as the binary is a target for flyby in2033 by the Lucy spacecraft (e.g., Levison et al. 2017). Sim-ilar smaller-scale systems may well exist among the Trojanpopulation , but their abundance is uncertain. Observation-ally, such small-scale binaries remain beyond our detection,and theoretical models of their survival depend on a num-ber of unknown parameters (e.g., Nesvorn´y et al. 2018, 2020;Nesvorn´y & Vokrouhlick´y 2019). The existence of Trojan bi-naries is interesting by itself, but in the context of our work,it is worth noting that, if such binaries exist, they likely serveas a feeding cradle for a population of Trojan pairs.Following this logic, then if the population of pairsamong the Trojans can become known and well character-ized, such that their dominant formation process is under-stood, that would in turn prove to be a source of new in-formation about Trojan binaries. Milani (1993) in his pi-oneering work on Jovian Trojan orbital architecture noteda case of L4-swarm objects (1583) Antilochus and (3801)Thrasymedes. Their suspicious orbital proximity led the au-thor to suggest that they may constitute a genetically re-lated couple of bodies. A viable formation process would bethrough the instability and dissociation of a former binary(Milani and Farinella, personal communication). Unfortu-nately, the Antilochus–Thrasymedes interesting configura-tion has not since been revisited, nor further studied in amore detail.This background information motivates us to conducta search for Jovian Trojan pairs. Unfortunately, even nowthe problem is not simple, and we consider our work to bean initial attempt, rather than providing a definitive solu-tion. In section 2, we explain our strategy, and describe thedifficulties in Trojan pair identification. This strategy led uspreliminarily identify the Jovian Trojans (258656) 2002 ES and 2013 CC as a potential pair. To test this hypothesis,we attempted to prove that these two bodies could be geneti- MNRAS000
The well-documented cases of pairs among asteroidsidentified to date all feature separation ages of less thana million years. Vokrouhlick´y & Nesvorn´y (2008) speculatedabout three processes that could have led to the formation ofthose paris: (i) collisional break-up of a single parent object,(ii) rotational fission of such an object driven by radiationtorques, and (iii) instability and separation of the compo-nents of a binary system. Whilst each of these possibilitiescan explain the origin of asteroid pairs, with some beingmore likely than others for individual pair cases, evidencehas been found that the majority of currently identified pairswere probably formed through the rotational fission of theirparent object (e.g., Pravec et al. 2010, 2019). It is worthnoting that Main belt binaries in the same size category(i.e., with primary diameters of one to a few kilometers), arealso believed to be primarily formed through the rotationalfission of their parent body (e.g., Pravec & Harris 2007;Margot et al. 2015a). This is an interesting population-scaleresult that informs us about a leading dynamical process forfew-km size asteroids in the Main belt. It would certainly be © T. R. Holt et al. desirable to extend this knowledge to other populations ofsmall Solar system bodies.Attempts to detect orbital pairs in other populationshave, to date, either failed or were not strictly convincing.For instance, the orbital evolution of bodies in the near-Earth population is very fast and chaotic and, at the sametime, the number of known objects is limited (see, e.g.,Moskovitz et al. 2019, and references therein). Searches inpopulations beyond the Main belt were not successful fordifferent reasons. Whilst dynamical chaos could also be rel-evant, a more important factor concerns the smallest sizeof bodies found at larger distance from the Sun. The small-est bodies found in Cybele zone, and amongst the Hildasor Jovian Trojans, are about an order of magnitude largerthan the smallest known asteroids in the inner Main beltor the Hungarias (e.g. Emery et al. 2015). The proposedpair-formation processes have a characteristic timescale thatrapidly increases as a function of parent body size. For thatreason, it is no surprise that, to date, no recently formed( ≤ My) traditional pairs sharing the same heliocentric or-bit have been detected beyond the Main belt. If any pairs doexist in these distant small-body populations, they should berevealed by their tight configuration in proper element spaceand long-term backward orbital propagation, if the stabilityin that particular zone of orbital phase space allows. Withthat guideline in mind, we focus here on the Jovian Trojanpopulation. The leap to the Trojan population might appearto contradict the logical steps of gradually extending ourknowledge of Main belt pairs by searches among the Cybeleor Hilda populations first. However, we argue that the caseof possible Jovian Trojan pairs is actually more interestingbecause of that population’s entirely different origin.The Jovian Trojan population consists of two swarmsof objects, librating on tadpole trajectories about the Jo-vian L4 and L5 Lagrange points. Indeed, 588 Achilles Wolf(1907) was the first discovered object to serve as an exampleof a solution to the restricted three body problem (Lagrange1772). Whilst originally considered to be just an extension ofthe main belt, and particularly the Hilda and Thule popula-tions, towards the orbit of Jupiter, the Jovian Trojans weresoon realised to be a totally distinct group of objects, witha unique history(see Emery et al. 2015, for a review). Mostimportantly, the majority of the Jovian Trojans are thoughtto have formed in a vast trans-Neptunian disk of planetesi-mals, at a heliocentric distance beyond ≃ au, and becamecaptured onto their current orbits during the planetesimal-driven instability of giant planets (see Nesvorn´y 2018, for re-view). The physical properties of the Trojans, such as theirmaterial strength or bulk density, are therefore most likelydifferent from most of the asteroidal populations, resem-bling rather those of comets and Centaurs with which theyshare the birth-zone. Though relatively stable, the JovianTrojans can escape their stable region (e.g. Di Sisto et al.2014; Holt et al. 2020, and references therein), and con-tribute to other populations, most notably the Centaurs (seeDi Sisto et al. 2019, and references therein). An example ofthis, (1173) Anchises, exhibits significant dynamical insta-bility on timescales of hundreds of millions of years, with theresult that it will likely one day escape the Jovian Trojanpopulation and become a Centaur before being ejected fromthe Solar system, disintegrating, or colliding with one of theplanets (Horner et al. 2012). Despite their importance as a source of information onthe Solar system’s past evolution, fact that the Jovian Tro-jans are markedly farther from Earth than the Main Belt hasmade them significantly more challenging targets for study.As a result, our knowledge of the collisional history, binarity,and the presence/absence of pairs in the Trojan populationremains far smaller than our knowledge of the main Asteroidbelt (e.g. Margot et al. 2015b). In fact, to date, no confirmedTrojan pairs have been discovered, and the true level of bi-narity in the population remains to be uncovered. The mostfamous confirmed binary in the Trojan population is (617)Patroclus, accompanied by a nearly equal size satellite Me-noetius (both in the 100 km range; e.g., Marchis et al. 2006;Buie et al. 2015). The Patroclus-Menoetius system is fullyevolved into a doubly synchronous spin-orbit configuration(see Davis & Scheeres 2020, and references therein), and rep-resents an example of the kind of binary systems which areexpected to be common among Trojans. A number of suchbinaries, comprising two components of almost equal size,have been found amongst the large trans-Neptunian objects(e.g., Noll et al. 2020). This comparison is of particular in-terest, given that the Patroclus system was, in all likelihood,implanted to the Trojan region from the trans-Neptunian re-gion source zone (e.g., Nesvorn´y et al. 2018). It seems likelythat the Patroclus system represents the closest example ofan Edgeworth-Kuiper belt binary system. Further informa-tion on the Patroclus system will become available in thein coming decades, as the binary is a target for flyby in2033 by the Lucy spacecraft (e.g., Levison et al. 2017). Sim-ilar smaller-scale systems may well exist among the Trojanpopulation , but their abundance is uncertain. Observation-ally, such small-scale binaries remain beyond our detection,and theoretical models of their survival depend on a num-ber of unknown parameters (e.g., Nesvorn´y et al. 2018, 2020;Nesvorn´y & Vokrouhlick´y 2019). The existence of Trojan bi-naries is interesting by itself, but in the context of our work,it is worth noting that, if such binaries exist, they likely serveas a feeding cradle for a population of Trojan pairs.Following this logic, then if the population of pairsamong the Trojans can become known and well character-ized, such that their dominant formation process is under-stood, that would in turn prove to be a source of new in-formation about Trojan binaries. Milani (1993) in his pi-oneering work on Jovian Trojan orbital architecture noteda case of L4-swarm objects (1583) Antilochus and (3801)Thrasymedes. Their suspicious orbital proximity led the au-thor to suggest that they may constitute a genetically re-lated couple of bodies. A viable formation process would bethrough the instability and dissociation of a former binary(Milani and Farinella, personal communication). Unfortu-nately, the Antilochus–Thrasymedes interesting configura-tion has not since been revisited, nor further studied in amore detail.This background information motivates us to conducta search for Jovian Trojan pairs. Unfortunately, even nowthe problem is not simple, and we consider our work to bean initial attempt, rather than providing a definitive solu-tion. In section 2, we explain our strategy, and describe thedifficulties in Trojan pair identification. This strategy led uspreliminarily identify the Jovian Trojans (258656) 2002 ES and 2013 CC as a potential pair. To test this hypothesis,we attempted to prove that these two bodies could be geneti- MNRAS000 , 1–22 (2020) pair of Jovian Trojans cally related using backward orbital integration, as describedin section 3. In section 4, we discuss potential formation pro-cesses for the pair, before presenting our concluding remarksand a call for observations in section 5. The Appendix A de-scribes our methods for the construction of Jovian Trojanproper elements. An up-to-date catalogue of those elements,which we have made publicly available online, is actually afruitful by-product of our work that may prove useful forfuture studies. We discuss some additional candidate pairsin Appendix B. The discovery of asteroid pairs was a direct by-product of a search for very young asteroid families(see Nesvorn´y et al. 2006; Nesvorn´y & Vokrouhlick´y 2006;Vokrouhlick´y & Nesvorn´y 2008). As a result, the primaryambition was to find pairs that formed recently, within thelast Myr, amongst the Main belt and Hungaria populations.In fact, the necessity for proven pairs to be young is essen-tially related to the method that allows their identification.Just like collisional families, asteroid pairs are identifiedas a result of the similarity of their heliocentric orbits. Thesearch for classical collisional families has traditionally beenperformed using clustering techniques in proper orbital ele-ment space, examining the proper semi-major axis a P , eccen-tricity e P and the sine of proper inclination sin I P (see, e.g.,Benjoya & Zappal`a 2002; Nesvorn´y et al. 2015, for reviews).The use of the proper elements allows us, with some care, tosearch for both young and old families. This is because theproper elements are believed to be stable over much longertimescales than other types of orbital elements, such as os-culating or mean, ideally on a timescale reaching hundredsof Myrs or Gyrs.There are, however, limitations to this method. In thecase of very old families, problems arise from instability ofthe proper orbital elements and the incompleteness of thedynamical model used to derive the proper elements. A dif-ferent problem occurs for very young families. The issue hasto do with the huge increase in the number of small-body ob-jects discovered over the past decades. Despite the fact thatthe very young families and asteroid pairs must have veryclose values of the proper orbital elements, it is difficult tostatistically discern them from random fluctuations of back-ground asteroids. Both occur at the same orbital distance inproper element space.This fundamental obstacle arises due to the low di-mensionality of proper element space, which consists ofjust three independent variables. In order to separate veryyoung asteroid families and asteroid pairs from the randomfluctuations of the background population, Nesvorn´y et al.(2006) and Vokrouhlick´y & Nesvorn´y (2008) realized thatthis problem can be overcome if the search is conductedin a higher-dimensional space. As a result, they used thefive-dimensional space of the osculating orbital elements, ne-glecting just the mean longitude. The mean orbital elementsare also suitable alternative parameters for such an analysis(e.g., Ro˙zek et al. 2011). In order to effectively use the twoextra dimensions, the searched structures must also be clus-tered in secular angles, the longitudes of ascending node andperihelion. This is perfectly justified for very young families and pairs that are expected to have separated at very lowvelocities.Previous searches for these young structures in thespace of osculating or mean orbital elements proved the use-fulness of the method, provided the age of the pair was lessthan about one Myr. Asteroid pairs will clearly exist thatformed earlier than this limit, but a differential precessionof their secular angles will result in them becoming effec-tively randomized, which will, in turn, render the identifi-cation procedure described above ineffective. A key pointhere is that the population of Main belt asteroids is cur-rently known to very small sizes, with objects detected withdiameters of one kilometer, or even smaller. The proposedformation processes for very young families and pairs areexpected to generate enough pairs within the last Myr that,even after accounting for discovery biases, we still have someof them in our catalogs.The situation is, however, different in the case of theJovian Trojan swarms. First, the characteristic size of thesmallest Trojans is ≃ km, with few objects being discov-ered that are smaller than this limit. Second, the formationprocesses of putative Trojan pairs, such as a rotational fis-sion or collisions, are significantly less efficient than in themain belt. As a result, no identifiable pairs among Trojansare expected to have been formed in the last − Myr,over which time, one would expect secular angles of any suchpairs to diverge from each other. We conducted a traditionalsearch for pairs in the five dimensional space of osculatingorbital elements (as in Vokrouhlick´y & Nesvorn´y 2008), butdid not find any candidates. If pairs do exist amongst theknown Trojans, their ages must be larger. In that case, how-ever, their secular angles would be randomized, as is thecase for old pairs in the main belt. Our candidate selec-tion method then returns back to the analysis of the Trojanproper elements, with further considerations based on addi-tional criteria.
The
AstDyS website, founded at the University ofPisa, and currently run by SpaceDys company (see https://newton.spacedys.com/astdys/ ), is a worldrenowned storehouse of proper orbital elements for Solarsystem minor bodies. It also contains data on the JovianTrojans, namely synthetic proper elements based on mathe-matical methods presented in the pioneering work of Milani(1993). We also note the work of Beaug´e & Roig (2001),which discusses an alternative approach to the calculationof Trojan proper elements, but these authors neither maketheir results readily available online, nor update them ona regular basis. For that reason, one possibility for thisstudy would be to use the
AstDyS data. However, thosedata have at least two drawbacks for our application. First,their last update occurred in June 2017. As a result, theyprovide information for a total of 5553 numbered andmulti-opposition Jovian Trojans. Given the efficiency ofall-sky surveys, this number has increased significantly inthe years since that update, with more than 7000 JovianTrojans now known for which observations span multipleoppositions. Second, the proper elements provided at
AstDyS are given to a precision of just four decimal places, which is not sufficient for our work. The
AstDyS database
MNRAS , 1–22 (2020)
T. R. Holt et al. (a) (b)
Figure 1.
Proper orbital elements of Jovian Trojans: semimajor axis da P vs sine of inclination sin I P (top panels), and semimajor axis da P vs eccentricity e P (bottom panels). the left panels show the orbits of 4607 objects at the L4 libration zone, whilst the right panels showthe orbits of 2721 objects at the L5 libration zone. These data were computed using the method described in the Appendix A and displaynumbered and multi-opposition orbits as of April 2020. Blue triangles indicate the largest objects of the previously identified Trojanfamilies (e.g., Rozehnal et al. 2016): (a) Eurybates, Hektor, Arkesilaos and 1996 RJ in the L4 zone, and (b) Ennomos and 2001 UV inthe L5 zone. The two green circles denote position of Jovian Trojans (1583) and (3801) that were previously identified as constituting asuspiciously close pair (see Milani 1993). The two overlapping red circles denote the location of our proposed pair candidate of (258656)2002 ES and 2013 CC . would, as a result, allow the determination of the orbitaldistance in the proper element space –Eq. (1)– with only ≃ to m s − accuracy, which is insufficient to characterizethe low velocity tail. For both of these reasons, in this work,we decided to determine our own synthetic proper elements.Details of the approach are given in Appendix A. Here,we only mention that our proper element definition andmathematical methods follow the work of Milani (1993),with substantial differences only for those orbits with verysmall libration amplitudes. Previous applications using thistechnique may be found in Broˇz & Rozehnal (2011) andRozehnal et al. (2016).Figure 1 shows our results, namely proper elementscomputed for 7328 Jovian Trojans (numbered and multi-opposition objects as of April 2020) projected onto the ( da P , sin I P ) and ( da P , e P ) planes for the L4 swarm (“Greeks”leading Jupiter on its orbit; left panels) and the L5 swarm(“Trojans” trailing behind Jupiter; right panels). The L4swarm is more numerous, partly as a result of four majorcollisional families that have been recognised in recent years(e.g., Rozehnal et al. 2016), and contains 4607 objects. The smaller L5 swarm contains only 2721 known objects, includ-ing the 2001 UV and Ennomos collisional families. Toproceed with an investigation of the orbital similarity be-tween members of the Trojan population, the basis of thepair and family recognition process, one must introduce ametric function in the space of the proper orbital elements.Several choices have been discussed by Milani (1993). Weopt for the d metric, also favoured by the author of thatwork, though we slightly adjust that metric , such that theorbital distance is given in velocity units. Given two orbits inthe Trojan L4 or L5 proper element space, obviously with-out mixing the two swarms, we define their distance δ V P asa quadratic form using the differences δ a P , δ e P and δ sin I P as δ V P = V J s (cid:18) δ da P a J (cid:19) + ( δ e P ) + ( δ sin I P ) , (1)where V J ≃ m s − and a J ≃ . au are mean orbitalvelocity and semimajor axis of Jupiter. Milani (1993) arguedthat this particular choice of the coefficients – ( . , , ) – MNRAS , 1–22 (2020) pair of Jovian Trojans helps to equally weight contributions from all three dimen-sions. Given the metric shown in Eq. (1), we computed distances ofall possible pairs in the L4 and L5 Trojans swarms, and orga-nized them in the form of a cumulative distribution N ( < δ V P ) (see also Vokrouhlick´y & Nesvorn´y 2008, for context). Theresults of this process are shown in Fig. 2. Whilst the largest δ V P values of approximately V J are set by the maximum ex-tension of the stable phase space of tadpole orbits associ-ated with Jupiter (Fig. 1), the smallest δ V P values of order ∼ − m s − are determined by a combination of several fac-tors. The number of known Jovian Trojans filling the stableorbital space is the first factor, compared with the typicalsmallest values δ V P ≃ m s − found by Milani (1993), whostudied just 80 and 94 Trojans in the L4 and L5 swarms,respectively. Additionally, small velocity differences occurwhen bodies become organized in structures like families.Last, the inevitable uncertainty of the proper elements con-tributes to the noise in δ V P . We determine the uncertaintiesof δ V P by a propagation of the proper element uncertain-ties described in the Appendix A. This effect is obviouslynot uniform, but organized in a complicated structure of achaotic web, generally increasing toward the border of thestable tadpole zone (see, e.g., Robutel & Gabern 2006). In-terestingly, the characteristic noise level from such determin-istic chaos is of the order of a few meters per second, aboutthe same as minimum distances between the orbits, as canbe seen in Fig. 2, where we show uncertainty intervals of δ V P for the low-velocity tail.It is also worth noting that for reasonably small val-ues of δ V P (hundred m s − or so), one would expect N ( <δ V P ) ∝ ( δ V P ) provided that: (i) Trojans fill the availablestable phase space at random, and (ii) the weighting coeffi-cients in the metric function (1) truly express isotropy, theexponent 3 is then a measure of the proper element spacedimension. For large δ V P values the cumulative distributions N ( < δ V P ) become shallower because of the finite extent ofthe stable orbital region. We also note that N ( < δ V P ) holdsglobal information about the whole L4 and L5 populations,while local structures, such as families and clusters, are al-most not seen in this distribution.We find it interesting that N ( < δ V P ) are broadly similarfor the L4 and L5 swarms, but they also differ in some im-portant characteristics, in particular, the smallest and thelargest δ V P values. This is due to the directly comparablepopulations of the two swarms and basically identical vol-umes of their stable phase space. However, the δ V P < m s − parts of the distributions have a different behaviourwhen approximated with a power-law N ( < δ V P ) ∝ ( δ V P ) α :(i) the L4 swarm has the canonical value α ≃ , while (ii)the L5 swarm is shallower, with approximately α ≃ / . Wehypothesize that this difference is caused by a presence ofthe prominent Trojan families in the L4 population. Fam-ily members efficiently contribute to the low- δ V P part of thedistribution. Given their small extent, it is also conceivablethat the mutual orbital distribution in families is approx-imately isotropic. The L5 population is less influenced byTrojan families, and, as a result, N ( < δ V P ) may reflect theparameters of the background Trojan population. This is af- fected both by the resonances that sculpt the stable orbitalzone in a complicated way and, perhaps , the initial fillingof the Trojan region by planetesimals. Finally, the weight-ing coefficients of the metric function (1), that express howdifferences in semimajor axis, eccentricity and inclinationcontribute to the whole, may also slightly affect the result(though our experiments with small changes in those valuesdid not yield significant differences). If combined altogether,the α value may be slightly shallower than , such as / wefound for the L5 population.Seeking details that could explain the difference in thepopulation exponents α in further detail, we analysed distri-butions of the proper elements. The most significant differ-ence concerns proper inclination I P . Figure 4 shows L4 andL5 Trojan distributions of I P for all bodies. The dashed linesare simple approximations with a function I P exp (− I P / C ) ,where the adjustable constant C characterizes width of thedistribution (the prominent families, such as Eurybates at ≃ ◦ among L4 or Ennomos at ≃ ◦ among L5, were ex-cluded from the fit). We found C ≃ . ◦ for L4 and C ≃ . ◦ for L5, implying the inclination distribution at L5 is slightlybroader. This confirms results in Di Sisto et al. (2014) . Itis not clear, whether this is due to the details of the captureprocess, or whether the escapees from the prominent Eu-rybates and Arkesilaos families in the L4 swarm contributeto the difference, and how it may affect the exponent α ofthe N ( < δ V P ) distribution discussed above. A full analysis ofthese interesting findings is beyond the aims of our work.Regarding the smallest δ V P values, neither of the two dis-tribution functions N ( < δ V P ) show a change in behaviour.In the context of our work, this implies no hint of a statis-tically significant population of very close orbits, a tracerof a possible Trojan pair population. In fact, given the lowdimensionality of the proper element space, this was not un-expected , given that the asteroid pairs in the Main beltwould not manifest themselves using a similar analysis. Theslight deviation of N ( < δ V P ) below ≃ m s − velocity to ashallower trend for the L5 swarm is interesting, but likelynot statistically robust enough to allow firm conclusions tobe drawn at the current time.We paid some attention to the smallest-distance couple(215110) 1997 NO –2011 PU , and could not conclusivelyprove that it represents a real pair of related objects (Ap-pendix B). A closer analysis of the second to sixth closestcouples in the L5 population indicates the possibility of avery compact cluster about Trojan (381148) 2007 GZ , butits status needs to be confirmed with more data in the future.In any case, because our interest here focuses on Trojans inthe low-velocity tail of the N ( < δ V P ) distribution, seeking pu-tative pairs, we also show in Fig. 3 location of couples thathave δ V P < m s − in both Trojan swarms. These wouldbe the most logical candidates for further inspection.A full frontal approach to this data would be to analysethe results from backward orbital integrations for these littlemore than 200 putative couples using the methods describedin section 3. However, this would require a significant com-putational effort, and thus we chose to adopt further cri-teria for candidate selection. For instance, data in the L4swarm show that the lowest δ V P couples are strongly con-centrated in the recognized families. The locally increaseddensity of Trojans in these regions obviously imply small dis-tances δ V P , but this also means such couples are most likely MNRAS , 1–22 (2020)
T. R. Holt et al. (a) (b)
Figure 2.
Cumulative distribution N ( < δ V P ) of Trojans with velocity distance δ V P in the proper elements space using the metricdescribed in Eq. (1): the left panel presents the results for the 4607 objects of the L4 Trojan swarm, with the right-hand panel showingthe 2721 members of the L5 Trojan swarm. The light gray solid lines indicate the N ( < δ V P ) ∝ ( δ V P ) relationship, for reference; curiously,the L5 distribution is better matched with N ( < δ V P ) ∝ ( δ V P ) / , shown with a dashed gray line. Blue symbols denote the population withthe smallest δ V P values, namely δ V P ≤ m s − for both the L4 and L5 Trojans. For sake of interest, we also show uncertainties in thedetermination of δ V P for these low-velocity couples with grey horizontal intervals. Position of three couples of interest is highlighted bylabels. These are the (1583)-(3801) couple with δ V P = ( . ± . ) m s − and (258656) 2002 ES -2013 CC couple δ V P = ( . ± . ) m s − among L4 Trojans, and (215110) 1997 NO -2011 PU couple with δ V P = ( . ± . ) m s − among L5 Trojans. not the objects that we seek. The correlation with Trojanfamilies is somewhat weaker in the L5 swarm, though severalof the small-distance couples are found in both the Ennomosand 2001 UV families. Other constitute compact clustersscattered in the background population, like that around(381148) 2007 GZ , as mentioned above.Sifting the δ V P < m s − couples unrelated to fami-lies would still leave us with too many candidates to pursuewith backward n -body simulations. Having experimentedwith several cases, we adopted the strategy of focusing onthose low- δ V P couples characterized by (i) the least popu-lated background, and (ii) located in the most dynamicallystable zones of the orbital phase space. The former conditionincreases the likelihood that the candidate couple is a realpair, and not just a fluke, whilst the latter condition wouldallow us to investigate the past orbital configuration of theputative pair across as lengthy a timescale as possible. Thisis particularly important for pairs in the Jovian Trojan pop-ulation, since no recently-formed pairs are to be expected, asdescribed above. Moreover, the expected large ages of pos-sible Trojan pairs do not allow us to seek their past orbitalconvergence in full six-dimensional Cartesian space of posi-tions and velocities. Even the most stable Trojan orbits havean estimated Lyapunov timescale of about − Myr. Inthis situation, our convergence scheme should rely on the be-haviour of secular angles, the longitudes of node and perihe-lion, and the related eccentricity and inclination (section 3).It is then advantageous to suppress the role of the last twoelements, the semimajor axis and the mean longitude , byletting them vary as little as possible . This favours loca-tions very near the tadpole libration center of either the L4or L5 swarms, where also the previous two conditions, low background population and maximum orbital stability, aresatisfied.
With all these criteria in mind, we found a candidate cou-ple of L4 objects, (258656) 2002 ES and 2013 CC .The proximity of these two objects to the libration cen-ter is reflected by the small values of all proper elements(see Fig. 1), namely da P ≃ ( . ± . ) × − au, e P ≃ ( . ± . )× − and sin I P ≃ ( . ± . )× − for (258656) 2002 ES , and da P ≃ ( . ± . )× − au, e P ≃ ( . ± . )× − and sin I P ≃ ( . ± . )× − for 2013 CC . The close proximity to L4 also indicates thatthe pair have been in stable orbits for the life of the Solarsystem (e.g., Holt et al. 2020). For reference, we also men-tion their libration amplitude, in the angular measure, whichis only about . ◦ , resp. . ◦ . There are only four otherL4 objects in our sample that have smaller libration ampli-tudes, and none among the known L5 objects, though thesehave generally larger proper eccentricity and/or inclinationvalues. The similarity of the two orbits is immediately ap-parent and quantitatively expressed with δ a P ≃ . × − auand δ e P ≃ . × − , both with negligible uncertainty, while δ sin I P ≃ . × − with a small uncertainty of . × − .This uncertainty amounts to about . ◦ difference in theproper inclination. All these values result in the velocity dif-ference δ V P ≃ . ± . m s − , using our adopted metric(1), dominated by the inclination contribution the contribu-tion from the difference in proper eccentricities is about %of the total, and the difference in proper semimajor axes is MNRAS000
With all these criteria in mind, we found a candidate cou-ple of L4 objects, (258656) 2002 ES and 2013 CC .The proximity of these two objects to the libration cen-ter is reflected by the small values of all proper elements(see Fig. 1), namely da P ≃ ( . ± . ) × − au, e P ≃ ( . ± . )× − and sin I P ≃ ( . ± . )× − for (258656) 2002 ES , and da P ≃ ( . ± . )× − au, e P ≃ ( . ± . )× − and sin I P ≃ ( . ± . )× − for 2013 CC . The close proximity to L4 also indicates thatthe pair have been in stable orbits for the life of the Solarsystem (e.g., Holt et al. 2020). For reference, we also men-tion their libration amplitude, in the angular measure, whichis only about . ◦ , resp. . ◦ . There are only four otherL4 objects in our sample that have smaller libration ampli-tudes, and none among the known L5 objects, though thesehave generally larger proper eccentricity and/or inclinationvalues. The similarity of the two orbits is immediately ap-parent and quantitatively expressed with δ a P ≃ . × − auand δ e P ≃ . × − , both with negligible uncertainty, while δ sin I P ≃ . × − with a small uncertainty of . × − .This uncertainty amounts to about . ◦ difference in theproper inclination. All these values result in the velocity dif-ference δ V P ≃ . ± . m s − , using our adopted metric(1), dominated by the inclination contribution the contribu-tion from the difference in proper eccentricities is about %of the total, and the difference in proper semimajor axes is MNRAS000 , 1–22 (2020) pair of Jovian Trojans (a) (b) Figure 3.
Proper orbital elements of Jovian Trojans as in Fig. 1. The red symbols in both L4 and L5 swarms show couples of Trojanswith δ V P < m s − , namely the lowest velocity tail in the distributions shown in Fig. 2: the primary component of each couple is shownusing a filled red circle, whilst the secondary is shown using a blue circle. In the L4 case their relation to the recognized families isapparent. In the L5 case their distribution is more scattered, though some are also associated with the 2001 UV and Ennomos families. Figure 4.
Number of Jovian Trojans with proper inclination I P (in degrees), showing the L4 (red) and L5 (blue) swarms. Thedashed lines represent an approximation I P exp (− I P / C ) for thebackground population (significant peaks due to Trojan familieseliminated), where we found C ≃ . ◦ for L4 and C ≃ . ◦ for L5. negligible). With that said, this couple would qualify among the closest in the population if it were not for the slightinclination offset of the two orbits.Not much physical information is available about thesetwo objects. Various databases providing orbital solutions(such as AstDyS , JPL or MPC ) yield an absolute magnitude for(258656) 2002 ES in the range . to . , and values inthe range . to . for 2013 CC . Given the mean albedo, p V ≃ . , for small Trojans (a value with an admittedlylarge scatter; e.g., Grav et al. 2011, 2012), we estimate theirsizes to be D ≃ . − . km for (258656) 2002 ES and D ≃ . − . km for 2013 CC . Unless the assumptionof similar albedoes is significantly in error, it is clear thatthe two bodies are similar in size, though not exactly thesame. No other physical parameters, such as the rotationperiod, thermal inertia and/or spectral colors, are known atthe present time. Further observational follow-up on theseobjects is therefore highly recommended. The small libration amplitude zone of the proper elementspace contains a relatively small number of bodies, as canbe seen in the left panel (a) in Fig. 5. Here, we used therange da P ≤ . au, expressing the proximity to the li- MNRAS , 1–22 (2020)
T. R. Holt et al. (a)
0 0.005 0.01 0 0.05 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 d a p [au] e p sin I p (b)
0 0.0005 0.001 0.0015 0.002 0 0.05 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 S [au ] e p sin I p Figure 5.
Left (a): The small-libration portion of the L4 stable orbital zone in the 3-dimensional proper element coordinates ( da P , e P , sin I P ) . The proximity to the libration was arbitrarily set by da P ≤ . au, whilst the extent of e P and sin I P is limited byorbital stability. We find objects (black symbols) in this zone for our dataset of Trojans. The candidate pair (258656) 2002 ES and2013 CC is highlighted with a red circle. The vertical intervals help to appreciate 3-dimensional nature of the display. Right (b): Thesame as on the left panel, but the da P was replaced with a surface area S = π ( da P ) . In this case, the small-amplitude Trojans aredistributed more uniformly. bration center, but left e P ≤ . and sin I P ≤ . , gener-ally capturing the width of the stable Trojan phase space(Levison et al. 1997; Nesvorn´y et al. 2002a; Tsiganis et al.2005; Di Sisto et al. 2014; Holt et al. 2020). We could havealso more strongly restricted the proper eccentricity andinclination values , but if this is done too aggressively, itwould result in the sample of observed Trojans available forour analysis becoming too small. With our limits, we find k = Trojans in the L4 space, including our candidate pair(258656) 2002 ES and 2013 CC .The proper element differences in the (258656)2002 ES and 2013 CC couple are δ a P = . × − au , δ e P = . , δ sin I P = . , much smaller than thescale of the chosen zone, assuming that all dimensions aretaken equally. In the first approximation, taking all di-mensions equally, and thus neglecting the weighting coef-ficients from Eq. 1 which are all of the order of unity, the ( δ a P , δ e P , δ sin I P ) differences in this couple define a small boxof which represents only a ≃ . × − fraction of the anal-ysed target zone. For statistical calculations , it is useful toimagine “numbered” boxes of the ( δ a P , δ e P , δ sin I P ) volumein the whole zone. Their total number of such boxes wouldthen be n ≃ . × .The simplest estimate of the statistical significant of the(258656) 2002 ES -2013 CC pair is based on the assump-tion that bodies were distributed in the analysed zone ran-domly/uniformly. We choose k numbers from n possibilities(i.e., one for each body from a set of “numbered” boxes). Or-dered, repeated selections are given as variations V ′ ( n , k ) = n k , while ordered, non-repeated as V ( n , k ) = n ! /( n − k ) ! . Thelikelihood that among the trials the box-numbers do not re-peat is simply the ratio V ( n , k )/ V ′ ( n , k ) , and we are interestedjust in the complementary probability: p = − V ( n , k ) V ′ ( n , k ) ≃ . × − . (2)We verified this result by directly running a Monte-Carlosimulation of the selection process. Thus, we find the prob-ability that the selected couple is only a random orbital co- incidence to be very low. Shrinking the width of the e P and sin I P to half the previously mentioned values did not changeour result significantly.As can be seen in the left panel (a) of Fig. 5, the as-sumption of a uniform distribution of background Trojansin the target zone is fair, but not exactly satisfied. This isthe result of the decreasing number of Trojans towards thelibration center (i.e., at very small values da P ). We there-fore repeated our analysis in a different system of coor-dinates. Keeping e P and sin I P , we now changed da P with S = π ( da P ) . The background reasoning is that the libra-tion point, da P = , represents a center about which thetadpole orbits move in 3-dimensions. In a Cartesian viewcentered at L4 the radial coordinate is to be replaced withthe surface area S = π ( da P ) . Re-mapping and re-binningour analysis in the ( S , e P , sin I P ) coordinate system, we ob-tained the situation shown in the right panel (b) of Fig. 5.Whilst still keeping the same number k = of Trojans inthe analysed zone, their distribution is now more uniform.Given the new box-definition by the (258656) 2002 ES and2013 CC couple, we now find the number of thus definedsmall boxes to be increased to n ≃ . × This is theresult of the candidate couple’s close proximity to the libra-tion center. As a result, the likelihood (Eq. 2) of the couplebeing just a fluke in a uniform distribution of objects nowbecomes smaller, namely p ≃ . × − .The probability p , defined and computed for the(258656) 2002 ES –2013 CC couple above, is appreciablysmall. It is both interesting and important to compare thisresult with the similarly defined quantity for other Trojancouples, especially amongst those that have a small δ V P dis-tance in the metrics (1). This will tell us whether the prob-ability p for (258656) 2002 ES –2013 CC is sufficientlysmall in absolute measure for the couple to be considereda true pair, whilst at the same time enabling our algorithmto better connect our p definition with the velocity metricsused above. Here we analyze the L4-swarm population , butthe same approach could equally be applied to the L5 case. MNRAS000
Left (a): The small-libration portion of the L4 stable orbital zone in the 3-dimensional proper element coordinates ( da P , e P , sin I P ) . The proximity to the libration was arbitrarily set by da P ≤ . au, whilst the extent of e P and sin I P is limited byorbital stability. We find objects (black symbols) in this zone for our dataset of Trojans. The candidate pair (258656) 2002 ES and2013 CC is highlighted with a red circle. The vertical intervals help to appreciate 3-dimensional nature of the display. Right (b): Thesame as on the left panel, but the da P was replaced with a surface area S = π ( da P ) . In this case, the small-amplitude Trojans aredistributed more uniformly. bration center, but left e P ≤ . and sin I P ≤ . , gener-ally capturing the width of the stable Trojan phase space(Levison et al. 1997; Nesvorn´y et al. 2002a; Tsiganis et al.2005; Di Sisto et al. 2014; Holt et al. 2020). We could havealso more strongly restricted the proper eccentricity andinclination values , but if this is done too aggressively, itwould result in the sample of observed Trojans available forour analysis becoming too small. With our limits, we find k = Trojans in the L4 space, including our candidate pair(258656) 2002 ES and 2013 CC .The proper element differences in the (258656)2002 ES and 2013 CC couple are δ a P = . × − au , δ e P = . , δ sin I P = . , much smaller than thescale of the chosen zone, assuming that all dimensions aretaken equally. In the first approximation, taking all di-mensions equally, and thus neglecting the weighting coef-ficients from Eq. 1 which are all of the order of unity, the ( δ a P , δ e P , δ sin I P ) differences in this couple define a small boxof which represents only a ≃ . × − fraction of the anal-ysed target zone. For statistical calculations , it is useful toimagine “numbered” boxes of the ( δ a P , δ e P , δ sin I P ) volumein the whole zone. Their total number of such boxes wouldthen be n ≃ . × .The simplest estimate of the statistical significant of the(258656) 2002 ES -2013 CC pair is based on the assump-tion that bodies were distributed in the analysed zone ran-domly/uniformly. We choose k numbers from n possibilities(i.e., one for each body from a set of “numbered” boxes). Or-dered, repeated selections are given as variations V ′ ( n , k ) = n k , while ordered, non-repeated as V ( n , k ) = n ! /( n − k ) ! . Thelikelihood that among the trials the box-numbers do not re-peat is simply the ratio V ( n , k )/ V ′ ( n , k ) , and we are interestedjust in the complementary probability: p = − V ( n , k ) V ′ ( n , k ) ≃ . × − . (2)We verified this result by directly running a Monte-Carlosimulation of the selection process. Thus, we find the prob-ability that the selected couple is only a random orbital co- incidence to be very low. Shrinking the width of the e P and sin I P to half the previously mentioned values did not changeour result significantly.As can be seen in the left panel (a) of Fig. 5, the as-sumption of a uniform distribution of background Trojansin the target zone is fair, but not exactly satisfied. This isthe result of the decreasing number of Trojans towards thelibration center (i.e., at very small values da P ). We there-fore repeated our analysis in a different system of coor-dinates. Keeping e P and sin I P , we now changed da P with S = π ( da P ) . The background reasoning is that the libra-tion point, da P = , represents a center about which thetadpole orbits move in 3-dimensions. In a Cartesian viewcentered at L4 the radial coordinate is to be replaced withthe surface area S = π ( da P ) . Re-mapping and re-binningour analysis in the ( S , e P , sin I P ) coordinate system, we ob-tained the situation shown in the right panel (b) of Fig. 5.Whilst still keeping the same number k = of Trojans inthe analysed zone, their distribution is now more uniform.Given the new box-definition by the (258656) 2002 ES and2013 CC couple, we now find the number of thus definedsmall boxes to be increased to n ≃ . × This is theresult of the candidate couple’s close proximity to the libra-tion center. As a result, the likelihood (Eq. 2) of the couplebeing just a fluke in a uniform distribution of objects nowbecomes smaller, namely p ≃ . × − .The probability p , defined and computed for the(258656) 2002 ES –2013 CC couple above, is appreciablysmall. It is both interesting and important to compare thisresult with the similarly defined quantity for other Trojancouples, especially amongst those that have a small δ V P dis-tance in the metrics (1). This will tell us whether the prob-ability p for (258656) 2002 ES –2013 CC is sufficientlysmall in absolute measure for the couple to be considereda true pair, whilst at the same time enabling our algorithmto better connect our p definition with the velocity metricsused above. Here we analyze the L4-swarm population , butthe same approach could equally be applied to the L5 case. MNRAS000 , 1–22 (2020) pair of Jovian Trojans -6 -5 -4 -3 -2 EurybatesHektor(9799) p [] δ V P [m/s] Figure 6.
The probability p that a pair is random fluke, com-puted using the method described in the text, versus its dis-tance δ V P , computed for all low-velocity pairs in the L4 zone usingEqs. (2) and (1). The pair (258656) 2002 ES and 2013 CC ishighlighted with a red circle. The pair (219902) 2002 EG and(432271) 2009 SH , discussed briefly in the Appendix B, is high-lighted with a green circle. The colored symbols denote pairs inthe identified L4 families: (i) Eurybates (blue), (ii) the core ofthe Hektor family (light blue), and (iii) (9799) 1996 RJ (cyan).The dashed line, p = − ( δ V P /(
10 m s − )) , is used to emphasizethat the candidate pair (258656) 2002 ES and 2013 CC is anoutlier in this population. The potentially complicated part of the procedureis that, for each selected couple, we have to (i) adaptthe box size ( δ a P , δ e P , δ sin I P ) , and (ii) the zone size ( ∆ a P , ∆ e P , ∆ sin I P ) , as well as the position to which the boxsize refers. The choice of the latter obviously varies becausethe local number density of bodies differs from place toplace. In order to prevent excessively small boxes in oneof the dimensions (as an example, due to an almost zerodifference δ e P (cid:17) ), we use the metric δ V P as a measure ofthe “diagonal” of the box and we define its respective vol-ume as ( δ V P ) /√ . Observing the typical spatial variationof the number density of Trojans, we use a fixed value for ∆ a P = . au, rejecting pairs with δ a P > . ∆ a P . In order toprevent a low number of bodies k in the zone, both ∆ e P and ∆ sin I P are then sequentially increased until k ≥ . Once weset the zone, we again define its volume as ( ∆ V P ) /√ , with ∆ V P the velocity distance of the corners connected with a di-agonal. The number of boxes n , as well as the probability p ,is then computed as before (Eq. (2)). Obviously, the wholealgorithm cannot be done manually, but an automated com-puter script was written to run the method.The statistical results of our analysis are shown inFig. 6. The pairs seem to be well organized in the ( p , δ V P ) plane, expressing an overall correlation between the twoquantities. As might be expected, the general trend is p ( δ V P ) ∝ ( δ V P ) , namely volume of the box. Nevertheless,the p vs δ V P values do not follow a single curve, due to thelocal number density being different for each of the couples.Those couples located within known families generally haverelatively high p values. This is to be expected, since thesurrounding zones are densely populated by Trojans, which causes the dimensions of the zone to be small. To illustratethis effect, we colored data for pairs in the largest familiesin the Fig. 6, identifying those in the (i) Eurybates family(blue), (ii) the core of the Hektor family (light blue), and(iii) the (9799) 1996 RJ family (cyan), after Nesvorn´y et al.(2015). The Eurybates family, the largest and most popu-lous in the Trojan population, has systematically the largest p values. This is because even a small zone quickly containsour threshold number of k = Trojans . We note that p ≃ ,or even formally larger, just indicates that a couple of Tro-jans in this zone is fully expected at their distance δ V P . Anexception to this rule is the (9799) 1996 RJ family, where wefind the smallest p values, which are clearly correlated with δ V P . This is because (9799) 1996 RJ is a very compact fam-ily located in isolation in a high-inclination portion of theTrojan phase space (see also Fig. 1). For each of the couplesselected in this family, the reference zone needs to be largeto contain the minimum required number of objects.Whilst the collisional families could clearly contain dy-namical pairs, their recognition is confused by the locallyhigh background of family members. We therefore excludeobjects located in families from our work. What remains isthen a diffuse background population of Trojans. For ev-ery fixed δ V P value, there are some background couples forwhich p extends to small values. The true Trojan dynamicalpairs, namely those objects genetically related to a commonparent, form the basis for our search among this populationof a low- p tail for sufficiently small δ V P values. There arepossibly a number of such cases, but amongst them, the onewhich is the most outlying from the p ( δ V P ) ∝ ( δ V P ) refer-ence level shown by the dashed curve in Fig. 6 is the case of(258656) 2002 ES –2013 CC (highlighted with red circle).Its p value is an order of magnitude lower when compared tocouples with similar δ V P values. This justifies the validity ofthe (258656) 2002 ES –2013 CC couple as a true asteroidpair, based on our statistical analysis alone. There are alsosome family-unrelated couples with p values comparable orsmaller, and these are briefly discussed in Appendix B.In the next section 3, we conduct a search for past or-bital convergence of the selected (258656) 2002 ES and2013 CC couple. If successful, this process add an impor-tant piece of evidence justifying the couple as a real pair ofgenetically related objects. We explain our methods in de-tail. These methods are also briefly applied to several othercandidate couples, with less success (Appendix B). The dynamics of the Jovian Trojans have been exten-sively studied using both analytical and numerical means(e.g., Milani 1993; Beaug´e & Roig 2001; Robutel & Gabern2006; Di Sisto et al. 2014; Holt et al. 2020, and referencestherein). Here, we confine ourselves to briefly recalling onlythe information necessary for understanding and interpret-ing our numerical simulations of the (258656) 2002 ES –2013 CC pair.As previously noted, the objects in this pair are nottypical, but are instead exceptional representatives of Tro-jan population. This is because they reside extremely closeto the L libration center. As a result, the evolution of theirsemimajor axis a and the resonant argument λ − λ J be charac- MNRAS , 1–22 (2020) T. R. Holt et al. terized by many small-amplitude and high-frequency terms.Those are, however, of the least importance for our analysis.More relevant is the behaviour of the eccentricity e , the in-clination I , the longitude of ascending node Ω , and the longi-tude of perihelion ̟ . Due to the small values of the eccentric-ity and inclination , it is also useful to think about complexnon-singular elements z = e exp ( ı̟ ) and ζ = sin I exp ( ı Ω ) .In linear perturbation theory, a fairly satisfactory zero ap-proximation, both z and ζ are represented by a finite num-ber of Fourier terms, namely the proper term and a fewforced planetary terms. A simpler description concerns ζ ,whose Fourier representation is dominated by the properterm with I P ≃ . ◦ , followed only by small contributionsfrom the s term, with I ≃ . ◦ , and a number of signifi-cantly smaller contributions. As a result, the osculating in-clination I is well represented by a constant I P and a periodicterm with amplitude I . Correspondingly, the osculating lon-gitude of the ascending node, Ω , steadily circulates with aperiod given by the proper s frequency, and experiences onlyvery small perturbation from the s term. The evolution of z is more complicated because it is represented by three termsof comparable amplitude. The largest-amplitude contribu-tion, ≃ . , is provided by the term with frequency g ,followed by proper g and g terms with comparable ampli-tudes of ≃ . and ≃ . . Whilst still very simple inthe Cartesian representation of z , the polar variables in thisplane (i.e. the eccentricity and especially longitude of peri-helion) exhibit a non-linear evolution, characteristic of manylow-eccentricity asteroid orbits. Equipped with this knowledge, we can now turn to in-vestigating the common origin of (258656) 2002 ES and 2013 CC . In studies of asteroid pairs, researchersseek to demonstrate a convergence of heliocentric orbitsof the proposed pair at some moment in the past (e.g.,Vokrouhlick´y & Nesvorn´y 2008). This is considered to bethe origin of the two objects from a common parent body,and the corresponding time in the past representative of theage of the pair. As typically achievable ages of the asteroidpairs in the Main belt are less then one Myr, with manyless than kyr, a convergence is often sought in Cartesianspace. This approach means to demonstrate that the two or-bits meet at the same point in space and have a very smallrelative velocity.The same condition can be expressed in heliocentric or-bital elements by making them basically equal at the forma-tion moment of the pair. For this work, we find it markedlymore useful to work with the orbital elements of our candi-date pair, as they can teach us more readily about the evo-lution of the orbits of the two objects. Therefore, in Fig. 7,we show the results of our initial numerical experiment. Weprovide the differences between the osculating heliocentricelements of the nominal orbits of (258656) 2002 ES and2013 CC over a short time interval of the past Myr. Weuse the swift_rmvs4 integrator (Levison & Duncan 1994)which allows us to efficiently include gravitational pertur-bations from all eight planets. The integration timstep usedwas days, and the state vectors of all propagated bodies,planets and the two Trojans, were output every years.We use a reference system defined by the invariable plane Figure 7.
Differences between the osculating orbital elementsof (258656) 2002 ES and 2013 CC from a Myr backwardintegration of their nominal orbits. Gravitational perturbationsfrom all planets were included and an invariable-plane referencesystem used. The differences of semimajor axis δ a , eccentricity δ e , longitude of pericenter δ̟ and longitude in orbit δλ (topfour panels) indicate a simultaneous collapse to near zero valuesat ≃ . Ma (gray vertical line). In contrast, the differences ofinclination δ I and longitude of ascending node δ Ω (the bottomtwo panels) do not converge at that epoch: the nodal longitudesof the two objects are still ≃ ◦ away from each other, and the in-clination difference shows steady oscillation about the mean valueof ≃ − . ◦ , namely a difference in their proper inclinations. Thesteady trend in δ Ω has a slope ≃ . arcsec yr − , very close tothe difference in proper frequencies s of (258656) 2002 ES and2013 CC . MNRAS , 1–22 (2020) pair of Jovian Trojans of the planetary system. The initial conditions of (258656)2002 ES and 2013 CC at MJD58800 epoch were obtainedfrom the AstDyS website.The differences in the orbital elements shown in Fig. 7oscillate with the dominant frequencies identified by theanalysis of z and ζ themselves. For instance, the principalperiodicity seen in δ I and δ Ω corresponds to the frequency s − s , whilst the principal periodicity seen in δ e and δ̟ corresponds to frequencies g and g − g . Differences δ a and δλ are characterized by higher frequencies, such as the plan-etary orbital frequencies, the libration frequency, and thenfollowed by a “forest” of lower frequencies starting with g .We also note a markedly different behavior of δ̟ and δ Ω , which can be understood from the above mentioned de-scription of the z and ζ non-singular elements of the two ob-jects. Observing the general behavior of the amplitude in the ( δ a , δ e , δ̟, δλ ) terms, we note a curious fact that those ampli-tudes become very small simultaneously for semimajor axis,eccentricity, longitude of perihelion and longitude in orbit ≃ . Myr ago (upper four panels in Fig. 7). However, anyhope for a clear orbital convergence at that epoch is removedby looking at behavior of the inclination and longitude of as-cending node differences (bottom two panels in Fig. 7). Wenote that δ I keeps steadily oscillating about a mean value of ≃ − . ◦ , namely a difference in the proper inclinations of(258656) 2002 ES and 2013 CC , without the amplitudeof those oscillations showing any tendency to shrink. At thesame time, the nodal difference stays large, and only slowlydecreases from ≃ − ◦ to ≃ − ◦ . This rate of decrease in δ Ω fits perfectly the difference in proper frequencies s of thetwo objects as to be expected. Hence some ≃ . Myr ago,the two orbits had basically identical ( a , e , ̟, λ ) values, butthe nodes were still offset by about ◦ . This is inconsis-tent with any believable low-velocity separation of the twoobjects from a common parent body at their origin. Whilstinconclusive about the origin of the (258656) 2002 ES and2013 CC couple, this Myr integration provides usefulhints for further analyses.
Extrapolating the trend seen in Fig. 7, we can estimate thatthe nodes of (258656) 2002 ES and 2013 CC became co-incident some Myr ago. Obviously, this is only the firstsuch configuration in the historical evolution of the two ob-jects. Assuming orbital stability, we also predict that theconfiguration will repeat with a ≃ Myr periodicity. Toprobe the long-term changes in the orbital architecture ofthe (258656) 2002 ES – 2013 CC couple, we extendedour previous simulation to Myr in the past. We note inpassing that the necessity to seek this pair’s age over such along timespan forces us to abandon any hopes of finding aconvergence in Cartesian coordinates. This is because of thesmall but non-negligible chaoticity of the integrated orbits,and principally results from an uncertainty in the thermalaccelerations that the objects would experience (as discussedbelow). Both of these factors would require a large numberof clones of (258656) 2002 ES and 2013 CC to investigatetheir past histories, and thus are computationally prohibitiveto pursue. We therefore choose to downsize the dimensional-ity of the space where a convergence is quantified, and focuson the behavior of secular evolution in just the non-singular elements z and ζ . Figure 8 shows the differences between theosculating δ Ω and δ̟ of the two objects, and pays specialattention to the time interval near δ Ω ≃ configurations.As expected, the first such configuration occurred about Myr ago. However, a closer look at the relevant panel ofFig. 8 indicates that suitable orbital convergence conditionsdid not occur at that time. Unlike ≃ . Mya, the orbitalplanes converge, but the perihelion longitudes are at themaximum of their oscillations. An even closer look at theepochs near nodal convergence shows that when δ̟ crosseszero, δ e is large, and vice versa. Once again, we thereforefind that the conditions of a low-velocity separation of thetwo orbits cannot be met at that epoch.Inspecting further epochs of nodal crossing, as shownin Fig. 8, we conclude that δ Ω ≃ in fact never exactlycoincides with δ̟ ≃ , a convergence pre-requisite. Here,however, we must revisit some of the assumptions madein our simulation. In particular, recall that (i) we usedonly nominal realizations of the orbits of both (258656)2002 ES and 2013 CC , and (ii)we included only gravita-tional perturbations from planets in our dynamical model.Both of these approximations are insufficient for a fullanalysis of our problem (see a similar discussion of theattempts to determine the origin of young asteroid clus-ters/families and pairs in Nesvorn´y & Vokrouhlick´y (2006),or Vokrouhlick´y & Nesvorn´y (2008)).First, the nominal orbital solution represents the best-fit of the available astrometric data. The inevitable uncer-tainties of the latter implies the uncertainty of the orbital fititself. Well-behaved orbital solutions are represented by fixedconfidence-level regions in the six-dimensional orbital space,using an ellipsoidal geometry, mathematically expressed byelements organized in the covariance matrix. Each orbitstarting in a high confidence-level zone ( ≥ − %, say)is statistically equivalent to the best-fit solution. whilst ini-tially very compact, these different solutions typically di-verge with time. We thus need to consider in our simulationnot only the best-fit orbits, but also a sample of those start-ing from the high-confidence zone. We call these“geometricalclones”.The second issue that needs to be considered is the va-lidity of the dynamical model used. the long-term dynamicsof small objects are known to be subject to perturbationsdue to the thermal acceleration known as the Yarkovskyeffect (e.g., Bottke et al. 2006; Vokrouhlick´y et al. 2015).Nominally, within the Trojan population, objects are onlyminimally affected by the Yarkovsky effect (Wang & Hou2017; Hellmich et al. 2019), which has the greatest influ-ence at smaller sizes. However, the two components in the(258656) 2002 ES -2013 CC couple are well within thissize range, and so it is warranted to see what dynamical ef-fects might be produced by Yarkovsky accelerations. Sincenone of the parameters needed for evaluation of the thermalaccelerations, such as the rotation state, the surface thermalinertia, and the bulk density, are known for either (258656)2002 ES or 2013 CC , we need to consider a suite of po-tential orbit histories, each generated by numerical integra-tion of test particles experiencing the a range of physicallyplausible thermal accelerations . These will be called theYarkovsky clones. We also note that the effect of thermalaccelerations was included in swift_rmvs4 using the samemethod as described in Vokrouhlick´y & Nesvorn´y (2008). MNRAS , 1–22 (2020) T. R. Holt et al.
Figure 8. the long-term behavior of the difference in osculating nodal and perihelion longitudes δ Ω (red) and δ̟ (blue) for the nominalorbits of (258656) 2002 ES and 2013 CC . The top panel shows the results from a backward integration in time to Myr. Thefour panels below show a zoom around the configurations where δ Ω becomes small, also indicated by the black rectangles in the toppanel. As inferred from data in Fig. 7, the first such situation occurs ≃ Myr in the past, and repeats with a period of ≃ Myr.The configuration of the nominal orbit becomes closest to true convergence at ≃ Myr and ≃ Myr in the past (right middle andbottom panels).
We conducted two sets of numerical simulations, one con-sidering only the geometrical clones (section 3.3.1), and theother considering only the Yarkovsky clones (section 3.3.2)of (258656) 2002 ES and 2013 CC . In each simulationset, we include the nominal orbit of the objects, comple-mented by a set of 20 clones. We ran a backward integra-tion of all orbits for . Gyr with an integration timestepof days. Every years, we evaluated the differences be-tween the osculating orbital elements of the 21 realizationsof (258656) 2002 ES with each of those of 2013 CC , and searched for the possibility of a convergent configuration.To quantify the latter, we used two conditions. First, as inNesvorn´y & Vokrouhlick´y (2006), we evaluated the targetfunction δ V = na q ( sin I δ Ω ) + . ( e δ̟ ) , (3)where ( n , a , e , I ) are the arithmetically-mean values of themean motion, semimajor axis, eccentricity and inclination ofthe two considered clones, and δ Ω and δ̟ are the differencesbetween the osculating longitude of the ascending node andperihelion for the two clones, respectively. This way, δ V hasthe dimension of velocity, and is constructed to provide, in MNRAS , 1–22 (2020) pair of Jovian Trojans a statistically mean sense, the necessary velocity perturba-tion required for a transfer between the secular angles of thetwo orbits. However, the analysis of the results presented inFig. 8 has shown that even a configuration with potential δ Ω ≃ and δ̟ ≃ , and therefore δ V ≃ , is not enough toguarantee a satisfactory orbital convergence, provided that δ e and δ I are simultaneously large. For that reason, we ad-mit as a potentially convergent configuration a case wherethe orbits of the two clones satisfy • δ V ≤ V lim , where V lim is some small value, we use typi-cally − m s − , and • δ e ≤ e lim and δ I ≤ I lim , where again we use suitablysmall values of e lim ≃ × − and I lim ≃ . ◦ namely dif-ferences in the corresponding proper elements of (258656)2002 ES and 2013 CC .We output information about these potentially convergingconfigurations for further analysis. In the next two sections,we comment on the results of our numerical experimentsthat use geometrical (section 3.3.1) and Yarkovsky clones(section 3.3.2) separately. Information about the orbit determination, needed for a con-struction of the geometrical clones, was taken from the
Ast-DyS database. the orbits of both (258656) 2002 ES and2013 CC are rather well constrained, reflecting numer-ous astrometric observations. Even the poorer of the two,2013 CC , was observed over seven oppositions, leading toa fractional accuracy of ≃ − in the semimajor axis, a ,and the Cartesian components of the non-singular elements, z and ζ . Only the mean longitude, λ , has a slightly worseaccuracy, namely ≃ × − degrees. These are the character-istic differences between the six orbital osculating elements E = ( a , z , ζ, λ ) of the clones in ≃ % confidence zone and thebest-fit solution E ⋆ . The solution is given at the initial epochMJD58800. Complete information about the parameters ofthe six-dimensional confidence zone ellipsoid in the space ofelements E is given by the covariance and normal matrices,also provided at the AstDyS website. Denoting Σ the normalmatrix, we may construct the initial orbital elements E ofthe geometric clones using E = T T ξ + E ⋆ , (4)where ξ is a six-dimensional vector whose componentsare random deviates of normal distribution (with vari-ance equal to unity), and the matrix T satisfies T T T = Σ (e.g., Milani & Gronchi 2010); T is obtained using theCholesky decomposition method. As mentioned above, weconstructed 20 geometric clones of both (258656) 2002 ES and 2013 CC at the initial epoch of our simulation.The bottom panel of Fig. 9 shows the maximum nodaldifference between the clones of (258656) 2002 ES and itsnominal orbit. Tiny differences between the orbital param-eters imply that the s frequency of the clone orbits is notexactly the same as that of the nominal orbit. However, thestability of this orbital zone ensures that the configuration ofthe clone orbits does not evolve, and thus initially the nodaldivergence is basically linear in time. Only beyond about . Gyr does the divergence become stronger than linear.
Figure 9.
The statistical distribution of convergent solutions forgeometric clones of (258656) 2002 ES and 2013 CC from sim-ulations of the nominal orbits of the two objects, plus 20 cloneseach, using the velocity cutoff δ V ≤ m s − , and eccentricity andinclination limits discussed in the text. Abscissa is time to thepast starting from Mya (there are no earlier solutions).Theleft ordinate in the upper two panels gives the number of recordedsolutions in kyr bins (red histogram). The gray line gives | δ Ω | of the nominal orbits of (258656) 2002 ES ’s and 2013 CC (seethe right ordinate and the red line on Fig. 8), aiming to aid in-terpretation of the results. The green line in the bottom panelshows the maximum difference in the longitude of the ascendingnode between the clones of (258656) 2002 ES and the longitudeof ascending node of its nominal orbit; up to about Myr thistrend is nearly linear, but becomes more complicated beyond thisepoch due to very weak orbital chaos.
This is an expression of a very weak instability that man-ifests itself in the behavior of the secular angle solely Gyrtimescales . The formal Lyapunov timescale of the orbits ofboth (258656) 2002 ES and 2013 CC is only ≃ Myr(see the
AstDyS database). This implies that a divergencein λ is dominant, whilst the divergence in the secular anglesis slower, as shown in Fig. 9. At Gyr, the nodal longi-tudes of clones of (258656) 2002 ES are thus spread over a ≃ ◦ range. A similar, and potentially slightly larger, effectis seen among the clones of 2013 CC , principally due totheir larger differences at the initial epoch. This divergencemay overcome the difficulties we experienced in attempt-ing to find an epoch at which the nominal orbits achieve aconverging configuration. For instance, in the bottom rightpanel of Fig. 8, we note that the nodal difference of thenominal orbits misses the epoch at which the difference ofpericenters basically shrinks to zero by about ◦ at ≃ Gyr.This may be compensated for if the orbits of suitable clones
MNRAS , 1–22 (2020) T. R. Holt et al. are used, instead of the nominal orbits. Obviously, a sat-isfactorily large nodal spread of the clone orbits must beattained.The top panel of Fig. 9 shows the statistical distribu-tion of the converging geometric clones of the two Trojans,organized in kyr wide bins. Obviously, the rather smallnumber of clones in our test run does not allow us to probethe convergence properties in great detail. For that reason,and with the rather tight limit δ V ≤ m s − chosen, thepossible solutions cluster only near the ≃ Myr epoch,though we note that, if a looser criterion δ V ≤ m s − was chosen, more solutions would also exist at ≃ Myr.Taken naively at a face value, we would conclude a possibleorigin of the (258656) 2002 ES – 2013 CC couple at thistime in the past, if the couple are not older than . Gyr,beyond which we did not continue our simulation. However,as is often in the case of a pair configuration which is notvery young, the so far neglected thermal accelerations in thedynamical model can prove to be a source of considerableuncertainty.This is analysed in section 3.3.2.
Our Yarkovsky clones all have the same initial conditionsas the nominal orbit, but they differ in the magnitude ofthermal accelerations used for their orbital propagation. Asin Vokrouhlick´y & Nesvorn´y (2008), we approximate ther-mal accelerations using a simple transverse component withthe magnitude inversely proportional to the square of theheliocentric distance. The magnitude of this acceleration isadjusted such that the resulting change in the semimajoraxis da / dt matches predictions from the theoretical formu-lation of Yarkovsky effect (see also Farnocchia et al. 2013,where a classical formalism used in cometary dynamics wasadopted). In order to estimate plausible da / dt values, weuse a simple approach describing the diurnal Yarkovsky ef-fect for a spherical body on a circular heliocentric orbit,presented in Vokrouhlick´y (1998). We use the following setof physical parameters: the surface thermal conductivity K ≃ . − . W m − K − , the surface thermal inertia Γ ≃ − [SI units] (for both see Delb´o et al. 2015),the bulk density ρ ≃ . g cm − (e.g., Carry 2012), rotationperiod P ≃ − hr, and size D ≃ km. The max-imum semimajor axis drift rate at zero obliquity is then ( da / dt ) max ≃ ( . ± . ) × − au Myr − . Our choice of aslow rotation period is tied to the working assumption that(258656) 2002 ES and 2013 CC are indeed a real Trojanpair. We argue in section 4.1 that the most plausible for-mation mechanism for such a pair is the destabilization ofa Trojan binary. If this is indeed the case, then before theirseparation, the two components were most likely spin-orbitsynchronized to periods of ≥ hr (e.g., Nesvorn´y et al.2020). If, however, the formation mechanism of the pair wasdifferent, such as the YORP-driven fission of a parent ob-ject (see section 4.2), the rotation periods P of (258656)2002 ES and 2013 CC could well be as short as a fewhours. In that case, ( da / dt ) max would be smaller by a fac-tor of to . Indeed, as a confirmation of our reasoning, wenote that scaling the value of the detected Yarkovsky signal × − au Myr − for the m size near-Earth asteroid101955 Bennu with P ≃ . hr (e.g., Chesley et al. 2014), wewould have ( da / dt ) max ≃ . × − au Myr − . In our simula- Figure 10. the statistical distribution of convergent solutionsfor the Yarkovsky clones (nominal orbits plus 20 clones each)of (258656) 2002 ES and 2013 CC , using the velocity cutoff δ V ≤ m s − , and eccentricity and inclination limits discussed inthe text. Abscissa is time to the past starting from Mya (thereare no earlier solutions). the left ordinate in the upper two panelsgives the number of recorded solutions in kyr bins. The toppanel (red histogram) gives the number of solutions for all pos-sible combinations of clones. The middle panel (blue histogram)for the case when only clones with the same sign of da / dt werecompared. The gray line gives | δ Ω | of the (258656) 2002 ES ’sand 2013 CC ’s nominal orbits (see the right ordinate and thered line on Fig. 8), aiming to aid interpretation of the results.The green line in the bottom panel shows the difference in thelongitude of ascending node between the Yarkovsky clone withmaximum positive drift rate ( da / dt ) max and the nominal orbit of(258656) 2002 ES . tion, we consider only the case of long rotation periods, andfix ( da / dt ) max ≃ . × − au Myr − . For each of the twoTrojans, (258656) 2002 ES and 2013 CC , we considerthe nominal orbit with da / dt = , and Yarkovsky clones.In both cases, clones have positive da / dt and cloneshave negative da / dt . Additionally, because in the case of thediurnal variant of the Yarkovsky effect da / dt ∝ cos γ , where γ is the spin axis obliquity, the positive/negative close da / dt values uniformly sample the interval to ( da / dt ) max , resp. −( da / dt ) max to . MNRAS000
Our Yarkovsky clones all have the same initial conditionsas the nominal orbit, but they differ in the magnitude ofthermal accelerations used for their orbital propagation. Asin Vokrouhlick´y & Nesvorn´y (2008), we approximate ther-mal accelerations using a simple transverse component withthe magnitude inversely proportional to the square of theheliocentric distance. The magnitude of this acceleration isadjusted such that the resulting change in the semimajoraxis da / dt matches predictions from the theoretical formu-lation of Yarkovsky effect (see also Farnocchia et al. 2013,where a classical formalism used in cometary dynamics wasadopted). In order to estimate plausible da / dt values, weuse a simple approach describing the diurnal Yarkovsky ef-fect for a spherical body on a circular heliocentric orbit,presented in Vokrouhlick´y (1998). We use the following setof physical parameters: the surface thermal conductivity K ≃ . − . W m − K − , the surface thermal inertia Γ ≃ − [SI units] (for both see Delb´o et al. 2015),the bulk density ρ ≃ . g cm − (e.g., Carry 2012), rotationperiod P ≃ − hr, and size D ≃ km. The max-imum semimajor axis drift rate at zero obliquity is then ( da / dt ) max ≃ ( . ± . ) × − au Myr − . Our choice of aslow rotation period is tied to the working assumption that(258656) 2002 ES and 2013 CC are indeed a real Trojanpair. We argue in section 4.1 that the most plausible for-mation mechanism for such a pair is the destabilization ofa Trojan binary. If this is indeed the case, then before theirseparation, the two components were most likely spin-orbitsynchronized to periods of ≥ hr (e.g., Nesvorn´y et al.2020). If, however, the formation mechanism of the pair wasdifferent, such as the YORP-driven fission of a parent ob-ject (see section 4.2), the rotation periods P of (258656)2002 ES and 2013 CC could well be as short as a fewhours. In that case, ( da / dt ) max would be smaller by a fac-tor of to . Indeed, as a confirmation of our reasoning, wenote that scaling the value of the detected Yarkovsky signal × − au Myr − for the m size near-Earth asteroid101955 Bennu with P ≃ . hr (e.g., Chesley et al. 2014), wewould have ( da / dt ) max ≃ . × − au Myr − . In our simula- Figure 10. the statistical distribution of convergent solutionsfor the Yarkovsky clones (nominal orbits plus 20 clones each)of (258656) 2002 ES and 2013 CC , using the velocity cutoff δ V ≤ m s − , and eccentricity and inclination limits discussed inthe text. Abscissa is time to the past starting from Mya (thereare no earlier solutions). the left ordinate in the upper two panelsgives the number of recorded solutions in kyr bins. The toppanel (red histogram) gives the number of solutions for all pos-sible combinations of clones. The middle panel (blue histogram)for the case when only clones with the same sign of da / dt werecompared. The gray line gives | δ Ω | of the (258656) 2002 ES ’sand 2013 CC ’s nominal orbits (see the right ordinate and thered line on Fig. 8), aiming to aid interpretation of the results.The green line in the bottom panel shows the difference in thelongitude of ascending node between the Yarkovsky clone withmaximum positive drift rate ( da / dt ) max and the nominal orbit of(258656) 2002 ES . tion, we consider only the case of long rotation periods, andfix ( da / dt ) max ≃ . × − au Myr − . For each of the twoTrojans, (258656) 2002 ES and 2013 CC , we considerthe nominal orbit with da / dt = , and Yarkovsky clones.In both cases, clones have positive da / dt and cloneshave negative da / dt . Additionally, because in the case of thediurnal variant of the Yarkovsky effect da / dt ∝ cos γ , where γ is the spin axis obliquity, the positive/negative close da / dt values uniformly sample the interval to ( da / dt ) max , resp. −( da / dt ) max to . MNRAS000 , 1–22 (2020) pair of Jovian Trojans Figure 11.
Two examples of converging solutions between Yarkovsky clones of (258656) 2002 ES and 2013 CC : left panels at ≃ . Mya, right panels at ≃ . Mya (gray vertical lines show the nominal convergence epochs). Each of the panels shows thedifferences between the osculating orbital elements of the clones: eccentricity (top), inclination (middle), and longitude of node (red) andperihelion (blue; bottom). The secular angles Ω and ̟ converge to better than . ◦ , corresponding to a negligible value of the targetfunction δ V ≤ . m s − (see Eq. 3). Differences in e and I are relatively larger, namely δ e ≃ . × − and δ I ≃ . ◦ (left), resp. δ I ≃ . ◦ (right). The dashed horizontal lines show the differences between the respective proper elements of (258656) 2002 ES and2013 CC . Note the ordinate of the middle panel (inclination) which is offset from zero. Figure 10 shows the results from our Yarkovsky clonesimulations. In contrast to the simulations where only thegeometrical clones were used (Fig. 9), there are many moreconvergent solutions, starting from
Mya. The reason isillustrated in the bottompanel of Fig. 10, which shows thedivergence of the osculating longitude of the ascending nodebetween the nominal orbit (no Yarkovsky effect) and theclone with the maximum positive drift-rate ( da / dt ) max of(258656) 2002 ES . Clones with smaller da / dt values havenodal differences smaller than the signal seen in Fig. 10,proportionally to their cos γ value.The nodal differences between various clones are nowmuch larger, reaching the maximum possible value of ◦ after at ≃ . Gya. The nodal difference to the nominal orbitof the clone with the maximum negative drift-rate value isabout the same but negative. This is because δ Ω now prop-agates nearly quadratically in time as opposed to the quasi-linear trend for the geometrical clones. Such a quadratictrend in node propagation is characteristic of Yarkovskystudies of asteroids (e.g., Vokrouhlick´y & Nesvorn´y 2008).In that case, the phenomenon was easily associated withthe principal dynamical perturbation produced by theYarkovsky effect, namely the secular drift in semimajor axis. As a result, the semimajor axis dependence of the s frequency produces, after a straightforward integration, aquadratic-in-time drift of the node. In our case of JovianTrojans, the effects are slightly subtler. This is because, inspite of a permanent transverse perturbing acceleration inorbits of the clones, their semimajor axis does not show anyconstant drift in time due to the resonant locking inherentto their presence in the Trojan population. However, otherelements –eccentricity and inclination– do display such a sec-ular drift, as previously found in Wang & Hou (2017) andHellmich et al. (2019). As the s frequency is also dependenton these values, it still displays a linear change as a functionof time, explaining the quadratic effect in node seen in theFig. 10.Returning to the pattern in the distribution of converg-ing solutions seen in Fig. 10, we note their clustering nearepochs when δ Ω of the (258656) 2002 ES and 2013 CC nominal orbits has been found to reach zero (the grey linein the top panels). This is to be expected, since the nodaldifference exhibits the most stable evolution in time. There-fore, when nominal orbits of the two Trojans have large δ Ω values, the clones will also follow the same pattern. This con-clusion will, however, weaken further into the past because MNRAS , 1–22 (2020) T. R. Holt et al. of the clone nodal divergence discussed above. As a result,beyond ∼ one Gyr into the past, the solution distributionspreads more in time. This is because specific clone combi-nations may now satisfy more easily our convergence con-ditions. Additionally, convergent solutions cluster in peaksseparated by about Myr, rather than exhibiting a contin-uous distribution about the δ Ω ≃ nodal conditions. This isdue to the δ̟ ≃ perihelion condition also facilitating theconvergence criteria we adopted.The middle panel in Fig. 10 shows the statistical dis-tribution of the number of converging solutions for a sub-sample of cases in which clones of (258656) 2002 ES and2013 CC both have the same sign of the associated da / dt drift. Translated using the diurnal Yarkovsky theory, thisalso implies that the two clones have the same sense of rota-tion: either both prograde, or both retrograde. The proposedformation mechanisms for this pair, namely a binary split orrotation fission, would both predict this property. There areobviously fewer solutions found, but the general pattern oftheir distribution is about the same as in the general casewhen all clones are taken into account.Figure 11 shows the conditions at convergence for twopairs of the Yarkovsky clones of (258656) 2002 ES and2013 CC : the left panels at the most recent possible clus-ter of solutions in the past (namely at ≃ . Mya), whilstthe right panel shows the cluster at an epoch which is moredistant in the past by two cycles of the differential motion oftheir orbital nodes (namely at ≃ . Mya). In general,the quality of the convergence is similar, including those so-lutions beyond Gya. In both cases, the formal convergenceof the secular angles is better than . ◦ .When inserted into Eq. (3), the equivalent velocity dif-ference is negligibly small δ V ≤ . m s − . At the con-vergence epoch, the osculating eccentricity values are alsosatisfactorily close to each other, namely δ e ≃ . × − .Using the Gauss equations (e.g., Nesvorn´y & Vokrouhlick´y2006), we estimate that this tiny eccentricity difference cor-responds to an orbital velocity change smaller than m s − in a statistical sense. This change is actually smaller than thedifference in proper eccentricity values of (258656) 2002 ES and 2013 CC . The inclination convergence turns out to bethe most troublesome element of the simulation: the per-sisting differences of ≃ . ◦ statistically correspond to avelocity change of ≃ m s − . Such a difference in the os-culating values of inclinations corresponds to the differenceof their proper values. In contrast, the acceptable true sep-aration velocity of the objects should be a fraction of theescape velocity from the effective parent body. With its sizeof ≃ km, the ideal condition of the separation in this pairwould require a velocity difference of ≤ m s − . The incli-nation difference at converging solutions is therefore nearlyan order of magnitude larger.One possibility to explain this mismatch may be relatedto our approximation of the Yarkovsky effect. By represent-ing it using the transverse acceleration only, the inclinationis not perturbed. In fact, a complete model of the thermalaccelerations may admit an out-of-plane component, pro-vided that the obliquities of the components of the pair arenot extreme (e.g., Vokrouhlick´y 1998). However, to fully usesuch a model, we would need to sample a multi-parametricspace of possible spin orientations and physical parameters for Yarkovsky clones, an effort which is postponed to furtherstudies.An alternative dynamical mechanism, that has not beenincluded in our simulations, consists of perturbations fromthe largest Trojans in the L swarm. As an example, weconsider the influence 624 Hektor, whose mass is estimatedto be ≃ kg (e.g., Carry 2012), about − of the massof dwarf-planet 1 Ceres. Nesvorn´y et al. (2002b) found that,statistically, the mean perturbation of the orbital inclina-tion produced by Ceres in the inner and middle parts of theMain belt is ≃ . ◦ in Gyr. Assuming the effect scales withthe square root of the perturber mass, we estimate that theapproximate effect of Hektor on small L Trojans would be ≃ . ◦ over Gyr, in a statistical sense. Therefore, at leasta part of the inclination mismatch reported above could wellbe due to the ongoing scattering influence of the most mas-sive Trojans.
We now briefly discuss possible formation processes for the(258656) 2002 ES –2013 CC pair. In principle, thesemechanisms coincide with the suggestions outlined in Sec. 6of Vokrouhlick´y & Nesvorn´y (2008). Building on that work,we will skip for now the possibility that these two Trojansare the two largest objects in a compact, collisionally-bornfamily. Given their comparable size, the collision requiredto form such a family must have been super-catastrophic,with many kilometer size fragments created and dominatingthe mass. Without information about them, it is hard tosay anything more about the putative collision conditions,including the probability of such a collision actually havingoccurred. The first possible origin for the (258656) 2002 ES –2013 CC pair consists of a model, in which the two objectswere formerly components in a binary system which under-went some kind of instability. We assume that the instabil-ity was not of a dynamical origin. Indeed, even if formedby gravitational collapse, the initial angular momentum ofthe binary would exceed that of a critically rotating sin-gle body of an equivalent mass by a factor of ≃ ( − ) (Nesvorn´y et al. 2019). This is not sufficient to drive tidalevolution, whilst conserving angular momentum, to the sta-bility limit at about half of the Hill sphere, even in theTrojan zone. The limiting configuration would require an-gular momentum at least twice as large. Additionally, timeconstraints may prevent evolution to such large separationswithin ≤ . Gyr. Therefore, the nature of the parent bi-nary instability must be different. We assume instead thatthis instability was triggered by a gentle-enough impact onone of the components. We leave aside other possibilities,such as binary instability produced during a close threebody encounter with a massive Trojan (Agnor & Hamilton2006; Nesvorn´y & Vokrouhlick´y 2019),for future investiga-tions, once the mechanisms are better understood in theJovian Trojan population.Let us start the likelihood analysis of the formation
MNRAS , 1–22 (2020) pair of Jovian Trojans of the (258656) 2002 ES –2013 CC pair via the sub-critical impact dissociation of a previously existing syn-chronous binary with a very simple, order-of-magnitude esti-mate. Assume that the needed imparted velocity by the im-pact onto a ≃ km size component in the binary is about m s − . Then, using the simple formulation in Nesvorn´y et al.(2011), a projectile of ≃ . km size is required. The char-acteristic impact velocity assumed was V imp ≃ . km s − (Davis et al. 2002). The Trojan population contains very ap-proximately N ≃ such objects (e.g., Wong & Brown2015; Emery et al. 2015, and Fig. 12).Using the mean impact probability p i ≃ × − km − yr − (e.g., Davis et al. 2002), we can therefore esti-mate the order-of-magnitude likelihood that such an eventwould occur within a timeframe of T ≃ . Gyr, namely p i R NT ≃ . (here R = . km is the radius of the tar-get body). This suggests that every such binary implantedto the Trojan population has a non-negligible (15%) chanceto be split via this process. Assuming that, initially, at leasthundreds of binaries were captured intact to the Trojan pop-ulation, a non-negligible number of Trojan pairs might havebeen created over the age of the Solar system. Obviously,in many cases, our ability to identify the pair produced inthis manner is low, due to unsuitable locations in the Tro-jan orbital phase space. Nonetheless, this result suggest thatsufficiently many such pairs could be produced that futurestudy might well reveal several more.We now substantiate this order-of-magnitude estimateusing a more involved numerical simulation. As outlinedabove, the mutual orbit of a binary can be affected by smallimpacts on to its components. The binary may become un-bound if the velocity change imparted by an impact ex-ceeds binary’s orbital speed ∼ . − m s − for bodies with D ≃ km (Petit & Mousis 2004).We investigate this process with the previously devel-oped collisional code (Morbidelli et al. 2009; Nesvorn´y et al.2011). The code, known as Boulder , employs a sta-tistical method to track the collisional fragmentationof planetesimal populations. A full description of the
Boulder code, tests, and various applications can befound in Morbidelli et al. (2009), Levison et al. (2009) andBottke et al. (2010). The binary module in
Boulder ac-counts for small, non-disruptive impacts on binary compo-nents, and computes the binary orbit change depending onthe linear momentum of impactors (see Nesvorn´y et al. 2011;Nesvorn´y & Vokrouhlick´y 2019).We account for impacts over the life of the Solar sys-tem, . Gyr. The captured population of Jovian Trojans isassumed to be similar to the present population, for ob-jects with large diameters. There are ≃ Trojans with D > km. The population is assumed to follow a powerlaw profile below km, with a cumulative index equal to − . (Fig. 12). The intrinsic impact probability and impactvelocity is the same as used for the order-of-magnitude esti-mate above. We adopt a standard disruption law for solid icefrom Benz & Asphaug (1999). Fragments are generated ac-cording to the method described in Morbidelli et al. (2009).These rules are implemented in the Boulder code, whichis then used to determine the collisional survival of Trojanbinaries (e.g., Nesvorn´y et al. 2018).Figure 12 shows the evolution of the size distributionfor the Jovian Trojans. The size distribution for D > km Figure 12.
The effects of collisional grinding on the cumulativesize distribution of Jovian Trojans. The upper bold line is theinitial distribution. The lower bold line is the size distribution at T = . Gyr. The gray lines show the changing size distributionin
Myr intervals. The dip in the final distribution near D = . km is produced by the strength-to-gravity transition of thedisruption law. Figure 13.
The survival probability of binaries with (258656)2002 ES –2013 CC components as a function of separation,here normalized to the sum of physical radii, R + R . The sur-vival probability decreases with separation because wide bina-ries have smaller orbital speeds and are easier to dissolve bya small impact. For reference, the Hill radius R H of a binarywith (258656) 2002 ES –2013 CC components, correspondingto mass ∼ × g (for g cm − density), is R H ≃ , km, ornearly a B /( R + R ) = . remains unchanged over . Gyr, but below D ≃ km theslope becomes shallower. This is consistent with Jovian Tro-jan observations that detect a shallower slope for D ≃ km(e.g., Wong & Brown 2015). If this interpretation is correct,the slope should become steeper below approximately m,for bodies that are too faint to be detected from the groundusing the current generation of observatories. The dip in thesize distribution is produced by the transition from strength-to-gravity dominated branches of the disruption law (e.g.,Nesvorn´y et al. 2018).We find that the survival chances of Trojan binariesare generally good, but drop significantly when the binaryseparation approaches . R H ( R H being the Hill sphere of MNRAS , 1–22 (2020) T. R. Holt et al. gravitational influence, see Fig. 13). This is expected be-cause binaries with semimajor axis a B > . R H are dynam-ically unstable (e.g., Porter & Grundy 2012). For a char-acteristic separation of a B /( R + R ) ≃ − , where a B is the binary semimajor axis and R and R are the bi-nary component radii, consistent with the pair (258656)2002 ES –2013 CC ( R + R = . km), which is quite com-mon among equal-size binaries in the Edgeworth-Kuiper belt(e.g., Noll et al. 2020), the survival probability is − %.There is plenty of room in this parameter space for Trojanpair formation by this mechanism. Assuming that the pair(258656) 2002 ES –2013 CC is an impact-dissolved bi-nary, we find that there should be . − . surviving binariesfor each pair such as (258656) 2002 ES –2013 CC . Giventhat the vast majority of Trojan pairs remain undetected(see the difficulties briefly outlined in the Appendix B),the obvious implication is that there should also be severalequal-size binaries among Jovian Trojans in this size range. An alternative formation mechanism that could explain theobserved properties of the (258656) 2002 ES –2013 CC pair is that they might be the result of the rotational fis-sion of their common parent object (this is indeed the fa-vorite mechanism for asteroid pair formation in the mainbelt; e.g., Pravec et al. 2010). The most probable drivingprocess for such a fission event is the Yarkovsky-O’Keefe-Radzievski-Paddack (YORP) effect, a radiative torque re-sulting from the combination of reflected and thermallyemitted radiation by the surface (being thus a complemen-tary phenomenon to the Yarkovsky effect; e.g., Bottke et al.2006; Vokrouhlick´y et al. 2015). The YORP effect is ableto constantly accelerate an asteroid’s rotation up to speedsthat meet the requisite conditions to cause the object tofission . The rotation frequency change Û ω satisfies generalscaling properties, such that Û ω ∝ /[ ρ ( aD ) ] , where ρ is thebulk density, a the orbital semimajor axis and D the size.However, the problematic part of the YORP effect, unlikethe Yarkovsky effect, is its large sensitivity to details of thesurface roughness. For that reason it is troublesome to de-termine the exact value of the strength of the YORP effectfor a given object, and we must satisfy ourselves with anorder-of-magnitude estimate in our case.If we were to determine the doubling timescale τ YORP = ω / Û ω (sometimes also the YORP cycle timescale; e.g.,Rubincam 2000), it would be reasonable to use the YORPdetection of the small near-Earth asteroid (101955) Bennuas a template, as we did above for the Yarkovsky effectin section 3.3.2. (101955) Bennu has τ YORP ≃ . Myr(e.g., Hergenrother et al. 2019). Adopting plainly the scal-ing τ YORP ∝ ρ ( aD ) / P (with P being the rotation period),we obtain τ YORP ≃ Gyr for a D ≃ km Trojan, theestimated size of a putative parent object of the (258656)2002 ES – 2013 CC pair. Note that τ YORP provides an es-timate of a timescale for doubling ω , as an example changingrotation period from hr to . hr, an approximate fissionlimit for a large internal strength Trojan model. Another τ YORP / ≃ . Gyr time would be needed if the initial ro-tation period of the parent object was hr. This shortertimescale would also be an appropriate estimate to reach thefission limit at a longer period of ≃ hr when the internal strength and bulk densities are low (e.g., French et al. 2015;Szab´o et al. 2017).If, however, we were to consider the results from numer-ical simulations of the YORP effect for a large statisticalsample of Gaussian-sphere shapes ˇCapek & Vokrouhlick´y(2004), which obtained τ YORP ≃ Myr for a typical mainbelt S-type asteroid of a km size, we would have τ YORP ≃ . Gyr for changing the parent object period from to . hr. Whilst these results are known to typically over-estimate the strength of the YORP effect by a factor of − , when compared to detections of the YORP effect forsmall near-Earth asteroids, we nonetheless get a timescaleshorter by a factor 2 to 3 than for the Bennu case. The take-away message is that the estimate of the YORP doublingtimescale prior the fission of the putative parent object of the(258656) 2002 ES and 2013 CC pair is very uncertain,with values ranging possibly from Gyr to some Gyr.Taken at a face value, the smaller values in this intervalare plausible as an explanation for the origin of the pair whencompared to the lifetime of the Solar system. It may not besurprising to find that some D ≃ km Jupiter Trojan objectsundergo a rotational fission during their lifetime. However,a more detailed inspection of the (258656) 2002 ES and2013 CC parameters speaks against this possibility. First,we note that the known rotation periods of Jovian Trojansrarely have values smaller than − hr (e.g., French et al.2015; Szab´o et al. 2017; Ryan et al. 2017), which suggests inturn that more than one τ YORP timescale would be needed toreach fission from a typical initial rotation state (though, ad-mittedly, these known data concern larger objects). More im-portantly, though, we note that the absolute magnitude dif-ference of (258656) 2002 ES and 2013 CC is ≃ ( . − . ) ,depending on the database used. This implies that the twoobjects are nearly of the same size. Pravec et al. (2010) ar-gued that the typical conditions of fission mechanics requireat least magnitude difference between the two compo-nents in pair. This is because some degree of size disparityis needed to make the two components separate onto dis-tinct heliocentric orbits. Whilst exceptions have been foundto this guideline (see e.g. Pravec et al. 2019), the majority ofthe known asteroid pairs, more than %, satisfy this con-dition of having a large enough magnitude disparity. Thecomponents in the (258656) 2002 ES – 2013 CC pair vi-olate this rule and would require special conditions for theirseparation to feasibly be the result of rotational fission. In this work, we identified the first potential dynamical pairin the Jovian Trojan population. In particular, we analysedthe distribution of Trojans in their proper orbital elementspace. Using information about the local density of objects,we also assessed the statistical significance of the proximityof potential couples. This procedure lead us to select a pair ofbodies, (258656) 2002 ES and 2013 CC , in the L4 swarmas a potential candidate pair. Interestingly, this suggestedpair is located very close to the L4 Lagrange point, withlow proper elements, semimajor axis ( da P ), eccentricity e P and sine of inclination ( sin I P ) values. Finally, as part of oureffort, we developed an up-to-date, highly accurate set of MNRAS , 1–22 (2020) pair of Jovian Trojans proper elements for the all Jovian Trojans, which we havemade publicly available (Appendix A).In order to further investigate the selected pair, we rana series of n -body simulations, which were used to look forpast convergences in the osculating nodal ( δ̟ ) and perihe-lion longitude ( δ Ω ) value for the two objects, whilst ensuringthat, at the time of such convergences the differences in theosculating eccentricity and inclination were also sufficientlysmall. Our simulations included both geometric clones, cre-ated from the uncertainties in the orbital elements of thebodies, and Yarkovsky clones, based on the estimated ther-mal accelerations that the two objects could experience, fora variety of realistic rotation rates. As a result, we obtained astatistical set of convergences, finding a larger pool of possi-bilities once the Yarkovsky clones were included. Our resultsreveal that the pair is at least ≃ Myr old, but are com-patible with the age being significantly older, potentially inthe Gyr time scale. By finding such possible convergences,we increase the confidence that the (258656) 2002 ES –2013 CC couple is a legitimate pair.We then considered the mechanisms by which the(258656) 2002 ES –2013 CC pair could have formed(compared with Vokrouhlick´y & Nesvorn´y 2008). The pairis not associated with any known collisional family, and assuch we do not favour the possibility of the pair having beenformed as a result of a catastrophic impact on a putativeparent body . The pair might have been formed throughthe rotational fission of their parent Trojan, since, for cer-tain initial conditions, the timescale for such an object tobe spun-up by the YORP effect to the point that it under-goes fission could be somewhat shorter than the age of theSolar system. However, this pair consists of two nearly equal-sized components, whilst the vast majority of observed pairsformed by rotational fission have a size ratio of at least . (see Pravec et al. 2010, 2019). For that reason, we considerthat the pair most likely formed as a result of the dissoci-ation of an equal-size binary. We can confirm that such ascenario is indeed feasible using an estimation of the binarysurvival rate in the size range of the (258656) 2002 ES –2013 CC pair, D ≃ km, over . Gyr, after implantationto the Trojan population early in Solar system’s history. Sta-tistically, this indicates that there should be many such pairswithin the Trojan population in this − km size range. Asthe Rubin Observatory’s Legacy Survey of Space and Time(LSST) comes online, it is expected to discover many Jo-vian Trojans in this size range (e.g., Schwamb et al. 2018).As new Trojans are discovered, our results suggest that fur-ther pairs should be revealed.The (258656) 2002 ES –2013 CC pair provides aninteresting clue to the past history of the Jovian Trojans,and the Solar system as a whole. So far, we know littlebeyond their dynamical properties and size estimations. Inparticular, lightcurve analysis could assist in constrainingthe formation mechanism, as this would provide an esti-mate of the rotational periods of the two objects. Due totheir small size, and dark albedo, the objects have rela-tively low apparent magnitudes, at best ≃ . magnitude invisible band. In order to further characterize these objects,observations using large Earth-based facilities, such as theSUBARU (Kashikawa et al. 2002) or Keck (Oke et al. 1995)telescopes, will be required. These objects would also ben-efit from future observations using the James Web (JWST, Rivkin et al. 2016) and Nancy Grace Roman Space Tele-scopes (RST, formerly known as WFIRST, Holler et al.2018). Time on these telescopes is competitive, but we rec-ommend proposals for observations of (258656) 2002 ES and 2013 CC be selected to further extend our understand-ing of this interesting pair of Trojans. DATA AVAILABILITY
The database of Jovian Trojanproper elements is accessible at https://sirrah.troja.mff.cuni.cz/~mira/mp/trojans_hildas/ ,and is available for community use. See Appendix A fordetails.
ACKNOWLEDGEMENTS
T.R.H was supported by the Australian Government Re-search Training Program Scholarship. The work of D.V. anM.B. was partially supported by the Czech Science Foun-dation (grant 18-06083S). This research has made use ofNASA Astrophysics Data System Bibliographic Services. Wethank Dr. Romina Di Sisto for their valuable review of thismanuscript.We dedicate this paper to the memory of Andrea Milaniand Paolo Farinella, who were the first to propose the idea ofa genetically connected pair of objects in the Jovian Trojanpopulation (Milani 1993).
REFERENCES
Agnor C. B., Hamilton D. P., 2006, Nature, 441, 192Beaug´e C., Roig F., 2001, Icarus, 153, 391Benjoya P., Zappal`a V., 2002, in Bottke W. F., Cellino A., Paolic-chi P., Binzel R. P., eds, Asteroids III. pp 613–618Benz W., Asphaug E., 1999, Icarus, 142, 5Bottke W. F., Vokrouhlick´y D., Rubincam D. P., Nesvorn´yD., 2006, Annual Review of Earth and Planetary Sciences,34, 157Bottke W. F., Nesvorn´y D., Vokrouhlick´y D., Morbidelli A., 2010,AJ, 139, 994Broˇz M., Rozehnal J., 2011, MNRAS, 414, 565Buie M. W., et al., 2015, AJ, 149, 113Carry B., 2012, Planet. Space Sci., 73, 98Chesley S. R., et al., 2014, Icarus, 235, 5Davis A. B., Scheeres D. J., 2020, Icarus, 341, 113439Davis D. R., Durda D. D., Marzari F., Campo Bagatin A., Gil-Hutton R., 2002, in Bottke W. F., Cellino A., Paolicchi P.,Binzel R. P., eds, Asteroids III. pp 545–558Delb´o M., Mueller M., Emery J. P., Rozitis B.,Capria M. T., 2015, in Michel P., DeMeo F. E.,Bottke W. F., eds, Asteroids IV. pp 107–128,doi:10.2458/azu uapress 9780816532131-ch006Di Sisto R. P., Ramos X. S., Beaug´e C., 2014, Icarus, 243, 287Di Sisto R. P., Ramos X. S., Gallardo T., 2019, Icarus, 319, 828Emery J. P., Marzari F., Morbidelli A., FrenchL. M., Grav T., 2015, in Michel P., DeMeo F. E.,Bottke W. F., eds, Asteroids IV. pp 203–220,doi:10.2458/azu uapress 9780816532131-ch011Farnocchia D., Chesley S. R., Vokrouhlick´y D., Milani A., SpotoF., Bottke W. F., 2013, Icarus, 224, 1French L. M., Stephens R. D., Coley D., Wasserman L. H., SiebenJ., 2015, Icarus, 254, 1MNRAS , 1–22 (2020) T. R. Holt et al.
Grav T., et al., 2011, ApJ, 742, 40Grav T., Mainzer A. K., Bauer J. M., Masiero J. R., Nugent C. R.,2012, ApJ, 759, 49Hellmich S., Mottola S., Hahn G., K¨uhrt E., de Niem D., 2019,A&A, 630, A148Hergenrother C. W., et al., 2019, Nature Communications,10, 1291Holler B. J., Milam S. N., Bauer J. M., AlcockC., Bannister M. T., Bjoraker G. L., 2018,J. Astron. Telesc. Instruments, Syst., 4, 1Holt T. R., et al., 2020, Mon. Not. R. Astron. Soc., 495, 4085Horner J., M¨uller T. G., Lykawka P. S., 2012,Mon. Not. R. Astron. Soc., 423, 2587Kashikawa N., et al., 2002, PASJ, 54, 819Lagrange J.-L., 1772, Essai sur le Probl`eme des Trois Corps, Prixde l’Acad´emie Royale des Sciences de Paris (printed in 1868,Œuvres de Lagrange, Tome VI, Gauthier-Villars, p. 229)Laskar J., Robutel P., 2001, Celestial Mechanics and DynamicalAstronomy, 80, 39Levison H. F., Duncan M. J., 1994, Icarus, 108, 18Levison H. F., Shoemaker E. M., Shoemaker C. S., 1997, Nature,385, 42Levison H. F., Bottke W. F., Gounelle M., Morbidelli A.,Nesvorn´y D., Tsiganis K., 2009, Nature, 460, 364Levison H. F., Olkin C. B., Noll K., Marchi S., LucyTeam 2017, in Lunar Planet. Sci. Conf.. p. 2025, http://adsabs.harvard.edu/abs/2017LPI....48.2025L
Marchis F., et al., 2006, Nature, 439, 565Margot J.-L., Pravec P., Taylor P., Carry B., Jacobson S., 2015a,in Michel P., DeMeo F. E., Bottke W. F., eds, Asteroids IV.pp 355–374, doi:10.2458/azu uapress 9780816532131-ch019Margot J.-L., Pravec P., Taylor P., Carry B., Jacobson S., 2015b,Asteroids IV, pp 355–373Milani A., 1993, Celestial Mechanics and Dynamical Astronomy,57, 59Milani A., Gronchi G. F., 2010, Theory of Orbital Determination.Cambridge University Press, CambridgeMorbidelli A., Bottke W. F., Nesvorn´y D., Levison H. F., 2009,Icarus, 204, 558Moskovitz N. A., et al., 2019, Icarus, 333, 165Nesvorn´y D., 2018, Annu. Rev. Astron. Astrophys., 56, 137Nesvorn´y D., Vokrouhlick´y D., 2006, AJ, 132, 1950Nesvorn´y D., Vokrouhlick´y D., 2019, Icarus, 331, 49Nesvorn´y D., Thomas F., Ferraz-Mello S., Morbidelli A., 2002a,Celest. Mech. Dyn. Astron., 82, 323Nesvorn´y D., Morbidelli A., Vokrouhlick´y D., Bottke W. F., BroˇzM., 2002b, Icarus, 157, 155Nesvorn´y D., Vokrouhlick´y D., Bottke W. F., 2006, Science,312, 1490Nesvorn´y D., Vokrouhlick´y D., Bottke W. F., Noll K., LevisonH. F., 2011, AJ, 141, 159Nesvorn´y D., Broˇz M., Carruba V., 2015, in Michel P., De-Meo F. E., Bottke W. F., eds, Asteroids IV. pp 297–321,doi:10.2458/azu uapress 9780816532131-ch016Nesvorn´y D., Vokrouhlick´y D., Bottke W. F., Levison H. F., 2018,Nature Astronomy, 2, 878Nesvorn´y D., Li R., Youdin A. N., Simon J. B., Grundy W. M.,2019, Nature Astronomy, 3, 808Nesvorn´y D., Vokrouhlick´y D., Bottke W. F., Levison H. F.,Grundy W. M., 2020, ApJ, 893, L16Noll K., Grundy W. M., Nesvorn´y D., Thirouin A.,2020, in Prialnik D., Barucci M. A., Young L., eds,The Trans-Neptunian Solar System. pp 201–224,doi:10.1016/B978-0-12-816490-7.00009-6Oke J. B., et al., 1995, Publ. Astron. Soc. Pacific, 107, 375Petit J. M., Mousis O., 2004, Icarus, 168, 409Porter S. B., Grundy W. M., 2012, Icarus, 220, 947Pravec P., Harris A. W., 2007, Icarus, 190, 250 Pravec P., Vokrouhlick´y D., 2009, Icarus, 204, 580Pravec P., et al., 2010, Nature, 466, 1085Pravec P., et al., 2019, Icarus, 333, 429Quinn T. R., Tremaine S., Duncan M., 1991, AJ, 101, 2287Rivkin A. S., Marchis F., Stansberry J. A., Takir D., Thomas C.,2016, Publ. Astron. Soc. Pacific, 128, 018003Robutel P., Gabern F., 2006, MNRAS, 372, 1463Rozehnal J., Broˇz M., Nesvorn´y D., Durda D. D., Walsh K.,Richardson D. C., Asphaug E., 2016, MNRAS, 462, 2319Ro˙zek A., Breiter S., Jopek T. J., 2011, MNRAS, 412, 987Rubincam D. P., 2000, Icarus, 148, 2Ryan E. L., Sharkey B. N. L., Woodward C. E., 2017, AJ, 153, 116Schwamb M. E., et al., 2018, preprint, 1802.01783( arXiv:1802.01783 )Szab´o G. M., et al., 2017, A&A, 599, A44Tsiganis K., Varvoglis H., Dvorak R., 2005,Celest. Mech. Dyn. Astron., 92, 71Vokrouhlick´y D., 1998, A&A, 335, 1093Vokrouhlick´y D., Nesvorn´y D., 2008, AJ, 136, 280Vokrouhlick´y D., Bottke W. F., Chesley S. R., ScheeresD. J., Statler T. S., 2015, in Michel P., DeMeoF. E., Bottke W. F., eds, Asteroids IV. pp 509–531,doi:10.2458/azu uapress 9780816532131-ch027Vokrouhlick´y D., et al., 2017, AJ, 153, 270Wang X., Hou X., 2017, MNRAS, 471, 243Wolf M., 1907, Astron. Nachrichten, 174, 47Wong I., Brown M. E., 2015, AJ, 150, 174ˇCapek D., Vokrouhlick´y D., 2004, Icarus, 172, 526ˇSidlichovsk´y M., Nesvorn´y D., 1996,Celestial Mechanics and Dynamical Astronomy, 65, 137
APPENDIX A: DETERMINATION OF THEJOVIAN TROJAN PROPER ELEMENTS
Here we briefly review our approach to compute syntheticproper elements for the currently known Jovian Trojan pop-ulation. The method is based on Milani (1993), see alsoBroˇz & Rozehnal (2011), though we needed several modi-fications of the digital filters in order to stabilize determi-nation of the proper elements for Trojans having very smalllibration amplitude. Our dynamical model included four gi-ant planets, with barycentric corrections to compensate forthe indirect perturbations for terrestrial planets. This ar-rangement suitably speeds up computations when dealingwith the whole population of many thousands of Trojans.Nevertheless, we also checked validity of our results usinga dynamical model including also the terrestrial planets ina full-fledged manner for a sub-sample of Trojans (notablythe low- δ V P that is of interest here). No significant differ-ences were observed. The initial planetary state vectors weretaken from the JPL ephemerides and those of the Trojansfrom the AstOrb catalogue as of April 28, 2020, from whichtheir population was also identified.We used well tested numerical package swift (e.g.,Levison & Duncan 1994), specifically the MVS2 symplecticintegrator (e.g., Laskar & Robutel 2001), that we adaptedfor our application in several ways. The most important wasan implementation of digital filters, helping us to eliminateshort-period and forced terms from osculating orbital ele-ments, necessary for identification of the proper terms. Dueto the absence of the direct perturbations from the terres-trial planets, we can allow a fixed integration timestep of . yr. The input sampling into the filtering routines was yr. We used a sequence of the convolution (Kaiser-window) MNRAS000
Here we briefly review our approach to compute syntheticproper elements for the currently known Jovian Trojan pop-ulation. The method is based on Milani (1993), see alsoBroˇz & Rozehnal (2011), though we needed several modi-fications of the digital filters in order to stabilize determi-nation of the proper elements for Trojans having very smalllibration amplitude. Our dynamical model included four gi-ant planets, with barycentric corrections to compensate forthe indirect perturbations for terrestrial planets. This ar-rangement suitably speeds up computations when dealingwith the whole population of many thousands of Trojans.Nevertheless, we also checked validity of our results usinga dynamical model including also the terrestrial planets ina full-fledged manner for a sub-sample of Trojans (notablythe low- δ V P that is of interest here). No significant differ-ences were observed. The initial planetary state vectors weretaken from the JPL ephemerides and those of the Trojansfrom the AstOrb catalogue as of April 28, 2020, from whichtheir population was also identified.We used well tested numerical package swift (e.g.,Levison & Duncan 1994), specifically the MVS2 symplecticintegrator (e.g., Laskar & Robutel 2001), that we adaptedfor our application in several ways. The most important wasan implementation of digital filters, helping us to eliminateshort-period and forced terms from osculating orbital ele-ments, necessary for identification of the proper terms. Dueto the absence of the direct perturbations from the terres-trial planets, we can allow a fixed integration timestep of . yr. The input sampling into the filtering routines was yr. We used a sequence of the convolution (Kaiser-window) MNRAS000 , 1–22 (2020) pair of Jovian Trojans filters A A B (e.g., Quinn et al. 1991) with decimation fac-tors 10 10 3, which were applied to the non-singular elements z = k + ı h = e exp ( ı̟ ) and ζ = q + ı p = sin I exp ( ı Ω ) . Theintermediate time window for this filtering procedure andoutput timestep was yr. At this stage, the short-periodterms with periods comparable to planetary orbital periodsor the libration period were efficiently suppressed from theresulting mean values ¯ z and ¯ ζ of eccentricity and inclina-tion variables. We then accumulated batches of 2048 valuesof ¯ z and ¯ ζ , and applied Fourier transformation (in particu-lar the FMFT method from ˇSidlichovsk´y & Nesvorn´y 1996),on the output. After rejecting signal associated with forcedplanetary frequencies (such as g , g or s to recall the prin-cipal ones), we were left with the proper values e P for theeccentricity and I P for the inclination as the amplitude ofthe remaining dominant terms. Our simulation spanned thetotal of Myr, and we computed proper elements in the ≃ kyr window described above many times over inter-vals with kyr shift in their origin. This way we had aseries of many tens of proper element realizations, allowingto access their stability and compute their mean and vari-ance. We also observed that the series of individual e P and I P still contained long-period signal (periods > Myr), whichin future studies may call for extension of integration win-dows. At this moment, we however, satisfied ourselves withour set-up. We also used the above outlined procedure forthe semimajor axis a , but instead of applying FMFT on itsmean values we simply computed its mean value ¯ a over a Myr interval. This helps us to determine semimajor axisvalue of the libration center for a given Trojan orbit.In order to obtain a reliable information about a stablelibration amplitude we need to apply a different method thathas been implemented in our code in parallel to computationof e P and I P . This is because the corresponding libration fre-quency is fast, f ≃ . deg yr − and ◦ / f ≃ yr, andmust not be under-sampled. A delicate issue consists of thefact that, at the same time, one has to suppress terms withperiod even shorter than the libration period, namely thosewhich are related to orbital periods of giant planets (prin-cipally Jupiter, ≃ . yr). We thus applied convolutionfilters B B, with decimation factors 3 3, to the osculatingvalues of the semimajor axis a and the longitude difference λ − λ ′ (the orbital elements labeled with prime correspond toJupiter), a resonant argument of the Trojan tadpole motion.These intermediate (mean) values of a and λ − λ ′ are com-puted with a yr cadence. In the next step, the intermediate a − a ′ were fitted by a straight line and the constant term a was subtracted. In the same way, the intermediate angle φ = λ − λ ′ − χ , where χ = ± ◦ depending on the L4 and L5libration points, was fitted by a straight line and the con-stant term φ was subtracted. Effectively, after subtractionsof the mean values was done, the tadpole motion around theLagrange point centers in these rescaled, zero-averages a − a ′ vs φ coordinates is centered at the origin. Consequently, thepolar angle ψ defined as (see, e.g., Milani 1993, a and a ′ inau) ψ = arctan (cid:18) a − a ′ . φ (cid:19) (A1)can be unfolded by ◦ , fitted by a straight line, with theslope defining the libration frequency f . The libration am-plitudes da P (in au) and D (in deg) are computed by the li b r a t i on a m p li t ude d a P [ au ] libration centre a [au] 0 0.05 0.1 0.15 0.2 e cc en t r i c i t y e p [] Figure A1.
Libration amplitude da P vs libration center ¯ a forTrojans in the L4 region. Colour corresponds to the proper eccen-tricity e P . The dependence of ¯ a ( da P , e P ) is systematic, indicatinga functional dependence. Fourier transform as amplitudes of spectral terms with fre-quency f . This second step uses a kyr cadence. Finally, weapply another averaging of da P and D values, defined on asimple running window with the output time step of Myr.Both da P and D may be considered as the third proper or-bital element alongside of e P and I P .We note that the value of libration center ¯ a is not uni-versal for all Trojans. Instead, its value functionally dependson the proper elements ( da P , e P , I P ) or ( D , e P , I P ) , see Fig. A1.Some authors (e.g., Broˇz & Rozehnal 2011; Rozehnal et al.2016) thus define an alternative set of proper elements ( a P = ¯ a + da P , e P , I P ) .We determined the above-introduced parameters,including different variants of orbital proper valuesand their uncertainty, for Jovian Trojans, popu-lation as of April 2020. These data can be found on https://sirrah.troja.mff.cuni.cz/~mira/mp/trojans_hildas/ . APPENDIX B: ARE THERE MORE LOW- δ V P COUPLES?
As also suggested by data in Fig. 6, the brief answer tothe topic of this Appendix is probably positive, but a fullanalysis if this issue is left to the future work. Here we onlyrestrict ourselves to illustrate difficulties one would quicklyface in attempting to prove the past orbital convergence ona Gyr timescales for most of the candidates.Let us consider another low- δ V P candidate couplecharacterized by small values of proper orbital elements ( da P , e P , sin I P ) , which helps to minimize the unrelated back-ground Trojan population (section 2). Staying near the L4libration point, we find 219902 (2002 EG ) and 432271(2009 SH ) at δ V P ≃ . m s − distance. This couple hasalso appreciably small probability p ≃ . × − to be arandom fluke and it has been highlighted by a green cir-cle in Fig. 6. The proper elements read da P ≃ ( . ± . ) × − au, e P ≃ ( . ± . ) × − and sin I P ≃( . ± . ) × − for (219902) 2002 EG , and da P ≃ MNRAS , 1–22 (2020) T. R. Holt et al.
Figure B1.
The same as Fig. 9 but for the (219902) 2002 EG –(432271) 2009 SH couple of Trojans: past orbital histories ofnominal orbits and 20 geometrical clones each compared every yr and convergent solutions within δ V ≤ m s − limit com-bined in kyr bins. Top panel gives the number of solutionsfor all possible combinations of clones (red histogram). The grayline gives | δ Ω | of the nominal orbit of (219902) 2002 EG and(432271) 2009 SH (see also the right ordinate). The green lineat the bottom panel shows the maximum difference in longitudeof ascending node between the clones of (219902) 2002 EG andthe longitude of ascending node of its nominal orbit, comparedwith the same information for (258656) 2002 ES given in Fig. 9). ( . ± . ) × − au, e P ≃ ( . ± . ) × − and sin I P ≃ ( . ± . )× − for (432271) 2009 SH (forreference, we again mention their quite small libration ampli-tudes . ◦ , resp. . ◦ ). This is a configuration reminiscentof the (258656) 2002 ES -2013 CC case, though each ofthe three proper elements is slightly larger now. The rela-tive velocity δ V P is again entirely dominated by the properinclination difference, this time somewhat smaller than inthe (258656) 2002 ES -2013 CC case (only ≃ . ◦ ).Assuming geometric albedo value . , we obtain sizes of ≃ . km and ≃ ( . − . ) km for (219902) 2002 EG and (432271) 2009 SH , considering absolute magnitudevalues from the major three small-body ephemerides sites asabove. While little larger, it still places this couple into thesame category of very small Trojans as (258656) 2002 ES -2013 CC .We repeated the convergence experiment using geo-metrical clones from section 3.3.1. In particular we consid-ered nominal (best-fit) orbits of (219902) 2002 EG and(432271) 2009 SH , and for each of them we constructed 20geometrical clone variants of the initial data at MJD58800 epoch. We again used information from the AstDyS web-site and noted that both initial orbits of components in thispossible couple have smaller uncertainties in all orbital ele-ments than the orbits of (258656) 2002 ES and 2013 CC .This is because their longer observation arcs and more dataavailable for the orbit determination. We propagated these42 (21+21) test bodies backward in time to . Gyr beforepresent. Perturbations from all planets were included andevery yr configuration of the nominal orbits and ac-companied clones for the two bodies compared. A criterionfor convergence included δ V ≤ m s − from Eq. (3), andsmall eccentricity and inclination differences. In particular,we required δ e ≤ − and δ I ≤ . ◦ . These values areonly slightly larger than the difference in the correspond-ing proper values and and each represent a few meters persecond contribution in (1).Results are shown in Fig. B1 which has the samestructure as the Fig. 9, previously given for the (258656)2002 ES and 2013 CC couple. The main take-away mes-sage is in the bottom panel, which shows maximum nodaldifference between clones of (219902) 2002 EG and itsnominal orbit as a function of time to the past. The slope ofthe initially linear trend (lasting approximately Myr) issimply given by maximum δ s proper frequency among clonesfrom the initial data difference. The non-linearity, which de-velops at later epochs, is due to orbital long-term chaoticity.While for the (258656) 2002 ES and 2013 CC couplethe chaotic effects were very minimum, the nodal differencebetween (258656) 2002 ES clones and the nominal orbitincreased to only ≃ ◦ in . Gyr. At the end of our run thenodal difference expanded to ≃ ◦ . Given the very limitednumber of clones we had, this works again identification ofconvergent solutions. Note that beyond ≃ Myr, wherewe would expect more convergent cases, we could satisfy theconvergence criteria of only few meters per second describedabove only rarely. CPU-demanding effort with many moreclones would be needed to achieve the desired convergencelimits.We repeated the same experiment for several other can-didate couples from the small- δ V P sample, including the caseof (215110) 1997 NO –2011 PU (see Fig. 2), but observedeven faster onset of the clone diffusion in the Trojan orbitalphase space. This was due to their large e P and/or sin I P val-ues, as well as larger libration amplitudes. Their systematicanalysis is beyond the scope of this paper. This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS000