Asteroid families in the first order resonances with Jupiter
aa r X i v : . [ a s t r o - ph . E P ] A p r Mon. Not. R. Astron. Soc. , 1–20 (2008) Printed 27 March 2018 (MN L A TEX style file v2.2)
Asteroid families in the first order resonances with Jupiter
M. Broˇz ⋆ and D. Vokrouhlick´y Institute of Astronomy, Charles University, Prague, V Holeˇsoviˇck´ach 2, 18000 Prague 8, Czech Republic
Accepted ???. Received ???; in original form ???
ABSTRACT
Asteroids residing in the first-order mean motion resonances with Jupiter hold impor-tant information about the processes that set the final architecture of giant planets.Here we revise current populations of objects in the J2/1 (Hecuba-gap group), J3/2(Hilda group) and J4/3 (Thule group) resonances. The number of multi-opposition as-teroids found is 274 for J2/1, 1197 for J3/2 and 3 for J4/3. By discovering a second andthird object in the J4/3 resonance, (186024) 2001 QG and (185290) 2006 UB ,this population becomes a real group rather than a single object. Using both hierarchi-cal clustering technique and colour identification we characterise a collisionally-bornasteroid family around the largest object (1911) Schubart in the J3/2 resonance. Thereis also a looser cluster around the largest asteroid (153) Hilda. Using N -body numeri-cal simulations we prove that the Yarkovsky effect (infrared thermal emission from thesurface of asteroids) causes a systematic drift in eccentricity for resonant asteroids,while their semimajor axis is almost fixed due to the strong coupling with Jupiter.This is a different mechanism from main belt families, where the Yarkovsky drift af-fects basically the semimajor axis. We use the eccentricity evolution to determine thefollowing ages: (1 . ± .
7) Gyr for the Schubart family and & ≃ Key words: celestial mechanics – minor planets, asteroids – methods: N -body sim-ulations. Populations of asteroids in the Jovian first order mean mo-tion resonances –J2/1, J3/2 and J4/3– are closely linked tothe orbital evolution of the giant planets. This is because oftheir orbital proximity to Jupiter. Stability or instability ofthese asteroid populations directly derives from the orbitalconfiguration of the giant planets. As such it is also sensitiveon the nature and amount of Jupiter’s migration and otherfiner details of its dynamics. As a result, the currently ob-served asteroids in the Jovian first order resonances contain ⋆ E-mail: [email protected]ff.cuni.cz (MB);[email protected] (DV). Interestingly, at their discovery (158) Hilda and (279) Thule,residing in the J3/2 and J4/3 resonances, immediately attractedattention of astronomers by vastly extending asteroid zone towardgiant planets and by their ability to apparently approach Jupiternear aphelia of their orbits (e.g., K¨uhnert 1876; Krueger 1889). valuable information about the early evolution of planetsand, if correctly understood and properly modelled, theymay help to constrain it.Apart from the Trojan clouds (not studied in this paper)the largest known population in the Jovian mean motionresonances occupies the J3/2 resonance, and is frequentlycalled the Hilda group. It was carefully studied in a par-allel series of papers by Schubart and Dahlgren and col-laborators during the past few decades. Schubart (1982a,b,1991) analysed short-term dynamics of Hilda-type orbitsand introduced quasi-constant orbital parameters that al-lowed their first classification. While pioneering, Schubart’swork had the disadvantage of having much smaller sample ofknown asteroids and computer power than today. Dahlgren& Lagerkvist (1995) and Dahlgren et al. (1997, 1998, 1999)conducted the first systematic spectroscopic and rotation-rate investigation of Hildas. They found about equal abun- c (cid:13) M. Broˇz and D. Vokrouhlick´y dance of D- and P-type asteroids and suggested spectral-size correlation such that P-types dominate large Hildas andD-type dominate smaller Hildas. They also suggested smallHildas have large lightcurve amplitudes, as an indication ofelongated or irregular shape, and that the distribution oftheir rotation rates is non-Maxwellian. Further analysis us-ing the Sloan Digital Sky Survey (SDSS) data however doesnot support significant dominance of either of the two spec-tral types for small sizes and indicates about equal mix ofthem (Gil-Hutton & Brunini 2008; see also below). Smallerpopulations of asteroids in the J2/1 and J4/3 received com-paratively less observational effort.Since the late 1990s powerful-enough computers alloweda more systematic analysis of fine details of the longer-termdynamics in the Jovian first order resonances. Ferraz-Mello& Michtchenko (1996) and Ferraz-Mello et al. (1998a,b) de-termined that asteroids in the J2/1 resonance can be verylong-lived, possibly primordial, yet their motion is compar-atively more chaotic than those in the J3/2 resonance. Thelatter paper showed that commensurability between the li-bration period and the period of Jupiter’s and Saturn’sGreat Inequality might have played an important role indepletion of the J2/1 resonance. This would have occurredwhen both giant planets were farther from their mutual 2:5mean motion configuration in the past. A still more completeanalysis was obtained by Nesvorn´y & Ferraz-Mello (1997)who also pointed out that the J4/3 resonance stable zone issurprisingly void of asteroids, containing only (279) Thule.Roig et al. (2001) and Broˇz et al. (2005) recently revised thepopulation of asteroids in the J2/1 resonance and classifiedthem into several groups according to their long-term orbitalstability. While the origin of the unstable resonant popula-tion was successfully interpreted using a model of a steady-state flow of main belt objects driven by the Yarkovsky semi-major axis drift, the origin of the long-lived asteroids in theJ2/1 remains elusive. Population of Hildas and Thule wasassumed primordial or captured by an adiabatic migrationof Jupiter (e.g., Franklin et al. 2004).It has been known for some time that the current con-figuration of giant planets does not correspond to that attheir birth. However, a new momentum to that hypothesiswas given by the so called Nice model (Tsiganis et al. 2005;Morbidelli et al. 2005; Gomes et al. 2005). The Nice modelpostulates the initial configuration of the giant planets wassuch that Jupiter and Saturn were interior of their mutual1:2 mean motion resonance (see also Morbidelli et al. 2007).The event of crossing this resonance had a major influenceon the final architecture of giant planets and strongly influ-enced structure of small-bodies populations in the Solar sys-tem. Morbidelli et al. (2005) showed that the population ofJupiters Trojan asteroids was destabilised and re-populatedduring this phase. In what follows we show that, within theNice model, the same most probably occurs for populationsof asteroids in the J3/2 and J4/3 resonances.The paper is organised as follows: in Section 2 we re-vise information about the current populations of asteroidsin the Jovian first order resonances. We use an up-to-date AstOrb database of asteroid orbits from the Lowell Observa- Note the former P-type objects were reclassified to X-type in anewer taxonomy by Bus and Binzel (e.g., Bus et al. 2002). tory ( ftp.lowell.edu ) as of September 2007 and eliminateonly single-opposition cases to assure accurate orbital infor-mation.In Section 3 we apply clustering techniques and extracttwo families of asteroids on similar orbits in the J3/2 reso-nance. We strengthen their case with an additional colouranalysis using the SDSS broadband data. We model thelong-term orbital evolution of these families and estimatetheir ages on the basis of Yarkovsky-driven dispersion in ec-centricity.In Section 4 we determine an orbital stability of the pu-tative primordial populations of planetesimals in the Jovianfirst-order resonances. We show that those in the J3/2 andJ4/3 are very efficiently eliminated when Jupiter and Sat-urn cross their mutual 1:2 mean motion resonance. We alsodetermine the removal rate of very small resonant asteroidsdue to the Yarkovsky/YORP effects.
Dynamics of asteroid motion in the Jovian first order reso-nances has been extensively studied by both analytical andnumerical methods in the past few decades (e.g., Murray1986; Sessin & Bressane 1988; Ferraz-Mello 1988; Lemaitre& Henrard 1990; Morbidelli & Moons 1993; Moons et al.1998; Nesvorn´y & Ferraz-Mello 1997; Roig et al. 2002;Schubart 2007). In what follows we review a minimum in-formation needed to understand our paper, referring an in-terested reader to the literature mentioned above for moreinsights.In the simplest framework of a circular restricted planarthree-body problem (Sun-Jupiter-asteroid) the fundamentaleffects of the resonant dynamics is reduced to a one-degreeof freedom problem defined by a pair of variables (Σ , σ ). ForJ( p + 1) /p resonance ( p = 1 , √ a (cid:16) − p − e (cid:17) , (1) σ = ( p + 1) λ ′ − p λ − ̟ , (2)where a is the semimajor axis, e the eccentricity, ̟ the lon-gitude of pericentre and λ the mean longitude in orbit of theasteroid, and λ ′ is the mean longitude in orbit of Jupiter.If the asteroid motion is not confined into the orbitalplane of the planet, we have an additional pair of resonantvariables (Σ z , σ z ) such thatΣ z = 2 p a (1 − e ) sin i , (3) σ z = ( p + 1) λ ′ − p λ − Ω , (4)where i denotes the inclination of asteroids orbit and Ω thelongitude of its node. Remaining still with the simple av-eraged model, orbital effects with shorter periods are ne-glected, the motion obeys an integral of motion N given by N = √ a (cid:18) p + 1 p − p − e cos i (cid:19) . (5)Because of this integral of motion, variations of Σ implyoscillations of both a and e .The two-degree of freedom character of the resonant c (cid:13) , 1–20 steroid families in the first order resonances motion prevents integrability. However, as an approxima-tion we may introduce a hierarchy by noting that pertur-bation described by the (Σ , σ ) variables is larger than thatdescribed by the (Σ z , σ z ) terms (e.g., Moons et al. 1998).This is usually true for real resonant asteroids of interest.Only the angle σ librates and σ z circulates with a very longperiod. The (Σ z , σ z ) dynamics thus produces a long-periodperturbation of the (Σ , σ ) motion.Within this model the minimum value of Σ in one reso-nant cycle (typically several hundreds of years) implies a isminimum and e is maximum. These values do not conserveexactly from one cycle to another because the (Σ z , σ z ) mo-tion produces small oscillations. Since Σ + Σ z − N = −√ a/p one needs to wait until Σ z reaches maximum over its cycleto attain ‘real’ minimum of a values and ‘real’ maximumof e values over a longer time interval. From (3) we notethe maximum of Σ z occurs for the maximum of i variations.This situations occurs typically once in a few thousands ofyears. In an ideal situation, these extremal values of ( a, e, i )would be constant and may serve as a set of proper orbitalelements.The motion of real asteroids in the Solar system is fur-ther complicated by Jupiter having non-zero and oscillat-ing value of eccentricity. This brings further perturbations(e.g., Ferraz-Mello 1988; Sessin & Bressane 1988 for a sim-ple analytic description) and sources of instability insidethe resonance. Despite the non-integrability, we follow Roiget al. (2001) and introduce pseudo-proper orbital elements ( a p , e p , sin i p ) as the osculating elements ( a, e, sin i ) at themoment, when the orbit satisfies the condition: σ = 0 ∧ dσdt < ∧ ̟ − ̟ ′ = 0 ∧ Ω − Ω ′ = 0 , (6)where ̟ ′ and Ω ′ denote the longitude of pericentre and thelongitude of node of Jupiter. As above, when (6) holds theosculating orbital elements are such that a attains minimum, e attains maximum and i attains maximum. Numerical ex-periments show that with a complete perturbation modeland a finite time-step it is difficult to satisfy all conditionsof (6) simultaneously. Following Roig et al. (2001) we thusrelax (6) to a more practical condition | σ | < ◦ ∧ ∆ σ ∆ t < ∧ | ̟ − ̟ ′ | < ◦ . (7)Because this condition is only approximate, we numericallyintegrate orbits of resonant asteroids for 1 Myr, over whichthe pseudo-proper orbital elements are recorded. We thencompute their mean value and standard deviation, which isan expression of the orbital stability over that interval oftime.In the case of the J3/2 and J4/3 resonances, we usethe condition (7) with a different sign ∆ σ ∆ t > σ ( t ). This intermediate stage serves to suppress oscillationsfaster than the libration period. The different sign of ∆ σ ∆ t just means our pseudo-proper orbital elements correspondto maximum value of a and minimum values of e and i , inorder to allow more direct comparison with previous analy-ses. Aside to this short-term integration we perform long-term runs to determine the stability of a particular reso-nant orbit. With this aim we conduct integrations spanning 4 Gyr for all resonant asteroids. Because of the inherent un-certainty in the initial conditions (orbital elements at thecurrent epoch), we perform such integration for the nomi-nal orbit and 10 clones that randomly span the uncertaintyellipsoid. We then define dynamical lifetime of the orbit asthe median of time intervals, for which the individual clonesstayed in the resonance.All integrations are performed using the Swift pack-age (Levison & Duncan 1994), slightly modified to includenecessary on-line digital filters and a second order symplec-tic integrator (Laskar & Robutel 2001). Most of numeri-cal simulations take into account gravitational interactionsonly, but in specific cases – and when explicitly mentioned– we include also Yarkovsky (thermal) accelerations. In thiscase we use an implementation described in detail by Broˇz(2006). Our simulations include 4 outer planets. We mod-ify the initial conditions of the planets and asteroids by abarycentric correction to partially account for the influenceof the terrestrial planets. The absence of the terrestrial plan-ets as perturbers is a reasonable approximation in the outerpart of the main belt and for orbits with e < . In order to determine, which objects are located in the J2/1mean motion resonance, we first extracted orbits from the
AstOrb database with osculating orbital elements in a broadbox around this resonance (see, e.g., Roig et al. 2001 for asimilar procedure). We obtained 7139 orbits, which we nu-merically integrated for 10 kyr. We recorded and analysedbehaviour of the resonance angle σ = 2 λ ′ − λ − ̟ fromEq. (2). Pericentric librators, for which σ oscillates about0 ◦ , were searched. We found 274 such cases; this extendsthe previous catalogue of Broˇz et al. (2005) almost twice.The newly identified resonant objects are mainly asteroidsdiscovered or recovered after 2005 with accurate enough or-bits. We disregard from our analysis asteroids at the borderof the resonance, for which σ ( t ) exhibits alternating peri-ods of libration and circulation, and also those asteroids forwhich σ oscillates but are not resonant anyway ( N . Among the short-lived objects we found 14 have dynamical lifetimes even less Our results for both J2/1 and J3/2 resonancesare summarised in tables available through a website http://sirrah.troja.mff.cuni.cz/yarko-site/ (those forJ4/3 bodies are given in Table 1). These contain listing of allresonant asteroids, their pseudo-proper orbital elements withc (cid:13) , 1–20
M. Broˇz and D. Vokrouhlick´y nu m be r o f a s t e r o i d s i n b i n N median lifetime in J2/1 t J2/1 (Myr)
122 Zhongguos ex t r e m e l y un s t a b l e unstable Griquas and Zhongguos
Figure 1.
A distribution of the median dynamical lifetimes forobjects in the J2/1 resonance. A division to several groups isdenoted: extremely unstable objects ( t J2 / t J2 /
70 Myr) and long-lived objects (Griquas andZhongguos, t J2 / >
70 Myr). than 2 Myr and we call them extremely unstable. Broˇz et al.(2005) suggested the unstable orbits in the J2/1 resonanceare resupplied from the adjacent main belt due to a per-manent flux driven by the Yarkovsky force, the extremelyunstable objects are most probably temporarily capturedJupiter-family comets. The origin of the long-lived popula-tion in this resonance is still not known.Figure 2 shows the pseudo-proper orbital elements ofthe J2/1 asteroids projected onto the ( a p , e p ) and ( a p , sin i p )planes. Our data confirm that the unstable population ofJ2/1 asteroids populates the resonance outskirts near itsseparatrix, where several secular resonances overlap andtrigger chaotic motion (e.g., Morbidelli & Moons 1993;Nesvorn´y & Ferraz-Mello 1997; Moons et al. 1998). At low-eccentricities the chaos is also caused by an overlap onnumerous secondary resonances (e.g., Lemaitre & Henrard1990). Two ‘islands’ of stability –A and B– harbour the long-lived population of bodies. The high-inclination island A,separated from the low-inclination island B by the ν sec-ular resonance, is much less populated. Our current searchidentifies 9 asteroids in the island A. The origin of the asym-metry in A/B islands is not known, but since the work ofMichtchenko & Ferraz-Mello (1997) and Ferraz-Mello et al.(1998a,b) it is suspected to be caused by instability due tothe libration period commensurability with the forcing termsproduced by the Great Inequality.The size-frequency distribution of objects of a popula-tion is an important property, complementing that of theorbital distribution. Figure 3 shows cumulative distribution N (
Pseudo-proper orbital elements for the 247 objectsin the J2/1 resonance projected onto the planes of semimajoraxis a p vs eccentricity e p (top) and semimajor axis a p vs sine ofinclination sin i p (bottom). Bars are standard deviations of theelements derived from 1 Myr numerical integration. Position ofseveral secular resonances embedded in J2/1 is shown in the upperpanel. The unstable population of asteroids (crosses) occupiesthe region of their overlap; the stable population (full circles)occupies two distinct zones –A and B– of low-eccentricity andlow-inclination orbits (e.g., Nesvorn´y & Ferraz-Mello 1997). Thepopulation of marginally stable asteroids (open squares) residesin region adjacent to the unstable borders of the resonance ornear the bridge over the stable regions associated with the ν secular resonance. ulation is steeper than it would correspond to a standardcollisionally evolved system (e.g., Dohnanyi 1969; O’Brien& Greenberg 2003) with γ = +0 .
5. The same result holdsfor both the short- and long-lived sub-populations in thisresonance separately.Albedos of J2/1 bodies are not known, except for(1362) Griqua for which Tedesco et al. (2002) give p V =0 . p V = 0 .
05. For sake of simplicity we convert absolutemagnitudes to sizes using this averaged value when needed.For instance in Fig. 4 we show a zoom on the long-lived pop-ulation of objects in the J2/1 resonance with symbol size c (cid:13) , 1–20 steroid families in the first order resonances nu m be r o f a s t e r o i d s N ( < H ) absolute magnitude H (mag)diameter D (km) for albedo p V = 0.05 +0.35 +0.70 J3/2J2/1J4/3
Figure 3.
Cumulative distributions N ( 15) mag for J2/1 and (10 . , . 5) mag forJ3/2; no such approximation is available for J4/3 where only threeobjects are currently known. The H values where the straight lineapproximations level off from the data roughly correspond to thecompleteness limit of the population (R. Jedicke, personal com-munication). For sake of a rough comparison, the upper abscissagives an estimate of sizes for the albedo value p V = 0 . 05, averageof the outer belt population. weighted by the estimated size of the body. We note largeobjects are located far from each other and they are quiteisolated — no small asteroids are in close surroundings. Boththese observations suggest that the long-lived J2/1 popula-tion does not contain recently-born collisional clusters. Because asteroids in the J3/2 constitute a rather isolatedgroup, it is easy to select their candidates: we simply ex-tracted from the AstOrb database those asteroids with semi-major axis in between 3 . . σ = 3 λ ′ − λ − ̟ . We obtained 1197cases for which σ librates about 0 ◦ and which have N > . p s eudo - p r ope r e cc en t r i c i t y e p pseudo-proper semimajor axis a p (AU) AB Griqua Zhongguo Kohman p s eudo - p r ope r i n c li na t i on s i n ( I p ) pseudo-proper semimajor axis a p (AU) AB Zhongguo Kohman Figure 4. A zoom on the ( a p , e p ) and ( a p , sin i p ) plots from Fig. 2with relative size of the resonant asteroids indicated by size of thecrosses. Note the large bodies, some of which are labelled, residefar from each other. largest Hilda members might be the corresponding diffusivemechanisms. Small enough members might be also suscep-tible to the resonant Yarkovsky effect (see Sec. 4.2 and theAppendix A).Data in Fig. 3 confirm earlier findings that the Hildagroup is characterized by an anomalously shallow size distri-bution. In between absolute magnitudes H = 10 . . N ( 02) only.Sub-populations among Hilda asteroids, namely twocollisional families, are studied in Section 3. In spite of a frequent terminology “Thule group”, asteroidsin the J4/3 resonance consisted of a single object (279) Thuleup to now. Nesvorn´y & Ferraz-Mello (1997) considered thissituation anomalous because the extent of the stable zone ofthis resonance is not much smaller than that of the J3/2 res-onance (see also Franklin et al. 2004). In the same way, ourknowledge about the low- e and low- i Thule-type stable or- c (cid:13) , 1–20 M. Broˇz and D. Vokrouhlick´y p s eudo - p r ope r e cc en t r i c i t y e p pseudo-proper semimajor axis a p (AU) 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 3.96 3.97 3.98 3.99 4 4.01 4.02 4.03 4.04 p s eudo - p r ope r i n c li na t i on s i n ( I p ) pseudo-proper semimajor axis a p (AU)(153) Hilda(1911) Schubart Figure 5. Pseudo-proper orbital elements for 1197 Hildas pro-jected onto the planes of semimajor axis a p vs eccentricity e p (top) and semimajor axis a p vs sine of inclination sin i p (bottom).Larger size of the symbol indicates larger physical size of the as-teroid. Because of Hildas orbital stability, the uncertainty in thepseudo-proper element values is typically smaller than the sym-bol size. Note a tight cluster around the proper inclination valuesin i p ≃ . i p ≃ . bits (e.g., a = 4 . 27 AU, e = 0 . i = 5 ◦ ) should be obser-vationally complete at about magnitudes H = 12 . H beyond this complete-ness limit the Thule population should be known at ∼ H = 9–13 mag are very likely missing in this reso-nance. Where does the existing population of small Thule-type asteroids begin?Our initial search in the broad box around the J4/3resonance detected only 13 objects. Six of them, includ-ing the well-known extinct comet (3552) Don Quixote (e.g.,Weissman et al. 2002), are on typical orbits of Jupiter-family e cc en t r i c i t y e [] semimajor axis a [AU] stableescaped Figure 6. Pseudo-proper semimajor axis vs eccentricity plot forasteroids in the J3/2 resonance. Bold crosses denote orbits, whichescaped during the 4 Gyr of evolution. Solid line is the positionof the libration centre (the outer separatrix is located further tothe right). nu m be r o f a s t e r o i d s N albedo p V (153) Hilda(1911) Schubart (279) Thule Hildas in J3/2 Figure 7. The histogram shows distribution of values of the ge-ometric albedo p V determined by Tedesco et al. (2002) for aster-oids located inside the J3/2 resonance. Three individual values –(153) Hilda, (1911) Schubart and (279) Thule (a J4/3 asteroid) –are indicated on top. comets that happen to reside near this resonance with veryhigh eccentricity and moderately high inclination. Two moreare single-opposition objects and one has only poorly con-strained orbit, leaving us with (279) Thule and three ad-ditional candidate objects: (52007) 2002 EQ47, (186024)2001 QG and (185290) 2006 UB .Figure 8 (top panels) shows short-term tracks of (279)Thule, (186024) 2001 QG and (185290) 2006 UB in res-onant variables √ 2Σ (cos σ, sin σ ) ≃ ( e cos σ, e sin σ ) of theJ4/3 resonance ( σ = 4 λ ′ − λ − ̟ in this case). In all casestheir orbits librate about the pericentric branch ( σ = 0 ◦ ) ofthis resonance, although this is complicated –mainly in thelow-eccentricity case of (279) Thule– by the forced terms dueto Jupiter’s eccentricity (see, e.g., Ferraz-Mello 1988; Sessin& Bressane 1988; Tsuchida 1990). The leftmost panel re-covers the 40 ◦ –50 ◦ libration of (279) Thule, determined pre-viously by Tsuchida (1990; Fig. 3). The other two smallerasteroids show librations with comparable amplitudes. Thelast object, (52007) 2002 EQ47, appears to reside on an un-stable orbit outside the J4/3 resonance. Our search thus c (cid:13)000 2Σ (cos σ, sin σ ) ≃ ( e cos σ, e sin σ ) of theJ4/3 resonance ( σ = 4 λ ′ − λ − ̟ in this case). In all casestheir orbits librate about the pericentric branch ( σ = 0 ◦ ) ofthis resonance, although this is complicated –mainly in thelow-eccentricity case of (279) Thule– by the forced terms dueto Jupiter’s eccentricity (see, e.g., Ferraz-Mello 1988; Sessin& Bressane 1988; Tsuchida 1990). The leftmost panel re-covers the 40 ◦ –50 ◦ libration of (279) Thule, determined pre-viously by Tsuchida (1990; Fig. 3). The other two smallerasteroids show librations with comparable amplitudes. Thelast object, (52007) 2002 EQ47, appears to reside on an un-stable orbit outside the J4/3 resonance. Our search thus c (cid:13)000 , 1–20 steroid families in the first order resonances Table 1. Data on presently known population of asteroids residing in the J4/3 Jovian mean motion reso-nance (Thule group). Pseudo-proper orbital elements ( a p , e p , sin i p ) are given together with their standarddeviations ( δa p , δe p , δ sin i p ) determined from a 1 Myr numerical integration. σ p, max is the maximum libra-tion amplitude in the Sessin’s ( K, H ) variables (see Fig. 8), H is the absolute magnitude from the AstOrb catalogue and D is the estimated size using p V = 0 . 04 geometric albedo (Tedesco et al. 2002).No. Name a p e p sin i p δa p δe p δ sin i p σ p, max H D [AU] [AU] [deg] [mag] [km]279 Thule 4.2855 0.119 0.024 0.0005 0.012 0.003 ∼ 50 8 . 57 126 . . 36 8 . . 75 11 . e s i n σ e cos σ e s i n σ e cos σ e s i n σ e cos σ e s i n σ e cos σ e s i n σ e cos σ e s i n σ e cos σ Figure 8. Top panels: orbits of (279) Thule [left], (186024) 2001 QG [middle] and (185290) 2006 UB [right] in resonant variables( e cos σ, e sin σ ) of the J4/3 resonance. Here e is the eccentricity and σ = 4 λ ′ − λ − ̟ , with λ and λ ′ the mean longitude in orbit of theasteroid and Jupiter and ̟ is the longitude of asteroid’s pericentre (in all cases the osculating orbital elements are used). Each panelshows results of a short-term, 10 kyr numerical integration. In all three cases the orbits librate about the pericentric branch ( σ = 0 ◦ ) ofthe resonance. Perturbations due to Jupiter’s eccentricity and its variations make the regular libration move in an epicyclic manner (e.g.,Ferraz-Mello 1988; Sessin & Brassane 1988). Bottom panels: filtered resonant variables ( e cos σ, e sin σ ), with short-period oscillationsforced by Jupiter removed by digital filtering. They are similar to Sessin’s ( K, H ) coordinates in Sessin & Brassane (1988) or Tsuchida(1990). In each of the cases the pericentric libration is clearly revealed. The maximum libration amplitude in these coordinates is denotedby σ p, max in Table 1. lead to the detection of two new asteroids in this resonance,increasing its population by a factor of three. Results of a long-term numerical integration of the nom-inal orbits plus 10 close clones, placed within an orbital un-certainty, reveal that the orbit of (279) Thule is stable over4 Gyr, but the orbits of (186024) 2001 QG and (185290)2006 UB are partially unstable. They are not ‘short-lived’ but 45 % and 60 % of clones, respectively, escaped be-fore 4 Gyr. Figure 9 shows pseudo-proper semimajor axis vstime for nominal orbits and their clones of all J4/3 objects;the escaping orbits leave the figure before the simulation While our initial search used AstOrb catalogue from Septem-ber 2007, we repeated it using the catalogue as of June 2008. Noadditional J4/3 objects were found. was ended at 4 Gyr. We suspect similar non-conservative ef-fects as mentioned above for the 20% fraction of long-term-unstable Hildas to bring these two Thule members onto theirmarginally stable orbits.We would like to point out that it took more than a cen-tury from the discovery of (279) Thule (Palisa 1888; Krueger1889) until further objects in this resonance were finally dis-covered. This is because there is an anomalously large gapin size of these bodies: (279) Thule is 127 km in size with p V = 0 . 04 AU (e.g., Tedesco et al. 2002), while the estimatedsizes of (186024) 2001 QG and (185290) 2006 UB forthe same value of albedo are 8 . . H = 13–15 mag absolute magnituderange using future-generation survey projects such as Pan-STARRS (e.g., Jedicke et al. 2007). Such a completed pop- c (cid:13) , 1–20 M. Broˇz and D. Vokrouhlick´y 0 1000 2000 3000 4000 4.28 4.29 4.3 a p [ A U ] t [Myr](279) Thule 2001 QG Figure 9. Pseudo-proper semimajor axis a p vs time t plot forthe 3 orbits of asteroids located inside the J4/3 resonance and10 close clones orbits (placed randomly within the uncertaintyellipsoids) for each of them. Thin lines denote the stable orbitsand thick lines unstable, which escaped from the J4/3 before thecompletion of the simulation at 4 Gyr. ulation may present an interesting constraint on the plan-etesimal size distribution 4 Gyr ago. Collisions and subsequent fragmentations are ubiquitousprocesses since planets formed in the Solar system. Becausethe characteristic dispersal velocities of the ejecta (as a ruleof thumb equal to the escape velocity of the parent body)are usually smaller than the orbital velocity, the resultingfragments initially reside on nearby orbits. If the orbitalchaoticity is not prohibitively large in the formation zone,we can recognise the outcome of such past fragmentations asdistinct clusters in the space of sufficiently stable orbital ele-ments. More than 30 collisional families are known and stud-ied in the main asteroid belt (e.g., Zappal`a et al. 2002) withimportant additions in the recent years (e.g., Nesvorn´y et al.2002; Nesvorn´y, Vokrouhlick´y & Bottke 2006). Similarly, col-lisional families have been found among the Trojan clouds ofJupiter (e.g., Milani 1993; Beaug´e & Roig 2001; Roig et al.2008), irregular satellites of Jupiter (e.g., Nesvorn´y et al.2003, 2004) and even trans-Neptunian objects (e.g., Brownet al. 2007). Mean motion resonances, other than the Trojanlibrators of Jupiter, are typically too chaotic to hold sta-ble asteroid populations, or the populations were too smallto enable search for families. The only remaining candidatepopulations are those in the Jovian first order resonances,with Hilda asteroids the most promising group. However,low expectations for an existence of collisional families likelyde-motivated systematic search. Note, the estimated intrin-sic collisional probability of Hilda asteroids is about a factor3 smaller than in the main asteroid belt (e.g., Dahlgren 1998;Dell’Oro et al. 2001) and the population is more than twoorders of magnitude smaller.In spite of the situation outlined above, Schubart(1982a, 1991) repeatedly noticed groups of Hilda-type as-teroids with very similar proper elements. For instance, inhis 1991 paper he lists 5 members of what we call Schubartfamily below and pointed out their nearly identical values of N m i n , N v cutoff [m/s] N min N max( N i ) Figure 10. The dependence of the minimum number of asteroids N min , to be considered a statistically significant cluster, on thecutoff velocity v cutoff (thick curve). The average number N andthe maximum number max( N i ) of asteroids, which are closer than v cutoff , is shown by dashed and thin curves. All quantities arevalid for the J3/2 population. The fact that max( N i ) is muchlarger than N and N min indicates a presence of a significantcluster (or clusters) among the Hilda group. the proper inclination. Already in his 1982 paper Schubartmentions a similarity of such clusters to Hirayama families,but later never got back to the topic to investigate this prob-lem with sufficient amount of data provided by the growingknowledge about the J3/2 population. Even a zero orderinspection of Fig. 5, in particular the bottom panel, impliesthe existence of two large clusters among the J3/2 popu-lation. In what follows we pay a closer analysis to both ofthem.We adopt an approach similar to the hierarchical-clustering method (HCM) frequently used for identificationof the asteroid families in the main belt (e.g., Zappal`a et al.1990, 1994, 2002). In the first step of our analysis, we com-pute the number of bodies N min which is assumed to con-stitute a statistically significant cluster for a given valueof the cutoff velocity v cutoff . We use a similar approach tothat of Beaug´e & Roig (2001): for all asteroids in the J3/2resonance we determine the number N i ( v cutoff ) of asteroidswhich are closer than v cutoff . Then we compute the averagevalue N = ¯ N i . According to Zappal`a et al. (1994), a clustermay be considered significant if N > N min = N + 2 √ N .The plots N ( v cutoff ) and the corresponding N min ( v cutoff ) forHilda population are shown in Fig. 10. We use a standardmetric ( d defined by Zappal`a et al., 1994), namely: δv = na p s (cid:18) δa p a p (cid:19) + 2 ( δe p ) + 2 ( δ sin i p ) , (8)where ( a p , e p , sin i p ) are 10 Myr averaged values of the res-onant pseudo-proper elements (we checked that our resultspractically do not depend on the width of this averaginginterval). Schubart lists 11 additional asteroids in the group on his web-site ∼ s24/hilda.htm ,but again he does not go into details of their putative collisionalorigin. c (cid:13)000 The dependence of the minimum number of asteroids N min , to be considered a statistically significant cluster, on thecutoff velocity v cutoff (thick curve). The average number N andthe maximum number max( N i ) of asteroids, which are closer than v cutoff , is shown by dashed and thin curves. All quantities arevalid for the J3/2 population. The fact that max( N i ) is muchlarger than N and N min indicates a presence of a significantcluster (or clusters) among the Hilda group. the proper inclination. Already in his 1982 paper Schubartmentions a similarity of such clusters to Hirayama families,but later never got back to the topic to investigate this prob-lem with sufficient amount of data provided by the growingknowledge about the J3/2 population. Even a zero orderinspection of Fig. 5, in particular the bottom panel, impliesthe existence of two large clusters among the J3/2 popu-lation. In what follows we pay a closer analysis to both ofthem.We adopt an approach similar to the hierarchical-clustering method (HCM) frequently used for identificationof the asteroid families in the main belt (e.g., Zappal`a et al.1990, 1994, 2002). In the first step of our analysis, we com-pute the number of bodies N min which is assumed to con-stitute a statistically significant cluster for a given valueof the cutoff velocity v cutoff . We use a similar approach tothat of Beaug´e & Roig (2001): for all asteroids in the J3/2resonance we determine the number N i ( v cutoff ) of asteroidswhich are closer than v cutoff . Then we compute the averagevalue N = ¯ N i . According to Zappal`a et al. (1994), a clustermay be considered significant if N > N min = N + 2 √ N .The plots N ( v cutoff ) and the corresponding N min ( v cutoff ) forHilda population are shown in Fig. 10. We use a standardmetric ( d defined by Zappal`a et al., 1994), namely: δv = na p s (cid:18) δa p a p (cid:19) + 2 ( δe p ) + 2 ( δ sin i p ) , (8)where ( a p , e p , sin i p ) are 10 Myr averaged values of the res-onant pseudo-proper elements (we checked that our resultspractically do not depend on the width of this averaginginterval). Schubart lists 11 additional asteroids in the group on his web-site ∼ s24/hilda.htm ,but again he does not go into details of their putative collisionalorigin. c (cid:13)000 , 1–20 steroid families in the first order resonances Figure 11. A stalactite diagram computed for the J3/2 popula-tion (Hildas). Two prominent groupings, the Schubart family andthe Hilda family, are indicated. Every group plotted here has atleast 5 or N min members, whichever is larger (see Fig. 10). Next, we construct a stalactite diagram for Hildas ina traditional way (e.g., Zappal`a et al. 1990): we start with(153) Hilda as the first central body and we find all bodiesassociated with it at v cutoff = 300 m / s, using a hierarchicalclustering method (HCM; Zappal`a et al. 1990, 1994). Thenwe select the asteroid with the lowest number (cataloguedesignation) from remaining (not associated) asteroids andrepeat the HCM association again and again, until no aster-oids are left. Then we repeat the whole procedure recursivelyfor all clusters detected at v cutoff = 300 m / s, but now for alower value, e.g., v cutoff = 299 m / s. We may continue until v cutoff = 0 m / s, but of course, for too low values of the cut-off velocity, no clusters can be detected and all asteroids aresingle. The resulting stalactite diagram at Fig. 11 is simplythe asteroid number (designation) vs v cutoff plot: a dot ata given place is plotted only if the asteroids belongs to acluster of at least max(5 , N min ( v cutoff )) bodies. We are notinterested in clusters with less than 5 members; they aremost probably random flukes.We can see two prominent clusters among Hildas: thefirst one around the asteroid (153) Hilda itself, and the sec-ond one around (1911) Schubart. In the remaining part ofthis Section we discuss each of them separately.The stalactite diagram constructed in the same way forZhongguos and Griquas is shown in Fig. 12. No groupingseem to be significant enough to be considered an impact-generated cluster. This is consistent with the discussion ofthe ( a p , e p , sin i p ) plots in Section 2.1. The Schubart family can be distinguished from the remain-ing population of Hildas on a large range of cutoff velocities:from 50 m/s to more than 100 m/s (Fig. 11). It merges withthe Hilda family at 200 m/s. For the purpose of our analysiswe selected v cutoff = 60 m / s as the nominal value. While thetotal number of Schubart family members is not too sensitiveto this cutoff value, we refrain from using too high v cutoff , forwhich we would expect and increasing number of interlop- Figure 12. A stalactite diagram computed for the long-lived J2/1population. There are no prominent groupings; 60 asteroids arenot associated with any others, even at v cutoff = 300 m / s. Everygroup plotted here has at least 5 members. nu m be r o f a s t e r o i d s N ( < H ) absolute magnitude H (mag)diameter D (km) for albedo p V = 0.044 J3/2HildaSchubart Figure 13. Cumulative distributions N ( 25) for Hilda and (12 . , . 75) for Schubartclusters. For sake of a rough comparison, the upper abscissa givesan estimate of the sizes for an albedo value p V = 0 . ers to be associated with the family, and the family wouldattain a rather peculiar shape in the ( a p , e p , sin i p ) space.Figure 13 shows the cumulative distribution of the ab-solute magnitudes for the Schubart family members, com-pared to the rest of the J3/2 population. Importantly, theslope γ = (+0 . ± . 02) of the N ( C/X D J3/2 nu m be r o f s pe c t r a Hilda spectral slope ( PC ) Schubart Figure 14. Distribution of spectral slopes (PC components ofthe 5 broad-band colours) of 153 asteroids in the J3/2 resonance(top), 21 Hilda family members (middle) and 4 Schubart familymembers (bottom). Data computed from the 3rd release of theSDSS catalogue of moving objects (ADR3; Ivezi´c et al. 2002).Note some objects have been observed multiple times by SDSSand the histograms show distribution of all observations (ratherthan averages for a given body). The slope values of Hildas rangefrom neutral to steep, with roughly two groups separated by thevalue PC = 0 . 3: (i) C/X-types with PC < . 3, and (ii) D-typeswith PC > . 3. Importantly, the Schubart family members arespectrally similar; the median value PC = 0 . 20 corresponds to aC- or X-type parent body. We also analysed the available SDSS catalogue of mov-ing objects (ADR3; Ivezi´c et al. 2002). We searched for theJ3/2 asteroids among the entries of this catalogue and com-puted the principal component PC of the spectrum in thevisible band. Note the PC value is an indicator of the spec-tral slope and allows thus to broadly distinguish principalspectral classes of asteroids (e.g., Bus et al. 2002). Figure 14shows our results. The top panel confirms the bimodal char-acter of the J3/2 population (see also Dahlgren et al. 1997,1998, 1999 and Gil-Hutton & Brunini 2008). More impor-tantly, though, the bottom panel indicates a spectral homo-geneity of the Schubart family, placing all members withinthe C/X taxonomy class branch. This finding strongly sup-ports collisional origin of the Schubart family.Tedesco et al. (2002) derive D = 80 km size for (1911)Schubart, corresponding to a very low albedo p V = 0 . D = 38 km size of (4230) vanden Bergh and exactly the same albedo; this asteroid isamong the five largest in the family. Assuming the samealbedo for all other family members, we can construct a size-frequency distribution (Fig. 15). The slope α ≃ ( − . ± . nu m be r o f a s t e r o i d s N ( > D ) diameter D / km − − − J3/2HildaSchubart Figure 15. Cumulative size distribution N ( >D ) for the wholeJ3/2 population (solid) and for members of the two suggestedcollisional clusters: Hilda (dashed) and Schubart (dotted). Weassumed the geometric albedo p V = 0 . 044 (and 0 . 025 in theSchubart-family case; see the text) for the conversion of absolutemagnitudes H to diameters D . Labels are the best-fitted values of α using N ( >D ) ∝ D α (straight lines) in the interval (11 , 20) kmfor the two families and (9 , 20) km for the J3/2 population. in the numerical simulations of disruptions (e.g., Durda et al.2007). If we sum the volumes of the observed members, we endup with a lower limit for the parent body size D PB = 110 km,provided there are no interlopers. We can also estimate thecontribution of small (unobserved) bodies using the follow-ing simple method: (i) we sum only the volumes of the ob-served bodies larger than an assumed completeness limit D complete = 10 km ( V complete = P i π D i ); (ii) we fit thecumulative size distribution by a power law (log N ( >D ) = α log[ D ] km + β ; α = − . β = 4 . 73 for the Schubart);(iii) we prolong this slope from D complete down to D min = 0and calculate the total volume of the parent body (provided α > − V PB = V complete + π β αα + 3 (cid:2) D a +3min − D a +3complete (cid:3) . (9)The result is D PB = q π V PB . = 130 km, some sort of an up-per limit. The volumetric ratio between the largest fragmentand the parent body is then V LF /V PB . = 0 . 2, a fairly typicalvalue for asteroid families in the main asteroid belt. Obvi-ously, the assumption of a single power-law extrapolationof the N ( >D ) at small sizes is only approximate and canlead to a result with a 10 % uncertainty. However, if we usean entirely different geometric method, developed by Tanga We also mention that so far asteroid disruption simulations didnot explore cases of weak-strength materials appropriate for thesuggested C/X spectral taxonomy of the Schubart family parentbody. c (cid:13) , 1–20 steroid families in the first order resonances et al. (1998), we obtain D PB ≃ d disrupt of a projectile nec-essary to disrupt the parent body of the Schubart family?Using Eq. (1) from Bottke et al. (2005): d disrupt = (cid:0) Q ∗ D /V (cid:1) D target . (10)and substituting Q ∗ D = 10 J / kg for the strength (somewhatlower than that of basaltic objects to accommodate the as-sumed C/X spectral type; e.g., Kenyon et al. 2008 and ref-erences therein), V imp = 4 . 78 km / s for the typical impactvelocity (see Dahlgren 1998) and D target ≃ 130 km, we ob-tain d disrupt ≃ 25 km. At this size the projectile populationis dominated by main belt bodies. Considering also differentintrinsic collisional probabilities between Hilda–Hilda aster-oids (2 . × − km − yr − ; Dahlgren 1998) and Hilda–mainbelt asteroids (0 . × − km − yr − ), we find it more likelythe Schubart family parent body was hit by a projectile orig-inating from the main belt. We repeated the same analysis as in Section 3.1 for theHilda family. The family remains statistically distinct fromthe whole J3/2 population in the range of cutoff velocities(130 , / s; we choose v cutoff = 150 m / s as the nominalvalue.The slope γ of the cumulative absolute magnitude dis-tribution N ( 02) (Fig. 13), again steeperthan for the total J3/2 population and comparable to thatof the Schubart family. The spectral slopes (PC ) are some-what spread from flat (C/X-compatible values; PC < . > . 3) — see Fig. 14.Overall, though, the C/X members prevail such that the D-type objects might be actually interlopers, at least accordingto a simple estimate based on the volume of the Hilda familyin the ( a p , e p , sin i p ) space, compared to the total volume ofthe J3/2 population.Tedesco et al. (2002) determine albedos for six familymembers. They range from 0.037 to 0.087, but three val-ues are close to the median albedo 0.044 of all J3/2 aster-oids. We thus consider this value to be representative of theHilda family. The corresponding cumulative size distributionis plotted in Fig. 15. Using the same method as in Section 3.1we estimate the size of the parent body D PB = 180–190 km,with V LF /V PB ≃ . 8. With the model of Tanga et al. (1999)we would obtain D PB ≃ 210 km and thus V LF /V PB ≃ . d disrupt = 50–55 km.While not so prominent as the Schubart family, we con-sider the group of asteroids around Hilda a fairly robust caseof a collisionally-born family too. In order to asses some limits for the age of the Schubart andHilda families, we perform a number of numerical tests. Inparticular, we simulate a disruption of a parent body insidethe resonance and numerically determine the long-term or-bital evolution of fragments. The evolved synthetic family e p a p / AU J3/2Schubart cluster initial osculating elementspseudo-proper elements Figure 16. The initial osculating elements of an impact-generated swarm of 139 fragments at the location of(1911) Schubart (bottom crosses), the corresponding pseudo-proper elements computed from the first My of evolution (uppercrosses) and the pseudo-proper elements of the observed Schubartfamily (circles). We show here projection onto the plane definedby semimajor axis and eccentricity. Dots are the pseudo-properelements of the background J3/2-population asteroids. The ini-tial synthetic swarm of asteroids poorly matches the observedfamily: it is both too extent in semimajor axis and too compactin eccentricity. e p a p / AU J3/2 Schubart cluster impact at f = 135 ° Figure 17. The same as Fig. 16, but for a different impact ge-ometry ( f = 135 ◦ , ω + f = 180 ◦ ). This choice of f maximizesthe initial spread of the synthetic family in proper eccentricity. Inthis case we do not show the initial osculating orbital elements. at different time steps is then compared with the observedfamily. Ideally, this approach should allow to constrain thetime elapsed since the family formed.As a first step, we need to create a synthetic family in-side the resonance. We use current orbital elements of thelargest family member, (1911) Schubart in this case, as rep-resentative to the parent body and only allow changes inthe true anomaly f and in the argument of pericentre ω atthe break-up event. By changing these two geometric pa- c (cid:13) , 1–20 M. Broˇz and D. Vokrouhlick´y rameters we can produce different initial positions of thefragments in the orbital element space. For sake of our test,fragments are assumed to be dispersed isotropically with re-spect to the parent body, with a velocity distribution givenby the model of Farinella et al. (1993, 1994). The numberof fragments dN ( v ) launched with relative velocities in theinterval ( v, v + dv ) is given by dN ( v ) = Cv ( v + v ) − ( κ +1) / dv , (11)with C a normalization constant, v esc the escape velocityfrom the parent body and κ = 3 . 25. To prevent excessiveescape velocities we introduce a maximum allowed value v max . Nominally, we set v max = 200 m/s, but in the nextSection 3.4 we also use restricted values of this parameterto test sensitivity of our results to initial conditions.To simulate an impact that might have created theSchubart family, we generated velocities randomly for139 fragments with v esc = 65 m / s (note the number of frag-ments in the synthetic family is equal to the number ofthe Schubart family members). The resulting swarm of frag-ments is shown in Fig. 16, for the impact geometry f = 0 ◦ and ω + f = 180 ◦ . We show both the initial osculating orbitalelements and the pseudo-proper elements.The synthetic family extends over significantly largerrange of the semimajor axis than the observed Schubartfamily, but all fragments still fall within the J3/2 resonance.The eccentricity distribution is, on the other hand, substan-tially more compact. Only the distribution of inclinations ofthe synthetic family roughly matches that of the observedfamily. We verified this holds also for other isotropic-impactgeometries (such as f = 135 ◦ and ω + f = 180 ◦ shown inFig. 17). The peculiar shape of the synthetic family in thepseudo-proper element space ( a p , e p ) is an outcome of theisotropic disruption, simply because some fragments fall tothe left from the libration centre of the J3/2 resonance (at3.97 AU) and they are ‘mapped’ to the right. This is becausethe pseudo-proper elements are the maxima and minima of a and e , respectively, over their resonant oscillations.The initial configuration of the synthetic family waspropagated for 4 Gyr, using the integrator described in Sec-tion 2. At this stage, we use only the gravitational perturba-tions from the 4 exterior giant planets. We performed suchsimulation for several impact geometries, as determined by f and ω , with similar results.Figure 18 shows the long-term evolution of the syn-thetic family. Because the family resides mostly in the sta-ble zone of the J3/2 resonance, only little evolution can beseen for most of the bodies. This is in accord with findingsof Nesvorn´y & Ferraz-Mello (1997) who concluded that thestable region in this resonance shows little or no diffusionover time scales comparable to the age of the Solar system.Only about 10 % of orbits that initially started at the out-skirts of the stable zone (with large libration amplitudes)escaped from the resonance during the 4 Gyr simulation.The removal of orbits with large semimajor axis a p helpsin part to reconcile the mismatch with the distribution ofthe observed Schubart family. However, the dispersion ineccentricity e p does not evolve much and it still shows largemismatch if compared to the observed family. Even in thecase f = 135 ◦ ( ω + f = 180 ◦ ; Fig. 17), which maximizes theinitial eccentricity dispersion of the synthetic fragments, thefinal value at 4 Gyr is about three times smaller than that Figure 18. The synthetic family from Fig. 16 evolved over 4 Gyr:the grey dots show evolutionary tracks of the fragments in thepseudo-proper orbital element space. Overall, stability of the J3/2resonance makes many fragments to stay very close to their ini-tial values. Only ∼ 10 % of fragments with the initial extremalvalues of a p (and thus the libration amplitude) escape from theresonance during the simulation. This helps in part to reduce themismatch with the observed family (circles) in semimajor axis,but is not sufficient to attain the Schubart-family full eccentricitydispersion. of the Schubart family. Clearly, our model is missing a keyelement to reproduce the current orbital configuration ofthis family.One possibility to resolve the problem could be to re-lease the assumption of an isotropic impact and exploreanisotropies in the initial velocity field. This is an obvioussuspect in all attempts to reconstruct orbital configurationsof the asteroid families, but we doubt it might help much inthis case. Exceedingly large relative velocities, compared tothe escape velocity of the estimated parent body, would berequired. Recall, the fragments located in the stable regionof the J3/2 resonance would hardly evolve over the age ofthe Solar system.A more radical solution is to complement the forcemodel, used for the long-term propagation, by additionaleffects. The only viable mechanism for the size range we aredealing with is the Yarkovsky effect. This tiny force, due toanisotropic thermal emission, has been proved to have deter-mining role in understanding fine structures of the asteroidfamilies in the main belt (e.g., Bottke et al. 2001; Vokrouh-lick´y et al. 2006a,b). In these applications the Yarkovskyeffect produces a steady drift of the semimajor axis, leavingother orbital elements basically constant.However, the situation is different for resonant orbits.The semimajor axis evolution is locked by the strong grav-itational influence of Jupiter. For that reason we first ransimplified simulations with the Yarkovsky forces — resultsof these tests are briefly described in the Appendix A. Wenext applied the model containing both gravitational andYarkovsky perturbations to the evolution of the syntheticfamily. Results of these experiments are described in thenext Section 3.4. c (cid:13) , 1–20 steroid families in the first order resonances Figure 19. The impact-generated swarm from Fig. 16 evolvedby planetary perturbations and the Yarkovsky forces, in the pro-jection on the pseudo-proper semimajor axis a p vs eccentricity e p plane. The grey dots indicate the evolutionary tracks over thewhole 4 Gyr time span and crosses denote the configuration at1.7 Gyr, when the eccentricity dispersion of the synthetic familyparticles roughly matches that of the observed Schubart family(circles). We ran our previous simulation of the long-term evolutionof the synthetic family with the Yarkovsky forces included.Our best guess of thermal parameters for bodies of the C/Xtype is: ρ s = ρ b = 1300 kg / m for the surface and bulk densi-ties, K = 0 . 01 W / m / K for the surface thermal conductivity, C = 680 J / kg for the heat capacity, A = 0 . 02 for the Bondalbedo and ǫ = 0 . 95 for the thermal emissivity parameter.Rotation periods are bound in the 2–12 hours range. Spinaxes orientations are assumed isotropic in space. Finally, weassign sizes to our test particles equal to the estimate of sizesfor Schubart family members, based on their reported abso-lute magnitudes and albedo p V = 0 . resonant Yarkovsky effect produces mainlysecular changes of eccentricity. We recall this systematicdrift in e must not be confused with the chaotic diffusionin e . We also note that inclination of the orbits remains sta-ble, in accord with a good match of the Schubart family bythe initial inclination distribution.Because the initial eccentricity dispersion of the syn-thetic family is much smaller than that of the observedone, its steady increase due to the combined effects of theYarkovsky forces and the resonant lock gives us a possibilityto date the origin of the family (see Vokrouhlick´y et al. 2006afor a similar method applied to families in the main belt).To proceed in a quantitative way, we use a 1-dimensional D KS ( K o l m ogo r o v - S m i r no v d i ff e r en c e ) time t / Myr p (> D KS ) = 0.01 v max = 200 m/s40, 60, 80 or 100 m/s Figure 20. The difference between the cumulative distributionof eccentricities for the observed Schubart family and our syn-thetic families expressed as the Kolmogorov–Smirnov (KS) dis-tance D KS vs time t . All models assume f = 0 ◦ and ω + f = 180 ◦ but they have different maximum values v max of the ejected par-ticles: v max = 40, 60, 80, 100 and 200 m/s (labels). For t < t > . p ( >D KS ) < . 01. Thus with a 99% probability level wecan exclude such age for the Schubart family. The best matchesare achieved for v max = 40 and 60 m/s, for which the most ex-treme particles have been eliminated from the synthetic family. Kolmogorov–Smirnov (KS) test to compare cumulative dis-tribution of pseudo-proper eccentricity values e p of the ob-served and synthetic families (e.g., Press et al. 2007).Figure 20 shows the KS distance D KS of the two ec-centricity distributions as a function of time. For sake ofa test, we also use smaller v max values of the initial veloc-ity field (essentially, this is like to start with a more com-pact synthetic family). Regardless of the v max value, ourmodel rejects Schubart family ages smaller than 1 Gyr andlarger than 2.4 Gyr with a 99 % confidence level. For agesin between 1.5 Gyr and 1.7 Gyr the KS-tested likelihood ofa similarity of the synthetic-family and the observed-family e p distributions can reach up to 50 %. We thus conclude themost likely age of the Schubart cluster is (1 . ± . 7) Gyr.We repeated the same analysis for Hilda family by cre-ating a synthetic family of 233 particles about (153) Hilda.In this case we used v esc = 110 m / s. The situation is actu-ally very similar to the Schubart — there is again a problemwith the small dispersion of eccentricities in case of a purelygravitational model. Using the model with the Yarkovskyeffect, we can eventually fit the spread of eccentricities andaccording to the KS test (Fig. 21) the age of the familymight be & c (cid:13) , 1–20 M. Broˇz and D. Vokrouhlick´y D KS ( K o l m ogo r o v - S m i r no v d i ff e r en c e ) time t / Myr p (> D KS ) = 0.01 v max = 200 m/s50, 100 or 150 m/s Figure 21. The same as Fig. 20 but for the Hilda family. TheKS distance D KS of the observed and modelled pseudo-propereccentricity e p distribution for several synthetic clusters with dif-ferent maximum velocities v max are plotted vs time t . In this case, t . family parent body is a very unlikely event during the last3.5 Gyr).We finally simulated a putative collision in the J2/1 res-onance, around the asteroid (3789) Zhongguo (the largestasteroid in the stable island B). There are two major differ-ences as compared to the J3/2 resonance: (i) the underlyingchaotic diffusion due to the gravitational perturbations islarger in the J2/1 resonance (e.g., Nesvorn´y & Ferraz-Mello1997), such that an initially compact cluster would fill thewhole stable region in 1–1.5 Gyr and consequently becomesunobservable; (ii) sizes of the observed asteroids are gener-ally smaller, which together with a slightly smaller helio-centric distance, accelerates the Yarkovsky drift in e p . Thelatter effect would likely shorten the time scale to 0.5–1 Gyr.Thus the non-existence of any significant orbital clusters inthe J2/1 resonance (Sec. 3 and Fig. 12) does not excludea collisional origin of the long-lived population by an eventolder than 1 Gyr. This would also solve the apparent prob-lem of the very steep size distribution of the stable J2/1population (Broˇz et al. 2005 and Fig. 3). Note the expectedcollisional lifetime of the smallest observed J2/1 asteroids isseveral Gyr (e.g., Bottke et al. 2005). We finally pay a brief attention to the overall stability of as-teroid populations in the first-order resonances with respectto different configurations of giant planets. We are focus-ing on the situations when the orbits of Jupiter and Saturnbecome resonant. This is motivated by currently adoptedviews about final stages of building planetary orbits archi-tecture, namely planet migration in a diluted planetesimaldisk (e.g., Malhotra 1995; Hahn & Malhotra 1999; Tsiga-nis et al. 2005; Morbidelli et al. 2005; Gomes et al. 2005 –these last three references are usually described as the Nicemodel). Morbidelli et al. (2005) proved that the primordial 385 390 395 400 405 410 415 9.525 9.526 9.527 9.528 9.529 9.53 pe r i od / y r semimajor axis of Saturn a S / AU Great Inequalitymean J2/1 libration(1362) Griqua Figure 22. Dependence of the period associated with circulationof the Great Inequality (GI) angle 5 λ S − λ J ( λ S and λ J aremean longitude in orbit for Saturn and Jupiter) on the osculatingsemimajor axis a S of Saturn (crosses); Jupiter’s orbit is fixed.Moving Saturn farther away from the location of 2:5 resonancewith Jupiter (at a S = 9 . 584 AU) makes the GI period shorter(with today’s a S = 9 . a S values shown in this figure the GI periodbecomes comparable to the average period ¯ P J2 / of librations inthe J2/1 resonance (solid line). Libration period for a particularcase of asteroid (1362) Griqua is shown by circles. Trojan asteroids were destabilized when Jupiter and Sat-urn crossed their mutual 1:2 mean motion resonance and,at the same time, the Trojan region was re-populated byparticles of the planetesimal disk. Since the mutual 1:2 res-onance of Jupiter and Saturn plays a central role in theNice model, and since these two planets had to cross other(weaker) mutual resonances such as 4:9 and 3:7 before theyacquired today’s orbits, one can naturally pose a questionabout the stability of primordial populations in the first-order mean motion resonances with Jupiter. Ferraz-Melloet al. (1998a,b) demonstrated that even subtler effects caninfluence the J2/1 population, namely the resonances be-tween the asteroid libration period in the J2/1 resonance andthe period of the Great Inequality (GI) terms in planetaryperturbations (i.e., those associated with Jupiter and Sat-urn proximity to their mutual 2:5 mean motion resonance;Fig. 22). A first glimpse to the stability of the first-orderresonance populations with respect to these effects is givenin Section 4.1.In Section 4.2 we also briefly estimate the change ofdynamical lifetimes for small J2/1 and J3/2 bodies causedby the Yarkovsky effect. In what follows we use a simple approach by only movingSaturn’s orbit into different resonance configurations withJupiter’s orbit. We do not let orbits of these planets migrate,but consider them static. With such a crude approach we canonly get a first hint about a relative role of depletion of theasteroid populations in the first-order resonances (note inreality the planets undergo steady, but likely not smooth,migration and exhibit jumps over different mutual resonantstates; e.g., supplementary materials of Tsiganis et al. 2005).The results are summarised in Fig. 23. (i) The Hilda c (cid:13) , 1–20 steroid families in the first order resonances nu m be r o f a s t e r o i d s i n b i n N lifetime in J2/1 (Myr)no MMR1:24:93:72:5G.I. vs J2/1 0 20 40 60 80 100 0.1 1 10 100 nu m be r o f a s t e r o i d s i n b i n N lifetime in J3/2 (Myr)no MMR1:24:93:72:5 0 1 2 3 4 5 6 7 8 0.1 1 10 100 nu m be r o f a s t e r o i d s i n b i n N lifetime in J4/3 (Myr)no MMR1:24:93:72:5 Figure 23. Histograms of dynamical lifetimes for asteroids in theJ2/1 [top], J3/2 [middle] and J4/3 [bottom] resonances, in caseJupiter and Saturn are at their current orbits, or in a mutual 1:2,4:9, 3:7, 2:5 resonance, or in a Great Inequality resonance (in caseof the J2/1 only). The histograms were computed for 106 long-lived asteroids in the J2/1, first 100 Hildas in the J3/2 and 8 inthe J4/3 (including short-lived). group in the J3/2 resonance is very unstable (on the timescale ∼ on the contrary J2/1 asteroids may survive several10 Myr in this configuration of Jupiter and Saturn, so thispopulation is not much affected by this phase by planetaryevolution (note, moreover, that Jupiter and Saturn wouldlikely cross the zone of other mutual 1:2 resonance in ∼ Jupiter Trojans, which are already known to be strongly unsta-ble (Morbidelli et al. 2005), would have the dynamical lifetime ofthe order 0 . a , q , Q / A U time t / Myr Saturn1:2 MMR 7 7.5 8 8.5 9 9.5 10 0.8 1 1.2 1.4 1.6 1.8 2 2.2 2.4 Figure 24. An N -body simulation of planetary migration drivenby dynamically cold planetesimal disk beyond Neptune, with 10 particles and total mass 50 M ⊕ , and including also 10 mass-lessparticles in the J3/2 resonance with Jupiter. Top: The semima-jor axis a S of Saturn vs time and the position of the 1:2 meanmotion resonance with Jupiter (estimated from the Kepler law(1 / − / a J ). Bottom: The same for asteroids initially locatedinside the J3/2 resonance with Jupiter. The J3/2 asteroids arestrongly destabilised at the very time of the 1:2 Jupiter–Saturnresonance crossing ( t = 1 . 25 Myr) and none of the 1000 test parti-cles survived in the J3/2 region after a mere 0.5 Myr. This meansmore than 99 . − . Initial conditions of planets were: a J = 5 . a S = 8 . 05 AU, a U = 12 . a N = 17 . − . We took the current orbitsof Hildas as the initial conditions for our test particles. Note thedestabilisation of the Hilda region is neither sensitive to preciseinitial conditions nor to the mass of the planetesimal disk; theonly relevant condition is that Jupiter and Saturn cross their mu-tual 1:2 resonance. the migration proceed very slowly, it may cause a significantdepletion of the primordial J2/1 population. In the exact 2:5resonance, the J2/1 population would not be affected at all.According to our preliminary tests with a more com-plete N -body model for planetary migration which includesa disk of 10 planetesimals beyond Neptune, the stronginstability of the J3/2 asteroids indeed occurs during theJupiter–Saturn 1:2 resonance crossing (see Fig. 24). A vastreservoir of planetesimals residing beyond the giant plan-ets, and to some extent also nearby regions of the outerasteroid belt, are probably capable to re-populate the J3/2resonance zone at the same time and form the currently ob-served populations. This is similar to the Trojan clouds ofJupiter (Morbidelli et al. 2005). c (cid:13) , 1–20 M. Broˇz and D. Vokrouhlick´y 50 100 150 200 250 300 350 400 1 10 100 1000 nu m be r o f a s t e r o i d s i n b i n N median lifetime in J2/1 t J2/1 (Myr)sizes multiplied by: 20.20.020.002 Figure 25. Logarithmic histograms of dynamical lifetimes forthe originally long-lived asteroids in the J2/1 resonance, in casethe Yarkovsky effect perturbs the orbits. The sizes of the objects(4–12 km) were multiplied by 2, 0.2, 0.02 and 0.002 for compari-son. Note a stronger instability starts to occur for the factor 0.02(i.e., for the sizes 0.08–0.24 km). The YORP reorientations arenot included in this model. In Section 3.4 we already discussed the influence of theYarkovsky effect on the families located inside the J3/2 res-onance. In course of time it modified eccentricities of theirmembers, but did not cause a large-scale instability; the fam-ilies remained inside the resonance all the time. Here weseek the size threshold for which the Yarkovsky would caseoverall instability by quickly removing the bodies from theresonance.We perform the following numerical test: we multiplysizes of the long-lived J2/1 objects by fudge factors 0.2,0.02 and 0.002, for which the Yarkovsky effect is stronger,and compare respective dynamical lifetimes with the originallong-lived objects. Results are summarized in Fig. 25. Wecan conclude that a significant destabilisation of the J2/1resonant population occurs for sizes ∼ . R is radius in kilo-metres and a orbital semimajor axis in AU): τ YORP ≃ 25 Myr ( R/ (2 . /a ) . (12)Since the Yarkovsky effect depends on the obliquity value,the systematic drift would be changed to a random walk forbodies whose spin axis undergo frequent re-orientations bythe YORP effect. We can thus expect that the YORP effectmight significantly prolong dynamical lifetimes of resonantobjects with sizes ∼ . τ YORP < . 25 Myr for them.We can also check, which orientation of the spin axis 10 20 30 40 50 60 70 80 90 100 1 10 100 1000 nu m be r o f a s t e r o i d s i n b i n N median lifetime in J2/1 t J2/1 (Myr)sizes multiplied by 0.02obliquities: 0 ° ° ° ° ° Figure 26. The same as Fig. 25, but now the Yarkovsky effect iscalculated only for small bodies (0.08–0.24 km) and varies due todifferent values of the obliquity γ = 0, 45 ◦ , 90 ◦ , 135 ◦ and 180 ◦ .Note the retrograde bodies (i.e., with a negative drift da/dt < makes the escape from the J2/1 resonance more likelyto happen. We consider 0.08–0.24 km bodies, clone them5 times and assign them different values of the obliquity γ = 0 ◦ , 45 ◦ , 90 ◦ , 135 ◦ and 180 ◦ . Figure 26 shows clearlythat the retrograde rotation increases the probability of theescape. This is consistent with the structure of the J2/1resonance, for which low- a separatrix does not continue to e = 0. We conclude the remaining very small (yet unobserv-able) asteroids inside the J2/1 may exhibit a preferentialprograde rotation.We perform a similar simulation for the J3/2 population(first 100 bodies with sizes 10–60 km), but now with YORPreorientations (Eq. 12) included. The results (Fig. 27) showthe J3/2 population is much less affected than the J2/1 bythe Yarkovsky/YORP perturbation.We conclude the Yarkovsky/YORP effect may destabi-lize the J2/1 and J3/2 bodies only partially on the 100 Myrtime scale and only for sizes smaller than ∼ . The main results of this paper can be summarised as fol-lows: (i) we provided an update of the observed J2/1, J3/2and J4/3 resonant populations; (ii) we discovered two newobjects in the J4/3 resonance; (iii) we described two as-teroid families located inside the J3/2 resonance (Schubartand Hilda) and provided an evidence that they are of a col-lisional origin; (iv) we reported a new mechanism how theYarkovsky effect systematically changes eccentricities of res-onant asteroids; we used this phenomenon to estimate theages of the Schubart and Hilda families ((1 . ± . 7) Gyr and & c (cid:13) , 1–20 steroid families in the first order resonances 50 100 150 200 250 300 350 400 450 500 1 10 100 1000 nu m be r o f a s t e r o i d s i n b i n N lifetime in J3/2 t J3/2 (Myr)sizes multiplied by: 20.20.020.002 Figure 27. The same as Fig. 25, but for the J3/2 resonance.Instability occurs for the sizes multiplied by 0.002 (i.e., 0.02–0.12 km). The YORP effect and corresponding reorientations ofthe spin axes are included in this case. a is transformedto a drift in e , due to a strong gravitational coupling withJupiter. We emphasize, this is is not a chaotic diffusion in e ,but a systematic drift in e .Another piece of information about the families in J3/2resonance is hidden in the eccentricity e p vs absolute magni-tude H plots (see Fig. 28 for the Hilda family). The triangu-lar shape (larger eccentricity dispersion of the family mem-bers for larger H ) is a well-known combination of two ef-fects: (i) larger ejection speed, and (ii) faster dispersal by theYarkovky forces for smaller fragments. Interestingly, thereis also a noticeable depletion of small bodies in the centreof the family and their concentration at the outskirts — aphenomenon known from ( a, H ) plots of main-belt families,which was interpreted as an interplay between the Yarkovskyand YORP effects (Vokrouhlick´y et al. 2006). Indeed theestimated ∼ D ≃ 10 km asteroids in the Hilda region(e.g., Vokrouhlick´y & ˇCapek 2002; ˇCapek & Vokrouhlick´y2004). Vokrouhlick´y et al. (2006a) pointed out that this cir-cumstance makes the uneven distribution of family membersmost pronounced.We postpone the following topics for the future work:(i) a more precise age determination for the resonant aster-oid families, based on the Yarkovsky/YORP evolution in the( e, H ) space; (ii) a more detailed modelling of analytic or N - The YORP effect tilts the spin axes of asteroids preferentiallytowards obliquity γ = 0 ◦ or 180 ◦ and this enhances the diurnalYarkovsky drift due to its cos γ dependence. ab s o l u t e m agn i t ude H pseudo-proper eccentricity e p (153) Hilda Figure 28. Pseudo-proper eccentricity e p vs absolute magnitude H plot for the members of the Hilda family. The triangular shapeand the depletion of asteroids in the centre of the family (around e p ≃ . 22) is discussed in the text. body migration of planets and its influence on the stabilityof resonant populations. ACKNOWLEDGEMENTS We thank David Nesvorn´y, Alessandro Morbidelli andWilliam F. Bottke for valuable discussions on the subjectand also Rodney Gomes for a constructive review. The workhas been supported by the Grant Agency of the Czech Re-public (grants 205/08/0064 and 205/08/P196) and the Re-search Program MSM0021620860 of the Czech Ministry ofEducation. We also acknowledge the usage of the Metacen-trum computing facilities ( http://meta.cesnet.cz/ ) andcomputers of the Observatory and Planetarium in HradecKr´alov´e ( ). REFERENCES Beaug´e C., Roig F., 2001, Icarus, 153, 391Bottke W. F., Vokrouhlick´y D., Broˇz M., Nesvorn´y D., MorbidelliA., 2001, Sci, 294, 1693Bottke W. F., Vokrouhlick´y D., Rubincam D. P., Broˇz M., 2002,in Bottke W.F., Cellino A., Paolicchi P., Binzel R.P., eds.,Asteroids III. University of Arizona Press, Tucson, p. 395Bottke W. F., Durda D. D., Nesvorn´y D., Jedicke R., MorbidelliA., Vokrouhlick´y D., Levison H. F., 2005, Icarus, 179, 63Bottke W. F., Vokrouhlick´y D., Rubincam D. P., Nesvorn´y D.,2006, Annu. Rev. Earth Planet. Sci., 34, 157Brown M.E., Barkume K.M., Ragozzine D., Schaller E.L., 2007,Nature, 446, 294Broˇz M., Vokrouhlick´y D., Roig F., Nesvorn´y D., Bottke W.F.,Morbidelli A., 2005, MNRAS, 359, 1437Broˇz M., 2006, PhD thesis, Charles University (see athttp://sirrah.troja.mff.cuni.cz/ ∼ mira/mp/)Bus S.J., Vilas F., Barucci M.A., 2002, in Bottke W.F., CellinoA., Paolicchi P., Binzel R.P., eds., Asteroids III. University ofArizona Press, Tucson, p. 169ˇCapek D., Vokrouhlick´y D., 2004, Icarus, 172, 526Chambers J.E., 1999, MNRAS, 304, 793Dahlgren M., 1998, A&A, 336, 1056Dahlgren M., Lagerkvist C.-I., 1995, A&A, 302, 907c (cid:13) , 1–20 M. Broˇz and D. Vokrouhlick´y Dahlgren M., Lagerkvist C.-I., Fitzsimmons A., Williams I. P.,Gordon M., 1997, A&A, 323, 606Dahlgren M., Lahulla J.F., Lagerkvist C.-I., Lagerros J., Mot-tola S., Erikson A., Gonano-Beurer M., di Martino M., 1998,Icarus, 133, 247Dahlgren M., Lahulla J.F., Lagerkvist C.-I., 1999, Icarus, 138, 259Davis 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. University of Arizona Press,Tucson, p. 545Dell’Oro A., Marzari F., Paolicchi P., Vanzani A., 2001, A&A,366, 1053Dohnanyi J.W., 1969, J. Geophys. Res., 74, 2531Durda D., Bottke W.F., Nesvorny D., Encke B., Merline W.J.,Asphaug E., Richardson D.C., 2007, Icarus, 186, 498Farinella P., Gonczi R., Froeschl´e Ch., Froeschl´e C., 1993, Icarus,101, 174Farinella P., Froeschl´e C., Gonczi R., 1994, in Milani A., Di Mar-tino M., Cellino A., eds., Asteroids, comets, meteors 1993.Kluwer Academic Publishers, Dordrecht, p. 205Ferraz-Mello S., 1988, AJ, 96, 400Ferraz-Mello S., Michtchenko T.A., 1996, A&A, 310, 1021Ferraz-Mello S., Michtchenko T.A., Nesvorn´y D., Roig F., Simula,A., 1998a, Planet. Sp. Sci., 46, 1425Ferraz-Mello S., Michtchenko T.A., Roig F., 1998b, AJ, 116, 1491Gil-Hutton R., Brunini A., 2008, Icarus, 193, 567Gomes R., Levison H. L., Tsiganis K., Morbidelli A., 2005, Na-ture, 435, 466Hahn J. M., Malhotra R., 1999, AJ, 117, 3041Henrard J., 1982, Celest. Mech., 27, 3Ivezi´c Z., Juri´c M., Lupton R.H., Tabachnik S., Quinn T., 2002,Survey and other telescope technologies and discoveries, (J.A.Tyson, S. Wolff, Eds.) Proceedings of SPIE Vol. 4836Jedicke R., Magnier E.A., Kaiser N., Chambers K.C., 2007, in Mi-lani A., Valsecchi G.B., Vokrouhlick´y D., eds., Near Earth As-teroids, our Celestial Neighbors. Cambridge University Press,Cambridge, p. 341Kenyon S. J., Bromley B. C., O’Brien D. P., Davis D. R., 2008, inBarucci M.A., Boehnhardt H., Cruikshank D.P., MorbidelliA., eds., The Solar System beyond Neptune. University ofArizona Press, Tucson, p. 293Krueger A., 1889, AN, 120, 221K¨uhnert F., 1876, AN, 87, 173Landau L. D., Lifschiz E. M., 1976, Mechanics, Pergamon Press,OxfordLaskar J., Robutel P., 2001, Celest. Mech. Dyn. Astron., 80, 39Lemaitre A., Henrard J., 1990, Icarus, 83, 391Levison H., Duncan M., 1994, Icarus, 108, 18Malhotra R., 1995, AJ, 110, 420Michtchenko T. A., Ferraz-Mello S., 1997, Planet. Sp. Sci., 45,1587Milani A., 1993, Celest. Mech. Dyn. Astron., 57, 59Milani A., Farinella P., 1995, Icarus, 115, 209Morbidelli A., Moons M., 1993, Icarus, 102, 316Morbidelli A., Levison H. L., Tsiganis K., Gomes R., 2005, Na-ture, 435, 462Morbidelli A., Tsiganis K., Crida A., Levison H.F., Gomes R.,2007, AJ, 134, 1790Moons M., Morbidelli A., Migliorini F., 1998, Icarus, 135, 458Murray C. D., 1986, Icarus, 65, 70Murray C. D., Dermott S. F., 1999, Solar System Dynamics, Cam-bridge University Press, CambridgeNesvorn´y D., Ferraz-Mello S., 1997, Icarus, 130, 247Nesvorn´y D., Beaug´e C., Dones L., 2004, AJ, 127, 1768Nesvorn´y D., Vokrouhlick´y D., Bottke W.F., 2006, Science, 312,1490Nesvorn´y D., Bottke W.F., Dones L., Levison H.F., 2002, Nature,417, 720 Nesvorn´y D., Alvarellos J.L.A., Dones L., Levison H.F., 2003, AJ,126, 398Nesvorn´y D., Jedicke R., Whiteley R.J., Ivezi´c Z., 2005, Icarus,173, 132O’Brien D.P., Greenberg R., 2003, Icarus, 164, 334Palisa J., 1888, AN, 120, 111Press W. R., Teukolsky S. A., Vetterling W., Flannery B. P.,2007, Numerical Recipes: The Art of Scientific Computing,Cambridge University Press, CambridgeQuinn T.R., Tremaine, S., Duncan, M. 1991, AJ, 101, 2287Roig F., Gil-Hutton R., 2006, Icarus, 183, 411Roig F., Nesvorn´y D., Ferraz-Mello S., 2002, MNRAS, 335, 417Roig F., Ribeiro A. O., Gil-Hutton R., 2008, A&A, in pressSchubart J., 1982a, A&A, 114, 200Schubart J., 1982b, Cel. Mech. Dyn. Astron., 28, 189Schubart J., 1991, A&A, 241, 297Schubart J., 2007, Icarus, 188, 189Sessin W., Bressane B.R., 1988, Icarus, 75, 97Tanga P., Cellino A., Michel P., Zappal`a V., Paolicchi P., Dell’OroA., 1998, Icarus, 141, 65Tedesco E.F., Noah P.V., Noah M., Price S.D., 2002, AJ, 123,1056Tsiganis K., Kneˇzevi´c Z., Varvoglis H., 2007, Icarus, 186, 484Tsiganis K., Gomes R., Morbidelli A., Levison H.L., 2005, Nature,435, 459Tsuchida M., 1990, Rev. Mexicana Astron. Astrophys., 21, 585Vokrouhlick´y D., Broˇz, 2002, in Celletti A., Ferraz-Mello S., Hen-rard J., eds., Modern Celestial Mechanics: from Theory toApplications. Kluwer Academic Publishers, Dordrecht, p. 467Vokrouhlick´y D., Broˇz M., Bottke W.F., Nesvorn´y D., MorbidelliA., 2006a, Icarus, 182, 118Vokrouhlick´y D., Broˇz M., Morbidelli A., Bottke W.F., Nesvorn´yD., Lazzaro D., Rivkin A. S., 2006b, Icarus, 182, 92Weissman P.R., Bottke W.F., Levison H.L., 2002, in Bottke W.F.,Cellino A., Paolicchi P., Binzel R.P., eds., Asteroids III. Uni-versity of Arizona Press, Tucson, p. 669Zappal`a V., Cellino A., Farinella P., Kneˇzevi´c Z., 1990, AJ, 100,2030Zappal`a V., Cellino A., Farinella P., Milani A., 1994, AJ, 107, 772Zappal`a V., Cellino A., Dell’Oro A., Paolicchi P., 2002, in BottkeW.F., Cellino A., Paolicchi P., Binzel R.P., eds., Asteroids III.University of Arizona Press, Tucson, p. 619 APPENDIX A: RESONANT YARKOVSKYEFFECT The effects of weak dissipative forces, such as the tidal force,gas-drag force and the Poynting-Robertson force, on bothnon-resonant and resonant orbits were extensively studiedin the past (e.g., Murray & Dermott 1999 and referencestherein). Interaction of the Yarkovsky drifting orbits withhigh-order, weak resonances was also numerically studied tosome extent (e.g., Vokrouhlick´y & Broˇz 2002) but no sys-tematic effort was paid to study Yarkovsky evolving orbitsin strong low-order resonances. Here we do not intend todevelop a detailed theory, rather give a numerical examplethat can both help to explain results presented in the maintext and motivate a more thorough analytical theory. The Yarkovsky effect outside the resonance. We con-structed the following simple numerical experiment: we tookthe current orbit of (1911) Schubart as a starting conditionand integrated the motion of two 0.1 km size objects withextreme obliquity values 0 ◦ and 180 ◦ . Their thermal param-eters were the same as in Section 3.4. Because the diurnalvariant of the Yarkovsky effect dominates the evolution, the c (cid:13) , 1–20 steroid families in the first order resonances o sc u l a t i ng e osculating a / AU N + = . N = . N – = . observed ∆ e expected ∆ a ( Σ ) / s i n σ (2 Σ ) cos σ a /d t > 0100 Myrd a /d t < 0 N time t / Myr N + N – Figure A1. Orbital evolution of two D = 0 . ◦ and 180 ◦ were assigned to the two bodies, such that outside the resonance they would migrate by theYarkovsky forces in opposite direction. The expected change ∆ a of the semimajor axis in 100 My is depicted by the arrow on top of theleft panel. The left and middle panels show 1 kyr orbital segments at the beginning and at the end of the simulation: (i) in the semimajoraxis a vs eccentricity e projection (left), and (ii) in the projection of Cartesian resonant variables √ 2Σ (cos σ, sin σ ) (see Eqs. 1 and 2;short-period variations have been removed for better visibility) (middle). The orbits slowly evolve from the initial N ≃ . 45 level-curveof the integral given by Eq. (5) to their final values of ≃ . 44 ( N + with da/dt > 0) and ≃ . 46 ( N − with da/dt < e as a changes, orbital evolution is characterised by a significant changeof the eccentricity ∆ e (also e p ) but only a small change in a (also a p ). extreme obliquities would mean the two test bodies wouldnormally (outside any resonances) drift in semimajor axisin two opposite directions (e.g., Bottke et al. 2002, 2006).The two orbits would secularly acquire ∆ a ≃ +0 . 25 AU or − . 25 AU in 100 My, about the extent shown by the arrowon top of the left panel of Fig. A1. Since the strength of theYarkovsky forces is inversely proportional to the size, we canreadily scale the results for larger bodies. The resonance without the Yarkovsky effect. If we in-clude gravitational perturbations by Jupiter only, within arestricted circular three body problem ( e J = 0), and removeshort-period oscillations by a digital filter, the parameter N from Eq. (5) would stay constant. The orbit would becharacterized by a stable libration in (Σ , σ ) variables withabout 30 ◦ amplitude in σ (see the curve labelled 0 Myr inthe middle panel of Fig. A1).While evolving, some parameters known as the adia-batic invariants are approximately conserved (see, e.g., Lan-dau & Lifschitz 1976; Henrard 1982; Murray & Dermott1999). One of the adiabatic invariants is N itself. Another,slightly more involved quantity, is the area J enclosed bythe resonant path in the √ 2Σ (cos σ, sin σ ) space: J = I √ dσ . (A1)We would thus expect these parameters be constant, exceptfor strong-enough perturbation or long-enough time scales(recall the adiabatic invariants are constant to the secondorder of the perturbing parameter only). Resonant Yarkovsky effect. Introducing the Yarkovskyforces makes the system to evolve slowly. The lock in the res-onance prevents the orbits to steadily drift away in the semi-major axis and the perturbation by the Yarkovsky forcesacts adiabatically. This is because (i) the time scale of theresonance oscillations is much shorter than the characteristictime scale of the orbital evolution driven by the Yarkovsky forces, and (ii) the strength of the resonant terms in theequations of motion are superior to the strength of theYarkovsky accelerations.We let the two J3/2 orbits evolve over 100 Myr(Fig. A1). At the end of our simulation the orbits movedfrom N ≃ . 45 to N + ≃ . 46 (for the outward migratingorbit) or to N − ≃ . 44 (for the inward migrating orbit), re-spectively. During this evolution, both orbits remained per-manently locked in the J2/1 resonance, librating about theperiodic orbit (dashed line in the left panel of Fig. A1). Be-cause the position of this centre has a steep progression inthe eccentricity and only small progression in the semima-jor axis, the evolution across different N -planes makes theorbital eccentricity evolve significantly more than the semi-major axis. This is also seen in the middle panel of Fig. A1,where the librating orbits significantly split farther/closerwith respect to origin of coordinates (note the polar dis-tance from the origin is basically a measure of the eccen-tricity). The shape of the librating orbit is modified suchthat the area J stays approximately constant. We have ver-ified that the relative change in both adiabatic invariants,acquired during the 100 Myr of evolution, is about the same: δN/N ∼ δJ/J ∼ × − . It is a direct expression of thestrength of the perturbation by the Yarkovsky forces.We can conclude the Yarkovsky effect results in a sig-nificantly different type of secular evolution for orbits ini-tially inside strong first-order mean motion resonances withJupiter. Instead of secularly pushing the orbital semimajoraxis inward or outward from the Sun, it drives the orbitaleccentricity to smaller or larger values, while leaving thesemimajor axis to follow the resonance centre.If we were to leave the orbital evolution continue inour simple model, the inward-migrating orbit would leavethe resonance toward the zone of low-eccentricity apocen-tric librators. Such bodies are observed just below the J2/1resonance. On the other hand, the outward-migrating orbit c (cid:13) , 1–20 M. Broˇz and D. Vokrouhlick´y would finally increase the eccentricity to the value when theorbit starts to cross the Jupiters orbit. Obviously, in a morecomplete model, with all planets included, the orbits wouldfirst encounter the unstable region surrounding the stableresonant zone. Such marginally stable populations exist inboth the J3/2 and J2/1 resonances.This paper has been typeset from a TEX/ L A TEX file preparedby the author. c (cid:13)000