Asteroid families - long-term dynamical evolution due to chaotic diffusion and Yarkovsky effect

Here we shall collect results from our long-term simulations of asteroid families with Yarkovsky effect and related material. Any comments and suggestions are welcome (you may also send us pieces of text you wished to include here); our idea is that we shall continuously update the information, and this page may later facilitate writing a paper on the material. At this moment, we gathered some information about Dora, Eos, Themis and Eunomia families that we have integrated here in Prague, and we also include a link to Nice for Flora, Maria and Eunomia (David N directory; all with the original version of the on-line filter).
  1. Preliminaries, initial data, etc.
    1. Integrator and implementation of the Yarkovsky effect
    2. Thermal parameters and other data related to the Yarkovsky effect
    3. Initial orbital data for TPs
  2. Numerical simulations - Dora, Eos, Themis, Eunomia (and Maria)
    1. Evolution of the mean semimajor axis a(t)
    2. Evolution of the proper elements: a-e projection (includes approximate position of the mean-motion and secular resonances)
    3. Evolution of the proper elements: a-sin(i) projection
    4. Animations with dt = 25 Myr
    5. Semimajor axis distribution of the family members N(a,t)
    6. Center and sigma of families, other histograms
    7. Interaction with high-order and three-body resonances
    8. Examples of interesting TPs (interaction with secular resonances)
    9. Real asteroids of Eos, Eunomia and Maria
  3. Problems and remote simulations
    1. Local directory (ls -l Yarko_fam/)
    2. Problems with dora-1f and Koronis runs
    3. Results on Flora, Eunomia, Maria and Koronis (D. Nesvorny and A. Morbidelli, Nice; B. Bottke, SWRI)
  4. Monte Carlo code including collisional evolution
    1. Motivation and a brief description of the model
    2. Comments on YORP
  5. References

Integrator and implementation of the Yarkovsky effect

The choice of integrator has been constrained by two restrictive demands following from the problem we want to solve: rapidity and confidence as far as the statistical properties averaged over an esemble of bodies are concerned. Symplectic integrators certainly most ideally meet these characteristics, so that we use the swift integrator by Levison and Duncan (???). We have supplemented the original code by the Yarkovsky effect (both diurnal and seasonal variants) subroutines; the reference sources for the corresponding analytic linearized solutions are Vokrouhlicky (1998) and Vokrouhlicky and Farinella (1999). As a technical remark we note that there are two possibilities how to incorporate the weakly dissipative effects in the symplectic integrators: either (i) include them as a perturbation of the free Keplerian propagation (e.g. Malhotra 1994), or (ii) couple them with gravitational perturbations of velocities (e.g. Codeiro et al. 1998, Mikkola 1998). We use the second option. Obviously, the Yarkovsky effect violates the excact symplecticity of the integrator, but we have verified (e.g. Broz 1999) that the degree of dissipation is small enough so that our integrator satisfies all test-bed checks (like reproducing the analytical drifts of the semimajor axis on circular orbits).

Thermal parameters and other data related to the Yarkovsky effect

The linearized modeling of the Yarkovsky effect is likely not to be the most significant source of error. The non-linear aspects of the underlying thermal model are generally small for low-eccentricity (~< 0.3) orbits (e.g. Vokrouhlicky and Farinella 1999, Spitale and Greenberg 2001). More important error source may be due to mismatched thermal parameters of the small asteroids, their irregular shapes and a possibility of an excited (non-principal axis) mode of rotation. Discussion of the last two effect is beyond the scope of the current text, but we shall briefly comment on our choice of the thermal parameters of the integrated asteroids.

For the S-type families (Flora, Maria, Eos, Eunomia and Koronis) we assumed a low-conductive regolith of the thermal conductivity K=0.001 W/m/K (the thermal capacity is 680 J/kg/K and the surface density is 1.5 g/cm^3). These values match reasonably well those of the lunar (and Phobos) surface, but - what more important - are consistent with observations of Eros by Near (???) and small asteroids by ISO (Muller and Lagerros 1999). Bond albedo 0.1 and thermal emissivity 0.95, fits also their mean for S-type asteroids. The assumed value of the bulk density - 2.5 g/cm^3 - corresponds to S-types.

For the C-type families (Themis and Dora) we also assumed low surface conductivity of K=0.01 W/m/K; here the somewhat higher value than given above for S-types should correspond to ice components in the surface regolith (ice is a good conductor and increases the effective value of K). Here the direct observational evidence is virtually absent, except of models of cometary activity. We may thus admit somewhat larger uncertainty in modeling the Yarkovsky effect for the C-type families. Bond albedo 0.02 and thermal emissivity 0.95, fits also their mean for S-type asteroids. The assumed lower value of the bulk density - 1.3 g/cm^3 - should follow from icy components and high porosity.

It may be worth to point out that despite of persisting uncertainty in the Yarkovsky effect modeling for individual orbits, the corresponding error of the statistical quantities (like the growth of the family dispersion in a) may be much smaller. Though we do not have an exact prove for this statement, this is a general wisdom from other studies.

Initial orbital data for TPs

Morby wrote a brief report from the last summer "get-together" where we specified the procedure how to generate the initial orbital data for TPs (planets are from JPL ephemerides). Here is the corresponding TEX source file or a compiled ps file.

Evolution of the mean semimajor axis a(t)

The mean semimajor axis is output from the simulation every 1.5 kyr; they are thus sampled in time densely enough to see the fine structures of the miniresonances. Each of the runs was split into two parts of 110 and 100 asteroids. In the first one (called *-1) we included larger objects so that this gets less dispersed; the secodn group (called *-2) includes smaller asteroids and the dispersion is larger. Integrations of Dora and Eos were completed to 300 Myr, Themis up to 500 Myr; remember also, that with exception of Flora we included only the four outer planets in our simulations. Here we go with specific comments on integrated families (by clicking on the figure you get a magnification).

Dora: integration time span 300 Myr, *-1 run contains asteroids of sizes in the range 4-19 km, while *-2 run contains asteroids of sizes from 4 km downto 2km (smallest objects in each of the simulations have D=2 km). Here is what Bill B wrote about the data of this family (with a plot in proper-element space). Dora is specific by the closeby J5/2 resonance. The four km-sized object in the run dora-2f, that crossed the J5/2 resonance, would likely be ejected by close encounters with inner planets. About 25% of small asteroids (and fast drifting) asteroids in dora-2f escaped through the J5/2 resonance. At the other side of the family, asteroids that reach the 3J-1S-1 resonance (at ~2.75 AU) are significantly scattered in eccentricity, so they would not be identified as family members.

Eos: integration time span 300 Myr, *-1 run contains asteroids of sizes in the range 4.3-65 km, while *-2 run contains asteroids of sizes from 4.3 km downto 2km. Here is what Bill B wrote about the data of this family (with two plots in proper(a) vs proper(e) space: R>15 km and 15 km>R>5 km). The strong mean motion resonance J9/4 is a significant feature in this family, and a number of asteroids may leake out by increasing its eccentricity. Another, but more subtle resonance, that may affect the long term dynamical evolution in this family is the secular resonance g+s-g6-s6. This cannot be noticed in the a(t) plot, but will clearly appear in the proper element space. Milani and Knezevic (1992) briefly analysed Eos members as regard to their possible location in this resonance, but limited computer power prevented them to perform a complete study (they did only 10Myr runs for several asteroids but definitely not all Eos members). Looking at Bill B' data for 5-15 km sized asteroids (a-e proper element plane), one may be tempted to "believe" that a non-negligible fraction of the Eos family members are located/close to this secular resonance, but it may be just an "optical" illusion. To be more precise at this point, we ran orbital evolution (without Yark-effect) of 520 Eos members (all currently known members according to Cellino's list) for something like 30 Myr. This should be enough to determine which of the asteroids reside in this peculiar resonance (our preliminary look at the data say that there are some 31 percent of Eos members interact with this secular resonance). We hope that this might be an interesting test to trace footprints of the ongoing mobility in the Eos family, since we clearly see that Yark-moving orbits get ``stick for a long time'' to this resonance. First results can be found below.

Additionally, I guess that David N did check that several larger members of the Eos family are clustered in the 8J-3S-3 three-body resonance. Indeed, we checked that several of the slowly drifting asteroids in eos-1f get temporarily captured in this resonance. Above the J9/4 resonance the motion seems to be quite irregular (other resonances?), which may add to the reason why we see so few Eos objects there.

Themis: integration time span 500 Myr, *-1 run contains asteroids of sizes in the range 4.4-110 km, while *-2 run contains asteroids of sizes from 4.4 km downto 2km. Here is what Bill B wrote about the data of this family (with two plots in proper(a) vs proper(e) space: R>20 km, 20 km >R>10 km and 5 km>R>10 km). The family is very broad and covers a portion of the outer belt with many miniresonances. Obviously, the principal resonance at the outer edge is J2/1. We centered the initial data around ~3.13 AU, so that only few (4) of the smallest asteroids reached the J2/1 resonance within the 500 Myr limit. All of these were ejected. We obviously keep track of the puzzling group of Zhongguos in our plans and we intend to perform an additional simulation by placing a significant number of test particles at the bottom part of the J2/1 resonances and let them drifting up in semimajor axis. Our simulation indicates that closeby (at larger eccentricity even overlapped) resonances 3J-2S-1 and J11/5 at the bottom part of the family are very efficient in dispersing bodies in the proper(a)-proper(e) plane, so that those that crossed this area may not be identified as the family members. Quite powerful are also the first-order resonance 5J-2S-2 (at ~3.173 AU) and J13/6 (at ~3.105 AU). There are also weaker resonances (7J-2S-3, J15/7, J17/8, etc.) with still large enough effect to be noticed in our simulation. As it will be seen below, these miniresonances are quite efficient in dispersing the family in eccentricity, which is "good" since Themis has also considerable extension in e.

See also alternative images (with different semimajor axis range).

Eunomia: here we give only preliminary figures of approximately 75 Myr integration (we plan to finish this integration to ~300 Myr as for the other families). The sizes of asteroids range from 45 km to 4 km in *-1 run and from 3 km to 2 km in *-2 run.

See Bill B. text and his (a, e) plot concerning this family. The evolution of Eunomia spanning 175 Myr was also calculated by D. Nesvorny and A. Morbidelli (with previous version of on-line filter :-(

Remarkable mean motion resonances in the Eunomia region seem to be J11/4 with Jupiter and weaker 3-body resonances 6J-1S-2, 9J-2S-3, 7J-4S-2. We also noted the role of s-s6-g5+g6 secular resonance (see below for examples).

Evolution of the proper elements (a-e projection; includes approximate position of the mean-motion and secular resonances)

In what follows we show evolutionary tracks of the integrated asteroids as projected in the proper(a)-proper(e) plane. The proper elements were computed for by applying a 5 Myr averaging window to the corresponding mean elements with a step of 0.1 Myr. All 210 particles are shown together, but you may also see separate plot for runs *-1 and *-2. Borders of secular resonaces are taken from Milani & Knezevic (1994); as mentioned above the remarkable secular resonances are g+s-g6-s6 (together with g+s-g5-s7, which is always very close to this resonance) in case of Eos family and s-s6-g5+g6 in the Eunomia region. Location of the two-body mean-motion resonances, corresponding to separatrix position of the given resonance within the restricted circular three-body problem, are taken from the data file provided by David N. The position of the three-body resonances are indicated by their centers only (for this moment).

The observed family members are superimposed over the simulation by crosses; data from Bill Bottke coincide well with those indicated by those produced by A Cellino. For this moment we used only the numbered members with numbers smaller than ~5400 (hope to extend to all currently numbered members soon). The grey crosses are fair members (Cellino's quality parameter 2 or 3), while red crosses are tentative members (Cellino's quality parameter 1).

In the right hand side plots you may compare the observed family (the orange cross then indicates the larges member of the family) with the initial conditions of our simulation: (i) osculating initial elements, and (ii) proper elements computed from the first 5 Myr window. The first 5 Myr evolution in the family typically does not change the elements much; the exception id Eos, where you may already notice the influence of the g+s-g6-s6 and g+s-g5-s7 resonances. Notice the typically large difference between the osculating and proper eccentricities, which is due to forced oscillations of the osculating elements due to planetary perturbations. The ``initial'' proper eccentricities are typically well-positioned in the observed family, with only a small misfit for Dora and little larger misfit for Koronis. Perhaps, we were not enough careful when producing the initial conditions, but I believe that this shift is not essential for our conclusions. What really matters is how the family disperses in time.

not yet available

Evolution of the proper elements (a-sin(i) projection; includes approximate position of the mean-motion and secular resonances)

Position of the secular resonance g+s-g6 -s6 depends sensitively on all proper elements (unlike the mean motion resonances). As a result, we plotted it projection for e=0.10 and e=0.06; the latter applies to the late evolution of the particles (see also a vs e plot above). The same is tru for the s-s6-g5 +g6 secular resonance in the Eunomia family. This later resonance seems to influence interestingly proper inclination (more than the proper eccentricity) resulting in large oscillations at the bottom part of the family. As for Eos, I have a ``feeling'' (by looking at the Cellino data file of Eunomia family) that proper inclinations do get anomalously perturbed around a ~ 2.6 AU -- is this really the ``footprint'' of this secular resonance or not? See the upper particle for clearly ``feeling'' librations of the corresponding resonance angle (period of ~10 Myr). Notice also the Themis case: though the weak resonances clearly influence the proper inclination the effect is rather small (see also discission below).

As above, the observed numbered (<5400) family members are superimposed with symbols (grey/red crosses depending on the quality factor of the identification) over the numerical simulation.

The right hand side plot indicate comparision of the family (in the appropriate proper-element projection), the initial osculating elements and the proper elements computed from the first 5 Myr of our simulation. Things are ``reasonably'' good except from Themis (see also below, where we plot the histograms of the proper element distribution). Anybody has an idea? small-inclination problem? after all the osculating elements typically oscillate between 0 - 3.5 degrees.

not yet available

Animations with dt = 25 Myr

Some of the features of our results may be better seen in the animations of the family evolution in the proper element space; here we give projection on the proper(a)-proper(e) plane for Dora, Eos, and Themis (you may also get individual snapshots at all times (EPS images here)). The same is available for proper(a)-proper(i) plane and Dora, Eos, and Themis family (you may also get individual snapshots at all times (EPS images here)).

Notice, for instance, how the Themis family becomes more dispersed in eccentricity for a > 3.15 AU due to interaction of the Yarkovsky drifting orbits with miniresonances. The real data of the family indicate that this is true; moreover, a number of real family members located in the 5J-2S-2 (and widely spread in proper eccentricity) are an indirect hint for capturing of the Yarkovsky-drifting orbits in this resonance (see the last slide of our simulation - at 500 Myr - with 8 particles captured in this resonance).

The animations above show one more feature: sizes of the asteroids (estimated), or TPs (modelled), are indicated by the size of the symbols (circles). Dora looks good and the family may be something like twice old than the timespan of our simulation (~ 0.5-0.8 Gyr say). As expected, Themis (for instance) is a more difficult case. Please do not be ``frustrated from scratch'' by seeing that only small (~ km-size) particles disperse enough in our simulation, while the observed family (gold symbols) indicate that also larger asteroids are well dispersed between ~3.07 - 3.22 AU. Take into account that the timespan of our simulation (0.5 Gyr) may be just a fraction of the family age. If the family were 2.5 Gyr old (i.e. 5 times more than our simulation) then ``multiply'' the size of the smallest TPs in our simualtion by this factor 5 (since Yark-mobility is inversely proportional to the size). We may thus get objects of the ~5 - 7.5 km radius (the smallest observed) well dispersed throughout the whole family in the 2.5 Gyr age. This conclusion is supported by the Monte-Carlo simulations discussed at the bottom of this page.

Semimajor axis distribution of the family members N(a,t)

Here we show how the (initially tight) distribution of semimajor axes of the family members diserses in course of time. For each of the family we somewhat arbitrarily give the distribution (i) at the initial time (black), (ii) at the middle time of the integration time span (gray) and (iii) at the end of the integration time span (red; this latter instant obviously depends on the family, being 300 Myr for Dora and Eos and 500 Myr for Themis).

The distributions are compared to the identified members (yellow curve); here we used data files by Cellino (provided to us by Morby) but later we would like to compare it also with data files of Bill B. Number of observed asteroids in each family (e.g. from 78 upto 550 in Cellino's files) is normalized to 210, to be able to directly compare with our simulation.

A complete list of histograms is N(a), N(e), N(i) is available.

family semimajor axis N(a) eccentricity N(e) inclination N(i)

2DO: comments on distributions of a, e, i of all other families.

Themis. You may notice that in eccentricity our simulation still "lags behind" the true family dispersion, but we should take into account that 500 Myr may represent only a 15-50% fraction of the family age (so that further dispersion in the modelled data may be expected).

Themis inclinations. We are a bit concerned about the Themis distribution of inclinations and its comparison with the data from A. Cellino's file. Though the eccentricity dispersion might eventually fit (after extrapolating evolution, doubling or tripling the family age), inclinations seems to present a problem. The family data have much broader distribution in i with a shallow maximum at ~ 1.3 degrees. This is in contrast with Themis proper inclination of 1.08 degrees (rather off the peak of the distribution); the miniresonances seem to apreciably affect the proper eccentricities, but much less the inclinations (see the above plotted evolution of the family in the proper a - proper i) projection). Anybody has some idea?

Additional informations on this issue: (i) a typical oscilation of osculating inclination is more than 3° (with mean value above 1.5°, which corresponds to the center of "our" histogram) - as could be seen in i(t) plot for one particular particle. (ii) The initial osculating inclination of Themis - parent body was set around 0.8° - see list of initial osculating elements (a,e,I,...) of themis-1 run and this (a, i) plot with comparison of initial osculating and proper elements with the real family. The fact that our proper inclination is centered around 1.8° means probably that we were not very careful in chossing the initial elements of the asteroids (they were ``too low in the inclination oscillation cycle''; ?). (iii) You can also compare proper elements from A. Cellino data-file with those by Milani and Knezevic (downloaded from AstDys) - they agree very well. (By the way - note a clear "trace" of J11/5 resonance at approx. 3.75 AU visible on this image - very few asteroids were identified as family members "behind" its border.)

Center and sigma of families, other histograms

Center of families in semimajor axis (mean value of all TP's), eccentricity or inclinations changes only little (see the following images), unless the family is located in the neighbourhood of a strong mean motion resonance. In this situation, a fair fraction of the asteroids may be ejected, and the family cannot "disperse symmetrically". This is the case of Dora due to J5/2 resonance. Here is how the family-center in a evolves in time. The fluctuations are likely associated with individual particles leaking out of the family through J5/2.

The initial and final standard deviations of the distribution in semimajor axes are as follows: Dora 0.005 AU and 0.025 AU, Eos 0.012 AU and 0.018 AU (both after 300 Myr), Themis 0.012 AU and 0.042 AU (after 500 Myr). As expected, in all cases the family dispersion increases (after a certain time the increase is nearly linear in time). Obviously, this is the Yarkovsky feature. As it will be discussed below, this trend should stop in reality because of collisional evolution in the family.

The corresponding values of dispersion for the eccentricity and inclination are also shown in the figures collected in the following table. Here the evolution is less reguliar, especially in the the case of Dora where only a few of miniresonances are located. The best seen is for Themis, where (on the contrary) a large number of miniresonances are located and are able to temporarily alter proper eccentricity or inclination. On contrary to the semimajor axis, the quasilinear trend here should be affected less by the collisional evolution of the family. As before, the fluctuations correspond to fast ejections of particles from the simulation.

To compute the standard deviation of the family dispersion, we should obviously ``somehow'' eliminate asteroids that leak out fo the family through principal resonances (e.g. J5/2 for Dora), since these are no more family members and would artificialy increase the computed dispersion. Ideally, we might contrast the family with a background population in this part of the asteroid belt and check at each timestep which of our integrated asteroids would be still identified as the family members. For this moment, we chose a somewhat simpler option and we use the following algorithm (applied to each element): (i) compute center of the family, (ii) compute standard deviation of the distribution and eliminate all asteroids that exceed 3-sigma distance from the center, (iii) iterate back to (i) until no more asteroids are eliminated by the procedure. We then report the final value of the family-center and standard deviation in the corresponding element.

The first column in the table shows evolution of the proper semimajor axis as a function of time in our simulations, computed center of the family (red curve), 3-sigma range (blue curves) and the standard deviation (orange curve and the right ordinate). You may see, that e.g. for Themis we have a few asteroids that do drift so fast that they escape from the 3-sigma range from the center and are ``eliminated'' from the family. The work is done similarly for the eccentricity and inclination.

family a(t) plot standard deviation in semimajor axis ... eccentricity ... inclination
(75 Myr)

Interaction of the Yarkovsky-drifting orbits with high-order and three-body resonances

Dependence of Delta e (ie. maximum - minimum eccentricity during the simulation) on the quality of the linear fit of a(t), characterized by its sigmafit (roughly a square root of the chi2 value), has been previously suggested by Morby as a good indicator to see how Yarkovsky replenishes miniresonances by particles. See here his mail.txt for detailed explanation. So for each particle we performed a linear fit of its proper semimajor axis as a function of time, computed the corresponding sigma and plotted against the difference of maximum and minimum values of e. Scales are logarithmic. The clear correlation of the plotted quantities shows, that particles which contribute most to the eccentricity dispersion of the family did undergo some Yarkovsky evolution. Note that if the body would originally sit in the miniresonance and leak slowly through in e will be situated in the left upper corner of this plot, since its stationary value of a might be well represented by a linear fit. A little caveat of the current version of the plot is that we did not select particles that remained in the simulation (and that would correspond to the "observable family today") but from all 210 particles. Thus the cloud of particles in the Dora case that have large value of Delta e are those that leaked through the J5/2 resonance. As expected, the most outstanding example of the Yarko-feeding of the miniresonances is Themis. In the case of Eos we have some particles with a good a(t) linear fit and large Delta e: quite likely these are those that were sliding along the g+s-g6-s6 secular resonance.

Resonance crossing analysis. You may also get a more direct feeling how the ``weak'' resonaces affect the mean semimajor axis by going here; for selected resonances (such as J13/5 (3J-1S-1), J9/4, J23/10 (5J-1S-2), J11/5 (3J-2S-1), J13/6, etc.) we zoom at their close vicinity. We also analyse in a quantitative manner the resonant captures/jumping for those resonance, that were encountered with ``statistically'' enough particles.

Examples of interesting TPs

  1. Secular resonance g+s-g6-s6 in Eos region. TPs from the Eos region might be captured in g+s-g6-s6 (or g+s-g5-s7, which is always very close). 110 particles (ie. 52 %) from total number of 210 interacted with these resonances at least for short time-interval - see the list of these TPs. You can also inspect separate proper (ae) and a(t) figures for each TP. (a(t) graph is "zoomed" to the surrounding of J9/4.)

    Critical angles. We produced plots for three particular test particles from eos-2f run: no. 1 - proper (a, e), a(t), e(t), g+s-g6-s6 (t) for 0-300 Myr; notice the first ~ 140 Myr phase, when the particle was ``sliding'' along this secular resonance (as a result of the Yark-mobility in the semimajor axis) until it hit the zone of J9/4 resonance and later was scattered out of the family; no. 3 - proper (a, e), e(t), sin i(t), g+s-g6-s6 (t), g+s-g5-s7 (t) (you also can view the evolution of node, peri and varpi for this asteroid and also for Jupiter and Saturn); and no. 62 - proper (a, e), (a, sin i), e(t), g+s-g6-s6 (t), g+s-g5-s7 (t); this is again an interesting case of a TP, that has been ejected from the g+s-g6-s6 resonance at ~ 125 Myr, likely as a result of the interaction with the 5J-1S-2 multiple resonance, but later again reach the secular resonance.

    We observe strong perturbations of e and i, evidently driven by the secular resonance - TPs "slide" exactly along expected resonant borders, and also corresponding librations of critical angles with the same periods. The g and s frequences for planets were taken from Milani and Knezevic (1997) software secre8. (The longitudes of nodes do not circulate in our ecliptic reference frame - David Nesvorny comment.)

  2. Eunomia secular resonances. Several asteroids from eunomia-2 run exhibit similar large oscilation of proper inclination with period upto 13 Myr. Proper (a, e) plot suggests s-s6-g5+g6 resonance (according to its analyticaly computed position). It is confirmed by critical angle plot for TP no. 52, which resembles above mentioned evolution of proper inclination.

  3. Suspicious increasing eccentricity of one TP from Themis. TP (id 110+34) from the themis-2, which is stable at a = 3.13 AU, exhibits a constant increase of eccentricity. This capture (in unknown resonance) persisted for more than 400 Myr, the particle started to drift in a again just before the end of simulation.

Real asteroids of Eos, Eunomia and Maria

520 Eos members (eos-a1 run)

We were interested to see which of the real members of the Eos family are located in the g+s-g6-s6 secular resonance. Our ``vague'' idea was to see whether analysis of the Eos-family shape in the proper element space, in particular a sort of a ``stream'' of asteroids in the bottom-left part, could support the necessity of the ongoing mobility in the family. Recall, that without mobility the asteroids would stay in this resonance at approximately fixed positions (since the resonance is anyway rather weak; one would notice only small oscillations of ``proper e'' and ``proper i'' with a typical period of ~ 5-10 Myr). In contrast, when Yark-effect moves the asteroids in the semimajor axis, the long-term orbital evolution follows/sticks to the resonance.

For purpose of this test, we integrated 520 known members of the Eos family (data from the Cellino's file; initial data from astrob.dat) for 25 Myr. Yarko-effect was not taken into account here, so that we could you the original version of swift; inner planets were not integrated. We looked at the results after 18Myr of this simulations have been completed and produced the graph on right. You can see Eos as projected on the proper (a,e) plane; 162 out of 520 asteroids interact with the g+s-g6-s6 resonance (ie. 31 %). We can say that because we individually looked at the evolution of the g+s-g6-s6 resonant angles. In particular, for asteroids indicated by red diamonds the corresponding resonant angle permanently librates, for asteroids denoted by blue squares we see alternation between librations and circulations (since they lie at the border of this resonance or some of them are in J9/4 resonance that temporarily drives them to this secular resonance), and the black crosses are out of the g+s-g6-s6 resonance. We also indicate a rough position of this resonance as given by the semianalytic approach by Milani and Knezevic; note, however, that this resonance is a really 3D structure in the proper element space (unlike the mean motion resonances); the indicated position corresponds to the mean inclination in the family, but we can also see some ``red-asteroids'' ouside, meaning that their inclination does not fit the family-mean (and vice-versa: there a few black crosses where we would expect the resonant asteroids, but the point is that their inclination is ``out''). Clearly, the bulk of the resonant asteroids is where we expect them and there is also the weak-feeling of a stream of asteroids at the left-bottom part. Obviously, we should not exaggerate, but we anyway found this test interesting. We would be happy to hear suggestions what more we could do?

A drawback of the previous plots stems from the fact, that we have projected all the Eos members onto a single planes of proper(a) vs. proper(e) or proper(a) vs. proper(I). This fact ``messed'' the 3D structures in the family. Note, that this is not important for the mean motion resonances, that are basically identified by the proper(a) value (and only weakly open with the increasing value of the proper(e) or proper(I)). However, the secular resonances, like the one we are studying here, are fundamentally 3D structures in the space of proper elements. In the following animation and plots we cut the Eos family onto slices in proper inclination (each of the 0.25 degrees thick) and show the data for each of the slices separately. This helps understanding the 3D aspects of the Eos structure. Few immediate conclusions: in each slice we plot (i) the Eos members and by red symbols we point out those that were found in the resonance (through libration of the corresponding resonance angle), and (ii) the approximate position of the resonance by using the software by Milani and Knezevic. We may observe, that the way how the ``red asteroids'' fit to the expected location of the resonance is perfect. More importantly, this analysis helps to really see how the Eos family is structed in the proper element space. It is nothing like a triaxial ellipsoid. In the low-inclination ``frames'' we may notice a ``jet'' of asteroids with low values of proper(a) and proper(e), a structure that may be formed by the resonance.

In our Yarko-simulation we do see how the TPs slide along this resonance. For sake of a direct comparison, we have plotted the final state of our 300 Myr evolution of the Eos family with the Yarkovsky effect using the same technique (i.e. in the frames of quasi-constant proper inclination). Here you can see the corresponding animation or the plots. Note, that 300 Myr may not be enough timescale to reach the position of the outermost asteroids in the zone of the z1 resonance, but we are able to reproduce the general features of the Eos structure here.

A still more involved analysis is summarized here:

NEW: Spectroscopic observations in the Eos region (here we discuss in a more detail the potentially new Eos members in the z1 region; special emphasis is given to selection of the objects to be observed).

Two more integrations were performed to check the behavior of z1 critical angle of each "candidate" asteroid for spectroscopic observation. (These are NOT mentioned on previous "public" spectroscopy web page!)

First we selected 56 bodies (eos-a2 run) with frequency z1 < 0.5 "/Myr: see their evolution of critical angle and eccentricity, list of non-librating and alternating (+) bodies (they were NOT dicarded from our list of candidates) and 3 slices in inclination (including also positions of largest asteroid (8340) Mumma and asteroids from eos-a3 run - see below).

In the next set of 20 asteroids (eos-a3) with 0.5 "/Myr < z1 < 1.0 "/Myr most of them also librate, surprisingly. Look at: critical angles and eccentricities and list of non-librating bodies (these WERE dicarded from the candidate list).

514 Eunomia members (eunomia-a1 run)

A similar simulation as in the case of Eos family - see s-s6-g5+g6 resonant angles, sections through the (a, e) plane: animation and individual figures.

81 Maria members (maria-a1 run)

We have carried out ~ 35Myr simulation of 81 members of the Maria family (as identified by A. Cellino); as before Yarko-effect not taken into account in this job. Here we show 5Myr running-box averaged mean elements (i.e. not yet properly defined resonant-proper elements) as projected onto the (a, e) and (a, i) planes (with borders of "all" 7 secular resonances in this region), and the evolution of z2 resonant angles for each TP and comparsion with their proper eccentricities; 6 out of the 81 asteroids are located in the z2 resonance, another 7 asteroids show a very slow circulation (i.e. period of circulation longer than ~ 3Myr). Typical examples of the long-term eccentricity evolution for particles in the z2 resonance are shown here: TP no. 31 or TP no. 43 (in both cases the secular resonance triggers eccentricity oscillation with periods of ~ 7.5Myr and ~ 5Myr, respectively). Obviously, this corresponds to the libration period of the corresponding resonant angle.

We are preparing (a, e) plots of the family for different slices of the ``proper'' inclination; check here later for the corresponding animation or the individual figures.

Problems with dora-1f and Koronis runs

In the run dora-1f the program had bad filter parameters between approximately 80 and 120 Myr. The proper eccentricity for this time interval was calculated as an 5 Myr average from osculating elements and thus could be more "noisy".

The integration of Koronis family is performed on a Windows NT computer, TPs sometimes perform unexpected jumps in semimajor axis during frequent restars of the runs.

Results on Flora, Eunomia, Maria and Koronis

Flora, Eunomia and Maria families were calculated by D. Nesvorny and A. Morbidelli from Nice observatory - check their web or its local copy for preliminary results.

See also SWRI web for new results on Flora, Eunomia, Maria and Koronis families:,

Motivation and a brief description of the model

The previously described simulations contain the major dynamical features of the system with many complex details. Interaction of the Yarkovsky drifting orbits with miniresonances is most outstanding of them. However, there are several important caveats to be reminded:
  1. number of integrated particles (210) is limited;
  2. the simulation timespan is typically smaller than the assumed age of the family (this may be a significant factor for Themis, Eunomia and Koronis families especially);
  3. collisional evolution of the family is entirely neglected (and this includes also the nondisruptive impacts that may affect orientation of the spin axes and thus strength and even direction of the Yarkovsky drift);
  4. other phenomena that may on a long-term affect orientation of the spin axes have been neglected (this includes mainly the YORP effect, that might be capable to despin rotation period and tilt the axis of a ~10km object on a time scale of ~100-500 Myr).
For these reasons, extrapolation of the numerically simulated results (e.g. increase of the semimajor axes dispersion given above) to the entire family age is far to be trivial. To get a first glimpse of the undergoing very-long-term processes in the families we have contructed a model, that significantly simplifies the dynamical side but incorporates collisional (and other so far neglected) aspects at least on an elementary level. The model resembles that of Vokrouhlicky and Farinella (2000), used for modeling the meterite delivery from the inner part of the main asteroid belt. Here are the principal features.
  1. Dynamical side of the evolution is entirely restricted to description of the Yarkovsky dispersion of the family in proper semimajor axes. We thus do not integrate orbits, but we just shift the proper semimajor axis as described by the linearized solution of both diurnal and seasonal variants of the Yarkovsky effect (consistent with the level of modeling of the corresponding phenomena in the numerical code above). Evolution of proper eccentricities and inclinations are excluded from the current version of our code.
  2. The porcedure by which we generate initial data overlaps to some extend with the one we used for the numerical simulations above. In particular, the initial distribution of the semimajor axes (following from the assumed velocity field of break-up fragments) is the same as above. The principal free, adjustable parameter here if f_KE, which indicates how much of the energy needed for disruption of the parent object is deposited to the kinetic energy of fragments. Smaller value means initially compact family; in the simulations of Themis we assume f_KE in the range 0.01 to 0.025, leading to initial values of the sigma(a) in the range 0.01 - 0.025 AU (the equivalent mean ejection velocities are 50-150 m/s). Initial spin axes of the asteroids are randomized in space. We assume initialy a power-law for size distribution of fragments/asteroids, so that the number dN of partiles in the log-bin d log R is ~ R^{3(1-q)} d log R (corresponding to mass distribution ~ m^-q d log m). There is a nice discrete model for generating such an ensemble of bodies by Kresak (1979). At large sizes, the families are steeper than the equilibrium (Dohnanyi) distribution with q ~ 1.83. Since these large bodies are likely to be primordial, we assume for Themis q ~ 1.89 initially; this means that Themis (largest body) represents about 20% of the parent body mass. By taking ~250 km size Themis as the largest body in our ensemble, we fix the absolute number of the family members. In our simulation we retain all bodies larger than 2 km in size (smaller objects are neglected). As a result, the initial population contains approximately 200000 objects. As before, the rotation periods are assumed to have Maxwellian distribution peaked at 8 hr and having cut-offs at 4 hr and 12 hr.
  3. We include an ad hoc model of the collisional evolution of the family; this means that our model is not equivalent to a self-consistent collisional model (e.g. Marzari ???). For asteroid of a given size we assume a given flux of impactors fixed throughout the evolution (i.e. not responding on the long-term changes of the size distribution in the family). Specificity of the Yarkovsky-collisional code is that it should consider also nondisruptive collisions with an impactor, that is able to alter rotational state of the target (change the rotation period and, in particular, reorient the spin axis). In principle, any projectile may to some extend influence the rotational parameters of the target; obviously smaller ones by a small amount but what matters is a cummulative affect of many small effects. Like the orbital semimajor axis, the orientation of the spin axis may "random walk" as a result of small effects. However, time from time a larger, but still nondisruptive impact, may significantly change the rotation state of the target. This was the simplified approach of Farinella et al. (1998) and we use it here as well. We thus assume a possibility of a random sequence of nondisruptive impacts that totaly reorient original direction of the spin axis (rotational effect is not modeled) obeying a Poissonian statistics with a characteristic time scale \tau_rot ~ cf1 * 15.0 * \sqrt{R} Myr; R is radius of target in meters and the cf1=1 result has been derived by Farinella and Vokrouhlicky (1999) for S-type objects and fluxes in the inner part of the asteroid belt. To allow for compositional modification, and also modificaion due to other location in the belt, we introduce an adjustable parameter cf1 (usually in the range 0.5-5). Similarly, the collisional disruptions are modelled as a Poissonian process with a different characteristic timescale of \tau_disr ~ cf2 * 16.8 * \sqrt{R} Myr (cf2 is again an adjustable free parameter). It is important to note that the ~ \sqrt{R} proportionality follows from the assumed Dohnanyi (equilibrium) size distribution of projectiles, which may not be exactly satisfied in reality. A slightly different power of R may then appear in the previous formulae for \tau_rot and \tau_disr.
  4. When one of the simulated objects is assumed to be disrupted, it is replaced by a swarm of fragments with the same power-law distribution as the initial population and mass equivalent to the disrupted body. The lower cut-off at 2 km size is again applied.
  5. Propagation of the system is done with steps of 0.2 Myr; at each of them we update proper semimajor axis of each of the particles (given estimated drift rate due to Yarkovsky effect) and test the possibility of the spin axis reorientation or disruption (by comparing a generated random with Poissonian probability that the corresponding event occured during the 0.2 Myr timestep). If the particle reaches some of the major resonances in the belt (e.g. J2/1 in the Themis case) it is assumed to be ejected from the system; the resonance border is approximated by a simple analytic formula.
  6. The current age of the family is a free, adjustable parameter in the simulation; for Themis it ranges from 2 Gyr to 3.5 Gyr.
Our goal is to contrast the initial and final parameters (semimajor and size distribution, etc.) of the simulated family and then contrast both with the corresponding parameters of the currently observed family. The first comparison shows the nature of simulated phenomena (e.g. rapidity with which the semimajor axis ditribution may increase due to the Yarkovsky effect). The latter comparison is of obvious interest.

To recall, the principal "free parameters" of the simulation are: (i) the age of the family, (ii) the f_KE factor that tunes initial compactness of the family (in semimajor axes), (iii) the empirical factors cf1 and cf2 of the axis-reorientation and disruption timescale. We may also test sensitivity of the results on changing the steepness q of the initial size distribution.

Which parameters we store in course and at the end of the simulation?

  1. Distribution of the semimajor-axes in the population for asteroids in four different size bins of radii: 1-5 km, 5-10 km, 10-15 km and 15-20 km;
  2. we output mean semimajor axis and standard deviation ("sigma") of the semimajor-axes distribution in each of the size bins every 100 Myr;
  3. size distribution of the population members;
  4. we keep trace of major disruption events (bodies with size > 40 km).
Here is an example of results that I get. In this particular case I took f_KE = 0.011 (fairly low initial dispersion of the family), age of 3.47 Gyr (large to somewhat compensate low value of f_KE), cf1 = cf2 = 4.4 (large). I have noticed that letting the objects "survive" longer (i.e. longer \tau_disr) helps expansion of the family, since larger objects near the center of the family do not typically disrupt within its age (and thus do not reduce sigma(a) at smaller sizes). Could we find justification for longer lifetime than expected by our simple formula at the outer part of the belt?

In the first figure you see initial (green) and final (black) distribution of the modeled Themis family in the above indicated four size bins. You may clearly see how the family expanded in course of time. Notice a decrease at ~ 3.07 AU (first two bins); this is my elementary attempt to include the influence of the J11/5 and 3J-2S-1 resonances. When the drifting objects spand longer time than ~ 50 Myr in the zone of these two resonances, I eject them from the simulation. Obviously, I would be happy for any suggestions how to improve my model in this respect. If Bill B would send me his data file of Themis I might plot here a comparison with the observed Themis a-distribution.

From the data above, I may compute how the standard deviation of the semimajor-axes distribution changes in time; here is the corresponding plot. In black you see the modeled family, green are the data from Bill B for the current family (as I discussed above, the real sigma(a) in the largest bin - 15-20 km - may be lower according to data file of A Cellino). Notice the expected features of the solution. The smallest bodies move fast so that their sigma(a) increases significantly in the first phase of the simulation. However, these bodies have also the shortest lifetime, so that they become to disappear from the family first (also falling in the J2/1 resonance). This produces flatenning of the sigma(a) "dynamics" in the smallest size bin for epochs later than ~ 0.5 Gyr. Notice, that this is by chance maximum time of our numerical simulation above. Small bump in the sigma(a) of the smallest size bin at ~ 1.75 Gyr is caused by disruption of one of the family large members (~100 km size in my simulation), that was obviously located near the center of the family. Dynamics of sigma(a) at larger size-bins is slower, but more steady.

Finally, let us have a look at the size distribution of the family members. The following plot shows my initial distribution (green), final distribution affected uniquely by the collisional evolution of the family (blue), "modeled" distribution (red), where I have applied some guessed completeness-factors in different size-bins, and the observed (black) size distribution of Themis members (here I use data file by A Cellino). The difference between green and blue curves indicates that family looses small members faster than replenishes from break-ups of larger members, so that the steepness of the distribution get shallower at small sizes. These small sizes are then affected by observational effect of our incomplete knowledge of the family members (blue -> red curves). Here I have assumed the following ad hoc completeness factor, by which I have multiplied the modeled final distribution. I would be very happy, if somebody could tell me more reasonable numbers. At the moment I looked at Fig. 5 of Jedicke and Metcalfe (1998), where they compare Spacewatch debiased population with observations performed by Palomar Leiden Survey; well, you see that this is already a sort of hybrid information. The magnitude bins in this figure were recomputed to the size bins using an approximate formula D ~ 5900 10^{-0.2 H} (in km); obviously I took take from the 3.0-3.5 AU range. Obviously, we cannot get the really good fit of the observed size distribution (black curve), since we do not have the truly self-consistent collisional model. For instance, I would believe that the hump at ~ 80-90 km sizes might have to do with changing the regime by which the bodies "keep together" - gravity vs. molecular cohesion - and this rules their strengthness with respect to disruptions. Things like that should be modeled by a dedicated study. I am basically satisfied, that we could obtain something "reasonable" without putting in some really crazy factors. Suggestions?

My next plan is to see what happens if I force the spin axes of the asteroids to evolve ``faster'', something that YORP suggests. One possibility would be to constrain the factor cf1 to be small, while still allowing cf2 large. However, this would not respect the right size dependence of the characteristic timescale of the YORP effect, which is much steeper than ~ \sqrt{R}, notably more like ~ R^2. Now, I am trying to work slowly on YORP with one of my students and we have some new findings, if compared to Dave Rubincam's Icarus paper, that I want to use here. First, let me explain what YORP does and what are the caveats of current modeling.

YORP secularly (``meaning periods longer than the orbital period'') changes (i) rotation rate and (ii) obliquity of the spin axis. The effect is nil for spherical bodies, but is generic for irregular shapes. The pioneering paper here in Rubincam (2000). Now you understand, that an infinite variety of shapes exist for asteroids (meteoroids) smaller than ~km across. See here for fun several well-determined shapes (axes are in kilometers):


Castalia1998 KY26


These data were obtained using radar imaging by the Ostro's group at JPL. In our simulations we have used all these cases, plus a couple of others, but we have also generated ``synthetic'' objects by the Gaussian approach developed by Kari Muinonen (1996, 1998). Now one of the - expected - troubles with YORP is that there is no generic result of YORP: YORP evolution for Golevka is different from that of 1998 KY26 and this is still different from the YORP evolution of Kleopatra (you should understand ``the shape of ...'' since the linear dimensions of the object do not matter for the qualitative result). My impression is that this feeling is missing in Dave Rubincam's paper: e.g. rotation of any object of the Kleopatra shape will only be decelerated by YORP and never accerated, Golevka's obliquity will be asymptotically approaching some ~40 degrees, while obliquity of 1998 KY26 will approach 90 degrees, etc. There is, however, one major feature: if nothing than YORP is assumed for the spin axis evolution, the obliquity will be always assymptotically approaching some value (most often 0 degrees or 90 degrees, but any intermediate value is also possible) and at this terminate state of the obliquity evolution the rotation rate is in most cases (~ 95 percent of cases) decelerated. Hence the likely long-term YORP evolution is that the axis is slowly tilted to some asymptotic value; as soon as this is approached the rotation rate is decreased enough so that an impact will reorient the axis (and speed the rotation) and the cycle begins again. What matters mostly for our application, I believe, is the right estimate of the timescale on which this cycle occurs. With my another PhD student, David Capek, I have tried to estimate the statistically mean value of the YORP cycle duration (ie. starting from a random orientation and mean rotation period of 5hr till the YORP despuns the rotation to 50hr). We used a statistical sample of 1000 synthetic asteroids whose shapes were generated by the Gaussian random sphere technique of Muinonen (1996, 1998) and Muinonen and Lagerros (1998). For R=1 km Themis asteroid we got 60 Myr, rather short if compared with our 500 Myr simulation of Themis members with fixed orientation of the rotation axes. The factor 60 Myr is comparable to ~ 30 Myr we would got by appropriate scaling of the Dave Rubincam's (2000) result for Eros. This means that the YORP effect is slightly amplified for the Eros bean-like shape than in average for the Gaussian random sphere of Muinonen. For larger objects the 60-Myr factor scales as ~ R^2 (as I have mentioned above). So for R=3 km Themis asteroid the YORP timescale becomes ~ 540 Myr, larger than our simulation. Note, however, that there is an large uncertainty in the numerical values mentioned above, maybe up to factor 2. From what is summarized above, I propose the following approach (to start the thing): since the asymptotic obliquity is not generic (0 degrees with about the same likelihood as 90 degrees), just use the previous Poissonian statistical approach by keeping the spin axis fixed until a a condition of impact-reorientation is randomly satisfied. The reorientation timescale is, however, \tau_rot ~60 * R^2 Myr (with R in km). At small sizes this timescale is quite shorter that of the collisional reorientations; however, thanks to steep dependence on R the collisional reorientations begin to dominate for larger asteroids with size > 3 km (see above). We have thus started simulation of the Themis evolution with this zero-level modeling of YORP; only smallest 100 members are included in our simulation, since the larger objects should not presumably get any, or just one, reorientation within the 500 Myr period.

As a remark, I just note here a few caveats of the current YORP modeling (we again are working to remove them):

  1. the previous model neglects thermal relaxation (ie. neglects finite thermal conductivity of the asteroid surface) and this means that it is applicable to larger asteroids with regolith; meteoroids that likely have more conductive surface need a thermophysical model of YORP; we have developed a first approximation and showed that the YORP-cycle timescale is prolonged by a factor of 10-100 by the finite thermal condutivity;
  2. we always assumed rotation around the principal axis of rotation; it is likely that the impact-event at the end of the YORP cycle may also induce a non-principal-axis rotation state that will be eventually damped by the internal dissipation; here again estimations of the damping timescale are rather uncertain but generally they are quite shorter than the YORP-cycle timescale for km-sized bodies but may become quite longer than the YORP-cycle timescale for smaller objects.
  3. other torques that may on a long-term affect the spin axis are neglected; most importantly we should include the gravitational torque due to the solar gravity which alone would cause a reguliar precession of the axis around the orbital angular momentum (hence keeping the obliquity constant); typical precession periods range 0.3 Myr to tens (or hundreds) of Myr, depending on sphericity of the body, its rotation rate etc. However, the situation might be more complicated because the orbital plane evolves due to the planetary perturbations. What matters then is the spectrum of the nonsingular p and q variables (combining node and inclination); if some of the orbital evolution lines are close to the forced precession frequency (as it may happen) then resonances may occur. Still more important is the fact that the p and q spectrum may contain clusters of non-isolated lines; if the precession frequency get into interaction with such a cluster, chaoticity may be onset and obliquity could undergo large variations (see Laskar and Robutel 1993, or some earlier works of Ward on Mars obliquity). I have been working on this topic over last few days and got some interesting results that you may find here. At the end, my conclusion is that these effects might not be fundamental (at the zero level of approximation) for our work and the YORP timescale is still the most relevant in our simulations.

Anyway here, is what happens if I include in my Monte Carlo code the above mentioned, somewhat extreme (and obviously simplified), modeling of YORP-induced effects on the spin axes (random reorientations after the characteristic timescale \tau_rot ~60 * R^2 Myr.

You see, that the result is not ``that excellent as before'', since the frequent reorientations of small-size asteroids prevent them to diffuse due to Yarkovsky orbital effect very far from their initial positions (hence the ``bell-shaped'' distribution of the final semimajor axes in the 1-5 km size bin). The sigma(a) increases more slowly than before. Obviously, the real bonus of this simulation (apart from the fact that it should be ``more realistic'' since it includes anyway unavoidable YORP effect) should appear in the eccentricity and inclination distributions that are not modelled here. Namely, random wandering in the semimajor axis ``up and down'' due to axis reorientations should help to disperse the family in these other proper elements thanks to interactions with miniresonances. I am still running some more tests to see if could get a better result, so that we shall see more next week.



  1. M. Broz, 1999, Long-term dynamics of small bodies in the solar system with gravitational and thermal perturbations (Thesis, Charles University, Prague).
  2. R.R. Codeiro, R.S. Gomes and R. Vieira Martins, 1997, Celest. Mech. Dyn. Astr. 65, 407.
  3. P. Farinella and D. Vokrouhlicky, 1999, Science 283, 1507.
  4. P. Farinella, D. Vokrouhlicky and W.K. Hartmann, 1998, Icarus 132, 378.
  5. R. Jedicke and T.S. Metcalfe, 1998, Icarus 131, 245.
  6. L. Kresak, 1979, Bull. Astron. Inst. Czechosl. 28, 65.
  7. J. Laskar and P. Robutel, 1993, Nature 361, 608.
  8. J. Laskar, F. Joutel and F. Boudin, 1993, Astron. Astrophys. 270, 522.
  9. R. Malhotra, 1994, Celest. Mech. Dyn. Astr. 60, 373.
  10. F. Marzari, P. Farinella and D. Davis, 1999, Icarus ???
  11. S. Mikkola, 1998, Celest. Mech. Dyn. Astr. 68, 249.
  12. A. Milani and Z. Knezevic, 1992, Icarus 98, 211
  13. T.G. Muller and J.S.V. Lagerros, 1999, Bull. Am. Astron. Soc. 31, 1075.
  14. K. Muinonen, 1996, Earth, Moon and Planets 72, 339.
  15. K. Muinonen, 1998, Astron. Astrophys. 332, 1087.
  16. K. Muinonen and J.S.V. Lagerros, 1998, Astron. Astrophys. 333, 753.
  17. O. Neron de Surgy and J. Laskar, 1997, Astron. Astrophys. 318, 975.
  18. D.P. Rubincam, 2000, Icarus, 148, 2.
  19. J.N. Spitale and R. Greenberg, 2001, Icarus 149, 222.
  20. D. Vokrouhlicky, 1998, Astron. Astrophys. 335, 1093.
  21. D. Vokrouhlicky and P. Farinella, 1999, Astron. J. 118, 3049.
  22. D. Vokrouhlicky and P. Farinella, 2000, Nature 407, 606.

David Vokrouhlicky, and Miroslav Broz,, Apr 18th 2001