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

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)
Dora
Eos
Themis

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

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
Dora
Eos
Themis
Eunomia (75 Myr)

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

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-g₆-s₆ in Eos region. TPs from the Eos region might be captured in g+s-g₆-s₆ (or g+s-g₅-s₇, 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 (a, e) 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-g₆-s₆ (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-g₆-s₆ (t), g+s-g₅-s₇ (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-g₆-s₆ (t), g+s-g₅-s₇ (t); this is again an interesting case of a TP, that has been ejected from the g+s-g₆-s₆ 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-s₆-g₅+g₆ 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-g₆-s₆ 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-g₆-s₆ resonance (ie. 31 %). We can say that because we individually looked at the evolution of the g+s-g₆-s₆ 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-g₆-s₆ 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 z₁ 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 z₁ region; special emphasis is given to selection of the objects to be observed).

Two more integrations were performed to check the behavior of z₁ 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 z₁ < 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 < z₁ < 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)

81 Maria members (maria-a1 run)

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

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

See also SWRI web for new results on Flora, Eunomia, Maria and Koronis families: http://www.boulder.swri.edu/~davidn/yarko/yarko.html, http://www.boulder.swri.edu/~bottke/Animation/.

Motivation and a brief description of the model

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).

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.

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).

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):

Golevka	Geographos

Castalia	1998 KY26

Kleopatra

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.

UNDER CONSTRUCTION (SEE LATER).

References

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.